Amesos2 - Direct Sparse Solver Interfaces Version of the Day
klu2_diagnostics.hpp
1/* ========================================================================== */
2/* === KLU_diagnostics ====================================================== */
3/* ========================================================================== */
4// @HEADER
5// ***********************************************************************
6//
7// KLU2: A Direct Linear Solver package
8// Copyright 2011 Sandia Corporation
9//
10// Under terms of Contract DE-AC04-94AL85000, with Sandia Corporation, the
11// U.S. Government retains certain rights in this software.
12//
13// This library is free software; you can redistribute it and/or modify
14// it under the terms of the GNU Lesser General Public License as
15// published by the Free Software Foundation; either version 2.1 of the
16// License, or (at your option) any later version.
17//
18// This library is distributed in the hope that it will be useful, but
19// WITHOUT ANY WARRANTY; without even the implied warranty of
20// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21// Lesser General Public License for more details.
22//
23// You should have received a copy of the GNU Lesser General Public
24// License along with this library; if not, write to the Free Software
25// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
26// USA
27// Questions? Contact Mike A. Heroux (maherou@sandia.gov)
28//
29// KLU2 is derived work from KLU, licensed under LGPL, and copyrighted by
30// University of Florida. The Authors of KLU are Timothy A. Davis and
31// Eka Palamadai. See Doc/KLU_README.txt for the licensing and copyright
32// information for KLU.
33//
34// ***********************************************************************
35// @HEADER
36
37/* Linear algebraic diagnostics:
38 * KLU_rgrowth: reciprocal pivot growth, takes O(|A|+|U|) time
39 * KLU_condest: condition number estimator, takes about O(|A|+5*(|L|+|U|)) time
40 * KLU_flops: compute # flops required to factorize A into L*U
41 * KLU_rcond: compute a really cheap estimate of the reciprocal of the
42 * condition number, min(abs(diag(U))) / max(abs(diag(U))).
43 * Takes O(n) time.
44 */
45
46#ifndef KLU2_DIAGNOSTICS_HPP
47#define KLU2_DIAGNOSTICS_HPP
48
49#include "klu2_internal.h"
50#include "klu2_tsolve.hpp"
51
52/* ========================================================================== */
53/* === KLU_rgrowth ========================================================== */
54/* ========================================================================== */
55
56/* Compute the reciprocal pivot growth factor. In MATLAB notation:
57 *
58 * rgrowth = min (max (abs ((R \ A (p,q)) - F))) ./ max (abs (U)))
59 */
60
61template <typename Entry, typename Int>
62Int KLU_rgrowth /* return TRUE if successful, FALSE otherwise */
63(
64 Int *Ap,
65 Int *Ai,
66 double *Ax,
67 KLU_symbolic<Entry, Int> *Symbolic,
68 KLU_numeric<Entry, Int> *Numeric,
69 KLU_common<Entry, Int> *Common
70)
71{
72 double temp, max_ai, max_ui, min_block_rgrowth ;
73 Entry aik ;
74 Int *Q, *Ui, *Uip, *Ulen, *Pinv ;
75 Unit *LU ;
76 Entry *Aentry, *Ux, *Ukk ;
77 double *Rs ;
78 Int i, newrow, oldrow, k1, k2, nk, j, oldcol, k, pend, len ;
79
80 /* ---------------------------------------------------------------------- */
81 /* check inputs */
82 /* ---------------------------------------------------------------------- */
83
84 if (Common == NULL)
85 {
86 return (FALSE) ;
87 }
88
89 if (Symbolic == NULL || Ap == NULL || Ai == NULL || Ax == NULL)
90 {
91 Common->status = KLU_INVALID ;
92 return (FALSE) ;
93 }
94
95 if (Numeric == NULL)
96 {
97 /* treat this as a singular matrix */
98 Common->rgrowth = 0 ;
99 Common->status = KLU_SINGULAR ;
100 return (TRUE) ;
101 }
102 Common->status = KLU_OK ;
103
104 /* ---------------------------------------------------------------------- */
105 /* compute the reciprocal pivot growth */
106 /* ---------------------------------------------------------------------- */
107
108 Aentry = (Entry *) Ax ;
109 Pinv = Numeric->Pinv ;
110 Rs = Numeric->Rs ;
111 Q = Symbolic->Q ;
112 Common->rgrowth = 1 ;
113
114 for (i = 0 ; i < Symbolic->nblocks ; i++)
115 {
116 k1 = Symbolic->R[i] ;
117 k2 = Symbolic->R[i+1] ;
118 nk = k2 - k1 ;
119 if (nk == 1)
120 {
121 continue ; /* skip singleton blocks */
122 }
123 LU = (Unit *) Numeric->LUbx[i] ;
124 Uip = Numeric->Uip + k1 ;
125 Ulen = Numeric->Ulen + k1 ;
126 Ukk = ((Entry *) Numeric->Udiag) + k1 ;
127 min_block_rgrowth = 1 ;
128 for (j = 0 ; j < nk ; j++)
129 {
130 max_ai = 0 ;
131 max_ui = 0 ;
132 oldcol = Q[j + k1] ;
133 pend = Ap [oldcol + 1] ;
134 for (k = Ap [oldcol] ; k < pend ; k++)
135 {
136 oldrow = Ai [k] ;
137 newrow = Pinv [oldrow] ;
138 if (newrow < k1)
139 {
140 continue ; /* skip entry outside the block */
141 }
142 ASSERT (newrow < k2) ;
143 if (Rs != NULL)
144 {
145 /* aik = Aentry [k] / Rs [oldrow] */
146 SCALE_DIV_ASSIGN (aik, Aentry [k], Rs [newrow]) ;
147 }
148 else
149 {
150 aik = Aentry [k] ;
151 }
152 /* temp = ABS (aik) */
153 KLU2_ABS (temp, aik) ;
154 if (temp > max_ai)
155 {
156 max_ai = temp ;
157 }
158 }
159
160 GET_POINTER (LU, Uip, Ulen, Ui, Ux, j, len) ;
161 for (k = 0 ; k < len ; k++)
162 {
163 /* temp = ABS (Ux [k]) */
164 KLU2_ABS (temp, Ux [k]) ;
165 if (temp > max_ui)
166 {
167 max_ui = temp ;
168 }
169 }
170 /* consider the diagonal element */
171 KLU2_ABS (temp, Ukk [j]) ;
172 if (temp > max_ui)
173 {
174 max_ui = temp ;
175 }
176
177 /* if max_ui is 0, skip the column */
178 if (SCALAR_IS_ZERO (max_ui))
179 {
180 continue ;
181 }
182 temp = max_ai / max_ui ;
183 if (temp < min_block_rgrowth)
184 {
185 min_block_rgrowth = temp ;
186 }
187 }
188
189 if (min_block_rgrowth < Common->rgrowth)
190 {
191 Common->rgrowth = min_block_rgrowth ;
192 }
193 }
194 return (TRUE) ;
195}
196
197
198/* ========================================================================== */
199/* === KLU_condest ========================================================== */
200/* ========================================================================== */
201
202/* Estimate the condition number. Uses Higham and Tisseur's algorithm
203 * (A block algorithm for matrix 1-norm estimation, with applications to
204 * 1-norm pseudospectra, SIAM J. Matrix Anal. Appl., 21(4):1185-1201, 2000.
205 */
206
207template <typename Entry, typename Int>
208Int KLU_condest /* return TRUE if successful, FALSE otherwise */
209(
210 Int Ap [ ],
211 double Ax [ ], /* TODO : Check !!! */
212 KLU_symbolic<Entry, Int> *Symbolic,
213 KLU_numeric<Entry, Int> *Numeric,
214 KLU_common<Entry, Int> *Common
215)
216{
217 double xj, Xmax, csum, anorm, ainv_norm, est_old, est_new, abs_value ;
218 Entry *Udiag, *Aentry, *X, *S ;
219 Int *R ;
220 Int nblocks, i, j, jmax, jnew, pend, n ;
221#ifndef COMPLEX
222 Int unchanged ;
223#endif
224
225 /* ---------------------------------------------------------------------- */
226 /* check inputs */
227 /* ---------------------------------------------------------------------- */
228
229 if (Common == NULL)
230 {
231 return (FALSE) ;
232 }
233 if (Symbolic == NULL || Ap == NULL || Ax == NULL)
234 {
235 Common->status = KLU_INVALID ;
236 return (FALSE) ;
237 }
238 abs_value = 0 ;
239 if (Numeric == NULL)
240 {
241 /* treat this as a singular matrix */
242 Common->condest = 1 / abs_value ;
243 Common->status = KLU_SINGULAR ;
244 return (TRUE) ;
245 }
246 Common->status = KLU_OK ;
247
248 /* ---------------------------------------------------------------------- */
249 /* get inputs */
250 /* ---------------------------------------------------------------------- */
251
252 n = Symbolic->n ;
253 nblocks = Symbolic->nblocks ;
254 R = Symbolic->R ;
255 Udiag = (Entry *) Numeric->Udiag ;
256
257 /* ---------------------------------------------------------------------- */
258 /* check if diagonal of U has a zero on it */
259 /* ---------------------------------------------------------------------- */
260
261 for (i = 0 ; i < n ; i++)
262 {
263 KLU2_ABS (abs_value, Udiag [i]) ;
264 if (SCALAR_IS_ZERO (abs_value))
265 {
266 Common->condest = 1 / abs_value ;
267 Common->status = KLU_SINGULAR ;
268 return (TRUE) ;
269 }
270 }
271
272 /* ---------------------------------------------------------------------- */
273 /* compute 1-norm (maximum column sum) of the matrix */
274 /* ---------------------------------------------------------------------- */
275
276 anorm = 0.0 ;
277 Aentry = (Entry *) Ax ;
278 for (i = 0 ; i < n ; i++)
279 {
280 pend = Ap [i + 1] ;
281 csum = 0.0 ;
282 for (j = Ap [i] ; j < pend ; j++)
283 {
284 KLU2_ABS (abs_value, Aentry [j]) ;
285 csum += abs_value ;
286 }
287 if (csum > anorm)
288 {
289 anorm = csum ;
290 }
291 }
292
293 /* ---------------------------------------------------------------------- */
294 /* compute estimate of 1-norm of inv (A) */
295 /* ---------------------------------------------------------------------- */
296
297 /* get workspace (size 2*n Entry's) */
298 X = (Entry *) Numeric->Xwork ; /* size n space used in KLU_solve, tsolve */
299 X += n ; /* X is size n */
300 S = X + n ; /* S is size n */
301
302 for (i = 0 ; i < n ; i++)
303 {
304 CLEAR (S [i]) ;
305 CLEAR (X [i]) ;
306 /* TODO : Check REAL(X[i]) -> X[i]*/
307 /*REAL (X [i]) = 1.0 / ((double) n) ;*/
308 X [i] = 1.0 / ((double) n) ;
309 }
310 jmax = 0 ;
311
312 ainv_norm = 0.0 ;
313 for (i = 0 ; i < 5 ; i++)
314 {
315 if (i > 0)
316 {
317 /* X [jmax] is the largest entry in X */
318 for (j = 0 ; j < n ; j++)
319 {
320 /* X [j] = 0 ;*/
321 CLEAR (X [j]) ;
322 }
323 /* TODO : Check REAL(X[i]) -> X[i]*/
324 /*REAL (X [jmax]) = 1 ;*/
325 X [jmax] = 1 ;
326 }
327
328 KLU_solve (Symbolic, Numeric, n, 1, (double *) X, Common) ;
329 est_old = ainv_norm ;
330 ainv_norm = 0.0 ;
331
332 for (j = 0 ; j < n ; j++)
333 {
334 /* ainv_norm += ABS (X [j]) ;*/
335 KLU2_ABS (abs_value, X [j]) ;
336 ainv_norm += abs_value ;
337 }
338
339#ifndef COMPLEX
340 unchanged = TRUE ;
341
342 for (j = 0 ; j < n ; j++)
343 {
344 double s = (X [j] >= 0) ? 1 : -1 ;
345 if (s != (Int) REAL (S [j]))
346 {
347 S [j] = s ;
348 unchanged = FALSE ;
349 }
350 }
351
352 if (i > 0 && (ainv_norm <= est_old || unchanged))
353 {
354 break ;
355 }
356#else
357 for (j = 0 ; j < n ; j++)
358 {
359 if (IS_NONZERO (X [j]))
360 {
361 KLU2_ABS (abs_value, X [j]) ;
362 SCALE_DIV_ASSIGN (S [j], X [j], abs_value) ;
363 }
364 else
365 {
366 CLEAR (S [j]) ;
367 /* TODO : Check REAL(X[i]) -> X[i]*/
368 /*REAL (S [j]) = 1 ; */
369 S [j] = 1 ;
370 }
371 }
372
373 if (i > 0 && ainv_norm <= est_old)
374 {
375 break ;
376 }
377#endif
378
379 for (j = 0 ; j < n ; j++)
380 {
381 X [j] = S [j] ;
382 }
383
384#ifndef COMPLEX
385 /* do a transpose solve */
386 KLU_tsolve (Symbolic, Numeric, n, 1, X, Common) ;
387#else
388 /* do a conjugate transpose solve */
389 KLU_tsolve (Symbolic, Numeric, n, 1, (double *) X, 1, Common) ;
390#endif
391
392 /* jnew = the position of the largest entry in X */
393 jnew = 0 ;
394 Xmax = 0 ;
395 for (j = 0 ; j < n ; j++)
396 {
397 /* xj = ABS (X [j]) ;*/
398 KLU2_ABS (xj, X [j]) ;
399 if (xj > Xmax)
400 {
401 Xmax = xj ;
402 jnew = j ;
403 }
404 }
405 if (i > 0 && jnew == jmax)
406 {
407 /* the position of the largest entry did not change
408 * from the previous iteration */
409 break ;
410 }
411 jmax = jnew ;
412 }
413
414 /* ---------------------------------------------------------------------- */
415 /* compute another estimate of norm(inv(A),1), and take the largest one */
416 /* ---------------------------------------------------------------------- */
417
418 for (j = 0 ; j < n ; j++)
419 {
420 CLEAR (X [j]) ;
421 if (j % 2)
422 {
423 /* TODO : Check REAL(X[i]) -> X[i]*/
424 /*REAL (X [j]) = 1 + ((double) j) / ((double) (n-1)) ;*/
425 X [j] = 1 + ((double) j) / ((double) (n-1)) ;
426 }
427 else
428 {
429 /* TODO : Check REAL(X[i]) -> X[i]*/
430 /*REAL (X [j]) = -1 - ((double) j) / ((double) (n-1)) ;*/
431 X [j] = -1 - ((double) j) / ((double) (n-1)) ;
432 }
433 }
434
435 KLU_solve (Symbolic, Numeric, n, 1, (double *) X, Common) ;
436
437 est_new = 0.0 ;
438 for (j = 0 ; j < n ; j++)
439 {
440 /* est_new += ABS (X [j]) ;*/
441 KLU2_ABS (abs_value, X [j]) ;
442 est_new += abs_value ;
443 }
444 est_new = 2 * est_new / (3 * n) ;
445 ainv_norm = MAX (est_new, ainv_norm) ;
446
447 /* ---------------------------------------------------------------------- */
448 /* compute estimate of condition number */
449 /* ---------------------------------------------------------------------- */
450
451 Common->condest = ainv_norm * anorm ;
452 return (TRUE) ;
453}
454
455
456/* ========================================================================== */
457/* === KLU_flops ============================================================ */
458/* ========================================================================== */
459
460/* Compute the flop count for the LU factorization (in Common->flops) */
461
462template <typename Entry, typename Int>
463Int KLU_flops /* return TRUE if successful, FALSE otherwise */
464(
465 KLU_symbolic<Entry, Int> *Symbolic,
466 KLU_numeric<Entry, Int> *Numeric,
467 KLU_common<Entry, Int> *Common
468)
469{
470 double flops = 0 ;
471 Int *R, *Ui, *Uip, *Llen, *Ulen ;
472 Unit **LUbx ;
473 Unit *LU ;
474 Int k, ulen, p, n, nk, block, nblocks, k1 ;
475
476 /* ---------------------------------------------------------------------- */
477 /* check inputs */
478 /* ---------------------------------------------------------------------- */
479
480 if (Common == NULL)
481 {
482 return (FALSE) ;
483 }
484 Common->flops = AMESOS2_KLU2_EMPTY ;
485 if (Numeric == NULL || Symbolic == NULL)
486 {
487 Common->status = KLU_INVALID ;
488 return (FALSE) ;
489 }
490 Common->status = KLU_OK ;
491
492 /* ---------------------------------------------------------------------- */
493 /* get the contents of the Symbolic object */
494 /* ---------------------------------------------------------------------- */
495
496 n = Symbolic->n ;
497 R = Symbolic->R ;
498 nblocks = Symbolic->nblocks ;
499
500 /* ---------------------------------------------------------------------- */
501 /* get the contents of the Numeric object */
502 /* ---------------------------------------------------------------------- */
503
504 LUbx = (Unit **) Numeric->LUbx ;
505
506 /* ---------------------------------------------------------------------- */
507 /* compute the flop count */
508 /* ---------------------------------------------------------------------- */
509
510 for (block = 0 ; block < nblocks ; block++)
511 {
512 k1 = R [block] ;
513 nk = R [block+1] - k1 ;
514 if (nk > 1)
515 {
516 Llen = Numeric->Llen + k1 ;
517 Uip = Numeric->Uip + k1 ;
518 Ulen = Numeric->Ulen + k1 ;
519 LU = LUbx [block] ;
520 for (k = 0 ; k < nk ; k++)
521 {
522 /* compute kth column of U, and update kth column of A */
523 GET_I_POINTER (LU, Uip, Ui, k) ;
524 ulen = Ulen [k] ;
525 for (p = 0 ; p < ulen ; p++)
526 {
527 flops += 2 * Llen [Ui [p]] ;
528 }
529 /* gather and divide by pivot to get kth column of L */
530 flops += Llen [k] ;
531 }
532 }
533 }
534 Common->flops = flops ;
535 return (TRUE) ;
536}
537
538
539/* ========================================================================== */
540/* === KLU_rcond ============================================================ */
541/* ========================================================================== */
542
543/* Compute a really cheap estimate of the reciprocal of the condition number,
544 * condition number, min(abs(diag(U))) / max(abs(diag(U))). If U has a zero
545 * pivot, or a NaN pivot, rcond will be zero. Takes O(n) time.
546 */
547
548template <typename Entry, typename Int>
549Int KLU_rcond /* return TRUE if successful, FALSE otherwise */
550(
551 KLU_symbolic<Entry, Int> *Symbolic, /* input, not modified */
552 KLU_numeric<Entry, Int> *Numeric, /* input, not modified */
553 KLU_common<Entry, Int> *Common /* result in Common->rcond */
554)
555{
556 double ukk, umin = 0, umax = 0 ;
557 Entry *Udiag ;
558 Int j, n ;
559
560 /* ---------------------------------------------------------------------- */
561 /* check inputs */
562 /* ---------------------------------------------------------------------- */
563
564 if (Common == NULL)
565 {
566 return (FALSE) ;
567 }
568 if (Symbolic == NULL)
569 {
570 Common->status = KLU_INVALID ;
571 return (FALSE) ;
572 }
573 if (Numeric == NULL)
574 {
575 Common->rcond = 0 ;
576 Common->status = KLU_SINGULAR ;
577 return (TRUE) ;
578 }
579 Common->status = KLU_OK ;
580
581 /* ---------------------------------------------------------------------- */
582 /* compute rcond */
583 /* ---------------------------------------------------------------------- */
584
585 n = Symbolic->n ;
586 Udiag = (Entry *) Numeric->Udiag ;
587 for (j = 0 ; j < n ; j++)
588 {
589 /* get the magnitude of the pivot */
590 KLU2_ABS (ukk, Udiag [j]) ;
591 if (SCALAR_IS_NAN (ukk) || SCALAR_IS_ZERO (ukk))
592 {
593 /* if NaN, or zero, the rcond is zero */
594 Common->rcond = 0 ;
595 Common->status = KLU_SINGULAR ;
596 return (TRUE) ;
597 }
598 if (j == 0)
599 {
600 /* first pivot entry */
601 umin = ukk ;
602 umax = ukk ;
603 }
604 else
605 {
606 /* subsequent pivots */
607 umin = MIN (umin, ukk) ;
608 umax = MAX (umax, ukk) ;
609 }
610 }
611
612 Common->rcond = umin / umax ;
613 if (SCALAR_IS_NAN (Common->rcond) || SCALAR_IS_ZERO (Common->rcond))
614 {
615 /* this can occur if umin or umax are Inf or NaN */
616 Common->rcond = 0 ;
617 Common->status = KLU_SINGULAR ;
618 }
619 return (TRUE) ;
620}
621
622#endif