C
C PAGE 220-223: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
C
C FILE: GAUSS.FOR
C
C GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING (GAUSS,SOLVE,TSTGAUS)
C
      DIMENSION A1(4,4),A2(4,4),A3(4,4),B1(4),B2(4),B3(4)
      DIMENSION L(4),S(4),X(4)
      DATA ((A1(I,J),I=1,4),J=1,4)/3.,1.,6.,0.,4.,5.,3.,0.,3.,-1.,7.,
     A     0.,0.,0.,0.,0./
      DATA (B1(I),I=1,4)/16.,-12.,102.,0./
      DATA ((A2(I,J),I=1,4),J=1,4)/3.,2.,1.,0.,2.,-3.,4.,0.,-5.,1.,-1.,
     A     0.,0.,0.,0.,0./
      DATA (B2(I),I=1,4)/4.,8.,3.,0./
      DATA ((A3(I,J),I=1,4),J=1,4)/1.,3.,5.,4.,-1.,2.,8.,2.,2.,1.,6.,
     A 5.,1.,4.,3.,3./
      DATA (B3(I),I=1,4)/5.,8.,10.,12./
C
      CALL TSTGAUS(3,A1,4,L,S,B1,X)
      CALL TSTGAUS(3,A2,4,L,S,B2,X)
      CALL TSTGAUS(4,A3,4,L,S,B3,X)
      END

      SUBROUTINE TSTGAUS(N,A,IA,L,S,B,X)
      DIMENSION  A(IA,N),B(N),X(N),S(N),L(N)
      PRINT 10,((A(I,J),J=1,N),I=1,N)
      PRINT 10,(B(I),I=1,N)
      CALL GAUSS(N,A,IA,L,S)
      CALL SOLVE(N,A,IA,L,B,X)
      PRINT 10,(X(I),I=1,N)
      RETURN
10    FORMAT(5X,3(F10.5,2X))
      END

      SUBROUTINE GAUSS(N,A,IA,L,S)
      DIMENSION  A(IA,N),L(N),S(N)
      DO 3 I = 1,N
        L(I) = I
        SMAX = 0.0
        DO 2 J = 1,N
          SMAX = AMAX1(SMAX,ABS(A(I,J)))
   2    CONTINUE
        S(I) = SMAX
   3  CONTINUE
      DO 7 K = 1,N-1
        RMAX = 0.0
        DO 4 I = K,N
          R = ABS(A(L(I),K))/S(L(I))
          IF(R .LE. RMAX)  GO TO 4
          J = I
          RMAX = R
   4    CONTINUE
        LK = L(J)
        L(J) = L(K)
        L(K) = LK
        DO 6 I = K+1,N
          XMULT = A(L(I),K)/A(LK,K)
          DO 5 J = K+1,N
            A(L(I),J) = A(L(I),J) - XMULT*A(LK,J)
   5      CONTINUE
          A(L(I),K) = XMULT
   6    CONTINUE
   7  CONTINUE
      RETURN
      END

      SUBROUTINE SOLVE(N,A,IA,L,B,X)
      DIMENSION  A(IA,N),L(N),B(N),X(N)
      DO 3 K = 1,N-1
        DO 2 I = K+1,N
          B(L(I)) = B(L(I)) - A(L(I),K)*B(L(K))
   2    CONTINUE
   3  CONTINUE
      X(N) = B(L(N))/A(L(N),N)
      DO 5 I = N-1,1,-1
        SUM = B(L(I))
        DO 4 J = I+1,N
          SUM = SUM - A(L(I),J)*X(J)
   4    CONTINUE
        X(I) = SUM/A(L(I),I)
   5  CONTINUE
      RETURN
      END