Blog

CUDA Deep Dive: 10B Galaxy Pairs on an RTX 5080

A Blackwell port of the galaxy-correlation solver, rebuilt as an ablation ladder: six binaries, one optimization apart. Two did not work.

21 August 2026 · HPC, GPU, CUDA, NVIDIA, Linux, Performance

An earlier post described the HIP version of this solver, which was tuned for RDNA 2 on an AMD RX 6900 XT. This post describes the NVIDIA side of the same problem: a native CUDA implementation written for the RTX 5080, which is a Blackwell part compiled here for sm_120.

The final runtime is not really the point of the exercise. What I wanted was the optimization history rebuilt as a measured ablation ladder, meaning six binaries that differ from one another by exactly one optimization, so that every claim made below has a measurement behind it rather than an argument. Two of the optimizations turned out not to work, and they are reported here alongside the ones that did.

The dataset itself can be inspected in a small 3D viewer: Galaxy Visualization.

What this is, and what it is not #

It is worth being clear about the goal before any numbers appear, because it determines which trade-offs count as acceptable.

This is an optimization exercise, not a scientific measurement. The aim is to take one well-defined problem and push it as hard as two different pieces of consumer hardware will go, first AMD RDNA 2 and now NVIDIA Blackwell. What matters here is total runtime and correctness: the program has to produce histograms that sum to exactly $N^2$ , and it has to do it as fast as the hardware allows. Those are the two criteria, and everything below is judged against them.

What is explicitly not being claimed is bit-exact numerical fidelity. Several of the steps in this post change the last few digits of the result, and one of them moves counts between neighbouring bins. That is a normal consequence of -ffast-math, of reassociating floating-point expressions, and of replacing a library function with a polynomial approximation. Where it happens it is measured and reported rather than glossed over, so that anyone reusing this code knows what they are getting. For studying the shape of a correlation function it is irrelevant. For quoting a figure off an individual bin it would not be, and such a use would call for the exact acosf path.

The problem #

The input consists of two catalogs of 100,000 galaxies each. One catalog holds real observations and the other a synthetic random distribution, and both are given as a right ascension $\alpha$ and a declination $\delta$ expressed in arcminutes.

For every pair of galaxies the angular separation has to be computed:

\theta_{12} = \arccos\left(\sin(\delta_1)\sin(\delta_2) + \cos(\delta_1)\cos(\delta_2)\cos(\alpha_1 - \alpha_2)\right)

The separations are binned into 360 * 4 = 1440 buckets of 0.25° each, for three histogram families, which are the real-real pairs (DD), the real-random pairs (DR) and the random-random pairs (RR). The three histograms are then combined with the Landy-Szalay estimator:

\omega = \frac{DD - 2DR + RR}{RR}

Since each family covers every pair, there are $N^2 = 10^{10}$ pairs per family. The correctness criterion is exact rather than statistical: each histogram has to sum to precisely 10000000000. Every run reported in this post asserts all three sums, and a variant that failed the assert would have been reported as a failure rather than timed.

One property of this particular dataset is that only 360 of the 1440 bins are ever populated, because the largest separation present falls in the bin starting at 89.75°.

Test rig #

All the measurements below were taken on a single machine, whose relevant configuration is given in the table:

Component	Value
GPU	NVIDIA GeForce RTX 5080 (GB203), compute capability 12.0
SMs	84, max 1536 threads/SM, 65536 regs/block
Memory	16.65 GB, 64 MB L2, 49152 B shared memory per block
Clock	2700 MHz (as reported by `cudaDevAttrClockRate`)
Driver	610.57.04 (CUDA UMD 13.3)
Toolkit	nvcc 13.3, V13.3.73
CPU	AMD Ryzen 9 3950X (16C/32T)
Kernel	7.2.0-1-cachyos
Governor	`performance`
Storage	Samsung 970 EVO Plus NVMe, LUKS, btrfs

The build flags are the ones used by the project’s own build script:

nvcc -O3 --use_fast_math -arch=sm_120 -lineinfo --ptxas-options=-v -DTILE=512

Measurement protocol #

Every variant is built from one LEVEL-guarded source file, so that the host code is byte-identical across the whole ladder and only the ablated feature differs between binaries.
Each variant was run 9 times. The first run is reported separately as the cold run, and the steady state figure is the median of runs 2 through 9, with the minimum quoted next to it.
Kernel time is measured with CUDA events rather than with gettimeofday, so that host-side overhead is excluded from it.
The desktop is idle, but a Wayland compositor (kwin_wayland) is live on the same GPU. With the CPU governor set to performance and nothing else using the GPU, most runs land within 0.02 ms of the median, with occasional excursions about 0.8 ms above it. Any difference of that order is therefore given a separate interleaved A/B, in which the two binaries alternate run by run, so that drift affects both arms equally.
Nsight Compute is not installed on this machine and RmProfilingAdminOnly is set to 1, so there are no hardware counters in this post. Occupancy is not measured; it is inferred from the ptxas register counts and the shared-memory footprint.

Before any of these measurements were counted, the scratchpad build was checked against the repository’s committed results/omega.out, which it reproduces byte for byte. That check is make verify in the source repository, so it can be repeated rather than taken on trust.

One caveat is worth stating explicitly. An earlier pass of the same suite, run with the governor set to powersave and with video playing on the same desktop, came out 3 to 6% slower across the whole ladder. That varies two things at once and is therefore not a controlled experiment, so it is not published here as a governor result. It is, however, a reasonable argument for fixing the state of the machine before benchmarking anything on it.

The ladder #

Level	What it adds	Kernel median	Kernel min	Registers	Speedup vs. previous
L0	one thread per $(i,j)$ , 32×32 blocks, per-pair `sincos`, `acosf`, no symmetry	188.653 ms	188.312 ms	32	-
L1	+ DD/RR symmetry	169.209 ms	168.696 ms	37	1.12×
L2	+ host-precomputed $\sin/\cos(\delta)$	163.435 ms	162.939 ms	40	1.04×
L3	+ block tiling with shared-memory $j$ -tile	34.027 ms	33.951 ms	30	4.80×
L4	+ `fast_acosf` polynomial	24.528 ms	24.514 ms	30	1.39×
L5	+ unroll, `__launch_bounds__`, diagonal block skip	22.747 ms	22.742 ms	34	1.08×

The ladder gives 8.29× end to end, and no level spills registers; ptxas reports 0 bytes spill stores, 0 bytes spill loads throughout. In wall-clock terms the program as a whole goes from 0.195 s to 0.032 s.

At L5 the kernel evaluates $10^{10}$ DR angles together with $N(N{+}1)/2$ angles each for DD and RR, which is 20,000,100,000 angle evaluations in 22.747 ms, or approximately 879 G pair-angles/s. L0, which has no symmetry, evaluates all 30,000,000,000 angles in 188.653 ms, which corresponds to 159 G/s.

L1: DD/RR symmetry #

DD and RR are symmetric in $(i, j)$ , so only the upper triangle needs to be computed, with the off-diagonal pairs weighted twice:

if (j >= i) {
  const unsigned int inc = 1U + static_cast<unsigned int>(j != i);
  hist_add(s_hist_DD, compute_histogram_index(dd_expr), inc);
  hist_add(s_hist_RR, compute_histogram_index(rr_expr), inc);
}

The weighting is exact rather than an approximation, since

\sum_{i=0}^{N-1}\sum_{j=i}^{N-1} \left(1 + [j \neq i]\right) = N + 2\cdot\frac{N(N-1)}{2} = N^2

Removing the redundant half of two of the three histogram families halves two thirds of the arithmetic, and it is worth 10.3%, from 188.653 ms to 169.209 ms. The size of that gap is the first useful signal in the ladder, because it indicates that the naive kernel is not arithmetic-bound. A thread in a block below the diagonal is still launched, still issues its loads and still evaluates the branch; only the arithmetic is skipped. The work is removed without the thread being removed, and the measured saving is correspondingly smaller than the arithmetic saving.

L2: precomputed sin/cos #

The kernel up to this point recomputes $\sin(\delta)$ and $\cos(\delta)$ for both galaxies of every pair, which amounts to eight transcendental evaluations per pair, repeated $10^{10}$ times. Those values depend only on the galaxy and not on the pair, so they can be hoisted into a host-side loop over the 100,000 elements of each catalog:

for (long int k = 0; k < N; ++k) {
  real_sin[k] = sinf(real_decl[k]);
  real_cos[k] = cosf(real_decl[k]);
  rand_sin[k] = sinf(rand_decl[k]);
  rand_cos[k] = cosf(rand_decl[k]);
}

The precomputation is worth 3.5% on the kernel, from 169.209 ms to 163.435 ms. It costs approximately 2.3 ms of host time, since the input phase moves from a median of 4.65 ms to one of 6.99 ms. Trading 2.3 ms of host time for 5.8 ms of kernel time is favourable wherever the kernel dominates the runtime, which it does here.

The register count is worth following through the ladder. It climbs from 32 at L0 to 37 at L1 and 40 at L2, which is to say that the first two optimizations buy speed with registers. At 40 registers per thread the kernel is approaching the point where occupancy would begin to suffer for it.

L3: block tiling, the one that mattered #

This is the largest single step in the ladder, at 4.80×, from 163.435 ms to 34.027 ms.

In the one-thread-per-pair arrangement there are $10^{10}$ threads, and each of them performs six global memory loads in order to compute a single angle. In the tiled arrangement each thread instead owns one $i$ galaxy, keeps its values in registers for the lifetime of the block, and iterates over a TILE-wide tile of $j$ galaxies that the block has cooperatively staged in shared memory:

const int tid = threadIdx.x;
const int i = blockIdx.x * TILE + tid;

// i-galaxy lives in registers for the whole block
const float ri_sin = d_real_sin[li];
const float ri_cos = d_real_cos[li];
const float ri_ra  = d_real_rasc[li];

// the j-tile is loaded cooperatively, once, and reused TILE times
if (tid < j_count) {
  sj_real_sin[tid] = d_real_sin[jload];
  sj_real_cos[tid] = d_real_cos[jload];
  sj_real_ra[tid]  = d_real_rasc[jload];
  // ... and the rand catalog
}
__syncthreads();

for (int jj = 0; jj < j_count; ++jj) { /* TILE pairs per loaded tile */ }

Six global loads per pair therefore become six global loads per TILE pairs. The register pressure also falls rather than rises, from 40 to 30, because the inner loop reads its $j$ values out of shared memory instead of keeping eight per-pair values live.

Two properties of this step are worth noting. It is the only step in the ladder that changes the memory access pattern rather than the arithmetic, and it is by a wide margin the largest, which together indicate what the kernel was bound by from the beginning. It is also numerically free: the tiled output is bit-identical to the output of L2, with 0 of the 360 populated bins differing. The 4.80× is obtained without moving a single count.

L4: fast acos, with an asterisk #

Replacing acosf with a minimax polynomial approximation is the second largest step, at 1.39×, from 34.027 ms to 24.528 ms:

__device__ __forceinline__ float fast_acosf(float x) {
  float negate = (float)(x < 0.0f);
  x = fabsf(x);
  float ret = -0.0187293f;
  ret = ret * x + 0.0742610f;
  ret = ret * x - 0.2121144f;
  ret = ret * x + 1.5707288f;
  ret = ret * sqrtf(1.0f - x);
  ret = ret - 2.0f * negate * ret;
  return negate * 3.14159265358979f + ret;
}

The comment above this function in the source claimed that the approximation error stays “well below the 0.25-degree bin width, so bin assignment matches acosf exactly here”. That claim is not correct, and the diff between the L3 and L4 outputs is what established it. The measured differences are the following:

353 of 360 populated bins differ
largest absolute shift in $\omega$ : 0.002257, in the bin starting at 89.75°
largest relative shift in $\omega$ : 17.19%, at 33.25° where $\omega$ is 0.000739, so that percentage is mostly a near-zero denominator talking
largest single-bin count movement: 12,332 counts of DD at 45.25°, which is 0.03% of that bin

The histogram sums are nevertheless still exactly $10^{10}$ , so the correctness assert does not fire. Counts near a bin boundary simply land in the neighbouring bin, and a redistribution of that kind is invisible to an assert on totals.

Measured against the criteria set out at the top, this is an acceptable trade. Runtime falls by 28% and every histogram still sums to exactly $10^{10}$ , which is what the exercise is optimizing for. The problem was never the approximation itself but the claim made about it: a comment that promises exact bin agreement invites the reader to skip the check, and the check is the only thing that would have caught this. The comment in the source has since been corrected to state the measured behaviour instead.

For calibration it is worth noting that the L2 precompute step is not bit-exact either. Reassociating the trigonometry changes the rounding, and all 360 bins differ, with a largest absolute $\omega$ shift of 0.008749 at the bin starting at 0.0° (0.37% of that bin’s value) and a largest relative shift of 0.60% at 31.0°, where $\omega$ is -0.0005. The fast acos is therefore not even the largest numerical perturbation in the ladder. It is only the one whose source comment promised otherwise.

L5: the last 8% #

This level combines three small changes: a #pragma unroll 8 on the inner loop, a __launch_bounds__(TILE) annotation on the kernel, and a skip of the DD and RR work for blocks that lie entirely below the diagonal:

// Whole block below the diagonal -> no j >= i anywhere -> skip DD/RR.
const bool do_sym = (j_base + j_count - 1) >= i_base;

The last of the three is the per-block form of the predicate that L1 introduced per thread, and it is the thing L1 could not do, since a per-thread predicate does not prevent the block from being scheduled in the first place. The difference between L4 and L5 is 1.8 ms, which is close enough to the noise floor of this machine to warrant an interleaved A/B of 20 alternating runs:

Variant	n	min	median	mean	max
L4	20	24.518 ms	24.534 ms	24.937 ms	25.942 ms
L5	20	22.740 ms	22.766 ms	23.070 ms	24.065 ms

The two distributions are completely separated, in that the worst L5 run is faster than the best L4 run. The register count rises from 30 to 34, which is the cost of the unroll, and there are still no spills. The output is bit-identical to that of L4.

Tile size sweep #

TILE sets both the block size and the width of the shared-memory tile, so changing it moves occupancy and data reuse at the same time. Five values were measured back to back:

TILE	Threads/block	Grid	Shared mem/block	Registers	Kernel median
64	64	1563²	19968 B	46	58.050 ms
128	128	782²	21504 B	34	30.857 ms
256	256	391²	24576 B	34	23.294 ms
512	512	196²	30720 B	34	22.749 ms
1024	1024	98²	43008 B	41	23.105 ms

Below 256 the fixed per-block cost dominates. Every block zeroes 3 × 1536 shared-memory bins and flushes them to global memory regardless of how many pairs it processes, and at TILE=64 that overhead is amortized over only 64 pairs per loaded tile. At the other end, TILE=1024 brings the shared-memory footprint to 43008 B against a per-block limit of 49152 B and raises the register count to 41, at which point the curve turns back upwards.

The best three values lie within 2.4% of one another, so 512 and 256 were compared directly in an interleaved A/B of 20 runs:

Variant	n	min	median	mean	max
TILE=512	20	22.736 ms	22.762 ms	23.046 ms	24.312 ms
TILE=256	20	23.275 ms	23.289 ms	23.367 ms	24.109 ms

TILE=512 is faster by 2.3%. The margin is small, but it is consistent, and the two distributions barely overlap.

The Linux side #

At this point the kernel accounts for approximately two thirds of the wall clock. The remainder is input and output, which is where the remaining savings are.

mmap + hand-rolled parser vs. fscanf #

The catalogs are stored as ASCII text: a count line followed by 100,000 rows of two floating-point values each. The straightforward way to read that is an fscanf loop, and it is very slow. Replacing it with mmap and a hand-written ASCII float parser gives the following:

int fd = open(file_path, O_RDONLY);
struct stat file_info;
fstat(fd, &file_info);
void *mapped_data = mmap(NULL, (size_t)file_info.st_size, PROT_READ, MAP_PRIVATE, fd, 0);
close(fd);

const char *cursor = (const char *)mapped_data;
const char *end = cursor + file_info.st_size;
// parse_int_fast / parse_float_fast walk the buffer directly

Interleaved A/B on the input phase, with a warm page cache and 15 runs each:

Reader	n	min	median	mean	max
`fscanf`	15	58.230 ms	59.310 ms	59.621 ms	61.414 ms
`mmap` + parser	15	6.314 ms	6.586 ms	6.735 ms	7.367 ms

That is 9.0× on the input phase. For comparison, the fscanf reader on its own costs more than 2.5 times the entire optimized GPU kernel. Once the kernel has been brought down to 23 ms, 59 ms spent interpreting scanf format strings is not a rounding error; it is the largest single item in the program.

The output is bit-identical as well. The hand-written parser accumulates in double before narrowing to float, and across the 200,000 parsed values it does not move a single histogram count relative to strtof.

Cold vs. warm page cache #

The gap between the first run and the rest has a specific cause, and it is worth measuring rather than attributing to warm-up in general terms. Clean page cache can be evicted without root privileges through posix_fadvise:

int fd = open(path, O_RDONLY);
posix_fadvise(fd, 0, 0, POSIX_FADV_DONTNEED);

The eviction was verified with mincore, which reports residency dropping from 409/409 and 601/601 pages to 0/0. Ten cold and warm pairs were then run alternately:

Reader	Cache	n	min	median	mean	max
`mmap`	warm	10	6.595 ms	7.242 ms	7.114 ms	7.453 ms
`mmap`	cold	10	10.323 ms	16.207 ms	14.840 ms	17.253 ms
`fscanf`	warm	10	58.680 ms	60.075 ms	60.372 ms	61.794 ms
`fscanf`	cold	10	61.875 ms	67.327 ms	66.826 ms	69.964 ms

The cold penalty is 9.0 ms for mmap and 7.3 ms for fscanf, and in both cases the cold distributions are wide, with a spread of approximately 7 ms between the minimum and the maximum. This is the cost of faulting 4 MB off an NVMe device through LUKS and btrfs, and it is largely independent of which parser runs afterwards.

It is also the reason the protocol above discards the first run. What that run measures is the page cache, not GPU warm-up.

Pinned host memory: no #

The standard recommendation for host-to-device transfers is to allocate the host buffers with cudaHostAlloc rather than calloc, so that the driver can DMA directly out of non-pageable memory. Six arrays are uploaded here, totalling approximately 2.4 MB, so the expectation was a small improvement.

Interleaved A/B, 15 runs each:

Buffers	Metric	n	min	median	mean	max
pageable	GPU phase	15	23.263 ms	23.301 ms	23.481 ms	24.093 ms
pinned	GPU phase	15	23.190 ms	23.211 ms	23.859 ms	29.916 ms
pageable	wall clock	15	0.031 s	0.031 s	0.031 s	0.032 s
pinned	wall clock	15	0.032 s	0.033 s	0.033 s	0.034 s

The GPU phase is unchanged for practical purposes. The pinned variant is 0.09 ms better on the median and 0.4 ms worse on the mean, and it produced a 29.9 ms outlier that the pageable variant never produced. The wall clock, however, is consistently 2 ms worse, because cudaHostAlloc is an expensive allocator and pinning the eight 400 KB host arrays costs more than the transfer saves at this size.

Pinned memory is the right answer when hundreds of megabytes are being streamed, or when copies are being overlapped with compute. For a single 2.4 MB upload it is a net loss, and that conclusion is only available because the wall clock was measured and not only the transfer.

What didn’t pay off #

Two of the things tried here did not pay off, and both are reported, because an ablation ladder in which every rung is a win has usually been edited after the fact.

1. Shared-memory bin padding. The histograms are 1440 bins padded to 1536, on the assumption that a power-of-two-aligned stride avoids bank conflicts across the 32 shared-memory banks. Interleaved A/B, 15 runs each:

Padding	n	min	median	mean	max
1536 (padded)	15	22.734 ms	22.754 ms	22.893 ms	23.493 ms
1440 (unpadded)	15	22.730 ms	22.752 ms	22.839 ms	23.500 ms

Every statistic agrees to within 0.06 ms, which is 0.3%. There is no measurable benefit on this GPU, and the padding costs 288 bytes of shared memory per block in exchange for it.

This is not an argument that padding is useless in general. It came over from the RDNA 2 tuning pass, where the LDS bank behaviour is different, and that side has not been re-measured here. On Blackwell, with this access pattern, the effect is simply not present. The atomics are scattered across bins by the data rather than by the thread index, so the conflict pattern that padding is designed to remove is not the pattern this kernel generates.

2. Pinned host memory, as described above. It is negative on wall clock rather than merely neutral.

Correctness guardrails #

Every run asserts all three sums:

if (histogramDRsum != 10000000000L) { /* bail */ }
if (histogramDDsum != 10000000000L) { /* bail */ }
if (histogramRRsum != 10000000000L) { /* bail */ }

This is a useful check. It catches races, dropped atomics and an off-by-one in the symmetry weighting immediately, because each of those changes the total.

The fast_acosf result above is a clear demonstration of its blind spot, however. Moving 12,332 counts from one bin into its neighbour does not change the sum at all, since a total is invariant under redistribution. The sum assert therefore has to be paired with a diff of the output against a known-good reference, which is what identified the discrepancy here, and which is the reason the first step in this work was to confirm that the scratchpad build reproduces the committed omega.out byte for byte.

Lessons #

Ablation is more informative than narration. Six binaries that differ by one change each turn a description of what was optimized into a table in which the 4.80× step is unmistakable and the 1.04× step is honestly labelled as 1.04×.
Memory layout mattered more than arithmetic, by a wide margin. Block tiling was worth more than the symmetry, the precomputation and the fast acos combined, and unlike those three it changed nothing numerically.
A correctness check has a blind spot, and the blind spot needs its own check. The sum asserts caught nothing about fast_acosf, whereas a diff against a reference output caught it in a single command.
A/B measurements should be interleaved. Even on an otherwise quiet machine the last two rungs of this ladder are within a couple of milliseconds of each other, and sequential measurement will manufacture or hide a difference of that size.
Input and output are part of the program. The fscanf reader was costing 2.5 times the entire optimized kernel.
Textbook optimizations are hypotheses rather than conclusions. Pinned memory and bin padding are both standard advice, and both measured neutral to negative here. The advice is not wrong; it is conditional, and the condition is the workload.

The whole history, Dione to Blackwell #

The table below collects every recorded stage of this program, from the original course implementation through the AMD tuning work to the current CUDA build. The AMD figures are the ones recorded in the project’s own README at each commit, on an RX 6900 XT. They have not been re-measured for this post, and they cannot be: there is no AMD GPU in the machine used here. The final row is the only one measured under the protocol described above.

Date	Stage	Hardware	Wall clock	Kernel
Course project	Original implementation	Dione cluster	8.7 s	not recorded
2025-06-19	CUDA to HIP port	RX 6900 XT	1.6 s	not recorded
2025-06-19	+ shared-memory histograms	RX 6900 XT	0.85 s	not recorded
(intermediate)	recorded as “previous record”	RX 6900 XT	0.6 s	not recorded
2026-02-03	Native HIP, RDNA 2 tuning	RX 6900 XT	0.33 s	274.685 ms
2026-02-05	+ further kernel tuning	RX 6900 XT	not updated	194.707 ms
2026-02-14	+ block-size tuning (32×32)	RX 6900 XT	0.21 s	145.535 ms
2026-02-14	+ `mmap` input, DD/RR symmetry	RX 6900 XT	0.152 s	114.500 ms
2026-02-14	Final native HIP	RX 6900 XT	0.154 s (best 0.152)	116.330 ms
2026-08-21	Native CUDA, Blackwell	RTX 5080	0.032 s	22.75 ms

Two cautions about reading this table. The hardware changes at the last row, so the step from 0.154 s to 0.032 s is not a software result: it is a different GPU, three years newer, running a differently structured kernel. Nothing here says CUDA is faster than HIP, or that NVIDIA is faster than AMD, because no such comparison was run. The only like-for-like measurement in this post is the ablation ladder, where every rung uses the same GPU, the same compiler and the same host code.

The second caution is that the AMD rows are historical records rather than measurements taken under a stated protocol. They come from different days, different driver versions and, in at least one case, a README that was not updated when the kernel time was. They are included because the question “where did this start” is a reasonable one, not because they would survive the scrutiny the last row was put through.

With that said, the trajectory is the point. The original course implementation took 8.7 seconds. The current one takes 0.032 seconds and still asserts, on every single run, that all three histograms sum to exactly ten billion. Total runtime and correctness were the two criteria stated at the top, and both of them held all the way down.

Appendix A: `src/galaxy_cuda.cu` #

The complete implementation is reproduced below, at the revision the measurements in this post were taken from. It is a single translation unit of 751 lines: the device kernel and its helpers first, then the memory-mapped catalog reader, and finally the host driver that allocates, launches, verifies the histogram sums and writes omega.out.

Three comments in this listing were corrected while the measurements were being taken, because the source claimed things the numbers did not support. The most significant is the one above fast_acosf, which previously asserted that the approximation left bin assignment identical to acosf. It does not, and the section on the fast acos above gives the measured extent of the difference.

One line of code was also corrected, after clang-tidy was pointed at the sources while the repository was being prepared for release. The dim3 members are unsigned int, so number_of_threads evaluated its whole product in 32 bits and only widened afterwards. It does not overflow at the tile size used here, but it did on the AMD configuration, where 3125 x 3125 x 1024 threads wrapped to the 1410065408 that the project’s README reported for years. The value is only ever printed, so no measurement in this post is affected.

A further round of corrections was made when the source was published as a standalone repository, and the listing below reflects them. A comment still pointed at a build script that had been deleted; parseargs_readinput was declared at block scope inside main; it and get_device were declared static but defined without it; a variable named threadblocks held a thread count; and the vestigial struct timezone argument to gettimeofday, which glibc ignores, was dropped. None of it touches the kernel: the binary still compiles to 34 registers with zero spills and reproduces results/omega.out byte for byte. That is what the four added lines are.

The build command is the one given under Test rig, with TILE set to 512.

  1// Native CUDA implementation tuned for NVIDIA RTX 5080 (Blackwell, sm_120).
  2// Same two-point angular correlation as galaxy_hip.cpp, restructured for
  3// NVIDIA:
  4//   - block-tiled pair loop (each thread owns one i, loops a shared-memory
  5//     j-tile) instead of one-thread-per-pair -> massive reuse / fewer loads
  6//   - precomputed per-galaxy sin/cos(decl) (no per-pair sincos)
  7//   - fast polynomial acos approximation (measured |err| <= 6.77e-5 rad;
  8//     this does shift a small number of counts between bins, see fast_acosf)
  9//   - DD/RR symmetry with whole-block diagonal skipping
 10//   - inner-loop unrolling; shared-memory atomic histograms (fast on Blackwell)
 11// Measured on an RTX 5080 (driver 610.57.04, CUDA 13.3, sm_120, performance
 12// governor), median of runs 2-9: 22.7 ms kernel, 0.032 s wall clock. The same
 13// code with the block tiling, the fast acos, the symmetry and the unrolling
 14// removed one at a time measures 188.7 ms, so the tuning is worth 8.3x. See
 15// the Makefile (`make cuda`, `make bench`) for the build and run protocol.
 16//
 17// The host side (catalog reader, timing, histogram assertions) is duplicated
 18// in galaxy_hip.cpp rather than shared through a header. That is deliberate:
 19// each backend is one self-contained translation unit that compiles with a
 20// single command and can be read end to end. The two share no device code.
 21#include <cuda_runtime.h>
 22#include <fcntl.h>
 23#include <inttypes.h>
 24#include <math.h>
 25#include <stdio.h>
 26#include <stdlib.h>
 27#include <sys/mman.h>
 28#include <sys/stat.h>
 29#include <sys/time.h>
 30#include <unistd.h>
 31
 32static float *real_rasc;
 33static float *real_decl;
 34static float *rand_rasc;
 35static float *rand_decl;
 36static float *real_sin; // precomputed sin(decl) for real catalog
 37static float *real_cos; // precomputed cos(decl) for real catalog
 38static float *rand_sin; // precomputed sin(decl) for rand catalog
 39static float *rand_cos; // precomputed cos(decl) for rand catalog
 40static constexpr long int N = 100000L;
 41static long *histogram_DR;
 42static long *histogram_DD;
 43static long *histogram_RR;
 44static constexpr float PI = 3.14159265358979323846f;
 45static long int CPUMemory = 0L;
 46static long int GPUMemory = 0L;
 47static constexpr int totaldegrees = 360;
 48static constexpr int binsperdegree = 4;
 49
 50// 1440 bins padded to 1536, carried over from the RDNA 2 tuning pass where the
 51// intent was to avoid LDS bank conflicts. Note 1536 is 3 * 512, not a power of
 52// two. On this GPU the padding measures as a no-op: an interleaved A/B of 1536
 53// against 1440 agrees to within 0.06 ms (0.3%) on every statistic, because the
 54// atomics scatter across bins by data rather than by thread index. It is kept
 55// only so the output stays comparable with the HIP build; 1440 is equally fine.
 56const int num_bins = binsperdegree * totaldegrees; // 1440
 57const int num_bins_padded = 1536;
 58
 59// Tile size = threads per block. Each block computes a TILE x TILE sub-block of
 60// the pair matrix: TILE threads each own one i and loop a shared-memory j-tile.
 61#ifndef TILE
 62#define TILE 512
 63#endif
 64
 65#define CUDA_ERR_CHECK(ans)                                                    \
 66  {                                                                            \
 67    gpuAssert((ans), __FILE__, __LINE__);                                      \
 68  }
 69static inline void gpuAssert(cudaError_t code, const char *file, int line,
 70                             bool abort = true) {
 71  if (code != cudaSuccess) {
 72    fprintf(stderr, "   GPUassert: %s %s %d\n", cudaGetErrorString(code), file,
 73            line);
 74    if (abort)
 75      exit(code);
 76  }
 77}
 78
 79__device__ __forceinline__ void hist_add(unsigned int *histogram, int bin_index,
 80                                         unsigned int increment = 1U) {
 81  atomicAdd(&histogram[bin_index], increment);
 82}
 83
 84__device__ __forceinline__ float fast_acosf(float x) {
 85  // Handbook-of-Math-Functions minimax approximation. Measured maximum error
 86  // against acosf is 6.77e-5 rad (0.0039 degrees) over a 2e8-point sweep of
 87  // [-1, 1], which is far below the 0.25-degree bin width.
 88  //
 89  // That bound does NOT make the binning identical, and an earlier version of
 90  // this comment wrongly claimed it did. A sample only has to sit within
 91  // 0.0039 degrees of a bin edge to move, and 0.65% of the sweep does exactly
 92  // that. Over the real catalogs it shifts counts in 353 of the 360 populated
 93  // bins, at most 0.03% of any one bin, with a largest omega change of
 94  // 0.002257 in the bin at 89.75 degrees. The histogram totals are unaffected,
 95  // so the N*N asserts below cannot see this: a total is invariant to
 96  // redistribution. Use acosf instead if you need exact bin agreement.
 97  float negate = (float)(x < 0.0f);
 98  x = fabsf(x);
 99  float ret = -0.0187293f;
100  ret = ret * x + 0.0742610f;
101  ret = ret * x - 0.2121144f;
102  ret = ret * x + 1.5707288f;
103  ret = ret * sqrtf(1.0f - x);
104  ret = ret - 2.0f * negate * ret;
105  return negate * 3.14159265358979f + ret;
106}
107
108__device__ __forceinline__ int compute_histogram_index(float expr) {
109  constexpr float angle = 57.29577951308232f;
110  expr = fminf(fmaxf(expr, -1.0f), 1.0f);
111
112  int histogram_index = int(fast_acosf(expr) * angle * binsperdegree);
113  return (histogram_index < 0)
114             ? 0
115             : ((histogram_index >= num_bins) ? num_bins - 1 : histogram_index);
116}
117
118// Block-tiled kernel with precomputed sin/cos(decl).
119// Grid is 2D over TILE-sized i-blocks (x) and j-blocks (y). Each of the TILE
120// threads owns one i; it keeps its real_i/rand_i data in registers and loops a
121// shared-memory tile of TILE j-galaxies. DD/RR use symmetry (j >= i); blocks
122// fully below the diagonal skip DD/RR entirely.
123__global__ void __launch_bounds__(TILE) fill_histograms(
124    const float *__restrict__ d_real_rasc, const float *__restrict__ d_real_sin,
125    const float *__restrict__ d_real_cos, const float *__restrict__ d_rand_rasc,
126    const float *__restrict__ d_rand_sin, const float *__restrict__ d_rand_cos,
127    unsigned long long int *d_histogram_DR,
128    unsigned long long int *d_histogram_DD,
129    unsigned long long int *d_histogram_RR) {
130  const int tid = threadIdx.x;
131  const int i_base = blockIdx.x * TILE;
132  const int j_base = blockIdx.y * TILE;
133  const int i = i_base + tid;
134
135  extern __shared__ unsigned int s_mem[];
136  unsigned int *s_hist_DR = s_mem;
137  unsigned int *s_hist_DD = s_hist_DR + num_bins_padded;
138  unsigned int *s_hist_RR = s_hist_DD + num_bins_padded;
139  // Coordinate tiles for the j-block (real and rand catalogs).
140  float *s_coord = (float *)(s_hist_RR + num_bins_padded);
141  float *sj_real_sin = s_coord;
142  float *sj_real_cos = sj_real_sin + TILE;
143  float *sj_real_ra = sj_real_cos + TILE;
144  float *sj_rand_sin = sj_real_ra + TILE;
145  float *sj_rand_cos = sj_rand_sin + TILE;
146  float *sj_rand_ra = sj_rand_cos + TILE;
147
148  for (int b = tid; b < num_bins_padded; b += TILE) {
149    s_hist_DR[b] = 0;
150    s_hist_DD[b] = 0;
151    s_hist_RR[b] = 0;
152  }
153
154  // Load this thread's i-galaxy data into registers.
155  const bool valid_i = i < N;
156  const int li = valid_i ? i : 0;
157  const float ri_sin = d_real_sin[li];
158  const float ri_cos = d_real_cos[li];
159  const float ri_ra = d_real_rasc[li];
160  const float di_sin = d_rand_sin[li];
161  const float di_cos = d_rand_cos[li];
162  const float di_ra = d_rand_rasc[li];
163
164  const int j_end = min(j_base + TILE, (int)N);
165  const int j_count = j_end - j_base;
166  // Whole block below the diagonal -> no j >= i anywhere -> skip DD/RR.
167  const bool do_sym = (j_base + j_count - 1) >= i_base;
168
169  // Cooperatively load the j-tile into shared memory.
170  const int jload = j_base + tid;
171  if (tid < j_count) {
172    sj_real_sin[tid] = d_real_sin[jload];
173    sj_real_cos[tid] = d_real_cos[jload];
174    sj_real_ra[tid] = d_real_rasc[jload];
175    sj_rand_sin[tid] = d_rand_sin[jload];
176    sj_rand_cos[tid] = d_rand_cos[jload];
177    sj_rand_ra[tid] = d_rand_rasc[jload];
178  }
179  __syncthreads();
180
181  if (valid_i) {
182#pragma unroll 8
183    for (int jj = 0; jj < j_count; ++jj) {
184      const int j = j_base + jj;
185
186      // DR: real_i vs rand_j (not symmetric, always counted).
187      const float dr_expr =
188          ri_sin * sj_rand_sin[jj] +
189          ri_cos * sj_rand_cos[jj] * cosf(ri_ra - sj_rand_ra[jj]);
190      hist_add(s_hist_DR, compute_histogram_index(dr_expr));
191
192      // DD/RR: symmetric, only j >= i.
193      if (do_sym && j >= i) {
194        const unsigned int inc = 1U + static_cast<unsigned int>(j != i);
195
196        const float dd_expr =
197            ri_sin * sj_real_sin[jj] +
198            ri_cos * sj_real_cos[jj] * cosf(ri_ra - sj_real_ra[jj]);
199        hist_add(s_hist_DD, compute_histogram_index(dd_expr), inc);
200
201        const float rr_expr =
202            di_sin * sj_rand_sin[jj] +
203            di_cos * sj_rand_cos[jj] * cosf(di_ra - sj_rand_ra[jj]);
204        hist_add(s_hist_RR, compute_histogram_index(rr_expr), inc);
205      }
206    }
207  }
208  __syncthreads();
209
210  for (int b = tid; b < num_bins; b += TILE) {
211    if (s_hist_DR[b] > 0)
212      atomicAdd(&d_histogram_DR[b], (unsigned long long int)s_hist_DR[b]);
213    if (s_hist_DD[b] > 0)
214      atomicAdd(&d_histogram_DD[b], (unsigned long long int)s_hist_DD[b]);
215    if (s_hist_RR[b] > 0)
216      atomicAdd(&d_histogram_RR[b], (unsigned long long int)s_hist_RR[b]);
217  }
218}
219
220// Forward declarations
221static int get_device();
222static int parseargs_readinput(int argc, char *argv[]);
223
224static inline bool is_ascii_whitespace(char c) {
225  return (c == ' ' || c == '\n' || c == '\r' || c == '\t' || c == '\v' ||
226          c == '\f');
227}
228
229static inline void skip_ascii_whitespace(const char *&cursor, const char *end) {
230  while (cursor < end && is_ascii_whitespace(*cursor))
231    ++cursor;
232}
233
234static bool parse_int_fast(const char *&cursor, const char *end, int *value) {
235  skip_ascii_whitespace(cursor, end);
236  if (cursor >= end)
237    return false;
238
239  int sign = 1;
240  if (*cursor == '+' || *cursor == '-') {
241    sign = (*cursor == '-') ? -1 : 1;
242    ++cursor;
243  }
244
245  if (cursor >= end || *cursor < '0' || *cursor > '9')
246    return false;
247
248  int parsed_value = 0;
249  while (cursor < end && *cursor >= '0' && *cursor <= '9') {
250    parsed_value = parsed_value * 10 + (*cursor - '0');
251    ++cursor;
252  }
253
254  *value = sign * parsed_value;
255  return true;
256}
257
258static bool parse_float_fast(const char *&cursor, const char *end,
259                             float *value) {
260  skip_ascii_whitespace(cursor, end);
261  if (cursor >= end)
262    return false;
263
264  int sign = 1;
265  if (*cursor == '+' || *cursor == '-') {
266    sign = (*cursor == '-') ? -1 : 1;
267    ++cursor;
268  }
269
270  double result = 0.0;
271  bool has_digits = false;
272
273  while (cursor < end && *cursor >= '0' && *cursor <= '9') {
274    has_digits = true;
275    result = result * 10.0 + (double)(*cursor - '0');
276    ++cursor;
277  }
278
279  if (cursor < end && *cursor == '.') {
280    ++cursor;
281    double place = 0.1;
282    while (cursor < end && *cursor >= '0' && *cursor <= '9') {
283      has_digits = true;
284      result += (double)(*cursor - '0') * place;
285      place *= 0.1;
286      ++cursor;
287    }
288  }
289
290  if (!has_digits)
291    return false;
292
293  if (cursor < end && (*cursor == 'e' || *cursor == 'E')) {
294    ++cursor;
295
296    int exp_sign = 1;
297    if (cursor < end && (*cursor == '+' || *cursor == '-')) {
298      exp_sign = (*cursor == '-') ? -1 : 1;
299      ++cursor;
300    }
301
302    if (cursor >= end || *cursor < '0' || *cursor > '9')
303      return false;
304
305    int exponent = 0;
306    while (cursor < end && *cursor >= '0' && *cursor <= '9') {
307      exponent = exponent * 10 + (*cursor - '0');
308      ++cursor;
309    }
310
311    exponent *= exp_sign;
312    if (exponent > 0) {
313      while (exponent--)
314        result *= 10.0;
315    } else if (exponent < 0) {
316      while (exponent++)
317        result *= 0.1;
318    }
319  }
320
321  *value = (float)(sign * result);
322  return true;
323}
324
325static int read_catalog_mmap(const char *file_path, float *output_rasc,
326                             float *output_decl, int expected_galaxies,
327                             float arcmin2rad) {
328  int fd = open(file_path, O_RDONLY);
329  if (fd < 0) {
330    printf("   ERROR: Cannot open data file %s\n", file_path);
331    return (EXIT_FAILURE);
332  }
333
334  struct stat file_info;
335  if (fstat(fd, &file_info) != 0 || file_info.st_size <= 0) {
336    printf("   ERROR: Cannot stat data file %s\n", file_path);
337    close(fd);
338    return (EXIT_FAILURE);
339  }
340
341  void *mapped_data =
342      mmap(NULL, (size_t)file_info.st_size, PROT_READ, MAP_PRIVATE, fd, 0);
343  close(fd);
344
345  if (mapped_data == MAP_FAILED) {
346    printf("   ERROR: Cannot memory-map data file %s\n", file_path);
347    return (EXIT_FAILURE);
348  }
349
350  const char *cursor = (const char *)mapped_data;
351  const char *end = cursor + file_info.st_size;
352
353  int number_of_galaxies = 0;
354  if (!parse_int_fast(cursor, end, &number_of_galaxies)) {
355    printf("   ERROR: Cannot read galaxy count in %s\n", file_path);
356    munmap(mapped_data, (size_t)file_info.st_size);
357    return (EXIT_FAILURE);
358  }
359
360  if (number_of_galaxies < expected_galaxies) {
361    printf("   ERROR: File %s has %d galaxies, expected at least %d\n",
362           file_path, number_of_galaxies, expected_galaxies);
363    munmap(mapped_data, (size_t)file_info.st_size);
364    return (EXIT_FAILURE);
365  }
366
367  for (int i = 0; i < expected_galaxies; ++i) {
368    float rasc = 0.0f;
369    float decl = 0.0f;
370
371    if (!parse_float_fast(cursor, end, &rasc) ||
372        !parse_float_fast(cursor, end, &decl)) {
373      printf("   ERROR: Cannot read line %d in data file %s\n", i + 1,
374             file_path);
375      munmap(mapped_data, (size_t)file_info.st_size);
376      return (EXIT_FAILURE);
377    }
378
379    output_rasc[i] = rasc * arcmin2rad;
380    output_decl[i] = decl * arcmin2rad;
381  }
382
383  munmap(mapped_data, (size_t)file_info.st_size);
384  return (EXIT_SUCCESS);
385}
386
387int main(int argc, char **argv) {
388  printf("        Native CUDA Galaxy Correlation - RTX 5080 (sm_120)\n");
389  long int histogramDRsum, histogramDDsum, histogramRRsum;
390  double walltime;
391  double inputReadTimeMs = 0.0;
392  double kernelExecutionTimeMs = 0.0;
393  double outputWriteTimeMs = 0.0;
394  struct timeval _ttime;
395  struct timeval inputStart, inputEnd;
396  get_device();
397
398  gettimeofday(&_ttime, NULL);
399  walltime = (double)_ttime.tv_sec + (double)_ttime.tv_usec / 1000000.;
400
401  // Allocate host memory
402  real_rasc = (float *)calloc(100000L, sizeof(float));
403  real_decl = (float *)calloc(100000L, sizeof(float));
404  rand_rasc = (float *)calloc(100000L, sizeof(float));
405  rand_decl = (float *)calloc(100000L, sizeof(float));
406  real_sin = (float *)calloc(100000L, sizeof(float));
407  real_cos = (float *)calloc(100000L, sizeof(float));
408  rand_sin = (float *)calloc(100000L, sizeof(float));
409  rand_cos = (float *)calloc(100000L, sizeof(float));
410  CPUMemory += 8L * 100000L * sizeof(float);
411
412  // Read input data from files
413  gettimeofday(&inputStart, NULL);
414  if (parseargs_readinput(argc, argv) != 0) {
415    printf("   Program stopped.\n");
416    return (EXIT_FAILURE);
417  }
418  // Precompute per-galaxy sin/cos(decl) once (removes redundant per-pair trig).
419  for (long int k = 0; k < N; ++k) {
420    real_sin[k] = sinf(real_decl[k]);
421    real_cos[k] = cosf(real_decl[k]);
422    rand_sin[k] = sinf(rand_decl[k]);
423    rand_cos[k] = cosf(rand_decl[k]);
424  }
425  gettimeofday(&inputEnd, NULL);
426  inputReadTimeMs = (inputEnd.tv_sec - inputStart.tv_sec) * 1000.0;
427  inputReadTimeMs += (inputEnd.tv_usec - inputStart.tv_usec) / 1000.0;
428
429  printf("   Input data read, now calculating histograms\n");
430
431  FILE *outfile;
432
433  if (argc != 4) {
434    printf("Usage: ./galaxy_cuda data_100k_arcmin.txt flat_100k_arcmin.txt "
435           "omega.out\n");
436    return (EXIT_FAILURE);
437  }
438
439  histogram_DR =
440      (long int *)calloc(totaldegrees * binsperdegree + 1ULL, sizeof(long int));
441  histogram_DD =
442      (long int *)calloc(totaldegrees * binsperdegree + 1ULL, sizeof(long int));
443  histogram_RR =
444      (long int *)calloc(totaldegrees * binsperdegree + 1ULL, sizeof(long int));
445  CPUMemory += 3L * (totaldegrees * binsperdegree + 1L) * sizeof(long int);
446
447  // Allocate GPU device memory (precomputed sin/cos + rasc per catalog)
448  float *d_real_sin, *d_real_cos, *d_real_rasc;
449  float *d_rand_sin, *d_rand_cos, *d_rand_rasc;
450
451  struct timeval t1, t2;
452  double gpuPhaseTimeMs;
453  gettimeofday(&t1, NULL);
454
455  int deviceCount = 0;
456  CUDA_ERR_CHECK(cudaGetDeviceCount(&deviceCount));
457  printf("   \nRunning on %d GPU(s)\n", deviceCount);
458
459  const int tiles_per_dim = (N + TILE - 1) / TILE;
460  dim3 threadsInBlock(TILE, 1, 1);
461  dim3 threadBlocks(tiles_per_dim, tiles_per_dim, 1);
462
463  // Widen before multiplying: dim3 members are unsigned int, so evaluating the
464  // whole product in 32 bits and casting afterwards overflows for large grids.
465  const long int number_of_threads =
466      (long int)threadsInBlock.x * threadsInBlock.y * threadsInBlock.z *
467      threadBlocks.x * threadBlocks.y * threadBlocks.z;
468  const long int threads_per_block = threadsInBlock.x;
469
470  CUDA_ERR_CHECK(cudaSetDevice(0));
471
472  printf(
473      "====================================================================\n");
474
475  printf("    Using GPU 0 (RTX 5080, sm_120)\n");
476  printf("    Tile: %d (%ld threads/block), grid %dx%d blocks\n", TILE,
477         threads_per_block, threadBlocks.x, threadBlocks.y);
478
479  // Allocate device memory
480  size_t inputDataArrayBytes = N * sizeof(float);
481  CUDA_ERR_CHECK(cudaMalloc((void **)&d_real_sin, inputDataArrayBytes));
482  CUDA_ERR_CHECK(cudaMalloc((void **)&d_real_cos, inputDataArrayBytes));
483  CUDA_ERR_CHECK(cudaMalloc((void **)&d_real_rasc, inputDataArrayBytes));
484  CUDA_ERR_CHECK(cudaMalloc((void **)&d_rand_sin, inputDataArrayBytes));
485  CUDA_ERR_CHECK(cudaMalloc((void **)&d_rand_cos, inputDataArrayBytes));
486  CUDA_ERR_CHECK(cudaMalloc((void **)&d_rand_rasc, inputDataArrayBytes));
487  GPUMemory += 6L * inputDataArrayBytes;
488
489  // Allocate histogram arrays
490  unsigned long long int *d_histogram_DR, *d_histogram_DD, *d_histogram_RR;
491  size_t histogramArrayBytes =
492      (totaldegrees * binsperdegree + 1ULL) * sizeof(unsigned long long);
493
494  CUDA_ERR_CHECK(cudaMalloc((void **)&d_histogram_DR, histogramArrayBytes));
495  CUDA_ERR_CHECK(cudaMalloc((void **)&d_histogram_DD, histogramArrayBytes));
496  CUDA_ERR_CHECK(cudaMalloc((void **)&d_histogram_RR, histogramArrayBytes));
497  CUDA_ERR_CHECK(cudaMemset(d_histogram_DR, 0, histogramArrayBytes));
498  CUDA_ERR_CHECK(cudaMemset(d_histogram_DD, 0, histogramArrayBytes));
499  CUDA_ERR_CHECK(cudaMemset(d_histogram_RR, 0, histogramArrayBytes));
500  GPUMemory += 3L * histogramArrayBytes;
501
502  // Copy input data to device
503  CUDA_ERR_CHECK(cudaMemcpy(d_real_sin, real_sin, inputDataArrayBytes,
504                            cudaMemcpyHostToDevice));
505  CUDA_ERR_CHECK(cudaMemcpy(d_real_cos, real_cos, inputDataArrayBytes,
506                            cudaMemcpyHostToDevice));
507  CUDA_ERR_CHECK(cudaMemcpy(d_real_rasc, real_rasc, inputDataArrayBytes,
508                            cudaMemcpyHostToDevice));
509  CUDA_ERR_CHECK(cudaMemcpy(d_rand_sin, rand_sin, inputDataArrayBytes,
510                            cudaMemcpyHostToDevice));
511  CUDA_ERR_CHECK(cudaMemcpy(d_rand_cos, rand_cos, inputDataArrayBytes,
512                            cudaMemcpyHostToDevice));
513  CUDA_ERR_CHECK(cudaMemcpy(d_rand_rasc, rand_rasc, inputDataArrayBytes,
514                            cudaMemcpyHostToDevice));
515
516  printf("    threadBlocks:\t\t{%d, %d, %d} blocks.\n    threadsInBlock:\t\t%d "
517         "threads.\n",
518         threadBlocks.x, threadBlocks.y, threadBlocks.z,
519         threadsInBlock.x * threadsInBlock.y * threadsInBlock.z);
520  printf("    Total number of threads:\t%ld\n", number_of_threads);
521
522  // Launch kernel with padded shared histograms + j-tile coordinate buffers.
523  size_t sharedMemSize =
524      3 * num_bins_padded * sizeof(unsigned int) + 6 * TILE * sizeof(float);
525  printf("    Shared memory per block:\t%zu bytes (padded: %d bins, tile %d)\n",
526         sharedMemSize, num_bins_padded, TILE);
527
528  // CUDA-event kernel timing (finer than gettimeofday).
529  cudaEvent_t kstart, kstop;
530  CUDA_ERR_CHECK(cudaEventCreate(&kstart));
531  CUDA_ERR_CHECK(cudaEventCreate(&kstop));
532  CUDA_ERR_CHECK(cudaEventRecord(kstart));
533
534  fill_histograms<<<threadBlocks, threadsInBlock, sharedMemSize, 0>>>(
535      d_real_rasc, d_real_sin, d_real_cos, d_rand_rasc, d_rand_sin, d_rand_cos,
536      d_histogram_DR, d_histogram_DD, d_histogram_RR);
537
538  CUDA_ERR_CHECK(cudaGetLastError());
539  CUDA_ERR_CHECK(cudaEventRecord(kstop));
540  CUDA_ERR_CHECK(cudaEventSynchronize(kstop));
541  float kernelMs = 0.0f;
542  CUDA_ERR_CHECK(cudaEventElapsedTime(&kernelMs, kstart, kstop));
543  kernelExecutionTimeMs = (double)kernelMs;
544  CUDA_ERR_CHECK(cudaEventDestroy(kstart));
545  CUDA_ERR_CHECK(cudaEventDestroy(kstop));
546
547  // Copy results back to host
548  CUDA_ERR_CHECK(cudaMemcpy(histogram_DR, d_histogram_DR, histogramArrayBytes,
549                            cudaMemcpyDeviceToHost));
550  CUDA_ERR_CHECK(cudaMemcpy(histogram_DD, d_histogram_DD, histogramArrayBytes,
551                            cudaMemcpyDeviceToHost));
552  CUDA_ERR_CHECK(cudaMemcpy(histogram_RR, d_histogram_RR, histogramArrayBytes,
553                            cudaMemcpyDeviceToHost));
554
555  // Free device memory
556  CUDA_ERR_CHECK(cudaFree(d_real_rasc));
557  CUDA_ERR_CHECK(cudaFree(d_real_sin));
558  CUDA_ERR_CHECK(cudaFree(d_real_cos));
559  CUDA_ERR_CHECK(cudaFree(d_rand_rasc));
560  CUDA_ERR_CHECK(cudaFree(d_rand_sin));
561  CUDA_ERR_CHECK(cudaFree(d_rand_cos));
562  CUDA_ERR_CHECK(cudaFree(d_histogram_DR));
563  CUDA_ERR_CHECK(cudaFree(d_histogram_DD));
564  CUDA_ERR_CHECK(cudaFree(d_histogram_RR));
565
566  gettimeofday(&t2, NULL);
567  gpuPhaseTimeMs = (t2.tv_sec - t1.tv_sec) * 1000.0;
568  gpuPhaseTimeMs += (t2.tv_usec - t1.tv_usec) / 1000.0;
569  printf("Kernel execution time: %f ms.\n", kernelExecutionTimeMs);
570
571  // Free host memory
572  free(real_rasc);
573  free(real_decl);
574  free(rand_rasc);
575  free(rand_decl);
576  free(real_sin);
577  free(real_cos);
578  free(rand_sin);
579  free(rand_cos);
580
581  // Verify histogram sums
582  histogramDRsum = 0L;
583  for (int i = 0; i < binsperdegree * totaldegrees; ++i)
584    histogramDRsum += histogram_DR[i];
585  printf("results:\n");
586  printf("   DR histogram sum = %ld\n", histogramDRsum);
587
588  if (histogramDRsum != 10000000000L) {
589    printf("   Incorrect histogram sum, exiting.. histogramDRsum: %ld\t\n   "
590           "percentage of target: %15f\n",
591           histogramDRsum, ((float)histogramDRsum / (float)(N * N)));
592    return (EXIT_FAILURE);
593  }
594
595  histogramDDsum = 0L;
596  for (int i = 0; i < binsperdegree * totaldegrees; ++i)
597    histogramDDsum += histogram_DD[i];
598  printf("   DD histogram sum = %ld\n", histogramDDsum);
599  if (histogramDDsum != 10000000000L) {
600    printf("   Incorrect histogram sum, exiting.. histogramDDsum: %ld\n",
601           histogramDDsum);
602    return (EXIT_FAILURE);
603  }
604
605  histogramRRsum = 0L;
606  for (int i = 0; i < binsperdegree * totaldegrees; ++i)
607    histogramRRsum += histogram_RR[i];
608  printf("   RR histogram sum = %ld\n", histogramRRsum);
609  if (histogramRRsum != 10000000000L) {
610    printf("   Incorrect histogram sum, exiting..histogramRRsum: %ld\n",
611           histogramRRsum);
612    return (EXIT_FAILURE);
613  }
614
615  struct timeval outputStart, outputEnd;
616  gettimeofday(&outputStart, NULL);
617
618  // Write results to output file
619  outfile = fopen(argv[3], "w");
620  if (outfile == NULL) {
621    printf("Cannot open output file %s\n", argv[3]);
622    return (-1);
623  }
624
625  fprintf(outfile, "bin start\t\tomega\t        hist_DD\t        hist_DR\t     "
626                   "   hist_RR\n");
627
628  for (int i = 0; i < binsperdegree * totaldegrees; ++i) {
629    if (histogram_RR[i] > 0) {
630      float omega = (histogram_DD[i] - 2 * histogram_DR[i] + histogram_RR[i]) /
631                    ((float)(histogram_RR[i]));
632      fprintf(outfile, "%6.3f\t%15f\t%15ld\t%15ld\t%15ld\n",
633              ((float)i) / binsperdegree, omega, histogram_DD[i],
634              histogram_DR[i], histogram_RR[i]);
635      if (i < 5)
636        printf("   %6.4f", omega);
637    } else {
638      if (i < 5)
639        printf("         ");
640    }
641  }
642  printf("\n");
643
644  fclose(outfile);
645
646  gettimeofday(&outputEnd, NULL);
647  outputWriteTimeMs = (outputEnd.tv_sec - outputStart.tv_sec) * 1000.0;
648  outputWriteTimeMs += (outputEnd.tv_usec - outputStart.tv_usec) / 1000.0;
649
650  // Free host memory
651  free(histogram_DR);
652  free(histogram_DD);
653  free(histogram_RR);
654
655  printf("   Results written to file %s\n", argv[3]);
656  printf("   CPU memory allocated  = %.2lf MB\n", CPUMemory / 1000000.0);
657  printf("   GPU memory allocated  = %.2lf MB\n", GPUMemory / 1000000.0);
658  printf("   Timing breakdown      = input %.3f ms | kernel %.3f ms | output "
659         "%.3f ms\n",
660         inputReadTimeMs, kernelExecutionTimeMs, outputWriteTimeMs);
661  printf("   GPU phase time        = %.3f ms\n", gpuPhaseTimeMs);
662
663  gettimeofday(&_ttime, NULL);
664  walltime =
665      (double)(_ttime.tv_sec) + (double)(_ttime.tv_usec / 1000000.0) - walltime;
666
667  printf("   Total wall clock time = %.3lf s\n", walltime);
668
669  return (EXIT_SUCCESS);
670}
671
672static int get_device() {
673  int deviceCount;
674  CUDA_ERR_CHECK(cudaGetDeviceCount(&deviceCount));
675
676  printf("   Found %d CUDA devices\n", deviceCount);
677  if (deviceCount < 0 || deviceCount > 128)
678    return (EXIT_FAILURE);
679
680  int device;
681  for (device = 0; device < deviceCount; ++device) {
682    cudaDeviceProp deviceProp;
683    CUDA_ERR_CHECK(cudaGetDeviceProperties(&deviceProp, device));
684    printf("      Device %s | device %d\n", deviceProp.name, device);
685    printf("         compute capability           =         %d.%d\n",
686           deviceProp.major, deviceProp.minor);
687    printf("         totalGlobalMemory            =        %.2lf GB\n",
688           deviceProp.totalGlobalMem / 1000000000.0);
689    printf("         l2CacheSize                  =    %8d B\n",
690           deviceProp.l2CacheSize);
691    printf("         regsPerBlock                 =    %8d\n",
692           deviceProp.regsPerBlock);
693    printf("         multiProcessorCount          =    %8d\n",
694           deviceProp.multiProcessorCount);
695    printf("         maxThreadsPerMultiprocessor  =    %8d\n",
696           deviceProp.maxThreadsPerMultiProcessor);
697    printf("         sharedMemPerBlock            =    %8d B\n",
698           (int)deviceProp.sharedMemPerBlock);
699    printf("         warpSize                     =    %8d\n",
700           deviceProp.warpSize);
701    int clockRateKHz = 0;
702    cudaDeviceGetAttribute(&clockRateKHz, cudaDevAttrClockRate, device);
703    printf("         clockRate                    =    %8.2lf MHz\n",
704           clockRateKHz / 1000.0);
705    printf("         maxThreadsPerBlock           =    %8d\n",
706           deviceProp.maxThreadsPerBlock);
707  }
708
709  CUDA_ERR_CHECK(cudaSetDevice(0));
710  CUDA_ERR_CHECK(cudaGetDevice(&device));
711  if (device != 0)
712    printf("   Unable to set device 0, using %d instead", device);
713  else
714    printf("   Using CUDA device %d\n\n", device);
715
716  return (EXIT_SUCCESS);
717}
718
719static int parseargs_readinput(int argc, char *argv[]) {
720  FILE *out_file;
721  constexpr int expected_galaxies = 100000;
722  const float arcmin2rad = 1.0f / 60.0f / 180.0f * PI;
723
724  if (argc != 4) {
725    printf("   Usage: galaxy real_data random_data output_file\n");
726    return (EXIT_FAILURE);
727  }
728
729  printf("   Running galaxy_cuda %s %s %s\n", argv[1], argv[2], argv[3]);
730
731  if (read_catalog_mmap(argv[1], real_rasc, real_decl, expected_galaxies,
732                        arcmin2rad) != 0) {
733    return (EXIT_FAILURE);
734  }
735  printf("   Successfully read %d lines from %s\n", expected_galaxies, argv[1]);
736
737  if (read_catalog_mmap(argv[2], rand_rasc, rand_decl, expected_galaxies,
738                        arcmin2rad) != 0) {
739    return (EXIT_FAILURE);
740  }
741  printf("   Successfully read %d lines from %s\n", expected_galaxies, argv[2]);
742
743  out_file = fopen(argv[3], "w");
744  if (out_file == NULL) {
745    printf("   ERROR: Cannot open output file %s\n", argv[3]);
746    return (EXIT_FAILURE);
747  }
748  fclose(out_file);
749
750  return (EXIT_SUCCESS);
751}

Source