CUDA矩阵乘法的优化

本文以 深入浅出谈cuda 中矩阵乘法优化的部分为主线，体会cuda的并行执行过程。

1 实验环境

个人笔记本电脑，显卡NIVIDA GTX950M，操作系统win10，IDE是VS2015，显卡的详细信息如图。

2 实验过程

2.1 初始版本

加载了4000个block，每个block 256个线程来做并行的矩阵乘法运算。

2.2 第一次改良——提高精度

由于GPU上只能用32位浮点数进行运算，大量累加时误差会累积。修改核函数中矩阵相乘计算的部分，利用 Kahan’s Summation Formula 来提高精确度。这样做提升了计算精度，但是没有改变运算时间。

2.3 第二次改良——使用共享内存

加载了1000个blcok，每个block 256个线程，共享内存大小32bit*1000，利用shared memory来实现同一block中各线程的数据共享。每个row的数据由一个block进行计算，最后同步，求和。

2.4 第三次改良——使用cudaMallocPitch()

考虑对齐问题，修改cudamalloc的部分。使用cudaMallocPitch() 申请逻辑二维数组，自动以最佳的倍数来配置内存，并把配置的宽度传回。

2.5 第四次改良——矩阵分块

将1000*1000的矩阵切分成8*8个16*16的矩阵，AB矩阵分块后两两相乘。加载64*64个block，每个block 16*16个线程（二维线程组织），从而将内存存取量降低到原来的1/8。

2.6 第五次改良——去掉if条件判断

在分块的基础上，利用矩阵乘法的特点，在申请完显存后将其置0，这样在核函数中做矩阵乘时就可以拿掉if条件判断语句，减少指令数。

3 遇到的问题和解决

1）每次核函数执行完，将矩阵c拷贝到内存时，经常不能正常获取到c。在cudaMemcpyDeviceToHost后调用cudaGetLastError，发现是从device拷贝到host时失败：CUDA error:unspecified launch failure,

这个错误在2.1和2.2中偶发，之后的改良中就没有再见过了，未解决。

2）在2.3及之后的改良中，发现每一次改进都达不到预期的加速效果。尤其是2.4使用pitch后，与之前几乎没有提升。

我想到使用cuda提供的方法，更细粒度地测量核函数运行时长，发现了关于kernel运行时间的讨论，其中提到了debug和release编译运行时间可能存在区别。

默认情况我们总是debug调试用，release分析性能用的。

于是我在release模式下重新配置了工程并编译运行，改进效果才有了第二节中明显的变化。

4 最后版本代码

#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <device_functions.h>
using namespace std;
#define NUM_THREADS 256 
#define BLOCK_SIZE 16

//初始化cuda
bool InitCUDA();
//生成a,b矩阵
void matgen(float *a, int Ida, int n);
//进行CUDA配置的函数，包括内存分配、传数、调用核函数计算等
clock_t matmultCUDA(const float* a, int lda, const float* b, int ldb, float* c, int ldc, int n);
//进行GPU计算的核函数
__global__ static void matMultCUDA(const float* a, size_t lda, const float* b, size_t ldb, float* c, size_t ldc, int n);
//以 CPU 进行矩阵乘法，用来进行验证答案正确与否
void matmult(const float* a, int Ida, const float* b, int Idb, float* c, int Idc, int n);
//计算两个矩阵的最大相对误差和平均相对误差，并把结果印出来
void compare_mat(const float*a, int Ida, const float* b, int Idb, int n);




int main() {
   
	
	int n = 1000;
	if (!InitCUDA()) return 0;

	//定义矩阵
	float *a, *b, *c, *d;
	a = (float*)malloc(sizeof(float) * n * n)