1 // =============================================================================
2 // === spqr_rcount =============================================================
3 // =============================================================================
4
5 // Count the number of explicit nonzeros in each column of R. Exact zero
6 // entries are excluded.
7
8 #include "spqr.hpp"
9
spqr_rcount(spqr_symbolic * QRsym,spqr_numeric<Entry> * QRnum,Long n1rows,Long econ,Long n2,int getT,Long * Ra,Long * Rb,Long * H2p,Long * p_nh)10 template <typename Entry> void spqr_rcount
11 (
12 // inputs, not modified
13 spqr_symbolic *QRsym,
14 spqr_numeric <Entry> *QRnum,
15
16 Long n1rows, // added to each row index of Ra and Rb
17 Long econ, // only get entries in rows n1rows to econ-1
18 Long n2, // Ra = R (:,0:n2-1), Rb = R (:,n2:n-1)
19 int getT, // if true, count Rb' instead of Rb
20
21 // input/output
22 // FUTURE : make Ra, Rb, H2 cholmod_sparse
23 Long *Ra, // size n2; Ra [j] += nnz (R (:,j)) if j < n2
24 Long *Rb, // If getT is false: size n-n2 and
25 // Rb [j-n2] += nnz (R (:,j)) if j >= n2.
26 // If getT is true: size econ, and
27 // Rb [i] += nnz (R (i, n2:n-1))
28 Long *H2p, // size rjsize+1. Column pointers for H.
29 // Only computed if H was kept during factorization.
30 // Only H2p [0..nh] is used.
31 Long *p_nh // number of Householder vectors (nh <= rjsize)
32 )
33 {
34 Entry **Rblock, *R, *Tau, *HTau ;
35 Long *Rp, *Rj, *Super, *HStair, *Stair, *Hm ;
36 char *Rdead ;
37 Long nf, j, f, col1, fp, pr, fn, rm, k, i, t, fm, h, getRa, getRb, nh,
38 row1, keepH, getH, hnz ;
39
40 // -------------------------------------------------------------------------
41 // get the contents of the QRsym and QRnum objects
42 // -------------------------------------------------------------------------
43
44 keepH = QRnum->keepH ;
45
46 getRa = (Ra != NULL) ;
47 getRb = (Rb != NULL) ;
48 getH = (H2p != NULL && p_nh != NULL) && keepH ;
49 if (!(getRa || getRb || getH))
50 {
51 // nothing to do
52 return ;
53 }
54
55 nf = QRsym->nf ;
56 // n = QRsym->n ;
57 Rblock = QRnum->Rblock ;
58 Rp = QRsym->Rp ;
59 Rj = QRsym->Rj ;
60 Super = QRsym->Super ;
61 Rdead = QRnum->Rdead ;
62
63 HStair = QRnum->HStair ;
64 HTau = QRnum->HTau ;
65 Hm = QRnum->Hm ;
66 Stair = NULL ;
67 Tau = NULL ;
68 fm = 0 ;
69 h = 0 ;
70 t = 0 ;
71 nh = 0 ;
72 hnz = 0 ;
73
74 // -------------------------------------------------------------------------
75 // examine the packed block for each front F
76 // -------------------------------------------------------------------------
77
78 row1 = n1rows ;
79 for (f = 0 ; f < nf ; f ++)
80 {
81 R = Rblock [f] ;
82 col1 = Super [f] ; // first pivot column in front F
83 fp = Super [f+1] - col1 ; // number of pivots in front F
84 pr = Rp [f] ; // pointer to row indices for F
85 fn = Rp [f+1] - pr ; // # of columns in front F
86
87 if (keepH)
88 {
89 Stair = HStair + pr ; // staircase of front F
90 Tau = HTau + pr ; // Householder coeff. for front F
91 fm = Hm [f] ; // # of rows in front F
92 h = 0 ; // H vbector starts in row h
93 }
94
95 rm = 0 ; // number of rows in R block
96 for (k = 0 ; k < fn ; k++)
97 {
98 // -----------------------------------------------------------------
99 // get the column and its staircase
100 // -----------------------------------------------------------------
101
102 if (k < fp)
103 {
104 // a pivotal column of front F
105 j = col1 + k ;
106 ASSERT (Rj [pr + k] == j) ;
107 if (keepH)
108 {
109 t = Stair [k] ; // length of R+H vector
110 ASSERT (t >= 0 && t <= fm) ;
111 if (t == 0)
112 {
113 t = rm ; // dead col, R only, no H
114 }
115 else if (rm < fm)
116 {
117 rm++ ; // column k is not dead
118 }
119 h = rm ; // H vector starts in row h
120 }
121 else
122 {
123 if (!Rdead [j])
124 {
125 rm++ ; // column k is not dead
126 }
127 }
128 }
129 else
130 {
131 // a non-pivotal column of front F
132 j = Rj [pr + k] ;
133 ASSERT (j >= Super [f+1] && j < QRsym->n) ;
134 if (keepH)
135 {
136 t = Stair [k] ; // length of R+H vector
137 ASSERT (t >= rm && t <= fm) ;
138 h = MIN (h+1, fm) ; // one more row of C to skip
139 }
140 }
141
142 // -----------------------------------------------------------------
143 // count nnz (R (0:econ-1,j)) for this R block
144 // -----------------------------------------------------------------
145
146 for (i = 0 ; i < rm ; i++)
147 {
148 // R (i,j) is nonzero (i local)
149 Entry rij = *(R++) ;
150 if (rij != (Entry) 0)
151 {
152 if (j < n2)
153 {
154 if (getRa && row1 + i < econ)
155 {
156 Ra [j]++ ;
157 }
158 }
159 else
160 {
161 if (getRb && row1 + i < econ)
162 {
163 if (getT)
164 {
165 Rb [row1+i]++ ;
166 }
167 else
168 {
169 Rb [j-n2]++ ;
170 }
171 }
172 }
173 }
174 }
175
176 // -----------------------------------------------------------------
177 // count nnz (H (:,pr+k))
178 // -----------------------------------------------------------------
179
180 if (keepH && t >= h)
181 {
182 // the Householder reflection is not empty
183 PR (("count terms in H k %ld, t %ld h %ld\n", k, t, h)) ;
184 if (getH && Tau [k] != (Entry) 0)
185 {
186 H2p [nh++] = hnz++ ; // count the implicit identity
187 for (i = h ; i < t ; i++)
188 {
189 Entry hij = *(R++) ;
190 if (hij != (Entry) 0)
191 {
192 hnz++ ; // H (i,pr+k) is nonzero
193 }
194 }
195 }
196 else
197 {
198 R += (t-h) ; // skip over the column of H
199 }
200 }
201 }
202 ASSERT (IMPLIES (keepH, QRnum->Hr [f] == rm)) ;
203 row1 += rm ; // count the squeezed rows of R
204 }
205
206 // -------------------------------------------------------------------------
207 // finalize the H column pointers
208 // -------------------------------------------------------------------------
209
210 if (getH)
211 {
212 H2p [nh] = hnz ;
213 *p_nh = nh ;
214 }
215 }
216
217
218 // =============================================================================
219
220 template void spqr_rcount <double>
221 (
222 // inputs, not modified
223 spqr_symbolic *QRsym,
224 spqr_numeric <double> *QRnum,
225
226 Long n1rows, // added to each row index of Ra and Rb
227 Long econ, // only get entries in rows n1rows to econ-1
228 Long n2, // Ra = R (:,0:n2-1), Rb = R (:,n2:n-1)
229 int getT, // if true, count Rb' instead of Rb
230
231 // input/output
232 Long *Ra, // size n2; Ra [j] += nnz (R (:,j)) if j < n2
233 Long *Rb, // If getT is false: size n-n2 and
234 // Rb [j-n2] += nnz (R (:,j)) if j >= n2.
235 // If getT is true: size econ, and
236 // Rb [i] += nnz (R (i, n2:n-1))
237 Long *H2p, // size rjsize+1. Column pointers for H.
238 // Only computed if H was kept during factorization.
239 // Only H2p [0..nh] is used.
240 Long *p_nh // number of Householder vectors (nh <= rjsize)
241 ) ;
242
243 // =============================================================================
244
245 template void spqr_rcount <Complex>
246 (
247 // inputs, not modified
248 spqr_symbolic *QRsym,
249 spqr_numeric <Complex> *QRnum,
250
251 Long n1rows, // added to each row index of Ra and Rb
252 Long econ, // only get entries in rows n1rows to econ-1
253 Long n2, // Ra = R (:,0:n2-1), Rb = R (:,n2:n-1)
254 int getT, // if true, count Rb' instead of Rb
255
256 // input/output
257 Long *Ra, // size n2; Ra [j] += nnz (R (:,j)) if j < n2
258 Long *Rb, // If getT is false: size n-n2 and
259 // Rb [j-n2] += nnz (R (:,j)) if j >= n2.
260 // If getT is true: size econ, and
261 // Rb [i] += nnz (R (i, n2:n-1))
262 Long *H2p, // size rjsize+1. Column pointers for H.
263 // Only computed if H was kept during factorization.
264 // Only H2p [0..nh] is used.
265 Long *p_nh // number of Householder vectors (nh <= rjsize)
266 ) ;
267