1 /* ========================================================================== */
2 /* === MatrixOps/cholmod_submatrix ========================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MatrixOps Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* C = A (rset,cset), where C becomes length(rset)-by-length(cset) in dimension.
11 * rset and cset can have duplicate entries. A and C must be unsymmetric. C
12 * is packed. If the sorted flag is TRUE on input, or rset is sorted and A is
13 * sorted, then C is sorted; otherwise C is unsorted.
14 *
15 * A NULL rset or cset means "[ ]" in MATLAB notation.
16 * If the length of rset or cset is negative, it denotes ":" in MATLAB notation.
17 *
18 * For permuting a matrix, this routine is an alternative to cholmod_ptranspose
19 * (which permutes and transposes a matrix and can work on symmetric matrices).
20 *
21 * The time taken by this routine is O(A->nrow) if the Common workspace needs
22 * to be initialized, plus O(C->nrow + C->ncol + nnz (A (:,cset))). Thus, if C
23 * is small and the workspace is not initialized, the time can be dominated by
24 * the call to cholmod_allocate_work. However, once the workspace is
25 * allocated, subsequent calls take less time.
26 *
27 * workspace: Iwork (max (A->nrow + length (rset), length (cset))).
28 * allocates temporary copy of C if it is to be returned sorted.
29 *
30 * Future work: A common case occurs where A has sorted columns, and rset is in
31 * the form lo:hi in MATLAB notation. This routine could exploit that case
32 * to run even faster when the matrix is sorted, particularly when lo is small.
33 *
34 * Only pattern and real matrices are supported. Complex and zomplex matrices
35 * are supported only when "values" is FALSE.
36 */
37
38 #ifndef NGPL
39 #ifndef NMATRIXOPS
40
41 #include "cholmod_internal.h"
42 #include "cholmod_matrixops.h"
43
44 /* ========================================================================== */
45 /* === check_subset ========================================================= */
46 /* ========================================================================== */
47
48 /* Check the rset or cset, and return TRUE if valid, FALSE if invalid */
49
check_subset(Int * set,Int len,Int n)50 static int check_subset (Int *set, Int len, Int n)
51 {
52 Int k ;
53 if (set == NULL)
54 {
55 return (TRUE) ;
56 }
57 for (k = 0 ; k < len ; k++)
58 {
59 if (set [k] < 0 || set [k] >= n)
60 {
61 return (FALSE) ;
62 }
63 }
64 return (TRUE) ;
65 }
66
67
68 /* ========================================================================== */
69 /* === cholmod_submatrix ==================================================== */
70 /* ========================================================================== */
71
CHOLMOD(submatrix)72 cholmod_sparse *CHOLMOD(submatrix)
73 (
74 /* ---- input ---- */
75 cholmod_sparse *A, /* matrix to subreference */
76 Int *rset, /* set of row indices, duplicates OK */
77 SuiteSparse_long rsize, /* size of rset, or -1 for ":" */
78 Int *cset, /* set of column indices, duplicates OK */
79 SuiteSparse_long csize, /* size of cset, or -1 for ":" */
80 int values, /* if TRUE compute the numerical values of C */
81 int sorted, /* if TRUE then return C with sorted columns */
82 /* --------------- */
83 cholmod_common *Common
84 )
85 {
86 double aij = 0 ;
87 double *Ax, *Cx ;
88 Int *Ap, *Ai, *Anz, *Ci, *Cp, *Head, *Rlen, *Rnext, *Iwork ;
89 cholmod_sparse *C ;
90 Int packed, ancol, anrow, cnrow, cncol, nnz, i, j, csorted, ilast, p,
91 pend, pdest, ci, cj, head, nr, nc ;
92 size_t s ;
93 int ok = TRUE ;
94
95 /* ---------------------------------------------------------------------- */
96 /* check inputs */
97 /* ---------------------------------------------------------------------- */
98
99 RETURN_IF_NULL_COMMON (NULL) ;
100 RETURN_IF_NULL (A, NULL) ;
101 values = (values && (A->xtype != CHOLMOD_PATTERN)) ;
102 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
103 values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
104 if (A->stype != 0)
105 {
106 /* A must be unsymmetric */
107 ERROR (CHOLMOD_INVALID, "symmetric upper or lower case not supported") ;
108 return (NULL) ;
109 }
110 Common->status = CHOLMOD_OK ;
111
112 /* ---------------------------------------------------------------------- */
113 /* allocate workspace */
114 /* ---------------------------------------------------------------------- */
115
116 ancol = A->ncol ;
117 anrow = A->nrow ;
118 nr = rsize ;
119 nc = csize ;
120 if (rset == NULL)
121 {
122 /* nr = 0 denotes rset = [ ], nr < 0 denotes rset = 0:anrow-1 */
123 nr = (nr < 0) ? (-1) : 0 ;
124 }
125 if (cset == NULL)
126 {
127 /* nr = 0 denotes cset = [ ], nr < 0 denotes cset = 0:ancol-1 */
128 nc = (nc < 0) ? (-1) : 0 ;
129 }
130 cnrow = (nr < 0) ? anrow : nr ; /* negative rset means rset = 0:anrow-1 */
131 cncol = (nc < 0) ? ancol : nc ; /* negative cset means cset = 0:ancol-1 */
132
133 if (nr < 0 && nc < 0)
134 {
135
136 /* ------------------------------------------------------------------ */
137 /* C = A (:,:), use cholmod_copy instead */
138 /* ------------------------------------------------------------------ */
139
140 /* workspace: Iwork (max (C->nrow,C->ncol)) */
141 PRINT1 (("submatrix C = A (:,:)\n")) ;
142 C = CHOLMOD(copy) (A, 0, values, Common) ;
143 if (Common->status < CHOLMOD_OK)
144 {
145 /* out of memory */
146 return (NULL) ;
147 }
148 return (C) ;
149 }
150 PRINT1 (("submatrix nr "ID" nc "ID" Cnrow "ID" Cncol "ID""
151 " Anrow "ID" Ancol "ID"\n", nr, nc, cnrow, cncol, anrow, ancol)) ;
152
153 /* s = MAX3 (anrow+MAX(0,nr), cncol, cnrow) ; */
154 s = CHOLMOD(add_size_t) (anrow, MAX (0,nr), &ok) ;
155 if (!ok)
156 {
157 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
158 return (NULL) ;
159 }
160 s = MAX3 (s, ((size_t) cncol), ((size_t) cnrow)) ;
161
162 CHOLMOD(allocate_work) (anrow, s, 0, Common) ;
163 if (Common->status < CHOLMOD_OK)
164 {
165 /* out of memory */
166 return (NULL) ;
167 }
168
169 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
170
171 /* ---------------------------------------------------------------------- */
172 /* get inputs */
173 /* ---------------------------------------------------------------------- */
174
175 Ap = A->p ;
176 Anz = A->nz ;
177 Ai = A->i ;
178 Ax = A->x ;
179 packed = A->packed ;
180
181 /* ---------------------------------------------------------------------- */
182 /* get workspace */
183 /* ---------------------------------------------------------------------- */
184
185 Head = Common->Head ; /* size anrow */
186 Iwork = Common->Iwork ;
187 Rlen = Iwork ; /* size anrow (i/i/l) */
188 Rnext = Iwork + anrow ; /* size nr (i/i/l), not used if nr < 0 */
189
190 /* ---------------------------------------------------------------------- */
191 /* construct inverse of rset and compute nnz (C) */
192 /* ---------------------------------------------------------------------- */
193
194 PRINT1 (("nr "ID" nc "ID"\n", nr, nc)) ;
195 PRINT1 (("anrow "ID" ancol "ID"\n", anrow, ancol)) ;
196 PRINT1 (("cnrow "ID" cncol "ID"\n", cnrow, cncol)) ;
197 DEBUG (for (i = 0 ; i < nr ; i++) PRINT2 (("rset ["ID"] = "ID"\n",
198 i, rset [i])));
199 DEBUG (for (i = 0 ; i < nc ; i++) PRINT2 (("cset ["ID"] = "ID"\n",
200 i, cset [i])));
201
202 /* C is sorted if A and rset are sorted, or if C has one row or less */
203 csorted = A->sorted || (cnrow <= 1) ;
204
205 if (!check_subset (rset, nr, anrow))
206 {
207 ERROR (CHOLMOD_INVALID, "invalid rset") ;
208 return (NULL) ;
209 }
210
211 if (!check_subset (cset, nc, ancol))
212 {
213 ERROR (CHOLMOD_INVALID, "invalid cset") ;
214 return (NULL) ;
215 }
216
217 nnz = 0 ;
218 if (nr < 0)
219 {
220 /* C = A (:,cset) where cset = [ ] or cset is not empty */
221 ASSERT (IMPLIES (cncol > 0, cset != NULL)) ;
222 for (cj = 0 ; cj < cncol ; cj++)
223 {
224 /* construct column cj of C, which is column j of A */
225 j = cset [cj] ;
226 nnz += (packed) ? (Ap [j+1] - Ap [j]) : MAX (0, Anz [j]) ;
227 }
228 }
229 else
230 {
231 /* C = A (rset,cset), where rset is not empty but cset might be empty */
232 /* create link lists in reverse order to preserve natural order */
233 ilast = anrow ;
234 for (ci = nr-1 ; ci >= 0 ; ci--)
235 {
236 /* row i of A becomes row ci of C; add ci to ith link list */
237 i = rset [ci] ;
238 head = Head [i] ;
239 Rlen [i] = (head == EMPTY) ? 1 : (Rlen [i] + 1) ;
240 Rnext [ci] = head ;
241 Head [i] = ci ;
242 if (i > ilast)
243 {
244 /* row indices in columns of C will not be sorted */
245 csorted = FALSE ;
246 }
247 ilast = i ;
248 }
249
250 #ifndef NDEBUG
251 for (i = 0 ; i < anrow ; i++)
252 {
253 Int k = 0 ;
254 Int rlen = (Head [i] != EMPTY) ? Rlen [i] : -1 ;
255 PRINT1 (("Row "ID" Rlen "ID": ", i, rlen)) ;
256 for (ci = Head [i] ; ci != EMPTY ; ci = Rnext [ci])
257 {
258 k++ ;
259 PRINT2 ((""ID" ", ci)) ;
260 }
261 PRINT1 (("\n")) ;
262 ASSERT (IMPLIES (Head [i] != EMPTY, k == Rlen [i])) ;
263 }
264 #endif
265
266 /* count nonzeros in C */
267 for (cj = 0 ; cj < cncol ; cj++)
268 {
269 /* count rows in column cj of C, which is column j of A */
270 j = (nc < 0) ? cj : (cset [cj]) ;
271 p = Ap [j] ;
272 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
273 for ( ; p < pend ; p++)
274 {
275 /* row i of A becomes multiple rows (ci) of C */
276 i = Ai [p] ;
277 ASSERT (i >= 0 && i < anrow) ;
278 if (Head [i] != EMPTY)
279 {
280 nnz += Rlen [i] ;
281 }
282 }
283 }
284 }
285 PRINT1 (("nnz (C) "ID"\n", nnz)) ;
286
287 /* rset and cset are now valid */
288 DEBUG (CHOLMOD(dump_subset) (rset, rsize, anrow, "rset", Common)) ;
289 DEBUG (CHOLMOD(dump_subset) (cset, csize, ancol, "cset", Common)) ;
290
291 /* ---------------------------------------------------------------------- */
292 /* allocate C */
293 /* ---------------------------------------------------------------------- */
294
295 C = CHOLMOD(allocate_sparse) (cnrow, cncol, nnz, csorted, TRUE, 0,
296 values ? A->xtype : CHOLMOD_PATTERN, Common) ;
297 if (Common->status < CHOLMOD_OK)
298 {
299 /* out of memory */
300 for (i = 0 ; i < anrow ; i++)
301 {
302 Head [i] = EMPTY ;
303 }
304 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
305 return (NULL) ;
306 }
307
308 Cp = C->p ;
309 Ci = C->i ;
310 Cx = C->x ;
311
312 /* ---------------------------------------------------------------------- */
313 /* C = A (rset,cset) */
314 /* ---------------------------------------------------------------------- */
315
316 pdest = 0 ;
317 if (nnz == 0)
318 {
319 /* C has no nonzeros */
320 for (cj = 0 ; cj <= cncol ; cj++)
321 {
322 Cp [cj] = 0 ;
323 }
324 }
325 else if (nr < 0)
326 {
327 /* C = A (:,cset), where cset is not empty */
328 for (cj = 0 ; cj < cncol ; cj++)
329 {
330 /* construct column cj of C, which is column j of A */
331 PRINT1 (("construct cj = j = "ID"\n", cj)) ;
332 j = cset [cj] ;
333 Cp [cj] = pdest ;
334 p = Ap [j] ;
335 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
336 for ( ; p < pend ; p++)
337 {
338 Ci [pdest] = Ai [p] ;
339 if (values)
340 {
341 Cx [pdest] = Ax [p] ;
342 }
343 pdest++ ;
344 ASSERT (pdest <= nnz) ;
345 }
346 }
347 }
348 else
349 {
350 /* C = A (rset,cset), where rset is not empty but cset might be empty */
351 for (cj = 0 ; cj < cncol ; cj++)
352 {
353 /* construct column cj of C, which is column j of A */
354 PRINT1 (("construct cj = "ID"\n", cj)) ;
355 j = (nc < 0) ? cj : (cset [cj]) ;
356 PRINT1 (("cj = "ID"\n", j)) ;
357 Cp [cj] = pdest ;
358 p = Ap [j] ;
359 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
360 for ( ; p < pend ; p++)
361 {
362 /* row (Ai [p]) of A becomes multiple rows (ci) of C */
363 PRINT2 (("i: "ID" becomes: ", Ai [p])) ;
364 if (values)
365 {
366 aij = Ax [p] ;
367 }
368 for (ci = Head [Ai [p]] ; ci != EMPTY ; ci = Rnext [ci])
369 {
370 PRINT3 ((""ID" ", ci)) ;
371 Ci [pdest] = ci ;
372 if (values)
373 {
374 Cx [pdest] = aij ;
375 }
376 pdest++ ;
377 ASSERT (pdest <= nnz) ;
378 }
379 PRINT2 (("\n")) ;
380 }
381 }
382 }
383 Cp [cncol] = pdest ;
384 ASSERT (nnz == pdest) ;
385
386 /* ---------------------------------------------------------------------- */
387 /* clear workspace */
388 /* ---------------------------------------------------------------------- */
389
390 for (ci = 0 ; ci < nr ; ci++)
391 {
392 Head [rset [ci]] = EMPTY ;
393 }
394
395 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
396
397 /* ---------------------------------------------------------------------- */
398 /* sort C, if requested */
399 /* ---------------------------------------------------------------------- */
400
401 ASSERT (CHOLMOD(dump_sparse) (C , "C before sort", Common) >= 0) ;
402
403 if (sorted && !csorted)
404 {
405 /* workspace: Iwork (max (C->nrow,C->ncol)) */
406 if (!CHOLMOD(sort) (C, Common))
407 {
408 /* out of memory */
409 CHOLMOD(free_sparse) (&C, Common) ;
410 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
411 return (NULL) ;
412 }
413 }
414
415 /* ---------------------------------------------------------------------- */
416 /* return result */
417 /* ---------------------------------------------------------------------- */
418
419 ASSERT (CHOLMOD(dump_sparse) (C , "Final C", Common) >= 0) ;
420 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
421 return (C) ;
422 }
423 #endif
424 #endif
425