1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/symbfact2 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 /* Usage:
12 *
13 * count = symbfact2 (A) returns row counts of R=chol(A)
14 * count = symbfact2 (A,'col') returns row counts of R=chol(A'*A)
15 * count = symbfact2 (A,'sym') same as symbfact2(A)
16 * count = symbfact2 (A,'lo') same as symbfact2(A'), uses tril(A)
17 * count = symbfact2 (A,'row') returns row counts of R=chol(A*A')
18 *
19 * [count, h, parent, post, R] = symbfact2 (...)
20 *
21 * h: height of the elimination tree
22 * parent: the elimination tree itself
23 * post: postordering of the elimination tree
24 * R: a 0-1 matrix whose structure is that of chol(A) or chol(A'*A)
25 * for the 'col' case
26 *
27 * symbfact2(A) and symbfact2(A,'sym') uses triu(A).
28 * symbfact2(A,'lo') uses tril(A).
29 *
30 * These forms return L = R' instead of R. They are faster and take less
31 * memory. They return the same count, h, parent, and post outputs.
32 *
33 * [count, h, parent, post, L] = symbfact2 (A,'col','L')
34 * [count, h, parent, post, L] = symbfact2 (A,'sym','L')
35 * [count, h, parent, post, L] = symbfact2 (A,'lo', 'L')
36 * [count, h, parent, post, L] = symbfact2 (A,'row','L')
37 */
38
39 #include "cholmod_matlab.h"
40
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])41 void mexFunction
42 (
43 int nargout,
44 mxArray *pargout [ ],
45 int nargin,
46 const mxArray *pargin [ ]
47 )
48 {
49 double dummy = 0 ;
50 double *Lx, *px ;
51 Long *Parent, *Post, *ColCount, *First, *Level, *Rp, *Ri, *Lp, *Li, *W ;
52 cholmod_sparse *A, Amatrix, *F, *Aup, *Alo, *R, *A1, *A2, *L, *S ;
53 cholmod_common Common, *cm ;
54 Long n, i, coletree, j, lnz, p, k, height, c ;
55 char buf [LEN] ;
56
57 /* ---------------------------------------------------------------------- */
58 /* start CHOLMOD and set defaults */
59 /* ---------------------------------------------------------------------- */
60
61 cm = &Common ;
62 cholmod_l_start (cm) ;
63 sputil_config (SPUMONI, cm) ;
64
65 /* ---------------------------------------------------------------------- */
66 /* get inputs */
67 /* ---------------------------------------------------------------------- */
68
69 if (nargout > 5 || nargin < 1 || nargin > 3)
70 {
71 mexErrMsgTxt (
72 "Usage: [count h parent post R] = symbfact2 (A, mode, Lmode)") ;
73 }
74
75 /* ---------------------------------------------------------------------- */
76 /* get input matrix A */
77 /* ---------------------------------------------------------------------- */
78
79 A = sputil_get_sparse_pattern (pargin [0], &Amatrix, &dummy, cm) ;
80 S = (A == &Amatrix) ? NULL : A ;
81
82 /* ---------------------------------------------------------------------- */
83 /* get A->stype, default is to use triu(A) */
84 /* ---------------------------------------------------------------------- */
85
86 A->stype = 1 ;
87 n = A->nrow ;
88 coletree = FALSE ;
89 if (nargin > 1)
90 {
91 buf [0] = '\0' ;
92 if (mxIsChar (pargin [1]))
93 {
94 mxGetString (pargin [1], buf, LEN) ;
95 }
96 c = buf [0] ;
97 if (tolower (c) == 'r')
98 {
99 /* unsymmetric case (A*A') if string starts with 'r' */
100 A->stype = 0 ;
101 }
102 else if (tolower (c) == 'c')
103 {
104 /* unsymmetric case (A'*A) if string starts with 'c' */
105 n = A->ncol ;
106 coletree = TRUE ;
107 A->stype = 0 ;
108 }
109 else if (tolower (c) == 's')
110 {
111 /* symmetric upper case (A) if string starts with 's' */
112 A->stype = 1 ;
113 }
114 else if (tolower (c) == 'l')
115 {
116 /* symmetric lower case (A) if string starts with 'l' */
117 A->stype = -1 ;
118 }
119 else
120 {
121 mexErrMsgTxt ("symbfact2: unrecognized mode") ;
122 }
123 }
124
125 if (A->stype && A->nrow != A->ncol)
126 {
127 mexErrMsgTxt ("symbfact2: A must be square") ;
128 }
129
130 /* ---------------------------------------------------------------------- */
131 /* compute the etree, its postorder, and the row/column counts */
132 /* ---------------------------------------------------------------------- */
133
134 Parent = cholmod_l_malloc (n, sizeof (Long), cm) ;
135 Post = cholmod_l_malloc (n, sizeof (Long), cm) ;
136 ColCount = cholmod_l_malloc (n, sizeof (Long), cm) ;
137 First = cholmod_l_malloc (n, sizeof (Long), cm) ;
138 Level = cholmod_l_malloc (n, sizeof (Long), cm) ;
139
140 /* F = A' */
141 F = cholmod_l_transpose (A, 0, cm) ;
142
143 if (A->stype == 1 || coletree)
144 {
145 /* symmetric upper case: find etree of A, using triu(A) */
146 /* column case: find column etree of A, which is etree of A'*A */
147 Aup = A ;
148 Alo = F ;
149 }
150 else
151 {
152 /* symmetric lower case: find etree of A, using tril(A) */
153 /* row case: find row etree of A, which is etree of A*A' */
154 Aup = F ;
155 Alo = A ;
156 }
157
158 cholmod_l_etree (Aup, Parent, cm) ;
159
160 if (cm->status < CHOLMOD_OK)
161 {
162 /* out of memory or matrix invalid */
163 mexErrMsgTxt ("symbfact2 failed: matrix corrupted!") ;
164 }
165
166 if (cholmod_l_postorder (Parent, n, NULL, Post, cm) != n)
167 {
168 /* out of memory or Parent invalid */
169 mexErrMsgTxt ("symbfact2 postorder failed!") ;
170 }
171
172 /* symmetric upper case: analyze tril(F), which is triu(A) */
173 /* column case: analyze F*F', which is A'*A */
174 /* symmetric lower case: analyze tril(A) */
175 /* row case: analyze A*A' */
176 cholmod_l_rowcolcounts (Alo, NULL, 0, Parent, Post, NULL, ColCount,
177 First, Level, cm) ;
178
179 if (cm->status < CHOLMOD_OK)
180 {
181 /* out of memory or matrix invalid */
182 mexErrMsgTxt ("symbfact2 failed: matrix corrupted!") ;
183 }
184
185 /* ---------------------------------------------------------------------- */
186 /* return results to MATLAB: count, h, parent, and post */
187 /* ---------------------------------------------------------------------- */
188
189 pargout [0] = sputil_put_int (ColCount, n, 0) ;
190 if (nargout > 1)
191 {
192 /* compute the elimination tree height */
193 height = 0 ;
194 for (i = 0 ; i < n ; i++)
195 {
196 height = MAX (height, Level [i]) ;
197 }
198 height++ ;
199 pargout [1] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
200 px = mxGetPr (pargout [1]) ;
201 px [0] = height ;
202 }
203 if (nargout > 2)
204 {
205 pargout [2] = sputil_put_int (Parent, n, 1) ;
206 }
207 if (nargout > 3)
208 {
209 pargout [3] = sputil_put_int (Post, n, 1) ;
210 }
211
212 /* ---------------------------------------------------------------------- */
213 /* construct L, if requested */
214 /* ---------------------------------------------------------------------- */
215
216 if (nargout > 4)
217 {
218
219 if (A->stype == 1)
220 {
221 /* symmetric upper case: use triu(A) only, A2 not needed */
222 A1 = A ;
223 A2 = NULL ;
224 }
225 else if (A->stype == -1)
226 {
227 /* symmetric lower case: use tril(A) only, A2 not needed */
228 A1 = F ;
229 A2 = NULL ;
230 }
231 else if (coletree)
232 {
233 /* column case: analyze F*F' */
234 A1 = F ;
235 A2 = A ;
236 }
237 else
238 {
239 /* row case: analyze A*A' */
240 A1 = A ;
241 A2 = F ;
242 }
243
244 /* count the total number of entries in L */
245 lnz = 0 ;
246 for (j = 0 ; j < n ; j++)
247 {
248 lnz += ColCount [j] ;
249 }
250
251 /* allocate the output matrix L (pattern-only) */
252 L = cholmod_l_allocate_sparse (n, n, lnz, TRUE, TRUE, 0,
253 CHOLMOD_PATTERN, cm) ;
254 Lp = L->p ;
255 Li = L->i ;
256
257 /* initialize column pointers */
258 lnz = 0 ;
259 for (j = 0 ; j < n ; j++)
260 {
261 Lp [j] = lnz ;
262 lnz += ColCount [j] ;
263 }
264 Lp [j] = lnz ;
265
266 /* create a copy of the column pointers */
267 W = First ;
268 for (j = 0 ; j < n ; j++)
269 {
270 W [j] = Lp [j] ;
271 }
272
273 /* get workspace for computing one row of L */
274 R = cholmod_l_allocate_sparse (n, 1, n, FALSE, TRUE, 0, CHOLMOD_PATTERN,
275 cm) ;
276 Rp = R->p ;
277 Ri = R->i ;
278
279 /* compute L one row at a time */
280 for (k = 0 ; k < n ; k++)
281 {
282 /* get the kth row of L and store in the columns of L */
283 cholmod_l_row_subtree (A1, A2, k, Parent, R, cm) ;
284 for (p = 0 ; p < Rp [1] ; p++)
285 {
286 Li [W [Ri [p]]++] = k ;
287 }
288 /* add the diagonal entry */
289 Li [W [k]++] = k ;
290 }
291
292 /* free workspace */
293 cholmod_l_free_sparse (&R, cm) ;
294
295 /* transpose L to get R, or leave as is */
296 if (nargin < 3)
297 {
298 /* R = L' */
299 R = cholmod_l_transpose (L, 0, cm) ;
300 cholmod_l_free_sparse (&L, cm) ;
301 L = R ;
302 }
303
304 /* fill numerical values of L with one's (only MATLAB needs this...) */
305 L->x = cholmod_l_malloc (lnz, sizeof (double), cm) ;
306 Lx = L->x ;
307 for (p = 0 ; p < lnz ; p++)
308 {
309 Lx [p] = 1 ;
310 }
311 L->xtype = CHOLMOD_REAL ;
312
313 /* return L (or R) to MATLAB */
314 pargout [4] = sputil_put_sparse (&L, cm) ;
315 }
316
317 /* ---------------------------------------------------------------------- */
318 /* free workspace */
319 /* ---------------------------------------------------------------------- */
320
321 cholmod_l_free (n, sizeof (Long), Parent, cm) ;
322 cholmod_l_free (n, sizeof (Long), Post, cm) ;
323 cholmod_l_free (n, sizeof (Long), ColCount, cm) ;
324 cholmod_l_free (n, sizeof (Long), First, cm) ;
325 cholmod_l_free (n, sizeof (Long), Level, cm) ;
326 cholmod_l_free_sparse (&F, cm) ;
327 cholmod_l_free_sparse (&S, cm) ;
328 cholmod_l_finish (cm) ;
329 cholmod_l_print_common (" ", cm) ;
330 /*
331 if (cm->malloc_count != ((nargout == 5) ? 3:0)) mexErrMsgTxt ("!") ;
332 */
333 }
334