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