      SUBROUTINE LINSYS(A,B,X,N,IORDER)
C
C   SOLVES THE LINEAR SYSTEM AX=B USING LU DECOMPOSITION
C   THIS ROUTINE ASSUMES THAT A IS STORED IN THE MAIN PROGRAM AS A
C   ONE-DIMENSIONAL ARRAY OF DIMENSION N*N, SO THE TRANSPOSE OF A
C   MUST BE TAKEN BEFORE CALLING LUDCMP AND SOLVLU.
C   RETURNS:
C   X       THE SOLUTION TO THE LINEAR SYSTEM
C   A       REPLACED BY ITS LU DECOMPOSITION IN COMPACT FORM
C   B       REORDERED TO CORRESPOND TO ROW INTERCHANGES PERFORMED ON A
C   IORDER  THE ORDER OF THE ELEMENTS OF B
C
      DIMENSION A(N,N),B(N),X(N)
      DIMENSION IORDER(N)
C
C   FIND TRANSPOSE OF A SO THAT THE 2 DIMENSIONAL REFERENCES TO A
C   CORRESPOND TO THE 1 DIMENSIONAL DECLARATION IN THE MAIN PROGRAM
C
      DO 10 I=1,N-1
         DO 20 J=I+1,N
            A(I,J)=A(J,I)
   20    CONTINUE
   10 CONTINUE
C
      CALL LUDCMP(A,N,N,IORDER)
      CALL SOLVLU(A,B,X,N,N,IORDER)
      RETURN
      END
C
      SUBROUTINE LUDCMP(A,N,NDIM,IORDER)
C
C   FROM "APPLIED NUMERICAL ANALYSIS" BY GERALD AND WHEATLEY,
C   THIRD EDITION (ADDISON-WESLEY, 1984).
C
C ---------------------------------------------------------------
C
C   THIS SUBROUTINE COMPUTES THE L AND U TRIANGULAR MATRICES
C   EQUIVALENT TO THE A MATRIX, SUCH THAT LU=A.  THESE MATRICES
C   ARE RETURNED IN THE SPACE OF A, IN COMPACT FORM.  THE U
C   MATRIX HAS ONES ON ITS DIAGONAL.  PARTIAL PIVOTING IS USED TO
C   GIVE MAXIMUM VALUED ELEMENTS ON THE DIAGONAL OF L.  THE ORDER
C   OF THE ROWS AFTER PIVOTING IS RETURNED IN THE INTEGER VECTOR
C   IORDER.  THIS SHOULD BE USED TO REORDER R.H.S. VECTORS BEFORE
C   SOLVING THE SYSTEM AX=B.
C
C ---------------------------------------------------------------
C
C   PARAMETERS:
C
C   A         THE N X N MATRIX OF COEFFICIENTS
C   N         THE NUMBER OF EQUATIONS
C   NDIM      THE FIRST DIMENSION OF A IN THE CALLING PROGRAM
C   IORDER    INTEGER VECTOR HOLDING ROW ORDER AFTER PIVOTING
C
C   THIS ROUTINE CALLS A SUBROUTINE APVT TO LOCATE THE PIVOT ROW AND
C   MAKE INTERCHANGES.
C
C ---------------------------------------------------------------
C
      DIMENSION A(NDIM,N)
      DIMENSION IORDER(N)
C
C ---------------------------------------------------------------
C
C   ESTABLISH INITIAL ORDERING IN ORDER VECTOR
C
      DO 10 I=1,N
         IORDER(I)=I
   10 CONTINUE
C
C   DO PIVOTING FOR FIRST COLUMN BY CALL TO SUBROUTINE APVT
C
      CALL APVT(A,N,NDIM,IORDER,1)
C
C ---------------------------------------------------------------
C
C   IF PIVOT ELEMENT VERY SMALL, PRINT ERROR MESSAGE AND RETURN
C
      IF (ABS(A(1,1)).LT.1.0E-5) THEN
         PRINT 100
         RETURN
      END IF
C
C   NOW COMPUTE ELEMENTS FOR FIRST ROW OF U
C
      DO 20 KCOL=2,N
         A(1,KCOL)=A(1,KCOL)/A(1,1)
   20 CONTINUE
C
C ---------------------------------------------------------------
C
C   COMPLETE THE COMPUTING OF L AND U ELEMENTS.  THE GENERAL PLAN IS TO
C   COMPUTE A COLUMN OF L'S, THEN CALL APVT TO INTERCHANGE ROWS, AND THEN
C   GET A ROW OF U'S.
C
      NM1=N-1
      DO 80 JCOL=2,NM1
C
C   FIRST COMPUTE A COLUMN OF L'S
C
      JM1=JCOL-1
      DO 50 IROW=JCOL,N
         SUM=0
         DO 40 KCOL=1,JM1
            SUM=SUM+A(IROW,KCOL)*A(KCOL,JCOL)
   40    CONTINUE
         A(IROW,JCOL)=A(IROW,JCOL)-SUM
   50 CONTINUE
C
C ---------------------------------------------------------------
C
C   NOW INTERCHANGE ROWS IF NEED BE, THEN TEST FOR TOO SMALL PIVOT
C
      CALL APVT(A,N,NDIM,IORDER,JCOL)
      IF (ABS(A(JCOL,JCOL)).LT.1.0E-5) THEN
         PRINT 100
         RETURN
      END IF
C
C   NOW WE GET A ROW OF U'S
C
      JP1=JCOL+1
      DO 70 KCOL=JP1,N
         SUM=0
         DO 60 IROW=1,JM1
            SUM=SUM+A(JCOL,IROW)*A(IROW,KCOL)
   60    CONTINUE
         A(JCOL,KCOL)=(A(JCOL,KCOL)-SUM)/A(JCOL,JCOL)
   70 CONTINUE         
   80 CONTINUE
C
C ---------------------------------------------------------------
C
C   STILL NEED TO GET LAST ELEMENT IN L MATRIX
C
      SUM=0
      DO 90 KCOL=1,NM1
         SUM=SUM+A(N,KCOL)*A(KCOL,N)
   90 CONTINUE
      A(N,N)=A(N,N)-SUM
      RETURN
C
  100 FORMAT(/,' VERY SMALL PIVOT ELEMENT INDICATES A NEARLY',
     +       ' SINGULAR MATRIX.')
      END
C
      SUBROUTINE APVT(A,N,NDIM,IORDER,JCOL)
C
C   FROM "APPLIED NUMERICAL ANALYSIS" BY GERALD AND WHEATLEY,
C   THIRD EDITION (ADDISON-WESLEY, 1984).
C
C ---------------------------------------------------------------
C
C   THIS SUBROUTINE FINDS THE LARGEST ELEMENT FOR PIVOT IN JCOL OF
C   MATRIX A, PERFORMS INTERCHANGES OF ELEMENTS IN A AND ALSO
C   INTERCHANGES THE ELEMENTS IN THE ORDER VECTOR.
C
C ---------------------------------------------------------------
C
C   PARAMETERS:
C
C   A         MATRIX OF COEFFICIENTS WHOSE ROWS ARE TO BE INTERCHANGED
C   N         NUMBER OF EQUATIONS
C   NDIM      FIRST DIMENSION OF A IN THE MAIN PROGRAM
C   IORDER    INTEGER VECTOR TO HOLD ROW ORDERING
C
C ---------------------------------------------------------------
C
      DIMENSION A(NDIM,N)
      DIMENSION IORDER(N)
C
C ---------------------------------------------------------------
C
C  FIND PIVOT ROW, CONSIDERING ONLY THE ELEMENTS ON AND BELOW DIAGONAL
C
      IPVT=JCOL
      BIG=ABS(A(JCOL,JCOL))
      JP1=JCOL+1
      DO 10 IROW=JP1,N
         ANEXT=ABS(A(IROW,JCOL))
         IF (ANEXT.GT.BIG) IPVT=IROW
   10 CONTINUE
C
C ---------------------------------------------------------------
C
C   NOW INTERCHANGE ROW ELEMENTS IN THE ROW WHOSE NUMBER EQUALS JCOL WITH
C   THE PIVOT ROW UNLESS PIVOT ROW IS JCOL.
C
      IF (IPVT.EQ.JCOL) THEN
         RETURN
      END IF
      DO 20 KCOL=1,N
         SAVE=A(JCOL,KCOL)
         A(JCOL,KCOL)=A(IPVT,KCOL)
         A(IPVT,KCOL)=SAVE
   20 CONTINUE
C
C ---------------------------------------------------------------
C
C   NOW SWITCH ELEMENTS IN THE ORDER VECTOR
C
      ISAVE=IORDER(JCOL)
      IORDER(JCOL)=IORDER(IPVT)
      IORDER(IPVT)=ISAVE
      RETURN
      END
C
      SUBROUTINE SOLVLU(RLU,B,X,N,NDIM,IORDER)
C
C   FROM "APPLIED NUMERICAL ANALYSIS" BY GERALD AND WHEATLEY,
C   THIRD EDITION (ADDISON-WESLEY, 1984).
C
C ---------------------------------------------------------------
C
C   THIS SUBROUTINE IS USED TO FIND THE SOLUTION TO A SYSTEM OF EQUATIONS,
C   AX=B, AFTER THE LU EQUIVALENT OF A HAS BEEN FOUND.  BEFORE USING THIS
C   ROUTINE, THE VECTOR B SHOULD BE SCALED IF MATRIX A WAS SCALED, USING
C   THE SAME SCALE FACTORS.  WITHIN THIS ROUTINE, THE ELEMENTS OF B ARE
C   REARRANGED IN THE SAME WAY THAT THE ROWS OF A WERE INTERCHANGED, USING
C   THE ORDER VECTOR WHICH HOLDS THE ROW ORDERINGS.  THE SOLUTION IS
C   RETURNED IN X.
C
C ---------------------------------------------------------------
C
C   PARAMETERS:
C
C   RLU       THE LU EQUIVALENT OF THE COEFFICIENT MATRIX
C   B         THE VECTOR OF RIGHT HAND SIDES
C   X         SOLUTION VECTOR
C   N         NUMBER OF EQUATIONS
C   NDIM      FIRST DIMENSION OF A IN THE MAIN PROGRAM
C   IORDER    INTEGER ARRAY OF ROW ORDER AS ARRANGED DURING PIVOTING
C
C ---------------------------------------------------------------
C
      DIMENSION RLU(NDIM,N),B(N),X(N)
      DIMENSION IORDER(N)
C
C ---------------------------------------------------------------
C
C   REARRANGE THE ELEMENTS OF THE B VECTOR.  X IS USED TO HOLD THEM.
C
      DO 10 I=1,N
         J=IORDER(I)
         X(I)=B(J)
   10 CONTINUE
C
C ---------------------------------------------------------------
C
C   COMPUTE THE B VECTOR, STORING BACK IN X
C
      X(1)=X(1)/RLU(1,1)
      DO 50 IROW=2,N
         IM1=IROW-1
         SUM=0
         DO 40 JCOL=1,IM1
            SUM=SUM+RLU(IROW,JCOL)*X(JCOL)
   40    CONTINUE
         X(IROW)=(X(IROW)-SUM)/RLU(IROW,IROW)
   50 CONTINUE
C
C ---------------------------------------------------------------
C
C   NOW GET THE SOLUTION VECTOR.  X(N)=B(N) ALREADY.
C
      DO 70 IROW=2,N
         NVBL=N-IROW+1
         SUM=0
         NP1=NVBL+1
         DO 60 JCOL=NP1,N
            SUM=SUM+RLU(NVBL,JCOL)*X(JCOL)
   60    CONTINUE
         X(NVBL)=X(NVBL)-SUM
   70 CONTINUE
C
      RETURN
      END
