Distributed & GPU Sample Sort, and API Selection Strategy

Having mastered shared-memory sample sort, we now extend this powerful algorithm to distributed-memory clusters using MPI and massively parallel GPUs using CUDA.

To conclude the course, we synthesize all concepts learned across Chapters 1 through 5, formulating a rigorous decision framework to answer the central question of parallel computing: Which parallel programming API is best suited for a given problem?

---

7.19 Distributed Sample Sort with MPI

In distributed memory, each of the $p$ processes initially owns a sublist of $n / p$ keys. The final sorted list is either gathered onto Process 0 or left distributed in global order across the processes.

Version 1: Standard MPI All-to-All Redistribution

The first implementation utilizes built-in MPI collective communication primitives:

/* Program 7.5: First MPI Sample Sort Outline */
loc_s = s / p;
Gen_sample(my_rank, loc_list, loc_n, loc_samp, loc_s);
Gather_to_0(loc_samp, global_sample);

if (my_rank == 0) {
    Find_splitters(global_sample, splitters);
}
Broadcast_from_0(splitters);

/* Count elements destined for each process */
Count_elts_going_to_procs(loc_list, splitters, my_to_counts);

/* Exchange send counts so each process knows incoming counts */
MPI_Alltoall(my_to_counts, 1, MPI_INT, my_fr_counts, 1, MPI_INT, comm);

/* Calculate send and receive displacement offsets */
Excl_prefix_sums(my_to_counts, my_to_offsets, p);
Excl_prefix_sums(my_fr_counts, my_fr_offsets, p);

/* Global data redistribution */
MPI_Alltoallv(loc_list, my_to_counts, my_to_offsets, MPI_INT,
              tlist, my_fr_counts, my_fr_offsets, MPI_INT, comm);

/* Local sort of received bucket */
Sort(tlist, my_new_count);

/* Gather variable-length buckets onto Process 0 */
MPI_Gather(&my_new_count, 1, MPI_INT, bkt_counts, 1, MPI_INT, 0, comm);
if (my_rank == 0) Excl_prefix_sums(bkt_counts, bkt_offsets, p);
MPI_Gatherv(tlist, my_new_count, MPI_INT, list, bkt_counts, bkt_offsets, MPI_INT, 0, comm);

Table 7.14: MPI Version 1 Run-times ($n = 2^{22} = 4,194,304$ integers, $s = 16,384$)

Processes ($p$)	MPI Version 1 Run-time	Speedup vs. Serial qsort ($0.739\,\text{s}$)	Parallel Efficiency ($E$)
1	$0.994\,\text{s}$	$0.74\times$	$1.00$
2	$0.477\,\text{s}$	$1.55\times$	$0.96$
4	$0.253\,\text{s}$	$2.92\times$	$0.91$
8	$0.132\,\text{s}$	$5.60\times$	$0.88$
16	$0.0704\,\text{s}$	$10.50\times$	$0.88$
32	$0.0598\,\text{s}$	$12.36\times$	$0.52$

---

Version 2: Custom Butterfly Communication & Nonblocking Memory Optimization

To eliminate the scaling bottleneck at $p = 32$ processes, Version 2 introduces two key innovations:

1. Parallel Odd-Even Sample Sort & Splitter Exchange: Processes sort subsamples in parallel and compute splitters by exchanging boundary keys with neighboring ranks (prev_max and my_min), avoiding serial bottlenecks on rank 0.
2. Logarithmic Butterfly Exchange: Replaces MPI_Alltoallv with a custom $\log_2(p)$-stage butterfly network, where partner ranks exchange upper/lower halves of their keys using bitwise XOR masks (partner = my_rank ^ bitmask).

flowchart TD
  subgraph ButterflyAlltoall["Figure 7.11: Butterfly Exchange across 4 Processes"]
      direction TB
      subgraph Stage1["Stage 1: Middle Splitter e2 = 6 (bitmask = 2)"]
          P0["Proc 0: [3, 5, 6]"] <-->|Exchange| P2["Proc 2: [7, 8, 12]"]
          P1["Proc 1: [2, 9, 10]"] <-->|Exchange| P3["Proc 3: [1, 4, 11]"]
      end
      subgraph Stage2["Stage 2: Splitters e1 = 4, e3 = 9 (bitmask = 1)"]
          P0b["Proc 0: [3, 5]"] <-->|Exchange| P1b["Proc 1: [1, 2, 4]"]
          P2b["Proc 2: [6, 7, 8, 12]"] <-->|Exchange| P3b["Proc 3: [9, 10, 11]"]
      end
      Stage1 --> Stage2
  end

Memory Optimization via MPI_Probe and MPI_Isend

Rather than pre-allocating an oversized $n$-element buffer on every rank, processes query incoming message sizes dynamically with MPI_Probe and MPI_Get_count, resizing buffers with realloc() only when necessary:

/* Program 7.7: Memory-Efficient Nonblocking Send/Receive */
void Send_recv(int snd_buf[], int count, int** rcv_buf_p,
               int* rcv_buf_sz_p, int partner, MPI_Comm comm) {
    MPI_Request req;
    MPI_Status status;
    int rcv_count;

    /* Start nonblocking transmission */
    MPI_Isend(snd_buf, count, MPI_INT, partner, 0, comm, &req);

    /* Inspect incoming message metadata without receiving */
    MPI_Probe(partner, 0, comm, &status);
    MPI_Get_count(&status, MPI_INT, &rcv_count);

    if (rcv_count > *rcv_buf_sz_p) {
        *rcv_buf_p = realloc(*rcv_buf_p, rcv_count * sizeof(int));
        *rcv_buf_sz_p = rcv_count;
    }

    /* Receive directly into sized buffer */
    MPI_Recv(*rcv_buf_p, rcv_count, MPI_INT, partner, 0, comm, MPI_STATUS_IGNORE);
    MPI_Wait(&req, MPI_STATUS_IGNORE);
}

Table 7.15: Run-time Comparison: Alltoallv vs. Butterfly ($n = 2^{22}$)

Processes	Version 1 (MPI_Alltoallv)	Version 2 (Butterfly)	Performance Gain
1	$0.994\,\text{s}$	$0.839\,\text{s}$	$+15.6\%$
2	$0.477\,\text{s}$	$0.403\,\text{s}$	$+15.5\%$
4	$0.253\,\text{s}$	$0.210\,\text{s}$	$+17.0\%$
8	$0.132\,\text{s}$	$0.110\,\text{s}$	$+16.7\%$
16	$0.0704\,\text{s}$	$0.0636\,\text{s}$	$+9.7\%$
32	$0.0598\,\text{s}$	$0.0467\,\text{s}$	$+21.9\%$

---

7.20 GPU Sample Sort with CUDA

Sorting on GPUs requires mapping bucket operations to thousands of SIMD/SIMT execution units.

Version 1: Single-Block Shared Memory Solver

In the basic CUDA implementation, sorting is confined to a single thread block (up to $1024$ threads):

1. The CPU host selects sample keys and computes splitters.
2. In kernel Dev_ssort, threads evaluate key counts and perform a parallel exclusive prefix sum in __shared__ memory.
3. Keys are routed into shared memory buckets, and each thread sorts its bucket with a single-threaded heap sort.

Table 7.16: Hillis-Steele Parallel Exclusive Prefix Sum Trace in Shared Memory

Step / Shift	Thread 0	Thread 1	Thread 2	Thread 3	Thread 4	Thread 5
Input ($x$)	$3$	$6$	$7$	$5$	$3$	$5$
Init (Right Shift 1)	$0$	$3$	$6$	$7$	$5$	$3$
Shift = 1	$0$	$3$	$9$	$13$	$12$	$8$
Shift = 2	$0$	$3$	$9$	$16$	$21$	$21$
Shift = 4	$0$	$3$	$9$	$16$	$21$	$24$

Warp Divergence in Serial Sorts

When 32 threads in the same warp run single-threaded heap sort on different buckets, divergent data branches force threads to execute disjoint code paths sequentially. Benchmarks show thread divergence increases execution time by over $50\%$ compared to independent blocks!

---

Version 2: Scalable Multi-Block Grid Sample Sort

To sort millions of elements without warp divergence, Version 2 maps one thread block per bucket:

1. Block Bitonic Pre-sort: Each block copies its sublist into shared memory and sorts it using bitonic sort.
2. Global Bitonic Sample Sort: The sample is sorted across the grid using multi-block bitonic sort.
3. Parallel Splitter Extraction: Threads extract splitters concurrently.
4. Binary Search & Atomic Counting: Threads search splitters in parallel, using atomicAdd into mat_counts[i][j].
5. Cache-Friendly Prefix Sums: Row-wise prefix sums establish bucket starting indices.
6. Bucket Padding & Bitonic Final Sort: Each bucket is padded with $\infty$ to a power of 2 and sorted in parallel using bitonic sort.

Table 7.20: Performance of Multi-Block CUDA Sample Sort ($s = n / 8$)

Blocks (blk_ct)	Threads/Block	List Size ($n$)	CUDA Run-time	Serial CPU qsort	GPU Speedup
256	256	$131,072$	$0.00158\,\text{s}$	$0.0290\,\text{s}$	$18.3\times$
256	512	$262,144$	$0.00214\,\text{s}$	$0.0549\,\text{s}$	$25.6\times$
256	1024	$524,288$	$0.00328\,\text{s}$	$0.122\,\text{s}$	$37.2\times$
512	1024	$1,048,576$	$0.00625\,\text{s}$	$0.246\,\text{s}$	$39.4\times$
1024	1024	$2,097,152$	$0.0142\,\text{s}$	$0.520\,\text{s}$	$36.6\times$

At 2 million elements, multi-block CUDA sample sort finishes in $14.2$ milliseconds, outperforming optimized sequential C qsort by over $36\times$!

---

7.21 Which API Should You Choose?

Selecting the optimal parallel programming API is one of the most critical architecture decisions in software engineering. Consider the following decision framework:

flowchart TD
  Start["Start: Parallel Architecture Decision"] --> Q1{"Is problem data-parallel<br/>with regular SIMD control flow?"}
  Q1 -- Yes --> GPU["NVIDIA CUDA<br/>Best for compute-intensive,<br/>regular data parallelism"]
  Q1 -- No --> Q2{"Does dataset exceed<br/>single-node physical RAM?"}
  Q2 -- Yes --> MPI["MPI (Distributed Memory)<br/>Massive scalability, clusters,<br/>high aggregate cache"]
  Q2 -- No --> Q3{"Is there an existing<br/>complex serial C/C++ codebase?"}
  Q3 -- Yes --> OMP["OpenMP (Shared Memory)<br/>Incremental directive insertion,<br/>minimal refactoring"]
  Q3 -- No --> Q4{"Does code require fine-grained<br/>thread signaling or custom locks?"}
  Q4 -- Yes --> PTH["POSIX Threads (Pthreads)<br/>Explicit low-level thread lifecycle<br/>and synchronization"]
  Q4 -- No --> OMP

Comprehensive API Comparison Matrix

Feature / Criterion	OpenMP	Pthreads	MPI	CUDA
Hardware Target	Shared-memory multicore	Shared-memory multicore	Distributed clusters / supercomputers	Massively parallel GPUs
Address Space	Unified Shared Memory	Unified Shared Memory	Private Address Spaces	Separate Host & Device Memory
Programming Paradigm	Compiler Pragmas & Directives	Explicit Library Functions	Explicit Message Passing	Kernel Extensions & Host API
Development Effort	Low (incremental parallelization)	Medium (explicit thread management)	High (complete algorithmic redesign)	Medium–High (memory transfers & kernel tuning)
Scalability Limit	Single-node cores ($8–128$)	Single-node cores ($8–128$)	Thousands of nodes ($10^3–10^6$)	Massive on-chip parallelism ($10^4–10^6$ threads)
Best Used For	Loop parallelization, tasks, quick porting	Complex thread lifecycles, background workers	Large-scale scientific HPC, big data simulations	High-throughput data parallelism, linear algebra

---

7.22 Chapter Summary & Advanced MPI Reference

Key Concepts Mastered:

1. Foster’s Methodology: Partitioning, Communication, Agglomeration, and Mapping provide a rigorous framework for building parallel algorithms from scratch.
2. Algorithmic Reductions: Exploiting physical symmetries (like Newton’s third law) halves arithmetic operations but requires careful synchronization (two-phase accumulation in OpenMP, the ring pass in MPI).
3. Programmer-Managed Caching: When GPU global memory bandwidth saturates, shared-memory tiling reduces DRAM traffic by $1000\times$.
4. Adaptive Sorting: Sample sort adapts to arbitrary key distributions through statistical splitters, scaling seamlessly across threads, processes, and GPU blocks.

---

Advanced MPI Functions Reference

/* 1. Variable-length gather onto root */
int MPI_Gatherv(
    const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    void* recvbuf, const int recvcounts[], const int displs[],
    MPI_Datatype recvtype, int root, MPI_Comm comm
);

/* 2. Equal-length all-to-all scatter-gather */
int MPI_Alltoall(
    const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    void* recvbuf, int recvcount, MPI_Datatype recvtype,
    MPI_Comm comm
);

/* 3. Variable-length all-to-all scatter-gather */
int MPI_Alltoallv(
    const void* sendbuf, const int sendcounts[], const int sdispls[], MPI_Datatype sendtype,
    void* recvbuf, const int recvcounts[], const int rdispls[], MPI_Datatype recvtype,
    MPI_Comm comm
);

/* 4. Nonblocking asynchronous send */
int MPI_Isend(
    const void* buf, int count, MPI_Datatype datatype,
    int dest, int tag, MPI_Comm comm, MPI_Request* request
);

/* 5. Complete nonblocking communication */
int MPI_Wait(MPI_Request* request, MPI_Status* status);

/* 6. Non-destructive incoming message inspection */
int MPI_Probe(int source, int tag, MPI_Comm comm, MPI_Status* status);

/* 7. Emergency communicator termination */
int MPI_Abort(MPI_Comm comm, int errorcode);