program lu1

      implicit none

      integer n
      parameter (n=4)
      double precision A(n,n), b(n)
      integer i, j

      A(1,1) = 1.0d0
      A(1,2) = 4.0d0
      A(1,3) = 5.0d0
      A(1,4) = 6.0d0
      A(2,1) = 1.0d0
      A(2,2) = 1.0d0
      A(2,3) = 1.0d0
      A(2,4) = 1.0d0
      A(3,1) = 3.0d0
      A(3,2) = 4.0d0
      A(3,3) = -1.0d0
      A(3,4) = 4.0d0
      A(4,1) = 6.0d0
      A(4,2) = 4.0d0
      A(4,3) = 3.0d0
      A(4,4) = 3.0d0

      b(1) = 1.0d0
      b(2) = 1.0d0
      b(3) = 1.0d0
      b(4) = 1.0d0

c     ---------------------- factor, solve
      call LUfactor1 ( A, n )
      call LUsolve1 ( A, b, n )

c     ---------------------- print LU, b solution
      print*, "A = "
      do i = 1, n
         print*, (A(i,j),j=1,n)
      end do
      print*, "b = ", b

      stop
      end


c     replaces matrix A with LU factorization

      subroutine LUfactor1 ( A, n )  ! factor A=LU, but not pivoting

      implicit none

c     input variables
      integer n
      double precision A(n,n)

c     local variables
      integer i, j, k
      double precision Ajj, Lij
      
      do j = 1, n
         Ajj = A(j,j)
         if (abs(Ajj) < 1.0d-15) then
            print *, "Zero Diagonal element -- LU fails"
            return
         end if
         do i = j+1, n
            Lij = A(i,j)/Ajj
            A(i,j) = Lij
            do k = j+1, n
               A(i,k) = A(i,k) - Lij*A(j,k)
            end do
         end do
      end do

      return
      end

c     replaces b vector with solution x in Ax=b

      subroutine LUsolve1 ( A, b, n )

      implicit none

c     input variables
      integer n
      double precision A(n,n)
      double precision b(n)

c     local variables
      double precision y(n)
      double precision rowsum
      integer i, j

      y(1) = b(1)
      do i = 2, n
         y(i) = b(i)
         rowsum = 0.0d0
         do j = 1, i-1
            rowsum = rowsum + A(i,j)*y(j)
         end do
         y(i) = y(i) - rowsum
      end do

      b(n) = y(n)/A(n,n)
      do i=n-1,1,-1
         b(i) = y(i)
         rowsum = 0.0d0
         do j = n,i+1,-1
            rowsum = rowsum + A(i,j)*b(j)
         end do
         b(i) = (b(i)-rowsum)/A(i,i)
      end do

      return
      end