      SUBROUTINE RSLV (M0,N,A,IA,JA,B,R,C,MAX,X,IWK,WK,IERR)
C ----------------------------------------------------------------------
C     RSLV EMPLOYS GAUSSIAN ELIMINATION WITH COLUMN INTERCHANGES TO
C     SOLVE THE NXN LINEAR SYSTEM AX = B. THE ARGUMENT M0 SPECIFIES
C     IF RSLV IS BEING CALLED FOR THE FIRST TIME, OR IF IT IS BEING
C     RECALLED WHERE A IS THE SAME MATRIX BUT B HAS BEEN MODIFIED.
C     ON AN INITIAL CALL TO THE ROUTINE (WHEN M0=0) THE LU DECOMPO-
C     SITION OF A IS OBTAINED WHERE U IS A UNIT UPPER TRIANGULAR
C     MATRIX. THEN THE EQUATIONS ARE SOLVED. ON SUBSEQUENT CALLS
C     (WHEN M0.NE.0) THE EQUATIONS ARE SOLVED USING THE DECOMPOSITION
C     OBTAINED ON THE INITIAL CALL TO RSLV.
C
C
C     INPUT ARGUMENTS WHEN M0=0 ---
C
C     N        NUMBER OF EQUATIONS AND UNKNOWNS.
C
C     A,IA,JA  THE MATRIX A STORED IN SPARSE FORM.
C
C     B        ARRAY OF N ENTRIES CONTAINING THE RIGHT HAND SIDE DATA.
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     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 WHEN M0=0 ---
C
C     C        INTEGER ARRAY OF N ENTRIES SPECIFYING THE ORDER OF
C              THE COLUMNS THAT WAS SELECTED BY THE ROUTINE.
C
C     X        REAL ARRAY OF N ENTRIES CONTAINING THE SOLUTION.
C              B AND X MAY SHARE THR SAME STORAGE AREA.
C
C     IERR     INTEGER SPECIFYING THE STATUS OF THE RESULTS. IF THE
C              SOLUTION OF AX = B IS OBTAINED THEN IERR = MAX(1,M)
C              WHERE M IS THE TOTAL NUMBER OF OFF-DIAGONAL NONZERO
C              ENTRIES OF L AND U. OTHERWISE IERR.LE.0.
C
C
C     GENERAL STORAGE AREAS ---
C
C     IWK      INTEGER ARRAY OF DIMENSION 4*N + MAX + 2.
C
C     WK       REAL ARRAY OF DIMENSION 2*N + MAX.
C
C
C     AFTER AN INITIAL CALL TO RSLV, THE ROUTINE MAY BE RECALLED WITH
C     M0.NE.0 FOR A NEW B. WHEN M0.NE.0 IT IS ASSUMED THAT N,A,IA,JA,
C     R,C,IWK,WK HAVE NOT BEEN MODIFIED. THE ROUTINE RETRIEVES THE LU
C     DECOMPOSITION WHICH WAS OBTAINED ON THE INITIAL CALL TO RSLV
C     AND SOLVES THE NEW EQUATIONS AX = B. IN THIS CASE A,IA,JA,MAX,
C     AND IERR ARE NOT REFERENCED.
C ----------------------------------------------------------------------
      DIMENSION A(1), IA(1), JA(1), B(N), X(N)
      INTEGER   R(N), C(N), Y, T, P
      DIMENSION IWK(1), WK(1)
C
C     SET INDICES TO DIVIDE TEMPORARY STORAGE
C
      Y = N + 1
      T = Y + N
      P = N + 1
      IT = P + N + 1
      IU = IT + N + 1
      JT = IU + N
      IF (M0 .NE. 0) GO TO 20
C
C     COMPUTE THE INVERSE PERMUTATION OF C
C
      IERR = 0
      IF (N .LE. 0) RETURN
      DO 10 K = 1,N
      L = C(K)
   10 IWK(L) = K
C
C     OBTAIN THE LU DECOMPOSITION OF A
C
      CALL SPLU (A,IA,JA,R,C,IWK(1),N,MAX,WK(1),WK(T),IWK(IT),IWK(JT),
     *             IWK(IU),WK(Y),IWK(P),IERR)
      IF (IERR .LT. 0) RETURN
      IERR = MAX0(1,IERR)
C
C     SOLVE THE SYSTEM OF EQUATIONS
C
   20 CALL RSLV1 (N,R,C,IWK(1),WK(1),WK(T),IWK(IT),IWK(JT),IWK(IU),
     *                  B,X,WK(Y))
      RETURN
      END
