*
      SUBROUTINE GLU1( A, N, N1, EPS, WK, IP, IER )
*
*        GLU1
*               COPYRIGHT : H.HASEGAWA, OCT.  4 1991 V.1
*
*               MODIFIED  : S.TAKEHIRO, JUL.  9 1993 (R*4)
*
*               SOLVES SIMULTANEOUS LINEAR EQUATIONS
*               BY GAUSSIAN ELIMINATION METHOD.
*
*  [INPUT]
      INTEGER N             !" ORDER OF MATRIX.
      INTEGER N1            !" SIZE OF ARRAY A.
*
      DOUBLE PRECISION  EPS !" PARAMETER TO CHECK SINGULARITY OF THE
*                           !" MATRIX. ( STANDARD VALUE 3.52D-15 )
*
*  [MODIFY]
      DOUBLE PRECISION  A(N1,*)  !" 2-DIM. ARRAY CONTAINING THE COEFFICIENTS.
*                                !" RESULT OF GAUSSIAN ELIMINATION.
*  [OUTPUT]
      INTEGER IP(*)              ! PIVOT NUMBER.
      INTEGER IER                !" = 0,  FOR NORMAL EXECUTION.
*                                !" = 1,  FOR SINGULAR MATRIX.
*                                !" = 3,  FOR INVALID ARGUEMENT.
*  [WORK]
      DOUBLE PRECISION  WK(*)    !"1-DIM. ARRAY.
*
*  [INTERNAL WORK]
      INTEGER I, J, K, IPK
      DOUBLE PRECISION    AMAX, AIK, W, T
*
*             LEFT HAND SIDE
      IF( EPS.LT.0.0D0 )  EPS = 3.52D-15
      IF( ( N.GT.N1 ).OR.( N.LE.0 ) )  THEN
         IER = 3
         WRITE(*,*) '  (SUBR. GLU1)  INVALID ARGUMENT.  N1, N =', N1, N
         RETURN
      END IF
*
      IER = 0
      DO 100 K = 1, N
*             FIND MAXIMUM ELEMENT IN THE K-TH COLUMN.
        AMAX = ABS(A(K,K))
        IPK = K
        DO 110 I = K+1, N
          AIK = ABS(A(I,K))
          IF( AIK.GT.AMAX )  THEN
             IPK = I
             AMAX = AIK
          END IF
  110   CONTINUE
        IP(K) = IPK
*
        IF( AMAX.GT.EPS )  THEN
           IF( IPK.NE.K )  THEN
              W = A(IPK,K)
              A(IPK,K) = A(K,K)
              A(K,K) = W
           END IF
*             COMPUTE ALFA
           DO 120 I = K+1, N
             A(I,K) = -A(I,K)/A(K,K)
  120        WK(I) = A(I,K)
*
           DO 130 J = K+1, N
             IF( IPK.NE.K )  THEN
                W = A(IPK,J)
                A(IPK,J) = A(K,J)
                A(K,J) = W
             END IF
*             GAUSSIAN ELIMINATION
             T = A(K,J)
             DO 140 I = K+1, N
  140          A(I,J) = A(I,J)+WK(I)*T
  130      CONTINUE
*             MATRIX IS SINGULAR.
        ELSE
           IER = 1
           IP(K) = K
           DO 150 I = K+1, N
  150        A(I,K) = 0.0D0
           WRITE(*,*) '  (SUBR. GLU1)  MATRIX IS SINGULAR AT K =', K
           RETURN
        END IF
  100 CONTINUE
      RETURN
      END
***********************************************************************
      SUBROUTINE GSLV1( A, N, N1, B, IP )
*
*        GSLV1
*               COPYRIGHT : H.HASEGAWA, OCT.  4 1991 V.1
*
*               MODIFIED  : S.TAKEHIRO, JUL.  9 1993 (R*4)
*
*               SOLVES SIMULTANEOUS LINEAR EQUATIONS
*               BY GAUSSIAN ELIMINATION METHOD.
*
*  [INPUT]
      INTEGER N                  !" ORDER OF MATRIX.
      INTEGER N1                 !" SIZE OF ARRAY A.
      DOUBLE PRECISION  A(N1,*)  !" RESULT OF GAUSSIAN ELIMINATION.
      INTEGER IP(*)              !" PIVOT NUMBER.
*
*  [MODIFY]
      DOUBLE PRECISION  B(*)     !" 1-DIM. ARRAY CONTAINING THE RIGHT HAND
*                                !" SIDE VECTOR. OUTPUT IS SOLUTION.
*
*  [INTERNAL WORK]
      INTEGER I, K
      DOUBLE PRECISION  T, W
*
*             FORWARD ELIMINATION PROCESS
      DO 100 K = 1, N
        IF( IP(K).NE.K ) THEN
           W = B(IP(K))
           B(IP(K)) = B(K)
           B(K) = W
        END IF
*             GAUSSIAN ELIMINATION
        T = B(K)
        DO 110 I = K+1, N
  110     B(I) = B(I)+A(I,K)*T
  100 CONTINUE
*             BACKWARD SUBSTITUTION PROCESS
      B(N) = B(N)/A(N,N)
      DO 200 K = N-1, 1, -1
        T = B(K+1)
        DO 210 I = 1, K
  210     B(I) = B(I)-A(I,K+1)*T
        B(K) = B(K)/A(K,K)
  200 CONTINUE
      RETURN
      END
