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