【CUDA】矩阵乘向量 matVecMul

Matrix - Vector Multiplication

矩阵-向量乘法是线性代数中的基本操作。它用于将一个矩阵与一个向量相乘。乘法的结果是与输入向量大小相同的向量。

矩阵和向量的乘法如图1所示。

图1

基础kernel与共享内存kernel

执行矩阵-向量乘法的基础kernel是使用单个线程执行输出向量的单个元素的乘法。这意味着所需的线程数等于输出向量中的元素数。线程以一维网格排列，每个线程被分配一个唯一的索引。线程的索引用于访问输入矩阵的对应行。对行和向量进行乘法操作，结果存储在输出向量的对应元素中。

在这种基础kernel中，每个线程将加载整个输入向量。这并不高效，因为输入向量被多次加载。为了避免这种情况，我们可以使用一个tile将输入向量存储在共享内存中。tile是一个大小等于块中线程数的一维数组。向量被加载到tile中，每个线程使用tile执行乘法。

Code

Host代码初始化输入矩阵和向量，随机赋值并调用kernel执行乘法。输入矩阵以线性化格式存储。

#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>#include <cuda_runtime.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include "image2gray.cuh"
#include "helper_cuda.h"
#include "error.cuh"const int FORTIME = 50;int main(int argc, char* argv[])
{int rows, cols, t_x;rows = 4096;cols = 4096;t_x = 1024;thrust::host_vector<float> h_in_mat(rows * cols);thrust::host_vector<float> h_in_vec(cols);srand(time(NULL));for (int i = 0; i < rows * cols; i++) {h_in_mat[i] = rand() / (float)RAND_MAX;if (i < cols)h_in_vec[i] = rand() / (float)RAND_MAX;}thrust::host_vector<float> h_out(rows);for (int i = 0; i < rows; ++i)for (int j = 0; j < cols; ++j)h_out[i] += h_in_mat[i * cols + j] * h_in_vec[j];thrust::device_vector<float> d_in_mat = h_in_mat;thrust::device_vector<float> d_in_vec = h_in_vec;thrust::device_vector<float> d_out(rows);dim3 block(t_x);dim3 grid((rows + block.x - 1) / block.x);cudaEvent_t start, stop;checkCudaErrors(cudaEventCreate(&start));checkCudaErrors(cudaEventCreate(&stop));checkCudaErrors(cudaEventRecord(start));for (int i = 0; i < FORTIME; i++)mat_vec_mul << <grid, block >> > (thrust::raw_pointer_cast(d_out.data()),thrust::raw_pointer_cast(d_in_mat.data()),thrust::raw_pointer_cast(d_in_vec.data()),rows, cols);checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));float elapsed_time;checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));std::cout << "elapsed_time:" << elapsed_time / FORTIME << std::endl;bool success = true;for (int i = 0; i < rows; ++i)if (abs(h_out[i] - d_out[i]) >= 0.001) {std::cout << "Error at " << i << ": " << h_out[i] << " != " << d_out[i] << " (Mat Vec Mul)" << std::endl;success = false;break;}if (success)std::cout << "Success (Mat Vec Mul)" << std::endl;checkCudaErrors(cudaEventRecord(start));for (int i = 0; i < FORTIME; i++)mat_vec_mul_tiles << <grid, block, block.x * sizeof(float) >> > (thrust::raw_pointer_cast(d_out.data()),thrust::raw_pointer_cast(d_in_mat.data()),thrust::raw_pointer_cast(d_in_vec.data()),rows, cols);checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));std::cout << "elapsed_time:" << elapsed_time / FORTIME << std::endl;success = true;for (int i = 0; i < rows; ++i)if (abs(h_out[i] - d_out[i]) >= 0.001) {std::cout << "Error at " << i << ": " << h_out[i] << " != " << d_out[i] << " (Mat Vec Mul Tiles)" << std::endl;success = false;break;}if (success)std::cout << "Success (Mat Vec Mul Tiles)" << std::endl;return 0;
}

Note：

helper_cuda.h 与error.cuh头文件为错误查验工具。后续会发布到Github。

去除checkCudaErrors等错误查验函数不影响程序运行。

下面显示了两个kernel。

template<typename T> __global__
void mat_vec_mul(T *out, T *in_mat, T *in_vec, int m, int n) {int i = blockIdx.x * blockDim.x + threadIdx.x;if(i >= m) return;T temp = 0;for (int k = 0; k < n; ++k)temp += in_mat[i * n + k] * in_vec[k];out[i] = temp;
}template<typename T> __global__
void mat_vec_mul_tiles(T *out, T *in_mat, T *in_vec, int m, int n) {extern __shared__ uint8_t tiles_uint8[];T *tiles = reinterpret_cast<T*>(tiles_uint8);int i = blockIdx.x * blockDim.x + threadIdx.x;T res = 0;int n_phases = (n + blockDim.x - 1) / blockDim.x;for (int phase = 0; phase < n_phases; ++phase){int elem_idx = phase * blockDim.x + threadIdx.x;tiles[threadIdx.x] = elem_idx < n ? in_vec[elem_idx] : 0;__syncthreads();if(i < m)for (int k = 0; k < blockDim.x && phase * blockDim.x + k < n; ++k)res += in_mat[i * n + phase * blockDim.x + k] * tiles[k];__syncthreads();}out[i] = res;
}

基础kernel

kernel首先计算线程的索引。如果索引在输出向量的范围内，则执行乘法。

int i = blockIdx.x * blockDim.x + threadIdx.x;
if(i >= m) return;

将矩阵的每一列与向量相应元素相乘，并将结果存储在输出向量的相应元素中。

T temp = 0;
for (int k = 0; k < n; ++k)temp += in_mat[i * n + k] * in_vec[k];
out[i] = temp;

共享内存kernel

kernel首先计算线程的索引。kernel不会检查边界条件，因为需要加载向量到tile中的超出输出向量范围的线程。

int i = blockIdx.x * blockDim.x + threadIdx.x;

现在，乘法被分成几个阶段。在每个阶段，线程将把向量的单个或零个元素（如果它们超出范围）加载到tile中。然后线程将执行加载到tile中的向量和输入矩阵行的相应元素的乘法。结果存储在输出向量的相应元素中。

在加载向量到tile中和执行乘法时，线程都会执行边界检查，以避免超出范围的访问。

T res = 0;
int n_phases = (n + blockDim.x - 1) / blockDim.x;
for (int phase = 0; phase < n_phases; ++phase){int elem_idx = phase * blockDim.x + threadIdx.x;tiles[threadIdx.x] = elem_idx < n ? in_vec[elem_idx] : 0;__syncthreads();if(i < m)for (int k = 0; k < blockDim.x && phase * blockDim.x + k < n; ++k)res += in_mat[i * n + phase * blockDim.x + k] * tiles[k];__syncthreads();
}
out[i] = res;