1 // =============================================================================
2 // === SuiteSparseQR ===========================================================
3 // =============================================================================
4
5 // QR factorization of a sparse matrix, and optionally solve a least squares
6 // problem, Q*R = A*E where A is m-by-n, E is a permutation matrix, R is
7 // upper triangular if A is full rank, and Q is an orthogonal matrix.
8 // R is upper trapezoidal or a "squeezed" upper trapezoidal matrix if A is
9 // found to be rank deficient.
10 //
11 // All output arrays are optional. If you pass NULL for the C, R, E, X,
12 // and/or H array handles, those arrays will not be returned.
13 //
14 // The Z output (either sparse or dense) is either C, C', or X where
15 // C = Q'*B and X = E*(R\(Q'*B)). The latter gives the result of X=A\B,
16 // which is the least-squares solution if m > n.
17 //
18 // To return full-sized results, set econ = m. Then C and R will have m rows,
19 // and C' will have m columns.
20 //
21 // To return economy-sized results, set econ = n. Then C and R will have k
22 // rows and C' will have k columns, where k = min(m,n).
23 //
24 // To return rank-sized results, set econ = 0. Then C and R will have k rows
25 // and C' will have k columns, where k = r = the estimated rank of A.
26 //
27 // In all three cases, k = max (min (m, econ), r).
28 //
29 // To compute Q, pass in B = speye (m), and then set getCTX to 1
30 // (Q is then returned as Zsparse).
31 //
32 // The Householder representation of Q is returned in H, HTau, and HPinv.
33 // E is a permutation vector represention of the column permutation, for
34 // the factorization Q*R=A*E.
35
36
37 #include "spqr.hpp"
38
39 #define XCHUNK 4 // FUTURE: make this a parameter
40
41 #define FREE_ALL \
42 spqr_freefac (&QR, cc) ; \
43 cholmod_l_free (rjsize+1, sizeof (Long), H2p, cc) ; \
44 cholmod_l_free_dense (&HTau, cc) ; \
45 cholmod_l_free_sparse (&H, cc) ; \
46 cholmod_l_free_sparse (&R, cc) ; \
47 cholmod_l_free_sparse (&Xsparse, cc) ; \
48 cholmod_l_free_sparse (&Zsparse, cc) ; \
49 cholmod_l_free_dense (&Zdense, cc) ; \
50 cholmod_l_free (xsize, sizeof (Entry), Xwork, cc) ; \
51 cholmod_l_free (csize, sizeof (Entry), C, cc) ; \
52 cholmod_l_free (wsize, sizeof (Entry), W, cc) ; \
53 cholmod_l_free (maxfrank, sizeof (Long), Rlive, cc) ; \
54 cholmod_l_free (maxfrank, sizeof (Entry *), Rcolp, cc) ; \
55 cholmod_l_free (n+bncols, sizeof (Long), E, cc) ; \
56 cholmod_l_free (m, sizeof (Long), HP1inv, cc) ;
57
58 // returns rank(A) estimate if successful, EMPTY otherwise
SuiteSparseQR(int ordering,double tol,Long econ,int getCTX,cholmod_sparse * A,cholmod_sparse * Bsparse,cholmod_dense * Bdense,cholmod_sparse ** p_Zsparse,cholmod_dense ** p_Zdense,cholmod_sparse ** p_R,Long ** p_E,cholmod_sparse ** p_H,Long ** p_HPinv,cholmod_dense ** p_HTau,cholmod_common * cc)59 template <typename Entry> Long SuiteSparseQR
60 (
61 // inputs, not modified
62 int ordering, // all, except 3:given treated as 0:fixed
63 double tol, // columns with 2-norm <= tol are treated as zero
64 Long econ, // number of rows of C and R to return; a value
65 // less than the rank r of A is treated as r, and
66 // a value greater than m is treated as m.
67 // That is, e = max(min(m,econ),rank(A)) gives the
68 // number of rows of C and R, and columns of C'.
69
70 int getCTX, // if 0: return Z = C of size econ-by-bncols
71 // if 1: return Z = C' of size bncols-by-econ
72 // if 2: return Z = X of size econ-by-bncols
73
74 cholmod_sparse *A, // m-by-n sparse matrix
75
76 // B is either sparse or dense. If Bsparse is non-NULL, B is sparse and
77 // Bdense is ignored. If Bsparse is NULL and Bdense is non-NULL, then B is
78 // dense. B is not present if both are NULL.
79 cholmod_sparse *Bsparse, // B is m-by-bncols
80 cholmod_dense *Bdense,
81
82 // output arrays, neither allocated nor defined on input.
83
84 // Z is the matrix C, C', or X, either sparse or dense. If p_Zsparse is
85 // non-NULL, then Z is returned as a sparse matrix. If p_Zsparse is NULL
86 // and p_Zdense is non-NULL, then Z is returned as a dense matrix. Neither
87 // are returned if both arguments are NULL.
88 cholmod_sparse **p_Zsparse,
89 cholmod_dense **p_Zdense,
90
91 cholmod_sparse **p_R, // the R factor
92 Long **p_E, // size n; fill-reducing ordering of A.
93 cholmod_sparse **p_H, // the Householder vectors (m-by-nh)
94 Long **p_HPinv, // size m; row permutation for H
95 cholmod_dense **p_HTau, // size 1-by-nh, Householder coefficients
96
97 // workspace and parameters
98 cholmod_common *cc
99 )
100 {
101 Long *Q1fill, *R1p, *R1j, *P1inv, *Zp, *Zi, *Rp, *Ri, *Rap, *Hp, *H2p,
102 *Hi, *HP1inv, // *Ap, *Ai,
103 *Rlive, *E, *Bp, *Bi ;
104 Entry *R1x, *B, *Zx, *Rx, *X2, *C, *Hx, *C1, *X1, *Xwork, // *Ax,
105 *W, **Rcolp, *Bx ;
106 Long i, j, k, p, p2, n1cols, n1rows, B_is_sparse, Z_is_sparse, getC,
107 getR, getH, getE, getX, getZ, iold, inew, rank, n2, rnz, znz, zn, zm,
108 pr, xsize, getT, nh, k1, k2, xchunk, xncol, m, n, csize, wsize, ldb,
109 maxfrank, rjsize, bncols, bnrows ;
110 spqr_symbolic *QRsym ;
111 spqr_numeric <Entry> *QRnum ;
112 SuiteSparseQR_factorization <Entry> *QR ;
113 cholmod_sparse *Xsparse, *Zsparse, *R, *H ;
114 cholmod_dense *Zdense, *HTau ;
115
116 double t0 = SuiteSparse_time ( ) ;
117
118 int ok = TRUE ;
119
120 // -------------------------------------------------------------------------
121 // get inputs
122 // -------------------------------------------------------------------------
123
124 if (p_Zsparse != NULL) *p_Zsparse = NULL ;
125 if (p_Zdense != NULL) *p_Zdense = NULL ;
126 if (p_R != NULL) *p_R = NULL ;
127 if (p_E != NULL) *p_E = NULL ;
128 if (p_H != NULL) *p_H = NULL ;
129 if (p_HPinv != NULL) *p_HPinv = NULL ;
130 if (p_HTau != NULL) *p_HTau = NULL ;
131
132 RETURN_IF_NULL_COMMON (EMPTY) ;
133 RETURN_IF_NULL (A, EMPTY) ;
134 Long xtype = spqr_type <Entry> ( ) ;
135 RETURN_IF_XTYPE_INVALID (A, EMPTY) ;
136 if (Bsparse != NULL) RETURN_IF_XTYPE_INVALID (Bsparse, EMPTY) ;
137 if (Bdense != NULL) RETURN_IF_XTYPE_INVALID (Bdense, EMPTY) ;
138 cc->status = CHOLMOD_OK ;
139
140 n = 0 ;
141 m = 0 ;
142 rnz = 0 ;
143 nh = 0 ;
144
145 Rap = NULL ;
146
147 R = NULL ;
148 Rp = NULL ;
149 Ri = NULL ;
150 Rx = NULL ;
151
152 Xsparse = NULL ;
153 Zsparse = NULL ;
154 Zdense = NULL ;
155
156 Zp = NULL ;
157 Zi = NULL ;
158 Zx = NULL ;
159
160 H = NULL ;
161 Hp = NULL ;
162 Hi = NULL ;
163 Hx = NULL ;
164 H2p = NULL ;
165
166 HTau = NULL ;
167 HP1inv = NULL ;
168
169 maxfrank = 0 ;
170 xsize = 0 ;
171 csize = 0 ;
172 wsize = 0 ;
173
174 Xwork = NULL ;
175 C = NULL ;
176 Rcolp = NULL ;
177 Rlive = NULL ;
178
179 E = NULL ;
180 W = NULL ;
181
182 m = A->nrow ;
183 n = A->ncol ;
184 // Ap = (Long *) A->p ;
185 // Ai = (Long *) A->i ;
186 // Ax = (Entry *) A->x ;
187
188 // B is an optional input. It can be sparse or dense
189 B_is_sparse = (Bsparse != NULL) ;
190 if (B_is_sparse)
191 {
192 // B is sparse
193 bncols = Bsparse->ncol ;
194 bnrows = Bsparse->nrow ;
195 Bp = (Long *) Bsparse->p ;
196 Bi = (Long *) Bsparse->i ;
197 Bx = (Entry *) Bsparse->x ;
198 ldb = 0 ; // unused
199 }
200 else if (Bdense != NULL)
201 {
202 // B is dense
203 bncols = Bdense->ncol ;
204 bnrows = Bdense->nrow ;
205 Bp = NULL ;
206 Bi = NULL ;
207 Bx = (Entry *) Bdense->x ;
208 ldb = Bdense->d ;
209 }
210 else
211 {
212 // B is not present
213 bncols = 0 ;
214 bnrows = m ;
215 Bp = NULL ;
216 Bi = NULL ;
217 Bx = NULL ;
218 ldb = 0 ;
219 }
220
221 if (bnrows != m)
222 {
223 ERROR (CHOLMOD_INVALID, "A and B must have the same # of rows") ;
224 return (EMPTY) ;
225 }
226 PR (("bncols : %ld\n", bncols)) ;
227
228 // Z = C, C', or X, can be sparse or dense
229 Z_is_sparse = (p_Zsparse != NULL) ;
230 getZ = (Z_is_sparse || p_Zdense != NULL) && (getCTX >= 0 && getCTX <= 2) ;
231
232 // Z = C is an optional output of size econ-by-bncols, sparse or dense
233 getC = getZ && (getCTX == 0) ;
234
235 // Z = C' is an optional output of size bncols-by-econ, sparse or dense
236 getT = getZ && (getCTX == 1) ;
237
238 // Z = X is an optional output of size n-by-bncols, sparse or dense
239 getX = getZ && (getCTX == 2) ;
240
241 // R is an optional output of size econ-by-n. It is always sparse
242 getR = (p_R != NULL) ;
243
244 // E is an optional output; it is a permutation vector of length n+bncols
245 getE = (p_E != NULL) ;
246
247 // H is an optional output, of size m-by-nh. It is always sparse
248 getH = (p_H != NULL && p_HPinv != NULL && p_HTau != NULL) ;
249
250 // at most one of C, C', or X will be returned as the Z matrix
251 ASSERT (getC + getT + getX <= 1) ;
252
253 PR (("getCTX %d getZ %ld getC %ld getT %ld getX %ld getR %ld getE %ld"
254 "getH %ld\n", getCTX, getZ, getC, getT, getX, getR, getE, getH)) ;
255
256 // -------------------------------------------------------------------------
257 // symbolic and numeric QR factorization of [A B], exploiting singletons
258 // -------------------------------------------------------------------------
259
260 QR = spqr_1factor (ordering, tol, bncols, getH, A, ldb, Bp, Bi, Bx, cc) ;
261
262 if (cc->status < CHOLMOD_OK)
263 {
264 // out of memory
265 return (EMPTY) ;
266 }
267
268 R1p = QR->R1p ;
269 R1j = QR->R1j ;
270 R1x = QR->R1x ;
271 P1inv = QR->P1inv ;
272 Q1fill = QR->Q1fill ;
273 QRsym = QR->QRsym ;
274 QRnum = QR->QRnum ;
275 n1rows = QR->n1rows ;
276 n1cols = QR->n1cols ;
277 n2 = n - n1cols ;
278 rank = QR->rank ;
279 tol = QR->tol ;
280
281 PR (("Singletons: n1cols %ld n1rows %ld\n", n1cols, n1rows)) ;
282 PR (("n2 %ld n %ld QRsym->n %ld\n", n2, n, QRsym->n)) ;
283
284 // -------------------------------------------------------------------------
285 // determine the economy size
286 // -------------------------------------------------------------------------
287
288 // econ gives the number of rows of C and R to return. It must be at least
289 // as large as the estimated rank of R (all entries in R are always
290 // returned). It can be no larger than m (all of C and R are returned).
291 // If X is requested but not C nor R, then only rows 0 to rank-1 of C is
292 // extracted from the QR factorization.
293
294 if (getX && !getR)
295 {
296 econ = rank ;
297 }
298 else
299 {
300 // constrain econ to be between rank and m
301 econ = MIN (m, econ) ;
302 econ = MAX (econ, rank) ;
303 }
304 PR (("getX %ld econ %ld\n", getX, econ)) ;
305
306 zm = getT ? bncols : econ ; // number of rows of Z, if present
307 zn = getT ? econ : bncols ; // number of columns of Z, if present
308
309 ASSERT (rank <= econ && econ <= m) ;
310
311 // -------------------------------------------------------------------------
312 // count the entries in the R factor of Y
313 // -------------------------------------------------------------------------
314
315 if (getZ)
316 {
317 // Zsparse is zm-by-zn, but with no entries so far
318 Zsparse = cholmod_l_allocate_sparse (zm, zn, 0, TRUE, TRUE, 0, xtype,
319 cc) ;
320 Zp = Zsparse ? ((Long *) Zsparse->p) : NULL ;
321 PR (("Z is zm %ld by zn %ld\n", zm, zn)) ;
322 }
323
324 if (getR)
325 {
326 // R is econ-by-n, but with no entries so far
327 R = cholmod_l_allocate_sparse (econ, n, 0, TRUE, TRUE, 0, xtype, cc) ;
328 Rp = R ? ((Long *) R->p) : NULL ;
329 Rap = Rp + n1cols ;
330 }
331
332 rjsize = QRsym->rjsize ; // rjsize is an upper bound on nh
333
334 if (getH)
335 {
336 H2p = (Long *) cholmod_l_malloc (rjsize+1, sizeof (Long), cc) ;
337 }
338
339 if (cc->status < CHOLMOD_OK)
340 {
341 // out of memory
342 FREE_ALL ;
343 return (EMPTY) ;
344 }
345
346 // count the number of nonzeros in each column of the R factor of Y,
347 // and count the entries in the Householder vectors
348 spqr_rcount (QRsym, QRnum, n1rows, econ, n2, getT, Rap, Zp, H2p, &nh) ;
349
350 #ifndef NDEBUG
351 if (getZ) for (k = 0 ; k < zn ; k++)
352 {
353 PR (("Zp [%ld] %ld count\n", k, Zp [k])) ;
354 }
355 if (getR) for (k = 0 ; k < n ; k++)
356 {
357 PR (("R factor count Rp [%ld] = %ld\n", k, Rp [k])) ;
358 }
359 #endif
360
361 // -------------------------------------------------------------------------
362 // count the entries in the singleton rows, R1
363 // -------------------------------------------------------------------------
364
365 if (getR)
366 {
367 // add the entries in the singleton rows R11 and R12
368 for (k = 0 ; k < n1rows ; k++)
369 {
370 for (p = R1p [k] ; p < R1p [k+1] ; p++)
371 {
372 j = R1j [p] ;
373 ASSERT (j >= k && j < n) ;
374 Rp [j]++ ;
375 }
376 }
377 }
378
379 #ifndef NDEBUG
380 if (getR) for (k = 0 ; k < n ; k++)
381 {
382 PR (("R factor count Rp [%ld] = %ld after R1\n", k, Rp [k])) ;
383 }
384 #endif
385
386 // -------------------------------------------------------------------------
387 // count the entries in B1 and add them to the counts of Z
388 // -------------------------------------------------------------------------
389
390 if (getZ && n1rows > 0)
391 {
392 // add the entries in the singleton rows of B1
393 if (B_is_sparse)
394 {
395 // B is sparse
396 for (k = 0 ; k < bncols ; k++)
397 {
398 // count the nonzero entries in column k of B1
399 for (p = Bp [k] ; p < Bp [k+1] ; p++)
400 {
401 iold = Bi [p] ;
402 inew = P1inv [iold] ;
403 if (inew < n1rows)
404 {
405 Zp [getT ? inew : k]++ ;
406 }
407 }
408 }
409 }
410 else
411 {
412 // B is dense
413 B = Bx ;
414 for (k = 0 ; k < bncols ; k++)
415 {
416 // count the nonzero entries in column k of B1
417 for (iold = 0 ; iold < m ; iold++)
418 {
419 inew = P1inv [iold] ;
420 if (inew < n1rows && B [iold] != (Entry) 0)
421 {
422 Zp [getT ? inew : k]++ ;
423 }
424 }
425 B += ldb ;
426 }
427 }
428 #ifndef NDEBUG
429 for (k = 0 ; k < zn ; k++)
430 {
431 PR (("after adding singletons Zp [%ld] %ld count\n", k, Zp [k])) ;
432 }
433 #endif
434 }
435
436 // -------------------------------------------------------------------------
437 // compute the Rp and Zp column pointers (skip if NULL)
438 // -------------------------------------------------------------------------
439
440 // no Long overflow can occur
441 rnz = spqr_cumsum (n, Rp) ; // Rp = cumsum ([0 Rp])
442 znz = spqr_cumsum (zn, Zp) ; // Zp = cumsum ([0 Zp])
443
444 // -------------------------------------------------------------------------
445 // allocate R, Z, and H
446 // -------------------------------------------------------------------------
447
448 if (getR)
449 {
450 // R now has space for rnz entries
451 cholmod_l_reallocate_sparse (rnz, R, cc) ;
452 Ri = (Long *) R->i ;
453 Rx = (Entry *) R->x ;
454 }
455
456 if (getZ)
457 {
458 // Zsparse now has space for znz entries
459 cholmod_l_reallocate_sparse (znz, Zsparse, cc) ;
460 Zi = (Long *) Zsparse->i ;
461 Zx = (Entry *) Zsparse->x ;
462 }
463
464 if (getH)
465 {
466 // H is m-by-nh with hnz entries, where nh <= rjsize
467 Long hnz = H2p [nh] ;
468 H = cholmod_l_allocate_sparse (m, nh, hnz, TRUE, TRUE, 0, xtype, cc) ;
469 // copy the column pointers from H2p to Hp
470 if (cc->status == CHOLMOD_OK)
471 {
472 Hp = (Long *) H->p ;
473 Hi = (Long *) H->i ;
474 Hx = (Entry *) H->x ;
475 for (k = 0 ; k <= nh ; k++)
476 {
477 Hp [k] = H2p [k] ;
478 }
479 }
480 // free the H2p workspace, and allocate HTau
481 cholmod_l_free (rjsize+1, sizeof (Long), H2p, cc) ;
482 H2p = NULL ;
483 HTau = cholmod_l_allocate_dense (1, nh, 1, xtype, cc) ;
484 }
485
486 if (cc->status < CHOLMOD_OK)
487 {
488 // out of memory
489 FREE_ALL ;
490 return (EMPTY) ;
491 }
492
493 // -------------------------------------------------------------------------
494 // place the singleton rows in R
495 // -------------------------------------------------------------------------
496
497 if (getR)
498 {
499 // add the entries in the singleton rows R11 and R12
500 for (k = 0 ; k < n1rows ; k++)
501 {
502 for (p = R1p [k] ; p < R1p [k+1] ; p++)
503 {
504 j = R1j [p] ;
505 ASSERT (j >= k && j < n) ;
506 pr = Rp [j]++ ;
507 Ri [pr] = k ;
508 Rx [pr] = R1x [p] ;
509 }
510 }
511 }
512
513 // -------------------------------------------------------------------------
514 // place the B1 entries in C or C'
515 // -------------------------------------------------------------------------
516
517 if (getZ && n1rows > 0)
518 {
519 if (B_is_sparse)
520 {
521 // B is sparse
522 PR (("B1 sparse:\n")) ;
523 for (k = 0 ; k < bncols ; k++)
524 {
525 // copy the nonzero entries in column k of B1 into C or C'
526 for (p = Bp [k] ; p < Bp [k+1] ; p++)
527 {
528 iold = Bi [p] ;
529 inew = P1inv [iold] ;
530 if (inew < n1rows)
531 {
532 if (getT)
533 {
534 p2 = Zp [inew]++ ;
535 Zi [p2] = k ;
536 Zx [p2] = spqr_conj (Bx [p]) ;
537 }
538 else
539 {
540 p2 = Zp [k]++ ;
541 Zi [p2] = inew ;
542 Zx [p2] = Bx [p] ;
543 }
544 }
545 }
546 }
547 }
548 else
549 {
550 // B is dense
551 B = Bx ;
552 PR (("B1 dense:\n")) ;
553 for (k = 0 ; k < bncols ; k++)
554 {
555 // copy the nonzero entries in column k of B1 into C or C'
556 for (iold = 0 ; iold < m ; iold++)
557 {
558 inew = P1inv [iold] ;
559 if (inew < n1rows && B [iold] != (Entry) 0)
560 {
561 if (getT)
562 {
563 p2 = Zp [inew]++ ;
564 Zi [p2] = k ;
565 Zx [p2] = spqr_conj (B [iold]) ;
566 }
567 else
568 {
569 p2 = Zp [k]++ ;
570 Zi [p2] = inew ;
571 Zx [p2] = B [iold] ;
572 }
573 }
574 }
575 B += ldb ;
576 }
577 }
578 }
579
580 // -------------------------------------------------------------------------
581 // place the R factor entries in Z and R, and get H
582 // -------------------------------------------------------------------------
583
584 spqr_rconvert (QRsym, QRnum, n1rows, econ, n2, getT, Rap, Ri, Rx,
585 Zp, Zi, Zx, Hp, Hi, Hx, HTau ? ((Entry *) HTau->x) : NULL) ;
586
587 // -------------------------------------------------------------------------
588 // finalize the column pointers of Z and R
589 // -------------------------------------------------------------------------
590
591 spqr_shift (n, Rp) ;
592 spqr_shift (zn, Zp) ;
593
594 // -------------------------------------------------------------------------
595 // Finalize Z: either Z = E*(R\Z) or Z = full (Z)
596 // -------------------------------------------------------------------------
597
598 if (getX)
599 {
600
601 // ---------------------------------------------------------------------
602 // solve the least squares problem: X = E*(R\C)
603 // ---------------------------------------------------------------------
604
605 // result X is n-by-bncols with leading dimension n, sparse or dense
606 maxfrank = QRnum->maxfrank ;
607
608 if (Z_is_sparse)
609 {
610
611 // -----------------------------------------------------------------
612 // create the sparse result X
613 // -----------------------------------------------------------------
614
615 // ask for space for n+1 entries; this is increased later if needed
616 Xsparse = cholmod_l_allocate_sparse (n, bncols, n+1, TRUE, TRUE, 0,
617 xtype, cc) ;
618 xncol = 0 ;
619
620 if (cc->status < CHOLMOD_OK)
621 {
622 // out of memory
623 FREE_ALL ;
624 return (EMPTY) ;
625 }
626
627 // Xwork will be a temporary dense result
628 xchunk = MIN (bncols, XCHUNK) ;
629
630 // xsize = n * xchunk ;
631 xsize = spqr_mult (n, xchunk, &ok) ;
632
633 // csize = rank * xchunk ;
634 csize = spqr_mult (rank, xchunk, &ok) ;
635
636 // wsize = maxfrank * xchunk ;
637 wsize = spqr_mult (maxfrank, xchunk, &ok) ;
638
639 if (ok)
640 {
641 Xwork = (Entry *) cholmod_l_malloc (xsize, sizeof (Entry), cc) ;
642 C = (Entry *) cholmod_l_calloc (csize, sizeof (Entry), cc) ;
643 W = (Entry *) cholmod_l_malloc (wsize, sizeof (Entry), cc) ;
644 }
645
646 // -----------------------------------------------------------------
647 // punt: use less workspace if we ran out of memory
648 // -----------------------------------------------------------------
649
650 if (!ok || (cc->status < CHOLMOD_OK && xchunk > 1))
651 {
652 // PUNT: out of memory; try again with xchunk = 1
653 cc->status = CHOLMOD_OK ;
654 ok = TRUE ;
655 cholmod_l_free (xsize, sizeof (Entry), Xwork, cc) ;
656 cholmod_l_free (csize, sizeof (Entry), C, cc) ;
657 cholmod_l_free (wsize, sizeof (Entry), W, cc) ;
658 xchunk = 1 ;
659 xsize = n ;
660 csize = rank ;
661 wsize = maxfrank ;
662 Xwork = (Entry *) cholmod_l_malloc (xsize, sizeof (Entry), cc) ;
663 C = (Entry *) cholmod_l_calloc (csize, sizeof (Entry), cc) ;
664 W = (Entry *) cholmod_l_malloc (wsize, sizeof (Entry), cc) ;
665 }
666
667 // -----------------------------------------------------------------
668 // use Xwork for the solve
669 // -----------------------------------------------------------------
670
671 X2 = Xwork ;
672
673 }
674 else
675 {
676
677 // -----------------------------------------------------------------
678 // construct the dense X
679 // -----------------------------------------------------------------
680
681 xchunk = bncols ;
682
683 // csize = rank * xchunk ;
684 csize = spqr_mult (rank, xchunk, &ok) ;
685
686 // wsize = maxfrank * xchunk ;
687 wsize = spqr_mult (maxfrank, xchunk, &ok) ;
688
689 if (ok)
690 {
691 C = (Entry *) cholmod_l_calloc (csize, sizeof (Entry), cc) ;
692 W = (Entry *) cholmod_l_malloc (wsize, sizeof (Entry), cc) ;
693 }
694
695 // allocate the dense X and use it for the solve
696 Zdense = cholmod_l_allocate_dense (n, bncols, n, xtype, cc) ;
697 X2 = Zdense ? ((Entry *) Zdense->x) : NULL ;
698 }
699
700 Rlive = (Long *) cholmod_l_malloc (maxfrank, sizeof (Long), cc) ;
701 Rcolp = (Entry **) cholmod_l_malloc (maxfrank, sizeof (Entry *), cc) ;
702
703 if (!ok || cc->status < CHOLMOD_OK)
704 {
705 // still out of memory (or problem too large)
706 FREE_ALL ;
707 ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
708 return (EMPTY) ;
709 }
710
711 for (k1 = 0 ; k1 < bncols ; k1 += xchunk)
712 {
713
714 // -----------------------------------------------------------------
715 // this iteration computes X (:,k1:k2-1)
716 // -----------------------------------------------------------------
717
718 k2 = MIN (bncols, k1+xchunk) ;
719 PR (("k1 %ld k2 %ld bncols %ld n1rows %ld\n",
720 k1, k2, bncols, n1rows)) ;
721
722 #ifndef NDEBUG
723 for (i = 0 ; i < rank*(k2-k1) ; i++)
724 {
725 ASSERT (C [i] == (Entry) 0) ;
726 }
727 #endif
728
729 // -----------------------------------------------------------------
730 // scatter C
731 // -----------------------------------------------------------------
732
733 C1 = C ;
734 for (k = k1 ; k < k2 ; k++)
735 {
736 // scatter C (0:rank-1,k1:k2-1) into workspace
737 for (p = Zp [k] ; p < Zp [k+1] ; p++)
738 {
739 i = Zi [p] ;
740 if (i < rank)
741 {
742 C1 [i] = Zx [p] ;
743 }
744 }
745 C1 += rank ;
746 }
747
748 // -----------------------------------------------------------------
749 // X = E*(R\C)
750 // -----------------------------------------------------------------
751
752 spqr_rsolve (QR, TRUE, k2-k1, rank, C, X2, Rcolp, Rlive, W, cc) ;
753
754 // -----------------------------------------------------------------
755 // clear workspace C ; skip if this is the last k (C is freed)
756 // -----------------------------------------------------------------
757
758 if (k2 < bncols)
759 {
760 C1 = C ;
761 for (k = k1 ; k < k2 ; k++)
762 {
763 for (p = Zp [k] ; p < Zp [k+1] ; p++)
764 {
765 i = Zi [p] ;
766 if (i < rank)
767 {
768 C1 [i] = 0 ;
769 }
770 }
771 C1 += rank ;
772 }
773 }
774
775 // -----------------------------------------------------------------
776 // save the result, either sparse or dense
777 // -----------------------------------------------------------------
778
779 if (Z_is_sparse)
780 {
781 // append the columns of X2 onto the sparse X
782 X1 = X2 ;
783 for (k = 0 ; k < k2-k1 ; k++)
784 {
785 spqr_append (X1, NULL, Xsparse, &xncol, cc) ;
786 X1 += n ;
787 if (cc->status < CHOLMOD_OK)
788 {
789 // out of memory
790 FREE_ALL ;
791 return (EMPTY) ;
792 }
793 }
794 }
795 }
796
797 // ---------------------------------------------------------------------
798 // free workspace
799 // ---------------------------------------------------------------------
800
801 C = (Entry *) cholmod_l_free (csize, sizeof (Entry), C, cc) ;
802 W = (Entry *) cholmod_l_free (wsize, sizeof (Entry), W, cc) ;
803 Rlive = (Long *) cholmod_l_free (maxfrank, sizeof (Long), Rlive, cc);
804 Rcolp = (Entry **) cholmod_l_free (maxfrank, sizeof (Entry *), Rcolp,
805 cc) ;
806
807 // ---------------------------------------------------------------------
808 // free the sparse Z
809 // ---------------------------------------------------------------------
810
811 cholmod_l_free_sparse (&Zsparse, cc) ;
812
813 // ---------------------------------------------------------------------
814 // finalize the sparse X
815 // ---------------------------------------------------------------------
816
817 if (Z_is_sparse)
818 {
819
820 // -----------------------------------------------------------------
821 // reduce X in size so that nnz(X) == nzmax(X)
822 // -----------------------------------------------------------------
823
824 znz = cholmod_l_nnz (Xsparse, cc) ;
825 cholmod_l_reallocate_sparse (znz, Xsparse, cc) ;
826 ASSERT (cc->status == CHOLMOD_OK) ;
827
828 // -----------------------------------------------------------------
829 // Z becomes the sparse X
830 // -----------------------------------------------------------------
831
832 Zsparse = Xsparse ;
833 Xsparse = NULL ;
834
835 // -----------------------------------------------------------------
836 // free the dense Xwork
837 // -----------------------------------------------------------------
838
839 cholmod_l_free (xsize, sizeof (Entry), Xwork, cc) ;
840 Xwork = NULL ;
841 xsize = 0 ;
842 }
843
844 }
845 else if ((getC || getT) && !Z_is_sparse)
846 {
847
848 // ---------------------------------------------------------------------
849 // convert C or C' to full
850 // ---------------------------------------------------------------------
851
852 Zdense = cholmod_l_sparse_to_dense (Zsparse, cc) ;
853 cholmod_l_free_sparse (&Zsparse, cc) ;
854
855 if (cc->status < CHOLMOD_OK)
856 {
857 // out of memory
858 FREE_ALL ;
859 return (EMPTY) ;
860 }
861 }
862
863 // -------------------------------------------------------------------------
864 // extract permutations, if requested, from the QR factors
865 // -------------------------------------------------------------------------
866
867 if (getE)
868 {
869 // preserve Q1fill if E has been requested
870 E = QR->Q1fill ;
871 QR->Q1fill = NULL ;
872 }
873
874 if (getH)
875 {
876 // preserve HP1inv if E has been requested. If there are no singletons,
877 // the QR->HP1inv is NULL, and QR->QRnum->HPinv serves in its place.
878 if (n1cols > 0)
879 {
880 // singletons are present, get E from QR->HP1inv
881 HP1inv = QR->HP1inv ;
882 QR->HP1inv = NULL ;
883 }
884 else
885 {
886 // no singletons are present, get E from QR->QRnum->HPinv
887 ASSERT (QR->HP1inv == NULL) ;
888 HP1inv = QRnum->HPinv ;
889 QRnum->HPinv = NULL ;
890 }
891 }
892
893 // -------------------------------------------------------------------------
894 // free the QR factors
895 // -------------------------------------------------------------------------
896
897 spqr_freefac (&QR, cc) ;
898
899 // -------------------------------------------------------------------------
900 // convert R to trapezoidal, if rank < n and ordering is not fixed
901 // -------------------------------------------------------------------------
902
903 if (getR && ordering != SPQR_ORDERING_FIXED && rank < n && tol >= 0)
904 {
905 Long *Rtrapp, *Rtrapi, *Qtrap ;
906 Entry *Rtrapx ;
907
908 // find Rtrap and Qtrap. This may fail if tol < 0 and the matrix
909 // is structurally rank deficient; in that case, k2 = EMPTY and
910 // Rtrap is returned NULL.
911 k2 = spqr_trapezoidal (n, Rp, Ri, Rx, bncols, Q1fill, TRUE,
912 &Rtrapp, &Rtrapi, &Rtrapx, &Qtrap, cc) ;
913
914 if (cc->status < CHOLMOD_OK)
915 {
916 // out of memory
917 // This is very unlikely to fail since the QR factorization object
918 // has just been freed
919 FREE_ALL ;
920 return (EMPTY) ;
921 }
922
923 ASSERT (k2 == EMPTY || rank == k2) ;
924
925 if (Rtrapp != NULL)
926 {
927 // Rtrap was created (it was skipped if R was already in upper
928 // trapezoidal form)
929
930 // free the old R and Q1fill
931 cholmod_l_free (n+1, sizeof (Long), Rp, cc) ;
932 cholmod_l_free (rnz, sizeof (Long), Ri, cc) ;
933 cholmod_l_free (rnz, sizeof (Entry), Rx, cc) ;
934 cholmod_l_free (n+bncols, sizeof (Long), E, cc) ;
935
936 // replace R and Q1fill with Rtrap and Qtrap
937 R->p = Rtrapp ;
938 R->i = Rtrapi ;
939 R->x = Rtrapx ;
940 E = Qtrap ;
941 }
942 }
943
944 // -------------------------------------------------------------------------
945 // return results
946 // -------------------------------------------------------------------------
947
948 if (getZ)
949 {
950 // return C, C', or X
951 if (Z_is_sparse)
952 {
953 *p_Zsparse = Zsparse ;
954 }
955 else
956 {
957 *p_Zdense = Zdense ;
958 }
959 }
960
961 if (getH)
962 {
963 // return H, of size m-by-nh in sparse column-oriented form
964 *p_H = H ;
965 *p_HTau = HTau ;
966 *p_HPinv = HP1inv ; // this is QR->HP1inv or QR->QRnum->HPinv
967 }
968
969 if (getR)
970 {
971 // return R, of size econ-by-n in sparse column-oriented form
972 *p_R = R ;
973 }
974
975 if (getE)
976 {
977 // return E of size n+bncols (NULL if ordering == SPQR_ORDERING_FIXED)
978 *p_E = E ;
979 }
980
981 double t3 = SuiteSparse_time ( ) ;
982 double total_time = t3 - t0 ;
983 cc->SPQR_solve_time =
984 total_time - cc->SPQR_analyze_time - cc->SPQR_factorize_time ;
985
986 return (rank) ;
987 }
988
989
990 template Long SuiteSparseQR <double>
991 (
992 // inputs, not modified
993 int ordering, // all, except 3:given treated as 0:fixed
994 double tol, // columns with 2-norm <= tol are treated as zero
995 Long econ, // number of rows of C and R to return; a value
996 // less than the rank r of A is treated as r, and
997 // a value greater than m is treated as m.
998
999 int getCTX, // if 0: return Z = C of size econ-by-bncols
1000 // if 1: return Z = C' of size bncols-by-econ
1001 // if 2: return Z = X of size econ-by-bncols
1002
1003 cholmod_sparse *A, // m-by-n sparse matrix
1004 cholmod_sparse *Bsparse,
1005 cholmod_dense *Bdense,
1006
1007 // output arrays, neither allocated nor defined on input.
1008
1009 // Z is the matrix C, C', or X
1010 cholmod_sparse **p_Zsparse,
1011 cholmod_dense **p_Zdense,
1012 cholmod_sparse **p_R, // the R factor
1013 Long **p_E, // size n; fill-reducing ordering of A.
1014 cholmod_sparse **p_H, // the Householder vectors (m-by-nh)
1015 Long **p_HPinv, // size m; row permutation for H
1016 cholmod_dense **p_HTau, // size 1-by-nh, Householder coefficients
1017
1018 // workspace and parameters
1019 cholmod_common *cc
1020 ) ;
1021
1022 template Long SuiteSparseQR <Complex>
1023 (
1024 // inputs, not modified
1025 int ordering, // all, except 3:given treated as 0:fixed
1026 double tol, // columns with 2-norm <= tol are treated as zero
1027 Long econ, // number of rows of C and R to return; a value
1028 // less than the rank r of A is treated as r, and
1029 // a value greater than m is treated as m.
1030
1031 int getCTX, // if 0: return Z = C of size econ-by-bncols
1032 // if 1: return Z = C' of size bncols-by-econ
1033 // if 2: return Z = X of size econ-by-bncols
1034
1035 cholmod_sparse *A, // m-by-n sparse matrix
1036 cholmod_sparse *Bsparse,
1037 cholmod_dense *Bdense,
1038
1039 // output arrays, neither allocated nor defined on input.
1040
1041 // Z is the matrix C, C', or X
1042 cholmod_sparse **p_Zsparse,
1043 cholmod_dense **p_Zdense,
1044 cholmod_sparse **p_R, // the R factor
1045 Long **p_E, // size n; fill-reducing ordering of A.
1046 cholmod_sparse **p_H, // the Householder vectors (m-by-nh)
1047 Long **p_HPinv, // size m; row permutation for H
1048 cholmod_dense **p_HTau, // size 1-by-nh, Householder coefficients
1049
1050 // workspace and parameters
1051 cholmod_common *cc
1052 ) ;
1053
1054
1055 // =============================================================================
1056 // === SuiteSparseQR overloaded functions ======================================
1057 // =============================================================================
1058
1059 // -----------------------------------------------------------------------------
1060 // X=A\B where X and B are dense
1061 // -----------------------------------------------------------------------------
1062
SuiteSparseQR(int ordering,double tol,cholmod_sparse * A,cholmod_dense * B,cholmod_common * cc)1063 template <typename Entry> cholmod_dense *SuiteSparseQR
1064 (
1065 int ordering, // all, except 3:given treated as 0:fixed
1066 double tol, // columns with 2-norm <= tol are treated as zero
1067 cholmod_sparse *A, // m-by-n sparse matrix
1068 cholmod_dense *B, // m-by-nrhs
1069 cholmod_common *cc // workspace and parameters
1070 )
1071 {
1072 cholmod_dense *X ;
1073 SuiteSparseQR <Entry> (ordering, tol, 0, 2, A,
1074 NULL, B, NULL, &X, NULL, NULL, NULL, NULL, NULL, cc) ;
1075 return (X) ;
1076 }
1077
1078 template cholmod_dense *SuiteSparseQR <double>
1079 (
1080 int ordering, // all, except 3:given treated as 0:fixed
1081 double tol, // columns with 2-norm <= tol are treated as zero
1082 cholmod_sparse *A, // m-by-n sparse matrix
1083 cholmod_dense *B, // m-by-nrhs
1084 cholmod_common *cc // workspace and parameters
1085 ) ;
1086
1087 template cholmod_dense *SuiteSparseQR <Complex>
1088 (
1089 int ordering, // all, except 3:given treated as 0:fixed
1090 double tol, // columns with 2-norm <= tol are treated as zero
1091 cholmod_sparse *A, // m-by-n sparse matrix
1092 cholmod_dense *B, // m-by-nrhs
1093 cholmod_common *cc // workspace and parameters
1094 ) ;
1095
1096 // -----------------------------------------------------------------------------
1097 // X=A\B where X and B are dense, default ordering and tol
1098 // -----------------------------------------------------------------------------
1099
SuiteSparseQR(cholmod_sparse * A,cholmod_dense * B,cholmod_common * cc)1100 template <typename Entry> cholmod_dense *SuiteSparseQR
1101 (
1102 cholmod_sparse *A, // m-by-n sparse matrix
1103 cholmod_dense *B, // m-by-nrhs
1104 cholmod_common *cc // workspace and parameters
1105 )
1106 {
1107 cholmod_dense *X ;
1108 SuiteSparseQR <Entry> (SPQR_ORDERING_DEFAULT, SPQR_DEFAULT_TOL, 0, 2, A,
1109 NULL, B, NULL, &X, NULL, NULL, NULL, NULL, NULL, cc) ;
1110 return (X) ;
1111 }
1112
1113 template cholmod_dense *SuiteSparseQR <double>
1114 (
1115 cholmod_sparse *A, // m-by-n sparse matrix
1116 cholmod_dense *B, // m-by-nrhs
1117 cholmod_common *cc // workspace and parameters
1118 ) ;
1119
1120 template cholmod_dense *SuiteSparseQR <Complex>
1121 (
1122 cholmod_sparse *A, // m-by-n sparse matrix
1123 cholmod_dense *B, // m-by-nrhs
1124 cholmod_common *cc // workspace and parameters
1125 ) ;
1126
1127
1128 // -----------------------------------------------------------------------------
1129 // X=A\B where X and B are sparse
1130 // -----------------------------------------------------------------------------
1131
SuiteSparseQR(int ordering,double tol,cholmod_sparse * A,cholmod_sparse * B,cholmod_common * cc)1132 template <typename Entry> cholmod_sparse *SuiteSparseQR
1133 (
1134 int ordering, // all, except 3:given treated as 0:fixed
1135 double tol, // columns with 2-norm <= tol are treated as zero
1136 cholmod_sparse *A, // m-by-n sparse matrix
1137 cholmod_sparse *B, // m-by-nrhs
1138 cholmod_common *cc // workspace and parameters
1139 )
1140 {
1141
1142 cholmod_sparse *X ;
1143 SuiteSparseQR <Entry> (ordering, tol, 0, 2, A,
1144 B, NULL, &X, NULL, NULL, NULL, NULL, NULL, NULL, cc) ;
1145
1146 return (X) ;
1147 }
1148
1149 template cholmod_sparse *SuiteSparseQR <double>
1150 (
1151 int ordering, // all, except 3:given treated as 0:fixed
1152 double tol, // columns with 2-norm <= tol are treated as zero
1153 cholmod_sparse *A, // m-by-n sparse matrix
1154 cholmod_sparse *B, // m-by-nrhs
1155 cholmod_common *cc // workspace and parameters
1156 ) ;
1157
1158 template cholmod_sparse *SuiteSparseQR <Complex>
1159 (
1160 int ordering, // all, except 3:given treated as 0:fixed
1161 double tol, // columns with 2-norm <= tol are treated as zero
1162 cholmod_sparse *A, // m-by-n sparse matrix
1163 cholmod_sparse *B, // m-by-nrhs
1164 cholmod_common *cc // workspace and parameters
1165 ) ;
1166
1167
1168 // -----------------------------------------------------------------------------
1169 // [Q,R,E] = qr(A), returning Q as a sparse matrix
1170 // -----------------------------------------------------------------------------
1171
SuiteSparseQR(int ordering,double tol,Long econ,cholmod_sparse * A,cholmod_sparse ** Q,cholmod_sparse ** R,Long ** E,cholmod_common * cc)1172 template <typename Entry> Long SuiteSparseQR // returns rank(A) estimate
1173 (
1174 int ordering, // all, except 3:given treated as 0:fixed
1175 double tol, // columns with 2-norm <= tol are treated as zero
1176 Long econ, // e = max(min(m,econ),rank(A))
1177 cholmod_sparse *A, // m-by-n sparse matrix
1178 // outputs
1179 cholmod_sparse **Q, // m-by-e sparse matrix
1180 cholmod_sparse **R, // e-by-n sparse matrix
1181 Long **E, // permutation of 0:n-1, NULL if identity
1182 cholmod_common *cc // workspace and parameters
1183 )
1184 {
1185 cholmod_sparse *I ;
1186 Long xtype = spqr_type <Entry> ( ) ;
1187 RETURN_IF_NULL_COMMON (EMPTY) ;
1188 RETURN_IF_NULL (A, EMPTY) ;
1189 Long m = A->nrow ;
1190 I = cholmod_l_speye (m, m, xtype, cc) ;
1191 Long rank = (I == NULL) ? EMPTY : SuiteSparseQR <Entry> (ordering, tol,
1192 econ, 1, A, I, NULL, Q, NULL, R, E, NULL, NULL, NULL, cc) ;
1193 cholmod_l_free_sparse (&I, cc) ;
1194 return (rank) ;
1195 }
1196
1197 template Long SuiteSparseQR <Complex> // returns rank(A) estimate
1198 (
1199 int ordering, // all, except 3:given treated as 0:fixed
1200 double tol, // columns with 2-norm <= tol are treated as zero
1201 Long econ, // e = max(min(m,econ),rank(A))
1202 cholmod_sparse *A, // m-by-n sparse matrix
1203 // outputs
1204 cholmod_sparse **Q, // m-by-e sparse matrix
1205 cholmod_sparse **R, // e-by-n sparse matrix
1206 Long **E, // permutation of 0:n-1
1207 cholmod_common *cc // workspace and parameters
1208 ) ;
1209
1210 template Long SuiteSparseQR <double> // returns rank(A) estimate
1211 (
1212 int ordering, // all, except 3:given treated as 0:fixed
1213 double tol, // columns with 2-norm <= tol are treated as zero
1214 Long econ, // e = max(min(m,econ),rank(A))
1215 cholmod_sparse *A, // m-by-n sparse matrix
1216 // outputs
1217 cholmod_sparse **Q, // m-by-e sparse matrix where e=max(econ,rank(A))
1218 cholmod_sparse **R, // e-by-n sparse matrix
1219 Long **E, // permutation of 0:n-1
1220 cholmod_common *cc // workspace and parameters
1221 ) ;
1222
1223 // -----------------------------------------------------------------------------
1224 // [Q,R,E] = qr(A), discarding Q
1225 // -----------------------------------------------------------------------------
1226
SuiteSparseQR(int ordering,double tol,Long econ,cholmod_sparse * A,cholmod_sparse ** R,Long ** E,cholmod_common * cc)1227 template <typename Entry> Long SuiteSparseQR // returns rank(A) estimate
1228 (
1229 int ordering, // all, except 3:given treated as 0:fixed
1230 double tol, // columns with 2-norm <= tol are treated as zero
1231 Long econ, // e = max(min(m,econ),rank(A))
1232 cholmod_sparse *A, // m-by-n sparse matrix
1233 // outputs
1234 cholmod_sparse **R, // e-by-n sparse matrix
1235 Long **E, // permutation of 0:n-1, NULL if identity
1236 cholmod_common *cc // workspace and parameters
1237 )
1238 {
1239 return (SuiteSparseQR <Entry> (ordering, tol, econ, 1, A,
1240 NULL, NULL, NULL, NULL, R, E, NULL, NULL, NULL, cc)) ;
1241 }
1242
1243 template Long SuiteSparseQR <Complex> // returns rank(A) estimate
1244 (
1245 int ordering, // all, except 3:given treated as 0:fixed
1246 double tol, // columns with 2-norm <= tol are treated as zero
1247 Long econ, // e = max(min(m,econ),rank(A))
1248 cholmod_sparse *A, // m-by-n sparse matrix
1249 // outputs
1250 cholmod_sparse **R, // e-by-n sparse matrix
1251 Long **E, // permutation of 0:n-1, NULL if identity
1252 cholmod_common *cc // workspace and parameters
1253 ) ;
1254
1255 template Long SuiteSparseQR <double> // returns rank(A) estimate
1256 (
1257 int ordering, // all, except 3:given treated as 0:fixed
1258 double tol, // columns with 2-norm <= tol are treated as zero
1259 Long econ, // e = max(min(m,econ),rank(A))
1260 cholmod_sparse *A, // m-by-n sparse matrix
1261 // outputs
1262 cholmod_sparse **R, // e-by-n sparse matrix
1263 Long **E, // permutation of 0:n-1, NULL if identity
1264 cholmod_common *cc // workspace and parameters
1265 ) ;
1266
1267 // -----------------------------------------------------------------------------
1268 // [C,R,E] = qr(A,B) where C and B are both dense
1269 // -----------------------------------------------------------------------------
1270
1271 // returns rank(A) estimate if successful, EMPTY otherwise
SuiteSparseQR(int ordering,double tol,Long econ,cholmod_sparse * A,cholmod_dense * B,cholmod_dense ** C,cholmod_sparse ** R,Long ** E,cholmod_common * cc)1272 template <typename Entry> Long SuiteSparseQR
1273 (
1274 // inputs, not modified
1275 int ordering, // all, except 3:given treated as 0:fixed
1276 double tol, // columns with 2-norm <= tol are treated as zero
1277 Long econ, // e = max(min(m,econ),rank(A))
1278 cholmod_sparse *A, // m-by-n sparse matrix
1279 cholmod_dense *B, // m-by-nrhs dense matrix
1280 // outputs
1281 cholmod_dense **C, // C = Q'*B, an e-by-nrhs dense matrix
1282 cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
1283 Long **E, // permutation of 0:n-1, NULL if identity
1284 cholmod_common *cc // workspace and parameters
1285 )
1286 {
1287 return (SuiteSparseQR <Entry> (ordering, tol, econ, 0, A, NULL, B,
1288 NULL, C, R, E, NULL, NULL, NULL, cc)) ;
1289 }
1290
1291 template Long SuiteSparseQR <double>
1292 (
1293 // inputs, not modified
1294 int ordering, // all, except 3:given treated as 0:fixed
1295 double tol, // columns with 2-norm <= tol are treated as zero
1296 Long econ, // e = max(min(m,econ),rank(A))
1297 cholmod_sparse *A, // m-by-n sparse matrix
1298 cholmod_dense *B, // m-by-nrhs dense matrix
1299 // outputs
1300 cholmod_dense **C, // C = Q'*B, an e-by-nrhs dense matrix
1301 cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
1302 Long **E, // permutation of 0:n-1, NULL if identity
1303 cholmod_common *cc // workspace and parameters
1304 ) ;
1305
1306
1307 template Long SuiteSparseQR <Complex>
1308 (
1309 // inputs, not modified
1310 int ordering, // all, except 3:given treated as 0:fixed
1311 double tol, // columns with 2-norm <= tol are treated as zero
1312 Long econ, // e = max(min(m,econ),rank(A))
1313 cholmod_sparse *A, // m-by-n sparse matrix
1314 cholmod_dense *B, // m-by-nrhs dense matrix
1315 // outputs
1316 cholmod_dense **C, // C = Q'*B, an e-by-nrhs dense matrix
1317 cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
1318 Long **E, // permutation of 0:n-1, NULL if identity
1319 cholmod_common *cc // workspace and parameters
1320 ) ;
1321
1322
1323 // -----------------------------------------------------------------------------
1324 // [C,R,E] = qr(A,B) where C and B are both sparse
1325 // -----------------------------------------------------------------------------
1326
1327 // returns rank(A) estimate if successful, EMPTY otherwise
SuiteSparseQR(int ordering,double tol,Long econ,cholmod_sparse * A,cholmod_sparse * B,cholmod_sparse ** C,cholmod_sparse ** R,Long ** E,cholmod_common * cc)1328 template <typename Entry> Long SuiteSparseQR
1329 (
1330 // inputs, not modified
1331 int ordering, // all, except 3:given treated as 0:fixed
1332 double tol, // columns with 2-norm <= tol are treated as zero
1333 Long econ, // e = max(min(m,econ),rank(A))
1334 cholmod_sparse *A, // m-by-n sparse matrix
1335 cholmod_sparse *B, // m-by-nrhs sparse matrix
1336 // outputs
1337 cholmod_sparse **C, // C = Q'*B, an e-by-nrhs sparse matrix
1338 cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
1339 Long **E, // permutation of 0:n-1, NULL if identity
1340 cholmod_common *cc // workspace and parameters
1341 )
1342 {
1343 return (SuiteSparseQR <Entry> (ordering, tol, econ, 0, A, B, NULL,
1344 C, NULL, R, E, NULL, NULL, NULL, cc)) ;
1345 }
1346
1347 template Long SuiteSparseQR <double>
1348 (
1349 // inputs, not modified
1350 int ordering, // all, except 3:given treated as 0:fixed
1351 double tol, // columns with 2-norm <= tol are treated as zero
1352 Long econ, // e = max(min(m,econ),rank(A))
1353 cholmod_sparse *A, // m-by-n sparse matrix
1354 cholmod_sparse *B, // m-by-nrhs sparse matrix
1355 // outputs
1356 cholmod_sparse **C, // C = Q'*B, an e-by-nrhs sparse matrix
1357 cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
1358 Long **E, // permutation of 0:n-1, NULL if identity
1359 cholmod_common *cc // workspace and parameters
1360 ) ;
1361
1362
1363 template Long SuiteSparseQR <Complex>
1364 (
1365 // inputs, not modified
1366 int ordering, // all, except 3:given treated as 0:fixed
1367 double tol, // columns with 2-norm <= tol are treated as zero
1368 Long econ, // e = max(min(m,econ),rank(A))
1369 cholmod_sparse *A, // m-by-n sparse matrix
1370 cholmod_sparse *B, // m-by-nrhs sparse matrix
1371 // outputs
1372 cholmod_sparse **C, // C = Q'*B, an e-by-nrhs sparse matrix
1373 cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
1374 Long **E, // permutation of 0:n-1, NULL if identity
1375 cholmod_common *cc // workspace and parameters
1376 ) ;
1377
1378 // -----------------------------------------------------------------------------
1379 // [Q,R,E] = qr(A) where Q is returned in Householder form
1380 // -----------------------------------------------------------------------------
1381
1382 // returns rank(A) estimate if successful, EMPTY otherwise
SuiteSparseQR(int ordering,double tol,Long econ,cholmod_sparse * A,cholmod_sparse ** R,Long ** E,cholmod_sparse ** H,Long ** HPinv,cholmod_dense ** HTau,cholmod_common * cc)1383 template <typename Entry> Long SuiteSparseQR
1384 (
1385 // inputs, not modified
1386 int ordering, // all, except 3:given treated as 0:fixed
1387 double tol, // columns with 2-norm <= tol are treated as zero
1388 Long econ, // e = max(min(m,econ),rank(A))
1389 cholmod_sparse *A, // m-by-n sparse matrix
1390 // outputs
1391 cholmod_sparse **R, // the R factor
1392 Long **E, // permutation of 0:n-1, NULL if identity
1393 cholmod_sparse **H, // the Householder vectors (m-by-nh)
1394 Long **HPinv, // size m; row permutation for H
1395 cholmod_dense **HTau, // size 1-by-nh, Householder coefficients
1396 cholmod_common *cc // workspace and parameters
1397 )
1398 {
1399 return (SuiteSparseQR <Entry> (ordering, tol, econ, EMPTY, A,
1400 NULL, NULL, NULL, NULL, R, E, H, HPinv, HTau, cc)) ;
1401 }
1402
1403 template Long SuiteSparseQR <double>
1404 (
1405 // inputs, not modified
1406 int ordering, // all, except 3:given treated as 0:fixed
1407 double tol, // columns with 2-norm <= tol are treated as zero
1408 Long econ, // e = max(min(m,econ),rank(A))
1409 cholmod_sparse *A, // m-by-n sparse matrix
1410 // outputs
1411 cholmod_sparse **R, // the R factor
1412 Long **E, // permutation of 0:n-1, NULL if identity
1413 cholmod_sparse **H, // the Householder vectors (m-by-nh)
1414 Long **HPinv, // size m; row permutation for H
1415 cholmod_dense **HTau, // size 1-by-nh, Householder coefficients
1416 cholmod_common *cc // workspace and parameters
1417 ) ;
1418
1419 template Long SuiteSparseQR <Complex>
1420 (
1421 // inputs, not modified
1422 int ordering, // all, except 3:given treated as 0:fixed
1423 double tol, // columns with 2-norm <= tol are treated as zero
1424 Long econ, // e = max(min(m,econ),rank(A))
1425 cholmod_sparse *A, // m-by-n sparse matrix
1426 // outputs
1427 cholmod_sparse **R, // the R factor
1428 Long **E, // permutation of 0:n-1, NULL if identity
1429 cholmod_sparse **H, // the Householder vectors (m-by-nh)
1430 Long **HPinv, // size m; row permutation for H
1431 cholmod_dense **HTau, // size 1-by-nh, Householder coefficients
1432 cholmod_common *cc // workspace and parameters
1433 ) ;
1434