Matriz exponencial de uma matriz assimétrica real com Fortran 95 e LAPACK

10

Recentemente, fiz uma pergunta na mesma linha para matrizes enviesadas-hermitianas. Inspirado pelo sucesso dessa pergunta, e depois de bater minha cabeça contra uma parede por algumas horas, estou olhando para o exponencial da matriz de matrizes assimétricas reais. O caminho para encontrar os autovalores e autovetores parece bastante complicado, e receio ter me perdido.

Antecedentes: Há algum tempo, fiz essa pergunta sobre a física teórica SE. O resultado me permite formular equações mestres como matrizes assimétricas reais. No caso independente do tempo, a equação principal é resolvida exponenciando essa matriz. No caso dependente do tempo, será necessária integração. Só estou preocupado com a independência do tempo no momento.

Ao olhar para as várias sub-rotinas eu acho que deveria ser chamado ( ? Gehrd , ? Orghr , ? Hseqr ...) não está claro se ele seria mais simples para lançar a matriz de real*8a complex*16e prossiga com as versões duplas complexos destas rotinas, ou fique com real*8e tente dobrar o número de minhas matrizes e criar uma matriz complexa delas mais tarde.

Então, quais rotinas devo chamar (e em que ordem) e devo usar as versões duplas reais ou as versões duplas complexas? Abaixo está uma tentativa de fazer isso com versões duplas reais. Eu fiquei preso encontrando os autovalores e autovetores de L*t.

function time_indep_master(s,L,t)
  ! s is the length of a side of L, which is square.
  ! L is a real*8, asymmetric square matrix.
  ! t is a real*8 value corresponding to time.
  ! This function (will) compute expm(L*t).

  integer, intent(in)    :: s
  real*8,  intent(in)    :: L(s,s), t
  real*8                 :: tau(s-1), work(s), wr(s), wi(s), vl
  real*8, dimension(s,s) :: time_indep_master, A, H, vr
  integer                :: info, m, ifaill(2*s), ifailr(2*s)
  logical                :: sel(s)

  A = L*t
  sel = .true.

  call dgehrd(s,1,s,A,s,tau,work,s,info)
  H = A
  call dorghr(s,1,s,A,s,tau,work,s,info)
  call dhseqr('e','v',s,1,s,H,s,wr,wi,A,s,work,s,info)
  call dhsein('r','q','n',sel,H,s,wr,wi,vl,1,vr,s,2*s,m,work,ifaill,ifailr,info)

  ! Confused now...

end function
qubyte
fonte

Respostas:

8

Gostaria de pensar seriamente se a matriz é realmente completamente arbitrária: existe alguma transformação que a tornaria hermitiana? A física garante que a matriz seja diagonalizável (com uma matriz autovetor razoavelmente condicionada)?

Se realmente não houver simetria a ser explorada, comece lendo Dezenove maneiras duvidosas de calcular o exponencial da matriz , que é a referência padrão (e é escrita pelo autor do MATLAB e co-autor do G & vL) .

Jack Poulson
fonte
11
O melhor que posso fazer é transformá-lo em uma matriz diagonal de blocos com blocos assimétricos. Isso por si só é muito interessante. A maioria desses blocos é , e eu posso resolvê-los analiticamente. Existe um bloco restante sem simetrias para explorar. 2×24×4
qubyte
11
Eu gosto desta resposta; o caso assimétrico possui armadilhas suficientes que vale a pena considerar se pode haver uma formulação do seu problema que leve a matrizes simétricas em vez de matrizes assimétricas.
JM
@ MarkS.Everitt: Você parece estar quase lá ... qual o tamanho das matrizes? ~ 36 x 36 de novo?
21412 Jack Poulson
Nesse caso, , mas existe a possibilidade de subir para . 16×1636×36
qubyte
2
@ MarkS.Everitt: Então, o seu problema agora é exatamente como exponenciar matrizes 4x4. Isso é suficientemente pequeno para que a análise assintótica seja irrelevante; portanto, a resposta dependerá inteiramente dos valores. Eu realmente não posso dizer mais, a menos que você traduza sua publicação de física vinculada em álgebra linear (o que é um superoperador?!?).
21412 Jack Poulson
7

Para desenvolver o que Jack disse, a abordagem padrão que parece ser usada em software (como EXPOKIT, mencionado em sua pergunta anterior) é escalar e esquadrar, seguida pela aproximação de Padé (métodos 2 e 3) ou pelos métodos do subespaço de Krylov (método 20) Em particular, se você estiver olhando para integradores exponenciais, considere os métodos do subespaço de Krylov e veja artigos sobre integradores exponenciais (algumas referências são mencionadas juntamente com o Método 20 no artigo de Moler e van Loan).

Se você está decidido a usar vetores próprios, considere o uso de sistemas triangulares de vetores próprios (método 15); como sua matriz pode ser não-diagonalizável, essa abordagem pode não ser a melhor opção, mas é melhor do que tentar calcular os vetores próprios e os valores próprios diretamente (ou seja, método 14).

Reduzir a forma de Hessenberg não é uma boa ideia (método 13).

Não está claro para mim se você seria melhor atendido com aritmética real ou complexa, uma vez que a aritmética complexa de Fortran é rápida, mas pode estourar / estourar (consulte "Qual é realmente melhor os compiladores Fortran?" ).

Você pode ignorar com segurança os Métodos 5-7 (os métodos baseados no solucionador de EDO são ineficientes), os Métodos 8 a 13 (caro), o Método 14 (o cálculo de vetores próprios de matrizes grandes é difícil, sem estrutura especial e propenso a erros numéricos em casos mal condicionados). e Método 16 (calcular a decomposição de Jordan de uma matriz é numericamente instável). Os métodos 17-19 são mais difíceis de implementar; em particular, os métodos 17 e 18 exigiriam mais leitura. O método 1 é uma opção de retorno para escalar e esquadrar se as aproximações de Padé não funcionarem bem.

Editar # 1 : Com base nos comentários em resposta à resposta de Jack, a diagonalização de blocos parece ser uma opção; nesse caso, algo como o Método 18 (diagonalização triangular de blocos) é um método muito bom de usar. Hesitei em recomendá-lo no início porque sua pergunta não especificou essa estrutura, mas se você tiver uma transformação que bloqueie a diagonalização de sua matriz, isso tira a maior parte da complexidade da abordagem. Você só quer ter certeza de usar o truque de GW Stewart para decompor cada matriz diagonal de bloco emBj

Bj=γjI+Ej,

onde é a média dos valores próprios da matriz diagonal do ésimo bloco. Essa decomposição tornará quase nilpotente, o que aumentará a precisão dos seus cálculos exponenciais da matriz. Esse truque é discutido na página 26 da versão do artigo de Moler & van Loan ao qual Jack vinculou.γjjEj

Geoff Oxberry
fonte
11
Votado, mas o pessoal que implementou o LAPACK não é ingênuo em relação à aritmética complexa, especialmente considerando o tempo todo que Kahan dedicou a analisá-la. No entanto, ele ainda deve computar na aritmética real se ela salvar flops; o custo da conversão posteriormente é apenas versus toda a álgebra linear densa de . O(n2)O(n3)
Jack Poulson
Sem dúvida, eles sabem o que estão fazendo; Não estou preocupado com a implementação do LAPACK. Estou mais surpreso com o comportamento do compilador Fortran.
Geoff Oxberry
2
Sim, o compilador pode ser mais um problema do que o LAPACK bem escrito. Ele pode ser desconcertante para descobrir que seu programa só foi falhando porque as implementações de valor absoluto e divisão utilizada pelo compilador foram mal feita ...
JM
-1

Eu tenho uma sub-rotina Fortran simples que calcula o expoente de uma matriz arbitrária. Eu verifiquei o comando Matlab e está tudo bem. Baseia-se em escala e quadratura. Eu escrevi alguns anos atrás.

Eu gostaria de definir outra sub-rotina, como as que eu baixo de gams.nist.gov. Mas sem sorte ainda.

freejack
fonte