我在 Fortran90 中有一段代码，其中我必须解决非线性（使用 Newton-Raphson 方法，为此我必须反转雅可比矩阵）和线性方程组。当我说非常小的时候，我的意思是两个操作的 n 个未知数，使用 。不幸的是， 不是先验的。你认为最快的选择是什么？的情况编写明确的公式， 使用其他方法 （例如英特尔 MKL 库的某些函数）。这是明智的还是我应该为 写逆矩阵的显式公式？
该代码将在有限元方法分析中被调用很多次，因此我一直在寻找最快的解决方案。
在 Fortran 中有效地反转微小矩阵
计算科学
有限元
矩阵
正则
2021-12-19 23:13:43
2个回答
这是我们在 QMC=Chem 中用于小矩阵的例程。QMC=Chem 是在 GPL 许可下，所以如果你的代码也在 GPL 下，你可以使用它们。
编辑：要获得相反的结果，您需要将结果乘以1.d0/det_l.
subroutine invert2(a,LDA,na,det_l)
implicit none
double precision :: a (LDA,na)
integer :: LDA
integer :: na
double precision :: det_l
double precision :: b(2,2)
double precision :: f
b(1,1) = a(1,1)
b(2,1) = a(2,1)
b(1,2) = a(1,2)
b(2,2) = a(2,2)
det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1)
a(1,1) = b(2,2)
a(2,1) = -b(2,1)
a(1,2) = -b(1,2)
a(2,2) = b(1,1)
end
subroutine invert3(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
integer, intent(in) :: LDA
integer, intent(in) :: na
double precision, intent(inout) :: det_l
double precision :: b(4,3)
integer :: i
double precision :: f
det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
do i=1,4
b(i,1) = a(i,1)
b(i,2) = a(i,2)
b(i,3) = a(i,3)
enddo
a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2)
a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3)
a(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
a(1,2) = b(1,3)*b(3,2) - b(1,2)*b(3,3)
a(2,2) = b(1,1)*b(3,3) - b(1,3)*b(3,1)
a(3,2) = b(1,2)*b(3,1) - b(1,1)*b(3,2)
a(1,3) = b(1,2)*b(2,3) - b(1,3)*b(2,2)
a(2,3) = b(1,3)*b(2,1) - b(1,1)*b(2,3)
a(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
end
subroutine invert4(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
integer, intent(in) :: LDA
integer, intent(in) :: na
double precision, intent(inout) :: det_l
double precision :: b(4,4)
integer :: i,j
double precision :: f
det_l = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
-a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) &
+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) &
+a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) &
+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) &
-a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)) &
-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) &
+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))
do i=1,4
b(1,i) = a(1,i)
b(2,i) = a(2,i)
b(3,i) = a(3,i)
b(4,i) = a(4,i)
enddo
a(1,1) = b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
a(3,1) = b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(1,2) = -b(1,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(1,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(1,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
a(2,2) = b(1,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(1,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(1,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
a(3,2) = -b(1,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(1,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(1,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(4,2) = b(1,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(1,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(1,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(1,3) = b(1,2)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))-b(1,3)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))+b(1,4)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))
a(2,3) = -b(1,1)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))+b(1,3)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))-b(1,4)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))
a(3,3) = b(1,1)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))-b(1,2)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))+b(1,4)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
a(4,3) = -b(1,1)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))+b(1,2)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))-b(1,3)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
a(1,4) = -b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))+b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))-b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))
a(2,4) = b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))-b(1,3)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))+b(1,4)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))
a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
a(4,4) = b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
end
对于如此小的矩阵，我只需调用 LAPACK ...
! Computes inverse of a matrix
subroutine MatInv(A,Ainv)
implicit none
real(8) :: A(:,:),Ainv(:,:)
real(8),allocatable :: work(:)
integer :: n,info
integer,allocatable :: ipiv(:)
Ainv=A
n=size(Ainv,1)
allocate(ipiv(n),work(n))
call dgetrf(n,n,Ainv,n,ipiv,info)
call dgetri(n,Ainv,n,ipiv,work,n,info)
deallocate(ipiv,work)
end subroutine MatInv
由于潜在的错误和代码膨胀，不值得推出自己的东西。