module LU

!----------------------------------------module coment
!  Version     :  V1.0    
!  Coded by    :  lh
!  Date        :  2012-8-2
!-----------------------------------------------------
!  Description :  LU 分解解方程
!    
!-----------------------------------------------------

contains
subroutine solve(A,b,x,N)
implicit none
integer :: N
real :: A(N,N),b(N),x(N)
real :: L(N,N),U(N,N)
real :: y(N)
call doolittle(A,L,U,N)
call downtri(L,b,y,N)
call uptri(U,y,x,N)
end subroutine solve

subroutine doolittle(A,L,U,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  lh
!  Date      :  
!-----------------------------------------------------
!  Purpose   :  LU分解之Doolittle函数
!              A=LU
!-----------------------------------------------------
!  Input  parameters  :
!       1.    A  方阵
!       2.    N  阶数
!  Output parameters  :
!       1.   L
!       2.   U
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.
!       2.
!----------------------------------------------------

implicit none
integer::N,i,k,r,m,s,j
real::A(N,N),L(N,N),U(N,N)
!U的第一行
U(1,:)=A(1,:)
!L的第一列
L(:,1)=A(:,1)/U(1,1)
do k=2,N  
    L(k,k)=1
    do j=k,N
       s=0
       do m=1,k-1
        s=s+L(k,m)*U(m,j)
       end do
       u(k,j)=A(k,j)-s
   end do
   do i=k+1,n
     s=0
     do m=1,k-1
      s=s+L(i,m)*U(m,k)
     end do
     L(i,k)=(A(i,k)-s)/U(k,k)
   end do
end do
end subroutine doolittle

subroutine uptri(A,b,x,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  lh
!  Date      :  2012-8-2
!-----------------------------------------------------
!  Purpose   :  上三角方程组的回带方法
!                 Ax=b
!-----------------------------------------------------
!  Input  parameters  :
!       1.   A(N,N)系数矩阵
!       2.   b(N)右向量
!       3.   N方程维数
!  Output parameters  :
!       1.  x  方程的根
!       2.
!  Common parameters  :
!
!----------------------------------------------------
implicit none
  integer::i,j,N
  real::A(N,N),b(N),x(N)
  x(N)=b(N)/A(N,N)
!回带部分
   do i=n-1,1,-1
    x(i)=b(i)
      do j=i+1,N
       x(i)=x(i)-A(i,j)*x(j)
       end do
    x(i)=x(i)/A(i,i)
end do
end subroutine uptri


subroutine downtri(A,b,x,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  lh
!  Date      :  2012-8-2
!-----------------------------------------------------
!  Purpose   :  下三角方程组的回带方法
!                 Ax=b
!-----------------------------------------------------
!  Input  parameters  :
!       1.   A(N,N)系数矩阵
!       2.   b(N)右向量
!       3.   N方程维数
!  Output parameters  :
!       1.  x  方程的根
!       2.
!  Common parameters  :
!
!----------------------------------------------------

implicit none
integer::i,j,N,k
real::A(N,N),b(N),x(N)
x(1)=b(1)/A(1,1)
   do k=2,N
    x(k)=b(k)
    do i=1,k-1
      x(k)=x(k)-A(k,i)*x(i)
   end do
   x(k)=x(k)/A(k,k)
end do
end subroutine downtri
end module LU

program main
!--------------------------------------program comment
!  Version   :  V1.0    
!  Coded by  :  lh
!  Date      :  2012-8-2
!-----------------------------------------------------
!  Purpose   :   LU分解计算线性方程组
!              Ax=b
!-----------------------------------------------------
!  In put data  files :
!       1.    A,b
!       2.
!  Output data files  :
!       1.    x
!       2.
!-----------------------------------------------------
use LU
implicit none
integer,parameter::N=4
integer :: i,j
real::A(N,N),L(N,N),U(N,N)
real::b(N),x(N)
open(unit=11,file='fin.txt')
open(unit=12,file='fout.txt')
read(11,*)
read(11,*)((A(i,j),j=1,N),i=1,N)
read(11,*) b
call solve(A,b,x,N)
write(12,101)x
101 format(T5,'LU分解计算线性方程组计算结果',//,4(/,F10.6))
end program main