      SUBROUTINE SPLU (A,IA,JA,R,C,IC,N,MAX,D,T,IT,JT,IU,W,P,IERR)
C ----------------------------------------------------------------------
C     SPLU EMPLOYS GAUSSIAN ELIMINATION WITH COLUMN INTERCHANGES TO
C     PERFORM THE LU DECOMPOSITION OF A SPARSE MATRIX. U IS A UNIT
C     UPPER TRIANGULAR MATRIX.
C
C
C     INPUT ARGUMENTS ---
C
C     A,IA,JA  THE SPARSE MATRIX TO BE DECOMPOSED.
C
C     R        INTEGER ARRAY OF N ENTRIES SPECIFYING THE ORDER OF
C              THE ROWS OF A.
C
C     C        INTEGER ARRAY OF N ENTRIES SPECIFYING A SUGGESTED
C              ORDER OF THE COLUMNS. C IS ALSO AN OUTPUT ARGUMENT.
C
C     IC       INTEGER ARRAY OF N ENTRIES WHICH IS THE INVERSE OF C
C              (I.E., IC(C(I)) = I). IC IS ALSO AN OUTPUT ARGUMENT.
C
C     N        ORDER OF THE MATRIX A.
C
C     MAX      INTEGER SPECIFYING THE MAXIMUM NUMBER OF OFF-DIAGONAL
C              NONZERO ENTRIES OF L AND U WHICH MAY BE STORED.
C
C
C     OUTPUT ARGUMENTS ---
C
C     C        INTEGER ARRAY OF N ENTRIES SPECIFYING THE ORDER OF
C              THE COLUMNS THAT WAS SELECTED BY THE ROUTINE.
C
C     IC       INTEGER ARRAY OF N ENTRIES WHICH IS THE INVERSE OF C.
C
C     D        REAL ARRAY CONTAINING THE N DIAGONAL ELEMENTS OF L.
C
C     T,IT,IU  T CONTAINS THE OFF-DIAGONAL NONZERO ELEMENTS OF L AND
C              U. FOR I = 1,...,N THE OFF-DIAGONAL NONZERO ELEMENTS
C              OF THE I-TH ROW OF L ARE STORED IN LOCATIONS
C              IT(I),...,IU(I)-1 OF T, AND THE OFF-DIAGONAL NONZERO
C              ELEMENTS OF THE I-TH ROW OF U ARE STORED IN LOCATIONS
C              IU(I),...,IT(I+1)-1 OF T.
C
C     JT       INTEGER ARRAY CONTAINING THE COLUMN INDICES (ACCORDING
C              TO THE ORGINAL COLUMN ORDERING) OF THE ELEMENTS OF T
C              (I.E., FOR EACH L(I,J) AND U(I,J) IN T, C(J) IS THE
C              CORRESPONDING COLUMN INDEX IN JT).
C
C     IERR     INTEGER SPECIFYING THE STATUS OF THE RESULTS. IF THE
C              LU DECOMPOSITION IS OBTAINED THEN IERR = THE NUMBER
C              OF OFF-DIAGONAL ENTRIES OF L AND U WHICH WERE STORED
C              IN T. OTHERWISE IERR IS ASSIGNED A NEGATIVE VALUE.
C
C
C     WORK SPACES ---
C
C     W        REAL ARRAY OF DIMENSION N.
C
C     P        INTEGER ARRAY OF DIMENSION N+1.
C ----------------------------------------------------------------------
      DIMENSION A(1), IA(1), JA(1)
      INTEGER   R(N), C(N), IC(N)
      DIMENSION T(MAX), IT(1), JT(MAX), IU(N)
      DIMENSION D(N), W(N)
      INTEGER   P(1), PM
      DATA ZERO/0.0/
C
      JPTR = 0
      IT(1) = 1
      DO 10 J = 1,N
   10 W(J) = ZERO
C
C     PERFORM THE LU FACTORIZATION OF THE R(K)-TH ROW OF A
C
      DO 100 K = 1,N
      LK = R(K)
      JMIN = IA(LK)
      JMAX = IA(LK+1) - 1
      IF (JMIN .GT. JMAX) GO TO 200
C
C     SET P TO THE REORDERED ROW OF A
C
      P(N+1) = N + 1
      JJ = JMAX
   20 LJ = JA(JJ)
      J = IC(LJ)
      W(J) = A(JJ)
      PM = N + 1
   21 M = PM
      PM = P(M)
      IF (PM - J) 21,210,22
   22 P(M) = J
      P(J) = PM
      JJ = JJ - 1
      IF (JJ .GE. JMIN) GO TO 20
C
C     PROCESS THE ENTRIES IN THE LOWER TRIANGLE OF A
C
      I = N + 1
   30 I = P(I)
      IF (I .GE. K) GO TO 50
      IF (W(I) .EQ. ZERO) GO TO 30
C
C     L(K,I) IS NONZERO. THEREFORE STORE IT IN L.
C
      JPTR = JPTR + 1
      IF (JPTR .GT. MAX) GO TO 230
      CONST = W(I)
      T(JPTR) = CONST
      JT(JPTR) = C(I)
      W(I) = ZERO
C
C     PERFORM ELIMINATION USING THE I-TH ROW OF U
C
      JMIN = IU(I)
      JMAX = IT(I+1) - 1
      IF (JMIN .GT. JMAX) GO TO 30
      PM = I
      DO 43 JJ = JMIN,JMAX
      LJ = JT(JJ)
      J = IC(LJ)
      IF (W(J) .NE. ZERO) GO TO 43
      IF (J - PM) 40,43,41
   40 PM = I
   41 M = PM
      PM = P(M)
      IF (PM - J) 41,43,42
   42 P(M) = J
      P(J) = PM
      PM = J
   43 W(J) = W(J) - CONST*T(JJ)
      GO TO 30
C
C     SEARCH FOR THE K-TH PIVOT ELEMENT
C
   50 IF (I .GT. N) GO TO 220
      WMAX = ABS(W(I))
      MAXI = I
      PM = I
   51 M = PM
      PM = P(M)
      IF (PM .GT. N) GO TO 60
      WI = ABS(W(PM))
      IF (WI .LE. WMAX) GO TO 51
      WMAX = WI
      MAXI = PM
      MAXIL = M
      GO TO 51
C
C     STORE THE PIVOT IN D
C
   60 IF (WMAX .EQ. ZERO) GO TO 220
      D(K) = W(MAXI)
C
C     PERFORM THE COLUMN INTERCHANGE
C
      IF (I .EQ. K) GO TO 70
      IF (I .EQ. MAXI) GO TO 70
      P(MAXIL) = P(MAXI)
      GO TO 80
   70 I = P(I)
C
   80 W(MAXI) = W(K)
      W(K) = ZERO
      LK = C(K)
      LL = C(MAXI)
      C(K) = LL
      C(MAXI) = LK
      IC(LK) = MAXI
      IC(LL) = K
C
C     THE REMAINING ELEMENTS OF P FORM THE K-TH ROW OF U
C
      IU(K) = JPTR + 1
   90 IF (I .GT. N) GO TO 100
      IF (W(I) .EQ. ZERO) GO TO 91
      JPTR = JPTR + 1
      IF (JPTR .GT. MAX) GO TO 230
      T(JPTR) = W(I)/D(K)
      JT(JPTR) = C(I)
      W(I) = ZERO
   91 I = P(I)
      GO TO 90
C
C     PREPARE FOR THE NEXT ROW
C
  100 IT(K+1) = JPTR + 1
C
      IERR = JPTR
      RETURN
C
C     -------------------- ERROR RETURN --------------------
C
C     ROW R(K) IS NULL
C
  200 IERR = -K
      RETURN
C
C     ROW R(K) HAS A DUPLICATE ENTRY
C
  210 IERR = -(N + K)
      RETURN
C
C     ZERO PIVOT IN ROW R(K)
C
  220 IERR = -(2*N + K)
      RETURN
C
C     STORAGE FOR L AND U EXCEEDED ON ROW R(K)
C
  230 IERR = -(3*N + K)
      RETURN
      END
