!此法利用LU分解 然后再利用上三角矩阵求逆方法 其中利用矩阵乘法 以及转置 估计精度不敢保证，不过对于低阶矩阵，精度可以达到要求。
program inverse_matrix !利用上三角矩阵逆矩阵的求法 构造两个三角矩阵 再分开求逆 再求乘积
implicit none
integer i,n
real*8,allocatable::a(:,:),l(:,:),u(:,:),lt(:,:),inverse_lt(:,:),inverse_l(:,:),inverse_u(:,:),inverse_a(:,:)
write(*,*) 'Please input the order of Matrix A:'
read(*,*) n
allocate(a(n,n),l(n,n),u(n,n),lt(n,n),inverse_lt(n,n),inverse_l(n,n),inverse_u(n,n),inverse_a(n,n))
do i=1,n
write(*,\"('The ',i2,' row of the matrix a is:')\") i
read(*,*) a(i,:)
end do
call LU_break(a,n,l,u)
call inverse_uptri_matrix(u,n,inverse_u)
call transpose_matrix(l,n,n,lt)
call inverse_uptri_matrix(lt,n,inverse_lt)
call transpose_matrix(inverse_lt,n,n,inverse_l)
call multiply_matrix(inverse_u,n,n,inverse_l,n,inverse_a)
do i=1,n
write(*,\"('The ',i2,' row of the inverse of matrix a is:')\") i
write(*,*) inverse_a(i,:)
end do
stop
end
subroutine LU_break(A,n,L,U) !普通LU分解，目的是为了便于求解逆矩阵
implicit none
integer k,i,j,t,n
real*8::A(n,n),L(n,n),U(n,n)
real*8 sum1,sum2
do k=1,n
do j=k,n
sum1=0
do t=1,k-1
sum1=sum1+L(k,t)*U(t,j)
end do
U(k,j)=A(k,j)-sum1
end do
do i=k+1,n
sum2=0
do t=1,k-1
sum2=sum2+L(i,t)*U(t,k)
end do
L(i,k)=(A(i,k)-sum2)/U(k,k)
end do
end do
do j=1,n
L(j,j)=1
end do
return
end
subroutine transpose_matrix(A,n,m,B) !求矩阵转置
implicit none
integer i,j,n,m
real*8 A(n,m),B(m,n),t
do i=1,m
do j=1,n
B(i,j)=A(j,i)
end do
end do
return
end
subroutine multiply_matrix(A,n,m,B,l,C) !矩阵间乘法
implicit none
integer i,j,k,l,n,m
real*8 sum
real*8 A(n,m),B(m,l),C(n,l)
do i=1,n
do j=1,l
sum=0
do k=1,m
sum=sum+A(i,k)*B(k,j)
end do
C(i,j)=sum
end do
end do
return
end
subroutine inverse_uptri_matrix(A,n,inverse_A) !上三角矩阵的逆矩阵函数
implicit none
integer i,j,k,n
real*8 sum
real*8 A(n,n),inverse_A(n,n)
do i=1,n
inverse_A(i,i)=1/A(i,i)
end do
do i=1,n
do j=1,n
if(i do k=i,j-1 sum=sum+inverse_A(i,k)*A(k,j) end do inverse_A(i,j)=-sum/A(j,j) end if end do end do return end 因篇幅问题不能全部显示，请点此查看更多更全内容