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