1 // =============================================================================
2 // === spqr_rmap ===============================================================
3 // =============================================================================
4
5 // R is squeezed, find the mapping that permutes it to trapezoidal form
6
7 // Rmap is a permutation that converts R from squeezed to upper trapezoidal.
8 // If R is already in upper triangular form, then Rmap is NULL (an implicit
9 // identity; Rmap [0:n-1] = 0:n-1), and this function is not used.
10
11 // If R is rank deficient, Rmap [j] = k if column j of the squeezed R is the
12 // kth column in the trapezoidal R. If j is a live column then
13 // k = Rmap [j] < QR->rank; otherwise k = Rmap [j] > QR->rank. RmapInv is
14 // the inverse of Rmap.
15
16 // Example: Suppose R has the following format:
17 //
18 // 0 1 2 3 4 5 6
19 // X x x x x x x
20 // . X x x x x x
21 // . . . X x x x
22 // . . . . . X x
23 // . . . . . . X
24 //
25 // Then Rmap is [0 1 5 2 6 3 4] and RmapInv is [0 1 3 5 6 2 4]. The rank of R
26 // is 5, and thus columns 2 and 4 (with Rmap [2] = 5 and Rmap [4] = 6) are both
27 // dead.
28
29 #include "spqr.hpp"
30
spqr_rmap(SuiteSparseQR_factorization<Entry> * QR,cholmod_common * cc)31 template <typename Entry> int spqr_rmap
32 (
33 SuiteSparseQR_factorization <Entry> *QR,
34 cholmod_common *cc
35 )
36 {
37 Long n, j, i, p, n1rows, n1cols ;
38 Long *Rmap, *RmapInv, *R1p, *R1j ;
39
40 n = QR->nacols ;
41 Rmap = QR->Rmap ;
42 RmapInv = QR->RmapInv ;
43
44 if (Rmap == NULL)
45 {
46 ASSERT (RmapInv == NULL) ;
47 QR->Rmap = Rmap = (Long *) cholmod_l_malloc (n, sizeof(Long), cc);
48 QR->RmapInv = RmapInv = (Long *) cholmod_l_malloc (n, sizeof(Long), cc);
49 if (cc->status < CHOLMOD_OK)
50 {
51 // out of memory
52 return (FALSE) ;
53 }
54 }
55
56 for (j = 0 ; j < n ; j++)
57 {
58 Rmap [j] = EMPTY ;
59 }
60
61 R1p = QR->R1p ;
62 R1j = QR->R1j ;
63 n1rows = QR->n1rows ;
64 n1cols = QR->n1cols ;
65
66 // find the mapping for the singleton rows
67 for (i = 0 ; i < n1rows ; i++)
68 {
69 // The ith row of R is a singleton row; find its corresponding
70 // pivotal column.
71 p = R1p [i] ;
72 ASSERT (R1p [i] < R1p [i+1]) ;
73 j = R1j [p] ;
74 ASSERT (j >= 0 && j < n1cols) ;
75 Rmap [j] = i ;
76 }
77
78 // find the mapping for the pivotal rows of the multifrontal R
79 char *Rdead = QR->QRnum->Rdead ;
80 for (j = n1cols ; j < n ; j++)
81 {
82 if (!Rdead [j-n1cols])
83 {
84 Rmap [j] = i++ ;
85 }
86 }
87 ASSERT (i == QR->rank) ;
88
89 // finish the mapping with the dead columns of R, both in the singleton
90 // part of R and the multifrontal part of R
91 for (j = 0 ; j < n ; j++)
92 {
93 if (Rmap [j] == EMPTY)
94 {
95 Rmap [j] = i++ ;
96 }
97 PR (("Rmap [%ld] = %ld (%d), rank = %ld\n",
98 j, Rmap [j], Rmap [j] >= QR->rank, QR->rank)) ;
99 }
100 ASSERT (i == n) ;
101
102 // construct the inverse of Rmap
103 for (j = 0 ; j < n ; j++)
104 {
105 i = Rmap [j] ;
106 RmapInv [i] = j ;
107 }
108 return (TRUE) ;
109 }
110
111 template int spqr_rmap <double>
112 (
113 SuiteSparseQR_factorization <double> *QR,
114 cholmod_common *cc
115 ) ;
116
117 template int spqr_rmap <Complex>
118 (
119 SuiteSparseQR_factorization <Complex> *QR,
120 cholmod_common *cc
121 ) ;
122