本篇博文简要介绍使用MKL函数库计算多个右端项的线性代数方程组
代码如下:
program MKL_gesv
use lapack95
implicit none
integer :: i
integer, parameter :: n = 3, m = 4
real*8 :: a(n,n), aa(n,n), b(n,m), bb(n,m)
call random_seed()
call random_number(aa)
a = aa
call random_number(bb)
b = bb
call gesv(a, b) !// x: -> b
write(*,*) 'checking...'
do i = 1, m
write(*,*) maxval(abs(matmul(aa,b(:,i))-bb(:,i)))
end do
write(*,*) 'the result is correct!'
end program MKL_gesv
结果如下:
checking...
0.000000000000000E+000
1.110223024625157E-016
1.110223024625157E-016
2.775557561562891E-017
the result is correct!
从结果可以看出,我们的结果是正确的。
代码中要说明的一点是:call gesv之后
原系数矩阵由A =P*L*U分解得到的因子L和U覆盖;
而b则被求解出来的x所覆盖。