        SUBROUTINE  CDRV
     *     (N, R,C,IC, IA,JA,A, B, Z, NSP,ISP,RSP,ESP, PATH, FLAG)        CD  12
C                                                                         CD  13
C    PARAMETERS                                                           CD  14
C    CLASS ABBREVIATIONS ARE--                                            CD  15
C       N - INTEGER VARIABLE                                              CD  16
C       F - REAL VARIABLE                                                 CD  17
C       V - SUPPLIES A VALUE TO THE DRIVER                                CD  18
C       R - RETURNS A RESULT FROM THE DRIVER                              CD  19
C       I - USED INTERNALLY BY THE DRIVER                                 CD  20
C       A - ARRAY                                                         CD  21
C                                                                         CD  22
C CLASS \ PARAMETER                                                       CD  23
C ------+----------                                                       CD  24
C       \                                                                 CD  25
C         THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED      CD  26
C    ROW-BY-ROW IN THE ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO       CD  27
C    ENTRIES IN EACH ROW, WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY      CD  28
C    LIES.  THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES    CD  29
C    OF M ARE STORED IN THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN     CD  30
C    JA(K) = J.  IN ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND   CD  31
C    HOW LONG IT IS.  THE INDEX POSITIONS IN JA AND A WHERE THE ROWS OF   CD  32
C    M BEGIN ARE STORED IN THE ARRAY IA;  I.E., IF M(I,J) IS THE FIRST    CD  33
C    NONZERO ENTRY (STORED) IN THE I-TH ROW AND A(K) = M(I,J),  THEN      CD  34
C    IA(I) = K.  MOREOVER, THE INDEX IN JA AND A OF THE FIRST LOCATION    CD  35
C    FOLLOWING THE LAST ELEMENT IN THE LAST ROW IS STORED IN IA(N+1).     CD  36
C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS GIVEN BY              CD  37
C    IA(I+1) - IA(I),  THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED     CD  38
C    CONSECUTIVELY IN                                                     CD  39
C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),                   CD  40
C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN     CD  41
C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).                  CD  42
C    FOR EXAMPLE, THE 5 BY 5 MATRIX                                       CD  43
C                ( 1. 0. 2. 0. 0.)                                        CD  44
C                ( 0. 3. 0. 0. 0.)                                        CD  45
C            M = ( 0. 4. 5. 6. 0.)                                        CD  46
C                ( 0. 0. 0. 7. 0.)                                        CD  47
C                ( 0. 0. 0. 8. 9.)                                        CD  48
C    WOULD BE STORED AS                                                   CD  49
C               \ 1  2  3  4  5  6  7  8  9                               CD  50
C            ---+--------------------------                               CD  51
C            IA \ 1  3  4  7  8 10                                        CD  52
C            JA \ 1  3  2  2  3  4  4  4  5                               CD  53
C             A \ 1. 2. 3. 4. 5. 6. 7. 8. 9.         .                    CD  54
C                                                                         CD  55
C NV    \ N     - NUMBER OF VARIABLES/EQUATIONS.                          CD  56
C FVA   \ A     - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED     CD  57
C       \           BY ROWS.                                              CD  58
C       \           SIZE = NUMBER OF NONZERO ENTRIES IN M.                CD  59
C NVA   \ IA    - POINTERS TO DELIMIT THE ROWS IN A.                      CD  60
C       \           SIZE = N+1.                                           CD  61
C NVA   \ JA    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A.      CD  62
C       \           SIZE = SIZE OF A.                                     CD  63
C FVA   \ B     - RIGHT-HAND SIDE B;  B AND Z CAN THE SAME ARRAY.         CD  64
C       \           SIZE = N.                                             CD  65
C FRA   \ Z     - SOLUTION X;  B AND Z CAN BE THE SAME ARRAY.             CD  66
C       \           SIZE = N.                                             CD  67
C                                                                         CD  68
C         THE ROWS AND COLUMNS OF THE ORIGINAL MATRIX M CAN BE            CD  69
C    REORDERED (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY)     CD  70
C    BEFORE CALLING THE DRIVER.  IF NO REORDERING IS DONE, THEN SET       CD  71
C    R(I) = C(I) = IC(I) = I  FOR I=1,...,N.  THE SOLUTION Z IS RETURNED  CD  72
C    IN THE ORIGINAL ORDER.                                               CD  73
C         IF THE COLUMNS HAVE BEEN REORDERED (I.E.,  C(I).NE.I  FOR SOME  CD  74
C    I), THEN THE DRIVER WILL CALL A SUBROUTINE (NROC) WHICH REARRANGES   CD  75
C    EACH ROW OF JA AND A, LEAVING THE ROWS IN THE ORIGINAL ORDER, BUT    CD  76
C    PLACING THE ELEMENTS OF EACH ROW IN INCREASING ORDER WITH RESPECT    CD  77
C    TO THE NEW ORDERING.  IF  PATH.NE.1,  THEN NROC IS ASSUMED TO HAVE   CD  78
C    BEEN CALLED ALREADY.                                                 CD  79
C                                                                         CD  80
C NVA   \ R     - ORDERING OF THE ROWS OF M.                              CD  81
C       \           SIZE = N.                                             CD  82
C NVA   \ C     - ORDERING OF THE COLUMNS OF M.                           CD  83
C       \           SIZE = N.                                             CD  84
C NVA   \ IC    - INVERSE OF THE ORDERING OF THE COLUMNS OF M;  I.E.,     CD  85
C       \           IC(C(I)) = I  FOR I=1,...,N.                          CD  86
C       \           SIZE = N.                                             CD  87
C                                                                         CD  88
C         THE SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS IS DIVIDED INTO  CD  89
C    THREE STAGES --                                                      CD  90
C      NSFC -- THE MATRIX M IS PROCESSED SYMBOLICALLY TO DETERMINE WHERE  CD  91
C               FILLIN WILL OCCUR DURING THE NUMERIC FACTORIZATION.       CD  92
C      NNFC -- THE MATRIX M IS FACTORED NUMERICALLY INTO THE PRODUCT LDU  CD  93
C               OF A UNIT LOWER TRIANGULAR MATRIX L, A DIAGONAL MATRIX    CD  94
C               D, AND A UNIT UPPER TRIANGULAR MATRIX U, AND THE SYSTEM   CD  95
C               MX = B  IS SOLVED.                                        CD  96
C      NNSC -- THE LINEAR SYSTEM  MX = B  IS SOLVED USING THE LDU         CD  97
C  OR           FACTORIZATION FROM NNFC.                                  CD  98
C      NNTC -- THE TRANSPOSED LINEAR SYSTEM  MT X = B  IS SOLVED USING    CD  99
C               THE LDU FACTORIZATION FROM NNF.                           CD 100
C    FOR SEVERAL SYSTEMS WHOSE COEFFICIENT MATRICES HAVE THE SAME         CD 101
C    NONZERO STRUCTURE, NSFC NEED BE DONE ONLY ONCE (FOR THE FIRST        CD 102
C    SYSTEM);  THEN NNFC IS DONE ONCE FOR EACH ADDITIONAL SYSTEM.  FOR    CD 103
C    SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX, NSFC AND NNFC      CD 104
C    NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM);  THEN NNSC OR NNTC    CD 105
C    IS DONE ONCE FOR EACH ADDITIONAL RIGHT-HAND SIDE.                    CD 106
C                                                                         CD 107
C NV    \ PATH  - PATH SPECIFICATION;  VALUES AND THEIR MEANINGS ARE --   CD 108
C       \           1  PERFORM NROC, NSFC, AND NNFC.                      CD 109
C       \           2  PERFORM NNFC ONLY  (NSFC IS ASSUMED TO HAVE BEEN   CD 110
C       \               DONE IN A MANNER COMPATIBLE WITH THE STORAGE      CD 111
C       \               ALLOCATION USED IN THE DRIVER).                   CD 112
C       \           3  PERFORM NNSC ONLY  (NSFC AND NNFC ARE ASSUMED TO   CD 113
C       \               HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE    CD 114
C       \               STORAGE ALLOCATION USED IN THE DRIVER).           CD 115
C       \           4  PERFORM NNTC ONLY  (NSFC AND NNFC ARE ASSUMED TO   CD 116
C       \               HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE    CD 117
C       \               STORAGE ALLOCATION USED IN THE DRIVER).           CD 118
C       \           5  PERFORM NROC AND NSFC.                             CD 119
C                                                                         CD 120
C         VARIOUS ERRORS ARE DETECTED BY THE DRIVER AND THE INDIVIDUAL    CD 121
C    SUBROUTINES.                                                         CD 122
C                                                                         CD 123
C NR    \ FLAG  - ERROR FLAG;  VALUES AND THEIR MEANINGS ARE --           CD 124
C       \             0     NO ERRORS DETECTED                            CD 125
C       \             N+K   NULL ROW IN A  --  ROW = K                    CD 126
C       \            2N+K   DUPLICATE ENTRY IN A  --  ROW = K             CD 127
C       \            3N+K   INSUFFICIENT STORAGE IN NSFC  --  ROW = K     CD 128
C       \            4N+1   INSUFFICIENT STORAGE IN NNFC                  CD 129
C       \            5N+K   NULL PIVOT  --  ROW = K                       CD 130
C       \            6N+K   INSUFFICIENT STORAGE IN NSFC  --  ROW = K     CD 131
C       \            7N+1   INSUFFICIENT STORAGE IN NNFC                  CD 132
C       \            8N+K   ZERO PIVOT  --  ROW = K                       CD 133
C       \           10N+1   INSUFFICIENT STORAGE IN CDRV                  CD 134
C       \           11N+1   ILLEGAL PATH SPECIFICATION                    CD 135
C                                                                         CD 136
C         WORKING STORAGE IS NEEDED FOR THE FACTORED FORM OF THE MATRIX   CD 137
C    M PLUS VARIOUS TEMPORARY VECTORS.  THE ARRAYS ISP AND RSP SHOULD BE  CD 138
C    EQUIVALENCED;  INTEGER STORAGE IS ALLOCATED FROM THE BEGINNING OF    CD 139
C    ISP AND REAL STORAGE FROM THE END OF RSP.                            CD 140
C                                                                         CD 141
C NV    \ NSP   - DECLARED DIMENSION OF RSP;  NSP GENERALLY MUST          CD 142
C       \           BE LARGER THAN  8N+2 + 2K  (WHERE  K = (NUMBER OF     CD 143
C       \           NONZERO ENTRIES IN M)).                               CD 144
C NVIRA \ ISP   - INTEGER WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS  CD 145
C       \           NEEDED BY THE SUBROUTINES;  ISP AND RSP SHOULD BE     CD 146
C       \           EQUIVALENCED.                                         CD 147
C       \           SIZE = LRATIO*NSP.                                    CD 148
C FVIRA \ RSP   - REAL WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS     CD 149
C       \           NEEDED BY THE SUBROUTINES;  ISP AND RSP SHOULD BE     CD 150
C       \           EQUIVALENCED.                                         CD 151
C       \           SIZE = NSP.                                           CD 152
C NR    \ ESP   - IF SUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE      CD 153
C       \           SYMBOLIC FACTORIZATION (NSFC), THEN ESP IS SET TO     CD 154
C       \           THE AMOUNT OF EXCESS STORAGE PROVIDED (NEGATIVE IF    CD 155
C       \           INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE     CD 156
C       \           NUMERIC FACTORIZATION (NNFC)).                        CD 157
C                                                                         CD 158
C                                                                         CD 159
C  CONVERSION TO DOUBLE PRECISION                                         CD 160
C                                                                         CD 161
C    TO CONVERT THESE ROUTINES FOR DOUBLE PRECISION ARRAYS, SIMPLY USE    CD 162
C    THE DOUBLE PRECISION DECLARATIONS IN PLACE OF THE REAL DECLARATIONS  CD 163
C    IN EACH SUBPROGRAM;  IN ADDITION, THE DATA VALUE OF THE INTEGER      CD 164
C    VARIABLE LRATIO MUST BE SET AS INDICATED IN SUBROUTINE CDRV          CD 165
C                                                                         CD 166
        INTEGER  R(1), C(1), IC(1),  IA(1), JA(1),  ISP(1), ESP,  PATH,   CD 167
     *     FLAG,  D, U, Q, ROW, TMP, AR,  UMAX                            CD 168
C       REAL  A(1), B(1), Z(1), RSP(1)                                    CD 169
        DOUBLE PRECISION  A(1), B(1), Z(1), RSP(1)                        CD 170
C                                                                         CD 171
C  SET LRATIO EQUAL TO THE RATIO BETWEEN THE LENGTH OF FLOATING POINT     CD 172
C  AND INTEGER ARRAY DATA;  E. G., LRATIO = 1 FOR (REAL, INTEGER),        CD 173
C  LRATIO = 2 FOR (DOUBLE PRECISION, INTEGER)                             CD 174
C                                                                         CD 175
        DATA LRATIO/1/                                                    CD 176
C                                                                         CD 177
        IF (PATH.LT.1 .OR. 5.LT.PATH)  GO TO 111                          CD 178
C******  INITIALIZE AND DIVIDE UP TEMPORARY STORAGE  *******************  CD 179
        IL   = 1                                                          CD 180
        IJL  = IL  + (N+1)                                                CD 181
        IU   = IJL +   N                                                  CD 182
        IJU  = IU  + (N+1)                                                CD 183
        IRL  = IJU +   N                                                  CD 184
        JRL  = IRL +   N                                                  CD 185
        JL   = JRL +   N                                                  CD 186
C                                                                         CD 187
C  ******  REORDER A IF NECESSARY, CALL NSFC IF FLAG IS SET  ***********  CD 188
        IF ((PATH-1) * (PATH-5) .NE. 0)  GO TO 5                          CD 189
          MAX = (LRATIO*NSP + 1 - JL) - (N+1) - 5*N                       CD 190
          JLMAX = MAX/2                                                   CD 191
          Q     = JL   + JLMAX                                            CD 192
          IRA   = Q    + (N+1)                                            CD 193
          JRA   = IRA  +   N                                              CD 194
          IRAC  = JRA  +   N                                              CD 195
          IRU   = IRAC +   N                                              CD 196
          JRU   = IRU  +   N                                              CD 197
          JUTMP = JRU  +   N                                              CD 198
          JUMAX = LRATIO*NSP  + 1 - JUTMP                                 CD 199
          ESP = MAX/LRATIO                                                CD 200
          IF (JLMAX.LE.0 .OR. JUMAX.LE.0)  GO TO 110                      CD 201
C                                                                         CD 202
          DO 1 I=1,N                                                      CD 203
            IF (C(I).NE.I)  GO TO 2                                       CD 204
   1        CONTINUE                                                      CD 205
          GO TO 3                                                         CD 206
   2      AR = NSP + 1 - N                                                CD 207
          CALL  NROC                                                      CD 208
     *       (N, IC, IA,JA,A, ISP(IL), RSP(AR), ISP(IU), FLAG)            CD 209
          IF (FLAG.NE.0)  GO TO 100                                       CD 210
C                                                                         CD 211
   3      CALL  NSFC                                                      CD 212
     *       (N, R, IC, IA,JA,                                            CD 213
     *        JLMAX, ISP(IL), ISP(JL), ISP(IJL),                          CD 214
     *        JUMAX, ISP(IU), ISP(JUTMP), ISP(IJU),                       CD 215
     *        ISP(Q), ISP(IRA), ISP(JRA), ISP(IRAC),                      CD 216
     *        ISP(IRL), ISP(JRL), ISP(IRU), ISP(JRU),  FLAG)              CD 217
          IF(FLAG .NE. 0)  GO TO 100                                      CD 218
C  ******  MOVE JU NEXT TO JL  *****************************************  CD 219
          JLMAX = ISP(IJL+N-1)                                            CD 220
          JU    = JL + JLMAX                                              CD 221
          JUMAX = ISP(IJU+N-1)                                            CD 222
          IF (JUMAX.LE.0)  GO TO 5                                        CD 223
          DO 4 J=1,JUMAX                                                  CD 224
   4        ISP(JU+J-1) = ISP(JUTMP+J-1)                                  CD 225
C                                                                         CD 226
C  ******  CALL REMAINING SUBROUTINES  *********************************  CD 227
   5    JLMAX = ISP(IJL+N-1)                                              CD 228
        JU    = JL  + JLMAX                                               CD 229
        JUMAX = ISP(IJU+N-1)                                              CD 230
        L     = (JU + JUMAX - 2 + LRATIO)  /  LRATIO    +    1            CD 231
        LMAX  = ISP(IL+N) - 1                                             CD 232
        D     = L   + LMAX                                                CD 233
        U     = D   + N                                                   CD 234
        ROW   = NSP + 1 - N                                               CD 235
        TMP   = ROW - N                                                   CD 236
        UMAX  = TMP - U                                                   CD 237
        ESP   = UMAX - (ISP(IU+N) - 1)                                    CD 238
C                                                                         CD 239
        IF ((PATH-1) * (PATH-2) .NE. 0)  GO TO 6                          CD 240
          IF (UMAX.LE.0)  GO TO 110                                       CD 241
          CALL NNFC                                                       CD 242
     *       (N,  R, C, IC,  IA, JA, A, Z, B,                             CD 243
     *        LMAX, ISP(IL), ISP(JL), ISP(IJL), RSP(L),  RSP(D),          CD 244
     *        UMAX, ISP(IU), ISP(JU), ISP(IJU), RSP(U),                   CD 245
     *        RSP(ROW), RSP(TMP),  ISP(IRL), ISP(JRL),  FLAG)             CD 246
          IF(FLAG .NE. 0)  GO TO 100                                      CD 247
C                                                                         CD 248
   6    IF ((PATH-3) .NE. 0)  GO TO 7                                     CD 249
          CALL NNSC                                                       CD 250
     *       (N,  R, C,  ISP(IL), ISP(JL), ISP(IJL), RSP(L),              CD 251
     *        RSP(D),    ISP(IU), ISP(JU), ISP(IJU), RSP(U),              CD 252
     *        Z, B,  RSP(TMP))                                            CD 253
C                                                                         CD 254
   7    IF ((PATH-4) .NE. 0)  GO TO 8                                     CD 255
          CALL NNTC                                                       CD 256
     *       (N,  R, C,  ISP(IL), ISP(JL), ISP(IJL), RSP(L),              CD 257
     *        RSP(D),    ISP(IU), ISP(JU), ISP(IJU), RSP(U),              CD 258
     *        Z, B,  RSP(TMP))                                            CD 259
   8    RETURN                                                            CD 260
C                                                                         CD 261
C ** ERROR:  ERROR DETECTED IN NROC, NSFC, NNFC, OR NNSC                  CD 262
 100    RETURN                                                            CD 263
C ** ERROR:  INSUFFICIENT STORAGE                                         CD 264
 110    FLAG = 10*N + 1                                                   CD 265
        RETURN                                                            CD 266
C ** ERROR:  ILLEGAL PATH SPECIFICATION                                   CD 267
 111    FLAG = 11*N + 1                                                   CD 268
        RETURN                                                            CD 269
        END                                                               CD 270
