1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
|
subroutine cdrv
* (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)
clll. optimize
c*** subroutine cdrv
c*** driver for subroutines for solving sparse nonsymmetric systems of
c linear equations (compressed pointer storage)
c
c
c parameters
c class abbreviations are--
c n - integer variable
c f - real variable
c v - supplies a value to the driver
c r - returns a result from the driver
c i - used internally by the driver
c a - array
c
c class - parameter
c ------+----------
c -
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 nonzero entry (stored) 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 nv - n - number of variables/equations.
c fva - a - nonzero entries of the coefficient matrix m, stored
c - by rows.
c - size = number of nonzero entries in m.
c nva - ia - pointers to delimit the rows in a.
c - size = n+1.
c nva - ja - column numbers corresponding to the elements of a.
c - size = size of a.
c fva - b - right-hand side b. b and z can the same array.
c - size = n.
c fra - z - solution x. b and z can be the same array.
c - size = n.
c
c the rows and columns of the original matrix m can be
c reordered (e.g., to reduce fillin or ensure numerical stability)
c before calling the driver. if no reordering is done, then set
c r(i) = c(i) = ic(i) = i for i=1,...,n. the solution z is returned
c in the original order.
c if the columns have been reordered (i.e., c(i).ne.i for some
c i), then the driver will call a subroutine (nroc) which rearranges
c each row of ja and a, leaving the rows in the original order, but
c placing the elements of each row in increasing order with respect
c to the new ordering. if path.ne.1, then nroc is assumed to have
c been called already.
c
c nva - r - ordering of the rows of m.
c - size = n.
c nva - c - ordering of the columns of m.
c - size = n.
c nva - ic - inverse of the ordering of the columns of m. i.e.,
c - ic(c(i)) = i for i=1,...,n.
c - size = n.
c
c the solution of the system of linear equations is divided into
c three stages --
c nsfc -- the matrix m is processed symbolically to determine where
c fillin will occur during the numeric factorization.
c nnfc -- the matrix m is factored numerically into the product ldu
c of a unit lower triangular matrix l, a diagonal matrix
c d, and a unit upper triangular matrix u, and the system
c mx = b is solved.
c nnsc -- the linear system mx = b is solved using the ldu
c or factorization from nnfc.
c nntc -- the transposed linear system mt x = b is solved using
c the ldu factorization from nnf.
c for several systems whose coefficient matrices have the same
c nonzero structure, nsfc need be done only once (for the first
c system). then nnfc is done once for each additional system. for
c several systems with the same coefficient matrix, nsfc and nnfc
c need be done only once (for the first system). then nnsc or nntc
c is done once for each additional right-hand side.
c
c nv - path - path specification. values and their meanings are --
c - 1 perform nroc, nsfc, and nnfc.
c - 2 perform nnfc only (nsfc is assumed to have been
c - done in a manner compatible with the storage
c - allocation used in the driver).
c - 3 perform nnsc only (nsfc and nnfc are assumed to
c - have been done in a manner compatible with the
c - storage allocation used in the driver).
c - 4 perform nntc only (nsfc and nnfc are assumed to
c - have been done in a manner compatible with the
c - storage allocation used in the driver).
c - 5 perform nroc and nsfc.
c
c various errors are detected by the driver and the individual
c subroutines.
c
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 in nsfc -- row = k
c - 4n+1 insufficient storage in nnfc
c - 5n+k null pivot -- row = k
c - 6n+k insufficient storage in nsfc -- row = k
c - 7n+1 insufficient storage in nnfc
c - 8n+k zero pivot -- row = k
c - 10n+1 insufficient storage in cdrv
c - 11n+1 illegal path specification
c
c working storage is needed for the factored form of the matrix
c m plus various temporary vectors. the arrays isp and rsp should be
c equivalenced. integer storage is allocated from the beginning of
c isp and real storage from the end of rsp.
c
c nv - nsp - declared dimension of rsp. nsp generally must
c - be larger than 8n+2 + 2k (where k = (number of
c - nonzero entries in m)).
c nvira - isp - integer working storage divided up into various arrays
c - needed by the subroutines. isp and rsp should be
c - equivalenced.
c - size = lratio*nsp.
c fvira - rsp - real working storage divided up into various arrays
c - needed by the subroutines. isp and rsp should be
c - equivalenced.
c - size = nsp.
c nr - esp - if sufficient storage was available to perform the
c - symbolic factorization (nsfc), then esp is set to
c - the amount of excess storage provided (negative if
c - insufficient storage was available to perform the
c - numeric factorization (nnfc)).
c
c
c conversion to double precision
c
c to convert these routines for double precision arrays..
c (1) use the double precision declarations in place of the real
c declarations in each subprogram, as given in comment cards.
c (2) change the data-loaded value of the integer lratio
c in subroutine cdrv, as indicated below.
c (3) change e0 to d0 in the constants in statement number 10
c in subroutine nnfc and the line following that.
c
integer r(1), c(1), ic(1), ia(1), ja(1), isp(1), esp, path,
* flag, d, u, q, row, tmp, ar, umax
double precision a(1), b(1), z(1), rsp(1)
c
c set lratio equal to the ratio between the length of floating point
c and integer array data. e. g., lratio = 1 for (real, integer),
c lratio = 2 for (double precision, integer)
c
data lratio/2/
c
if (path.lt.1 .or. 5.lt.path) go to 111
c******initialize and divide up temporary storage *******************
il = 1
ijl = il + (n+1)
iu = ijl + n
iju = iu + (n+1)
irl = iju + n
jrl = irl + n
jl = jrl + n
c
c ****** reorder a if necessary, call nsfc if flag is set ***********
if ((path-1) * (path-5) .ne. 0) go to 5
max = (lratio*nsp + 1 - jl) - (n+1) - 5*n
jlmax = max/2
q = jl + jlmax
ira = q + (n+1)
jra = ira + n
irac = jra + n
iru = irac + n
jru = iru + n
jutmp = jru + n
jumax = lratio*nsp + 1 - jutmp
esp = max/lratio
if (jlmax.le.0 .or. jumax.le.0) go to 110
c
do 1 i=1,n
if (c(i).ne.i) go to 2
1 continue
go to 3
2 ar = nsp + 1 - n
call nroc
* (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)
if (flag.ne.0) go to 100
c
3 call nsfc
* (n, r, ic, ia,ja,
* jlmax, isp(il), isp(jl), isp(ijl),
* jumax, isp(iu), isp(jutmp), isp(iju),
* isp(q), isp(ira), isp(jra), isp(irac),
* isp(irl), isp(jrl), isp(iru), isp(jru), flag)
if(flag .ne. 0) go to 100
c ****** move ju next to jl *****************************************
jlmax = isp(ijl+n-1)
ju = jl + jlmax
jumax = isp(iju+n-1)
if (jumax.le.0) go to 5
do 4 j=1,jumax
4 isp(ju+j-1) = isp(jutmp+j-1)
c
c ****** call remaining subroutines *********************************
5 jlmax = isp(ijl+n-1)
ju = jl + jlmax
jumax = isp(iju+n-1)
l = (ju + jumax - 2 + lratio) / lratio + 1
lmax = isp(il+n) - 1
d = l + lmax
u = d + n
row = nsp + 1 - n
tmp = row - n
umax = tmp - u
esp = umax - (isp(iu+n) - 1)
c
if ((path-1) * (path-2) .ne. 0) go to 6
if (umax.lt.0) go to 110
call nnfc
* (n, r, c, ic, ia, ja, a, z, b,
* lmax, isp(il), isp(jl), isp(ijl), rsp(l), rsp(d),
* umax, isp(iu), isp(ju), isp(iju), rsp(u),
* rsp(row), rsp(tmp), isp(irl), isp(jrl), flag)
if(flag .ne. 0) go to 100
c
6 if ((path-3) .ne. 0) go to 7
call nnsc
* (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
* rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
* z, b, rsp(tmp))
c
7 if ((path-4) .ne. 0) go to 8
call nntc
* (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
* rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
* z, b, rsp(tmp))
8 return
c
c ** error.. error detected in nroc, nsfc, nnfc, or nnsc
100 return
c ** error.. insufficient storage
110 flag = 10*n + 1
return
c ** error.. illegal path specification
111 flag = 11*n + 1
return
end
|
