HIP Deep Dive: 10B Galaxy-Pair Angles on AMD (~0.15s)

This is the full, spoiler-allowed deep dive of my HIP solution for the galaxy-correlation assignment.

No abstractions, no hand-waving: this is the exact strategy that took me from an early ~1.60s implementation to ~0.154s wall-clock on later optimized runs (best observed ~0.152s in the latest phase).

If you want to see the data, I also built a small 3D viewer for it: Galaxy Visualization

Problem statement and target metric #

Input:

• 100,000 real galaxies
• 100,000 synthetic (random) galaxies

Each row contains right ascension $\alpha$ and declination $\delta$ in arcminutes, converted to radians during ingest.

Core workload:

• Build histograms for DR, DD, and RR
• Total pair workload is $N^2 = 100000^2 = 10^{10}$ comparisons per histogram family
• Angular bins: 360 * 4 = 1440 bins

The spherical-angle computation can be written as:

// Pseudo-C: angular separation for one pair
float dot = sinf(delta1) * sinf(delta2)
          + cosf(delta1) * cosf(delta2) * cosf(alpha1 - alpha2);
dot = fminf(fmaxf(dot, -1.0f), 1.0f);
float theta12 = acosf(dot);

The final omega value per bin is:

// Pseudo-C: Landy-Szalay estimator at one histogram bin
float omega = (DD[bin] - 2.0f * DR[bin] + RR[bin]) / RR[bin];

Acronym guide #

• HIP: Heterogeneous-computing Interface for Portability (AMD’s CUDA-like GPU programming model)
• LDS: Local Data Share (AMD on-chip shared memory, similar in role to CUDA shared memory)
• RDNA 2: AMD GPU architecture generation used for this tuning pass
• DR/DD/RR: Data-Random, Data-Data, Random-Random pair histograms
• RA/DEC: Right Ascension / Declination sky coordinates
• mmap: Memory mapping file I/O directly into process address space

End-to-end architecture #

The implementation is intentionally simple at a high level:

1. Read both catalogs with mmap + custom ASCII parsers.
2. Copy RA/DEC arrays to GPU.
3. Launch one HIP kernel over a 2D $(i, j)$ space.
4. Accumulate DR/DD/RR in shared-memory histograms.
5. Flush per-block histograms into global 64-bit histograms.
6. Copy histograms back, validate sums, compute omega, write output.

Why this maps well to GPU #

Pair generation is embarrassingly parallel. Histogram updates are not.

That means the real challenge is minimizing contention and global atomic pressure while preserving exact counts.

I map one thread to one $(i, j)$ pair:

const int i = blockIdx.x * blockDim.x + threadIdx.x;
const int j = blockIdx.y * blockDim.y + threadIdx.y;

And use a 2D launch grid:

dim3 threadsInBlock(BLOCK_SIZE_X, BLOCK_SIZE_Y);
dim3 threadBlocks(
    (N + threadsInBlock.x - 1) / threadsInBlock.x,
    (N + threadsInBlock.y - 1) / threadsInBlock.y);

Kernel deep dive #

1) Shared-memory histogram staging #

Each block allocates three shared histograms:

• s_hist_DR
• s_hist_DD
• s_hist_RR

Instead of 1440 bins, I allocate 1536 padded bins per histogram to reduce LDS bank conflict pressure on RDNA-class hardware.

const int num_bins = 1440;
const int num_bins_padded = 1536;
extern __shared__ unsigned int s_hist[];

unsigned int *s_hist_DR = s_hist;
unsigned int *s_hist_DD = s_hist_DR + num_bins_padded;
unsigned int *s_hist_RR = s_hist_DD + num_bins_padded;

Shared memory footprint per block is:

// Pseudo-C: shared-memory budget per block
int histograms = 3;                 // DR, DD, RR
int padded_bins = 1536;
int bytes_per_counter = 4;          // uint32
int shared_bytes = histograms * padded_bins * bytes_per_counter; // 18432

That is a good fit for occupancy while still giving enough local accumulation capacity.

2) Numerical safety in binning #

Floating-point dot products can drift slightly outside $[-1, 1]$, so I clamp before acosf:

expr = fminf(fmaxf(expr, -1.0f), 1.0f);
int histogram_index = int(acosf(expr) * 57.29577951308232f * binsperdegree);
histogram_index = min(max(histogram_index, 0), num_bins - 1);

This removes NaNs at the domain boundary and avoids out-of-range bin writes.

3) DR, DD, RR computation path #

Per valid thread:

• Always compute DR for $(i, j)$.
• Compute DD and RR only for upper triangle (j >= i).
• Use a symmetric increment:
  • diagonal ($i = j$): +1
  • off-diagonal ($i \neq j$): +2

if (j >= i)
{
    const unsigned int symmetric_increment = 1U + static_cast<unsigned int>(j != i);
    // DD update with symmetric_increment
    // RR update with symmetric_increment
}

This cuts DD/RR math work roughly in half while preserving exact totals.

4) Why symmetry is exact #

For DD/RR, the full matrix has $N^2$ entries. Upper triangle with mirrored weighting gives:

// Pseudo-C: exact symmetry-preserving accumulation
for (int i = 0; i < N; ++i) {
    for (int j = i; j < N; ++j) {
        int w = (j == i) ? 1 : 2;
        DD[bin_dd(i, j)] += w;
        RR[bin_rr(i, j)] += w;
    }
}

So we keep mathematically exact counts and avoid duplicate trigonometric work.

5) Atomic strategy #

• Shared memory atomics: 32-bit atomicAdd to LDS arrays.
• Global memory atomics: 64-bit atomicAdd when flushing block totals.

if (s_hist_DR[b] > 0)
    atomicAdd(&d_histogram_DR[b], (unsigned long long int)s_hist_DR[b]);

Rationale:

• LDS updates are much cheaper than global updates.
• 32-bit in LDS is enough per block (bounded per-kernel contribution).
• 64-bit global counters safely hold totals up to $10^{10}$.

Input pipeline: why mmap mattered #

I replaced fscanf loops with:

• open + fstat
• mmap whole file
• custom fast ASCII parsers (parse_int_fast, parse_float_fast)

This improved startup cost and made timing more stable.

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

I still keep strict validation:

• verify line count availability
• fail on malformed rows
• convert arcminutes to radians at ingest

Build and hardware tuning choices #

From the build script and runtime tuning:

• -O3
• -ffast-math
• -munsafe-fp-atomics
• --offload-arch=gfx1030
• -mwavefrontsize64

I tested multiple block shapes (16x32, 32x16, 16x16, and 32x32). The right answer depends on occupancy, LDS pressure, and register usage on your exact GPU/driver stack.

Correctness guardrails #

I fail hard if histogram sums are wrong:

if (histogramDRsum != 10000000000L) exit(EXIT_FAILURE);
if (histogramDDsum != 10000000000L) exit(EXIT_FAILURE);
if (histogramRRsum != 10000000000L) exit(EXIT_FAILURE);

This catches race/logic regressions immediately.

On Linux x86-64, long is 64-bit, so 10^10 fits on host-side sum checks.

Timing protocol and measured progression #

My benchmark protocol is 5 back-to-back runs:

• Run 1: warm-up/cache effects
• Runs 2-5: steady-state behavior

I also print explicit phase timing:

• input
• kernel
• output
• end-to-end GPU phase

Runtime progression for the galaxy pair angle computation across different implementation stages.

Best observed run in the latest phase: 0.152s.

Key utility pattern: fail-fast HIP errors #

#define HIP_ERR_CHECK(ans)                    \
    {                                         \
        gpuAssert((ans), __FILE__, __LINE__); \
    }

static inline void
gpuAssert(hipError_t code, const char *file, int line, bool abort = true)
{
    if (code != hipSuccess)
    {
        fprintf(stderr, "   GPUassert: %s %s %d\n", hipGetErrorString(code), file, line);
        if (abort)
            exit(code);
    }
}

Practical lessons from this implementation #

1. Histogram contention dominates long before floating-point math does.
2. Shared-memory staging + padded bins is worth real time on AMD GPUs.
3. Symmetry in DD/RR is one of the highest-impact algorithmic wins.
4. Input path (mmap + parser) matters once kernels are fast.
5. Hard correctness checks save days of benchmarking bad results.