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