1 //------------------------------------------------------------------------------
2 // GB_emult_01_phase0: find vectors of C to compute for C=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 // The eWise multiply of two matrices, C=A.*B, C<M>=A.*B, or C<!M>=A.*B starts
11 // with this phase, which determines which vectors of C need to be computed.
12
13 // On input, A and B are the two matrices being ewise multiplied, and M is the
14 // optional mask matrix. If present, it is not complemented.
15
16 // The M, A, and B matrices are sparse or hypersparse. C will be sparse
17 // (if Ch is returned NULL) or hypersparse (if Ch is returned non-NULL).
18
19 // Ch: the vectors to compute in C. Not allocated, but equal to either
20 // A->h, B->h, or M->h, or NULL if C is not hypersparse.
21
22 // C_to_A: if A is hypersparse, and Ch is not A->h, then C_to_A [k] = kA
23 // if the kth vector j = Ch [k] is equal to Ah [kA]. If j does not appear
24 // in A, then C_to_A [k] = -1. Otherwise, C_to_A is returned as NULL.
25 // C is always hypersparse in this case.
26
27 // C_to_B: if B is hypersparse, and Ch is not B->h, then C_to_B [k] = kB
28 // if the kth vector j = Ch [k] is equal to Bh [kB]. If j does not appear
29 // in B, then C_to_B [k] = -1. Otherwise, C_to_B is returned as NULL.
30 // C is always hypersparse in this case.
31
32 // C_to_M: if M is hypersparse, and Ch is not M->h, then C_to_M [k] = kM
33 // if the kth vector j = GBH (Ch, k) is equal to Mh [kM].
34 // If j does not appear in M, then C_to_M [k] = -1. Otherwise, C_to_M is
35 // returned as NULL. C is always hypersparse in this case.
36
37 // FUTURE:: exploit A==M, B==M, and A==B aliases
38
39 #include "GB_emult.h"
40
GB_emult_01_phase0(int64_t * p_Cnvec,const int64_t * restrict * Ch_handle,size_t * Ch_size_handle,int64_t * restrict * C_to_M_handle,size_t * C_to_M_size_handle,int64_t * restrict * C_to_A_handle,size_t * C_to_A_size_handle,int64_t * restrict * C_to_B_handle,size_t * C_to_B_size_handle,int * C_sparsity,const GrB_Matrix M,const GrB_Matrix A,const GrB_Matrix B,GB_Context Context)41 GrB_Info GB_emult_01_phase0 // find vectors in C for C=A.*B or C<M>=A.*B
42 (
43 int64_t *p_Cnvec, // # of vectors to compute in C
44 const int64_t *restrict *Ch_handle, // Ch is M->h, A->h, B->h, or NULL
45 size_t *Ch_size_handle,
46 int64_t *restrict *C_to_M_handle, // C_to_M: size Cnvec, or NULL
47 size_t *C_to_M_size_handle,
48 int64_t *restrict *C_to_A_handle, // C_to_A: size Cnvec, or NULL
49 size_t *C_to_A_size_handle,
50 int64_t *restrict *C_to_B_handle, // C_to_B: size Cnvec, or NULL
51 size_t *C_to_B_size_handle,
52 int *C_sparsity, // sparsity structure of C
53 // original input:
54 const GrB_Matrix M, // optional mask, may be NULL
55 const GrB_Matrix A,
56 const GrB_Matrix B,
57 GB_Context Context
58 )
59 {
60
61 //--------------------------------------------------------------------------
62 // check inputs
63 //--------------------------------------------------------------------------
64
65 // M, A, and B can be jumbled for this phase
66
67 ASSERT (p_Cnvec != NULL) ;
68 ASSERT (Ch_handle != NULL) ;
69 ASSERT (Ch_size_handle != NULL) ;
70 ASSERT (C_to_A_handle != NULL) ;
71 ASSERT (C_to_B_handle != NULL) ;
72
73 ASSERT_MATRIX_OK_OR_NULL (M, "M for emult phase0", GB0) ;
74 ASSERT (!GB_ZOMBIES (M)) ;
75 ASSERT (GB_JUMBLED_OK (M)) ; // pattern not accessed
76 ASSERT (!GB_PENDING (M)) ;
77
78 ASSERT_MATRIX_OK (A, "A for emult phase0", GB0) ;
79 ASSERT (!GB_ZOMBIES (A)) ;
80 ASSERT (GB_JUMBLED_OK (B)) ; // pattern not accessed
81 ASSERT (!GB_PENDING (A)) ;
82
83 ASSERT_MATRIX_OK (B, "B for emult phase0", GB0) ;
84 ASSERT (!GB_ZOMBIES (B)) ;
85 ASSERT (GB_JUMBLED_OK (A)) ; // pattern not accessed
86 ASSERT (!GB_PENDING (B)) ;
87
88 ASSERT (A->vdim == B->vdim) ;
89 ASSERT (A->vlen == B->vlen) ;
90 ASSERT (GB_IMPLIES (M != NULL, A->vdim == M->vdim)) ;
91 ASSERT (GB_IMPLIES (M != NULL, A->vlen == M->vlen)) ;
92
93 //--------------------------------------------------------------------------
94 // initializations
95 //--------------------------------------------------------------------------
96
97 (*p_Cnvec) = 0 ;
98 (*Ch_handle) = NULL ;
99 (*Ch_size_handle) = 0 ;
100 if (C_to_M_handle != NULL)
101 {
102 (*C_to_M_handle) = NULL ;
103 }
104 (*C_to_A_handle) = NULL ;
105 (*C_to_B_handle) = NULL ;
106
107 ASSERT ((*C_sparsity) == GxB_SPARSE || (*C_sparsity) == GxB_HYPERSPARSE) ;
108
109 const int64_t *restrict Ch = NULL ; size_t Ch_size = 0 ;
110 int64_t *restrict C_to_M = NULL ; size_t C_to_M_size = 0 ;
111 int64_t *restrict C_to_A = NULL ; size_t C_to_A_size = 0 ;
112 int64_t *restrict C_to_B = NULL ; size_t C_to_B_size = 0 ;
113
114 //--------------------------------------------------------------------------
115 // get content of M, A, and B
116 //--------------------------------------------------------------------------
117
118 int64_t n = A->vdim ;
119
120 int64_t Anvec = A->nvec ;
121 int64_t vlen = A->vlen ;
122 const int64_t *restrict Ah = A->h ;
123 bool A_is_hyper = (Ah != NULL) ;
124
125 int64_t Bnvec = B->nvec ;
126 const int64_t *restrict Bh = B->h ;
127 bool B_is_hyper = (Bh != NULL) ;
128
129 int64_t Mnvec = 0 ;
130 const int64_t *restrict Mh = NULL ;
131 bool M_is_hyper = false ;
132
133 if (M != NULL)
134 {
135 Mnvec = M->nvec ;
136 Mh = M->h ;
137 M_is_hyper = (Mh != NULL) ;
138 }
139
140 //--------------------------------------------------------------------------
141 // determine how to construct the vectors of C
142 //--------------------------------------------------------------------------
143
144 if (M != NULL)
145 {
146
147 //----------------------------------------------------------------------
148 // 8 cases to consider: A, B, M can each be hyper or sparse
149 //----------------------------------------------------------------------
150
151 // Mask is present and not complemented
152
153 if (A_is_hyper)
154 {
155 if (B_is_hyper)
156 {
157 if (M_is_hyper)
158 {
159
160 //----------------------------------------------------------
161 // (1) A hyper, B hyper, M hyper: C hyper
162 //----------------------------------------------------------
163
164 // Ch = smaller of Mh, Bh, Ah
165
166 int64_t nvec = GB_IMIN (Anvec, Bnvec) ;
167 nvec = GB_IMIN (nvec, Mnvec) ;
168 if (nvec == Anvec)
169 {
170 Ch = Ah ; Ch_size = A->h_size ;
171 }
172 else if (nvec == Bnvec)
173 {
174 Ch = Bh ; Ch_size = B->h_size ;
175 }
176 else // (nvec == Mnvec)
177 {
178 Ch = Mh ; Ch_size = M->h_size ;
179 }
180
181 }
182 else
183 {
184
185 //----------------------------------------------------------
186 // (2) A hyper, B hyper, M sparse: C hyper
187 //----------------------------------------------------------
188
189 // Ch = smaller of Ah, Bh
190 if (Anvec <= Bnvec)
191 {
192 Ch = Ah ; Ch_size = A->h_size ;
193 }
194 else
195 {
196 Ch = Bh ; Ch_size = B->h_size ;
197 }
198 }
199
200 }
201 else
202 {
203
204 if (M_is_hyper)
205 {
206
207 //----------------------------------------------------------
208 // (3) A hyper, B sparse, M hyper: C hyper
209 //----------------------------------------------------------
210
211 // Ch = smaller of Mh, Ah
212 if (Anvec <= Mnvec)
213 {
214 Ch = Ah ; Ch_size = A->h_size ;
215 }
216 else
217 {
218 Ch = Mh ; Ch_size = M->h_size ;
219 }
220
221 }
222 else
223 {
224
225 //----------------------------------------------------------
226 // (4) A hyper, B sparse, M sparse: C hyper
227 //----------------------------------------------------------
228
229 Ch = Ah ; Ch_size = A->h_size ;
230 }
231 }
232
233 }
234 else
235 {
236
237 if (B_is_hyper)
238 {
239 if (M_is_hyper)
240 {
241
242 //----------------------------------------------------------
243 // (5) A sparse, B hyper, M hyper: C hyper
244 //----------------------------------------------------------
245
246 // Ch = smaller of Mh, Bh
247
248 if (Bnvec <= Mnvec)
249 {
250 Ch = Bh ; Ch_size = B->h_size ;
251 }
252 else
253 {
254 Ch = Mh ; Ch_size = M->h_size ;
255 }
256
257 }
258 else
259 {
260
261 //----------------------------------------------------------
262 // (6) A sparse, B hyper, M sparse: C hyper
263 //----------------------------------------------------------
264
265 Ch = Bh ; Ch_size = B->h_size ;
266
267 }
268 }
269 else
270 {
271
272 if (M_is_hyper)
273 {
274
275 //----------------------------------------------------------
276 // (7) A sparse, B sparse, M hyper: C hyper
277 //----------------------------------------------------------
278
279 Ch = Mh ; Ch_size = M->h_size ;
280
281 }
282 else
283 {
284
285 //----------------------------------------------------------
286 // (8) A sparse, B sparse, M sparse: C sparse
287 //----------------------------------------------------------
288
289 Ch = NULL ;
290 }
291 }
292 }
293
294 }
295 else
296 {
297
298 //----------------------------------------------------------------------
299 // 4 cases to consider: A, B can be hyper or sparse
300 //----------------------------------------------------------------------
301
302 // Mask is not present, or present and complemented.
303
304 if (A_is_hyper)
305 {
306 if (B_is_hyper)
307 {
308
309 //--------------------------------------------------------------
310 // (1) A hyper, B hyper: C hyper
311 //--------------------------------------------------------------
312
313 // Ch = smaller of Ah, Bh
314 if (Anvec <= Bnvec)
315 {
316 Ch = Ah ; Ch_size = A->h_size ;
317 }
318 else
319 {
320 Ch = Bh ; Ch_size = B->h_size ;
321 }
322 }
323 else
324 {
325
326 //--------------------------------------------------------------
327 // (2) A hyper, B sparse: C hyper
328 //--------------------------------------------------------------
329
330 Ch = Ah ; Ch_size = A->h_size ;
331
332 }
333
334 }
335 else
336 {
337
338 if (B_is_hyper)
339 {
340
341 //--------------------------------------------------------------
342 // (3) A sparse, B hyper: C hyper
343 //--------------------------------------------------------------
344
345 Ch = Bh ; Ch_size = B->h_size ;
346
347 }
348 else
349 {
350
351 //--------------------------------------------------------------
352 // (4) A sparse, B sparse: C sparse
353 //--------------------------------------------------------------
354
355 Ch = NULL ;
356 }
357 }
358 }
359
360 //--------------------------------------------------------------------------
361 // find Cnvec
362 //--------------------------------------------------------------------------
363
364 int64_t Cnvec ;
365
366 if (Ch == NULL)
367 {
368 // C is sparse
369 (*C_sparsity) = GxB_SPARSE ;
370 Cnvec = n ;
371 }
372 else
373 {
374 // C is hypersparse; one of A, B, or M are hypersparse
375 ASSERT (A_is_hyper || B_is_hyper || M_is_hyper) ;
376 (*C_sparsity) = GxB_HYPERSPARSE ;
377 if (Ch == Ah)
378 {
379 Cnvec = Anvec ;
380 }
381 else if (Ch == Bh)
382 {
383 Cnvec = Bnvec ;
384 }
385 else // (Ch == Mh)
386 {
387 Cnvec = Mnvec ;
388 }
389 }
390
391 //--------------------------------------------------------------------------
392 // determine the number of threads to use
393 //--------------------------------------------------------------------------
394
395 GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
396 int nthreads = GB_nthreads (Cnvec, chunk, nthreads_max) ;
397
398 //--------------------------------------------------------------------------
399 // construct C_to_M mapping
400 //--------------------------------------------------------------------------
401
402 if (M_is_hyper && Ch != Mh)
403 {
404 // allocate C_to_M
405 C_to_M = GB_MALLOC_WERK (Cnvec, int64_t, &C_to_M_size) ;
406 if (C_to_M == NULL)
407 {
408 // out of memory
409 return (GrB_OUT_OF_MEMORY) ;
410 }
411
412 // compute C_to_M
413 ASSERT (Ch != NULL) ;
414
415 const int64_t *restrict Mp = M->p ;
416
417 int64_t k ;
418 #pragma omp parallel for num_threads(nthreads) schedule(static)
419 for (k = 0 ; k < Cnvec ; k++)
420 {
421 int64_t pM, pM_end, kM = 0 ;
422 int64_t j = Ch [k] ;
423 GB_lookup (true, Mh, Mp, vlen, &kM, Mnvec-1, j, &pM, &pM_end) ;
424 C_to_M [k] = (pM < pM_end) ? kM : -1 ;
425 }
426 }
427
428 //--------------------------------------------------------------------------
429 // construct C_to_A mapping
430 //--------------------------------------------------------------------------
431
432 if (A_is_hyper && Ch != Ah)
433 {
434 // allocate C_to_A
435 C_to_A = GB_MALLOC_WERK (Cnvec, int64_t, &C_to_A_size) ;
436 if (C_to_A == NULL)
437 {
438 // out of memory
439 GB_FREE_WERK (&C_to_M, C_to_M_size) ;
440 return (GrB_OUT_OF_MEMORY) ;
441 }
442
443 // compute C_to_A
444 ASSERT (Ch != NULL) ;
445 const int64_t *restrict Ap = A->p ;
446
447 int64_t k ;
448 #pragma omp parallel for num_threads(nthreads) schedule(static)
449 for (k = 0 ; k < Cnvec ; k++)
450 {
451 int64_t pA, pA_end, kA = 0 ;
452 int64_t j = Ch [k] ;
453 GB_lookup (true, Ah, Ap, vlen, &kA, Anvec-1, j, &pA, &pA_end) ;
454 C_to_A [k] = (pA < pA_end) ? kA : -1 ;
455 }
456 }
457
458 //--------------------------------------------------------------------------
459 // construct C_to_B mapping
460 //--------------------------------------------------------------------------
461
462 if (B_is_hyper && Ch != Bh)
463 {
464 // allocate C_to_B
465 C_to_B = GB_MALLOC_WERK (Cnvec, int64_t, &C_to_B_size) ;
466 if (C_to_B == NULL)
467 {
468 // out of memory
469 GB_FREE_WERK (&C_to_M, C_to_M_size) ;
470 GB_FREE_WERK (&C_to_A, C_to_A_size) ;
471 return (GrB_OUT_OF_MEMORY) ;
472 }
473
474 // compute C_to_B
475 ASSERT (Ch != NULL) ;
476 const int64_t *restrict Bp = B->p ;
477
478 int64_t k ;
479 #pragma omp parallel for num_threads(nthreads) schedule(static)
480 for (k = 0 ; k < Cnvec ; k++)
481 {
482 int64_t pB, pB_end, kB = 0 ;
483 int64_t j = Ch [k] ;
484 GB_lookup (true, Bh, Bp, vlen, &kB, Bnvec-1, j, &pB, &pB_end) ;
485 C_to_B [k] = (pB < pB_end) ? kB : -1 ;
486 }
487 }
488
489 //--------------------------------------------------------------------------
490 // return result
491 //--------------------------------------------------------------------------
492
493 (*p_Cnvec) = Cnvec ;
494 (*Ch_handle) = Ch ;
495 (*Ch_size_handle) = Ch_size ;
496 if (C_to_M_handle != NULL)
497 {
498 (*C_to_M_handle) = C_to_M ;
499 (*C_to_M_size_handle) = C_to_M_size ;
500 }
501 (*C_to_A_handle) = C_to_A ; (*C_to_A_size_handle) = C_to_A_size ;
502 (*C_to_B_handle) = C_to_B ; (*C_to_B_size_handle) = C_to_B_size ;
503
504 //--------------------------------------------------------------------------
505 // The code below describes what the output contains:
506 //--------------------------------------------------------------------------
507
508 #ifdef GB_DEBUG
509 ASSERT (A != NULL) ; // A and B are always present
510 ASSERT (B != NULL) ;
511 int64_t jlast = -1 ;
512 for (int64_t k = 0 ; k < Cnvec ; k++)
513 {
514
515 // C(:,j) is in the list, as the kth vector
516 int64_t j ;
517 if (Ch == NULL)
518 {
519 // C will be constructed as sparse
520 j = k ;
521 }
522 else
523 {
524 // C will be constructed as hypersparse
525 j = Ch [k] ;
526 }
527
528 // vectors j in Ch are sorted, and in the range 0:n-1
529 ASSERT (j >= 0 && j < n) ;
530 ASSERT (j > jlast) ;
531 jlast = j ;
532
533 // see if A (:,j) exists
534 if (C_to_A != NULL)
535 {
536 // A is hypersparse
537 ASSERT (A_is_hyper)
538 int64_t kA = C_to_A [k] ;
539 ASSERT (kA >= -1 && kA < A->nvec) ;
540 if (kA >= 0)
541 {
542 int64_t jA = A->h [kA] ;
543 ASSERT (j == jA) ;
544 }
545 }
546 else if (A_is_hyper)
547 {
548 // A is hypersparse, and Ch is a shallow copy of A->h
549 ASSERT (Ch == A->h) ;
550 }
551
552 // see if B (:,j) exists
553 if (C_to_B != NULL)
554 {
555 // B is hypersparse
556 ASSERT (B_is_hyper)
557 int64_t kB = C_to_B [k] ;
558 ASSERT (kB >= -1 && kB < B->nvec) ;
559 if (kB >= 0)
560 {
561 int64_t jB = B->h [kB] ;
562 ASSERT (j == jB) ;
563 }
564 }
565 else if (B_is_hyper)
566 {
567 // A is hypersparse, and Ch is a shallow copy of A->h
568 ASSERT (Ch == B->h) ;
569 }
570
571 // see if M (:,j) exists
572 if (Ch != NULL && M != NULL && Ch == M->h)
573 {
574 // Ch is the same as Mh
575 ASSERT (M != NULL) ;
576 ASSERT (M->h != NULL) ;
577 ASSERT (Ch != NULL && M->h != NULL && Ch [k] == M->h [k]) ;
578 ASSERT (C_to_M == NULL) ;
579 }
580 else if (C_to_M != NULL)
581 {
582 // M is present and hypersparse
583 ASSERT (M != NULL) ;
584 ASSERT (M->h != NULL) ;
585 int64_t kM = C_to_M [k] ;
586 ASSERT (kM >= -1 && kM < M->nvec) ;
587 if (kM >= 0)
588 {
589 int64_t jM = M->h [kM] ;
590 ASSERT (j == jM) ;
591 }
592 }
593 else
594 {
595 // M is not present, or in sparse form
596 ASSERT (M == NULL || M->h == NULL) ;
597 }
598 }
599
600 #endif
601
602 return (GrB_SUCCESS) ;
603 }
604
605