Coding the forward-backward substitution

The routine accepts the LU factorization of the coeficient matrix, stored compactly in , the vector of pivots and a right hand side vector . It returns the solution of the system stored in the vector (i.e. is overwritten with the solution).

A sample substitution routine might be

! Check if the arguments match
IF ( size(A,1).NE.size(b,1) ) THEN
  PRINT*,"Error in solve: A and b must match"
  RETURN
END IF
n=SIZE(A,1) ! the dimension of the matrix
!
! Permute the vector b,   b <- Pb
!
DO j=1,n-1  ! columns 1:n-1
!
  IF (Pivot(j).NE.j) THEN     ! If pivot is non-diagonal
     Tmp         = b(j)       ! permute elements j and pivot(j)
     b(j)        = b(Pivot(j))
     b(Pivot(j)) = Tmp
  END IF
!
END DO

!
!  Forward Substitution,  y <- L\b
!
DO i=2,n  
  b(i) = b(i) - DOT_PRODUCT( A(i,1:i-1), b(1:i-1) )
END DO           
!
!
!  Backward Substitution, x <- U\y
!
b(n) = b(n)/A(n,n)
DO i=n-1,1,-1                ! elements to be processed
  b(i) = (b(i) - DOT_PRODUCT(A(i,i+1:n),b(i+1:n)))/A(i,i)
END DO