减少CUDA中的矩阵行或列

Question

我正在使用CUDA和cuBLAS来执行矩阵操作。

我需要总结矩阵的行（或列）。目前我正在通过将矩阵乘以一个矢量来实现，但这看起来效率不高。

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

谢谢。

Answer 1

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

您可以轻松地分析cublas_gemv()的存储带宽。这与仅仅读取整个矩阵数据一次非常接近，这可以看作矩阵行/列求和的理论峰值性能。

额外的操作“×1.0”不会导致太大的业绩减少的原因是：

1. cublas_gemv()基本上是一个MEM带宽的限制操作，额外的运算指令不会成为瓶颈;
2. FMA指令进一步降低指令吞吐量;
3. mem的一个向量通常比矩阵的要小得多，并且可以通过GPU轻松缓存以降低到mem带宽。

cublas_gemv()还可以帮助您处理矩阵布局问题。它适用于row/col-major和任意填充。

我也问过a similar question这个。我的实验显示cublas_gemv()比使用Thrust::reduce_by_key的分段缩减好，这是另一种矩阵行总和的方法。

Answer 2

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

Reduce matrix rows with CUDA

和

Reduce matrix columns with CUDA。

这里我只想指出如何通过乘以同一矩阵的行来减少矩阵的列的方法可以推广到执行矢量集合的线性组合。换句话说，如果一个人想计算以下矢量基础膨胀

其中f(x_m)是函数f(x)的样品，而\psi_n的是基函数和c_n的是膨胀系数，那么可以将\psi_n的组织成一个N x M矩阵并将系数c_n排成一行，然后使用cublas<t>gemv计算向量x矩阵乘法。

下面，我要检举完全样例：

#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; 
}