      SUBROUTINE NROC (N, IC, IA, JA, A, JAR, AR, P, FLAG)
CLLL. OPTIMIZE
C
C       ----------------------------------------------------------------
C
C               YALE SPARSE MATRIX PACKAGE - NONSYMMETRIC CODES
C                    SOLVING THE SYSTEM OF EQUATIONS MX = B 
C
C    I.   CALLING SEQUENCES
C         THE COEFFICIENT MATRIX CAN BE PROCESSED BY AN ORDERING ROUTINE
C    (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY) BEFORE USING
C    THE REMAINING SUBROUTINES.  IF NO REORDERING IS DONE, THEN SET
C    R(I) = C(I) = IC(I) = I  FOR I=1,...,N.  IF AN ORDERING SUBROUTINE
C    IS USED, THEN NROC SHOULD BE USED TO REORDER THE COEFFICIENT MATRIX
C    THE CALLING SEQUENCE IS --
C        (       (MATRIX ORDERING))
C        (NROC   (MATRIX REORDERING))
C         NSFC   (SYMBOLIC FACTORIZATION TO DETERMINE WHERE FILLIN WILL
C                  OCCUR DURING NUMERIC FACTORIZATION)
C         NNFC   (NUMERIC FACTORIZATION INTO PRODUCT LDU OF UNIT LOWER
C                  TRIANGULAR MATRIX L, DIAGONAL MATRIX D, AND UNIT
C                  UPPER TRIANGULAR MATRIX U, AND SOLUTION OF LINEAR
C                  SYSTEM)
C         NNSC   (SOLUTION OF LINEAR SYSTEM FOR ADDITIONAL RIGHT-HAND 
C                  SIDE USING LDU FACTORIZATION FROM NNFC)
C    (IF ONLY ONE SYSTEM OF EQUATIONS IS TO BE SOLVED, THEN THE
C    SUBROUTINE TRK SHOULD BE USED.)
C
C    II.  STORAGE OF SPARSE MATRICES
C         THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED
C    ROW-BY-ROW IN THE ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO
C    ENTRIES IN EACH ROW, WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY
C    LIES.  THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES
C    OF M ARE STORED IN THE ARRAY JA.  I.E., IF  A(K) = M(I,J),  THEN 
C    JA(K) = J.  IN ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND
C    HOW LONG IT IS.  THE INDEX POSITIONS IN JA AND A WHERE THE ROWS OF
C    M BEGIN ARE STORED IN THE ARRAY IA.  I.E., IF M(I,J) IS THE FIRST
C    (LEFTMOST) ENTRY IN THE I-TH ROW AND  A(K) = M(I,J),  THEN
C    IA(I) = K.  MOREOVER, THE INDEX IN JA AND A OF THE FIRST LOCATION
C    FOLLOWING THE LAST ELEMENT IN THE LAST ROW IS STORED IN IA(N+1). 
C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS GIVEN BY
C    IA(I+1) - IA(I),  THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED 
C    CONSECUTIVELY IN
C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),
C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN 
C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).
C    FOR EXAMPLE, THE 5 BY 5 MATRIX
C                ( 1. 0. 2. 0. 0.)
C                ( 0. 3. 0. 0. 0.)
C            M = ( 0. 4. 5. 6. 0.)
C                ( 0. 0. 0. 7. 0.)
C                ( 0. 0. 0. 8. 9.)
C    WOULD BE STORED AS
C               - 1  2  3  4  5  6  7  8  9
C            ---+--------------------------
C            IA - 1  3  4  7  8 10
C            JA - 1  3  2  2  3  4  4  4  5
C             A - 1. 2. 3. 4. 5. 6. 7. 8. 9.         .
C
C         THE STRICT UPPER (LOWER) TRIANGULAR PORTION OF THE MATRIX
C    U (L) IS STORED IN A SIMILAR FASHION USING THE ARRAYS  IU, JU, U 
C    (IL, JL, L)  EXCEPT THAT AN ADDITIONAL ARRAY IJU (IJL) IS USED TO
C    COMPRESS STORAGE OF JU (JL) BY ALLOWING SOME SEQUENCES OF COLUMN 
C    (ROW) INDICES TO USED FOR MORE THAN ONE ROW (COLUMN)  (N.B., L IS
C    STORED BY COLUMNS).  IJU(K) (IJL(K)) POINTS TO THE STARTING
C    LOCATION IN JU (JL) OF ENTRIES FOR THE KTH ROW (COLUMN).
C    COMPRESSION IN JU (JL) OCCURS IN TWO WAYS.  FIRST, IF A ROW
C    (COLUMN) I WAS MERGED INTO THE CURRENT ROW (COLUMN) K, AND THE
C    NUMBER OF ELEMENTS MERGED IN FROM (THE TAIL PORTION OF) ROW
C    (COLUMN) I IS THE SAME AS THE FINAL LENGTH OF ROW (COLUMN) K, THEN
C    THE KTH ROW (COLUMN) AND THE TAIL OF ROW (COLUMN) I ARE IDENTICAL
C    AND IJU(K) (IJL(K)) POINTS TO THE START OF THE TAIL.  SECOND, IF 
C    SOME TAIL PORTION OF THE (K-1)ST ROW (COLUMN) IS IDENTICAL TO THE
C    HEAD OF THE KTH ROW (COLUMN), THEN IJU(K) (IJL(K)) POINTS TO THE 
C    START OF THAT TAIL PORTION.  FOR EXAMPLE, THE NONZERO STRUCTURE OF
C    THE STRICT UPPER TRIANGULAR PART OF THE MATRIX
C            D 0 X X X
C            0 D 0 X X
C            0 0 D X 0
C            0 0 0 D X
C            0 0 0 0 D
C    WOULD BE REPRESENTED AS
C                - 1 2 3 4 5 6
C            ----+------------
C             IU - 1 4 6 7 8 8
C             JU - 3 4 5 4
C            IJU - 1 2 4 3           .
C    THE DIAGONAL ENTRIES OF L AND U ARE ASSUMED TO BE EQUAL TO ONE AND
C    ARE NOT STORED.  THE ARRAY D CONTAINS THE RECIPROCALS OF THE
C    DIAGONAL ENTRIES OF THE MATRIX D.
C
C    III. ADDITIONAL STORAGE SAVINGS
C         IN NSFC, R AND IC CAN BE THE SAME ARRAY IN THE CALLING
C    SEQUENCE IF NO REORDERING OF THE COEFFICIENT MATRIX HAS BEEN DONE.
C         IN NNFC, R, C, AND IC CAN ALL BE THE SAME ARRAY IF NO
C    REORDERING HAS BEEN DONE.  IF ONLY THE ROWS HAVE BEEN REORDERED, 
C    THEN C AND IC CAN BE THE SAME ARRAY.  IF THE ROW AND COLUMN
C    ORDERINGS ARE THE SAME, THEN R AND C CAN BE THE SAME ARRAY.  Z AND
C    ROW CAN BE THE SAME ARRAY.
C         IN NNSC OR NNTC, R AND C CAN BE THE SAME ARRAY IF NO
C    REORDERING HAS BEEN DONE OR IF THE ROW AND COLUMN ORDERINGS ARE THE
C    SAME.  Z AND B CAN BE THE SAME ARRAY.  HOWEVER, THEN B WILL BE
C    DESTROYED.
C
C    IV.  PARAMETERS
C         FOLLOWING IS A LIST OF PARAMETERS TO THE PROGRAMS.  NAMES ARE
C    UNIFORM AMONG THE VARIOUS SUBROUTINES.  CLASS ABBREVIATIONS ARE --
C       N - INTEGER VARIABLE
C       F - REAL VARIABLE
C       V - SUPPLIES A VALUE TO A SUBROUTINE
C       R - RETURNS A RESULT FROM A SUBROUTINE
C       I - USED INTERNALLY BY A SUBROUTINE
C       A - ARRAY
C
C CLASS - PARAMETER 
C ------+---------- 
C FVA   - A     - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED 
C       -           BY ROWS.
C       -           SIZE = NUMBER OF NONZERO ENTRIES IN M.
C FVA   - B     - RIGHT-HAND SIDE B.
C       -           SIZE = N. 
C NVA   - C     - ORDERING OF THE COLUMNS OF M.
C       -           SIZE = N. 
C FVRA  - D     - RECIPROCALS OF THE DIAGONAL ENTRIES OF THE MATRIX D.
C       -           SIZE = N. 
C NR    - FLAG  - ERROR FLAG.  VALUES AND THEIR MEANINGS ARE --
C       -            0     NO ERRORS DETECTED
C       -            N+K   NULL ROW IN A  --  ROW = K
C       -           2N+K   DUPLICATE ENTRY IN A  --  ROW = K
C       -           3N+K   INSUFFICIENT STORAGE FOR JL  --  ROW = K
C       -           4N+1   INSUFFICIENT STORAGE FOR L
C       -           5N+K   NULL PIVOT  --  ROW = K
C       -           6N+K   INSUFFICIENT STORAGE FOR JU  --  ROW = K
C       -           7N+1   INSUFFICIENT STORAGE FOR U
C       -           8N+K   ZERO PIVOT  --  ROW = K
C NVA   - IA    - POINTERS TO DELIMIT THE ROWS OF A.
C       -           SIZE = N+1.
C NVRA  - IJL   - POINTERS TO THE FIRST ELEMENT IN EACH COLUMN IN JL, 
C       -           USED TO COMPRESS STORAGE IN JL.
C       -           SIZE = N. 
C NVRA  - IJU   - POINTERS TO THE FIRST ELEMENT IN EACH ROW IN JU, USED
C       -           TO COMPRESS STORAGE IN JU.
C       -           SIZE = N. 
C NVRA  - IL    - POINTERS TO DELIMIT THE COLUMNS OF L.
C       -           SIZE = N+1.
C NVRA  - IU    - POINTERS TO DELIMIT THE ROWS OF U.
C       -           SIZE = N+1.
C NVA   - JA    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A.
C       -           SIZE = SIZE OF A.
C NVRA  - JL    - ROW NUMBERS CORRESPONDING TO THE ELEMENTS OF L.
C       -           SIZE = JLMAX.
C NV    - JLMAX - DECLARED DIMENSION OF JL.  JLMAX MUST BE LARGER THAN
C       -           THE NUMBER OF NONZEROS IN THE STRICT LOWER TRIANGLE
C       -           OF M PLUS FILLIN MINUS COMPRESSION.
C NVRA  - JU    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF U.
C       -           SIZE = JUMAX.
C NV    - JUMAX - DECLARED DIMENSION OF JU.  JUMAX MUST BE LARGER THAN
C       -           THE NUMBER OF NONZEROS IN THE STRICT UPPER TRIANGLE
C       -           OF M PLUS FILLIN MINUS COMPRESSION.
C FVRA  - L     - NONZERO ENTRIES IN THE STRICT LOWER TRIANGULAR PORTION
C       -           OF THE MATRIX L, STORED BY COLUMNS.
C       -           SIZE = LMAX.
C NV    - LMAX  - DECLARED DIMENSION OF L.  LMAX MUST BE LARGER THAN
C       -           THE NUMBER OF NONZEROS IN THE STRICT LOWER TRIANGLE
C       -           OF M PLUS FILLIN  (IL(N+1)-1 AFTER NSFC).
C NV    - N     - NUMBER OF VARIABLES/EQUATIONS.
C NVA   - R     - ORDERING OF THE ROWS OF M.
C       -           SIZE = N. 
C FVRA  - U     - NONZERO ENTRIES IN THE STRICT UPPER TRIANGULAR PORTION
C       -           OF THE MATRIX U, STORED BY ROWS.
C       -           SIZE = UMAX.
C NV    - UMAX  - DECLARED DIMENSION OF U.  UMAX MUST BE LARGER THAN
C       -           THE NUMBER OF NONZEROS IN THE STRICT UPPER TRIANGLE
C       -           OF M PLUS FILLIN  (IU(N+1)-1 AFTER NSFC).
C FRA   - Z     - SOLUTION X. 
C       -           SIZE = N. 
C
C       ----------------------------------------------------------------
C
C*** SUBROUTINE NROC
C*** REORDERS ROWS OF A, LEAVING ROW ORDER UNCHANGED
C
C
C       INPUT PARAMETERS.. N, IC, IA, JA, A
C       OUTPUT PARAMETERS.. JA, A, FLAG 
C
C       PARAMETERS USED INTERNALLY..
C NIA   - P     - AT THE KTH STEP, P IS A LINKED LIST OF THE REORDERED
C       -           COLUMN INDICES OF THE KTH ROW OF A.  P(N+1) POINTS
C       -           TO THE FIRST ENTRY IN THE LIST.
C       -           SIZE = N+1.
C NIA   - JAR   - AT THE KTH STEP,JAR CONTAINS THE ELEMENTS OF THE
C       -           REORDERED COLUMN INDICES OF A.
C       -           SIZE = N. 
C FIA   - AR    - AT THE KTH STEP, AR CONTAINS THE ELEMENTS OF THE
C       -           REORDERED ROW OF A. 
C       -           SIZE = N. 
C
      INTEGER  IC(1), IA(1), JA(1), JAR(1), P(1), FLAG
      DOUBLE PRECISION  A(1), AR(1)
C
C  ******  FOR EACH NONEMPTY ROW  *******************************
      DO 5 K=1,N
        JMIN = IA(K)
        JMAX = IA(K+1) - 1
        IF(JMIN .GT. JMAX) GO TO 5
        P(N+1) = N + 1
C  ******  INSERT EACH ELEMENT IN THE LIST  *********************
        DO 3 J=JMIN,JMAX
          NEWJ = IC(JA(J))
          I = N + 1 
   1      IF(P(I) .GE. NEWJ) GO TO 2
            I = P(I)
            GO TO 1 
   2      IF(P(I) .EQ. NEWJ) GO TO 102
          P(NEWJ) = P(I)
          P(I) = NEWJ
          JAR(NEWJ) = JA(J)
          AR(NEWJ) = A(J)
   3      CONTINUE
C  ******  REPLACE OLD ROW IN JA AND A  *************************
        I = N + 1
        DO 4 J=JMIN,JMAX
          I = P(I)
          JA(J) = JAR(I)
   4      A(J) = AR(I)
   5    CONTINUE
      FLAG = 0
      RETURN
C
C ** ERROR.. DUPLICATE ENTRY IN A
 102  FLAG = N + K
      RETURN
      END 
