1 // =============================================================================
2 // === spqr_solve mexFunction ==================================================
3 // =============================================================================
4
5 #include "spqr_mx.hpp"
6
7 /* x = A\b using a sparse QR factorization.
8
9 x = spqr_solve (A,b,opts) ;
10
11 A 2nd optional output gives statistics, in a self-explainable struct.
12
13 Let [m n] = size (A).
14
15 opts is an optional struct with the following fields. Fields not
16 present are set to their defaults. The defaults are:
17
18 opts.tol = 20 * (m+n) * eps * sqrt(max(diag(A'*A)))
19 opts.ordering = "colamd"
20
21 opts.tol: Default 20 * (m+n) * eps * sqrt(max(diag(A'*A))),
22 that is, 20 * (m+n) * eps * (max 2-norm of the columns of A).
23
24 If tol >= 0, then an approximate rank-revealing QR factorization
25 is used. A diagonal entry of R with magnitude < tol is treated
26 as zero; the matrix A is considered rank deficient.
27
28 If tol = -1, then no rank detection is attempted. This allows
29 the sparse QR to use less memory than when attempting to detect
30 rank, even for full-rank matrices. Thus, if you know your
31 matrix has full rank, tol = -1 will reduce memory usage. Zeros
32 may appear on the diagonal of R if A is rank deficient.
33
34 If tol <= -2, then the default tolerance is used.
35
36 opts.ordering: Default is "default"
37 See spqr.m for info.
38
39 opts.solution: default "basic"
40 Defines how an underdetermined system should be solved (m < n).
41 "basic": compute a basic solution, using
42 [C,R,E]=spqr(A,B) ; X=E*(R\C) ; where C=Q'*B
43 "min2norm": compute a minimum 2-norm solution, using
44 [C,R,E]=spqr(A') ; X = Q*(R'\(E'*B))
45 A min2norm solution is more costly to compute, in time
46 and memory, than a basic solution.
47 This option is ignored if m >= n.
48 */
49
50 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
51
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])52 void mexFunction
53 (
54 int nargout,
55 mxArray *pargout [ ],
56 int nargin,
57 const mxArray *pargin [ ]
58 )
59 {
60 Long *Bp, *Bi ;
61 double *Ax, *Bx, dummy ;
62 Long m, n, k, bncols, p, i, rank, A_complex, B_complex, is_complex,
63 anz, bnz ;
64 spqr_mx_options opts ;
65 cholmod_sparse *A, Amatrix, *Xsparse ;
66 cholmod_dense *Xdense ;
67 cholmod_common Common, *cc ;
68 char msg [LEN+1] ;
69
70 double t0 = (nargout > 1) ? SuiteSparse_time ( ) : 0 ;
71
72 // -------------------------------------------------------------------------
73 // start CHOLMOD and set parameters
74 // -------------------------------------------------------------------------
75
76 cc = &Common ;
77 cholmod_l_start (cc) ;
78 spqr_mx_config (SPUMONI, cc) ;
79
80 // -------------------------------------------------------------------------
81 // check inputs
82 // -------------------------------------------------------------------------
83
84 if (nargout > 2)
85 {
86 mexErrMsgIdAndTxt ("MATLAB:maxlhs", "Too many output arguments") ;
87 }
88 if (nargin < 2)
89 {
90 mexErrMsgIdAndTxt ("MATLAB:minrhs", "Not enough input arguments") ;
91 }
92 if (nargin > 3)
93 {
94 mexErrMsgIdAndTxt ("MATLAB:maxrhs", "Too many input arguments") ;
95 }
96
97 // -------------------------------------------------------------------------
98 // get the input matrix A (must be sparse)
99 // -------------------------------------------------------------------------
100
101 if (!mxIsSparse (pargin [0]))
102 {
103 mexErrMsgIdAndTxt ("QR:invalidInput", "A must be sparse") ;
104 }
105
106 A = spqr_mx_get_sparse (pargin [0], &Amatrix, &dummy) ;
107 m = A->nrow ;
108 n = A->ncol ;
109 A_complex = mxIsComplex (pargin [0]) ;
110
111 B_complex = mxIsComplex (pargin [1]) ;
112 is_complex = (A_complex || B_complex) ;
113 Ax = spqr_mx_merge_if_complex (pargin [0], is_complex, &anz, cc) ;
114 if (is_complex)
115 {
116 // A has been converted from real or zomplex to complex
117 A->x = Ax ;
118 A->z = NULL ;
119 A->xtype = CHOLMOD_COMPLEX ;
120 }
121
122 // -------------------------------------------------------------------------
123 // determine usage and parameters
124 // -------------------------------------------------------------------------
125
126 spqr_mx_get_options ((nargin < 3) ? NULL : pargin [2], &opts, m, 3, cc) ;
127
128 opts.Qformat = SPQR_Q_DISCARD ;
129 opts.econ = 0 ;
130 opts.permvector = TRUE ;
131 opts.haveB = TRUE ;
132
133 // -------------------------------------------------------------------------
134 // get the input matrix B (sparse or dense)
135 // -------------------------------------------------------------------------
136
137 if (!mxIsNumeric (pargin [1]))
138 {
139 mexErrMsgIdAndTxt ("QR:invalidInput", "invalid non-numeric B") ;
140 }
141 if (mxGetM (pargin [1]) != m)
142 {
143 mexErrMsgIdAndTxt ("QR:invalidInput",
144 "A and B must have the same number of rows") ;
145 }
146
147 cholmod_sparse Bsmatrix, *Bsparse ;
148 cholmod_dense Bdmatrix, *Bdense ;
149
150 // convert from real or zomplex to complex
151 Bx = spqr_mx_merge_if_complex (pargin [1], is_complex, &bnz, cc) ;
152
153 int B_is_sparse = mxIsSparse (pargin [1]) ;
154 if (B_is_sparse)
155 {
156 Bsparse = spqr_mx_get_sparse (pargin [1], &Bsmatrix, &dummy) ;
157 Bdense = NULL ;
158 if (is_complex)
159 {
160 // Bsparse has been converted from real or zomplex to complex
161 Bsparse->x = Bx ;
162 Bsparse->z = NULL ;
163 Bsparse->xtype = CHOLMOD_COMPLEX ;
164 }
165 }
166 else
167 {
168 Bsparse = NULL ;
169 Bdense = spqr_mx_get_dense (pargin [1], &Bdmatrix, &dummy) ;
170 if (is_complex)
171 {
172 // Bdense has been converted from real or zomplex to complex
173 Bdense->x = Bx ;
174 Bdense->z = NULL ;
175 Bdense->xtype = CHOLMOD_COMPLEX ;
176 }
177 }
178
179 // -------------------------------------------------------------------------
180 // X = A\B
181 // -------------------------------------------------------------------------
182
183 if (opts.min2norm && m < n)
184 {
185 #ifndef NEXPERT
186 // This requires SuiteSparseQR_expert.cpp
187 if (is_complex)
188 {
189 if (B_is_sparse)
190 {
191 // X and B are both sparse and complex
192 Xsparse = SuiteSparseQR_min2norm <Complex> (opts.ordering,
193 opts.tol, A, Bsparse, cc) ;
194 pargout [0] = spqr_mx_put_sparse (&Xsparse, cc) ;
195 }
196 else
197 {
198 // X and B are both dense and complex
199 Xdense = SuiteSparseQR_min2norm <Complex> (opts.ordering,
200 opts.tol, A, Bdense, cc) ;
201 pargout [0] = spqr_mx_put_dense (&Xdense, cc) ;
202 }
203 }
204 else
205 {
206 if (B_is_sparse)
207 {
208 // X and B are both sparse and real
209 Xsparse = SuiteSparseQR_min2norm <double> (opts.ordering,
210 opts.tol, A, Bsparse, cc) ;
211 pargout [0] = spqr_mx_put_sparse (&Xsparse, cc) ;
212 }
213 else
214 {
215 // X and B are both dense and real
216 Xdense = SuiteSparseQR_min2norm <double> (opts.ordering,
217 opts.tol, A, Bdense, cc) ;
218 pargout [0] = spqr_mx_put_dense (&Xdense, cc) ;
219 }
220 }
221 #else
222 mexErrMsgIdAndTxt ("QR:notInstalled", "min2norm method not installed") ;
223 #endif
224
225 }
226 else
227 {
228
229 if (is_complex)
230 {
231 if (B_is_sparse)
232 {
233 // X and B are both sparse and complex
234 Xsparse = SuiteSparseQR <Complex> (opts.ordering, opts.tol,
235 A, Bsparse, cc) ;
236 pargout [0] = spqr_mx_put_sparse (&Xsparse, cc) ;
237 }
238 else
239 {
240 // X and B are both dense and complex
241 Xdense = SuiteSparseQR <Complex> (opts.ordering, opts.tol,
242 A, Bdense, cc) ;
243 pargout [0] = spqr_mx_put_dense (&Xdense, cc) ;
244 }
245 }
246 else
247 {
248 if (B_is_sparse)
249 {
250 // X and B are both sparse and real
251 Xsparse = SuiteSparseQR <double> (opts.ordering, opts.tol,
252 A, Bsparse, cc) ;
253 pargout [0] = spqr_mx_put_sparse (&Xsparse, cc) ;
254 }
255 else
256 {
257 // X and B are both dense and real
258 Xdense = SuiteSparseQR <double> (opts.ordering, opts.tol,
259 A, Bdense, cc) ;
260 pargout [0] = spqr_mx_put_dense (&Xdense, cc) ;
261 }
262 }
263 }
264
265 // -------------------------------------------------------------------------
266 // info output
267 // -------------------------------------------------------------------------
268
269 if (nargout > 1)
270 {
271 double flops = cc->SPQR_flopcount ;
272 double t = SuiteSparse_time ( ) - t0 ;
273 pargout [1] = spqr_mx_info (cc, t, flops) ;
274 }
275
276 // -------------------------------------------------------------------------
277 // warn if rank deficient
278 // -------------------------------------------------------------------------
279
280 rank = cc->SPQR_istat [4] ;
281 if (rank < MIN (m,n))
282 {
283 // snprintf would be safer, but Windows is oblivious to safety ...
284 // (Visual Studio C++ 2008 does not recognize snprintf!)
285 sprintf (msg, "rank deficient. rank = %ld tol = %g\n", rank,
286 cc->SPQR_tol_used) ;
287 mexWarnMsgIdAndTxt ("MATLAB:rankDeficientMatrix", msg) ;
288 }
289
290 if (is_complex)
291 {
292 // free the merged complex copies of A and B
293 cholmod_l_free (anz, sizeof (Complex), Ax, cc) ;
294 cholmod_l_free (bnz, sizeof (Complex), Bx, cc) ;
295 }
296
297 cholmod_l_finish (cc) ;
298 if (opts.spumoni > 0) spqr_mx_spumoni (&opts, is_complex, cc) ;
299 }
300