ccccccccccccccccccccccccc     Program 4.1     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE ELGS(A,N,IR,INDX)
C
C Subroutine for the partial pivoting Gaussian elimination.
C A is the coefficient matrix in the input and the transformed
C matrix in the output.  INDX records the pivoting order.  IR
C is the indicator for the rescaling.
C
      PARAMETER (NMAX=100)
      DIMENSION A(N,N),INDX(N),C(NMAX)
      LOGICAL IR
C
      IF (N.GT.NMAX) STOP 'The matrix dimension is too large.'
C
C Initialize the index
C
      DO     50    I = 1, N
        INDX(I) = I
   50 CONTINUE
C
C Rescale the coefficients
C
        DO     100   I = 1, N
          C1= 0.0
          DO    90   J = 1, N
            C1 = AMAX1(C1,ABS(A(I,J)))
   90     CONTINUE
          C(I) = C1
  100   CONTINUE
        DO     140   I = 1, N
          DO   130   J = 1, N
            A(I,J) = A(I,J)/C(I)
  130     CONTINUE
  140   CONTINUE
C
C Search the pivoting elements
C
      DO     200   J = 1, N-1
        PI1 = 0.0
        DO   150   I = J, N
          PI = ABS(A(INDX(I),J))
          IF (PI.GT.PI1) THEN
            PI1 = PI
            K   = I
          ELSE
          ENDIF
  150   CONTINUE
C
C Eliminate the elements below the diagonal
C
        ITMP    = INDX(J)
        INDX(J) = INDX(K)
        INDX(K) = ITMP
        DO   170   I = J+1, N
          PJ  = A(INDX(I),J)/A(INDX(J),J)
          A(INDX(I),J) = PJ
          DO 160   K = J+1, N
            A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
  160     CONTINUE
  170   CONTINUE
  200 CONTINUE
C
C Rescale the coefficients back to the original magnitudes
C
      IF(IR) THEN
        DO     240   I = 1, N
          DO   230   J = 1, N
            A(I,J) = C(I)*A(I,J)
  230     CONTINUE
  240   CONTINUE
      ELSE
      ENDIF
C
      RETURN
      END
