cuda - 减少CUDA中的矩阵行或列

Question

我正在使用 CUDAcuBLAS来执行矩阵运算。

我需要对矩阵的行（或列）求和。目前我通过将矩阵与一个向量相乘来做到这一点，但这似乎并不那么有效。

有没有更好的办法？中找不到任何东西cuBLAS。

Answer 1

实际上，将矩阵与一个向量相乘cublas_gemv()是一种非常有效的方法，除非您正在考虑手动编写自己的内核。

您可以轻松分析cublas_gemv(). 非常接近于一次简单地读取整个矩阵数据，可以看作是矩阵行/列求和的理论峰值性能。

额外的操作“x1.0”不会导致性能大幅下降，因为：

1. cublas_gemv()基本上是内存带宽限制操作，额外的算术指令不会成为瓶颈；
2. FMA指令进一步降低指令吞吐量；
3. 个向量的 mem 通常比矩阵的要小得多，并且可以很容易地被 GPU 缓存以减少到 mem 带宽。

cublas_gemv()还可以帮助您处理矩阵布局问题。它适用于行/列主要和任意填充。

我也问过类似的问题。我的实验表明cublas_gemv()比分段 reduce using 更好Thrust::reduce_by_key，这是矩阵行求和的另一种方法。

Answer 2

与此相关的帖子包含有关同一主题的有用答案，可在

使用 CUDA 减少矩阵行

和

使用 CUDA 减少矩阵列。

在这里，我只想指出如何通过将行乘以同一矩阵来减少矩阵列的方法可以推广到执行向量集合的线性组合。换句话说，如果要计算以下向量基展开

其中f(x_m)是函数的样本f(x)，而\psi_n是基函数，c_n是扩展系数，那么\psi_n可以在N x M矩阵中组织，系数c_n在行向量中，然后使用计算向量 x 矩阵乘法cublas<t>gemv.

下面，我报告了一个完整的示例：

#include <cublas_v2.h>

#include <thrust/device_vector.h>
#include <thrust/random.h>

#include <stdio.h>
#include <iostream>

#include "Utilities.cuh"

/********************************************/
/* LINEAR COMBINATION FUNCTION - FLOAT CASE */
/********************************************/
void linearCombination(const float * __restrict__ d_coeff, const float * __restrict__ d_basis_functions_real, float * __restrict__ d_linear_combination,
                       const int N_basis_functions, const int N_sampling_points, const cublasHandle_t handle) {

    float alpha = 1.f;
    float beta  = 0.f;
    cublasSafeCall(cublasSgemv(handle, CUBLAS_OP_N, N_sampling_points, N_basis_functions, &alpha, d_basis_functions_real, N_sampling_points, 
                               d_coeff, 1, &beta, d_linear_combination, 1));

}

void linearCombination(const double * __restrict__ d_coeff, const double * __restrict__ d_basis_functions_real, double * __restrict__ d_linear_combination,
                       const int N_basis_functions, const int N_sampling_points, const cublasHandle_t handle) {

    double alpha = 1.;
    double beta  = 0.;
    cublasSafeCall(cublasDgemv(handle, CUBLAS_OP_N, N_sampling_points, N_basis_functions, &alpha, d_basis_functions_real, N_sampling_points, 
                               d_coeff, 1, &beta, d_linear_combination, 1));

}

/********/
/* MAIN */
/********/
int main()
{
    const int N_basis_functions = 5;     // --- Number of rows                  -> Number of basis functions
    const int N_sampling_points = 8;     // --- Number of columns               -> Number of sampling points of the basis functions

    // --- Random uniform integer distribution between 10 and 99
    thrust::default_random_engine rng;
    thrust::uniform_int_distribution<int> dist(10, 99);

    // --- Matrix allocation and initialization
    thrust::device_vector<float> d_basis_functions_real(N_basis_functions * N_sampling_points);
    for (size_t i = 0; i < d_basis_functions_real.size(); i++) d_basis_functions_real[i] = (float)dist(rng);

    thrust::device_vector<double> d_basis_functions_double_real(N_basis_functions * N_sampling_points);
    for (size_t i = 0; i < d_basis_functions_double_real.size(); i++) d_basis_functions_double_real[i] = (double)dist(rng);

    /************************************/
    /* COMPUTING THE LINEAR COMBINATION */
    /************************************/
    cublasHandle_t handle;
    cublasSafeCall(cublasCreate(&handle));

    thrust::device_vector<float>  d_linear_combination_real(N_sampling_points);
    thrust::device_vector<double> d_linear_combination_double_real(N_sampling_points);
    thrust::device_vector<float>  d_coeff_real(N_basis_functions, 1.f);
    thrust::device_vector<double> d_coeff_double_real(N_basis_functions, 1.);

    linearCombination(thrust::raw_pointer_cast(d_coeff_real.data()), thrust::raw_pointer_cast(d_basis_functions_real.data()), thrust::raw_pointer_cast(d_linear_combination_real.data()),
                      N_basis_functions, N_sampling_points, handle);
    linearCombination(thrust::raw_pointer_cast(d_coeff_double_real.data()), thrust::raw_pointer_cast(d_basis_functions_double_real.data()), thrust::raw_pointer_cast(d_linear_combination_double_real.data()),
                      N_basis_functions, N_sampling_points, handle);

    /*************************/
    /* DISPLAYING THE RESULT */
    /*************************/
    std::cout << "Real case \n\n";
    for(int j = 0; j < N_sampling_points; j++) {
        std::cout << "Column " << j << " - [ ";
        for(int i = 0; i < N_basis_functions; i++)
            std::cout << d_basis_functions_real[i * N_sampling_points + j] << " ";
        std::cout << "] = " << d_linear_combination_real[j] << "\n";
    }

    std::cout << "\n\nDouble real case \n\n";
    for(int j = 0; j < N_sampling_points; j++) {
        std::cout << "Column " << j << " - [ ";
        for(int i = 0; i < N_basis_functions; i++)
            std::cout << d_basis_functions_double_real[i * N_sampling_points + j] << " ";
        std::cout << "] = " << d_linear_combination_double_real[j] << "\n";
    }

    return 0;
}