1 // =============================================================================
2 // === spqr_1factor ============================================================
3 // =============================================================================
4
5 /* Compute the QR factorization of a sparse matrix, both symbolic and numeric.
6
7 This function exploits column singletons, and thus cannot be
8 split into symbolic and numeric phases.
9
10 Outline:
11
12 function [C,R,E,X,H] = sparseQR (A,B)
13 E = fill reducing ordering of A, with singletons coming first
14 S = A*E
15 Y = [S B] or [S2 B2]
16 QR factorization of Y, obtaining factor R, Householder vectors H,
17 and applying H to B2 to get C2
18 X = E*(R\C)
19
20 ordering options:
21 0 or 3: fixed
22 1: natural (only look for singletons)
23 2: colamd after finding singletons
24 4: CHOLMOD fter finding singletons
25 5: amd(A'*A) after finding singletons
26 6: metis(A'*A) after finding singletons
27 7: SuiteSparseQR default (selects COLAMD, AMD, or METIS)
28
29 First, the column singletons of the sparse matrix A are found. If the
30 column singletons and their corresponding rows are permuted to the left and
31 top, respectively, of a matrix, and the matrix is not rank deficient, the
32 result looks like:
33
34 c x x x x x x x
35 . c x x x x x x
36 . . c x x x x x
37 . . . c x x x x
38 . . . . a a a a
39 . . . . a a a a
40 . . . . a a a a
41 . . . . a a a a
42 . . . . a a a a
43
44 where c denotes a column singleton j and its corresponding row i, x denotes
45 an entry in a singleton row, and a denotes the rest of the matrix. No
46 column in the S2 submatrix (the "a" entries) has just one entry in it
47 unless that entry falls below tol.
48
49 A singleton entry c is not accepted if its magnitude falls below tol. If
50 tol is negative, any entry present in the data structure of A is
51 acceptable, even an explicitly-stored zero entry).
52
53 If the matrix is structurally rank-deficient, then the permuted matrix can
54 look like the following, where in this case the 3rd column is "dead":
55
56 c x x x x x x x
57 . c x x x x x x
58 . . . c x x x x
59 . . . . a a a a
60 . . . . a a a a
61 . . . . a a a a
62 . . . . a a a a
63 . . . . a a a a
64 . . . . a a a a
65
66 After the column singletons are found, the remaining S2 matrix (and R12) is
67 ordered with a fill-reducing ordering, resulting in the ordering E. The
68 columns of B are not permuted.
69
70 The matrix [A*E B] is split into a block 2-by-3 form:
71
72 R11 R12 B1
73 0 S2 B2
74
75 where R1 = [R11 R12] are the singleton rows of A, and [R11 ; 0] are the
76 singleton columns. R1 is stored in row-oriented form (with sorted rows)
77 and Y = [S2 B2] is stored in column-oriented form (with sorted columns if
78 the input matrix A has sorted columns). S2 contains no empty columns, and
79 any columns in S2 with just one entry have a 2-norm less than tol. Note
80 that the rows of Y are permuted according to the singleton rows.
81 */
82
83 #include "spqr.hpp"
84
spqr_1factor(int ordering,double tol,Long bncols,int keepH,cholmod_sparse * A,Long ldb,Long * Bp,Long * Bi,Entry * Bx,cholmod_common * cc)85 template <typename Entry> SuiteSparseQR_factorization <Entry> *spqr_1factor
86 (
87 // inputs, not modified
88 int ordering, // all, except 3:given treated as 0:fixed
89 double tol, // only accept singletons above tol. If tol <= -2,
90 // then use the default tolerance
91 Long bncols, // number of columns of B
92 int keepH, // if TRUE, keep the Householder vectors
93 cholmod_sparse *A, // m-by-n sparse matrix
94 Long ldb, // if dense, the leading dimension of B
95 Long *Bp, // size bncols+1, column pointers of B
96 Long *Bi, // size bnz = Bp [bncols], row indices of B
97 Entry *Bx, // size bnz, numerical values of B
98
99 // workspace and parameters
100 cholmod_common *cc
101 )
102 {
103 spqr_symbolic *QRsym ;
104 spqr_numeric <Entry> *QRnum ;
105 SuiteSparseQR_factorization <Entry> *QR ;
106 Long *Yp, *Yi, *Q1fill, *R1p, *R1j, *P1inv, *Ap, *Ai ;
107 Entry *Yx, *R1x, *Ax ;
108 Long noY, anz, a2nz, r1nz, ynz, i, j, k, p, p2, bnz, py, n1rows,
109 n1cols, n2, Bsparse, d, iold, inew, m, n ;
110 cholmod_sparse *Y = NULL ;
111
112 double t0 = SuiteSparse_time ( ) ;
113 double t1, t2 ;
114
115 // -------------------------------------------------------------------------
116 // get inputs and allocate result
117 // -------------------------------------------------------------------------
118
119 m = A->nrow ;
120 n = A->ncol ;
121 Ap = (Long *) A->p ;
122 Ai = (Long *) A->i ;
123 Ax = (Entry *) A->x ;
124
125 QR = (SuiteSparseQR_factorization <Entry> *)
126 cholmod_l_malloc (1, sizeof (SuiteSparseQR_factorization <Entry>), cc) ;
127
128 if (cc->status < CHOLMOD_OK)
129 {
130 // out of memory
131 return (NULL) ;
132 }
133
134 QR->QRsym = NULL ;
135 QR->QRnum = NULL ;
136
137 QR->R1p = NULL ;
138 QR->R1j = NULL ;
139 QR->R1x = NULL ;
140 QR->P1inv = NULL ;
141 QR->Q1fill = NULL ;
142 QR->Rmap = NULL ;
143 QR->RmapInv = NULL ;
144 QR->HP1inv = NULL ;
145
146 QR->narows = m ;
147 QR->nacols = n ;
148 QR->n1rows = 0 ;
149 QR->n1cols = 0 ;
150
151 QR->r1nz = 0 ;
152 r1nz = 0 ;
153
154 // B is an optional input. It can be sparse or dense
155 Bsparse = (Bp != NULL && Bi != NULL) ;
156 if (Bx == NULL)
157 {
158 // B is not present; force bncols to be zero
159 bncols = 0 ;
160 }
161
162 QR->bncols = bncols ;
163
164 // -------------------------------------------------------------------------
165 // find the default tol, if requested
166 // -------------------------------------------------------------------------
167
168 if (tol <= SPQR_DEFAULT_TOL)
169 {
170 tol = spqr_tol <Entry> (A, cc) ;
171 }
172 if (tol < 0)
173 {
174 // no rank detection will be performed
175 QR->allow_tol = FALSE ;
176 tol = EMPTY ;
177 }
178 else
179 {
180 QR->allow_tol = TRUE ;
181 }
182 QR->tol = tol ;
183
184 // -------------------------------------------------------------------------
185 // find singletons and construct column pointers for the A part of Y
186 // -------------------------------------------------------------------------
187
188 // These return R1p, P1inv, and Y; but they are all NULL if out of memory.
189 // Note that only Y->p is allocated (Y->i and Y->x are dummy placeholders
190 // of one Long and one Entry, each, actually). The entries of Y are
191 // allocated later, below.
192
193 if (ordering == SPQR_ORDERING_GIVEN)
194 {
195 ordering = SPQR_ORDERING_FIXED ;
196 }
197
198 if (ordering == SPQR_ORDERING_FIXED)
199 {
200 // fixed ordering: find column singletons without permuting columns
201 Q1fill = NULL ;
202 spqr_1fixed <Entry> (tol, bncols, A,
203 &R1p, &P1inv, &Y, &n1cols, &n1rows, cc) ;
204 }
205 else
206 {
207 // natural or fill-reducing ordering: find column singletons with
208 // column permutations allowed, then permute the pruned submatrix with
209 // a fill-reducing ordering if ordering is not SPQR_ORDERING_NATURAL.
210 spqr_1colamd <Entry> (ordering, tol, bncols, A, &Q1fill,
211 &R1p, &P1inv, &Y, &n1cols, &n1rows, cc) ;
212 ordering = cc->SPQR_istat [7] ;
213 }
214
215 if (cc->status < CHOLMOD_OK)
216 {
217 // out of memory
218 spqr_freefac (&QR, cc) ;
219 return (NULL) ;
220 }
221
222 QR->R1p = R1p ;
223 QR->P1inv = P1inv ;
224 QR->Q1fill = Q1fill ;
225 QR->n1rows = n1rows ;
226 QR->n1cols = n1cols ;
227
228 noY = (Y == NULL) ; // A will be factorized, not Y
229 ASSERT (noY == (n1cols == 0 && bncols == 0)) ;
230 Yp = noY ? NULL : (Long *) Y->p ;
231 anz = Ap [n] ; // nonzeros in A
232 a2nz = noY ? anz : Yp [n-n1cols] ; // nonzeros in S2
233 n2 = n - n1cols ; // number of columns of S2
234
235 // Y is NULL, or of size (m-n1rows)-by-(n-n1cols+bncols)
236 ASSERT (IMPLIES (Y != NULL, ((Long) Y->nrow == m-n1rows))) ;
237 ASSERT (IMPLIES (Y != NULL, ((Long) Y->ncol == n-n1cols+bncols))) ;
238
239 // Y, if allocated, has no space for any entries yet
240 ynz = 0 ;
241
242 // -------------------------------------------------------------------------
243 // construct the column pointers for the B or B2 part of Y
244 // -------------------------------------------------------------------------
245
246 if (noY)
247 {
248
249 // A will be factorized instead of Y. There is no B. C or X can exist
250 // as empty matrices with rows but no columns.
251 // There are no singletons.
252 ASSERT (Yp == NULL) ;
253 ASSERT (R1p == NULL) ;
254 ASSERT (P1inv == NULL) ;
255 ASSERT (n1rows == 0) ;
256 ASSERT (n1cols == 0) ;
257 ASSERT (a2nz == Ap [n]) ;
258 ASSERT (bncols == 0) ;
259
260 }
261 else if (n1cols == 0)
262 {
263
264 // ---------------------------------------------------------------------
265 // construct the column pointers for the B part of Y = [S B]
266 // ---------------------------------------------------------------------
267
268 ASSERT (R1p == NULL) ;
269 ASSERT (P1inv == NULL) ;
270 ASSERT (n1rows == 0) ;
271 ASSERT (a2nz == Ap [n]) ;
272
273 ynz = a2nz ;
274 if (Bsparse)
275 {
276 // B is sparse
277 for (k = 0 ; k < bncols ; k++)
278 {
279 Yp [(n-n1cols)+k] = ynz ;
280 d = Bp [k+1] - Bp [k] ;
281 ynz += d ;
282 }
283 }
284 else
285 {
286 // B is dense
287 Entry *B1 = Bx ;
288 for (k = 0 ; k < bncols ; k++)
289 {
290 // count the nonzero entries in column k of B
291 Yp [(n-n1cols)+k] = ynz ;
292 d = 0 ;
293 for (i = 0 ; i < m ; i++)
294 {
295 if (B1 [i] != (Entry) 0)
296 {
297 d++ ;
298 }
299 }
300 B1 += ldb ;
301 ynz += d ;
302 }
303 }
304 Yp [(n-n1cols)+bncols] = ynz ;
305
306 }
307 else
308 {
309
310 // ---------------------------------------------------------------------
311 // construct the column pointers for the B2 part of Y = [S2 B2]
312 // ---------------------------------------------------------------------
313
314 ynz = a2nz ;
315 if (Bsparse)
316 {
317 // B is sparse
318 for (k = 0 ; k < bncols ; k++)
319 {
320 // count the nonzero entries in column k of B2
321 Yp [(n-n1cols)+k] = ynz ;
322 d = 0 ;
323 for (p = Bp [k] ; p < Bp [k+1] ; p++)
324 {
325 iold = Bi [p] ;
326 inew = P1inv [iold] ;
327 if (inew >= n1rows)
328 {
329 d++ ;
330 }
331 }
332 ynz += d ;
333 }
334 }
335 else
336 {
337 // B is dense
338 Entry *B1 = Bx ;
339 for (k = 0 ; k < bncols ; k++)
340 {
341 // count the nonzero entries in column k of B2
342 Yp [(n-n1cols)+k] = ynz ;
343 d = 0 ;
344 for (iold = 0 ; iold < m ; iold++)
345 {
346 inew = P1inv [iold] ;
347 if (inew >= n1rows && B1 [iold] != (Entry) 0)
348 {
349 d++ ;
350 }
351 }
352 B1 += ldb ;
353 ynz += d ;
354 }
355 }
356 Yp [(n-n1cols)+bncols] = ynz ;
357 }
358
359
360 // -------------------------------------------------------------------------
361 // allocate the nonzeros for Y
362 // -------------------------------------------------------------------------
363
364 if (noY)
365 {
366 // no singletons found, and B is empty. pass Y=A to QR factorization,
367 // and pass in Q1fill as the "user-provided" ordering
368 ASSERT (Yp == NULL) ;
369 Yi = NULL ;
370 Yx = NULL ;
371 }
372 else
373 {
374 cholmod_l_reallocate_sparse (ynz, Y, cc) ;
375 Yi = (Long *) Y->i ;
376 Yx = (Entry *) Y->x ;
377 }
378
379 if (cc->status < CHOLMOD_OK)
380 {
381 // out of memory
382 spqr_freefac (&QR, cc) ;
383 cholmod_l_free_sparse (&Y, cc) ;
384 return (NULL) ;
385 }
386
387 // -------------------------------------------------------------------------
388 // create the pattern and values of Y and R1
389 // -------------------------------------------------------------------------
390
391 if (noY)
392 {
393
394 // ---------------------------------------------------------------------
395 // R1 does not exist
396 // ---------------------------------------------------------------------
397
398 ASSERT (R1p == NULL) ;
399 R1j = NULL ;
400 R1x = NULL ;
401
402 }
403 else if (n1cols == 0)
404 {
405
406 // ---------------------------------------------------------------------
407 // R1 does not exist
408 // ---------------------------------------------------------------------
409
410 ASSERT (R1p == NULL) ;
411 R1j = NULL ;
412 R1x = NULL ;
413
414 // ---------------------------------------------------------------------
415 // construct the A part of Y = [S B]
416 // ---------------------------------------------------------------------
417
418 ASSERT (anz == a2nz) ;
419 py = 0 ;
420 for (k = 0 ; k < n ; k++)
421 {
422 j = Q1fill ? Q1fill [k] : k ;
423 ASSERT (py == Yp [k]) ;
424 for (p = Ap [j] ; p < Ap [j+1] ; p++)
425 {
426 Yi [py] = Ai [p] ;
427 Yx [py] = Ax [p] ;
428 py++ ;
429 }
430 }
431 ASSERT (py == anz) ;
432 ASSERT (py == Yp [n]) ;
433
434 // ---------------------------------------------------------------------
435 // construct the B part of Y = [S B]
436 // ---------------------------------------------------------------------
437
438 if (Bsparse)
439 {
440 // B is sparse
441 bnz = Bp [bncols] ;
442 for (p = 0 ; p < bnz ; p++)
443 {
444 Yi [py++] = Bi [p] ;
445 }
446 py = anz ;
447 for (p = 0 ; p < bnz ; p++)
448 {
449 Yx [py++] = Bx [p] ;
450 }
451 }
452 else
453 {
454 // B is dense
455 Entry *B1 = Bx ;
456 for (k = 0 ; k < bncols ; k++)
457 {
458 ASSERT (py == Yp [n+k]) ;
459 for (i = 0 ; i < m ; i++)
460 {
461 Entry bij = B1 [i] ;
462 if (bij != (Entry) 0)
463 {
464 Yi [py] = i ;
465 Yx [py] = bij ;
466 py++ ;
467 }
468 }
469 B1 += ldb ;
470 }
471 }
472 ASSERT (py == ynz) ;
473
474 }
475 else
476 {
477
478 // ---------------------------------------------------------------------
479 // R1p = cumsum ([0 R1p])
480 // ---------------------------------------------------------------------
481
482 r1nz = spqr_cumsum (n1rows, R1p) ; // Long overflow cannot occur
483 PR (("total nonzeros in R1: %ld\n", r1nz)) ;
484
485 // ---------------------------------------------------------------------
486 // allocate R1
487 // ---------------------------------------------------------------------
488
489 R1j = (Long *) cholmod_l_malloc (r1nz, sizeof (Long ), cc) ;
490 R1x = (Entry *) cholmod_l_malloc (r1nz, sizeof (Entry), cc) ;
491 QR->R1j = R1j ;
492 QR->R1x = R1x ;
493 QR->r1nz = r1nz ;
494
495 if (cc->status < CHOLMOD_OK)
496 {
497 // out of memory
498 spqr_freefac (&QR, cc) ;
499 cholmod_l_free_sparse (&Y, cc) ;
500 return (NULL) ;
501 }
502
503 // ---------------------------------------------------------------------
504 // scan A and construct R11
505 // ---------------------------------------------------------------------
506
507 // At this point, R1p [i] points to the start of row i:
508 // for (Long t = 0 ; t <= n1rows ; t++) Rsave [t] = R1p [t] ;
509
510 for (k = 0 ; k < n1cols ; k++)
511 {
512 j = Q1fill ? Q1fill [k] : k ;
513 for (p = Ap [j] ; p < Ap [j+1] ; p++)
514 {
515 // row i of A is row inew after singleton permutation
516 i = Ai [p] ;
517 inew = P1inv [i] ;
518 ASSERT (inew < n1rows) ;
519 // A (i,j) is in a singleton row. It becomes R1 (inew,k)
520 p2 = R1p [inew]++ ;
521 ASSERT (p2 < R1p [inew+1]) ;
522 R1j [p2] = k ;
523 R1x [p2] = Ax [p] ;
524 }
525 }
526
527 // ---------------------------------------------------------------------
528 // scan A and construct R12 and the S2 part of Y = [S2 B2]
529 // ---------------------------------------------------------------------
530
531 py = 0 ;
532 for ( ; k < n ; k++)
533 {
534 j = Q1fill ? Q1fill [k] : k ;
535 ASSERT (py == Yp [k-n1cols]) ;
536 for (p = Ap [j] ; p < Ap [j+1] ; p++)
537 {
538 // row i of A is row inew after singleton permutation
539 i = Ai [p] ;
540 inew = P1inv [i] ;
541 if (inew < n1rows)
542 {
543 // A (i,j) is in a singleton row. It becomes R1 (inew,k)
544 p2 = R1p [inew]++ ;
545 ASSERT (p2 < R1p [inew+1]) ;
546 R1j [p2] = k ;
547 R1x [p2] = Ax [p] ;
548 }
549 else
550 {
551 // A (i,j) is not in a singleton row. Place it in
552 // Y (inew-n1rows, k-n1cols)
553 Yi [py] = inew - n1rows ;
554 Yx [py] = Ax [p] ;
555 py++ ;
556 }
557 }
558 }
559 ASSERT (py == Yp [n-n1cols]) ;
560
561 // ---------------------------------------------------------------------
562 // restore the row pointers for R1
563 // ---------------------------------------------------------------------
564
565 spqr_shift (n1rows, R1p) ;
566
567 // the row pointers are back to what they were:
568 // for (Long t = 0 ; t <= n1rows ; t++) ASSERT (Rsave [t] == R1p [t]) ;
569
570 // ---------------------------------------------------------------------
571 // construct the B2 part of Y = [S2 B2]
572 // ---------------------------------------------------------------------
573
574 if (Bsparse)
575 {
576 // B is sparse
577 for (k = 0 ; k < bncols ; k++)
578 {
579 // construct the nonzero entries in column k of B2
580 ASSERT (py == Yp [k+(n-n1cols)]) ;
581 for (p = Bp [k] ; p < Bp [k+1] ; p++)
582 {
583 iold = Bi [p] ;
584 inew = P1inv [iold] ;
585 if (inew >= n1rows)
586 {
587 Yi [py] = inew - n1rows ;
588 Yx [py] = Bx [p] ;
589 py++ ;
590 }
591 }
592 }
593 }
594 else
595 {
596 // B is dense
597 Entry *B1 = Bx ;
598 for (k = 0 ; k < bncols ; k++)
599 {
600 // construct the nonzero entries in column k of B2
601 ASSERT (py == Yp [k+(n-n1cols)]) ;
602 for (iold = 0 ; iold < m ; iold++)
603 {
604 inew = P1inv [iold] ;
605 if (inew >= n1rows)
606 {
607 Entry bij = B1 [iold] ;
608 if (bij != (Entry) 0)
609 {
610 Yi [py] = inew - n1rows ;
611 Yx [py] = bij ;
612 py++ ;
613 }
614 }
615 }
616 B1 += ldb ;
617 }
618 }
619 ASSERT (py == ynz) ;
620 }
621
622 // -------------------------------------------------------------------------
623 // QR factorization of A or Y
624 // -------------------------------------------------------------------------
625
626 if (noY)
627 {
628 // factorize A, with fill-reducing ordering already given in Q1fill
629 QRsym = spqr_analyze (A, SPQR_ORDERING_GIVEN, Q1fill,
630 tol >= 0, keepH, cc) ;
631 t1 = SuiteSparse_time ( ) ;
632 QRnum = spqr_factorize <Entry> (&A, FALSE, tol, n, QRsym, cc) ;
633 }
634 else
635 {
636 // fill-reducing ordering is already applied to Y; free Y when loaded
637 QRsym = spqr_analyze (Y, SPQR_ORDERING_FIXED, NULL,
638 tol >= 0, keepH, cc) ;
639 t1 = SuiteSparse_time ( ) ;
640 QRnum = spqr_factorize <Entry> (&Y, TRUE, tol, n2, QRsym, cc) ;
641 // Y has been freed
642 ASSERT (Y == NULL) ;
643 }
644
645 // record the actual ordering used (this will have been changed to GIVEN
646 // or FIXED, in spqr_analyze, but change it back to the ordering used by
647 // spqr_1fixed or spqr_1colamd.
648 cc->SPQR_istat [7] = ordering ;
649
650 QR->QRsym = QRsym ;
651 QR->QRnum = QRnum ;
652
653 if (cc->status < CHOLMOD_OK)
654 {
655 // out of memory
656 spqr_freefac (&QR, cc) ;
657 return (NULL) ;
658 }
659
660 // singletons do not take part in the symbolic analysis, and any columns
661 // of B are tacked on and do take part.
662 ASSERT (QRsym->n == n - n1cols + bncols) ;
663
664 cc->SPQR_istat [0] += r1nz ; // nnz (R)
665
666 // rank estimate of A, including singletons but excluding the columns of
667 // of B, in case [A B] was factorized.
668 QR->rank = n1rows + QRnum->rank1 ;
669 PR (("rank estimate of A: QR->rank = %ld = %ld + %ld\n",
670 QR->rank, n1rows, QRnum->rank1)) ;
671
672 // -------------------------------------------------------------------------
673 // construct global row permutation if H is kept and singletons exist
674 // -------------------------------------------------------------------------
675
676 // If there are no singletons, then HP1inv [0:m-1] and HPinv [0:m-1] would
677 // be identical, so HP1inv is not needed.
678
679 ASSERT ((n1cols == 0) == (P1inv == NULL)) ;
680 ASSERT (IMPLIES (n1cols == 0, n1rows == 0)) ;
681 ASSERT (n1cols >= n1rows) ;
682
683 if (keepH && n1cols > 0)
684 {
685 // construct the global row permutation. Currently, the row indices
686 // in H reflect the global R. P1inv is the singleton permutation,
687 // where a row index of Y = (P1inv (row of A) - n1rows), and
688 // row of R2 = QRnum->HPinv (row of Y). Combine these two into
689 // HP1inv, where a global row of R = HP1inv (a row of A)
690
691 Long kk ;
692 Long *HP1inv, *HPinv ;
693 QR->HP1inv = HP1inv = (Long *) cholmod_l_malloc (m, sizeof (Long), cc) ;
694 HPinv = QRnum->HPinv ;
695
696 if (cc->status < CHOLMOD_OK)
697 {
698 // out of memory
699 spqr_freefac (&QR, cc) ;
700 return (NULL) ;
701 }
702
703 for (i = 0 ; i < m ; i++)
704 {
705 // i is a row of A, k is a row index after row singletons are
706 // permuted. Then kk is a row index of the global R.
707 k = P1inv ? P1inv [i] : i ;
708 ASSERT (k >= 0 && k < m) ;
709 if (k < n1rows)
710 {
711 kk = k ;
712 }
713 else
714 {
715 // k-n1rows is a row index of Y, the matrix factorized by
716 // the QR factorization kernels (in QRsym and QRnum).
717 // HPinv [k-n1rows] gives a row index of R2, to which n1rows
718 // must be added to give a row of the global R.
719 kk = HPinv [k - n1rows] + n1rows ;
720 }
721 ASSERT (kk >= 0 && kk < m) ;
722 HP1inv [i] = kk ;
723 }
724 }
725
726 // -------------------------------------------------------------------------
727 // find the mapping for the squeezed R, if A is rank deficient
728 // -------------------------------------------------------------------------
729
730 if (QR->rank < n && !spqr_rmap <Entry> (QR, cc))
731 {
732 // out of memory
733 spqr_freefac (&QR, cc) ;
734 return (NULL) ;
735 }
736
737 // -------------------------------------------------------------------------
738 // output statistics
739 // -------------------------------------------------------------------------
740
741 cc->SPQR_istat [4] = QR->rank ; // estimated rank of A
742 cc->SPQR_istat [5] = n1cols ; // number of columns singletons
743 cc->SPQR_istat [6] = n1rows ; // number of singleton rows
744 cc->SPQR_tol_used = tol ; // tol used
745
746 t2 = SuiteSparse_time ( ) ;
747 cc->SPQR_analyze_time = t1 - t0 ; // analyze time, including singletons
748 cc->SPQR_factorize_time = t2 - t1 ; // factorize time
749
750 return (QR) ;
751 }
752
753
754 // =============================================================================
755
756 template SuiteSparseQR_factorization <double> *spqr_1factor <double>
757 (
758 // inputs, not modified
759 int ordering, // all, except 3:given treated as 0:fixed
760 double tol, // only accept singletons above tol
761 Long bncols, // number of columns of B
762 int keepH, // if TRUE, keep the Householder vectors
763 cholmod_sparse *A, // m-by-n sparse matrix
764 Long ldb, // if dense, the leading dimension of B
765 Long *Bp, // size bncols+1, column pointers of B
766 Long *Bi, // size bnz = Bp [bncols], row indices of B
767 double *Bx, // size bnz, numerical values of B
768
769 // workspace and parameters
770 cholmod_common *cc
771 ) ;
772
773 // =============================================================================
774
775 template SuiteSparseQR_factorization <Complex> *spqr_1factor <Complex>
776 (
777 // inputs, not modified
778 int ordering, // all, except 3:given treated as 0:fixed
779 double tol, // only accept singletons above tol
780 Long bncols, // number of columns of B
781 int keepH, // if TRUE, keep the Householder vectors
782 cholmod_sparse *A, // m-by-n sparse matrix
783 Long ldb, // if dense, the leading dimension of B
784 Long *Bp, // size bncols+1, column pointers of B
785 Long *Bi, // size bnz = Bp [bncols], row indices of B
786 Complex *Bx, // size bnz, numerical values of B
787
788 // workspace and parameters
789 cholmod_common *cc
790 ) ;
791