Code: CROUT

       program crout

       interface
         subroutine croutludecomp(c,n)
           integer n
           real, dimension(1:n, 1:n) :: c
         end subroutine
         subroutine lubkfor(a,b,n)
            integer n
            real, dimension (1:n,1:n) :: a
            real, dimension (1:n)     :: b
          end subroutine
       end interface

       parameter(n =3)
       integer i
       real, dimension(n*(n))::vector
       real, dimension(n, n)::a
       real, dimension(n) :: b,x

       vector = [ 3.0, 5.0, 2.0, 0.0, 8.0, 2.0, 6.0, 2.0, 8.0]
       a = reshape( vector, [3,3], order= [2, 1])
       b = [ 8.0, -7.0, 26.0]
       x = [0.0,0.0,0.0]

       print *, ' Original matrix is : '
       do i = 1, n
         print *, a(i, 1:n)
       end do
       print *, ' right hand side vector is :'
       print *, b

       call croutludecomp(a,n)
       call lubkfor(a,b,n)

       print *, '   '

       print *, ' Final LU decomposition matrix  '
       do i = 1, n
         print *, a(i, 1:n)
       end do

       print *, ' Solution vector : '
       print *, b


       end

!**********************************************************************
!  LUBKFOR
!**********************************************************************
!  This routine takes as input the output matrix of CROUTLUDECOMP
!  and a right hand side vector, and returns the solution of the
!  original set of equations Ax=b.  The solution vector x is
!  returned in b
!**********************************************************************
   subroutine lubkfor(a,b,n)
   integer n
   real, dimension (1:n,1:n) :: a
   real, dimension (1:n)     :: b
   ! local var
   real, dimension (1:n)     :: y
   integer i

   ! Forward substitution: find y in the equations Ly=b

   y(1) = b(1)
   do i=2,n
      y(i) = b(i) - sum( a(i,1:i-1) * y(1:i-1) )
   enddo

   ! Backward substitution: find x in the equations Ux=y
   !     (put solution x in b)

   b(n) = y(n) / a(n,n)
   do i=n-1,1,-1
      b(i) = (y(i) - sum( a(i,i+1:n) * b(i+1:n) ) )/a(i,i)
   enddo
end


!**********************************************************************
!  CROUTLUDECOMP
!**********************************************************************
!  This routine employs Crout's algorithm to perform an LU decompostion
!  of the input matrix a.  The result is returned in a.
!**********************************************************************
   subroutine croutludecomp(a,n)
   integer n
   real, dimension (1:n, 1:n)  :: a
   ! local var
   integer i,j

!  main loop for crout's algorithm

   do j = 1,n
      do i=2,j
         a(i,j) = a(i,j) - sum( a(i,1:i-1) * a(1:i-1,j) )
      enddo
      if(abs(a(j,j)) .lt. 1.0e-6) then
         write(*,*) 'error: croutludecomp: input matrix is singular'
         write(*,*) '       to machine precision.  Need pivoting'
         stop
      endif
      do i=j+1,n
         if(j .eq. 1) then
            a(i,j) = a(i,j)/a(j,j)
         else
            a(i,j) = ( a(i,j) - sum( a(i,1:j-1) * a(1:j-1,j) ) )/a(j,j)
         endif
      enddo
   enddo
end