fortran使用MKL函数库计算矩阵指数函数

计算矩阵指数，matlab直接可以用expm，但是fortran笔者了解到MKL好像还没有子程序直接可以使用
闲暇之余，给出一个求矩阵指数的fortran代码仅供参考
代码如下，如有必要，读者为了以后的方便也可将下面代码直接封装
program test_expm
        use lapack95
        implicit none
        integer, parameter :: n = 3
        integer :: i, info, ipiv(n)
        real :: a(n,n), wr(n), wi(n), vr(n,n), vl(n,n)
        real :: temp(n,n)

        a = reshape([1, 1, 0, 0, 0, 2, 0, 0, -1],shape(a))
        a = transpose(a)

        call geev(a, wr, wi, vl, vr, info)
        
        write(*,*) '特征值为:'
        write(*,'(*(f10.4))') wr

        write(*,*) '特征向量为:'
        do i = 1, n
                write(*,'(*(f10.4))') vr(i,:)
        end do
        
        if (info == 0) then
                write(*,*) '特征值与特征向量计算成功!'
        else
                write(*,*) '特征值与特征向量计算失败!'
        end if
        
        write(*,*)
        write(*,*) '开始计算矩阵指数......'
        temp = 0.
        forall( i = 1:n ) temp(i,i) = exp(wr(i))

        vr = vl
        call getrf(vr, ipiv)
        call getri(vr, ipiv)  !// 求特征向量的逆矩阵

        !// 计算矩阵指数
        a = matmul(matmul(vl,temp),vr)
        write(*,*) '矩阵A的矩阵指数为:'
        do i = 1, n
                write(*,'(*(f10.4))') a(:,i)
        end do
        
end program test_expm

执行结果如下
 特征值为:
    1.0000    0.0000   -1.0000
 特征向量为:
    1.0000   -0.7071    0.4082
    0.0000    0.7071   -0.8165
    0.0000    0.0000    0.4082
 特征值与特征向量计算成功!

 开始计算矩阵指数......
 矩阵A的矩阵指数为:
    2.7183    1.7183    1.0862
    0.0000    1.0000    1.2642
    0.0000    0.0000    0.3679

如有必要，读者可与matlab结果检验一下。

要注意的地方：如果所求矩阵阶数较大，比如几千阶，程序中的matmul函数建议更换为MKl中的gemm计算
具体可参考https://blog.csdn.net/chd_lkl/article/details/94437455

还有就是阶数过大，在求逆矩阵那一块如有必要，也可进行适当修改