|
5164
|
1 c======================================================================= |
|
|
2 c== umf4hb ============================================================= |
|
|
3 c======================================================================= |
|
|
4 |
|
|
5 c----------------------------------------------------------------------- |
|
|
6 c UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE |
|
|
7 c Dept, Univ. of Florida. All Rights Reserved. See ../Doc/License for |
|
|
8 c License. web: http://www.cise.ufl.edu/research/sparse/umfpack |
|
|
9 c----------------------------------------------------------------------- |
|
|
10 |
|
|
11 c umf4hb64: |
|
|
12 c read a sparse matrix in the Harwell/Boeing format, factorizes |
|
|
13 c it, and solves Ax=b. Also saves and loads the factors to/from a |
|
|
14 c file. Saving to a file is not required, it's just here to |
|
|
15 c demonstrate how to use this feature of UMFPACK. This program |
|
|
16 c only works on square RUA-type matrices. |
|
|
17 c |
|
|
18 c This is HIGHLY non-portable. It may not work with your C and |
|
|
19 c FORTRAN compilers. See umf4_f77wrapper.c for more details. |
|
|
20 c |
|
|
21 c usage (for example): |
|
|
22 c |
|
|
23 c in a Unix shell: |
|
|
24 c umf4hb64 < HB/arc130.rua |
|
|
25 |
|
|
26 integer*8 |
|
|
27 $ nzmax, nmax |
|
|
28 parameter (nzmax = 5000000, nmax = 160000) |
|
|
29 integer*8 |
|
|
30 $ Ap (nmax), Ai (nzmax), n, nz, totcrd, ptrcrd, i, j, p, |
|
|
31 $ indcrd, valcrd, rhscrd, ncol, nrow, nrhs, nzrhs, nel, |
|
|
32 $ numeric, symbolic, status, sys, filenum |
|
|
33 |
|
|
34 character title*72, key*30, type*3, ptrfmt*16, |
|
|
35 $ indfmt*16, valfmt*20, rhsfmt*20 |
|
|
36 double precision Ax (nzmax), x (nmax), b (nmax), aij, xj, |
|
|
37 $ r (nmax), control (20), info (90) |
|
|
38 character rhstyp*3 |
|
|
39 |
|
|
40 c ---------------------------------------------------------------- |
|
|
41 c read the Harwell/Boeing matrix |
|
|
42 c ---------------------------------------------------------------- |
|
|
43 |
|
|
44 read (5, 10, err = 998) |
|
|
45 $ title, key, |
|
|
46 $ totcrd, ptrcrd, indcrd, valcrd, rhscrd, |
|
|
47 $ type, nrow, ncol, nz, nel, |
|
|
48 $ ptrfmt, indfmt, valfmt, rhsfmt |
|
|
49 if (rhscrd .gt. 0) then |
|
|
50 c new Harwell/Boeing format: |
|
|
51 read (5, 20, err = 998) rhstyp, nrhs, nzrhs |
|
|
52 endif |
|
|
53 10 format (a72, a8 / 5i14 / a3, 11x, 4i14 / 2a16, 2a20) |
|
|
54 20 format (a3, 11x, 2i14) |
|
|
55 |
|
|
56 print *, 'Matrix key: ', key |
|
|
57 |
|
|
58 n = nrow |
|
|
59 if (type .ne. 'RUA' .or. nrow .ne. ncol) then |
|
|
60 print *, 'Error: can only handle square RUA matrices' |
|
|
61 stop |
|
|
62 endif |
|
|
63 if (n .ge. nmax .or. nz .gt. nzmax) then |
|
|
64 print *, ' Matrix too big!' |
|
|
65 stop |
|
|
66 endif |
|
|
67 |
|
|
68 c read the matrix (1-based) |
|
|
69 read (5, ptrfmt, err = 998) (Ap (p), p = 1, ncol+1) |
|
|
70 read (5, indfmt, err = 998) (Ai (p), p = 1, nz) |
|
|
71 read (5, valfmt, err = 998) (Ax (p), p = 1, nz) |
|
|
72 |
|
|
73 c ---------------------------------------------------------------- |
|
|
74 c create the right-hand-side, assume x (i) = 1 + i/n |
|
|
75 c ---------------------------------------------------------------- |
|
|
76 |
|
|
77 do 30 i = 1,n |
|
|
78 b (i) = 0 |
|
|
79 30 continue |
|
|
80 c b = A*x |
|
|
81 do 50 j = 1,n |
|
|
82 xj = j |
|
|
83 xj = 1 + xj / n |
|
|
84 do 40 p = Ap (j), Ap (j+1)-1 |
|
|
85 i = Ai (p) |
|
|
86 aij = Ax (p) |
|
|
87 b (i) = b (i) + aij * xj |
|
|
88 40 continue |
|
|
89 50 continue |
|
|
90 |
|
|
91 c ---------------------------------------------------------------- |
|
|
92 c convert from 1-based to 0-based |
|
|
93 c ---------------------------------------------------------------- |
|
|
94 |
|
|
95 do 60 j = 1, n+1 |
|
|
96 Ap (j) = Ap (j) - 1 |
|
|
97 60 continue |
|
|
98 do 70 p = 1, nz |
|
|
99 Ai (p) = Ai (p) - 1 |
|
|
100 70 continue |
|
|
101 |
|
|
102 c ---------------------------------------------------------------- |
|
|
103 c factor the matrix and save to a file |
|
|
104 c ---------------------------------------------------------------- |
|
|
105 |
|
|
106 c set default parameters |
|
|
107 call umf4def (control) |
|
|
108 |
|
|
109 c print control parameters. set control (1) to 1 to print |
|
|
110 c error messages only |
|
|
111 control (1) = 2 |
|
|
112 call umf4pcon (control) |
|
|
113 |
|
|
114 c pre-order and symbolic analysis |
|
|
115 call umf4sym (n, n, Ap, Ai, Ax, symbolic, control, info) |
|
|
116 |
|
|
117 c print statistics computed so far |
|
|
118 c call umf4pinf (control, info) could also be done. |
|
|
119 print 80, info (1), info (16), |
|
|
120 $ (info (21) * info (4)) / 2**20, |
|
|
121 $ (info (22) * info (4)) / 2**20, |
|
|
122 $ info (23), info (24), info (25) |
|
|
123 80 format ('symbolic analysis:',/, |
|
|
124 $ ' status: ', f5.0, /, |
|
|
125 $ ' time: ', e10.2, ' (sec)'/, |
|
|
126 $ ' estimates (upper bound) for numeric LU:', /, |
|
|
127 $ ' size of LU: ', f10.2, ' (MB)', /, |
|
|
128 $ ' memory needed: ', f10.2, ' (MB)', /, |
|
|
129 $ ' flop count: ', e10.2, / |
|
|
130 $ ' nnz (L): ', f10.0, / |
|
|
131 $ ' nnz (U): ', f10.0) |
|
|
132 |
|
|
133 c check umf4sym error condition |
|
|
134 if (info (1) .lt. 0) then |
|
|
135 print *, 'Error occurred in umf4sym: ', info (1) |
|
|
136 stop |
|
|
137 endif |
|
|
138 |
|
|
139 c numeric factorization |
|
|
140 call umf4num (Ap, Ai, Ax, symbolic, numeric, control, info) |
|
|
141 |
|
|
142 c print statistics for the numeric factorization |
|
|
143 c call umf4pinf (control, info) could also be done. |
|
|
144 print 90, info (1), info (66), |
|
|
145 $ (info (41) * info (4)) / 2**20, |
|
|
146 $ (info (42) * info (4)) / 2**20, |
|
|
147 $ info (43), info (44), info (45) |
|
|
148 90 format ('numeric factorization:',/, |
|
|
149 $ ' status: ', f5.0, /, |
|
|
150 $ ' time: ', e10.2, /, |
|
|
151 $ ' actual numeric LU statistics:', /, |
|
|
152 $ ' size of LU: ', f10.2, ' (MB)', /, |
|
|
153 $ ' memory needed: ', f10.2, ' (MB)', /, |
|
|
154 $ ' flop count: ', e10.2, / |
|
|
155 $ ' nnz (L): ', f10.0, / |
|
|
156 $ ' nnz (U): ', f10.0) |
|
|
157 |
|
|
158 c check umf4num error condition |
|
|
159 if (info (1) .lt. 0) then |
|
|
160 print *, 'Error occurred in umf4num: ', info (1) |
|
|
161 stop |
|
|
162 endif |
|
|
163 |
|
|
164 c save the symbolic analysis to the file s0.umf |
|
|
165 c note that this is not needed until another matrix is |
|
|
166 c factorized, below. |
|
|
167 filenum = 0 |
|
|
168 call umf4ssym (symbolic, filenum, status) |
|
|
169 if (status .lt. 0) then |
|
|
170 print *, 'Error occurred in umf4ssym: ', status |
|
|
171 stop |
|
|
172 endif |
|
|
173 |
|
|
174 c save the LU factors to the file n0.umf |
|
|
175 call umf4snum (numeric, filenum, status) |
|
|
176 if (status .lt. 0) then |
|
|
177 print *, 'Error occurred in umf4snum: ', status |
|
|
178 stop |
|
|
179 endif |
|
|
180 |
|
|
181 c free the symbolic analysis |
|
|
182 call umf4fsym (symbolic) |
|
|
183 |
|
|
184 c free the numeric factorization |
|
|
185 call umf4fnum (numeric) |
|
|
186 |
|
|
187 c No LU factors (symbolic or numeric) are in memory at this point. |
|
|
188 |
|
|
189 c ---------------------------------------------------------------- |
|
|
190 c load the LU factors back in, and solve the system |
|
|
191 c ---------------------------------------------------------------- |
|
|
192 |
|
|
193 c At this point the program could terminate and load the LU |
|
|
194 C factors (numeric) from the n0.umf file, and solve the |
|
|
195 c system (see below). Note that the symbolic object is not |
|
|
196 c required. |
|
|
197 |
|
|
198 c load the numeric factorization back in (filename: n0.umf) |
|
|
199 call umf4lnum (numeric, filenum, status) |
|
|
200 if (status .lt. 0) then |
|
|
201 print *, 'Error occurred in umf4lnum: ', status |
|
|
202 stop |
|
|
203 endif |
|
|
204 |
|
|
205 c solve Ax=b, without iterative refinement |
|
|
206 sys = 0 |
|
|
207 call umf4sol (sys, x, b, numeric, control, info) |
|
|
208 if (info (1) .lt. 0) then |
|
|
209 print *, 'Error occurred in umf4sol: ', info (1) |
|
|
210 stop |
|
|
211 endif |
|
|
212 |
|
|
213 c free the numeric factorization |
|
|
214 call umf4fnum (numeric) |
|
|
215 |
|
|
216 c No LU factors (symbolic or numeric) are in memory at this point. |
|
|
217 |
|
|
218 c print final statistics |
|
|
219 call umf4pinf (control, info) |
|
|
220 |
|
|
221 c print the residual. x (i) should be 1 + i/n |
|
|
222 call resid (n, nz, Ap, Ai, Ax, x, b, r) |
|
|
223 |
|
|
224 c ---------------------------------------------------------------- |
|
|
225 c load the symbolic analysis back in, and factorize a new matrix |
|
|
226 c ---------------------------------------------------------------- |
|
|
227 |
|
|
228 c Again, the program could terminate here, recreate the matrix, |
|
|
229 c and refactorize. Note that umf4sym is not called. |
|
|
230 |
|
|
231 c load the symbolic factorization back in (filename: s0.umf) |
|
|
232 call umf4lsym (symbolic, filenum, status) |
|
|
233 if (status .lt. 0) then |
|
|
234 print *, 'Error occurred in umf4lsym: ', status |
|
|
235 stop |
|
|
236 endif |
|
|
237 |
|
|
238 c arbitrarily change the values of the matrix but not the pattern |
|
|
239 do 100 p = 1, nz |
|
|
240 Ax (p) = Ax (p) + 3.14159 / 100.0 |
|
|
241 100 continue |
|
|
242 |
|
|
243 c numeric factorization of the modified matrix |
|
|
244 call umf4num (Ap, Ai, Ax, symbolic, numeric, control, info) |
|
|
245 if (info (1) .lt. 0) then |
|
|
246 print *, 'Error occurred in umf4num: ', info (1) |
|
|
247 stop |
|
|
248 endif |
|
|
249 |
|
|
250 c free the symbolic analysis |
|
|
251 call umf4fsym (symbolic) |
|
|
252 |
|
|
253 c create a new right-hand-side, assume x (i) = 7 - i/n |
|
|
254 do 110 i = 1,n |
|
|
255 b (i) = 0 |
|
|
256 110 continue |
|
|
257 c b = A*x, with the modified matrix A (note that A is now 0-based) |
|
|
258 do 130 j = 1,n |
|
|
259 xj = j |
|
|
260 xj = 7 - xj / n |
|
|
261 do 120 p = Ap (j) + 1, Ap (j+1) |
|
|
262 i = Ai (p) + 1 |
|
|
263 aij = Ax (p) |
|
|
264 b (i) = b (i) + aij * xj |
|
|
265 120 continue |
|
|
266 130 continue |
|
|
267 |
|
|
268 c ---------------------------------------------------------------- |
|
|
269 c solve Ax=b, with iterative refinement |
|
|
270 c ---------------------------------------------------------------- |
|
|
271 |
|
|
272 sys = 0 |
|
|
273 call umf4solr (sys, Ap, Ai, Ax, x, b, numeric, control, info) |
|
|
274 if (info (1) .lt. 0) then |
|
|
275 print *, 'Error occurred in umf4solr: ', info (1) |
|
|
276 stop |
|
|
277 endif |
|
|
278 |
|
|
279 c print the residual. x (i) should be 7 - i/n |
|
|
280 call resid (n, nz, Ap, Ai, Ax, x, b, r) |
|
|
281 |
|
|
282 c ---------------------------------------------------------------- |
|
|
283 c solve Ax=b, without iterative refinement, broken into steps |
|
|
284 c ---------------------------------------------------------------- |
|
|
285 |
|
|
286 c the factorization is PAQ=LU, PRAQ=LU, or P(R\A)Q=LU. |
|
|
287 |
|
|
288 c x = R*b (or x=R\b, or x=b, as appropriate) |
|
|
289 call umf4scal (x, b, numeric, status) |
|
|
290 if (status .lt. 0) then |
|
|
291 print *, 'Error occurred in umf4scal: ', status |
|
|
292 stop |
|
|
293 endif |
|
|
294 |
|
|
295 c solve P'Lr=x for r (using r as workspace) |
|
|
296 sys = 3 |
|
|
297 call umf4sol (sys, r, x, numeric, control, info) |
|
|
298 if (info (1) .lt. 0) then |
|
|
299 print *, 'Error occurred in umf4sol: ', info (1) |
|
|
300 stop |
|
|
301 endif |
|
|
302 |
|
|
303 c solve UQ'x=r for x |
|
|
304 sys = 9 |
|
|
305 call umf4sol (sys, x, r, numeric, control, info) |
|
|
306 if (info (1) .lt. 0) then |
|
|
307 print *, 'Error occurred in umf4sol: ', info (1) |
|
|
308 stop |
|
|
309 endif |
|
|
310 |
|
|
311 c free the numeric factorization |
|
|
312 call umf4fnum (numeric) |
|
|
313 |
|
|
314 c print the residual. x (i) should be 7 - i/n |
|
|
315 call resid (n, nz, Ap, Ai, Ax, x, b, r) |
|
|
316 |
|
|
317 stop |
|
|
318 998 print *, 'Read error: Harwell/Boeing matrix' |
|
|
319 stop |
|
|
320 end |
|
|
321 |
|
|
322 c======================================================================= |
|
|
323 c== resid ============================================================== |
|
|
324 c======================================================================= |
|
|
325 |
|
|
326 c Compute the residual, r = Ax-b, its max-norm, and print the max-norm |
|
|
327 C Note that A is zero-based. |
|
|
328 |
|
|
329 subroutine resid (n, nz, Ap, Ai, Ax, x, b, r) |
|
|
330 integer*8 |
|
|
331 $ n, nz, Ap (n+1), Ai (n), j, i, p |
|
|
332 double precision Ax (nz), x (n), b (n), r (n), rmax, aij |
|
|
333 |
|
|
334 do 10 i = 1, n |
|
|
335 r (i) = -b (i) |
|
|
336 10 continue |
|
|
337 |
|
|
338 do 30 j = 1,n |
|
|
339 do 20 p = Ap (j) + 1, Ap (j+1) |
|
|
340 i = Ai (p) + 1 |
|
|
341 aij = Ax (p) |
|
|
342 r (i) = r (i) + aij * x (j) |
|
|
343 20 continue |
|
|
344 30 continue |
|
|
345 |
|
|
346 rmax = 0 |
|
|
347 do 40 i = 1, n |
|
|
348 rmax = max (rmax, r (i)) |
|
|
349 40 continue |
|
|
350 |
|
|
351 print *, 'norm (A*x-b): ', rmax |
|
|
352 return |
|
|
353 end |
