Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/Demo/umf4hb.f @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
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 umf4hb:
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 umf4hb < HB/arc130.rua
25
26 integer
27 $ nzmax, nmax
28 parameter (nzmax = 5000000, nmax = 160000)
29 integer
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
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