global void mulsum(int* a, int* b, int* out, int n) { int tid = threadIdx.x; extern shared int tmp[]; /* ▲

• │ ┌───────────────────────────────┐
• │ │External and shared, allocated │
• └──┤by the cuda runtime when kernel│
• │ is launched │
• └───────────────────────────────┘ */ if (tid >= n) return;

tmp[tid] = a[tid] * b[tid];

__syncthreads();

if (tid == 0) { int sum = 0; for (int i=0; i < n; i++) sum += tmp[i]; *out = sum; } }

int main() { #define N 10

int a[N]; for (int i=0; i < N; i++) a[i] = i;

int b[N]; for (int i=0; i < N; i++) b[i] = i;

int* d_a; cudaMalloc(&d_a, sizeof(int[N])); cudaMemcpy(d_a, a, sizeof(int[N]), cudaMemcpyHostToDevice);

int* d_b; cudaMalloc(&d_b, sizeof(int[N])); cudaMemcpy(d_b, b, sizeof(int[N]), cudaMemcpyHostToDevice);

int* d_out; cudaMalloc(&d_out, sizeof(int));

mulsum<<<1, N, sizeof(int[N])>>>(d_a, d_b, d_out, N); /* ▲

•         │
  

• ┌────────┴────────────────────────────┐
• │Size of shared memory to be allocated│
• │ for kernel launch │
• └─────────────────────────────────────┘ */

int out; cudaMemcpy(&out, d_out, sizeof(int), cudaMemcpyDeviceToHost); printf(“%d\n”, out);

cudaFree(d_a); cudaFree(d_b); cudaFree(d_out); return 0; }

```console
$ nvcc mul-sum-reduce.cu && ./a.out
285

Notice how our shared memory is declared:

  extern __shared__ int tmp[];

It is external because we are not allocating the memory in our kernel; it’s allocated by the cuda runtime when we pass the third parameter to the kernel launch parameters:

mulsum<<<1, N, sizeof(int)*N>>>(d_a, d_b, d_out, N);
               ▲
               │
 ┌─────────────┴───────────────────────┐
 │Size of shared memory to be allocated│
 │         for kernel launch           │
 └─────────────────────────────────────┘

There can only be one segment of shared memory in a kernel launch, so the shared memory segment will be interpreted as whatever type we declare our shared memory with. In this case, it’s an array of ints. Although there is strictly one segment of shared memory in a kernel launch, you can still declare multiple variables as __shared__, so long as they all fit in the allocated shared memroy.

We also introduced another CUDA extension to the host language: __syncthreads(). __syncthreads() is a fence or barrier, a point which no thread in that block can cross until all threads have reached it. __syncthreads() There are many other CUDA primitives for atomic, and synchronization operations, such as atomicAdd.

CUDA Mat-Vec Multiply

We again return to our matvecmul example, this time armed with some knowledge about the CUDA runtime and some software and hardware abstractions.

#include <cstdio>
#include <cuda.h>
#include <cuda_runtime.h>

__global__ void matvecmul(int* mat, int* vec, int* outv,
                          int m, int n) {

  int rowidx = blockIdx.x;
  int colidx = threadIdx.x;

  extern __shared__ int tmp[];

  if (colidx < n && rowidx < m) {
    tmp[colidx] = mat[colidx + (rowidx * n)] * vec[colidx];

    __syncthreads();

    if (colidx == 0) {
      int sum = 0;
      for (int i=0; i < n; i++)
        sum += tmp[i];
      outv[rowidx] = sum;
    }
  }
}

int main() {
  #define M 10
  #define N 15

  int a[M*N];
  for (int i=0; i < M*N; i++) a[i] = i;

  int b[N];
  for (int i=0; i < N; i++) b[i] = i;

  int* d_a;
  cudaMalloc(&d_a, sizeof(int[M*N]));
  cudaMemcpy(d_a, a, sizeof(int[M*N]), cudaMemcpyHostToDevice);

  int* d_b;
  cudaMalloc(&d_b, sizeof(int[N]));
  cudaMemcpy(d_b, b, sizeof(int[N]), cudaMemcpyHostToDevice);

  int* d_c;
  cudaMalloc(&d_c, sizeof(int[M]));

  matvecmul<<<M, N, sizeof(int[N])>>>(d_a, d_b, d_c, M, N);

  int c[M];
  cudaMemcpy(c, d_c, sizeof(int[M]), cudaMemcpyDeviceToHost);

  cudaFree(d_a);
  cudaFree(d_b);
  cudaFree(d_c);
  return 0;
}

After adding some printing code to our example above, we get the following:

$ nvcc matvecmul.cu && ./a.out
Matrix:
   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14
  15   16   17   18   19   20   21   22   23   24   25   26   27   28   29
  30   31   32   33   34   35   36   37   38   39   40   41   42   43   44
  45   46   47   48   49   50   51   52   53   54   55   56   57   58   59
  60   61   62   63   64   65   66   67   68   69   70   71   72   73   74
  75   76   77   78   79   80   81   82   83   84   85   86   87   88   89
  90   91   92   93   94   95   96   97   98   99  100  101  102  103  104
 105  106  107  108  109  110  111  112  113  114  115  116  117  118  119
 120  121  122  123  124  125  126  127  128  129  130  131  132  133  134
 135  136  137  138  139  140  141  142  143  144  145  146  147  148  149

Vector:
   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14

Output:
1015 2590 4165 5740 7315 8890 10465 12040 13615 15190

In our CUDA C example, we launch a block for each row of our matrix. This way, we can share memory between each thread operating on a given row of the matrix. A single thread per row can then perform the sum reduction and assign the value to the index in the output vector.

CUDA C++

Naive Approach

In our previous CUDA C example, we weren’t really using CUDA C perse, but CUDA C++. If you read the introductory chapter in the official NVIDIA CUDA Programming guide, you’ll see that CUDA C is really just CUDA mixed with the common language subset between C and C++ on the host. We were using CUDA C++ the whole time, we just restricted ourselves to the C subset of C++ for simplicity.

#include <cstdio>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/sequence.h>
#include <cuda.h>
#include <cuda_runtime.h>

/* ┌───────────────────────────┐
   │Using thrust's pointer type│
   │     instead of int*       │
   └───────────────────────────┘ */
__global__ void matvecmul(
    thrust::device_ptr<int> mat,
    thrust::device_ptr<int> vec,
    thrust::device_ptr<int> outv,
    int m, int n) {

  int rowidx = blockIdx.x;
  int colidx = threadIdx.x;

  extern __shared__ int tmp[];

  if (colidx < n && rowidx < m)
  {
    tmp[colidx] = mat[colidx + (rowidx * n)] * vec[colidx];

    __syncthreads();

    if (colidx == 0) {
      int sum = 0;
      for (int i=0; i < n; i++)
        sum += tmp[i];
      outv[rowidx] = sum;
    }
  }
}

int main() {
  #define M 10
  #define N 15

  /* ┌────────────────────────────────┐
     │ Using thrust's host and device │
     │    vectors over raw arrays.    │
     │No longer need to use cudaMemcpy│
     │        or cudaMalloc!          │
     └────────────────────────────────┘ */
  thrust::device_vector<int> a(M*N);
  thrust::sequence(a.begin(), a.end(), 0);

  thrust::device_vector<int> b(N);
  thrust::sequence(b.begin(), b.end(), 0);

  thrust::device_vector<int> c(M, 0);

  matvecmul<<<M, N, sizeof(int[N])>>>(a.data(), b.data(), c.data(), M, N);

  /* ┌────────────────────────────┐
     │The assignment operator will│
     │   perform the cudaMemcpy   │
     └────────────────────────────┘ */
  thrust::host_vector<int> out = c;

  puts("Output:");
  for (int i=0; i < M; i++)
    printf("%d ", out[i]);
  puts("");
  return 0;
}

$ nvcc thrust-ex.cu && ./a.out
Output:
1015 2590 4165 5740 7315 8890 10465 12040 13615 15190

As you can see, the code looks quite similar, except for the lack of memory management. This is hiding a few extra details as well. In our original CUDA example, we first allocate and assign to memory on the host before copying it to the device. In this example, we allocate memory on the device first, and perform assignments on the device.

In this line, we allocate device memory for our matrix, and assign values to it on the device.

  thrust::device_vector<int> a(M*N);
  thrust::sequence(a.begin(), a.end(), 0);

thrust::sequence is almost identical to std::iota or for (int i=0; i < N; i++) vec[i] = i;, except that it may execute on the device. In this new example, we launch three kernels instead of one: one for each call to thrust::sequence, and one for our manual kernel launch. You can look at the details of the ptx assembly in Compiler Explorer here.

More Sophisticated Approach

Remember all that fuss about fundamental algorithms in the earlier sections? How our fundamental algorithm here is a transform-reduce?

Well, in our first-pass CUDA implementation, we don’t really use this to our advantage. Our kernel contains the following lines:

    if (colidx == 0) {
      int sum = 0;
      for (int i=0; i < n; i++)
        sum += tmp[i];
      outv[rowidx] = sum;
    }

While we used a raw loop at least once per block, an optimized reduction will perform something like the following:

Extremely performant reductions are actually quite hard to get right - it’s easy to get some parallelism in a reduction, but it takes significant effort to truly maximize the speed you can get from a GPU. This slidedeck from the NVIDIA developer blog details various approaches to optimizing a reduction operation on a GPU. Thrust’s reduce implementation will even select different strategies and launch parameters based on the sizes of the data it operates on.

The point of this Thrust discussion is not to dissuade you from writing raw CUDA kernels - it’s to dissuage you from doing it too early. In the majority of cases, it’s likely that using a library around raw CUDA kernels will result in faster code and less development time. Once you have already written your code using known algorithms, once you have tested your code to demonstrate its correctness, once you have profiled your code to demonstrate where the performance bottlenecks are on the target architectures you care about, then it makes sense to write raw CUDA kernels.

Use sane defaults, only optimize after profiling and testing

So let’s try again, using Thrust’s parallel algorithms to compute the reductions for each row of the matrix-vector multiplication (godbolt link here):

#include <cstdio>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/tuple.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/sequence.h>

__global__
void broadcast_to_matrix(thrust::device_ptr<int> mat,
    thrust::device_ptr<int> vec,
    int m, int n) {
  const auto col = blockIdx.x;
  const auto row = threadIdx.x;

  if (row < m and col < n)
    mat[col+(row*n)] = vec[col];
}

int main() {
  #define M 10
  #define N 15

  thrust::device_vector<int> a(M*N);
  thrust::sequence(a.begin(), a.end(), 0);

  thrust::device_vector<int> b(N);
  thrust::sequence(b.begin(), b.end(), 0);

  thrust::device_vector<int> broadcasted_b(M*N);
  broadcast_to_matrix<<<N, M>>>(broadcasted_b.data(), b.data(), M, N);

  thrust::host_vector<int> c(M);

  thrust::zip_iterator iter(thrust::make_tuple(a.begin(), broadcasted_b.begin()));

  for (int i=0; i < M; i++)
    c[i] = thrust::transform_reduce(iter+(i*N), iter+(i*N)+N,
        [] __device__ (auto tup) -> int {
          return thrust::get<0>(tup) * thrust::get<1>(tup);
        },
        0,
        thrust::plus<int>()
        );

  puts("Output:");
  for (int i=0; i < M; i++)
    printf("%d ", c[i]);
  puts("");
  return 0;
}

Since we’re using some more fancy C++ features, we have to pass some extra flags to the compiler:

$ nvcc -std=c++17 --extended-lambda fancy.cu && ./a.out
Output:
1015 2590 4165 5740 7315 8890 10465 12040 13615 15190

The first kernel is launched manually, since we’re just broadcasting values from the vector to another vector to match the shape of our matrix. We perform this step so we can create a zip iterator from the matrix and the broadcasted vector:

  thrust::zip_iterator iter(thrust::make_tuple(a.begin(), broadcasted_b.begin()));

This means we can feed the single iterator into our transform_reduce operation. Elements obtained by the zip iterator are passed into our lambda function, and we simply multiply the two values together, before using the plus functor to reduce the vector of intermediate values for a given row into a scalar:

  for (int i=0; i < M; i++)
    c[i] = thrust::transform_reduce(iter+(i*N), iter+(i*N)+N,
        [] __device__ (auto tup) -> int {
          return thrust::get<0>(tup) * thrust::get<1>(tup);
        },
        0,
        thrust::plus<int>()
        );

If we want to use threading on the host as well, we can even use an OpenMP directive:

#pragma openmp parallel for
  for (int i=0; i < M; i++)
    c[i] = thrust::transform_reduce(iter+(i*N), iter+(i*N)+N,
        [] __device__ (auto tup) -> int {
          return thrust::get<0>(tup) * thrust::get<1>(tup);
        },
        0,
        thrust::plus<int>()
        );

We’ll have to tell nvcc to pass the -fopenmp flag to the host compiler:

$ nvcc -std=c++17 --extended-lambda -Xcompiler -fopenmp fancy.cu

The Best Tool for the Job

We have by now hopefully learned that we should use the most specialized tool for the job, and we should write kernels by hand only when we’re sure we can do better than your libraries of choice. We can take this principle one step further with a little extra knowledge of our problem.

A matrix-vector product is a very common linear algebra operation, and a member of the Basic Linear Algebra Subroutines interface, which CUDA provides a library for (CUBLAS). Because this is such a common operation, NVIDIA provides an extremely fast implementation - far more optimized than anything we would write by hand.

This knowledge of our problem leads us to using the most appropriate library, and likely to the fastest solution.

#include <cstdio>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/sequence.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>

int main() {
  #define M 10
  #define N 15

  cublasHandle_t ch;
  cublasCreate(&ch);

  thrust::device_vector<double> a(M*N);
  thrust::sequence(a.begin(), a.end(), 0);

  const double alpha = 1.0;
  const double beta = 0.0;

  thrust::device_vector<double> b(N);
  thrust::sequence(b.begin(), b.end(), 0);

  thrust::device_vector<double> c(M);

  #define PTR(x) thrust::raw_pointer_cast(x.data())
  cublasDgemv(
      ch,
      CUBLAS_OP_T,
      N, M,
      &alpha,
      PTR(a), N,
      PTR(b), 1,
      &beta,
      PTR(c), 1
      );
  #undef PTR

  thrust::host_vector<double> hc = c;

  puts("Output:");
  for (int i=0; i < M; i++)
    printf("%.1f ", hc[i]);
  puts("");

  cublasDestroy(ch);

  return 0;
}

$ nvcc cublas.cu -lcublas && ./a.out
Output:
1015.0 2590.0 4165.0 5740.0 7315.0 8890.0 10465.0 12040.0 13615.0 15190.0

Conclusion

BQN Example

Personally, I use BQN to prototype solutions to problems and to better understand the fundamental algorithms at play; you don’t have to know an APL in order to understand this, but it might be helpful. Feel free to skip this section; it is not critical to understanding the concepts.

Here’s a permalink to the BQN snippet.

   # Same matrix as in our C example
   mat ← 3‿3⥊↕10
┌─       
╵ 0 1 2  
  3 4 5  
  6 7 8  
        ┘
   # Same vector as in our C example
   vec ← 3⥊2
⟨ 2 2 2 ⟩

   +˝⎉1 mat×vec
⟨ 6 24 42 ⟩

The core algorithm is seen in the final expression:

+˝⎉1 mat×vec
▲    ▲
│    │     ┌───────────────────────────┐
│    └─────┤Multiply rows of mat by vec│
│          │        element-wise       │
│          └───────────────────────────┘
│     ┌─────────────────────────┐
│     │Sum-reduce rows of matrix│
└─────┤ resulting from mat×vec  │
      └─────────────────────────┘

Alternatively:

Links, References, Additional Reading

• BQN matvecmul example
• Matrix-Vector Product image
• UT Austin slides on loop-carried dependencies and parallelism
• How to Write Parallel Programs: A First Course
• Thrust parallel algorithms library
• ADSP podcast episode from the lead HPC architect at NVIDIA discussing speed vs efficiency
• Bryce Adelstein Lelbach’s talk on C++ standard parallelism
• Kokkos mdspan single header
• CUDA C Introduction Slides
• More sophisticated CUDA matrix-vector product implementations
• Slides on CUDA reduction operation