减少CUDA中的矩阵行或列

Question

我正在使用 CUDA 和 cuBLAS 来执行矩阵运算。

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

有没有更好的方法？在cuBLAS 中找不到任何内容。

Answer 1

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

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

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

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

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

我也问过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 可以组织在@ 987654331@ 矩阵和系数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;
}