一步步做程序优化-讲一个用于OpenACC优化的程序（转载）

一步步做程序优化【1】讲一个用于OpenACC优化的程序

分析下A,B,C为三个矩阵，A为m*n维，B为n*k维，C为m*k维，用A和B来计算C，计算方法是：C = alpha*A*B + beta*C。它的程序如下：

// 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;
    for(j = 0; j < m; j++) 
    {
        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_Tutorial_Labs_CUDA中的lab0。

C中的每个元素的计算是独立的，完全可以并行化。后面的系列文章将会讲述优化这个程序的整个过程。

一步步做程序优化【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, "
**** error : memory allocation failed ****

");
    return 1;
  }

  fprintf( stdout, "---- Initialization of the Matrices ----

");
  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 ----
");


  // run sgemm on GPU (NB_RUNS iterations)
  printf("Run on GPU
");

  // 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("

Run on CPU
");

  // 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, "

The error is %e with index %d:%d @ %e (CPU) / %e (GPU)
", error, index_i, index_j, valueCPU, valueGPU);
    fprintf( stdout, "The error is is too big!
");
    return -1;
  }


  fprintf( stdout, "

---- Results ----");
  fprintf( stdout, "
");
  fprintf( stdout, "Sizes of matrices: M:%i  N:%i  K:%i

", m, n, k);
  fprintf( stdout, "Best HWA time    : %f ms
", best_measure_GPU / 1e3 );
  fprintf( stdout, "Mean HWA time    : %f ms
", sum_measure_GPU / NB_RUNS / 1e3);
  fprintf( stdout, "
");
  fprintf( stdout, "Best CPU time    : %f ms
", best_measure_CPU / 1e3 );
  fprintf( stdout, "Mean CPU time    : %f ms
", sum_measure_CPU / NB_RUNS / 1e3);
  fprintf( stdout, "
");
  fprintf( stdout, "Global HWA time  : %f ms
", global_GPU_time / 1e3 );
  fprintf( stdout, "Global CPU time  : %f ms
", global_CPU_time / 1e3 );
  fprintf( stdout, "
");
  fprintf( stdout, "Speed-up         : %f (computed on the best time)",
           ((float)best_measure_CPU)/best_measure_GPU);


  fprintf( stdout, "
");

  free(my_a);
  free(my_b);
  free(c_hwa);
  free(c_cpu);

  return 0;
}

注意上述代码中，测试了两次统一个函数的执行结果，下面加入两句简单的指令，然后编译执行下，看一下加速比情况。

在第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的计算资源了。种种好处 只有在你试用之后才能知道的奥。

后面的博客我还会讲解更多的指令，还有一些有意思的细节。欢迎大家关注奥！