1 //------------------------------------------------------------------------------ 2 // GB_AxB_dot3_template: C<M>=A'*B via dot products, where C is sparse/hyper 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 and M are both sparse or hyper, and C->h is a copy of M->h. 11 // M is present, and not complemented. It may be valued or structural. 12 13 { 14 15 int tid ; 16 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \ 17 reduction(+:nzombies) 18 for (tid = 0 ; tid < ntasks ; tid++) 19 { 20 21 //---------------------------------------------------------------------- 22 // get the task descriptor 23 //---------------------------------------------------------------------- 24 25 int64_t kfirst = TaskList [tid].kfirst ; 26 int64_t klast = TaskList [tid].klast ; 27 int64_t pC_first = TaskList [tid].pC ; 28 int64_t pC_last = TaskList [tid].pC_end ; 29 int64_t bpleft = 0 ; // Ch is not jumbled 30 int64_t task_nzombies = 0 ; // # of zombies found by this task 31 32 //---------------------------------------------------------------------- 33 // compute all vectors in this task 34 //---------------------------------------------------------------------- 35 36 for (int64_t k = kfirst ; k <= klast ; k++) 37 { 38 39 //------------------------------------------------------------------ 40 // get C(:,k) and M(:k) 41 //------------------------------------------------------------------ 42 43 #if defined ( GB_MASK_SPARSE_AND_STRUCTURAL ) 44 // M and C are sparse 45 const int64_t j = k ; 46 #else 47 // M and C are either both sparse or both hypersparse 48 const int64_t j = GBH (Ch, k) ; 49 #endif 50 51 int64_t pC_start = Cp [k] ; 52 int64_t pC_end = Cp [k+1] ; 53 if (k == kfirst) 54 { 55 // First vector for task; may only be partially owned. 56 pC_start = pC_first ; 57 pC_end = GB_IMIN (pC_end, pC_last) ; 58 } 59 else if (k == klast) 60 { 61 // Last vector for task; may only be partially owned. 62 pC_end = pC_last ; 63 } 64 else 65 { 66 // task completely owns this vector C(:,k). 67 } 68 69 //------------------------------------------------------------------ 70 // get B(:,j) 71 //------------------------------------------------------------------ 72 73 #if GB_B_IS_HYPER 74 // B is hyper 75 int64_t pB_start, pB_end ; 76 GB_lookup (true, Bh, Bp, vlen, &bpleft, bnvec-1, j, 77 &pB_start, &pB_end) ; 78 #elif GB_B_IS_SPARSE 79 // B is sparse 80 const int64_t pB_start = Bp [j] ; 81 const int64_t pB_end = Bp [j+1] ; 82 #else 83 // B is bitmap or full 84 const int64_t pB_start = j * vlen ; 85 #endif 86 87 #if (GB_B_IS_SPARSE || GB_B_IS_HYPER) 88 const int64_t bjnz = pB_end - pB_start ; 89 if (bjnz == 0) 90 { 91 // no work to do if B(:,j) is empty, except for zombies 92 task_nzombies += (pC_end - pC_start) ; 93 for (int64_t pC = pC_start ; pC < pC_end ; pC++) 94 { 95 // C(i,j) is a zombie 96 int64_t i = Mi [pC] ; 97 Ci [pC] = GB_FLIP (i) ; 98 } 99 continue ; 100 } 101 #if (GB_A_IS_SPARSE || GB_A_IS_HYPER) 102 // Both A and B are sparse; get first and last in B(:,j) 103 const int64_t ib_first = Bi [pB_start] ; 104 const int64_t ib_last = Bi [pB_end-1] ; 105 #endif 106 #endif 107 108 //------------------------------------------------------------------ 109 // C(:,j)<M(:,j)> = A(:,i)'*B(:,j) 110 //------------------------------------------------------------------ 111 112 for (int64_t pC = pC_start ; pC < pC_end ; pC++) 113 { 114 115 //-------------------------------------------------------------- 116 // get C(i,j) and M(i,j) 117 //-------------------------------------------------------------- 118 119 bool cij_exists = false ; 120 GB_CIJ_DECLARE (cij) ; 121 122 // get the value of M(i,j) 123 int64_t i = Mi [pC] ; 124 #if !defined ( GB_MASK_SPARSE_AND_STRUCTURAL ) 125 // if M is structural, no need to check its values 126 if (GB_mcast (Mx, pC, msize)) 127 #endif 128 { 129 130 //---------------------------------------------------------- 131 // the mask allows C(i,j) to be computed 132 //---------------------------------------------------------- 133 134 #if GB_A_IS_HYPER 135 // A is hyper 136 int64_t pA, pA_end ; 137 int64_t apleft = 0 ; // M might be jumbled 138 GB_lookup (true, Ah, Ap, vlen, &apleft, anvec-1, i, 139 &pA, &pA_end) ; 140 const int64_t ainz = pA_end - pA ; 141 if (ainz > 0) 142 #elif GB_A_IS_SPARSE 143 // A is sparse 144 int64_t pA = Ap [i] ; 145 const int64_t pA_end = Ap [i+1] ; 146 const int64_t ainz = pA_end - pA ; 147 if (ainz > 0) 148 #else 149 // A is bitmap or full 150 const int64_t pA = i * vlen ; 151 #endif 152 { 153 // C(i,j) = A(:,i)'*B(:,j) 154 #include "GB_AxB_dot_cij.c" 155 } 156 } 157 158 if (!GB_CIJ_EXISTS) 159 { 160 // C(i,j) is a zombie 161 task_nzombies++ ; 162 Ci [pC] = GB_FLIP (i) ; 163 } 164 } 165 } 166 nzombies += task_nzombies ; 167 } 168 } 169 170 #undef GB_A_IS_SPARSE 171 #undef GB_A_IS_HYPER 172 #undef GB_A_IS_BITMAP 173 #undef GB_A_IS_FULL 174 #undef GB_B_IS_SPARSE 175 #undef GB_B_IS_HYPER 176 #undef GB_B_IS_BITMAP 177 #undef GB_B_IS_FULL 178 179