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