program Doolittle

    implicit none
    integer i,j,n,t,k
    real,dimension (:,:),allocatable::a,l,u,w
    real,dimension (:,:),allocatable::b,x,y,v
    print*,"輸入方程階數：n"
    read *,n
    allocate(a(n,n),l(n,n),u(n,n),w(n,n),b(n,1),x(n,1),y(n,1),v(n,1))
    print*,"輸入系數矩陣A:n×n階："
    read*,((a(i,j),j=1,n),i=1,n)
     print*, "輸入矩陣B：n×1階："
    read*,b
    do i=2,n
       do j=1,i-1
          u(i,j)=0
       end do
    end do
    do i=1,n-1
      do j=i+1,n
          l(i,j)=0
       end do
    end do
    do i=1,n
      l(i,i)=1
    end do
    do j=1,n
       u(1,j)=a(1,j)
    end do
    do i=2,n
       l(i,1)=a(i,1)/u(1,1)
    end do
    do i=2,n
       do j=i,n
          w=0
          do t=1,i-1
            w(i,j)=w(i,j)+l(i,t)*u(t,j)
          end do
          u(i,j)=a(i,j)-w(i,j)
       end do
       j=i
       if(u(j,j)/=0)then
           do k=j+1,n
              w=0
              do t=1,j-1
                 w=w+l(k,t)*u(t,j)
              end do
              l(k,j)=(a(k,j)-w(i,j))/u(j,j)
           end do
        end if
     end do
     y(1,1)=b(1,1)
     do i=2,n
        v=0
        do j=1,i-1
           v=v+l(i,j)*y(j,1)
        end do
        y(i,1)=b(i,1)-v(i,1)
     end do
     x(n,1)=y(n,1)/u(n,n)
     do i=n-1,1,-1
        v=0
        do j=i+1,n
          v(i,1)=v(i,1)+u(i,j)*x(j,1)
        end do
        x(i,1)=(y(i,1)-v(i,1))/u(i,i)
     end do
     print*,"方程的解為："
     do i=1,n 
        write (*,*) "x(",i,")=",x(i,1)
     end do
    end program Doolittle