HPC - Lecture 14 - CUDA libraries

Overview: The cuBLAS Library

Cuda-X suite: https://developer.nvidia.com/cuda/cuda-x-libraries

cuBLAS introduction

The cuBLAS library provides basic linear algebra operations in the CUDA environment, based on the BLAS (Basic Linear Algebra Subprograms) library.

It provides basic routine with optimized memory access.

Today, all scientific computing libraries rely on the BLAS library for matrix and vector processing!

Since cuBLAS was originally written for Fortran to maintain compatibility with it, cuBLAS stores matrices by columns and uses 1-based indexing.

PBLAS (Parallel BLAS)

Various parallel versions of basic linear algebra algorithms have been developed over years.

PBLAS is a parallel library for basic linear algebra operations in a MIMD distributed-memory environment.

PBLAS uses:

BLAS for computing
BLACS for communications

BLACS (Basic Linear Algebra Communication Subprograms) is a communication library developed for different environments (including MPI).

It is a tool that facilitates the development of parallel codes for linear algebra operations in a distributed-memory environment. It is based on communication of 2D data structures → MPI-like routines for sending data and collecting results, specifically defined for vectors and matrices.

PLASMA (Parallel Linear Algebra for Scalable Multi-core Architectures)

Is a library introduced in 2006 for multicore architectures.

The PLASMA library provides a toolkit of routines for basic linear algebra operations (like BLAS). In addition, it provides some routines for solving systems of linear equations, for example, by LU factorization. → Scientific Computing

cuBLAS architecture

The cuBLAS library provides a considerable number of functions that allow matrice and vectors of different data type to be processed .

The library is organized into three levels:

Level 1: Operations between vectors (e.g., cublas<t>dot).
Level 2: Operations between matrices and vectors (e.g., cublas<t>gemv).
Level 3: Operations between matrices (e.g., cublas<t>gemm).

For each function, there are 4 different versions based on data type: single precision, double precision, single precision complex, and double precision complex.

Moving data between devices

Copy of vector from host to device - host → device

cublasSetVector(n, sizeof(type), x, incx, y, incy);

It copies n values from a host vector x to a device vector y.

type is the data type of the allocated vector.
incx, incy is the strides between consecutive values read from x and written to y.

N.B. arrays x and y must be allocated in advance.

Copy of vector from device to host - device → host

cublasGetVector(n, sizeof(type), x, incx, y, incy);

It copies n values from a device vector x to a host vector y.

type is the data type of the allocated vector.
incx, incy is the strides between consecutive values read from x and written to y.

Copy of matrix from host to device - host → device

cublasSetMatrix(rows, cols, sizeof(type), A, lda, B, ldb);

It copies rows by cols matrix A, stored in the host memory, to matrix B stored in the device memory.

lda, ldb are the leading dimensions of matrices

It’s like the previous functions but for 2D data structures

Copy of matrix from device to host - device → host

cublasGetMatrix(rows, cols, sizeof(type), A, lda, B, ldb);

It copies rows by cols matrix A, stored in the device memory, to matrix B stored in the host memory.

lda, ldb: the leading dimensions of matrices

N.B. This function is very similar to the previous one, but performs the inverse operation.

Matrix management in cuBLAS

To ensure a full compatibility with the fortran language, when using the cuBLAS library the matrix entries can be stored by columns and starting the indexing by 1, instead of 0.

By storing matrices in this way, some inconveniences can be avoided.

To use the cuBLAS library, the header file cublas_v2.h must be included.

cuBLAS handle

To use the cuBLAS library, a handle must be initialized to manage library resources, especially when using multiple host threads or multiple GPUs.

Initialization: cublasCreate(&name_handle)
Destruction: cublasDestroy(&name_handle)

The handle is a pointer, which can be declared as: cublasHandle_t name_handle In this way a special data structure is initialized for the cuBLAS library.

The handle must be passed as a parameter to every cuBLAS function.

Error Management

Developers can implement an error capture mechanism using the cublasStatus_t status type. A simple check on the variables allows us to obtain the errors that occurred.

There is a set of predefined values:

Success: CUBLAS_STATUS_SUCCESS
Errors: CUBLAS_STATUS_NON_INITIALIZED, CUBLAS_STATUS_ALLOC_FAILED, CUBLAS_STATUS_EXECUTION_FAILED, CUBLAS_STATUS_INVALID_VALUE, etc.

cuBLAS level 1

Using level 1 functions, operations between vectors and on a single vector can be performed. There are 13 functions available, of which the most suitable can be chosen based on the data type to be processed.

For each function, there are 4 different versions, one for each type of data: single precision, double precision, single precision complex and double precision complex.

dot product between vectors

The function for the dot product is: cublas<t>dot() where the parameter <t> indicates the data type.

cublasStatus_t cublasSdot(cublasHandle_t handle,
                          int n, const float *x,
                          int incx, const float *y,
                          int incy, float *result);

cublasHandle_t handle is the cuBLAS handle
n is the number of elements in the vectors x and y
incx is the stride between consecutive elements of the vector x
incy is the stride between consecutive elements of the vector y
return cublasStatus_t → Error Management

If i want to work with doubles instead of floats, i’d use:

cublasStatus_t cublasDdot(cublasHandle_t handle,
                          int n, const double *x,
                          int incx, const double *y,
                          int incy, double *result);

cuBLAS supports also cuComplex, that is a built-in datatype for single-precision complex numbers (that contains both a real and imaginary component).

cublasStatus_t cublasCdot(cublasHandle_t handle,
                          int n, const cuComplex *x,
                          int incx,
                          const cuComplex *y,
                          int incy, cuComplex *result);

use cublasZdot for the double complex numbers.

cuBLAS level 2

Using level 2 functions, operations between matrices and vectors can be performed. There are 23 functions available, from which the most suitable can be choosen, based on the data type to be processed.

Matrix-Vector Multiplication

The function for matrix-vector multiplcation is:

cublas<t>gemv

For example, for data type float:

cublasStatus_t cublasSgemv(cublasHandle_t handle,
                           cublasOperation_t trans,
                           int m, int n,
                           const float *alpha,
                           const float *A, int lda,
                           const float *y, int incy,
                           const float *beta,
                           float *r, int incr);

change the “S” in cublasStatus_t cublasSgemv with D for double, C for cuComplex and Z for cuDoubleComplex.

For each function, the parameters are:

The sizes n and m of the matrix and vector
The matrix A and the vector y
The leading dimension lda
The handle required to use the cuBLAS library
The type of operation CUBLAS_OP_X r=alpha*op(A)*y+beta*r
The scalar values alpha and beta (for a simple matrixvector multiplication, they are usually 1 and 0, respectively)
The result vector r
The values incy and incr, which indicate the stride between consecutive elements in vectors y and r (generally equal to 1).

The type of operation that the function cublas<t>gemv performs is: r = alpha * op(A) * y + beta * r

cublasOperation_t trans: specifies if matrix A is transposed (CUBLAS_OP_T) or not (CUBLAS_OP_N)
Parameters include sizes n, m, matrix A, vector y, scalar alpha/beta, and result vector r.

cuBLAS level 3

Using level 3 functions, operations between matrices can be performed.

There are 9 functions available, from which the most suitable can be chosen based on the data type to be processed. More precisely, for each function, there are 4 different versions, one for each type of data: single precision, double precision, single precision complex or double precision complex.

Matrix-Matrix Multiplication

The function for matrix-matrix multiplication is: cublas<t>gemm().

For example, for floats:

cublasStatus_t cublasSgemm(cublasHandle_t handle,
                           cublasOperation_t transa,
                           cublasOperation_t transb,
                           int m, int n, int k,
                           const float *alpha,
                           const float *A, int lda,
                           const float *B, int ldb,
                           const float *beta,
                           float *C, int ldc);

and similarly double, single-precision and double-precision complex are available changing the letter.

For each function, the parameters are:

The sizes n, m, k of matrices
The matrices A and B
The leading dimension lda
The leading dimension ldb
The handle required to use the cuBLAS library
The type of operation CUBLAS_OP_X C=alpha*op(A)*op(B)+beta*C
The scalar values alpha and beta (for a simple matrixmatrix multiplication, they are usually 1 and 0, respectively)
The matrix C
The leading dimension ldc

Performs: C = alpha * op(A) * op(B) + beta * C

transa, transb: specify if op(A)/op(B) are A/B or their transposes.

Steps to Use cuBLAS

Create a CUBLAS handle.
Initialize vectors/matrices in host memory.
Transfer data from host to device using cublasSetVector(n, sizeof(type), x, incx, y, incy) or cublasSetMatrix(rows, cols, sizeof(type), A, lda, B, ldb).
Perform the operation using the appropriate cuBLAS function.
Transfer results from device to host using cublasGetVector or cublasGetMatrix.

Example of program

Main:

#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 
void initializeArray(double *, int);
void initializeArrayRandom(double *, int);
void initializeMatrix(double *, int, int);
void initializeMatrixRandom(double *, int, int);
void printArray(double *, int);
void printMatrix(double *, int, int);
 
int main(int argn, char *argv[])
{
    int n, m;
    int random = 0;
    // Arrays on the host
    double *vector_host, *final_host, *matrix_host;
 
    // Arrays on the device
    double *vector_device, *vector_final, *matrix_device;
    // Scalar values for the product
    double alpha = 1.0f;
    double beta = 0.0f;
    // error-check variable
    cublasStatus_t status;
    // cuBLAS handle declaration
    cublasHandle_t handle;
    // cuBLAS handle creation
    status = cublasCreate(&handle);
    // Error check
    if (status != CUBLAS_STATUS_SUCCESS)
    {
        printf("cublasCreate failed\n");
        return 2;
    }
    printf("computation of matrix-vector multiplication\n"
           "using the cuBLAS library\n");
    if (argn < 4)
    {
        printf("Default settings\n");
        n = 8;
        m = 8;
    }
    else
    {
        m = atoi(argv[1]);
        n = atoi(argv[2]);
        random = atoi(argv[3]);
    }
 
    // Vector and matrix allocation on the host
    vector_host = (double *)malloc(n * sizeof(double));
    final_host = (double *)malloc(m * sizeof(double));
    matrix_host = (double *)malloc((m * n + 1) * sizeof(double));
    // allocation on the device
    cudaMalloc((void **)&vector_device, n * sizeof(double));
    cudaMalloc((void **)&vector_final, m * sizeof(double));
    cudaMalloc((void **)&matrix_device, m * n * sizeof(double));
 
    // vector initialization
    if (random == 0)
        initializeArray(vector_host, n);
    else
        initializeArrayRandom(vector_host, n);
    // matrix initialization
    if (random == 0)
        initializeMatrix(matrix_host, m, n);
    else
        initializeMatrixRandom(matrix_host, m, n);
    // print the generated vector and matrix
    if (m <= 10 && n <= 10)
    {
        printArray(vector_host, n);
        printMatrix(matrix_host, m, n);
    }
    else
        printf("too many values to print...\n");
 
    // Copy the vector from host to device
    status = cublasSetVector(n, sizeof(double),
                             vector_host, 1, vector_device, 1);
    if (status != CUBLAS_STATUS_SUCCESS)
    {
        printf("cublasSetVector failed\n");
        return 2;
    }
 
    // Copy of matrix from host to device
    status = cublasSetMatrix(m, n, sizeof(double),
                             matrix_host, m, matrix_device, m);
    if (status != CUBLAS_STATUS_SUCCESS)
    {
        printf("cublasSetMatrix failed\n");
        return 2;
    }
 
    // Finally, the matrix-vector multiplication
    status = cublasDgemv(handle, CUBLAS_OP_N, m, n, &alpha,
                         matrix_device, m, vector_device, 1,
                         &beta, vector_final, 1);
    if (status != CUBLAS_STATUS_SUCCESS)
    {
        printf("cublasDgemv failed\n");
        return 2;
    }
    // Copy the result from device to host
    status = cublasGetVector(m, sizeof(double), vector_final,
                             1, final_host, 1);
    if (status != CUBLAS_STATUS_SUCCESS)
    {
        printf("cublasGetVector failed\n");
        return 2;
    }
    // Final output
    printf("matXvet final result\n");
    printArray(final_host, m);
    // Destroy the handle
    cublasDestroy(handle);
    // Free memory
    cudaFree(vector_device);
    cudaFree(vector_final);
    cudaFree(matrix_device);
    free(vector_host);
    free(final_host);
    free(matrix_host);
}

Functions used

void initializeArray(double *array, int n)
{
    int i;
    for (i = 0; i < n; i++)
        array[i] = i;
}
 
void initializeArrayRandom(double *array, int n)
{
    int i;
    double random_number;
    srand((unsigned)time(0));
    for (i = 0; i < n; i++)
    {
        random_number = 200 +
            (double)rand() / ((double)RAND_MAX / (500 - 200));
        array[i] = random_number;
    }
}
 
void initializeMatrix(double *matrix, int M, int N)
{
    int k = 0;
    int i, j;
    // column-wise generation (starting from index 1)
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < M; j++)
        {
            matrix[k++] = j * N + i;
            printf("k=%d value ", k);
            printf("Element = %f\n", (double)(j * N + i));
        }
    }
}
 
void initializeMatrixRandom(double *matrix, int M, int N)
{
    int i, j, k = 0;
    double random_number;
    srand((unsigned)time(0));
    // column-wise generation (starting from index 1)
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < M; j++)
        {
            random_number = 200 +
                (double)rand() / ((double)RAND_MAX / (500 - 200));
            matrix[k++] = random_number;
        }
    }
}
 
void printArray(double *array, int n)
{
    int i;
    printf("vector\n");
    for (i = 0; i < n; i++)
        printf("%f ", array[i]);
    printf("\n\n");
}
 
void printMatrix(double *matrix, int M, int N)
{
    int i, j;
    printf("matrix\n");
    for (i = 0; i < M; i++)
    {
        for (j = 0; j < N; j++)
            printf("%f ", matrix[j * M + i]);
        printf("\n");
    }
}

Obsidian + 🪴 Quartz 4.0

Explorer

HPC - Lecture 14 - CUDA libraries

HPC - Lecture 14 - CUDA libraries

Overview: The cuBLAS Library

cuBLAS introduction

PBLAS (Parallel BLAS)

PLASMA (Parallel Linear Algebra for Scalable Multi-core Architectures)

cuBLAS architecture

Moving data between devices

Matrix management in cuBLAS

cuBLAS handle

Error Management

cuBLAS level 1

dot product between vectors

cuBLAS level 2

Matrix-Vector Multiplication

cuBLAS level 3

Matrix-Matrix Multiplication

Steps to Use cuBLAS

Example of program

Matrix Management Note

Graph View

Table of Contents