-
!
此法利用
LU
分解<
/p>
然后再利用上三角矩
阵
求逆方法
其中利用矩
p>
阵
乘法
以及
转
置
估
计
精度不
敢保
证
,不
过对
于低
阶
矩
阵
,精度可以
达
到要求。
progra
m
inverse_
m
atrix
!
利用上三角矩
阵
逆矩
阵
的求法
< br>
构
造
两
个
三角
矩
阵
再分
开
求逆
再求乘
积
im
plicit none
integer i,n
real*8,allocatab
l
e::a(:,:),l(:,:),u(:,:),lt(:,:),invers
e_lt(:,:),inverse
_l(:,:),inverse_u(:,:)
,inverse_a(:,:)
write(*,*) 'Please
input the order of M
atrix
A
:'
read(*,*) n <
/p>
allocate(a(n,n),l(n,n),u(n,n),lt(n,n),in
verse_lt(n,n),inverse_l(n,
n),inverse_u(
n,n),inverse_a(n,n))
do i=1,n
write(*,
',i2,'
row of
the m
atrix a is:')
read(*,*) a(i,:)
end do
call LU_break(a,n,l,u)
call
inverse_uptri_
m
atrix(u,n,inv
erse_u)
call
transpose_
m
atrix(l,n,n,lt)
call inverse_uptri_
m
< br>atrix(lt,n,inverse_lt)
call transpos
e_
m
atrix(inverse_lt,n,n,inve
rse_l)
call m
ultiply_
m
atrix(inverse_u,n,n,inverse_l,n,in
verse_a)
do i=1,n
write(*,
',i2,'
row
of
the
inverse
of
m
atrix
a
is:')
write(*,*) inverse_a(i,:)
end do
stop
end
subroutine LU_break(A,n,L,U)
!
普通
LU
分解,
目的是
为
了便于求解逆矩
阵
im
plicit none
integer
k,i,j,t,n
real*8::A(n,n),L(n,n),U(n,n)
real*8
sum
1,sum
2
do k=1,n
do
j=k,n
sum
1=0
do
t=1,k-
1
sum
1=sum
< br>1+L(k,t)*U(t,j)
end do
U(k,j)=A(k,j)
-
p>
sum
1
end do
do
i=k+1,n
sum
2=0
do
t=1,k-
1
sum
2=sum
< br>2+L(i,t)*
U(t,k)
end do
L(i
,k)=(A(i,k)
-
sum
2)
/
U(k,k)
end do
end do
do j=1,n
L(j,j)=1
end do
return
end
subroutine transpose_
m
atrix(A,n,m
,B)
!
求矩
阵转
置
im
plicit none
integer i,j,n,m
-
-
-
-
-
-
-
-
-
上一篇:调酒用具的知识
下一篇:最全最实用的自然拼读法教学讲义