Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_rcond.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_rcond =============================================== */
00003 /* ========================================================================== */
00004
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
00007  * The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
00008  * Lesser General Public License.  See lesser.txt for a text of the license.
00009  * CHOLMOD is also available under other licenses; contact authors for details.
00010  * http://www.cise.ufl.edu/research/sparse
00011  * -------------------------------------------------------------------------- */
00012
00013 /* Return a rough estimate of the reciprocal of the condition number:
00014  * the minimum entry on the diagonal of L (or absolute entry of D for an LDL'
00015  * factorization) divided by the maximum entry (squared for LL').  L can be
00016  * real, complex, or zomplex.  Returns -1 on error, 0 if the matrix is singular
00017  * or has a zero entry on the diagonal of L, 1 if the matrix is 0-by-0, or
00018  * min(diag(L))/max(diag(L)) otherwise.  Never returns NaN; if L has a NaN on
00019  * the diagonal it returns zero instead.
00020  *
00021  * For an LL' factorization,  (min(diag(L))/max(diag(L)))^2 is returned.
00022  * For an LDL' factorization, (min(diag(D))/max(diag(D))) is returned.
00023  */
00024
00025 #ifndef NCHOLESKY
00026
00027 #include "amesos_cholmod_internal.h"
00028 #include "amesos_cholmod_cholesky.h"
00029
00030 /* ========================================================================== */
00031 /* === LMINMAX ============================================================== */
00032 /* ========================================================================== */
00033
00034 /* Update lmin and lmax for one entry L(j,j) */
00035
00036 #define FIRST_LMINMAX(Ljj,lmin,lmax) \
00037 { \
00038     double ljj = Ljj ; \
00039     if (IS_NAN (ljj)) \
00040     { \
00041   return (0) ; \
00042     } \
00043     lmin = ljj ; \
00044     lmax = ljj ; \
00045 }
00046
00047 #define LMINMAX(Ljj,lmin,lmax) \
00048 { \
00049     double ljj = Ljj ; \
00050     if (IS_NAN (ljj)) \
00051     { \
00052   return (0) ; \
00053     } \
00054     if (ljj < lmin) \
00055     { \
00056   lmin = ljj ; \
00057     } \
00058     else if (ljj > lmax) \
00059     { \
00060   lmax = ljj ; \
00061     } \
00062 }
00063
00064 /* ========================================================================== */
00065 /* === cholmod_rcond ======================================================== */
00066 /* ========================================================================== */
00067
00068 double CHOLMOD(rcond)     /* return min(diag(L)) / max(diag(L)) */
00069 (
00070     /* ---- input ---- */
00071     cholmod_factor *L,
00072     /* --------------- */
00073     cholmod_common *Common
00074 )
00075 {
00076     double lmin, lmax, rcond ;
00077     double *Lx ;
00078     Int *Lpi, *Lpx, *Super, *Lp ;
00079     Int n, e, nsuper, s, k1, k2, psi, psend, psx, nsrow, nscol, jj, j ;
00080
00081     /* ---------------------------------------------------------------------- */
00082     /* check inputs */
00083     /* ---------------------------------------------------------------------- */
00084
00085     RETURN_IF_NULL_COMMON (EMPTY) ;
00086     RETURN_IF_NULL (L, EMPTY) ;
00087     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, EMPTY) ;
00088     Common->status = CHOLMOD_OK ;
00089
00090     /* ---------------------------------------------------------------------- */
00091     /* get inputs */
00092     /* ---------------------------------------------------------------------- */
00093
00094     n = L->n ;
00095     if (n == 0)
00096     {
00097   return (1) ;
00098     }
00099     if (L->minor < L->n)
00100     {
00101   return (0) ;
00102     }
00103
00104     e = (L->xtype == CHOLMOD_COMPLEX) ? 2 : 1 ;
00105
00106     if (L->is_super)
00107     {
00108   /* L is supernodal */
00109   nsuper = L->nsuper ;  /* number of supernodes in L */
00110   Lpi = L->pi ;   /* column pointers for integer pattern */
00111   Lpx = L->px ;   /* column pointers for numeric values */
00112   Super = L->super ;  /* supernode sizes */
00113   Lx = L->x ;   /* numeric values */
00114   FIRST_LMINMAX (Lx [0], lmin, lmax) ;  /* first diagonal entry of L */
00115   for (s = 0 ; s < nsuper ; s++)
00116   {
00117       k1 = Super [s] ;    /* first column in supernode s */
00118       k2 = Super [s+1] ;    /* last column in supernode is k2-1 */
00119       psi = Lpi [s] ;   /* first row index is L->s [psi] */
00120       psend = Lpi [s+1] ;   /* last row index is L->s [psend-1] */
00121       psx = Lpx [s] ;   /* first numeric entry is Lx [psx] */
00122       nsrow = psend - psi ; /* supernode is nsrow-by-nscol */
00123       nscol = k2 - k1 ;
00124       for (jj = 0 ; jj < nscol ; jj++)
00125       {
00126     LMINMAX (Lx [e * (psx + jj + jj*nsrow)], lmin, lmax) ;
00127       }
00128   }
00129     }
00130     else
00131     {
00132   /* L is simplicial */
00133   Lp = L->p ;
00134   Lx = L->x ;
00135   if (L->is_ll)
00136   {
00137       /* LL' factorization */
00138       FIRST_LMINMAX (Lx [Lp [0]], lmin, lmax) ;
00139       for (j = 1 ; j < n ; j++)
00140       {
00141     LMINMAX (Lx [e * Lp [j]], lmin, lmax) ;
00142       }
00143   }
00144   else
00145   {
00146       /* LDL' factorization, the diagonal might be negative */
00147       FIRST_LMINMAX (fabs (Lx [Lp [0]]), lmin, lmax) ;
00148       for (j = 1 ; j < n ; j++)
00149       {
00150     LMINMAX (fabs (Lx [e * Lp [j]]), lmin, lmax) ;
00151       }
00152   }
00153     }
00154     rcond = lmin / lmax ;
00155     if (L->is_ll)
00156     {
00157   rcond = rcond*rcond ;
00158     }
00159     return (rcond) ;
00160 }
00161 #endif
```