1 //------------------------------------------------------------------------------ 2 // GB_AxB_dot2_template: C=A'B, C<!M>=A'*B, or C<M>=A'*B via dot products 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 // A and B are sparse, bitmap, or full; never hypersparse. If the input 11 // matrices A and/or B are hypersparse, they are packed into sparse matrices, 12 // and C is unpacked from bitmap to sparse/hypersparse when done. 13 14 #if ( !GB_A_IS_HYPER && !GB_B_IS_HYPER ) 15 { 16 17 //-------------------------------------------------------------------------- 18 // C=A'*B, C<M>=A'*B, or C<!M>=A'*B where C is bitmap 19 //-------------------------------------------------------------------------- 20 21 int tid ; 22 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \ 23 reduction(+:cnvals) 24 for (tid = 0 ; tid < ntasks ; tid++) 25 { 26 27 //---------------------------------------------------------------------- 28 // get the task descriptor 29 //---------------------------------------------------------------------- 30 31 const int a_tid = tid / nbslice ; 32 const int b_tid = tid % nbslice ; 33 const int64_t kA_start = A_slice [a_tid] ; 34 const int64_t kA_end = A_slice [a_tid+1] ; 35 const int64_t kB_start = B_slice [b_tid] ; 36 const int64_t kB_end = B_slice [b_tid+1] ; 37 int64_t task_cnvals = 0 ; 38 39 //---------------------------------------------------------------------- 40 // C=A'*B, C<M>=A'*B, or C<!M>=A'*B via dot products 41 //---------------------------------------------------------------------- 42 43 for (int64_t j = kB_start ; j < kB_end ; j++) 44 { 45 46 //------------------------------------------------------------------ 47 // get B(:,j) and C(:,j) 48 //------------------------------------------------------------------ 49 50 const int64_t pC_start = j * cvlen ; 51 52 #if GB_B_IS_SPARSE 53 // B is sparse (never hypersparse) 54 const int64_t pB_start = Bp [j] ; 55 const int64_t pB_end = Bp [j+1] ; 56 const int64_t bjnz = pB_end - pB_start ; 57 if (bjnz == 0) 58 { 59 // no work to do if B(:,j) is empty, except to clear Cb 60 memset (&Cb [pC_start + kA_start], 0, kA_end - kA_start) ; 61 continue ; 62 } 63 #if GB_A_IS_SPARSE 64 // Both A and B are sparse; get first and last in B(:,j) 65 const int64_t ib_first = Bi [pB_start] ; 66 const int64_t ib_last = Bi [pB_end-1] ; 67 #endif 68 #else 69 // B is bitmap or full 70 const int64_t pB_start = j * vlen ; 71 #endif 72 73 //------------------------------------------------------------------ 74 // C(:,j)<#M(:,j)> = A'*B(:,j), or C(:,j) = A'*B(:,j) if no mask 75 //------------------------------------------------------------------ 76 77 for (int64_t i = kA_start ; i < kA_end ; i++) 78 { 79 80 //-------------------------------------------------------------- 81 // get C(i,j), M(i,j), and clear the C(i,j) bitmap 82 //-------------------------------------------------------------- 83 84 int64_t pC = pC_start + i ; // C is bitmap 85 86 #if defined ( GB_ANY_SPECIALIZED ) 87 // M is bitmap and structural; Mask_comp true 88 Cb [pC] = 0 ; 89 if (!Mb [pC]) 90 #elif defined ( GB_MASK_IS_PRESENT ) 91 bool mij ; 92 if (M_is_bitmap) 93 { 94 // M is bitmap 95 mij = Mb [pC] && GB_mcast (Mx, pC, msize) ; 96 } 97 else if (M_is_full) 98 { 99 // M is full 100 mij = GB_mcast (Mx, pC, msize) ; 101 } 102 else // M is sparse or hyper 103 { 104 // M has been scattered into the C bitmap 105 mij = (Cb [pC] > 1) ; 106 } 107 Cb [pC] = 0 ; 108 if (mij ^ Mask_comp) 109 #else 110 // M is not present 111 Cb [pC] = 0 ; 112 #endif 113 { 114 115 //---------------------------------------------------------- 116 // the mask allows C(i,j) to be computed 117 //---------------------------------------------------------- 118 119 #if GB_A_IS_SPARSE 120 // A is sparse 121 int64_t pA = Ap [i] ; 122 const int64_t pA_end = Ap [i+1] ; 123 const int64_t ainz = pA_end - pA ; 124 if (ainz > 0) 125 #else 126 // A is bitmap or full 127 const int64_t pA = i * vlen ; 128 #endif 129 { 130 // C(i,j) = A(:,i)'*B(:,j) 131 bool cij_exists = false ; 132 GB_CIJ_DECLARE (cij) ; 133 #include "GB_AxB_dot_cij.c" 134 } 135 } 136 } 137 } 138 cnvals += task_cnvals ; 139 } 140 } 141 #endif 142 143 #undef GB_A_IS_SPARSE 144 #undef GB_A_IS_HYPER 145 #undef GB_A_IS_BITMAP 146 #undef GB_A_IS_FULL 147 #undef GB_B_IS_SPARSE 148 #undef GB_B_IS_HYPER 149 #undef GB_B_IS_BITMAP 150 #undef GB_B_IS_FULL 151 152