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

Question

I have to compute the logarithm of a 3x3 matrix (A) in a Fortran 77 program.

In Python I used

scipy.linalg.logm(A)

for this job but I don´t find a solution for this task in Fortran.

So far I figured out that this operation is not possible with the built-in functions of Fortran 77. Also I searched in the documentation of the Intel Math Kernel Library, but haven´t found a suitable subroutine. I found the subroutine F01EJ in the NAG Fortran Library, but unfortunatly I have no access to this commercial library. I know the papers of Higham et al., but I would like to avoid to implement the algorithms by myself since I think that this must be a solved problem.

Can anyone provide me with a hint to a subroutine that calculates the matrix logarithm?

Solution:

I have implemented a subroutine for the logarithm of a 3x3 matrix according to the way proposed by @kvantour and it worked out fine for me. Maybe the (not very pretty) quick and dirty code of the subroutine is useful for others with the same problem:

  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

In order to compute the logarithm of a matrix, you need to compute the eigenvalues λ_{1,2,3} of your original matrix. If all your eigenvalues are bigger then ZERO, then you can compute the logarithm.

So the steps I would take are :

1. Compute the eigenvalues λ_{1,2,3} and its corresponding eigenvectors v_{1,2,3}.
2. Check if all eigenvalues are bigger then ZERO.
3. If condition (2) is satisfied, compute the Log of the eigenvalues.
4. Using the logarithmic values, reconstruct the matrix with the eigenvectors.

Simple methods to compute the eigenvector and eigenvalues of a matrix can be found in numerical recipes (check the old book section). More advanced methods can be based on blas and lapack - you are interested in the routine dgeev.

The best way would be to attempt to implement the eigendecomposition yourself as it is just a 3x3 matrix.