目次

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

はじめに

プログラムの速度的なボトルネックのほとんどはループにあり、多くの場合そこを改善することで高速化を図ることができます。
ここでは、Nvidiaから無料で使えるコンパイラが提供されており、それを使った高速化手法を紹介いたします。NVIDIA HPC SDKから入手できます。これまでPGIが販売していましたがNvidiaに移管されたもので、CPUの自動並列化機能やopenACCなどを備えており、また最適化情報も出力してくれますので、プログラムの高速化に役立つものと思われます。

ここではC++で書いたプログラムを示していますが、Cでもほぼ同じ方法が可能です。
今回はSIMDに限定した話ですが、パソコンで使われているSIMDは適用できる条件が厳しく、通常は最内ループで使われます。
それを使える条件は概ね次の通りです。

1. 連続したメモリ領域に書き込まれたデータを連続的にアクセスすること

1次元の配列の場合には端から順番にアクセスすれば良いのですが、
x[i][j]のように2次元の場合はiを固定して　j＝0,1,2...のようにアクセスするとメモリ領域を連続して順番にアクセスできます。逆方向はダメなようです。

2.ループ実行の順番が変わっても結果が変わらないこと

例えばループの最初の結果を次の計算で使うようなことがあると、最初の結果が出ないと次の計算ができないことになります。SIMDでは複数回分の計算を一度にやりますので、困ったことになります。
すなわちループ内の計算が独立している必要があります。ただし、合計や最大値、最小値を求める場合には、前の結果を使いますが、実行順が入れ替わっても同じ結果になりますので問題はありません。(並列計算reductionで調べてみてください)

3. ループに入る時に回る回数が分かっていること

ループをforやwhileなど違う書き方をしてみました。プログラムはtest1.cppをご覧ください。コンパイルした時のSIMD使用の有無と実行時間(100万回実行)を表にしています

	ループ	SIMD使用	実行時間 (msec)
sum1	範囲ベースfor	X	785
sum2	従来型for	O	207
sum3	while	O	207
sum4	while break	X	786

これらの関数は、連続したメモリ領域に保存されていることは全部に共通していますが、SIMDが使えるかはループに入る前に何回回るか事前にコンパイラが分かっているかで決まっているようです。sum1は回数が分かりそうなものですが、このコンパイラには分かってもらえないようです。コンパイラの癖がかなりありますので、コンパイラに理解しやすいように書く必要があります。ともかくSIMDを使うかで4倍近くの差ができます。(AVXはdoubleを4つ同時に計算できる)

4. ループ内で計算が完結すること

ループの中で二乗の合計計算しています。プログラムはtest2.cppをご覧ください。

	二乗関数	SIMD使用	実行時間 (msec)
sqsum1	ループ内に sum+=x[i]*x[i];	O	209
sqsum2	ループ内に double x1=x[i]; sum+=x1*x1; // 一旦変数の代入	O	208
sqsum3	ループ内に sum+=x[i]*x[i]; // pragmaでベクトル化とunrollを抑制	X	849
sqsum4	インライン関数	O	208
sqsum5	外部	X	3959

いろいろな書き方をしていますが、SIMDが使えると同じように高速化されています。sqsum3はpragmaでSIMDを使わないように指定しさらにunroll最適化も使わないようにしています。SIMD使用の有無で上と同じように4倍程度の差が出ています。sqsum4はインライン関数を使っており、ループ内に直接書いたのと同じ効果があります。インライン関数にすると関数コールのオーバーヘッドがなくなるだけではなく、ベクトル化の効果もあります。

5. ループ内に分岐を入れない

ループの中はSIMDの命令だけで処理できるように書く必要があります。if文などの分岐はSIMDでは処理できませんので遅くなってしまいます。できるだけコンパイラが提供している関数を使ってみましょう。CPUコード上で分岐にならないことが多いようです。だだし、絶対値や最大値の計算などをif文で書いてもコンパイラが気を利かせてくれる場合もあります。

Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz AVX
CentOS 8.2
nvc++ 20.7

使用したプログラム

ループに入る時に回る回数が分かっていること

test1f.cpp

#include <vector>
using namespace std;

double sum1(const vector<double>& x){
    double sum=0.0;
    for(auto& xi : x){
        sum+=xi;
    }
    return(sum);
}

double sum2(const vector<double>& x){
    double sum=0.0;
    for(int i=0; i<x.size(); i++){
        sum+=x[i];
    }
    return(sum);
}

double sum3(const vector<double>& x){
    double sum=0.0;
    int i=0;
    while(i < x.size()){
        sum+=x[i];
        i++;
    }
    return(sum);
}

double sum4(const vector<double>& x){
    double sum=0.0;
    int i=0;
    while(true){
        if(i >= x.size()) break;
        sum+=x[i];
        i++;
    }
    return(sum);
}

test1.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <chrono>

using namespace std;

double sum1(const vector<double>& x);
double sum2(const vector<double>& x);
double sum3(const vector<double>& x);
double sum4(const vector<double>& x);

int main(){
    vector<double> data(1000); // 入力データの作成
    for(int i=0; i<1000; i++){
        data[i]=sqrt(i);
    }
    
    chrono::system_clock::time_point start, end;
    double time;
    
    volatile double result;
    result=0.0;
    start=chrono::system_clock::now();
    for(int i=0; i<1000000; i++){
        result+=sum1(data);
    }
    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;
/* 以下省略
sum2(data)
sum3(data)
sum4(data)
についても同様に
*/

コマンド

nvc++ -Minfo=vect -fast -c test1f.cpp
nvc++ -O0 test.cpp test1f.o
./a.out

ループ内で計算が完結すること

test2f.cpp

#include <vector>
using namespace std;

double sqsum1(const vector<double>& x){
    double sum=0.0;
    for(int i=0; i<x.size(); i++){
        sum+=x[i]*x[i];
    }
    return(sum);
}

double sqsum2(const vector<double>& x){
    double sum=0.0;
    for(int i=0; i<x.size(); i++){
        double x1=x[i];
        sum+=x1*x1;
    }
    return(sum);
}

double sqsum3(const vector<double>& x){
    double sum=0.0;
    #pragma loop novector
    #pragma loop nounroll
    for(int i=0; i<x.size(); i++){
        sum+=x[i]*x[i];
    }
    return(sum);
}

inline double sq4(double x){return(x*x);}

double sqsum4(const vector<double>& x){
    double sum=0.0;
    for(int i=0; i<x.size(); i++){
        sum+=sq4(x[i]);
    }
    return(sum);
}

double sq5(double x){return(x*x);} // 外部から使う

test2fn.cpp

#include <vector>
using namespace std;

double sq5(double x);

double sqsum5(const vector<double>& x){
    double sum=0.0;
    for(int i=0; i<x.size(); i++){
        sum+=sq5(x[i]);
    }
    return(sum);
}

test2.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <chrono>
using namespace std;

double sqsum1(const vector<double>& x);
double sqsum2(const vector<double>& x);
double sqsum3(const vector<double>& x);
double sqsum4(const vector<double>& x);
double sqsum5(const vector<double>& x);

int main(){
    vector<double> data(1000);
    for(int i=0; i<1000; i++){
        data[i]=sqrt(i);
    }
    
    chrono::system_clock::time_point start, end;
    double time;
    
    volatile double result;
    result=0.0;
    start=chrono::system_clock::now();
    for(int i=0; i<1000000; i++){
        result+=sqsum1(data);
    }
    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;
    
/* 以下省略
seqsum2(data)
seqsum3(data)
seqsum4(data)
seqsum5(data)
についても同様に
*/

nvc++ -Minfo=all -fast -c test2f.cpp
nvc++ -Minfo=all -fast -c test2fn.cpp -Mipa=inline
nvc++ -O0 test2.cpp test2f.o test2fn.o