1 /* ========================================================================== */
2 /* === Tcov/cm ============================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Tcov Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* A program for exhaustive statement-coverage for CHOLMOD, AMD, COLAMD, and
11 * CCOLAMD. It tests every line of code in all three packages.
12 *
13 * For a complete test, all CHOLMOD modules, AMD, COLAMD, CCOLAMD, METIS,
14 * LAPACK, and the BLAS are required. A partial test can be performed without
15 * the Supernodal and/or Partition modules. METIS is not required if the
16 * Partition module is not installed. LAPACK and the BLAS are not required
17 * if the Supernodal module is not installed.
18 *
19 * Usage:
20 *
21 * cm < input > output
22 *
23 * where "input" contains a sparse matrix in triplet form. The first line of
24 * the file contains four or five integers:
25 *
26 * nrow ncol nnz stype complex
27 *
28 * where the matrix is nrow-by-ncol. nnz is the number of (i,j,aij) triplets
29 * in the rest of the file, one triplet per line. stype is -1 (symmetric
30 * with entries in lower triangular part provided), 0 (unsymmetric), or 1
31 * (symmetric upper). If the 5th entry is missing, or 0, then the matrix is
32 * real; if 1 the matrix is complex, if 2 the matrix is zomplex. Each
33 * subsequent line contains
34 *
35 * i j aij
36 *
37 * for the row index, column index, and value of A(i,j). Duplicate entries
38 * are summed. If stype is 2 or 3, the rest of the file is ignored, and a
39 * special test matrix is constructed (2: arrowhead, 3: tridiagonal plus a
40 * dense row). Test matrices are located in the Matrix/ subdirectory.
41 *
42 * For complex matrices, each line consists of
43 *
44 * i j xij zij
45 *
46 * where xij is the real part of A(i,j) and zij is the imaginary part.
47 *
48 * cm takes one optional parameter. If present (it does not matter what the
49 * argument is, actually) then extension memory-failure tests are performed.
50 */
51
52 #include "cm.h"
53
54 /* ========================================================================== */
55 /* === global variables ===================================================== */
56 /* ========================================================================== */
57
58 double zero [2], one [2], minusone [2] ;
59 cholmod_common Common, *cm ;
60 cholmod_dense *M1 ;
61 Int dot = 0 ;
62 double Zero [2] ;
63
64
65 /* ========================================================================== */
66 /* === my_rand ============================================================== */
67 /* ========================================================================== */
68
69 /* The POSIX example of rand, duplicated here so that the same sequence will
70 * be generated on different machines. */
71
72 static unsigned long next = 1 ;
73
74 /* RAND_MAX assumed to be 32767 */
my_rand(void)75 Int my_rand (void)
76 {
77 next = next * 1103515245 + 12345 ;
78 return ((unsigned)(next/65536) % /* 32768 */ (MY_RAND_MAX + 1)) ;
79 }
80
my_srand(unsigned seed)81 void my_srand (unsigned seed)
82 {
83 next = seed ;
84 }
85
my_seed(void)86 unsigned long my_seed (void)
87 {
88 return (next) ;
89 }
90
91
92 /* ========================================================================== */
93 /* === progress ============================================================= */
94 /* ========================================================================== */
95
96 /* print a "." on stderr to indicate progress */
97
98 #define LINEWIDTH 70
99 #define COUNT 100
100
progress(Int force,char s)101 void progress (Int force, char s)
102 {
103 dot++ ;
104 if (force)
105 {
106 dot += (COUNT - (dot % COUNT)) ;
107 }
108 if (dot % COUNT == 0)
109 {
110 fprintf (stderr, "%c", s) ;
111 }
112 if (dot % (COUNT*LINEWIDTH) == 0)
113 {
114 fprintf (stderr, "\n") ;
115 }
116 fflush (stdout) ;
117 fflush (stderr) ;
118 }
119
120
121 /* ========================================================================== */
122 /* === my_handler =========================================================== */
123 /* ========================================================================== */
124
125 /* An error occurred that should not have occurred */
126
my_handler(int status,const char * file,int line,const char * msg)127 void my_handler (int status, const char *file, int line, const char *msg)
128 {
129 printf ("Error handler: file %s line %d status %d: %s\n",
130 file, line, status, msg) ;
131 if (status < CHOLMOD_OK || status > CHOLMOD_DSMALL)
132 {
133 fprintf (stderr, "\n\n************************************************"
134 "********************************\n"
135 "*** Test failure: file: %s line: %d\n"
136 "*** status: %d message: %s\n"
137 "***********************************************************"
138 "*********************\n", file, line, status, msg);
139 fflush (stderr) ;
140 fflush (stdout) ;
141 abort ( ) ;
142 }
143 }
144
145
146 /* ========================================================================== */
147 /* === Assert =============================================================== */
148 /* ========================================================================== */
149
Assert(int truth,char * file,int line)150 void Assert (int truth, char *file, int line)
151 {
152 if (!truth)
153 {
154 my_handler (-1, file, line, "") ;
155 }
156 }
157
158
159 /* ========================================================================== */
160 /* === nrand ================================================================ */
161 /* ========================================================================== */
162
163 /* return a random Int between 0 and n-1 */
164
nrand(Int n)165 Int nrand (Int n)
166 {
167 return ((n <= 0) ? (0) : (rand ( ) % n)) ;
168 }
169
170 /* ========================================================================== */
171 /* === xrand ================================================================ */
172 /* ========================================================================== */
173
174 /* return a random double between 0 and x */
175
xrand(double range)176 double xrand (double range)
177 {
178 return ((range * (double) (my_rand ( ))) / MY_RAND_MAX) ;
179 }
180
181
182 /* ========================================================================== */
183 /* === prand ================================================================ */
184 /* ========================================================================== */
185
186 /* allocate and construct a random permutation of 0:n-1 */
187
prand(Int n)188 Int *prand (Int n)
189 {
190 Int *P ;
191 Int t, j, k ;
192 P = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
193 if (P == NULL)
194 {
195 ERROR (CHOLMOD_INVALID, "cannot create random perm") ;
196 return (NULL) ;
197 }
198 for (k = 0 ; k < n ; k++)
199 {
200 P [k] = k ;
201 }
202 for (k = 0 ; k < n-1 ; k++)
203 {
204 j = k + nrand (n-k) ;
205 t = P [j] ;
206 P [j] = P [k] ;
207 P [k] = t ;
208 }
209 CHOLMOD(print_perm) (P, n, n, "Prandom", cm) ;
210 return (P) ;
211 }
212
213
214 /* ========================================================================== */
215 /* === rand_set ============================================================= */
216 /* ========================================================================== */
217
218 /* allocate and construct a random set of 0:n-1, possibly with duplicates */
219
rand_set(Int len,Int n)220 Int *rand_set (Int len, Int n)
221 {
222 Int *cset ;
223 Int k ;
224 cset = CHOLMOD(malloc) (len, sizeof (Int), cm) ;
225 if (cset == NULL)
226 {
227 ERROR (CHOLMOD_INVALID, "cannot create random set") ;
228 return (NULL) ;
229 }
230 for (k = 0 ; k < len ; k++)
231 {
232 cset [k] = nrand (n) ;
233 }
234 CHOLMOD(print_subset) (cset, len, n, "random cset", cm) ;
235 return (cset) ;
236 }
237
238
239 /* ========================================================================== */
240 /* === read_triplet ========================================================= */
241 /* ========================================================================== */
242
243 /* Read a triplet matrix from a file. */
244
245 #define MAXLINE 1024
246
read_triplet(FILE * f)247 cholmod_triplet *read_triplet
248 (
249 FILE *f
250 )
251 {
252 cholmod_triplet *T ;
253 double *Tx, *Tz ;
254 long long x1, x2, x3, x4, x5 ;
255 Int *Ti, *Tj ;
256 Int n, j, k, nrow, ncol, nz, stype, arrowhead, tridiag_plus_denserow,
257 xtype, is_complex ;
258 char s [MAXLINE] ;
259
260 /* ---------------------------------------------------------------------- */
261 /* read in a triplet matrix from a file */
262 /* ---------------------------------------------------------------------- */
263
264 dot = 0 ;
265 xtype = 0 ;
266 if (fgets (s, MAXLINE, f) == NULL)
267 {
268 return (NULL) ;
269 }
270
271 x1 = 0 ;
272 x2 = 0 ;
273 x3 = 0 ;
274 x4 = 0 ;
275 x5 = 0 ;
276 k = sscanf (s, "%lld %lld %lld %lld %lld\n", &x1, &x2, &x3, &x4, &x5) ;
277 nrow = x1 ;
278 ncol = x2 ;
279 nz = x3 ;
280 stype = x4 ;
281 xtype = x5 ;
282
283 xtype++ ;
284 is_complex = (xtype != CHOLMOD_REAL) ;
285
286 printf ("read_triplet: nrow "ID" ncol "ID" nz "ID" stype "ID" xtype "ID"\n",
287 nrow, ncol, nz, stype, xtype) ;
288
289 arrowhead = FALSE ;
290 tridiag_plus_denserow = FALSE ;
291
292 n = MAX (nrow, ncol) ;
293 if (stype == 2)
294 {
295 /* ignore nz and the rest of the file, and create an arrowhead matrix */
296 arrowhead = TRUE ;
297 nz = nrow + ncol + 1 ;
298 stype = (nrow == ncol) ? (1) : (0) ;
299 }
300 else if (stype == 3)
301 {
302 tridiag_plus_denserow = TRUE ;
303 nrow = n ;
304 ncol = n ;
305 nz = 4*n + 4 ;
306 stype = 0 ;
307 }
308
309 T = CHOLMOD(allocate_triplet) (nrow, ncol, nz, stype,
310 is_complex ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL, cm) ;
311 if (T == NULL)
312 {
313 ERROR (CHOLMOD_INVALID, "cannot create triplet matrix") ;
314 return (NULL) ;
315 }
316 Ti = T->i ;
317 Tj = T->j ;
318 Tx = T->x ;
319 Tz = T->z ;
320
321 if (arrowhead)
322 {
323 for (k = 0 ; k < MIN (nrow,ncol) ; k++)
324 {
325 Ti [k] = k ;
326 Tj [k] = k ;
327 Tx [k] = nrow + xrand (1) ; /* RAND */
328 if (is_complex)
329 {
330 Tz [k] = nrow + xrand (1) ; /* RAND */
331 }
332 }
333 for (j = 0 ; j < ncol ; j++)
334 {
335 Ti [k] = 0 ;
336 Tj [k] = j ;
337 Tx [k] = - xrand (1) ; /* RAND */
338 if (is_complex)
339 {
340 Tz [k] = - xrand (1) ; /* RAND */
341 }
342 k++ ;
343 }
344 T->nnz = k ;
345 }
346 else if (tridiag_plus_denserow)
347 {
348 /* dense row, except for the last column */
349 for (k = 0 ; k < n-1 ; k++)
350 {
351 Ti [k] = 0 ;
352 Tj [k] = k ;
353 Tx [k] = xrand (1) ; /* RAND */
354 if (is_complex)
355 {
356 Tz [k] = xrand (1) ; /* RAND */
357 }
358 }
359
360 /* diagonal */
361 for (j = 0 ; j < n ; j++)
362 {
363 Ti [k] = j ;
364 Tj [k] = j ;
365 Tx [k] = nrow + xrand (1) ; /* RAND */
366 if (is_complex)
367 {
368 Tz [k] = nrow + xrand (1) ; /* RAND */
369 }
370 k++ ;
371 }
372
373 /* superdiagonal */
374 for (j = 1 ; j < n ; j++)
375 {
376 Ti [k] = j-1 ;
377 Tj [k] = j ;
378 Tx [k] = xrand (1) ; /* RAND */
379 if (is_complex)
380 {
381 Tz [k] = xrand (1) ; /* RAND */
382 }
383 k++ ;
384 }
385
386 /* subdiagonal */
387 for (j = 0 ; j < n-1 ; j++)
388 {
389 Ti [k] = j+1 ;
390 Tj [k] = j ;
391 Tx [k] = xrand (1) ; /* RAND */
392 if (is_complex)
393 {
394 Tz [k] = xrand (1) ; /* RAND */
395 }
396 k++ ;
397 }
398
399 /* a few extra terms in the last column */
400 Ti [k] = MAX (0, n-3) ;
401 Tj [k] = n-1 ;
402 Tx [k] = xrand (1) ; /* RAND */
403 if (is_complex)
404 {
405 Tz [k] = xrand (1) ; /* RAND */
406 }
407 k++ ;
408
409 Ti [k] = MAX (0, n-4) ;
410 Tj [k] = n-1 ;
411 Tx [k] = xrand (1) ; /* RAND */
412 if (is_complex)
413 {
414 Tz [k] = xrand (1) ; /* RAND */
415 }
416 k++ ;
417
418 Ti [k] = MAX (0, n-5) ;
419 Tj [k] = n-1 ;
420 Tx [k] = xrand (1) ; /* RAND */
421 if (is_complex)
422 {
423 Tz [k] = xrand (1) ; /* RAND */
424 }
425 k++ ;
426
427 T->nnz = k ;
428 }
429 else
430 {
431 if (is_complex)
432 {
433 for (k = 0 ; k < nz ; k++)
434 {
435 if (fscanf (f,""ID" "ID" %lg %lg\n", Ti+k, Tj+k, Tx+k, Tz+k)
436 == EOF)
437 {
438 ERROR (CHOLMOD_INVALID, "Error reading triplet matrix\n") ;
439 }
440 }
441 }
442 else
443 {
444 for (k = 0 ; k < nz ; k++)
445 {
446 if (fscanf (f, ""ID" "ID" %lg\n", Ti+k, Tj+k, Tx+k) == EOF)
447 {
448 ERROR (CHOLMOD_INVALID, "Error reading triplet matrix\n") ;
449 }
450 }
451 }
452 T->nnz = nz ;
453 }
454
455 CHOLMOD(triplet_xtype) (xtype, T, cm) ;
456
457 /* ---------------------------------------------------------------------- */
458 /* print the triplet matrix */
459 /* ---------------------------------------------------------------------- */
460
461 cm->print = 4 ;
462 CHOLMOD(print_triplet) (T, "T input", cm) ;
463 cm->print = 1 ;
464 fprintf (stderr, "Test matrix: "ID"-by-"ID" with "ID" entries, stype: "ID
465 "\n",
466 (Int) T->nrow, (Int) T->ncol, (Int) T->nnz, (Int) T->stype) ;
467 printf ("\n\n======================================================\n"
468 "Test matrix: "ID"-by-"ID" with "ID" entries, stype: "ID"\n",
469 (Int) T->nrow, (Int) T->ncol, (Int) T->nnz, (Int) T->stype) ;
470
471 if (MAX (nrow, ncol) > NLARGE)
472 {
473 fprintf (stderr, "Please wait, this will take a while ...") ;
474 dot = 39*LINEWIDTH ;
475 }
476 return (T) ;
477 }
478
479
480 /* ========================================================================== */
481 /* === zeros ================================================================ */
482 /* ========================================================================== */
483
484 /* Same as cholmod_zeros or cholmod_l_zeros, except it allows for a leading
485 * dimension that is different than nrow */
486
zeros(Int nrow,Int ncol,Int d,Int xtype)487 cholmod_dense *zeros (Int nrow, Int ncol, Int d, Int xtype)
488 {
489 cholmod_dense *X ;
490 double *Xx, *Xz ;
491 Int i, nz ;
492 X = CHOLMOD(allocate_dense) (nrow, ncol, d, xtype, cm) ;
493 if (X == NULL)
494 {
495 return (NULL) ;
496 }
497 Xx = X->x ;
498 Xz = X->z ;
499 nz = MAX (1, X->nzmax) ;
500 switch (X->xtype)
501 {
502 case CHOLMOD_REAL:
503 for (i = 0 ; i < nz ; i++)
504 {
505 Xx [i] = 0 ;
506 }
507 break ;
508 case CHOLMOD_COMPLEX:
509 for (i = 0 ; i < 2*nz ; i++)
510 {
511 Xx [i] = 0 ;
512 }
513 break ;
514 case CHOLMOD_ZOMPLEX:
515 for (i = 0 ; i < nz ; i++)
516 {
517 Xx [i] = 0 ;
518 }
519 for (i = 0 ; i < nz ; i++)
520 {
521 Xz [i] = 0 ;
522 }
523 break ;
524 }
525 return (X) ;
526 }
527
528
529 /* ========================================================================== */
530 /* === xtrue ================================================================ */
531 /* ========================================================================== */
532
533 /* Allocate and construct a dense matrix, X(i,j) = i+j*d+1 */
534
xtrue(Int nrow,Int ncol,Int d,Int xtype)535 cholmod_dense *xtrue (Int nrow, Int ncol, Int d, Int xtype)
536 {
537 double *x, *z ;
538 cholmod_dense *X ;
539 Int j, i ;
540 X = zeros (nrow, ncol, d, xtype) ;
541 if (X == NULL)
542 {
543 ERROR (CHOLMOD_INVALID, "cannot create dense matrix") ;
544 return (NULL) ;
545 }
546 x = X->x ;
547 z = X->z ;
548
549 if (xtype == CHOLMOD_REAL)
550 {
551 for (j = 0 ; j < ncol ; j++)
552 {
553 for (i = 0 ; i < nrow ; i++)
554 {
555 x [i+j*d] = i+j*d + 1 ;
556 }
557 }
558 }
559 else if (xtype == CHOLMOD_COMPLEX)
560 {
561 for (j = 0 ; j < ncol ; j++)
562 {
563 for (i = 0 ; i < nrow ; i++)
564 {
565 x [2*(i+j*d) ] = i+j*d + 1 ;
566 x [2*(i+j*d)+1] = ((double) (j+i*d + 1))/10 ;
567 }
568 }
569 }
570 else
571 {
572 for (j = 0 ; j < ncol ; j++)
573 {
574 for (i = 0 ; i < nrow ; i++)
575 {
576 x [i+j*d] = i+j*d + 1 ;
577 z [i+j*d] = ((double) (j+i*d + 1))/10 ;
578 }
579 }
580 }
581 return (X) ;
582 }
583
584
585 /* ========================================================================== */
586 /* === rhs ================================================================== */
587 /* ========================================================================== */
588
589 /* Create a right-hand-side, b = A*x, where x is a known solution */
590
rhs(cholmod_sparse * A,Int nrhs,Int d)591 cholmod_dense *rhs (cholmod_sparse *A, Int nrhs, Int d)
592 {
593 Int n ;
594 cholmod_dense *W, *Z, *B ;
595
596 if (A == NULL)
597 {
598 ERROR (CHOLMOD_INVALID, "cannot compute rhs") ;
599 return (NULL) ;
600 }
601
602 n = A->nrow ;
603
604 /* B = zeros (n,rhs) but with leading dimension d */
605 B = zeros (n, nrhs, d, A->xtype) ;
606
607 /* ---------------------------------------------------------------------- */
608 /* create a known solution */
609 /* ---------------------------------------------------------------------- */
610
611 Z = xtrue (n, nrhs, d, A->xtype) ;
612
613 /* ---------------------------------------------------------------------- */
614 /* compute B = A*Z or A*A'*Z */
615 /* ---------------------------------------------------------------------- */
616
617 if (A->stype == 0)
618 {
619 /* W = A'*Z */
620 W = CHOLMOD(zeros) (A->ncol, nrhs, A->xtype, cm) ;
621 CHOLMOD(sdmult) (A, TRUE, one, zero, Z, W, cm) ;
622 /* B = A*W */
623 CHOLMOD(sdmult) (A, FALSE, one, zero, W, B, cm) ;
624 CHOLMOD(free_dense) (&W, cm) ;
625 }
626 else
627 {
628 /* B = A*Z */
629 CHOLMOD(sdmult) (A, FALSE, one, zero, Z, B, cm) ;
630 }
631 CHOLMOD(free_dense) (&Z, cm) ;
632 return (B) ;
633 }
634
635
636 /* ========================================================================== */
637 /* === resid ================================================================ */
638 /* ========================================================================== */
639
640 /* compute r = norm (A*x-b)/norm(b) or r = norm (A*A'*x-b)/norm(b) */
641
resid(cholmod_sparse * A,cholmod_dense * X,cholmod_dense * B)642 double resid (cholmod_sparse *A, cholmod_dense *X, cholmod_dense *B)
643 {
644 double r, bnorm ;
645 cholmod_dense *R, *X2, *B2 ;
646 cholmod_sparse *C, *A2 ;
647 Int d, n, nrhs, xtype ;
648
649 if (A == NULL || X == NULL || B == NULL)
650 {
651 ERROR (CHOLMOD_INVALID, "cannot compute resid") ;
652 return (1) ;
653 }
654
655 cm->status = CHOLMOD_OK ;
656 n = B->nrow ;
657
658 /* ---------------------------------------------------------------------- */
659 /* convert all inputs to an identical xtype */
660 /* ---------------------------------------------------------------------- */
661
662 xtype = MAX (A->xtype, X->xtype) ;
663 xtype = MAX (xtype, B->xtype) ;
664 A2 = NULL ;
665 B2 = NULL ;
666 X2 = NULL ;
667 if (A->xtype != xtype)
668 {
669 A2 = CHOLMOD(copy_sparse) (A, cm) ;
670 CHOLMOD(sparse_xtype) (xtype, A2, cm) ;
671 A = A2 ;
672 }
673 if (X->xtype != xtype)
674 {
675 X2 = CHOLMOD(copy_dense) (X, cm) ;
676 CHOLMOD(dense_xtype) (xtype, X2, cm) ;
677 X = X2 ;
678 }
679 if (B->xtype != xtype)
680 {
681 B2 = CHOLMOD(copy_dense) (B, cm) ;
682 CHOLMOD(dense_xtype) (xtype, B2, cm) ;
683 B = B2 ;
684 }
685
686 if (cm->status < CHOLMOD_OK)
687 {
688 ERROR (CHOLMOD_INVALID, "cannot compute resid") ;
689 CHOLMOD(free_sparse) (&A2, cm) ;
690 CHOLMOD(free_dense) (&B2, cm) ;
691 CHOLMOD(free_dense) (&X2, cm) ;
692 return (1) ;
693 }
694
695 /* ---------------------------------------------------------------------- */
696 /* get the right-hand-side, B, and its norm */
697 /* ---------------------------------------------------------------------- */
698
699 nrhs = B->ncol ;
700 d = B->d ;
701 if (nrhs == 1)
702 {
703 /* inf-norm, 1-norm, or 2-norm (random choice) */
704 bnorm = CHOLMOD(norm_dense) (B, nrand (2), cm) ;
705 }
706 else
707 {
708 /* inf-norm or 1-norm (random choice) */
709 bnorm = CHOLMOD(norm_dense) (B, nrand (1), cm) ;
710 }
711
712 /* ---------------------------------------------------------------------- */
713 /* compute the residual */
714 /* ---------------------------------------------------------------------- */
715
716 if (A->stype == 0)
717 {
718 if (n < 10 && A->xtype == CHOLMOD_REAL)
719 {
720 /* test cholmod_aat, C = A*A' */
721 C = CHOLMOD(aat) (A, NULL, 0, 1, cm) ;
722
723 /* R = B */
724 R = CHOLMOD(copy_dense) (B, cm) ;
725 /* R = C*X - R */
726 CHOLMOD(sdmult) (C, FALSE, one, minusone, X, R, cm) ;
727 CHOLMOD(free_sparse) (&C, cm) ;
728
729 }
730 else
731 {
732 /* W = A'*X */
733 cholmod_dense *W ;
734 W = CHOLMOD(zeros) (A->ncol, nrhs, A->xtype, cm) ;
735 CHOLMOD(sdmult) (A, TRUE, one, zero, X, W, cm) ;
736 /* R = B */
737 R = CHOLMOD(copy_dense) (B, cm) ;
738 /* R = A*W - R */
739 CHOLMOD(sdmult) (A, FALSE, one, minusone, W, R, cm) ;
740 CHOLMOD(free_dense) (&W, cm) ;
741 }
742 }
743 else
744 {
745 /* R = B */
746 R = CHOLMOD(copy_dense) (B, cm) ;
747 /* R = A*X - R */
748 CHOLMOD(sdmult) (A, FALSE, one, minusone, X, R, cm) ;
749 }
750
751 /* ---------------------------------------------------------------------- */
752 /* r = norm (R) / norm (B) */
753 /* ---------------------------------------------------------------------- */
754
755 r = CHOLMOD(norm_dense) (R, 1, cm) ;
756
757 CHOLMOD(free_dense) (&R, cm) ;
758 CHOLMOD(free_sparse) (&A2, cm) ;
759 CHOLMOD(free_dense) (&B2, cm) ;
760 CHOLMOD(free_dense) (&X2, cm) ;
761
762 if (bnorm > 0)
763 {
764 r /= bnorm ;
765 }
766 return (r) ;
767 }
768
769
770 /* ========================================================================== */
771 /* === resid_sparse ========================================================= */
772 /* ========================================================================== */
773
774 /* compute r = norm (A*x-b)/norm(b) or r = norm (A*A'*x-b)/norm(b) */
775
resid_sparse(cholmod_sparse * A,cholmod_sparse * X,cholmod_sparse * B)776 double resid_sparse (cholmod_sparse *A, cholmod_sparse *X, cholmod_sparse *B)
777 {
778 double r, bnorm ;
779 cholmod_sparse *R, *W, *AT, *W2 ;
780 cholmod_dense *X2, *B2 ;
781 Int n, nrhs, xtype ;
782
783 if (A == NULL || X == NULL || B == NULL)
784 {
785 ERROR (CHOLMOD_INVALID, "cannot compute resid") ;
786 return (1) ;
787 }
788
789 /* ---------------------------------------------------------------------- */
790 /* compute the residual */
791 /* ---------------------------------------------------------------------- */
792
793 xtype = MAX (A->xtype, X->xtype) ;
794 xtype = MAX (xtype, B->xtype) ;
795
796 if (xtype > CHOLMOD_REAL)
797 {
798
799 /* ------------------------------------------------------------------ */
800 /* convert X and B to dense if any is complex or zomplex */
801 /* ------------------------------------------------------------------ */
802
803 X2 = CHOLMOD(sparse_to_dense) (X, cm) ;
804 B2 = CHOLMOD(sparse_to_dense) (B, cm) ;
805 r = resid (A, X2, B2) ;
806 CHOLMOD(free_dense) (&X2, cm) ;
807 CHOLMOD(free_dense) (&B2, cm) ;
808
809 }
810 else
811 {
812
813 /* ------------------------------------------------------------------ */
814 /* all inputs are real */
815 /* ------------------------------------------------------------------ */
816
817 n = B->nrow ;
818 nrhs = B->ncol ;
819 /* inf-norm or 1-norm (random choice) */
820 bnorm = CHOLMOD(norm_sparse) (B, nrand (1), cm) ;
821
822 if (A->stype == 0)
823 {
824 /* W = A'*X */
825 AT = CHOLMOD(transpose) (A, 1, cm) ;
826 W = CHOLMOD(ssmult) (AT, X, 0, TRUE, FALSE, cm) ;
827 CHOLMOD(free_sparse) (&AT, cm) ;
828 /* W2 = A*W */
829 W2 = CHOLMOD(ssmult) (A, W, 0, TRUE, FALSE, cm) ;
830 CHOLMOD(free_sparse) (&W, cm) ;
831 /* R = W2 - B */
832 R = CHOLMOD(add) (W2, B, one, minusone, TRUE, FALSE, cm) ;
833 CHOLMOD(free_sparse) (&W2, cm) ;
834 }
835 else
836 {
837 /* W = A*X */
838 W = CHOLMOD(ssmult) (A, X, 0, TRUE, FALSE, cm) ;
839 /* R = W - B */
840 R = CHOLMOD(add) (W, B, one, minusone, TRUE, FALSE, cm) ;
841 CHOLMOD(free_sparse) (&W, cm) ;
842 }
843
844 r = CHOLMOD(norm_sparse) (R, 1, cm) ;
845 CHOLMOD(free_sparse) (&R, cm) ;
846 if (bnorm > 0)
847 {
848 r /= bnorm ;
849 }
850 }
851
852 return (r) ;
853 }
854
855
856 /* ========================================================================== */
857 /* === resid3 =============================================================== */
858 /* ========================================================================== */
859
860 /* r = norm (A1*A2*A3*x - b) / norm (b) */
861
resid3(cholmod_sparse * A1,cholmod_sparse * A2,cholmod_sparse * A3,cholmod_dense * X,cholmod_dense * B)862 double resid3 (cholmod_sparse *A1, cholmod_sparse *A2, cholmod_sparse *A3,
863 cholmod_dense *X, cholmod_dense *B)
864 {
865 double r, bnorm ;
866 cholmod_dense *R, *W1, *W2, *X2, *B2 ;
867 cholmod_sparse *C1, *C2, *C3 ;
868 Int n, nrhs, d, xtype ;
869
870 if (A1 == NULL || X == NULL || B == NULL)
871 {
872 ERROR (CHOLMOD_INVALID, "cannot compute resid3") ;
873 return (1) ;
874 }
875
876 cm->status = CHOLMOD_OK ;
877
878 n = B->nrow ;
879
880 /* ---------------------------------------------------------------------- */
881 /* convert all inputs to an identical xtype */
882 /* ---------------------------------------------------------------------- */
883
884 xtype = MAX (A1->xtype, X->xtype) ;
885 xtype = MAX (xtype, B->xtype) ;
886 if (A2 != NULL)
887 {
888 xtype = MAX (xtype, A2->xtype) ;
889 }
890 if (A3 != NULL)
891 {
892 xtype = MAX (xtype, A3->xtype) ;
893 }
894
895 C1 = NULL ;
896 C2 = NULL ;
897 C3 = NULL ;
898 B2 = NULL ;
899 X2 = NULL ;
900
901 if (A1->xtype != xtype)
902 {
903 C1 = CHOLMOD(copy_sparse) (A1, cm) ;
904 CHOLMOD(sparse_xtype) (xtype, C1, cm) ;
905 A1 = C1 ;
906 }
907
908 if (A2 != NULL && A2->xtype != xtype)
909 {
910 C2 = CHOLMOD(copy_sparse) (A2, cm) ;
911 CHOLMOD(sparse_xtype) (xtype, C2, cm) ;
912 A2 = C2 ;
913 }
914
915 if (A3 != NULL && A3->xtype != xtype)
916 {
917 C3 = CHOLMOD(copy_sparse) (A3, cm) ;
918 CHOLMOD(sparse_xtype) (xtype, C3, cm) ;
919 A3 = C3 ;
920 }
921
922 if (X->xtype != xtype)
923 {
924 X2 = CHOLMOD(copy_dense) (X, cm) ;
925 CHOLMOD(dense_xtype) (xtype, X2, cm) ;
926 X = X2 ;
927 }
928
929 if (B->xtype != xtype)
930 {
931 B2 = CHOLMOD(copy_dense) (B, cm) ;
932 CHOLMOD(dense_xtype) (xtype, B2, cm) ;
933 B = B2 ;
934 }
935
936 if (cm->status < CHOLMOD_OK)
937 {
938 ERROR (CHOLMOD_INVALID, "cannot compute resid3") ;
939 CHOLMOD(free_sparse) (&C1, cm) ;
940 CHOLMOD(free_sparse) (&C2, cm) ;
941 CHOLMOD(free_sparse) (&C3, cm) ;
942 CHOLMOD(free_dense) (&B2, cm) ;
943 CHOLMOD(free_dense) (&X2, cm) ;
944 return (1) ;
945 }
946
947 /* ---------------------------------------------------------------------- */
948 /* get B and its norm */
949 /* ---------------------------------------------------------------------- */
950
951 nrhs = B->ncol ;
952 d = B->d ;
953 bnorm = CHOLMOD(norm_dense) (B, 1, cm) ;
954
955 /* ---------------------------------------------------------------------- */
956 /* compute the residual */
957 /* ---------------------------------------------------------------------- */
958
959 if (A3 != NULL)
960 {
961 /* W1 = A3*X */
962 W1 = CHOLMOD(zeros) (n, nrhs, xtype, cm) ;
963 CHOLMOD(sdmult) (A3, FALSE, one, zero, X, W1, cm) ;
964 }
965 else
966 {
967 W1 = X ;
968 }
969
970 if (A2 != NULL)
971 {
972 /* W2 = A2*W1 */
973 W2 = CHOLMOD(eye) (n, nrhs, xtype, cm) ;
974 CHOLMOD(sdmult) (A2, FALSE, one, zero, W1, W2, cm) ;
975 }
976 else
977 {
978 W2 = W1 ;
979 }
980
981 /* R = B */
982 R = CHOLMOD(copy_dense) (B, cm) ;
983
984 /* R = A1*W2 - R */
985 CHOLMOD(sdmult) (A1, FALSE, one, minusone, W2, R, cm) ;
986
987 /* ---------------------------------------------------------------------- */
988 /* r = norm (R) / norm (B) */
989 /* ---------------------------------------------------------------------- */
990
991 r = CHOLMOD(norm_dense) (R, 1, cm) ;
992 CHOLMOD(free_dense) (&R, cm) ;
993
994 CHOLMOD(free_sparse) (&C1, cm) ;
995 CHOLMOD(free_sparse) (&C2, cm) ;
996 CHOLMOD(free_sparse) (&C3, cm) ;
997 CHOLMOD(free_dense) (&B2, cm) ;
998 CHOLMOD(free_dense) (&X2, cm) ;
999
1000 if (A3 != NULL)
1001 {
1002 CHOLMOD(free_dense) (&W1, cm) ;
1003 }
1004 if (A2 != NULL)
1005 {
1006 CHOLMOD(free_dense) (&W2, cm) ;
1007 }
1008 if (bnorm > 0)
1009 {
1010 r /= bnorm ;
1011 }
1012 return (r) ;
1013 }
1014
1015
1016 /* ========================================================================== */
1017 /* === pnorm ================================================================ */
1018 /* ========================================================================== */
1019
1020 /* r = norm (x-Pb) or r = norm(x-P'b). This is lengthy because CHOLMOD does
1021 * not provide any operations on dense matrices, and because it used to test
1022 * the sparse-to-dense conversion routine. Multiple methods are used to compute
1023 * the same thing.
1024 */
1025
pnorm(cholmod_dense * X,Int * P,cholmod_dense * B,Int inv)1026 double pnorm (cholmod_dense *X, Int *P, cholmod_dense *B, Int inv)
1027 {
1028 cholmod_dense *R, *X2, *B2 ;
1029 cholmod_factor *L ;
1030 double *xx, *xz, *bx, *bz, *rx, *rz ;
1031 Int *Pinv, *Perm ;
1032 double rnorm, r ;
1033 Int i, j, k, n, nrhs, xtype, ok, save, lxtype ;
1034
1035 if (X == NULL || P == NULL || B == NULL)
1036 {
1037 ERROR (CHOLMOD_INVALID, "cannot compute pnorm") ;
1038 return (1) ;
1039 }
1040
1041 save = cm->prefer_zomplex ;
1042 n = X->nrow ;
1043 nrhs = X->ncol ;
1044 rnorm = 0 ;
1045
1046 Pinv = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
1047 if (Pinv != NULL)
1048 {
1049 for (k = 0 ; k < n ; k++)
1050 {
1051 Pinv [P [k]] = k ;
1052 }
1053 }
1054
1055 xtype = MAX (X->xtype, B->xtype) ;
1056
1057 R = CHOLMOD(zeros) (n, nrhs, CHOLMOD_ZOMPLEX, cm) ;
1058 B2 = CHOLMOD(copy_dense) (B, cm) ;
1059 ok = R != NULL && B2 != NULL ;
1060 ok = ok && CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, B2, cm) ;
1061
1062 for (lxtype = CHOLMOD_REAL ; ok && lxtype <= CHOLMOD_ZOMPLEX ; lxtype++)
1063 {
1064 /* create a fake factor object */
1065 L = CHOLMOD(allocate_factor) (n, cm) ;
1066 CHOLMOD(change_factor) (lxtype, TRUE, FALSE, TRUE, TRUE, L, cm) ;
1067 ok = ok && (L != NULL && L->Perm != NULL && Pinv != NULL) ;
1068 if (ok)
1069 {
1070 L->ordering = CHOLMOD_GIVEN ;
1071 Perm = L->Perm ;
1072 for (k = 0 ; k < n ; k++)
1073 {
1074 Perm [k] = Pinv [k] ;
1075 }
1076 }
1077 for (k = 0 ; k <= 1 ; k++)
1078 {
1079 /* solve the inverse permutation system, X2 = P*X or X2 = P'*X */
1080 cm->prefer_zomplex = k ;
1081 X2 = CHOLMOD(solve) (inv ? CHOLMOD_Pt : CHOLMOD_P, L, X, cm) ;
1082
1083 ok = ok && CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, X2, cm) ;
1084 if (ok && X2 != NULL)
1085 {
1086 rx = R->x ;
1087 rz = R->z ;
1088 xx = X2->x ;
1089 xz = X2->z ;
1090 bx = B2->x ;
1091 bz = B2->z ;
1092 for (j = 0 ; j < nrhs ; j++)
1093 {
1094 for (i = 0 ; i < n ; i++)
1095 {
1096 rx [i+j*n] = xx [i+j*n] - bx [i+j*n] ;
1097 rz [i+j*n] = xz [i+j*n] - bz [i+j*n] ;
1098 }
1099 }
1100 }
1101 r = CHOLMOD(norm_dense) (R, 0, cm) ;
1102 rnorm = MAX (r, rnorm) ;
1103 CHOLMOD(free_dense) (&X2, cm) ;
1104 }
1105 CHOLMOD(free_factor) (&L, cm) ;
1106 }
1107
1108 CHOLMOD(free_dense) (&B2, cm) ;
1109 CHOLMOD(free_dense) (&R, cm) ;
1110
1111 if (xtype == CHOLMOD_REAL)
1112 {
1113 /* X and B are both real */
1114 cholmod_sparse *Bs, *Pb ;
1115 Bs = CHOLMOD(dense_to_sparse) (B, TRUE, cm) ;
1116 Pb = CHOLMOD(submatrix) (Bs, inv ? Pinv : P, n, NULL, -1, TRUE, TRUE,cm);
1117 X2 = CHOLMOD(sparse_to_dense) (Pb, cm) ;
1118 R = CHOLMOD(zeros) (n, nrhs, CHOLMOD_REAL, cm) ;
1119 if (R != NULL && X != NULL && X2 != NULL)
1120 {
1121 rx = R->x ;
1122 xx = X->x ;
1123 bx = X2->x ;
1124 for (j = 0 ; j < nrhs ; j++)
1125 {
1126 for (i = 0 ; i < n ; i++)
1127 {
1128 rx [i+j*n] = xx [i+j*n] - bx [i+j*n] ;
1129 }
1130 }
1131 }
1132 CHOLMOD(free_sparse) (&Bs, cm) ;
1133 CHOLMOD(free_sparse) (&Pb, cm) ;
1134 CHOLMOD(free_dense) (&X2, cm) ;
1135 r = CHOLMOD(norm_dense) (R, 1, cm) ;
1136 rnorm = MAX (rnorm, r) ;
1137 CHOLMOD(free_dense) (&R, cm) ;
1138 }
1139
1140 CHOLMOD(free) (n, sizeof (Int), Pinv, cm) ;
1141 cm->prefer_zomplex = save ;
1142 return (rnorm) ;
1143 }
1144
1145
1146 /* ========================================================================== */
1147 /* === prune_row ============================================================ */
1148 /* ========================================================================== */
1149
1150 /* Set row k and column k of a packed matrix A to zero. Set A(k,k) to 1
1151 * if space is available. */
1152
prune_row(cholmod_sparse * A,Int k)1153 void prune_row (cholmod_sparse *A, Int k)
1154 {
1155 double *Ax ;
1156 Int *Ap, *Ai ;
1157 Int ncol, p, i, j, nz ;
1158
1159 if (A == NULL)
1160 {
1161 ERROR (CHOLMOD_INVALID, "nothing to prune") ;
1162 return ;
1163 }
1164
1165 Ap = A->p ;
1166 Ai = A->i ;
1167 Ax = A->x ;
1168 nz = 0 ;
1169 ncol = A->ncol ;
1170
1171 for (j = 0 ; j < ncol ; j++)
1172 {
1173 p = Ap [j] ;
1174 Ap [j] = nz ;
1175 if (j == k && nz < Ap [j+1])
1176 {
1177 Ai [nz] = k ;
1178 Ax [nz] = 1 ;
1179 nz++ ;
1180 }
1181 else
1182 {
1183 for ( ; p < Ap [j+1] ; p++)
1184 {
1185 i = Ai [p] ;
1186 if (i != k)
1187 {
1188 Ai [nz] = i ;
1189 Ax [nz] = Ax [p] ;
1190 nz++ ;
1191 }
1192 }
1193 }
1194 }
1195 Ap [ncol] = nz ;
1196 }
1197
1198
1199 /* ========================================================================== */
1200 /* === do_matrix =========================================================== */
1201 /* ========================================================================== */
1202
do_matrix(cholmod_sparse * A)1203 double do_matrix (cholmod_sparse *A)
1204 {
1205 double err, maxerr = 0 ;
1206 Int print, precise, maxprint, minprint, nmethods ;
1207
1208 if (A == NULL)
1209 {
1210 ERROR (CHOLMOD_INVALID, "no matrix") ;
1211 return (1) ;
1212 }
1213
1214 /* ---------------------------------------------------------------------- */
1215 /* determine print level, based on matrix size */
1216 /* ---------------------------------------------------------------------- */
1217
1218 if (A->nrow <= 4)
1219 {
1220 minprint = 5 ;
1221 maxprint = 5 ;
1222 }
1223 else if (A->nrow <= 8)
1224 {
1225 minprint = 4 ;
1226 maxprint = 4 ;
1227 }
1228 else
1229 {
1230 minprint = 1 ;
1231 maxprint = 1 ;
1232 }
1233
1234 /* ---------------------------------------------------------------------- */
1235 /* for all print levels and precisions, do: */
1236 /* ---------------------------------------------------------------------- */
1237
1238 for (print = minprint ; print <= maxprint ; print++)
1239 {
1240 for (precise = 0 ; precise <= (print >= 4) ? 1 : 0 ; precise++)
1241 {
1242 Int save1, save2 ;
1243
1244 maxerr = 0 ;
1245 my_srand (42) ; /* RAND reset */
1246 cm->print = print ;
1247 cm->precise = precise ;
1248 printf ("\n----------Print level %d precise: %d\n",
1249 cm->print, cm->precise) ;
1250
1251 save1 = cm->final_asis ;
1252 save2 = cm->final_super ;
1253 cm->final_asis = FALSE ;
1254 cm->final_super = TRUE ;
1255 OK (CHOLMOD(print_common) ("cm", cm)) ;
1256 cm->final_asis = save1 ;
1257 cm->final_super = save2 ;
1258
1259 /* -------------------------------------------------------------- */
1260 /* test various matrix operations */
1261 /* -------------------------------------------------------------- */
1262
1263 err = test_ops (A) ; /* RAND */
1264 MAXERR (maxerr, err, 1) ;
1265
1266 /* -------------------------------------------------------------- */
1267 /* solve the augmented system */
1268 /* -------------------------------------------------------------- */
1269
1270 err = aug (A) ; /* no random number use */
1271 MAXERR (maxerr, err, 1) ;
1272
1273 /* -------------------------------------------------------------- */
1274 /* solve using different methods */
1275 /* -------------------------------------------------------------- */
1276
1277 printf ("test_solver (1)\n") ;
1278 cm->nmethods = 9 ;
1279 cm->final_asis = TRUE ;
1280 err = test_solver (A) ; /* RAND reset */
1281 MAXERR (maxerr, err, 1) ;
1282
1283 printf ("test_solver (2)\n") ;
1284 cm->final_asis = TRUE ;
1285 for (nmethods = 0 ; nmethods < 7 ; nmethods++)
1286 {
1287 cm->nmethods = nmethods ;
1288 cm->method [0].ordering = CHOLMOD_NATURAL ;
1289 err = test_solver (A) ; /* RAND reset */
1290 MAXERR (maxerr, err, 1) ;
1291 }
1292
1293 printf ("test_solver (3)\n") ;
1294 cm->nmethods = 1 ;
1295 cm->method [0].ordering = CHOLMOD_NESDIS ;
1296 err = test_solver (A) ; /* RAND reset */
1297 MAXERR (maxerr, err, 1) ;
1298
1299 printf ("test_solver (3b)\n") ;
1300 cm->nmethods = 1 ;
1301 cm->method [0].ordering = CHOLMOD_NESDIS ;
1302 cm->method [0].nd_camd = 2 ;
1303 err = test_solver (A) ; /* RAND reset */
1304 MAXERR (maxerr, err, 1) ;
1305
1306 printf ("test_solver (3c)\n") ;
1307 cm->nmethods = 1 ;
1308 cm->method [0].ordering = CHOLMOD_NATURAL ;
1309 err = test_solver (A) ; /* RAND reset */
1310 MAXERR (maxerr, err, 1) ;
1311
1312 printf ("test_solver (4)\n") ;
1313 cm->nmethods = 1 ;
1314 cm->method[0].ordering = CHOLMOD_METIS ;
1315 err = test_solver (A) ; /* RAND reset */
1316 MAXERR (maxerr, err, 1) ;
1317
1318 printf ("test_solver (5)\n") ;
1319 cm->nmethods = 1 ;
1320 cm->method [0].ordering = CHOLMOD_AMD ;
1321 CHOLMOD(free_work) (cm) ;
1322 err = test_solver (A) ; /* RAND reset */
1323 MAXERR (maxerr, err, 1) ;
1324
1325 printf ("test_solver (6)\n") ;
1326 cm->nmethods = 1 ;
1327 cm->method[0].ordering = CHOLMOD_COLAMD ;
1328 err = test_solver (A) ; /* RAND reset */
1329 MAXERR (maxerr, err, 1) ;
1330
1331 /* -------------------------------------------------------------- */
1332 /* restore default control parameters */
1333 /* -------------------------------------------------------------- */
1334
1335 OK (CHOLMOD(print_common) ("cm", cm)) ;
1336 CHOLMOD(defaults) (cm) ; cm->useGPU = 0 ;
1337 }
1338 }
1339
1340 printf ("do_matrix max error %.1g\n", maxerr) ;
1341
1342 return (maxerr) ;
1343 }
1344
1345
1346 /* ========================================================================== */
1347 /* === main ================================================================= */
1348 /* ========================================================================== */
1349
1350 /* Usage:
1351 * cm < matrix do not perform intensive memory-failure tests
1352 * cm -m < matrix do perform memory tests
1353 * cm -s < matrix matrix is singular, nan error expected
1354 *
1355 * (The memory tests are performed if any argument is given to cm).
1356 */
1357
main(int argc,char ** argv)1358 int main (int argc, char **argv)
1359 {
1360 cholmod_triplet *T ;
1361 cholmod_sparse *A, *C, *AT ;
1362 char *s ;
1363 double err = 0, maxerr = 0 ;
1364 Int n = 0, nmin = 0, nrow = 0, ncol = 0, save ;
1365 int singular, do_memory, i, do_nantests, ok ;
1366 double v = CHOLMOD_VERSION, tic [2], t ;
1367 int version [3] ;
1368 char *p ;
1369 const char* env_use_gpu;
1370
1371 SuiteSparse_start ( ) ;
1372 SuiteSparse_tic (tic) ;
1373 printf ("Testing CHOLMOD (%g): %d ", v, CHOLMOD(version) (version)) ;
1374 printf ("(%d.%d.%d)\n", version [0], version [1], version [2]) ;
1375 v = SUITESPARSE_VERSION ;
1376 printf ("for SuiteSparse (%g): %d ", v, SuiteSparse_version (version)) ;
1377 printf ("(%d.%d.%d)\n", version [0], version [1], version [2]) ;
1378 printf ("%s: argc: %d\n", argv [0], argc) ;
1379 my_srand (42) ; /* RAND */
1380
1381 /* Ignore floating point exceptions. Infs and NaNs are generated
1382 on purpose. */
1383 signal (SIGFPE, SIG_IGN) ;
1384
1385 /* query the CHOLMOD_USE_GPU environment variable */
1386 env_use_gpu = getenv ("CHOLMOD_USE_GPU") ;
1387 if ( env_use_gpu )
1388 {
1389 /* CHOLMOD_USE_GPU environment variable is set to something */
1390 if ( atoi ( env_use_gpu ) == 0 )
1391 {
1392 printf ("CHOLMOD_USE_GPU 0\n") ;
1393 }
1394 else
1395 {
1396 printf ("CHOLMOD_USE_GPU 1 (ignored for this test)\n") ;
1397 }
1398 }
1399 else
1400 {
1401 printf ("CHOLMOD_USE_GPU not present\n") ;
1402 }
1403
1404 fflush (stdout) ;
1405
1406 singular = FALSE ;
1407 do_memory = FALSE ;
1408 do_nantests = FALSE ;
1409 for (i = 1 ; i < argc ; i++)
1410 {
1411 s = argv [i] ;
1412 if (s [0] == '-' && s [1] == 'm') do_memory = TRUE ;
1413 if (s [0] == '-' && s [1] == 's') singular = TRUE ;
1414 if (s [0] == '-' && s [1] == 'n') do_nantests = TRUE ;
1415 }
1416
1417 printf ("do_memory: %d singular: %d\n", do_memory, singular) ;
1418
1419 /* ---------------------------------------------------------------------- */
1420 /* test SuiteSparse malloc functions */
1421 /* ---------------------------------------------------------------------- */
1422
1423 p = SuiteSparse_malloc (0, 0) ;
1424 OKP (p) ;
1425 p [0] = 'a' ;
1426 SuiteSparse_free (p) ;
1427 p = SuiteSparse_calloc (0, 0) ;
1428 OKP (p) ;
1429 p [0] = 'a' ;
1430 p = SuiteSparse_realloc (0, 0, 0, p, &ok) ;
1431 OK (ok) ;
1432 OKP (p) ;
1433 p [0] = 'a' ;
1434 SuiteSparse_free (p) ;
1435 p = SuiteSparse_malloc (SuiteSparse_long_max, 1024) ;
1436 NOP (p) ;
1437 p = SuiteSparse_calloc (SuiteSparse_long_max, 1024) ;
1438 NOP (p) ;
1439 p = SuiteSparse_realloc (0, 0, 0, NULL, &ok) ;
1440 OK (ok) ;
1441 OKP (p) ;
1442 p [0] = 'a' ;
1443 SuiteSparse_free (p) ;
1444 p = SuiteSparse_realloc (SuiteSparse_long_max, 0, 1024, NULL, &ok) ;
1445 NOP (p) ;
1446 NOT (ok) ;
1447
1448 /* ---------------------------------------------------------------------- */
1449 /* initialize CHOLMOD */
1450 /* ---------------------------------------------------------------------- */
1451
1452 cm = &Common ;
1453 OK (CHOLMOD(start) (cm)) ; cm->useGPU = 0 ;
1454
1455 /* ---------------------------------------------------------------------- */
1456 /* test all methods with NULL common */
1457 /* ---------------------------------------------------------------------- */
1458
1459 /* no user error handler, since lots of errors will be raised here */
1460 cm->error_handler = NULL ;
1461 null_test (NULL) ;
1462 save = cm->itype ;
1463 cm->itype = -999 ;
1464 null_test (cm) ;
1465 cm->itype = save ;
1466 null_test2 ( ) ;
1467 CHOLMOD(finish) (cm) ;
1468 OK (cm->malloc_count == 0) ;
1469 OK (CHOLMOD(start) (cm)) ; cm->useGPU = 0 ;
1470
1471 /* ---------------------------------------------------------------------- */
1472 /* create basic scalars */
1473 /* ---------------------------------------------------------------------- */
1474
1475 Zero [0] = 0 ;
1476 Zero [1] = 0 ;
1477
1478 zero [0] = 0 ;
1479 zero [1] = 0 ;
1480 one [0] = 1 ;
1481 one [1] = 0 ;
1482 minusone [0] = -1 ;
1483 minusone [1] = 0 ;
1484 M1 = CHOLMOD(ones) (1, 1, CHOLMOD_REAL, cm) ;
1485
1486 if (M1 != NULL)
1487 {
1488 ((double *) (M1->x)) [0] = -1 ;
1489 }
1490
1491 /* ---------------------------------------------------------------------- */
1492 /* read in a triplet matrix and use it to test CHOLMOD */
1493 /* ---------------------------------------------------------------------- */
1494
1495 for ( ; (T = read_triplet (stdin)) != NULL ; ) /* RAND */
1496 {
1497
1498 if (T->nrow > 1000000)
1499 {
1500 /* do huge-problem tests only, but only for 32-bit systems */
1501 if (sizeof (Int) == sizeof (int))
1502 {
1503 huge ( ) ;
1504 }
1505 CHOLMOD(free_triplet) (&T, cm) ;
1506 continue ;
1507 }
1508
1509 maxerr = 0 ;
1510 CHOLMOD(defaults) (cm) ; cm->useGPU = 0 ;
1511 cm->error_handler = my_handler ;
1512 cm->print = 4 ;
1513 cm->precise = FALSE ;
1514 cm->nmethods = 8 ;
1515
1516 /* ------------------------------------------------------------------ */
1517 /* convert triplet to sparse matrix */
1518 /* ------------------------------------------------------------------ */
1519
1520 A = CHOLMOD(triplet_to_sparse) (T, 0, cm) ;
1521 AT = CHOLMOD(transpose) (A, 0, cm) ;
1522 OK (CHOLMOD(print_sparse) (A, "Test matrix, A", cm)) ;
1523 C = unpack (A) ; /* RAND */
1524 OK (CHOLMOD(print_sparse) (C, "Unpacked/unsorted version of A", cm)) ;
1525 cm->print = 1 ;
1526
1527 if (T != NULL)
1528 {
1529 nrow = T->nrow ;
1530 ncol = T->ncol ;
1531 n = MAX (nrow, ncol) ;
1532 nmin = MIN (nrow, ncol) ;
1533 }
1534
1535 /* ------------------------------------------------------------------ */
1536 /* basic error tests */
1537 /* ------------------------------------------------------------------ */
1538
1539 null2 (T, do_nantests) ; /* RAND */
1540 printf ("Null2 OK : no error\n") ;
1541 if (do_nantests)
1542 {
1543 maxerr = 0 ;
1544 goto done ; /* yes, this is ugly */
1545 }
1546
1547 /* ------------------------------------------------------------------ */
1548 /* raw factorization tests */
1549 /* ------------------------------------------------------------------ */
1550
1551 cm->error_handler = NULL ;
1552 err = raw_factor (A, 2) ; /* RAND */
1553 cm->error_handler = my_handler ;
1554 MAXERR (maxerr, err, 1) ;
1555 printf ("raw factorization error %.1g\n", err) ;
1556
1557 err = raw_factor2 (A, 0., 0) ;
1558 MAXERR (maxerr, err, 1) ;
1559 printf ("raw factorization error2 %.1g\n", err) ;
1560
1561 err = raw_factor2 (A, 1e-16, 0) ;
1562 MAXERR (maxerr, err, 1) ;
1563 printf ("raw factorization error3 %.1g\n", err) ;
1564
1565 if (n < 1000 && A && T && A->stype == 1)
1566 {
1567 /* factorize a symmetric matrix (upper part stored), possibly
1568 * with ignored entries in lower triangular part. */
1569 cholmod_sparse *F ;
1570 int save = T->stype ;
1571
1572 T->stype = 0 ;
1573 F = CHOLMOD(triplet_to_sparse) (T, 0, cm) ;
1574 T->stype = save ;
1575
1576 /*
1577 ET = CHOLMOD(transpose) (E, 2, cm) ;
1578 if (E) E->stype = 0 ;
1579 if (ET) ET->stype = 0 ;
1580 F = CHOLMOD(add) (E, ET, one, one, 1, 1, cm) ;
1581 */
1582
1583 if (F) F->stype = 1 ;
1584
1585 err = raw_factor2 (F, 0., 0) ;
1586 MAXERR (maxerr, err, 1) ;
1587 printf ("raw factorization error4 %.1g\n", err) ;
1588
1589 err = raw_factor2 (F, 1e-16, 0) ;
1590 MAXERR (maxerr, err, 1) ;
1591 printf ("raw factorization error5 %.1g\n", err) ;
1592
1593 cm->dbound = 1e-15 ;
1594
1595 err = raw_factor2 (F, 0., 0) ;
1596 MAXERR (maxerr, err, 1) ;
1597 printf ("raw factorization error6 %.1g\n", err) ;
1598
1599 err = raw_factor2 (F, 1e-16, 0) ;
1600 MAXERR (maxerr, err, 1) ;
1601 printf ("raw factorization error7 %.1g\n", err) ;
1602
1603 err = raw_factor2 (F, 1e-16, 1) ;
1604 MAXERR (maxerr, err, 1) ;
1605 printf ("raw factorization error8 %.1g\n", err) ;
1606
1607 cm->dbound = 0 ;
1608
1609 /*
1610 CHOLMOD(free_sparse) (&E, cm) ;
1611 CHOLMOD(free_sparse) (&ET, cm) ;
1612 */
1613
1614 CHOLMOD(free_sparse) (&F, cm) ;
1615 }
1616
1617 /* ------------------------------------------------------------------ */
1618 /* matrix ops */
1619 /* ------------------------------------------------------------------ */
1620
1621 err = test_ops (A) ; /* RAND */
1622 MAXERR (maxerr, err, 1) ;
1623 printf ("initial testops error %.1g\n", err) ;
1624
1625 /* ------------------------------------------------------------------ */
1626 /* analyze, factorize, solve */
1627 /* ------------------------------------------------------------------ */
1628
1629 err = solve (A) ; /* RAND */
1630 MAXERR (maxerr, err, 1) ;
1631 printf ("initial solve error %.1g\n", err) ;
1632
1633 /* ------------------------------------------------------------------ */
1634 /* CCOLAMD tests */
1635 /* ------------------------------------------------------------------ */
1636
1637 cctest (A) ; /* RAND reset */
1638 cctest (AT) ; /* RAND reset */
1639
1640 /* ------------------------------------------------------------------ */
1641 /* COLAMD tests */
1642 /* ------------------------------------------------------------------ */
1643
1644 ctest (A) ;
1645 ctest (AT) ;
1646
1647 /* ------------------------------------------------------------------ */
1648 /* AMD tests */
1649 /* ------------------------------------------------------------------ */
1650
1651 if (n < NLARGE || A->stype)
1652 {
1653 /* for unsymmetric matrices, this forms A*A' and A'*A, which can
1654 * fail if A has a dense row or column. So it is only done for
1655 * modest-sized unsymmetric matrices. */
1656 amdtest (A) ;
1657 amdtest (AT) ;
1658 camdtest (A) ; /* RAND */
1659 camdtest (AT) ; /* RAND */
1660 }
1661
1662 if (n < NSMALL)
1663 {
1664
1665 /* -------------------------------------------------------------- */
1666 /* do_matrix with an unpacked matrix */
1667 /* -------------------------------------------------------------- */
1668
1669 /* try with an unpacked matrix, and a non-default dbound */
1670 cm->dbound = 1e-15 ;
1671 err = do_matrix (C) ; /* RAND reset */
1672 MAXERR (maxerr, err, 1) ;
1673 cm->dbound = 0 ;
1674
1675 /* -------------------------------------------------------------- */
1676 /* do_matrix: analyze, factorize, and solve, with many options */
1677 /* -------------------------------------------------------------- */
1678
1679 err = do_matrix (A) ; /* RAND reset */
1680 MAXERR (maxerr, err, 1) ;
1681
1682 /* -------------------------------------------------------------- */
1683 /* pretend to solve an LP */
1684 /* -------------------------------------------------------------- */
1685
1686 if (nrow != ncol)
1687 {
1688 cm->print = 2 ;
1689 err = lpdemo (T) ; /* RAND */
1690 cm->print = 1 ;
1691 MAXERR (maxerr, err, 1) ;
1692 cm->print = 5; CHOLMOD(print_common) ("Common", cm);cm->print=1;
1693 cm->nmethods = 1 ;
1694 cm->method [0].ordering = CHOLMOD_COLAMD ;
1695 err = lpdemo (T) ; /* RAND */
1696 MAXERR (maxerr, err, 1) ;
1697 printf ("initial lp error %.1g, dbound %g\n", err, cm->dbound) ;
1698 cm->nmethods = 0 ;
1699 cm->method [0].ordering = CHOLMOD_GIVEN ;
1700 }
1701 }
1702
1703 progress (1, '|') ;
1704
1705 if (n < NSMALL && do_memory)
1706 {
1707
1708 /* -------------------------------------------------------------- */
1709 /* Exhaustive memory-error handling */
1710 /* -------------------------------------------------------------- */
1711
1712 memory_tests (T) ; /* RAND */
1713 }
1714
1715 /* ------------------------------------------------------------------ */
1716 /* free matrices and print results */
1717 /* ------------------------------------------------------------------ */
1718
1719 done: /* an ugly "goto" target; added to minimize code changes
1720 * when added "do_nantests", for version 1.0.2 */
1721
1722 CHOLMOD(free_sparse) (&C, cm) ;
1723 CHOLMOD(free_sparse) (&AT, cm) ;
1724 CHOLMOD(free_sparse) (&A, cm) ;
1725 CHOLMOD(free_triplet) (&T, cm) ;
1726
1727 fprintf (stderr, "\n "
1728 " Test OK") ;
1729 if (nrow <= ncol && !singular)
1730 {
1731 /* maxerr should be NaN if nrow > ncol, so don't print it */
1732 fprintf (stderr, ", maxerr %.1g", maxerr) ;
1733 }
1734 fprintf (stderr, "\n") ;
1735 my_srand (42) ; /* RAND reset */
1736 }
1737
1738 /* ---------------------------------------------------------------------- */
1739 /* finalize CHOLMOD */
1740 /* ---------------------------------------------------------------------- */
1741
1742 CHOLMOD(free_dense) (&M1, cm) ;
1743 CHOLMOD(finish) (cm) ;
1744
1745 cm->print = 5 ; OK (CHOLMOD(print_common) ("cm", cm)) ;
1746 printf ("malloc count "ID" inuse "ID"\n",
1747 (Int) cm->malloc_count,
1748 (Int) cm->memory_inuse) ;
1749 OK (cm->malloc_count == 0) ;
1750 OK (cm->memory_inuse == 0) ;
1751 t = SuiteSparse_toc (tic) ;
1752 if (nrow > ncol)
1753 {
1754 /* maxerr should be NaN, so don't print it */
1755 printf ("All tests successful: time %.1g\n", t) ;
1756 }
1757 else
1758 {
1759 printf ("All tests successful: max error %.1g time: %.1g\n", maxerr, t);
1760 }
1761
1762 SuiteSparse_finish ( ) ;
1763 return (0) ;
1764 }
1765