| blob | 7b7a53e182c5577b6a6d701b2296fa673f228fb7 |
1 *DECK CDRV
2 subroutine cdrv
4 c*** subroutine cdrv
5 c*** driver for subroutines for solving sparse nonsymmetric systems of
6 c linear equations (compressed pointer storage)
7 c
8 c
9 c parameters
10 c class abbreviations are--
11 c n - integer variable
12 c f - real variable
13 c v - supplies a value to the driver
14 c r - returns a result from the driver
15 c i - used internally by the driver
16 c a - array
17 c
18 c class - parameter
19 c ------+----------
20 c -
21 c the nonzero entries of the coefficient matrix m are stored
22 c row-by-row in the array a. to identify the individual nonzero
23 c entries in each row, we need to know in which column each entry
24 c lies. the column indices which correspond to the nonzero entries
25 c of m are stored in the array ja. i.e., if a(k) = m(i,j), then
26 c ja(k) = j. in addition, we need to know where each row starts and
27 c how long it is. the index positions in ja and a where the rows of
28 c m begin are stored in the array ia. i.e., if m(i,j) is the first
29 c nonzero entry (stored) in the i-th row and a(k) = m(i,j), then
30 c ia(i) = k. moreover, the index in ja and a of the first location
31 c following the last element in the last row is stored in ia(n+1).
32 c thus, the number of entries in the i-th row is given by
33 c ia(i+1) - ia(i), the nonzero entries of the i-th row are stored
34 c consecutively in
35 c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
36 c and the corresponding column indices are stored consecutively in
37 c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
38 c for example, the 5 by 5 matrix
39 c ( 1. 0. 2. 0. 0.)
40 c ( 0. 3. 0. 0. 0.)
41 c m = ( 0. 4. 5. 6. 0.)
42 c ( 0. 0. 0. 7. 0.)
43 c ( 0. 0. 0. 8. 9.)
44 c would be stored as
45 c - 1 2 3 4 5 6 7 8 9
46 c ---+--------------------------
47 c ia - 1 3 4 7 8 10
48 c ja - 1 3 2 2 3 4 4 4 5
49 c a - 1. 2. 3. 4. 5. 6. 7. 8. 9. .
50 c
51 c nv - n - number of variables/equations.
52 c fva - a - nonzero entries of the coefficient matrix m, stored
53 c - by rows.
54 c - size = number of nonzero entries in m.
55 c nva - ia - pointers to delimit the rows in a.
56 c - size = n+1.
57 c nva - ja - column numbers corresponding to the elements of a.
58 c - size = size of a.
59 c fva - b - right-hand side b. b and z can the same array.
60 c - size = n.
61 c fra - z - solution x. b and z can be the same array.
62 c - size = n.
63 c
64 c the rows and columns of the original matrix m can be
65 c reordered (e.g., to reduce fillin or ensure numerical stability)
66 c before calling the driver. if no reordering is done, then set
67 c r(i) = c(i) = ic(i) = i for i=1,...,n. the solution z is returned
68 c in the original order.
69 c if the columns have been reordered (i.e., c(i).ne.i for some
70 c i), then the driver will call a subroutine (nroc) which rearranges
71 c each row of ja and a, leaving the rows in the original order, but
72 c placing the elements of each row in increasing order with respect
73 c to the new ordering. if path.ne.1, then nroc is assumed to have
74 c been called already.
75 c
76 c nva - r - ordering of the rows of m.
77 c - size = n.
78 c nva - c - ordering of the columns of m.
79 c - size = n.
80 c nva - ic - inverse of the ordering of the columns of m. i.e.,
81 c - ic(c(i)) = i for i=1,...,n.
82 c - size = n.
83 c
84 c the solution of the system of linear equations is divided into
85 c three stages --
86 c nsfc -- the matrix m is processed symbolically to determine where
87 c fillin will occur during the numeric factorization.
88 c nnfc -- the matrix m is factored numerically into the product ldu
89 c of a unit lower triangular matrix l, a diagonal matrix
90 c d, and a unit upper triangular matrix u, and the system
91 c mx = b is solved.
92 c nnsc -- the linear system mx = b is solved using the ldu
93 c or factorization from nnfc.
94 c nntc -- the transposed linear system mt x = b is solved using
95 c the ldu factorization from nnf.
96 c for several systems whose coefficient matrices have the same
97 c nonzero structure, nsfc need be done only once (for the first
98 c system). then nnfc is done once for each additional system. for
99 c several systems with the same coefficient matrix, nsfc and nnfc
100 c need be done only once (for the first system). then nnsc or nntc
101 c is done once for each additional right-hand side.
102 c
103 c nv - path - path specification. values and their meanings are --
104 c - 1 perform nroc, nsfc, and nnfc.
105 c - 2 perform nnfc only (nsfc is assumed to have been
106 c - done in a manner compatible with the storage
107 c - allocation used in the driver).
108 c - 3 perform nnsc only (nsfc and nnfc are assumed to
109 c - have been done in a manner compatible with the
110 c - storage allocation used in the driver).
111 c - 4 perform nntc only (nsfc and nnfc are assumed to
112 c - have been done in a manner compatible with the
113 c - storage allocation used in the driver).
114 c - 5 perform nroc and nsfc.
115 c
116 c various errors are detected by the driver and the individual
117 c subroutines.
118 c
119 c nr - flag - error flag. values and their meanings are --
120 c - 0 no errors detected
121 c - n+k null row in a -- row = k
122 c - 2n+k duplicate entry in a -- row = k
123 c - 3n+k insufficient storage in nsfc -- row = k
124 c - 4n+1 insufficient storage in nnfc
125 c - 5n+k null pivot -- row = k
126 c - 6n+k insufficient storage in nsfc -- row = k
127 c - 7n+1 insufficient storage in nnfc
128 c - 8n+k zero pivot -- row = k
129 c - 10n+1 insufficient storage in cdrv
130 c - 11n+1 illegal path specification
131 c
132 c working storage is needed for the factored form of the matrix
133 c m plus various temporary vectors. the arrays isp and rsp should be
134 c equivalenced. integer storage is allocated from the beginning of
135 c isp and real storage from the end of rsp.
136 c
137 c nv - nsp - declared dimension of rsp. nsp generally must
138 c - be larger than 8n+2 + 2k (where k = (number of
139 c - nonzero entries in m)).
140 c nvira - isp - integer working storage divided up into various arrays
141 c - needed by the subroutines. isp and rsp should be
142 c - equivalenced.
143 c - size = lratio*nsp.
144 c fvira - rsp - real working storage divided up into various arrays
145 c - needed by the subroutines. isp and rsp should be
146 c - equivalenced.
147 c - size = nsp.
148 c nr - esp - if sufficient storage was available to perform the
149 c - symbolic factorization (nsfc), then esp is set to
150 c - the amount of excess storage provided (negative if
151 c - insufficient storage was available to perform the
152 c - numeric factorization (nnfc)).
153 c
154 c
155 c conversion to double precision
156 c
157 c to convert these routines for double precision arrays..
158 c (1) use the double precision declarations in place of the real
159 c declarations in each subprogram, as given in comment cards.
160 c (2) change the data-loaded value of the integer lratio
161 c in subroutine cdrv, as indicated below.
162 c (3) change e0 to d0 in the constants in statement number 10
163 c in subroutine nnfc and the line following that.
164 c
167 c real a(*), b(*), z(*), rsp(*)
169 c
170 c set lratio equal to the ratio between the length of floating point
171 c and integer array data. e. g., lratio = 1 for (real, integer),
172 c lratio = 2 for (double precision, integer)
173 c
175 c
177 c******initialize and divide up temporary storage *******************
185 c
186 c ****** reorder a if necessary, call nsfc if flag is set ***********
200 c
206 call nroc
209 c
217 c ****** move ju next to jl *****************************************
224 c
225 c ****** call remaining subroutines *********************************
237 c
240 call nnfc
246 c
248 call nnsc
252 c
254 call nntc
259 c
260 c ** error.. error detected in nroc, nsfc, nnfc, or nnsc
262 c ** error.. insufficient storage
264 return
265 c ** error.. illegal path specification
267 return
268 end