        SUBROUTINE SDRV
     *     (N, P,IP, IA,JA,A, B, Z, NSP,ISP,RSP,ESP, PATH, FLAG)          SD   6
C                                                                         SD   7
C  DESCRIPTION                                                            SD   8
C                                                                         SD   9
C    SDRV SOLVES SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEMS OF LINEAR     SD  10
C    EQUATIONS.  THE SOLUTION PROCESS IS DIVIDED INTO THREE STAGES --     SD  11
C                                                                         SD  12
C      SSF - THE COEFFICIENT MATRIX M IS FACTORED SYMBOLICALLY TO         SD  13
C            DETERMINE WHERE FILLIN WILL OCCUR DURING THE NUMERIC         SD  14
C            FACTORIZATION.                                               SD  15
C                                                                         SD  16
C      SNF - M IS FACTORED NUMERICALLY INTO THE PRODUCT UT-D-U, WHERE     SD  17
C            D IS DIAGONAL AND U IS UNIT UPPER TRIANGULAR.                SD  18
C                                                                         SD  19
C      SNS - THE LINEAR SYSTEM  MX = B  IS SOLVED USING THE UT-D-U        SD  20
C            FACTORIZATION FROM SNF.                                      SD  21
C                                                                         SD  22
C    FOR SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX, SSF AND SNF    SD  23
C    NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM);  THEN SNS IS DONE     SD  24
C    ONCE FOR EACH ADDITIONAL RIGHT-HAND SIDE.  FOR SEVERAL SYSTEMS       SD  25
C    WHOSE COEFFICIENT MATRICES HAVE THE SAME NONZERO STRUCTURE, SSF      SD  26
C    NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM);  THEN SNF AND SNS     SD  27
C    ARE DONE ONCE FOR EACH ADDITIONAL SYSTEM.                            SD  28
C                                                                         SD  29
C                                                                         SD  30
C  STORAGE OF SPARSE MATRICES                                             SD  31
C                                                                         SD  32
C    THE NONZERO ENTRIES OF THE MATRIX M ARE STORED ROW-BY-ROW IN THE     SD  33
C    ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO ENTRIES IN EACH ROW,    SD  34
C    WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY LIES.  THESE COLUMN       SD  35
C    INDICES ARE STORED IN THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN  SD  36
C    JA(K) = J.  TO IDENTIFY THE INDIVIDUAL ROWS, WE NEED TO KNOW WHERE   SD  37
C    EACH ROW STARTS.  THESE ROW POINTERS ARE STORED IN THE ARRAY IA;     SD  38
C    I.E., IF M(I,J) IS THE FIRST NONZERO ENTRY (STORED) IN THE I-TH ROW  SD  39
C    AND  A(K) = M(I,J),  THEN  IA(I) = K.  MOREOVER, IA(N+1) POINTS TO   SD  40
C    THE FIRST LOCATION FOLLOWING THE LAST ELEMENT IN THE LAST ROW.       SD  41
C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS  IA(I+1) - IA(I),     SD  42
C    THE NONZERO ENTRIES IN THE I-TH ROW ARE STORED CONSECUTIVELY IN      SD  43
C                                                                         SD  44
C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),                   SD  45
C                                                                         SD  46
C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN     SD  47
C                                                                         SD  48
C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).                  SD  49
C                                                                         SD  50
C    SINCE THE COEFFICIENT MATRIX IS SYMMETRIC, ONLY THE NONZERO ENTRIES  SD  51
C    IN THE UPPER TRIANGLE NEED BE STORED, FOR EXAMPLE, THE MATRIX        SD  52
C                                                                         SD  53
C             ( 1  0  2  3  0 )                                           SD  54
C             ( 0  4  0  0  0 )                                           SD  55
C         M = ( 2  0  5  6  0 )                                           SD  56
C             ( 3  0  6  7  8 )                                           SD  57
C             ( 0  0  0  8  9 )                                           SD  58
C                                                                         SD  59
C    COULD BE STORED AS                                                   SD  60
C                                                                         SD  61
C            \ 1  2  3  4  5  6  7  8  9 10 11 12 13                      SD  62
C         ---+--------------------------------------                      SD  63
C         IA \ 1  4  5  8 12 14                                           SD  64
C         JA \ 1  3  4  2  1  3  4  1  3  4  5  4  5                      SD  65
C          A \ 1  2  3  4  2  5  6  3  6  7  8  8  9                      SD  66
C                                                                         SD  67
C    OR (SYMMETRICALLY) AS                                                SD  68
C                                                                         SD  69
C            \ 1  2  3  4  5  6  7  8  9                                  SD  70
C         ---+--------------------------                                  SD  71
C         IA \ 1  4  5  7  9 10                                           SD  72
C         JA \ 1  3  4  2  3  4  4  5  5                                  SD  73
C          A \ 1  2  3  4  5  6  7  8  9          .                       SD  74
C                                                                         SD  75
C                                                                         SD  76
C  REORDERING THE ROWS AND COLUMNS OF M                                   SD  77
C                                                                         SD  78
C    A SYMMETRIC PERMUTATION OF THE ROWS AND COLUMNS OF THE COEFFICIENT   SD  79
C    MATRIX M (E.G., WHICH REDUCES FILLIN OR ENHANCES NUMERICAL           SD  80
C    STABILITY) MUST BE SPECIFIED.  THE SOLUTION Z IS RETURNED IN THE     SD  81
C    ORIGINAL ORDER.                                                      SD  82
C                                                                         SD  83
C    TO SPECIFY THE TRIVIAL ORDERING (I.E., THE IDENTITY PERMUTATION),    SD  84
C    SET  P(I) = IP(I) = I,  I=1,...,N.  IN THIS CASE, P AND IP CAN BE    SD  85
C    THE SAME ARRAY.                                                      SD  86
C                                                                         SD  87
C    IF A NONTRIVIAL ORDERING (I.E., NOT THE IDENTITY PERMUTATION) IS     SD  88
C    SPECIFIED AND M IS STORED SYMMETRICALLY (I.E., NOT BOTH M(I,J) AND   SD  89
C    M(J,I) ARE STORED FOR I NE J), THEN ODRV SHOULD BE CALLED (WITH      SD  90
C    PATH = 3 OR 5) TO SYMMETRICALLY REORDER (IA,JA,A) BEFORE CALLING     SD  91
C    SDRV.  THIS IS TO ENSURE THAT IF M(I,J) WILL BE IN THE UPPER         SD  92
C    TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN M(I,J) IS       SD  93
C    STORED IN ROW I (AND THUS M(J,I) IS NOT STORED);  WHEREAS IF M(I,J)  SD  94
C    WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS STORED IN  SD  95
C    ROW J (AND THUS M(I,J) IS NOT STORED).                               SD  96
C                                                                         SD  97
C                                                                         SD  98
C  PARAMETERS                                                             SD  99
C                                                                         SD 100
C    N    - NUMBER OF VARIABLES/EQUATIONS                                 SD 101
C                                                                         SD 102
C    P    - INTEGER ONE-DIMENSIONAL ARRAY SPECIFYING A PERMUTATION OF     SD 103
C           THE ROWS AND COLUMNS OF M;  DIMENSION = N                     SD 104
C                                                                         SD 105
C    IP   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE INVERSE OF THE   SD 106
C           PERMUTATION SPECIFIED IN P;  I.E., IP(P(I)) = I, I=1,...,N;   SD 107
C           DIMENSION = N                                                 SD 108
C                                                                         SD 109
C    IA   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO DELIMIT  SD 110
C           ROWS IN JA AND A;  DIMENSION = N+1                            SD 111
C                                                                         SD 112
C    JA   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE COLUMN INDICES   SD 113
C           CORRESPONDING TO THE ELEMENTS OF A;  DIMENSION = NUMBER OF    SD 114
C           NONZERO ENTRIES IN M STORED                                   SD 115
C                                                                         SD 116
C    A    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE NONZERO ENTRIES IN  SD 117
C           THE COEFFICIENT MATRIX M, STORED BY ROWS;  DIMENSION =        SD 118
C           NUMBER OF NONZERO ENTRIES IN M STORED                         SD 119
C                                                                         SD 120
C    B    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE RIGHT-HAND SIDE B;  SD 121
C           B AND Z CAN BE THE SAME ARRAY;  DIMENSION = N                 SD 122
C                                                                         SD 123
C    Z    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE SOLUTION X;  Z AND  SD 124
C           B CAN BE THE SAME ARRAY;  DIMENSION = N                       SD 125
C                                                                         SD 126
C    NSP  - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAYS ISP AND      SD 127
C           RSP;  NSP MUST BE (SUBSTANTIALLY) LARGER THAN  3N+2K,  WHERE  SD 128
C           K = NUMBER OF NONZERO ENTRIES IN THE UPPER TRIANGLE OF M      SD 129
C                                                                         SD 130
C    ISP  - INTEGER ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE;  ISP  SD 131
C           AND RSP SHOULD BE EQUIVALENCED;  DIMENSION = NSP              SD 132
C                                                                         SD 133
C    RSP  - REAL ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE;  RSP     SD 134
C           AND ISP SHOULD BE EQUIVALENCED;  DIMENSION = NSP              SD 135
C                                                                         SD 136
C    ESP  - INTEGER VARIABLE;  IF SUFFICIENT STORAGE WAS AVAILABLE TO     SD 137
C           PERFORM THE SYMBOLIC FACTORIZATION (SSF), THEN ESP IS SET TO  SD 138
C           THE AMOUNT OF EXCESS STORAGE PROVIDED (NEGATIVE IF            SD 139
C           INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE NUMERIC     SD 140
C           FACTORIZATION (SNF))                                          SD 141
C                                                                         SD 142
C    PATH - INTEGER PATH SPECIFICATION;  VALUES AND THEIR MEANINGS ARE -  SD 143
C             1  PERFORM SSF, SNF, AND SNS                                SD 144
C             2  PERFORM SNF AND SNS (ISP/RSP IS ASSUMED TO HAVE BEEN     SD 145
C                  SET UP IN AN EARLIER CALL TO SDRV (FOR SSF))           SD 146
C             3  PERFORM SNS ONLY (ISP/RSP IS ASSUMED TO HAVE BEEN SET    SD 147
C                  UP IN AN EARLIER CALL TO SDRV (FOR SSF AND SNF))       SD 148
C             4  PERFORM SSF                                              SD 149
C             5  PERFORM SSF AND SNF                                      SD 150
C             6  PERFORM SNF ONLY (ISP/RSP IS ASSUMED TO HAVE BEEN SET    SD 151
C                  UP IN AN EARLIER CALL TO SDRV (FOR SSF))               SD 152
C                                                                         SD 153
C    FLAG - INTEGER ERROR FLAG;  VALUES AND THEIR MEANINGS ARE -          SD 154
C               0     NO ERRORS DETECTED                                  SD 155
C              2N+K   DUPLICATE ENTRY IN A  --  ROW = K                   SD 156
C              6N+K   INSUFFICIENT STORAGE IN SSF  --  ROW = K            SD 157
C              7N+1   INSUFFICIENT STORAGE IN SNF                         SD 158
C              8N+K   ZERO PIVOT  --  ROW = K                             SD 159
C             10N+1   INSUFFICIENT STORAGE IN SDRV                        SD 160
C             11N+1   ILLEGAL PATH SPECIFICATION                          SD 161
C                                                                         SD 162
C                                                                         SD 163
C  CONVERSION FROM REAL TO DOUBLE PRECISION                               SD 164
C                                                                         SD 165
C    CHANGE THE REAL DECLARATIONS IN SDRV, SNF, AND SNS TO DOUBLE         SD 166
C    PRECISION DECLARATIONS;  AND CHANGE THE VALUE IN THE DATA STATEMENT  SD 167
C    FOR THE INTEGER VARIABLE RATIO (IN SDRV) FROM 1 TO 2.                SD 168
C                                                                         SD 169
C-----------------------------------------------------------------------  SD 170
C                                                                         SD 171
        INTEGER  P(1), IP(1),  IA(1), JA(1),  ISP(1), ESP,  PATH,  FLAG,  SD 172
     *     RATIO,  Q, MARK, D, U, TMP,  UMAX                              SD 173
C       REAL  A(1),  B(1),  Z(1),  RSP(1)                                 SD 174
C       DATA  RATIO/1/                                                    SD 175
        DOUBLE PRECISION  A(1),  B(1),  Z(1),  RSP(1)                     SD 176
        DATA  RATIO/2/                                                    SD 177
C                                                                         SD 178
C----VALIDATE PATH SPECIFICATION                                          SD 179
        IF (PATH.LT.1 .OR. 6.LT.PATH)  GO TO 111                          SD 180
C                                                                         SD 181
C----ALLOCATE STORAGE AND FACTOR M SYMBOLICALLY TO DETERMINE FILL-IN      SD 182
        IJU   = 1                                                         SD 183
        IU    = IJU     +  N                                              SD 184
        JL    = IU      + N+1                                             SD 185
        JU    = JL      +  N                                              SD 186
        Q     = (NSP+1) -  N                                              SD 187
        MARK  = Q       -  N                                              SD 188
        JUMAX = MARK    - JU                                              SD 189
C                                                                         SD 190
        IF ((PATH-1) * (PATH-4) * (PATH-5) .NE. 0)  GO TO 1               SD 191
          IF (JUMAX.LE.0)  GO TO 110                                      SD 192
          CALL SSF                                                        SD 193
     *       (N,  P, IP,  IA, JA,  ISP(IJU), ISP(JU), ISP(IU), JUMAX,     SD 194
     *        ISP(Q), ISP(MARK), ISP(JL), FLAG)                           SD 195
          IF (FLAG.NE.0)  GO TO 100                                       SD 196
C                                                                         SD 197
C----ALLOCATE STORAGE AND FACTOR M NUMERICALLY                            SD 198
   1    IL   = JU      + ISP(IJU+(N-1))                                   SD 199
        TMP  = ((IL-1)+(RATIO-1)) / RATIO  +  1                           SD 200
        D    = TMP     + N                                                SD 201
        U    = D       + N                                                SD 202
        UMAX = (NSP+1) - U                                                SD 203
        ESP  = UMAX    - (ISP(IU+N)-1)                                    SD 204
C                                                                         SD 205
        IF ((PATH-1) * (PATH-2) * (PATH-5) * (PATH-6) .NE. 0)  GO TO 2    SD 206
          IF (UMAX.LE.0)  GO TO 110                                       SD 207
          CALL SNF                                                        SD 208
     *       (N,  P, IP,  IA, JA, A,                                      SD 209
     *        RSP(D),  ISP(IJU), ISP(JU), ISP(IU), RSP(U), UMAX,          SD 210
     *        ISP(IL),  ISP(JL),  FLAG)                                   SD 211
          IF (FLAG.NE.0)  GO TO 100                                       SD 212
C                                                                         SD 213
C----SOLVE SYSTEM OF LINEAR EQUATIONS  MX = B                             SD 214
   2    IF ((PATH-1) * (PATH-2) * (PATH-3) .NE. 0)  GO TO 3               SD 215
          IF (UMAX.LE.0)  GO TO 110                                       SD 216
          CALL SNS                                                        SD 217
     *       (N,  P,  RSP(D), ISP(IJU), ISP(JU), ISP(IU), RSP(U),  Z, B,  SD 218
     *        RSP(TMP))                                                   SD 219
C                                                                         SD 220
   3    RETURN                                                            SD 221
C                                                                         SD 222
C ** ERROR -- ERROR DETECTED IN SSF, SNF, OR SNS                          SD 223
 100    RETURN                                                            SD 224
C ** ERROR -- INSUFFICIENT STORAGE                                        SD 225
 110    FLAG = 10*N + 1                                                   SD 226
        RETURN                                                            SD 227
C ** ERROR -- ILLEGAL PATH SPECIFICATION                                  SD 228
 111    FLAG = 11*N + 1                                                   SD 229
        RETURN                                                            SD 230
        END                                                               SD 231
