分析下A,B,C为三个矩阵，A为m*n维，B为n*k维，C为m*k维，用A和B来计算C，计算方法是：C = alpha*A*B + beta*C。它的程序如下：
1 // C = alpha*A*B + beta*C
2 void mySgemm(int m, int n, int k, float alpha, float beta,
3 float *A, float *B, float *C)
4 {
5 int i, j, l;
6 float ab;
7 for(j = 0; j < m; j++)
8 {
9 for(i = 0 ;i < k ;i++)
10 {
11 ab = 0.0f;
12 for(l = 0 ;l < n ;l++)
13 {
14 ab += A[j*n+l] * B[l*k+i];
15 }
16 C[j*k+i] = alpha*ab + beta*C[j*k+i];
17 }
18 }
19 }
这个程序修改自HMPP_Tutorial_Labs_CUDA中的lab0。
C中的每个元素的计算是独立的，完全可以并行化。后面的系列文章将会讲述优化这个程序的整个过程。
一步步做程序优化【2】OpenACC指令
这个写了很长时间了，但是一直没有顾上额。把这个版本稍微修改一下，只需要加上一个指令，我们就可以得到不错的效率奥。
看代码吧：
1 // C = alpha*A*B + beta*C
2 void mySgemm(int m, int n, int k, float alpha, float beta,
3 float *A, float *B, float *C)
4 {
5 int i, j, l;
6 float ab;
7 #pragma acc kernels copy(A[0:m*n],B[0:m*n],C[0:m*n])
8 #pragma acc loop independent
9 for(j = 0; j < m; j++)
10 {
11 #pragma acc loop independent
12 for(i = 0 ;i < k ;i++)
13 {
14 ab = 0.0f;
15 for(l = 0 ;l < n ;l++)
16 {
17 ab += A[j*n+l] * B[l*k+i];
18 }
19 C[j*k+i] = alpha*ab + beta*C[j*k+i];
20 }
21 }
22 }
这样，我们只是加入了几个指导语句，剩下的事是编译器帮我们做的奥，你原先的测试程序并不需要任何改变奥。
我之前讲过HMPP编译器的安装和使用，http://blog.csdn.net/bendanban/article/details/7662583大家可以使用HMPP编译器编译这段代码，在Linux下（安装好CUDA，HMPP之后）我们可以使用一下命令编译：
$hmpp --codelet-required gcc your_program.c
执行一下，你会发现速度相当的快了（你要有支持CUDA的显卡才行奥）
大家可以写一个测试程序来调用这个函数，随便你用什么编译器，只要你可以在你的测试程序里找到本文中提供的程序，你完全可以使用高效的函数奥。
为了得到更高的效率，我修改一下这个代码：
1 // C = alpha*A*B + beta*C
2 void mySgemm(int m, int n, int k, float alpha, float beta,
3 float *A, float *B, float *C)
4 {
5 int i, j, l;
6 float ab;
7 #pragma acc kernels copyin(A[0:m*n],B[0:m*n]) copy(C[0:m*n])
8 #pragma acc loop independent
9 for(j = 0; j < m; j++)
10 {
11 #pragma acc loop independent
12 for(i = 0 ;i < k ;i++)
13 {
14 ab = 0.0f;
15 for(l = 0 ;l < n ;l++)
16 {
17 ab += A[j*n+l] * B[l*k+i];
18 }
19 C[j*k+i] = alpha*ab + beta*C[j*k+i];
20 }
21 }
22 }
这样A和B两个矩阵就可只是传输到GPU上，而C传到GPU，计算结束后会倍传回来。
在copy（）中，A[0:m*n]，表示从第0个元素一共计算m*n个元素，第一个是起始位置，第二个量表示数据长度。
大家把代码拷贝走，去试试吧！！！
一步步做程序优化【3】OpenHMPP指令(更加灵活的使用异构计算)
1、简介下HMPP
HMPP指的是（Hybrid Multicore Parallel Programming），他是由CAPS(http://www.caps-entreprise.com(英文)；www.caps-entreprise.com.cn(中文)) 发起的一种异构计算的标准，他的出现可以大大减少我们的程序优化时间。大家可以参考我之前的几篇讲解HMPP的文章去获得HMPP的试用版。
HMPP是一种基于编译指导语句（类似与OpenMP）的标准，它与OpenMP的区别是：OMP是基于CPU的并行标准，HMPP是基于异构平台的标准（例如CPU+GPU，CPU+MIC），它支持C和Fortran两种语言。另外HMPP编译器可以根据你的#pragma指令产生CUDA代码，也可以直接编译CUDA代码！
总之，HMPP编译器非常强大！
2、使用HMPP以及OpenACC的一个推荐原则。
使用HMPP是为了尽可能不改变原有代码的基础上只需要添加少量的#pragma 语句就可一获得几十甚至几千倍的加速比。当然前提是你原有的代码要可以正确的按照算法设计的目的执行才行。
3、继续优化矩阵相乘的那段代码
1）重新贴一边需要优化的代码：（特别注意这段代码来值CAPS，这是原始代码，我没有做实质性的修改）
1 /*
2 * Copyright 2008 - 2012 CAPS entreprise. All rights reserved.
3 */
4
5
6 #include <getopt.h>
7 #include <sys/time.h>
8 #include <stdlib.h>
9 #include <stdio.h>
10 #include <string.h>
11 #include <math.h>
12
13 // Number of execution
14 #define NB_RUNS 5
15
16 // Size of the matrix
17 #define SIZE 256
18
19 // Initialization random value
20 #define SRAND_VALUE 5347
21
22 // Use to initialize the matrix
23 float randFloat(float low, float high)
24 {
25 float t = (float)rand() / (float)RAND_MAX;
26 return (1.0f - t) * low + t * high;
27 }
28
29 ////////////////////////////////////////////////////////////////////////////////
30 // sgemm_codelet
31 ////////////////////////////////////////////////////////////////////////////////
32 void mySgemm( int m, int n, int k, float alpha, float beta,
33 float a[m][n], float b[n][k], float c[m][k] )
34 {
35 int i,j,l; // Induction variables
36 float ab; // Temporary result
37
38 for( j = 0 ; j < m ; j++ ) {
39 for( i = 0 ; i < k ; i++ ) {
40 ab=0.0f;
41 for( l = 0 ; l < n ; l++ ){
42 ab += a[j][l] * b[l][i];
43 }
44 c[j][i] = alpha * ab + beta * c[j][i];
45 }
46 }
47 }
48
49
50 ////////////////////////////////////////////////////////////////////////////////
51 // Main program
52 ////////////////////////////////////////////////////////////////////////////////
53 int main(int argc, char **argv)
54 {
55
56 int m=SIZE, n=SIZE, k = SIZE;
57
58 float *my_a=NULL, *b=NULL, *c_hwa=NULL, *c_cpu=NULL;
59 int i, j, ii;
60 // For timer measures
61 struct timeval tv_global_begin, tv_global_end; // global timer (all iterations)
62 struct timeval tv_begin, tv_end; // local timer (1 iteration)
63
64 unsigned long long int best_measure_GPU = 0;
65 unsigned long long int sum_measure_GPU = 0;
66
67 unsigned long long int best_measure_CPU = 0;
68 unsigned long long int sum_measure_CPU = 0;
69
70 unsigned long long int global_CPU_time = 0;
71 unsigned long long int global_GPU_time = 0;
72
73 unsigned long long int current;
74
75 float alpha, beta;
76
77 double error = 0.0;
78 int index_i = 0.0;
79 int index_j = 0.0;
80 double valueCPU = 0.0;
81 double valueGPU = 0.0;
82
83
84
85 // Allocating CPU memory
86 my_a = (float *)malloc(m* n * sizeof(float));
87 my_b = (float *)malloc(n * k * sizeof(float));
88 c_hwa = (float *)malloc(m * k * sizeof(float));
89 c_cpu = (float *)malloc(m * k * sizeof(float));
90
91 if((my_a == NULL) || (my_b == NULL) || (c_hwa == NULL) || (c_cpu == NULL)) {
92 fprintf( stderr, "
**** error : memory allocation failed ****
");
93 return 1;
94 }
95
96 fprintf( stdout, "---- Initialization of the Matrices ----
");
97 srand(SRAND_VALUE);
98
99 //Generate options set
100
101 for(i = 0; i < m; i++){
102 for(j = 0; j < n; j++){
103 my_a[i*n+j] = randFloat(0.0001f, 1.0f);
104 }
105 }
106
107
108 for(i = 0; i < n; i++){
109 for(j = 0; j < k; j++){
110 my_b[i*k+j] = randFloat(0.0001f, 1.0f);
111 }
112 }
113
114
115 for(i = 0; i < m; i++){
116 for(j = 0; j < k; j++) {
117 c_cpu[i*k+j] = randFloat(1.0, 20.0f);
118 c_hwa[i*k+j] = c_cpu[i*k+j];
119 }
120 }
121
122 alpha = 0.5;
123 beta = randFloat(1.0, 2.0);
124
125 fprintf( stdout, "---- Running calculations ----
");
126
127
128 // run sgemm on GPU (NB_RUNS iterations)
129 printf("Run on GPU
");
130
131 // Start timer
132 gettimeofday(&tv_global_begin, NULL);
133
134
135 for( i=0; i<NB_RUNS; i++ ) {
136 printf("%d ",i);
137 gettimeofday(&tv_begin, NULL);
138
139 mySgemm( m, n, k, alpha, beta, my_a, my_b, c_hwa );
140 gettimeofday(&tv_end, NULL);
141
142 current = (tv_end.tv_sec-tv_begin.tv_sec)*1e6 + tv_end.tv_usec-tv_begin.tv_usec;
143
144 if( ( best_measure_GPU == 0 ) || ( best_measure_GPU > current ) ){
145 best_measure_GPU = current;
146 }
147 sum_measure_GPU += current;
148 }
149
150
151
152 gettimeofday(&tv_global_end, NULL);
153 global_GPU_time = (tv_global_end.tv_sec-tv_global_begin.tv_sec)*1e6 + tv_global_end.tv_usec-tv_global_begin.tv_usec;
154 // run sgemm on CPU (NB_RUNS iterations)
155 printf("
Run on CPU
");
156
157 // Start timer
158 gettimeofday(&tv_global_begin, NULL);
159
160 for( i=0; i<NB_RUNS; i++ ) {
161 printf("%d ",i);
162 gettimeofday(&tv_begin, NULL);
163
164 mySgemm( m, n, k, alpha, beta, my_a, my_b, c_cpu );
165
166 gettimeofday(&tv_end, NULL);
167 current = (tv_end.tv_sec-tv_begin.tv_sec)*1e6 + tv_end.tv_usec-tv_begin.tv_usec;
168
169 if( ( best_measure_CPU == 0 ) || ( best_measure_CPU > current ) ){
170 best_measure_CPU = current;
171 }
172 sum_measure_CPU += current;
173 }
174
175 gettimeofday(&tv_global_end, NULL);
176 global_CPU_time = (tv_global_end.tv_sec-tv_global_begin.tv_sec)*1e6 + tv_global_end.tv_usec-tv_global_begin.tv_usec;
177
178
179 // Compute error between GPU and CPU
180 for( ii = 0; ii < m; ii++){
181 for(j = 0; j < k; j++){
182 double lerror = fabs((c_hwa[ii*k+j]-c_cpu[ii*k+j])/c_cpu[ii*k+j]);
183 if (lerror > error) {
184 error = lerror;
185 valueCPU = c_cpu[ii*k+j];
186 valueGPU = c_hwa[ii*k+j];
187 index_i = ii;
188 index_j = j;
189 }
190 }
191 }
192
193 if (error > 2e-06) {
194 fprintf( stdout, "
The error is %e with index %d:%d @ %e (CPU) / %e (GPU)
", error, index_i, index_j, valueCPU, valueGPU);
195 fprintf( stdout, "The error is is too big!
");
196 return -1;
197 }
198
199
200 fprintf( stdout, "
---- Results ----");
201 fprintf( stdout, "
");
202 fprintf( stdout, "Sizes of matrices: M:%i N:%i K:%i
", m, n, k);
203 fprintf( stdout, "Best HWA time : %f ms
", best_measure_GPU / 1e3 );
204 fprintf( stdout, "Mean HWA time : %f ms
", sum_measure_GPU / NB_RUNS / 1e3);
205 fprintf( stdout, "
");
206 fprintf( stdout, "Best CPU time : %f ms
", best_measure_CPU / 1e3 );
207 fprintf( stdout, "Mean CPU time : %f ms
", sum_measure_CPU / NB_RUNS / 1e3);
208 fprintf( stdout, "
");
209 fprintf( stdout, "Global HWA time : %f ms
", global_GPU_time / 1e3 );
210 fprintf( stdout, "Global CPU time : %f ms
", global_CPU_time / 1e3 );
211 fprintf( stdout, "
");
212 fprintf( stdout, "Speed-up : %f (computed on the best time)",
213 ((float)best_measure_CPU)/best_measure_GPU);
214
215
216 fprintf( stdout, "
");
217
218 free(my_a);
219 free(my_b);
220 free(c_hwa);
221 free(c_cpu);
222
223 return 0;
224 }
注意上述代码中，测试了两次统一个函数的执行结果，下面加入两句简单的指令，然后编译执行下，看一下加速比情况。
在第31与第32行插入一下语句：
#pragma hmpp mylab codelet, target=CUDA, args[*].transfer=atcall
在138行插入：
#pragma hmpp mylab callsite
编译执行：
[]$hmpp --codelet-required gcc source.c
执行结果：
---- Initialization of the Matrices ----
---- Running calculations ---- Run on GPU 0 1 2 3 4 Run on CPU 0 1 2 3 4 ---- Results ---- Sizes of matrices: M:256 N:256 K:256 Best HWA time : 1.436000 ms Mean HWA time : 21.837000 ms Best CPU time : 86.995000 ms Mean CPU time : 87.583000 ms Global HWA time : 109.192000 ms Global CPU time : 437.922000 ms Speed-up : 60.581478 (computed on the best time)
加速比是60倍多！我只是键入了两行指令而已！！
当然HMPP并没有到这里这么简单，它提供了很多指令，指令学习并不难，也就是说我们不用直接学习CUDA或者OpenCL就可以很方便的使用GPU的计算资源了。种种好处 只有在你试用之后才能知道的奥。
后面的博客我还会讲解更多的指令，还有一些有意思的细节。欢迎大家关注奥！