1 //------------------------------------------------------------------------------ 2 // GB_meta16_definitions.h: methods that depend on the sparsity of A and 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 // Define macros that depend on the sparsity of A and B for GB_meta16_factory. 11 12 //------------------------------------------------------------------------------ 13 // GB_GET_B_j: prepare to iterate over B(:,j) 14 //------------------------------------------------------------------------------ 15 16 #undef GB_GET_B_j 17 18 #if defined ( GB_META16 ) 19 20 // this method appears in the GB_meta16_factory 21 #if ( GB_B_IS_HYPER || GB_A_IS_HYPER ) 22 23 // A or B are hyper 24 #define GB_GET_B_j \ 25 GB_GET_B_j_FOR_ALL_FORMATS (GB_A_IS_HYPER,GB_B_IS_SPARSE,GB_B_IS_HYPER) 26 27 #else 28 29 #if ( GB_B_IS_SPARSE ) 30 31 // B is sparse 32 #define GB_GET_B_j \ 33 int64_t j = kk ; \ 34 int64_t pB = Bp [kk] ; \ 35 int64_t pB_end = Bp [kk+1] ; \ 36 int64_t bjnz = pB_end - pB ; \ 37 GB_GET_T_FOR_SECONDJ 38 39 #else 40 41 // B is bitmap or full 42 #define GB_GET_B_j \ 43 int64_t j = kk ; \ 44 int64_t pB = kk * bvlen ; \ 45 int64_t pB_end = pB + bvlen ; \ 46 int64_t bjnz = bvlen ; \ 47 GB_GET_T_FOR_SECONDJ 48 49 #endif 50 51 #endif 52 53 #else 54 55 // define GB_GET_B_j for all sparsity formats 56 #define GB_GET_B_j \ 57 GB_GET_B_j_FOR_ALL_FORMATS (A_is_hyper, B_is_sparse, B_is_hyper) 58 59 #endif 60 61 //------------------------------------------------------------------------------ 62 // GB_GET_B_kj_INDEX: get the index k of the entry B(k,j) 63 //------------------------------------------------------------------------------ 64 65 #undef GB_GET_B_kj_INDEX 66 67 #if defined ( GB_META16 ) 68 69 #if ( GB_B_IS_HYPER || GB_B_IS_SPARSE ) 70 71 // B is hyper or sparse 72 #define GB_GET_B_kj_INDEX \ 73 int64_t k = Bi [pB] 74 75 #elif ( GB_B_IS_BITMAP ) 76 77 // B is bitmap 78 #define GB_GET_B_kj_INDEX \ 79 if (!Bb [pB]) continue ; \ 80 int64_t k = pB % bvlen 81 82 #else 83 84 // B is full 85 #define GB_GET_B_kj_INDEX \ 86 int64_t k = pB % bvlen 87 88 #endif 89 90 #else 91 92 // for any format of B 93 #define GB_GET_B_kj_INDEX \ 94 if (!GBB (Bb, pB)) continue ; \ 95 int64_t k = GBI (Bi, pB, bvlen) 96 97 #endif 98 99 //------------------------------------------------------------------------------ 100 // GB_GET_A_k: prepare to iterate over the vector A(:,k) 101 //------------------------------------------------------------------------------ 102 103 #undef GB_GET_A_k 104 105 #if defined ( GB_META16 ) 106 107 #if ( GB_A_IS_HYPER ) 108 109 // A is hyper 110 #define GB_GET_A_k GB_GET_A_k_FOR_ALL_FORMATS (true) 111 112 #elif ( GB_A_IS_SPARSE ) 113 114 // A is sparse 115 #define GB_GET_A_k \ 116 int64_t pA_start = Ap [k] ; \ 117 int64_t pA_end = Ap [k+1] ; \ 118 int64_t aknz = pA_end - pA_start 119 120 #else 121 122 // A is bitmap or full 123 #define GB_GET_A_k \ 124 int64_t pA_start = k * avlen ; \ 125 int64_t pA_end = pA_start + avlen ; \ 126 int64_t aknz = avlen 127 128 #endif 129 130 #else 131 132 // define GB_GET_A_k for all sparsity formats 133 #define GB_GET_A_k GB_GET_A_k_FOR_ALL_FORMATS (A_is_hyper) 134 135 #endif 136 137 //------------------------------------------------------------------------------ 138 // GB_GET_A_ik_INDEX: get the index i of the entry A(i,k) 139 //------------------------------------------------------------------------------ 140 141 #undef GB_GET_A_ik_INDEX 142 143 #if defined ( GB_META16 ) 144 145 #if ( GB_A_IS_HYPER || GB_A_IS_SPARSE ) 146 147 // A is hyper or sparse 148 #define GB_GET_A_ik_INDEX \ 149 int64_t i = Ai [pA] 150 151 #elif ( GB_A_IS_BITMAP ) 152 153 // A is bitmap 154 #define GB_GET_A_ik_INDEX \ 155 if (!Ab [pA]) continue ; \ 156 int64_t i = pA % avlen 157 158 #else 159 160 // A is full 161 #define GB_GET_A_ik_INDEX \ 162 int64_t i = pA % avlen 163 164 #endif 165 166 #else 167 168 // for any format of A 169 #define GB_GET_A_ik_INDEX \ 170 if (!GBB (Ab, pA)) continue ; \ 171 int64_t i = GBI (Ai, pA, avlen) 172 173 #endif 174 175 //------------------------------------------------------------------------------ 176 // GB_COMPUTE_C_j_WHEN_NNZ_B_j_IS_ONE: compute C(:,j) when nnz(B(:,j)) == 1 177 //------------------------------------------------------------------------------ 178 179 // C(:,j) = A(:,k)*B(k,j) when there is a single entry in B(:,j) 180 // The mask must not be present. A must be sparse or hypersparse. 181 182 #undef GB_COMPUTE_C_j_WHEN_NNZ_B_j_IS_ONE 183 184 #if GB_IS_ANY_PAIR_SEMIRING 185 186 // ANY_PAIR: result is purely symbolic; no numeric work to do 187 #define GB_COMPUTE_C_j_WHEN_NNZ_B_j_IS_ONE \ 188 ASSERT (A_is_sparse || A_is_hyper) ; \ 189 GB_GET_B_kj_INDEX ; /* get index k of B(k,j) */ \ 190 GB_GET_A_k ; /* get A(:,k) */ \ 191 memcpy (Ci + pC, Ai + pA_start, aknz * sizeof (int64_t)) ; \ 192 /* C becomes jumbled if A is jumbled */ \ 193 task_C_jumbled = task_C_jumbled || A_jumbled ; 194 195 #else 196 197 // typical semiring 198 #define GB_COMPUTE_C_j_WHEN_NNZ_B_j_IS_ONE \ 199 ASSERT (A_is_sparse || A_is_hyper) ; \ 200 GB_GET_B_kj_INDEX ; /* get index k of B(k,j) */ \ 201 GB_GET_A_k ; /* get A(:,k) */ \ 202 GB_GET_B_kj ; /* bkj = B(k,j) */ \ 203 /* scan A(:,k) */ \ 204 for (int64_t pA = pA_start ; pA < pA_end ; pA++) \ 205 { \ 206 GB_GET_A_ik_INDEX ; /* get index i of A(i,k) */ \ 207 GB_MULT_A_ik_B_kj ; /* t = A(i,k)*B(k,j) */ \ 208 GB_CIJ_WRITE (pC, t) ; /* Cx [pC] = t */ \ 209 Ci [pC++] = i ; \ 210 } \ 211 /* C becomes jumbled if A is jumbled */ \ 212 task_C_jumbled = task_C_jumbled || A_jumbled ; 213 214 #endif 215 216 //------------------------------------------------------------------------------ 217 // GB_COMPUTE_DENSE_C_j: compute C(:,j)=A*B(:,j) when C(:,j) is completely dense 218 //------------------------------------------------------------------------------ 219 220 // This method is not used for the saxpy3 generic method. 221 222 #undef GB_COMPUTE_DENSE_C_j 223 224 #if GB_IS_ANY_PAIR_SEMIRING 225 226 // ANY_PAIR: result is purely symbolic; no numeric work to do 227 #define GB_COMPUTE_DENSE_C_j \ 228 for (int64_t i = 0 ; i < cvlen ; i++) \ 229 { \ 230 Ci [pC + i] = i ; \ 231 } 232 233 #else 234 235 // typical semiring 236 #define GB_COMPUTE_DENSE_C_j \ 237 for (int64_t i = 0 ; i < cvlen ; i++) \ 238 { \ 239 Ci [pC + i] = i ; \ 240 GB_CIJ_WRITE (pC + i, GB_IDENTITY) ; /* C(i,j)=0 */ \ 241 } \ 242 for ( ; pB < pB_end ; pB++) /* scan B(:,j) */ \ 243 { \ 244 GB_GET_B_kj_INDEX ; /* get index k of B(k,j) */ \ 245 GB_GET_A_k ; /* get A(:,k) */ \ 246 if (aknz == 0) continue ; /* skip if A(:,k) empty */ \ 247 GB_GET_B_kj ; /* bkj = B(k,j) */ \ 248 /* scan A(:,k) */ \ 249 for (int64_t pA = pA_start ; pA < pA_end ; pA++) \ 250 { \ 251 GB_GET_A_ik_INDEX ; /* get index i of A(i,k) */ \ 252 GB_MULT_A_ik_B_kj ; /* t = A(i,k)*B(k,j) */ \ 253 GB_CIJ_UPDATE (pC + i, t) ; /* Cx [pC+i]+=t */ \ 254 } \ 255 } 256 257 #endif 258 259 //------------------------------------------------------------------------------ 260 // GB_SCAN_M_j_OR_A_k: compute C(:,j) using linear scan or binary search 261 //------------------------------------------------------------------------------ 262 263 // C(:,j)<M(:,j)>=A(:,k)*B(k,j) using one of two methods 264 #undef GB_SCAN_M_j_OR_A_k 265 #define GB_SCAN_M_j_OR_A_k(A_ok_for_binary_search) \ 266 { \ 267 if (A_ok_for_binary_search && aknz > 256 && mjnz_much < aknz && \ 268 mjnz < mvlen && aknz < avlen) \ 269 { \ 270 /* M and A are both sparse, and nnz(M(:,j)) is much less than */ \ 271 /* nnz(A(:,k)); scan M(:,j), and do binary search for A(i,k).*/ \ 272 /* This requires that A is not jumbled. */ \ 273 int64_t pA = pA_start ; \ 274 for (int64_t pM = pM_start ; pM < pM_end ; pM++) \ 275 { \ 276 GB_GET_M_ij (pM) ; /* get M(i,j) */ \ 277 if (!mij) continue ; /* skip if M(i,j)=0 */ \ 278 int64_t i = Mi [pM] ; \ 279 bool found ; /* search for A(i,k) */ \ 280 if (M_jumbled) pA = pA_start ; \ 281 int64_t apright = pA_end - 1 ; \ 282 GB_BINARY_SEARCH (i, Ai, pA, apright, found) ; \ 283 if (found) \ 284 { \ 285 /* C(i,j)<M(i,j)> += A(i,k) * B(k,j) for this method. */ \ 286 /* M(i,j) is always 1, as given in the hash table */ \ 287 GB_IKJ ; \ 288 } \ 289 } \ 290 } \ 291 else \ 292 { \ 293 /* A(:,j) is sparse enough relative to M(:,j) */ \ 294 /* M and/or A can dense, and either can be jumbled. */ \ 295 /* scan A(:,k), and lookup M(i,j) (in the hash table) */ \ 296 for (int64_t pA = pA_start ; pA < pA_end ; pA++) \ 297 { \ 298 GB_GET_A_ik_INDEX ; /* get index i of A(i,j) */ \ 299 /* do C(i,j)<M(i,j)> += A(i,k) * B(k,j) for this method */ \ 300 /* M(i,j) may be 0 or 1, as given in the hash table */ \ 301 GB_IKJ ; \ 302 } \ 303 } \ 304 } 305 306