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
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
|
subroutine nsfc
* (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju,
* q, ira,jra, irac, irl,jrl, iru,jru, flag)
clll. optimize
c*** subroutine nsfc
c*** symbolic ldu-factorization of nonsymmetric sparse matrix
c (compressed pointer storage)
c
c
c input variables.. n, r, ic, ia, ja, jlmax, jumax.
c output variables.. il, jl, ijl, iu, ju, iju, flag.
c
c parameters used internally..
c nia - q - suppose m* is the result of reordering m. if
c - processing of the ith row of m* (hence the ith
c - row of u) is being done, q(j) is initially
c - nonzero if m*(i,j) is nonzero (j.ge.i). since
c - values need not be stored, each entry points to the
c - next nonzero and q(n+1) points to the first. n+1
c - indicates the end of the list. for example, if n=9
c - and the 5th row of m* is
c - 0 x x 0 x 0 0 x 0
c - then q will initially be
c - a a a a 8 a a 10 5 (a - arbitrary).
c - as the algorithm proceeds, other elements of q
c - are inserted in the list because of fillin.
c - q is used in an analogous manner to compute the
c - ith column of l.
c - size = n+1.
c nia - ira, - vectors used to find the columns of m. at the kth
c nia - jra, step of the factorization, irac(k) points to the
c nia - irac head of a linked list in jra of row indices i
c - such that i .ge. k and m(i,k) is nonzero. zero
c - indicates the end of the list. ira(i) (i.ge.k)
c - points to the smallest j such that j .ge. k and
c - m(i,j) is nonzero.
c - size of each = n.
c nia - irl, - vectors used to find the rows of l. at the kth step
c nia - jrl of the factorization, jrl(k) points to the head
c - of a linked list in jrl of column indices j
c - such j .lt. k and l(k,j) is nonzero. zero
c - indicates the end of the list. irl(j) (j.lt.k)
c - points to the smallest i such that i .ge. k and
c - l(i,j) is nonzero.
c - size of each = n.
c nia - iru, - vectors used in a manner analogous to irl and jrl
c nia - jru to find the columns of u.
c - size of each = n.
c
c internal variables..
c jlptr - points to the last position used in jl.
c juptr - points to the last position used in ju.
c jmin,jmax - are the indices in a or u of the first and last
c elements to be examined in a given row.
c for example, jmin=ia(k), jmax=ia(k+1)-1.
c
integer cend, qm, rend, rk, vj
integer ia(1), ja(1), ira(1), jra(1), il(1), jl(1), ijl(1)
integer iu(1), ju(1), iju(1), irl(1), jrl(1), iru(1), jru(1)
integer r(1), ic(1), q(1), irac(1), flag
c
c ****** initialize pointers ****************************************
np1 = n + 1
jlmin = 1
jlptr = 0
il(1) = 1
jumin = 1
juptr = 0
iu(1) = 1
do 1 k=1,n
irac(k) = 0
jra(k) = 0
jrl(k) = 0
1 jru(k) = 0
c ****** initialize column pointers for a ***************************
do 2 k=1,n
rk = r(k)
iak = ia(rk)
if (iak .ge. ia(rk+1)) go to 101
jaiak = ic(ja(iak))
if (jaiak .gt. k) go to 105
jra(k) = irac(jaiak)
irac(jaiak) = k
2 ira(k) = iak
c
c ****** for each column of l and row of u **************************
do 41 k=1,n
c
c ****** initialize q for computing kth column of l *****************
q(np1) = np1
luk = -1
c ****** by filling in kth column of a ******************************
vj = irac(k)
if (vj .eq. 0) go to 5
3 qm = np1
4 m = qm
qm = q(m)
if (qm .lt. vj) go to 4
if (qm .eq. vj) go to 102
luk = luk + 1
q(m) = vj
q(vj) = qm
vj = jra(vj)
if (vj .ne. 0) go to 3
c ****** link through jru *******************************************
5 lastid = 0
lasti = 0
ijl(k) = jlptr
i = k
6 i = jru(i)
if (i .eq. 0) go to 10
qm = np1
jmin = irl(i)
jmax = ijl(i) + il(i+1) - il(i) - 1
long = jmax - jmin
if (long .lt. 0) go to 6
jtmp = jl(jmin)
if (jtmp .ne. k) long = long + 1
if (jtmp .eq. k) r(i) = -r(i)
if (lastid .ge. long) go to 7
lasti = i
lastid = long
c ****** and merge the corresponding columns into the kth column ****
7 do 9 j=jmin,jmax
vj = jl(j)
8 m = qm
qm = q(m)
if (qm .lt. vj) go to 8
if (qm .eq. vj) go to 9
luk = luk + 1
q(m) = vj
q(vj) = qm
qm = vj
9 continue
go to 6
c ****** lasti is the longest column merged into the kth ************
c ****** see if it equals the entire kth column *********************
10 qm = q(np1)
if (qm .ne. k) go to 105
if (luk .eq. 0) go to 17
if (lastid .ne. luk) go to 11
c ****** if so, jl can be compressed ********************************
irll = irl(lasti)
ijl(k) = irll + 1
if (jl(irll) .ne. k) ijl(k) = ijl(k) - 1
go to 17
c ****** if not, see if kth column can overlap the previous one *****
11 if (jlmin .gt. jlptr) go to 15
qm = q(qm)
do 12 j=jlmin,jlptr
if (jl(j).lt.qm) go to 12
if (jl(j).eq.qm) go to 13
go to 15
12 continue
go to 15
13 ijl(k) = j
do 14 i=j,jlptr
if (jl(i) .ne. qm) go to 15
qm = q(qm)
if (qm .gt. n) go to 17
14 continue
jlptr = j - 1
c ****** move column indices from q to jl, update vectors ***********
15 jlmin = jlptr + 1
ijl(k) = jlmin
if (luk .eq. 0) go to 17
jlptr = jlptr + luk
if (jlptr .gt. jlmax) go to 103
qm = q(np1)
do 16 j=jlmin,jlptr
qm = q(qm)
16 jl(j) = qm
17 irl(k) = ijl(k)
il(k+1) = il(k) + luk
c
c ****** initialize q for computing kth row of u ********************
q(np1) = np1
luk = -1
c ****** by filling in kth row of reordered a ***********************
rk = r(k)
jmin = ira(k)
jmax = ia(rk+1) - 1
if (jmin .gt. jmax) go to 20
do 19 j=jmin,jmax
vj = ic(ja(j))
qm = np1
18 m = qm
qm = q(m)
if (qm .lt. vj) go to 18
if (qm .eq. vj) go to 102
luk = luk + 1
q(m) = vj
q(vj) = qm
19 continue
c ****** link through jrl, ******************************************
20 lastid = 0
lasti = 0
iju(k) = juptr
i = k
i1 = jrl(k)
21 i = i1
if (i .eq. 0) go to 26
i1 = jrl(i)
qm = np1
jmin = iru(i)
jmax = iju(i) + iu(i+1) - iu(i) - 1
long = jmax - jmin
if (long .lt. 0) go to 21
jtmp = ju(jmin)
if (jtmp .eq. k) go to 22
c ****** update irl and jrl, *****************************************
long = long + 1
cend = ijl(i) + il(i+1) - il(i)
irl(i) = irl(i) + 1
if (irl(i) .ge. cend) go to 22
j = jl(irl(i))
jrl(i) = jrl(j)
jrl(j) = i
22 if (lastid .ge. long) go to 23
lasti = i
lastid = long
c ****** and merge the corresponding rows into the kth row **********
23 do 25 j=jmin,jmax
vj = ju(j)
24 m = qm
qm = q(m)
if (qm .lt. vj) go to 24
if (qm .eq. vj) go to 25
luk = luk + 1
q(m) = vj
q(vj) = qm
qm = vj
25 continue
go to 21
c ****** update jrl(k) and irl(k) ***********************************
26 if (il(k+1) .le. il(k)) go to 27
j = jl(irl(k))
jrl(k) = jrl(j)
jrl(j) = k
c ****** lasti is the longest row merged into the kth ***************
c ****** see if it equals the entire kth row ************************
27 qm = q(np1)
if (qm .ne. k) go to 105
if (luk .eq. 0) go to 34
if (lastid .ne. luk) go to 28
c ****** if so, ju can be compressed ********************************
irul = iru(lasti)
iju(k) = irul + 1
if (ju(irul) .ne. k) iju(k) = iju(k) - 1
go to 34
c ****** if not, see if kth row can overlap the previous one ********
28 if (jumin .gt. juptr) go to 32
qm = q(qm)
do 29 j=jumin,juptr
if (ju(j).lt.qm) go to 29
if (ju(j).eq.qm) go to 30
go to 32
29 continue
go to 32
30 iju(k) = j
do 31 i=j,juptr
if (ju(i) .ne. qm) go to 32
qm = q(qm)
if (qm .gt. n) go to 34
31 continue
juptr = j - 1
c ****** move row indices from q to ju, update vectors **************
32 jumin = juptr + 1
iju(k) = jumin
if (luk .eq. 0) go to 34
juptr = juptr + luk
if (juptr .gt. jumax) go to 106
qm = q(np1)
do 33 j=jumin,juptr
qm = q(qm)
33 ju(j) = qm
34 iru(k) = iju(k)
iu(k+1) = iu(k) + luk
c
c ****** update iru, jru ********************************************
i = k
35 i1 = jru(i)
if (r(i) .lt. 0) go to 36
rend = iju(i) + iu(i+1) - iu(i)
if (iru(i) .ge. rend) go to 37
j = ju(iru(i))
jru(i) = jru(j)
jru(j) = i
go to 37
36 r(i) = -r(i)
37 i = i1
if (i .eq. 0) go to 38
iru(i) = iru(i) + 1
go to 35
c
c ****** update ira, jra, irac **************************************
38 i = irac(k)
if (i .eq. 0) go to 41
39 i1 = jra(i)
ira(i) = ira(i) + 1
if (ira(i) .ge. ia(r(i)+1)) go to 40
irai = ira(i)
jairai = ic(ja(irai))
if (jairai .gt. i) go to 40
jra(i) = irac(jairai)
irac(jairai) = i
40 i = i1
if (i .ne. 0) go to 39
41 continue
c
ijl(n) = jlptr
iju(n) = juptr
flag = 0
return
c
c ** error.. null row in a
101 flag = n + rk
return
c ** error.. duplicate entry in a
102 flag = 2*n + rk
return
c ** error.. insufficient storage for jl
103 flag = 3*n + k
return
c ** error.. null pivot
105 flag = 5*n + k
return
c ** error.. insufficient storage for ju
106 flag = 6*n + k
return
end
|
