1 /* ========================================================================== */
2 /* === Cholesky/cholmod_rcond =============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Return a rough estimate of the reciprocal of the condition number:
10 * the minimum entry on the diagonal of L (or absolute entry of D for an LDL'
11 * factorization) divided by the maximum entry (squared for LL'). L can be
12 * real, complex, or zomplex. Returns -1 on error, 0 if the matrix is singular
13 * or has a zero entry on the diagonal of L, 1 if the matrix is 0-by-0, or
14 * min(diag(L))/max(diag(L)) otherwise. Never returns NaN; if L has a NaN on
15 * the diagonal it returns zero instead.
16 *
17 * For an LL' factorization, (min(diag(L))/max(diag(L)))^2 is returned.
18 * For an LDL' factorization, (min(diag(D))/max(diag(D))) is returned.
19 */
20
21 #ifndef NCHOLESKY
22
23 #include "cholmod_internal.h"
24 #include "cholmod_cholesky.h"
25
26 /* ========================================================================== */
27 /* === LMINMAX ============================================================== */
28 /* ========================================================================== */
29
30 /* Update lmin and lmax for one entry L(j,j) */
31
32 #define FIRST_LMINMAX(Ljj,lmin,lmax) \
33 { \
34 double ljj = Ljj ; \
35 if (IS_NAN (ljj)) \
36 { \
37 return (0) ; \
38 } \
39 lmin = ljj ; \
40 lmax = ljj ; \
41 }
42
43 #define LMINMAX(Ljj,lmin,lmax) \
44 { \
45 double ljj = Ljj ; \
46 if (IS_NAN (ljj)) \
47 { \
48 return (0) ; \
49 } \
50 if (ljj < lmin) \
51 { \
52 lmin = ljj ; \
53 } \
54 else if (ljj > lmax) \
55 { \
56 lmax = ljj ; \
57 } \
58 }
59
60 /* ========================================================================== */
61 /* === cholmod_rcond ======================================================== */
62 /* ========================================================================== */
63
CHOLMOD(rcond)64 double CHOLMOD(rcond) /* return min(diag(L)) / max(diag(L)) */
65 (
66 /* ---- input ---- */
67 cholmod_factor *L,
68 /* --------------- */
69 cholmod_common *Common
70 )
71 {
72 double lmin, lmax, rcond ;
73 double *Lx ;
74 Int *Lpi, *Lpx, *Super, *Lp ;
75 Int n, e, nsuper, s, k1, k2, psi, psend, psx, nsrow, nscol, jj, j ;
76
77 /* ---------------------------------------------------------------------- */
78 /* check inputs */
79 /* ---------------------------------------------------------------------- */
80
81 RETURN_IF_NULL_COMMON (EMPTY) ;
82 RETURN_IF_NULL (L, EMPTY) ;
83 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, EMPTY) ;
84 Common->status = CHOLMOD_OK ;
85
86 /* ---------------------------------------------------------------------- */
87 /* get inputs */
88 /* ---------------------------------------------------------------------- */
89
90 n = L->n ;
91 if (n == 0)
92 {
93 return (1) ;
94 }
95 if (L->minor < L->n)
96 {
97 return (0) ;
98 }
99
100 e = (L->xtype == CHOLMOD_COMPLEX) ? 2 : 1 ;
101
102 if (L->is_super)
103 {
104 /* L is supernodal */
105 nsuper = L->nsuper ; /* number of supernodes in L */
106 Lpi = L->pi ; /* column pointers for integer pattern */
107 Lpx = L->px ; /* column pointers for numeric values */
108 Super = L->super ; /* supernode sizes */
109 Lx = L->x ; /* numeric values */
110 FIRST_LMINMAX (Lx [0], lmin, lmax) ; /* first diagonal entry of L */
111 for (s = 0 ; s < nsuper ; s++)
112 {
113 k1 = Super [s] ; /* first column in supernode s */
114 k2 = Super [s+1] ; /* last column in supernode is k2-1 */
115 psi = Lpi [s] ; /* first row index is L->s [psi] */
116 psend = Lpi [s+1] ; /* last row index is L->s [psend-1] */
117 psx = Lpx [s] ; /* first numeric entry is Lx [psx] */
118 nsrow = psend - psi ; /* supernode is nsrow-by-nscol */
119 nscol = k2 - k1 ;
120 for (jj = 0 ; jj < nscol ; jj++)
121 {
122 LMINMAX (Lx [e * (psx + jj + jj*nsrow)], lmin, lmax) ;
123 }
124 }
125 }
126 else
127 {
128 /* L is simplicial */
129 Lp = L->p ;
130 Lx = L->x ;
131 if (L->is_ll)
132 {
133 /* LL' factorization */
134 FIRST_LMINMAX (Lx [Lp [0]], lmin, lmax) ;
135 for (j = 1 ; j < n ; j++)
136 {
137 LMINMAX (Lx [e * Lp [j]], lmin, lmax) ;
138 }
139 }
140 else
141 {
142 /* LDL' factorization, the diagonal might be negative */
143 FIRST_LMINMAX (fabs (Lx [Lp [0]]), lmin, lmax) ;
144 for (j = 1 ; j < n ; j++)
145 {
146 LMINMAX (fabs (Lx [e * Lp [j]]), lmin, lmax) ;
147 }
148 }
149 }
150 rcond = lmin / lmax ;
151 if (L->is_ll)
152 {
153 rcond = rcond*rcond ;
154 }
155 return (rcond) ;
156 }
157 #endif
158