如何使用 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. 计算特征值λ_{1,2,3}及其对应的特征向量v_{{1,2 ,3}}.
2. 检查是否所有特征值都大于ZERO。
3. 如果满足条件 (2)，则计算特征值的Log。
4. 使用对数值，通过特征向量重构矩阵。

计算矩阵的特征向量和特征值的简单方法可以在numerical recipes 中找到（查看旧书部分）。更高级的方法可以基于blas and lapack——有兴趣的套路dgeev。

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