Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_diagnostics.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_diagnostics ====================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Linear algebraic diagnostics:
00006  * KLU_rgrowth: reciprocal pivot growth, takes O(|A|+|U|) time
00007  * KLU_condest: condition number estimator, takes about O(|A|+5*(|L|+|U|)) time
00008  * KLU_flops: compute # flops required to factorize A into L*U
00009  * KLU_rcond: compute a really cheap estimate of the reciprocal of the
00010  *    condition number, min(abs(diag(U))) / max(abs(diag(U))).
00011  *    Takes O(n) time.
00012  */
00013 
00014 #include "amesos_klu_internal.h"
00015 
00016 /* ========================================================================== */
00017 /* === KLU_rgrowth ========================================================== */
00018 /* ========================================================================== */
00019 
00020 /* Compute the reciprocal pivot growth factor.  In MATLAB notation:
00021  *
00022  *   rgrowth = min (max (abs ((R \ A (p,q)) - F))) ./ max (abs (U)))
00023  */
00024 
00025 Int KLU_rgrowth   /* return TRUE if successful, FALSE otherwise */
00026 (
00027     Int *Ap,
00028     Int *Ai,
00029     double *Ax,
00030     KLU_symbolic *Symbolic,
00031     KLU_numeric *Numeric,
00032     KLU_common *Common
00033 )
00034 {
00035     double temp, max_ai, max_ui, min_block_rgrowth ;
00036     Entry aik ;
00037     Int *Q, *Ui, *Uip, *Ulen, *Pinv ;
00038     Unit *LU ;
00039     Entry *Aentry, *Ux, *Ukk ;
00040     double *Rs ;
00041     Int i, newrow, oldrow, k1, k2, nk, j, oldcol, k, pend, len ;
00042 
00043     /* ---------------------------------------------------------------------- */
00044     /* check inputs */
00045     /* ---------------------------------------------------------------------- */
00046 
00047     if (Common == NULL)
00048     {
00049   return (FALSE) ;
00050     }
00051 
00052     if (Symbolic == NULL || Ap == NULL || Ai == NULL || Ax == NULL)
00053     {
00054   Common->status = KLU_INVALID ;
00055   return (FALSE) ;
00056     }
00057 
00058     if (Numeric == NULL)
00059     {
00060   /* treat this as a singular matrix */
00061   Common->rgrowth = 0 ;
00062   Common->status = KLU_SINGULAR ;
00063   return (TRUE) ;
00064     }
00065     Common->status = KLU_OK ;
00066 
00067     /* ---------------------------------------------------------------------- */
00068     /* compute the reciprocal pivot growth */
00069     /* ---------------------------------------------------------------------- */
00070 
00071     Aentry = (Entry *) Ax ;
00072     Pinv = Numeric->Pinv ;
00073     Rs = Numeric->Rs ;
00074     Q = Symbolic->Q ;
00075     Common->rgrowth = 1 ;
00076 
00077     for (i = 0 ; i < Symbolic->nblocks ; i++)
00078     {
00079   k1 = Symbolic->R[i] ;
00080   k2 = Symbolic->R[i+1] ;
00081   nk = k2 - k1 ;
00082   if (nk == 1)
00083   {
00084       continue ;      /* skip singleton blocks */
00085   }
00086   LU = (Unit *) Numeric->LUbx[i] ;
00087   Uip = Numeric->Uip + k1 ;
00088   Ulen = Numeric->Ulen + k1 ;
00089   Ukk = ((Entry *) Numeric->Udiag) + k1 ;
00090   min_block_rgrowth = 1 ;
00091   for (j = 0 ; j < nk ; j++)
00092   {
00093       max_ai = 0 ;
00094       max_ui = 0 ;
00095       oldcol = Q[j + k1] ;
00096       pend = Ap [oldcol + 1] ;
00097       for (k = Ap [oldcol] ; k < pend ; k++)
00098       {
00099     oldrow = Ai [k] ;
00100     newrow = Pinv [oldrow] ;
00101                 if (newrow < k1)
00102     {
00103         continue ;  /* skip entry outside the block */
00104     }
00105     ASSERT (newrow < k2) ;
00106     if (Rs != NULL)
00107     {
00108         /* aik = Aentry [k] / Rs [oldrow] */
00109         SCALE_DIV_ASSIGN (aik, Aentry [k], Rs [newrow]) ;
00110     }
00111     else
00112     {
00113         aik = Aentry [k] ;
00114     }
00115     /* temp = ABS (aik) */
00116     ABS (temp, aik) ;
00117     if (temp > max_ai)
00118     {
00119         max_ai = temp ;
00120     }
00121       }
00122 
00123       GET_POINTER (LU, Uip, Ulen, Ui, Ux, j, len) ;
00124       for (k = 0 ; k < len ; k++)
00125       {
00126           /* temp = ABS (Ux [k]) */
00127           ABS (temp, Ux [k]) ;
00128     if (temp > max_ui)
00129     {
00130         max_ui = temp ;
00131     }
00132       }
00133       /* consider the diagonal element */
00134       ABS (temp, Ukk [j]) ;
00135       if (temp > max_ui)
00136       {
00137     max_ui = temp ;
00138       }
00139 
00140       /* if max_ui is 0, skip the column */
00141       if (SCALAR_IS_ZERO (max_ui))
00142       {
00143     continue ;
00144       }
00145       temp = max_ai / max_ui ;
00146       if (temp < min_block_rgrowth)
00147       {
00148     min_block_rgrowth = temp ;
00149       }
00150   }
00151 
00152   if (min_block_rgrowth < Common->rgrowth)
00153   {
00154       Common->rgrowth = min_block_rgrowth ;
00155   }
00156     }
00157     return (TRUE) ;
00158 }
00159 
00160 
00161 /* ========================================================================== */
00162 /* === KLU_condest ========================================================== */
00163 /* ========================================================================== */
00164 
00165 /* Estimate the condition number.  Uses Higham and Tisseur's algorithm
00166  * (A block algorithm for matrix 1-norm estimation, with applications to
00167  * 1-norm pseudospectra, SIAM J. Matrix Anal. Appl., 21(4):1185-1201, 2000.
00168  */
00169 
00170 Int KLU_condest   /* return TRUE if successful, FALSE otherwise */
00171 (
00172     Int Ap [ ],
00173     double Ax [ ],
00174     KLU_symbolic *Symbolic,
00175     KLU_numeric *Numeric,
00176     KLU_common *Common
00177 )
00178 {
00179     double xj, Xmax, csum, anorm, ainv_norm, est_old, est_new, abs_value ;
00180     Entry *Udiag, *Aentry, *X, *S ;
00181     Int *R ;
00182     Int nblocks, i, j, jmax, jnew, pend, n ;
00183 #ifndef COMPLEX
00184     Int unchanged ;
00185 #endif
00186 
00187     /* ---------------------------------------------------------------------- */
00188     /* check inputs */
00189     /* ---------------------------------------------------------------------- */
00190 
00191     if (Common == NULL)
00192     {
00193   return (FALSE) ;
00194     }
00195     if (Symbolic == NULL || Ap == NULL || Ax == NULL)
00196     {
00197   Common->status = KLU_INVALID ;
00198   return (FALSE) ;
00199     }
00200     abs_value = 0 ;
00201     if (Numeric == NULL)
00202     {
00203   /* treat this as a singular matrix */
00204   Common->condest = 1 / abs_value ;
00205   Common->status = KLU_SINGULAR ;
00206   return (TRUE) ;
00207     }
00208     Common->status = KLU_OK ;
00209 
00210     /* ---------------------------------------------------------------------- */
00211     /* get inputs */
00212     /* ---------------------------------------------------------------------- */
00213 
00214     n = Symbolic->n ;
00215     nblocks = Symbolic->nblocks ;
00216     R = Symbolic->R ;
00217     Udiag = Numeric->Udiag ;
00218 
00219     /* ---------------------------------------------------------------------- */
00220     /* check if diagonal of U has a zero on it */
00221     /* ---------------------------------------------------------------------- */
00222 
00223     for (i = 0 ; i < n ; i++)
00224     {
00225   ABS (abs_value, Udiag [i]) ;
00226   if (SCALAR_IS_ZERO (abs_value))
00227   {
00228       Common->condest = 1 / abs_value ;
00229       Common->status = KLU_SINGULAR ;
00230       return (TRUE) ;
00231   }
00232     }
00233 
00234     /* ---------------------------------------------------------------------- */
00235     /* compute 1-norm (maximum column sum) of the matrix */
00236     /* ---------------------------------------------------------------------- */
00237 
00238     anorm =  0.0 ;
00239     Aentry = (Entry *) Ax ;
00240     for (i = 0 ; i < n ; i++)
00241     {
00242   pend = Ap [i + 1] ;
00243   csum = 0.0 ;
00244   for (j = Ap [i] ; j < pend ; j++)
00245   {
00246       ABS (abs_value, Aentry [j]) ;
00247       csum += abs_value ;
00248   }
00249   if (csum > anorm)
00250   {
00251       anorm = csum ;
00252   }
00253     }
00254 
00255     /* ---------------------------------------------------------------------- */
00256     /* compute estimate of 1-norm of inv (A) */
00257     /* ---------------------------------------------------------------------- */
00258 
00259     /* get workspace (size 2*n Entry's) */
00260     X = Numeric->Xwork ;      /* size n space used in KLU_solve, tsolve */
00261     X += n ;          /* X is size n */
00262     S = X + n ;         /* S is size n */
00263 
00264     for (i = 0 ; i < n ; i++)
00265     {
00266   CLEAR (S [i]) ;
00267   CLEAR (X [i]) ;
00268   REAL (X [i]) = 1.0 / ((double) n) ;
00269     }
00270     jmax = 0 ;
00271 
00272     ainv_norm = 0.0 ;
00273     for (i = 0 ; i < 5 ; i++)
00274     {
00275   if (i > 0)
00276   {
00277       /* X [jmax] is the largest entry in X */
00278       for (j = 0 ; j < n ; j++)
00279       {
00280     /* X [j] = 0 ;*/
00281     CLEAR (X [j]) ;
00282       }
00283       REAL (X [jmax]) = 1 ;
00284   }
00285 
00286   KLU_solve (Symbolic, Numeric, n, 1, (double *) X, Common) ;
00287   est_old = ainv_norm ;
00288   ainv_norm = 0.0 ;
00289 
00290   for (j = 0 ; j < n ; j++)
00291   {
00292       /* ainv_norm += ABS (X [j]) ;*/
00293       ABS (abs_value, X [j]) ;
00294       ainv_norm += abs_value ;
00295   }
00296 
00297 #ifndef COMPLEX
00298   unchanged = TRUE ;
00299 
00300   for (j = 0 ; j < n ; j++)
00301   {
00302       double s = (X [j] >= 0) ? 1 : -1 ;
00303       if (s != (Int) REAL (S [j]))
00304       {
00305     S [j] = s ;
00306     unchanged = FALSE ;
00307       }
00308   }
00309 
00310   if (i > 0 && (ainv_norm <= est_old || unchanged))
00311   {
00312       break ;
00313   }
00314 #else
00315         for (j = 0 ; j < n ; j++)
00316   {
00317       if (IS_NONZERO (X [j]))
00318       {
00319     ABS (abs_value, X [j]) ;
00320     SCALE_DIV_ASSIGN (S [j], X [j], abs_value) ;
00321       }
00322       else
00323       {
00324     CLEAR (S [j]) ;
00325     REAL (S [j]) = 1 ;
00326       }
00327   }
00328 
00329   if (i > 0 && ainv_norm <= est_old)
00330         {
00331             break ;
00332         }
00333 #endif
00334 
00335   for (j = 0 ; j < n ; j++)
00336   {
00337       X [j] = S [j] ;
00338   }
00339 
00340 #ifndef COMPLEX
00341   /* do a transpose solve */
00342   KLU_tsolve (Symbolic, Numeric, n, 1, X, Common) ;
00343 #else
00344   /* do a conjugate transpose solve */
00345   KLU_tsolve (Symbolic, Numeric, n, 1, (double *) X, 1, Common) ;
00346 #endif
00347 
00348   /* jnew = the position of the largest entry in X */
00349   jnew = 0 ;
00350   Xmax = 0 ;
00351   for (j = 0 ; j < n ; j++)
00352   {
00353       /* xj = ABS (X [j]) ;*/
00354       ABS (xj, X [j]) ;
00355       if (xj > Xmax)
00356       {
00357     Xmax = xj ;
00358     jnew = j ;
00359       }
00360   }
00361   if (i > 0 && jnew == jmax)
00362   {
00363       /* the position of the largest entry did not change
00364        * from the previous iteration */
00365       break ;
00366   }
00367   jmax = jnew ;
00368     }
00369 
00370     /* ---------------------------------------------------------------------- */
00371     /* compute another estimate of norm(inv(A),1), and take the largest one */
00372     /* ---------------------------------------------------------------------- */
00373 
00374     for (j = 0 ; j < n ; j++)
00375     {
00376   CLEAR (X [j]) ;
00377   if (j % 2)
00378   {
00379       REAL (X [j]) = 1 + ((double) j) / ((double) (n-1)) ;
00380   }
00381   else
00382   {
00383       REAL (X [j]) = -1 - ((double) j) / ((double) (n-1)) ;
00384   }
00385     }
00386 
00387     KLU_solve (Symbolic, Numeric, n, 1, (double *) X, Common) ;
00388 
00389     est_new = 0.0 ;
00390     for (j = 0 ; j < n ; j++)
00391     {
00392   /* est_new += ABS (X [j]) ;*/
00393   ABS (abs_value, X [j]) ;
00394   est_new += abs_value ;
00395     }
00396     est_new = 2 * est_new / (3 * n) ;
00397     ainv_norm = MAX (est_new, ainv_norm) ;
00398 
00399     /* ---------------------------------------------------------------------- */
00400     /* compute estimate of condition number */
00401     /* ---------------------------------------------------------------------- */
00402 
00403     Common->condest = ainv_norm * anorm ;
00404     return (TRUE) ;
00405 }
00406 
00407 
00408 /* ========================================================================== */
00409 /* === KLU_flops ============================================================ */
00410 /* ========================================================================== */
00411 
00412 /* Compute the flop count for the LU factorization (in Common->flops) */
00413 
00414 Int KLU_flops   /* return TRUE if successful, FALSE otherwise */
00415 (
00416     KLU_symbolic *Symbolic,
00417     KLU_numeric *Numeric,
00418     KLU_common *Common
00419 )
00420 {
00421     double flops = 0 ;
00422     Int *R, *Ui, *Uip, *Llen, *Ulen ;
00423     Unit **LUbx ;
00424     Unit *LU ;
00425     Int k, ulen, p, n, nk, block, nblocks, k1 ;
00426 
00427     /* ---------------------------------------------------------------------- */
00428     /* check inputs */
00429     /* ---------------------------------------------------------------------- */
00430 
00431     if (Common == NULL)
00432     {
00433   return (FALSE) ;
00434     }
00435     Common->flops = EMPTY ;
00436     if (Numeric == NULL || Symbolic == NULL)
00437     {
00438   Common->status = KLU_INVALID ;
00439   return (FALSE) ;
00440     }
00441     Common->status = KLU_OK ;
00442 
00443     /* ---------------------------------------------------------------------- */
00444     /* get the contents of the Symbolic object */
00445     /* ---------------------------------------------------------------------- */
00446 
00447     n = Symbolic->n ;
00448     R = Symbolic->R ;
00449     nblocks = Symbolic->nblocks ;
00450 
00451     /* ---------------------------------------------------------------------- */
00452     /* get the contents of the Numeric object */
00453     /* ---------------------------------------------------------------------- */
00454 
00455     LUbx = (Unit **) Numeric->LUbx ;
00456 
00457     /* ---------------------------------------------------------------------- */
00458     /* compute the flop count */
00459     /* ---------------------------------------------------------------------- */
00460 
00461     for (block = 0 ; block < nblocks ; block++)
00462     {
00463   k1 = R [block] ;
00464   nk = R [block+1] - k1 ;
00465   if (nk > 1)
00466   {
00467       Llen = Numeric->Llen + k1 ;
00468       Uip  = Numeric->Uip  + k1 ;
00469       Ulen = Numeric->Ulen + k1 ;
00470       LU = LUbx [block] ;
00471       for (k = 0 ; k < nk ; k++)
00472       {
00473     /* compute kth column of U, and update kth column of A */
00474     GET_I_POINTER (LU, Uip, Ui, k) ;
00475     ulen = Ulen [k] ;
00476     for (p = 0 ; p < ulen ; p++)
00477     {
00478         flops += 2 * Llen [Ui [p]] ;
00479     }
00480     /* gather and divide by pivot to get kth column of L */
00481     flops += Llen [k] ;
00482       }
00483   }
00484     }
00485     Common->flops = flops ;
00486     return (TRUE) ;
00487 }
00488 
00489 
00490 /* ========================================================================== */
00491 /* === KLU_rcond ============================================================ */
00492 /* ========================================================================== */
00493 
00494 /* Compute a really cheap estimate of the reciprocal of the condition number,
00495  * condition number, min(abs(diag(U))) / max(abs(diag(U))).  If U has a zero
00496  * pivot, or a NaN pivot, rcond will be zero.  Takes O(n) time.
00497  */   
00498 
00499 Int KLU_rcond   /* return TRUE if successful, FALSE otherwise */
00500 (
00501     KLU_symbolic *Symbolic, /* input, not modified */
00502     KLU_numeric *Numeric, /* input, not modified */
00503     KLU_common *Common    /* result in Common->rcond */
00504 )
00505 {
00506     double ukk, umin, umax ;
00507     Entry *Udiag ;
00508     Int j, n ;
00509 
00510     /* ---------------------------------------------------------------------- */
00511     /* check inputs */
00512     /* ---------------------------------------------------------------------- */
00513 
00514     if (Common == NULL)
00515     {
00516   return (FALSE) ;
00517     }
00518     if (Symbolic == NULL)
00519     {
00520   Common->status = KLU_INVALID ;
00521   return (FALSE) ;
00522     }
00523     if (Numeric == NULL)
00524     {
00525   Common->rcond = 0 ;
00526   Common->status = KLU_SINGULAR ;
00527   return (TRUE) ;
00528     }
00529     Common->status = KLU_OK ;
00530 
00531     /* ---------------------------------------------------------------------- */
00532     /* compute rcond */
00533     /* ---------------------------------------------------------------------- */
00534 
00535     n = Symbolic->n ;
00536     Udiag = Numeric->Udiag ;
00537     for (j = 0 ; j < n ; j++)
00538     {
00539   /* get the magnitude of the pivot */
00540   ABS (ukk, Udiag [j]) ;
00541   if (SCALAR_IS_NAN (ukk) || SCALAR_IS_ZERO (ukk))
00542   {
00543       /* if NaN, or zero, the rcond is zero */
00544       Common->rcond = 0 ;
00545       Common->status = KLU_SINGULAR ;
00546       return (TRUE) ;
00547   }
00548   if (j == 0)
00549   {
00550       /* first pivot entry */
00551       umin = ukk ;
00552       umax = ukk ;
00553   }
00554   else
00555   {
00556       /* subsequent pivots */
00557       umin = MIN (umin, ukk) ;
00558       umax = MAX (umax, ukk) ;
00559   }
00560     }
00561 
00562     Common->rcond = umin / umax ;
00563     if (SCALAR_IS_NAN (Common->rcond) || SCALAR_IS_ZERO (Common->rcond))
00564     {
00565   /* this can occur if umin or umax are Inf or NaN */
00566   Common->rcond = 0 ;
00567   Common->status = KLU_SINGULAR ;
00568     }
00569     return (TRUE) ;
00570 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines