c Updated 10/24/2001.
c
ccccccccccccccccccccccccc     Program 4.3     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
      PROGRAM  EX43
C
C An example of solving linear equation set A(N,N)*X(N) = B(N)
C with the partial-pivoting Gaussian elimination scheme.  The
C numerical values are for the Wheatstone bridge example discussed
C in Section 4.1 in the book with all resistances being 100 ohms
C and the voltage 200 volts.
C
      PARAMETER (N=3)
      DIMENSION X(N),B(N),A(N,N),INDX(N)
      DATA B/200.0,0.0,0.0/,
     *     ((A(I,J), J=1,N),I=1,N) /100.0,100.0,100.0,-100.0,
     *                   300.0,-100.0,-100.0,-100.0, 300.0/
C
      CALL LEGS (A,N,B,X,INDX)
C 
      WRITE (6, 999) (X(I), I=1,N)
      STOP
  999 FORMAT (F16.8)
      END
C
      SUBROUTINE LEGS(A,N,B,X,INDX)
C
C Subroutine to solve the equation A(N,N)*X(N) = B(N) with the
C partial-pivoting Gaussian elimination scheme.
C
      DIMENSION A(N,N),B(N),X(N),INDX(N)
C
      CALL ELGS(A,N,INDX)
C
      DO     100 I = 1, N-1
        DO    90 J = I+1, N
            B(INDX(J)) = B(INDX(J))
     *                  -A(INDX(J),I)*B(INDX(I))
   90   CONTINUE
  100 CONTINUE
C
      X(N) = B(INDX(N))/A(INDX(N),N)
      DO     200 I = N-1, 1, -1
        X(I) = B(INDX(I))
        DO   190 J = I+1, N
          X(I) = X(I)-A(INDX(I),J)*X(J)
  190   CONTINUE
          X(I) =  X(I)/A(INDX(I),I)
  200 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE ELGS(A,N,INDX)
C
C Subroutine to perform the partial-pivoting Gaussian elimination.
C A(N,N) is the original matrix in the input and transformed
C matrix plus the pivoting element ratios below the diagonal in
C the output.  INDX(N) records the pivoting order.

C
      DIMENSION A(N,N),INDX(N),C(N)
C
C Initialize the index
C
      DO     50    I = 1, N
        INDX(I) = I
   50 CONTINUE
C
C Find the rescaling factors, one from each row
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
C
C Search the pivoting (largest) element from each column
C
      DO     200   J = 1, N-1
        PI1 = 0.0
        DO   150   I = J, N
          PI = ABS(A(INDX(I),J))/C(INDX(I))
          IF (PI.GT.PI1) THEN
            PI1 = PI
            K   = I
          ELSE
          ENDIF
  150   CONTINUE
C
C Interchange the rows via INDX(N) to record pivoting order
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)
C
C Record pivoting ratios below the diagonal
C
          A(INDX(I),J) = PJ
C
C Modify other elements accordingly
C
          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
      RETURN
      END