Mercurial > octave-nkf
comparison liboctave/UMFPACK/UMFPACK/Demo/umfpack_di_demo.c @ 5164:57077d0ddc8e
[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 /* ========================================================================== */ | |
2 /* === umfpack_di_demo ====================================================== */ | |
3 /* ========================================================================== */ | |
4 | |
5 | |
6 /* -------------------------------------------------------------------------- */ | |
7 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */ | |
8 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */ | |
9 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */ | |
10 /* -------------------------------------------------------------------------- */ | |
11 | |
12 /* | |
13 A demo of UMFPACK: umfpack_di_* version. | |
14 | |
15 First, factor and solve a 5-by-5 system, Ax=b, using default parameters. | |
16 Then solve A'x=b using the factors of A. Modify one entry (A (1,4) = 0, | |
17 where the row and column indices range from 0 to 4. The pattern of A | |
18 has not changed (it has explicitly zero entry), so a reanalysis with | |
19 umfpack_di_symbolic does not need to be done. Refactorize (with | |
20 umfpack_di_numeric), and solve Ax=b. Note that the pivot ordering has | |
21 changed. Next, change all of the entries in A, but not the pattern. | |
22 | |
23 Finally, compute C = A', and do the symbolic and numeric factorization of C. | |
24 Factorizing A' can sometimes be better than factorizing A itself (less work | |
25 and memory usage). Solve C'x=b twice; the solution is the same as the | |
26 solution to Ax=b. | |
27 | |
28 A note about zero-sized arrays: UMFPACK uses many user-provided arrays of | |
29 size n (order of the matrix), and of size nz (the number of nonzeros in a | |
30 matrix). n cannot be zero; UMFPACK does not handle zero-dimensioned arrays. | |
31 However, nz can be zero. If you attempt to malloc an array of size nz = 0, | |
32 however, malloc will return a null pointer which UMFPACK will report as a | |
33 "missing argument." Thus, nz1 in this code is set to MAX (nz,1), and | |
34 similarly for lnz and unz. Lnz can never be zero, however, since L is always | |
35 unit diagonal. | |
36 */ | |
37 | |
38 /* -------------------------------------------------------------------------- */ | |
39 /* definitions */ | |
40 /* -------------------------------------------------------------------------- */ | |
41 | |
42 #include <stdio.h> | |
43 #include <stdlib.h> | |
44 #include "umfpack.h" | |
45 | |
46 #define ABS(x) ((x) >= 0 ? (x) : -(x)) | |
47 | |
48 #define MAX(a,b) (((a) > (b)) ? (a) : (b)) | |
49 #ifndef TRUE | |
50 #define TRUE (1) | |
51 #endif | |
52 #ifndef FALSE | |
53 #define FALSE (0) | |
54 #endif | |
55 | |
56 /* -------------------------------------------------------------------------- */ | |
57 /* triplet form of the matrix. The triplets can be in any order. */ | |
58 /* -------------------------------------------------------------------------- */ | |
59 | |
60 static int n = 5, nz = 12 ; | |
61 static int Arow [ ] = { 0, 4, 1, 1, 2, 2, 0, 1, 2, 3, 4, 4} ; | |
62 static int Acol [ ] = { 0, 4, 0, 2, 1, 2, 1, 4, 3, 2, 1, 2} ; | |
63 static double Aval [ ] = {2., 1., 3., 4., -1., -3., 3., 6., 2., 1., 4., 2.} ; | |
64 static double b [ ] = {8., 45., -3., 3., 19.}, x [5], r [5] ; | |
65 | |
66 | |
67 /* -------------------------------------------------------------------------- */ | |
68 /* error: print a message and exit */ | |
69 /* -------------------------------------------------------------------------- */ | |
70 | |
71 static void error | |
72 ( | |
73 char *message | |
74 ) | |
75 { | |
76 printf ("\n\n====== error: %s =====\n\n", message) ; | |
77 exit (1) ; | |
78 } | |
79 | |
80 | |
81 /* -------------------------------------------------------------------------- */ | |
82 /* resid: compute the residual, r = Ax-b or r = A'x=b and return maxnorm (r) */ | |
83 /* -------------------------------------------------------------------------- */ | |
84 | |
85 static double resid | |
86 ( | |
87 int transpose, | |
88 int Ap [ ], | |
89 int Ai [ ], | |
90 double Ax [ ] | |
91 ) | |
92 { | |
93 int i, j, p ; | |
94 double norm ; | |
95 | |
96 for (i = 0 ; i < n ; i++) | |
97 { | |
98 r [i] = -b [i] ; | |
99 } | |
100 if (transpose) | |
101 { | |
102 for (j = 0 ; j < n ; j++) | |
103 { | |
104 for (p = Ap [j] ; p < Ap [j+1] ; p++) | |
105 { | |
106 i = Ai [p] ; | |
107 r [j] += Ax [p] * x [i] ; | |
108 } | |
109 } | |
110 } | |
111 else | |
112 { | |
113 for (j = 0 ; j < n ; j++) | |
114 { | |
115 for (p = Ap [j] ; p < Ap [j+1] ; p++) | |
116 { | |
117 i = Ai [p] ; | |
118 r [i] += Ax [p] * x [j] ; | |
119 } | |
120 } | |
121 } | |
122 norm = 0. ; | |
123 for (i = 0 ; i < n ; i++) | |
124 { | |
125 norm = MAX (ABS (r [i]), norm) ; | |
126 } | |
127 return (norm) ; | |
128 } | |
129 | |
130 | |
131 /* -------------------------------------------------------------------------- */ | |
132 /* main program */ | |
133 /* -------------------------------------------------------------------------- */ | |
134 | |
135 int main (int argc, char **argv) | |
136 { | |
137 double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL], *Ax, *Cx, *Lx, *Ux, | |
138 *W, t [2], *Dx, rnorm, *Rb, *y, *Rs ; | |
139 int *Ap, *Ai, *Cp, *Ci, row, col, p, lnz, unz, nr, nc, *Lp, *Li, *Ui, *Up, | |
140 *P, *Q, *Lj, i, j, k, anz, nfr, nchains, *Qinit, fnpiv, lnz1, unz1, nz1, | |
141 status, *Front_npivcol, *Front_parent, *Chain_start, *Wi, *Pinit, n1, | |
142 *Chain_maxrows, *Chain_maxcols, *Front_1strow, *Front_leftmostdesc, | |
143 nzud, do_recip ; | |
144 void *Symbolic, *Numeric ; | |
145 | |
146 /* ---------------------------------------------------------------------- */ | |
147 /* initializations */ | |
148 /* ---------------------------------------------------------------------- */ | |
149 | |
150 umfpack_tic (t) ; | |
151 | |
152 printf ("\n%s demo: _di_ version\n", UMFPACK_VERSION) ; | |
153 | |
154 /* get the default control parameters */ | |
155 umfpack_di_defaults (Control) ; | |
156 | |
157 /* change the default print level for this demo */ | |
158 /* (otherwise, nothing will print) */ | |
159 Control [UMFPACK_PRL] = 6 ; | |
160 | |
161 /* print the license agreement */ | |
162 umfpack_di_report_status (Control, UMFPACK_OK) ; | |
163 Control [UMFPACK_PRL] = 5 ; | |
164 | |
165 /* print the control parameters */ | |
166 umfpack_di_report_control (Control) ; | |
167 | |
168 /* ---------------------------------------------------------------------- */ | |
169 /* print A and b, and convert A to column-form */ | |
170 /* ---------------------------------------------------------------------- */ | |
171 | |
172 /* print the right-hand-side */ | |
173 printf ("\nb: ") ; | |
174 (void) umfpack_di_report_vector (n, b, Control) ; | |
175 | |
176 /* print the triplet form of the matrix */ | |
177 printf ("\nA: ") ; | |
178 (void) umfpack_di_report_triplet (n, n, nz, Arow, Acol, Aval, | |
179 Control) ; | |
180 | |
181 /* convert to column form */ | |
182 nz1 = MAX (nz,1) ; /* ensure arrays are not of size zero. */ | |
183 Ap = (int *) malloc ((n+1) * sizeof (int)) ; | |
184 Ai = (int *) malloc (nz1 * sizeof (int)) ; | |
185 Ax = (double *) malloc (nz1 * sizeof (double)) ; | |
186 if (!Ap || !Ai || !Ax) | |
187 { | |
188 error ("out of memory") ; | |
189 } | |
190 | |
191 status = umfpack_di_triplet_to_col (n, n, nz, Arow, Acol, Aval, | |
192 Ap, Ai, Ax, (int *) NULL) ; | |
193 | |
194 if (status < 0) | |
195 { | |
196 umfpack_di_report_status (Control, status) ; | |
197 error ("umfpack_di_triplet_to_col failed") ; | |
198 } | |
199 | |
200 /* print the column-form of A */ | |
201 printf ("\nA: ") ; | |
202 (void) umfpack_di_report_matrix (n, n, Ap, Ai, Ax, 1, Control) ; | |
203 | |
204 /* ---------------------------------------------------------------------- */ | |
205 /* symbolic factorization */ | |
206 /* ---------------------------------------------------------------------- */ | |
207 | |
208 status = umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, | |
209 Control, Info) ; | |
210 if (status < 0) | |
211 { | |
212 umfpack_di_report_info (Control, Info) ; | |
213 umfpack_di_report_status (Control, status) ; | |
214 error ("umfpack_di_symbolic failed") ; | |
215 } | |
216 | |
217 /* print the symbolic factorization */ | |
218 | |
219 printf ("\nSymbolic factorization of A: ") ; | |
220 (void) umfpack_di_report_symbolic (Symbolic, Control) ; | |
221 | |
222 /* ---------------------------------------------------------------------- */ | |
223 /* numeric factorization */ | |
224 /* ---------------------------------------------------------------------- */ | |
225 | |
226 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, | |
227 Control, Info) ; | |
228 if (status < 0) | |
229 { | |
230 umfpack_di_report_info (Control, Info) ; | |
231 umfpack_di_report_status (Control, status) ; | |
232 error ("umfpack_di_numeric failed") ; | |
233 } | |
234 | |
235 /* print the numeric factorization */ | |
236 printf ("\nNumeric factorization of A: ") ; | |
237 (void) umfpack_di_report_numeric (Numeric, Control) ; | |
238 | |
239 /* ---------------------------------------------------------------------- */ | |
240 /* solve Ax=b */ | |
241 /* ---------------------------------------------------------------------- */ | |
242 | |
243 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, | |
244 Numeric, Control, Info) ; | |
245 umfpack_di_report_info (Control, Info) ; | |
246 umfpack_di_report_status (Control, status) ; | |
247 if (status < 0) | |
248 { | |
249 error ("umfpack_di_solve failed") ; | |
250 } | |
251 printf ("\nx (solution of Ax=b): ") ; | |
252 (void) umfpack_di_report_vector (n, x, Control) ; | |
253 rnorm = resid (FALSE, Ap, Ai, Ax) ; | |
254 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
255 | |
256 /* ---------------------------------------------------------------------- */ | |
257 /* compute the determinant */ | |
258 /* ---------------------------------------------------------------------- */ | |
259 | |
260 status = umfpack_di_get_determinant (x, r, Numeric, Info) ; | |
261 umfpack_di_report_status (Control, status) ; | |
262 if (status < 0) | |
263 { | |
264 error ("umfpack_di_get_determinant failed") ; | |
265 } | |
266 printf ("determinant: (%g", x [0]) ; | |
267 printf (") * 10^(%g)\n", r [0]) ; | |
268 | |
269 /* ---------------------------------------------------------------------- */ | |
270 /* solve Ax=b, broken down into steps */ | |
271 /* ---------------------------------------------------------------------- */ | |
272 | |
273 /* Rb = R*b */ | |
274 Rb = (double *) malloc (n * sizeof (double)) ; | |
275 y = (double *) malloc (n * sizeof (double)) ; | |
276 if (!Rb || !y) error ("out of memory") ; | |
277 | |
278 status = umfpack_di_scale (Rb, b, Numeric) ; | |
279 if (status < 0) error ("umfpack_di_scale failed") ; | |
280 /* solve Ly = P*(Rb) */ | |
281 status = umfpack_di_solve (UMFPACK_Pt_L, Ap, Ai, Ax, y, Rb, | |
282 Numeric, Control, Info) ; | |
283 if (status < 0) error ("umfpack_di_solve failed") ; | |
284 /* solve UQ'x=y */ | |
285 status = umfpack_di_solve (UMFPACK_U_Qt, Ap, Ai, Ax, x, y, | |
286 Numeric, Control, Info) ; | |
287 if (status < 0) error ("umfpack_di_solve failed") ; | |
288 printf ("\nx (solution of Ax=b, solve is split into 3 steps): ") ; | |
289 (void) umfpack_di_report_vector (n, x, Control) ; | |
290 rnorm = resid (FALSE, Ap, Ai, Ax) ; | |
291 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
292 | |
293 free (Rb) ; | |
294 free (y) ; | |
295 | |
296 /* ---------------------------------------------------------------------- */ | |
297 /* solve A'x=b */ | |
298 /* ---------------------------------------------------------------------- */ | |
299 | |
300 status = umfpack_di_solve (UMFPACK_At, Ap, Ai, Ax, x, b, | |
301 Numeric, Control, Info) ; | |
302 umfpack_di_report_info (Control, Info) ; | |
303 if (status < 0) | |
304 { | |
305 error ("umfpack_di_solve failed") ; | |
306 } | |
307 printf ("\nx (solution of A'x=b): ") ; | |
308 (void) umfpack_di_report_vector (n, x, Control) ; | |
309 rnorm = resid (TRUE, Ap, Ai, Ax) ; | |
310 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
311 | |
312 /* ---------------------------------------------------------------------- */ | |
313 /* modify one numerical value in the column-form of A */ | |
314 /* ---------------------------------------------------------------------- */ | |
315 | |
316 /* change A (1,4), look for row index 1 in column 4. */ | |
317 row = 1 ; | |
318 col = 4 ; | |
319 for (p = Ap [col] ; p < Ap [col+1] ; p++) | |
320 { | |
321 if (row == Ai [p]) | |
322 { | |
323 printf ("\nchanging A (%d,%d) to zero\n", row, col) ; | |
324 Ax [p] = 0.0 ; | |
325 break ; | |
326 } | |
327 } | |
328 printf ("\nmodified A: ") ; | |
329 (void) umfpack_di_report_matrix (n, n, Ap, Ai, Ax, 1, Control) ; | |
330 | |
331 /* ---------------------------------------------------------------------- */ | |
332 /* redo the numeric factorization */ | |
333 /* ---------------------------------------------------------------------- */ | |
334 | |
335 /* The pattern (Ap and Ai) hasn't changed, so the symbolic factorization */ | |
336 /* doesn't have to be redone, no matter how much we change Ax. */ | |
337 | |
338 /* We don't need the Numeric object any more, so free it. */ | |
339 umfpack_di_free_numeric (&Numeric) ; | |
340 | |
341 /* Note that a memory leak would have occurred if the old Numeric */ | |
342 /* had not been free'd with umfpack_di_free_numeric above. */ | |
343 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, | |
344 Control, Info) ; | |
345 if (status < 0) | |
346 { | |
347 umfpack_di_report_info (Control, Info) ; | |
348 umfpack_di_report_status (Control, status) ; | |
349 error ("umfpack_di_numeric failed") ; | |
350 } | |
351 printf ("\nNumeric factorization of modified A: ") ; | |
352 (void) umfpack_di_report_numeric (Numeric, Control) ; | |
353 | |
354 /* ---------------------------------------------------------------------- */ | |
355 /* solve Ax=b, with the modified A */ | |
356 /* ---------------------------------------------------------------------- */ | |
357 | |
358 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, | |
359 Numeric, Control, Info) ; | |
360 umfpack_di_report_info (Control, Info) ; | |
361 if (status < 0) | |
362 { | |
363 umfpack_di_report_status (Control, status) ; | |
364 error ("umfpack_di_solve failed") ; | |
365 } | |
366 printf ("\nx (with modified A): ") ; | |
367 (void) umfpack_di_report_vector (n, x, Control) ; | |
368 rnorm = resid (FALSE, Ap, Ai, Ax) ; | |
369 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
370 | |
371 /* ---------------------------------------------------------------------- */ | |
372 /* modify all of the numerical values of A, but not the pattern */ | |
373 /* ---------------------------------------------------------------------- */ | |
374 | |
375 for (col = 0 ; col < n ; col++) | |
376 { | |
377 for (p = Ap [col] ; p < Ap [col+1] ; p++) | |
378 { | |
379 row = Ai [p] ; | |
380 printf ("changing ") ; | |
381 printf ("A (%d,%d) from %g", row, col, Ax [p]) ; | |
382 Ax [p] = Ax [p] + col*10 - row ; | |
383 printf (" to %g\n", Ax [p]) ; | |
384 } | |
385 } | |
386 printf ("\ncompletely modified A (same pattern): ") ; | |
387 (void) umfpack_di_report_matrix (n, n, Ap, Ai, Ax, 1, Control) ; | |
388 | |
389 /* ---------------------------------------------------------------------- */ | |
390 /* save the Symbolic object to file, free it, and load it back in */ | |
391 /* ---------------------------------------------------------------------- */ | |
392 | |
393 /* use the default filename, "symbolic.umf" */ | |
394 printf ("\nSaving symbolic object:\n") ; | |
395 status = umfpack_di_save_symbolic (Symbolic, (char *) NULL) ; | |
396 if (status < 0) | |
397 { | |
398 umfpack_di_report_status (Control, status) ; | |
399 error ("umfpack_di_save_symbolic failed") ; | |
400 } | |
401 printf ("\nFreeing symbolic object:\n") ; | |
402 umfpack_di_free_symbolic (&Symbolic) ; | |
403 printf ("\nLoading symbolic object:\n") ; | |
404 status = umfpack_di_load_symbolic (&Symbolic, (char *) NULL) ; | |
405 if (status < 0) | |
406 { | |
407 umfpack_di_report_status (Control, status) ; | |
408 error ("umfpack_di_load_symbolic failed") ; | |
409 } | |
410 printf ("\nDone loading symbolic object\n") ; | |
411 | |
412 /* ---------------------------------------------------------------------- */ | |
413 /* redo the numeric factorization */ | |
414 /* ---------------------------------------------------------------------- */ | |
415 | |
416 umfpack_di_free_numeric (&Numeric) ; | |
417 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, | |
418 Control, Info) ; | |
419 if (status < 0) | |
420 { | |
421 umfpack_di_report_info (Control, Info) ; | |
422 umfpack_di_report_status (Control, status) ; | |
423 error ("umfpack_di_numeric failed") ; | |
424 } | |
425 printf ("\nNumeric factorization of completely modified A: ") ; | |
426 (void) umfpack_di_report_numeric (Numeric, Control) ; | |
427 | |
428 /* ---------------------------------------------------------------------- */ | |
429 /* solve Ax=b, with the modified A */ | |
430 /* ---------------------------------------------------------------------- */ | |
431 | |
432 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, | |
433 Numeric, Control, Info) ; | |
434 umfpack_di_report_info (Control, Info) ; | |
435 if (status < 0) | |
436 { | |
437 umfpack_di_report_status (Control, status) ; | |
438 error ("umfpack_di_solve failed") ; | |
439 } | |
440 printf ("\nx (with completely modified A): ") ; | |
441 (void) umfpack_di_report_vector (n, x, Control) ; | |
442 rnorm = resid (FALSE, Ap, Ai, Ax) ; | |
443 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
444 | |
445 /* ---------------------------------------------------------------------- */ | |
446 /* free the symbolic and numeric factorization */ | |
447 /* ---------------------------------------------------------------------- */ | |
448 | |
449 umfpack_di_free_symbolic (&Symbolic) ; | |
450 umfpack_di_free_numeric (&Numeric) ; | |
451 | |
452 /* ---------------------------------------------------------------------- */ | |
453 /* C = transpose of A */ | |
454 /* ---------------------------------------------------------------------- */ | |
455 | |
456 Cp = (int *) malloc ((n+1) * sizeof (int)) ; | |
457 Ci = (int *) malloc (nz1 * sizeof (int)) ; | |
458 Cx = (double *) malloc (nz1 * sizeof (double)) ; | |
459 if (!Cp || !Ci || !Cx) | |
460 { | |
461 error ("out of memory") ; | |
462 } | |
463 status = umfpack_di_transpose (n, n, Ap, Ai, Ax, | |
464 (int *) NULL, (int *) NULL, Cp, Ci, Cx) ; | |
465 if (status < 0) | |
466 { | |
467 umfpack_di_report_status (Control, status) ; | |
468 error ("umfpack_di_transpose failed: ") ; | |
469 } | |
470 printf ("\nC (transpose of A): ") ; | |
471 (void) umfpack_di_report_matrix (n, n, Cp, Ci, Cx, 1, Control) ; | |
472 | |
473 /* ---------------------------------------------------------------------- */ | |
474 /* symbolic factorization of C */ | |
475 /* ---------------------------------------------------------------------- */ | |
476 | |
477 status = umfpack_di_symbolic (n, n, Cp, Ci, Cx, &Symbolic, | |
478 Control, Info) ; | |
479 if (status < 0) | |
480 { | |
481 umfpack_di_report_info (Control, Info) ; | |
482 umfpack_di_report_status (Control, status) ; | |
483 error ("umfpack_di_symbolic failed") ; | |
484 } | |
485 printf ("\nSymbolic factorization of C: ") ; | |
486 (void) umfpack_di_report_symbolic (Symbolic, Control) ; | |
487 | |
488 /* ---------------------------------------------------------------------- */ | |
489 /* copy the contents of Symbolic into user arrays print them */ | |
490 /* ---------------------------------------------------------------------- */ | |
491 | |
492 printf ("\nGet the contents of the Symbolic object for C:\n") ; | |
493 printf ("(compare with umfpack_di_report_symbolic output, above)\n") ; | |
494 Pinit = (int *) malloc ((n+1) * sizeof (int)) ; | |
495 Qinit = (int *) malloc ((n+1) * sizeof (int)) ; | |
496 Front_npivcol = (int *) malloc ((n+1) * sizeof (int)) ; | |
497 Front_1strow = (int *) malloc ((n+1) * sizeof (int)) ; | |
498 Front_leftmostdesc = (int *) malloc ((n+1) * sizeof (int)) ; | |
499 Front_parent = (int *) malloc ((n+1) * sizeof (int)) ; | |
500 Chain_start = (int *) malloc ((n+1) * sizeof (int)) ; | |
501 Chain_maxrows = (int *) malloc ((n+1) * sizeof (int)) ; | |
502 Chain_maxcols = (int *) malloc ((n+1) * sizeof (int)) ; | |
503 if (!Pinit || !Qinit || !Front_npivcol || !Front_parent || !Chain_start || | |
504 !Chain_maxrows || !Chain_maxcols || !Front_1strow || | |
505 !Front_leftmostdesc) | |
506 { | |
507 error ("out of memory") ; | |
508 } | |
509 | |
510 status = umfpack_di_get_symbolic (&nr, &nc, &n1, &anz, &nfr, &nchains, | |
511 Pinit, Qinit, Front_npivcol, Front_parent, Front_1strow, | |
512 Front_leftmostdesc, Chain_start, Chain_maxrows, Chain_maxcols, | |
513 Symbolic) ; | |
514 | |
515 if (status < 0) | |
516 { | |
517 error ("symbolic factorization invalid") ; | |
518 } | |
519 | |
520 printf ("From the Symbolic object, C is of dimension %d-by-%d\n", nr, nc); | |
521 printf (" with nz = %d, number of fronts = %d,\n", nz, nfr) ; | |
522 printf (" number of frontal matrix chains = %d\n", nchains) ; | |
523 | |
524 printf ("\nPivot columns in each front, and parent of each front:\n") ; | |
525 k = 0 ; | |
526 for (i = 0 ; i < nfr ; i++) | |
527 { | |
528 fnpiv = Front_npivcol [i] ; | |
529 printf (" Front %d: parent front: %d number of pivot cols: %d\n", | |
530 i, Front_parent [i], fnpiv) ; | |
531 for (j = 0 ; j < fnpiv ; j++) | |
532 { | |
533 col = Qinit [k] ; | |
534 printf ( | |
535 " %d-th pivot column is column %d in original matrix\n", | |
536 k, col) ; | |
537 k++ ; | |
538 } | |
539 } | |
540 | |
541 printf ("\nNote that the column ordering, above, will be refined\n") ; | |
542 printf ("in the numeric factorization below. The assignment of pivot\n") ; | |
543 printf ("columns to frontal matrices will always remain unchanged.\n") ; | |
544 | |
545 printf ("\nTotal number of pivot columns in frontal matrices: %d\n", k) ; | |
546 | |
547 printf ("\nFrontal matrix chains:\n") ; | |
548 for (j = 0 ; j < nchains ; j++) | |
549 { | |
550 printf (" Frontal matrices %d to %d are factorized in a single\n", | |
551 Chain_start [j], Chain_start [j+1] - 1) ; | |
552 printf (" working array of size %d-by-%d\n", | |
553 Chain_maxrows [j], Chain_maxcols [j]) ; | |
554 } | |
555 | |
556 /* ---------------------------------------------------------------------- */ | |
557 /* numeric factorization of C */ | |
558 /* ---------------------------------------------------------------------- */ | |
559 | |
560 status = umfpack_di_numeric (Cp, Ci, Cx, Symbolic, &Numeric, | |
561 Control, Info) ; | |
562 if (status < 0) | |
563 { | |
564 error ("umfpack_di_numeric failed") ; | |
565 } | |
566 printf ("\nNumeric factorization of C: ") ; | |
567 (void) umfpack_di_report_numeric (Numeric, Control) ; | |
568 | |
569 /* ---------------------------------------------------------------------- */ | |
570 /* extract the LU factors of C and print them */ | |
571 /* ---------------------------------------------------------------------- */ | |
572 | |
573 if (umfpack_di_get_lunz (&lnz, &unz, &nr, &nc, &nzud, Numeric) < 0) | |
574 { | |
575 error ("umfpack_di_get_lunz failed") ; | |
576 } | |
577 /* ensure arrays are not of zero size */ | |
578 lnz1 = MAX (lnz,1) ; | |
579 unz1 = MAX (unz,1) ; | |
580 Lp = (int *) malloc ((n+1) * sizeof (int)) ; | |
581 Lj = (int *) malloc (lnz1 * sizeof (int)) ; | |
582 Lx = (double *) malloc (lnz1 * sizeof (double)) ; | |
583 Up = (int *) malloc ((n+1) * sizeof (int)) ; | |
584 Ui = (int *) malloc (unz1 * sizeof (int)) ; | |
585 Ux = (double *) malloc (unz1 * sizeof (double)) ; | |
586 P = (int *) malloc (n * sizeof (int)) ; | |
587 Q = (int *) malloc (n * sizeof (int)) ; | |
588 Dx = (double *) NULL ; /* D vector not requested */ | |
589 Rs = (double *) malloc (n * sizeof (double)) ; | |
590 if (!Lp || !Lj || !Lx || !Up || !Ui || !Ux || !P || !Q || !Rs) | |
591 { | |
592 error ("out of memory") ; | |
593 } | |
594 status = umfpack_di_get_numeric (Lp, Lj, Lx, Up, Ui, Ux, | |
595 P, Q, Dx, &do_recip, Rs, Numeric) ; | |
596 if (status < 0) | |
597 { | |
598 error ("umfpack_di_get_numeric failed") ; | |
599 } | |
600 | |
601 printf ("\nL (lower triangular factor of C): ") ; | |
602 (void) umfpack_di_report_matrix (n, n, Lp, Lj, Lx, 0, Control) ; | |
603 printf ("\nU (upper triangular factor of C): ") ; | |
604 (void) umfpack_di_report_matrix (n, n, Up, Ui, Ux, 1, Control) ; | |
605 printf ("\nP: ") ; | |
606 (void) umfpack_di_report_perm (n, P, Control) ; | |
607 printf ("\nQ: ") ; | |
608 (void) umfpack_di_report_perm (n, Q, Control) ; | |
609 printf ("\nScale factors: row i of A is to be ") ; | |
610 if (do_recip) | |
611 { | |
612 printf ("multiplied by the ith scale factor\n") ; | |
613 } | |
614 else | |
615 { | |
616 printf ("divided by the ith scale factor\n") ; | |
617 } | |
618 for (i = 0 ; i < n ; i++) printf ("%d: %g\n", i, Rs [i]) ; | |
619 | |
620 /* ---------------------------------------------------------------------- */ | |
621 /* convert L to triplet form and print it */ | |
622 /* ---------------------------------------------------------------------- */ | |
623 | |
624 /* Note that L is in row-form, so it is the row indices that are created */ | |
625 /* by umfpack_di_col_to_triplet. */ | |
626 | |
627 printf ("\nConverting L to triplet form, and printing it:\n") ; | |
628 Li = (int *) malloc (lnz1 * sizeof (int)) ; | |
629 if (!Li) | |
630 { | |
631 error ("out of memory") ; | |
632 } | |
633 if (umfpack_di_col_to_triplet (n, Lp, Li) < 0) | |
634 { | |
635 error ("umfpack_di_col_to_triplet failed") ; | |
636 } | |
637 printf ("\nL, in triplet form: ") ; | |
638 (void) umfpack_di_report_triplet (n, n, lnz, Li, Lj, Lx, Control) ; | |
639 | |
640 /* ---------------------------------------------------------------------- */ | |
641 /* save the Numeric object to file, free it, and load it back in */ | |
642 /* ---------------------------------------------------------------------- */ | |
643 | |
644 /* use the default filename, "numeric.umf" */ | |
645 printf ("\nSaving numeric object:\n") ; | |
646 status = umfpack_di_save_numeric (Numeric, (char *) NULL) ; | |
647 if (status < 0) | |
648 { | |
649 umfpack_di_report_status (Control, status) ; | |
650 error ("umfpack_di_save_numeric failed") ; | |
651 } | |
652 printf ("\nFreeing numeric object:\n") ; | |
653 umfpack_di_free_numeric (&Numeric) ; | |
654 printf ("\nLoading numeric object:\n") ; | |
655 status = umfpack_di_load_numeric (&Numeric, (char *) NULL) ; | |
656 if (status < 0) | |
657 { | |
658 umfpack_di_report_status (Control, status) ; | |
659 error ("umfpack_di_load_numeric failed") ; | |
660 } | |
661 printf ("\nDone loading numeric object\n") ; | |
662 | |
663 /* ---------------------------------------------------------------------- */ | |
664 /* solve C'x=b */ | |
665 /* ---------------------------------------------------------------------- */ | |
666 | |
667 status = umfpack_di_solve (UMFPACK_At, Cp, Ci, Cx, x, b, | |
668 Numeric, Control, Info) ; | |
669 umfpack_di_report_info (Control, Info) ; | |
670 if (status < 0) | |
671 { | |
672 umfpack_di_report_status (Control, status) ; | |
673 error ("umfpack_di_solve failed") ; | |
674 } | |
675 printf ("\nx (solution of C'x=b): ") ; | |
676 (void) umfpack_di_report_vector (n, x, Control) ; | |
677 rnorm = resid (TRUE, Cp, Ci, Cx) ; | |
678 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
679 | |
680 /* ---------------------------------------------------------------------- */ | |
681 /* solve C'x=b again, using umfpack_di_wsolve instead */ | |
682 /* ---------------------------------------------------------------------- */ | |
683 | |
684 printf ("\nSolving C'x=b again, using umfpack_di_wsolve instead:\n") ; | |
685 Wi = (int *) malloc (n * sizeof (int)) ; | |
686 W = (double *) malloc (5*n * sizeof (double)) ; | |
687 if (!Wi || !W) | |
688 { | |
689 error ("out of memory") ; | |
690 } | |
691 | |
692 status = umfpack_di_wsolve (UMFPACK_At, Cp, Ci, Cx, x, b, | |
693 Numeric, Control, Info, Wi, W) ; | |
694 umfpack_di_report_info (Control, Info) ; | |
695 if (status < 0) | |
696 { | |
697 umfpack_di_report_status (Control, status) ; | |
698 error ("umfpack_di_wsolve failed") ; | |
699 } | |
700 printf ("\nx (solution of C'x=b): ") ; | |
701 (void) umfpack_di_report_vector (n, x, Control) ; | |
702 rnorm = resid (TRUE, Cp, Ci, Cx) ; | |
703 printf ("maxnorm of residual: %g\n\n", rnorm) ; | |
704 | |
705 /* ---------------------------------------------------------------------- */ | |
706 /* free everything */ | |
707 /* ---------------------------------------------------------------------- */ | |
708 | |
709 /* This is not strictly required since the process is exiting and the */ | |
710 /* system will reclaim the memory anyway. It's useful, though, just as */ | |
711 /* a list of what is currently malloc'ed by this program. Plus, it's */ | |
712 /* always a good habit to explicitly free whatever you malloc. */ | |
713 | |
714 free (Ap) ; | |
715 free (Ai) ; | |
716 free (Ax) ; | |
717 | |
718 free (Cp) ; | |
719 free (Ci) ; | |
720 free (Cx) ; | |
721 | |
722 free (Pinit) ; | |
723 free (Qinit) ; | |
724 free (Front_npivcol) ; | |
725 free (Front_1strow) ; | |
726 free (Front_leftmostdesc) ; | |
727 free (Front_parent) ; | |
728 free (Chain_start) ; | |
729 free (Chain_maxrows) ; | |
730 free (Chain_maxcols) ; | |
731 | |
732 free (Lp) ; | |
733 free (Lj) ; | |
734 free (Lx) ; | |
735 | |
736 free (Up) ; | |
737 free (Ui) ; | |
738 free (Ux) ; | |
739 | |
740 free (P) ; | |
741 free (Q) ; | |
742 | |
743 free (Li) ; | |
744 | |
745 free (Wi) ; | |
746 free (W) ; | |
747 | |
748 umfpack_di_free_symbolic (&Symbolic) ; | |
749 umfpack_di_free_numeric (&Numeric) ; | |
750 | |
751 /* ---------------------------------------------------------------------- */ | |
752 /* print the total time spent in this demo */ | |
753 /* ---------------------------------------------------------------------- */ | |
754 | |
755 umfpack_toc (t) ; | |
756 printf ("\numfpack_di_demo complete.\nTotal time: %5.2f seconds" | |
757 " (CPU time), %5.2f seconds (wallclock time)\n", t [1], t [0]) ; | |
758 return (0) ; | |
759 } |