[Fortran] 纯文本查看 复制代码
program test
implicit none
real(kind=8),dimension(3,3) :: A, L, U
real(kind=8),dimension(3) ::B, ANS, C
integer :: n
n=3
data A/1,2,3,2,5,1,3,2,5/
B=(/14,18,20/)
call LU_method(n,A,B,L,U,C,ANS)
print*,ANS
end program test
subroutine LU_method(n,A,B,L,U,C,ANS)
implicit none
integer ::n,i,j,k,r
real(kind=8) ::m1,m2,sum
real(kind=8),dimension(n,n) :: A, L, U
real(kind=8),dimension(n) :: B, ANS, C
do i=2,n
L(i,1)=A(i,1)/U(1,1)
end do
do j=1,n
U(1,j)=A(1,j)
end do
do r=2,n
if(A(r,r)==0)then
print*,'解不出来'
stop
else
if(r < j)then
U(r,j)=0
else
do k=1,r-1
m1=0
do j=r,n
m1=m1+L(r,k)*U(k,j)
U(r,j)=A(r,j)-m1
end do
end do
end if
if(r > i)then
L(i,r)=0
else
do k=1,r-1
m2=0
do i=r+1,n
m2=m2+L(i,k)*U(k,r)
L(i,r)=(A(i,r)-m2)/U(r,r)
end do
end do
end if
end if
!求解L*C=B
C(1)=B(1)
do i=2,n
sum=0
do k=1,i-1
sum=sum+L(i,k)*C(k)
end do
C(i)=B(i)-sum
end do
!求解U*ANS=C
ANS(n)=C(n)/U(n,n)
do i=n-1,1,-1
sum=0
do k=n-1,1,-1
sum=sum+U(i,k)*ANS(k)
end do
ANS(i)=(C(i)-sum)/U(i,i)
end do
end do
end subroutine LU_method