目次

C/C++で書いたループを高速化する(その1) ベクトル化
C/C++で書いたループを高速化する(その2) openACC <- ここ

はじめに

NvidiaコパイラでopenACCを使ってみました。コンパイラの入手方法はその1を見てください。そこそこの計算量がないと差が出にくいので、XYZの3次元座標1024点の相互距離の平均値を計算するプログラムを作ってみました。average_distance1とaverage_distance2の二つの関数を下に書いてあります。データの並びが違うだけです。その1に触れていますが、後者はCPUで実行するとSIMDを使える並びです。

CPUで実行する場合には
nvc++ -fast -Minfo=all distance1.cpp average_distance1.cpp
でコンパイルしています。関数の中にpragmaを記入していますがここでは無視されます。SIMDの使用の有無で速度の違いがでています。

次にGPUで実行する場合です。
nvc++ -acc -Minfo=all distance1.cpp average_distance1.cpp
両関数とものGPUが使われます。
CPUもGPUもだいぶ古いものですが、速くなっていることがわかります。データの並び方で差はほとんどないようです。
このようにループの直前に#pragma~の1行加えただけでGPUが使えるようになります。

結果

	CPU SIMD	CPU 実行時間(ms)	GPU(ms)
average_distance1	X	2.909	0.907
average_distance2	O	1.458	0.963

(実行時間はループ一回分に換算しています)

CPU Intel Core i7-2600 CPU 3.40GHz
GPU GTX 750
CentOS 8.2

使用したプログラム

average_distance1.cpp

# include <cmath>

inline double distance(double x1, double y1, double z1, double x2, double y2, double z2){
    return(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)));
}

double average_distance1(int n, const double xyz[][3]){
    double sum=0.0;
    #pragma acc kernels
    for(int i=0; i<n; i++)
        for(int j=i+1; j<n; j++)
            sum+=distance(xyz[i][0], xyz[i][1], xyz[i][2], xyz[j][0], xyz[j][1], xyz[j][2]);
    return(sum/(n*(n-1.0)/2.0));
}

double average_distance2(int n, const double xyz2[][1024]){
    double sum=0.0;
    #pragma acc kernels
    for(int i=0; i<n; i++)
        for(int j=i+1; j<n; j++)
            sum+=distance(xyz2[0][i], xyz2[1][i], xyz2[2][i], xyz2[0][j], xyz2[1][j], xyz2[2][j]);
    return(sum/(n*(n-1.0)/2.0));
}

distance1.cpp

# include <iostream>
# include <chrono>
# include <random>

using namespace std;

double average_distance1(int n, const double xyz[][3]);
double average_distance2(int n, const double xyz2[][1024]);

int main(){
    std::mt19937 mt(123);
    std::uniform_real_distribution<double> u_rand(-10.0, 10.0);
    
    double data1[1024][3];
    for(int i=0; i<1024; i++){
         data1[i][0]=u_rand(mt);
         data1[i][1]=u_rand(mt);
         data1[i][2]=u_rand(mt);
    }
    
    double data2[3][1024];
    for(int i=0; i<1024; i++){
         data2[0][i]=data1[i][0];
         data2[1][i]=data1[i][1];
         data2[2][i]=data1[i][2];
    }
    
    chrono::system_clock::time_point start, end;
    double time;
    volatile double result;
    
    start=chrono::system_clock::now();
    for(int i=0; i<1000; i++){
        result=average_distance1(1024, data1);
    }
    end=chrono::system_clock::now();
    time=static_cast<double>(chrono::duration_cast<chrono::microseconds>(end - start).count()/1000.0);
    cout << result << endl;
    cout << time << endl;
    
    start=chrono::system_clock::now();
    for(int i=0; i<1000; i++){
        result=average_distance2(1024, data2);
    }
    end=chrono::system_clock::now();
    time=static_cast<double>(chrono::duration_cast<chrono::microseconds>(end - start).count()/1000.0);
    cout << result << endl;
    cout << time << endl;
}