Fortran求逆矩阵子程序
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
! 求逆矩阵
subroutine inverse(A,IA)
implicit none
real :: A(:,:), IA(:,:)
real, allocatable :: B(:,:)
integer :: i,j,N
N = size(A,1)
allocate(B(N,N))
! 先把IA设定成单位矩阵
forall(i=1:N,j=1:N,i==j) IA(i,j)=1.0
forall(i=1:N,j=1:N,i/=j) IA(i,j)=0.0
! 保存原先的矩阵A, 使用B来计算
B=A
! 把B化成对角线矩阵(除了对角线外,都为0) call Upper(B,IA,N) ! 先把B化成上三角矩阵call Lower(B,IA,N) ! 再把B化成下三角矩阵! 求解
forall(i=1:N) IA(i,:)=IA(i,:)/B(i,i)
return
end subroutine
! 求上三角矩阵的子程序
subroutine Upper(M,S,N)
implicit none
integer :: N
real :: M(N,N)
real :: S(N,N)
integer :: I,J
real :: E
do I=1,N-1
do J=I+1,N
E=M(J,I)/M(I,I)
M(J,I:N)=M(J,I:N)-M(I,I:N)*E
S(J,:)=S(J,:)-S(I,:)*E
end do
end do
return
end subroutine Upper
! 求下三角矩阵的子程序subroutine Lower(M,S,N)
implicit none
integer :: N
real :: M(N,N)
real :: S(N,N)
integer :: I,J
real :: E
do I=N,2,-1
do J=I-1,1,-1
E=M(J,I)/M(I,I)
M(J,1:N)=M(J,1:N)-M(I,1:N)*E
S(J,:)=S(J,:)-S(I,:)*E
end do
end do
return
end subroutine Lower