exponential of a matrix in Fortran using Taylor Series expansion

I am having a matrix of n x n and I want to calculate exponential_of_matrix(matrix_name) in Fortran. Is there anyone who knows to calculate the exponential of a matrix using Taylor Series Expansion?

Taylor Series Expansion of e^x = 1 + x + x^2/2! + x^3/3! + x^4/4! + .....

I tried to write a matrix multiplication subroutine which may be useful during writing matrix exponential subroutine.

! mat_mul: matrix_a(n,m) matrix_b(m,l) product(n,l)

subroutine mat_mul(matrix_a,matrix_b, product,n,m,l)
    real*8 matrix_a(n,m), matrix_b(m,l), product(n,l)
    integer i,j,k
    do i=1,n
        do j=1,l
            product(i,j) = 0
            do k=1,m
                product(i,j) = product(i,j) + matrix_a(i,k) * matrix_b(k,j)
            end do
        end do
    end do
end subroutine mat_mul

Solution 1:

First, if you want to multiply matrices you should really use matmul¹ rather than doing it manually, so

do i=1,n
  do j=1,l
    product(i,j) = 0
    do k=1,m
      product(i,j) = product(i,j) + matrix_a(i,k) * matrix_b(k,j)
    end do
  end do
end do

becomes

product = matmul(matrix_a, matrix_b)

This code is clearer, and significantly faster.

Second, are you sure you want to use Taylor series? If the matrix is diagonalisable then matrix exponentiation is typically calculated via diagonalisation:

1. Diagonalise the matrix to get its eigenvectors v_i and eigenvalues e_i.
2. Taking the exponential of the eigenvalues e^e_i.
3. Construct the matrix with eigenvectors v_i and eigenvalues e^e_i.

¹ or even a BLAS/LAPACK distribution if you need more performance.