矩阵乘法

▶ 计算矩阵矩阵乘法 A_m×n B_n×p == C_m×p 过程。

▶ 原始矩阵乘法，一个线程计算结果矩阵中的一个元素。

  1 #include <stdio.h>
  2 #include <stdlib.h>
  3 #include <malloc.h>
  4 #include <time.h>
  5 #include "cuda_runtime.h"
  6 #include "device_launch_parameters.h"
  7 
  8 #define M           (1024 + 13)
  9 #define N           (1024 + 31)
 10 #define P           (1024 + 7)
 11 #define THREAD      16
 12 #define SEED        (unsigned int)clock()
 13 #define CEIL(x,y)   (( x - 1 ) /  y + 1)
 14 
 15 void checkNULL(void *input)
 16 {
 17     if (input == NULL)
 18     {
 19         printf("
	find a NULL!");
 20         exit(1);
 21     }
 22     return;
 23 }
 24 
 25 void checkCudaError(cudaError input)
 26 {
 27     if (input != cudaSuccess)
 28     {
 29         printf("
	find a cudaError!");
 30         exit(1);
 31     }
 32     return;
 33 }
 34 
 35 __global__ void times(float *a, float *b, float *c)// A(mn)×B(np)=C(np)
 36 {
 37     int k;
 38     float temp;
 39     int idx = blockIdx.x * blockDim.x + threadIdx.x;// B、C中精确列号
 40     int idy = blockIdx.y * blockDim.y + threadIdx.y;// A、C中精确行号
 41 
 42     for (; idy < M; idy += gridDim.y*blockDim.y)
 43     {
 44         for (; idx < P; idx += gridDim.x*blockDim.x)
 45         {
 46             for (k = 0, temp = 0.0f; k < N; k++)
 47                 temp += a[idy * N + k] * b[k * P + idx];
 48             c[idy * P + idx] = temp;
 49         }
 50     }
 51     return;
 52 }
 53 
 54 int main()
 55 {
 56     int i, j, k, flag;
 57     float *a, *b, *c, *dev_a, *dev_b, *dev_c, temp, elapsedTime;
 58     clock_t time, cputime;
 59     cudaEvent_t start, stop;
 60     cudaEventCreate(&start);
 61     cudaEventCreate(&stop);
 62 
 63     printf("
	Matrix size: A(%d, %d) × B(%d, %d) == C(%d, %d)", M, N, N, P, M, P);
 64     printf("
	Start!");
 65 
 66     a = (float *)malloc(sizeof(float) * M * N);
 67     b = (float *)malloc(sizeof(float) * N * P);
 68     c = (float *)malloc(sizeof(float) * M * P);
 69     checkNULL(a);
 70     checkNULL(b);
 71     checkNULL(c);
 72     checkCudaError(cudaMalloc((float **)&dev_a, sizeof(float) * M * N));
 73     checkCudaError(cudaMalloc((float **)&dev_b, sizeof(float) * N * P));
 74     checkCudaError(cudaMalloc((float **)&dev_c, sizeof(float) * M * P));
 75 
 76     time = clock();
 77     srand(SEED);
 78     for (i = 0; i < M * N; a[i] = (float)rand() / RAND_MAX, i++);
 79     for (i = 0; i < N * P; b[i] = (float)rand() / RAND_MAX, i++);
 80     printf("
	Initialized! Time:	%6.2f ms", (float)(clock() - time));
 81 
 82     // 用CPU计算
 83     time = clock();
 84     for (i = 0; i < M; i++)
 85     {
 86         for (j = 0; j < P; j++)
 87         {
 88             for (k = 0, temp = 0.0f; k < N; k++)
 89                 c[i * P + j] += a[i * N + k] * b[k * P + j];
 90         }
 91     }
 92     cputime = clock() - time;
 93     printf("
	CPU cost time:	%6.2f ms", (float)cputime);
 94 
 95     // 用GPU计算
 96     time = clock();
 97     cudaEventRecord(start, 0);
 98 
 99     cudaMemcpy(dev_a, a, sizeof(float) * M * N, cudaMemcpyHostToDevice);
100     cudaMemcpy(dev_b, b, sizeof(float) * N * P, cudaMemcpyHostToDevice);
101 
102     times << < dim3(CEIL(P, THREAD), CEIL(M, THREAD), 1), dim3(THREAD, THREAD) >> > (dev_a, dev_b, dev_c);
103 
104     cudaMemcpy(c, dev_c, sizeof(float) * M * P, cudaMemcpyDeviceToHost);
105     cudaDeviceSynchronize();
106 
107     time = clock() - time;
108     cudaEventRecord(stop, 0);
109 
110     cudaEventSynchronize(stop);
111     cudaEventElapsedTime(&elapsedTime, start, stop);
112     printf("
	GPU cost time by CPU API:	%6.2f ms
	GPU cost time by GPU API:	%6.2f ms

", (float)time, elapsedTime);
113     printf("
	AccRatio =	%3.1f

", (float)cputime / time);
114 
115     // 检验结果（第一次CPU计算的结果被覆盖掉了，这里相当于再算一次，逐个检查）
116     for (i = 0, flag = 1; i < M; i++)
117     {
118         for (j = 0; j < P; j++)
119         {
120             for (k = 0, temp = 0.0f; k < N; k++)
121                 temp += a[i * N + k] * b[k * P + j];
122             if (temp != c[i * P + j])
123             {
124                 printf("
	Calculate error at (%d,%d),
		c[i] == %f,CPU == %f", i, j, c[i * P + j], temp);
125                 flag = 0;
126             }
127         }
128     }
129     if (flag == 1)
130         printf("
	Calculate correctly!");
131 
132     cudaEventDestroy(start);
133     cudaEventDestroy(stop);
134     cudaFree(dev_a);
135     cudaFree(dev_b);
136     cudaFree(dev_c);
137     getchar();
138     return 0;
139 }

● 输出结果

        Matrix size: A(1037, 1055) × B(1055, 1031) == C(1037, 1031)
        Start!
        Initialized! Time:      127.00 ms
        CPU cost time:  4049.00 ms
        GPU cost time by CPU API:       103.00 ms
        GPU cost time by GPU API:       103.02 ms


        AccRatio =      39.3


        Calculate correctly!

● 矩阵初始化消耗的时间，在计算超大型矩阵的过程中初始化需要占用时间（因为矩阵元素是在CPU中逐个初始化的）。

● 比较了GPU运算过程的CPU和GPU计时的结果，在如上代码中（在内存拷贝前开始计时，在内存回拷完成后结束计时），两种计时方式差距很小，基本上认为结果相同。

● 下方为相同的矩阵乘法在CPU中的计算耗时，可见目前GPU效率约为CPU的39倍。

● 若去除内存拷贝时间，则GPU耗时可减少到100ms左右，注意此时若想利用利用CPU计时，则必须在第二次掐表以前加入 cudaDeviceSynchronize(); 。因为核函数调用的同时，主机代码会继续往下执行，只到同步语句或者亟需GPU预算结果的语句才会停止，在上面的例子中，若不加入同步语句，则计时结果总是为0 ms。

● 可以进一步细致调节线程粒度 THREAD，耗时在小范围内有变化。

▶ 改进一，使用共享内存减少全局内存访问量

 1 __global__ void times(float *a, float *b, float *c)// A(mn)×B(np)=C(mp)
 2 {
 3     __shared__ float shareA[THREAD * THREAD];
 4     __shared__ float shareB[THREAD * THREAD];
 5 
 6     int k, t;
 7     float temp;
 8     int idx = blockIdx.x * blockDim.x + threadIdx.x;// B、C中精确列号
 9     int idy = blockIdx.y * blockDim.y + threadIdx.y;// A、C中精确行号  
10 
11     for (k = 0, temp = 0.0f; k < CEIL(N, THREAD); k++)// N为行程长，k为行程步数
12     {
13         shareA[threadIdx.y * THREAD + threadIdx.x] = (idy < M && k * THREAD + threadIdx.x < N) ? a[idy * N + k * THREAD + threadIdx.x] : 0.0f;
14         shareB[threadIdx.y * THREAD + threadIdx.x] = (idx < P && k * THREAD + threadIdx.y < N) ? b[(k * THREAD + threadIdx.y) * P + idx] : 0.0f;
15         __syncthreads();
16 
17         for (t = 0; t < THREAD; t++)// THREAD为行程长，t为行程步数
18             temp += shareA[threadIdx.y * THREAD + t] * shareB[t * THREAD + threadIdx.x];
19         __syncthreads();
20     }
21     if (idx < P && idy < M && idy * P + idx < M * P)
22         c[idy * P + idx] = temp;
23     return;
24 }

● 输出结果

        Matrix size: A(1037, 1055) × B(1055, 1031) == C(1037, 1031)
        Start!
        Initialized! Time:      123.00 ms
        CPU cost time:  4044.00 ms
        GPU cost time by CPU API:       147.00 ms
        GPU cost time by GPU API:       146.75 ms


        AccRatio =      27.5


        Calculate correctly!

▶改进二，共享内存存储方式调整

 1 __global__ void times(float *a, float *b, float *c)// A(mn)×B(np)=C(mp)
 2 {
 3     __shared__ float shareA[THREAD * THREAD];
 4     __shared__ float shareB[THREAD * THREAD];
 5 
 6     int k, t;
 7     float temp;
 8     int idx = blockIdx.x * blockDim.x + threadIdx.x;// B、C中精确列号
 9     int idy = blockIdx.y * blockDim.y + threadIdx.y;// A、C中精确行号  
10 
11     for (k = 0, temp = 0.0f; k < CEIL(N, THREAD); k++)// N为行程长，k为行程步数
12     {
13         // shareA 改为转置（原始方法读取 a，写入 shareMemory 时挑适当位置
14         shareA[threadIdx.x * THREAD + threadIdx.y] = (idy < M && k * THREAD + threadIdx.x < N) ? a[idy * N + k * THREAD + threadIdx.x] : 0.0f;
15         shareB[threadIdx.y * THREAD + threadIdx.x] = (idx < P && k * THREAD + threadIdx.y < N) ? b[(k * THREAD + threadIdx.y) * P + idx] : 0.0f;
16         __syncthreads();
17 
18         for (t = 0; t < THREAD; t++)// THREAD为行程长，t为行程步数
19             // 注意A的行列发生了交换
20             temp += shareA[t * THREAD + threadIdx.y] * shareB[t * THREAD + threadIdx.x];
21         __syncthreads();
22     }
23     if (idx < P && idy < M && idy * P + idx < M * P)
24         c[idy * P + idx] = temp;
25     return;
26 }

● 输出结果

        Matrix size: A(1037, 1055) × B(1055, 1031) == C(1037, 1031)
        Start!
        Initialized! Time:      121.00 ms
        CPU cost time:  4059.00 ms
        GPU cost time by CPU API:       135.00 ms
        GPU cost time by GPU API:       134.80 ms


        AccRatio =      30.1


        Calculate correctly!