1 /* ========================================================================== */
2 /* === Demo/cholmod_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_demo matrixfile
19 * cholmod_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 * cholmod_demo is the same as cholmod_l_demo, except for the size of the
31 * basic integer (int vs SuiteSparse_long)
32 */
33
34 #include "cholmod_demo.h"
35 #define NTRIALS 100
36
37 /* ff is a global variable so that it can be closed by my_handler */
38 FILE *ff ;
39
40 /* halt if an error occurs */
my_handler(int status,const char * file,int line,const char * message)41 static void my_handler (int status, const char *file, int line,
42 const char *message)
43 {
44 printf ("cholmod error: file: %s line: %d status: %d: %s\n",
45 file, line, status, message) ;
46 if (status < 0)
47 {
48 if (ff != NULL) fclose (ff) ;
49 exit (0) ;
50 }
51 }
52
main(int argc,char ** argv)53 int main (int argc, char **argv)
54 {
55 double resid [4], t, ta, tf, ts [3], tot, bnorm, xnorm, anorm, rnorm, fl,
56 anz, axbnorm, rnorm2, resid2, rcond ;
57 FILE *f ;
58 cholmod_sparse *A ;
59 cholmod_dense *X = NULL, *B, *W, *R = NULL ;
60 double one [2], zero [2], minusone [2], beta [2], xlnz ;
61 cholmod_common Common, *cm ;
62 cholmod_factor *L ;
63 double *Bx, *Rx, *Xx ;
64 int i, n, isize, xsize, ordering, xtype, s, ss, lnz ;
65 int trial, method, L_is_super ;
66 int ver [3] ;
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 if (argc > 1)
78 {
79 if ((f = fopen (argv [1], "r")) == NULL)
80 {
81 my_handler (CHOLMOD_INVALID, __FILE__, __LINE__,
82 "unable to open file") ;
83 }
84 ff = f ;
85 }
86 else
87 {
88 f = stdin ;
89 }
90
91 /* ---------------------------------------------------------------------- */
92 /* start CHOLMOD and set parameters */
93 /* ---------------------------------------------------------------------- */
94
95 cm = &Common ;
96 cholmod_start (cm) ;
97 CHOLMOD_FUNCTION_DEFAULTS ; /* just for testing (not required) */
98
99 /* use default parameter settings, except for the error handler. This
100 * demo program terminates if an error occurs (out of memory, not positive
101 * definite, ...). It makes the demo program simpler (no need to check
102 * CHOLMOD error conditions). This non-default parameter setting has no
103 * effect on performance. */
104 cm->error_handler = my_handler ;
105
106 /* Note that CHOLMOD will do a supernodal LL' or a simplicial LDL' by
107 * default, automatically selecting the latter if flop/nnz(L) < 40. */
108
109 /* ---------------------------------------------------------------------- */
110 /* create basic scalars */
111 /* ---------------------------------------------------------------------- */
112
113 zero [0] = 0 ;
114 zero [1] = 0 ;
115 one [0] = 1 ;
116 one [1] = 0 ;
117 minusone [0] = -1 ;
118 minusone [1] = 0 ;
119 beta [0] = 1e-6 ;
120 beta [1] = 0 ;
121
122 /* ---------------------------------------------------------------------- */
123 /* read in a matrix */
124 /* ---------------------------------------------------------------------- */
125
126 printf ("\n---------------------------------- cholmod_demo:\n") ;
127 cholmod_version (ver) ;
128 printf ("cholmod version %d.%d.%d\n", ver [0], ver [1], ver [2]) ;
129 SuiteSparse_version (ver) ;
130 printf ("SuiteSparse version %d.%d.%d\n", ver [0], ver [1], ver [2]) ;
131 A = cholmod_read_sparse (f, cm) ;
132 if (ff != NULL)
133 {
134 fclose (ff) ;
135 ff = NULL ;
136 }
137 xtype = A->xtype ;
138 anorm = 1 ;
139 #ifndef NMATRIXOPS
140 anorm = cholmod_norm_sparse (A, 0, cm) ;
141 printf ("norm (A,inf) = %g\n", anorm) ;
142 printf ("norm (A,1) = %g\n", cholmod_norm_sparse (A, 1, cm)) ;
143 #endif
144 cholmod_print_sparse (A, "A", cm) ;
145
146 if (A->nrow > A->ncol)
147 {
148 /* Transpose A so that A'A+beta*I will be factorized instead */
149 cholmod_sparse *C = cholmod_transpose (A, 2, cm) ;
150 cholmod_free_sparse (&A, cm) ;
151 A = C ;
152 printf ("transposing input matrix\n") ;
153 }
154
155 /* ---------------------------------------------------------------------- */
156 /* create an arbitrary right-hand-side */
157 /* ---------------------------------------------------------------------- */
158
159 n = A->nrow ;
160 B = cholmod_zeros (n, 1, xtype, cm) ;
161 Bx = B->x ;
162
163 #if GHS
164 {
165 /* b = A*ones(n,1), used by Gould, Hu, and Scott in their experiments */
166 cholmod_dense *X0 ;
167 X0 = cholmod_ones (A->ncol, 1, xtype, cm) ;
168 cholmod_sdmult (A, 0, one, zero, X0, B, cm) ;
169 cholmod_free_dense (&X0, cm) ;
170 }
171 #else
172 if (xtype == CHOLMOD_REAL)
173 {
174 /* real case */
175 for (i = 0 ; i < n ; i++)
176 {
177 double x = n ;
178 Bx [i] = 1 + i / x ;
179 }
180 }
181 else
182 {
183 /* complex case */
184 for (i = 0 ; i < n ; i++)
185 {
186 double x = n ;
187 Bx [2*i ] = 1 + i / x ; /* real part of B(i) */
188 Bx [2*i+1] = (x/2 - i) / (3*x) ; /* imag part of B(i) */
189 }
190 }
191 #endif
192
193 cholmod_print_dense (B, "B", cm) ;
194 bnorm = 1 ;
195 #ifndef NMATRIXOPS
196 bnorm = cholmod_norm_dense (B, 0, cm) ; /* max norm */
197 printf ("bnorm %g\n", bnorm) ;
198 #endif
199
200 /* ---------------------------------------------------------------------- */
201 /* analyze and factorize */
202 /* ---------------------------------------------------------------------- */
203
204 t = CPUTIME ;
205 L = cholmod_analyze (A, cm) ;
206 ta = CPUTIME - t ;
207 ta = MAX (ta, 0) ;
208
209 printf ("Analyze: flop %g lnz %g\n", cm->fl, cm->lnz) ;
210
211 if (A->stype == 0)
212 {
213 printf ("Factorizing A*A'+beta*I\n") ;
214 t = CPUTIME ;
215 cholmod_factorize_p (A, beta, NULL, 0, L, cm) ;
216 tf = CPUTIME - t ;
217 tf = MAX (tf, 0) ;
218 }
219 else
220 {
221 printf ("Factorizing A\n") ;
222 t = CPUTIME ;
223 cholmod_factorize (A, L, cm) ;
224 tf = CPUTIME - t ;
225 tf = MAX (tf, 0) ;
226 }
227
228 cholmod_print_factor (L, "L", cm) ;
229
230 /* determine the # of integers's and reals's in L. See cholmod_free */
231 if (L->is_super)
232 {
233 s = L->nsuper + 1 ;
234 xsize = L->xsize ;
235 ss = L->ssize ;
236 isize =
237 n /* L->Perm */
238 + n /* L->ColCount, nz in each column of 'pure' L */
239 + s /* L->pi, column pointers for L->s */
240 + s /* L->px, column pointers for L->x */
241 + s /* L->super, starting column index of each supernode */
242 + ss ; /* L->s, the pattern of the supernodes */
243 }
244 else
245 {
246 /* this space can increase if you change parameters to their non-
247 * default values (cm->final_pack, for example). */
248 lnz = L->nzmax ;
249 xsize = lnz ;
250 isize =
251 n /* L->Perm */
252 + n /* L->ColCount, nz in each column of 'pure' L */
253 + n+1 /* L->p, column pointers */
254 + lnz /* L->i, integer row indices */
255 + n /* L->nz, nz in each column of L */
256 + n+2 /* L->next, link list */
257 + n+2 ; /* L->prev, link list */
258 }
259
260 /* solve with Bset will change L from simplicial to supernodal */
261 rcond = cholmod_rcond (L, cm) ;
262 L_is_super = L->is_super ;
263
264 /* ---------------------------------------------------------------------- */
265 /* solve */
266 /* ---------------------------------------------------------------------- */
267
268 for (method = 0 ; method <= 3 ; method++)
269 {
270 double x = n ;
271 resid [method] = -1 ; /* not yet computed */
272
273 if (method == 0)
274 {
275 /* basic solve, just once */
276 t = CPUTIME ;
277 X = cholmod_solve (CHOLMOD_A, L, B, cm) ;
278 ts [0] = CPUTIME - t ;
279 ts [0] = MAX (ts [0], 0) ;
280 }
281 else if (method == 1)
282 {
283 /* basic solve, many times, but keep the last one */
284 t = CPUTIME ;
285 for (trial = 0 ; trial < NTRIALS ; trial++)
286 {
287 cholmod_free_dense (&X, cm) ;
288 Bx [0] = 1 + trial / x ; /* tweak B each iteration */
289 X = cholmod_solve (CHOLMOD_A, L, B, cm) ;
290 }
291 ts [1] = CPUTIME - t ;
292 ts [1] = MAX (ts [1], 0) / NTRIALS ;
293 }
294 else if (method == 2)
295 {
296 /* solve with reused workspace */
297 cholmod_dense *Ywork = NULL, *Ework = NULL ;
298 cholmod_free_dense (&X, cm) ;
299
300 t = CPUTIME ;
301 for (trial = 0 ; trial < NTRIALS ; trial++)
302 {
303 Bx [0] = 1 + trial / x ; /* tweak B each iteration */
304 cholmod_solve2 (CHOLMOD_A, L, B, NULL, &X, NULL,
305 &Ywork, &Ework, cm) ;
306 }
307 cholmod_free_dense (&Ywork, cm) ;
308 cholmod_free_dense (&Ework, cm) ;
309 ts [2] = CPUTIME - t ;
310 ts [2] = MAX (ts [2], 0) / NTRIALS ;
311 }
312 else
313 {
314 /* solve with reused workspace and sparse Bset */
315 cholmod_dense *Ywork = NULL, *Ework = NULL ;
316 cholmod_dense *X2 = NULL, *B2 = NULL ;
317 cholmod_sparse *Bset, *Xset = NULL ;
318 int *Bsetp, *Bseti, *Xsetp, *Xseti, xlen, j, k, *Lnz ;
319 double *X1x, *X2x, *B2x, err ;
320 FILE *timelog = fopen ("timelog.m", "w") ;
321 if (timelog) fprintf (timelog, "results = [\n") ;
322
323 B2 = cholmod_zeros (n, 1, xtype, cm) ;
324 B2x = B2->x ;
325
326 Bset = cholmod_allocate_sparse (n, 1, 1, FALSE, TRUE, 0,
327 CHOLMOD_PATTERN, cm) ;
328 Bsetp = Bset->p ;
329 Bseti = Bset->i ;
330 Bsetp [0] = 0 ; /* nnz(B) is 1 (it can be anything) */
331 Bsetp [1] = 1 ;
332 resid [3] = 0 ;
333
334 for (i = 0 ; i < MIN (100,n) ; i++)
335 {
336 /* B (i) is nonzero, all other entries are ignored
337 (implied to be zero) */
338 Bseti [0] = i ;
339 if (xtype == CHOLMOD_REAL)
340 {
341 B2x [i] = 3.1 * i + 0.9 ;
342 }
343 else
344 {
345 B2x [2*i ] = i + 0.042 ;
346 B2x [2*i+1] = i - 92.7 ;
347 }
348
349 /* first get the entire solution, to compare against */
350 cholmod_solve2 (CHOLMOD_A, L, B2, NULL, &X, NULL,
351 &Ywork, &Ework, cm) ;
352
353 /* now get the sparse solutions; this will change L from
354 supernodal to simplicial */
355
356 if (i == 0)
357 {
358 /* first solve can be slower because it has to allocate
359 space for X2, Xset, etc, and change L.
360 So don't time it */
361 cholmod_solve2 (CHOLMOD_A, L, B2, Bset, &X2, &Xset,
362 &Ywork, &Ework, cm) ;
363 }
364
365 t = CPUTIME ;
366 for (trial = 0 ; trial < NTRIALS ; trial++)
367 {
368 /* solve Ax=b but only to get x(i).
369 b is all zero except for b(i).
370 This takes O(xlen) time */
371 cholmod_solve2 (CHOLMOD_A, L, B2, Bset, &X2, &Xset,
372 &Ywork, &Ework, cm) ;
373 }
374 t = CPUTIME - t ;
375 t = MAX (t, 0) / NTRIALS ;
376
377 /* check the solution and log the time */
378 Xsetp = Xset->p ;
379 Xseti = Xset->i ;
380 xlen = Xsetp [1] ;
381 X1x = X->x ;
382 X2x = X2->x ;
383 Lnz = L->nz ;
384
385 /*
386 printf ("\ni %d xlen %d (%p %p)\n", i, xlen, X1x, X2x) ;
387 */
388
389 if (xtype == CHOLMOD_REAL)
390 {
391 fl = 2 * xlen ;
392 for (k = 0 ; k < xlen ; k++)
393 {
394 j = Xseti [k] ;
395 fl += 4 * Lnz [j] ;
396 err = X1x [j] - X2x [j] ;
397 err = ABS (err) ;
398 resid [3] = MAX (resid [3], err) ;
399 }
400 }
401 else
402 {
403 fl = 16 * xlen ;
404 for (k = 0 ; k < xlen ; k++)
405 {
406 j = Xseti [k] ;
407 fl += 16 * Lnz [j] ;
408 err = X1x [2*j ] - X2x [2*j ] ;
409 err = ABS (err) ;
410 resid [3] = MAX (resid [3], err) ;
411 err = X1x [2*j+1] - X2x [2*j+1] ;
412 err = ABS (err) ;
413 resid [3] = MAX (resid [3], err) ;
414 }
415 }
416 if (timelog) fprintf (timelog, "%g %g %g %g\n",
417 (double) i, (double) xlen, fl, t);
418
419 /* clear B for the next test */
420 if (xtype == CHOLMOD_REAL)
421 {
422 B2x [i] = 0 ;
423 }
424 else
425 {
426 B2x [2*i ] = 0 ;
427 B2x [2*i+1] = 0 ;
428 }
429
430 }
431
432 if (timelog)
433 {
434 fprintf (timelog, "] ; resid = %g ;\n", resid [3]) ;
435 fprintf (timelog, "lnz = %g ;\n", cm->lnz) ;
436 fprintf (timelog, "t = %g ; %% dense solve time\n", ts [2]) ;
437 fclose (timelog) ;
438 }
439
440 #ifndef NMATRIXOPS
441 resid [3] = resid [3] / cholmod_norm_dense (X, 1, cm) ;
442 #endif
443
444 cholmod_free_dense (&Ywork, cm) ;
445 cholmod_free_dense (&Ework, cm) ;
446 cholmod_free_dense (&X2, cm) ;
447 cholmod_free_dense (&B2, cm) ;
448 cholmod_free_sparse (&Xset, cm) ;
449 cholmod_free_sparse (&Bset, cm) ;
450 }
451
452 /* ------------------------------------------------------------------ */
453 /* compute the residual */
454 /* ------------------------------------------------------------------ */
455
456 if (method < 3)
457 {
458 #ifndef NMATRIXOPS
459
460 if (A->stype == 0)
461 {
462 /* (AA'+beta*I)x=b is the linear system that was solved */
463 /* W = A'*X */
464 W = cholmod_allocate_dense (A->ncol, 1, A->ncol, xtype, cm) ;
465 cholmod_sdmult (A, 2, one, zero, X, W, cm) ;
466 /* R = B - beta*X */
467 cholmod_free_dense (&R, cm) ;
468 R = cholmod_zeros (n, 1, xtype, cm) ;
469 Rx = R->x ;
470 Xx = X->x ;
471 if (xtype == CHOLMOD_REAL)
472 {
473 for (i = 0 ; i < n ; i++)
474 {
475 Rx [i] = Bx [i] - beta [0] * Xx [i] ;
476 }
477 }
478 else
479 {
480 /* complex case */
481 for (i = 0 ; i < n ; i++)
482 {
483 Rx [2*i ] = Bx [2*i ] - beta [0] * Xx [2*i ] ;
484 Rx [2*i+1] = Bx [2*i+1] - beta [0] * Xx [2*i+1] ;
485 }
486 }
487 /* R = A*W - R */
488 cholmod_sdmult (A, 0, one, minusone, W, R, cm) ;
489 cholmod_free_dense (&W, cm) ;
490 }
491 else
492 {
493 /* Ax=b was factorized and solved, R = B-A*X */
494 cholmod_free_dense (&R, cm) ;
495 R = cholmod_copy_dense (B, cm) ;
496 cholmod_sdmult (A, 0, minusone, one, X, R, cm) ;
497 }
498 rnorm = -1 ;
499 xnorm = 1 ;
500 rnorm = cholmod_norm_dense (R, 0, cm) ; /* max abs. entry */
501 xnorm = cholmod_norm_dense (X, 0, cm) ; /* max abs. entry */
502 axbnorm = (anorm * xnorm + bnorm + ((n == 0) ? 1 : 0)) ;
503 resid [method] = rnorm / axbnorm ;
504 #else
505 printf ("residual not computed (requires CHOLMOD/MatrixOps)\n") ;
506 #endif
507 }
508 }
509
510 tot = ta + tf + ts [0] ;
511
512 /* ---------------------------------------------------------------------- */
513 /* iterative refinement (real symmetric case only) */
514 /* ---------------------------------------------------------------------- */
515
516 resid2 = -1 ;
517 #ifndef NMATRIXOPS
518 if (A->stype != 0 && A->xtype == CHOLMOD_REAL)
519 {
520 cholmod_dense *R2 ;
521
522 /* R2 = A\(B-A*X) */
523 R2 = cholmod_solve (CHOLMOD_A, L, R, cm) ;
524 /* compute X = X + A\(B-A*X) */
525 Xx = X->x ;
526 Rx = R2->x ;
527 for (i = 0 ; i < n ; i++)
528 {
529 Xx [i] = Xx [i] + Rx [i] ;
530 }
531 cholmod_free_dense (&R2, cm) ;
532 cholmod_free_dense (&R, cm) ;
533
534 /* compute the new residual, R = B-A*X */
535 cholmod_free_dense (&R, cm) ;
536 R = cholmod_copy_dense (B, cm) ;
537 cholmod_sdmult (A, 0, minusone, one, X, R, cm) ;
538 rnorm2 = cholmod_norm_dense (R, 0, cm) ;
539 resid2 = rnorm2 / axbnorm ;
540 }
541 #endif
542
543 cholmod_free_dense (&R, cm) ;
544
545 /* ---------------------------------------------------------------------- */
546 /* print results */
547 /* ---------------------------------------------------------------------- */
548
549 anz = cm->anz ;
550 for (i = 0 ; i < CHOLMOD_MAXMETHODS ; i++)
551 {
552 fl = cm->method [i].fl ;
553 xlnz = cm->method [i].lnz ;
554 cm->method [i].fl = -1 ;
555 cm->method [i].lnz = -1 ;
556 ordering = cm->method [i].ordering ;
557 if (fl >= 0)
558 {
559 printf ("Ordering: ") ;
560 if (ordering == CHOLMOD_POSTORDERED) printf ("postordered ") ;
561 if (ordering == CHOLMOD_NATURAL) printf ("natural ") ;
562 if (ordering == CHOLMOD_GIVEN) printf ("user ") ;
563 if (ordering == CHOLMOD_AMD) printf ("AMD ") ;
564 if (ordering == CHOLMOD_METIS) printf ("METIS ") ;
565 if (ordering == CHOLMOD_NESDIS) printf ("NESDIS ") ;
566 if (xlnz > 0)
567 {
568 printf ("fl/lnz %10.1f", fl / xlnz) ;
569 }
570 if (anz > 0)
571 {
572 printf (" lnz/anz %10.1f", xlnz / anz) ;
573 }
574 printf ("\n") ;
575 }
576 }
577
578 printf ("ints in L: %15.0f, doubles in L: %15.0f\n",
579 (double) isize, (double) xsize) ;
580 printf ("factor flops %g nnz(L) %15.0f (w/no amalgamation)\n",
581 cm->fl, cm->lnz) ;
582 if (A->stype == 0)
583 {
584 printf ("nnz(A): %15.0f\n", cm->anz) ;
585 }
586 else
587 {
588 printf ("nnz(A*A'): %15.0f\n", cm->anz) ;
589 }
590 if (cm->lnz > 0)
591 {
592 printf ("flops / nnz(L): %8.1f\n", cm->fl / cm->lnz) ;
593 }
594 if (anz > 0)
595 {
596 printf ("nnz(L) / nnz(A): %8.1f\n", cm->lnz / cm->anz) ;
597 }
598 printf ("analyze cputime: %12.4f\n", ta) ;
599 printf ("factor cputime: %12.4f mflop: %8.1f\n", tf,
600 (tf == 0) ? 0 : (1e-6*cm->fl / tf)) ;
601 printf ("solve cputime: %12.4f mflop: %8.1f\n", ts [0],
602 (ts [0] == 0) ? 0 : (1e-6*4*cm->lnz / ts [0])) ;
603 printf ("overall cputime: %12.4f mflop: %8.1f\n",
604 tot, (tot == 0) ? 0 : (1e-6 * (cm->fl + 4 * cm->lnz) / tot)) ;
605 printf ("solve cputime: %12.4f mflop: %8.1f (%d trials)\n", ts [1],
606 (ts [1] == 0) ? 0 : (1e-6*4*cm->lnz / ts [1]), NTRIALS) ;
607 printf ("solve2 cputime: %12.4f mflop: %8.1f (%d trials)\n", ts [2],
608 (ts [2] == 0) ? 0 : (1e-6*4*cm->lnz / ts [2]), NTRIALS) ;
609 printf ("peak memory usage: %12.0f (MB)\n",
610 (double) (cm->memory_usage) / 1048576.) ;
611 printf ("residual (|Ax-b|/(|A||x|+|b|)): ") ;
612 for (method = 0 ; method <= 3 ; method++)
613 {
614 printf ("%8.2e ", resid [method]) ;
615 }
616 printf ("\n") ;
617 if (resid2 >= 0)
618 {
619 printf ("residual %8.1e (|Ax-b|/(|A||x|+|b|))"
620 " after iterative refinement\n", resid2) ;
621 }
622
623 printf ("rcond %8.1e\n\n", rcond) ;
624
625 if (L_is_super)
626 {
627 cholmod_gpu_stats (cm) ;
628 }
629
630 cholmod_free_factor (&L, cm) ;
631 cholmod_free_dense (&X, cm) ;
632
633 /* ---------------------------------------------------------------------- */
634 /* free matrices and finish CHOLMOD */
635 /* ---------------------------------------------------------------------- */
636
637 cholmod_free_sparse (&A, cm) ;
638 cholmod_free_dense (&B, cm) ;
639 cholmod_finish (cm) ;
640
641 return (0) ;
642 }
643