Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/Demo/umfpack_dl_demo.c @ 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 /* ========================================================================== */
2 /* === umfpack_dl_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_dl_* 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_dl_symbolic does not need to be done. Refactorize (with
20 umfpack_dl_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 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 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 long transpose,
88 long Ap [ ],
89 long Ai [ ],
90 double Ax [ ]
91 )
92 {
93 long 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 long *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: _dl_ version\n", UMFPACK_VERSION) ;
153
154 /* get the default control parameters */
155 umfpack_dl_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_dl_report_status (Control, UMFPACK_OK) ;
163 Control [UMFPACK_PRL] = 5 ;
164
165 /* print the control parameters */
166 umfpack_dl_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_dl_report_vector (n, b, Control) ;
175
176 /* print the triplet form of the matrix */
177 printf ("\nA: ") ;
178 (void) umfpack_dl_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 = (long *) malloc ((n+1) * sizeof (long)) ;
184 Ai = (long *) malloc (nz1 * sizeof (long)) ;
185 Ax = (double *) malloc (nz1 * sizeof (double)) ;
186 if (!Ap || !Ai || !Ax)
187 {
188 error ("out of memory") ;
189 }
190
191 status = umfpack_dl_triplet_to_col (n, n, nz, Arow, Acol, Aval,
192 Ap, Ai, Ax, (long *) NULL) ;
193
194 if (status < 0)
195 {
196 umfpack_dl_report_status (Control, status) ;
197 error ("umfpack_dl_triplet_to_col failed") ;
198 }
199
200 /* print the column-form of A */
201 printf ("\nA: ") ;
202 (void) umfpack_dl_report_matrix (n, n, Ap, Ai, Ax, 1, Control) ;
203
204 /* ---------------------------------------------------------------------- */
205 /* symbolic factorization */
206 /* ---------------------------------------------------------------------- */
207
208 status = umfpack_dl_symbolic (n, n, Ap, Ai, Ax, &Symbolic,
209 Control, Info) ;
210 if (status < 0)
211 {
212 umfpack_dl_report_info (Control, Info) ;
213 umfpack_dl_report_status (Control, status) ;
214 error ("umfpack_dl_symbolic failed") ;
215 }
216
217 /* print the symbolic factorization */
218
219 printf ("\nSymbolic factorization of A: ") ;
220 (void) umfpack_dl_report_symbolic (Symbolic, Control) ;
221
222 /* ---------------------------------------------------------------------- */
223 /* numeric factorization */
224 /* ---------------------------------------------------------------------- */
225
226 status = umfpack_dl_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
227 Control, Info) ;
228 if (status < 0)
229 {
230 umfpack_dl_report_info (Control, Info) ;
231 umfpack_dl_report_status (Control, status) ;
232 error ("umfpack_dl_numeric failed") ;
233 }
234
235 /* print the numeric factorization */
236 printf ("\nNumeric factorization of A: ") ;
237 (void) umfpack_dl_report_numeric (Numeric, Control) ;
238
239 /* ---------------------------------------------------------------------- */
240 /* solve Ax=b */
241 /* ---------------------------------------------------------------------- */
242
243 status = umfpack_dl_solve (UMFPACK_A, Ap, Ai, Ax, x, b,
244 Numeric, Control, Info) ;
245 umfpack_dl_report_info (Control, Info) ;
246 umfpack_dl_report_status (Control, status) ;
247 if (status < 0)
248 {
249 error ("umfpack_dl_solve failed") ;
250 }
251 printf ("\nx (solution of Ax=b): ") ;
252 (void) umfpack_dl_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_dl_get_determinant (x, r, Numeric, Info) ;
261 umfpack_dl_report_status (Control, status) ;
262 if (status < 0)
263 {
264 error ("umfpack_dl_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_dl_scale (Rb, b, Numeric) ;
279 if (status < 0) error ("umfpack_dl_scale failed") ;
280 /* solve Ly = P*(Rb) */
281 status = umfpack_dl_solve (UMFPACK_Pt_L, Ap, Ai, Ax, y, Rb,
282 Numeric, Control, Info) ;
283 if (status < 0) error ("umfpack_dl_solve failed") ;
284 /* solve UQ'x=y */
285 status = umfpack_dl_solve (UMFPACK_U_Qt, Ap, Ai, Ax, x, y,
286 Numeric, Control, Info) ;
287 if (status < 0) error ("umfpack_dl_solve failed") ;
288 printf ("\nx (solution of Ax=b, solve is split into 3 steps): ") ;
289 (void) umfpack_dl_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_dl_solve (UMFPACK_At, Ap, Ai, Ax, x, b,
301 Numeric, Control, Info) ;
302 umfpack_dl_report_info (Control, Info) ;
303 if (status < 0)
304 {
305 error ("umfpack_dl_solve failed") ;
306 }
307 printf ("\nx (solution of A'x=b): ") ;
308 (void) umfpack_dl_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 (%ld,%ld) to zero\n", row, col) ;
324 Ax [p] = 0.0 ;
325 break ;
326 }
327 }
328 printf ("\nmodified A: ") ;
329 (void) umfpack_dl_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_dl_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_dl_free_numeric above. */
343 status = umfpack_dl_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
344 Control, Info) ;
345 if (status < 0)
346 {
347 umfpack_dl_report_info (Control, Info) ;
348 umfpack_dl_report_status (Control, status) ;
349 error ("umfpack_dl_numeric failed") ;
350 }
351 printf ("\nNumeric factorization of modified A: ") ;
352 (void) umfpack_dl_report_numeric (Numeric, Control) ;
353
354 /* ---------------------------------------------------------------------- */
355 /* solve Ax=b, with the modified A */
356 /* ---------------------------------------------------------------------- */
357
358 status = umfpack_dl_solve (UMFPACK_A, Ap, Ai, Ax, x, b,
359 Numeric, Control, Info) ;
360 umfpack_dl_report_info (Control, Info) ;
361 if (status < 0)
362 {
363 umfpack_dl_report_status (Control, status) ;
364 error ("umfpack_dl_solve failed") ;
365 }
366 printf ("\nx (with modified A): ") ;
367 (void) umfpack_dl_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 (%ld,%ld) 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_dl_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_dl_save_symbolic (Symbolic, (char *) NULL) ;
396 if (status < 0)
397 {
398 umfpack_dl_report_status (Control, status) ;
399 error ("umfpack_dl_save_symbolic failed") ;
400 }
401 printf ("\nFreeing symbolic object:\n") ;
402 umfpack_dl_free_symbolic (&Symbolic) ;
403 printf ("\nLoading symbolic object:\n") ;
404 status = umfpack_dl_load_symbolic (&Symbolic, (char *) NULL) ;
405 if (status < 0)
406 {
407 umfpack_dl_report_status (Control, status) ;
408 error ("umfpack_dl_load_symbolic failed") ;
409 }
410 printf ("\nDone loading symbolic object\n") ;
411
412 /* ---------------------------------------------------------------------- */
413 /* redo the numeric factorization */
414 /* ---------------------------------------------------------------------- */
415
416 umfpack_dl_free_numeric (&Numeric) ;
417 status = umfpack_dl_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
418 Control, Info) ;
419 if (status < 0)
420 {
421 umfpack_dl_report_info (Control, Info) ;
422 umfpack_dl_report_status (Control, status) ;
423 error ("umfpack_dl_numeric failed") ;
424 }
425 printf ("\nNumeric factorization of completely modified A: ") ;
426 (void) umfpack_dl_report_numeric (Numeric, Control) ;
427
428 /* ---------------------------------------------------------------------- */
429 /* solve Ax=b, with the modified A */
430 /* ---------------------------------------------------------------------- */
431
432 status = umfpack_dl_solve (UMFPACK_A, Ap, Ai, Ax, x, b,
433 Numeric, Control, Info) ;
434 umfpack_dl_report_info (Control, Info) ;
435 if (status < 0)
436 {
437 umfpack_dl_report_status (Control, status) ;
438 error ("umfpack_dl_solve failed") ;
439 }
440 printf ("\nx (with completely modified A): ") ;
441 (void) umfpack_dl_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_dl_free_symbolic (&Symbolic) ;
450 umfpack_dl_free_numeric (&Numeric) ;
451
452 /* ---------------------------------------------------------------------- */
453 /* C = transpose of A */
454 /* ---------------------------------------------------------------------- */
455
456 Cp = (long *) malloc ((n+1) * sizeof (long)) ;
457 Ci = (long *) malloc (nz1 * sizeof (long)) ;
458 Cx = (double *) malloc (nz1 * sizeof (double)) ;
459 if (!Cp || !Ci || !Cx)
460 {
461 error ("out of memory") ;
462 }
463 status = umfpack_dl_transpose (n, n, Ap, Ai, Ax,
464 (long *) NULL, (long *) NULL, Cp, Ci, Cx) ;
465 if (status < 0)
466 {
467 umfpack_dl_report_status (Control, status) ;
468 error ("umfpack_dl_transpose failed: ") ;
469 }
470 printf ("\nC (transpose of A): ") ;
471 (void) umfpack_dl_report_matrix (n, n, Cp, Ci, Cx, 1, Control) ;
472
473 /* ---------------------------------------------------------------------- */
474 /* symbolic factorization of C */
475 /* ---------------------------------------------------------------------- */
476
477 status = umfpack_dl_symbolic (n, n, Cp, Ci, Cx, &Symbolic,
478 Control, Info) ;
479 if (status < 0)
480 {
481 umfpack_dl_report_info (Control, Info) ;
482 umfpack_dl_report_status (Control, status) ;
483 error ("umfpack_dl_symbolic failed") ;
484 }
485 printf ("\nSymbolic factorization of C: ") ;
486 (void) umfpack_dl_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_dl_report_symbolic output, above)\n") ;
494 Pinit = (long *) malloc ((n+1) * sizeof (long)) ;
495 Qinit = (long *) malloc ((n+1) * sizeof (long)) ;
496 Front_npivcol = (long *) malloc ((n+1) * sizeof (long)) ;
497 Front_1strow = (long *) malloc ((n+1) * sizeof (long)) ;
498 Front_leftmostdesc = (long *) malloc ((n+1) * sizeof (long)) ;
499 Front_parent = (long *) malloc ((n+1) * sizeof (long)) ;
500 Chain_start = (long *) malloc ((n+1) * sizeof (long)) ;
501 Chain_maxrows = (long *) malloc ((n+1) * sizeof (long)) ;
502 Chain_maxcols = (long *) malloc ((n+1) * sizeof (long)) ;
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_dl_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 %ld-by-%ld\n", nr, nc);
521 printf (" with nz = %ld, number of fronts = %ld,\n", nz, nfr) ;
522 printf (" number of frontal matrix chains = %ld\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 %ld: parent front: %ld number of pivot cols: %ld\n",
530 i, Front_parent [i], fnpiv) ;
531 for (j = 0 ; j < fnpiv ; j++)
532 {
533 col = Qinit [k] ;
534 printf (
535 " %ld-th pivot column is column %ld 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: %ld\n", k) ;
546
547 printf ("\nFrontal matrix chains:\n") ;
548 for (j = 0 ; j < nchains ; j++)
549 {
550 printf (" Frontal matrices %ld to %ld are factorized in a single\n",
551 Chain_start [j], Chain_start [j+1] - 1) ;
552 printf (" working array of size %ld-by-%ld\n",
553 Chain_maxrows [j], Chain_maxcols [j]) ;
554 }
555
556 /* ---------------------------------------------------------------------- */
557 /* numeric factorization of C */
558 /* ---------------------------------------------------------------------- */
559
560 status = umfpack_dl_numeric (Cp, Ci, Cx, Symbolic, &Numeric,
561 Control, Info) ;
562 if (status < 0)
563 {
564 error ("umfpack_dl_numeric failed") ;
565 }
566 printf ("\nNumeric factorization of C: ") ;
567 (void) umfpack_dl_report_numeric (Numeric, Control) ;
568
569 /* ---------------------------------------------------------------------- */
570 /* extract the LU factors of C and print them */
571 /* ---------------------------------------------------------------------- */
572
573 if (umfpack_dl_get_lunz (&lnz, &unz, &nr, &nc, &nzud, Numeric) < 0)
574 {
575 error ("umfpack_dl_get_lunz failed") ;
576 }
577 /* ensure arrays are not of zero size */
578 lnz1 = MAX (lnz,1) ;
579 unz1 = MAX (unz,1) ;
580 Lp = (long *) malloc ((n+1) * sizeof (long)) ;
581 Lj = (long *) malloc (lnz1 * sizeof (long)) ;
582 Lx = (double *) malloc (lnz1 * sizeof (double)) ;
583 Up = (long *) malloc ((n+1) * sizeof (long)) ;
584 Ui = (long *) malloc (unz1 * sizeof (long)) ;
585 Ux = (double *) malloc (unz1 * sizeof (double)) ;
586 P = (long *) malloc (n * sizeof (long)) ;
587 Q = (long *) malloc (n * sizeof (long)) ;
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_dl_get_numeric (Lp, Lj, Lx, Up, Ui, Ux,
595 P, Q, Dx, &do_recip, Rs, Numeric) ;
596 if (status < 0)
597 {
598 error ("umfpack_dl_get_numeric failed") ;
599 }
600
601 printf ("\nL (lower triangular factor of C): ") ;
602 (void) umfpack_dl_report_matrix (n, n, Lp, Lj, Lx, 0, Control) ;
603 printf ("\nU (upper triangular factor of C): ") ;
604 (void) umfpack_dl_report_matrix (n, n, Up, Ui, Ux, 1, Control) ;
605 printf ("\nP: ") ;
606 (void) umfpack_dl_report_perm (n, P, Control) ;
607 printf ("\nQ: ") ;
608 (void) umfpack_dl_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 ("%ld: %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_dl_col_to_triplet. */
626
627 printf ("\nConverting L to triplet form, and printing it:\n") ;
628 Li = (long *) malloc (lnz1 * sizeof (long)) ;
629 if (!Li)
630 {
631 error ("out of memory") ;
632 }
633 if (umfpack_dl_col_to_triplet (n, Lp, Li) < 0)
634 {
635 error ("umfpack_dl_col_to_triplet failed") ;
636 }
637 printf ("\nL, in triplet form: ") ;
638 (void) umfpack_dl_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_dl_save_numeric (Numeric, (char *) NULL) ;
647 if (status < 0)
648 {
649 umfpack_dl_report_status (Control, status) ;
650 error ("umfpack_dl_save_numeric failed") ;
651 }
652 printf ("\nFreeing numeric object:\n") ;
653 umfpack_dl_free_numeric (&Numeric) ;
654 printf ("\nLoading numeric object:\n") ;
655 status = umfpack_dl_load_numeric (&Numeric, (char *) NULL) ;
656 if (status < 0)
657 {
658 umfpack_dl_report_status (Control, status) ;
659 error ("umfpack_dl_load_numeric failed") ;
660 }
661 printf ("\nDone loading numeric object\n") ;
662
663 /* ---------------------------------------------------------------------- */
664 /* solve C'x=b */
665 /* ---------------------------------------------------------------------- */
666
667 status = umfpack_dl_solve (UMFPACK_At, Cp, Ci, Cx, x, b,
668 Numeric, Control, Info) ;
669 umfpack_dl_report_info (Control, Info) ;
670 if (status < 0)
671 {
672 umfpack_dl_report_status (Control, status) ;
673 error ("umfpack_dl_solve failed") ;
674 }
675 printf ("\nx (solution of C'x=b): ") ;
676 (void) umfpack_dl_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_dl_wsolve instead */
682 /* ---------------------------------------------------------------------- */
683
684 printf ("\nSolving C'x=b again, using umfpack_dl_wsolve instead:\n") ;
685 Wi = (long *) malloc (n * sizeof (long)) ;
686 W = (double *) malloc (5*n * sizeof (double)) ;
687 if (!Wi || !W)
688 {
689 error ("out of memory") ;
690 }
691
692 status = umfpack_dl_wsolve (UMFPACK_At, Cp, Ci, Cx, x, b,
693 Numeric, Control, Info, Wi, W) ;
694 umfpack_dl_report_info (Control, Info) ;
695 if (status < 0)
696 {
697 umfpack_dl_report_status (Control, status) ;
698 error ("umfpack_dl_wsolve failed") ;
699 }
700 printf ("\nx (solution of C'x=b): ") ;
701 (void) umfpack_dl_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_dl_free_symbolic (&Symbolic) ;
749 umfpack_dl_free_numeric (&Numeric) ;
750
751 /* ---------------------------------------------------------------------- */
752 /* print the total time spent in this demo */
753 /* ---------------------------------------------------------------------- */
754
755 umfpack_toc (t) ;
756 printf ("\numfpack_dl_demo complete.\nTotal time: %5.2f seconds"
757 " (CPU time), %5.2f seconds (wallclock time)\n", t [1], t [0]) ;
758 return (0) ;
759 }