1 /* ========================================================================== */
2 /* === Check/cholmod_check ================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Check Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Routines to check and print the contents of the 5 CHOLMOD objects:
10 *
11 * No CHOLMOD routine calls the check or print routines. If a user wants to
12 * check CHOLMOD's input parameters, a separate call to the appropriate check
13 * routine should be used before calling other CHOLMOD routines.
14 *
15 * cholmod_check_common check statistics and workspace in Common
16 * cholmod_check_sparse check sparse matrix in compressed column form
17 * cholmod_check_dense check dense matrix
18 * cholmod_check_factor check factorization
19 * cholmod_check_triplet check sparse matrix in triplet form
20 *
21 * cholmod_print_common print statistics in Common
22 * cholmod_print_sparse print sparse matrix in compressed column form
23 * cholmod_print_dense print dense matrix
24 * cholmod_print_factor print factorization
25 * cholmod_print_triplet print sparse matrix in triplet form
26 *
27 * In addition, this file contains routines to check and print three types of
28 * integer vectors:
29 *
30 * cholmod_check_perm check a permutation of 0:n-1 (no duplicates)
31 * cholmod_check_subset check a subset of 0:n-1 (duplicates OK)
32 * cholmod_check_parent check an elimination tree
33 *
34 * cholmod_print_perm print a permutation
35 * cholmod_print_subset print a subset
36 * cholmod_print_parent print an elimination tree
37 *
38 * Each Common->print level prints the items at or below the given level:
39 *
40 * 0: print nothing; just check the data structures and return TRUE/FALSE
41 * 1: error messages
42 * 2: warning messages
43 * 3: one-line summary of each object printed
44 * 4: short summary of each object (first and last few entries)
45 * 5: entire contents of the object
46 *
47 * No CHOLMOD routine calls these routines, so no printing occurs unless
48 * the user specifically calls a cholmod_print_* routine. Thus, the default
49 * print level is 3.
50 *
51 * Common->precise controls the # of digits printed for numerical entries
52 * (5 if FALSE, 15 if TRUE).
53 *
54 * If SuiteSparse_config.printf_func is NULL, then no printing occurs. The
55 * cholmod_check_* and cholmod_print_* routines still check their inputs and
56 * return TRUE/FALSE if the object is valid or not.
57 *
58 * This file also includes debugging routines that are enabled only when
59 * NDEBUG is defined in cholmod_internal.h (cholmod_dump_*).
60 */
61
62 #ifndef NCHECK
63
64 #include "cholmod_internal.h"
65 #include "cholmod_check.h"
66
67 /* ========================================================================== */
68 /* === printing definitions ================================================= */
69 /* ========================================================================== */
70
71 #ifdef LONG
72 #define I8 "%8ld"
73 #define I_8 "%-8ld"
74 #else
75 #define I8 "%8d"
76 #define I_8 "%-8d"
77 #endif
78
79 #define PR(i,format,arg) \
80 { \
81 if (print >= i && SuiteSparse_config.printf_func != NULL) \
82 { \
83 SuiteSparse_config.printf_func (format, arg) ; \
84 } \
85 }
86
87 #define P1(format,arg) PR(1,format,arg)
88 #define P2(format,arg) PR(2,format,arg)
89 #define P3(format,arg) PR(3,format,arg)
90 #define P4(format,arg) PR(4,format,arg)
91
92 #define ERR(msg) \
93 { \
94 P1 ("\nCHOLMOD ERROR: %s: ", type) ; \
95 if (name != NULL) \
96 { \
97 P1 ("%s", name) ; \
98 } \
99 P1 (": %s\n", msg) ; \
100 ERROR (CHOLMOD_INVALID, "invalid") ; \
101 return (FALSE) ; \
102 }
103
104 /* print a numerical value */
105 #define PRINTVALUE(value) \
106 { \
107 if (Common->precise) \
108 { \
109 P4 (" %23.15e", value) ; \
110 } \
111 else \
112 { \
113 P4 (" %.5g", value) ; \
114 } \
115 }
116
117 /* start printing */
118 #define ETC_START(count,limit) \
119 { \
120 count = (init_print == 4) ? (limit) : (-1) ; \
121 }
122
123 /* re-enable printing if condition is met */
124 #define ETC_ENABLE(condition,count,limit) \
125 { \
126 if ((condition) && init_print == 4) \
127 { \
128 count = limit ; \
129 print = 4 ; \
130 } \
131 }
132
133 /* turn off printing if limit is reached */
134 #define ETC_DISABLE(count) \
135 { \
136 if ((count >= 0) && (count-- == 0) && print == 4) \
137 { \
138 P4 ("%s", " ...\n") ; \
139 print = 3 ; \
140 } \
141 }
142
143 /* re-enable printing, or turn if off after limit is reached */
144 #define ETC(condition,count,limit) \
145 { \
146 ETC_ENABLE (condition, count, limit) ; \
147 ETC_DISABLE (count) ; \
148 }
149
150 #define BOOLSTR(x) ((x) ? "true " : "false")
151
152 /* ========================================================================== */
153 /* === print_value ========================================================== */
154 /* ========================================================================== */
155
print_value(Int print,Int xtype,double * Xx,double * Xz,Int p,cholmod_common * Common)156 static void print_value
157 (
158 Int print,
159 Int xtype,
160 double *Xx,
161 double *Xz,
162 Int p,
163 cholmod_common *Common)
164 {
165 if (xtype == CHOLMOD_REAL)
166 {
167 PRINTVALUE (Xx [p]) ;
168 }
169 else if (xtype == CHOLMOD_COMPLEX)
170 {
171 P4 ("%s", "(") ;
172 PRINTVALUE (Xx [2*p ]) ;
173 P4 ("%s", " , ") ;
174 PRINTVALUE (Xx [2*p+1]) ;
175 P4 ("%s", ")") ;
176 }
177 else if (xtype == CHOLMOD_ZOMPLEX)
178 {
179 P4 ("%s", "(") ;
180 PRINTVALUE (Xx [p]) ;
181 P4 ("%s", " , ") ;
182 PRINTVALUE (Xz [p]) ;
183 P4 ("%s", ")") ;
184 }
185 }
186
187 /* ========================================================================== */
188 /* === cholmod_check_common ================================================= */
189 /* ========================================================================== */
190
191 /* Print and verify the contents of Common */
192
check_common(Int print,const char * name,cholmod_common * Common)193 static int check_common
194 (
195 Int print,
196 const char *name,
197 cholmod_common *Common
198 )
199 {
200 double fl, lnz ;
201 double *Xwork ;
202 Int *Flag, *Head ;
203 SuiteSparse_long mark ;
204 Int i, nrow, nmethods, ordering, xworksize, amd_backup, init_print ;
205 const char *type = "common" ;
206
207 /* ---------------------------------------------------------------------- */
208 /* print control parameters and statistics */
209 /* ---------------------------------------------------------------------- */
210
211 RETURN_IF_NULL_COMMON (FALSE) ;
212 init_print = print ;
213
214 P2 ("%s", "\n") ;
215
216 P1 ("CHOLMOD version %d", CHOLMOD_MAIN_VERSION) ;
217 P1 (".%d", CHOLMOD_SUB_VERSION) ;
218 P1 (".%d", CHOLMOD_SUBSUB_VERSION) ;
219 P1 (", %s: ", CHOLMOD_DATE) ;
220
221 if (name != NULL)
222 {
223 P1 ("%s: ", name) ;
224 }
225 switch (Common->status)
226 {
227
228 case CHOLMOD_OK:
229 P1 ("%s", "status: OK\n") ;
230 break ;
231
232 case CHOLMOD_OUT_OF_MEMORY:
233 P1 ("%s", "status: ERROR, out of memory\n") ;
234 break ;
235
236 case CHOLMOD_INVALID:
237 P1 ("%s", "status: ERROR, invalid parameter\n") ;
238 break ;
239
240 case CHOLMOD_TOO_LARGE:
241 P1 ("%s", "status: ERROR, problem too large\n") ;
242 break ;
243
244 case CHOLMOD_NOT_INSTALLED:
245 P1 ("%s", "status: ERROR, method not installed\n") ;
246 break ;
247
248 case CHOLMOD_GPU_PROBLEM:
249 P1 ("%s", "status: ERROR, GPU had a fatal error\n") ;
250 break ;
251
252 case CHOLMOD_NOT_POSDEF:
253 P1 ("%s", "status: warning, matrix not positive definite\n") ;
254 break ;
255
256 case CHOLMOD_DSMALL:
257 P1 ("%s", "status: warning, diagonal entry has tiny abs. value\n") ;
258 break ;
259
260 default:
261 ERR ("unknown status") ;
262 }
263
264 P2 (" Architecture: %s\n", CHOLMOD_ARCHITECTURE) ;
265 P3 (" sizeof(int): %d\n", (int) sizeof (int)) ;
266 P3 (" sizeof(SuiteSparse_long): %d\n", (int) sizeof (SuiteSparse_long));
267 P3 (" sizeof(void *): %d\n", (int) sizeof (void *)) ;
268 P3 (" sizeof(double): %d\n", (int) sizeof (double)) ;
269 P3 (" sizeof(Int): %d (CHOLMOD's basic integer)\n", (int) sizeof (Int)) ;
270 P3 (" sizeof(BLAS_INT): %d (integer used in the BLAS)\n",
271 (int) sizeof (BLAS_INT)) ;
272
273 if (Common->fl != EMPTY)
274 {
275 P2 ("%s", " Results from most recent analysis:\n") ;
276 P2 (" Cholesky flop count: %.5g\n", Common->fl) ;
277 P2 (" Nonzeros in L: %.5g\n", Common->lnz) ;
278 }
279 if (Common->modfl != EMPTY)
280 {
281 P2 (" Update/downdate flop count: %.5g\n", Common->modfl) ;
282 }
283
284 P2 (" memory blocks in use: %8.0f\n", (double) (Common->malloc_count)) ;
285 P2 (" memory in use (MB): %8.1f\n",
286 (double) (Common->memory_inuse) / 1048576.) ;
287 P2 (" peak memory usage (MB): %8.1f\n",
288 (double) (Common->memory_usage) / 1048576.) ;
289
290 /* ---------------------------------------------------------------------- */
291 /* primary control parameters and related ordering statistics */
292 /* ---------------------------------------------------------------------- */
293
294 P3 (" maxrank: update/downdate rank: "ID"\n",
295 (Int) CHOLMOD(maxrank) (0, Common)) ;
296 P3 (" supernodal control: %d", Common->supernodal) ;
297 P3 (" %g ", Common->supernodal_switch) ;
298 if (Common->supernodal <= CHOLMOD_SIMPLICIAL)
299 {
300 P3 ("%s", "(always do simplicial)\n") ;
301 }
302 else if (Common->supernodal == CHOLMOD_AUTO)
303 {
304 P3 ("(supernodal if flops/lnz >= %g)\n", Common->supernodal_switch) ;
305 }
306 else
307 {
308 P3 ("%s", "(always do supernodal)\n") ;
309 }
310
311 nmethods = MIN (Common->nmethods, CHOLMOD_MAXMETHODS) ;
312 nmethods = MAX (0, nmethods) ;
313
314 if (nmethods > 0)
315 {
316 P3 ("%s", " nmethods: number of ordering methods to try: ") ;
317 P3 (""ID"\n", nmethods) ;
318 amd_backup = (nmethods > 1) || (nmethods == 1 &&
319 (Common->method [0].ordering == CHOLMOD_METIS ||
320 Common->method [0].ordering == CHOLMOD_NESDIS)) ;
321 }
322 else
323 {
324 P3 ("%s", " nmethods=0: default strategy: Try user permutation if "
325 "given. Try AMD.\n") ;
326 #ifndef NPARTITION
327 if (Common->default_nesdis)
328 {
329 P3 ("%s", " Try NESDIS if AMD reports flops/nnz(L) >= 500 and "
330 "nnz(L)/nnz(A) >= 5.\n") ;
331 }
332 else
333 {
334 P3 ("%s", " Try METIS if AMD reports flops/nnz(L) >= 500 and "
335 "nnz(L)/nnz(A) >= 5.\n") ;
336 }
337 #endif
338 P3 ("%s", " Select best ordering tried.\n") ;
339 Common->method [0].ordering = CHOLMOD_GIVEN ;
340 Common->method [1].ordering = CHOLMOD_AMD ;
341 Common->method [2].ordering =
342 (Common->default_nesdis ? CHOLMOD_NESDIS : CHOLMOD_METIS) ;
343 amd_backup = FALSE ;
344 #ifndef NPARTITION
345 nmethods = 3 ;
346 #else
347 nmethods = 2 ;
348 #endif
349 }
350
351 for (i = 0 ; i < nmethods ; i++)
352 {
353 P3 (" method "ID": ", i) ;
354 ordering = Common->method [i].ordering ;
355 fl = Common->method [i].fl ;
356 lnz = Common->method [i].lnz ;
357 switch (ordering)
358 {
359
360 case CHOLMOD_NATURAL:
361 P3 ("%s", "natural\n") ;
362 break ;
363
364 case CHOLMOD_GIVEN:
365 P3 ("%s", "user permutation (if given)\n") ;
366 break ;
367
368 case CHOLMOD_AMD:
369 P3 ("%s", "AMD (or COLAMD if factorizing AA')\n") ;
370 amd_backup = FALSE ;
371 break ;
372
373 case CHOLMOD_COLAMD:
374 P3 ("%s", "AMD if factorizing A, COLAMD if factorizing AA')\n");
375 amd_backup = FALSE ;
376 break ;
377
378 case CHOLMOD_METIS:
379 P3 ("%s", "METIS_NodeND nested dissection\n") ;
380 break ;
381
382 case CHOLMOD_NESDIS:
383 P3 ("%s", "CHOLMOD nested dissection\n") ;
384
385 P3 (" nd_small: # nodes in uncut subgraph: "ID"\n",
386 (Int) (Common->method [i].nd_small)) ;
387 P3 (" nd_compress: compress the graph: %s\n",
388 BOOLSTR (Common->method [i].nd_compress)) ;
389 P3 (" nd_camd: use constrained min degree: %s\n",
390 BOOLSTR (Common->method [i].nd_camd)) ;
391 break ;
392
393 default:
394 P3 (ID, ordering) ;
395 ERR ("unknown ordering method") ;
396 break ;
397
398 }
399
400 if (!(ordering == CHOLMOD_NATURAL || ordering == CHOLMOD_GIVEN))
401 {
402 if (Common->method [i].prune_dense < 0)
403 {
404 P3 (" prune_dense: for pruning dense nodes: %s\n",
405 " none pruned") ;
406 }
407 else
408 {
409 P3 (" prune_dense: for pruning dense nodes: "
410 "%.5g\n",
411 Common->method [i].prune_dense) ;
412 P3 (" a dense node has degree "
413 ">= max(16,(%.5g)*sqrt(n))\n",
414 Common->method [i].prune_dense) ;
415 }
416 }
417
418 if (ordering == CHOLMOD_COLAMD || ordering == CHOLMOD_NESDIS)
419 {
420 if (Common->method [i].prune_dense2 < 0)
421 {
422 P3 (" prune_dense2: for pruning dense rows for AA':"
423 " %s\n", " none pruned") ;
424 }
425 else
426 {
427 P3 (" prune_dense2: for pruning dense rows for AA':"
428 " %.5g\n", Common->method [i].prune_dense2) ;
429 P3 (" a dense row has degree "
430 ">= max(16,(%.5g)*sqrt(ncol))\n",
431 Common->method [i].prune_dense2) ;
432 }
433 }
434
435 if (fl != EMPTY) P3 (" flop count: %.5g\n", fl) ;
436 if (lnz != EMPTY) P3 (" nnz(L): %.5g\n", lnz) ;
437 }
438
439 /* backup AMD results, if any */
440 if (amd_backup)
441 {
442 P3 ("%s", " backup method: ") ;
443 P3 ("%s", "AMD (or COLAMD if factorizing AA')\n") ;
444 fl = Common->method [nmethods].fl ;
445 lnz = Common->method [nmethods].lnz ;
446 if (fl != EMPTY) P3 (" AMD flop count: %.5g\n", fl) ;
447 if (lnz != EMPTY) P3 (" AMD nnz(L): %.5g\n", lnz) ;
448 }
449
450 /* ---------------------------------------------------------------------- */
451 /* arcane control parameters */
452 /* ---------------------------------------------------------------------- */
453
454 if (Common->final_asis)
455 {
456 P4 ("%s", " final_asis: TRUE, leave as is\n") ;
457 }
458 else
459 {
460 P4 ("%s", " final_asis: FALSE, convert when done\n") ;
461 if (Common->final_super)
462 {
463 P4 ("%s", " final_super: TRUE, leave in supernodal form\n") ;
464 }
465 else
466 {
467 P4 ("%s", " final_super: FALSE, convert to simplicial form\n") ;
468 }
469 if (Common->final_ll)
470 {
471 P4 ("%s", " final_ll: TRUE, convert to LL' form\n") ;
472 }
473 else
474 {
475 P4 ("%s", " final_ll: FALSE, convert to LDL' form\n") ;
476 }
477 if (Common->final_pack)
478 {
479 P4 ("%s", " final_pack: TRUE, pack when done\n") ;
480 }
481 else
482 {
483 P4 ("%s", " final_pack: FALSE, do not pack when done\n") ;
484 }
485 if (Common->final_monotonic)
486 {
487 P4 ("%s", " final_monotonic: TRUE, ensure L is monotonic\n") ;
488 }
489 else
490 {
491 P4 ("%s",
492 " final_monotonic: FALSE, do not ensure L is monotonic\n") ;
493 }
494 P4 (" final_resymbol: remove zeros from amalgamation: %s\n",
495 BOOLSTR (Common->final_resymbol)) ;
496 }
497
498 P4 (" dbound: LDL' diagonal threshold: % .5g\n Entries with abs. value"
499 " less than dbound are replaced with +/- dbound.\n",
500 Common->dbound) ;
501
502 P4 (" grow0: memory reallocation: % .5g\n", Common->grow0) ;
503 P4 (" grow1: memory reallocation: % .5g\n", Common->grow1) ;
504 P4 (" grow2: memory reallocation: %g\n", (double) (Common->grow2)) ;
505
506 P4 ("%s", " nrelax, zrelax: supernodal amalgamation rule:\n") ;
507 P4 ("%s", " s = # columns in two adjacent supernodes\n") ;
508 P4 ("%s", " z = % of zeros in new supernode if they are merged.\n") ;
509 P4 ("%s", " Two supernodes are merged if") ;
510 P4 (" (s <= %g) or (no new zero entries) or\n",
511 (double) (Common->nrelax [0])) ;
512 P4 (" (s <= %g and ", (double) (Common->nrelax [1])) ;
513 P4 ("z < %.5g%%) or", Common->zrelax [0] * 100) ;
514 P4 (" (s <= %g and ", (double) (Common->nrelax [2])) ;
515 P4 ("z < %.5g%%) or", Common->zrelax [1] * 100) ;
516 P4 (" (z < %.5g%%)\n", Common->zrelax [2] * 100) ;
517
518 /* ---------------------------------------------------------------------- */
519 /* check workspace */
520 /* ---------------------------------------------------------------------- */
521
522 mark = Common->mark ;
523 nrow = Common->nrow ;
524 Flag = Common->Flag ;
525 Head = Common->Head ;
526 if (nrow > 0)
527 {
528 if (mark < 0 || Flag == NULL || Head == NULL)
529 {
530 ERR ("workspace corrupted (Flag and/or Head missing)") ;
531 }
532 for (i = 0 ; i < nrow ; i++)
533 {
534 if (Flag [i] >= mark)
535 {
536 PRINT0 (("Flag ["ID"]="ID", mark = %ld\n", i, Flag [i], mark)) ;
537 ERR ("workspace corrupted (Flag)") ;
538 }
539 }
540 for (i = 0 ; i <= nrow ; i++)
541 {
542 if (Head [i] != EMPTY)
543 {
544 PRINT0 (("Head ["ID"] = "ID",\n", i, Head [i])) ;
545 ERR ("workspace corrupted (Head)") ;
546 }
547 }
548 }
549 xworksize = Common->xworksize ;
550 Xwork = Common->Xwork ;
551 if (xworksize > 0)
552 {
553 if (Xwork == NULL)
554 {
555 ERR ("workspace corrupted (Xwork missing)") ;
556 }
557 for (i = 0 ; i < xworksize ; i++)
558 {
559 if (Xwork [i] != 0.)
560 {
561 PRINT0 (("Xwork ["ID"] = %g\n", i, Xwork [i])) ;
562 ERR ("workspace corrupted (Xwork)") ;
563 }
564 }
565 }
566
567 /* workspace and parameters are valid */
568 P3 ("%s", " OK\n") ;
569 P4 ("%s", "\n") ;
570 return (TRUE) ;
571 }
572
573
CHOLMOD(check_common)574 int CHOLMOD(check_common)
575 (
576 cholmod_common *Common
577 )
578 {
579 return (check_common (0, NULL, Common)) ;
580 }
581
582
CHOLMOD(print_common)583 int CHOLMOD(print_common)
584 (
585 /* ---- input ---- */
586 const char *name, /* printed name of Common object */
587 /* --------------- */
588 cholmod_common *Common
589 )
590 {
591 Int print = (Common == NULL) ? 3 : (Common->print) ;
592 return (check_common (print, name, Common)) ;
593 }
594
595
596 /* ========================================================================== */
597 /* === cholmod_gpu_stats ==================================================== */
598 /* ========================================================================== */
599
600 /* Print CPU / GPU statistics. If the timer is not installed, the times are
601 reported as zero, but this function still works. Likewise, the function
602 still works if the GPU BLAS is not installed. */
603
CHOLMOD(gpu_stats)604 int CHOLMOD(gpu_stats)
605 (
606 cholmod_common *Common /* input */
607 )
608 {
609 double cpu_time, gpu_time ;
610 int print ;
611
612 RETURN_IF_NULL_COMMON (FALSE) ;
613 print = Common->print ;
614
615 P2 ("%s", "\nCHOLMOD GPU/CPU statistics:\n") ;
616 P2 ("SYRK CPU calls %12.0f", (double) Common->CHOLMOD_CPU_SYRK_CALLS) ;
617 P2 (" time %12.4e\n", Common->CHOLMOD_CPU_SYRK_TIME) ;
618 P2 (" GPU calls %12.0f", (double) Common->CHOLMOD_GPU_SYRK_CALLS) ;
619 P2 (" time %12.4e\n", Common->CHOLMOD_GPU_SYRK_TIME) ;
620 P2 ("GEMM CPU calls %12.0f", (double) Common->CHOLMOD_CPU_GEMM_CALLS) ;
621 P2 (" time %12.4e\n", Common->CHOLMOD_CPU_GEMM_TIME) ;
622 P2 (" GPU calls %12.0f", (double) Common->CHOLMOD_GPU_GEMM_CALLS) ;
623 P2 (" time %12.4e\n", Common->CHOLMOD_GPU_GEMM_TIME) ;
624 P2 ("POTRF CPU calls %12.0f", (double) Common->CHOLMOD_CPU_POTRF_CALLS) ;
625 P2 (" time %12.4e\n", Common->CHOLMOD_CPU_POTRF_TIME) ;
626 P2 (" GPU calls %12.0f", (double) Common->CHOLMOD_GPU_POTRF_CALLS) ;
627 P2 (" time %12.4e\n", Common->CHOLMOD_GPU_POTRF_TIME) ;
628 P2 ("TRSM CPU calls %12.0f", (double) Common->CHOLMOD_CPU_TRSM_CALLS) ;
629 P2 (" time %12.4e\n", Common->CHOLMOD_CPU_TRSM_TIME) ;
630 P2 (" GPU calls %12.0f", (double) Common->CHOLMOD_GPU_TRSM_CALLS) ;
631 P2 (" time %12.4e\n", Common->CHOLMOD_GPU_TRSM_TIME) ;
632
633 cpu_time = Common->CHOLMOD_CPU_SYRK_TIME + Common->CHOLMOD_CPU_TRSM_TIME +
634 Common->CHOLMOD_CPU_GEMM_TIME + Common->CHOLMOD_CPU_POTRF_TIME ;
635
636 gpu_time = Common->CHOLMOD_GPU_SYRK_TIME + Common->CHOLMOD_GPU_TRSM_TIME +
637 Common->CHOLMOD_GPU_GEMM_TIME + Common->CHOLMOD_GPU_POTRF_TIME ;
638
639 P2 ("time in the BLAS: CPU %12.4e", cpu_time) ;
640 P2 (" GPU %12.4e", gpu_time) ;
641 P2 (" total: %12.4e\n", cpu_time + gpu_time) ;
642
643 P2 ("assembly time %12.4e", Common->CHOLMOD_ASSEMBLE_TIME) ;
644 P2 (" %12.4e\n", Common->CHOLMOD_ASSEMBLE_TIME2) ;
645 return (TRUE) ;
646 }
647
648
649 /* ========================================================================== */
650 /* === cholmod_check_sparse ================================================= */
651 /* ========================================================================== */
652
653 /* Ensure that a sparse matrix in column-oriented form is valid, and optionally
654 * print it. Returns the number of entries on the diagonal or -1 if error.
655 *
656 * workspace: Iwork (nrow)
657 */
658
check_sparse(Int * Wi,Int print,const char * name,cholmod_sparse * A,SuiteSparse_long * nnzdiag,cholmod_common * Common)659 static SuiteSparse_long check_sparse
660 (
661 Int *Wi,
662 Int print,
663 const char *name,
664 cholmod_sparse *A,
665 SuiteSparse_long *nnzdiag,
666 cholmod_common *Common
667 )
668 {
669 double *Ax, *Az ;
670 Int *Ap, *Ai, *Anz ;
671 Int nrow, ncol, nzmax, sorted, packed, j, p, pend, i, nz, ilast,
672 init_print, dnz, count, xtype ;
673 const char *type = "sparse" ;
674
675 /* ---------------------------------------------------------------------- */
676 /* print header information */
677 /* ---------------------------------------------------------------------- */
678
679 P4 ("%s", "\n") ;
680 P3 ("%s", "CHOLMOD sparse: ") ;
681 if (name != NULL)
682 {
683 P3 ("%s: ", name) ;
684 }
685
686 if (A == NULL)
687 {
688 ERR ("null") ;
689 }
690
691 nrow = A->nrow ;
692 ncol = A->ncol ;
693 nzmax = A->nzmax ;
694 sorted = A->sorted ;
695 packed = A->packed ;
696 xtype = A->xtype ;
697 Ap = A->p ;
698 Ai = A->i ;
699 Ax = A->x ;
700 Az = A->z ;
701 Anz = A->nz ;
702 nz = CHOLMOD(nnz) (A, Common) ;
703
704 P3 (" "ID"", nrow) ;
705 P3 ("-by-"ID", ", ncol) ;
706 P3 ("nz "ID",", nz) ;
707 if (A->stype > 0)
708 {
709 P3 ("%s", " upper.") ;
710 }
711 else if (A->stype < 0)
712 {
713 P3 ("%s", " lower.") ;
714 }
715 else
716 {
717 P3 ("%s", " up/lo.") ;
718 }
719
720 P4 ("\n nzmax "ID", ", nzmax) ;
721 if (nz > nzmax)
722 {
723 ERR ("nzmax too small") ;
724 }
725 if (!sorted)
726 {
727 P4 ("%s", "un") ;
728 }
729 P4 ("%s", "sorted, ") ;
730 if (!packed)
731 {
732 P4 ("%s", "un") ;
733 }
734 P4 ("%s", "packed, ") ;
735
736 switch (A->itype)
737 {
738 case CHOLMOD_INT: P4 ("%s", "\n scalar types: int, ") ; break ;
739 case CHOLMOD_INTLONG: ERR ("mixed int/long type unsupported") ;
740 case CHOLMOD_LONG: P4 ("%s", "\n scalar types: SuiteSparse_long, ");
741 break ;
742 default: ERR ("unknown itype") ;
743 }
744
745 switch (A->xtype)
746 {
747 case CHOLMOD_PATTERN: P4 ("%s", "pattern") ; break ;
748 case CHOLMOD_REAL: P4 ("%s", "real") ; break ;
749 case CHOLMOD_COMPLEX: P4 ("%s", "complex") ; break ;
750 case CHOLMOD_ZOMPLEX: P4 ("%s", "zomplex") ; break ;
751 default: ERR ("unknown xtype") ;
752 }
753
754 switch (A->dtype)
755 {
756 case CHOLMOD_DOUBLE: P4 ("%s", ", double\n") ; break ;
757 case CHOLMOD_SINGLE: ERR ("float unsupported") ;
758 default: ERR ("unknown dtype") ;
759 }
760
761 if (A->itype != ITYPE || A->dtype != DTYPE)
762 {
763 ERR ("integer and real type must match routine") ;
764 }
765
766 if (A->stype && nrow != ncol)
767 {
768 ERR ("symmetric but not square") ;
769 }
770
771 /* check for existence of Ap, Ai, Anz, Ax, and Az arrays */
772 if (Ap == NULL)
773 {
774 ERR ("p array not present") ;
775 }
776 if (Ai == NULL)
777 {
778 ERR ("i array not present") ;
779 }
780 if (!packed && Anz == NULL)
781 {
782 ERR ("nz array not present") ;
783 }
784 if (xtype != CHOLMOD_PATTERN && Ax == NULL)
785 {
786 ERR ("x array not present") ;
787 }
788 if (xtype == CHOLMOD_ZOMPLEX && Az == NULL)
789 {
790 ERR ("z array not present") ;
791 }
792
793 /* packed matrices must start at Ap [0] = 0 */
794 if (packed && Ap [0] != 0)
795 {
796 ERR ("p [0] must be zero") ;
797 }
798 if (packed && (Ap [ncol] < Ap [0] || Ap [ncol] > nzmax))
799 {
800 ERR ("p [ncol] invalid") ;
801 }
802
803 /* ---------------------------------------------------------------------- */
804 /* allocate workspace if needed */
805 /* ---------------------------------------------------------------------- */
806
807 if (!sorted)
808 {
809 if (Wi == NULL)
810 {
811 CHOLMOD(allocate_work) (0, nrow, 0, Common) ;
812 Wi = Common->Iwork ; /* size nrow, (i/i/l) */
813 }
814 if (Common->status < CHOLMOD_OK)
815 {
816 return (FALSE) ; /* out of memory */
817 }
818 for (i = 0 ; i < nrow ; i++)
819 {
820 Wi [i] = EMPTY ;
821 }
822 }
823
824 /* ---------------------------------------------------------------------- */
825 /* check and print each column */
826 /* ---------------------------------------------------------------------- */
827
828 init_print = print ;
829 dnz = 0 ;
830 ETC_START (count, 8) ;
831
832 for (j = 0 ; j < ncol ; j++)
833 {
834 ETC (j == ncol-1, count, 4) ;
835 p = Ap [j] ;
836 if (packed)
837 {
838 pend = Ap [j+1] ;
839 nz = pend - p ;
840 }
841 else
842 {
843 /* Note that Anz [j] < 0 is treated as zero */
844 nz = MAX (0, Anz [j]) ;
845 pend = p + nz ;
846 }
847 P4 (" col "ID":", j) ;
848 P4 (" nz "ID"", nz) ;
849 P4 (" start "ID"", p) ;
850 P4 (" end "ID"", pend) ;
851 P4 ("%s", ":\n") ;
852 if (p < 0 || pend > nzmax)
853 {
854 ERR ("pointer invalid") ;
855 }
856 if (nz < 0 || nz > nrow)
857 {
858 ERR ("nz invalid") ;
859 }
860 ilast = EMPTY ;
861
862 for ( ; p < pend ; p++)
863 {
864 ETC (j == ncol-1 && p >= pend-4, count, -1) ;
865 i = Ai [p] ;
866 P4 (" "I8":", i) ;
867
868 print_value (print, xtype, Ax, Az, p, Common) ;
869
870 if (i == j)
871 {
872 dnz++ ;
873 }
874 if (i < 0 || i >= nrow)
875 {
876 ERR ("row index out of range") ;
877 }
878 if (sorted && i <= ilast)
879 {
880 ERR ("row indices out of order") ;
881 }
882 if (!sorted && Wi [i] == j)
883 {
884 ERR ("duplicate row index") ;
885 }
886 P4 ("%s", "\n") ;
887 ilast = i ;
888 if (!sorted)
889 {
890 Wi [i] = j ;
891 }
892 }
893 }
894
895 /* matrix is valid */
896 P4 (" nnz on diagonal: "ID"\n", dnz) ;
897 P3 ("%s", " OK\n") ;
898 P4 ("%s", "\n") ;
899 *nnzdiag = dnz ;
900 return (TRUE) ;
901 }
902
903
CHOLMOD(check_sparse)904 int CHOLMOD(check_sparse)
905 (
906 /* ---- input ---- */
907 cholmod_sparse *A, /* sparse matrix to check */
908 /* --------------- */
909 cholmod_common *Common
910 )
911 {
912 SuiteSparse_long nnzdiag ;
913 RETURN_IF_NULL_COMMON (FALSE) ;
914 Common->status = CHOLMOD_OK ;
915 return (check_sparse (NULL, 0, NULL, A, &nnzdiag, Common)) ;
916 }
917
918
CHOLMOD(print_sparse)919 int CHOLMOD(print_sparse)
920 (
921 /* ---- input ---- */
922 cholmod_sparse *A, /* sparse matrix to print */
923 const char *name, /* printed name of sparse matrix */
924 /* --------------- */
925 cholmod_common *Common
926 )
927 {
928 SuiteSparse_long nnzdiag ;
929 RETURN_IF_NULL_COMMON (FALSE) ;
930 Common->status = CHOLMOD_OK ;
931 return (check_sparse (NULL, Common->print, name, A, &nnzdiag, Common)) ;
932 }
933
934
935 /* ========================================================================== */
936 /* === cholmod_check_dense ================================================== */
937 /* ========================================================================== */
938
939 /* Ensure a dense matrix is valid, and optionally print it. */
940
check_dense(Int print,const char * name,cholmod_dense * X,cholmod_common * Common)941 static int check_dense
942 (
943 Int print,
944 const char *name,
945 cholmod_dense *X,
946 cholmod_common *Common
947 )
948 {
949 double *Xx, *Xz ;
950 Int i, j, d, nrow, ncol, nzmax, nz, init_print, count, xtype ;
951 const char *type = "dense" ;
952
953 /* ---------------------------------------------------------------------- */
954 /* print header information */
955 /* ---------------------------------------------------------------------- */
956
957 P4 ("%s", "\n") ;
958 P3 ("%s", "CHOLMOD dense: ") ;
959 if (name != NULL)
960 {
961 P3 ("%s: ", name) ;
962 }
963
964 if (X == NULL)
965 {
966 ERR ("null") ;
967 }
968
969 nrow = X->nrow ;
970 ncol = X->ncol ;
971 nzmax = X->nzmax ;
972 d = X->d ;
973 Xx = X->x ;
974 Xz = X->z ;
975 xtype = X->xtype ;
976
977 P3 (" "ID"", nrow) ;
978 P3 ("-by-"ID", ", ncol) ;
979 P4 ("\n leading dimension "ID", ", d) ;
980 P4 ("nzmax "ID", ", nzmax) ;
981 if (d * ncol > nzmax)
982 {
983 ERR ("nzmax too small") ;
984 }
985 if (d < nrow)
986 {
987 ERR ("leading dimension must be >= # of rows") ;
988 }
989 if (Xx == NULL)
990 {
991 ERR ("null") ;
992 }
993
994 switch (X->xtype)
995 {
996 case CHOLMOD_PATTERN: ERR ("pattern unsupported") ; break ;
997 case CHOLMOD_REAL: P4 ("%s", "real") ; break ;
998 case CHOLMOD_COMPLEX: P4 ("%s", "complex") ; break ;
999 case CHOLMOD_ZOMPLEX: P4 ("%s", "zomplex") ; break ;
1000 default: ERR ("unknown xtype") ;
1001 }
1002
1003 switch (X->dtype)
1004 {
1005 case CHOLMOD_DOUBLE: P4 ("%s", ", double\n") ; break ;
1006 case CHOLMOD_SINGLE: ERR ("single unsupported") ;
1007 default: ERR ("unknown dtype") ;
1008 }
1009
1010 /* ---------------------------------------------------------------------- */
1011 /* check and print each entry */
1012 /* ---------------------------------------------------------------------- */
1013
1014 if (print >= 4)
1015 {
1016 init_print = print ;
1017 ETC_START (count, 9) ;
1018 nz = nrow * ncol ;
1019 for (j = 0 ; j < ncol ; j++)
1020 {
1021 ETC (j == ncol-1, count, 5) ;
1022 P4 (" col "ID":\n", j) ;
1023 for (i = 0 ; i < nrow ; i++)
1024 {
1025 ETC (j == ncol-1 && i >= nrow-4, count, -1) ;
1026 P4 (" "I8":", i) ;
1027
1028 print_value (print, xtype, Xx, Xz, i+j*d, Common) ;
1029
1030 P4 ("%s", "\n") ;
1031 }
1032 }
1033 }
1034
1035 /* dense is valid */
1036 P3 ("%s", " OK\n") ;
1037 P4 ("%s", "\n") ;
1038 return (TRUE) ;
1039 }
1040
1041
CHOLMOD(check_dense)1042 int CHOLMOD(check_dense)
1043 (
1044 /* ---- input ---- */
1045 cholmod_dense *X, /* dense matrix to check */
1046 /* --------------- */
1047 cholmod_common *Common
1048 )
1049 {
1050 RETURN_IF_NULL_COMMON (FALSE) ;
1051 Common->status = CHOLMOD_OK ;
1052 return (check_dense (0, NULL, X, Common)) ;
1053 }
1054
1055
CHOLMOD(print_dense)1056 int CHOLMOD(print_dense)
1057 (
1058 /* ---- input ---- */
1059 cholmod_dense *X, /* dense matrix to print */
1060 const char *name, /* printed name of dense matrix */
1061 /* --------------- */
1062 cholmod_common *Common
1063 )
1064 {
1065 RETURN_IF_NULL_COMMON (FALSE) ;
1066 Common->status = CHOLMOD_OK ;
1067 return (check_dense (Common->print, name, X, Common)) ;
1068 }
1069
1070
1071 /* ========================================================================== */
1072 /* === cholmod_check_subset ================================================= */
1073 /* ========================================================================== */
1074
1075 /* Ensure S (0:len-1) is a subset of 0:n-1. Duplicates are allowed. S may be
1076 * NULL. A negative len denotes the set 0:n-1.
1077 *
1078 * To check the rset and cset for A(rset,cset), where nc and nr are the length
1079 * of cset and rset respectively:
1080 *
1081 * cholmod_check_subset (cset, nc, A->ncol, Common) ;
1082 * cholmod_check_subset (rset, nr, A->nrow, Common) ;
1083 *
1084 * workspace: none
1085 */
1086
check_subset(Int * S,SuiteSparse_long len,size_t n,Int print,const char * name,cholmod_common * Common)1087 static int check_subset
1088 (
1089 Int *S,
1090 SuiteSparse_long len,
1091 size_t n,
1092 Int print,
1093 const char *name,
1094 cholmod_common *Common
1095 )
1096 {
1097 Int i, k, init_print, count ;
1098 const char *type = "subset" ;
1099
1100 init_print = print ;
1101
1102 if (S == NULL)
1103 {
1104 /* zero len denotes S = [ ], negative len denotes S = 0:n-1 */
1105 len = (len < 0) ? (-1) : 0 ;
1106 }
1107
1108 P4 ("%s", "\n") ;
1109 P3 ("%s", "CHOLMOD subset: ") ;
1110 if (name != NULL)
1111 {
1112 P3 ("%s: ", name) ;
1113 }
1114
1115 P3 (" len: %ld ", len) ;
1116 if (len < 0)
1117 {
1118 P3 ("%s", "(denotes 0:n-1) ") ;
1119 }
1120 P3 ("n: "ID"", (Int) n) ;
1121 P4 ("%s", "\n") ;
1122
1123 if (len <= 0 || S == NULL)
1124 {
1125 P3 ("%s", " OK\n") ;
1126 P4 ("%s", "\n") ;
1127 return (TRUE) ;
1128 }
1129
1130 if (print >= 4)
1131 {
1132 ETC_START (count, 8) ;
1133 for (k = 0 ; k < ((Int) len) ; k++)
1134 {
1135 ETC (k == ((Int) len) - 4, count, -1) ;
1136 i = S [k] ;
1137 P4 (" "I8":", k) ;
1138 P4 (" "ID"\n", i) ;
1139 if (i < 0 || i >= ((Int) n))
1140 {
1141 ERR ("entry out range") ;
1142 }
1143 }
1144 }
1145 else
1146 {
1147 for (k = 0 ; k < ((Int) len) ; k++)
1148 {
1149 i = S [k] ;
1150 if (i < 0 || i >= ((Int) n))
1151 {
1152 ERR ("entry out range") ;
1153 }
1154 }
1155 }
1156 P3 ("%s", " OK\n") ;
1157 P4 ("%s", "\n") ;
1158 return (TRUE) ;
1159 }
1160
1161
CHOLMOD(check_subset)1162 int CHOLMOD(check_subset)
1163 (
1164 /* ---- input ---- */
1165 Int *Set, /* Set [0:len-1] is a subset of 0:n-1. Duplicates OK */
1166 SuiteSparse_long len, /* size of Set (an integer array), or < 0 if 0:n-1 */
1167 size_t n, /* 0:n-1 is valid range */
1168 /* --------------- */
1169 cholmod_common *Common
1170 )
1171 {
1172 RETURN_IF_NULL_COMMON (FALSE) ;
1173 Common->status = CHOLMOD_OK ;
1174 return (check_subset (Set, len, n, 0, NULL, Common)) ;
1175 }
1176
1177
CHOLMOD(print_subset)1178 int CHOLMOD(print_subset)
1179 (
1180 /* ---- input ---- */
1181 Int *Set, /* Set [0:len-1] is a subset of 0:n-1. Duplicates OK */
1182 SuiteSparse_long len, /* size of Set (an integer array), or < 0 if 0:n-1 */
1183 size_t n, /* 0:n-1 is valid range */
1184 const char *name, /* printed name of Set */
1185 /* --------------- */
1186 cholmod_common *Common
1187 )
1188 {
1189 RETURN_IF_NULL_COMMON (FALSE) ;
1190 Common->status = CHOLMOD_OK ;
1191 return (check_subset (Set, len, n, Common->print, name, Common)) ;
1192 }
1193
1194
1195 /* ========================================================================== */
1196 /* === cholmod_check_perm =================================================== */
1197 /* ========================================================================== */
1198
1199 /* Ensure that Perm [0..len-1] is a permutation of a subset of 0:n-1. Perm
1200 * may be NULL, which is interpreted as the identity permutation. There can
1201 * be no duplicate entries (len must be <= n).
1202 *
1203 * If n <= Common->nrow, then this routine takes O(len) time and does not
1204 * allocate any memory, by using Common->Flag. Otherwise, it takes O(n) time
1205 * and ensures that Common->Iwork is at least n*sizeof(Int) in size.
1206 *
1207 * To check the fset: cholmod_check_perm (fset, fsize, ncol, Common) ;
1208 * To check a permutation: cholmod_check_perm (Perm, n, n, Common) ;
1209 *
1210 * workspace: Flag (n) if n <= Common->nrow, Iwork (n) otherwise.
1211 */
1212
check_perm(Int * Wi,Int print,const char * name,Int * Perm,size_t len,size_t n,cholmod_common * Common)1213 static int check_perm
1214 (
1215 Int *Wi,
1216 Int print,
1217 const char *name,
1218 Int *Perm,
1219 size_t len,
1220 size_t n,
1221 cholmod_common *Common
1222 )
1223 {
1224 Int *Flag ;
1225 Int i, k, mark, init_print, count ;
1226 const char *type = "perm" ;
1227
1228 /* ---------------------------------------------------------------------- */
1229 /* checks that take O(1) time */
1230 /* ---------------------------------------------------------------------- */
1231
1232 if (Perm == NULL || n == 0)
1233 {
1234 /* Perm is valid implicit identity, or empty */
1235 return (TRUE) ;
1236 }
1237
1238 /* ---------------------------------------------------------------------- */
1239 /* checks that take O(n) time or require memory allocation */
1240 /* ---------------------------------------------------------------------- */
1241
1242 init_print = print ;
1243 ETC_START (count, 8) ;
1244
1245 if (Wi == NULL && n <= Common->nrow)
1246 {
1247 /* use the Common->Flag array if it's big enough */
1248 mark = CHOLMOD(clear_flag) (Common) ;
1249 Flag = Common->Flag ;
1250 ASSERT (CHOLMOD(dump_work) (TRUE, FALSE, 0, Common)) ;
1251 if (print >= 4)
1252 {
1253 for (k = 0 ; k < ((Int) len) ; k++)
1254 {
1255 ETC (k >= ((Int) len) - 4, count, -1) ;
1256 i = Perm [k] ;
1257 P4 (" "I8":", k) ;
1258 P4 (""ID"\n", i) ;
1259 if (i < 0 || i >= ((Int) n) || Flag [i] == mark)
1260 {
1261 CHOLMOD(clear_flag) (Common) ;
1262 ERR ("invalid permutation") ;
1263 }
1264 Flag [i] = mark ;
1265 }
1266 }
1267 else
1268 {
1269 for (k = 0 ; k < ((Int) len) ; k++)
1270 {
1271 i = Perm [k] ;
1272 if (i < 0 || i >= ((Int) n) || Flag [i] == mark)
1273 {
1274 CHOLMOD(clear_flag) (Common) ;
1275 ERR ("invalid permutation") ;
1276 }
1277 Flag [i] = mark ;
1278 }
1279 }
1280 CHOLMOD(clear_flag) (Common) ;
1281 ASSERT (CHOLMOD(dump_work) (TRUE, FALSE, 0, Common)) ;
1282 }
1283 else
1284 {
1285 if (Wi == NULL)
1286 {
1287 /* use Common->Iwork instead, but initialize it first */
1288 CHOLMOD(allocate_work) (0, n, 0, Common) ;
1289 Wi = Common->Iwork ; /* size n, (i/i/i) is OK */
1290 }
1291 if (Common->status < CHOLMOD_OK)
1292 {
1293 return (FALSE) ; /* out of memory */
1294 }
1295 for (i = 0 ; i < ((Int) n) ; i++)
1296 {
1297 Wi [i] = FALSE ;
1298 }
1299 if (print >= 4)
1300 {
1301 for (k = 0 ; k < ((Int) len) ; k++)
1302 {
1303 ETC (k >= ((Int) len) - 4, count, -1) ;
1304 i = Perm [k] ;
1305 P4 (" "I8":", k) ;
1306 P4 (""ID"\n", i) ;
1307 if (i < 0 || i >= ((Int) n) || Wi [i])
1308 {
1309 ERR ("invalid permutation") ;
1310 }
1311 Wi [i] = TRUE ;
1312 }
1313 }
1314 else
1315 {
1316 for (k = 0 ; k < ((Int) len) ; k++)
1317 {
1318 i = Perm [k] ;
1319 if (i < 0 || i >= ((Int) n) || Wi [i])
1320 {
1321 ERR ("invalid permutation") ;
1322 }
1323 Wi [i] = TRUE ;
1324 }
1325 }
1326 }
1327
1328 /* perm is valid */
1329 return (TRUE) ;
1330 }
1331
1332
CHOLMOD(check_perm)1333 int CHOLMOD(check_perm)
1334 (
1335 /* ---- input ---- */
1336 Int *Perm, /* Perm [0:len-1] is a permutation of subset of 0:n-1 */
1337 size_t len, /* size of Perm (an integer array) */
1338 size_t n, /* 0:n-1 is valid range */
1339 /* --------------- */
1340 cholmod_common *Common
1341 )
1342 {
1343 RETURN_IF_NULL_COMMON (FALSE) ;
1344 Common->status = CHOLMOD_OK ;
1345 return (check_perm (NULL, 0, NULL, Perm, len, n, Common)) ;
1346 }
1347
1348
CHOLMOD(print_perm)1349 int CHOLMOD(print_perm)
1350 (
1351 /* ---- input ---- */
1352 Int *Perm, /* Perm [0:len-1] is a permutation of subset of 0:n-1 */
1353 size_t len, /* size of Perm (an integer array) */
1354 size_t n, /* 0:n-1 is valid range */
1355 const char *name, /* printed name of Perm */
1356 /* --------------- */
1357 cholmod_common *Common
1358 )
1359 {
1360 Int ok, print ;
1361 RETURN_IF_NULL_COMMON (FALSE) ;
1362 Common->status = CHOLMOD_OK ;
1363 print = Common->print ;
1364 P4 ("%s", "\n") ;
1365 P3 ("%s", "CHOLMOD perm: ") ;
1366 if (name != NULL)
1367 {
1368 P3 ("%s: ", name) ;
1369 }
1370 P3 (" len: "ID"", (Int) len) ;
1371 P3 (" n: "ID"", (Int) n) ;
1372 P4 ("%s", "\n") ;
1373 ok = check_perm (NULL, print, name, Perm, len, n, Common) ;
1374 if (ok)
1375 {
1376 P3 ("%s", " OK\n") ;
1377 P4 ("%s", "\n") ;
1378 }
1379 return (ok) ;
1380 }
1381
1382
1383 /* ========================================================================== */
1384 /* === cholmod_check_parent ================================================= */
1385 /* ========================================================================== */
1386
1387 /* Ensure that Parent is a valid elimination tree of nodes 0 to n-1.
1388 * If j is a root of the tree then Parent [j] is EMPTY (-1).
1389 *
1390 * NOTE: this check will fail if applied to the component tree (CParent) in
1391 * cholmod_nested_dissection, unless it has been postordered and renumbered.
1392 *
1393 * workspace: none
1394 */
1395
check_parent(Int * Parent,size_t n,Int print,const char * name,cholmod_common * Common)1396 static int check_parent
1397 (
1398 Int *Parent,
1399 size_t n,
1400 Int print,
1401 const char *name,
1402 cholmod_common *Common
1403 )
1404 {
1405 Int j, p, init_print, count ;
1406 const char *type = "parent" ;
1407
1408 init_print = print ;
1409
1410 P4 ("%s", "\n") ;
1411 P3 ("%s", "CHOLMOD parent: ") ;
1412 if (name != NULL)
1413 {
1414 P3 ("%s: ", name) ;
1415 }
1416
1417 P3 (" n: "ID"", (Int) n) ;
1418 P4 ("%s", "\n") ;
1419
1420 if (Parent == NULL)
1421 {
1422 ERR ("null") ;
1423 }
1424
1425 /* ---------------------------------------------------------------------- */
1426 /* checks that take O(n) time */
1427 /* ---------------------------------------------------------------------- */
1428
1429 ETC_START (count, 8) ;
1430 for (j = 0 ; j < ((Int) n) ; j++)
1431 {
1432 ETC (j == ((Int) n) - 4, count, -1) ;
1433 p = Parent [j] ;
1434 P4 (" "I8":", j) ;
1435 P4 (" "ID"\n", p) ;
1436 if (!(p == EMPTY || p > j))
1437 {
1438 ERR ("invalid") ;
1439 }
1440 }
1441 P3 ("%s", " OK\n") ;
1442 P4 ("%s", "\n") ;
1443 return (TRUE) ;
1444 }
1445
1446
CHOLMOD(check_parent)1447 int CHOLMOD(check_parent)
1448 (
1449 /* ---- input ---- */
1450 Int *Parent, /* Parent [0:n-1] is an elimination tree */
1451 size_t n, /* size of Parent */
1452 /* --------------- */
1453 cholmod_common *Common
1454 )
1455 {
1456 RETURN_IF_NULL_COMMON (FALSE) ;
1457 Common->status = CHOLMOD_OK ;
1458 return (check_parent (Parent, n, 0, NULL, Common)) ;
1459 }
1460
1461
CHOLMOD(print_parent)1462 int CHOLMOD(print_parent)
1463 (
1464 /* ---- input ---- */
1465 Int *Parent, /* Parent [0:n-1] is an elimination tree */
1466 size_t n, /* size of Parent */
1467 const char *name, /* printed name of Parent */
1468 /* --------------- */
1469 cholmod_common *Common
1470 )
1471 {
1472 RETURN_IF_NULL_COMMON (FALSE) ;
1473 Common->status = CHOLMOD_OK ;
1474 return (check_parent (Parent, n, Common->print, name, Common)) ;
1475 }
1476
1477
1478 /* ========================================================================== */
1479 /* === cholmod_check_factor ================================================= */
1480 /* ========================================================================== */
1481
check_factor(Int * Wi,Int print,const char * name,cholmod_factor * L,cholmod_common * Common)1482 static int check_factor
1483 (
1484 Int *Wi,
1485 Int print,
1486 const char *name,
1487 cholmod_factor *L,
1488 cholmod_common *Common
1489 )
1490 {
1491 double *Lx, *Lz ;
1492 Int *Lp, *Li, *Lnz, *Lnext, *Lprev, *Perm, *ColCount, *Lpi, *Lpx, *Super,
1493 *Ls ;
1494 Int n, nzmax, j, p, pend, i, nz, ordering, space, is_monotonic, minor,
1495 count, precise, init_print, ilast, lnz, head, tail, jprev, plast,
1496 jnext, examine_super, nsuper, s, k1, k2, psi, psend, psx, nsrow, nscol,
1497 ps2, psxend, ssize, xsize, maxcsize, maxesize, nsrow2, jj, ii, xtype ;
1498 Int check_Lpx ;
1499 const char *type = "factor" ;
1500
1501 /* ---------------------------------------------------------------------- */
1502 /* print header information */
1503 /* ---------------------------------------------------------------------- */
1504
1505 P4 ("%s", "\n") ;
1506 P3 ("%s", "CHOLMOD factor: ") ;
1507 if (name != NULL)
1508 {
1509 P3 ("%s: ", name) ;
1510 }
1511
1512 if (L == NULL)
1513 {
1514 ERR ("null") ;
1515 }
1516
1517 n = L->n ;
1518 minor = L->minor ;
1519 ordering = L->ordering ;
1520 xtype = L->xtype ;
1521
1522 Perm = L->Perm ;
1523 ColCount = L->ColCount ;
1524 lnz = 0 ;
1525
1526 precise = Common->precise ;
1527
1528 P3 (" "ID"", n) ;
1529 P3 ("-by-"ID"", n) ;
1530
1531 if (minor < n)
1532 {
1533 P3 (" not positive definite (column "ID")", minor) ;
1534 }
1535
1536 switch (L->itype)
1537 {
1538 case CHOLMOD_INT: P4 ("%s", "\n scalar types: int, ") ; break ;
1539 case CHOLMOD_INTLONG: ERR ("mixed int/long type unsupported") ;
1540 case CHOLMOD_LONG: P4 ("%s", "\n scalar types: SuiteSparse_long, ");
1541 break ;
1542 default: ERR ("unknown itype") ;
1543 }
1544
1545 switch (L->xtype)
1546 {
1547 case CHOLMOD_PATTERN: P4 ("%s", "pattern") ; break ;
1548 case CHOLMOD_REAL: P4 ("%s", "real") ; break ;
1549 case CHOLMOD_COMPLEX: P4 ("%s", "complex") ; break ;
1550 case CHOLMOD_ZOMPLEX: P4 ("%s", "zomplex") ; break ;
1551 default: ERR ("unknown xtype") ;
1552 }
1553
1554 switch (L->dtype)
1555 {
1556 case CHOLMOD_DOUBLE: P4 ("%s", ", double\n") ; break ;
1557 case CHOLMOD_SINGLE: ERR ("single unsupported") ;
1558 default: ERR ("unknown dtype") ;
1559 }
1560
1561 if (L->itype != ITYPE || L->dtype != DTYPE)
1562 {
1563 ERR ("integer and real type must match routine") ;
1564 }
1565
1566 if (L->is_super)
1567 {
1568 P3 ("%s", " supernodal") ;
1569 }
1570 else
1571 {
1572 P3 ("%s", " simplicial") ;
1573 }
1574
1575 if (L->is_ll)
1576 {
1577 P3 ("%s", ", LL'.") ;
1578 }
1579 else
1580 {
1581 P3 ("%s", ", LDL'.") ;
1582 }
1583
1584 P4 ("%s", "\n ordering method used: ") ;
1585 switch (L->ordering)
1586 {
1587 case CHOLMOD_POSTORDERED:P4 ("%s", "natural (postordered)") ; break ;
1588 case CHOLMOD_NATURAL: P4 ("%s", "natural") ; break ;
1589 case CHOLMOD_GIVEN: P4 ("%s", "user-provided") ; break ;
1590 case CHOLMOD_AMD: P4 ("%s", "AMD") ; break ;
1591 case CHOLMOD_COLAMD: P4 ("%s", "AMD for A, COLAMD for A*A'") ;break ;
1592 #ifndef NPARTITION
1593 case CHOLMOD_METIS: P4 ("%s", "METIS NodeND") ; break ;
1594 case CHOLMOD_NESDIS: P4 ("%s", "CHOLMOD nested dissection") ; break ;
1595 #endif
1596 default: ERR ("unknown ordering") ;
1597 }
1598
1599 P4 ("%s", "\n") ;
1600
1601 init_print = print ;
1602
1603 if (L->is_super && L->xtype == CHOLMOD_ZOMPLEX)
1604 {
1605 ERR ("Supernodal zomplex L not supported") ;
1606 }
1607
1608 /* ---------------------------------------------------------------------- */
1609 /* check L->Perm */
1610 /* ---------------------------------------------------------------------- */
1611
1612 if (!check_perm (Wi, print, name, Perm, n, n, Common))
1613 {
1614 return (FALSE) ;
1615 }
1616
1617 /* ---------------------------------------------------------------------- */
1618 /* check L->ColCount */
1619 /* ---------------------------------------------------------------------- */
1620
1621 if (ColCount == NULL)
1622 {
1623 ERR ("ColCount vector invalid") ;
1624 }
1625
1626 ETC_START (count, 8) ;
1627 for (j = 0 ; j < n ; j++)
1628 {
1629 ETC (j >= n-4, count, -1) ;
1630 P4 (" col: "ID" ", j) ;
1631 nz = ColCount [j] ;
1632 P4 ("colcount: "ID"\n", nz) ;
1633 if (nz < 0 || nz > n-j)
1634 {
1635 ERR ("ColCount out of range") ;
1636 }
1637 }
1638
1639 /* ---------------------------------------------------------------------- */
1640 /* check factor */
1641 /* ---------------------------------------------------------------------- */
1642
1643 if (L->xtype == CHOLMOD_PATTERN && !(L->is_super))
1644 {
1645
1646 /* ------------------------------------------------------------------ */
1647 /* check simplicial symbolic factor */
1648 /* ------------------------------------------------------------------ */
1649
1650 /* nothing else to do */ ;
1651
1652 }
1653 else if (L->xtype != CHOLMOD_PATTERN && !(L->is_super))
1654 {
1655
1656 /* ------------------------------------------------------------------ */
1657 /* check simplicial numerical factor */
1658 /* ------------------------------------------------------------------ */
1659
1660 P4 ("monotonic: %d\n", L->is_monotonic) ;
1661 nzmax = L->nzmax ;
1662 P3 (" nzmax "ID".", nzmax) ;
1663 P4 ("%s", "\n") ;
1664 Lp = L->p ;
1665 Li = L->i ;
1666 Lx = L->x ;
1667 Lz = L->z ;
1668 Lnz = L->nz ;
1669 Lnext = L->next ;
1670 Lprev = L->prev ;
1671
1672 /* check for existence of Lp, Li, Lnz, Lnext, Lprev, and Lx arrays */
1673 if (Lp == NULL)
1674 {
1675 ERR ("p array not present") ;
1676 }
1677 if (Li == NULL)
1678 {
1679 ERR ("i array not present") ;
1680 }
1681 if (Lnz == NULL)
1682 {
1683 ERR ("nz array not present") ;
1684 }
1685 if (Lx == NULL)
1686 {
1687 ERR ("x array not present") ;
1688 }
1689 if (xtype == CHOLMOD_ZOMPLEX && Lz == NULL)
1690 {
1691 ERR ("z array not present") ;
1692 }
1693 if (Lnext == NULL)
1694 {
1695 ERR ("next array not present") ;
1696 }
1697 if (Lprev == NULL)
1698 {
1699 ERR ("prev array not present") ;
1700 }
1701
1702 ETC_START (count, 8) ;
1703
1704 /* check each column of L */
1705 plast = 0 ;
1706 is_monotonic = TRUE ;
1707 for (j = 0 ; j < n ; j++)
1708 {
1709 ETC (j >= n-3, count, -1) ;
1710 p = Lp [j] ;
1711 nz = Lnz [j] ;
1712 pend = p + nz ;
1713 lnz += nz ;
1714
1715 P4 (" col "ID":", j) ;
1716 P4 (" nz "ID"", nz) ;
1717 P4 (" start "ID"", p) ;
1718 P4 (" end "ID"", pend) ;
1719
1720 if (Lnext [j] < 0 || Lnext [j] > n)
1721 {
1722 ERR ("invalid link list") ;
1723 }
1724 space = Lp [Lnext [j]] - p ;
1725
1726 P4 (" space "ID"", space) ;
1727 P4 (" free "ID":\n", space - nz) ;
1728
1729 if (p < 0 || pend > nzmax || space < 1)
1730 {
1731 ERR ("pointer invalid") ;
1732 }
1733 if (nz < 1 || nz > (n-j) || nz > space)
1734 {
1735 ERR ("nz invalid") ;
1736 }
1737 ilast = j-1 ;
1738
1739 if (p < plast)
1740 {
1741 is_monotonic = FALSE ;
1742 }
1743 plast = p ;
1744
1745 i = Li [p] ;
1746 P4 (" "I8":", i) ;
1747 if (i != j)
1748 {
1749 ERR ("diagonal missing") ;
1750 }
1751
1752 print_value (print, xtype, Lx, Lz, p, Common) ;
1753
1754 P4 ("%s", "\n") ;
1755 ilast = j ;
1756 for (p++ ; p < pend ; p++)
1757 {
1758 ETC_DISABLE (count) ;
1759 i = Li [p] ;
1760 P4 (" "I8":", i) ;
1761 if (i < j || i >= n)
1762 {
1763 ERR ("row index out of range") ;
1764 }
1765 if (i <= ilast)
1766 {
1767 ERR ("row indices out of order") ;
1768 }
1769
1770 print_value (print, xtype, Lx, Lz, p, Common) ;
1771
1772 P4 ("%s", "\n") ;
1773 ilast = i ;
1774 }
1775 }
1776
1777 if (L->is_monotonic && !is_monotonic)
1778 {
1779 ERR ("columns not monotonic") ;
1780 }
1781
1782 /* check the link list */
1783 head = n+1 ;
1784 tail = n ;
1785 j = head ;
1786 jprev = EMPTY ;
1787 count = 0 ;
1788 for ( ; ; )
1789 {
1790 if (j < 0 || j > n+1 || count > n+2)
1791 {
1792 ERR ("invalid link list") ;
1793 }
1794 jnext = Lnext [j] ;
1795 if (j >= 0 && j < n)
1796 {
1797 if (jprev != Lprev [j])
1798 {
1799 ERR ("invalid link list") ;
1800 }
1801 }
1802 count++ ;
1803 if (j == tail)
1804 {
1805 break ;
1806 }
1807 jprev = j ;
1808 j = jnext ;
1809 }
1810 if (Lnext [tail] != EMPTY || count != n+2)
1811 {
1812 ERR ("invalid link list") ;
1813 }
1814
1815 }
1816 else
1817 {
1818
1819 /* ------------------------------------------------------------------ */
1820 /* check supernodal numeric or symbolic factor */
1821 /* ------------------------------------------------------------------ */
1822
1823 nsuper = L->nsuper ;
1824 ssize = L->ssize ;
1825 xsize = L->xsize ;
1826 maxcsize = L->maxcsize ;
1827 maxesize = L->maxesize ;
1828 Ls = L->s ;
1829 Lpi = L->pi ;
1830 Lpx = L->px ;
1831 Super = L->super ;
1832 Lx = L->x ;
1833 ETC_START (count, 8) ;
1834
1835 P4 (" ssize "ID" ", ssize) ;
1836 P4 ("xsize "ID" ", xsize) ;
1837 P4 ("maxcsize "ID" ", maxcsize) ;
1838 P4 ("maxesize "ID"\n", maxesize) ;
1839
1840 if (Ls == NULL)
1841 {
1842 ERR ("invalid: L->s missing") ;
1843 }
1844 if (Lpi == NULL)
1845 {
1846 ERR ("invalid: L->pi missing") ;
1847 }
1848 if (Lpx == NULL)
1849 {
1850 ERR ("invalid: L->px missing") ;
1851 }
1852 if (Super == NULL)
1853 {
1854 ERR ("invalid: L->super missing") ;
1855 }
1856
1857 if (L->xtype != CHOLMOD_PATTERN)
1858 {
1859 /* numerical supernodal factor */
1860 if (Lx == NULL)
1861 {
1862 ERR ("invalid: L->x missing") ;
1863 }
1864 if (Ls [0] == EMPTY)
1865 {
1866 ERR ("invalid: L->s not defined") ;
1867 }
1868 examine_super = TRUE ;
1869 }
1870 else
1871 {
1872 /* symbolic supernodal factor, but only if it has been computed */
1873 examine_super = (Ls [0] != EMPTY) ;
1874 }
1875
1876 if (examine_super)
1877 {
1878 if (Lpi [0] != 0 || MAX (1, Lpi [nsuper]) != ssize)
1879 {
1880 PRINT0 (("Lpi [0] "ID", Lpi [nsuper = "ID"] = "ID"\n",
1881 Lpi [0], nsuper, Lpi [nsuper])) ;
1882 ERR ("invalid: L->pi invalid") ;
1883 }
1884
1885 /* If Lpx [0] is 123456, then supernodes are present but
1886 Lpx [0...nsuper] is not defined, so don't check it. This is
1887 used in the non-GPU accelerated SPQR */
1888 check_Lpx = (Lpx [0] != 123456) ;
1889 if (check_Lpx && (Lpx [0] != 0 || MAX (1, Lpx[nsuper]) != xsize))
1890 {
1891 ERR ("invalid: L->px invalid") ;
1892 }
1893
1894 /* check and print each supernode */
1895 for (s = 0 ; s < nsuper ; s++)
1896 {
1897 k1 = Super [s] ;
1898 k2 = Super [s+1] ;
1899 psi = Lpi [s] ;
1900 psend = Lpi [s+1] ;
1901 nsrow = psend - psi ;
1902 nscol = k2 - k1 ;
1903 nsrow2 = nsrow - nscol ;
1904 ps2 = psi + nscol ;
1905
1906 if (check_Lpx)
1907 {
1908 psx = Lpx [s] ;
1909 psxend = Lpx [s+1] ;
1910 }
1911
1912 ETC (s == nsuper-1, count, 4) ;
1913
1914 P4 (" supernode "ID", ", s) ;
1915 P4 ("col "ID" ", k1) ;
1916 P4 ("to "ID". ", k2-1) ;
1917 P4 ("nz in first col: "ID".\n", nsrow) ;
1918
1919 if (check_Lpx)
1920 {
1921 P4 (" values start "ID", ", psx) ;
1922 P4 ("end "ID"\n", psxend) ;
1923 }
1924
1925 if (k1 > k2 || k1 < 0 || k2 > n || nsrow < nscol || nsrow2 < 0
1926 || (check_Lpx && psxend - psx != nsrow * nscol))
1927 {
1928 ERR ("invalid supernode") ;
1929 }
1930
1931 lnz += nscol * nsrow - (nscol*nscol - nscol)/2 ;
1932
1933 if (L->xtype != CHOLMOD_PATTERN)
1934 {
1935 /* print each column of the supernode */
1936 for (jj = 0 ; jj < nscol ; jj++)
1937 {
1938 ETC_ENABLE (s == nsuper-1 && jj >= nscol-3, count, -1) ;
1939 j = k1 + jj ;
1940 P4 (" col "ID"\n", j) ;
1941 ilast = j ;
1942 i = Ls [psi + jj] ;
1943 P4 (" "I8":", i) ;
1944 if (i != j)
1945 {
1946 ERR ("row index invalid") ;
1947 }
1948
1949 /* PRINTVALUE (Lx [psx + jj + jj*nsrow]) ; */
1950 print_value (print, xtype, Lx, NULL,
1951 psx + jj + jj*nsrow, Common) ;
1952
1953 P4 ("%s", "\n") ;
1954 for (ii = jj + 1 ; ii < nsrow ; ii++)
1955 {
1956 ETC_DISABLE (count) ;
1957 i = Ls [psi + ii] ;
1958 P4 (" "I8":", i) ;
1959 if (i <= ilast || i > n)
1960 {
1961 ERR ("row index out of range") ;
1962 }
1963
1964 /* PRINTVALUE (Lx [psx + ii + jj*nsrow]) ; */
1965 print_value (print, xtype, Lx, NULL,
1966 psx + ii + jj*nsrow, Common) ;
1967
1968 P4 ("%s", "\n") ;
1969 ilast = i ;
1970 }
1971 }
1972 }
1973 else
1974 {
1975 /* just print the leading column of the supernode */
1976 P4 (" col "ID"\n", k1) ;
1977 for (jj = 0 ; jj < nscol ; jj++)
1978 {
1979 ETC (s == nsuper-1 && jj >= nscol-3, count, -1) ;
1980 j = k1 + jj ;
1981 i = Ls [psi + jj] ;
1982 P4 (" "I8"", i) ;
1983 if (i != j)
1984 {
1985 ERR ("row index invalid") ;
1986 }
1987 P4 ("%s", "\n") ;
1988 }
1989 ilast = j ;
1990 for (ii = nscol ; ii < nsrow ; ii++)
1991 {
1992 ETC_DISABLE (count) ;
1993 i = Ls [psi + ii] ;
1994 P4 (" "I8"", i) ;
1995 if (i <= ilast || i > n)
1996 {
1997 ERR ("row index out of range") ;
1998 }
1999 P4 ("%s", "\n") ;
2000 ilast = i ;
2001 }
2002 }
2003 }
2004 }
2005 }
2006
2007 /* factor is valid */
2008 P3 (" nz "ID"", lnz) ;
2009 P3 ("%s", " OK\n") ;
2010 P4 ("%s", "\n") ;
2011 return (TRUE) ;
2012 }
2013
2014
CHOLMOD(check_factor)2015 int CHOLMOD(check_factor)
2016 (
2017 /* ---- input ---- */
2018 cholmod_factor *L, /* factor to check */
2019 /* --------------- */
2020 cholmod_common *Common
2021 )
2022 {
2023 RETURN_IF_NULL_COMMON (FALSE) ;
2024 Common->status = CHOLMOD_OK ;
2025 return (check_factor (NULL, 0, NULL, L, Common)) ;
2026 }
2027
2028
CHOLMOD(print_factor)2029 int CHOLMOD(print_factor)
2030 (
2031 /* ---- input ---- */
2032 cholmod_factor *L, /* factor to print */
2033 const char *name, /* printed name of factor */
2034 /* --------------- */
2035 cholmod_common *Common
2036 )
2037 {
2038 RETURN_IF_NULL_COMMON (FALSE) ;
2039 Common->status = CHOLMOD_OK ;
2040 return (check_factor (NULL, Common->print, name, L, Common)) ;
2041 }
2042
2043
2044 /* ========================================================================== */
2045 /* === cholmod_check_triplet ================================================ */
2046 /* ========================================================================== */
2047
2048 /* Ensure a triplet matrix is valid, and optionally print it. */
2049
check_triplet(Int print,const char * name,cholmod_triplet * T,cholmod_common * Common)2050 static int check_triplet
2051 (
2052 Int print,
2053 const char *name,
2054 cholmod_triplet *T,
2055 cholmod_common *Common
2056 )
2057 {
2058 double *Tx, *Tz ;
2059 Int *Ti, *Tj ;
2060 Int i, j, p, nrow, ncol, nzmax, nz, xtype, init_print, count ;
2061 const char *type = "triplet" ;
2062
2063 /* ---------------------------------------------------------------------- */
2064 /* print header information */
2065 /* ---------------------------------------------------------------------- */
2066
2067 P4 ("%s", "\n") ;
2068 P3 ("%s", "CHOLMOD triplet: ") ;
2069 if (name != NULL)
2070 {
2071 P3 ("%s: ", name) ;
2072 }
2073
2074 if (T == NULL)
2075 {
2076 ERR ("null") ;
2077 }
2078
2079 nrow = T->nrow ;
2080 ncol = T->ncol ;
2081 nzmax = T->nzmax ;
2082 nz = T->nnz ;
2083 Ti = T->i ;
2084 Tj = T->j ;
2085 Tx = T->x ;
2086 Tz = T->z ;
2087 xtype = T->xtype ;
2088
2089
2090 P3 (" "ID"", nrow) ;
2091 P3 ("-by-"ID", ", ncol) ;
2092 P3 ("nz "ID",", nz) ;
2093 if (T->stype > 0)
2094 {
2095 P3 ("%s", " upper.") ;
2096 }
2097 else if (T->stype < 0)
2098 {
2099 P3 ("%s", " lower.") ;
2100 }
2101 else
2102 {
2103 P3 ("%s", " up/lo.") ;
2104 }
2105
2106 P4 ("\n nzmax "ID", ", nzmax) ;
2107 if (nz > nzmax)
2108 {
2109 ERR ("nzmax too small") ;
2110 }
2111
2112 switch (T->itype)
2113 {
2114 case CHOLMOD_INT: P4 ("%s", "\n scalar types: int, ") ; break ;
2115 case CHOLMOD_INTLONG: ERR ("mixed int/long type unsupported") ;
2116 case CHOLMOD_LONG: P4 ("%s", "\n scalar types: SuiteSparse_long, ");
2117 break ;
2118 default: ERR ("unknown itype") ;
2119 }
2120
2121 switch (T->xtype)
2122 {
2123 case CHOLMOD_PATTERN: P4 ("%s", "pattern") ; break ;
2124 case CHOLMOD_REAL: P4 ("%s", "real") ; break ;
2125 case CHOLMOD_COMPLEX: P4 ("%s", "complex") ; break ;
2126 case CHOLMOD_ZOMPLEX: P4 ("%s", "zomplex") ; break ;
2127 default: ERR ("unknown xtype") ;
2128 }
2129
2130 switch (T->dtype)
2131 {
2132 case CHOLMOD_DOUBLE: P4 ("%s", ", double\n") ; break ;
2133 case CHOLMOD_SINGLE: ERR ("single unsupported") ;
2134 default: ERR ("unknown dtype") ;
2135 }
2136
2137 if (T->itype != ITYPE || T->dtype != DTYPE)
2138 {
2139 ERR ("integer and real type must match routine") ;
2140 }
2141
2142 if (T->stype && nrow != ncol)
2143 {
2144 ERR ("symmetric but not square") ;
2145 }
2146
2147 /* check for existence of Ti, Tj, Tx arrays */
2148 if (Tj == NULL)
2149 {
2150 ERR ("j array not present") ;
2151 }
2152 if (Ti == NULL)
2153 {
2154 ERR ("i array not present") ;
2155 }
2156
2157 if (xtype != CHOLMOD_PATTERN && Tx == NULL)
2158 {
2159 ERR ("x array not present") ;
2160 }
2161 if (xtype == CHOLMOD_ZOMPLEX && Tz == NULL)
2162 {
2163 ERR ("z array not present") ;
2164 }
2165
2166 /* ---------------------------------------------------------------------- */
2167 /* check and print each entry */
2168 /* ---------------------------------------------------------------------- */
2169
2170 init_print = print ;
2171 ETC_START (count, 8) ;
2172
2173 for (p = 0 ; p < nz ; p++)
2174 {
2175 ETC (p >= nz-4, count, -1) ;
2176 i = Ti [p] ;
2177 P4 (" "I8":", p) ;
2178 P4 (" "I_8"", i) ;
2179 if (i < 0 || i >= nrow)
2180 {
2181 ERR ("row index out of range") ;
2182 }
2183 j = Tj [p] ;
2184 P4 (" "I_8"", j) ;
2185 if (j < 0 || j >= ncol)
2186 {
2187 ERR ("column index out of range") ;
2188 }
2189
2190 print_value (print, xtype, Tx, Tz, p, Common) ;
2191
2192 P4 ("%s", "\n") ;
2193 }
2194
2195 /* triplet matrix is valid */
2196 P3 ("%s", " OK\n") ;
2197 P4 ("%s", "\n") ;
2198 return (TRUE) ;
2199 }
2200
2201
2202
CHOLMOD(check_triplet)2203 int CHOLMOD(check_triplet)
2204 (
2205 /* ---- input ---- */
2206 cholmod_triplet *T, /* triplet matrix to check */
2207 /* --------------- */
2208 cholmod_common *Common
2209 )
2210 {
2211 RETURN_IF_NULL_COMMON (FALSE) ;
2212 Common->status = CHOLMOD_OK ;
2213 return (check_triplet (0, NULL, T, Common)) ;
2214 }
2215
2216
CHOLMOD(print_triplet)2217 int CHOLMOD(print_triplet)
2218 (
2219 /* ---- input ---- */
2220 cholmod_triplet *T, /* triplet matrix to print */
2221 const char *name, /* printed name of triplet matrix */
2222 /* --------------- */
2223 cholmod_common *Common
2224 )
2225 {
2226 RETURN_IF_NULL_COMMON (FALSE) ;
2227 Common->status = CHOLMOD_OK ;
2228 return (check_triplet (Common->print, name, T, Common)) ;
2229 }
2230
2231
2232
2233 /* ========================================================================== */
2234 /* === CHOLMOD debugging routines =========================================== */
2235 /* ========================================================================== */
2236
2237 #ifndef NDEBUG
2238
2239 /* The global variables present only when debugging enabled. */
2240 int CHOLMOD(dump) = 0 ;
2241 int CHOLMOD(dump_malloc) = -1 ;
2242
2243 /* workspace: no debug routines use workspace in Common */
2244
2245 /* ========================================================================== */
2246 /* === cholmod_dump_init ==================================================== */
2247 /* ========================================================================== */
2248
CHOLMOD(dump_init)2249 void CHOLMOD(dump_init) (const char *s, cholmod_common *Common)
2250 {
2251 int i = 0 ;
2252 FILE *f ;
2253 f = fopen ("debug", "r") ;
2254 CHOLMOD(dump) = 0 ;
2255 if (f != NULL)
2256 {
2257 i = fscanf (f, "%d", &CHOLMOD(dump)) ;
2258 fclose (f) ;
2259 }
2260 PRINT1 (("%s: cholmod_dump_init, D = %d\n", s, CHOLMOD(dump))) ;
2261 }
2262
2263
2264 /* ========================================================================== */
2265 /* === cholmod_dump_sparse ================================================== */
2266 /* ========================================================================== */
2267
2268 /* returns nnz (diag (A)) or EMPTY if error */
2269
CHOLMOD(dump_sparse)2270 SuiteSparse_long CHOLMOD(dump_sparse)
2271 (
2272 cholmod_sparse *A,
2273 const char *name,
2274 cholmod_common *Common
2275 )
2276 {
2277 Int *Wi ;
2278 SuiteSparse_long nnzdiag ;
2279 Int ok ;
2280
2281 if (CHOLMOD(dump) < -1)
2282 {
2283 /* no checks if debug level is -2 or less */
2284 return (0) ;
2285 }
2286
2287 RETURN_IF_NULL_COMMON (FALSE) ;
2288 RETURN_IF_NULL (A, FALSE) ;
2289 Wi = malloc (MAX (1, A->nrow) * sizeof (Int)) ;
2290 ok = check_sparse (Wi, CHOLMOD(dump), name, A, &nnzdiag, Common) ;
2291 if (Wi != NULL) free (Wi) ;
2292 return (ok ? nnzdiag : EMPTY) ;
2293 }
2294
2295
2296 /* ========================================================================== */
2297 /* === cholmod_dump_factor ================================================== */
2298 /* ========================================================================== */
2299
CHOLMOD(dump_factor)2300 int CHOLMOD(dump_factor)
2301 (
2302 cholmod_factor *L,
2303 const char *name,
2304 cholmod_common *Common
2305 )
2306 {
2307 Int *Wi ;
2308 int ok ;
2309
2310 if (CHOLMOD(dump) < -1)
2311 {
2312 /* no checks if debug level is -2 or less */
2313 return (TRUE) ;
2314 }
2315 RETURN_IF_NULL_COMMON (FALSE) ;
2316 RETURN_IF_NULL (L, FALSE) ;
2317 Wi = malloc (MAX (1, L->n) * sizeof (Int)) ;
2318 ok = check_factor (Wi, CHOLMOD(dump), name, L, Common) ;
2319 if (Wi != NULL) free (Wi) ;
2320 return (ok) ;
2321 }
2322
2323
2324 /* ========================================================================== */
2325 /* === cholmod_dump_perm ==================================================== */
2326 /* ========================================================================== */
2327
CHOLMOD(dump_perm)2328 int CHOLMOD(dump_perm)
2329 (
2330 Int *Perm,
2331 size_t len,
2332 size_t n,
2333 const char *name,
2334 cholmod_common *Common
2335 )
2336 {
2337 Int *Wi ;
2338 int ok ;
2339
2340 if (CHOLMOD(dump) < -1)
2341 {
2342 /* no checks if debug level is -2 or less */
2343 return (TRUE) ;
2344 }
2345 RETURN_IF_NULL_COMMON (FALSE) ;
2346 Wi = malloc (MAX (1, n) * sizeof (Int)) ;
2347 ok = check_perm (Wi, CHOLMOD(dump), name, Perm, len, n,Common) ;
2348 if (Wi != NULL) free (Wi) ;
2349 return (ok) ;
2350 }
2351
2352
2353 /* ========================================================================== */
2354 /* === cholmod_dump_dense =================================================== */
2355 /* ========================================================================== */
2356
CHOLMOD(dump_dense)2357 int CHOLMOD(dump_dense)
2358 (
2359 cholmod_dense *X,
2360 const char *name,
2361 cholmod_common *Common
2362 )
2363 {
2364 if (CHOLMOD(dump) < -1)
2365 {
2366 /* no checks if debug level is -2 or less */
2367 return (TRUE) ;
2368 }
2369 RETURN_IF_NULL_COMMON (FALSE) ;
2370 return (check_dense (CHOLMOD(dump), name, X, Common)) ;
2371 }
2372
2373
2374 /* ========================================================================== */
2375 /* === cholmod_dump_triplet ================================================= */
2376 /* ========================================================================== */
2377
CHOLMOD(dump_triplet)2378 int CHOLMOD(dump_triplet)
2379 (
2380 cholmod_triplet *T,
2381 const char *name,
2382 cholmod_common *Common
2383 )
2384 {
2385 if (CHOLMOD(dump) < -1)
2386 {
2387 /* no checks if debug level is -2 or less */
2388 return (TRUE) ;
2389 }
2390 RETURN_IF_NULL_COMMON (FALSE) ;
2391 return (check_triplet (CHOLMOD(dump), name, T, Common)) ;
2392 }
2393
2394
2395 /* ========================================================================== */
2396 /* === cholmod_dump_subset ================================================== */
2397 /* ========================================================================== */
2398
CHOLMOD(dump_subset)2399 int CHOLMOD(dump_subset)
2400 (
2401 Int *S,
2402 size_t len,
2403 size_t n,
2404 const char *name,
2405 cholmod_common *Common
2406 )
2407 {
2408 if (CHOLMOD(dump) < -1)
2409 {
2410 /* no checks if debug level is -2 or less */
2411 return (TRUE) ;
2412 }
2413 RETURN_IF_NULL_COMMON (FALSE) ;
2414 return (check_subset (S, len, n, CHOLMOD(dump), name, Common)) ;
2415 }
2416
2417
2418 /* ========================================================================== */
2419 /* === cholmod_dump_parent ================================================== */
2420 /* ========================================================================== */
2421
CHOLMOD(dump_parent)2422 int CHOLMOD(dump_parent)
2423 (
2424 Int *Parent,
2425 size_t n,
2426 const char *name,
2427 cholmod_common *Common
2428 )
2429 {
2430 if (CHOLMOD(dump) < -1)
2431 {
2432 /* no checks if debug level is -2 or less */
2433 return (TRUE) ;
2434 }
2435 RETURN_IF_NULL_COMMON (FALSE) ;
2436 return (check_parent (Parent, n, CHOLMOD(dump), name, Common)) ;
2437 }
2438
2439
2440 /* ========================================================================== */
2441 /* === cholmod_dump_real ==================================================== */
2442 /* ========================================================================== */
2443
CHOLMOD(dump_real)2444 void CHOLMOD(dump_real)
2445 (
2446 const char *name,
2447 Real *X, SuiteSparse_long nrow, SuiteSparse_long ncol, int lower,
2448 int xentry, cholmod_common *Common
2449 )
2450 {
2451 /* dump an nrow-by-ncol real dense matrix */
2452 SuiteSparse_long i, j ;
2453 double x, z ;
2454 if (CHOLMOD(dump) < -1)
2455 {
2456 /* no checks if debug level is -2 or less */
2457 return ;
2458 }
2459 PRINT1 (("%s: dump_real, nrow: %ld ncol: %ld lower: %d\n",
2460 name, nrow, ncol, lower)) ;
2461 for (j = 0 ; j < ncol ; j++)
2462 {
2463 PRINT2 ((" col %ld\n", j)) ;
2464 for (i = 0 ; i < nrow ; i++)
2465 {
2466 /* X is stored in column-major form */
2467 if (lower && i < j)
2468 {
2469 PRINT2 ((" %5ld: -", i)) ;
2470 }
2471 else
2472 {
2473 x = *X ;
2474 PRINT2 ((" %5ld: %e", i, x)) ;
2475 if (xentry == 2)
2476 {
2477 z = *(X+1) ;
2478 PRINT2 ((", %e", z)) ;
2479 }
2480 }
2481 PRINT2 (("\n")) ;
2482 X += xentry ;
2483 }
2484 }
2485 }
2486
2487
2488 /* ========================================================================== */
2489 /* === cholmod_dump_super =================================================== */
2490 /* ========================================================================== */
2491
CHOLMOD(dump_super)2492 void CHOLMOD(dump_super)
2493 (
2494 SuiteSparse_long s,
2495 Int *Super, Int *Lpi, Int *Ls, Int *Lpx, double *Lx,
2496 int xentry,
2497 cholmod_common *Common
2498 )
2499 {
2500 Int k1, k2, do_values, psi, psx, nsrow, nscol, psend, ilast, p, i ;
2501 if (CHOLMOD(dump) < -1)
2502 {
2503 /* no checks if debug level is -2 or less */
2504 return ;
2505 }
2506 k1 = Super [s] ;
2507 k2 = Super [s+1] ;
2508 nscol = k2 - k1 ;
2509 do_values = (Lpx != NULL) && (Lx != NULL) ;
2510 psi = Lpi [s] ;
2511 psend = Lpi [s+1] ;
2512 nsrow = psend - psi ;
2513 PRINT1 (("\nSuper %ld, columns "ID" to "ID", "ID" rows "ID" cols\n",
2514 s, k1, k2-1, nsrow, nscol)) ;
2515 ilast = -1 ;
2516 for (p = psi ; p < psend ; p++)
2517 {
2518 i = Ls [p] ;
2519 PRINT2 ((" "ID" : p-psi "ID"\n", i, p-psi)) ;
2520 ASSERT (IMPLIES (p-psi < nscol, i == k1 + (p-psi))) ;
2521 if (p-psi == nscol-1) PRINT2 (("------\n")) ;
2522 ASSERT (i > ilast) ;
2523 ilast = i ;
2524 }
2525 if (do_values)
2526 {
2527 psx = Lpx [s] ;
2528 CHOLMOD(dump_real) ("Supernode", Lx + xentry*psx, nsrow, nscol, TRUE,
2529 xentry, Common) ;
2530 }
2531 }
2532
2533
2534 /* ========================================================================== */
2535 /* === cholmod_dump_mem ===================================================== */
2536 /* ========================================================================== */
2537
CHOLMOD(dump_mem)2538 int CHOLMOD(dump_mem)
2539 (
2540 const char *where,
2541 SuiteSparse_long should,
2542 cholmod_common *Common
2543 )
2544 {
2545 SuiteSparse_long diff = should - Common->memory_inuse ;
2546 if (diff != 0)
2547 {
2548 PRINT0 (("mem: %-15s peak %10g inuse %10g should %10g\n",
2549 where, (double) Common->memory_usage, (double) Common->memory_inuse,
2550 (double) should)) ;
2551 PRINT0 (("mem: %s diff %ld !\n", where, diff)) ;
2552 }
2553 return (diff == 0) ;
2554 }
2555
2556
2557 /* ========================================================================== */
2558 /* === cholmod_dump_partition =============================================== */
2559 /* ========================================================================== */
2560
2561 /* make sure we have a proper separator (for debugging only)
2562 *
2563 * workspace: none
2564 */
2565
CHOLMOD(dump_partition)2566 int CHOLMOD(dump_partition)
2567 (
2568 SuiteSparse_long n,
2569 Int *Cp,
2570 Int *Ci,
2571 Int *Cnw, /* can be NULL */
2572 Int *Part,
2573 SuiteSparse_long sepsize,
2574 cholmod_common *Common
2575 )
2576 {
2577 Int chek [3], which, ok, i, j, p ;
2578 PRINT1 (("bisect sepsize %ld\n", sepsize)) ;
2579 ok = TRUE ;
2580 chek [0] = 0 ;
2581 chek [1] = 0 ;
2582 chek [2] = 0 ;
2583 for (j = 0 ; j < n ; j++)
2584 {
2585 PRINT2 (("--------j "ID" in part "ID" nw "ID"\n", j, Part [j],
2586 Cnw ? (Cnw[j]):1));
2587 which = Part [j] ;
2588 for (p = Cp [j] ; p < Cp [j+1] ; p++)
2589 {
2590 i = Ci [p] ;
2591 PRINT3 (("i "ID", part "ID"\n", i, Part [i])) ;
2592 if (which == 0)
2593 {
2594 if (Part [i] == 1)
2595 {
2596 PRINT0 (("Error! "ID" "ID"\n", i, j)) ;
2597 ok = FALSE ;
2598 }
2599 }
2600 else if (which == 1)
2601 {
2602 if (Part [i] == 0)
2603 {
2604 PRINT0 (("Error! "ID" "ID"\n", i, j)) ;
2605 ok = FALSE ;
2606 }
2607 }
2608 }
2609 if (which < 0 || which > 2)
2610 {
2611 PRINT0 (("Part out of range\n")) ;
2612 ok = FALSE ;
2613 }
2614 chek [which] += (Cnw ? (Cnw [j]) : 1) ;
2615 }
2616 PRINT1 (("sepsize %ld check "ID" "ID" "ID"\n",
2617 sepsize, chek[0], chek[1],chek[2]));
2618 if (sepsize != chek[2])
2619 {
2620 PRINT0 (("mismatch!\n")) ;
2621 ok = FALSE ;
2622 }
2623 return (ok) ;
2624 }
2625
2626
2627 /* ========================================================================== */
2628 /* === cholmod_dump_work ==================================================== */
2629 /* ========================================================================== */
2630
CHOLMOD(dump_work)2631 int CHOLMOD(dump_work) (int flag, int head, SuiteSparse_long wsize,
2632 cholmod_common *Common)
2633 {
2634 double *W ;
2635 Int *Flag, *Head ;
2636 Int k, nrow, mark ;
2637
2638 if (CHOLMOD(dump) < -1)
2639 {
2640 /* no checks if debug level is -2 or less */
2641 return (TRUE) ;
2642 }
2643
2644 RETURN_IF_NULL_COMMON (FALSE) ;
2645 nrow = Common->nrow ;
2646 Flag = Common->Flag ;
2647 Head = Common->Head ;
2648 W = Common->Xwork ;
2649 mark = Common->mark ;
2650
2651 if (wsize < 0)
2652 {
2653 /* check all of Xwork */
2654 wsize = Common->xworksize ;
2655 }
2656 else
2657 {
2658 /* check on the first wsize doubles in Xwork */
2659 wsize = MIN (wsize, (Int) (Common->xworksize)) ;
2660 }
2661
2662 if (flag)
2663 {
2664 for (k = 0 ; k < nrow ; k++)
2665 {
2666 if (Flag [k] >= mark)
2667 {
2668 PRINT0 (("Flag invalid, Flag ["ID"] = "ID", mark = "ID"\n",
2669 k, Flag [k], mark)) ;
2670 ASSERT (0) ;
2671 return (FALSE) ;
2672 }
2673 }
2674 }
2675
2676 if (head)
2677 {
2678 for (k = 0 ; k < nrow ; k++)
2679 {
2680 if (Head [k] != EMPTY)
2681 {
2682 PRINT0 (("Head invalid, Head ["ID"] = "ID"\n", k, Head [k])) ;
2683 ASSERT (0) ;
2684 return (FALSE) ;
2685 }
2686 }
2687 }
2688
2689 for (k = 0 ; k < wsize ; k++)
2690 {
2691 if (W [k] != 0.)
2692 {
2693 PRINT0 (("W invalid, W ["ID"] = %g\n", k, W [k])) ;
2694 ASSERT (0) ;
2695 return (FALSE) ;
2696 }
2697 }
2698
2699 return (TRUE) ;
2700 }
2701 #endif
2702 #endif
2703