1 //------------------------------------------------------------------------------
2 // GB_bitmap_emult_template: C = A.*B, C<M>=A.*B, and C<!M>=A.*B, C bitmap
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 // C is bitmap.  A and B are bitmap or full.  M depends on the method
11 
12 {
13 
14     //--------------------------------------------------------------------------
15     // get C, A, and B
16     //--------------------------------------------------------------------------
17 
18     const int8_t  *restrict Ab = A->b ;
19     const int8_t  *restrict Bb = B->b ;
20     const int64_t vlen = A->vlen ;
21 
22     ASSERT (GB_IS_BITMAP (A) || GB_IS_FULL (A) || GB_as_if_full (A)) ;
23     ASSERT (GB_IS_BITMAP (B) || GB_IS_FULL (A) || GB_as_if_full (B)) ;
24 
25     const GB_ATYPE *restrict Ax = (GB_ATYPE *) A->x ;
26     const GB_BTYPE *restrict Bx = (GB_BTYPE *) B->x ;
27           int8_t   *restrict Cb = C->b ;
28           GB_CTYPE *restrict Cx = (GB_CTYPE *) C->x ;
29     const int64_t cnz = GB_NNZ_HELD (C) ;
30 
31     //--------------------------------------------------------------------------
32     // C=A.*B, C<M>=A.*B, or C<!M>=A.*B: C is bitmap
33     //--------------------------------------------------------------------------
34 
35     // TODO modify this method so it can modify C in-place, and also use the
36     // accum operator.
37     int64_t cnvals = 0 ;
38 
39     if (ewise_method == GB_EMULT_METHOD_05)
40     {
41 
42         //----------------------------------------------------------------------
43         // C is bitmap, M is not present
44         //----------------------------------------------------------------------
45 
46         //      ------------------------------------------
47         //      C       =           A       .*      B
48         //      ------------------------------------------
49         //      bitmap  .           bitmap          bitmap  (method: 05)
50         //      bitmap  .           bitmap          full    (method: 05)
51         //      bitmap  .           full            bitmap  (method: 05)
52 
53         int tid ;
54         #pragma omp parallel for num_threads(C_nthreads) schedule(static) \
55             reduction(+:cnvals)
56         for (tid = 0 ; tid < C_nthreads ; tid++)
57         {
58             int64_t pstart, pend, task_cnvals = 0 ;
59             GB_PARTITION (pstart, pend, cnz, tid, C_nthreads) ;
60             for (int64_t p = pstart ; p < pend ; p++)
61             {
62                 if (GBB (Ab, p) && GBB (Bb,p))
63                 {
64                     // C (i,j) = A (i,j) + B (i,j)
65                     GB_GETA (aij, Ax, p) ;
66                     GB_GETB (bij, Bx, p) ;
67                     GB_BINOP (GB_CX (p), aij, bij, p % vlen, p / vlen) ;
68                     Cb [p] = 1 ;
69                     task_cnvals++ ;
70                 }
71             }
72             cnvals += task_cnvals ;
73         }
74 
75     }
76     else if (ewise_method == GB_EMULT_METHOD_06)
77     {
78 
79         //----------------------------------------------------------------------
80         // C is bitmap, !M is sparse or hyper
81         //----------------------------------------------------------------------
82 
83         //      ------------------------------------------
84         //      C       <!M>=       A       .*      B
85         //      ------------------------------------------
86         //      bitmap  sparse      bitmap          bitmap  (method: 06)
87         //      bitmap  sparse      bitmap          full    (method: 06)
88         //      bitmap  sparse      full            bitmap  (method: 06)
89 
90         // M is sparse and complemented.  If M is sparse and not
91         // complemented, then C is constructed as sparse, not bitmap.
92         ASSERT (M != NULL) ;
93         ASSERT (Mask_comp) ;
94         ASSERT (GB_IS_SPARSE (M) || GB_IS_HYPERSPARSE (M)) ;
95 
96         // C(i,j) = A(i,j) .* B(i,j) can only be computed where M(i,j) is
97         // not present in the sparse pattern of M, and where it is present
98         // but equal to zero.
99 
100         //----------------------------------------------------------------------
101         // scatter M into the C bitmap
102         //----------------------------------------------------------------------
103 
104         GB_bitmap_M_scatter_whole (C, M, Mask_struct, GB_BITMAP_M_SCATTER_SET_2,
105             M_ek_slicing, M_ntasks, M_nthreads, Context) ;
106 
107         // C(i,j) has been marked, in Cb, with the value 2 where M(i,j)=1.
108         // These positions will not be computed in C(i,j).  C(i,j) can only
109         // be modified where Cb [p] is zero.
110 
111         int tid ;
112         #pragma omp parallel for num_threads(C_nthreads) schedule(static) \
113             reduction(+:cnvals)
114         for (tid = 0 ; tid < C_nthreads ; tid++)
115         {
116             int64_t pstart, pend, task_cnvals = 0 ;
117             GB_PARTITION (pstart, pend, cnz, tid, C_nthreads) ;
118             for (int64_t p = pstart ; p < pend ; p++)
119             {
120                 if (Cb [p] == 0)
121                 {
122                     // M(i,j) is zero, so C(i,j) can be computed
123                     if (GBB (Ab, p) && GBB (Bb, p))
124                     {
125                         // C (i,j) = A (i,j) + B (i,j)
126                         GB_GETA (aij, Ax, p) ;
127                         GB_GETB (bij, Bx, p) ;
128                         GB_BINOP (GB_CX (p), aij, bij, p % vlen, p / vlen) ;
129                         Cb [p] = 1 ;
130                         task_cnvals++ ;
131                     }
132                 }
133                 else
134                 {
135                     // M(i,j) == 1, so C(i,j) is not computed
136                     Cb [p] = 0 ;
137                 }
138             }
139             cnvals += task_cnvals ;
140         }
141 
142     }
143     else // if (ewise_method == GB_EMULT_METHOD_07)
144     {
145 
146         //----------------------------------------------------------------------
147         // C is bitmap; M is bitmap or full
148         //----------------------------------------------------------------------
149 
150         //      ------------------------------------------
151         //      C      <M> =        A       .*      B
152         //      ------------------------------------------
153         //      bitmap  bitmap      bitmap          bitmap  (method: 07)
154         //      bitmap  bitmap      bitmap          full    (method: 07)
155         //      bitmap  bitmap      full            bitmap  (method: 07)
156 
157         //      ------------------------------------------
158         //      C      <M> =        A       .*      B
159         //      ------------------------------------------
160         //      bitmap  full        bitmap          bitmap  (method: 07)
161         //      bitmap  full        bitmap          full    (method: 07)
162         //      bitmap  full        full            bitmap  (method: 07)
163 
164         //      ------------------------------------------
165         //      C      <!M> =       A       .*      B
166         //      ------------------------------------------
167         //      bitmap  bitmap      bitmap          bitmap  (method: 07)
168         //      bitmap  bitmap      bitmap          full    (method: 07)
169         //      bitmap  bitmap      full            bitmap  (method: 07)
170 
171         //      ------------------------------------------
172         //      C      <!M> =       A       .*      B
173         //      ------------------------------------------
174         //      bitmap  full        bitmap          bitmap  (method: 07)
175         //      bitmap  full        bitmap          full    (method: 07)
176         //      bitmap  full        full            bitmap  (method: 07)
177 
178         ASSERT (GB_IS_BITMAP (M) || GB_IS_FULL (M)) ;
179 
180         const int8_t  *restrict Mb = M->b ;
181         const GB_void *restrict Mx = (GB_void *) (Mask_struct ? NULL : (M->x)) ;
182         size_t msize = M->type->size ;
183 
184         #undef  GB_GET_MIJ
185         #define GB_GET_MIJ(p)                                           \
186             bool mij = GBB (Mb, p) && GB_mcast (Mx, p, msize) ;         \
187             if (Mask_comp) mij = !mij ; /* TODO: use ^ */
188 
189         int tid ;
190         #pragma omp parallel for num_threads(C_nthreads) schedule(static) \
191             reduction(+:cnvals)
192         for (tid = 0 ; tid < C_nthreads ; tid++)
193         {
194             int64_t pstart, pend, task_cnvals = 0 ;
195             GB_PARTITION (pstart, pend, cnz, tid, C_nthreads) ;
196             for (int64_t p = pstart ; p < pend ; p++)
197             {
198                 GB_GET_MIJ (p) ;
199                 if (mij)
200                 {
201                     // M(i,j) is true, so C(i,j) can be computed
202                     if (GBB (Ab, p) && GBB (Bb, p))
203                     {
204                         // C (i,j) = A (i,j) + B (i,j)
205                         GB_GETA (aij, Ax, p) ;
206                         GB_GETB (bij, Bx, p) ;
207                         GB_BINOP (GB_CX (p), aij, bij, p % vlen, p / vlen) ;
208                         Cb [p] = 1 ;
209                         task_cnvals++ ;
210                     }
211                 }
212                 else
213                 {
214                     // M(i,j) == 1, so C(i,j) is not computed
215                     Cb [p] = 0 ;
216                 }
217             }
218             cnvals += task_cnvals ;
219         }
220     }
221 
222     C->nvals = cnvals ;
223 }
224 
225