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