      SUBROUTINE CROUT1(A,KA,N,IEND,INDEX,TEMP,JP,IERR)
C     ******************************************************************
C     CROUT PROCEDURE FOR INVERTING MATRICES
C     ******************************************************************
C     A IS A MATRIX OF ORDER N WHERE N IS GREATER THAN OR EQUAL TO 1.
C     THE INVERSE OF A IS COMPUTED AND STORED IN A.
C
C     KA = LENGTH OF THE COLUMNS OF THE ARRAY A
C     JP = THE NUMBER OF THE COLUMN THAT CONTAINS THE SMALLEST PIVOT
C
C     IEND MAY BE 0,1,...,N-1.  IT IS ASSUMED THAT EACH OF THE FIRST
C     IEND COLUMNS OF THE MATRIX A CONTAINS ONLY ONE NONZERO ELEMENT
C     AND THAT THE NONZERO ELEMENT IS 1 OR -1.
C
C     INDEX IS AN ARRAY OF DIMENSION N-1 OR LARGER THAT IS USED BY THE
C     ROUTINE FOR KEEPING TRACK OF THE ROW INTERCHANGES THAT ARE MADE.
C
C     TEMP IS A TEMPORARY STORAGE ARRAY.
C
C     IERR REPORTS THE STATUS OF THE RESULTS. IF A IS NONSINGULAR THEN
C     THE INVERSE OF A IS COMPUTED AND IERR=0. OTHERWISE IF A IS FOUND
C     TO BE SINGULAR THEN IERR=1 AND THE ROUTINE TERMINATES.
C     --------------------
      DIMENSION A(1),INDEX(1),TEMP(N)
      DOUBLE PRECISION DSUM
      MAX = KA*N
      MCOL = IEND*KA
      IF (IEND .EQ. 0) GO TO 100
C
C           PROCESS THE FIRST IEND COLUMNS OF A
C
      KCOL = 0
      DO 32 K = 1,IEND
      KK = KCOL + K
      NK = KCOL + N
      DO 10 LK = KK,NK
      IF (A(LK)) 20,10,30
   10 CONTINUE
      JP = K
      GO TO 300
C
   20 L = LK - KCOL
      LJ0 = MCOL + L
      DO 21 LJ = LJ0,MAX,KA
   21 A(LJ) = -A(LJ)
C
   30 L = LK - KCOL
      INDEX(K) = L
      IF (K .EQ. L) GO TO 32
      LJ = LK
      DO 31 KJ = KK,MAX,KA
      C = A(KJ)
      A(KJ) = A(LJ)
      A(LJ) = C
   31 LJ = LJ + KA
   32 KCOL = KCOL + KA
C
C           PROCESS THE REMAINING COLUMNS OF A
C
  100 NM1 = N - 1
      JP = 1
      IERR = 0
      PMIN = 0.0
      IBEG = IEND + 1
      IF (IBEG .EQ. N) GO TO 190
C
      K = IBEG
      KM1 = IEND
      KP1 = K + 1
      KCOL = MCOL
      KK = KCOL + K
      DO 172 KCOUNT = IBEG,NM1
C
C     SEARCH FOR THE K-TH PIVOT ELEMENT (K=IBEG,...,N-1)
C
      L = K
      S = ABS(A(KK))
      DO 110 I = KP1,N
      IK = KCOL + I
      C = ABS(A(IK))
      IF (S .GE. C) GO TO 110
      L = I
      S = C
  110 CONTINUE
C
      IF (K.GT.IBEG .AND. S.GE.PMIN) GO TO 120
      JP = K
      PMIN = S
      IF (S .EQ. 0.0) GO TO 300
C
C              INTERCHANGING ROWS K AND L
C
  120 INDEX(K) = L
      IF (K .EQ. L) GO TO 130
      KJ0 = MCOL + K
      LJ  = MCOL + L
      DO 121 KJ = KJ0,MAX,KA
      C = A(KJ)
      A(KJ) = A(LJ)
      A(LJ) = C
  121 LJ = LJ + KA
C
C       COMPUTE THE K-TH ROW OF U (K=IBEG,...,N-1)
C
  130 C = A(KK)
      IF (K .GT. IBEG) GO TO 140
      KJ0 = KK + KA
      DO 131 KJ = KJ0,MAX,KA
  131 A(KJ) = A(KJ)/C
      GO TO 160
C
  140 KL = MCOL + K
      DO 141 L = IBEG,KM1
      TEMP(L) = A(KL)
  141 KL = KL + KA
C
      KJ0 = KK + KA
      DO 151 KJ = KJ0,MAX,KA
      JCOL = KJ - K
      DSUM = -A(KJ)
        DO 150 L = IBEG,KM1
        LJ = JCOL + L
  150   DSUM = DSUM + DBLE(TEMP(L))*DBLE(A(LJ))
  151 A(KJ) = SNGL(-DSUM)/C
C
C      COMPUTE THE K-TH COLUMN OF L (K=IBEG+1,...,N)
C
  160 KM1 = K
      K = KP1
      KP1 = K + 1
      KCOL = KCOL + KA
      KK = KCOL + K
      DO 161 L = IBEG,KM1
      LK = KCOL + L
  161 TEMP(L) = A(LK)
C
      DO 171 I = K,N
      IL = MCOL + I
      DSUM = 0.D0
        DO 170 L = IBEG,KM1
        DSUM = DSUM + DBLE(A(IL))*DBLE(TEMP(L))
  170   IL = IL + KA
  171 A(IL) = DBLE(A(IL)) - DSUM
  172 CONTINUE
C
C           CHECK THE N-TH PIVOT ELEMENT
C
  190 NCOL = MAX - KA
      NN = NCOL + N
      C = ABS(A(NN))
      IF (C .GT. PMIN) GO TO 200
      JP = N
      IF (C .EQ. 0.0) GO TO 300
C
C          REPLACE L WITH THE INVERSE OF L
C
  200 IF (IBEG .EQ. N) GO TO 213
      JJ = MCOL + IBEG
      I = KA + 1
      DO 212 J = IBEG,NM1
      A(JJ) = 1.0/A(JJ)
      TEMP(J) = A(JJ)
      KJ = JJ
        DO 211 KM1 = J,NM1
        K = KM1 + 1
        KJ = KJ + 1
        DSUM = 0.D0
        KL = KJ
          DO 210 L = J,KM1
          DSUM = DSUM + DBLE(A(KL)*TEMP(L))
  210     KL = KL + KA
        A(KJ) = SNGL(-DSUM)/A(KL)
  211   TEMP(K) = A(KJ)
  212 JJ = JJ + I
  213 A(NN) = 1.0/A(NN)
      IF (N .EQ. 1) RETURN
C
C       SOLVE UX = Y WHERE Y IS THE INVERSE OF L
C
      DO 242 NMK = 1,NM1
      K = N - NMK
      LMIN = MAX0(IBEG,K+1)
      KL = (LMIN-1)*KA + K
        DO 230 L = LMIN,N
        TEMP(L) = A(KL)
        A(KL) = 0.0
  230   KL = KL + KA
C
      KJ0 = MCOL + K
        DO 241 KJ = KJ0,MAX,KA
        DSUM = -A(KJ)
        LJ = (KJ - K) + LMIN
          DO 240 L = LMIN,N
          DSUM = DSUM + DBLE(TEMP(L)*A(LJ))
  240     LJ = LJ + 1
  241   A(KJ) = -DSUM
  242 CONTINUE
C
C                 COLUMN INTERCHANGES
C
      JCOL = NCOL - KA
      DO 251 NMJ = 1,NM1
      J = N - NMJ
      K = INDEX(J)
      IF (J .EQ. K) GO TO 251
      IJ = JCOL
      IK = (K-1)*KA
        DO 250 I = 1,N
        IJ = IJ + 1
        IK = IK + 1
        C = A(IJ)
        A(IJ) = A(IK)
  250   A(IK) = C
  251 JCOL = JCOL - KA
      RETURN
C
C                    ERROR RETURN
C
  300 IERR = 1
      RETURN
      END
