significant load imbalance across tiles, we use it here as an example to
illustrate how load balancing is defined within our framework.

The inputs are the three iterators from the last stage plus an atom and tile
count. The load-balance algorithm developer, then, implements tiles() and
atoms() procedure calls, which return the C++ range of tiles and atoms to be
processed by the current thread, effectively creating a map between assigned
processor ids and segments of the workload. Listing 3 shows a complete example
of the thread-mapped schedule. Although a simple algorithm, it can deliver
high performance for well-balanced workloads with coarse-grained parallelism
(a small number of atoms per tile), such as multiplying a sparse vector by a
dense vector. Furthermore, our abstraction is not limited to only simple
scheduling algorithms, as Section 2 provides examples of more complex load
balancing algorithms.

Figure 3: A thread-mapped load-balancing algorithm expressed as C++ ranges,
incorporating the atoms and tiles defined as iterators from Listing 2. Each
tile is mapped to a thread, where the thread id corresponds to the index of
the tile in the tile set. All atoms within a tile are sequentially processed
by the thread. After a tile is processed, a thread is mapped to the next tile,
obtained by striding the index by the grid size of the kernel.

⬇

class schedule_t {

// Construct a thread-mapped schedule.

__host__ __device__ schedule_t(atoms_it_t atoms_it,

tiles_it_t tiles_it,

atoms_it_t atoms_per_tile_it,

size_t num_atoms, size_t num_tiles) :

m_atoms_it(atoms_it), m_tiles_it(tiles_it),

m_atoms_per_tile_it(atoms_per_tile_it),

m_num_atoms(num_atoms),

m_num_tiles(num_tiles) {}

// Returns a range of tiles to process in "this" thread.

// Stride by grid dimension.

__host__ __device__ auto tiles() {

auto begin = m_tiles_it(blockDim.x * blockIdx.x + threadIdx.x);

auto end = m_tiles_it(m_num_tiles);

return range(begin, end).step(gridDim.x * blockDim.x);

}

// Returns a range of atoms to process in "this" thread.

__host__ __device__ auto atoms(

const std::size_t& tile) {

auto begin = m_atoms_per_tile_it[tile];

auto end = m_atoms_per_tile_it[tile + 1];

return range(begin, end).step(1);

}

};

using schedule_t = thread_mapped_schedule_t;

####  3 Implementing Work Execution

Our framework is designed to explicitly let the user _own_ the kernel launch
boundary. Owning a CUDA kernel boundary means that the user is responsible for
maintaining and configuring launch parameters and implementing the CUDA kernel
used to define the application. Although this design decision comes at a cost
of convenience and simplicity, it offers significant flexibility in what users
can express through our abstraction. This design decision is motivated by the
following reasons. (1) Users are not required to add a complex dependency to
their existing workflow/libraries, therefore making code maintenance simpler
and more scalable as they do not have to rely on our framework to incorporate
new CUDA constructs and features. (2) Users are free to express anything and
everything CUDA allows within their kernels while consuming our load-balanced
C++ ranges. This allows for versatility in what can be expressed, as the users
can now specify multiple load-balanced work domains, range-based for loops,
and even fusing multiple computations to build more complex algorithms within
a single kernel. (3) Higher-level APIs can be used to build simpler higher-
level abstractions that _do_ own the kernel boundary and provide simpler APIs
at the cost of flexibility.

As an input to this stage, users consume the load-balanced C++ ranges to
implement their computation. This can be done in multiple ways, but one of the
most common patterns is a nested range-based for loop that loops over all the
assigned tiles and atoms ranges. Listing 4 shows a simple example of a CUDA
kernel that implements the SpMV algorithm using CSR format and thread-mapped
load-balancing algorithm constructed in Listings 2 and 3. In this example, the
outer for loop within each thread iterates over the assigned rows of the
sparse matrix (tiles), and the inner loop sequentially processes the assigned
nonzeros (atoms) within each row. In Section 3 we implement and discuss more
complex kernels and computations.

Figure 4: Sparse-Vector Matrix Multiplication (SpMV) implemented within our
load-balancing abstraction using range-based nested for loops. The sparse
matrix is represented using a CSR-based format, where $x$ is the dense input
vector and $y$ is the dense output vector ($y=Ax$). Lines 11–14 use the load-
balancing schedule implemented in Listing 3 and the iterators defined in
Listing 2 to construct the load-balanced work to be processed. Lines 16 and 19
show the for loops within each thread, which iterate over the assigned rows of
the sparse matrix and sequentially process the assigned atoms within each row.
Line 20 shows the actual computation performed on each work atom (nonzero),
and Line 21 writes the result to the dense output vector $y$.

⬇

// Implements load-balanced SpMV kernel.

__global__ void spmv(const size_t rows,

const size_t cols,

const size_t nnz,

const int* offsets,

const int* indices,

const float* values,

const float* x, float* y) {

// Configure load-balancing.

// Input: iterators defined for CSR format.

schedule_t config(

atoms_iter, tile_iter,

atoms_per_tile_it,

nnz, rows);

// Consume rows using a range-based for loop.

for (auto row : config.tiles()) {

type_t sum = 0;

// Consume atoms using a range-based for loop.

for (auto nz : config.atoms(row))

sum += values[nz] * x[indices[nz]];

y[row] = sum;

}

}

// Launches SpMV kernel.

constexpr size_t blocks = 256;

size_t grid = (rows + blocks - 1) / blocks;

spmv<<<grid, blocks>>>(rows, cols, nnz, offsets, indices, values, x, y);

###  4 Implementation Details

####  1 Flexible, Composable CUDA-enabled Ranges

Composability of load-balanced primitives and applications using our API is a
conscious design choice within our framework supported through the use of
CUDA-enabled C++ ranges. Our framework does not _own_ the kernel boundary
(kernel launch), which forces our APIs to be focused and contained within the
kernels. This allows programmers to build and maintain their own kernels while
still benefiting from our framework’s load-balancing capabilities. This is
largely implemented using device-wide C++ functions and classes tagged with
CUDA’s __device__ keyword.111A method decorated with the __device__ keyword
allows the CUDA compiler to generate a device-callable entry point. This
allows the code to be called from within kernels [69]. We implemented and
expose several different types of specialized ranges that were particularly
useful in implementing load-balanced schedules:

  * •

step_range: A range that iterates from begin to end in steps of step. Useful
for defining load balancing schedules that require a custom stepping range or
process a constant number of work items per thread (which can be defined using
step).

  * •

infinite_range: A range that iterates from begin to infinity. Useful for
defining load balancing schedules in persistent kernel mode [90], where the
kernel persistently runs until all work is consumed or an algorithm has
converged.

  * •

grid_stride_range: A specialized case of step_range that iterates from begin
to end in steps of step using the CUDA kernel’s grid size. Also supports
block_stride and warp_stride variants that iterate in steps of the block or
warp size, respectively.

####  2 Implementing Non-Trivial Load-Balancing

As we describe in Section 1, we can decouple and express existing load-
balancing techniques as a set of C++ ranges. To illustrate the potential of
this abstraction, we begin by decoupling and expressing a state-of-the-art
load-balancing algorithm known as merge-path [35] previously used for
balancing CSR-based SpMV and SpMM [64, 88], and implement three additional
load balancing algorithms (warp-, block- and group-mapped), all of which are
available in our library for programmers to use. Our new group-mapped
algorithm is a tile-per-group-based schedule, where a group is defined as a
collection of threads of any arbitrary size (not limited to a warp or block
size). Our group-mapped schedule is a generalization of the tile-per-thread,
-warp or -block schedules [65, 14] using CUDA’s Cooperative Groups programming
model [40].

##### Merge-path load balancing

In the language of a sparse matrix, merge-path assumes that each non-zero in
the matrix and each new row in the matrix are an equivalent amount of work,
then evenly divides $\text{nnzs}+\text{rows}$ work across the set of worker
threads. Each thread then performs a 2-D binary search within the nonzero
indices and row offsets of a CSR matrix to find the starting position of the
row and nonzero it needs to process. Threads then sequentially process the
rows and nonzeros from the starting position until they reach the end of their
assigned work [64].

We implement this algorithm as a load-balancing schedule in our abstraction by
expressing it in two steps: (1) Setup: The initialization step of the C++
schedule class computes the number of work units per thread, conducts a binary
search as described above, and stores the starting position of each tile and
atom in a thread-local variable. (2) Ranges: The second step of the algorithm
builds the ranges for each thread to process as “complete” tiles and “partial”
tiles [64]. If a thread’s atom range lies entirely within one tile, it is
“complete”, and is processed in a simple nested loop. If a thread’s range
crosses a tile boundary, the thread processes its work in a separate nested
loop.

Because we decouple the load-balancing method (Section 2, and above) from the
work execution (Section 3), we can use this merge-path implementation to
implement not only SpMV but also any other algorithm whose work can be divided
into tiles and atoms, e.g., a graph neighborhood-traversal algorithm used to
implement breadth-first search [84]. Just as importantly, the merge-path
schedule is now no longer limited to a CSR-based sparse format. Supporting
other formats only requires building the necessary slightly more complex
iterators that are able to count atoms per tile (the computation that the CSR
implementation achieves with the row offsets array in Listing 2).

##### Warp- and block-level load balancing

The goal of a warp- or block-level load-balancing schedule is to assign an
equal share of tiles to each warp or block, which are then sequentially
processed. The work atoms within each tile will be processed in parallel by
the available threads within a warp or a block. Each thread strides by the
size of the warp or block to process a new work atom until the end of work is
reached.

The imbalance across different processing units is left for the hardware
scheduler to handle. This scheduler depends on the oversubscription model of
CUDA, where the programmer can launch a larger number of warps or blocks than
the GPU can physically schedule at any given time. As the processing units
finish processing their work, new ones are scheduled from the oversubscribed
pool [65, 14].

##### Group-level load balancing

Group-level load balancing generalizes warp- and block-level schedules.
Instead of requiring that group sizes are the size of a warp or block, as
above, this method leverages CUDA’s Cooperative Groups (CG) programming model
[40] to allow programmer-specified dynamically sized groups of arbitrary size.
Within these groups, the CG model permits detailed control of the group’s
synchronization behaviors as well as simple parallel group-level collectives
such as reduce or scan. We leverage this powerful tool to implement a
generalized group-level load balancing schedule, effectively giving us the
warp- and block-level schedules above for free when the group size equals that
of a warp or a block.

Our schedule assigns work tiles to a group, and each group looks at its equal
share of tiles and computes the number of atoms for each tile and stores it in
a scratchpad memory (CUDA’s _shared memory_). The group then performs a
parallel prefix-sum, a widely used parallel algorithm that inputs an array and
produces a new array where the element at any position is a sum of all
previous elements [9]. We use this prefix-sum array for two purposes: (1) the
last element of a prefix-sum array indicates the aggregated number of work
atoms that a group has to process, and (2) the position of each sum in the
prefix-sum array corresponds to the work tile to which those atoms belong. The
setup phase of the schedule builds the prefix-sum arrays per group in the
scratchpad memory, and the ranged-loop of the schedule returns the atom to
process in each thread. The corresponding tile, if needed, is obtained by a
simple get_tile(atom_id) operation, which executes a binary search within the
prefix-sum array to find the tile corresponding to the atom being processed.

Relying on the CG model for this load-balancing schedule has a unique
advantage of configuring the group size (effectively software constructs that
directly map onto the hardware) per the shape of the problem and the
underlying hardware architecture. For example, targeting GPUs where the warp
size is not 32 threads (AMD’s GPU architecture supports a warp size of 64 [3])
is now possible with a simple compile-time constant, or configuring the group
size to perfectly align with the structure of the problem.

####  3 Application Space

Our work definition (Section 1), composable APIs (Section 1), and multiple
sophisticated, high-performance load-balancing schedules (Section 2) together
provide for a versatile and extensible framework with plenty of room for
application-specific optimizations. In Listing 4 we already demonstrated how
to implement the SpMV algorithm using our framework. A simple and natural
extension is to implement Sparse-Matrix Matrix Multiplication (SpMM). Listing
5 shows the minor change necessary, which adds another loop over the columns
of the B matrix around the existing code from Listing 4 to implement SpMM.
This implementation could also be extended to support Gustavson’s General
Sparse Matrix-Matrix Multiplication (SpGEMM), using two kernels and an
allocation stage; the first kernel would compute the size of the output rows
used to allocate the memory for the output sparse matrix and the second kernel
would perform the multiply-accumulation.

Figure 5: A simple loop wrapped around SpMV introduced in Listing 4 allows us
to represent the slightly more complex SpMM load-balanced computation.

⬇

// ... Inside the CUDA kernel.

// Loop over all the assigned rows.

for (auto row : config.tiles()) {

// Loop over all the columns of Matrix B.

for (auto col : range(size_t(0), B.cols)

.stride(size_t(1))) { /// < New Loop

float sum = 0;

// Loop over all the assigned nonzeros.

for (auto nz : config.atoms(row))

sum += values[nz] * B(nz, col);

// Output the sum to Matrix-C.

C(row, col) = sum;

}

}

Beyond sparse linear algebra, we can use our framework to address applications
in other domains. Listing 6 implements the graph primitive Single-Source
Shortest Path (SSSP) using our group-level load-balancing schedule. SSSP’s
performance on GPUs is largely gated by good load balancing [84, 14], but if
the programmer chooses a load-balancing schedule from our library, the details
of load balancing are completely hidden. Moreover, the same schedules that
were used in one application domain (e.g., sparse linear algebra) are easily
reusable in this different application domain.

Figure 6: The parallel single-source shortest path (SSSP) graph primitive
expressed using our load-balanced schedule.

⬇

// ... Inside the CUDA kernel.

// Loop over all the assigned edges to process.

for (auto edge : config.atoms()) {

auto source = config.get_tile(edge);

// G is the graph data structure

auto neighbor = G.get_neighbor(source, edge);

auto weight = G.get_edge_weight(edge);

float source_dist = dist[source];

float neighbor_dist = source_dist + weight;

// Check if the destination node has been

// claimed as someone’s child.

float recover_distance =

atomicMin(&(dist[neighbor]), neighbor_dist);

// Add the neighbor to the frontier.

if (neighbor_dist < recover_distance)

out_frontier[neighbor] = true;

}

\par// ... Outside the CUDA kernel.

// Loop until the frontier is empty.

###  5 Evaluation

We aim to show that our framework, built on our load balancing abstraction,
enables both high performance and better programmability for sparse-irregular
problems. Our evaluation below uses our SpMV implementation as a benchmark
against state-of-the-art implementations provided within NVIDIA’s (open-
source) CUB library and production (closed-source) cuSparse library. We
considered (and implemented) several additional applications for evaluation,
including SSSP, BFS, and SpMM. We found they led to similar high-level
conclusions. Thus our evaluation here focuses on SpMV. Our test corpus
consists of approximately the _entire_ SuiteSparse Matrix Collection [29] with
a broad scope of sparse matrices from many different high-performance
computing domains. We ran all experiments on a Ubuntu 20.04 LTS-based
workstation with an NVIDIA Tesla V100 GPU and CUDA 11.7.

####  1 Performance Overhead

Figure 7: Load balancing abstraction performance overhead analysis using our
Merge-path SpMV vs. CUB across all SuiteSparse datasets. Our runtimes almost
perfectly match CUB’s for all datasets. The small number of datasets where CUB
is faster is due to a simple heuristic that CUB uses for sparse matrices with
the number of columns equal to 1 (i.e., a vector).

Our first and foremost goal is to ensure that the elements within our
abstraction do not add any additional performance overhead to the existing
load balancing techniques and algorithms developed using them. To verify this,
we compare the runtime performance of our SpMV implementation using the merge-
path schedule to the implementation provided by NVIDIA’s CUB library [63]
(also used for Merrill and Garland’s merge-path SpMV paper [64]) on the
SuiteSparse collection. As previously mentioned, and in contrast to our
design, CUB contains a hardwired implementation of the merge-path scheduling
algorithms and does not decouple workload balancing from the actual SpMV
computation. CUB’s approach is not reusable for any other irregular parallel
problem without significant changes to the implementation.

Figure 7 plots the number of nonzeros (i.e., the total work) vs. runtime for
our work vs. CUB’s implementation. Our implementation has minimal performance
overhead when using our abstraction: a geomean slowdown of 2.5% vs. CUB, with
92% of datasets achieving at least 90% of CUB’s performance. Figure 7 shows
our implementation almost perfectly matches CUB for all datasets, except for
some datasets with fewer than 100,000 nonzeros. Upon further investigation, we
identify that CUB uses a simple heuristic to launch a thread-mapped SpMV
kernel where the number of columns of a given input matrix equals 1 (i.e., a
sparse vector). Unlike our more general implementation, CUB’s simple (but
specialized) thread-mapped SpMV kernel has no load-balancing overhead for a
perfectly balanced workload such as SpVV computation.

####  2 Improved Performance Response

Figure 8: Complete performance landscape of SpMV across all SuiteSparse
datasets using 3 load balancing schedules vs. NVIDIA’s cuSparse library. Our 3
different SpMV implementations are made possible with very little code change.

We also compare our work to NVIDIA’s vendor library for sparse computations,
cuSparse. Figure 8 shows the performance response of our SpMV implementation
using each of our scheduling algorithms individually vs. cuSparse’s state-of-
the-art implementation. Switching between any of our implementations requires
very little code change; in the case of merge-path and thread-mapped, we need
only update a single C++ enum (identifier) to select the desired load
balancing schedule.

Figure 9: Speedup of our framework’s SpMV vs. cuSparse’s SpMV across
SuiteSparse using a heuristic (Section 2) to choose the appropriate load-
balancing schedule.

We then combine our scheduling algorithms into one implementation for SpMV,
demonstrating noticeable performance improvements over cuSparse (Figure 9).
This is primarily possible due to our ability to quickly experiment with
different heuristic schemes with a variety of available load-balancing
schedules. Here, we use merge-path unless either the number of rows or columns
are less than the threshold $\alpha$ and the nonzeros of a given matrix are
less than threshold $\beta$ (we choose $\alpha=500$ and $\beta=10000$ for
SuiteSparse). In this case, we use thread-mapped or group-mapped load
balancing instead of merge-path. Our system shows a peak performance speedup
of 39$\times$ and a geomean performance speedup of 2.7$\times$ vs. cuSparse.

Our framework not only allows programmers to express computations efficiently
and simply (i.e. without worrying about the load-balancing algorithms), but
also quickly optimize a given application using a range of scheduling
algorithms, both with minor code changes.

####  3 Lines of Code (LOC)

We are able to achieve these performance gains with minimal code complexity.
Table 1 shows lines of code (LOC) for our framework when compared to the
state-of-the-art open-source implementation of merge-path and thread-mapped
within NVIDIA’s CUB library. We deliver the same performance results as
highlighted in the previous section with 14$\times$ and 1$\times$ fewer lines
of code for merge-path and thread-mapped scheduling algorithms, respectively.
Using our merge-path implementation only requires $\sim$15 additional LoC to
the trivial thread-mapped schedule.

Furthermore, we extend the same SpMV computation to our novel _group-mapped_
load balancing schedule (that can also be specialized to perform _block-_ and
_warp-mapped_ load balancing) within the same 30 LoC.

Load Balancing Algorithm | NVIDIA/CUB | Our Work
---|---|---
Merge-Path | 503 | 36
Thread-Mapped | 22 | 21
Group-Mapped | N/A | 30
Warp-Mapped | N/A | 30 (free)
Block-Mapped | N/A | 30 (free)

Table 1: Lines of code (LoC) comparison for NVIDIA’s CUB library versus our
work for SpMV application implemented using merge-path, thread-mapped and
group-mapped (warp- and block-mapped use the _exact_ same code for group-
mapped) load balancing algorithms. We report only non-commented lines of code,
formatted using clang-format tool with the Chromium style guide that
contributes to the kernel implementation [34].

###  6 Related Work

Load balancing is the key to achieving high performance on GPUs for sparse,
irregular parallel problems. Several high-performance computing applications
deploy sophisticated load balancing algorithms on the GPUs. For instance,
high-performance sparse-matrix vector multiplication (SpMV) leverages merge-
path [64] (discussed in detail in this paper) and nonzero splitting algorithm,
which partitions the number of non-zeros in a sparse-matrix evenly across the
number of threads [8, 26, 78]. Sparse-matrix matrix multiplication (SpMM) and
sparse matricized tensor times Khatri-Rao product (SpMTTKRP) use binning and
bundling algorithms [88, 33, 67], which attempt to bin like-length work
together such that they are processed together.

While some applications actively perform work to load balance a given input,
others store the input in more efficient, already-load-balanced/partitioned
formats. These include the F-COO format (a variant of coordinate format) used
for SpMTTKRP and Sparse-Tensor Tensor Multiplication (SpTTM), where each
thread gets the same number of nonzeros to process [56].

Many of the above GPU load-balancing algorithms, along with other novel
techniques, were first described in the graph analytics domain. Davidson et
al. and Merrill and Garland were the first to present Warp, Block-level and
Thread-Warp-CTA dynamic load balancing techniques for Single-Source Shortest
Path (SSSP) and Breadth-First Search (BFS) respectively [28, 65]. Logarithmic
Radix Binning (LRB) is a particularly effective technique for binning work
based on a logarithmic work estimate, used for the Triangle Counting graph
algorithm and more [36, 32]. Gunrock, GraphIT, and GraphBLAST are graph
analytics libraries that implement several different graph algorithms such as
BFS, SSSP, PageRank, Graph Coloring, and more, built on these previously
mentioned load-balancing techniques [84, 14, 89]. Although many of these are
effective load balancing techniques with high-performance implementations,
they all tightly couple workload scheduling with the application itself. Our
framework is designed to separate these two concerns, allowing the application
to be independent of the load-balancing algorithm, and therefore be expressed
simply. Our approach also allows these previously proposed techniques to be
implemented within our framework, and be used for applications beyond those
originally targeted.

Relatively few GPU works target generalized load balancing for irregular
workloads. Most of these are focused on providing a singular, dynamic load-
balancing solution centered on task parallelism, often using a GPU queue-based
data structure. Cederman and Tsigas proposed a task-based approach to load
balancing an octree partitioning workload using lock-free and lock-based
methods [18]. Two Tzeng works provide task-management frameworks that
implement load balancing of tasks using a single monolithic task queue and
distributed queues with task stealing and donation [82, 83]. CUIRRE, a
framework for load balancing and characterizing irregular applications on
GPUs, also uses a task-pool approach [92], and more recently, Atos, a task-
parallel GPU dynamic scheduling framework, targets asynchronous algorithms
[22]. All of these works deploy either a centralized or a distributed queue-
like data structure on the GPUs, each making design decisions on how the queue
is to be partitioned and updated. Except for the most recent Atos work, most
earlier works focus on a coarse-grained parallelism approach of effectively
distributing tasks to the GPU. Our work takes advantage of more modern GPU
architectures, which are more effectively utilized by a fine-grained
parallelism approach (parallelizing over work atoms instead of work tiles).
Unlike our abstraction, these aforementioned works also rely on a singular
load balancing solution, whereas our abstraction flexibly adapts to many
different load balancing techniques, static and dynamic, and allows for new
schedules to be implemented within our framework.

###  7 Conclusion

In this paper, we present a programming model for GPU load balancing for
sparse irregular parallel problems. Our model is built on the idea of
separation of concerns between workload mapping and work execution. In the
future, we are interested in expanding our model to a multi-GPU environment,
and implementing load balancing schedules that span across the GPU boundary
covering multiple devices and nodes for massive parallel problems. Our current
work focuses solely on load balancing, but we also identify locality to be
another key factor for high performance. We are interested in identifying an
orthogonal model that builds an abstraction for caching and locality into our
existing load balancing framework.

In part, this chapter has shown the importance of load-balancing to the
performance of irregular problems. While regular problems typically map in a
more straightforward way to GPUs, and thus suffer from fewer load-balancing
issues, the next chapter shows that even dense, arithmetically intense,
operations such as matrix-matrix product (GEMM) can suffer from poor resource
utilization due to load imbalance. However, whereas load imbalance for sparse
problems typically arises from variable-length data and unpredictable access
patterns, load imbalance in GEMM arises from work decompositions that are a
poor fit to the physical resources of the architecture. In the following
chapter we introduce _Stream-K_ , a work-centric parallel decomposition for
dense operations that addresses the resource imbalance issues.

##  Chapter 4 Work-centric Parallel Decomposition for Dense Matrix-Matrix
Multiplication on the GPU

General matrix-matrix product (GEMM), convolution, and other similar
computations constitute the dominant workloads in many deep learning and
scientific computing applications. High-performance processors such as GPUs,
for example, are designed to achieve nearly 100% of their theoretical peak
math throughput when computing GEMM. Doing so, however, requires a work
decomposition that perfectly occupies the underlying physical cores. As we
show, attaining such high levels of processor utilization across a broad
landscape of problems shapes and sizes can be challenging.

Classically, GEMM implementations block their computation using a _data-
parallel_ tiling of the output matrix, assigning the independent production of
output tiles among concurrent threads (or thread groups) [2, 49, 66]. The work
per output tile is regular, and tile production tends to dispatch across idle
physical cores in “waves”. The overall workload is well-balanced and processor
utilization is highest when there are many waves, i.e., the number of output
tiles greatly oversubscribes the number of cores.

However, such oversubscription has shrunk considerably as processors have
grown in size. An increased core count will require fewer waves to produce a
given tile count. Bigger cores will compel larger matrix blocking factors,
leading to fewer waves of larger tiles. In general, execution schedules with
fewer waves are much more likely to suffer from _quantization inefficiency_ ,
i.e., the processor underutilization that occurs when the number of output
tiles is not an even multiple of the number of processor cores. When the last
wave is partially full, the unused cores must wait for the remaining threads
to execute millions (if not billions) of multiply-accumulate (MAC)
instructions before they are able to execute any dependent work.

Figure 1(a) illustrates such a scenario on a hypothetical GPU with four
streaming multiprocessor cores (SMs). If we block a $384\times 384\times 128$
GEMM computation into nine $128\times 128$ output tiles, a _data-parallel_
decomposition cannot achieve more than 75% of the processor’s rated
throughput. This theoretical utilization ceiling can be improved to 90% by
halving the tile size as shown in Figure 1(b). However, the finer-grained
blocking factor will be less cache and scratchpad efficient, and may preclude
any practical performance improvement.

(a) Data parallel decomposition with grid size $g$=9 CTAs,
large $128\times 128\times 128$ CTA work volumes,
and 75% processor utilization ceiling

(b) _Data parallel_ decomposition with grid size $g$=18 CTAs,
smaller $128\times 64\times 128$ CTA work volumes,
and 90% processor utilization ceiling

Figure 1: _Data-parallel_ execution schedules for $384\times 384\times 128$
GEMM across a hypothetical four-SM GPU.

(a) _Fixed-split_ decomposition with splitting factor $s$=2,
grid size $g$=18 CTAs, smaller $128\times 128\times 64$ CTA work volumes,
and 90% quantization efficiency

(b) Basic _Stream-K_ decomposition with grid size $g$=4 CTAs,
larger $128\times 128\times 288$ CTA work volumes,
and nearly 100% quantization efficiency

Figure 2: Tile-splitting execution schedules for $384\times 384\times 128$
GEMM across a hypothetical four-SM GPU.

Quantization inefficiency is a concern for increasingly wide processors such
as GPUs, where ALUs-per-core and cores-per-processor both currently number in
the hundreds. Consequently, many common GEMM-like workloads now exhibit a
final, partially full wave that comprises a significant fraction of the total
computation time.

The current remedy employed by GPU-based math and deep learning libraries is
to deploy an ensemble of tiling configurations. When the ideal blocking factor
does not quantize well, the library chooses among tiling alternatives with
smaller concurrent work volumes, such as those illustrated in Figure 1(b) and
Figure 2(a).

Tile-based ensembles, however, present performance and logistical challenges
for math libraries seeking to deliver the best-achievable performance across
diverse problem sizes and shapes. Distributable code size can be problematic
for large ensembles. For example, NVIDIA’s cuBLAS library [72] is hundreds of
megabytes, often providing more than twenty pre-compiled kernel
specializations per architecture for a given API entry point. Large ensembles
also require sophisticated selection heuristics. In our evaluation, we show
these heuristics can struggle to consistently identify the optimal
configuration for arbitrary problems.

Unlike these tile-based methods, our _Stream-K_ decomposition always
distributes an even share (within one) of the aggregate multiply-accumulate
loop iterations required by the GEMM computation across SMs. Because the
instruction workload of a single MAC-loop iteration is far smaller than that
of an entire output tile, any variance in core workload is practically
negligible. _Stream-K_ uses the ideal blocking factor regardless of problem
shape, has communication overheads that scale with processor width (rather
than output tiles), and compiles to a single kernel.

We use an enormous corpus of 32,824 GEMM shapes and sizes to evaluate _Stream-
K_ , which we implemented within NVIDIA’s CUTLASS library [49]. In comparison
with CUTLASS’s _data-parallel_ implementation of the same blocking factor,
_Stream-K_ provides a substantially higher performance response across our
landscape of GEMM problems, demonstrating up to 14$\times$ speedup on NVIDIA
A100 GPUs.

To highlight the practical challenges of ensemble-based solutions, we also
evaluate NVIDIA’s cuBLAS library as well as an oracle-driven ensemble of
_data-parallel_ CUTLASS tilings. Relative to both ensembles, we show that our
single-kernel _Stream-K_ achieves both (1) higher average performance, and (2)
higher performance consistency. Versus cuBLAS, _Stream-K_ demonstrates up to
6.7$\times$ speedup and virtually no instances of slowdown for compute-bound
problems.

###  1 Background

General Matrix Multiplication (GEMM) is defined as the product
$\textbf{C}=\alpha\textbf{AB}+\beta\textbf{C}$ where $\alpha$ and $\beta$ are
scalar values and A, B, and C are matrices. (For simplicity, we assume
$\alpha=1$, $\beta=0$ throughout this paper.) We refer to the _shape_ of a
given GEMM problem by the volumetric extents of its computation. For example,
a $m\times n\times k$ GEMM consumes $m\times k$ and $k\times n$ input matrices
A and B, respectively, performs $m\times n\times k$ multiply-accumulate
operations, and produces an $m\times n$ output matrix C.

GEMM is a performance-critical subroutine in many large-scale engineering and
scientific applications. It plays an important role in matrix factorization
methods such as LU, QR, and Cholesky decomposition. High-performance modeling
and simulation applications in engineering, climate simulation, cosmology,
quantum chemistry, and other scientific domains rely on these factorization
methods.

Matrix multiplication is also the fundamental building block of modern deep
learning (DL) methods. The training of deep neural networks (DNNs) is often
performed on massive datasets across large distributed systems [61]. Many DL
training and inference operations are cast as matrix multiplications. For
example, image recognition and computer vision models rely on convolution,
which can be implemented directly as the product of filter and image datasets
[23]. Transformer architectures, which have come to dominate natural language
processing and other applications, are almost entirely limited by the
performance of large matrix products.

Early work on GPU matrix-matrix multiplication from Larsen and McAllister
framed the computation as a multi-texture multiplication and blending
operation [53]. The user programmable shared memory provided by subsequent GPU
architectures enabled higher-performing _data parallel_ schemes with two
levels of blocking (shared memory and registers) with tile sizes informed via
extensive microbenchmarking analysis [7, 66, 79, 81] and auto-tuning [25, 48,
55].

The MAGMA GPU math library was perhaps the first to optimize for diverse GEMM
problem shapes [51]. Their solution applied a constrained set of tiling
parameters to a templated CUDA C++ code stencil, generating several hundred
_data-parallel_ variants per API primitive (e.g., hgemm_tt() for half-
precision transpose-transpose GEMM). They evaluated these variants to distill
a small ensemble of typically three to five kernels that collectively perform
well across a diversity of problem shapes. Kernel selection and dispatch for a
given problem was governed by size thresholds expressed via simple handwritten
rules.

Subsequent GPU math libraries have employed more sophisticated code-generation
and kernel-selection components. For example, the ISAAC project uses machine
learning techniques to predict an optimal tiling and/or splitting
parameterization for a given GEMM shape, which can then be instantiated either
online or offline via a PTX-level code generator [80].

NVIDIA’s cuBLAS [72] library has provided an extended `cublasGemmEx` interface
that allows the caller to select from among 24 different GEMM “algorithms”.
Carefully trained heuristics choose between this large space of alternatives
when using the default interface. These algorithms implement a variety of
different _data-parallel_ and _fixed-split_ variants, and it is common for
cuBLAS to have assembled each variant into its own architecture-specific
kernel program for code optimization purposes. The cross product of GEMM API
functionality, strategic variants, and microarchitecture has resulted in
distributions that are increasingly enormous, exceeding hundreds of megabytes
of executable code.

Given the fast-paced and rapidly changing nature of contemporary deep
learning, recent work has focused on programming models for simplifying the
expression and construction high performance kernels that alter or supplement
the GEMM computation. The CUTLASS C++ library provides data-movement and
multiply-accumulation classes for composing custom GEMM-like computations at
all levels of the GPU thread hierarchy [49]. Triton [81] is a domain-specific
language for tensor programming centered on the expression, transformation,
and optimization of block/tile concepts. Other domain-specific programming
languages such as Halide [76] and TVM [21] separate the expression of
pointwise operators from that of loop scheduling. Fireiron [39] further adds
data movement constructs into the scheduling grammar.

###  2 Work Decomposition Strategies

Modern processors typically store A, B, and C in a large, slow, distant memory
and have access to a small, fast, scratchpad or cache memory. A primary goal
for any GEMM implementation is to leverage these local storage resources so
that the resulting implementation is computation-bound.

####  1 Sequential Cache-Blocked

The classic cache-blocked formulation of GEMM divides its computational volume
into blocks and chooses a traversal order that exposes memory locality.
Algorithm 6 presents a simplified implementation comprising six loops. The
innermost three loops iterate within the blocking factors BLK_M, BLK_N, and
BLK_K, while the outermost three iterate across them. If the cache can capture
one block from each of the three matrices, the resulting data reuse among
those elements will significantly reduce the number of last-level memory
accesses [52].

Algorithm 6 Sequential cache-blocked GEMM.

inconsolata

1:$\rhd$ tile-processing outer loops

2:for mm $\leftarrow$ 0 to m step BLK_M do

3: for nn $\leftarrow$ 0 to n step BLK_N do

4: $\rhd$ zero-initialize output tile

5: for mmm $\leftarrow$ mm to (mm + BLK_M) do

6: for nnn $\leftarrow$ nn to (nn + BLK_N) do

7: C[mmm,nnn] $\leftarrow$ 0

8: end for

9: end for

10: $\rhd$ perform the MAC iterations for this tile

11: for kk $\leftarrow$ 0 to k step BLK_K do

12: $\rhd$ MAC iteration (fully unrolled)

13: for mmm $\leftarrow$ mm to (mm + BLK_M) do

14: for nnn $\leftarrow$ nn to (nn + BLK_N) do

15: for kkk $\leftarrow$ kk to (kk + BLK_K) do

16: C[mmm,nnn] $\leftarrow$ C[mmm, nnn] +

17: (A[mmm,kkk] $\times$ B[kkk,nnn])

18: end for

19: end for

20: end for

21: end for

22: end for

23:end forinconsolatait

####  2 Data-parallel

As shown in Algorithm 7, the _data-parallel_ GPU formulation of GEMM is
decomposed across a grid of parallel thread blocks, or _cooperative thread
arrays_ (CTAs)111Blocks of GPU threads are coscheduled in CTAs, which
virtualize the hardware’s streaming multiprocessor cores (SMs).. The grid is
sized such that each CTA produces its own (BLK_M $\times$ BLK_N) output tile.

For exposition, the MacLoop() subroutine of Algorithm 8 encapsulates the
multiply-accumulate (MAC) workloads that compute the values of the CTA’s
output tile. It performs a sequence of _MAC-loop_ iterations in the
accumulation domain, e.g., the _k_ -axis for GEMM. Each _MAC-loop_ iteration
comprises a per-thread volume of (BLK_M $\times$ BLK_N $\times$ BLK_K) $/$
CTA_THREADS MAC operations. As the computation proceeds, fragments of the
input matrices are staged through the SM’s shared memory for local reuse among
individual threads.

Although this particular presentation of MacLoop() deploys one thread per
output tile element, the sophisticated implementations in CUTLASS [49] and
cuBLAS [49] will: (1) fully unroll the per-thread MAC-loop iteration; (2)
implement additional blocking at the warp and/or thread levels; and (3)
orchestrate a software pipeline of shared memory data movement across MAC-loop
iterations.

Unfortunately, this classic _data-parallel_ decomposition is liable to suffer
from quantization inefficiency on modern GPUs, as illustrated in Figure 1.
Although an ensemble of diverse blocking factors may uncover opportunities for
greater processor utilization, it is unlikely to facilitate perfect
quantizations for arbitrary problem sizes. Furthermore, smaller blocking
factors have two drawbacks: (1) fewer instructions per MAC-loop iteration for
covering the latencies of global and shared memory transfers in pipelined
implementations; and (2) a higher proportion of memory operations relative to
MAC instructions, which may prevent them from being computation-bound.

Algorithm 7 Data-parallel GPU GEMM.

inconsolata

1:_shared_ accum[BLK_M,BLK_N]

2:iters_per_tile $\leftarrow\lceil$k/BLK_K$\rceil$

3:$\rhd$ instantiate one CTA per output tile

4:fork CTA[x] in [ $\lceil$m/BLK_M$\rceil$ $\times$ $\lceil$n/BLK_N$\rceil$ ]
do

5: $\rhd$ perform the MAC iterations for this tile

6: accum $\leftarrow$ MacLoop(x, 0, iters_per_tile)

7: $\rhd$ store accumulators to output tile

8: StoreTile(C, x, accum)

9:joininconsolatait

Algorithm 8 CTA-wide MacLoop() subroutine for performing a sequence of MAC-
loop iterations.

inconsolata

1:procedure MacLoop(tile_idx, iter_begin, iter_end)

2: _shared_ accum[BLK_M,BLK_N]

3: _shared_ frag_a[BLK_M,BLK_K]

4: _shared_ frag_b[BLK_K,BLK_N]

5: $\rhd$ determine output tile coordinates

6: mm $\leftarrow$ BLK_M $\times$ (tile_idx / $\lceil$m/BLK_M$\rceil$)

7: nn $\leftarrow$ BLK_N $\times$ (tile_idx % $\lceil$m/BLK_M$\rceil$)

8: $\rhd$ zero-initialize local accumulator storage

9: accum $\leftarrow$ 0

10: $\rhd$ perform the specified range of MAC iters for this tile

11: for iter $\leftarrow$ iter_begin to iter_end do

12: kk $\leftarrow$ iter $\times$ BLK_K

13: $\rhd$ copy global matrix fragments to local storage

14: frag_a $\leftarrow$ LoadFragment(A, mm, kk)

15: frag_b $\leftarrow$ LoadFragment(B, kk, nn)

16: fork THREAD[mmm,nnn] in [BLK_M, BLK_N] do

17: $\rhd$ MAC iteration per thread (fully unrolled)

18: for kkk $\leftarrow$ 0 to BLK_K do

19: accum[mmm, nnn] $\leftarrow$ accum[mmm,nnn] +

20: (frag_a[mmm,kkk] $\times$ frag_b[kkk,nnn])

21: end for

22: join

23: end for

24: return accum

25:end procedureinconsolatait

####  3 Fixed-split

Alternatively, the granularity of work assigned to each CTA can be reduced via
parallelization across the accumulation dimension. For a given output tile,
the associativity of addition allows the iteration domain to be split among
multiple concurrent CTAs, followed by a dependent “fixup” step to reduce the
partial sums computed by each CTA. We highlight this _fixed-split_ approach in
Algorithm 9, where each output tile is cooperatively produced by $s$ CTAs.
Notably, it functions identically to the _data-parallel_ decomposition when
the splitting factor $s=1$.

The _fixed-split_ decomposition is also featured in CUTLASS and cuBLAS. The
splitting factor is implemented as a runtime parameter, allowing a single
kernel executable to support multiple work volumes while retaining the ideal
blocking factors for optimal data sharing and latency hiding. However, as
illustrated in Figure 2(a), the prospect of achieving a perfect quantization
from a uniform tile-splitting is unlikely. Furthermore, the extra overheads of
communication and synchronization scale with both the overall problem size as
well as the splitting factor.

Algorithm 9 Fixed-split GPU GEMM with splitting factor $s$.

inconsolata

1:_shared_ accum[BLK_M,BLK_N]

2:iters_per_tile $\leftarrow\lceil$k/BLK_K$\rceil$

3:iters_per_split $\leftarrow\lceil$iters_per_tile/$s$$\rceil$

4:$\rhd$ instantiate $s$ CTAs per output tile

5:fork CTA[x,y] in [ $\lceil$m/BLK_M$\rceil$ $\times$ $\lceil$n/BLK_N$\rceil$,
$s$] do

6: $\rhd$ perform the range of MAC iterations for this split

7: iter $\leftarrow$ y $\times$ iters_per_split

8: iter_end $\leftarrow$ min(iters_per_tile, iter + iters_per_split)

9: accum $\leftarrow$ MacLoop(x, iter, iter_end)

10: $\rhd$ consolidate partial-sums across CTAs

11: if y $\neq$ 0 then

12: $\rhd$ store accumulators to temporary global storage

13: StorePartials(partials[x,y], accum)

14: Signal(flags[x,y])

15: else

16: $\rhd$ accumulate partial sums from other CTAs contributing to this tile

17: for cta $\leftarrow$ 1 to $s$ do

18: Wait(flags[x,cta])

19: accum $\leftarrow$ accum + LoadPartials(partials[x,cta])

20: end for

21: $\rhd$ store accumulators to output tile

22: StoreTile(C, tile_id, accum)

23: end if

24:joininconsolatait

####  4 Stream-K

Our _Stream-K_ decomposition is a tile-splitting parallelization in which the
splitting seams are completely dissociated from the tiling structure itself.
Although we employ familiar blocking and tiling strategies for data reuse, we
instead quantize the GEMM computation into MAC-loop iterations, i.e., small
volumes of CTA-wide BLK_M $\times$ BLK_N $\times$ BLK_K work. As presented in
Algorithm 10, _Stream-K_ evenly partitions the GEMM’s aggregate workload of
MAC-loop iterations across a constant-sized grid of $g$ CTAs. Each CTA’s range
of MAC-loop iterations is mapped contiguously into the $m\rightarrow
n\rightarrow k$ linearization of the GEMM shape, crossing output-tile
boundaries as it may.

Algorithm 10 Basic Stream-K GPU GEMM with grid size $g$.

inconsolata

1:_shared_ accum[BLK_M,BLK_N]

2:iters_per_tile $\leftarrow\lceil$k/BLK_K$\rceil$

3:total_iters $\leftarrow\lceil$m/BLK_M$\rceil\times\lceil$n/ BLK_N$\rceil$
$\times$ iters_per_tile

4:iters_per_cta $\leftarrow\lceil$total_iters / g$\rceil$

5:$\rhd$ instantiate g CTAs

6:fork CTA[x] in [g] do

7: iter $\leftarrow$ x $\times$ iters_per_cta

8: iter_end $\leftarrow$ iter + iters_per_cta

9: $\rhd$ iteration-processing outer loop

10: while iter $<$ iter_end do

11: tile_idx $\leftarrow$ iter / iters_per_tile

12: tile_iter $\leftarrow$ tile_idx $\times$ iters_per_tile

13: tile_iter_end $\leftarrow$ tile_iter + iters_per_tile

14: $\rhd$ perform the range of MAC iterations for this tile

15: local_iter $\leftarrow$ iter - tile_iter

16: local_iter_end $\leftarrow$ min(iter_end, tile_iter_end) - tile_iter

17: accum $\leftarrow$ MacLoop(tile_id, local_iter, local_iter_end)

18: $\rhd$ consolidate partial-sums across CTAs

19: tile_started $\leftarrow$ iter = tile_iter

20: tile_ended $\leftarrow$ (iter_end $\geq$ tile_iter_end)

21: if $\neg$tile_started then

22: $\rhd$ store accum to temporary global storage

23: StorePartials(partials[x], accum)

24: Signal(flags[x])

25: else

26: $\rhd$ store accumulators to output tile

27: if $\neg$tile_ended then

28: $\rhd$ accumulate partial sums from other CTA

29: $\rhd$ contributing to this tile

30: cta_end $\leftarrow$ tile_iter_end / iters_per_tile

31: for cta $\leftarrow$ (x+1) in cta_end do

32: Wait(flags[cta])

33: accum $\leftarrow$ accum + LoadPartials(partials[cta])

34: end for

35: end if

36: StoreTile(C, tile_id, accum)

37: end if

38: iter $\leftarrow$ tile_iter_end

39: end while

40:joininconsolatait

Should a given CTA’s starting and/or ending iterations not coincide with tile
boundaries (as is expected to be the common case), it must consolidate its
partial results with those of the other CTA(s) also covering that tile. In
this basic implementation, each output tile in C is written by the CTA that
performed that tile’s $k=0$ MAC-loop iteration. Before it can do so, however,
it must accumulate any partial sums shared from other CTAs in temporary global
storage. Notably, _Stream-K’_ s communication, synchronization, and global
storage overheads are independent of problem size, scaling instead with the
number of CTAs.

A secondary benefit of _Stream-K_ is that synchronization-waiting is likely
negligible when the number of output tiles is greater than the number of CTAs.
In this regime, each output tile is covered by at most two CTAs, and the tile-
processing skew ensures that the accumulating CTA will not need its peer
contributions until well after those collaborators have finished producing
them.

Continuing our earlier example, Figure 2(b) illustrates the basic _Stream-K_
execution schedule of the $384\times 384\times 128$ GEMM problem on a
hypothetical four-SM GPU. To fully occupy the GPU, we launch $g=4$ CTAs.
Assuming BLK_M $=128$, BLK_N $=128$, and BLK_K $=4$, each CTA is tasked with a
$128\times 128\times 288$ work volume comprising 72 MAC-loop iterations. This
results in a 100% quantization efficiency, as all four SMs will execute the
same number of MAC instructions.

Additionally, the work volume of a single MAC-loop iteration is 32$\times$
smaller than that of an entire output tile. Consequently, a 32-way _fixed-
split_ decomposition would also provide a 100% quantization efficiency, but at
the expense of an 8$\times$ larger “fixup” overhead. Furthermore, _Stream-K_
is better able to hide the latency of inter-CTA synchronization due to the
temporal skew between writers and readers when sharing partial sums.

_Stream-K_ also generalizes to both _fixed-split_ and _data-parallel_
decompositions. When the grid size $g$ is an even multiple of the number of
output tiles, _Stream-K_ functions exactly as the _fixed-split_ decomposition.
Similarly, when $g$ equals the number of output tiles, _Stream-K_ behaves
identically to the _data-parallel_ decomposition. We take advantage of this
generalization to create an optimized hybridization of the _Stream-K_
decomposition in following section (2).

###  3 Implementation Details

(a) Basic Stream-K

(b) DP + one-tile SK

(c) Two-tile SK + DP

Figure 3: Basic Stream-K vs. hybrid execution schedules for $896\times
384\times 128$ GEMM across a hypothetical four-SM GPU.

The work decomposition we introduced in the last section can be instantiated
in a number of different ways to suit the needs of different hardware
architectures and software library designs. Our implementation targets NVIDIA
GPUs and is designed to be integrated into existing libraries like cuBLAS and
CUTLASS. In this section, we describe how we configure the kernels we launch
and introduce a hybridization scheme that helps ensure users achieve maximum
GEMM performance across the widest possible range of problem shapes.

We also emphasize that these are truly internal implementation details. They
are completely transparent to the user of a BLAS-like library and do not alter
the library’s interface. The only observable impact is the improved
performance characteristics that we analyze in Section 4.

####  1 Kernel Configuration

The tile size chosen for blocking the GEMM computation is, of course, a
critical parameter controlling the performance of the GEMM kernel. For modern
NVIDIA GPUs, appropriate tile sizes are determined by the shape of matrices
supported by the GPU’s Tensor Cores. Based on extensive empirical experience,
we selected the smallest CTA-wide tile size capable of achieving 99% of the
GPU’s peak TFLOP/s for very large GEMM volumes for each supported precision.
For the NVIDIA A100 GPU used in our experiments, these sizes are
64$\times$64$\times$16 for FP64 problems and 128$\times$128$\times$32 for
FP16$\rightarrow$32 problems.

Achieving maximal GEMM performance from _Stream-K_ parallelization also
requires some degree of dynamic problem-specific configuration. Before
launching a kernel we choose a grid size likely to yield the best performance
on the specific problem shape at hand. This is in contrast to ensemble-based
approaches which accommodate diverse problem shapes through the static
generation of many kernel variants based on workload decomposition and
blocking factor.

Our grid size selection heuristic is based on a simple analytical model that
minimizes the cost of reading, writing, and accumulating partial sums while
equally distributing the MAC-loop iterations per CTA. Details of this
analytical model are provided in the section below. Parameters to the model
are trivially chosen with empirical measurements and need only be done once
per target architecture. The resulting parameters can then be compiled
statically into the library. Again, this is in contrast to ensemble-based
approaches that rely on potentially complex heuristics and machine learning
models for kernel selection at run time.

##### Analytical Modeling for Stream-K Configuration

In practice, it is not always advantageous to invoke the _Stream-K_
decomposition with as many CTAs as can be actively resident on the GPU.
Because it is a tile-splitting approach, it incurs fixup costs above and
beyond the simple _data-parallel_ decomposition. Consequently, the fundamental
proposition is one of strong scaling: how much additional parallelism can be
expressed before the extra overhead causes a negative return on investment.
Depending on the problem shape, the optimal splitting could be enough to fill
the entire processor (i.e., $g\leftarrow p$), no splitting at all (i.e.,
$g\leftarrow t$), or somewhere in between.

To predict this inflection point, we present a simple approach to model the
runtime of _Stream-K_ as a function of grid size $g$. In the absence of other
work on the GPU, the runtime of the entire _Stream-K_ schedule will be the
same as that of one of its tile-outputting CTAs, which we formulate as
follows:

$\displaystyle time_{CTA}(g)\leftarrow$
$\displaystyle\mathpzc{a}+\mathpzc{b}(FixupPeers(g)>1)$
$\displaystyle+\mathpzc{c}(ItersPerCta(g))+\mathpzc{d}(FixupPeers(g)-1)$

where:

$\displaystyle ItersPerCta(g)\leftarrow$
$\displaystyle\left\lceil\frac{\lceil\frac{m}{\text{BLK\\_M}}\rceil\times\lceil\frac{n}{\text{BLK\\_N}}\rceil\times\lceil\frac{k}{\text{BLK\\_K}}\rceil}{g}\right\rceil$
$\displaystyle FixupPeers(g)\leftarrow$
$\displaystyle\left\lceil\frac{\left\lceil\frac{k}{\text{BLK\\_K}}\right\rceil}{IterationsPerCta(g)}\right\rceil$

This CTA runtime model comprises four components. The $\mathpzc{a}$ workload
encompasses the one-time, fixed-size costs incurred by each CTA, e.g., the
grid launch latency, the initial compulsory cache misses, the cost of writing
the final output tile to C, etc. The second component $\mathpzc{b}$
incorporates the conditional costs of outputting temporary partial sums for
scenarios where the number of output tiles does not quantize perfectly across
the processor. The third—the per-iteration workload $\mathpzc{c}$—represents
the instruction and stall workload of each MAC-iteration. The final, per-
collaborator workload $\mathpzc{d}$ is the cost of reading and accumulating
the partial sums from another CTA covering the same tile. The set of workload
constants {$\mathpzc{a}$, $\mathpzc{b}$, $\mathpzc{c}$, $\mathpzc{d}$} will be
unique to each combination of blocking factors, matrix data type, and GPU
microarchitecture, and can be determined empirically via microbenchmarks.

(a) GEMM $256\times 3584\times 8192$
56 output tiles, 256 iterations per tile
$g_{best}\leftarrow 108$ CTAs, 132/133 iterations per CTA

(b) GEMM $128\times 128\times 16384$
1 output tile, 512 iterations per tile
$g_{best}\leftarrow 8$ CTAs, 64 iterations per CTA

(c) GEMM $1024\times 1024\times 1024$
64 output tiles, 32 iterations per tile
$g_{best}\leftarrow 64$ CTAs, 32 iterations per CTA

Figure 4: Modeled Stream-K performance on NVIDIA A100 (108 SMs) for BLK_M=128,
BLK_N=128, BLK_K=32

Figure 4 illustrates the behavior of our grid size selection model as
parameterized for fp16-precision GEMM on NVIDIA’s A100 GPU using blocking
factors BLK_M $=128$, BLK_N $=128$, and BLK_K $=32$. Specifically, we
highlight three strong-scaling GEMM scenarios where the number of output tiles
is insufficient to produce a single full wave across the processor’s 108 SM
cores.

The first GEMM shape accumulates through a large-sized $k$-dimension to
produce a short, wide output matrix. In this scenario, the reduction in MAC-
loop time relative to the increasing costs of seam fixup is monotonically
improving. Consequently, the optimal grid size coincides with maximal
parallelism at $g=108$ CTAs.

The second shape accumulates through a medium-sized $k$-dimension to produce a
square matrix with 64 output tiles. In this case, the fixup costs of
$\mathpzc{b}$ and $\mathpzc{d}$ outweigh any reduction in MAC-loop iteration
count, as seen by the global minima “dip” at $g=64$ CTAs.

Figure 5: Strong-scaling comparison of _data-parallel_ and _Stream-K_
execution schedules for $128\times 128\times 384$ GEMM across a hypothetical
four-SM GPU. _Data-parallel_ causes the enormous $k$-dimension to be
sequentially processed within single CTA, whereas _Stream-K_ is able to take
advantage of the parallelism available across the $k$-dimension.

The third shape produces a single output tile after accumulating through an
enormous $k$-dimension, analogous to the execution schedule in Figure 5.
Although the opportunity for strong scaling is quite large, the per-peer cost
of serial reduction is entirely incurred by a single CTA. These accumulation
costs begin to outweigh any further reductions in iteration count for grid
sizes $g>8$.

####  2 Data-parallel Hybridization

The basic _Stream-K_ decomposition can, in certain cases, exhibit tile-
processing skew that leads to potentially adverse effects on cache
performance. When the number of output tiles $t$ is not an even multiple of
the grid size $g$, the starting $k$-offset for the first MAC-loop iteration in
each CTA will be different. Depending on the sizes and shapes of the input
matrices and blocking factors, this skew may preclude these fragments from
seeing reuse across CTAs in the GPU’s cache structure. In Figure 3(a), for
example, the initial $k$-axis fragment offsets for each of the four CTAs will
be $k=0$, $k=32$, $k=64$, and $k=96$, respectively. Furthermore, this
32-element skew between CTAs will persist for the duration of the GEMM
computation.

Tile-processing skew is a direct consequence of _Stream-K_ ’s workload
balancing strategy. However, we can take measures to limit its duration by
applying _Stream-K_ ’s iteration balancing to a smaller, tile-aligned region
of the total iteration domain such that the remaining tiles can be produced in
full, temporally aligned waves.

The simplest hybrid scheme is the “ _data-parallel_ \+ one-tile _Stream-K_ ”
schedule illustrated in Figure 3(b). It applies iteration balancing only among
the tiles otherwise remaining for a final, partially full _data-parallel_
wave. The total number of full waves is $w=\lfloor t/p\rfloor$, where $t$ is
the number of output tiles and $p$ is the number of SM cores in the GPU.
Consequently, each _Stream-K_ CTA receives an even share of iterations that is
less than one tile’s worth. Unfortunately, this strategy has little ability to
hide the synchronization latency for the exchange of partial sums when three
or more CTAs cover the same tile. In these scenarios, the accumulating CTA may
be forced to wait for the contributions of other CTAs to become visible, as
all but the last will be completing their final iterations at roughly the same
time. Furthermore, the basic version of our scheme for aggregating partials is
serialized within a single CTA, and thus will likely cause SM workload
imbalance when the number of contributing CTAs per tile is large.

We address these problems with our “two-tile _Stream-K_ \+ _data-parallel_ ”
hybrid schedule, illustrated in Figure 3(c). It performs one fewer full data-
parallel wave in exchange for each _Stream-K_ CTA receiving more than one
tile’s worth of iterations (but fewer than two). This provides much better
latency hiding when $w\geq 2$, and each accumulating CTA will only need to
receive partials from one other contributing CTA. Otherwise, it behaves
identically to the “ _DP + one tile SK_ ” schedule. This hybrid approach
results in both improved memory access patterns and latency hiding. It also
shows the versatility of the generic _Stream-K_ looping structure to implement
different scheduling policies within the same kernel instance.

###  4 Performance Evaluation

We have implemented our _Stream-K_ decomposition using NVIDIA’s CUTLASS
library of CUDA C++ template abstractions for authoring GEMM-like
computations. CUTLASS provides the optimized equivalent of the CTA-wide
MacLoop() subroutine in Algorithm 8, which performs blocking, tiling, and
software-pipelined data movement that is analogous to the closed-source cuBLAS
and cuDNN implementations. Our evaluation encompasses both (1) double-
precision FP64 GEMM, and (2) mixed-precision FP16$\rightarrow$32 GEMM. For the
latter, the input matrices A and B comprise half-precision FP16 values, yet
the internal accumulation and output matrix C values are single-precision
FP32.

###### Hardware environment.

Our test GPU is the NVIDIA A100, which contains 108 SM cores and can issue
13,824 data-parallel instructions per cycle. It provides 40 MB of L2 cache and
1,555 GB/s of global memory bandwidth. For measurement stability, we lock the
power envelope at 400 W and SM clocks at 1005 MHz ($\sim$71% of their dynamic
peak). This establishes FP64 tensor-core peak throughput of 13.9 TFLOP/s, and
mixed FP16$\rightarrow$32 tensor-core peak throughput of 222.3 TFLOP/s.

Figure 6: The test domain of 32,824 GEMM problem shapes and sizes used for
performance evaluation. $\\{m\\}=\\{128\ldots 8192\\}$, $\\{n\\}=\\{128\ldots
8192\\}$, $\\{k\\}=\\{128\ldots 8192\\}$

###### Dataset.

Our test corpus intends to approximate the enormous breadth and scope of
device-wide GEMM problems that GPU math kernel libraries are designed to
accommodate. As shown in Figure 6, we evaluate 32,768 different problem sizes
and shapes, log-sampled at random within a domain of $m$, $n$, and $k$ matrix
dimensions whose volume spans six orders of magnitude.

###### Methodology.

For both GEMM precisions, we build a single _Stream-K_ kernel that has been
specialized per the guidelines in the Section 3. Furthermore, these kernels
implement our “two-tile _Stream-K_ \+ _data-parallel_ ” hybrid decomposition.
Our evaluation compares each _Stream-K_ kernel with:

  1. 1.

the default _data-parallel_ CUTLASS kernel of the same blocking factor;

  2. 2.

the cuBLAS ensemble for that precision (CUDA 11.6); and

  3. 3.

an idealized oracle that will always select the highest performing _data-
parallel_ CUTLASS blocking factor to execute for a given GEMM instance.

For FP64 problems, this oracle selects among the ensemble of
(32$\times$32$\times$16), (32$\times$64$\times$16), (64$\times$64$\times$16),
(64$\times$128$\times$16), and (128$\times$128$\times$16) blocking factor
specializations. For FP16 $\rightarrow$ 32, it selects among the ensemble of
(64$\times$64$\times$64), (64$\times$128$\times$32),
(128$\times$128$\times$32), and (128$\times$256$\times$32) blocking factor
specializations. These specific specializations are an open-sourced strict
subsets alternative of the corresponding cuBLAS GEMM kernel ensembles.

(a) CUTLASS FP16$\rightarrow$32 _data-parallel_ “roofline”
performance ($\text{blocking factors}=$ 128$\times$128$\times$32).

(b) cuBLAS (ensemble)

(c) Idealized CUTLASS oracle (ensemble)

(d) _Stream-K_ ($\text{blocking factors}=$ 128$\times$128$\times$32)

Figure 7: FP16$\rightarrow$FP32 GEMM “roofline” performance utilization
landscapes on NVIDIA A100 across 32K GEMM problem shapes and sizes.

(a) CUTLASS _data-parallel_
($\text{blocking factors}=$ 64$\times$64$\times$16)

(b) cuBLAS (ensemble)

(c) Idealized CUTLASS oracle (ensemble)

(d) _Stream-K_ ($\text{blocking factors}=$ 64$\times$64$\times$16)

Figure 8: FP64 GEMM “roofline” performance utilization landscapes on NVIDIA
A100 across 32K problem shapes and sizes.

|  vs. CUTLASS $64\times 64\times 16$ |  vs. cuBLAS |  vs. cuBLAS $>150$ ops/B |  vs. CUTLASS oracle
---|---|---|---|---
Average |  $1.23\times$ |  $1.06\times$ |  $1.03\times$ |  $1.05\times$
StdDev |  $0.45$ |  $0.10$ |  $0.03$ |  $0.09$
Min |  $0.77\times$ |  $0.68\times$ |  $0.99\times$ |  $0.70\times$
Max |  $5.63\times$ |  $2.55\times$ |  $1.24\times$ |  $1.64\times$

Table 1: _Stream-K_ FP64 Relative Performance

|  vs. CUTLASS $128\times 128\times 32$ |  vs. cuBLAS |  vs. cuBLAS $>150$ ops/B |  vs. CUTLASS oracle
---|---|---|---|---
Average |  $1.63\times$ |  $1.13\times$ |  $1.15\times$ |  $1.12\times$
StdDev |  $1.46$ |  $0.45$ |  $0.12$ |  $0.37$
Min |  $0.80\times$ |  $0.64\times$ |  $0.98\times$ |  $0.61\times$
Max |  $14.7\times$ |  $6.74\times$ |  $1.85\times$ |  $4.63\times$

Table 2: _Stream-K_ FP16$\rightarrow$32 Relative Performance

The “roofline” plots of Figure 8(a) and Figure 7(a) highlight the spread of
performance produced by the singleton _data-parallel_ CUTLASS kernels. They
plot the percentage of FP64 and FP16$\rightarrow$32 processor utilization as a
function of computational intensity. Ideally, a GEMM implementation’s
performance response would manifest as a narrow band that adheres tightly to
the machine’s bandwidth- and compute-bound performance ceilings. Here, the
_data-parallel_ kernels exhibit a fairly large dynamic range for any given
regime of arithmetic intensity. In contrast, the performance responses from
the equivalent _Stream-K_ kernels in Figure 8(d) and Figure 7(d) are much
tighter. These observations are corroborated by Table 2 and Table 2, which
show the _Stream-K_ kernels outperforming their _data-parallel_ FP64 and
FP16$\rightarrow$32 equivalents by an average of 1.23$\times$ and
1.63$\times$, respectively. For extreme strong-scaling scenarios where
$m\times n$ is small and $k$ is large, our _Stream-K_ kernels demonstrate up
to 5.63$\times$ and 14.7$\times$ speedup, respectively.

The second columns of Table 2 and Table 2 compare our _Stream-K_ performance
with that of cuBLAS. On average, our FP64 and FP16$\rightarrow$32 _Stream-K_
GEMM kernels respectively deliver 6% and 13% greater throughput than their
corresponding cuBLAS ensembles, with peak improvement of 2.55$\times$ and
6.74$\times$. This is a significant improvement over the breadth of 32K GEMM
problem shapes and sizes with 20$\times$ _less_ executable code (a single
kernel for each precision) than NVIDIA’s vendor GEMM library, cuBLAS.

Furthermore, the contrast between the FP64 and FP16$\rightarrow$32 cuBLAS
performance responses (Figure 8(b) and Figure 7(b)) versus those of our
hypothetical CUTLASS oracle ensembles (Figure 8(c) and Figure 7(c)) reveal the
difficulties of designing kernel selection heuristics that deliver
consistently good performance. Despite having access to the same blocking
factor specializations, cuBLAS exhibits substantially wider dynamic ranges
than the idealized _data-parallel_ CUTLASS oracle. The performance spreads of
our _Stream-K_ kernels are narrower still, achieving up to 4.6$\times$ the
idealized oracle performance and underscoring their ability to achieve
utilization levels that are simply not possible from tile-centric work
decompositions.

(a) FP64 _Stream-K_ speedup vs. cuBLAS.

(b) FP16$\rightarrow$32 _Stream-K_ speedup vs. cuBLAS.

Figure 9: _Stream-K_ speedup vs. cuBLAS.

Finally, we observe regimes of small, bandwidth-bound problem shapes where our
largish blocking factors do not compete well against cuBLAS. However, if we
restrict our scope to the domain of compute-bound problems (i.e., FP64
problems having compute intensity $>$ 150 ops/byte and FP16$\rightarrow$32
problems $>$ 400 ops/byte), Figure 9(a) and Figure 9(b) demonstrate that our
singleton _Stream-K_ kernels achieve unilaterally higher performance than the
cuBLAS ensembles. The “noisy” relative performance in the regimes below these
thresholds is not surprising, as _Stream-K_ is attempting to make memory-bound
computations run faster by adding more memory workload. This suggests a few
avenues for future work, namely separate cost-modeling for the memory-bound
regime and/or the bundling of a second _Stream-K_ kernel having smaller tile
size into a two-kernel ensemble.

###  5 Conclusion

We presented _Stream-K_ , a novel parallel workload decomposition technique
for scheduling general matrix multiplication (GEMM) and similar computations
on wide architectures such as GPUs. Unlike other tile-splitting techniques,
the MAC-loop iteration is our unit of workload quantization across processor
cores. This affords excellent strong scaling and workload balancing because
its cost is (1) a constant with respect to the problem shape, and (2)
substantially smaller than that of an entire output tile.

Furthermore, _Stream-K_ produces an $O(p)$ number of splitting seams that are
bound by the number of processor cores. Consequently, the overheads of strong
scaling and workload balancing scale with processor width rather than problem
size. This is a welcome feature for many applications that cannot afford to
allocate large amounts of temporary storage equivalent to the problem output.

Finally, we evaluated our _Stream-K_ approach across a broad spectrum of GEMM
shapes and sizes. We showed that a single blocking configuration of _Stream-K_
can (1) achieve levels of absolute performance that match and/or exceed that
of NVIDIA’s cuBLAS library, even when the latter is operating at near-peak
processor utilization, and (2) do so with much higher levels of performance
consistency. Additionally, _Stream-K_ is an attractive option for library
construction and maintenance, as it presents an opportunity to reduce
distribution sizes by an order of magnitude and removes the need for complex
handcoded heuristics or machine learning models for kernel selection without
compromising performance.

For future works, we identify cache-aware, tile-access patterns such as Morton
Order, an avenue for optimization. We also believe that _Stream-K_
decomposition could provide a similar improved performance response for other
GEMM-like workloads that struggle with the same quantization inefficiencies.

##  Chapter 5 Conclusion

This dissertation identifies load balancing as the key component to achieve
high-performance for both regular and irregular computations on the GPU. In
efforts to address challenges associated with load balancing, we propose an
abstraction that provides a new way to program irregular applications.
Expressed using our simple, yet powerful abstraction, we are able to extend
previously application-specific load-balancing techniques to other workloads
and parallel domains. We also highlight the need to consider resource
imbalance for regular fine-grained problems, and show that our general
_Stream-K_ schedule rids of any complicated heuristics and provides state-of-
the art performance on an already well-optimized problem like General Matrix-
Multiplication. Built on the insights from this work, we propose the following
future research directions.

###  1 Future Research Directions

####  1 Multi-GPU Load Balancing

In Chapter 3 we built an abstraction for load-balancing centered around
_atoms_ , _tiles_ and _sets_ as fundamental building blocks targetted for
single GPU implementations. Some of the underlying load-balancing algorithms
can easily extend to support multi-GPU problems. This is going to be the key
advancement going forward, as sparse-irregular problems in the real world can
be extremely large (for example, graphs with billions of vertices and edges).
These large problems do not fit in a single GPU’s memory, and to accelerate
computation on these problems, the abstraction must be extended to support
multi-GPU and multi-node systems.

Our insights in this dissertation inspire a new research direction, where
specific modular algorithms (such as merge-path 3 SpMV) are well suited for a
multi GPU enviroment and can be easily exposed with the existing building
blocks with an additional levels of hierarchy for _sets_ of _sets_ , where
each _set_ is large enough to fit in a single GPU’s memory. We call workloads
such as merge-path SpMV, isomorphic problems. Isomorphism here implies that
regardless of where the computation occurs, a thread, block, or a separate
GPU, the structure of the algorithm stays the same.

A much more challenging research direction is further extending the multi GPU
programming model to then support amorphous workloads, where different levels
of compute hierarchy will require different kinds of algorithms to complete a
workload. For example, the Radix Sorting algorithm scales nicely for a single
GPU compute hierarchy as implemented in the CUB library [63], but it may not
be the best choice at a multi-GPU boundary as communicating bit-level
information among different GPUs may not be enough to fully utilize the memory
bandwidth. Therefore, we may require different sorting algorithms at different
parts of the compute hierarchy—and similarly, may require different load-
balancing schedules at the compute boundaries. This will make for an
interesting challenge to explore abstractions and APIs for hybrid load-
balancing schemes with heuristics to context switch depending on the amount of
work and the compute boundary one is working at.

####  2 Heuristics using Roofline Model

Figure 1: Speedup of an oracle heuristics SpMV vs. cuSparse’s SpMV across
SuiteSparse. Selects the best load-balancing schedule from the available ones
for a given problem using exhaustive search.

One optimization we highlighted in our results was the use of simple
heuristics to determine which load-balancing algorithm to use. However, our
results show that our simple heuristics are far from an oracle schedule (see
Figure 1) that will pick the right load-balancing algorithm for all sparse
datasets and workloads. One direction for research is to look at the
traditional roofline analytical model and expand it to capture the
irregularity within sparse-linear algebra operations like SpMV or SpMM. Using
this model, we can then build more robust heuristics to select the among the
different supported load-balancing schedules.

####  3 Standardization of Sparse-Matrix Formats

In Chapter 3 we highlight the diversity and importance of sparse-matrix
formats for sparse-linear algebra and graph analytics. Our APIs for the load-
balancing abstraction were designed to accomodate multiple formats for
different workloads at both high- and low-level of the abstraction. Based on
these insights, we propose the standardization of sparse-matrix formats as new
C++ containers for both sequential and parallel programming use cases. The
research challenge here is to find a general and high-performant way to build
and use these containers for various different domains (graphs, sparse-linear
algebra, machine learning, etc.) One potential way to build a general, high-
performant sparse-matrix container will be to use specialized C++ iterators to
define the “sparse-view” of the container. For example, instead of a sparse-
matrix accessed using the three arrays of a CSR (offsets, indices and values),
a sparse-matrix will be accessed using an iterator over the matrix that
provides the tuple: row, column, value (see Listing 2).

Figure 2: Standardized sparse matrix accessed within a generic SpMV. The
underlying sparse format can be changed to a CSC, COO, DIA, but the high-level
loop remains the same.

⬇

// Implements a "format" unaware SpMV.

template<typename type_t>

void spmv(

matrix_t<type_t, std::size_t, csr_view_t>& A,

std::vector<type_t>& x,

std::vector<type_t>& y) {

// Internally a csr type,

// a sparse matrix simply accessed by

// specialized iterator that return row,

// column, value.

for (auto& [i, j, v] : A) {

y[i] += v * x[j];

}

}

####  4 Programming Model for Locality

We believe that along with efficient load-balancing algorithms for irregular
problems, there is also a need for further research in an orthogonal model for
better leveraging locality within these problems. A dynamic tiling approach
for grouping similar data for better caching and a light-weight reordering
approach for reordering sparse-matrices are two possible directions that can
be exploited within this programming model. The challenge is the tradeoffs
that arises when discussing between the impact load-balancing might have on
locality and vice-versa, and this needs to be better quantified/researched.

####  5 _Stream-K_ Abstraction

In our work in Chapter 4, we presented a work-centric approach for parallel
decomposition targeted at GEMM. Impactful future-work will be using a further
generalized version of _Stream-K_ that can distill any regular computation
(such as tensor algebra, convolutions and other machine learning workloads) to
“amount of work” each thread group needs to perform, and a loop that consumes
the assigned work within the thread group until all work is completed (using a
persistent-CTA). The goal of this future work will be to make the performance
benefits attained from using _Stream-K_ more accessible to a wide breadth of
regular problems that suffer from wave quantization and require a
sophisticated heuristics driven approach (e.g. kernels within cuDNN).

####  6 Beyond Today’s CUDA

One thrust we see going forward is looking at new CUDA programming model
additions and support that will better allow for addressing scheduling fine-
grained parallel computations on the GPUs.

###### CUDA and C++

Our load-balancing abstraction and APIs are fully realized with the help of
tools available within modern C++. However, current CUDA (major version 11)
has limited support for the latest C++ standard. We look forward to future
CUDA/C++ with native support for the C++20 Ranges library that extends and
generalizes the standard algorithms and iterator libraries making composable
APIs even more powerful and efficient.

###### Reconfigurable CUDA Grids

A natural way to think about the load-balancing or scheduling problem within
CUDA is through CUDA’s dynamic parallelism, where one can launch CUDA kernels
from within a kernel. Irregular problems can then be load-balanced by first a
kernel looking at the amount of work available for a given problem and then
launching a CUDA kernel with the exact amount of resources required to process
that part of the problem. However, dynamic parallelism within CUDA is not
performant due to the work required by the driver and the APIs to configure a
fresh kernel launch per CUDA thread. One approach to addressing poor dynamic
parallelism performance is by utilizing CUDA’s compute hierarchy. A single
grid can be oversubscribed with the number of threads and/or resources, and
then reconfigured to fit the problem size; effectively “freeing” up resources
for other threads to consume. Since no new work is being spawned, a
reconfiguration on-the-fly can be achieved at a much cheaper cost than
launching new kernels from within a kernel. We propose this as an interesting
direction to improve CUDA for irregular parallel problems.

## References

  * CPP [2017] ISO international standard ISO/IEC 14882:2017(E) - programming language C++.  Technical report, International Organization for Standardization (ISO), 2017.  https://isocpp.org/std/the-standard.
  * Abdelfattah et al. [2016] Ahmad Abdelfattah, David Keyes, and Hatem Ltaief.  KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators.  _ACM Transactions on Mathematical Software_ , 42(3):1–31, June 2016.  ISSN 1557-7295.  10.1145/2818311.
  * Advanced Micro Devices, Inc [2022] Advanced Micro Devices, Inc.  HIP programming guide v5.2.  June 2022.  URL https://docs.amd.com/bundle/HIP-Programming-Guide-v5.2/page/Introduction.html.
  * Aila and Laine [2009] Timo Aila and Samuli Laine.  Understanding the efficiency of ray traversal on GPUs.  In _Proceedings of High Performance Graphics_ , HPG ’09, pages 145–149, August 2009.  10.1145/1572769.1572792.
  * Arifuzzaman et al. [2013] Shaikh Arifuzzaman, Maleq Khan, and Madhav V. Marathe.  PATRIC: a parallel algorithm for counting triangles in massive networks.  In Qi He, Arun Iyengar, Wolfgang Nejdl, Jian Pei, and Rajeev Rastogi, editors, _22nd ACM International Conference on Information and Knowledge Management_ , CIKM ’13, pages 529–538. ACM, 2013.  10.1145/2505515.2505545.
  * Ashari et al. [2014] Arash Ashari, Naser Sedaghati, John Eisenlohr, Srinivasan Parthasarathy, and P. Sadayappan.  Fast sparse matrix-vector multiplication on GPUs for graph applications.  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC ’14, pages 781–792, November 2014.  10.1109/SC.2014.69.
  * Barrachina et al. [2008] Sergio Barrachina, Maribel Castillo, Francisco D. Igual, Rafael Mayo, and Enrique S. Quintana-Orti.  Evaluation and tuning of the level 3 CUBLAS for graphics processors.  In _2008 IEEE International Symposium on Parallel and Distributed Processing_. IEEE, April 2008.  10.1109/ipdps.2008.4536485.
  * Baxter [2013–2016] Sean Baxter.  Moderngpu: Patterns and behaviors for GPU computing.  http://moderngpu.github.io/moderngpu, 2013–2016.
  * Blelloch [1990] Guy E. Blelloch.  Prefix sums and their applications.  Technical Report CMU-CS-90-190, School of Computer Science, Carnegie Mellon University, November 1990.
  * Blelloch et al. [1993] Guy E. Blelloch, Jonathan C. Hardwick, Siddhartha Chatterjee, Jay Sipelstein, and Marco Zagha.  Implementation of a portable nested data-parallel language.  In _Proceedings of the Fourth ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’93, pages 102–111, May 1993.  ISBN 0-89791-589-5.  10.1145/155332.155343.
  * Blumofe and Leiserson [1999] Robert D. Blumofe and Charles E. Leiserson.  Scheduling multithreaded computations by work stealing.  _Journal of the ACM_ , 46(5):720–748, September 1999.  10.1145/324133.324234.
  * Blumofe et al. [1995] Robert D. Blumofe, Christopher F. Joerg, Bradley C. Kuszmaul, Charles E. Leiserson, Keith H. Randall, and Yuli Zhou.  Cilk: An efficient multithreaded runtime system.  In _Proceedings of the Fifth ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’95, pages 207–216, July 1995.  10.1145/209936.209958.
  * Bolz et al. [2003] Jeff Bolz, Ian Farmer, Eitan Grinspun, and Peter Schröder.  Sparse matrix solvers on the GPU: Conjugate gradients and multigrid.  _ACM Transactions on Graphics_ , 22(3):917–924, July 2003.
  * Brahmakshatriya et al. [2021] Ajay Brahmakshatriya, Yunming Zhang, Changwan Hong, Shoaib Kamil, Julian Shun, and Saman Amarasinghe.  Compiling graph applications for GPUs with GraphIt.  In _IEEE/ACM International Symposium on Code Generation and Optimization_ , CGO 2021, pages 248–261, 2021.  10.1109/CGO51591.2021.9370321.
  * Busato and Bombieri [2015] F. Busato and N. Bombieri.  BFS-4K: An efficient implementation of BFS for kepler GPU architectures.  _IEEE Transactions on Parallel and Distributed Systems_ , 26(7):1826–1838, July 2015.  ISSN 1045-9219.  10.1109/TPDS.2014.2330597.
  * Busato and Bombieri [2017] Federico Busato and Nicola Bombieri.  A performance, power, and energy efficiency analysis of load balancing techniques for GPUs.  In _12th IEEE International Symposium on Industrial Embedded Systems_ , SIES 2017, pages 1–8, 2017.  10.1109/SIES.2017.7993387.
  * Busato and Bombieri [2019] Federico Busato and Nicola Bombieri.  Configuring graph traversal applications for GPUs: Analysis of implementation strategies and their correlation with graph characteristics.  In _International Conference on High Performance Computing Simulation_ , HPCS 2019, pages 145–151, 2019.  10.1109/HPCS48598.2019.9188204.
  * Cederman and Tsigas [2008] Daniel Cederman and Philippas Tsigas.  On dynamic load-balancing on graphics processors.  In _Graphics Hardware_ , GH ’08, pages 57–64, June 2008.  10.2312/EGGH/EGGH08/057-064.
  * Cederman and Tsigas [2011] Daniel Cederman and Philippas Tsigas.  Dynamic load balancing using work-stealing.  In Wen-mei W. Hwu, editor, _GPU Computing Gems_ , volume 2, chapter 35, pages 485–499. Morgan Kaufmann, October 2011.  10.1016/B978-0-12-385963-1.00035-6.
  * Chen et al. [2010] Long Chen, Oreste Villa, Sriram Krishnamoorthy, and Guang R Gao.  Dynamic load balancing on single- and multi-GPU systems.  In _2010 IEEE International Symposium on Parallel & Distributed Processing_, IPDPS 2010. IEEE, April 2010.  10.1109/IPDPS.2010.5470413.
  * Chen et al. [2018] Tianqi Chen, Thierry Moreau, Ziheng Jiang, Haichen Shen, Eddie Q. Yan, Leyuan Wang, Yuwei Hu, Luis Ceze, Carlos Guestrin, and Arvind Krishnamurthy.  TVM: End-to-end optimization stack for deep learning.  _CoRR_ , (1802.04799), February 2018.
  * Chen et al. [2022] Yuxin Chen, Benjamin Brock, Serban Porumbescu, Aydın Buluç, Katherine Yelick, and John D. Owens.  Atos: A task-parallel GPU scheduler for graph analytics.  In _Proceedings of the International Conference on Parallel Processing_ , ICPP 2022, August/September 2022.  10.1145/3545008.3545056.
  * Chetlur et al. [2014] Sharan Chetlur, Cliff Woolley, Philippe Vandermersch, Jonathan Cohen, John Tran, Bryan Catanzaro, and Evan Shelhamer.  cuDNN: Efficient primitives for deep learning.  _CoRR_ , (1410.0759), October 2014.
  * Choquette et al. [2018] Jack Choquette, Oliver Giroux, and Denis Foley.  Volta: Performance and programmability.  _IEEE Micro_ , 38(2):42–52, April 2018.  10.1109/MM.2018.022071134.
  * Cui et al. [2010] Xiang Cui, Yifeng Chen, Changyou Zhang, and Hong Mei.  Auto-tuning dense matrix multiplication for GPGPU with cache.  In _Proceedings of the 16th International Conference on Parallel and Distributed Systems_ , ICPADS 2010, pages 237–242, December 2010.  10.1109/icpads.2010.64.
  * Dalton et al. [2015] Steven Dalton, Sean Baxter, Duane Merrill, Luke Olson, and Michael Garland.  Optimizing sparse matrix operations on GPUs using merge path.  In _IEEE International Parallel and Distributed Processing Symposium_ , pages 407–416, 2015.  10.1109/IPDPS.2015.98.
  * Dathathri et al. [2018] Roshan Dathathri, Gurbinder Gill, Loc Hoang, Hoang-Vu Dang, Alex Brooks, Nikoli Dryden, Marc Snir, and Keshav Pingali.  Gluon: a communication-optimizing substrate for distributed heterogeneous graph analytics.  In Jeffrey S. Foster and Dan Grossman, editors, _Proceedings of the 39th ACM SIGPLAN Conference on Programming Language Design and Implementation_ , PLDI 2018, pages 752–768. ACM, June 2018.  10.1145/3192366.3192404.
  * Davidson et al. [2014] Andrew Davidson, Sean Baxter, Michael Garland, and John D. Owens.  Work-efficient parallel GPU methods for single-source shortest paths.  In _Proceedings of the 28th IEEE International Parallel and Distributed Processing Symposium_ , IPDPS 2014, pages 349–359, May 2014.  10.1109/IPDPS.2014.45.  URL http://escholarship.org/uc/item/8qr166v2.
  * Davis and Hu [2011] Timothy A. Davis and Yifan Hu.  The University of Florida sparse matrix collection.  _ACM Trans. Math. Softw._ , 38(1), December 2011.  ISSN 0098-3500.  10.1145/2049662.2049663.  URL https://doi.org/10.1145/2049662.2049663.
  * Esfahani et al. [2022] Mohsen Koohi Esfahani, Peter Kilpatrick, and Hans Vandierendonck.  LOTUS: locality optimizing triangle counting.  In Jaejin Lee, Kunal Agrawal, and Michael F. Spear, editors, _27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’22, pages 219–233. ACM, April 2022.  10.1145/3503221.3508402.
  * Filippone et al. [2017] Salvatore Filippone, Valeria Cardellini, Davide Barbieri, and Alessandro Fanfarillo.  Sparse matrix-vector multiplication on GPGPUs.  _ACM Trans. Math. Softw._ , 43(4), January 2017.  ISSN 0098-3500.  10.1145/3017994.  URL https://doi.org/10.1145/3017994.
  * Fox et al. [2019] James Fox, Alok Tripathy, and Oded Green.  Improving scheduling for irregular applications with logarithmic radix binning.  _2019 IEEE High Performance Extreme Computing Conference_ , pages 1–7, 2019.  10.1109/HPEC.2019.8916333.
  * Gale et al. [2020] Trevor Gale, Matei Zaharia, Cliff Young, and Erich Elsen.  Sparse GPU kernels for deep learning.  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC ’20, November 2020\.  ISBN 9781728199986.  10.1109/SC41405.2020.00021.
  * [34] Google.  Chromium C++ style guide.  URL https://chromium.googlesource.com/chromium/src/+/HEAD/styleguide/c++/c++.md.
  * Green et al. [2012] Oded Green, Robert McColl, and David A. Bader.  GPU merge path: A GPU merging algorithm.  In _Proceedings of the 26th ACM International Conference on Supercomputing_ , ICS ’12, pages 331–340, June 2012.  ISBN 978-1-4503-1316-2.  10.1145/2304576.2304621.
  * Green et al. [2018] Oded Green, James Fox, Alex Watkins, Alok Tripathy, Kasimir Gabert, Euna Kim, Xiaojing An, Kumar Aatish, and David A. Bader.  Logarithmic radix binning and vectorized triangle counting.  In _IEEE High Performance Extreme Computing Conference_ , HPEC 2018, pages 1–7, 2018.  10.1109/HPEC.2018.8547581.
  * Gremse et al. [2015] Felix Gremse, Andreas Höfter, Lars Ole Schwen, Fabian Kiessling, and Uwe Naumann.  GPU-accelerated sparse matrix-matrix multiplication by iterative row merging.  _SIAM Journal on Scientific Computing_ , 37(1):C54–C71, 2015.  10.1137/130948811.
  * Gustavson [1978] Fred G. Gustavson.  Two fast algorithms for sparse matrices: Multiplication and permuted transposition.  _ACM Transactions on Mathematical Software_ , 4(3):250–269, September 1978.  10.1145/355791.355796.
  * Hagedorn et al. [2020] Bastian Hagedorn, Archibald Samuel Elliott, Henrik Barthels, Rastislav Bodik, and Vinod Grover.  Fireiron: A data-movement-aware scheduling language for GPUs.  In _Proceedings of the ACM International Conference on Parallel Architectures and Compilation Techniques_ , PACT ’20, September 2020.  10.1145/3410463.3414632.
  * Harris and Perelygin [2017] Mark Harris and Kyrylo Perelygin.  Cooperative Groups: Flexible CUDA thread programming, 2017.  URL https://developer.nvidia.com/blog/cooperative-groups/.
  * Harris et al. [2007] Mark Harris, Shubhabrata Sengupta, and John D. Owens.  Parallel prefix sum (scan) with CUDA.  In Hubert Nguyen, editor, _GPU Gems 3_ , chapter 39, pages 851–876. Addison Wesley, August 2007.  URL http://www.idav.ucdavis.edu/publications/print_pub?pub_id=916.
  * Hong et al. [2017] Changwan Hong, Aravind Sukumaran-Rajam, Jinsung Kim, and P. Sadayappan.  Multigraph: Efficient graph processing on gpus.  In _26th International Conference on Parallel Architectures and Compilation Techniques_ , PACT 2017, pages 27–40. IEEE Computer Society, September 2017.  10.1109/PACT.2017.48.
  * Hu et al. [2021] Lin Hu, Lei Zou, and Yu Liu.  Accelerating triangle counting on GPU.  In Guoliang Li, Zhanhuai Li, Stratos Idreos, and Divesh Srivastava, editors, _International Conference on Management of Data_ , SIGMOD ’21, pages 736–748. ACM, June 2021.  10.1145/3448016.3452815.
  * Hwang [1992] Kai Hwang.  _Advanced Computer Architecture: Parallelism, Scalability, Programmability_.  McGraw-Hill Higher Education, 1st edition, 1992.  ISBN 0070316228.
  * Jamshidi et al. [2020] Kasra Jamshidi, Rakesh Mahadasa, and Keval Vora.  Peregrine: a pattern-aware graph mining system.  In Angelos Bilas, Kostas Magoutis, Evangelos P. Markatos, Dejan Kostic, and Margo I. Seltzer, editors, _Fifteenth EuroSys Conference 2020_ , EuroSys ’20, pages 13:1–13:16. ACM, April 2020.  10.1145/3342195.3387548.
  * Jatala et al. [2020a] Vishwesh Jatala, Roshan Dathathri, Gurbinder Gill, Loc Hoang, V. Krishna Nandivada, and Keshav Pingali.  A study of graph analytics for massive datasets on distributed multi-GPUs.  In _IEEE International Parallel and Distributed Processing Symposium_ , IPDPS 2020, pages 84–94. IEEE, May 2020a.  10.1109/IPDPS47924.2020.00019.
  * Jatala et al. [2020b] Vishwesh Jatala, Loc Hoang, Roshan Dathathri, Gurbinder Gill, V. Krishna Nandivada, and Keshav Pingali.  An adaptive load balancer for graph analytics on GPUs.  _CoRR_ , abs/1911.09135., 2020b.  URL https://arxiv.org/pdf/1911.09135.
  * Jiang and Snir [2005] Changhao Jiang and Marc Snir.  Automatic tuning matrix multiplication performance on graphics hardware.  In _Proceedings of the 14th International Conference on Parallel Architectures and Compilation Techniques_ , PACT ’05, pages 185–194, 2005.  10.1109/pact.2005.10.
  * Kerr et al. [2017] Andrew Kerr, Duane Merrill, Julien Demouth, and John Tran.  CUTLASS: Fast linear algebra in CUDA C++, 2017.  URL https://devblogs.nvidia.com/cutlass-linear-algebra-cuda/.
  * Khorasani et al. [2015] Farzad Khorasani, Rajiv Gupta, and Laxmi N. Bhuyan.  Scalable SIMD-efficient graph processing on GPUs.  In _2015 International Conference on Parallel Architectures and Compilation_ , PACT 2015, pages 39–50. IEEE Computer Society, 2015.  10.1109/PACT.2015.15.
  * Kurzak et al. [2012] Jakub Kurzak, Stanimire Tomov, and Jack Dongarra.  Autotuning GEMM kernels for the fermi GPU.  _IEEE Transactions on Parallel and Distributed Systems_ , 23(11):2045–2057, November 2012.  10.1109/tpds.2011.311.
  * Lam et al. [1991] Monica D. Lam, Edward E. Rothberg, and Michael E. Wolf.  The cache performance and optimizations of blocked algorithms.  In _Proceedings of the Fourth International Conference on Architectural Support for Programming Languages and Operating Systems_ , ASPLOS IV, pages 63–74, 1991.  10.1145/106972.106981.
  * Larsen and McAllister [2001] E. Scott Larsen and David McAllister.  Fast matrix multiplies using graphics hardware.  In _Proceedings of the 2001 ACM/IEEE Conference on Supercomputing_ , SC ’01, pages 55:1–55:6, 2001.  ISBN 158113293X.  10.1145/582034.582089.
  * Lee et al. [2020] Jeongmyung Lee, Seokwon Kang, Yongseung Yu, Yong-Yeon Jo, Sang-Wook Kim, and Yongjun Park.  Optimization of GPU-based sparse matrix multiplication for large sparse networks.  In _IEEE 36th International Conference on Data Engineering_ , ICDE 2020, pages 925–936, 2020.  10.1109/ICDE48307.2020.00085.
  * Li et al. [2009] Yinan Li, Jack Dongarra, and Stanimire Tomov.  A note on auto-tuning GEMM for GPUs.  In _International Conference on Computational Science_ , ICCS 2009, pages 884–892, 2009.  ISBN 978-3-642-01970-8.  10.1007/978-3-642-01970-8_89.
  * Liu et al. [2017] Bangtian Liu, Chengyao Wen, Anand D. Sarwate, and Maryam Mehri Dehnavi.  A unified optimization approach for sparse tensor operations on GPUs.  _IEEE International Conference on Cluster Computing_ , pages 47–57, September 2017.  10.1109/CLUSTER.2017.75.
  * Liu et al. [2018] Junhong Liu, Xin He, Weifeng Liu, and Guangming Tan.  Register-based implementation of the sparse general matrix-matrix multiplication on GPUs.  In Andreas Krall and Thomas R. Gross, editors, _Proceedings of the 23rd ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’18, pages 407–408. ACM, 2018.  10.1145/3178487.3178529.
  * Liu et al. [2013] Xing Liu, Mikhail Smelyanskiy, Edmond Chow, and Pradeep Dubey.  Efficient sparse matrix-vector multiplication on x86-based many-core processors.  In _Proceedings of the 27th International ACM Conference on International Conference on Supercomputing_ , ICS ’13, pages 273–282, New York, NY, USA, 2013. Association for Computing Machinery.  ISBN 9781450321303.  10.1145/2464996.2465013.  URL https://doi.org/10.1145/2464996.2465013.
  * Liu et al. [2022] Zhuang Liu, Hanzi Mao, Chao-Yuan Wu, Christoph Feichtenhofer, Trevor Darrell, and Saining Xie.  A ConvNet for the 2020s.  _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2022.
  * Manca et al. [2016] Emanuele Manca, Andrea Manconi, Alessandro Orro, Giuliano Armano, and Luciano Milanesi.  CUDA-quicksort: an improved GPU-based implementation of quicksort.  _Concurrency and Computation Practice and Experience_ , 28(1):21–43, January 2016.  10.1002/cpe.3611.
  * Mattson et al. [2020] Peter Mattson, Christine Cheng, Gregory Diamos, Cody Coleman, Paulius Micikevicius, David Patterson, Hanlin Tang, Gu-Yeon Wei, Peter Bailis, Victor Bittorf, David Brooks, Dehao Chen, Debo Dutta, Udit Gupta, Kim Hazelwood, Andy Hock, Xinyuan Huang, Daniel Kang, David Kanter, Naveen Kumar, Jeffery Liao, Deepak Narayanan, Tayo Oguntebi, Gennady Pekhimenko, Lillian Pentecost, Vijay Janapa Reddi, Taylor Robie, Tom St John, Carole-Jean Wu, Lingjie Xu, Cliff Young, and Matei Zaharia.  MLPerf training benchmark.  In I. Dhillon, D. Papailiopoulos, and V. Sze, editors, _Proceedings of Machine Learning and Systems_ , volume 2, pages 336–349, 2020\.
  * Mehrabi et al. [2021] Atefeh Mehrabi, Donghyuk Lee, Niladrish Chatterjee, Daniel J. Sorin, Benjamin C. Lee, and Mike O’Connor.  Learning sparse matrix row permutations for efficient spmm on GPU architectures.  In _International Symposium on Performance Analysis of Systems and Software_ , ISPASS 2021, pages 48–58. IEEE, March 2021.  10.1109/ISPASS51385.2021.00016.
  * Merrill [2013–2016] Duane Merrill.  Cub: Flexible library of cooperative threadblock primitives and other utilities for CUDA kernel programming.  https://nvlabs.github.io/cub/, 2013–2016.
  * Merrill and Garland [2016] Duane Merrill and Michael Garland.  Merge-based parallel sparse matrix-vector multiplication.  In _International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC ’16, pages 678–689, November 2016.  ISBN 9781467388153.  10.1109/SC.2016.57.
  * Merrill et al. [2012] Duane Merrill, Michael Garland, and Andrew Grimshaw.  Scalable GPU graph traversal.  In _Proceedings of the 17th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’12, pages 117–128, February 2012\.  10.1145/2145816.2145832.
  * Nath et al. [2010] Rajib Nath, Stanimire Tomov, and Jack Dongarra.  An improved Magma Gemm for Fermi graphics processing units.  _The International Journal of High Performance Computing Applications_ , 24(4):511–515, November 2010.  10.1177/1094342010385729.
  * Nisa et al. [2019] Israt Nisa, Jiajia Li, Aravind Sukumaran-Rajam, Richard Vuduc, and P. Sadayappan.  Load-balanced sparse MTTKRP on GPUs.  In _IEEE International Parallel and Distributed Processing Symposium_ , IPDPS 2019, pages 123–133, May 2019.  10.1109/IPDPS.2019.00023.
  * Niu et al. [2022] Yuyao Niu, Zhengyang Lu, Haonan Ji, Shuhui Song, Zhou Jin, and Weifeng Liu.  TileSpGEMM: a tiled algorithm for parallel sparse general matrix-matrix multiplication on GPUs.  In Jaejin Lee, Kunal Agrawal, and Michael F. Spear, editors, _27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , pages 90–106. ACM, 2022.  10.1145/3503221.3508431.
  * NVIDIA Corporation [2016] NVIDIA Corporation.  NVIDIA CUDA C programming guide.  PG-02829-001_v8.0, September 2016.
  * NVIDIA Corporation [2019] NVIDIA Corporation.  Tips for optimizing GPU performance using tensor cores, 2019.  https://developer.nvidia.com/blog/optimizing-gpu-performance-tensor-cores/.
  * NVIDIA Corporation [2020a] NVIDIA Corporation.  NVIDIA A100 tensor core GPU architecture, 2020a.  https://images.nvidia.com/aem-dam/en-zz/Solutions/data-center/nvidia-ampere-architecture-whitepaper.pdf.
  * NVIDIA Corporation [2020b] NVIDIA Corporation.  CUDA cuBLAS library (v9.2), 2020b.  http://developer.nvidia.com/cublas.
  * NVIDIA Corporation [2021] NVIDIA Corporation.  Deep learning perfomance documentation, 2021.  https://docs.nvidia.com/deeplearning/performance/dl-performance-matrix-multiplication/index.html.
  * Pai and Pingali [2016] Sreepathi Pai and Keshav Pingali.  A compiler for throughput optimization of graph algorithms on GPUs.  In _Proceedings of the 2016 ACM SIGPLAN International Conference on Object-Oriented Programming, Systems, Languages, and Applications_ , OOPSLA 2016, pages 1–19, November 2016.  ISBN 978-1-4503-4444-9.  10.1145/2983990.2984015.
  * Parger et al. [2020] Mathias Parger, Martin Winter, Daniel Mlakar, and Markus Steinberger.  spECK: accelerating GPU sparse matrix-matrix multiplication through lightweight analysis.  In _ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , pages 362–375, February 2020.  10.1145/3332466.3374521.
  * Ragan-Kelley et al. [2013] Jonathan Ragan-Kelley, Connelly Barnes, Andrew Adams, Sylvain Paris, Frédo Durand, and Saman Amarasinghe.  Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines.  In _Proceedings of the 34th ACM SIGPLAN Conference on Programming Language Design and Implementation_ , PLDI ’13, pages 519–530, 2013\.  ISBN 9781450320146.  10.1145/2491956.2462176.
  * Steinberger et al. [2014] Markus Steinberger, Michael Kenzel, Pedro Boechat, Bernhard Kerbl, Mark Dokter, and Dieter Schmalstieg.  Whippletree: Task-based scheduling of dynamic workloads on the GPU.  _ACM Transactions on Graphics_ , 33(6):228:1–228:11, November 2014.  ISSN 0730-0301.  10.1145/2661229.2661250.
  * Steinberger et al. [2017] Markus Steinberger, Rhaleb Zayer, and Hans-Peter Seidel.  Globally homogeneous, locally adaptive sparse matrix-vector multiplication on the GPU.  In _Proceedings of the International Conference on Supercomputing_ , ICS 2017, pages 13:1–13:11. ACM, 2017.  10.1145/3079079.3079086.
  * Tan et al. [2011] Guangming Tan, Linchuan Li, Sean Treichle, Everett Phillips, Yungang Bao, and Ninghui Sun.  Fast implementation of DGEMM on fermi GPU.  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC11, pages 35:1–35:11, Seattle, Washington, 2011.  10.1145/2063384.2063431.
  * Tillet and Cox [2017] Philippe Tillet and David Cox.  Input-aware auto-tuning of compute-bound HPC kernels.  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC17, pages 43:1–43:12, November 2017.  ISBN 9781450351140.  10.1145/3126908.3126939.
  * Tillet et al. [2019] Philippe Tillet, H. T. Kung, and David Cox.  Triton: An intermediate language and compiler for tiled neural network computations.  In _Proceedings of the 3rd ACM SIGPLAN International Workshop on Machine Learning and Programming Languages_ , MAPL 2019, pages 10–19, 2019.  ISBN 9781450367196.  10.1145/3315508.3329973.
  * Tzeng et al. [2010] Stanley Tzeng, Anjul Patney, and John D. Owens.  Task management for irregular-parallel workloads on the GPU.  In _Proceedings of High Performance Graphics_ , HPG ’10, pages 29–37, June 2010.  10.2312/EGGH/HPG10/029-037.  URL https://escholarship.org/uc/item/9r15d4zk.
  * Tzeng et al. [2012] Stanley Tzeng, Brandon Lloyd, and John D. Owens.  A GPU task-parallel model with dependency resolution.  _IEEE Computer_ , 45(8):34–41, August 2012.  10.1109/MC.2012.255.  URL https://escholarship.org/uc/item/4956q122.
  * Wang et al. [2017] Yangzihao Wang, Yuechao Pan, Andrew Davidson, Yuduo Wu, Carl Yang, Leyuan Wang, Muhammad Osama, Chenshan Yuan, Weitang Liu, Andy T. Riffel, and John D. Owens.  Gunrock: GPU graph analytics.  _ACM Transactions on Parallel Computing_ , 4(1):3:1–3:49, August 2017.  10.1145/3108140.  URL http://escholarship.org/uc/item/9gj6r1dj.
  * Winter et al. [2019] Martin Winter, Daniel Mlakar, Rhaleb Zayer, Hans-Peter Seidel, and Markus Steinberger.  Adaptive sparse matrix-matrix multiplication on the GPU.  In Jeffrey K. Hollingsworth and Idit Keidar, editors, _Proceedings of the 24th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’19, pages 68–81. ACM, February 2019\.  10.1145/3293883.3295701.
  * Wu et al. [2015] Yuduo Wu, Yangzihao Wang, Yuechao Pan, Carl Yang, and John D. Owens.  Performance characterization of high-level programming models for GPU graph analytics.  In _IEEE International Symposium on Workload Characterization_ , IISWC-2015, pages 66–75, October 2015.  10.1109/IISWC.2015.13.  URL http://escholarship.org/uc/item/2t69m5ht.
  * Yan et al. [2014] Shengen Yan, Chao Li, Yunquan Zhang, and Huiyang Zhou.  yaSpMV: yet another SpMV framework on GPUs.  In José E. Moreira and James R. Larus, editors, _ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’14, pages 107–118. ACM, 2014.  10.1145/2555243.2555255.
  * Yang et al. [2018] Carl Yang, Aydın Buluç, and John D. Owens.  Design principles for sparse matrix multiplication on the GPU.  In Marco Aldinucci, Luca Padovani, and Massimo Torquati, editors, _Euro-Par 2018: Proceedings of the 24th International European Conference on Parallel and Distributed Computing_ , pages 672–687, August 2018\.  10.1007/978-3-319-96983-1_48.  URL https://escholarship.org/uc/item/5h35w3b7.
  * Yang et al. [2022] Carl Yang, Aydın Buluç, and John D. Owens.  GraphBLAST: A high-performance linear algebra-based graph framework on the GPU.  _ACM Transactions on Mathematical Software_ , 48(1):1:1–1:51, February 2022.  10.1145/3466795.  URL https://escholarship.org/uc/item/292901ks.
  * Zhang et al. [2022] Lingqi Zhang, Mohamed Wahib, Peng Chen, Jintao Meng, Xiao Wang, and Satoshi Matsuoka.  Persistent kernels for iterative memory-bound gpu applications.  _CoRR_ , April 2022.
  * Zhang et al. [2015] Peter Zhang, Eric Holk, John Matty, Samantha Misurda, Marcin Zalewski, Jonathan Chu, Scott McMillan, and Andrew Lumsdaine.  Dynamic parallelism for simple and efficient GPU graph algorithms.  In Antonino Tumeo, John Feo, and Oreste Villa, editors, _Proceedings of the 5th Workshop on Irregular Applications - Architectures and Algorithms_ , IA3 2015, pages 11:1–11:4. ACM, 9 2015.  10.1145/2833179.2833189.
  * Zhang et al. [2014] Tao Zhang, Wei Shu, and Min-You Wu.  CUIRRE: An open-source library for load balancing and characterizing irregular applications on GPUs.  _Journal of Parallel and Distributed Computing_ , 74(10):2951–2966, October 2014.  10.1016/j.jpdc.2014.07.004.