Как я могу вычислить логарифм матрицы, используя Fortran 77? - PullRequest
Новая
       67

Как я могу вычислить логарифм матрицы, используя Fortran 77?

0 голосов
/ 07 июня 2018

Мне нужно вычислить логарифм матрицы 3х3 (A) в программе на Фортране 77.

В Python я использовал

scipy.linalg.logm(A)

для этой работы, но я не нашел решения этой задачи в Фортране.

До сих пор я выяснил, что эта операция невозможна со встроенными функциями Fortran 77. Также я искал в документации библиотеки Intel Math Kernel, но не нашел подходящей подпрограммы.Я нашел подпрограмму F01EJ в библиотеке NAG Fortran, но, к сожалению, у меня нет доступа к этой коммерческой библиотеке.Я знаком с работами Хайама и др., Но я бы не хотел реализовывать алгоритмы самостоятельно, так как считаю, что это должно быть решенной проблемой.

Может ли кто-нибудь дать мне подсказку для подпрограммы, которая вычисляет матричный логарифм?

Решение:

Я реализовал подпрограмму для логарифмаматрицы 3х3 по способу, предложенному @kvantour, и у меня все получилось.Возможно (не очень красивый) быстрый и грязный код подпрограммы полезен для других с такой же проблемой: 

  subroutine calclogM(M,logM)

  implicit none

  double precision,dimension(3,3)::M,logM,VL,VR,logMapo,VRinv
  integer::n,INFO,LWORK,I,J

  double precision,dimension(3)::WR,WI,logWR,ipiv
  double precision,dimension(24)::WORK

  n=3
  LWORK=24


  call DGEEV( 'N', 'V', n, M, n, WR, WI, VL, n, VR,
 1                  n, WORK, LWORK, INFO )

C Check if all eigenvalues are greater than zero

  if (WR(1) .le. 0.D0) then
    write(*,*) 'Unable to compute matrix logarithm!'
    GOTO 111
  end if

  if (WR(2) .le. 0.D0) then
    write(*,*) 'Unable to compute matrix logarithm!'
    GOTO 111
  end if

  if (WR(3) .le. 0.D0) then
    write(*,*) 'Unable to compute matrix logarithm!'
    GOTO 111
  end if

  DO I = 1, 3
    DO J = 1, 3
        logMapo(I,J) = 0.D0
    END DO
  END DO

C Then Mapo will be a diagonal matrix whose diagonal elements 
C are eigenvalues of M. Replace each diagonal element of Mapo by its 
C (natural) logarithm in order to obtain logMapo.

  DO I = 1, 3
    LogMapo(I,I)=log(WR(I))
  END DO


C Calculate inverse of V with LU Factorisation
C Copy VR to VRinv
  DO I = 1, 3
    DO J = 1, 3
        VRinv(I,J) = VR(I,J)
    END DO
  END DO

  call dgetrf( n, n, VRinv, n, ipiv, info )
  write(*,*) 'INFO',INFO

  call dgetri( n, VRinv, n, ipiv, WORK, LWORK, INFO )
  write(*,*) 'INFO',INFO

C Build the logM Matrix
  logM = matmul(matmul(VR,logMapo),VRinv)      

111  end subroutine calclogM

1 Ответ

0 голосов
/ 07 июня 2018

Чтобы вычислить логарифм матрицы, вам необходимо вычислить собственные значения λ {1,2,3} вашей исходной матрицы.Если все ваши собственные значения больше, чем ZERO, то вы можете вычислить логарифм.

Итак, я бы предпринял следующие шаги:

  1. Вычислим собственные значения λ {1,2,3} и соответствующие ему собственные векторы v {1,2,3} .
  2. Проверьте, все ли собственные значения больше, чем ZERO.
  3. Если условие (2) выполнено, вычислите Log из собственных значений.
  4. Используя логарифмические значения, восстановите матрицу с помощьюсобственные векторы.

Простые методы вычисления собственных векторов и собственных значений матрицы можно найти в числовых рецептах (см. старый раздел книги).Более продвинутые методы могут быть основаны на blas и lapack - вас интересует процедура dgeev .

Лучший способ - попытаться реализовать собственную декомпозицию самостоятельно какэто всего лишь матрица 3х3.

...