        SUBROUTINE YSOLVE(A,N,NP,B,M,MP,ISTAT)
*
*  Module number: FOS/HRS utility
*
*  Module name: ysolve
*
*  Keyphrase:
*  ----------
*       Solution of Linear System
*  Description:
*  ------------
*       Numerical Recipe Routine.
*       This routine solves a linear system using Gauss-Jordan
*       elimination.
*
*  FORTRAN name: ysolve.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*       none
*
*  Subroutines Called:
*  -------------------
*  SDAS:
*
*  History:
*  --------
*  Version      Date        Author          Description
*       1       Aug 87  D. Lindler      Conversion of Numerical Recipe routine
*-------------------------------------------------------------------------------
C
C INPUT PARAMETERS
C
C       A - NP x NP matrix (double precision)
C       N - number of equations and unknowns (integer)
C       NP - dimensions of A (integer)
C       B - array of vectors to solve for (NP x MP) (double precision)
C       MP - number of vectors to solve
C
C OUTPUT PARAMETERS
C
C       A - replaced with its inverse
C       B - replaced with solution vectors
C
C-----------------------------------------------------------------------
        DOUBLE PRECISION A,B,BIG,DUM,PIVINV
        INTEGER IPIV,INDXR,INDXC,J,I,IROW,ICOL,K,ISTAT,NP,L,N,LL,MP,M
        DIMENSION A(NP,NP),B(NP,MP),IPIV(100),INDXR(100),INDXC(100)

        DO 10 J=1,N
                IPIV(J)=0
10      CONTINUE
        DO 110 I=1,N
                BIG=0.
                DO 30 J=1,N
                        IF(IPIV(J).NE.1)THEN
                                DO 20 K=1,N
                                        IF (IPIV(K).EQ.0) THEN
                                          IF (DABS(A(J,K)).GE.BIG)THEN
                                                  BIG=DABS(A(J,K))
                                                  IROW=J
                                                  ICOL=K
                                          ENDIF
                                        ELSE IF (IPIV(K).GT.1) THEN
                                                  ISTAT=1
                                                GO TO 999
                                        ENDIF
20                              CONTINUE
                        ENDIF
30              CONTINUE
                IPIV(ICOL)=IPIV(ICOL)+1
                IF (IROW.NE.ICOL) THEN
                        DO 40 L=1,N
                                DUM=A(IROW,L)
                                A(IROW,L)=A(ICOL,L)
                                A(ICOL,L)=DUM
40                      CONTINUE
                        DO 50 L=1,M
                                DUM=B(IROW,L)
                                B(IROW,L)=B(ICOL,L)
                                B(ICOL,L)=DUM
50                      CONTINUE
                ENDIF
                INDXR(I)=IROW
                INDXC(I)=ICOL
                IF (A(ICOL,ICOL).EQ.0.) THEN
                        ISTAT=1
                        GO TO 999
                ENDIF
                PIVINV=1./A(ICOL,ICOL)
                A(ICOL,ICOL)=1.
                DO 60 L=1,N
                        A(ICOL,L)=A(ICOL,L)*PIVINV
60              CONTINUE
                DO 70 L=1,M
                        B(ICOL,L)=B(ICOL,L)*PIVINV
70              CONTINUE
                DO 100 LL=1,N
                        IF(LL.NE.ICOL)THEN
                                DUM=A(LL,ICOL)
                                A(LL,ICOL)=0.
                                DO 80 L=1,N
                                        A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
80                              CONTINUE
                                DO 90 L=1,M
                                        B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
90                              CONTINUE
                        ENDIF
100             CONTINUE
110     CONTINUE
        DO 130 L=N,1,-1
                IF(INDXR(L).NE.INDXC(L))THEN
                        DO 120 K=1,N
                                DUM=A(K,INDXR(L))
                                A(K,INDXR(L))=A(K,INDXC(L))
                                A(K,INDXC(L))=DUM
120                     CONTINUE
                ENDIF
130     CONTINUE
        ISTAT=0
999     RETURN
        END
