如何使用 Fortran 77 计算矩阵的对数？

Question

我必须在 Fortran 77 程序中计算 3x3 矩阵 (A) 的对数。

在Python中我使用了

scipy.linalg.logm(A)

对于这项工作，但我没有在 Fortran 中找到该任务的解决方案。

到目前为止，我发现Fortran 77的内置函数无法实现此操作。另外我在Intel Math Kernel Library的文档中进行了搜索，但没有找到合适的子程序。我在 NAG Fortran 库中找到了子程序 F01EJ，但不幸的是我无法访问这个商业库。我知道 Higham 等人的论文，但我想避免自己实现这些算法，因为我认为这一定是一个已解决的问题。

谁能给我一个计算矩阵对数的子程序的提示？

解法：

我已经根据@kvantour 提出的方法为 3x3 矩阵的对数实现了一个子程序，对我来说效果很好。也许子例程的（不是很漂亮）快速而肮脏的代码对其他有同样问题的人有用：

  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

Answer 1

为了计算矩阵的对数，您需要计算特征值 λ_{1,2,3}你的原始矩阵。如果你所有的特征值都大于 ZERO，那么你可以计算对数。

所以我要采取的步骤是：

计算特征值 λ_{1,2,3} 及其对应的特征向量 v_{1,2,3}.
检查是否所有特征值都大于 ZERO。
如果满足条件(2)，则计算特征值的Log。
使用对数值，用特征向量重建矩阵。

可以在 numerical recipes (check the old book section). More advanced methods can be based on blas and lapack - you are interested in the routine dgeev 中找到计算矩阵的特征向量和特征值的简单方法。

最好的方法是尝试自己实现特征分解，因为它只是一个 3x3 矩阵。

如何使用 Fortran 77 计算矩阵的对数？

How can I calculate the logarithm of a matrix using Fortran 77?

fortran

logarithm

matrix

linear-algebra

intel-mkl