1 /* ========================================================================== */
2 /* === Core/cholmod_aat ===================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 * Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* C = A*A' or C = A(:,f)*A(:,f)'
11 *
12 * A can be packed or unpacked, sorted or unsorted, but must be stored with
13 * both upper and lower parts (A->stype of zero). C is returned as packed,
14 * C->stype of zero (both upper and lower parts present), and unsorted. See
15 * cholmod_ssmult in the MatrixOps Module for a more general matrix-matrix
16 * multiply.
17 *
18 * You can trivially convert C into a symmetric upper/lower matrix by
19 * changing C->stype = 1 or -1 after calling this routine.
20 *
21 * workspace:
22 * Flag (A->nrow),
23 * Iwork (max (A->nrow, A->ncol)) if fset present,
24 * Iwork (A->nrow) if no fset,
25 * W (A->nrow) if mode > 0,
26 * allocates temporary copy for A'.
27 *
28 * A can be pattern or real. Complex or zomplex cases are supported only
29 * if the mode is <= 0 (in which case the numerical values are ignored).
30 */
31
32 #include "cholmod_internal.h"
33 #include "cholmod_core.h"
34
CHOLMOD(aat)35 cholmod_sparse *CHOLMOD(aat)
36 (
37 /* ---- input ---- */
38 cholmod_sparse *A, /* input matrix; C=A*A' is constructed */
39 Int *fset, /* subset of 0:(A->ncol)-1 */
40 size_t fsize, /* size of fset */
41 int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag)
42 * -2: pattern only, no diagonal, add 50% + n extra
43 * space to C */
44 /* --------------- */
45 cholmod_common *Common
46 )
47 {
48 double fjt ;
49 double *Ax, *Fx, *Cx, *W ;
50 Int *Ap, *Anz, *Ai, *Fp, *Fi, *Cp, *Ci, *Flag ;
51 cholmod_sparse *C, *F ;
52 Int packed, j, i, pa, paend, pf, pfend, n, mark, cnz, t, p, values, diag,
53 extra ;
54
55 /* ---------------------------------------------------------------------- */
56 /* check inputs */
57 /* ---------------------------------------------------------------------- */
58
59 RETURN_IF_NULL_COMMON (NULL) ;
60 RETURN_IF_NULL (A, NULL) ;
61 values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
62 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
63 values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
64 if (A->stype)
65 {
66 ERROR (CHOLMOD_INVALID, "matrix cannot be symmetric") ;
67 return (NULL) ;
68 }
69 Common->status = CHOLMOD_OK ;
70
71 /* ---------------------------------------------------------------------- */
72 /* allocate workspace */
73 /* ---------------------------------------------------------------------- */
74
75 diag = (mode >= 0) ;
76 n = A->nrow ;
77 CHOLMOD(allocate_work) (n, MAX (A->ncol, A->nrow), values ? n : 0, Common) ;
78 if (Common->status < CHOLMOD_OK)
79 {
80 return (NULL) ; /* out of memory */
81 }
82 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n : 0, Common)) ;
83
84 /* ---------------------------------------------------------------------- */
85 /* get inputs */
86 /* ---------------------------------------------------------------------- */
87
88 ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;
89
90 /* get the A matrix */
91 Ap = A->p ;
92 Anz = A->nz ;
93 Ai = A->i ;
94 Ax = A->x ;
95 packed = A->packed ;
96
97 /* get workspace */
98 W = Common->Xwork ; /* size n, unused if values is FALSE */
99 Flag = Common->Flag ; /* size n, Flag [0..n-1] < mark on input*/
100
101 /* ---------------------------------------------------------------------- */
102 /* F = A' or A(:,f)' */
103 /* ---------------------------------------------------------------------- */
104
105 /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
106 F = CHOLMOD(ptranspose) (A, values, NULL, fset, fsize, Common) ;
107 if (Common->status < CHOLMOD_OK)
108 {
109 return (NULL) ; /* out of memory */
110 }
111
112 Fp = F->p ;
113 Fi = F->i ;
114 Fx = F->x ;
115
116 /* ---------------------------------------------------------------------- */
117 /* count the number of entries in the result C */
118 /* ---------------------------------------------------------------------- */
119
120 cnz = 0 ;
121 for (j = 0 ; j < n ; j++)
122 {
123 /* clear the Flag array */
124 /* mark = CHOLMOD(clear_flag) (Common) ; */
125 CHOLMOD_CLEAR_FLAG (Common) ;
126 mark = Common->mark ;
127
128 /* exclude the diagonal, if requested */
129 if (!diag)
130 {
131 Flag [j] = mark ;
132 }
133
134 /* for each nonzero F(t,j) in column j, do: */
135 pfend = Fp [j+1] ;
136 for (pf = Fp [j] ; pf < pfend ; pf++)
137 {
138 /* F(t,j) is nonzero */
139 t = Fi [pf] ;
140
141 /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
142 pa = Ap [t] ;
143 paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
144 for ( ; pa < paend ; pa++)
145 {
146 i = Ai [pa] ;
147 if (Flag [i] != mark)
148 {
149 Flag [i] = mark ;
150 cnz++ ;
151 }
152 }
153 }
154 if (cnz < 0)
155 {
156 break ; /* integer overflow case */
157 }
158 }
159
160 extra = (mode == -2) ? (cnz/2 + n) : 0 ;
161
162 mark = CHOLMOD(clear_flag) (Common) ;
163
164 /* ---------------------------------------------------------------------- */
165 /* check for integer overflow */
166 /* ---------------------------------------------------------------------- */
167
168 if (cnz < 0 || (cnz + extra) < 0)
169 {
170 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
171 CHOLMOD(clear_flag) (Common) ;
172 CHOLMOD(free_sparse) (&F, Common) ;
173 return (NULL) ; /* problem too large */
174 }
175
176 /* ---------------------------------------------------------------------- */
177 /* allocate C */
178 /* ---------------------------------------------------------------------- */
179
180 C = CHOLMOD(allocate_sparse) (n, n, cnz + extra, FALSE, TRUE, 0,
181 values ? A->xtype : CHOLMOD_PATTERN, Common) ;
182 if (Common->status < CHOLMOD_OK)
183 {
184 CHOLMOD(free_sparse) (&F, Common) ;
185 return (NULL) ; /* out of memory */
186 }
187
188 Cp = C->p ;
189 Ci = C->i ;
190 Cx = C->x ;
191
192 /* ---------------------------------------------------------------------- */
193 /* C = A*A' */
194 /* ---------------------------------------------------------------------- */
195
196 cnz = 0 ;
197
198 if (values)
199 {
200
201 /* pattern and values */
202 for (j = 0 ; j < n ; j++)
203 {
204 /* clear the Flag array */
205 mark = CHOLMOD(clear_flag) (Common) ;
206
207 /* start column j of C */
208 Cp [j] = cnz ;
209
210 /* for each nonzero F(t,j) in column j, do: */
211 pfend = Fp [j+1] ;
212 for (pf = Fp [j] ; pf < pfend ; pf++)
213 {
214 /* F(t,j) is nonzero */
215 t = Fi [pf] ;
216 fjt = Fx [pf] ;
217
218 /* add the nonzero pattern of A(:,t) to the pattern of C(:,j)
219 * and scatter the values into W */
220 pa = Ap [t] ;
221 paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
222 for ( ; pa < paend ; pa++)
223 {
224 i = Ai [pa] ;
225 if (Flag [i] != mark)
226 {
227 Flag [i] = mark ;
228 Ci [cnz++] = i ;
229 }
230 W [i] += Ax [pa] * fjt ;
231 }
232 }
233
234 /* gather the values into C(:,j) */
235 for (p = Cp [j] ; p < cnz ; p++)
236 {
237 i = Ci [p] ;
238 Cx [p] = W [i] ;
239 W [i] = 0 ;
240 }
241 }
242
243 }
244 else
245 {
246
247 /* pattern only */
248 for (j = 0 ; j < n ; j++)
249 {
250 /* clear the Flag array */
251 mark = CHOLMOD(clear_flag) (Common) ;
252
253 /* exclude the diagonal, if requested */
254 if (!diag)
255 {
256 Flag [j] = mark ;
257 }
258
259 /* start column j of C */
260 Cp [j] = cnz ;
261
262 /* for each nonzero F(t,j) in column j, do: */
263 pfend = Fp [j+1] ;
264 for (pf = Fp [j] ; pf < pfend ; pf++)
265 {
266 /* F(t,j) is nonzero */
267 t = Fi [pf] ;
268
269 /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
270 pa = Ap [t] ;
271 paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
272 for ( ; pa < paend ; pa++)
273 {
274 i = Ai [pa] ;
275 if (Flag [i] != mark)
276 {
277 Flag [i] = mark ;
278 Ci [cnz++] = i ;
279 }
280 }
281 }
282 }
283 }
284
285 Cp [n] = cnz ;
286 ASSERT (IMPLIES (mode != -2, MAX (1,cnz) == C->nzmax)) ;
287
288 /* ---------------------------------------------------------------------- */
289 /* clear workspace and free temporary matrices and return result */
290 /* ---------------------------------------------------------------------- */
291
292 CHOLMOD(free_sparse) (&F, Common) ;
293 CHOLMOD(clear_flag) (Common) ;
294 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n : 0, Common)) ;
295 DEBUG (i = CHOLMOD(dump_sparse) (C, "aat", Common)) ;
296 ASSERT (IMPLIES (mode < 0, i == 0)) ;
297 return (C) ;
298 }
299