1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/ldlchol mexFunction =================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MATLAB Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * MATLAB(tm) is a Trademark of The MathWorks, Inc.
9 * -------------------------------------------------------------------------- */
10
11 /* Numeric LDL' factorization. Note that LL' and LDL' are faster than R'R
12 * and use less memory. The LDL' factorization methods use tril(A).
13 * The unit diagonal L can be obtained with tril(LD,-1)+speye(n) and the
14 * diagonal D can be obtained with D = diag(diag(LD)) ;
15 *
16 * LD = ldlchol (A) return the LDL' factorization of A
17 * [LD,p] = ldlchol (A) like [R,p] = chol(A), except LD is always square
18 * [LD,p,q] = ldlchol (A) factorizes A(q,q) into L*D*L'
19 *
20 * LD = ldlchol (A,beta) return the LDL' factorization of A*A'+beta*I
21 * [LD,p] = ldlchol (A,beta) like [R,p] = chol(A*A'+beta+I)
22 * [LD,p,q] = ldlchol (A,beta) factorizes A(q,:)*A(q,:)'+beta*I into L*D*L'
23 *
24 * Explicit zeros may appear in the LD matrix. The pattern of LD matches the
25 * pattern of L as computed by symbfact2, even if some entries in LD are
26 * explicitly zero. This is to ensure that ldlupdate works properly. You must
27 * NOT modify LD in MATLAB itself and then use ldlupdate if LD contains
28 * explicit zero entries; ldlupdate will fail catastrophically in this case.
29 *
30 * You MAY modify LD in MATLAB if you do not pass it back to ldlupdate. Just be
31 * aware that LD contains explicit zero entries, contrary to the standard
32 * practice in MATLAB of removing those entries from all sparse matrices.
33 *
34 * Note that CHOLMOD uses a supernodal LL' factorization and then converts it
35 * to LDL' for large matrices, and a simplicial LDL' factorization otherwise.
36 * Two-by-two block pivoting is not performed, in any case. Thus, ldlchol
37 * will not be able to perform an LDL' factorization of an arbitrary indefinite
38 * matrix. The matrix
39 *
40 * 0 1
41 * 1 0
42 *
43 * will fail, for example. You can tell CHOLMOD to always use its simplicial
44 * LDL' factorization by adding the statement
45 *
46 * cm->supernodal = CHOLMOD_SIMPLICIAL ;
47 *
48 * but ldlchol will be much slower for large matrices. It still will not be
49 * able to handle the above matrix, but it can then handle matrices such as
50 *
51 * -2 1
52 * 1 -2
53 *
54 * or any other symmetric matrix for which all leading minors are
55 * well-conditioned.
56 */
57
58 #include "cholmod_matlab.h"
59
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])60 void mexFunction
61 (
62 int nargout,
63 mxArray *pargout [ ],
64 int nargin,
65 const mxArray *pargin [ ]
66 )
67 {
68 double dummy = 0, beta [2], *px ;
69 cholmod_sparse Amatrix, *A, *Lsparse ;
70 cholmod_factor *L ;
71 cholmod_common Common, *cm ;
72 Long n, minor ;
73
74 /* ---------------------------------------------------------------------- */
75 /* start CHOLMOD and set parameters */
76 /* ---------------------------------------------------------------------- */
77
78 cm = &Common ;
79 cholmod_l_start (cm) ;
80 sputil_config (SPUMONI, cm) ;
81
82 /* convert to packed LDL' when done */
83 cm->final_asis = FALSE ;
84 cm->final_super = FALSE ;
85 cm->final_ll = FALSE ;
86 cm->final_pack = TRUE ;
87 cm->final_monotonic = TRUE ;
88
89 /* since numerically zero entries are NOT dropped from the symbolic
90 * pattern, we DO need to drop entries that result from supernodal
91 * amalgamation. */
92 cm->final_resymbol = TRUE ;
93
94 cm->quick_return_if_not_posdef = (nargout < 2) ;
95
96 /* This will disable the supernodal LL', which will be slow. */
97 /* cm->supernodal = CHOLMOD_SIMPLICIAL ; */
98
99 /* ---------------------------------------------------------------------- */
100 /* get inputs */
101 /* ---------------------------------------------------------------------- */
102
103 if (nargin < 1 || nargin > 2 || nargout > 3)
104 {
105 mexErrMsgTxt ("usage: [L,p,q] = ldlchol (A,beta)") ;
106 }
107
108 n = mxGetM (pargin [0]) ;
109
110 if (!mxIsSparse (pargin [0]))
111 {
112 mexErrMsgTxt ("A must be sparse") ;
113 }
114 if (nargin == 1 && n != mxGetN (pargin [0]))
115 {
116 mexErrMsgTxt ("A must be square") ;
117 }
118
119 /* get sparse matrix A, use tril(A) */
120 A = sputil_get_sparse (pargin [0], &Amatrix, &dummy, -1) ;
121
122 if (nargin == 1)
123 {
124 A->stype = -1 ; /* use lower part of A */
125 beta [0] = 0 ;
126 beta [1] = 0 ;
127 }
128 else
129 {
130 A->stype = 0 ; /* use all of A, factorizing A*A' */
131 beta [0] = mxGetScalar (pargin [1]) ;
132 beta [1] = 0 ;
133 }
134
135 /* use natural ordering if no q output parameter */
136 if (nargout < 3)
137 {
138 cm->nmethods = 1 ;
139 cm->method [0].ordering = CHOLMOD_NATURAL ;
140 cm->postorder = FALSE ;
141 }
142
143 /* ---------------------------------------------------------------------- */
144 /* analyze and factorize */
145 /* ---------------------------------------------------------------------- */
146
147 L = cholmod_l_analyze (A, cm) ;
148 cholmod_l_factorize_p (A, beta, NULL, 0, L, cm) ;
149
150 if (nargout < 2 && cm->status != CHOLMOD_OK)
151 {
152 mexErrMsgTxt ("matrix is not positive definite") ;
153 }
154
155 /* ---------------------------------------------------------------------- */
156 /* convert L to a sparse matrix */
157 /* ---------------------------------------------------------------------- */
158
159 /* the conversion sets L->minor back to n, so get a copy of it first */
160 minor = L->minor ;
161 Lsparse = cholmod_l_factor_to_sparse (L, cm) ;
162 if (Lsparse->xtype == CHOLMOD_COMPLEX)
163 {
164 /* convert Lsparse from complex to zomplex */
165 cholmod_l_sparse_xtype (CHOLMOD_ZOMPLEX, Lsparse, cm) ;
166 }
167
168 /* ---------------------------------------------------------------------- */
169 /* return results to MATLAB */
170 /* ---------------------------------------------------------------------- */
171
172 /* return L as a sparse matrix (it may contain numerically zero entries) */
173 pargout [0] = sputil_put_sparse (&Lsparse, cm) ;
174
175 /* return minor (translate to MATLAB convention) */
176 if (nargout > 1)
177 {
178 pargout [1] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
179 px = mxGetPr (pargout [1]) ;
180 px [0] = ((minor == n) ? 0 : (minor+1)) ;
181 }
182
183 /* return permutation */
184 if (nargout > 2)
185 {
186 pargout [2] = sputil_put_int (L->Perm, n, 1) ;
187 }
188
189 /* ---------------------------------------------------------------------- */
190 /* free workspace and the CHOLMOD L, except for what is copied to MATLAB */
191 /* ---------------------------------------------------------------------- */
192
193 cholmod_l_free_factor (&L, cm) ;
194 cholmod_l_finish (cm) ;
195 cholmod_l_print_common (" ", cm) ;
196 /*
197 if (cm->malloc_count != 3 + mxIsComplex (pargout[0])) mexErrMsgTxt ("!") ;
198 */
199 }
200