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