!此 法 利 用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
