C -----------------------------------------------------------------------
C        Subroutine LUFBSUB.FOR by www.numerical-methods.com             |
C -----------------------------------------------------------------------
C
C Subroutine LUFBSUB returns the solution of a real matrix-vector system
C  of equations
C                               A x = b
C
C  in which the original matrix A has been replaced by an LU
C  factorisation, that is the matrices P, L and U
C
C                               PA  = LU
C
C  where L and U have over-written the original matrix A [ The
C  lower triangle of the returned matrix with a set of 1s down the
C  diagonal is the matrix L. The upper triangle and the diagonal of the
C  returned matrix is the matrx U. perm records the row interchanges
C  that result from pivotting; perm(i)=j means that the P(i,j)=1.]
C  The input matrix components to LUFBSUB.FOR are typically the output
C  from a LU factorisation routine such as LUFAC.FOR.
C
C Sub LUFBSUB(A, N, PERM, B)
C  real       A: L and U stored in A
C  integer    N: the dimension of the matrix/vector
C  integer PERM: an n-vector, the column index of the permutation matrix P
C  real       B: the vector b
C
C As a result of the subroutine, the input vector b is over-written
C and then contains the solution x.
C
C Source of the code: http://www.numerical-methods.com/fortran/LUFBSUB.for
C Source of the user-guide: http://www.numerical-methods.com/fortran/
C  LUfbsub.htm
C
C Licence: This is 'open source'; the software may be used and applied
C  within other systems as long as its provenance is appropriately
C  acknowledged. See the GNU Licence http://www.gnu.org/licenses/lgpl.txt
C  for more information or contact webmaster@numerical-methods.com

      SUBROUTINE LUFBSUB(MAXN, N, A, PERM, B, PB)
C Input parameters
C ----------------
C The limit on the dimension of the matrices A and B 
      INTEGER    MAXN
C The dimension of the matrices
      INTEGER    N
C The matrix A
      REAL*8     A(MAXN,MAXN)
C The vector PERM
      INTEGER    PERM(MAXN)
C The vector b
      REAL*8     B(MAXN)
C Workspace
      REAL*8     PB(MAXN)

      INTEGER I,J
      REAL*8 SUM

C Initialise pb as b with rows exchanged as indicated by perm
      DO 100 I=1,N
        PB(I)=B(PERM(I))
100   CONTINUE

C Solve Ly=b by forward substitution (pb stores y on completion)
      DO 200 I=1,N
        SUM=PB(I)
        DO 210 J=1,I-1
          SUM=SUM-A(I,J)*PB(J)
210     CONTINUE
      PB(I)=SUM
200   CONTINUE

C Solve Ux=y by back substitution (pb stores the solution 'x' on
C  completion)
      DO 300 I=N,1,-1
        SUM=PB(I)
        DO 310 J=I+1,N
          SUM=SUM-A(I,J)*PB(J)
310     CONTINUE
        PB(I)=SUM/A(I,I)
300   CONTINUE


C b is over-written by pb, the solution 'x'.
      DO 400 I=1,N
        B(I)=PB(I)
400   CONTINUE


      END