1 // =============================================================================
2 // === spqr_trapezoidal ========================================================
3 // =============================================================================
4
5 // Permute the columns of a "squeezed" R, possibly rank-deficient, into
6 // upper trapezoidal form. On input, Qfill gives the column permutation of
7 // A that gave the factor R (Q*R = A(:,Qfill). On output, T is upper
8 // trapezoidal, and Qtrap is the adjusted column ordering so that
9 // Q*T = A(:,Qtrap). T is permuted so that T = [T1 T2], where T1 is square
10 // and upper triangular (or upper trapezoidal), with nnz (diag (T1)) == rank.
11 // On input, the row indices of the matrix R must be sorted.
12
13 #include "spqr.hpp"
14
spqr_trapezoidal(Long n,Long * Rp,Long * Ri,Entry * Rx,Long bncols,Long * Qfill,int skip_if_trapezoidal,Long ** p_Tp,Long ** p_Ti,Entry ** p_Tx,Long ** p_Qtrap,cholmod_common * cc)15 template <typename Entry> Long spqr_trapezoidal // rank of R; EMPTY on failure
16 (
17 // inputs, not modified
18
19 // FUTURE : make R and T cholmod_sparse:
20 Long n, // R is m-by-n (m is not needed here; can be economy R)
21 Long *Rp, // size n+1, column pointers of R
22 Long *Ri, // size rnz = Rp [n], row indices of R
23 Entry *Rx, // size rnz, numerical values of R
24
25 Long bncols, // number of columns of B
26
27 Long *Qfill, // size n+bncols, fill-reducing ordering. Qfill [k] = j if
28 // the jth column of A is the kth column of R. If Qfill is
29 // NULL, then it is assumed to be the identity
30 // permutation.
31
32 int skip_if_trapezoidal, // if R is already in trapezoidal form,
33 // and skip_if_trapezoidal is TRUE, then
34 // the matrix T is not created.
35
36 // outputs, not allocated on input
37 Long **p_Tp, // size n+1, column pointers of T
38 Long **p_Ti, // size rnz, row indices of T
39 Entry **p_Tx, // size rnz, numerical values of T
40
41 Long **p_Qtrap, // size n+bncols, modified Qfill
42
43 // workspace and parameters
44 cholmod_common *cc
45 )
46 {
47 Entry *Tx ;
48 Long *Tp, *Ti, *Qtrap ;
49 Long rnz, i, rank, k, p, pend, len, t1nz, t2nz, k1, k2, p1, p2, found_dead,
50 is_trapezoidal ;
51
52 // -------------------------------------------------------------------------
53 // find the rank of R, nnz(T1), and nnz(T2)
54 // -------------------------------------------------------------------------
55
56 rank = 0 ; // rank of R
57 t1nz = 0 ; // number of nonzeros in T1
58 t2nz = 0 ; // number of nonzeros in T2
59 found_dead = FALSE ; // true when first dead column is found
60 is_trapezoidal = TRUE ; // becomes false if is R not in trapezoidal form
61
62 *p_Tp = NULL ;
63 *p_Ti = NULL ;
64 *p_Tx = NULL ;
65 *p_Qtrap = NULL ;
66
67 for (k = 0 ; k < n ; k++)
68 {
69 // get the row index of the last entry in R (:,k)
70 p = Rp [k] ;
71 pend = Rp [k+1] ;
72 len = pend - p ;
73 i = (len > 0) ? Ri [pend - 1] : EMPTY ;
74
75 // determine status of column k
76 if (i > rank)
77 {
78 // R is not upper triangular, squeezed or otherwise. Do not create
79 // T and Qtrap; R is left in its original non-trapezoidal form.
80 PR (("R not upper, k = %ld\n", k)) ;
81 return (EMPTY) ;
82 }
83 else if (i == rank)
84 {
85 // column k is live
86 rank++ ;
87 t1nz += len ;
88 if (found_dead)
89 {
90 // If all live columns appear first (if any), then all dead
91 // columns, then the matrix is already in upper trapezoidal
92 // form. We just found a live column after one or more dead
93 // columns, so the matrix R is not upper trapezoidal.
94 is_trapezoidal = FALSE ;
95 }
96 }
97 else
98 {
99 // column k is dead
100 found_dead = TRUE ;
101 t2nz += len ;
102 }
103 }
104
105 // -------------------------------------------------------------------------
106 // quick return if already trapezoidal
107 // -------------------------------------------------------------------------
108
109 if (is_trapezoidal)
110 {
111 PR (("already trapezoidal\n")) ;
112 if (skip_if_trapezoidal)
113 {
114 // do not create T
115 return (rank) ;
116 }
117 }
118
119 // -------------------------------------------------------------------------
120 // allocate the results (T and Qtrap)
121 // -------------------------------------------------------------------------
122
123 rnz = Rp [n] ;
124
125 Tp = (Long *) cholmod_l_malloc (n+1, sizeof (Long), cc) ;
126 Ti = (Long *) cholmod_l_malloc (rnz, sizeof (Long), cc) ;
127 Tx = (Entry *) cholmod_l_malloc (rnz, sizeof (Entry), cc) ;
128 Qtrap = (Long *) cholmod_l_malloc (n+bncols, sizeof (Long), cc) ;
129
130 if (cc->status < CHOLMOD_OK)
131 {
132 // out of memory
133 cholmod_l_free (n+1, sizeof (Long), Tp, cc) ;
134 cholmod_l_free (rnz, sizeof (Long), Ti, cc) ;
135 cholmod_l_free (rnz, sizeof (Entry), Tx, cc) ;
136 cholmod_l_free (n+bncols, sizeof (Long), Qtrap, cc) ;
137 return (EMPTY) ;
138 }
139
140 PR (("rank %ld of T, nnz(T1) %ld nnz(T2) %ld\n", rank, t1nz, t2nz)) ;
141
142 // -------------------------------------------------------------------------
143 // find the column pointers Tp and permutation Qtrap
144 // -------------------------------------------------------------------------
145
146 k1 = 0 ; // first column of T1
147 k2 = rank ; // first column of T2
148 p1 = 0 ; // T1 starts at Ti,Tx [0]
149 p2 = t1nz ; // T2 starts at Ti,Tx [p2]
150 rank = 0 ; // rank of R
151
152 for (k = 0 ; k < n ; k++)
153 {
154 // get the row index of the last entry in R (:,k)
155 p = Rp [k] ;
156 pend = Rp [k+1] ;
157 len = pend - p ;
158 i = (len > 0) ? Ri [pend - 1] : EMPTY ;
159
160 // column k is live if i is a new row index
161 ASSERT (! (i > rank)) ;
162 if (i == rank)
163 {
164 // column k is live; place it in T1
165 rank++ ;
166 Tp [k1] = p1 ;
167 Qtrap [k1] = Qfill ? Qfill [k] : k ;
168 k1++ ;
169 for ( ; p < pend ; p++)
170 {
171 Ti [p1] = Ri [p] ;
172 Tx [p1] = Rx [p] ;
173 p1++ ;
174 }
175 }
176 else
177 {
178 // column k is dead; place it in T2
179 Tp [k2] = p2 ;
180 Qtrap [k2] = Qfill ? Qfill [k] : k ;
181 k2++ ;
182 for ( ; p < pend ; p++)
183 {
184 Ti [p2] = Ri [p] ;
185 Tx [p2] = Rx [p] ;
186 p2++ ;
187 }
188 }
189 }
190
191 for ( ; k < n+bncols ; k++)
192 {
193 Qtrap [k] = Qfill ? Qfill [k] : k ;
194 }
195
196 // -------------------------------------------------------------------------
197 // finalize the column pointers and return the results
198 // -------------------------------------------------------------------------
199
200 ASSERT (k1 == rank) ;
201 ASSERT (k2 == n) ;
202 ASSERT (p1 == t1nz) ;
203 ASSERT (p2 == rnz) ;
204 Tp [n] = rnz ;
205 *p_Tp = Tp ;
206 *p_Ti = Ti ;
207 *p_Tx = Tx ;
208 *p_Qtrap = Qtrap ;
209 return (rank) ;
210 }
211
212
213 // =============================================================================
214
215 template Long spqr_trapezoidal <double> // rank of R, or EMPTY on failure
216 (
217 // inputs, not modified
218
219 Long n, // R is m-by-n (m is not needed here; can be economy R)
220 Long *Rp, // size n+1, column pointers of R
221 Long *Ri, // size rnz = Rp [n], row indices of R
222 double *Rx, // size rnz, numerical values of R
223
224 Long bncols, // number of columns of B
225
226 Long *Qfill, // size n+bncols, fill-reducing ordering. Qfill [k] = j if
227 // the jth column of A is the kth column of R. If Qfill is
228 // NULL, then it is assumed to be the identity
229 // permutation.
230
231 int skip_if_trapezoidal, // if R is already in trapezoidal form,
232 // and skip_if_trapezoidal is TRUE, then
233 // the matrix T is not created.
234
235 // outputs, not allocated on input
236 Long **p_Tp, // size n+1, column pointers of T
237 Long **p_Ti, // size rnz, row indices of T
238 double **p_Tx, // size rnz, numerical values of T
239
240 Long **p_Qtrap, // size n+bncols, modified Qfill
241
242 // workspace and parameters
243 cholmod_common *cc
244 ) ;
245
246 // =============================================================================
247
248 template Long spqr_trapezoidal <Complex> // rank of R, or EMPTY on failure
249 (
250 // inputs, not modified
251
252 Long n, // R is m-by-n (m is not needed here; can be economy R)
253 Long *Rp, // size n+1, column pointers of R
254 Long *Ri, // size rnz = Rp [n], row indices of R
255 Complex *Rx, // size rnz, numerical values of R
256
257 Long bncols, // number of columns of B
258
259 Long *Qfill, // size n+bncols, fill-reducing ordering. Qfill [k] = j if
260 // the jth column of A is the kth column of R. If Qfill is
261 // NULL, then it is assumed to be the identity
262 // permutation.
263
264 int skip_if_trapezoidal, // if R is already in trapezoidal form,
265 // and skip_if_trapezoidal is TRUE, then
266 // the matrix T is not created.
267
268 // outputs, not allocated on input
269 Long **p_Tp, // size n+1, column pointers of T
270 Long **p_Ti, // size rnz, row indices of T
271 Complex **p_Tx, // size rnz, numerical values of T
272
273 Long **p_Qtrap, // size n+bncols, modified Qfill
274
275 // workspace and parameters
276 cholmod_common *cc
277 ) ;
278