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