1 /* ========================================================================== */
2 /* === Core/cholmod_add ===================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 * Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* C = alpha*A + beta*B, or spones(A+B). Result is packed, with sorted or
11 * unsorted columns. This routine is much faster and takes less memory if C
12 * is allowed to have unsorted columns.
13 *
14 * If A and B are both symmetric (in upper form) then C is the same. Likewise,
15 * if A and B are both symmetric (in lower form) then C is the same.
16 * Otherwise, C is unsymmetric. A and B must have the same dimension.
17 *
18 * workspace: Flag (nrow), W (nrow) if values, Iwork (max (nrow,ncol)).
19 * allocates temporary copies for A and B if they are symmetric.
20 * allocates temporary copy of C if it is to be returned sorted.
21 *
22 * A and B can have an xtype of pattern or real. Complex or zomplex cases
23 * are supported only if the "values" input parameter is FALSE.
24 */
25
26 #include "cholmod_internal.h"
27 #include "cholmod_core.h"
28
CHOLMOD(add)29 cholmod_sparse *CHOLMOD(add)
30 (
31 /* ---- input ---- */
32 cholmod_sparse *A, /* matrix to add */
33 cholmod_sparse *B, /* matrix to add */
34 double alpha [2], /* scale factor for A */
35 double beta [2], /* scale factor for B */
36 int values, /* if TRUE compute the numerical values of C */
37 int sorted, /* if TRUE, sort columns of C */
38 /* --------------- */
39 cholmod_common *Common
40 )
41 {
42 double *Ax, *Bx, *Cx, *W ;
43 Int apacked, up, lo, nrow, ncol, bpacked, nzmax, pa, paend, pb, pbend, i,
44 j, p, mark, nz ;
45 Int *Ap, *Ai, *Anz, *Bp, *Bi, *Bnz, *Flag, *Cp, *Ci ;
46 cholmod_sparse *A2, *B2, *C ;
47
48 /* ---------------------------------------------------------------------- */
49 /* check inputs */
50 /* ---------------------------------------------------------------------- */
51
52 RETURN_IF_NULL_COMMON (NULL) ;
53 RETURN_IF_NULL (A, NULL) ;
54 RETURN_IF_NULL (B, NULL) ;
55 values = values &&
56 (A->xtype != CHOLMOD_PATTERN) && (B->xtype != CHOLMOD_PATTERN) ;
57 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
58 values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
59 RETURN_IF_XTYPE_INVALID (B, CHOLMOD_PATTERN,
60 values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
61 if (A->nrow != B->nrow || A->ncol != B->ncol)
62 {
63 /* A and B must have the same dimensions */
64 ERROR (CHOLMOD_INVALID, "A and B dimesions do not match") ;
65 return (NULL) ;
66 }
67 /* A and B must have the same numerical type if values is TRUE (both must
68 * be CHOLMOD_REAL, this is implicitly checked above) */
69
70 Common->status = CHOLMOD_OK ;
71
72 /* ---------------------------------------------------------------------- */
73 /* allocate workspace */
74 /* ---------------------------------------------------------------------- */
75
76 nrow = A->nrow ;
77 ncol = A->ncol ;
78 CHOLMOD(allocate_work) (nrow, MAX (nrow,ncol), values ? nrow : 0, Common) ;
79 if (Common->status < CHOLMOD_OK)
80 {
81 return (NULL) ; /* out of memory */
82 }
83
84 /* ---------------------------------------------------------------------- */
85 /* get inputs */
86 /* ---------------------------------------------------------------------- */
87
88 if (nrow <= 1)
89 {
90 /* C will be implicitly sorted, so no need to sort it here */
91 sorted = FALSE ;
92 }
93
94 /* convert A or B to unsymmetric, if necessary */
95 A2 = NULL ;
96 B2 = NULL ;
97
98 if (A->stype != B->stype)
99 {
100 if (A->stype)
101 {
102 /* workspace: Iwork (max (nrow,ncol)) */
103 A2 = CHOLMOD(copy) (A, 0, values, Common) ;
104 if (Common->status < CHOLMOD_OK)
105 {
106 return (NULL) ; /* out of memory */
107 }
108 A = A2 ;
109 }
110 if (B->stype)
111 {
112 /* workspace: Iwork (max (nrow,ncol)) */
113 B2 = CHOLMOD(copy) (B, 0, values, Common) ;
114 if (Common->status < CHOLMOD_OK)
115 {
116 CHOLMOD(free_sparse) (&A2, Common) ;
117 return (NULL) ; /* out of memory */
118 }
119 B = B2 ;
120 }
121 }
122
123 /* get the A matrix */
124 ASSERT (A->stype == B->stype) ;
125 up = (A->stype > 0) ;
126 lo = (A->stype < 0) ;
127
128 Ap = A->p ;
129 Anz = A->nz ;
130 Ai = A->i ;
131 Ax = A->x ;
132 apacked = A->packed ;
133
134 /* get the B matrix */
135 Bp = B->p ;
136 Bnz = B->nz ;
137 Bi = B->i ;
138 Bx = B->x ;
139 bpacked = B->packed ;
140
141 /* get workspace */
142 W = Common->Xwork ; /* size nrow, used if values is TRUE */
143 Flag = Common->Flag ; /* size nrow, Flag [0..nrow-1] < mark on input */
144
145 /* ---------------------------------------------------------------------- */
146 /* allocate the result C */
147 /* ---------------------------------------------------------------------- */
148
149 /* If integer overflow occurs, nzmax < 0 and the allocate fails properly
150 * (likewise in most other matrix manipulation routines). */
151
152 nzmax = CHOLMOD(nnz) (A, Common) + CHOLMOD(nnz) (B, Common) ;
153
154 C = CHOLMOD(allocate_sparse) (nrow, ncol, nzmax, FALSE, TRUE,
155 SIGN (A->stype), values ? A->xtype : CHOLMOD_PATTERN, Common) ;
156 if (Common->status < CHOLMOD_OK)
157 {
158 CHOLMOD(free_sparse) (&A2, Common) ;
159 CHOLMOD(free_sparse) (&B2, Common) ;
160 return (NULL) ; /* out of memory */
161 }
162
163 Cp = C->p ;
164 Ci = C->i ;
165 Cx = C->x ;
166
167 /* ---------------------------------------------------------------------- */
168 /* compute C = alpha*A + beta*B */
169 /* ---------------------------------------------------------------------- */
170
171 nz = 0 ;
172 for (j = 0 ; j < ncol ; j++)
173 {
174 Cp [j] = nz ;
175
176 /* clear the Flag array */
177 /* mark = CHOLMOD(clear_flag) (Common) ; */
178 CHOLMOD_CLEAR_FLAG (Common) ;
179 mark = Common->mark ;
180
181 /* scatter B into W */
182 pb = Bp [j] ;
183 pbend = (bpacked) ? (Bp [j+1]) : (pb + Bnz [j]) ;
184 for (p = pb ; p < pbend ; p++)
185 {
186 i = Bi [p] ;
187 if ((up && i > j) || (lo && i < j))
188 {
189 continue ;
190 }
191 Flag [i] = mark ;
192 if (values)
193 {
194 W [i] = beta [0] * Bx [p] ;
195 }
196 }
197
198 /* add A and gather from W into C(:,j) */
199 pa = Ap [j] ;
200 paend = (apacked) ? (Ap [j+1]) : (pa + Anz [j]) ;
201 for (p = pa ; p < paend ; p++)
202 {
203 i = Ai [p] ;
204 if ((up && i > j) || (lo && i < j))
205 {
206 continue ;
207 }
208 Flag [i] = EMPTY ;
209 Ci [nz] = i ;
210 if (values)
211 {
212 Cx [nz] = W [i] + alpha [0] * Ax [p] ;
213 W [i] = 0 ;
214 }
215 nz++ ;
216 }
217
218 /* gather remaining entries into C(:,j), using pattern of B */
219 for (p = pb ; p < pbend ; p++)
220 {
221 i = Bi [p] ;
222 if ((up && i > j) || (lo && i < j))
223 {
224 continue ;
225 }
226 if (Flag [i] == mark)
227 {
228 Ci [nz] = i ;
229 if (values)
230 {
231 Cx [nz] = W [i] ;
232 W [i] = 0 ;
233 }
234 nz++ ;
235 }
236 }
237 }
238
239 Cp [ncol] = nz ;
240
241 /* ---------------------------------------------------------------------- */
242 /* reduce C in size and free temporary matrices */
243 /* ---------------------------------------------------------------------- */
244
245 ASSERT (MAX (1,nz) <= C->nzmax) ;
246 CHOLMOD(reallocate_sparse) (nz, C, Common) ;
247 ASSERT (Common->status >= CHOLMOD_OK) ;
248
249 /* clear the Flag array */
250 mark = CHOLMOD(clear_flag) (Common) ;
251
252 CHOLMOD(free_sparse) (&A2, Common) ;
253 CHOLMOD(free_sparse) (&B2, Common) ;
254
255 /* ---------------------------------------------------------------------- */
256 /* sort C, if requested */
257 /* ---------------------------------------------------------------------- */
258
259 if (sorted)
260 {
261 /* workspace: Iwork (max (nrow,ncol)) */
262 if (!CHOLMOD(sort) (C, Common))
263 {
264 CHOLMOD(free_sparse) (&C, Common) ;
265 if (Common->status < CHOLMOD_OK)
266 {
267 return (NULL) ; /* out of memory */
268 }
269 }
270 }
271
272 /* ---------------------------------------------------------------------- */
273 /* return result */
274 /* ---------------------------------------------------------------------- */
275
276 ASSERT (CHOLMOD(dump_sparse) (C, "add", Common) >= 0) ;
277 return (C) ;
278 }
279