Como o LAPACK resolve sistemas tridiagonais e por quê?

9

No meu projeto, tenho que resolver duas matrizes tridiagonais a cada passo do tempo, por isso é crucial ter um bom solucionador para elas. Fiz minha própria implementação, exatamente da maneira clássica de fazê-lo, descrita na Wikipedia. Então tentei usar o Lapack e, para minha surpresa, foi mais lento!

Agora, dentro de Lapack, parece que ele resolve a fatoração da LU e me pergunto por que, não é mais complexo do que poderia ser?

Além disso, encontrei um algoritmo no livro "Receitas Numéricas" do nr.com que divide recursivamente o sistema em problemas tridiagonais menores. Parecia promissor. Existem outras guloseimas por aí?

Atualização: o tamanho do problema é de cerca de 1000 x 1000. Eu usei o GotoBLAS, ele também fornece uma biblioteca Lapack 3.1.1. O problema não é simétrico. Usei a rotina Lapack para matrizes tridiagonais gerais.

banded-matrix lapack tiam
fonte
2
Você precisará indicar quais rotinas LAPACK foram usadas para isso. Observe que o dgtsv executa rotação parcial, mas seu código pode não fazer isso. Indique também com qual implementação do LAPACK você testou e com que tamanho de problema você comparou. Além disso, seu problema é simétrico positivo definido?
precisa
Eu adicionei algumas informações na formulação da pergunta.
tiam
Seu aplicativo tem algo a ver com métodos de volume finito?
Inquérito
São diferenças finitas, mas nessa perspectiva é mais ou menos o mesmo, eu acho.
tiam

Respostas:

6

Você está usando uma implementação de referência que faz rotação parcial. As soluções tridiagonais fazem muito pouco trabalho e não entram no BLAS. Provavelmente é mais lento que o seu código, porque faz rotação parcial. O código fonte do dgtsv é direto.

Se você resolver com a mesma matriz várias vezes, convém armazenar os fatores usando dgttrf e dgttrs . É possível que as implementações em um LAPACK otimizado, como MKL, ACML ou ESSL, tenham melhor desempenho.

Jed Brown
fonte
Estou um pouco curioso. O Gaussian Elim com PP funcionaria para todas as matrizes, incluindo TriDiagonal. No CFD, usamos um método especial para casos de FVM 1D chamado TDMA . O que você acha que seria mais rápido para o caso que ele está discutindo? Embora eu não tenha certeza absoluta de que suas matrizes são diagonalmente dominantes.
Inquérito
O TDMA é o que eu implementei no meu código. A questão é por que o Lapack super-rápido usaria o procedimento de rotação parcial em uma matriz específica, que é resolvida mais rapidamente por um método tão fácil como o TDMA.
tiam
É exatamente o mesmo algoritmo (eliminação gaussiana especializada para uma matriz tridiagonal), mas sua implementação não faz pivotamento parcial, portanto, pode ser numericamente instável. Essa rotação não é gratuita e você está comparando com a implementação de referência. A implementação de referência não é otimizada para desempenho e a rotação parcial não é livre.
precisa
Entendo o que você quer dizer, aproveito meu conhecimento sobre os sistemas que estou resolvendo. Outras implementações do LAPACK aumentam o desempenho devido à adaptação a arquitetura específica ou vão além disso?
tiam