Problem 2

I wrote a subroutine that solves a system of linear equations by LU decomposition using Crout's algorithm. This routine is of no practical use, as Crout's method alone is unconditionally unstable. The algorithm takes an n x n matrix A, and returns an n x n matrix whose interpretation is A = L U, where L and U are lower and upper triangular matrices, respectively. The diagonal elements have been chosen to have a value of 1.0, so both the L and U matrices can be packed into the original A matrix.

The method is summarized as follows:

• For each do the following two procedures:
  • First, for , assign the matrix element as

    

    (when i=1 in the above equation, the sum is taken to be zero).

  • Second, for assign the matrix element as

    

This algorithm is employed in my routine CROUTLUDECOMP. The first step in the main loop of the algorithm is executed by the code

      do i=2,j
         a(i,j) = a(i,j) - sum( a(i,1:i-1) * a(1:i-1,j) )
      enddo

The second step is executed by the code

      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

Once Crout's algorithm has been employed, it is a simple matter of frontwards and backwards substitution to find the solution:

   ! 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

• Output from CROUT
• Code: CROUT