      SUBROUTINE DEC (N, NDIM, A, IP, IER)
C-----------------------------------------------------------------------
C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C INPUT..
C    N = ORDER OF MATRIX.
C    NDIM = DECLARED DIMENSION OF ARRAY  A .
C    A = MATRIX TO BE TRIANGULARIZED.
C OUTPUT..
C    A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
C    A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
C    IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
C    IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
C    IER = 0 IF A NONSINGULAR, OR K IF A FOUND TO BE
C          SINGULAR AT STAGE K.
C USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
C DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
C IF IP(N)=0, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
C INTERCHANGES FINISHED IN U , ONLY PARTLY IN L .
C
C REFERENCE..
C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
C COMM. ASSOC. COMPUT. MACH., 15 (1972), P. 274.
C-----------------------------------------------------------------------
      INTEGER IER, IP, N, NDIM
      INTEGER I, J, K, KP1, M, NM1
      REAL A
      REAL T
      REAL ONE, ZERO
      DIMENSION A(NDIM,N),IP(N)
      DATA ONE /1.0E0/, ZERO /0.0E0/
C
      IER = 0
      IP(N) = 1
      IF (N .EQ. 1) GO TO 70
      NM1 = N - 1
      DO 60 K = 1,NM1
        KP1 = K + 1
C  FIND THE PIVOT IN COLUMN K.  SEARCH ROWS K TO N. --------------------
        M = K
        DO 10 I = KP1,N
  10      IF (ABS(A(I,K)) .GT. ABS(A(M,K))) M = I
        IP(K) = M
C  INTERCHANGE ELEMENTS IN ROWS K AND M. -------------------------------
        T = A(M,K)
        IF (M .EQ. K) GO TO 20
        IP(N) = -IP(N)
        A(M,K) = A(K,K)
        A(K,K) = T
 20     IF (T .EQ. ZERO) GO TO 80
C  STORE MULTIPLIERS IN A(I,K), I = K+1,...,N. -------------------------
        T = ONE/T
        DO 30 I = KP1,N
 30       A(I,K) = -A(I,K)*T
C  APPLY MULTIPLIERS TO OTHER COLUMNS OF A. ----------------------------
        DO 50 J = KP1,N
          T = A(M,J)
          A(M,J) = A(K,J)
          A(K,J) = T
          IF (T .EQ. ZERO) GO TO 50
          DO 40 I = KP1,N
 40         A(I,J) = A(I,J) + A(I,K)*T
 50       CONTINUE
 60     CONTINUE
 70   K = N
      IF (A(N,N) .EQ. ZERO) GO TO 80
      RETURN
 80   IER = K
      IP(N) = 0
      RETURN
C-----------------------  END OF SUBROUTINE DEC  -----------------------
      END
