|
| 1 | +// bench_nvwarp_l1.cu — warp-scope vendor A/B: glass::warp:: (SIMT shuffles) vs |
| 2 | +// glass::nvidia::warp:: (cub::WarpReduce) on the three ops the vendor warp tier |
| 3 | +// ships (reduce / dot / nrm2). Same launch shape for both legs: 8 warps per |
| 4 | +// block, ONE problem per warp (the warp packing model), full warps only. |
| 5 | +// |
| 6 | +// Compilation: nvcc -std=c++17 -arch=native -O3 --expt-relaxed-constexpr -I.. bench_nvwarp_l1.cu -o bench_nvwarp_l1 |
| 7 | +// Usage: ./bench_nvwarp_l1 |
| 8 | +// |
| 9 | +// Every leg is correctness-gated against a host reference before it is timed — |
| 10 | +// a wrong result or failed launch prints FAIL and exits nonzero rather than |
| 11 | +// timing garbage. |
| 12 | + |
| 13 | +#include <cstdio> |
| 14 | +#include <cstdlib> |
| 15 | +#include <cstring> |
| 16 | +#include <ctime> |
| 17 | +#include <cmath> |
| 18 | +#include <vector> |
| 19 | +#include "../glass.cuh" |
| 20 | +#include "../glass-nvidia.cuh" // glass::nvidia::warp::{reduce,dot,nrm2} |
| 21 | + |
| 22 | +static const int WARPS = 8; // warps (= problems) per block |
| 23 | +static int NPROB = 8192; // set per section in main |
| 24 | + |
| 25 | +static double elapsed_ms(struct timespec a, struct timespec b) { |
| 26 | + return (b.tv_sec - a.tv_sec) * 1e3 + (b.tv_nsec - a.tv_nsec) * 1e-6; |
| 27 | +} |
| 28 | + |
| 29 | +// ─── kernels: p = global warp index = one problem per warp ─────────────────── |
| 30 | +template<typename T, uint32_t N> |
| 31 | +__global__ void kw_reduce(T* x, int np) { |
| 32 | + int p = blockIdx.x * blockDim.y + threadIdx.y; if (p >= np) return; |
| 33 | + glass::warp::reduce<T, N>(x + (size_t)p * N); |
| 34 | +} |
| 35 | +template<typename T, uint32_t N> |
| 36 | +__global__ void kw_dot(T* x, T* y, T* out, int np) { |
| 37 | + int p = blockIdx.x * blockDim.y + threadIdx.y; if (p >= np) return; |
| 38 | + T r = glass::warp::dot<T, N>(x + (size_t)p * N, y + (size_t)p * N); |
| 39 | + if ((threadIdx.x & 31) == 0) out[p] = r; |
| 40 | +} |
| 41 | +template<typename T, uint32_t N> |
| 42 | +__global__ void kw_nrm2(T* x, T* out, int np) { |
| 43 | + int p = blockIdx.x * blockDim.y + threadIdx.y; if (p >= np) return; |
| 44 | + T r = glass::warp::nrm2<T, N>(x + (size_t)p * N); |
| 45 | + if ((threadIdx.x & 31) == 0) out[p] = r; |
| 46 | +} |
| 47 | + |
| 48 | +template<typename T> |
| 49 | +__device__ T* warp_scratch(char* s) { |
| 50 | + constexpr size_t B = glass::nvidia::warp::warp_reduce_scratch_bytes<T>(); |
| 51 | + return reinterpret_cast<T*>(s + threadIdx.y * (B ? B : sizeof(T))); |
| 52 | +} |
| 53 | +template<typename T, uint32_t N> |
| 54 | +__global__ void kn_reduce(T* x, int np) { |
| 55 | + extern __shared__ char s[]; |
| 56 | + int p = blockIdx.x * blockDim.y + threadIdx.y; if (p >= np) return; |
| 57 | + glass::nvidia::warp::reduce<T, N>(x + (size_t)p * N, warp_scratch<T>(s)); |
| 58 | +} |
| 59 | +template<typename T, uint32_t N> |
| 60 | +__global__ void kn_dot(T* x, T* y, T* out, int np) { |
| 61 | + extern __shared__ char s[]; |
| 62 | + int p = blockIdx.x * blockDim.y + threadIdx.y; if (p >= np) return; |
| 63 | + glass::nvidia::warp::dot<T, N>(x + (size_t)p * N, y + (size_t)p * N, |
| 64 | + out + p, warp_scratch<T>(s)); |
| 65 | +} |
| 66 | +template<typename T, uint32_t N> |
| 67 | +__global__ void kn_nrm2(T* x, T* out, int np) { |
| 68 | + extern __shared__ char s[]; |
| 69 | + int p = blockIdx.x * blockDim.y + threadIdx.y; if (p >= np) return; |
| 70 | + glass::nvidia::warp::nrm2<T, N>(x + (size_t)p * N, out + p, warp_scratch<T>(s)); |
| 71 | +} |
| 72 | + |
| 73 | +// ─── harness ───────────────────────────────────────────────────────────────── |
| 74 | +template<typename F> |
| 75 | +static double time_ns_per_prob(F launch, int reps) { |
| 76 | + cudaGetLastError(); |
| 77 | + launch(); cudaDeviceSynchronize(); |
| 78 | + if (cudaGetLastError() != cudaSuccess) return -1.0; |
| 79 | + double best = 1e30; |
| 80 | + for (int t = 0; t < 3; t++) { |
| 81 | + struct timespec t0, t1; |
| 82 | + clock_gettime(CLOCK_MONOTONIC, &t0); |
| 83 | + for (int r = 0; r < reps; r++) launch(); |
| 84 | + cudaDeviceSynchronize(); |
| 85 | + clock_gettime(CLOCK_MONOTONIC, &t1); |
| 86 | + double ns = elapsed_ms(t0, t1) * 1e6 / ((double)reps * NPROB); |
| 87 | + if (ns < best) best = ns; |
| 88 | + } |
| 89 | + return best; |
| 90 | +} |
| 91 | + |
| 92 | +template<typename T> |
| 93 | +static void fill(std::vector<T>& h) { |
| 94 | + for (size_t i = 0; i < h.size(); i++) |
| 95 | + h[i] = (T)0.25 * (T)((int)(i % 7) - 3); |
| 96 | +} |
| 97 | + |
| 98 | +template<typename T> |
| 99 | +static bool close(T a, T b) { |
| 100 | + double tol = sizeof(T) == 4 ? 1e-4 : 1e-10; |
| 101 | + double d = std::fabs((double)a - (double)b); |
| 102 | + return d <= tol * (1.0 + std::fabs((double)b)); |
| 103 | +} |
| 104 | + |
| 105 | +template<typename T, uint32_t N> |
| 106 | +static void run_op(const char* dt) { |
| 107 | + const size_t elems = (size_t)NPROB * N; |
| 108 | + std::vector<T> hx(elems), hy(elems), hout(NPROB); |
| 109 | + fill(hx); fill(hy); |
| 110 | + for (size_t i = 0; i < elems; i++) hy[i] = hy[(i * 31 + 7) % elems]; |
| 111 | + T *dx, *dy, *dout; |
| 112 | + cudaMalloc(&dx, elems * sizeof(T)); cudaMalloc(&dy, elems * sizeof(T)); |
| 113 | + cudaMalloc(&dout, NPROB * sizeof(T)); |
| 114 | + const dim3 blk(32, WARPS); |
| 115 | + const int grid = (NPROB + WARPS - 1) / WARPS; |
| 116 | + constexpr size_t SB = glass::nvidia::warp::warp_reduce_scratch_bytes<T>(); |
| 117 | + const size_t smem = WARPS * (SB ? SB : sizeof(T)); |
| 118 | + const int reps = NPROB >= 8192 ? 250 : NPROB >= 1024 ? 500 : 2000; |
| 119 | + |
| 120 | + // host references (problem 0 and NPROB-1 spot-checked per problem below) |
| 121 | + auto refs = [&](int p, double& rsum, double& rdot, double& rnrm) { |
| 122 | + rsum = rdot = rnrm = 0.0; |
| 123 | + for (uint32_t i = 0; i < N; i++) { |
| 124 | + double xi = (double)hx[(size_t)p * N + i], yi = (double)hy[(size_t)p * N + i]; |
| 125 | + rsum += xi; rdot += xi * yi; rnrm += xi * xi; |
| 126 | + } |
| 127 | + rnrm = std::sqrt(rnrm); |
| 128 | + }; |
| 129 | + auto check = [&](const char* leg, int which) -> bool { // 0=reduce 1=dot 2=nrm2 |
| 130 | + cudaMemcpy(dx, hx.data(), elems * sizeof(T), cudaMemcpyHostToDevice); |
| 131 | + cudaMemcpy(dy, hy.data(), elems * sizeof(T), cudaMemcpyHostToDevice); |
| 132 | + cudaMemset(dout, 0, NPROB * sizeof(T)); |
| 133 | + switch (which) { |
| 134 | + case 0: strcmp(leg, "cub") ? kw_reduce<T, N><<<grid, blk>>>(dx, NPROB) |
| 135 | + : kn_reduce<T, N><<<grid, blk, smem>>>(dx, NPROB); break; |
| 136 | + case 1: strcmp(leg, "cub") ? kw_dot<T, N><<<grid, blk>>>(dx, dy, dout, NPROB) |
| 137 | + : kn_dot<T, N><<<grid, blk, smem>>>(dx, dy, dout, NPROB); break; |
| 138 | + case 2: strcmp(leg, "cub") ? kw_nrm2<T, N><<<grid, blk>>>(dx, dout, NPROB) |
| 139 | + : kn_nrm2<T, N><<<grid, blk, smem>>>(dx, dout, NPROB); break; |
| 140 | + } |
| 141 | + if (cudaDeviceSynchronize() != cudaSuccess) { printf("FAIL launch %s\n", leg); return false; } |
| 142 | + std::vector<T> gx(elems), gout(NPROB); |
| 143 | + cudaMemcpy(gx.data(), dx, elems * sizeof(T), cudaMemcpyDeviceToHost); |
| 144 | + cudaMemcpy(gout.data(), dout, NPROB * sizeof(T), cudaMemcpyDeviceToHost); |
| 145 | + for (int p : {0, NPROB - 1}) { |
| 146 | + double rs, rd, rn; refs(p, rs, rd, rn); |
| 147 | + T got = which == 0 ? gx[(size_t)p * N] : gout[p]; |
| 148 | + double want = which == 0 ? rs : which == 1 ? rd : rn; |
| 149 | + if (!close(got, (T)want)) { |
| 150 | + printf("FAIL check op=%d leg=%s N=%u p=%d got=%g want=%g\n", |
| 151 | + which, leg, N, p, (double)got, want); |
| 152 | + return false; |
| 153 | + } |
| 154 | + } |
| 155 | + return true; |
| 156 | + }; |
| 157 | + |
| 158 | + const char* names[3] = {"reduce", "dot", "nrm2"}; |
| 159 | + for (int which = 0; which < 3; which++) { |
| 160 | + if (!check("simt", which) || !check("cub", which)) exit(1); |
| 161 | + cudaMemcpy(dx, hx.data(), elems * sizeof(T), cudaMemcpyHostToDevice); |
| 162 | + cudaMemcpy(dy, hy.data(), elems * sizeof(T), cudaMemcpyHostToDevice); |
| 163 | + double simt = -1, cub = -1; |
| 164 | + switch (which) { |
| 165 | + case 0: |
| 166 | + simt = time_ns_per_prob([&]{ kw_reduce<T, N><<<grid, blk>>>(dx, NPROB); }, reps); |
| 167 | + cub = time_ns_per_prob([&]{ kn_reduce<T, N><<<grid, blk, smem>>>(dx, NPROB); }, reps); |
| 168 | + break; |
| 169 | + case 1: |
| 170 | + simt = time_ns_per_prob([&]{ kw_dot<T, N><<<grid, blk>>>(dx, dy, dout, NPROB); }, reps); |
| 171 | + cub = time_ns_per_prob([&]{ kn_dot<T, N><<<grid, blk, smem>>>(dx, dy, dout, NPROB); }, reps); |
| 172 | + break; |
| 173 | + case 2: |
| 174 | + simt = time_ns_per_prob([&]{ kw_nrm2<T, N><<<grid, blk>>>(dx, dout, NPROB); }, reps); |
| 175 | + cub = time_ns_per_prob([&]{ kn_nrm2<T, N><<<grid, blk, smem>>>(dx, dout, NPROB); }, reps); |
| 176 | + break; |
| 177 | + } |
| 178 | + printf("%-6s N=%-3u | simt=%8.3f cub=%8.3f | cub/simt=%.3f -> %s\n", |
| 179 | + names[which], N, simt, cub, cub / simt, |
| 180 | + cub < simt * 0.98 ? "CUB" : simt < cub * 0.98 ? "SIMT" : "tie"); |
| 181 | + } |
| 182 | + cudaFree(dx); cudaFree(dy); cudaFree(dout); |
| 183 | +} |
| 184 | + |
| 185 | +template<typename T> |
| 186 | +static void run_dtype(const char* dt) { |
| 187 | + printf("#### NPROB=%d dtype=%s (8 warps/block, 1 problem/warp) ####\n", NPROB, dt); |
| 188 | + run_op<T, 4>(dt); run_op<T, 8>(dt); run_op<T, 16>(dt); run_op<T, 32>(dt); |
| 189 | + run_op<T, 64>(dt); run_op<T, 128>(dt); run_op<T, 256>(dt); |
| 190 | +} |
| 191 | + |
| 192 | +int main() { |
| 193 | + cudaDeviceProp prop; cudaGetDeviceProperties(&prop, 0); |
| 194 | + printf("device: %s (sm_%d%d)\n", prop.name, prop.major, prop.minor); |
| 195 | + for (int np : {64, 1024, 8192}) { |
| 196 | + NPROB = np; |
| 197 | + run_dtype<float>("f32"); |
| 198 | + run_dtype<double>("f64"); |
| 199 | + } |
| 200 | + return 0; |
| 201 | +} |
0 commit comments