1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/lsubsolve mexFunction ================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MATLAB Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * http://www.suitesparse.com
8 * MATLAB(tm) is a Trademark of The MathWorks, Inc.
9 * -------------------------------------------------------------------------- */
10
11 /* [x xset] = lsubsolve (L,kind,P,b,system)
12 *
13 * L is a sparse lower triangular matrix, of size n-by-n. kind = 0 if
14 * L is from an LL' factorization, 1 if from LDL' (in which case L contains
15 * L in the lower part and D on the diagonal). P is a permutation vector, or
16 * empty (which means P is identity). b is a sparse column vector.
17 * system is a number between 0 and 8:
18 *
19 * Given L or LD, a permutation P, and a sparse right-hand size b,
20 * solve one of the following systems:
21 *
22 * Ax=b 0: CHOLMOD_A also applies the permutation L->Perm
23 * LDL'x=b 1: CHOLMOD_LDLt does not apply L->Perm
24 * LDx=b 2: CHOLMOD_LD
25 * DL'x=b 3: CHOLMOD_DLt
26 * Lx=b 4: CHOLMOD_L
27 * L'x=b 5: CHOLMOD_Lt
28 * Dx=b 6: CHOLMOD_D
29 * x=Pb 7: CHOLMOD_P apply a permutation (P is L->Perm)
30 * x=P'b 8: CHOLMOD_Pt apply an inverse permutation
31 *
32 * The solution x is a dense vector, but it is a subset of the entire solution,
33 * x is zero except where xset is 1. xset is reach of b (or P*b) in the graph
34 * of L. If P is empty then it is treated as the identity permutation.
35 *
36 * No zeros can be dropped from the stucture of the Cholesky factorization
37 * because this function gets its elimination tree from L itself. See ldlchol
38 * for a method of constructing a sparse L with explicit zeros that are
39 * normally dropped by MATLAB.
40 *
41 * This function is only meant for testing the cholmod_solve2 function. The
42 * cholmod_solve2 function takes O(flops) time, but the setup time in this
43 * mexFunction wrapper can dominate that time.
44 */
45
46 #include "cholmod_matlab.h"
47
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])48 void mexFunction
49 (
50 int nargout,
51 mxArray *pargout [ ],
52 int nargin,
53 const mxArray *pargin [ ]
54 )
55 {
56 double dummy = 0, *Px, *Xsetx ;
57 Long *Lp, *Lnz, *Xp, *Xi, xnz, *Perm, *Lprev, *Lnext, *Xsetp ;
58 cholmod_sparse *Bset, Bmatrix, *Xset ;
59 cholmod_dense *Bdense, *X, *Y, *E ;
60 cholmod_factor *L ;
61 cholmod_common Common, *cm ;
62 Long k, j, n, head, tail, xsetlen ;
63 int sys, kind ;
64
65 /* ---------------------------------------------------------------------- */
66 /* start CHOLMOD and set parameters */
67 /* ---------------------------------------------------------------------- */
68
69 cm = &Common ;
70 cholmod_l_start (cm) ;
71 sputil_config (SPUMONI, cm) ;
72
73 /* ---------------------------------------------------------------------- */
74 /* check inputs */
75 /* ---------------------------------------------------------------------- */
76
77 if (nargin != 5 || nargout > 2)
78 {
79 mexErrMsgTxt ("usage: [x xset] = lsubsolve (L,kind,P,b,system)") ;
80 }
81
82 n = mxGetN (pargin [0]) ;
83 if (!mxIsSparse (pargin [0]) || n != mxGetM (pargin [0]))
84 {
85 mexErrMsgTxt ("lsubsolve: L must be sparse and square") ;
86 }
87 if (mxGetNumberOfElements (pargin [1]) != 1)
88 {
89 mexErrMsgTxt ("lsubsolve: kind must be a scalar") ;
90 }
91
92 if (mxIsSparse (pargin [2]) ||
93 !(mxIsEmpty (pargin [2]) || mxGetNumberOfElements (pargin [2]) == n))
94 {
95 mexErrMsgTxt ("lsubsolve: P must be size n, or empty") ;
96 }
97
98 if (mxGetM (pargin [3]) != n || mxGetN (pargin [3]) != 1)
99 {
100 mexErrMsgTxt ("lsubsolve: b wrong dimension") ;
101 }
102 if (!mxIsSparse (pargin [3]))
103 {
104 mexErrMsgTxt ("lxbpattern: b must be sparse") ;
105 }
106 if (mxGetNumberOfElements (pargin [4]) != 1)
107 {
108 mexErrMsgTxt ("lsubsolve: system must be a scalar") ;
109 }
110
111 /* ---------------------------------------------------------------------- */
112 /* get the inputs */
113 /* ---------------------------------------------------------------------- */
114
115 kind = (int) sputil_get_integer (pargin [1], FALSE, 0) ;
116 sys = (int) sputil_get_integer (pargin [4], FALSE, 0) ;
117
118 /* ---------------------------------------------------------------------- */
119 /* get the sparse b */
120 /* ---------------------------------------------------------------------- */
121
122 /* get sparse matrix B (unsymmetric) */
123 Bset = sputil_get_sparse (pargin [3], &Bmatrix, &dummy, 0) ;
124 Bdense = cholmod_l_sparse_to_dense (Bset, cm) ;
125 Bset->x = NULL ;
126 Bset->z = NULL ;
127 Bset->xtype = CHOLMOD_PATTERN ;
128
129 /* ---------------------------------------------------------------------- */
130 /* construct a shallow copy of the input sparse matrix L */
131 /* ---------------------------------------------------------------------- */
132
133 /* the construction of the CHOLMOD takes O(n) time and memory */
134
135 /* allocate the CHOLMOD symbolic L */
136 L = cholmod_l_allocate_factor (n, cm) ;
137 L->ordering = CHOLMOD_NATURAL ;
138
139 /* get the MATLAB L */
140 L->p = mxGetJc (pargin [0]) ;
141 L->i = mxGetIr (pargin [0]) ;
142 L->x = mxGetPr (pargin [0]) ;
143 L->z = mxGetPi (pargin [0]) ;
144
145 /* allocate and initialize the rest of L */
146 L->nz = cholmod_l_malloc (n, sizeof (Long), cm) ;
147 Lp = L->p ;
148 Lnz = L->nz ;
149 for (j = 0 ; j < n ; j++)
150 {
151 Lnz [j] = Lp [j+1] - Lp [j] ;
152 }
153
154 /* these pointers are not accessed in cholmod_solve2 */
155 L->prev = cholmod_l_malloc (n+2, sizeof (Long), cm) ;
156 L->next = cholmod_l_malloc (n+2, sizeof (Long), cm) ;
157 Lprev = L->prev ;
158 Lnext = L->next ;
159
160 head = n+1 ;
161 tail = n ;
162 Lnext [head] = 0 ;
163 Lprev [head] = -1 ;
164 Lnext [tail] = -1 ;
165 Lprev [tail] = n-1 ;
166 for (j = 0 ; j < n ; j++)
167 {
168 Lnext [j] = j+1 ;
169 Lprev [j] = j-1 ;
170 }
171 Lprev [0] = head ;
172
173 L->xtype = (mxIsComplex (pargin [0])) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
174 L->nzmax = Lp [n] ;
175
176 /* get the permutation */
177 if (mxIsEmpty (pargin [2]))
178 {
179 L->Perm = NULL ;
180 Perm = NULL ;
181 }
182 else
183 {
184 L->ordering = CHOLMOD_GIVEN ;
185 L->Perm = cholmod_l_malloc (n, sizeof (Long), cm) ;
186 Perm = L->Perm ;
187 Px = mxGetPr (pargin [2]) ;
188 for (k = 0 ; k < n ; k++)
189 {
190 Perm [k] = ((Long) Px [k]) - 1 ;
191 }
192 }
193
194 /* set the kind, LL' or LDL' */
195 L->is_ll = (kind == 0) ;
196 /*
197 cholmod_l_print_factor (L, "L", cm) ;
198 */
199
200 /* ---------------------------------------------------------------------- */
201 /* solve the system */
202 /* ---------------------------------------------------------------------- */
203
204 X = cholmod_l_zeros (n, 1, L->xtype, cm) ;
205 Xset = NULL ;
206 Y = NULL ;
207 E = NULL ;
208
209 cholmod_l_solve2 (sys, L, Bdense, Bset, &X, &Xset, &Y, &E, cm) ;
210
211 cholmod_l_free_dense (&Y, cm) ;
212 cholmod_l_free_dense (&E, cm) ;
213
214 /* ---------------------------------------------------------------------- */
215 /* return result */
216 /* ---------------------------------------------------------------------- */
217
218 pargout [0] = sputil_put_dense (&X, cm) ;
219
220 /* fill numerical values of Xset with one's */
221 Xsetp = Xset->p ;
222 xsetlen = Xsetp [1] ;
223 Xset->x = cholmod_l_malloc (xsetlen, sizeof (double), cm) ;
224 Xsetx = Xset->x ;
225 for (k = 0 ; k < xsetlen ; k++)
226 {
227 Xsetx [k] = 1 ;
228 }
229 Xset->xtype = CHOLMOD_REAL ;
230
231 pargout [1] = sputil_put_sparse (&Xset, cm) ;
232
233 /* ---------------------------------------------------------------------- */
234 /* free workspace and the CHOLMOD L, except for what is copied to MATLAB */
235 /* ---------------------------------------------------------------------- */
236
237 L->p = NULL ;
238 L->i = NULL ;
239 L->x = NULL ;
240 L->z = NULL ;
241 cholmod_l_free_factor (&L, cm) ;
242 cholmod_l_finish (cm) ;
243 cholmod_l_print_common (" ", cm) ;
244 }
245