        SUBROUTINE GAUSSJ(A,N,NP,B,M,MP,ISTAT)
*
*  Module number: FOS/HRS utility
*
*  Module name: gaussj
*
*  Keyphrase:
*  ----------
*       Solution of Linear System
*  Description:
*  ------------
*       Numerical Recipe Routine.
*       This routine solves a linear system using Gauss-Jordan
*       elimination.
*
*  FORTRAN name: gaussj.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*       none
*
*  Subroutines Called:
*  -------------------
*  SDAS:
*	UMSPUT
*
*  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-----------------------------------------------------------------------
      IMPLICIT NONE
C
C     INCLUDE FILE FOR THE IRAF77 FORTRAN INTERFACE TO THE IRAF VOS
C
C
C     FILE I/O ACCESS MODES
C
      INTEGER   RDONLY
      PARAMETER (RDONLY = 1)
      INTEGER   RDWRIT
      PARAMETER (RDWRIT = 2)
      INTEGER   WRONLY
      PARAMETER (WRONLY = 3)
      INTEGER   APPEND
      PARAMETER (APPEND = 4)
C
C     CODES FOR DATA TYPES
C
      INTEGER   TYBOOL
      PARAMETER (TYBOOL = 1)
      INTEGER   TYCHAR
      PARAMETER (TYCHAR = 2)
      INTEGER   TYINT
      PARAMETER (TYINT = 4)
      INTEGER   TYREAL
      PARAMETER (TYREAL = 6)
      INTEGER   TYDOUB
      PARAMETER (TYDOUB = 7)
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
      INTEGER USRLOG
      PARAMETER (USRLOG = 4)
C
C     UHDAS HEADER PARM TYPES -- CB, DAO, 5-SEP-87
C
      INTEGER GENHDR
      PARAMETER (GENHDR = 0)
      INTEGER IMSPEC
      PARAMETER (IMSPEC = 1)
C
C     THIS SECTION IS FOR PARAMETERS RELEVANT TO TABLE I/O.
C
C     THESE MAY BE SET BY UTPPTI AND/OR READ BY UTPGTI:
C
C                                       LENGTH OF ROW (UNIT = SIZE OF REAL)
      INTEGER   TBRLEN
      PARAMETER (TBRLEN = 1)
C                                       INCREASE ROW LENGTH
      INTEGER   TBIRLN
      PARAMETER (TBIRLN = 2)
C                                       NUMBER OF ROWS TO ALLOCATE
      INTEGER   TBALLR
      PARAMETER (TBALLR = 3)
C                                       INCREASE ALLOC NUM OF ROWS
      INTEGER   TBIALR
      PARAMETER (TBIALR = 4)
C                                       WHICH TYPE OF TABLE? (ROW OR COLUMN)
      INTEGER   TBWTYP
      PARAMETER (TBWTYP = 5)
C                                       MAXIMUM NUMBER OF USER PARAMETERS
      INTEGER   TBMXPR
      PARAMETER (TBMXPR = 6)
C                                       MAXIMUM NUMBER OF COLUMNS
      INTEGER   TBMXCL
      PARAMETER (TBMXCL = 7)
C                                       TYPE = ROW-ORDERED TABLE
      INTEGER   TBTYPR
      PARAMETER (TBTYPR = 11)
C                                       TYPE = COLUMN-ORDERED TABLE
      INTEGER   TBTYPC
      PARAMETER (TBTYPC = 12)
C
C     THESE MAY BE READ BY UTPGTI BUT MAY NOT BE SET:
C
C                                       NUMBER OF ROWS WRITTEN TO
      INTEGER   TBNROW
      PARAMETER (TBNROW = 21)
C
C     END IRAF77.INC
        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)
        CHARACTER*80 MESS
C
C  DATA DECLARATIONS
C
        DATA MESS/'Attempt to invert a singular matrix'/
        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
                                                CALL UMSPUT(MESS,STDOUT+
     *                                            STDERR+USRLOG,0,ISTAT)
                                                  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
                        CALL UMSPUT(MESS,STDOUT+USRLOG+STDERR,0,ISTAT)
                        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
