1 // =============================================================================
2 // === SuiteSparseQR_C =========================================================
3 // =============================================================================
4
5 // This C++ file provides a set of C-callable wrappers so that a C program can
6 // call SuiteSparseQR.
7
8 #include "spqr.hpp"
9 #include "SuiteSparseQR_C.h"
10
11 extern "C" {
12
13 // =============================================================================
14 // === SuiteSparseQR_C =========================================================
15 // =============================================================================
16
17 // Primary sparse QR function, with all inputs/outputs available. The primary
18 // uses of this function are to perform any one of the the MATLAB equivalent
19 // statements:
20 //
21 // X = A\B % where B is sparse or dense
22 // [C,R,E] = qr (A,B) % where Q*R=A*E and C=Q'*B
23 // [Q,R,E] = qr (A) % with Q in Householder form (H, HPinv, HTau)
24 // [Q,R,E] = qr (A) % where Q is discarded (the "Q-less" QR)
25 // R = qr (A) % as above, but with E=I
26 //
27 // To obtain the factor Q in sparse matrix form, use SuiteSparseQR_C_QR instead.
28
SuiteSparseQR_C(int ordering,double tol,Long econ,int getCTX,cholmod_sparse * A,cholmod_sparse * Bsparse,cholmod_dense * Bdense,cholmod_sparse ** Zsparse,cholmod_dense ** Zdense,cholmod_sparse ** R,Long ** E,cholmod_sparse ** H,Long ** HPinv,cholmod_dense ** HTau,cholmod_common * cc)29 Long SuiteSparseQR_C // returns rank(A) estimate, (-1) if failure
30 (
31 // inputs:
32 int ordering, // all, except 3:given treated as 0:fixed
33 double tol, // columns with 2-norm <= tol are treated as 0
34 Long econ, // e = max(min(m,econ),rank(A))
35 int getCTX, // 0: Z=C (e-by-k), 1: Z=C', 2: Z=X (e-by-k)
36 cholmod_sparse *A, // m-by-n sparse matrix to factorize
37 cholmod_sparse *Bsparse,// sparse m-by-k B
38 cholmod_dense *Bdense, // dense m-by-k B
39 // outputs:
40 cholmod_sparse **Zsparse, // sparse Z
41 cholmod_dense **Zdense, // dense Z
42 cholmod_sparse **R, // R factor, e-by-n
43 Long **E, // size n column permutation, NULL if identity
44 cholmod_sparse **H, // m-by-nh Householder vectors
45 Long **HPinv, // size m row permutation
46 cholmod_dense **HTau, // 1-by-nh Householder coefficients
47 cholmod_common *cc // workspace and parameters
48 )
49 {
50 RETURN_IF_NULL_COMMON (EMPTY) ;
51 RETURN_IF_NULL (A, EMPTY) ;
52 cc->status = CHOLMOD_OK ;
53
54 return ((A->xtype == CHOLMOD_REAL) ?
55 SuiteSparseQR <double> (ordering, tol, econ, getCTX, A, Bsparse,
56 Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc) :
57 SuiteSparseQR <Complex> (ordering, tol, econ, getCTX, A, Bsparse,
58 Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc)) ;
59 }
60
61 // =============================================================================
62 // === SuiteSparseQR_C_QR ======================================================
63 // =============================================================================
64
65 // [Q,R,E] = qr(A), returning Q as a sparse matrix
66
SuiteSparseQR_C_QR(int ordering,double tol,Long econ,cholmod_sparse * A,cholmod_sparse ** Q,cholmod_sparse ** R,Long ** E,cholmod_common * cc)67 Long SuiteSparseQR_C_QR // returns rank(A) estimate, (-1) if failure
68 (
69 // inputs:
70 int ordering, // all, except 3:given treated as 0:fixed
71 double tol, // columns with 2-norm <= tol are treated as 0
72 Long econ, // e = max(min(m,econ),rank(A))
73 cholmod_sparse *A, // m-by-n sparse matrix to factorize
74 // outputs:
75 cholmod_sparse **Q, // m-by-e sparse matrix
76 cholmod_sparse **R, // e-by-n sparse matrix
77 Long **E, // size n column permutation, NULL if identity
78 cholmod_common *cc // workspace and parameters
79 )
80 {
81 RETURN_IF_NULL_COMMON (EMPTY) ;
82 RETURN_IF_NULL (A, EMPTY) ;
83 cc->status = CHOLMOD_OK ;
84
85 return ((A->xtype == CHOLMOD_REAL) ?
86 SuiteSparseQR <double> (ordering, tol, econ, A, Q, R, E, cc) :
87 SuiteSparseQR <Complex> (ordering, tol, econ, A, Q, R, E, cc)) ;
88 }
89
90 // =============================================================================
91 // === SuiteSparseQR_C_backslash ===============================================
92 // =============================================================================
93
94 // X = A\B where B is dense
SuiteSparseQR_C_backslash(int ordering,double tol,cholmod_sparse * A,cholmod_dense * B,cholmod_common * cc)95 cholmod_dense *SuiteSparseQR_C_backslash // returns X, NULL if failure
96 (
97 // inputs:
98 int ordering, // all, except 3:given treated as 0:fixed
99 double tol, // columns with 2-norm <= tol are treated as 0
100 cholmod_sparse *A, // m-by-n sparse matrix
101 cholmod_dense *B, // m-by-k
102 cholmod_common *cc // workspace and parameters
103 )
104 {
105 RETURN_IF_NULL_COMMON (NULL) ;
106 RETURN_IF_NULL (A, NULL) ;
107 RETURN_IF_NULL (B, NULL) ;
108 cc->status = CHOLMOD_OK ;
109
110 return ((A->xtype == CHOLMOD_REAL) ?
111 SuiteSparseQR <double> (ordering, tol, A, B, cc) :
112 SuiteSparseQR <Complex> (ordering, tol, A, B, cc)) ;
113 }
114
115 // =============================================================================
116 // === SuiteSparseQR_C_backslash_default =======================================
117 // =============================================================================
118
119 // X = A\B where B is dense, using default ordering and tol
SuiteSparseQR_C_backslash_default(cholmod_sparse * A,cholmod_dense * B,cholmod_common * cc)120 cholmod_dense *SuiteSparseQR_C_backslash_default // returns X, NULL if failure
121 (
122 // inputs:
123 cholmod_sparse *A, // m-by-n sparse matrix
124 cholmod_dense *B, // m-by-k
125 cholmod_common *cc // workspace and parameters
126 )
127 {
128 return (SuiteSparseQR_C_backslash (SPQR_ORDERING_DEFAULT, SPQR_DEFAULT_TOL,
129 A, B, cc)) ;
130 }
131
132 // =============================================================================
133 // === SuiteSparseQR_C_backslash_sparse ========================================
134 // =============================================================================
135
136 // X = A\B where B is sparse
SuiteSparseQR_C_backslash_sparse(int ordering,double tol,cholmod_sparse * A,cholmod_sparse * B,cholmod_common * cc)137 cholmod_sparse *SuiteSparseQR_C_backslash_sparse // returns X, NULL if failure
138 (
139 // inputs:
140 int ordering, // all, except 3:given treated as 0:fixed
141 double tol, // columns with 2-norm <= tol are treated as 0
142 cholmod_sparse *A, // m-by-n sparse matrix
143 cholmod_sparse *B, // m-by-k
144 cholmod_common *cc // workspace and parameters
145 )
146 {
147 RETURN_IF_NULL_COMMON (NULL) ;
148 RETURN_IF_NULL (A, NULL) ;
149 RETURN_IF_NULL (B, NULL) ;
150 cc->status = CHOLMOD_OK ;
151
152 return ((A->xtype == CHOLMOD_REAL) ?
153 SuiteSparseQR <double> (ordering, tol, A, B, cc) :
154 SuiteSparseQR <Complex> (ordering, tol, A, B, cc)) ;
155 }
156
157 #ifndef NEXPERT
158
159 // =============================================================================
160 // === C wrappers for expert routines ==========================================
161 // =============================================================================
162
163 // =============================================================================
164 // === SuiteSparseQR_C_factorize ===============================================
165 // =============================================================================
166
SuiteSparseQR_C_factorize(int ordering,double tol,cholmod_sparse * A,cholmod_common * cc)167 SuiteSparseQR_C_factorization *SuiteSparseQR_C_factorize
168 (
169 // inputs:
170 int ordering, // all, except 3:given treated as 0:fixed
171 double tol, // columns with 2-norm <= tol are treated as 0
172 cholmod_sparse *A, // m-by-n sparse matrix
173 cholmod_common *cc // workspace and parameters
174 )
175 {
176 RETURN_IF_NULL_COMMON (FALSE) ;
177 RETURN_IF_NULL (A, NULL) ;
178 cc->status = CHOLMOD_OK ;
179
180 SuiteSparseQR_C_factorization *QR ;
181 QR = (SuiteSparseQR_C_factorization *)
182 cholmod_l_malloc (1, sizeof (SuiteSparseQR_C_factorization), cc) ;
183 if (cc->status < CHOLMOD_OK)
184 {
185 return (NULL) ;
186 }
187 QR->xtype = A->xtype ;
188 QR->factors = (A->xtype == CHOLMOD_REAL) ?
189 ((void *) SuiteSparseQR_factorize <double> (ordering, tol, A, cc)) :
190 ((void *) SuiteSparseQR_factorize <Complex> (ordering, tol, A, cc)) ;
191 if (cc->status < CHOLMOD_OK)
192 {
193 SuiteSparseQR_C_free (&QR, cc) ;
194 }
195 return (QR) ;
196 }
197
198 // =============================================================================
199 // === SuiteSparseQR_C_symbolic ================================================
200 // =============================================================================
201
SuiteSparseQR_C_symbolic(int ordering,int allow_tol,cholmod_sparse * A,cholmod_common * cc)202 SuiteSparseQR_C_factorization *SuiteSparseQR_C_symbolic
203 (
204 // inputs:
205 int ordering, // all, except 3:given treated as 0:fixed
206 int allow_tol, // if FALSE, tol is ignored by the numeric
207 // factorization, and no rank detection is performed
208 cholmod_sparse *A, // sparse matrix to factorize (A->x ignored)
209 cholmod_common *cc // workspace and parameters
210 )
211 {
212 RETURN_IF_NULL_COMMON (NULL) ;
213 RETURN_IF_NULL (A, NULL) ;
214 cc->status = CHOLMOD_OK ;
215
216 SuiteSparseQR_C_factorization *QR ;
217 QR = (SuiteSparseQR_C_factorization *)
218 cholmod_l_malloc (1, sizeof (SuiteSparseQR_C_factorization), cc) ;
219 if (cc->status < CHOLMOD_OK)
220 {
221 // out of memory
222 return (NULL) ;
223 }
224 QR->xtype = A->xtype ;
225 QR->factors = (A->xtype == CHOLMOD_REAL) ?
226 ((void *) SuiteSparseQR_symbolic <double> (ordering, allow_tol, A,cc)):
227 ((void *) SuiteSparseQR_symbolic <Complex> (ordering, allow_tol, A,cc));
228 if (cc->status < CHOLMOD_OK)
229 {
230 // out of memory
231 SuiteSparseQR_C_free (&QR, cc) ;
232 }
233 return (QR) ;
234 }
235
236 // =============================================================================
237 // === SuiteSparseQR_C_numeric =================================================
238 // =============================================================================
239
240 // numeric QR factorization; must be preceded by a call to
241 // SuiteSparseQR_C_symbolic.
242
SuiteSparseQR_C_numeric(double tol,cholmod_sparse * A,SuiteSparseQR_C_factorization * QR,cholmod_common * cc)243 int SuiteSparseQR_C_numeric // returns TRUE if successful, FALSE otherwise
244 (
245 // inputs:
246 double tol, // treat columns with 2-norm <= tol as zero
247 cholmod_sparse *A, // sparse matrix to factorize
248 // input/output:
249 SuiteSparseQR_C_factorization *QR,
250 cholmod_common *cc // workspace and parameters
251 )
252 {
253 RETURN_IF_NULL_COMMON (FALSE) ;
254 RETURN_IF_NULL (A, FALSE) ;
255 RETURN_IF_NULL (QR, FALSE) ;
256 cc->status = CHOLMOD_OK ;
257
258 if (QR->xtype == CHOLMOD_REAL)
259 {
260 SuiteSparseQR_factorization <double> *QR2 ;
261 QR2 = (SuiteSparseQR_factorization <double> *) (QR->factors) ;
262 SuiteSparseQR_numeric <double> (tol, A, QR2, cc) ;
263 }
264 else
265 {
266 SuiteSparseQR_factorization <Complex> *QR2 ;
267 QR2 = (SuiteSparseQR_factorization <Complex> *) (QR->factors) ;
268 SuiteSparseQR_numeric <Complex> (tol, A, QR2, cc) ;
269 }
270 return (TRUE) ;
271 }
272
273 // =============================================================================
274 // === SuiteSparseQR_C_free ====================================================
275 // =============================================================================
276
SuiteSparseQR_C_free(SuiteSparseQR_C_factorization ** QR_handle,cholmod_common * cc)277 int SuiteSparseQR_C_free
278 (
279 // input/output:
280 SuiteSparseQR_C_factorization **QR_handle,
281 cholmod_common *cc
282 )
283 {
284 RETURN_IF_NULL_COMMON (FALSE) ;
285
286 SuiteSparseQR_C_factorization *QR ;
287 if (QR_handle == NULL || *QR_handle == NULL) return (TRUE) ;
288 QR = *QR_handle ;
289 if (QR->xtype == CHOLMOD_REAL)
290 {
291 SuiteSparseQR_factorization <double> *QR2 ;
292 QR2 = (SuiteSparseQR_factorization <double> *) (QR->factors) ;
293 spqr_freefac <double> (&QR2, cc) ;
294 }
295 else
296 {
297 SuiteSparseQR_factorization <Complex> *QR2 ;
298 QR2 = (SuiteSparseQR_factorization <Complex> *) (QR->factors) ;
299 spqr_freefac <Complex> (&QR2, cc) ;
300 }
301 cholmod_l_free (1, sizeof (SuiteSparseQR_C_factorization), QR, cc) ;
302 *QR_handle = NULL ;
303 return (TRUE) ;
304 }
305
306 // =============================================================================
307 // === SuiteSparseQR_C_solve ===================================================
308 // =============================================================================
309
310 // Solve an upper or lower triangular system using R from the QR factorization
311 //
312 // system=SPQR_RX_EQUALS_B (0): X = R\B B is m-by-k and X is n-by-k
313 // system=SPQR_RETX_EQUALS_B (1): X = E*(R\B) as above, E is a permutation
314 // system=SPQR_RTX_EQUALS_B (2): X = R'\B B is n-by-k and X is m-by-k
315 // system=SPQR_RTX_EQUALS_ETB (3): X = R'\(E'*B) as above, E is a permutation
316
SuiteSparseQR_C_solve(int system,SuiteSparseQR_C_factorization * QR,cholmod_dense * B,cholmod_common * cc)317 cholmod_dense* SuiteSparseQR_C_solve // returnx X, or NULL if failure
318 (
319 // inputs:
320 int system, // which system to solve
321 SuiteSparseQR_C_factorization *QR, // of an m-by-n sparse matrix A
322 cholmod_dense *B, // right-hand-side, m-by-k or n-by-k
323 cholmod_common *cc // workspace and parameters
324 )
325 {
326 RETURN_IF_NULL (QR, NULL) ;
327 return ((QR->xtype == CHOLMOD_REAL) ?
328 SuiteSparseQR_solve <double> (system,
329 (SuiteSparseQR_factorization <double> *) QR->factors, B, cc) :
330 SuiteSparseQR_solve <Complex> (system,
331 (SuiteSparseQR_factorization <Complex> *) QR->factors, B, cc)) ;
332 }
333
334 // =============================================================================
335 // === SuiteSparseQR_C_qmult ===================================================
336 // =============================================================================
337
338 // Applies Q in Householder form (as stored in the QR factorization object
339 // returned by SuiteSparseQR_C_factorize) to a dense matrix X.
340 //
341 // method SPQR_QTX (0): Y = Q'*X
342 // method SPQR_QX (1): Y = Q*X
343 // method SPQR_XQT (2): Y = X*Q'
344 // method SPQR_XQ (3): Y = X*Q
345
346 // returns Y of size m-by-n, or NULL on failure
SuiteSparseQR_C_qmult(int method,SuiteSparseQR_C_factorization * QR,cholmod_dense * X,cholmod_common * cc)347 cholmod_dense *SuiteSparseQR_C_qmult
348 (
349 // inputs:
350 int method, // 0,1,2,3
351 SuiteSparseQR_C_factorization *QR, // of an m-by-n sparse matrix A
352 cholmod_dense *X, // size m-by-n with leading dimension ldx
353 cholmod_common *cc // workspace and parameters
354 )
355 {
356 RETURN_IF_NULL (QR, NULL) ;
357 return ((QR->xtype == CHOLMOD_REAL) ?
358 SuiteSparseQR_qmult <double> (method,
359 (SuiteSparseQR_factorization <double> *) QR->factors, X, cc) :
360 SuiteSparseQR_qmult <Complex> (method,
361 (SuiteSparseQR_factorization <Complex> *) QR->factors, X, cc)) ;
362 }
363
364 #endif
365
366 // =============================================================================
367 }
368