cuBLAS и cuFFT: готовые GPU-библиотеки вместо своих ядер

Введение

Для широкого класса задач — линейной алгебры (умножение матриц, решение систем уравнений) и преобразования Фурье — NVIDIA предоставляет высокооптимизированные библиотеки cuBLAS и cuFFT, реализующие эти операции на GPU значительно эффективнее, чем большинство самописных ядер, особенно с учётом многих лет оптимизации под конкретные архитектуры GPU. В статье покажем, как подключить и использовать обе библиотеки на практических примерах.

Концепция

cuBLAS — GPU-реализация BLAS (Basic Linear Algebra Subprograms), стандартного интерфейса для операций над векторами и матрицами (умножение матриц, скалярное произведение, норма вектора и так далее). Использование начинается с создания контекста (cublasHandle_t) через cublasCreate(), который передаётся в каждый вызов операции (например, cublasSgemm для умножения матриц одинарной точности).

cuFFT — GPU-реализация быстрого преобразования Фурье, используется для анализа сигналов, обработки изображений (например, свёртки через теорему о свёртке), научных вычислений. Аналогично cuBLAS, работа начинается с создания «плана» преобразования (cufftPlan1d/cufftPlan2d и так далее), который описывает размер и тип преобразования, а затем выполняется через cufftExecC2C или аналогичные функции в зависимости от типа данных.

Пример кода

#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <iostream>
#include <vector>

void matrixMultiplyWithCublas()
{
    const int n = 4; // размер квадратных матриц n x n

    std::vector<float> hostA(n * n, 1.0f);
    std::vector<float> hostB(n * n, 2.0f);
    std::vector<float> hostC(n * n, 0.0f);

    float *deviceA, *deviceB, *deviceC;
    cudaMalloc(&deviceA, n * n * sizeof(float));
    cudaMalloc(&deviceB, n * n * sizeof(float));
    cudaMalloc(&deviceC, n * n * sizeof(float));

    cudaMemcpy(deviceA, hostA.data(), n * n * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(deviceB, hostB.data(), n * n * sizeof(float), cudaMemcpyHostToDevice);

    cublasHandle_t handle;
    cublasCreate(&handle);

    const float alpha = 1.0f, beta = 0.0f;
    // C = alpha * A * B + beta * C  (стандартная сигнатура GEMM в BLAS)
    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                n, n, n,
                &alpha,
                deviceA, n,
                deviceB, n,
                &beta,
                deviceC, n);

    cudaMemcpy(hostC.data(), deviceC, n * n * sizeof(float), cudaMemcpyDeviceToHost);

    std::cout << "Результат умножения матриц (cuBLAS): C[0] = " << hostC[0] << std::endl;

    cublasDestroy(handle);
    cudaFree(deviceA);
    cudaFree(deviceB);
    cudaFree(deviceC);
}

#include <cufft.h>
#include <cuda_runtime.h>
#include <iostream>
#include <vector>

void performFftWithCufft()
{
    const int signalSize = 256;

    std::vector<cufftComplex> hostSignal(signalSize);
    for (int i = 0; i < signalSize; ++i) {
        hostSignal[i].x = static_cast<float>(i % 16); // вещественная часть
        hostSignal[i].y = 0.0f;                        // мнимая часть
    }

    cufftComplex *deviceSignal;
    cudaMalloc(&deviceSignal, signalSize * sizeof(cufftComplex));
    cudaMemcpy(deviceSignal, hostSignal.data(),
               signalSize * sizeof(cufftComplex), cudaMemcpyHostToDevice);

    cufftHandle plan;
    cufftPlan1d(&plan, signalSize, CUFFT_C2C, 1);
    cufftExecC2C(plan, deviceSignal, deviceSignal, CUFFT_FORWARD);
    cudaDeviceSynchronize();

    cudaMemcpy(hostSignal.data(), deviceSignal,
               signalSize * sizeof(cufftComplex), cudaMemcpyDeviceToHost);

    std::cout << "Первая гармоника после FFT: " << hostSignal[0].x
              << " + " << hostSignal[0].y << "i" << std::endl;

    cufftDestroy(plan);
    cudaFree(deviceSignal);
}

Пояснения к коду

В примере с cuBLAS cublasSgemm выполняет умножение матриц одинарной точности (S в названии означает single precision) по стандартной формуле GEMM из BLAS: C = alpha * op(A) * op(B) + beta * C, где CUBLAS_OP_N указывает, что матрицы используются без транспонирования. Все операции работают над данными, уже находящимися на GPU (deviceA, deviceB, deviceC) — сама библиотека не выполняет копирование память, об этом заботится вызывающий код.

В примере с cuFFT создаётся «план» одномерного комплексного преобразования Фурье через cufftPlan1d, который затем выполняется через cufftExecC2C с указанием направления CUFFT_FORWARD (прямое преобразование, в противоположность CUFFT_INVERSE). План создаётся один раз и может быть переиспользован для множества преобразований одного и того же размера — это важно, поскольку создание плана включает предварительный анализ оптимальной стратегии вычисления, что не бесплатно по времени.

Подводные камни

Column-major порядок хранения матриц в cuBLAS (унаследован от оригинального Fortran BLAS), в отличие от привычного row-major в C/C++ — при работе с матрицами, хранящимися в обычном C-порядке, нужно либо транспонировать логику вызова (передавая операции в «обратном» порядке через трюк с переключением A и B местами), либо явно транспонировать данные, что легко перепутать и получить неверный результат без явной ошибки выполнения.
Повторное создание cufftHandle/cublasHandle_t на каждый вызов вместо переиспользования — создание контекста и плана не бесплатно, и при частых вызовах операций в цикле обработки (например, потоковой обработки сигнала) стоит создавать handle/plan один раз заранее.
Точность вычислений (single vs double precision). Использование cublasSgemm (float) там, где требуется точность cublasDgemm (double), может приводить к численным ошибкам в задачах с большим динамическим диапазоном значений или длинными цепочками вычислений — выбор точности должен быть осознанным, а не «по умолчанию».
Отсутствие проверки кодов ошибок специфичных для библиотек. cublasStatus_t и cufftResult — отдельные типы ошибок, отличные от cudaError_t, и их также нужно явно проверять (аналогично паттерну CUDA_CHECK из статьи про обработку ошибок) — игнорирование ошибок инициализации handle или плана может привести к работе с невалидным контекстом, дающим непредсказуемые результаты.