Matriz exponencial de una matriz asimétrica real con Fortran 95 y LAPACK

Hace poco hice una pregunta en la misma línea para las matrices ermitas sesgadas. Inspirado por el éxito de esa pregunta, y después de golpearme la cabeza contra la pared durante un par de horas, estoy mirando la matriz exponencial de las matrices asimétricas reales. La ruta para encontrar los valores propios y los vectores propios parece bastante complicada, y me temo que me he perdido.

Antecedentes: Hace algún tiempo hice esta pregunta sobre la física teórica SE. El resultado me permite formular ecuaciones maestras como matrices asimétricas reales. En el caso independiente del tiempo, la ecuación maestra se resuelve exponiendo esta matriz. En el caso dependiente del tiempo, requerirá integración. Solo me preocupa la independencia temporal en este momento.

Al mirar las diversas subrutinas creo que debería estar llamando a ( ? Gehrd , ? Orghr , ? Hseqr ...) no está claro si sería más sencillo para emitir la matriz a partir real*8de complex*16y proceder con las versiones dobles complejas de estas rutinas, o quédate real*8y toma el éxito de duplicar el número de mis matrices y hacer una matriz compleja de ellas más tarde.

Entonces, ¿a qué rutinas debería llamar (y en qué orden), y debería usar las versiones dobles reales o las versiones dobles complejas? A continuación hay un intento de hacer esto con versiones dobles reales. Me he quedado atascado encontrando los valores propios y los vectores propios 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

linear-algebra fortran lapack qubyte
fuente

Respuestas:

Primero pensaría mucho si la matriz es completamente arbitraria o no: ¿hay alguna transformación que la convierta en Hermitiana? ¿La física garantiza que la matriz debe ser diagonalizable (con una matriz de vectores propios razonablemente condicionada)?

Si resulta que realmente no hay ninguna simetría para explotar, entonces debe comenzar leyendo Diecinueve formas dudosas de calcular la matriz exponencial , que es la referencia estándar (y está escrita por el autor de MATLAB y el coautor de G & vL) .

Jack Poulson
fuente

2 \times 2

$2\times2$

4 \times 4

$4\times4$

Me gusta esta respuesta el caso asimétrico posee suficientes dificultades que vale la pena considerar si puede haber una formulación de su problema que conduzca a matrices simétricas en lugar de otras asimétricas.

@ MarkS.Everitt: Pareces estar casi allí ... ¿qué tan grandes son las matrices? ~ 36 x 36 otra vez?

Jack Poulson

16 \times 16

$16\times16$

36 \times 36

$36\times36$

@ MarkS.Everitt: Entonces, su problema ahora es efectivamente cómo exponer matrices 4x4. Esto es lo suficientemente pequeño para que el análisis asintótico sea irrelevante, por lo que la respuesta dependerá completamente de los valores. Realmente no puedo decir nada a menos que traduzcas tu publicación física vinculada en álgebra lineal (¿qué es un superoperador?!?).

Jack Poulson

Para construir sobre lo que Jack ha dicho, el enfoque estándar que parece usarse en software (como EXPOKIT, mencionado en su pregunta anterior) es escalar y cuadrar seguido de aproximación de Padé (Métodos 2 y 3) o métodos de subespacio de Krylov (Método 20) En particular, si está buscando integradores exponenciales, querrá considerar los métodos del subespacio de Krylov y mirar los documentos sobre integradores exponenciales (se mencionan algunas referencias junto con el Método 20 en el documento Moler & van Loan).

Si está empeñado en usar vectores propios, considere usar sistemas triangulares de vectores propios (Método 15); Como su matriz puede no ser desglosable, este enfoque podría no ser el mejor, pero es mejor que tratar de calcular los vectores y valores propios directamente (es decir, Método 14).

La reducción a la forma de Hessenberg no es una buena idea (Método 13).

No es evidente para mí si estaría mejor servido con aritmética real o compleja, ya que la aritmética compleja de Fortran es rápida, pero podría desbordarse / subdesbordarse (consulte "¿Cuánto mejor son realmente los compiladores de Fortran?" ).

Puede ignorar los Métodos 5-7 (los métodos basados en solucionadores de EDO son ineficientes), los Métodos 8-13 (caros), el Método 14 (calcular vectores propios de matrices grandes es difícil sin una estructura especial y propenso a errores numéricos en casos mal condicionados) y Método 16 (calcular la descomposición de Jordan de una matriz es numéricamente inestable). Los métodos 17-19 son más difíciles de implementar; en particular, los métodos 17 y 18 requerirían más lectura. El método 1 es una opción alternativa para escalar y cuadrar si las aproximaciones de Padé no funcionan bien.

$B_{j}$

B_{j} = γ_{j} I + E_{j},

$B_j = \gamma_j I + E_j,$

$\gamma_{j}$ $j$ $E_{j}$

Geoff Oxberry
fuente

O (n^{2})

$O(n^2)$

O (n^{3})

$O(n^3)$

Sin duda saben lo que están haciendo; No me preocupa la implementación de LAPACK. Estoy más sorprendido por el comportamiento del compilador Fortran.

Geoff Oxberry

Sí, el compilador podría ser más problemático que un LAPACK bien escrito. Puede ser desconcertante descubrir que su programa estaba fallando solo porque las implementaciones para el valor absoluto y la división utilizada por el compilador fueron fallidas ...

-1

Tengo una subrutina Fortran simple que calcula el exponente de una matriz arbitraria. Lo comprobé con el comando de Matlab y está bien. Se basa en escalar y cuadrar. Lo escribí hace unos años.

Me hubiera gustado encontrar otra subrutina, como las que descargo de gams.nist.gov. Pero aún no hay suerte.

freejack
fuente