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 
61 template <typename Entry, typename Int>
62 Int 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 
207 template <typename Entry, typename Int>
208 Int 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 
462 template <typename Entry, typename Int>
463 Int 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 = 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 
548 template <typename Entry, typename Int>
549 Int 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