一步步做程序優化【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;
}