1 // =============================================================================
2 // === spqr_rhpack =============================================================
3 // =============================================================================
4
5 /* Copy the R matrix and optionally the H matrix in the frontal matrix F and
6 store them in packed form. The columns of R and H are interleaved. F can
7 be overwritten with R+H (the pack can occur in-place), but in that case,
8 the C matrix is destroyed. If H is not kept and the pack is done in-place,
9 then H is destroyed. H is not stored in packed form if it is not kept.
10
11 If dead columns appear in F, the leading part of R is a squeezed upper
12 triangular matrix. In the example below, m = 9, n = 6, npiv = 4,
13 and the 3rd column is dead (Stair [2] == 0).
14
15 0: r r r r r r <- Stair [2] = 0, denotes dead column
16 1: h r r r r r
17 2: h h . r r r
18 3: h h . h c c
19 4: - h . h h c <- Stair [0] = 4
20 5: . h . h h h
21 6: . - . h h h <- Stair [1] = 6
22 7: . . . - - h <- Stair [3] = Stair [4] = 6
23 8: . . . . . h
24 - <- Stair [5] = 9
25
26 The number of rows of R is equal to the number of leading npiv columns that
27 are not dead; this number must be <= m. In this example, rm = 3.
28
29 The staircase defines the number of entries in each column, or is equal
30 to zero to denote a dead column.
31 */
32
33 #include "spqr.hpp"
34
spqr_rhpack(int keepH,Long m,Long n,Long npiv,Long * Stair,Entry * F,Entry * R,Long * p_rm)35 template <typename Entry> Long spqr_rhpack // returns # of entries in R+H
36 (
37 // input, not modified
38 int keepH, // if true, then H is packed
39 Long m, // # of rows in F
40 Long n, // # of columns in F
41 Long npiv, // number of pivotal columns in F
42 Long *Stair, // size npiv; column j is dead if Stair [j] == 0.
43 // Only the first npiv columns can be dead.
44
45 // input, not modified (unless the pack occurs in-place)
46 Entry *F, // m-by-n frontal matrix in column-major order
47
48 // output, contents not defined on input
49 Entry *R, // packed columns of R+H
50 Long *p_rm // number of rows in R block
51 )
52 {
53 Entry *R0 = R ;
54 Long i, k, h, t, rm ;
55
56 // -------------------------------------------------------------------------
57 // get inputs
58 // -------------------------------------------------------------------------
59
60 ASSERT (m >= 0 && n >= 0 && npiv <= n && npiv >= 0) ;
61
62 if (m <= 0 || n <= 0)
63 {
64 *p_rm = 0 ; // no rows in R block
65 return (0) ; // nothing to do
66 }
67
68 ASSERT (R != NULL && F != NULL) ;
69 ASSERT (R <= F // can be packed in-place, in F
70 || R >= F + m*n) ; // or must appear after F
71
72 // -------------------------------------------------------------------------
73 // pack the squeezed upper triangular part of R
74 // -------------------------------------------------------------------------
75
76 rm = 0 ; // number of rows in R block
77 for (k = 0 ; k < npiv ; k++)
78 {
79 // get the staircase
80 t = Stair [k] ; // F (0:t-1,k) contains R and H
81 ASSERT (t >= 0 && t <= m) ;
82 if (t == 0)
83 {
84 t = rm ; // dead col, R (0:rm-1,k) only, no H
85 }
86 else if (rm < m)
87 {
88 rm++ ; // column k is not dead
89 }
90 if (keepH)
91 {
92 // pack R (0:rm-1,k) and H (rm:t-1,k)
93 for (i = 0 ; i < t ; i++)
94 {
95 *(R++) = F [i] ;
96 }
97 }
98 else
99 {
100 // pack R (0:rm-1,k), discarding H
101 for (i = 0 ; i < rm ; i++)
102 {
103 *(R++) = F [i] ;
104 }
105 }
106 F += m ; // advance to the next column of F
107 }
108
109 // -------------------------------------------------------------------------
110 // pack the rectangular part of R
111 // -------------------------------------------------------------------------
112
113 h = rm ; // the column of H starts in row h
114 for ( ; k < n ; k++)
115 {
116
117 // pack R (0:rm-1,k)
118 for (i = 0 ; i < rm ; i++)
119 {
120 *(R++) = F [i] ;
121 }
122
123 if (keepH)
124 {
125 // pack H (h:t-1,k)
126 t = Stair [k] ; // get the staircase
127 ASSERT (t >= rm && t <= m) ;
128 h = MIN (h+1, m) ; // one more row of C to skip over
129 for (i = h ; i < t ; i++)
130 {
131 *(R++) = F [i] ;
132 }
133 }
134
135 F += m ; // advance to the next column of F
136 }
137
138 *p_rm = rm ; // return # of rows in R block
139 return (R-R0) ; // return # of packed entries
140 }
141
142
143 // =============================================================================
144
145 template Long spqr_rhpack <double> // returns # of entries in R+H
146 (
147 // input, not modified
148 int keepH, // if true, then H is packed
149 Long m, // # of rows in F
150 Long n, // # of columns in F
151 Long npiv, // number of pivotal columns in F
152 Long *Stair, // size npiv; column j is dead if Stair [j] == 0.
153 // Only the first npiv columns can be dead.
154
155 // input, not modified (unless the pack occurs in-place)
156 double *F, // m-by-n frontal matrix in column-major order
157
158 // output, contents not defined on input
159 double *R, // packed columns of R+H
160 Long *p_rm // number of rows in R block
161 ) ;
162
163 // =============================================================================
164
165 template Long spqr_rhpack <Complex> // returns # of entries in R+H
166 (
167 // input, not modified
168 int keepH, // if true, then H is packed
169 Long m, // # of rows in F
170 Long n, // # of columns in F
171 Long npiv, // number of pivotal columns in F
172 Long *Stair, // size npiv; column j is dead if Stair [j] == 0.
173 // Only the first npiv columns can be dead.
174
175 // input, not modified (unless the pack occurs in-place)
176 Complex *F, // m-by-n frontal matrix in column-major order
177
178 // output, contents not defined on input
179 Complex *R, // packed columns of R+H
180 Long *p_rm // number of rows in R block
181 ) ;
182