1 //------------------------------------------------------------------------------ 2 // GB_emult_02_template: C = A.*B when A is sparse/hyper and B is bitmap/full 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 sparse, with the same sparsity structure as A. No mask is present, or 11 // M is bitmap/full. A is sparse/hyper, and B is bitmap/full. This method 12 // also handles the case when the original input A is bitmap/full and B is 13 // sparse/hyper, by computing B.*A with the operator flipped. 14 15 { 16 17 //-------------------------------------------------------------------------- 18 // get A, B, and C 19 //-------------------------------------------------------------------------- 20 21 const int64_t *restrict Ap = A->p ; 22 const int64_t *restrict Ah = A->h ; 23 const int64_t *restrict Ai = A->i ; 24 const int64_t vlen = A->vlen ; 25 26 const int8_t *restrict Bb = B->b ; 27 28 const int64_t *restrict kfirst_Aslice = A_ek_slicing ; 29 const int64_t *restrict klast_Aslice = A_ek_slicing + A_ntasks ; 30 const int64_t *restrict pstart_Aslice = A_ek_slicing + A_ntasks * 2 ; 31 32 #if GB_FLIPPED 33 const GB_BTYPE *restrict Ax = (GB_BTYPE *) A->x ; 34 const GB_ATYPE *restrict Bx = (GB_ATYPE *) B->x ; 35 #else 36 const GB_ATYPE *restrict Ax = (GB_ATYPE *) A->x ; 37 const GB_BTYPE *restrict Bx = (GB_BTYPE *) B->x ; 38 #endif 39 40 const int64_t *restrict Cp = C->p ; 41 int64_t *restrict Ci = C->i ; 42 GB_CTYPE *restrict Cx = (GB_CTYPE *) C->x ; 43 44 //-------------------------------------------------------------------------- 45 // C=A.*B or C<#M>=A.*B 46 //-------------------------------------------------------------------------- 47 48 if (M == NULL) 49 { 50 51 //---------------------------------------------------------------------- 52 // C = A.*B 53 //---------------------------------------------------------------------- 54 55 if (GB_IS_BITMAP (B)) 56 { 57 58 //------------------------------------------------------------------ 59 // C=A.*B where A is sparse/hyper and B is bitmap 60 //------------------------------------------------------------------ 61 62 int tid ; 63 #pragma omp parallel for num_threads(A_nthreads) schedule(dynamic,1) 64 for (tid = 0 ; tid < A_ntasks ; tid++) 65 { 66 int64_t kfirst = kfirst_Aslice [tid] ; 67 int64_t klast = klast_Aslice [tid] ; 68 for (int64_t k = kfirst ; k <= klast ; k++) 69 { 70 int64_t j = GBH (Ah, k) ; 71 int64_t pB_start = j * vlen ; 72 int64_t pA, pA_end, pC ; 73 GB_get_pA_and_pC (&pA, &pA_end, &pC, tid, k, kfirst, klast, 74 pstart_Aslice, Cp_kfirst, Cp, vlen, Ap, vlen) ; 75 for ( ; pA < pA_end ; pA++) 76 { 77 int64_t i = Ai [pA] ; 78 int64_t pB = pB_start + i ; 79 if (!Bb [pB]) continue ; 80 // C (i,j) = A (i,j) .* B (i,j) 81 Ci [pC] = i ; 82 GB_GETA (aij, Ax, pA) ; 83 GB_GETB (bij, Bx, pB) ; 84 #if GB_FLIPPED 85 GB_BINOP (GB_CX (pC), bij, aij, i, j) ; 86 #else 87 GB_BINOP (GB_CX (pC), aij, bij, i, j) ; 88 #endif 89 pC++ ; 90 } 91 } 92 } 93 94 } 95 else 96 { 97 98 //------------------------------------------------------------------ 99 // C=A.*B where A is sparse/hyper and B is full 100 //------------------------------------------------------------------ 101 102 int tid ; 103 #pragma omp parallel for num_threads(A_nthreads) schedule(dynamic,1) 104 for (tid = 0 ; tid < A_ntasks ; tid++) 105 { 106 int64_t kfirst = kfirst_Aslice [tid] ; 107 int64_t klast = klast_Aslice [tid] ; 108 for (int64_t k = kfirst ; k <= klast ; k++) 109 { 110 int64_t j = GBH (Ah, k) ; 111 int64_t pB_start = j * vlen ; 112 int64_t pA, pA_end ; 113 GB_get_pA (&pA, &pA_end, tid, k, kfirst, klast, 114 pstart_Aslice, Ap, vlen) ; 115 for ( ; pA < pA_end ; pA++) 116 { 117 // C (i,j) = A (i,j) .* B (i,j) 118 int64_t i = Ai [pA] ; 119 int64_t pB = pB_start + i ; 120 // Ci [pA] = i ; already defined 121 GB_GETA (aij, Ax, pA) ; 122 GB_GETB (bij, Bx, pB) ; 123 #if GB_FLIPPED 124 GB_BINOP (GB_CX (pA), bij, aij, i, j) ; 125 #else 126 GB_BINOP (GB_CX (pA), aij, bij, i, j) ; 127 #endif 128 } 129 } 130 } 131 } 132 133 } 134 else 135 { 136 137 //---------------------------------------------------------------------- 138 // C<#M>=A.*B where A is sparse/hyper, M and B are bitmap/full 139 //---------------------------------------------------------------------- 140 141 const int8_t *restrict Mb = M->b ; 142 const GB_void *restrict Mx = (Mask_struct) ? NULL : ((GB_void *) M->x) ; 143 const size_t msize = M->type->size ; 144 145 int tid ; 146 #pragma omp parallel for num_threads(A_nthreads) schedule(dynamic,1) 147 for (tid = 0 ; tid < A_ntasks ; tid++) 148 { 149 int64_t kfirst = kfirst_Aslice [tid] ; 150 int64_t klast = klast_Aslice [tid] ; 151 for (int64_t k = kfirst ; k <= klast ; k++) 152 { 153 int64_t j = GBH (Ah, k) ; 154 int64_t pB_start = j * vlen ; 155 int64_t pA, pA_end, pC ; 156 GB_get_pA_and_pC (&pA, &pA_end, &pC, tid, k, kfirst, klast, 157 pstart_Aslice, Cp_kfirst, Cp, vlen, Ap, vlen) ; 158 for ( ; pA < pA_end ; pA++) 159 { 160 int64_t i = Ai [pA] ; 161 int64_t pB = pB_start + i ; 162 if (!GBB (Bb, pB)) continue ; 163 bool mij = GBB (Mb, pB) && GB_mcast (Mx, pB, msize) ; 164 mij = mij ^ Mask_comp ; 165 if (!mij) continue ; 166 // C (i,j) = A (i,j) .* B (i,j) 167 Ci [pC] = i ; 168 GB_GETA (aij, Ax, pA) ; 169 GB_GETB (bij, Bx, pB) ; 170 #if GB_FLIPPED 171 GB_BINOP (GB_CX (pC), bij, aij, i, j) ; 172 #else 173 GB_BINOP (GB_CX (pC), aij, bij, i, j) ; 174 #endif 175 pC++ ; 176 } 177 } 178 } 179 } 180 } 181 182