Matrix multiplication on GPU using CUDA with CUBLAS, CURAND and Thrust
Posted on May 31, 2012 by Sol

The code for this tutorial is on GitHub: https://github.com/sol-prog/cuda_cublas_curand_thrust.

Matrix multiplication is an essential building block for numerous numerical algorithms, for this reason most numerical libraries implements matrix multiplication. One of the oldest and most used matrix multiplication implementation GEMM is found in the BLAS library. While the reference BLAS implementation is not particularly fast there are a number of third party optimized BLAS implementations like MKL from Intel, ACML from AMD or CUBLAS from NVIDIA.

In this post I’m going to show you how you can multiply two arrays on a CUDA device with CUBLAS. A typical approach to this will be to create three arrays on CPU (the host in CUDA terminology), initialize them, copy the arrays on GPU (the device on CUDA terminology), do the actual matrix multiplication on GPU and finally copy the result on CPU. Our first example will follow the above suggested algorithm, in a second example we are going to significantly simplify the low level memory manipulation required by CUDA using Thrust which aims to be a replacement for the C++ STL on GPU.

Let’s start by allocating space for our three arrays on CPU and GPU:

 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
#include <cstdlib>

int main() {
    // Allocate 3 arrays on CPU
    int nr_rows_A, nr_cols_A, nr_rows_B, nr_cols_B, nr_rows_C, nr_cols_C;

    // for simplicity we are going to use square arrays
    nr_rows_A = nr_cols_A = nr_rows_B = nr_cols_B = nr_rows_C = nr_cols_C = 3;
    
    float *h_A = (float *)malloc(nr_rows_A * nr_cols_A * sizeof(float));
    float *h_B = (float *)malloc(nr_rows_B * nr_cols_B * sizeof(float));
    float *h_C = (float *)malloc(nr_rows_C * nr_cols_C * sizeof(float));

    // Allocate 3 arrays on GPU
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A,nr_rows_A * nr_cols_A * sizeof(float));
    cudaMalloc(&d_B,nr_rows_B * nr_cols_B * sizeof(float));
    cudaMalloc(&d_C,nr_rows_C * nr_cols_C * sizeof(float));

    // ....

    //Free GPU memory
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);  

    // Free CPU memory
    free(h_A);
    free(h_B);
    free(h_C);    

    return 0;
}

Please note the way in which we allocate memory for data on CPU using C’s malloc (line 10) and GPU using CUDA’s cudaMalloc (line 16), at the end of the main function we can free the device memory with cudaFree. The above code uses d_ for prefixing the arrays located on GPU and h_ for the host (CPU) arrays. We use 3x3 arrays in this example for simplicity, in a real application you should use much larger arrays for using the device efficiently.

Next step is to initialize the arrays A and B with some values, we are going to fill them with random numbers on the GPU using the CURAND library. In order to do this we could write a function that receives as input a device array, after the function is executed the array will be filled with random nubers:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
...
#include <curand.h>
...

// Fill the array A(nr_rows_A, nr_cols_A) with random numbers on GPU
void GPU_fill_rand(float *A, int nr_rows_A, int nr_cols_A) {
    // Create a pseudo-random number generator
    curandGenerator_t prng;
    curandCreateGenerator(&prng, CURAND_RNG_PSEUDO_DEFAULT);

    // Set the seed for the random number generator using the system clock
    curandSetPseudoRandomGeneratorSeed(prng, (unsigned long long) clock());

    // Fill the array with random numbers on the device
    curandGenerateUniform(prng, A, nr_rows_A * nr_cols_A);
}

...

If the arrays A and B are already filled with meaningful data, on CPU, or you simply want to fill them with random data on CPU and avoid the use of CURAND, you will need to transfer the host arrays to GPU:

1
2
    cudaMemcpy(d_A,h_A,nr_rows_A * nr_cols_A * sizeof(float),cudaMemcpyHostToDevice);
    cudaMemcpy(d_B,h_B,nr_rows_B * nr_cols_B * sizeof(float),cudaMemcpyHostToDevice);

As a side note, BLAS uses internally column-major order storage (Fortran order) and not the typical row-major storage used in C or C++ for multidimensional arrays.

Now that the arrays A and B are initialized and transfered on GPU we could write a function that will do the actual multiplication:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
...
#include <cublas_v2.h>
...

// Multiply the arrays A and B on GPU and save the result in C
// C(m,n) = A(m,k) * B(k,n)
void gpu_blas_mmul(const float *A, const float *B, float *C, const int m, const int k, const int n) {
    int lda=m,ldb=k,ldc=m;
    const float alf = 1;
    const float bet = 0;
    const float *alpha = &alf;
    const float *beta = &bet;

    // Create a handle for CUBLAS
    cublasHandle_t handle;
    cublasCreate(&handle);

    // Do the actual multiplication
    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);

    // Destroy the handle
    cublasDestroy(handle);
}

Observation - If you need to do more than one matrix multiplication in your code it is advisable to move the create/destroy handle code (lines 15 - 16 and 22) from the above function in the main function, and use the same handle for all multiplications. In this case the gpu_blas_mmul function became:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
void gpu_blas_mmul(cublasHandle_t &handle, const float *A, const float *B, float *C, const int m, const int k, const int n) {
    int lda=m,ldb=k,ldc=m;
    const float alf = 1;
    const float bet = 0;
    const float *alpha = &alf;
    const float *beta = &bet;

    // Do the actual multiplication
    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
}

We can also implement an utility function for printing the result of the multiplication on the standard output:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
//Print matrix A(nr_rows_A, nr_cols_A) storage in column-major format
void print_matrix(const float *A, int nr_rows_A, int nr_cols_A) {

    for(int i = 0; i < nr_rows_A; ++i){
        for(int j = 0; j < nr_cols_A; ++j){
            std::cout << A[j * nr_rows_A + i] << " ";
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;
}

Using the above pieces we present here for reference the complete main function:

 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
int main() {
    // Allocate 3 arrays on CPU
    int nr_rows_A, nr_cols_A, nr_rows_B, nr_cols_B, nr_rows_C, nr_cols_C;

    // for simplicity we are going to use square arrays
    nr_rows_A = nr_cols_A = nr_rows_B = nr_cols_B = nr_rows_C = nr_cols_C = 3;
    
    float *h_A = (float *)malloc(nr_rows_A * nr_cols_A * sizeof(float));
    float *h_B = (float *)malloc(nr_rows_B * nr_cols_B * sizeof(float));
    float *h_C = (float *)malloc(nr_rows_C * nr_cols_C * sizeof(float));

    // Allocate 3 arrays on GPU
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A,nr_rows_A * nr_cols_A * sizeof(float));
    cudaMalloc(&d_B,nr_rows_B * nr_cols_B * sizeof(float));
    cudaMalloc(&d_C,nr_rows_C * nr_cols_C * sizeof(float));

    // Fill the arrays A and B on GPU with random numbers
    GPU_fill_rand(d_A, nr_rows_A, nr_cols_A);
    GPU_fill_rand(d_B, nr_rows_B, nr_cols_B);

    // Optionally we can copy the data back on CPU and print the arrays
    cudaMemcpy(h_A,d_A,nr_rows_A * nr_cols_A * sizeof(float),cudaMemcpyDeviceToHost);
    cudaMemcpy(h_B,d_B,nr_rows_B * nr_cols_B * sizeof(float),cudaMemcpyDeviceToHost);
    std::cout << "A =" << std::endl;
    print_matrix(h_A, nr_rows_A, nr_cols_A);
    std::cout << "B =" << std::endl;
    print_matrix(h_B, nr_rows_B, nr_cols_B);

    // Multiply A and B on GPU
    gpu_blas_mmul(d_A, d_B, d_C, nr_rows_A, nr_cols_A, nr_cols_B);

    // Copy (and print) the result on host memory
    cudaMemcpy(h_C,d_C,nr_rows_C * nr_cols_C * sizeof(float),cudaMemcpyDeviceToHost);
    std::cout << "C =" << std::endl;
    print_matrix(h_C, nr_rows_C, nr_cols_C);

    //Free GPU memory
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);  

    // Free CPU memory
    free(h_A);
    free(h_B);
    free(h_C);

    return 0;
}

For compiling the above example open a Terminal (Mac and Linux) or a Visual Studio Command Prompt window (for Windows) and write:

1
nvcc mmul_1.cu -lcublas -lcurand -o mmul_1

This is the result of running the code on my Mac:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
sols-MacBook-Pro:Desktop sol$ nvcc mmul_1.cu -lcublas -lcurand -o mmul_1
sols-MacBook-Pro:Desktop sol$ ./mmul_1
A =
0.365931 0.887901 0.48295 
0.610939 0.453606 0.978245 
0.626738 0.275708 0.804269 

B =
0.975249 0.504392 0.520455 
0.423145 0.328515 0.128649 
0.767074 0.0733948 0.217057 

C =
1.10304 0.511707 0.409506 
1.53815 0.528967 0.588657 
1.34482 0.465725 0.53623 

sols-MacBook-Pro:Desktop sol$ 

As mentioned in the introduction, we could use Thrust to further simplify the above code. For example we could avoid completely the need to manually manage memory on the host and device using a Thrust vector for storing our data. Reimplementing the above example with Thrust will halve the number of lines of code from the main function:

 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
int main() {
    // Allocate 3 arrays on CPU
    int nr_rows_A, nr_cols_A, nr_rows_B, nr_cols_B, nr_rows_C, nr_cols_C;

    // for simplicity we are going to use square arrays
    nr_rows_A = nr_cols_A = nr_rows_B = nr_cols_B = nr_rows_C = nr_cols_C = 3;

    thrust::device_vector<float> d_A(nr_rows_A * nr_cols_A), d_B(nr_rows_B * nr_cols_B), d_C(nr_rows_C * nr_cols_C);

    // Fill the arrays A and B on GPU with random numbers
    GPU_fill_rand(thrust::raw_pointer_cast(&d_A[0]), nr_rows_A, nr_cols_A);
    GPU_fill_rand(thrust::raw_pointer_cast(&d_B[0]), nr_rows_B, nr_cols_B);

    // Optionally we can print the data
    std::cout << "A =" << std::endl;
    print_matrix(d_A, nr_rows_A, nr_cols_A);
    std::cout << "B =" << std::endl;
    print_matrix(d_B, nr_rows_B, nr_cols_B);

    // Multiply A and B on GPU
    gpu_blas_mmul(thrust::raw_pointer_cast(&d_A[0]), thrust::raw_pointer_cast(&d_B[0]), thrust::raw_pointer_cast(&d_C[0]), nr_rows_A, nr_cols_A, nr_cols_B);

    //Print the result
    std::cout << "C =" << std::endl;
    print_matrix(d_C, nr_rows_C, nr_cols_C);

    return 0;
}

Please note in the above highlighted lines the way in which we convert a Thrust device_vector to a CUDA device pointer.

If you are interested in learning CUDA, I would recommend reading CUDA Application Design and Development by Rob Farber.

or CUDA by Example: An Introduction to General-Purpose GPU Programming by J. Sanders and E. Kandrot.

blog comments powered by Disqus