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

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. Right ascension and declination are Earth’s longitude and latitude projected onto the sky: $\alpha$ sweeps around the celestial equator and $\delta$ climbs from it toward the pole, so together they pin one galaxy to one point on the celestial sphere, with Earth at the centre.

A celestial sphere around Earth, with a galaxy located by a right-ascension angle swept along the equator and a declination angle climbing from the equator to the galaxy — A galaxy's position is two angles measured from Earth: right ascension α along the celestial equator, and declination δ climbing from it.

For every pair of galaxies the angular separation is the angle between their two radius vectors, as seen from Earth at the centre of that sphere:

Two galaxies on the celestial sphere, each connected to Earth by a radius vector, with the angle theta-one-two marked between the two vectors at Earth — The angular separation θ₁₂ between two galaxies is the angle their radius vectors make at Earth, which is exactly arccos of the vectors' dot product.

Expressed in right ascension and declination rather than in the underlying Cartesian vectors, that dot product becomes:

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

This is the one line of trigonometry the whole exercise optimizes: it runs once per pair, $N^2$ times per histogram family, and everything from L0 to L5 below is about computing it, and sorting its result into a bin, as fast as the hardware allows.

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