1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/cholmod 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 /* Supernodal sparse Cholesky backslash, x = A\b. Factorizes PAP' in LL' then
12 * solves a sparse linear system. Uses the diagonal and upper triangular part
13 * of A only. A must be sparse. b can be sparse or dense.
14 *
15 * Usage:
16 *
17 * x = cholmod2 (A, b)
18 * [x stats] = cholmod2 (A, b, ordering) % a scalar: 0,-1,-2, or -3
19 * [x stats] = cholmod2 (A, b, p) % a permutation vector
20 *
21 * The 3rd argument select the ordering method to use. If not present or -1,
22 * the default ordering strategy is used (AMD, and then try METIS if AMD finds
23 * an ordering with high fill-in, and use the best method tried).
24 *
25 * Other options for the ordering parameter:
26 *
27 * 0 natural (no etree postordering)
28 * -1 use CHOLMOD's default ordering strategy (AMD, then try METIS)
29 * -2 AMD, and then try NESDIS (not METIS) if AMD has high fill-in
30 * -3 use AMD only
31 * -4 use METIS only
32 * -5 use NESDIS only
33 * -6 natural, but with etree postordering
34 * p user permutation (vector of size n, with a permutation of 1:n)
35 *
36 * stats(1) estimate of the reciprocal of the condition number
37 * stats(2) ordering used:
38 * 0: natural, 1: given, 2:amd, 3:metis, 4:nesdis, 5:colamd,
39 * 6: natural but postordered.
40 * stats(3) nnz(L)
41 * stats(4) flop count in Cholesky factorization. Excludes solution
42 * of upper/lower triangular systems, which can be easily
43 * computed from stats(3) (roughly 4*nnz(L)*size(b,2)).
44 * stats(5) memory usage in MB.
45 */
46
47 #include "cholmod_matlab.h"
48
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])49 void mexFunction
50 (
51 int nargout,
52 mxArray *pargout [ ],
53 int nargin,
54 const mxArray *pargin [ ]
55 )
56 {
57 double dummy = 0, rcond, *p ;
58 cholmod_sparse Amatrix, Bspmatrix, *A, *Bs, *Xs ;
59 cholmod_dense Bmatrix, *X, *B ;
60 cholmod_factor *L ;
61 cholmod_common Common, *cm ;
62 Long n, B_is_sparse, ordering, k, *Perm ;
63
64 /* ---------------------------------------------------------------------- */
65 /* start CHOLMOD and set parameters */
66 /* ---------------------------------------------------------------------- */
67
68 cm = &Common ;
69 cholmod_l_start (cm) ;
70 sputil_config (SPUMONI, cm) ;
71
72 /* There is no supernodal LDL'. If cm->final_ll = FALSE (the default), then
73 * this mexFunction will use a simplicial LDL' when flops/lnz < 40, and a
74 * supernodal LL' otherwise. This may give suprising results to the MATLAB
75 * user, so always perform an LL' factorization by setting cm->final_ll
76 * to TRUE. */
77
78 cm->final_ll = TRUE ;
79 cm->quick_return_if_not_posdef = TRUE ;
80
81 /* ---------------------------------------------------------------------- */
82 /* get inputs */
83 /* ---------------------------------------------------------------------- */
84
85 if (nargout > 2 || nargin < 2 || nargin > 3)
86 {
87 mexErrMsgTxt ("usage: [x,rcond] = cholmod2 (A,b,ordering)") ;
88 }
89 n = mxGetM (pargin [0]) ;
90 if (!mxIsSparse (pargin [0]) || (n != mxGetN (pargin [0])))
91 {
92 mexErrMsgTxt ("A must be square and sparse") ;
93 }
94 if (n != mxGetM (pargin [1]))
95 {
96 mexErrMsgTxt ("# of rows of A and B must match") ;
97 }
98
99 /* get sparse matrix A. Use triu(A) only. */
100 A = sputil_get_sparse (pargin [0], &Amatrix, &dummy, 1) ;
101
102 /* get sparse or dense matrix B */
103 B = NULL ;
104 Bs = NULL ;
105 B_is_sparse = mxIsSparse (pargin [1]) ;
106 if (B_is_sparse)
107 {
108 /* get sparse matrix B (unsymmetric) */
109 Bs = sputil_get_sparse (pargin [1], &Bspmatrix, &dummy, 0) ;
110 }
111 else
112 {
113 /* get dense matrix B */
114 B = sputil_get_dense (pargin [1], &Bmatrix, &dummy) ;
115 }
116
117 /* get the ordering option */
118 if (nargin < 3)
119 {
120 /* use default ordering */
121 ordering = -1 ;
122 }
123 else
124 {
125 /* use a non-default option */
126 ordering = mxGetScalar (pargin [2]) ;
127 }
128
129 p = NULL ;
130 Perm = NULL ;
131
132 if (ordering == 0)
133 {
134 /* natural ordering */
135 cm->nmethods = 1 ;
136 cm->method [0].ordering = CHOLMOD_NATURAL ;
137 cm->postorder = FALSE ;
138 }
139 else if (ordering == -1)
140 {
141 /* default strategy ... nothing to change */
142 }
143 else if (ordering == -2)
144 {
145 /* default strategy, but with NESDIS in place of METIS */
146 cm->default_nesdis = TRUE ;
147 }
148 else if (ordering == -3)
149 {
150 /* use AMD only */
151 cm->nmethods = 1 ;
152 cm->method [0].ordering = CHOLMOD_AMD ;
153 cm->postorder = TRUE ;
154 }
155 else if (ordering == -4)
156 {
157 /* use METIS only */
158 cm->nmethods = 1 ;
159 cm->method [0].ordering = CHOLMOD_METIS ;
160 cm->postorder = TRUE ;
161 }
162 else if (ordering == -5)
163 {
164 /* use NESDIS only */
165 cm->nmethods = 1 ;
166 cm->method [0].ordering = CHOLMOD_NESDIS ;
167 cm->postorder = TRUE ;
168 }
169 else if (ordering == -6)
170 {
171 /* natural ordering, but with etree postordering */
172 cm->nmethods = 1 ;
173 cm->method [0].ordering = CHOLMOD_NATURAL ;
174 cm->postorder = TRUE ;
175 }
176 else if (ordering == -7)
177 {
178 /* always try both AMD and METIS, and pick the best */
179 cm->nmethods = 2 ;
180 cm->method [0].ordering = CHOLMOD_AMD ;
181 cm->method [1].ordering = CHOLMOD_METIS ;
182 cm->postorder = TRUE ;
183 }
184 else if (ordering >= 1)
185 {
186 /* assume the 3rd argument is a user-provided permutation of 1:n */
187 if (mxGetNumberOfElements (pargin [2]) != n)
188 {
189 mexErrMsgTxt ("invalid input permutation") ;
190 }
191 /* copy from double to integer, and convert to 0-based */
192 p = mxGetPr (pargin [2]) ;
193 Perm = cholmod_l_malloc (n, sizeof (Long), cm) ;
194 for (k = 0 ; k < n ; k++)
195 {
196 Perm [k] = p [k] - 1 ;
197 }
198 /* check the permutation */
199 if (!cholmod_l_check_perm (Perm, n, n, cm))
200 {
201 mexErrMsgTxt ("invalid input permutation") ;
202 }
203 /* use only the given permutation */
204 cm->nmethods = 1 ;
205 cm->method [0].ordering = CHOLMOD_GIVEN ;
206 cm->postorder = FALSE ;
207 }
208 else
209 {
210 mexErrMsgTxt ("invalid ordering option") ;
211 }
212
213 /* ---------------------------------------------------------------------- */
214 /* analyze and factorize */
215 /* ---------------------------------------------------------------------- */
216
217 L = cholmod_l_analyze_p (A, Perm, NULL, 0, cm) ;
218 cholmod_l_free (n, sizeof (Long), Perm, cm) ;
219 cholmod_l_factorize (A, L, cm) ;
220
221 rcond = cholmod_l_rcond (L, cm) ;
222
223 if (rcond == 0)
224 {
225 mexWarnMsgTxt ("Matrix is indefinite or singular to working precision");
226 }
227 else if (rcond < DBL_EPSILON)
228 {
229 mexWarnMsgTxt ("Matrix is close to singular or badly scaled.") ;
230 mexPrintf (" Results may be inaccurate. RCOND = %g.\n", rcond) ;
231 }
232
233 /* ---------------------------------------------------------------------- */
234 /* solve and return solution to MATLAB */
235 /* ---------------------------------------------------------------------- */
236
237 if (B_is_sparse)
238 {
239 /* solve AX=B with sparse X and B; return sparse X to MATLAB */
240 Xs = cholmod_l_spsolve (CHOLMOD_A, L, Bs, cm) ;
241 pargout [0] = sputil_put_sparse (&Xs, cm) ;
242 }
243 else
244 {
245 /* solve AX=B with dense X and B; return dense X to MATLAB */
246 X = cholmod_l_solve (CHOLMOD_A, L, B, cm) ;
247 pargout [0] = sputil_put_dense (&X, cm) ;
248 }
249
250 /* return statistics, if requested */
251 if (nargout > 1)
252 {
253 pargout [1] = mxCreateDoubleMatrix (1, 5, mxREAL) ;
254 p = mxGetPr (pargout [1]) ;
255 p [0] = rcond ;
256 p [1] = L->ordering ;
257 p [2] = cm->lnz ;
258 p [3] = cm->fl ;
259 p [4] = cm->memory_usage / 1048576. ;
260 }
261
262 cholmod_l_free_factor (&L, cm) ;
263 cholmod_l_finish (cm) ;
264 cholmod_l_print_common (" ", cm) ;
265 /*
266 if (cm->malloc_count !=
267 (mxIsComplex (pargout [0]) + (mxIsSparse (pargout[0]) ? 3:1)))
268 mexErrMsgTxt ("memory leak!") ;
269 */
270 }
271