计算矩阵指数,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
还有就是阶数过大,在求逆矩阵那一块如有必要,也可进行适当修改