1 /* ========================================================================== */
2 /* === UMFPACK_report_numeric =============================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
9 /* -------------------------------------------------------------------------- */
10
11 /*
12 User-callable. Prints the Numeric object.
13 See umfpack_report_numeric.h for details.
14
15 Dynamic memory usage: Allocates a size n*sizeof(Int) workspace via a single
16 call to UMF_malloc and then frees all of it via UMF_free on return. The
17 workspace is not allocated if an early error return occurs before the
18 workspace is needed.
19 */
20
21 #include "umf_internal.h"
22 #include "umf_valid_numeric.h"
23 #include "umf_report_perm.h"
24 #include "umf_report_vector.h"
25 #include "umf_malloc.h"
26 #include "umf_free.h"
27
28
29 PRIVATE Int report_L
30 (
31 NumericType *Numeric,
32 Int Pattern [ ],
33 Int prl
34 ) ;
35
36
37 PRIVATE Int report_U
38 (
39 NumericType *Numeric,
40 Int Pattern [ ],
41 Int prl
42 ) ;
43
44 /* ========================================================================== */
45 /* === UMFPACK_report_numeric =============================================== */
46 /* ========================================================================== */
47
UMFPACK_report_numeric(void * NumericHandle,const double Control[UMFPACK_CONTROL])48 GLOBAL Int UMFPACK_report_numeric
49 (
50 void *NumericHandle,
51 const double Control [UMFPACK_CONTROL]
52 )
53 {
54 Int prl, *W, nn, n_row, n_col, n_inner, num_fixed_size, numeric_size,
55 npiv ;
56 NumericType *Numeric ;
57
58 prl = GET_CONTROL (UMFPACK_PRL, UMFPACK_DEFAULT_PRL) ;
59
60 if (prl <= 2)
61 {
62 return (UMFPACK_OK) ;
63 }
64
65 PRINTF (("Numeric object: ")) ;
66
67 Numeric = (NumericType *) NumericHandle ;
68 if (!UMF_valid_numeric (Numeric))
69 {
70 PRINTF (("ERROR: LU factors invalid\n\n")) ;
71 return (UMFPACK_ERROR_invalid_Numeric_object) ;
72 }
73
74 n_row = Numeric->n_row ;
75 n_col = Numeric->n_col ;
76 nn = MAX (n_row, n_col) ;
77 n_inner = MIN (n_row, n_col) ;
78 npiv = Numeric->npiv ;
79
80 DEBUG1 (("n_row " ID " n_col " ID " nn " ID " n_inner " ID " npiv " ID "\n",
81 n_row, n_col, nn, n_inner, npiv)) ;
82
83 /* size of Numeric object, except Numeric->Memory and Numeric->Upattern */
84 /* see also UMF_set_stats */
85 num_fixed_size =
86 UNITS (NumericType, 1) /* Numeric structure */
87 + UNITS (Entry, n_inner+1) /* D */
88 + UNITS (Int, n_row+1) /* Rperm */
89 + UNITS (Int, n_col+1) /* Cperm */
90 + 6 * UNITS (Int, npiv+1) /* Lpos, Uilen, Uip, Upos, Lilen, Lip */
91 + ((Numeric->scale != UMFPACK_SCALE_NONE) ?
92 UNITS (Entry, n_row) : 0) ; /* Rs */
93
94 DEBUG1 (("num fixed size: " ID "\n", num_fixed_size)) ;
95 DEBUG1 (("Numeric->size " ID "\n", Numeric->size)) ;
96 DEBUG1 (("ulen units " ID "\n", UNITS (Int, Numeric->ulen))) ;
97
98 /* size of Numeric->Memory is Numeric->size */
99 /* size of Numeric->Upattern is Numeric->ulen */
100 numeric_size = num_fixed_size + Numeric->size
101 + UNITS (Int, Numeric->ulen) ;
102
103 DEBUG1 (("numeric total size " ID "\n", numeric_size)) ;
104
105 if (prl >= 4)
106 {
107 PRINTF (("\n n_row: " ID " n_col: " ID "\n", n_row, n_col)) ;
108
109 PRINTF ((" relative pivot tolerance used: %g\n",
110 Numeric->relpt)) ;
111 PRINTF ((" relative symmetric pivot tolerance used: %g\n",
112 Numeric->relpt2)) ;
113
114 PRINTF ((" matrix scaled: ")) ;
115 if (Numeric->scale == UMFPACK_SCALE_NONE)
116 {
117 PRINTF (("no")) ;
118 }
119 else if (Numeric->scale == UMFPACK_SCALE_SUM)
120 {
121 PRINTF (("yes (divided each row by sum abs value in each row)\n")) ;
122 PRINTF ((" minimum sum (abs (rows of A)): %.5e\n",
123 Numeric->rsmin)) ;
124 PRINTF ((" maximum sum (abs (rows of A)): %.5e",
125 Numeric->rsmax)) ;
126 }
127 else if (Numeric->scale == UMFPACK_SCALE_MAX)
128 {
129 PRINTF (("yes (divided each row by max abs value in each row)\n")) ;
130 PRINTF ((" minimum max (abs (rows of A)): %.5e\n",
131 Numeric->rsmin)) ;
132 PRINTF ((" maximum max (abs (rows of A)): %.5e",
133 Numeric->rsmax)) ;
134 }
135 PRINTF (("\n")) ;
136
137 PRINTF ((" initial allocation parameter used: %g\n",
138 Numeric->alloc_init)) ;
139 PRINTF ((" frontal matrix allocation parameter used: %g\n",
140 Numeric->front_alloc_init)) ;
141 PRINTF ((" final total size of Numeric object (Units): " ID "\n",
142 numeric_size)) ;
143 PRINTF ((" final total size of Numeric object (MBytes): %.1f\n",
144 MBYTES (numeric_size))) ;
145 PRINTF ((" peak size of variable-size part (Units): " ID "\n",
146 Numeric->max_usage)) ;
147 PRINTF ((" peak size of variable-size part (MBytes): %.1f\n",
148 MBYTES (Numeric->max_usage))) ;
149 PRINTF ((" largest actual frontal matrix size: " ID "\n",
150 Numeric->maxfrsize)) ;
151 PRINTF ((" memory defragmentations: " ID "\n",
152 Numeric->ngarbage)) ;
153 PRINTF ((" memory reallocations: " ID "\n",
154 Numeric->nrealloc)) ;
155 PRINTF ((" costly memory reallocations: " ID "\n",
156 Numeric->ncostly)) ;
157 PRINTF ((" entries in compressed pattern (L and U): " ID "\n",
158 Numeric->isize)) ;
159 PRINTF ((" number of nonzeros in L (excl diag): " ID "\n",
160 Numeric->lnz)) ;
161 PRINTF ((" number of entries stored in L (excl diag): " ID "\n",
162 Numeric->nLentries)) ;
163 PRINTF ((" number of nonzeros in U (excl diag): " ID "\n",
164 Numeric->unz)) ;
165 PRINTF ((" number of entries stored in U (excl diag): " ID "\n",
166 Numeric->nUentries)) ;
167 PRINTF ((" factorization floating-point operations: %g\n",
168 Numeric->flops)) ;
169 PRINTF ((" number of nonzeros on diagonal of U: " ID "\n",
170 Numeric->nnzpiv)) ;
171 PRINTF ((" min abs. value on diagonal of U: %.5e\n",
172 Numeric->min_udiag)) ;
173 PRINTF ((" max abs. value on diagonal of U: %.5e\n",
174 Numeric->max_udiag)) ;
175 PRINTF ((" reciprocal condition number estimate: %.2e\n",
176 Numeric->rcond)) ;
177 }
178
179 W = (Int *) UMF_malloc (nn, sizeof (Int)) ;
180 if (!W)
181 {
182 PRINTF ((" ERROR: out of memory to check Numeric object\n\n")) ;
183 return (UMFPACK_ERROR_out_of_memory) ;
184 }
185
186 if (Numeric->Rs)
187 {
188 #ifndef NRECIPROCAL
189 if (Numeric->do_recip)
190 {
191 PRINTF4 (("\nScale factors applied via multiplication\n")) ;
192 }
193 else
194 #endif
195 {
196 PRINTF4 (("\nScale factors applied via division\n")) ;
197 }
198 PRINTF4 (("Scale factors, Rs: ")) ;
199 (void) UMF_report_vector (n_row, Numeric->Rs, (double *) NULL,
200 prl, FALSE, TRUE) ;
201 }
202 else
203 {
204 PRINTF4 (("Scale factors, Rs: (not present)\n")) ;
205 }
206
207 PRINTF4 (("\nP: row ")) ;
208 if (UMF_report_perm (n_row, Numeric->Rperm, W, prl, 0) != UMFPACK_OK)
209 {
210 (void) UMF_free ((void *) W) ;
211 return (UMFPACK_ERROR_invalid_Numeric_object) ;
212 }
213
214 PRINTF4 (("\nQ: column ")) ;
215 if (UMF_report_perm (n_col, Numeric->Cperm, W, prl, 0) != UMFPACK_OK)
216 {
217 (void) UMF_free ((void *) W) ;
218 return (UMFPACK_ERROR_invalid_Numeric_object) ;
219 }
220
221 if (!report_L (Numeric, W, prl))
222 {
223 (void) UMF_free ((void *) W) ;
224 PRINTF ((" ERROR: L factor invalid\n\n")) ;
225 return (UMFPACK_ERROR_invalid_Numeric_object) ;
226 }
227
228 if (!report_U (Numeric, W, prl))
229 {
230 (void) UMF_free ((void *) W) ;
231 PRINTF ((" ERROR: U factor invalid\n\n")) ;
232 return (UMFPACK_ERROR_invalid_Numeric_object) ;
233 }
234
235 /* The diagonal of U is in "merged" (Entry) form, not "split" form. */
236 PRINTF4 (("\ndiagonal of U: ")) ;
237 (void) UMF_report_vector (n_inner, (double *) Numeric->D, (double *) NULL,
238 prl, FALSE, FALSE) ;
239
240 (void) UMF_free ((void *) W) ;
241
242 PRINTF4 ((" Numeric object: ")) ;
243 PRINTF (("OK\n\n")) ;
244 return (UMFPACK_OK) ;
245 }
246
247
248 /* ========================================================================== */
249 /* === report_L ============================================================= */
250 /* ========================================================================== */
251
report_L(NumericType * Numeric,Int Pattern[],Int prl)252 PRIVATE Int report_L
253 (
254 NumericType *Numeric,
255 Int Pattern [ ],
256 Int prl
257 )
258 {
259 Int k, deg, *ip, j, row, n_row, *Lpos, *Lilen, valid, k1,
260 *Lip, newLchain, llen, prl1, pos, lp, p, npiv, n1, *Li ;
261 Entry *xp, *Lval ;
262
263 /* ---------------------------------------------------------------------- */
264
265 ASSERT (prl >= 3) ;
266
267 n_row = Numeric->n_row ;
268 npiv = Numeric->npiv ;
269 n1 = Numeric->n1 ;
270 Lpos = Numeric->Lpos ;
271 Lilen = Numeric->Lilen ;
272 Lip = Numeric->Lip ;
273 prl1 = prl ;
274 deg = 0 ;
275
276 PRINTF4 ((
277 "\nL in Numeric object, in column-oriented compressed-pattern form:\n"
278 " Diagonal entries are all equal to 1.0 (not stored)\n")) ;
279
280 ASSERT (Pattern != (Int *) NULL) ;
281
282 /* ---------------------------------------------------------------------- */
283 /* print L */
284 /* ---------------------------------------------------------------------- */
285
286 k1 = 12 ;
287
288 /* ---------------------------------------------------------------------- */
289 /* print the singleton columns of L */
290 /* ---------------------------------------------------------------------- */
291
292 for (k = 0 ; k < n1 ; k++)
293 {
294 if (k1 > 0)
295 {
296 prl = prl1 ;
297 }
298 lp = Lip [k] ;
299 deg = Lilen [k] ;
300 Li = (Int *) (Numeric->Memory + lp) ;
301 lp += UNITS (Int, deg) ;
302 Lval = (Entry *) (Numeric->Memory + lp) ;
303 if (k1-- > 0)
304 {
305 prl = prl1 ;
306 }
307 else if (prl == 4)
308 {
309 PRINTF ((" ...\n")) ;
310 prl-- ;
311 }
312 PRINTF4 (("\n column " ID ":", INDEX (k))) ;
313 PRINTF4 ((" length " ID ".\n", deg)) ;
314 for (j = 0 ; j < deg ; j++)
315 {
316 row = Li [j] ;
317 PRINTF4 (("\trow " ID " : ", INDEX (row))) ;
318 if (prl >= 4) PRINT_ENTRY (Lval [j]) ;
319 if (row <= k || row >= n_row)
320 {
321 return (FALSE) ;
322 }
323 PRINTF4 (("\n")) ;
324 /* truncate printout, but continue to check L */
325 if (prl == 4 && j == 9 && deg > 10)
326 {
327 PRINTF (("\t...\n")) ;
328 prl-- ;
329 }
330 }
331 }
332
333 /* ---------------------------------------------------------------------- */
334 /* print the regular columns of L */
335 /* ---------------------------------------------------------------------- */
336
337 for (k = n1 ; k < npiv ; k++)
338 {
339 /* if prl is 4, print the first 10 entries of the first 10 columns */
340 if (k1 > 0)
341 {
342 prl = prl1 ;
343 }
344
345 lp = Lip [k] ;
346 newLchain = (lp < 0) ;
347 if (newLchain)
348 {
349 lp = -lp ;
350 deg = 0 ;
351 }
352
353 if (k1-- > 0)
354 {
355 prl = prl1 ;
356 }
357 else if (prl == 4)
358 {
359 PRINTF ((" ...\n")) ;
360 prl-- ;
361 }
362
363 PRINTF4 (("\n column " ID ":", INDEX (k))) ;
364
365 /* ------------------------------------------------------------------ */
366 /* make column of L in Pattern [0..deg-1] */
367 /* ------------------------------------------------------------------ */
368
369 /* remove pivot row */
370 pos = Lpos [k] ;
371 if (pos != EMPTY)
372 {
373 PRINTF4 ((" remove row " ID " at position " ID ".",
374 INDEX (Pattern [pos]), INDEX (pos))) ;
375 valid = (!newLchain) && (deg > 0) && (pos < deg) && (pos >= 0)
376 && (Pattern [pos] == k) ;
377 if (!valid)
378 {
379 return (FALSE) ;
380 }
381 Pattern [pos] = Pattern [--deg] ;
382 }
383
384 /* concatenate the pattern */
385 llen = Lilen [k] ;
386 if (llen < 0)
387 {
388 return (FALSE) ;
389 }
390 p = lp + UNITS (Int, llen) ;
391 xp = (Entry *) (Numeric->Memory + p) ;
392 if ((llen > 0 || deg > 0)
393 && (p + (Int) UNITS (Entry, deg) > Numeric->size))
394 {
395 return (FALSE) ;
396 }
397 if (llen > 0)
398 {
399 PRINTF4 ((" add " ID " entries.", llen)) ;
400 ip = (Int *) (Numeric->Memory + lp) ;
401 for (j = 0 ; j < llen ; j++)
402 {
403 Pattern [deg++] = *ip++ ;
404 }
405 }
406
407 /* ------------------------------------------------------------------ */
408 /* print column k of L */
409 /* ------------------------------------------------------------------ */
410
411 PRINTF4 ((" length " ID ".", deg)) ;
412 if (newLchain)
413 {
414 PRINTF4 ((" Start of Lchain.")) ;
415 }
416 PRINTF4 (("\n")) ;
417
418 for (j = 0 ; j < deg ; j++)
419 {
420 row = Pattern [j] ;
421 PRINTF4 (("\trow " ID " : ", INDEX (row))) ;
422 if (prl >= 4) PRINT_ENTRY (*xp) ;
423 if (row <= k || row >= n_row)
424 {
425 return (FALSE) ;
426 }
427 PRINTF4 (("\n")) ;
428 xp++ ;
429 /* truncate printout, but continue to check L */
430 if (prl == 4 && j == 9 && deg > 10)
431 {
432 PRINTF (("\t...\n")) ;
433 prl-- ;
434 }
435 }
436 }
437
438 PRINTF4 (("\n")) ;
439 return (TRUE) ;
440 }
441
442
443 /* ========================================================================== */
444 /* === report_U ============================================================= */
445 /* ========================================================================== */
446
report_U(NumericType * Numeric,Int Pattern[],Int prl)447 PRIVATE Int report_U
448 (
449 NumericType *Numeric,
450 Int Pattern [ ],
451 Int prl
452 )
453 {
454 /* ---------------------------------------------------------------------- */
455
456 Int k, deg, j, *ip, col, *Upos, *Uilen, k1, prl1, pos,
457 *Uip, n_col, ulen, p, newUchain, up, npiv, n1, *Ui ;
458 Entry *xp, *Uval ;
459
460 /* ---------------------------------------------------------------------- */
461
462 ASSERT (prl >= 3) ;
463
464 n_col = Numeric->n_col ;
465 npiv = Numeric->npiv ;
466 n1 = Numeric->n1 ;
467 Upos = Numeric->Upos ;
468 Uilen = Numeric->Uilen ;
469 Uip = Numeric->Uip ;
470 prl1 = prl ;
471
472 PRINTF4 ((
473 "\nU in Numeric object, in row-oriented compressed-pattern form:\n"
474 " Diagonal is stored separately.\n")) ;
475
476 ASSERT (Pattern != (Int *) NULL) ;
477
478 k1 = 12 ;
479
480 /* ---------------------------------------------------------------------- */
481 /* print the sparse part of U */
482 /* ---------------------------------------------------------------------- */
483
484 deg = Numeric->ulen ;
485 if (deg > 0)
486 {
487 /* make last pivot row of U (singular matrices only) */
488 for (j = 0 ; j < deg ; j++)
489 {
490 Pattern [j] = Numeric->Upattern [j] ;
491 }
492 }
493
494 PRINTF4 (("\n row " ID ": length " ID ". End of Uchain.\n", INDEX (npiv-1),
495 deg)) ;
496
497 for (k = npiv-1 ; k >= n1 ; k--)
498 {
499
500 /* ------------------------------------------------------------------ */
501 /* print row k of U */
502 /* ------------------------------------------------------------------ */
503
504 /* if prl is 3, print the first 10 entries of the first 10 columns */
505 if (k1 > 0)
506 {
507 prl = prl1 ;
508 }
509
510 up = Uip [k] ;
511 ulen = Uilen [k] ;
512 if (ulen < 0)
513 {
514 return (FALSE) ;
515 }
516 newUchain = (up < 0) ;
517 if (newUchain)
518 {
519 up = -up ;
520 p = up + UNITS (Int, ulen) ;
521 }
522 else
523 {
524 p = up ;
525 }
526 xp = (Entry *) (Numeric->Memory + p) ;
527 if (deg > 0 && (p + (Int) UNITS (Entry, deg) > Numeric->size))
528 {
529 return (FALSE) ;
530 }
531 for (j = 0 ; j < deg ; j++)
532 {
533 col = Pattern [j] ;
534 PRINTF4 (("\tcol " ID " :", INDEX (col))) ;
535 if (prl >= 4) PRINT_ENTRY (*xp) ;
536 if (col <= k || col >= n_col)
537 {
538 return (FALSE) ;
539 }
540 PRINTF4 (("\n")) ;
541 xp++ ;
542 /* truncate printout, but continue to check U */
543 if (prl == 4 && j == 9 && deg > 10)
544 {
545 PRINTF (("\t...\n")) ;
546 prl-- ;
547 }
548 }
549
550 /* ------------------------------------------------------------------ */
551 /* make row k-1 of U in Pattern [0..deg-1] */
552 /* ------------------------------------------------------------------ */
553
554 if (k1-- > 0)
555 {
556 prl = prl1 ;
557 }
558 else if (prl == 4)
559 {
560 PRINTF ((" ...\n")) ;
561 prl-- ;
562 }
563
564 if (k > 0)
565 {
566 PRINTF4 (("\n row " ID ": ", INDEX (k-1))) ;
567 }
568
569 if (newUchain)
570 {
571 /* next row is a new Uchain */
572 if (k > 0)
573 {
574 deg = ulen ;
575 PRINTF4 (("length " ID ". End of Uchain.\n", deg)) ;
576 if (up + (Int) UNITS (Int, ulen) > Numeric->size)
577 {
578 return (FALSE) ;
579 }
580 ip = (Int *) (Numeric->Memory + up) ;
581 for (j = 0 ; j < deg ; j++)
582 {
583 Pattern [j] = *ip++ ;
584 }
585 }
586 }
587 else
588 {
589 if (ulen > 0)
590 {
591 PRINTF4 (("remove " ID " entries. ", ulen)) ;
592 }
593 deg -= ulen ;
594 if (deg < 0)
595 {
596 return (FALSE) ;
597 }
598 pos = Upos [k] ;
599 if (pos != EMPTY)
600 {
601 /* add the pivot column */
602 PRINTF4 (("add column " ID " at position " ID ". ",
603 INDEX (k), INDEX (pos))) ;
604 if (pos < 0 || pos > deg)
605 {
606 return (FALSE) ;
607 }
608 Pattern [deg++] = Pattern [pos] ;
609 Pattern [pos] = k ;
610 }
611 PRINTF4 (("length " ID ".\n", deg)) ;
612 }
613 }
614
615 /* ---------------------------------------------------------------------- */
616 /* print the singleton rows of U */
617 /* ---------------------------------------------------------------------- */
618
619 for (k = n1 - 1 ; k >= 0 ; k--)
620 {
621 if (k1 > 0)
622 {
623 prl = prl1 ;
624 }
625 up = Uip [k] ;
626 deg = Uilen [k] ;
627 Ui = (Int *) (Numeric->Memory + up) ;
628 up += UNITS (Int, deg) ;
629 Uval = (Entry *) (Numeric->Memory + up) ;
630 if (k1-- > 0)
631 {
632 prl = prl1 ;
633 }
634 else if (prl == 4)
635 {
636 PRINTF ((" ...\n")) ;
637 prl-- ;
638 }
639 PRINTF4 (("\n row " ID ":", INDEX (k))) ;
640 PRINTF4 ((" length " ID ".\n", deg)) ;
641 for (j = 0 ; j < deg ; j++)
642 {
643 col = Ui [j] ;
644 PRINTF4 (("\tcol " ID " : ", INDEX (col))) ;
645 if (prl >= 4) PRINT_ENTRY (Uval [j]) ;
646 if (col <= k || col >= n_col)
647 {
648 return (FALSE) ;
649 }
650 PRINTF4 (("\n")) ;
651 /* truncate printout, but continue to check U */
652 if (prl == 4 && j == 9 && deg > 10)
653 {
654 PRINTF (("\t...\n")) ;
655 prl-- ;
656 }
657 }
658 }
659
660 prl = prl1 ;
661 PRINTF4 (("\n")) ;
662 return (TRUE) ;
663 }
664