LU矩阵直接三角分解法

program zhijie !
implicit none
real(8),dimension(:,:),allocatable::a
real(8),dimension(:),allocatable::b,x
integer::i,j,n
print,"输入维度n"
read(,)n
allocate(a(n,n),b(n),x(n))
print,"输入nn的数组a"
read(,)((a(i,j),j=1,n),i=1,n)
print,"n的数组b"
read(,),b
call zhijie_s_j(n,a,b,x)
print,x
end program zhijie
!****************************************************
subroutine zhijie_s_j(n,a,b,x)
implicit none
integer::i,j,k,n
real(8)::sum1,sum2,sum3,sum4,d1,d2,d3,d4
real(8),dimension(n,n)::a,u,l
real(8),dimension(n)::b,c,x,y
do i=1,n
u(1,i)=a(1,i)
end do
do i=2,n
l(i,1)=a(i,1)/u(1,1)
end do
do j=2,n
do i=j,n
sum1=0
do k=1,j-1
d1=l(j,k)u(k,i)
sum1=sum1+d1
end do
u(j,i)=a(j,i)-sum1
end do
do i=j+1,n
if(j/=n)then
sum2=0
do k=1,j-1
d2=l(i,k)u(k,j)
sum2=sum2+d2
end do
l(i,j)=(a(i,j)-sum2)/u(j,j)
end if
end do
end do ! 计算U和L的数值
!***************************************************
y(1)=b(1)
do i=2,n
sum3=0
do k=1,i-1
d3=l(i,k)y(k)
sum3=sum3+d3
end do
y(i)=b(i)-sum3
end do !计算y
x(n)=y(n)/u(n,n)
do i=n-1,1,-1
sum4=0
do k=i+1,n
d4=u(i,k)*x(k)
sum4=sum4+d4
end do
x(i)=(y(i)-sum4)/u(i,i)
end do !计算x
end subroutine zhijie_s_j