数值计算 求矩阵逆 fortran源代码

!

此法利用

LU

分解< /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)

-

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