1 //------------------------------------------------------------------------------
2 // GB_emult: C = A.*B, C<M>=A.*B, or C<!M>=A.*B
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7
8 //------------------------------------------------------------------------------
9
10 // GB_emult, does C=A.*B, C<M>=A.*B, C<!M>=A.*B, using the given operator
11 // element-wise on the matrices A and B. The result is typecasted as needed.
12 // The pattern of C is the intersection of the pattern of A and B, intersection
13 // with the mask M or !M, if present.
14
15 // Let the op be z=f(x,y) where x, y, and z have type xtype, ytype, and ztype.
16 // If both A(i,j) and B(i,j) are present, then:
17
18 // C(i,j) = (ctype) op ((xtype) A(i,j), (ytype) B(i,j))
19
20 // If just A(i,j) is present but not B(i,j), or
21 // if just B(i,j) is present but not A(i,j), then C(i,j) is not present.
22
23 // ctype is the type of matrix C, and currently it is always op->ztype.
24
25 // The pattern of C is the intersection of A and B, and also intersection with
26 // M if present and not complemented.
27
28 // TODO: if C is bitmap on input and C_sparsity is GxB_BITMAP, then C=A.*B,
29 // C<M>=A.*B and C<M>+=A.*B can all be done in-place. Also, if C is bitmap
30 // but T<M>=A.*B is sparse (M sparse, with A and B bitmap), then it too can
31 // be done in place.
32
33 #include "GB_emult.h"
34 #include "GB_add.h"
35
36 #define GB_FREE_WORK \
37 { \
38 GB_FREE_WERK (&TaskList, TaskList_size) ; \
39 GB_FREE_WERK (&C_to_M, C_to_M_size) ; \
40 GB_FREE_WERK (&C_to_A, C_to_A_size) ; \
41 GB_FREE_WERK (&C_to_B, C_to_B_size) ; \
42 }
43
44 #define GB_FREE_ALL \
45 { \
46 GB_FREE_WORK ; \
47 GB_phbix_free (C) ; \
48 }
49
GB_emult(GrB_Matrix C,const GrB_Type ctype,const bool C_is_csc,const GrB_Matrix M,const bool Mask_struct,const bool Mask_comp,bool * mask_applied,const GrB_Matrix A,const GrB_Matrix B,const GrB_BinaryOp op,GB_Context Context)50 GrB_Info GB_emult // C=A.*B, C<M>=A.*B, or C<!M>=A.*B
51 (
52 GrB_Matrix C, // output matrix, static header
53 const GrB_Type ctype, // type of output matrix C
54 const bool C_is_csc, // format of output matrix C
55 const GrB_Matrix M, // optional mask, unused if NULL
56 const bool Mask_struct, // if true, use the only structure of M
57 const bool Mask_comp, // if true, use !M
58 bool *mask_applied, // if true, the mask was applied
59 const GrB_Matrix A, // input A matrix
60 const GrB_Matrix B, // input B matrix
61 const GrB_BinaryOp op, // op to perform C = op (A,B)
62 GB_Context Context
63 )
64 {
65
66 //--------------------------------------------------------------------------
67 // check inputs
68 //--------------------------------------------------------------------------
69
70 GrB_Info info ;
71 ASSERT (C != NULL && C->static_header) ;
72
73 ASSERT_MATRIX_OK (A, "A for emult phased", GB0) ;
74 ASSERT_MATRIX_OK (B, "B for emult phased", GB0) ;
75 ASSERT_MATRIX_OK_OR_NULL (M, "M for emult phased", GB0) ;
76 ASSERT_BINARYOP_OK (op, "op for emult phased", GB0) ;
77 ASSERT (A->vdim == B->vdim && A->vlen == B->vlen) ;
78 ASSERT (GB_IMPLIES (M != NULL, A->vdim == M->vdim && A->vlen == M->vlen)) ;
79
80 //--------------------------------------------------------------------------
81 // declare workspace
82 //--------------------------------------------------------------------------
83
84 GB_task_struct *TaskList = NULL ; size_t TaskList_size = 0 ;
85 int64_t *restrict C_to_M = NULL ; size_t C_to_M_size = 0 ;
86 int64_t *restrict C_to_A = NULL ; size_t C_to_A_size = 0 ;
87 int64_t *restrict C_to_B = NULL ; size_t C_to_B_size = 0 ;
88
89 //--------------------------------------------------------------------------
90 // delete any lingering zombies and assemble any pending tuples
91 //--------------------------------------------------------------------------
92
93 // some cases can allow M, A, and/or B to be jumbled
94 GB_MATRIX_WAIT_IF_PENDING_OR_ZOMBIES (M) ;
95 GB_MATRIX_WAIT_IF_PENDING_OR_ZOMBIES (A) ;
96 GB_MATRIX_WAIT_IF_PENDING_OR_ZOMBIES (B) ;
97
98 //--------------------------------------------------------------------------
99 // determine the sparsity of C and the method to use
100 //--------------------------------------------------------------------------
101
102 bool apply_mask ; // if true, mask is applied during emult
103 int ewise_method ; // method to use
104 int C_sparsity = GB_emult_sparsity (&apply_mask, &ewise_method,
105 M, Mask_comp, A, B) ;
106
107 //--------------------------------------------------------------------------
108 // C<M or !M> = A.*B
109 //--------------------------------------------------------------------------
110
111 switch (ewise_method)
112 {
113
114 case GB_EMULT_METHOD_ADD : // A and B both full (or as-if-full)
115
116 // ------------------------------------------
117 // C = A .* B
118 // ------------------------------------------
119 // full . full full (GB_add)
120 // ------------------------------------------
121 // C <M> = A .* B
122 // ------------------------------------------
123 // sparse sparse full full (GB_add or 03)
124 // bitmap bitmap full full (GB_add or 07)
125 // bitmap full full full (GB_add or 07)
126 // ------------------------------------------
127 // C <!M>= A .* B
128 // ------------------------------------------
129 // bitmap sparse full full (GB_add or 06)
130 // bitmap bitmap full full (GB_add or 07)
131 // bitmap full full full (GB_add or 07)
132
133 // A and B are both full (or as-if-full). The mask M may be
134 // anything. GB_add computes the same thing in this case, so it is
135 // used instead, to reduce the code needed for GB_emult. GB_add
136 // must be used for C=A.*B if all 3 matrices are full. Otherwise,
137 // GB_emult method can be used as well.
138
139 return (GB_add (C, ctype, C_is_csc, M, Mask_struct,
140 Mask_comp, mask_applied, A, B, op, Context)) ;
141
142 case GB_EMULT_METHOD_02A : // A sparse/hyper, B bitmap/full
143
144 // ------------------------------------------
145 // C = A .* B
146 // ------------------------------------------
147 // sparse . sparse bitmap (method: 02a)
148 // sparse . sparse full (method: 02a)
149 // ------------------------------------------
150 // C <M> = A .* B
151 // ------------------------------------------
152 // sparse bitmap sparse bitmap (method: 02a)
153 // sparse bitmap sparse full (method: 02a)
154 // sparse full sparse bitmap (method: 02a)
155 // sparse full sparse full (method: 02a)
156 // ------------------------------------------
157 // C <!M>= A .* B
158 // ------------------------------------------
159 // sparse sparse sparse bitmap (02a: M later)
160 // sparse sparse sparse full (02a: M later)
161 // ------------------------------------------
162 // C <!M> = A .* B
163 // ------------------------------------------
164 // sparse bitmap sparse bitmap (method: 02a)
165 // sparse bitmap sparse full (method: 02a)
166 // sparse full sparse bitmap (method: 02a)
167 // sparse full sparse full (method: 02a)
168
169 // A is sparse/hyper and B is bitmap/full. M is either not
170 // present, not applied (!M when sparse/hyper), or bitmap/full.
171 // This method does not handle the case when M is sparse/hyper,
172 // unless M is ignored and applied later.
173
174 return (GB_emult_02 (C, ctype, C_is_csc,
175 (apply_mask) ? M : NULL, Mask_struct, Mask_comp,
176 A, B, op, false, Context)) ;
177
178 case GB_EMULT_METHOD_02B : // A bitmap/full, B sparse/hyper
179
180 // ------------------------------------------
181 // C = A .* B
182 // ------------------------------------------
183 // sparse . bitmap sparse (method: 02b)
184 // sparse . full sparse (method: 02b)
185 // ------------------------------------------
186 // C <M> = A .* B
187 // ------------------------------------------
188 // sparse bitmap bitmap sparse (method: 02b)
189 // sparse bitmap full sparse (method: 02b)
190 // sparse full bitmap sparse (method: 02b)
191 // sparse full full sparse (method: 02b)
192 // ------------------------------------------
193 // C <!M>= A .* B
194 // ------------------------------------------
195 // sparse sparse bitmap sparse (02b: M later)
196 // sparse sparse full sparse (02b: M later)
197 // ------------------------------------------
198 // C <!M> = A .* B
199 // ------------------------------------------
200 // sparse bitmap bitmap sparse (method: 02b)
201 // sparse bitmap full sparse (method: 02b)
202 // sparse full bitmap sparse (method: 02b)
203 // sparse full full sparse (method: 02b)
204
205 // A is bitmap/full and B is sparse/hyper, with flipxy true.
206 // M is not present, not applied, or bitmap/full
207 // Note that A and B are flipped.
208
209 return (GB_emult_02 (C, ctype, C_is_csc,
210 (apply_mask) ? M : NULL, Mask_struct, Mask_comp,
211 B, A, op, true, Context)) ;
212
213 case GB_EMULT_METHOD_01 :
214
215 // ------------------------------------------
216 // C = A .* B
217 // ------------------------------------------
218 // sparse . sparse sparse (method: 01)
219 // ------------------------------------------
220 // C <M> = A .* B
221 // ------------------------------------------
222 // sparse sparse sparse sparse (method: 01)
223 // sparse bitmap sparse sparse (method: 01)
224 // sparse full sparse sparse (method: 01)
225 // ------------------------------------------
226 // C <!M>= A .* B
227 // ------------------------------------------
228 // sparse sparse sparse sparse (01: M later)
229 // sparse bitmap sparse sparse (method: 01)
230 // sparse full sparse sparse (method: 01)
231
232 // TODO: break method 01 into different methods
233 break ;
234
235 case GB_EMULT_METHOD_05 : // C is bitmap, M is not present
236
237 // ------------------------------------------
238 // C = A .* B
239 // ------------------------------------------
240 // bitmap . bitmap bitmap (method: 05)
241 // bitmap . bitmap full (method: 05)
242 // bitmap . full bitmap (method: 05)
243
244 case GB_EMULT_METHOD_06 : // C is bitmap, !M is sparse/hyper
245
246 // ------------------------------------------
247 // C <!M>= A .* B
248 // ------------------------------------------
249 // bitmap sparse bitmap bitmap (method: 06)
250 // bitmap sparse bitmap full (method: 06)
251 // bitmap sparse full bitmap (method: 06)
252 // bitmap sparse full full (GB_add or 06)
253
254 case GB_EMULT_METHOD_07 : // C is bitmap, M is bitmap/full
255
256 // ------------------------------------------
257 // C <M> = A .* B
258 // ------------------------------------------
259 // bitmap bitmap bitmap bitmap (method: 07)
260 // bitmap bitmap bitmap full (method: 07)
261 // bitmap bitmap full bitmap (method: 07)
262 // bitmap bitmap full full (GB_add or 07)
263 // bitmap full bitmap bitmap (method: 07)
264 // bitmap full bitmap full (method: 07)
265 // bitmap full full bitmap (method: 07)
266 // bitmap full full full (GB_add or 07)
267 // ------------------------------------------
268 // C <!M> = A .* B
269 // ------------------------------------------
270 // bitmap bitmap bitmap bitmap (method: 07)
271 // bitmap bitmap bitmap full (method: 07)
272 // bitmap bitmap full bitmap (method: 07)
273 // bitmap bitmap full full (GB_add or 07)
274 // bitmap full bitmap bitmap (method: 07)
275 // bitmap full bitmap full (method: 07)
276 // bitmap full full bitmap (method: 07)
277 // bitmap full full full (GB_add or 07)
278
279 // For methods 05, 06, and 07, C is constructed as bitmap.
280 // Both A and B are bitmap/full. M is either not present,
281 // complemented, or not complemented and bitmap/full. The
282 // case when M is not complemented and sparse/hyper is handled
283 // by method 03, which constructs C as sparse/hyper (the same
284 // structure as M), not bitmap.
285
286 return (GB_bitmap_emult (C, ewise_method, ctype, C_is_csc,
287 M, Mask_struct, Mask_comp, mask_applied, A, B,
288 op, Context)) ;
289
290 case GB_EMULT_METHOD_03 :
291
292 // ------------------------------------------
293 // C <M>= A .* B
294 // ------------------------------------------
295 // sparse sparse bitmap bitmap (method: 03)
296 // sparse sparse bitmap full (method: 03)
297 // sparse sparse full bitmap (method: 03)
298 // sparse sparse full full (GB_add or 03)
299
300 return (GB_emult_03 (C, ctype, C_is_csc, M, Mask_struct,
301 mask_applied, A, B, op, Context)) ;
302
303 case GB_EMULT_METHOD_04A : break ; // punt
304
305 // ------------------------------------------
306 // C <M>= A .* B
307 // ------------------------------------------
308 // sparse sparse sparse bitmap (method: 04a)
309 // sparse sparse sparse full (method: 04a)
310
311 // TODO: this will use 04 (M,A,B, flipxy=false)
312
313 // The method will compute the 2-way intersection of M and A,
314 // using the same parallization as C=A.*B when both A and B are
315 // both sparse. It will then lookup B in O(1) time.
316 // M and A must not be jumbled.
317
318 case GB_EMULT_METHOD_04B : break ; // punt
319
320 // ------------------------------------------
321 // C <M>= A .* B
322 // ------------------------------------------
323 // sparse sparse bitmap sparse (method: 04b)
324 // sparse sparse full sparse (method: 04b)
325
326 // TODO: this will use 04 (M,B,A, flipxy=true)
327 // M and B must not be jumbled.
328
329 default:;
330 }
331
332 //--------------------------------------------------------------------------
333 // method 01 (and for now, 04a and 04b)
334 //--------------------------------------------------------------------------
335
336 ASSERT (C_sparsity == GxB_SPARSE || C_sparsity == GxB_HYPERSPARSE) ;
337
338 GB_MATRIX_WAIT (M) ;
339 GB_MATRIX_WAIT (A) ;
340 GB_MATRIX_WAIT (B) ;
341
342 GBURBLE ("emult:(%s<%s>=%s.*%s) ",
343 GB_sparsity_char (C_sparsity),
344 GB_sparsity_char_matrix (M),
345 GB_sparsity_char_matrix (A),
346 GB_sparsity_char_matrix (B)) ;
347
348 //--------------------------------------------------------------------------
349 // initializations
350 //--------------------------------------------------------------------------
351
352 int64_t Cnvec, Cnvec_nonempty ;
353 int64_t *Cp = NULL ; size_t Cp_size = 0 ;
354 const int64_t *Ch = NULL ; size_t Ch_size = 0 ;
355 int C_ntasks = 0, C_nthreads ;
356
357 //--------------------------------------------------------------------------
358 // phase0: finalize the sparsity C and find the vectors in C
359 //--------------------------------------------------------------------------
360
361 GB_OK (GB_emult_01_phase0 (
362 // computed by phase0:
363 &Cnvec, &Ch, &Ch_size, &C_to_M, &C_to_M_size, &C_to_A, &C_to_A_size,
364 &C_to_B, &C_to_B_size,
365 // input/output to phase0:
366 &C_sparsity,
367 // original input:
368 (apply_mask) ? M : NULL, A, B, Context)) ;
369
370 // C is still sparse or hypersparse, not bitmap or full
371 ASSERT (C_sparsity == GxB_SPARSE || C_sparsity == GxB_HYPERSPARSE) ;
372
373 //--------------------------------------------------------------------------
374 // phase1: split C into tasks, and count entries in each vector of C
375 //--------------------------------------------------------------------------
376
377 // phase1a: split C into tasks
378 GB_OK (GB_ewise_slice (
379 // computed by phase1a:
380 &TaskList, &TaskList_size, &C_ntasks, &C_nthreads,
381 // computed by phase0:
382 Cnvec, Ch, C_to_M, C_to_A, C_to_B, false,
383 // original input:
384 (apply_mask) ? M : NULL, A, B, Context)) ;
385
386 // count the number of entries in each vector of C
387 GB_OK (GB_emult_01_phase1 (
388 // computed by phase1:
389 &Cp, &Cp_size, &Cnvec_nonempty,
390 // from phase1a:
391 TaskList, C_ntasks, C_nthreads,
392 // from phase0:
393 Cnvec, Ch, C_to_M, C_to_A, C_to_B,
394 // original input:
395 (apply_mask) ? M : NULL, Mask_struct, Mask_comp, A, B, Context)) ;
396
397 //--------------------------------------------------------------------------
398 // phase2: compute the entries (indices and values) in each vector of C
399 //--------------------------------------------------------------------------
400
401 // Cp is either freed by phase2, or transplanted into C.
402 // Either way, it is not freed here.
403
404 GB_OK (GB_emult_01_phase2 (
405 // computed or used by phase2:
406 C, ctype, C_is_csc, op,
407 // from phase1:
408 &Cp, Cp_size, Cnvec_nonempty,
409 // from phase1a:
410 TaskList, C_ntasks, C_nthreads,
411 // from phase0:
412 Cnvec, Ch, Ch_size, C_to_M, C_to_A, C_to_B, C_sparsity,
413 // from GB_emult_sparsity:
414 ewise_method,
415 // original input:
416 (apply_mask) ? M : NULL, Mask_struct, Mask_comp, A, B, Context)) ;
417
418 //--------------------------------------------------------------------------
419 // free workspace and return result
420 //--------------------------------------------------------------------------
421
422 GB_FREE_WORK ;
423 ASSERT_MATRIX_OK (C, "C output for emult phased", GB0) ;
424 (*mask_applied) = apply_mask ;
425 return (GrB_SUCCESS) ;
426 }
427
428