4

在 GPU 上的某些计算中，我需要缩放矩阵中的行，以便给定行中的所有元素总和为 1。

| 一个1,1一个1,2 ... 一个1,N | | 阿尔法1 *a 1,1阿尔法1 *a 1,2 ... 阿尔法1 *a 1,N |
| 一个2,1一个2,2 ... 一个2,N | => | 阿尔法2 *a 2,1阿尔法2 *a 2,2 ... 阿尔法2 *a 2,N |
| . . | | . . |
| 一个N,1一个N,2 ... 一个N,N | | alpha N *a N,1 alpha N *a N,2 ... alpha N *a N,N |

在哪里

α i = 1.0/(a i,1 + a i,2 + ... + a i,N )

我需要alpha's 的向量和缩放矩阵，我想在尽可能少的 blas 调用中做到这一点。该代码将在 nvidia CUDA 硬件上运行。有谁知道任何聪明的方法来做到这一点？

4

3 回答 3

6

Cublas 5.0 引入了一个类似于 blas 的例程，称为 cublas(Type)dgmm，它是将矩阵乘以对角矩阵（由向量表示）。

有一个左侧选项（将缩放行）或右侧选项将缩放列。

有关详细信息，请参阅 CUBLAS 5.0 文档。

因此，在您的问题中，您需要在 GPU 上创建一个包含所有 alpha 的向量，并将 cublasdgmm 与 left 选项一起使用。

于 2012-09-27T18:55:18.487 回答
2

如果您将 BLASgemv与单位向量一起使用，则结果将是您需要的比例因子倒数 (1/alpha) 的向量。那是容易的部分。

逐行应用因素有点困难，因为标准的 BLAS 没有您可以使用的 Hadamard 乘积运算符。另外因为您提到了 BLAS，我假设您正在为矩阵使用列主要顺序存储，这对于行明智的操作来说并不那么简单。真正缓慢的方法是scal在每一行上使用一个间距进行 BLAS，但这将需要每行一次 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;
    }
}

这只是有一堆块逐列逐行，随着它们的进行而缩放。应该适用于任何大小的列主要有序矩阵。内存访问应该合并，但根据您对性能的担心程度，您可以尝试一些优化。它至少给出了一个大致的想法。

于 2012-02-15T12:42:09.687 回答
2

我想用一个考虑使用 CUDA Thrust'sthrust::transform和 of cuBLAS's的例子来更新上面的答案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 版本似乎更方便。

于 2015-05-18T20:10:14.747 回答