1 /* ========================================================================== */
2 /* === Core/cholmod_band ==================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 * Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* C = tril (triu (A,k1), k2)
11 *
12 * C is a matrix consisting of the diagonals of A from k1 to k2.
13 *
14 * k=0 is the main diagonal of A, k=1 is the superdiagonal, k=-1 is the
15 * subdiagonal, and so on. If A is m-by-n, then:
16 *
17 * k1=-m C = tril (A,k2)
18 * k2=n C = triu (A,k1)
19 * k1=0 and k2=0 C = diag(A), except C is a matrix, not a vector
20 *
21 * Values of k1 and k2 less than -m are treated as -m, and values greater
22 * than n are treated as n.
23 *
24 * A can be of any symmetry (upper, lower, or unsymmetric); C is returned in
25 * the same form, and packed. If A->stype > 0, entries in the lower
26 * triangular part of A are ignored, and the opposite is true if
27 * A->stype < 0. If A has sorted columns, then so does C.
28 * C has the same size as A.
29 *
30 * If inplace is TRUE, then the matrix A is modified in place.
31 * Only packed matrices can be converted in place.
32 *
33 * C can be returned as a numerical valued matrix (if A has numerical values
34 * and mode > 0), as a pattern-only (mode == 0), or as a pattern-only but with
35 * the diagonal entries removed (mode < 0).
36 *
37 * workspace: none
38 *
39 * A can have an xtype of pattern or real. Complex and zomplex cases supported
40 * only if mode <= 0 (in which case the numerical values are ignored).
41 */
42
43 #include "cholmod_internal.h"
44 #include "cholmod_core.h"
45
band(cholmod_sparse * A,SuiteSparse_long k1,SuiteSparse_long k2,int mode,int inplace,cholmod_common * Common)46 static cholmod_sparse *band /* returns C, or NULL if failure */
47 (
48 /* ---- input or in/out if inplace is TRUE --- */
49 cholmod_sparse *A,
50 /* ---- input ---- */
51 SuiteSparse_long k1, /* ignore entries below the k1-st diagonal */
52 SuiteSparse_long k2, /* ignore entries above the k2-nd diagonal */
53 int mode, /* >0: numerical, 0: pattern, <0: pattern (no diagonal) */
54 int inplace, /* if TRUE, then convert A in place */
55 /* --------------- */
56 cholmod_common *Common
57 )
58 {
59 double *Ax, *Cx ;
60 Int packed, nz, j, p, pend, i, ncol, nrow, jlo, jhi, ilo, ihi, sorted,
61 values, diag ;
62 Int *Ap, *Anz, *Ai, *Cp, *Ci ;
63 cholmod_sparse *C ;
64
65 /* ---------------------------------------------------------------------- */
66 /* check inputs */
67 /* ---------------------------------------------------------------------- */
68
69 RETURN_IF_NULL_COMMON (NULL) ;
70 RETURN_IF_NULL (A, NULL) ;
71 values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
72 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
73 values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
74 packed = A->packed ;
75 diag = (mode >= 0) ;
76 if (inplace && !packed)
77 {
78 /* cannot operate on an unpacked matrix in place */
79 ERROR (CHOLMOD_INVALID, "cannot operate on unpacked matrix in-place") ;
80 return (NULL) ;
81 }
82 Common->status = CHOLMOD_OK ;
83
84 /* ---------------------------------------------------------------------- */
85 /* get inputs */
86 /* ---------------------------------------------------------------------- */
87
88
89 PRINT1 (("k1 %ld k2 %ld\n", k1, k2)) ;
90 Ap = A->p ;
91 Anz = A->nz ;
92 Ai = A->i ;
93 Ax = A->x ;
94 sorted = A->sorted ;
95
96
97 if (A->stype > 0)
98 {
99 /* ignore any entries in strictly lower triangular part of A */
100 k1 = MAX (k1, 0) ;
101 }
102 if (A->stype < 0)
103 {
104 /* ignore any entries in strictly upper triangular part of A */
105 k2 = MIN (k2, 0) ;
106 }
107 ncol = A->ncol ;
108 nrow = A->nrow ;
109
110 /* ensure k1 and k2 are in the range -nrow to +ncol to
111 * avoid possible integer overflow if k1 and k2 are huge */
112 k1 = MAX (-nrow, k1) ;
113 k1 = MIN (k1, ncol) ;
114 k2 = MAX (-nrow, k2) ;
115 k2 = MIN (k2, ncol) ;
116
117 /* consider columns jlo to jhi. columns outside this range are empty */
118 jlo = MAX (k1, 0) ;
119 jhi = MIN (k2+nrow, ncol) ;
120
121 if (k1 > k2)
122 {
123 /* nothing to do */
124 jlo = ncol ;
125 jhi = ncol ;
126 }
127
128 /* ---------------------------------------------------------------------- */
129 /* allocate C, or operate on A in place */
130 /* ---------------------------------------------------------------------- */
131
132 if (inplace)
133 {
134 /* convert A in place */
135 C = A ;
136 }
137 else
138 {
139 /* count the number of entries in the result C */
140 nz = 0 ;
141 if (sorted)
142 {
143 for (j = jlo ; j < jhi ; j++)
144 {
145 ilo = j-k2 ;
146 ihi = j-k1 ;
147 p = Ap [j] ;
148 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
149 for ( ; p < pend ; p++)
150 {
151 i = Ai [p] ;
152 if (i > ihi)
153 {
154 break ;
155 }
156 if (i >= ilo && (diag || i != j))
157 {
158 nz++ ;
159 }
160 }
161 }
162 }
163 else
164 {
165 for (j = jlo ; j < jhi ; j++)
166 {
167 ilo = j-k2 ;
168 ihi = j-k1 ;
169 p = Ap [j] ;
170 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
171 for ( ; p < pend ; p++)
172 {
173 i = Ai [p] ;
174 if (i >= ilo && i <= ihi && (diag || i != j))
175 {
176 nz++ ;
177 }
178 }
179 }
180 }
181 /* allocate C; A will not be modified. C is sorted if A is sorted */
182 C = CHOLMOD(allocate_sparse) (A->nrow, ncol, nz, sorted, TRUE,
183 A->stype, values ? A->xtype : CHOLMOD_PATTERN, Common) ;
184 if (Common->status < CHOLMOD_OK)
185 {
186 return (NULL) ; /* out of memory */
187 }
188 }
189
190 Cp = C->p ;
191 Ci = C->i ;
192 Cx = C->x ;
193
194 /* ---------------------------------------------------------------------- */
195 /* construct C */
196 /* ---------------------------------------------------------------------- */
197
198 /* columns 0 to jlo-1 are empty */
199 for (j = 0 ; j < jlo ; j++)
200 {
201 Cp [j] = 0 ;
202 }
203
204 nz = 0 ;
205 if (sorted)
206 {
207 if (values)
208 {
209 /* pattern and values */
210 ASSERT (diag) ;
211 for (j = jlo ; j < jhi ; j++)
212 {
213 ilo = j-k2 ;
214 ihi = j-k1 ;
215 p = Ap [j] ;
216 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
217 Cp [j] = nz ;
218 for ( ; p < pend ; p++)
219 {
220 i = Ai [p] ;
221 if (i > ihi)
222 {
223 break ;
224 }
225 if (i >= ilo)
226 {
227 Ci [nz] = i ;
228 Cx [nz] = Ax [p] ;
229 nz++ ;
230 }
231 }
232 }
233 }
234 else
235 {
236 /* pattern only, perhaps with no diagonal */
237 for (j = jlo ; j < jhi ; j++)
238 {
239 ilo = j-k2 ;
240 ihi = j-k1 ;
241 p = Ap [j] ;
242 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
243 Cp [j] = nz ;
244 for ( ; p < pend ; p++)
245 {
246 i = Ai [p] ;
247 if (i > ihi)
248 {
249 break ;
250 }
251 if (i >= ilo && (diag || i != j))
252 {
253 Ci [nz++] = i ;
254 }
255 }
256 }
257 }
258 }
259 else
260 {
261 if (values)
262 {
263 /* pattern and values */
264 ASSERT (diag) ;
265 for (j = jlo ; j < jhi ; j++)
266 {
267 ilo = j-k2 ;
268 ihi = j-k1 ;
269 p = Ap [j] ;
270 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
271 Cp [j] = nz ;
272 for ( ; p < pend ; p++)
273 {
274 i = Ai [p] ;
275 if (i >= ilo && i <= ihi)
276 {
277 Ci [nz] = i ;
278 Cx [nz] = Ax [p] ;
279 nz++ ;
280 }
281 }
282 }
283 }
284 else
285 {
286 /* pattern only, perhaps with no diagonal */
287 for (j = jlo ; j < jhi ; j++)
288 {
289 ilo = j-k2 ;
290 ihi = j-k1 ;
291 p = Ap [j] ;
292 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
293 Cp [j] = nz ;
294 for ( ; p < pend ; p++)
295 {
296 i = Ai [p] ;
297 if (i >= ilo && i <= ihi && (diag || i != j))
298 {
299 Ci [nz++] = i ;
300 }
301 }
302 }
303 }
304 }
305
306 /* columns jhi to ncol-1 are empty */
307 for (j = jhi ; j <= ncol ; j++)
308 {
309 Cp [j] = nz ;
310 }
311
312 /* ---------------------------------------------------------------------- */
313 /* reduce A in size if done in place */
314 /* ---------------------------------------------------------------------- */
315
316 if (inplace)
317 {
318 /* free the unused parts of A, and reduce A->i and A->x in size */
319 ASSERT (MAX (1,nz) <= A->nzmax) ;
320 CHOLMOD(reallocate_sparse) (nz, A, Common) ;
321 ASSERT (Common->status >= CHOLMOD_OK) ;
322 }
323
324 /* ---------------------------------------------------------------------- */
325 /* return the result C */
326 /* ---------------------------------------------------------------------- */
327
328 DEBUG (i = CHOLMOD(dump_sparse) (C, "band", Common)) ;
329 ASSERT (IMPLIES (mode < 0, i == 0)) ;
330 return (C) ;
331 }
332
333
334 /* ========================================================================== */
335 /* === cholmod_band ========================================================= */
336 /* ========================================================================== */
337
CHOLMOD(band)338 cholmod_sparse *CHOLMOD(band)
339 (
340 /* ---- input ---- */
341 cholmod_sparse *A, /* matrix to extract band matrix from */
342 SuiteSparse_long k1, /* ignore entries below the k1-st diagonal */
343 SuiteSparse_long k2, /* ignore entries above the k2-nd diagonal */
344 int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */
345 /* --------------- */
346 cholmod_common *Common
347 )
348 {
349 return (band (A, k1, k2, mode, FALSE, Common)) ;
350 }
351
352
353 /* ========================================================================== */
354 /* === cholmod_band_inplace ================================================= */
355 /* ========================================================================== */
356
CHOLMOD(band_inplace)357 int CHOLMOD(band_inplace)
358 (
359 /* ---- input ---- */
360 SuiteSparse_long k1, /* ignore entries below the k1-st diagonal */
361 SuiteSparse_long k2, /* ignore entries above the k2-nd diagonal */
362 int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */
363 /* ---- in/out --- */
364 cholmod_sparse *A, /* matrix from which entries not in band are removed */
365 /* --------------- */
366 cholmod_common *Common
367 )
368 {
369 return (band (A, k1, k2, mode, TRUE, Common) != NULL) ;
370 }
371