Parallel algorithms

A reduction combines a collection of elements into a single result using an operation such as addition, minimum, maximum, logical conjunction or a custom associative combination. The operation must have well-understood algebraic properties. Associativity permits regrouping, which is the foundation of parallel reduction trees. An identity value permits clean handling of empty or padded segments. Commutativity may allow arbitrary reordering, but not every associative operation is commutative: matrix multiplication and string concatenation provide familiar counterexamples. A correct implementation should preserve element order when the operator requires it, even if it still uses a balanced tree of contiguous segments.

To reduce sixteen integer values, eight workers can each read two and compute a local sum. Four sums are created in the next combination level, then two, then one. The tree has logarithmic depth in the number of input elements if enough execution resources exist. A sequential sum needs n-1 addition operations on a dependency chain of length n-1. A balanced reduction still performs n-1 additions overall but organizes them into roughly log2(n) dependent levels. This illustrates the distinction between work and span. A parallel algorithm is attractive when it preserves reasonable total work while reducing span substantially.

The simplest tree reduction assumes a power-of-two length, but real datasets rarely respect such a restriction. Padding with the operation's identity can make the logical length a power of two if storage and synchronization are handled safely. Alternatively, an implementation can use a ragged tree with correct bounds checks. Testing n equal to 0, 1, 3, 17 and values around a block boundary exposes many indexing defects that n=1,024 conceals. An empty reduction needs an explicit API contract: some operations return their identity, while others reject an empty input because no appropriate identity is defined.

Floating-point reduction introduces a more subtle concern. Floating-point addition rounds after each operation, and changing grouping can change results. For example, a large positive number, a small value and the matching large negative number may produce different rounded sums depending on order. A tree can improve numerical error properties compared with certain sequential orders, but it cannot guarantee bitwise identity with a conventional left-to-right sum. Numerically sensitive applications may use pairwise summation, compensated summation, superaccumulators or reproducible reduction libraries. A tolerance should depend on the problem's conditioning and desired accuracy rather than an arbitrary fixed count of decimal digits.

Contention is the main reason to avoid a shared global accumulator updated by every thread. Even if each update is atomic and correct, the single memory location may become a coherence hotspot. Private local reductions reduce synchronization frequency; a final hierarchy combines partial values. On GPUs, a block can reduce locally in shared memory or with lane-level operations, then write one partial result. A separate kernel combines block partials, or suitable atomic operations may be used when nondeterministic ordering is acceptable. The choice depends on payload size, hardware support and the desired reproducibility contract.

In distributed memory, reduction is a communication problem as well as an arithmetic problem. A local sum on each rank is cheap when data remain local, but combining across ranks requires a reduction tree or another collective algorithm. For small payloads, startup latency and tree depth dominate. For large arrays, pipelined or ring-based communication can use bandwidth more effectively. Data types and operations must match across ranks. If one rank contributes a mean calculated over a different number of samples, simply reducing the means may give the wrong global mean. Reduce sums and counts separately, then divide by the total count.

Reduction performance is best evaluated with work, span, memory traffic and communication depth together. A clever tree that minimizes arithmetic depth can still be slower than a simple loop on small arrays because of launch overhead. A tree that uses more temporary storage may displace useful data from cache. For large arrays, memory bandwidth may dominate long before arithmetic reduction time becomes significant. Measure a range of n and p, and compare to a well-tuned serial reference under equivalent numerical acceptance criteria.

A reduction is not merely a library primitive. It is an algebraic specification combined with a hardware mapping. State the operator, identity, data type, overflow policy and ordering expectations; then choose the tree and communication strategy. This careful framing also prepares the reader for parallel scans, which compute many related prefix reductions rather than only one final result.

A scan, or prefix operation, produces intermediate accumulations rather than only the final reduction. For an array [2,3,5,7], an inclusive prefix sum returns [2,5,10,17]; an exclusive prefix sum returns [0,2,5,10] when zero is the identity. Both forms are valuable, but confusing them causes off-by-one errors. Scans support stream compaction, radix sorting, histogram offsets, memory allocation, sparse matrix construction and parallel recurrence transformations. Their central property is that every output is a reduction of a prefix of the input sequence.

A sequential scan traverses the array once, maintaining a running accumulator, so it does O(n) work but has O(n) dependency depth. A naive parallel design might independently sum each prefix; that creates O(n squared) total work, wasting hardware despite abundant parallelism. Tree-based scans achieve much better work-depth tradeoffs. A common conceptual algorithm performs an upsweep to build partial sums, then a downsweep to distribute prefix values. For a power-of-two array, the two tree traversals have logarithmic depth and linear work. Each phase has dependencies that must be respected through barriers or kernel boundaries.

The upsweep-downsweep description is easiest to understand with eight items. Adjacent pairs first compute partial sums. Pairs of pairs combine, ultimately producing the total. The downsweep starts with the identity at the root and propagates prefix information: the left child inherits the parent's prefix, while the right child receives that prefix combined with the left child's aggregate. Leaves then contain exclusive prefix results. To obtain an inclusive result, combine each original input with its exclusive prefix. The logic preserves order for associative but noncommutative operators when traversal is arranged appropriately.

On GPUs, a scan is usually hierarchical because a single block cannot globally synchronize arbitrary other blocks in an ordinary launch. Each block computes local prefixes and a block total. A second scan computes offsets for the block totals. A final pass adds each block's offset to its local prefixes. This design incurs extra passes and temporary storage but scales to arrays much larger than one block. Advanced algorithms can fuse stages with carefully defined synchronization and progress guarantees, but these are harder to verify. A straightforward three-stage implementation is often the best teaching baseline.

Scan correctness requires an identity and a valid associative operator. For string concatenation, order matters: 'ab' followed by 'cd' must not be reordered as 'cdab'. For floating-point addition, regrouping changes rounding as in reductions, so mathematical associativity is only approximate in machine arithmetic. For integer counts, overflow behavior must be specified. If prefix counts exceed the maximum representable integer, the resulting memory offsets can become dangerous even when individual input values are valid. Use data types large enough for the total and check dimensions before allocating output arrays.

Stream compaction provides a worked design context. Suppose every input record has a predicate indicating whether it should survive filtering. Evaluate the predicate to obtain a zero-one flag for each record. An exclusive scan of flags yields the destination position of each surviving record. A thread with flag one writes its record to output at that position. The final output length is the total sum of flags. The procedure preserves input order and can run in parallel without output collisions because each surviving record has a unique prefix count. The proof uses the fact that prefix counts strictly increase at positions where the flag is one.

Scan performance involves more than theoretical depth. Each stage reads and writes arrays, and the algorithm may be bandwidth-bound. Temporary buffers, launch overhead and communication of block totals affect runtime. An optimized implementation may use warp-level primitives or fuse local scanning with later computation. Compare not only standalone scan throughput but also the surrounding pipeline: a filter-plus-scan-plus-scatter pipeline might benefit from processing tiles that stay in cache or shared memory. Saving one array pass can outweigh a marginal improvement in the scan's arithmetic operations.

Test with arrays of length zero, one, prime numbers, and lengths slightly below and above block boundaries. Check both inclusive and exclusive conventions with a serial reference. Add input patterns of all zeros, all ones and alternating flags to validate compaction. A scan is a foundational parallel algorithm precisely because its output structure exposes ordering and synchronization mistakes so clearly.

Sorting is a fundamental data transformation whose parallelization is more challenging than elementwise mapping. Comparison sorting must coordinate which values occupy which positions, and intermediate partitions may have very different sizes. Parallel merge sort divides the input into subarrays, sorts them concurrently and merges their results. Parallel quicksort partitions around pivots and recurses, but poor pivots can create highly uneven work. Sample sort selects splitters to partition records among processes, then locally sorts each partition. The best algorithm depends on whether keys are bounded integers or arbitrary comparable objects, on how large records are, and on the memory and communication hierarchy.

Consider parallel merge sort on n records. A serial merge-sort algorithm performs O(n log n) comparisons. Parallel recursive sorting exposes independent subproblems, but each merge can introduce a long dependency if implemented sequentially. A parallel merge uses binary searches or block partitions to identify independent output regions, reducing span at the cost of extra index work. The algorithm must preserve the ordering relation, including behavior when keys are equal. A stable sort retains input order among equal keys; parallel implementations sometimes sacrifice stability unless designed for it. If stability is required for a downstream database operation, this is a correctness requirement rather than an optional aesthetic feature.

Sorting distributed data usually involves communication. Suppose each of eight ranks holds a million records and the final order must be globally sorted by key. Sorting each rank's local array is insufficient, because a large key on rank zero may need to move after a small key on rank seven. A sample-sort strategy collects representative keys, selects splitters, partitions data by key ranges, exchanges partitions among ranks, and sorts locally. Skewed key distributions can overload a few destinations, so splitter selection and duplicate-key handling matter. If all keys are identical, a naive range partition might send nearly every record to one rank unless the algorithm deliberately balances equals.

Sorting network algorithms such as bitonic sort execute a fixed pattern of comparisons and swaps. Their predictable control flow can be attractive on hardware that benefits from regular operations, but they may perform more comparisons than asymptotically optimal comparison sorts. Bitonic networks naturally accommodate powers of two; a correct general implementation must pad safely or use a specialized irregular-size method. Network-based sorting is a good teaching example because every comparator has a clearly defined dependency level, making span and synchronization costs explicit.

Integer keys sometimes enable radix sorting. Rather than comparing arbitrary pairs, radix sorting processes a fixed number of bits or digits, building histograms and prefix offsets for each pass. A parallel radix sorter uses scans to calculate stable destination positions. Its complexity depends on key width, digit size and memory traffic. Larger digits reduce the number of passes but require larger histograms, creating different cache and shared-memory tradeoffs. A floating-point sort requires deliberate handling of sign representation, NaNs and signed zeros if bitwise ordering tricks are used.

Sorting benchmarks are particularly vulnerable to misleading input choices. Already sorted input, reverse order, many duplicates and uniform random keys can exercise very different paths. A quicksort implementation that performs well on random keys may degrade on adversarial pivots. A radix sort that excels on fixed-width integers may not apply to variable-length strings. Report key type, record size, distribution, stability requirements, device, total input length and memory limits. For a distributed sort, include bytes moved over the network and final output placement in the timing boundary.

Verification should go beyond checking that neighboring output keys are nondecreasing. A sorting algorithm could accidentally lose elements, duplicate values or corrupt associated payloads while still returning a sorted key sequence. Compare a multiset checksum or exact record counts, verify that the output is a permutation of the input and test equal-key stability when required. Small arrays with unique identifiers attached to duplicate keys are especially effective. On a multi-rank output, check global order across partition boundaries as well as order within each rank.

The parallel sorting lesson is about balancing comparison work, memory movement and global redistribution. An algorithm with an excellent asymptotic comparison bound can lose to a simpler method if it moves records inefficiently. The strongest design chooses its partitioning strategy from the data distribution and target memory hierarchy, then verifies both global order and complete preservation of records.

Dense matrix computation is important because it offers substantial arithmetic reuse when organized correctly. Multiplying an m-by-k matrix A by a k-by-n matrix B requires approximately 2mkn floating-point operations under the common convention that a multiply-add counts as two operations. The naive loop nest repeatedly loads elements whose values will be reused for many products. Blocking or tiling reorganizes the computation so a submatrix of A and a submatrix of B remain in faster memory while producing many output elements. This improves arithmetic intensity, the amount of useful computation per byte moved from slower memory.

A small matrix example makes the mapping explicit. For matrices A=[[1,2],[3,4]] and B=[[5,6],[7,8]], C=AB has entries 19,22,43,50. Each output element is a dot product of one row of A with one column of B, so four outputs can be computed independently if the inputs are immutable. For larger matrices, assign tiles of C to worker groups rather than one output scalar per expensive task. Each tile needs appropriate portions of A and B. A good implementation avoids two workers writing the same C tile unless a controlled reduction combines their partial results.

The classical three-loop multiplication has O(n cubed) arithmetic for square matrices, but performance depends strongly on loop order and storage layout. For row-major arrays, an inner loop that walks contiguous B elements can exploit cache and vectorization more effectively than one that repeatedly strides down columns. Packing blocks into contiguous scratch buffers can further improve locality. These transformations preserve the mathematical operation but may change floating-point reduction order. Numerical verification should account for that and include checks with known matrices, identities, zeros and non-square dimensions.

On a cluster, matrix multiplication requires distributed data placement. A two-dimensional process grid can assign blocks of C to ranks and communicate panels of A and B as needed. SUMMA-style algorithms use repeated broadcasts along process-grid rows and columns, with local matrix multiplications between communication phases. Performance depends on panel width, memory footprint and overlap. A naive strategy that sends all of A and B to every rank can exceed available memory and network bandwidth. Scalable algorithms strive to reuse communicated blocks, balance work and respect each rank's memory capacity.

GPU matrix kernels often use shared-memory tiles, register accumulators and specialized arithmetic units. Hardware matrix-multiply accelerators may support reduced-precision inputs with accumulation in a different precision. Their impressive nominal throughput does not guarantee acceptable accuracy for every scientific problem. An ill-conditioned linear algebra task may require higher precision or iterative refinement. The data type, accumulation method and tolerated error are part of the algorithm's contract, especially when results are used in numerical solvers rather than visualization.

Parallel matrix-vector multiplication offers a contrasting example. Multiplying a dense n-by-n matrix by a vector performs O(n squared) arithmetic and reads O(n squared) matrix entries. Unless the matrix or vector is repeatedly reused, the arithmetic intensity is much lower than matrix-matrix multiplication. Many optimized matrix-vector kernels are memory-bandwidth limited. This difference explains why simply placing both operations on a GPU yields very different scaling. The roofline framework predicts the sensitivity to bandwidth versus compute throughput by comparing operations to bytes moved.

Sparse matrices introduce irregular memory accesses. A compressed sparse row format stores nonzero values and column indices with row offsets. Different rows may contain very different numbers of nonzeros, producing load imbalance when one thread handles each row. Grouping long rows, using segmented reductions or reordering the sparse structure can improve balance, but preprocessing and storage overhead matter. For iterative solvers that reuse the matrix thousands of times, preprocessing may be worthwhile; for one multiplication it may not. Correctness checks must preserve the relation between stored values, column indices and row boundaries.

A robust matrix performance study separates transfer, packing, computation and numerical validation. Use dimensions not aligned to tile sizes, compare against a trusted small reference and measure sustained operations per second relative to relevant precision. A fast multiplication that silently discards small terms or uses an unreported approximation is not equivalent to a precise baseline. Parallel matrix algorithms succeed by coordinating mathematical decomposition with data movement, memory reuse and a transparent accuracy contract.