Масштабирование строк матрицы с помощью CUDA

Если вы используете BLAS gemv с единичным вектором, результатом будет вектор обратной величины коэффициентов масштабирования (1/альфа), которые вам нужны. Это легкая часть.

Применение факторов по строкам немного сложнее, потому что в стандартном BLAS нет ничего похожего на оператор продукта Адамара, который вы могли бы использовать. Кроме того, поскольку вы упоминаете BLAS, я предполагаю, что вы используете хранилище основного порядка столбцов для своих матриц, что не так просто для операций по строкам. очень медленный способ сделать это — выполнить BLAS scal для каждой строки с шагом, но для этого потребуется один вызов BLAS на строку, и доступ к памяти с шагом убьет производительность. из-за влияния на объединение и когерентность кэша L1.

Мое предложение состояло бы в том, чтобы использовать ваше собственное ядро ​​для второй операции. Это не обязательно должно быть так сложно, возможно, только что-то вроде этого:

template<typename T>
__global__ void rowscale(T * X, const int M, const int N, const int LDA,
                             const T * ralpha)
{
    for(int row=threadIdx.x; row<M; row+=gridDim.x) {
        const T rscale = 1./ralpha[row]; 
        for(int col=blockIdx.x; col<N; col+=blockDim.x) 
            X[row+col*LDA] *= rscale;
    }
}

Это просто набор блоков, которые шагают по строкам по столбцам, масштабируясь по мере их продвижения. Должен работать для любой упорядоченной матрицы столбца любого размера. Доступ к памяти должен быть объединен, но в зависимости от того, насколько вы беспокоитесь о производительности, вы можете попробовать ряд оптимизаций. По крайней мере, это дает общее представление о том, что делать.

В Cublas 5.0 появилась подобная blas процедура под названием cublas(Type)dgmm, которая представляет собой умножение матрицы на диагональную матрицу (представленную вектором).

Существует левый вариант (который будет масштабировать строки) или правый вариант, который будет масштабировать столбец.

Подробную информацию см. в документации CUBLAS 5.0.

Итак, в вашей задаче вам нужно создать вектор, содержащий всю альфу на графическом процессоре, и использовать cublasdgmm с левой опцией.

Я хочу обновить ответы выше примером, рассматривающим использование CUDA Thrust thrust::transform и cuBLAS cublas<t>dgmm. Я пропускаю вычисление коэффициентов масштабирования alpha, так как это уже рассматривалось в

Уменьшение строк матрицы с помощью CUDA

и

Уменьшение столбцов матрицы с помощью CUDA

Ниже приведен полный пример:

#include <thrust/device_vector.h>
#include <thrust/reduce.h>
#include <thrust/random.h>
#include <thrust/sort.h>
#include <thrust/unique.h>
#include <thrust/equal.h>

#include <cublas_v2.h>

#include "Utilities.cuh"
#include "TimingGPU.cuh"

/**************************************************************/
/* CONVERT LINEAR INDEX TO ROW INDEX - NEEDED FOR APPROACH #1 */
/**************************************************************/
template <typename T>
struct linear_index_to_row_index : public thrust::unary_function<T,T> {

    T Ncols; // --- Number of columns

    __host__ __device__ linear_index_to_row_index(T Ncols) : Ncols(Ncols) {}

    __host__ __device__ T operator()(T i) { return i / Ncols; }
};

/***********************/
/* RECIPROCAL OPERATOR */
/***********************/
struct Inv: public thrust::unary_function<float, float>
{
    __host__ __device__ float operator()(float x)
    {
        return 1.0f / x;
    }
};

/********/
/* MAIN */
/********/
int main()
{
    /**************************/
    /* SETTING UP THE PROBLEM */
    /**************************/

    const int Nrows = 10;           // --- Number of rows
    const int Ncols =  3;           // --- Number of columns  

    // --- Random uniform integer distribution between 0 and 100
    thrust::default_random_engine rng;
    thrust::uniform_int_distribution<int> dist1(0, 100);

    // --- Random uniform integer distribution between 1 and 4
    thrust::uniform_int_distribution<int> dist2(1, 4);

    // --- Matrix allocation and initialization
    thrust::device_vector<float> d_matrix(Nrows * Ncols);
    for (size_t i = 0; i < d_matrix.size(); i++) d_matrix[i] = (float)dist1(rng);

    // --- Normalization vector allocation and initialization
    thrust::device_vector<float> d_normalization(Nrows);
    for (size_t i = 0; i < d_normalization.size(); i++) d_normalization[i] = (float)dist2(rng);

    printf("\n\nOriginal matrix\n");
    for(int i = 0; i < Nrows; i++) {
        std::cout << "[ ";
        for(int j = 0; j < Ncols; j++)
            std::cout << d_matrix[i * Ncols + j] << " ";
        std::cout << "]\n";
    }

    printf("\n\nNormlization vector\n");
    for(int i = 0; i < Nrows; i++) std::cout << d_normalization[i] << "\n";

    TimingGPU timerGPU;

    /*********************************/
    /* ROW NORMALIZATION WITH THRUST */
    /*********************************/

    thrust::device_vector<float> d_matrix2(d_matrix);

    timerGPU.StartCounter();
    thrust::transform(d_matrix2.begin(), d_matrix2.end(),
                      thrust::make_permutation_iterator(
                                d_normalization.begin(),
                                thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols))),
                      d_matrix2.begin(),
                      thrust::divides<float>());
    std::cout << "Timing - Thrust = " << timerGPU.GetCounter() << "\n";

    printf("\n\nNormalized matrix - Thrust case\n");
    for(int i = 0; i < Nrows; i++) {
        std::cout << "[ ";
        for(int j = 0; j < Ncols; j++)
            std::cout << d_matrix2[i * Ncols + j] << " ";
        std::cout << "]\n";
    }

    /*********************************/
    /* ROW NORMALIZATION WITH CUBLAS */
    /*********************************/
    d_matrix2 = d_matrix;

    cublasHandle_t handle;
    cublasSafeCall(cublasCreate(&handle));

    timerGPU.StartCounter();
    thrust::transform(d_normalization.begin(), d_normalization.end(), d_normalization.begin(), Inv());
    cublasSafeCall(cublasSdgmm(handle, CUBLAS_SIDE_RIGHT, Ncols, Nrows, thrust::raw_pointer_cast(&d_matrix2[0]), Ncols, 
                   thrust::raw_pointer_cast(&d_normalization[0]), 1, thrust::raw_pointer_cast(&d_matrix2[0]), Ncols));
    std::cout << "Timing - cuBLAS = " << timerGPU.GetCounter() << "\n";

    printf("\n\nNormalized matrix - cuBLAS case\n");
    for(int i = 0; i < Nrows; i++) {
        std::cout << "[ ";
        for(int j = 0; j < Ncols; j++)
            std::cout << d_matrix2[i * Ncols + j] << " ";
        std::cout << "]\n";
    }

    return 0;
}

Файлы Utilities.cu и Utilities.cuh хранятся здесь и опущены здесь. TimingGPU.cu и TimingGPU.cuh поддерживаются здесь и также опускаются.

Я протестировал приведенный выше код на Kepler K20c, и вот результат:

                 Thrust      cuBLAS
2500 x 1250      0.20ms      0.25ms
5000 x 2500      0.77ms      0.83ms

Из времени cuBLAS я исключаю время cublasCreate. Даже с учетом этого версия CUDA Thrust кажется более удобной.