1 //------------------------------------------------------------------------------ 2 // GB_select_phase2: C=select(A,thunk) 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 { 11 //-------------------------------------------------------------------------- 12 // get A 13 //-------------------------------------------------------------------------- 14 15 const int64_t *restrict Ap = A->p ; 16 const int64_t *restrict Ah = A->h ; 17 const int64_t *restrict Ai = A->i ; 18 const GB_ATYPE *restrict Ax = (GB_ATYPE *) A->x ; 19 size_t asize = A->type->size ; 20 int64_t avlen = A->vlen ; 21 int64_t avdim = A->vdim ; 22 // if A is bitmap, the bitmap selector is always used instead 23 ASSERT (!GB_IS_BITMAP (A)) ; 24 #ifndef GB_DIAG_SELECTOR 25 // if A is full, all opcodes except DIAG use the bitmap selector instead 26 ASSERT (!GB_IS_FULL (A)) ; 27 #endif 28 29 const int64_t *restrict kfirst_Aslice = A_ek_slicing ; 30 const int64_t *restrict klast_Aslice = A_ek_slicing + A_ntasks ; 31 const int64_t *restrict pstart_Aslice = A_ek_slicing + A_ntasks * 2 ; 32 33 //-------------------------------------------------------------------------- 34 // C = select (A) 35 //-------------------------------------------------------------------------- 36 37 int tid ; 38 #pragma omp parallel for num_threads(A_nthreads) schedule(dynamic,1) 39 for (tid = 0 ; tid < A_ntasks ; tid++) 40 { 41 42 // if kfirst > klast then task tid does no work at all 43 int64_t kfirst = kfirst_Aslice [tid] ; 44 int64_t klast = klast_Aslice [tid] ; 45 46 //---------------------------------------------------------------------- 47 // selection from vectors kfirst to klast 48 //---------------------------------------------------------------------- 49 50 for (int64_t k = kfirst ; k <= klast ; k++) 51 { 52 53 //------------------------------------------------------------------ 54 // find the part of A(:,k) to be operated on by this task 55 //------------------------------------------------------------------ 56 57 int64_t pA_start, pA_end, pC ; 58 GB_get_pA_and_pC (&pA_start, &pA_end, &pC, tid, k, kfirst, klast, 59 pstart_Aslice, Cp_kfirst, Cp, avlen, Ap, avlen) ; 60 61 //------------------------------------------------------------------ 62 // compact Ai and Ax [pA_start ... pA_end-1] into Ci and Cx 63 //------------------------------------------------------------------ 64 65 #if defined ( GB_ENTRY_SELECTOR ) 66 67 GB_GET_J ; 68 for (int64_t pA = pA_start ; pA < pA_end ; pA++) 69 { 70 // A is never full; that case is now handled by the 71 // bitmap selector instead. 72 // int64_t i = GBI (Ai, pA, avlen) ; 73 ASSERT (Ai != NULL) ; 74 int64_t i = Ai [pA] ; 75 if (GB_TEST_VALUE_OF_ENTRY (pA)) 76 { 77 ASSERT (pC >= Cp [k] && pC < Cp [k+1]) ; 78 Ci [pC] = i ; 79 // Cx [pC] = Ax [pA] ; 80 GB_SELECT_ENTRY (Cx, pC, Ax, pA) ; 81 pC++ ; 82 } 83 } 84 85 #elif defined ( GB_TRIU_SELECTOR ) \ 86 || defined ( GB_RESIZE_SELECTOR ) 87 88 // keep pA_start to Zp[k]-1 89 int64_t p = GB_IMIN (Zp [k], pA_end) ; 90 int64_t mynz = p - pA_start ; 91 if (mynz > 0) 92 { 93 ASSERT (pC >= Cp [k] && pC + mynz <= Cp [k+1]) ; 94 ASSERT (Ai != NULL) ; 95 // if (Ai != NULL) 96 { 97 // A and C are both sparse or hypersparse 98 memcpy (Ci +pC, Ai +pA_start, mynz*sizeof (int64_t)) ; 99 } 100 #if 0 101 else 102 { 103 // A is full and C is sparse: for triu: the bitmap 104 // selector is used. For resize, A is converted to 105 // hypersparse first. 106 ASSERT (GB_DEAD_CODE) ; 107 int64_t i_start = pA_start % avlen ; 108 for (int64_t s = 0 ; s < mynz ; s++) 109 { 110 int64_t i = i_start + s ; 111 ASSERT (GBI (Ai, pA_start+s, avlen) == i) ; 112 Ci [pC+s] = i ; 113 } 114 } 115 #endif 116 memcpy (Cx +pC*asize, Ax +pA_start*asize, mynz*asize) ; 117 } 118 119 #elif defined ( GB_DIAG_SELECTOR ) 120 121 // task that owns the diagonal entry does this work 122 // A can be sparse or full, but not bitmap 123 int64_t p = Zp [k] ; 124 if (pA_start <= p && p < pA_end) 125 { 126 ASSERT (pC >= Cp [k] && pC + 1 <= Cp [k+1]) ; 127 Ci [pC] = GBI (Ai, p, avlen) ; 128 memcpy (Cx +pC*asize, Ax +p*asize, asize) ; 129 } 130 131 #elif defined ( GB_OFFDIAG_SELECTOR ) 132 133 // keep pA_start to Zp[k]-1 134 int64_t p = GB_IMIN (Zp [k], pA_end) ; 135 int64_t mynz = p - pA_start ; 136 if (mynz > 0) 137 { 138 ASSERT (pC >= Cp [k] && pC + mynz <= Cp [k+1]) ; 139 ASSERT (Ai != NULL) ; 140 // if (Ai != NULL) 141 { 142 // A and C are both sparse or hypersparse 143 memcpy (Ci +pC, Ai +pA_start, mynz*sizeof (int64_t)) ; 144 } 145 #if 0 146 else 147 { 148 // A is full and C is sparse or hypersparse: 149 // this is now always handled by the bitmap selector 150 ASSERT (GB_DEAD_CODE) ; 151 int64_t i_start = pA_start % avlen ; 152 for (int64_t s = 0 ; s < mynz ; s++) 153 { 154 int64_t i = i_start + s ; 155 ASSERT (GBI (Ai, pA_start+s, avlen) == i) ; 156 Ci [pC+s] = i ; 157 } 158 } 159 #endif 160 memcpy (Cx +pC*asize, Ax +pA_start*asize, mynz*asize) ; 161 pC += mynz ; 162 } 163 164 // keep Zp[k]+1 to pA_end-1 165 p = GB_IMAX (Zp [k]+1, pA_start) ; 166 mynz = pA_end - p ; 167 if (mynz > 0) 168 { 169 ASSERT (pA_start <= p && p < pA_end) ; 170 ASSERT (pC >= Cp [k] && pC + mynz <= Cp [k+1]) ; 171 ASSERT (Ai != NULL) ; 172 // if (Ai != NULL) 173 { 174 // A and C are both sparse or hypersparse 175 memcpy (Ci +pC, Ai +p, mynz*sizeof (int64_t)) ; 176 } 177 #if 0 178 else 179 { 180 // A is full and C is sparse or hypersparse 181 ASSERT (GB_DEAD_CODE) ; 182 int64_t i_start = p % avlen ; 183 for (int64_t s = 0 ; s < mynz ; s++) 184 { 185 int64_t i = i_start + s ; 186 ASSERT (GBI (Ai, p+s, avlen) == i) ; 187 Ci [pC+s] = i ; 188 } 189 } 190 #endif 191 memcpy (Cx +pC*asize, Ax +p*asize, mynz*asize) ; 192 } 193 194 #elif defined ( GB_TRIL_SELECTOR ) 195 196 // keep Zp [k] to pA_end-1 197 int64_t p = GB_IMAX (Zp [k], pA_start) ; 198 int64_t mynz = pA_end - p ; 199 if (mynz > 0) 200 { 201 ASSERT (pA_start <= p && p + mynz <= pA_end) ; 202 ASSERT (pC >= Cp [k] && pC + mynz <= Cp [k+1]) ; 203 ASSERT (Ai != NULL) ; 204 // if (Ai != NULL) 205 { 206 // A and C are both sparse or hypersparse 207 memcpy (Ci +pC, Ai +p, mynz*sizeof (int64_t)) ; 208 } 209 #if 0 210 else 211 { 212 // A is full and C is sparse or hypersparse: 213 // this is now always handled by the bitmap selector 214 ASSERT (GB_DEAD_CODE) ; 215 int64_t i_start = p % avlen ; 216 for (int64_t s = 0 ; s < mynz ; s++) 217 { 218 int64_t i = i_start + s ; 219 ASSERT (GBI (Ai, p+s, avlen) == i) ; 220 Ci [pC+s] = i ; 221 } 222 } 223 #endif 224 memcpy (Cx +pC*asize, Ax +p*asize, mynz*asize) ; 225 } 226 227 #endif 228 } 229 } 230 } 231 232