SUBROUTINE GAUSEL (MAX, N, A, NR, B, IERR)
C
C ****
C
C   FUNCTION              - COMPUTES SOLUTION TO A SET OF SIMULTANEOUS
C                             LINEAR EQUATIONS (DOES NOT GIVE PIVOT OR
C                             DETERMINANT DATA)
C
C   PARAMETERS  MAX       - MAXIMUM ROW DIMENSION OF B
C               N         - ORDER OF A
C               A(N,N)    - INPUT MATRIX OF COEFFICIENTS (DESTROYED)
C               NR        - NUMBER OF COLUMNS IN B
C               B(MAX,NR) - MATRIX OF CONSTANTS (REPLACED BY SOLUTIONS)
C               IERR      - INTEGER ERROR CODE
C                             = 0   NORMAL RETURN
C                             = 5   INPUT MATRIX IS SINGULAR
C
C   REQUIRED SUBPROGRAMS  - NONE
C
C   REQUIRED FUNCTIONS    - DABS
C
C   AUTHOR/IMPLEMENTER    - R.C. WARD / R.C. WARD
C
C   LANGUAGE              - FORTRAN (IBM DOUBLE PRECISION)
C
C   DATE RELEASED         - NOV. 14, 1975
C
C   LATEST VERSION        - sept. 1986 (minor modification:bb)
C
C ****
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N),B(MAX,NR)
      NM1 = N-1
      IF (NM1 .EQ. 0) GO TO 140
C ****
C     FIND LARGEST REMAINING ELEMENT IN I-TH COLUMN FOR PIVOT
C ****
      DO 100 I=1,NM1
         BIG = 0.D0
         DO 20 K=I,N
            TERM = DABS(A(K,I))
            IF (TERM - BIG) 20,20,10
  10        BIG = TERM
            L = K
  20     CONTINUE
         IF (BIG .EQ. 0.D0) GO TO 160
         IF (I-L) 50,80,50
C ****
C     PIVOT ROWS OF A AND B
C ****
  50     CONTINUE
         DO 60 J=1,N
            TEMP = A(I,J)
            A(I,J) = A(L,J)
            A(L,J) = TEMP
  60     CONTINUE
         DO 70 J=1,NR
            TEMP = B(I,J)
            B(I,J) = B(L,J)
            B(L,J) = TEMP
  70     CONTINUE
  80     CONTINUE
C ****
C     STORE PIVOT AND PERFORM COLUMN OPERATIONS ON A AND B
C ****
         IP1 = I+1
         DO 100 II=IP1,N
            A(II,I) = A(II,I)/A(I,I)
            X3 = A(II,I)
            DO 90 K=IP1,N
               A(II,K) = A(II,K) - X3*A(I,K)
  90        CONTINUE
            DO 100 K=1,NR
               B(II,K) = B(II,K) - X3*B(I,K)
 100  CONTINUE
C ****
C     PERFORM BACK SUBSTITUTION
C ****
      DO 110 IC=1,NR
         B(N,IC) = B(N,IC)/A(N,N)
 110  CONTINUE
      DO 130 KK=1,NM1
         I = N-KK
         IP1 = I+1
         DO 130 J=1,NR
            SUM = B(I,J)
            DO 120 K=IP1,N
               SUM = SUM - A(I,K)*B(K,J)
 120        CONTINUE
            B(I,J) = SUM/A(I,I)
 130  CONTINUE
      GO TO 170
 140  CONTINUE
      IF (A(1,1) .EQ. 0.D0) GO TO 160
      DO 150 J=1,NR
         B(1,J) = B(1,J)/A(1,1)
 150  CONTINUE
      GO TO 170
 160  CONTINUE
      IERR = 5
 170  CONTINUE
      RETURN
      END