      SUBROUTINE NSPIV1 (N,IA,JA,A,B,MAX,R,C,IC,X,Y,P,IU,JU,U,IERR)
C
C
C  NSPIV1 USES SPARSE GAUSSIAN ELIMINATION WITH
C  COLUMN INTERCHANGES TO SOLVE THE LINEAR SYSTEM A X = B.  THE
C  ELIMINATION PHASE PERFORMS ROW OPERATIONS ON A AND B TO OBTAIN
C  A UNIT UPPER TRIANGULAR MATRIX U AND A VECTOR Y.  THE SOLUTION
C  PHASE SOLVES U X = Y.
C
C
C  SEE SPSLV FOR DESCRIPTIONS OF ALL INPUT AND OUTPUT ARGUMENTS
C  OTHER THAN THOSE DESCRIBED BELOW
C
C  IC  INTEGER ARRAY OF N ENTRIES WHICH IS THE INVERSE OF C
C      (I.E., IC(C(I)) = I). IC IS BOTH AN INPUT AND OUTPUT
C      ARGUMENT.
C
C  INPUT ARGUMENTS (USED INTERNALLY ONLY)---
C
C  Y   REAL ARRAY OF N ENTRIES USED TO COMPUTE THE UPDATED
C      RIGHT HAND SIDE
C
C  P   INTEGER ARRAY OF N+1 ENTRIES USED FOR A LINKED LIST.
C      P(N+1) IS THE LIST HEADER, AND THE ENTRY FOLLOWING
C      P(K) IS IN P(P(K)).  THUS, P(N+1) IS THE FIRST DATA
C      ITEM, P(P(N+1)) IS THE SECOND, ETC.  A POINTER OF
C      N+1 MARKS THE END OF THE LIST
C
C  IU  INTEGER ARRAY OF N+1 ENTRIES USED FOR ROW POINTERS TO U
C      (SEE MATRIX STORAGE DESCRIPTION BELOW)
C
C  JU  INTEGER ARRAY OF MAX ENTRIES USED FOR COLUMN NUMBERS OF
C      THE NONZEROES IN THE STRICT UPPER TRIANGLE OF U.  (SEE
C      MATRIX STORAGE DESCRIPTION BELOW)
C
C  U   REAL ARRAY OF MAX ENTRIES USED FOR THE ACTUAL NONZEROES IN
C      THE STRICT UPPER TRIANGLE OF U.  (SEE MATRIX STORAGE
C      DESCRIPTION BELOW)
C
C
C  STORAGE OF SPARSE MATRICES---
C
C  THE SPARSE MATRIX A IS STORED USING THREE ARRAYS IA, JA, AND A.
C  THE ARRAY A CONTAINS THE NONZEROES OF THE MATRIX ROW-BY-ROW, NOT
C  NECESSARILY IN ORDER OF INCREASING COLUMN NUMBER.  THE ARRAY JA
C  CONTAINS THE COLUMN NUMBERS CORRESPONDING TO THE NONZEROES STORED
C  IN THE ARRAY A (I.E., IF THE NONZERO STORED IN A(K) IS IN
C  COLUMN J, THEN JA(K) = J).  THE ARRAY IA CONTAINS POINTERS TO THE
C  ROWS OF NONZEROES/COLUMN INDICES IN THE ARRAY A/JA (I.E.,
C  A(IA(I))/JA(IA(I)) IS THE FIRST ENTRY FOR ROW I IN THE ARRAY A/JA).
C  IA(N+1) IS SET SO THAT IA(N+1) - IA(1) = THE NUMBER OF NONZEROES IN
C  A.  IU, JU, AND U ARE USED IN A SIMILAR WAY TO STORE THE STRICT UPPER
C  TRIANGLE OF U, EXCEPT THAT JU ACTUALLY CONTAINS C(J) INSTEAD OF J
C
C
      REAL A(1),B(1),U(1),X(1),Y(1)
      REAL DK,LKI,ONE,XPV,XPVMAX,YK,ZERO
      INTEGER C(1),IA(1),IC(1),IU(1),JA(1),JU(1),P(1),R(1)
      INTEGER CK,PK,PPK,PV,V,VI,VJ,VK
C
C
      ONE = 1.0
      ZERO = 0.0
C
C  INITIALIZE WORK STORAGE AND POINTERS TO JU
C
      DO 10 J=1,N
        X(J) = ZERO
 10     CONTINUE
      IU(1) = 1
      JUPTR = 0
C
C  PERFORM SYMBOLIC AND NUMERIC FACTORIZATION ROW BY ROW
C  VK (VI,VJ) IS THE GRAPH VERTEX FOR ROW K (I,J) OF U
C
      DO 170 K=1,N
C
C  INITIALIZE LINKED LIST AND FREE STORAGE FOR THIS ROW
C  THE R(K)-TH ROW OF A BECOMES THE K-TH ROW OF U.
C
        P(N+1) = N+1
        VK = R(K)
C
C  SET UP ADJACENCY LIST FOR VK, ORDERED IN
C  CURRENT COLUMN ORDER OF U.  THE LOOP INDEX
C  GOES DOWNWARD TO EXPLOIT ANY COLUMNS
C  FROM A IN CORRECT RELATIVE ORDER
C
        JMIN = IA(VK)
        JMAX = IA(VK+1) - 1
        IF (JMIN .GT. JMAX)  GO TO 1002
        J = JMAX
 20       JAJ = JA(J)
          VJ = IC(JAJ)
C
C  STORE A(K,J) IN WORK VECTOR
C
          X(VJ) = A(J)
C  THIS CODE INSERTS VJ INTO ADJACENCY LIST OF VK
          PPK = N+1
 30       PK = PPK
          PPK = P(PK)
          IF (PPK - VJ)  30,1003,40
 40       P(VJ) = PPK
          P(PK) = VJ
          J = J - 1
          IF (J .GE. JMIN) GO TO 20
C
C  THE FOLLOWING CODE COMPUTES THE K-TH ROW OF U
C
        VI = N+1
        YK = B(VK)
 50     VI = P(VI)
        IF (VI .GE. K) GO TO 110
C
C  VI LT VK -- PROCESS THE L(K,I) ELEMENT AND MERGE THE
C  ADJACENCY OF VI WITH THE ORDERED ADJACENCY OF VK
C
        LKI = - X(VI)
        X(VI) = ZERO
C
C  ADJUST RIGHT HAND SIDE TO REFLECT ELIMINATION
C
        YK = YK + LKI * Y(VI)
        PPK = VI
        JMIN = IU(VI)
        JMAX = IU(VI+1) - 1
        IF (JMIN .GT. JMAX)  GO TO 50
        DO 100 J=JMIN,JMAX
          JUJ = JU(J)
          VJ = IC(JUJ)
C
C  IF VJ IS ALREADY IN THE ADJACENCY OF VK,
C  SKIP THE INSERTION
C
          IF (X(VJ) .NE. ZERO)  GO TO 90
C
C  INSERT VJ IN ADJACENCY LIST OF VK.
C  RESET PPK TO VI IF WE HAVE PASSED THE CORRECT
C  INSERTION SPOT.  (THIS HAPPENS WHEN THE ADJACENCY OF
C  VI IS NOT IN CURRENT COLUMN ORDER DUE TO PIVOTING.)
C
          IF (VJ - PPK) 60,90,70
 60       PPK = VI
 70       PK = PPK
          PPK = P(PK)
          IF (PPK - VJ)  70,90,80
 80       P(VJ) = PPK
          P(PK) = VJ
          PPK = VJ
C
C  COMPUTE L(K,J) = L(K,J) - L(K,I)*U(I,J) FOR L(K,I) NONZERO
C  COMPUTE U*(K,J) = U*(K,J) - L(K,I)*U(I,J) FOR U(K,J) NONZERO
C  (U*(K,J) = U(K,J)*D(K,K))
C
 90       X(VJ) = X(VJ) + LKI * U(J)
 100      CONTINUE
        GO TO 50
C
C  PIVOT--INTERCHANGE LARGEST ENTRY OF K-TH ROW OF U WITH
C  THE DIAGONAL ENTRY.
C
C  FIND LARGEST ENTRY, COUNTING OFF-DIAGONAL NONZEROES
C
 110    IF (VI .GT. N) GO TO 1004
        XPVMAX = ABS(X(VI))
        MAXC = VI
        NZCNT = 0
        PV = VI
 120      V = PV
          PV = P(PV)
          IF (PV .GT. N) GO TO 130
          NZCNT = NZCNT + 1
          XPV = ABS(X(PV))
          IF (XPV .LE. XPVMAX) GO TO 120
          XPVMAX = XPV
          MAXC = PV
          MAXCL = V
          GO TO 120
 130    IF (XPVMAX .EQ. ZERO) GO TO 1004
C
C  IF VI = K, THEN THERE IS AN ENTRY FOR DIAGONAL
C  WHICH MUST BE DELETED.  OTHERWISE, DELETE THE
C  ENTRY WHICH WILL BECOME THE DIAGONAL ENTRY
C
        IF (VI .EQ. K) GO TO 140
        IF (VI .EQ. MAXC) GO TO 140
        P(MAXCL) = P(MAXC)
        GO TO 150
 140    VI = P(VI)
C
C  COMPUTE D(K) = 1/L(K,K) AND PERFORM INTERCHANGE.
C
 150    DK = ONE / X(MAXC)
        X(MAXC) = X(K)
        I = C(K)
        C(K) = C(MAXC)
        C(MAXC) = I
        CK = C(K)
        IC(CK) = K
        IC(I) = MAXC
        X(K) = ZERO
C
C  UPDATE RIGHT HAND SIDE.
C
        Y(K) = YK * DK
C
C  COMPUTE VALUE FOR IU(K+1) AND CHECK FOR STORAGE OVERFLOW
C
        IU(K+1) = IU(K) + NZCNT
        IF (IU(K+1) .GT. MAX+1) GO TO 1005
C
C  MOVE COLUMN INDICES FROM LINKED LIST TO JU.
C  COLUMNS ARE STORED IN CURRENT ORDER WITH ORIGINAL
C  COLUMN NUMBER (C(J)) STORED FOR CURRENT COLUMN J
C
        IF (VI .GT. N)  GO TO 170
        J = VI
 160      JUPTR = JUPTR + 1
          JU(JUPTR) = C(J)
          U(JUPTR) = X(J) * DK
          X(J) = ZERO
          J = P(J)
          IF (J .LE. N) GO TO 160
 170    CONTINUE
C
C  BACKSOLVE U X = Y, AND REORDER X TO CORRESPOND WITH A
C
      K = N
      DO 200 I=1,N
        YK = Y(K)
        JMIN = IU(K)
        JMAX = IU(K+1) - 1
        IF (JMIN .GT. JMAX)  GO TO 190
        DO 180 J=JMIN,JMAX
          JUJ = JU(J)
          JUJ = IC(JUJ)
          YK = YK - U(J) * Y(JUJ)
 180      CONTINUE
 190    Y(K) = YK
        CK = C(K)
        X(CK) = YK
        K = K-1
 200    CONTINUE
C
C  RETURN WITH IERR = NUMBER OF OFF-DIAGONAL NONZEROES IN U
C
      IERR = IU(N+1) - IU(1)
      RETURN
C
C  ERROR RETURNS
C
C  ROW K OF A IS NULL
C
 1002 IERR = -K
      RETURN
C
C  ROW K OF A HAS A DUPLICATE ENTRY
C
 1003 IERR = -(N+K)
      RETURN
C
C  ZERO PIVOT IN ROW K
C
 1004 IERR = -(2*N+K)
      RETURN
C
C  STORAGE FOR U EXCEEDED ON ROW K
C
 1005 IERR = -(3*N+K)
      RETURN
      END
