CUTLASS: Fast Linear Algebra in CUDA C++

This article delves into NVIDIA’s open-source CUTLASS library, which is a collection of templates for implementing high-performance General Matrix Multiplication (GEMM) in CUDA C++. The article thoroughly analyzes the hierarchical structure of efficient GEMM computation on modern GPUs, from thread blocks, warps to thread-level blocking strategies, and covers key techniques such as mixed-precision computation, Tensor Core utilization, and software pipelining. Additionally, it demonstrates how to easily integrate custom element-wise operations (such as bias addition and ReLU activation functions) using CUTLASS, making it a powerful tool for compute-intensive applications like deep learning.

1. Introduction to CUTLASS
2. Efficient Matrix Multiplication on GPUs
3. Hierarchical GEMM Structure

• Thread Block Blocking
• Warp Blocking
• Thread Blocking

4. WMMA GEMM: Utilizing Tensor Core

5. Complete GEMM Process

6. Software Pipelining: Hiding Latency

7. Detailed Components of the CUTLASS Library

• CUTLASS GEMM Device Functions
• CUTLASS GEMM Strategies

8. Fusing Element-wise Operations: An Example with SGEMM

9. Performance of Tesla V100 (Volta)

10. Getting Started with CUTLASS

Introduction to CUTLASS

Matrix multiplication is a core computation in many scientific applications, especially in the field of deep learning. Many operations in modern deep neural networks are either defined as matrix multiplications or can be transformed into matrix multiplications.

For example, the NVIDIA cuDNN library uses various forms of matrix multiplication to implement convolutions in neural networks, such as the classic direct convolution formula, which is the matrix product between images and filter datasets. Convolution computations based on Fast Fourier Transform (FFT)^[1] or Winograd methods^[2] also have matrix multiplication as their core routine.

When building cuDNN, we started from the high-performance General Matrix Multiplication (GEMM) implementation in the cuBLAS library, supplementing and customizing it for efficient convolution computation. Today, adjusting these GEMM strategies and algorithms is crucial for providing optimal performance for various problems and applications in deep learning.

With CUTLASS, we aim to provide developers with the necessary techniques and structures to develop new algorithms in CUDA C++ using high-performance GEMM components. Efficient and flexible applications of dense linear algebra are vital in deep learning and the broader GPU computing ecosystem.

CUTLASS (CUDA Templates for Linear Algebra Subroutines) is a collection of CUDA C++ templates and abstractions for implementing high-performance GEMM computations at all levels and scales within CUDA kernels. Unlike other templated GPU libraries for dense linear algebra (such as the MAGMA library^[3]), CUTLASS aims to decompose the “active components” of GEMM into fundamental components abstracted by C++ template classes, allowing programmers to easily customize and specialize them in their own CUDA kernels.

Our CUTLASS primitives include extensive support for mixed-precision computations, providing specialized data movement and multiply-accumulate abstractions for handling 8-bit integers, half-precision floating-point (FP16), single-precision floating-point (FP32), and double-precision floating-point (FP64) types. One of the most exciting features of CUTLASS is the implementation of matrix multiplication running on the new Volta architecture Tensor Core using the WMMA API.

The Tensor Core of the Tesla V100 is a programmable matrix multiply-accumulate unit that can efficiently deliver performance of up to 125 Tensor TFLOP/s.

CUTLASS: Fast Linear Algebra in CUDA C++ — Matrix-matrix multiplication hierarchy used in CUTLASS

Efficient Matrix Multiplication on GPUs

The GEMM computation C = alpha * A * B + beta * C, where A, B, and C are matrices. A is a M×K matrix, B is a K×N matrix, and C is a M×N matrix. For simplicity, we assume scalars alpha=beta=1 in the following example. Later, we will show how to implement custom element-wise operations that support arbitrary scaling functions using CUTLASS.

The simplest implementation consists of three nested loops:

for (int i = 0; i < M; ++i)
    for (int j = 0; j < N; ++j)
       for (int k = 0; k < K; ++k) 
            C[i][j] += A[i][k] * B[k][j];

The element at position <code>(i, j) in C is the dot product of the i row of A and the j column of B with K elements. Ideally, performance should be limited by the arithmetic throughput of the processor. Indeed, for large square matrices where M=N=K, the number of mathematical operations in the matrix product is O(N³), while the amount of data required is O(N²), resulting in a computational intensity of the order of N. However, achieving this theoretical computational intensity requires reusing each element O(N) times. Unfortunately, the above “dot product” algorithm relies on keeping large data in fast on-chip caches, which leads to thrashing as M, N, and K grow.

A better form is to restructure the loop by making the K dimension the outermost loop. This form of computation loads a column of A and a row of B at once, computes their outer product, and accumulates the result of this outer product into the matrix C. After that, this column of A and this row of B will no longer be used.

for (int k = 0; k < K; ++k)     // K dimension is now the outermost loop
    for (int i = 0; i < M; ++i)
        for (int j = 0; j < N; ++j)
            C[i][j] += A[i][k] * B[k][j];

One issue with this method is that it requires all M×N elements of C to be active to store the result of each multiply-accumulate instruction, ideally, the write speed of the memory should match the compute speed of the multiply-accumulate instructions. We can reduce the data size of C by partitioning it into tiles of size Mtile×Ntile, ensuring that these tiles fit into on-chip memory. We then apply the “outer product” formula to each tile. This leads to the following nested loops.

for (int m = 0; m < M; m += Mtile)                // Traverse M dimension
    for (int n = 0; n < N; n += Ntile)            // Traverse N dimension
        for (int k = 0; k < K; ++k)
            for (int i = 0; i < Mtile; ++i)       // Compute a tile
                for (int j = 0; j < Ntile; ++j) {
                    int row = m + i;
                    int col = n + j;
                    C[row][col] += A[row][k] * B[k][col];
                }

For each tile of C, the tiles of A and B are fetched only once, achieving a computational intensity of O(N). The size of each C tile can be chosen to match the capacity of the L1 cache or registers of the target processor, and the outer loops of the nesting can be easily parallelized. This is a significant improvement!

Further restructuring provides more opportunities to exploit locality and parallelism. Instead of specifically accumulating the outer product of vectors, we accumulate the product of matrices by traversing the K dimension in tiles. We typically refer to this concept as accumulating matrix products.

Hierarchical GEMM Structure

CUTLASS efficiently implements GEMM for GPUs by decomposing the computation into a hierarchy of thread block blocking, warp blocking, and thread blocking, applying the strategy of accumulating matrix products. This hierarchy closely reflects the NVIDIA CUDA programming model, as shown in Figure 1. Here, you can see the movement of data from global memory to shared memory (matrix to thread block blocking), from shared memory to register files (thread block blocking to warp blocking), and from register files to CUDA cores for computation (warp blocking to thread blocking).

Thread Block Blocking

Each thread block computes its corresponding output GEMM portion by iteratively loading data blocks from the input matrices and calculating the accumulated matrix product (C += A * B). Figure 2 shows the computation performed by a single thread block, highlighting the data blocks used in one iteration of its main loop.

The CUDA thread block blocking structure is further divided into warps (groups of threads executing together in a SIMT manner). Warps provide a useful organization for GEMM computations and are an explicit part of the WMMA API that we will discuss later.

Figure 3 shows a detailed view of a block-level matrix multiplication structure. The tiles of A and B are loaded from global memory and stored in shared memory accessible to all warps. As shown in Figure 3, the output tiles of the thread block are spatially partitioned among the warps. We refer to the storage of this output tile as an accumulator, as it stores the results of the accumulated matrix products. Each accumulator is updated once per mathematical operation, so it needs to be located in the fastest memory in the SM: register files.

The parameter BlockItems{X,Y,K} is a compile-time constant specified by the programmer to adjust the aspect ratio of the GEMM computation for the target processor and specific GEMM configuration (e.g., M, N, K, data types, etc.). In the figure, we show a thread block with eight warps and 256 threads, which is a typical configuration for large SGEMM (FP32 GEMM) blocking implemented in CUTLASS.

Warp Blocking

Once the data is stored in shared memory, each warp computes a series of accumulated matrix products by traversing the K dimension of the thread block tiles, loading sub-matrices (or fragments) from shared memory and calculating the accumulated outer product. Figure 4 shows a detailed view. The size of the K dimension of the fragments is typically very small to maximize the computational intensity relative to the amount of data loaded from shared memory, thus avoiding shared memory bandwidth becoming a bottleneck.

We note that the warp-centric organization effectively implements efficient GEMM kernels but does not rely on implicit warp synchronization for synchronization. The CUTLASS GEMM kernels synchronize well by calling __syncthreads() at appropriate times.

Thread Blocking

The CUDA programming model is defined in terms of thread blocks and individual threads. Therefore, the warp structure is mapped to operations performed by individual threads. Threads cannot access each other’s registers, so we must choose an organization that allows values stored in registers to be reused across multiple mathematical instructions executed by the same thread. This leads to a 2D blocking structure within threads, as shown in the detailed view in Figure 5. Each thread issues a series of independent mathematical instructions to the CUDA core and computes an accumulated outer product.

In Figure 5, the upper left quadrant of the warp is shaded in gray. The 32 cells correspond to the 32 threads within a warp. This arrangement causes multiple threads in the same row or column to fetch the same elements of the fragments of A and B respectively. To maximize computational intensity, this basic structure can be replicated to form a complete warp-level accumulator block, resulting in a 8×8 overall thread block, computed from the outer product of 8×1 and 1×8 fragments. This is illustrated by the four green accumulator blocks in the figure.

WMMA GEMM: Utilizing Tensor Core

The warp blocking structure can be implemented using the CUDA Warp Matrix Multiply API (WMMA) introduced in CUDA 9, targeting the Tensor Core of the Volta V100 GPU. For more details on the WMMA API, refer to the article “Programming Tensor Cores with CUDA 9^[6]“.

Each Tensor Core provides a 4x4x4 matrix processing array that performs the operation D = A * B + C, where A, B, C, and D are 4×4 matrices, as shown in Figure 6. The matrix multiplication inputs A and B are FP16 matrices, while the accumulator matrices C and D can be either FP16 or FP32 matrices.

In practice, the WMMA API is an alternative to the thread blocking structure of warp-level matrix multiply operations described in the previous section. The WMMA API does not decompose the warp blocking structure into scalar and vector elements owned by individual threads, but rather provides an abstraction for matrix fragment loading/storing and multiply-accumulate mathematical operations for warp cooperation.

Figure 7 shows the warp blocking structure designed for the CUDA WMMA API. Calling wmma::load_matrix_sync loads fragments of A and B into instances of the nvcuda::wmma::fragment<> template, and the elements of the warp blocking accumulators are constructed as an array of nvcuda::wmma::fragment<accumulator> objects. These fragments store a 2D matrix distributed among the warp threads. Finally, calling nvcuda::wmma::mma_sync() for each accumulator fragment (as well as the corresponding fragments from A and B) performs the matrix multiply-accumulate operation over the warp using Tensor Core.

CUTLASS implements a WMMA API-based GEMM in the file block_task_wmma.h. The dimensions of the warp blocking must be multiples of the matrix multiplication shapes defined by the <code>nvcuda::wmma template for the target CUDA compute capability. In CUDA 9.0, the basic WMMA size is <code>16×16×16.

Complete GEMM Process

The complete GEMM structure can be represented as nested loops executed by the threads of a thread block, as shown in the code listing below. Except for the outermost “main” loop, all loops have constant iteration counts that can be fully unrolled by the compiler. For brevity, address and index calculations are omitted here but are detailed in the CUTLASS source code.

// Device function for computing the accumulated matrix product of a thread block
__device__ void block_matrix_product(int K_dim) {
    
    // Fragments for storing data fetched from SMEM
    value_t frag_a[ThreadItemsY];
    value_t frag_b[ThreadItemsX];

    // Accumulator storage
    accum_t accumulator[ThreadItemsX][ThreadItemsY];

    // GEMM main loop - iterating over the entire K dimension - not unrolled
    for (int kblock = 0; kblock < K_dim; kblock += BlockItemsK) {

        // Load A and B tiles from global memory and store to SMEM
        // (not shown for brevity - see CUTLASS source code)
        ...

        __syncthreads();
        
        // Warp blocking structure - iterating over thread block tiles
        #pragma unroll
        for (int warp_k = 0; warp_k < BlockItemsK; warp_k += WarpItemsK) {

            // Fetch frag_a and frag_b corresponding to k index from SMEM
            // (not shown for brevity - see CUTLASS source code)
            ...

            // Thread blocking structure - accumulate an outer product
            #pragma unroll
            for (int thread_x = 0; thread_x < ThreadItemsX; ++thread_x) {
                #pragma unroll
                for (int thread_y=0; thread_y < ThreadItemsY; ++thread_y) {
                    accumulator[thread_x][thread_y] += frag_a[y]*frag_b[x];
                }
            }
        }

        __syncthreads();
    }   
}

WarpItemsK refers to the size of the dot product for the target mathematical operation. For SGEMM (<code>FP32 GEMM), DGEMM (<code>FP64), and HGEMM (<code>FP16), the dot product length for scalar multiply-accumulate instructions is 1. For IGEMM (8-bit integer GEMM), CUTLASS targets the four-element integer dot product instruction (<code>IDP4A), setting <code>WarpItemsK=4. For WMMA-based GEMM, we choose the <code>wmma::fragment<> template's <code>K dimension. Currently, this is defined as <code>WarpItemsK=16.

Software Pipelining: Hiding Latency

Blocked matrix multiplication extensively uses register files to hold fragments and accumulator tiles, as well as a significant amount of shared memory allocation. The relatively high on-chip storage requirements limit occupancy, which is the maximum number of thread blocks that can be simultaneously executed on a single SM. Therefore, the number of warps and thread blocks that GEMM implementations can accommodate in each SM is far less than the typical number of GPU compute workloads.We use software pipelining techniques to hide data movement latency by concurrently executing all stages of the GEMM hierarchy within the loop and feeding the output of each stage to its dependent stage in the next iteration, as shown in Figure 8.

The GEMM CUDA kernel issues three concurrent operation streams within the pipeline, corresponding to the data flow stages in the GEMM hierarchy (Figure 1). The relative sizes of each stage in the illustration represent the lengths of operation latencies, with the orange arrows highlighting the data dependencies between each stream stage. After data is stored in shared memory, __syncthreads() is called to synchronize all warps so they can read from shared memory without contention. The final mathematical stage of the pipeline overlaps with the loading operation from shared memory that feeds data into the first mathematical stage of the next main loop iteration.

In practice, CUDA programmers achieve instruction-level concurrency between pipelined stages by interleaving CUDA statements for each stage in the program text and relying on the CUDA compiler to issue the correct instruction scheduling in the compiled code. Extensive use of #pragma unroll and compile-time constants enables the CUDA compiler to unroll loops and map array elements to registers, both of which are crucial for tunable and efficient implementations. For examples, see <code>block_task::consume_tile().

We use double buffering at each level of the GEMM hierarchy, allowing upstream pipelined stages to write data to shared memory or registers while dependent pipelined stages load data from their storage elements. Notably, this eliminates the need for a second __syncthreads(), as one shared memory buffer is written while another is read. The cost of double buffering is twice the shared memory capacity and twice the number of registers used to hold shared memory fetches.

The actual amount of latency hiding achievable depends on the sizes of thread blocks, warps, and thread blocking, as well as the throughput of the active mathematical functional units within the SM. While larger blocking can yield more opportunities for data reuse and potentially provide more latency hiding, the physical capacity of the SM register file and shared memory limits the maximum blocking size. Fortunately, NVIDIA GPUs have sufficient storage resources to execute sufficiently large GEMM blocks, making their performance limited by mathematical computations!

Detailed Components of the CUTLASS Library

CUTLASS is the implementation of the hierarchical GEMM structure as CUDA C++ template classes. We hope these templates can be included in existing device-side CUDA kernels and functions, but we also provide an example kernel and launch interface for quick start. Similar to CUB, extensive use of template parameters and compile-time constants makes CUTLASS tunable and flexible.

CUTLASS implements the abstractions required for efficient GEMM implementations. Specialized “block loaders” efficiently move data from global memory to shared memory, adapting to the layout of the source data while also achieving efficient, conflict-free loading into registers. For certain layouts, IGEMM requires some data reorganization to target CUDA’s four-element integer dot product instructions, which is done when data is stored to SMEM.

CUTLASS GEMM Device Functions

The following example from dispatch.h defines a <code>block_task type and instantiates a floating-point data GEMM for a hypothetical column-major input matrix. The block_task_policy_t defines the GEMM blocking size and is discussed in detail in the next section.

/// CUTLASS SGEMM Example
__global__ void gemm_kernel(
    float *C, 
    float const *A, 
    float const *B, 
    int M, 
    int N, 
    int K) {

    // Define GEMM blocking size - discussed in the next section
    typedef block_task_policy <
        128, // BlockItemsY: height of a tile
        32,  // BlockItemsX - width of a tile
        8,   // ThreadItemsY - height of a thread block
        4,   // ThreadItemsX - width of a thread block
        8,   // BlockItemsK - depth of a tile
        true, // UseDoubleScratchTiles - whether to use double buffering for SMEM
        block_raster_enum::Default // Block rasterization strategy
    > block_task_policy_t;

    // Define post-processing functor
    typedef gemm::blas_scaled_epilogue<float, float, float> epilogue_op_t ;

    // Define block_task type.
    typedef block_task < 
        block_task_policy_t, 
        float, 
        float, 
        matrix_transform_t::NonTranspose, 
        4, 
        matrix_transform_t::NonTranspose, 
        4, 
        epilogue_op_t, 
        4, 
        true 
    > block_task_t;

    // Declare statically allocated shared storage
    __shared__ block_task_t::scratch_storage_t smem;

    // Construct and run task
    block_task_t(
        reinterpret_cast<block_task_t::scratch_storage_t>(&smem),
        &smem,
        A,
        B,
        C,
        epilogue_op_t(1, 0),
        M,
        N,
        K).run();
}

The shared memory allocation smem is used by the block_task_t instance to store the thread block-level tiles in the matrix multiplication computation.

epilogue_op_t is a template parameter that specifies a functor used to update the output matrix after the matrix multiplication operation is complete. This allows you to easily combine matrix multiplication with custom element-wise operations, which we will describe in more detail later. CUTLASS provides the <code>gemm::blas_scaled_epilogue functor implementation for computing the familiar GEMM operation C = alpha * A * B + beta * C (defined in <code>epilogue_function.h).

CUTLASS GEMM Strategies

CUTLASS organizes compile-time constants specifying the blocking sizes at each level of the GEMM hierarchy into specializations of the gemm::block_task_policy template, which has the following declaration.

template <
    int BlockItemsY,            /// Height of a tile in matrix C
    int BlockItemsX,            /// Width of a tile in matrix C
    int ThreadItemsY,           /// Height of thread block in C
    int ThreadItemsX,           /// Width of thread block in C
    int BlockItemsK,            /// Number of K partition subgroups in a block
    bool UseDoubleScratchTiles, /// Whether to use double buffering for shared memory
    grid_raster_strategy::kind_t RasterStrategy  /// Grid rasterization strategy
> struct block_task_policy;

Several effective GEMM blocking structure strategies are defined in dispatch_policies.h, and we show one of the strategies below. This strategy decomposes the matrix multiplication operation into CUDA blocks, with each block spanning a <code>128×32 tile of the output matrix. The thread block sizes for storing <code>A and <code>B are <code>128×8 and <code>8×32, respectively. This strategy is optimized for GEMM computations where the <code>C matrix is relatively narrow in the <code>N dimension.

/// GEMM task policy specialization for tall skinny SGEMM
template <>
struct gemm_policy<float, float, problem_size_t::Tall> :
    block_task_policy<
        128,       // BlockItemsY - height of a tile
        32,        // BlockItemsX - width of a tile
        8,         // ThreadItemsY - height of a thread block
        4,         // ThreadItemsX - width of a thread block
        8,         // BlockItemsK - depth of a tile
        true,      // UseDoubleScratchTiles - whether to use double buffering for SMEM
        grid_raster_strategy::Default>   // Grid rasterization strategy
{};

The sizes of the thread block fragments are ThreadItemsY×1 and <code>ThreadItemsX×1. In the example above, these are the <code>8×1 vector from <code>A and the <code>4×1 vector from <code>B.

Once the policy type is defined, we can define the type of gemm::block_task, which is a CUTLASS GEMM. This template has the following parameter list.

template <
    /// Parameterization of block_task_policy
    typename block_task_policy_t,

    /// Value type of the multiplicands (matrices A and B)
    typename value_t,
    /// Value type of the accumulator (matrix C and scalars)
    typename accum_t,

    /// Layout enumeration for matrix A
    matrix_transform_t::kind_t  TransformA,
    
    /// Alignment of operand A (in bytes)
    int LdgAlignA,
 
    /// Layout enumeration for matrix B
    matrix_transform_t::kind_t  TransformB,

    /// Alignment of operand B (in bytes)
    int LdgAlignB,

    /// Post-processing functor applied to the matrix product
    typename epilogue_op_t,

    /// Alignment of operand C (in bytes)
    int LdgAlignC,

    /// Whether GEMM supports matrices of sizes not multiples of BlockItems{XY}
    bool Ragged              
> struct block_task;

value_t and <code>accum_t specify the types of the source operands and the accumulator matrix, respectively. <code>TransformA and <code>TransformB specify the layouts of operands <code>A and <code>B. While we do not discuss matrix layouts in detail, CUTLASS supports all combinations of row-major and column-major input matrices.

LdgAlignA and <code>LdgAlignB specify the guaranteed alignments, allowing CUTLASS device code to use vector memory operations. For example, 8-byte alignment allows CUTLASS to load elements of type <code>float in the form of double-element vectors. This reduces code size and improves performance by reducing the number of memory operations executed in the GPU. More critically, <code>Ragged handling indicates whether the matrix dimensions can be of arbitrary sizes (satisfying alignment guarantees). If this template parameter is <code>false, then the dimensions of matrices <code>A, <code>B, and <code>C must be multiples of the blocking parameters in <code>block_task_policy.

Fusing Element-wise Operations: An Example with SGEMM

Deep learning computations often perform simple element-wise operations after GEMM computations, such as computing activation functions. These bandwidth-limited layers can be fused at the end of the GEMM operation to eliminate additional kernel launches and avoid round trips to global memory.

The following example demonstrates a simple application of the GEMM template that adds a bias term to the scaled matrix product and then applies the ReLU function to clamp the result to non-negative values. By separating the post-processing into a functor, it becomes straightforward to pass parameters such as pointers to additional matrices and tensor parameters or additional scaling factors, without burdening the GEMM implementation.

First, we define a class that implements the gemm::epilogue_op concept. The constructor and other methods are not shown here, but the element-wise implementations of the bias and ReLU operations are shown in the implementation of the function call operator.

template <typename accum_t, typename scalar_t, typename output_t>
struct fused_bias_relu_epilogue {

    // Data members pass additional parameters to the post-processing
    scalar_t const *Bias;
    accum_t threshold;

    /// Host and device callable constructor, initializing data members
    inline __device__ __host__
    fused_bias_relu_epilogue(
        scalar_t const *Bias,
        accum_t threshold
    ): Bias(Bias), threshold(threshold) { }

    /// Apply bias + ReLU operation
    inline __device__ __host__
    output_t operator()(
        accum_t accumulator,  /// Element of the matrix product result
        output_t c,           /// Element of the source accumulator matrix C
        size_t idx            /// Index of the c element; can be used to load elements from other
                              /// matrices of the same structure
        ) const {

        // Compute result by scaling the matrix product, adding bias, and adding the scaled accumulator element.
        accum_t result = output_t(
            alpha * scalar_t(accumulator) +
            Bias[idx] +                    // Load and add bias
            beta * scalar_t(c)
        );

        // Apply clamping function
        return max(threshold, result);
    }
};

Then we apply this operator as a post-processing operation.

// New: Define the type of the custom post-processing functor
typedef fused_bias_relu_epilogue_t<float, float, float> 
    bias_relu_epilogue_t;

/// Compute GEMM fused with bias and ReLU operation
__global__ void gemm_bias_relu(
    ...,                                    /// GEMM parameters not shown
    bias_relu_epilogue_t bias_relu_op) {    /// bias_relu_op constructed by the caller

    // Define block_task type.
    typedef block_task<
        block_task_policy_t,          // Same policy as previous example
        float,
        float,
        matrix_transform_t::NonTranspose,
        4,
        matrix_transform_t::NonTranspose,
        4,
        bias_relu_epilogue_t,         // New: Custom post-processing functor type
        4,
        true
    > block_task_t ;

    // Declare statically allocated shared storage
    __shared__ block_task_t::scratch_storage_t smem;
    
    // Construct and run task
    block_task_t(
        reinterpret_cast<block_task_t::scratch_storage_t>(&smem),
        &smem,
        A,
        B,
        C,
        bias_relu_op,                 // New: Custom post-processing object
        M,
        N,
        K).run();
}

This simple example illustrates the value of combining general programming techniques with efficient GEMM implementations.

Performance of Tesla V100 (Volta)

CUTLASS is highly efficient, with performance for scalar GEMM computations comparable to cuBLAS. Figure 9 shows the performance of CUTLASS compiled with CUDA 9.0 on the NVIDIA Tesla V100 GPU for large matrix dimensions (M=10240, <code>N=K=4096). Figure 9 displays the relative performance of CUTLASS for each supported computation data type and all row-major and column-major layout arrangements of the input operands compared to cuBLAS.

In most cases, CUTLASS C++ achieves performance within a few percent of the hand-tuned assembly kernel performance in cuBLAS. For WMMA GEMM (WGEMM in Figure 9), CUTLASS has not yet reached the same performance as cuBLAS, but we are closely collaborating with the CUDA compiler and GPU architecture teams to develop techniques to achieve similar performance levels in CUDA code.

Getting Started with CUTLASS

In this blog post, we have not covered many interesting details, so we recommend checking out the CUTLASS repository^[7] and trying CUTLASS for yourself. The cutlass_test example program demonstrates how to invoke CUTLASS GEMM kernels, validate their results, and measure their performance. We look forward to our future updates and welcome your feedback or contact through the comments below!

Acknowledgments: Special thanks to Joel McCormack for his technical insights and elucidation, particularly regarding NVIDIA microarchitecture and the techniques employed in cuBLAS and cuDNN.

We hope this detailed analysis in English helps you gain a deeper understanding of CUTLASS!

This article is a translation and reconstruction of CUTLASS: Fast Linear Algebra in CUDA C++^[8]

References

[1] Matrix Multiplication: https://arxiv.org/abs/1410.0759[2] Fast Fourier Transform (FFT): https://arxiv.org/abs/1312.5851[3] Winograd Method: https://arxiv.org/abs/1509.09308[4] MAGMA Library: http://icl.cs.utk.edu/magma/index.html[5] Programming Tensor Cores with CUDA 9: https://developer.nvidia.com/zh-cn/blog/programming-tensor-cores-cuda-9/[6] CUTLASS Repository: https://github.com/NVIDIA/cutlass[7] CUTLASS: Fast Linear Algebra in CUDA C++: https://developer.nvidia.com/blog/cutlass-linear-algebra-cuda/

CUTLASS: Fast Linear Algebra in CUDA C++

Table of Contents

Introduction to CUTLASS

Efficient Matrix Multiplication on GPUs

Hierarchical GEMM Structure

Thread Block Blocking

Warp Blocking

Thread Blocking

WMMA GEMM: Utilizing Tensor Core

Complete GEMM Process

Software Pipelining: Hiding Latency

Detailed Components of the CUTLASS Library

CUTLASS GEMM Device Functions

CUTLASS GEMM Strategies

Fusing Element-wise Operations: An Example with SGEMM

Performance of Tesla V100 (Volta)

Getting Started with CUTLASS

References

Leave a Comment Cancel reply

Table of Contents

Introduction to CUTLASS

Efficient Matrix Multiplication on GPUs

Hierarchical GEMM Structure

Thread Block Blocking

Warp Blocking

Thread Blocking

WMMA GEMM: Utilizing Tensor Core

Complete GEMM Process

Software Pipelining: Hiding Latency

Detailed Components of the CUTLASS Library

CUTLASS GEMM Device Functions

CUTLASS GEMM Strategies

Fusing Element-wise Operations: An Example with SGEMM

Performance of Tesla V100 (Volta)

Getting Started with CUTLASS

References

Related posts

Leave a Comment Cancel reply