1 // =============================================================================
2 // === spqr mexFunction ========================================================
3 // =============================================================================
4
5 #include "spqr_mx.hpp"
6
7 // This function mimics the existing MATLAB QR syntax, with extensions.
8 // See spqr.m for details.
9
10 #define EMPTY (-1)
11 #define TRUE 1
12 #define FALSE 0
13
14 // =============================================================================
15 // === is_zero =================================================================
16 // =============================================================================
17
18 // Return TRUE if the MATLAB matrix s is a scalar equal to zero
19
is_zero(const mxArray * s)20 static int is_zero (const mxArray *s)
21 {
22 return (mxIsNumeric (s)
23 && (mxGetNumberOfElements (s) == 1)
24 && (mxGetScalar (s) == 0)) ;
25 }
26
27
28 // =============================================================================
29 // === spqr mexFunction ========================================================
30 // =============================================================================
31
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])32 void mexFunction
33 (
34 int nargout,
35 mxArray *pargout [ ],
36 int nargin,
37 const mxArray *pargin [ ]
38 )
39 {
40 Long *Ap, *Ai, *E, *Bp, *Bi, *HPinv ;
41 double *Ax, *Bx, dummy, tol ;
42 Long m, n, anz, bnz, is_complex, econ, A_complex, B_complex ;
43 spqr_mx_options opts ;
44 cholmod_sparse *A, Amatrix, *R, *Q, *Csparse, Bsmatrix, *Bsparse, *H ;
45 cholmod_dense *Cdense, Bdmatrix, *Bdense, *HTau ;
46 cholmod_common Common, *cc ;
47
48 double t0 = (nargout > 3) ? SuiteSparse_time ( ) : 0 ;
49
50 // -------------------------------------------------------------------------
51 // start CHOLMOD and set parameters
52 // -------------------------------------------------------------------------
53
54 cc = &Common ;
55 cholmod_l_start (cc) ;
56 spqr_mx_config (SPUMONI, cc) ;
57
58 // -------------------------------------------------------------------------
59 // check inputs
60 // -------------------------------------------------------------------------
61
62 // nargin can be 1, 2, or 3
63 // nargout can be 0, 1, 2, 3, or 4
64
65 if (nargout > 4)
66 {
67 mexErrMsgIdAndTxt ("MATLAB:maxlhs", "Too many output arguments") ;
68 }
69 if (nargin < 1)
70 {
71 mexErrMsgIdAndTxt ("MATLAB:minrhs", "Not enough input arguments") ;
72 }
73 if (nargin > 3)
74 {
75 mexErrMsgIdAndTxt ("MATLAB:maxrhs", "Too many input arguments") ;
76 }
77
78 // -------------------------------------------------------------------------
79 // get the input matrix A
80 // -------------------------------------------------------------------------
81
82 if (!mxIsSparse (pargin [0]))
83 {
84 mexErrMsgIdAndTxt ("QR:invalidInput", "A must be sparse") ;
85 }
86
87 A = spqr_mx_get_sparse (pargin [0], &Amatrix, &dummy) ;
88 Ap = (Long *) A->p ;
89 Ai = (Long *) A->i ;
90 m = A->nrow ;
91 n = A->ncol ;
92 A_complex = mxIsComplex (pargin [0]) ;
93
94 // -------------------------------------------------------------------------
95 // determine usage and parameters
96 // -------------------------------------------------------------------------
97
98 if (nargin == 1)
99 {
100
101 // ---------------------------------------------------------------------
102 // [ ] = qr (A)
103 // ---------------------------------------------------------------------
104
105 spqr_mx_get_options (NULL, &opts, m, nargout, cc) ;
106 // R = qr (A)
107 // [Q,R] = qr (A)
108 // [Q,R,E] = qr (A)
109
110 }
111 else if (nargin == 2)
112 {
113
114 // ---------------------------------------------------------------------
115 // [ ] = qr (A,0), [ ] = qr (A,opts), or [ ] = qr (A,B)
116 // ---------------------------------------------------------------------
117
118 if (is_zero (pargin [1]))
119 {
120
121 // -----------------------------------------------------------------
122 // [ ... ] = qr (A,0)
123 // -----------------------------------------------------------------
124
125 spqr_mx_get_options (NULL, &opts, m, nargout, cc) ;
126 opts.econ = n ;
127 opts.permvector = TRUE ;
128 // R = qr (A,0)
129 // [Q,R] = qr (A,0)
130 // [Q,R,E] = qr (A,0)
131
132 }
133 else if (mxIsEmpty (pargin [1]) || mxIsStruct (pargin [1]))
134 {
135
136 // -----------------------------------------------------------------
137 // [ ] = qr (A,opts)
138 // -----------------------------------------------------------------
139
140 spqr_mx_get_options (pargin [1], &opts, m, nargout, cc) ;
141 // R = qr (A,opts)
142 // [Q,R] = qr (A,opts)
143 // [Q,R,E] = qr (A,opts)
144
145 }
146 else
147 {
148
149 // -----------------------------------------------------------------
150 // [ ] = qr (A,B)
151 // -----------------------------------------------------------------
152
153 spqr_mx_get_options (NULL, &opts, m, nargout, cc) ;
154 opts.haveB = TRUE ;
155 opts.Qformat = SPQR_Q_DISCARD ;
156 if (nargout <= 1)
157 {
158 mexErrMsgIdAndTxt ("MATLAB:minlhs",
159 "Not enough output arguments") ;
160 }
161 // [C,R] = qr (A,B)
162 // [C,R,E] = qr (A,B)
163
164 }
165
166 }
167 else // if (nargin == 3)
168 {
169
170 // ---------------------------------------------------------------------
171 // [ ] = qr (A,B,opts)
172 // ---------------------------------------------------------------------
173
174 if (is_zero (pargin [2]))
175 {
176
177 // -----------------------------------------------------------------
178 // [ ] = qr (A,B,0)
179 // -----------------------------------------------------------------
180
181 spqr_mx_get_options (NULL, &opts, m, nargout, cc) ;
182 opts.econ = n ;
183 opts.permvector = TRUE ;
184 opts.haveB = TRUE ;
185 opts.Qformat = SPQR_Q_DISCARD ;
186 if (nargout <= 1)
187 {
188 mexErrMsgIdAndTxt ("MATLAB:minlhs",
189 "Not enough output arguments") ;
190 }
191 // [C,R] = qr (A,B,0)
192 // [C,R,E] = qr (A,B,0)
193
194 }
195 else if (mxIsEmpty (pargin [2]) || mxIsStruct (pargin [2]))
196 {
197
198 // -----------------------------------------------------------------
199 // [ ] = qr (A,B,opts)
200 // -----------------------------------------------------------------
201
202 spqr_mx_get_options (pargin [2], &opts, m, nargout, cc) ;
203 opts.haveB = TRUE ;
204 opts.Qformat = SPQR_Q_DISCARD ;
205 if (nargout <= 1)
206 {
207 mexErrMsgIdAndTxt ("MATLAB:minlhs",
208 "Not enough output arguments") ;
209 }
210 // [C,R] = qr (A,B,opts)
211 // [C,R,E] = qr (A,B,opts)
212
213 }
214 else
215 {
216 mexErrMsgIdAndTxt ("QR:invalidInput", "Invalid opts argument") ;
217 }
218 }
219
220 int order = opts.ordering ;
221 tol = opts.tol ;
222 econ = opts.econ ;
223
224 // -------------------------------------------------------------------------
225 // get A and convert to merged-complex if needed
226 // -------------------------------------------------------------------------
227
228 if (opts.haveB)
229 {
230 B_complex = mxIsComplex (pargin [1]) ;
231 }
232 else
233 {
234 B_complex = FALSE ;
235 }
236
237 is_complex = (A_complex || B_complex) ;
238 Ax = spqr_mx_merge_if_complex (pargin [0], is_complex, &anz, cc) ;
239 if (is_complex)
240 {
241 // A has been converted from real or zomplex to complex
242 A->x = Ax ;
243 A->z = NULL ;
244 A->xtype = CHOLMOD_COMPLEX ;
245 }
246
247 // -------------------------------------------------------------------------
248 // analyze, factorize, and get the results
249 // -------------------------------------------------------------------------
250
251 if (opts.haveB)
252 {
253
254 // ---------------------------------------------------------------------
255 // get B, and convert to complex if necessary
256 // ---------------------------------------------------------------------
257
258 if (!mxIsNumeric (pargin [1]))
259 {
260 mexErrMsgIdAndTxt ("QR:invalidInput", "invalid non-numeric B") ;
261 }
262 if (mxGetM (pargin [1]) != m)
263 {
264 mexErrMsgIdAndTxt ("QR:invalidInput",
265 "A and B must have the same number of rows") ;
266 }
267
268 // convert from real or zomplex to complex
269 Bx = spqr_mx_merge_if_complex (pargin [1], is_complex, &bnz, cc) ;
270
271 int B_is_sparse = mxIsSparse (pargin [1]) ;
272 if (B_is_sparse)
273 {
274 Bsparse = spqr_mx_get_sparse (pargin [1], &Bsmatrix, &dummy) ;
275 Bdense = NULL ;
276 if (is_complex)
277 {
278 // Bsparse has been converted from real or zomplex to complex
279 Bsparse->x = Bx ;
280 Bsparse->z = NULL ;
281 Bsparse->xtype = CHOLMOD_COMPLEX ;
282 }
283 }
284 else
285 {
286 Bsparse = NULL ;
287 Bdense = spqr_mx_get_dense (pargin [1], &Bdmatrix, &dummy) ;
288 if (is_complex)
289 {
290 // Bdense has been converted from real or zomplex to complex
291 Bdense->x = Bx ;
292 Bdense->z = NULL ;
293 Bdense->xtype = CHOLMOD_COMPLEX ;
294 }
295 }
296
297 // ---------------------------------------------------------------------
298 // [C,R,E] = qr (A,B,...) or [C,R] = qr (A,B,...)
299 // ---------------------------------------------------------------------
300
301 if (is_complex)
302 {
303
304 // -----------------------------------------------------------------
305 // [C,R,E] = qr (A,B): complex case
306 // -----------------------------------------------------------------
307
308 if (B_is_sparse)
309 {
310 // B and C are both sparse and complex
311 SuiteSparseQR <Complex> (order, tol, econ, A, Bsparse,
312 &Csparse, &R, &E, cc) ;
313 pargout [0] = spqr_mx_put_sparse (&Csparse, cc) ;
314 }
315 else
316 {
317 // B and C are both dense and complex
318 SuiteSparseQR <Complex> (order, tol, econ, A, Bdense,
319 &Cdense, &R, &E, cc) ;
320 pargout [0] = spqr_mx_put_dense (&Cdense, cc) ;
321 }
322
323 }
324 else
325 {
326
327 // -----------------------------------------------------------------
328 // [C,R,E] = qr (A,B): real case
329 // -----------------------------------------------------------------
330
331 if (B_is_sparse)
332 {
333 // B and C are both sparse and real
334 SuiteSparseQR <double> (order, tol, econ, A, Bsparse,
335 &Csparse, &R, &E, cc) ;
336 pargout [0] = spqr_mx_put_sparse (&Csparse, cc) ;
337 }
338 else
339 {
340 // B and C are both dense and real
341 SuiteSparseQR <double> (order, tol, econ, A, Bdense,
342 &Cdense, &R, &E, cc) ;
343 pargout [0] = spqr_mx_put_dense (&Cdense, cc) ;
344 }
345 }
346
347 pargout [1] = spqr_mx_put_sparse (&R, cc) ;
348
349 }
350 else if (nargout <= 1)
351 {
352
353 // ---------------------------------------------------------------------
354 // R = qr (A) or R = qr (A,opts)
355 // ---------------------------------------------------------------------
356
357 if (is_complex)
358 {
359 SuiteSparseQR <Complex> (0, tol, econ, A, &R, NULL, cc) ;
360 }
361 else
362 {
363 SuiteSparseQR <double> (0, tol, econ, A, &R, NULL, cc) ;
364 }
365 pargout [0] = spqr_mx_put_sparse (&R, cc) ;
366
367 }
368 else
369 {
370
371 // ---------------------------------------------------------------------
372 // [Q,R,E] = qr (A,...) or [Q,R] = qr (A,...)
373 // ---------------------------------------------------------------------
374
375 if (opts.Qformat == SPQR_Q_DISCARD)
376 {
377
378 // -----------------------------------------------------------------
379 // Q is discarded, and Q = [ ] is returned as a placeholder
380 // -----------------------------------------------------------------
381
382 if (is_complex)
383 {
384 SuiteSparseQR <Complex> (order, tol, econ, A, &R, &E, cc);
385 }
386 else
387 {
388 SuiteSparseQR <double> (order, tol, econ, A, &R, &E, cc) ;
389 }
390 pargout [0] = mxCreateDoubleMatrix (0, 0, mxREAL) ;
391
392 }
393 else if (opts.Qformat == SPQR_Q_MATRIX)
394 {
395
396 // -----------------------------------------------------------------
397 // Q is a sparse matrix
398 // -----------------------------------------------------------------
399
400 if (is_complex)
401 {
402 SuiteSparseQR <Complex> (order, tol, econ, A, &Q, &R, &E, cc) ;
403 }
404 else
405 {
406 SuiteSparseQR <double> (order, tol, econ, A, &Q, &R, &E, cc) ;
407 }
408 pargout [0] = spqr_mx_put_sparse (&Q, cc) ;
409
410 }
411 else
412 {
413
414 // -----------------------------------------------------------------
415 // H is kept, and Q is a struct containing H, Tau, and P
416 // -----------------------------------------------------------------
417
418 mxArray *Tau, *P, *Hmatlab ;
419 if (is_complex)
420 {
421 SuiteSparseQR <Complex> (order, tol, econ, A,
422 &R, &E, &H, &HPinv, &HTau, cc) ;
423 }
424 else
425 {
426 SuiteSparseQR <double> (order, tol, econ, A,
427 &R, &E, &H, &HPinv, &HTau, cc) ;
428 }
429
430 Tau = spqr_mx_put_dense (&HTau, cc) ;
431 Hmatlab = spqr_mx_put_sparse (&H, cc) ;
432
433 // Q.P contains the inverse row permutation
434 P = mxCreateDoubleMatrix (1, m, mxREAL) ;
435 double *Tx = mxGetPr (P) ;
436 for (Long i = 0 ; i < m ; i++)
437 {
438 Tx [i] = HPinv [i] + 1 ;
439 }
440
441 // return Q
442 const char *Qstruct [ ] = { "H", "Tau", "P" } ;
443 pargout [0] = mxCreateStructMatrix (1, 1, 3, Qstruct) ;
444 mxSetFieldByNumber (pargout [0], 0, 0, Hmatlab) ;
445 mxSetFieldByNumber (pargout [0], 0, 1, Tau) ;
446 mxSetFieldByNumber (pargout [0], 0, 2, P) ;
447
448 }
449 pargout [1] = spqr_mx_put_sparse (&R, cc) ;
450 }
451
452 // -------------------------------------------------------------------------
453 // return E
454 // -------------------------------------------------------------------------
455
456 if (nargout > 2)
457 {
458 pargout [2] = spqr_mx_put_permutation (E, n, opts.permvector, cc) ;
459 }
460
461 // -------------------------------------------------------------------------
462 // free copy of merged-complex, if needed
463 // -------------------------------------------------------------------------
464
465 if (is_complex)
466 {
467 // this was allocated by merge_if_complex
468 cholmod_l_free (anz, sizeof (Complex), Ax, cc) ;
469 if (opts.haveB)
470 {
471 cholmod_l_free (bnz, sizeof (Complex), Bx, cc) ;
472 }
473 }
474
475 // -------------------------------------------------------------------------
476 // info output
477 // -------------------------------------------------------------------------
478
479 if (nargout > 3)
480 {
481 double flops = cc->SPQR_flopcount ;
482 double t = SuiteSparse_time ( ) - t0 ;
483 pargout [3] = spqr_mx_info (cc, t, flops) ;
484 }
485
486 cholmod_l_finish (cc) ;
487 if (opts.spumoni > 0) spqr_mx_spumoni (&opts, is_complex, cc) ;
488 }
489