1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/chol2 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 R'R factorization. Note that LL' and LDL' are faster than R'R
12 * and use less memory. The R'R factorization methods use triu(A), just like
13 * MATLAB's built-in chol.
14 *
15 * R = chol2 (A) same as R = chol (A), just faster
16 * [R,p] = chol2 (A) save as [R,p] = chol(A), just faster
17 * [R,p,q] = chol2 (A) factorizes A(q,q) into R'*R
18 *
19 * A must be sparse. It can be complex or real.
20 *
21 * R is returned with no explicit zero entries. This means it might not be
22 * chordal, and R cannot be passed back to CHOLMOD for an update/downdate or
23 * for a fast simplicial solve. spones (R) will be equal to the R returned
24 * by symbfact2 if no numerically zero entries are dropped, or a subset
25 * otherwise.
26 */
27
28 #include "cholmod_matlab.h"
29
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])30 void mexFunction
31 (
32 int nargout,
33 mxArray *pargout [ ],
34 int nargin,
35 const mxArray *pargin [ ]
36 )
37 {
38 double dummy = 0, *px ;
39 cholmod_sparse Amatrix, *A, *Lsparse, *R ;
40 cholmod_factor *L ;
41 cholmod_common Common, *cm ;
42 Long n, minor ;
43
44 /* ---------------------------------------------------------------------- */
45 /* start CHOLMOD and set parameters */
46 /* ---------------------------------------------------------------------- */
47
48 cm = &Common ;
49 cholmod_l_start (cm) ;
50 sputil_config (SPUMONI, cm) ;
51
52 /* convert to packed LL' when done */
53 cm->final_asis = FALSE ;
54 cm->final_super = FALSE ;
55 cm->final_ll = TRUE ;
56 cm->final_pack = TRUE ;
57 cm->final_monotonic = TRUE ;
58
59 /* no need to prune entries due to relaxed supernodal amalgamation, since
60 * zeros are dropped with sputil_drop_zeros instead */
61 cm->final_resymbol = FALSE ;
62
63 cm->quick_return_if_not_posdef = (nargout < 2) ;
64
65 /* ---------------------------------------------------------------------- */
66 /* get inputs */
67 /* ---------------------------------------------------------------------- */
68
69 if (nargin != 1 || nargout > 3)
70 {
71 mexErrMsgTxt ("usage: [R,p,q] = chol2 (A)") ;
72 }
73
74 n = mxGetN (pargin [0]) ;
75
76 if (!mxIsSparse (pargin [0]) || n != mxGetM (pargin [0]))
77 {
78 mexErrMsgTxt ("A must be square and sparse") ;
79 }
80
81 /* get input sparse matrix A. Use triu(A) only */
82 A = sputil_get_sparse (pargin [0], &Amatrix, &dummy, 1) ;
83
84 /* use natural ordering if no q output parameter */
85 if (nargout < 3)
86 {
87 cm->nmethods = 1 ;
88 cm->method [0].ordering = CHOLMOD_NATURAL ;
89 cm->postorder = FALSE ;
90 }
91
92 /* ---------------------------------------------------------------------- */
93 /* analyze and factorize */
94 /* ---------------------------------------------------------------------- */
95
96 L = cholmod_l_analyze (A, cm) ;
97 cholmod_l_factorize (A, L, cm) ;
98
99 if (nargout < 2 && cm->status != CHOLMOD_OK)
100 {
101 mexErrMsgTxt ("matrix is not positive definite") ;
102 }
103
104 /* ---------------------------------------------------------------------- */
105 /* convert L to a sparse matrix */
106 /* ---------------------------------------------------------------------- */
107
108 /* the conversion sets L->minor back to n, so get a copy of it first */
109 minor = L->minor ;
110 Lsparse = cholmod_l_factor_to_sparse (L, cm) ;
111 if (Lsparse->xtype == CHOLMOD_COMPLEX)
112 {
113 /* convert Lsparse from complex to zomplex */
114 cholmod_l_sparse_xtype (CHOLMOD_ZOMPLEX, Lsparse, cm) ;
115 }
116
117 if (minor < n)
118 {
119 /* remove columns minor to n-1 from Lsparse */
120 sputil_trim (Lsparse, minor, cm) ;
121 }
122
123 /* drop zeros from Lsparse */
124 sputil_drop_zeros (Lsparse) ;
125
126 /* Lsparse is lower triangular; conjugate transpose to get R */
127 R = cholmod_l_transpose (Lsparse, 2, cm) ;
128 cholmod_l_free_sparse (&Lsparse, cm) ;
129
130 /* ---------------------------------------------------------------------- */
131 /* return results to MATLAB */
132 /* ---------------------------------------------------------------------- */
133
134 /* return R */
135 pargout [0] = sputil_put_sparse (&R, cm) ;
136
137 /* return minor (translate to MATLAB convention) */
138 if (nargout > 1)
139 {
140 pargout [1] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
141 px = mxGetPr (pargout [1]) ;
142 px [0] = ((minor == n) ? 0 : (minor+1)) ;
143 }
144
145 /* return permutation */
146 if (nargout > 2)
147 {
148 pargout [2] = sputil_put_int (L->Perm, n, 1) ;
149 }
150
151 /* ---------------------------------------------------------------------- */
152 /* free workspace and the CHOLMOD L, except for what is copied to MATLAB */
153 /* ---------------------------------------------------------------------- */
154
155 cholmod_l_free_factor (&L, cm) ;
156 cholmod_l_finish (cm) ;
157 cholmod_l_print_common (" ", cm) ;
158 /*
159 if (cm->malloc_count != (3 + mxIsComplex (pargout[0]))) mexErrMsgTxt ("!") ;
160 */
161 }
162