一步步做程序优化【2】OpenACC指令
这个写了很长时间了，但是一直没有顾上额。把这个版本稍微修改一下，只需要加上一个指令，我们就可以得到不错的效率奥。
看代码吧：
// C = alpha*A*B + beta*C
void mySgemm(int m, int n, int k, float alpha, float beta,\
float *A, float *B, float *C)
{
int i, j, l;
float ab;
#pragma acc kernels copy(A[0:m*n],B[0:m*n],C[0:m*n])
#pragma acc loop independent
for(j = 0; j < m; j++)
{
#pragma acc loop independent
for(i = 0 ;i < k ;i++)
{
ab = 0.0f;
for(l = 0 ;l < n ;l++)
{
ab += A[j*n+l] * B[l*k+i];
}
C[j*k+i] = alpha*ab + beta*C[j*k+i];
}
}
}
这样，我们只是加入了几个指导语句，剩下的事是编译器帮我们做的奥，你原先的测试程序并不需要任何改变奥。
我之前讲过HMPP编译器的安装和使用，http://blog.csdn.net/bendanban/article/details/7662583大家可以使用HMPP编译器编译这段代码，在Linux下（安装好CUDA，HMPP之后）我们可以使用一下命令编译：
$hmpp --codelet-required gcc your_program.c
执行一下，你会发现速度相当的快了（你要有支持CUDA的显卡才行奥）
大家可以写一个测试程序来调用这个函数，随便你用什么编译器，只要你可以在你的测试程序里找到本文中提供的程序，你完全可以使用高效的函数奥。
为了得到更高的效率，我修改一下这个代码：
// C = alpha*A*B + beta*C
void mySgemm(int m, int n, int k, float alpha, float beta,\
float *A, float *B, float *C)
{
int i, j, l;
float ab;
#pragma acc kernels copyin(A[0:m*n],B[0:m*n]) copy(C[0:m*n])
#pragma acc loop independent
for(j = 0; j < m; j++)
{
#pragma acc loop independent
for(i = 0 ;i < k ;i++)
{
ab = 0.0f;
for(l = 0 ;l < n ;l++)
{
ab += A[j*n+l] * B[l*k+i];
}
C[j*k+i] = alpha*ab + beta*C[j*k+i];
}
}
}
这样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，这是原始代码，我没有做实质性的修改）
/**//*
* Copyright 2008 - 2012 CAPS entreprise. All rights reserved.
*/
#include <getopt.h>
#include <sys/time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
// Number of execution
#define NB_RUNS 5
// Size of the matrix
#define SIZE 256
// Initialization random value
#define SRAND_VALUE 5347
// Use to initialize the matrix
float randFloat(float low, float high)
{
float t = (float)rand() / (float)RAND_MAX;
return (1.0f - t) * low + t * high;
}
/**/////////////////////////////////////////////////////////////////////////////////// sgemm_codelet
/**/////////////////////////////////////////////////////////////////////////////////void mySgemm( int m, int n, int k, float alpha, float beta,
float a[m][n], float b[n][k], float c[m][k] )
{
int i,j,l; // Induction variables
float ab; // Temporary result
for( j = 0 ; j < m ; j++ ) {
for( i = 0 ; i < k ; i++ ) {
ab=0.0f;
for( l = 0 ; l < n ; l++ ){
ab += a[j][l] * b[l][i];
}
c[j][i] = alpha * ab + beta * c[j][i];
}
}
}
/**/////////////////////////////////////////////////////////////////////////////////// Main program
/**/////////////////////////////////////////////////////////////////////////////////int main(int argc, char **argv)
{
int m=SIZE, n=SIZE, k = SIZE;
float *my_a=NULL, *b=NULL, *c_hwa=NULL, *c_cpu=NULL;
int i, j, ii;
// For timer measures
struct timeval tv_global_begin, tv_global_end; // global timer (all iterations)
struct timeval tv_begin, tv_end; // local timer (1 iteration)
unsigned long long int best_measure_GPU = 0;
unsigned long long int sum_measure_GPU = 0;
unsigned long long int best_measure_CPU = 0;
unsigned long long int sum_measure_CPU = 0;
unsigned long long int global_CPU_time = 0;
unsigned long long int global_GPU_time = 0;
unsigned long long int current;
float alpha, beta;
double error = 0.0;
int index_i = 0.0;
int index_j = 0.0;
double valueCPU = 0.0;
double valueGPU = 0.0;
// Allocating CPU memory
my_a = (float *)malloc(m* n * sizeof(float));
my_b = (float *)malloc(n * k * sizeof(float));
c_hwa = (float *)malloc(m * k * sizeof(float));
c_cpu = (float *)malloc(m * k * sizeof(float));
if((my_a == NULL) || (my_b == NULL) || (c_hwa == NULL) || (c_cpu == NULL)) {
fprintf( stderr, "\n**** error : memory allocation failed ****\n\n");
return 1;
}
fprintf( stdout, "---- Initialization of the Matrices ----\n\n");
srand(SRAND_VALUE);
//Generate options set
for(i = 0; i < m; i++){
for(j = 0; j < n; j++){
my_a[i*n+j] = randFloat(0.0001f, 1.0f);
}
}
for(i = 0; i < n; i++){
for(j = 0; j < k; j++){
my_b[i*k+j] = randFloat(0.0001f, 1.0f);
}
}
for(i = 0; i < m; i++){
for(j = 0; j < k; j++) {
c_cpu[i*k+j] = randFloat(1.0, 20.0f);
c_hwa[i*k+j] = c_cpu[i*k+j];
}
}
alpha = 0.5;
beta = randFloat(1.0, 2.0);
fprintf( stdout, "---- Running calculations ----\n");
// run sgemm on GPU (NB_RUNS iterations)
printf("Run on GPU\n");
// Start timer
gettimeofday(&tv_global_begin, NULL);
for( i=0; i<NB_RUNS; i++ ) {
printf("%d ",i);
gettimeofday(&tv_begin, NULL);
mySgemm( m, n, k, alpha, beta, my_a, my_b, c_hwa );
gettimeofday(&tv_end, NULL);
current = (tv_end.tv_sec-tv_begin.tv_sec)*1e6 + tv_end.tv_usec-tv_begin.tv_usec;
if( ( best_measure_GPU == 0 ) || ( best_measure_GPU > current ) ){
best_measure_GPU = current;
}
sum_measure_GPU += current;
}
gettimeofday(&tv_global_end, NULL);
global_GPU_time = (tv_global_end.tv_sec-tv_global_begin.tv_sec)*1e6 + tv_global_end.tv_usec-tv_global_begin.tv_usec;
// run sgemm on CPU (NB_RUNS iterations)
printf("\n\nRun on CPU\n");
// Start timer
gettimeofday(&tv_global_begin, NULL);
for( i=0; i<NB_RUNS; i++ ) {
printf("%d ",i);
gettimeofday(&tv_begin, NULL);
mySgemm( m, n, k, alpha, beta, my_a, my_b, c_cpu );
gettimeofday(&tv_end, NULL);
current = (tv_end.tv_sec-tv_begin.tv_sec)*1e6 + tv_end.tv_usec-tv_begin.tv_usec;
if( ( best_measure_CPU == 0 ) || ( best_measure_CPU > current ) ){
best_measure_CPU = current;
}
sum_measure_CPU += current;
}
gettimeofday(&tv_global_end, NULL);
global_CPU_time = (tv_global_end.tv_sec-tv_global_begin.tv_sec)*1e6 + tv_global_end.tv_usec-tv_global_begin.tv_usec;
// Compute error between GPU and CPU
for( ii = 0; ii < m; ii++){
for(j = 0; j < k; j++){
double lerror = fabs((c_hwa[ii*k+j]-c_cpu[ii*k+j])/c_cpu[ii*k+j]);
if (lerror > error) {
error = lerror;
valueCPU = c_cpu[ii*k+j];
valueGPU = c_hwa[ii*k+j];
index_i = ii;
index_j = j;
}
}
}
if (error > 2e-06) {
fprintf( stdout, "\n\nThe error is %e with index %d:%d @ %e (CPU) / %e (GPU)\n", error, index_i, index_j, valueCPU, valueGPU);
fprintf( stdout, "The error is is too big!\n");
return -1;
}
fprintf( stdout, "\n\n---- Results ----");
fprintf( stdout, "\n");
fprintf( stdout, "Sizes of matrices: M:%i N:%i K:%i\n\n", m, n, k);
fprintf( stdout, "Best HWA time : %f ms\n", best_measure_GPU / 1e3 );
fprintf( stdout, "Mean HWA time : %f ms\n", sum_measure_GPU / NB_RUNS / 1e3);
fprintf( stdout, "\n");
fprintf( stdout, "Best CPU time : %f ms\n", best_measure_CPU / 1e3 );
fprintf( stdout, "Mean CPU time : %f ms\n", sum_measure_CPU / NB_RUNS / 1e3);
fprintf( stdout, "\n");
fprintf( stdout, "Global HWA time : %f ms\n", global_GPU_time / 1e3 );
fprintf( stdout, "Global CPU time : %f ms\n", global_CPU_time / 1e3 );
fprintf( stdout, "\n");
fprintf( stdout, "Speed-up : %f (computed on the best time)",
((float)best_measure_CPU)/best_measure_GPU);
fprintf( stdout, "\n");
free(my_a);
free(my_b);
free(c_hwa);
free(c_cpu);
return 0;
}