1 /* ========================================================================== */
2 /* === Demo/cholmod_l_demo ================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Demo Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Read in a matrix from a file, and use CHOLMOD to solve Ax=b if A is
10 * symmetric, or (AA'+beta*I)x=b otherwise. The file format is a simple
11 * triplet format, compatible with most files in the Matrix Market format.
12 * See cholmod_read.c for more details. The readhb.f program reads a
13 * Harwell/Boeing matrix (excluding element-types) and converts it into the
14 * form needed by this program. reade.f reads a matrix in Harwell/Boeing
15 * finite-element form.
16 *
17 * Usage:
18 * cholmod_l_demo matrixfile
19 * cholmod_l_demo < matrixfile
20 *
21 * The matrix is assumed to be positive definite (a supernodal LL' or simplicial
22 * LDL' factorization is used).
23 *
24 * Requires the Core, Cholesky, MatrixOps, and Check Modules.
25 * Optionally uses the Partition and Supernodal Modules.
26 * Does not use the Modify Module.
27 *
28 * See cholmod_simple.c for a simpler demo program.
29 *
30 * SuiteSparse_long is normally defined as long, except for WIN64.
31 */
32
33 #include "cholmod_demo.h"
34 #define NTRIALS 100
35
36 /* ff is a global variable so that it can be closed by my_handler */
37 FILE *ff ;
38
39 /* halt if an error occurs */
my_handler(int status,const char * file,int line,const char * message)40 static void my_handler (int status, const char *file, int line,
41 const char *message)
42 {
43 printf ("cholmod error: file: %s line: %d status: %d: %s\n",
44 file, line, status, message) ;
45 if (status < 0)
46 {
47 if (ff != NULL) fclose (ff) ;
48 exit (0) ;
49 }
50 }
51
main(int argc,char ** argv)52 int main (int argc, char **argv)
53 {
54 double resid [4], t, ta, tf, ts [3], tot, bnorm, xnorm, anorm, rnorm, fl,
55 anz, axbnorm, rnorm2, resid2, rcond ;
56 FILE *f ;
57 cholmod_sparse *A ;
58 cholmod_dense *X = NULL, *B, *W, *R = NULL ;
59 double one [2], zero [2], minusone [2], beta [2], xlnz ;
60 cholmod_common Common, *cm ;
61 cholmod_factor *L ;
62 double *Bx, *Rx, *Xx, *Bz, *Xz, *Rz ;
63 SuiteSparse_long i, n, isize, xsize, ordering, xtype, s, ss, lnz ;
64 int trial, method, L_is_super ;
65 int ver [3] ;
66 int prefer_zomplex, nmethods ;
67
68 ts[0] = 0.;
69 ts[1] = 0.;
70 ts[2] = 0.;
71
72 /* ---------------------------------------------------------------------- */
73 /* get the file containing the input matrix */
74 /* ---------------------------------------------------------------------- */
75
76 ff = NULL ;
77 prefer_zomplex = 0 ;
78 if (argc > 1)
79 {
80 if ((f = fopen (argv [1], "r")) == NULL)
81 {
82 my_handler (CHOLMOD_INVALID, __FILE__, __LINE__,
83 "unable to open file") ;
84 }
85 ff = f ;
86 prefer_zomplex = (argc > 2) ;
87 }
88 else
89 {
90 f = stdin ;
91 }
92
93 /* ---------------------------------------------------------------------- */
94 /* start CHOLMOD and set parameters */
95 /* ---------------------------------------------------------------------- */
96
97 cm = &Common ;
98 cholmod_l_start (cm) ;
99 CHOLMOD_FUNCTION_DEFAULTS ; /* just for testing (not required) */
100
101 /* cm->useGPU = 1; */
102 cm->prefer_zomplex = prefer_zomplex ;
103
104 /* use default parameter settings, except for the error handler. This
105 * demo program terminates if an error occurs (out of memory, not positive
106 * definite, ...). It makes the demo program simpler (no need to check
107 * CHOLMOD error conditions). This non-default parameter setting has no
108 * effect on performance. */
109 cm->error_handler = my_handler ;
110
111 /* Note that CHOLMOD will do a supernodal LL' or a simplicial LDL' by
112 * default, automatically selecting the latter if flop/nnz(L) < 40. */
113
114 /* ---------------------------------------------------------------------- */
115 /* create basic scalars */
116 /* ---------------------------------------------------------------------- */
117
118 zero [0] = 0 ;
119 zero [1] = 0 ;
120 one [0] = 1 ;
121 one [1] = 0 ;
122 minusone [0] = -1 ;
123 minusone [1] = 0 ;
124 beta [0] = 1e-6 ;
125 beta [1] = 0 ;
126
127 /* ---------------------------------------------------------------------- */
128 /* read in a matrix */
129 /* ---------------------------------------------------------------------- */
130
131 printf ("\n---------------------------------- cholmod_l_demo:\n") ;
132 cholmod_l_version (ver) ;
133 printf ("cholmod version %d.%d.%d\n", ver [0], ver [1], ver [2]) ;
134 SuiteSparse_version (ver) ;
135 printf ("SuiteSparse version %d.%d.%d\n", ver [0], ver [1], ver [2]) ;
136 A = cholmod_l_read_sparse (f, cm) ;
137 if (ff != NULL)
138 {
139 fclose (ff) ;
140 ff = NULL ;
141 }
142 anorm = 1 ;
143 #ifndef NMATRIXOPS
144 anorm = cholmod_l_norm_sparse (A, 0, cm) ;
145 printf ("norm (A,inf) = %g\n", anorm) ;
146 printf ("norm (A,1) = %g\n", cholmod_l_norm_sparse (A, 1, cm)) ;
147 #endif
148
149 if (prefer_zomplex && A->xtype == CHOLMOD_COMPLEX)
150 {
151 /* Convert to zomplex, just for testing. In a zomplex matrix,
152 the real and imaginary parts are in separate arrays. MATLAB
153 uses zomplex matrix exclusively. */
154 double *Ax = A->x ;
155 SuiteSparse_long nz = cholmod_l_nnz (A, cm) ;
156 printf ("nz: %ld\n", nz) ;
157 double *Ax2 = cholmod_l_malloc (nz, sizeof (double), cm) ;
158 double *Az2 = cholmod_l_malloc (nz, sizeof (double), cm) ;
159 for (i = 0 ; i < nz ; i++)
160 {
161 Ax2 [i] = Ax [2*i ] ;
162 Az2 [i] = Ax [2*i+1] ;
163 }
164 cholmod_l_free (A->nzmax, 2*sizeof(double), Ax, cm) ;
165 A->x = Ax2 ;
166 A->z = Az2 ;
167 A->xtype = CHOLMOD_ZOMPLEX ;
168 /* cm->print = 5 ; */
169 }
170
171 xtype = A->xtype ;
172 cholmod_l_print_sparse (A, "A", cm) ;
173
174 #if 0
175 if ( 0 ) {
176 // scale diagonal
177 printf ("\n\n SCALING DIAGONAL \n\n");
178
179 // create diagonal
180 printf ("%ld,%ld,%d\n", A->nrow, A->ncol, A->xtype );
181
182 cholmod_sparse *D = cholmod_l_speye (A->nrow, A->ncol, A->xtype, cm );
183 printf ("sparse done \n");
184 cholmod_l_print_sparse (D, "D", cm);
185
186 D->stype = 1;
187 cholmod_l_print_sparse (D, "D", cm);
188
189 double alpha[2];
190 double beta[2];
191 alpha[0] = 1.0;
192 alpha[1] = 1.0;
193 beta[0] = 1.0e9; // 9 works, 467doesn't
194 beta[1] = 1.0e0;
195
196 cholmod_sparse *C = cholmod_l_add (A, D, alpha, beta, 1, 0, cm );
197 cholmod_l_print_sparse (C, "C", cm);
198
199 A = C;
200
201 }
202 #endif
203
204 if (A->nrow > A->ncol)
205 {
206 /* Transpose A so that A'A+beta*I will be factorized instead */
207 cholmod_sparse *C = cholmod_l_transpose (A, 2, cm) ;
208 cholmod_l_free_sparse (&A, cm) ;
209 A = C ;
210 printf ("transposing input matrix\n") ;
211 }
212
213 /* ---------------------------------------------------------------------- */
214 /* create an arbitrary right-hand-side */
215 /* ---------------------------------------------------------------------- */
216
217 n = A->nrow ;
218 B = cholmod_l_zeros (n, 1, xtype, cm) ;
219 Bx = B->x ;
220 Bz = B->z ;
221
222 #if GHS
223 {
224 /* b = A*ones(n,1), used by Gould, Hu, and Scott in their experiments */
225 cholmod_dense *X0 ;
226 X0 = cholmod_l_ones (A->ncol, 1, xtype, cm) ;
227 cholmod_l_sdmult (A, 0, one, zero, X0, B, cm) ;
228 cholmod_l_free_dense (&X0, cm) ;
229 }
230 #else
231 if (xtype == CHOLMOD_REAL)
232 {
233 /* real case */
234 for (i = 0 ; i < n ; i++)
235 {
236 double x = n ;
237 Bx [i] = 1 + i / x ;
238 }
239 }
240 else if (xtype == CHOLMOD_COMPLEX)
241 {
242 /* complex case */
243 for (i = 0 ; i < n ; i++)
244 {
245 double x = n ;
246 Bx [2*i ] = 1 + i / x ; /* real part of B(i) */
247 Bx [2*i+1] = (x/2 - i) / (3*x) ; /* imag part of B(i) */
248 }
249 }
250 else /* (xtype == CHOLMOD_ZOMPLEX) */
251 {
252 /* zomplex case */
253 for (i = 0 ; i < n ; i++)
254 {
255 double x = n ;
256 Bx [i] = 1 + i / x ; /* real part of B(i) */
257 Bz [i] = (x/2 - i) / (3*x) ; /* imag part of B(i) */
258 }
259 }
260
261 #endif
262
263 cholmod_l_print_dense (B, "B", cm) ;
264 bnorm = 1 ;
265 #ifndef NMATRIXOPS
266 bnorm = cholmod_l_norm_dense (B, 0, cm) ; /* max norm */
267 printf ("bnorm %g\n", bnorm) ;
268 #endif
269
270 /* ---------------------------------------------------------------------- */
271 /* analyze and factorize */
272 /* ---------------------------------------------------------------------- */
273
274 t = CPUTIME ;
275 L = cholmod_l_analyze (A, cm) ;
276 ta = CPUTIME - t ;
277 ta = MAX (ta, 0) ;
278
279 printf ("Analyze: flop %g lnz %g\n", cm->fl, cm->lnz) ;
280
281 if (A->stype == 0)
282 {
283 printf ("Factorizing A*A'+beta*I\n") ;
284 t = CPUTIME ;
285 cholmod_l_factorize_p (A, beta, NULL, 0, L, cm) ;
286 tf = CPUTIME - t ;
287 tf = MAX (tf, 0) ;
288 }
289 else
290 {
291 printf ("Factorizing A\n") ;
292 t = CPUTIME ;
293 cholmod_l_factorize (A, L, cm) ;
294 tf = CPUTIME - t ;
295 tf = MAX (tf, 0) ;
296 }
297
298 cholmod_l_print_factor (L, "L", cm) ;
299
300 /* determine the # of integers's and reals's in L. See cholmod_free */
301 if (L->is_super)
302 {
303 s = L->nsuper + 1 ;
304 xsize = L->xsize ;
305 ss = L->ssize ;
306 isize =
307 n /* L->Perm */
308 + n /* L->ColCount, nz in each column of 'pure' L */
309 + s /* L->pi, column pointers for L->s */
310 + s /* L->px, column pointers for L->x */
311 + s /* L->super, starting column index of each supernode */
312 + ss ; /* L->s, the pattern of the supernodes */
313 }
314 else
315 {
316 /* this space can increase if you change parameters to their non-
317 * default values (cm->final_pack, for example). */
318 lnz = L->nzmax ;
319 xsize = lnz ;
320 isize =
321 n /* L->Perm */
322 + n /* L->ColCount, nz in each column of 'pure' L */
323 + n+1 /* L->p, column pointers */
324 + lnz /* L->i, integer row indices */
325 + n /* L->nz, nz in each column of L */
326 + n+2 /* L->next, link list */
327 + n+2 ; /* L->prev, link list */
328 }
329
330 /* solve with Bset will change L from simplicial to supernodal */
331 rcond = cholmod_l_rcond (L, cm) ;
332 L_is_super = L->is_super ;
333
334 /* ---------------------------------------------------------------------- */
335 /* solve */
336 /* ---------------------------------------------------------------------- */
337
338 if (n >= 1000)
339 {
340 nmethods = 1 ;
341 }
342 else if (xtype == CHOLMOD_ZOMPLEX)
343 {
344 nmethods = 2 ;
345 }
346 else
347 {
348 nmethods = 3 ;
349 }
350 printf ("nmethods: %d\n", nmethods) ;
351
352 for (method = 0 ; method <= nmethods ; method++)
353 {
354 double x = n ;
355 resid [method] = -1 ; /* not yet computed */
356
357 if (method == 0)
358 {
359 /* basic solve, just once */
360 t = CPUTIME ;
361 X = cholmod_l_solve (CHOLMOD_A, L, B, cm) ;
362 ts [0] = CPUTIME - t ;
363 ts [0] = MAX (ts [0], 0) ;
364 }
365 else if (method == 1)
366 {
367 /* basic solve, many times, but keep the last one */
368 t = CPUTIME ;
369 for (trial = 0 ; trial < NTRIALS ; trial++)
370 {
371 cholmod_l_free_dense (&X, cm) ;
372 Bx [0] = 1 + trial / x ; /* tweak B each iteration */
373 X = cholmod_l_solve (CHOLMOD_A, L, B, cm) ;
374 }
375 ts [1] = CPUTIME - t ;
376 ts [1] = MAX (ts [1], 0) / NTRIALS ;
377 }
378 else if (method == 2)
379 {
380 /* solve with reused workspace */
381 cholmod_dense *Ywork = NULL, *Ework = NULL ;
382 cholmod_l_free_dense (&X, cm) ;
383
384 t = CPUTIME ;
385 for (trial = 0 ; trial < NTRIALS ; trial++)
386 {
387 Bx [0] = 1 + trial / x ; /* tweak B each iteration */
388 cholmod_l_solve2 (CHOLMOD_A, L, B, NULL, &X, NULL,
389 &Ywork, &Ework, cm) ;
390 }
391 cholmod_l_free_dense (&Ywork, cm) ;
392 cholmod_l_free_dense (&Ework, cm) ;
393 ts [2] = CPUTIME - t ;
394 ts [2] = MAX (ts [2], 0) / NTRIALS ;
395
396 }
397 else
398 {
399 /* solve with reused workspace and sparse Bset */
400 cholmod_dense *Ywork = NULL, *Ework = NULL ;
401 cholmod_dense *X2 = NULL, *B2 = NULL ;
402 cholmod_sparse *Bset, *Xset = NULL ;
403 SuiteSparse_long *Bsetp, *Bseti, *Xsetp, *Xseti, xlen, j, k, *Lnz ;
404 double *X1x, *X2x, *B2x, err ;
405 FILE *timelog = fopen ("timelog.m", "w") ;
406 if (timelog) fprintf (timelog, "results = [\n") ;
407
408 B2 = cholmod_l_zeros (n, 1, xtype, cm) ;
409 B2x = B2->x ;
410
411 Bset = cholmod_l_allocate_sparse (n, 1, 1, FALSE, TRUE, 0,
412 CHOLMOD_PATTERN, cm) ;
413 Bsetp = Bset->p ;
414 Bseti = Bset->i ;
415 Bsetp [0] = 0 ; /* nnz(B) is 1 (it can be anything) */
416 Bsetp [1] = 1 ;
417 resid [3] = 0 ;
418
419 for (i = 0 ; i < MIN (100,n) ; i++)
420 {
421 /* B (i) is nonzero, all other entries are ignored
422 (implied to be zero) */
423 Bseti [0] = i ;
424 if (xtype == CHOLMOD_REAL)
425 {
426 B2x [i] = 3.1 * i + 0.9 ;
427 }
428 else /* (xtype == CHOLMOD_COMPLEX) */
429 {
430 B2x [2*i ] = i + 0.042 ;
431 B2x [2*i+1] = i - 92.7 ;
432 }
433
434 /* first get the entire solution, to compare against */
435 cholmod_l_solve2 (CHOLMOD_A, L, B2, NULL, &X, NULL,
436 &Ywork, &Ework, cm) ;
437
438 /* now get the sparse solutions; this will change L from
439 supernodal to simplicial */
440
441 if (i == 0)
442 {
443 /* first solve can be slower because it has to allocate
444 space for X2, Xset, etc, and change L.
445 So don't time it */
446 cholmod_l_solve2 (CHOLMOD_A, L, B2, Bset, &X2, &Xset,
447 &Ywork, &Ework, cm) ;
448 }
449
450 t = CPUTIME ;
451 for (trial = 0 ; trial < NTRIALS ; trial++)
452 {
453 /* solve Ax=b but only to get x(i).
454 b is all zero except for b(i).
455 This takes O(xlen) time */
456 cholmod_l_solve2 (CHOLMOD_A, L, B2, Bset, &X2, &Xset,
457 &Ywork, &Ework, cm) ;
458 }
459 t = CPUTIME - t ;
460 t = MAX (t, 0) / NTRIALS ;
461
462 /* check the solution and log the time */
463 Xsetp = Xset->p ;
464 Xseti = Xset->i ;
465 xlen = Xsetp [1] ;
466 X1x = X->x ;
467 X2x = X2->x ;
468 Lnz = L->nz ;
469
470 if (xtype == CHOLMOD_REAL)
471 {
472 fl = 2 * xlen ;
473 for (k = 0 ; k < xlen ; k++)
474 {
475 j = Xseti [k] ;
476 fl += 4 * Lnz [j] ;
477 err = X1x [j] - X2x [j] ;
478 err = ABS (err) ;
479 resid [3] = MAX (resid [3], err) ;
480 }
481 }
482 else /* (xtype == CHOLMOD_COMPLEX) */
483 {
484 fl = 16 * xlen ;
485 for (k = 0 ; k < xlen ; k++)
486 {
487 j = Xseti [k] ;
488 fl += 16 * Lnz [j] ;
489 err = X1x [2*j ] - X2x [2*j ] ;
490 err = ABS (err) ;
491 resid [3] = MAX (resid [3], err) ;
492 err = X1x [2*j+1] - X2x [2*j+1] ;
493 err = ABS (err) ;
494 resid [3] = MAX (resid [3], err) ;
495 }
496 }
497
498 if (timelog) fprintf (timelog, "%g %g %g %g\n",
499 (double) i, (double) xlen, fl, t);
500
501 /* clear B for the next test */
502 if (xtype == CHOLMOD_REAL)
503 {
504 B2x [i] = 0 ;
505 }
506 else /* (xtype == CHOLMOD_COMPLEX) */
507 {
508 B2x [2*i ] = 0 ;
509 B2x [2*i+1] = 0 ;
510 }
511 }
512
513 if (timelog)
514 {
515 fprintf (timelog, "] ; resid = %g ;\n", resid [3]) ;
516 fprintf (timelog, "lnz = %g ;\n", cm->lnz) ;
517 fprintf (timelog, "t = %g ; %% dense solve time\n", ts [2]) ;
518 fclose (timelog) ;
519 }
520
521 #ifndef NMATRIXOPS
522 resid [3] = resid [3] / cholmod_l_norm_dense (X, 1, cm) ;
523 #endif
524
525 cholmod_l_free_dense (&Ywork, cm) ;
526 cholmod_l_free_dense (&Ework, cm) ;
527 cholmod_l_free_dense (&X2, cm) ;
528 cholmod_l_free_dense (&B2, cm) ;
529 cholmod_l_free_sparse (&Xset, cm) ;
530 cholmod_l_free_sparse (&Bset, cm) ;
531 }
532
533 /* ------------------------------------------------------------------ */
534 /* compute the residual */
535 /* ------------------------------------------------------------------ */
536
537 if (method < 3)
538 {
539 #ifndef NMATRIXOPS
540 if (A->stype == 0)
541 {
542 /* (AA'+beta*I)x=b is the linear system that was solved */
543 /* W = A'*X */
544 W = cholmod_l_allocate_dense (A->ncol, 1, A->ncol, xtype, cm) ;
545 cholmod_l_sdmult (A, 2, one, zero, X, W, cm) ;
546 /* R = B - beta*X */
547 cholmod_l_free_dense (&R, cm) ;
548 R = cholmod_l_zeros (n, 1, xtype, cm) ;
549 Rx = R->x ;
550 Rz = R->z ;
551 Xx = X->x ;
552 Xz = X->z ;
553 if (xtype == CHOLMOD_REAL)
554 {
555 for (i = 0 ; i < n ; i++)
556 {
557 Rx [i] = Bx [i] - beta [0] * Xx [i] ;
558 }
559 }
560 else if (xtype == CHOLMOD_COMPLEX)
561 {
562 /* complex case */
563 for (i = 0 ; i < n ; i++)
564 {
565 Rx [2*i ] = Bx [2*i ] - beta [0] * Xx [2*i ] ;
566 Rx [2*i+1] = Bx [2*i+1] - beta [1] * Xx [2*i+1] ;
567 }
568 }
569 else /* (xtype == CHOLMOD_ZOMPLEX) */
570 {
571 /* zomplex case */
572 for (i = 0 ; i < n ; i++)
573 {
574 Rx [i] = Bx [i] - beta [0] * Xx [i] ;
575 Rz [i] = Bz [i] - beta [1] * Xz [i] ;
576 }
577 }
578
579 /* R = A*W - R */
580 cholmod_l_sdmult (A, 0, one, minusone, W, R, cm) ;
581 cholmod_l_free_dense (&W, cm) ;
582 }
583 else
584 {
585 /* Ax=b was factorized and solved, R = B-A*X */
586 cholmod_l_free_dense (&R, cm) ;
587 R = cholmod_l_copy_dense (B, cm) ;
588 cholmod_l_sdmult (A, 0, minusone, one, X, R, cm) ;
589 }
590 rnorm = cholmod_l_norm_dense (R, 0, cm) ; /* max abs. entry */
591 xnorm = cholmod_l_norm_dense (X, 0, cm) ; /* max abs. entry */
592
593 axbnorm = (anorm * xnorm + bnorm + ((n == 0) ? 1 : 0)) ;
594 resid [method] = rnorm / axbnorm ;
595 #else
596 printf ("residual not computed (requires CHOLMOD/MatrixOps)\n") ;
597 #endif
598 }
599 }
600
601 tot = ta + tf + ts [0] ;
602
603 /* ---------------------------------------------------------------------- */
604 /* iterative refinement (real symmetric case only) */
605 /* ---------------------------------------------------------------------- */
606
607 resid2 = -1 ;
608 #ifndef NMATRIXOPS
609 if (A->stype != 0 && A->xtype == CHOLMOD_REAL)
610 {
611 cholmod_dense *R2 ;
612
613 /* R2 = A\(B-A*X) */
614 R2 = cholmod_l_solve (CHOLMOD_A, L, R, cm) ;
615 /* compute X = X + A\(B-A*X) */
616 Xx = X->x ;
617 Rx = R2->x ;
618 for (i = 0 ; i < n ; i++)
619 {
620 Xx [i] = Xx [i] + Rx [i] ;
621 }
622 cholmod_l_free_dense (&R2, cm) ;
623 cholmod_l_free_dense (&R, cm) ;
624
625 /* compute the new residual, R = B-A*X */
626 cholmod_l_free_dense (&R, cm) ;
627 R = cholmod_l_copy_dense (B, cm) ;
628 cholmod_l_sdmult (A, 0, minusone, one, X, R, cm) ;
629 rnorm2 = cholmod_l_norm_dense (R, 0, cm) ;
630 resid2 = rnorm2 / axbnorm ;
631 }
632 #endif
633
634 cholmod_l_free_dense (&R, cm) ;
635
636 /* ---------------------------------------------------------------------- */
637 /* print results */
638 /* ---------------------------------------------------------------------- */
639
640 anz = cm->anz ;
641 for (i = 0 ; i < CHOLMOD_MAXMETHODS ; i++)
642 /* for (i = 4 ; i < 3 ; i++) */
643 {
644 fl = cm->method [i].fl ;
645 xlnz = cm->method [i].lnz ;
646 cm->method [i].fl = -1 ;
647 cm->method [i].lnz = -1 ;
648 ordering = cm->method [i].ordering ;
649 if (fl >= 0)
650 {
651 printf ("Ordering: ") ;
652 if (ordering == CHOLMOD_POSTORDERED) printf ("postordered ") ;
653 if (ordering == CHOLMOD_NATURAL) printf ("natural ") ;
654 if (ordering == CHOLMOD_GIVEN) printf ("user ") ;
655 if (ordering == CHOLMOD_AMD) printf ("AMD ") ;
656 if (ordering == CHOLMOD_METIS) printf ("METIS ") ;
657 if (ordering == CHOLMOD_NESDIS) printf ("NESDIS ") ;
658 if (xlnz > 0)
659 {
660 printf ("fl/lnz %10.1f", fl / xlnz) ;
661 }
662 if (anz > 0)
663 {
664 printf (" lnz/anz %10.1f", xlnz / anz) ;
665 }
666 printf ("\n") ;
667 }
668 }
669
670 printf ("ints in L: %15.0f, doubles in L: %15.0f\n",
671 (double) isize, (double) xsize) ;
672 printf ("factor flops %g nnz(L) %15.0f (w/no amalgamation)\n",
673 cm->fl, cm->lnz) ;
674 if (A->stype == 0)
675 {
676 printf ("nnz(A): %15.0f\n", cm->anz) ;
677 }
678 else
679 {
680 printf ("nnz(A*A'): %15.0f\n", cm->anz) ;
681 }
682 if (cm->lnz > 0)
683 {
684 printf ("flops / nnz(L): %8.1f\n", cm->fl / cm->lnz) ;
685 }
686 if (anz > 0)
687 {
688 printf ("nnz(L) / nnz(A): %8.1f\n", cm->lnz / cm->anz) ;
689 }
690 printf ("analyze cputime: %12.4f\n", ta) ;
691 printf ("factor cputime: %12.4f mflop: %8.1f\n", tf,
692 (tf == 0) ? 0 : (1e-6*cm->fl / tf)) ;
693 printf ("solve cputime: %12.4f mflop: %8.1f\n", ts [0],
694 (ts [0] == 0) ? 0 : (1e-6*4*cm->lnz / ts [0])) ;
695 printf ("overall cputime: %12.4f mflop: %8.1f\n",
696 tot, (tot == 0) ? 0 : (1e-6 * (cm->fl + 4 * cm->lnz) / tot)) ;
697 printf ("solve cputime: %12.4f mflop: %8.1f (%d trials)\n", ts [1],
698 (ts [1] == 0) ? 0 : (1e-6*4*cm->lnz / ts [1]), NTRIALS) ;
699 printf ("solve2 cputime: %12.4f mflop: %8.1f (%d trials)\n", ts [2],
700 (ts [2] == 0) ? 0 : (1e-6*4*cm->lnz / ts [2]), NTRIALS) ;
701 printf ("peak memory usage: %12.0f (MB)\n",
702 (double) (cm->memory_usage) / 1048576.) ;
703 printf ("residual (|Ax-b|/(|A||x|+|b|)): ") ;
704 for (method = 0 ; method <= nmethods ; method++)
705 {
706 printf ("%8.2e ", resid [method]) ;
707 }
708 printf ("\n") ;
709 if (resid2 >= 0)
710 {
711 printf ("residual %8.1e (|Ax-b|/(|A||x|+|b|))"
712 " after iterative refinement\n", resid2) ;
713 }
714 printf ("rcond %8.1e\n\n", rcond) ;
715
716 if (L_is_super)
717 {
718 cholmod_l_gpu_stats (cm) ;
719 }
720
721 cholmod_l_free_factor (&L, cm) ;
722 cholmod_l_free_dense (&X, cm) ;
723
724 /* ---------------------------------------------------------------------- */
725 /* free matrices and finish CHOLMOD */
726 /* ---------------------------------------------------------------------- */
727
728 cholmod_l_free_sparse (&A, cm) ;
729 cholmod_l_free_dense (&B, cm) ;
730 cholmod_l_finish (cm) ;
731
732 return (0) ;
733 }
734