2

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
Ian Bush
  • 6,996
  • 1
  • 21
  • 27
Swapnal Shahil
  • 84
  • 1
  • 7

2 Answers2

1

First, if you want to multiply matrices you should really use matmul1 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.

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

veryreverie
  • 2,871
  • 2
  • 13
  • 26
0

I have written this subroutine using Taylor Series Expansion. Instead of mat_mul subroutine, its better to use MATMUL(matrix_A,matrix_B).

subroutine exponent_matrix(a,expo_mat,n)
    real a(n,n), prod(n,n), expo_mat(n,n), term(n,n)
    do i=1,n
        do j=1,n
            if(i .eq. j) then
                expo_mat(i,j) =1
                term(i,j) = 1
            else
                expo_mat(i,j) = 0
                term(i,j) = 0
            end if
        end do
    end do
    
    do i=1,5000
        call mat_mul(term,a,prod,n,n,n)
        term = prod/i
        expo_mat = expo_mat + term
    end do
end subroutine exponent_matrix
Swapnal Shahil
  • 84
  • 1
  • 7
  • 1
    Note real*8 is not and has never been part of the standard Fortran language, and as an extension may not be available on all platforms. Please learn about kind values - see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Ian Bush Jan 16 '22 at 10:48