参照条件
トピックの一部として、アーキテクチャ指向のソフトウェア最適化のための小さな競争を発表します。
要するに、コードはソースデータを一見して操作する一連の関数と、最適化されていないバージョンをn回実行し、次に最適化されたバージョンを実行し、桁数を比較し、それらが一致する場合にランタイム率を与えるドライバーローションで構成されます。 このように：
元のコードを実行しています...完了
最適化されたコードの実行...完了
結果を確認しています...合格
実行回数：3
元のコードの平均時間：11.954537秒。
最適化されたコード平均時間：1.052994秒。
スピードアップ：11.35
任意の最適化手法、任意のオプションを備えたGCCコンパイラ、および2つのクアッドコアIntel Xeon E5420 2.5 GHzプロセッサを搭載したサーバーを使用できます。
ところで、ここにコードがあります：
/*
* contestopt.c:
*/
#include "contest.h"
void function3( double a[][NUM], double b[][NUM], double c[][NUM]);
void function2( double a[][NUM]);
double function4(unsigned int second);
double ***p;
double ****x;
double ****y;
double ****z;
double ***wrk1;
double ***wrk2;
double ***bnd;
static double a[NUM][NUM], b[NUM][NUM];
void function1( double result[][NUM], double gos[2], unsigned int second)
{
function2(a);
function2(b);
function3(a, b, result);
gos[0] = function4(second);
}
void function3( double a[][NUM], double b[][NUM], double result[][NUM])
{
int i, j, k;
for (i = 0; i < NUM; i++) {
for (j = 0; j < NUM; j++) {
for (k = 0; k < NUM; k++) {
result[i][j] = result[i][j] + a[i][k] * b[k][j];
}
}
}
}
void function2( double a[][NUM])
{
int i, j;
double first = 0.0;
double second = 1.0;
a[0][0] = first;
a[0][1] = second;
for (i = 0; i < NUM; i++) {
for (j = 0; j < NUM; j++) {
if (!(i == 0 && (j == 0 || j == 1))) {
a[i][j] = first + second;
first = second;
second = a[i][j];
}
if (j % 20 == 0 && j != 0) {
first = first * (j + 1) / (NUM);
second = second * (j + 1) / (NUM);
}
}
first = ((i + 1) * NUM) / first;
second = ((i + 1) * NUM) / second;
}
}
void function5( double ***A, double val)
{
int i, j, k;
for (i = 0; i < XX; i++)
for (j = 0; j < YY; j++)
for (k = 0; k < ZZ; k++)
A[i][j][k] = val;
}
void function6(unsigned int second)
{
int mimax = XX;
int mjmax = YY;
int mkmax = ZZ;
int imax = mimax - 1;
int jmax = mjmax - 1;
int kmax = mkmax - 1;
int i, j, k, l;
double val;
srand((unsigned int ) second);
p = ( double ***) malloc( sizeof ( double **) * XX);
wrk1 = ( double ***) malloc( sizeof ( double **) * XX);
wrk2 = ( double ***) malloc( sizeof ( double **) * XX);
bnd = ( double ***) malloc( sizeof ( double **) * XX);
x = ( double ****) malloc( sizeof ( double ***) * XX);
y = ( double ****) malloc( sizeof ( double ***) * XX);
z = ( double ****) malloc( sizeof ( double ***) * XX);
for (i = 0; i < XX; i++) {
p[i] = ( double **) malloc( sizeof ( double *) * YY);
wrk1[i] = ( double **) malloc( sizeof ( double *) * YY);
wrk2[i] = ( double **) malloc( sizeof ( double *) * YY);
bnd[i] = ( double **) malloc( sizeof ( double *) * YY);
x[i] = ( double ***) malloc( sizeof ( double **) * YY);
y[i] = ( double ***) malloc( sizeof ( double **) * YY);
z[i] = ( double ***) malloc( sizeof ( double **) * YY);
for (j = 0; j < YY; j++) {
p[i][j] = ( double *) malloc( sizeof ( double ) * ZZ);
wrk1[i][j] = ( double *) malloc( sizeof ( double ) * ZZ);
wrk2[i][j] = ( double *) malloc( sizeof ( double ) * ZZ);
bnd[i][j] = ( double *) malloc( sizeof ( double ) * ZZ);
x[i][j] = ( double **) malloc( sizeof ( double *) * ZZ);
y[i][j] = ( double **) malloc( sizeof ( double *) * ZZ);
z[i][j] = ( double **) malloc( sizeof ( double *) * ZZ);
for (k = 0; k < ZZ; k++) {
x[i][j][k] = ( double *) malloc( sizeof ( double ) * 4);
y[i][j][k] = ( double *) malloc( sizeof ( double ) * 3);
z[i][j][k] = ( double *) malloc( sizeof ( double ) * 3);
}
}
}
val = ( double ) rand() / (10.0 * RAND_MAX);
for (i = 0; i < XX; i++)
for (j = 0; j < YY; j++)
for (k = 0; k < ZZ; k++) {
for (l = 0; l < 3; l++) {
x[i][j][k][l] = 1.0;
y[i][j][k][l] = val;
z[i][j][k][l] = 1.0;
}
x[i][j][k][3] = 1.0;
}
function5(bnd, 1.0);
function5(wrk1, 0.0);
function5(wrk2, 0.0);
for (i = 0; i < imax; i++)
for (j = 0; j < jmax; j++)
for (k = 0; k < kmax; k++)
p[i][j][k] = ( double ) (k * k) / ( double ) (kmax * kmax);
}
void function7()
{
int i, j, k;
for (i = 0; i < XX; i++) {
for (j = 0; j < YY; j++) {
free(p[i][j]);
free(wrk1[i][j]);
free(wrk2[i][j]);
free(bnd[i][j]);
for (k = 0; k < ZZ; k++) {
free(x[i][j][k]);
free(y[i][j][k]);
free(z[i][j][k]);
}
free(x[i][j]);
free(y[i][j]);
free(z[i][j]);
}
free(p[i]);
free(wrk1[i]);
free(wrk2[i]);
free(bnd[i]);
free(x[i]);
free(y[i]);
free(z[i]);
}
free(x);
free(y);
free(z);
free(p);
free(wrk1);
free(wrk2);
free(bnd);
}
double function4(unsigned int second)
{
int nn = NN;
double gosa, s0, ss;
int mimax = XX;
int mjmax = YY;
int mkmax = ZZ;
int imax = mimax - 1;
int jmax = mjmax - 1;
int kmax = mkmax - 1;
int iCount, jCount, kCount, loop, i, j, k;
function6(second);
for (loop = 0; loop < nn; loop++) {
gosa = 0.0;
for (k = 1; k < kmax - 1; k++)
for (j = 1; j < jmax - 1; j++)
for (i = 1; i < imax - 1; i++) {
s0 = x[i][j][k][0] * p[i + 1][j][k]
+ x[i][j][k][1] * p[i][j][k]
+ x[i][j][k][2] * p[i][j][k]
+ y[i][j][k][0] * (p[i + 1][j + 1][k] - p[i + 1][j - 1][k]
- p[i - 1][j + 1][k] + p[i - 1][j - 1][k])
+ y[i][j][k][1] * (p[i][j + 1][k + 1] - p[i][j - 1][k + 1]
- p[i][j + 1][k - 1] + p[i][j - 1][k - 1])
+ y[i][j][k][2] * (p[i + 1][j][k + 1] - p[i - 1][j][k + 1]
- p[i + 1][j][k - 1] + p[i - 1][j][k - 1])
+ z[i][j][k][0] * p[i - 1][j][k]
+ z[i][j][k][1] * p[i][j - 1][k]
+ z[i][j][k][2] * p[i][j][k - 1] + wrk1[i][j][k];
ss = (s0 * x[i][j][k][3] - p[i][j][k]) * bnd[i][j][k];
gosa = gosa + ss;
wrk2[i][j][k] = p[i][j][k] + 0.1 * ss;
}
for (iCount = 1; iCount < imax - 1; iCount++)
for (jCount = 1; jCount < jmax - 1; jCount++)
for (kCount = 1; kCount < kmax - 1; kCount++)
p[iCount][jCount][kCount] = wrk2[iCount][jCount][kCount];
}
function7();
return gosa;
}
double exponent( int i, double x)
{
double z = 1.0;
int j;
for (j = 0; j < i; j++) {
z = z * x;
}
return z;
}
double function8( float *a, float p)
{
int i, j;
double h, x;
double s = p;
h = 1.0 / VAL;
for (i = 1; i <= VAL; i++) {
x = h * (i - 0.5);
for (j = 0; j < ORDER; j++) {
s = (a[j] * exponent(j, x)) + s;
}
}
return s * h;
}
* This source code was highlighted with Source Code Highlighter .
値でテスト：
#NUM 600を定義
#define XX 65
#def YY 33
#define ZZ 33
#define NN 50
#define VAL 10000000
#define ORDER 30
#define TOLERANCE 1e-8
おおよその種類のドライバー：
/*
* driver.c: Contest driver.
*/
#include <sys/time.h>
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <time.h>
#include "contest.h"
#include "hpctimer.h"
double _RES[NUM][NUM] = { 0.0 };
double _GOS[2] = { 0.0, 0.0 };
double RES[NUM][NUM] = { 0.0 };
double GOS[2] = { 0.0, 0.0 };
float A[ORDER] = { 0.0 };
double function8( float *a, float p);
void function1( double result[][NUM], double gos[2], unsigned int second);
int is_results_eq( double _res[][NUM], double _gos[2], float _p,
double res[][NUM], double gos[2], float p)
{
int i, j;
if (fabsf(_p - p) > 1E-1) {
printf( "_p p\n" );
printf( "%f\n" , _p);
printf( "%f\n" , p);
return 0;
}
if (fabs(_gos[0] - gos[0]) > 1E-1) {
printf( "_gos gos\n" );
printf( "%f\n" , _gos[0]);
printf( "%f\n" , gos[0]);
return 0;
}
for (i = 0; i < NUM; i++)
for (j = 0; j < NUM; j++) {
if (fabs(_RES[i][j] - RES[i][j]) > TOLERANCE) {
printf( "[%d, %d] =\n%f\n%f\n" , i, j, _RES[i][j], RES[i][j]);
return 0;
}
}
return 1;
}
int main( int argc, char **argv)
{
unsigned int second;
unsigned int i, j;
float p = 9.0;
double _res, res;
int nruns = 3;
hpctimer_time_t t0, t1;
double torig = 0.0, topt = 0.0;
srand(time(NULL));
second = 1219304613 + (1000 - (2000.0 * ( double ) rand() / ( double ) (RAND_MAX + 1.0)));
for (i = 0; i < ORDER; i++) {
A[i] = ( float )i;
}
for (i = 0; i < NUM; i++) {
for (j = 0; j < NUM; j++) {
_RES[i][j] = 0.0;
RES[i][j] = 0.0;
}
}
/* Original code */
printf( "Executing original code ... " );
fflush(stdout);
t0 = hpctimer_gettime();
for (i = 0; i < nruns; i++) {
_function1(_RES, _GOS, second);
_res = _function8(A, p);
}
t1 = hpctimer_gettime();
torig = hpctimer_getdiff(t0, t1) / nruns;
printf( "done\n" );
fflush(stdout);
/* Optimized code */
printf( "Executing optimized code ... " );
fflush(stdout);
t0 = hpctimer_gettime();
for (i = 0; i < nruns; i++) {
function1(RES, GOS, second);
res = function8(A, p);
}
t1 = hpctimer_gettime();
topt = hpctimer_getdiff(t0, t1) / nruns;
printf( "done\n" );
fflush(stdout);
printf( "Checking results ... " );
fflush(stdout);
if (!is_results_eq(_RES, _GOS, _res, RES, GOS, res)) {
printf( "RESULTS ARE DIFFERENT!\n" );
} else {
printf( "PASSED\n" );
printf( "Number of runs: %d\n" , nruns);
printf( "Original code average time: %.6f sec.\n" , torig);
printf( "Optimized code average time: %.6f sec.\n" , topt);
printf( "Speedup: %.2f\n" , torig / topt);
}
return 0;
}
* This source code was highlighted with Source Code Highlighter .
加速、ルール1
ソースコードレベルでの最適化について言うことができます：「コードが何をするかは問題ではありません。どのようにそれを行うかが重要です。」 しかし、私が思いついた最初の最適化ルールは、この声明に多少反しています：
最初のルール
喫煙材料とアルゴリズム
コードを一見すると、彼が肛門を介していくつかの簡単なアクションを実行することが明らかになります。
そのため、大まかに一目見ただけでも、嫌な気持ちで指数関数を切り取ることができます。これにより、有用なジェスチャの数が大幅に設定されます（プロファイラを実行すると、合計カウント時間の80％を消費します）。 彼女が扱う博覧会は、function8で簡単に動的に実装されます。
function8自体では、積分の計算は無意識のうちに表示されますが、その最適化の範囲については後で説明します。
function4で使用される差分スキームは、コードを書き直しても最適化されていませんが、gosaの蓄積を別のループに取り込むことで、さらに高速化できます。
function3の行列乗算はさまざまな方法で加速されますが、最初の段階では、乗算された行列が完全に同一であり、この関数で消費されるメモリが半分になることに気付くだけで十分です。 さらに、行列を二乗するエレガントなアルゴリズムをグーグルで検索します。 Strassen、Coppersmith-Vinohrad、その他の変態については話していない。
加速、ルール2
コードの作成者がそれほど多くの冗長性を詰め込んでいるかどうかはわかりませんが、得られた高速化の主な部分は第2規則によるものです。
第二のルール
意味のないデータ構造を破壊する
ちなみに、このルールは美しく、最適化だけでなく、読みやすさとコンパイルの改善にも有効です。 ある素晴らしい女性を叙情的に引用して、「それが意味をなさない再帰を使う意味はない」。
この意味で、function4は下品なところまでロードされます。 4次元配列x、y、zがどのように初期化されるかに気づくのを恐れて、すぐにそれらを破壊します。 これらの後にwrk1とbndを続けることができます。 mimax、iCount、ssなどの些細なことについては話していない。 その結果、コードははるかに短くなり、見た目も美しくなります。
function6とfunction7は完全に削除され、有用なアクションの残りはfunction4を論理的または時間的にロードしません。したがって、それらは無痛で転送されます。 function5は突然自己破壊します。 function3では、既に述べたように、1つのソースマトリックスを使用でき、function2を1回呼び出すことができます。
すでに呼吸がずっと簡単になり、結果としてスピードアップが悪用を促し始めます。
加速、ルール3
第三のルール
変換予測エラーを最小限に抑える
今からすべてを説明します。 これはブランチ最適化手法です。 コンパイルされたループまたは分岐演算子は、次のようになります（大幅に簡略化されています）。
1 mov eax, ebx
2 jne $a
3 mov ecx, ebx
4 jmp b
5 a: ...
6 b: ...
パイプラインでこの時点で何が起こっているか想像してみてください。 命令はいくつかの段階で処理され（フェッチ-デコード-実行-ライトバック）、命令2を読み込んだ後、 遷移予測モジュールは動作するはずです。この操作の結果がまだわからず、次の命令を選択できないためです。
遷移予測モジュールが「間違っている」場合はどうなりますか？ そうです-コンベアが降ろされています。 そして、これはオーバーヘッドです。 遷移予測モジュールの呼び出し回数を減らすことで、遷移予測モジュールのエラー数を減らすことができます。
どうやって？
- 最も可能性の高いブランチをブランチの先頭に配置します
- （ループ内のネストに関して）不変分岐をより高くします
- ループ展開を使用する
for (i = 0; i < n; i++) {
sum += p[i];
}
書く
for (i = 0; i < n; i+=4) {
sum += p[i];
sum += p[i + 1];
sum += p[i + 2];
sum += p[i + 3];
}
遷移予測の数を4倍減らします。
この手法には、展開されたピースの反復回数の多様性、最適な展開の長さなどの落とし穴があります。 しかし、個々のケースでは、バランスを見つけて、この比較的簡単な方法で余分な加速をつかむことができます。
ループの展開の例は、function3の最適化バージョンで見ることができます。
function2で非常に無意味な不変ブランチのサイクルから削除を追加するのは良いことですが、そこからスクラッチするのはそれほど簡単ではありません。 誰が私より優れた実装をしているか-よくやった。
読み物：R.ガーバー、A。ビル「ソフトウェアの最適化：レシピ集」 2010、インテル
加速ルール4
第4規則
メモリを効率的に使用する
特に、キャッシュメモリを使用します。 リンクの局所性という特性があります -同じオブジェクトに繰り返し/頻繁にアクセスするプログラムの特性。 また、一時的な（一時的な）ローカリゼーションと空間的な（空間的な）ローカリゼーションがあります。 キャッシングはプリエンプティブであり、これに基づいて、キャッシュミスを避けるためにシリアルアドレスに連絡することは理にかなっています。
同じ行列乗算で、行はメモリ内に連続して配置され、列は分散され、2番目の行列の転置により既知の加速効果が得られます。
話をしたので、例を挙げます：構造（struct）がキャッシュに収まらず、その中に配列が含まれる場合、構造内の配列ではなく、構造の配列を使用する方が良いです。
また、移植可能なプログラムを作成するために、キャッシュのサイズにアタッチしないでください。
そして、他の多くの微妙な点。 読む-Chris Kaspersky 「ソフトウェア最適化テクニック。 メモリの効率的な使用 。 」 2003年、BHV-ピーターズバーグ
加速ルール5
どんなプロセッサアーキテクチャがありますか？
第5規則
可能な場合は計算をベクトル化します
幸いなことに、SSEを使用するときが来ました。幸いなことに、配列で機能するほとんどすべてのコードをベクトル化できます。
小さな教育プログラム。
SSEアーキテクチャ（x86）：
- MMX-整数ベクトル処理（64ビットレジスタ）
- SSE-実数ベクトルと整数ベクトル（128ビット）で動作します
mulss xmm1, xmm0 | mulps xmm1, xmm0 XMM0 4.0 3.0 2.0 1.0 | 4.0 3.0 2.0 1.0 * | * * * * XMM1 5.0 5.0 5.0 5.0 | 5.0 5.0 5.0 5.0 = | = = = = XMM1 4.0 3.0 2.0 5.0 | 20.0 15.0 10.0 5.0
SSE命令は、コピー、算術、比較、ビット演算、およびあらゆる種類のトリッキーな変換を実行できます。
それらの使用方法は？ アセンブラーベクトル命令があり、特別なC ++クラスがあり、結局のところ、自動ベクトル化があります-最も移植性があります。 組み込み関数を使用したかった。
組み込み関数は、SSE命令への高レベルのアクセスを提供するためにコンパイラーによってサポートされる関数とデータ型のセットです。
xmmintrin.hを接続し、ベクトル化されたコードを記述し、-msse3フラグでコンパイルします。 コンパイラは、XMM＃レジスタを割り当て、コマンドとアドレッシングモードをスケジュールします。 便利で直感的にシンプル。
たとえば、1つの命令で4つの浮動小数点数を追加します。
#include <xmmintrin.h>
void add( float *a, float *b, float *c)
{
__mm128 t0, t1;
t0 = _mm_load_ps(a);
t1 = _mm_load_ps(b);
t0 = _mm_add_ps(t0, t1);
_mm_store_ps(c, t0);
}
この特定のケースでは、行列の乗算が組み込み関数にコピーされます（同時にブロックになります）。 function8も原則としてベクトル化されますが、ここには特別なケースがありますので、今のところは忍耐に努めています。
この段階で、コードは夢を描き始め、各オペレーターは厚いサイクルから切り取られ、生産性の目に見える飛躍をもたらします。
加速、ルール6
境界条件を見てください。 可能な場合は、データ構造をより単純なものに置き換えます。 高速アルゴリズムが存在する問題に計算を減らします。
第六規則
モスクを使用
ドライバーでは、行列乗算の結果がTOLERANCEと正確に比較され、他の2つの数値は10分の1まで正確です。 科学的な突く方法を使用すると、function8では、間隔の3分の2が必要な精度を達成するのに十分であり、最初の300万（ちょうど）の反復は破棄できることがわかります。 さらに、6から7千回の反復ごとに考慮することで十分です。
取得します
double function8( float *a, float p)
{
int i, j, val, step;
double h, x, z;
double s = p;
h = 1.0 / VAL;
val = 0.33 * VAL;
step = val / 500;
for (i = val; i <= VAL; i += step) {
x = h * i;
z = 1.0;
for (j = 1; j < ORDER; j++) {
z *= x;
s += j * z;
}
}
return s * h * step;
}
そして喜ぶ。 好きですか？ しません あなたは完全にべき乗の操作を取り除くことができます。
double function8( float *a, float p)
{
unsigned int j;
double s = p / VAL;
for (j = 0; j < ORDER; j++) {
s += (a[j] / (j + 1));
}
return s;
}
元のO（VAL）* O（ORDER ^ 2）と比較して、O（ORDER）の複雑さがわかります。 同じ結果は、ベルヌーイ数を使用することで達成できます（それほど簡単ではありませんが、アルゴリズムを手動で作成することで確認できます）。
2乗行列の塗りつぶしの特性に注意を払い、偶数行のゼロを確認できます。
3次元配列を1次元配列に置き換えて、次のようにアドレス指定できます。
m1 =（YY-1）*（ZZ-1）; m2 =（ZZ-1）; p [i * m1 + j * m2 + k]。 この小さなトリックにより、memsetは多次元行列を埋めることができます。
小さな関数には__inlineを使用できますし、使用する必要さえあります。
その他。 この場合、8コアのマシンが利用できるため、OpenMPを使用できません。 function3、function4、およびfunction8の並列化には、多くの想像力さえ必要ありません。 もう少しわかりにくいのは、function1で並列セクションを使用する可能性です（function3とfunction4はデータに依存しないためです）。
デザートの場合-最適化されたプログラムのバージョン。 それは理想的なふりではありませんが、正直な問題の加速は数百回です。
/*
* contestopt.c:
*/
#include <sys/time.h>
#include <omp.h>
#include "contest.h"
void function3( double a[][NUM], double c[][NUM]);
void function2( double a[][NUM]);
double function4(unsigned int second);
static double p[XX][YY][ZZ] __attribute__ ((aligned(64)));
static double wrk2[XX][YY][ZZ] __attribute__ ((aligned(64)));
static double a[NUM][NUM] __attribute__ ((aligned(64)));
void function1( double result[][NUM], double gos[2], unsigned int second)
{
#pragma omp parallel sections
{
#pragma omp section
function3(a,result);
#pragma omp section
gos[0] = function4(second);
}
}
void function3( double a[][NUM], double result[][NUM])
{
int i, i2, j, j2, k, k2;
double *restrict rres;
double *restrict rmul1;
double *restrict rmul2;
unsigned int BS = 8;
int remainder = NUM % 8;
int limit = NUM - remainder;
function2(a);
if (!(NUM & 1)) {
#pragma omp for private (k, j, i, i2, j2, k2, rres, rmul1, rmul2)
for (i = 0; i < limit; i += BS)
for (j = 0; j < limit; j += BS)
for (k = 0; k < limit; k += BS)
for (i2 = 0, rres = &result[i][j], rmul1 = &a[i][k]; i2 < BS;
++i2, rres += NUM, rmul1 += NUM) {
_mm_prefetch (&rmul1[8], _MM_HINT_NTA);
if (!(NUM & 1))
for (k2 = 0, rmul2 = &a[k][j]; k2 < BS; ++k2, rmul2 += NUM) {
__m128d m1d = _mm_load_sd (&rmul1[k2]);
m1d = _mm_unpacklo_pd (m1d, m1d);
j2 = 0;
__m128d m2 = _mm_load_pd (&rmul2[j2]);
__m128d r2 = _mm_load_pd (&rres[j2]);
_mm_store_pd (&rres[j2],
_mm_add_pd (_mm_mul_pd (m2, m1d), r2));
j2 +=2;
m2 = _mm_load_pd (&rmul2[j2]);
r2 = _mm_load_pd (&rres[j2]);
_mm_store_pd (&rres[j2],
_mm_add_pd (_mm_mul_pd (m2, m1d), r2));
j2 +=2;
m2 = _mm_load_pd (&rmul2[j2]);
r2 = _mm_load_pd (&rres[j2]);
_mm_store_pd (&rres[j2],
_mm_add_pd (_mm_mul_pd (m2, m1d), r2));
j2 +=2;
m2 = _mm_load_pd (&rmul2[j2]);
r2 = _mm_load_pd (&rres[j2]);
_mm_store_pd (&rres[j2],
_mm_add_pd (_mm_mul_pd (m2, m1d), r2));
}
}
} else {
#pragma omp for private (k, j, i, i2, j2, k2, rres, rmul1, rmul2)
for (i = 0; i < limit; i += BS)
for (j = 0; j < limit; j += BS)
for (k = 0; k < limit; k += BS)
for (i2 = 0, rres = &result[i][j], rmul1 = &a[i][k]; i2 < BS;
++i2, rres += NUM, rmul1 += NUM) {
_mm_prefetch (&rmul1[8], _MM_HINT_NTA);
for (k2 = 0, rmul2 = &a[k][j]; k2 < BS; ++k2, rmul2 += NUM) {
for (j2 = 0; j2 < BS; j2++ ) {
rres[j2] += rmul1[k2] * rmul2[j2];
}
}
}
}
if (remainder) {
#pragma omp for private (i,j,k)
for (i = 0; i < limit; ++i)
for (k = NUM - remainder; k < NUM; ++k)
for (j = 0; j < limit; ++j)
result[i][j] += a[i][k] * a[k][j];
#pragma omp for private (i,j,k)
for (i = limit; i < NUM; ++i)
for (k = 0; k < NUM; ++k)
for (j = 0; j < NUM; ++j)
result[i][j] += a[i][k] * a[k][j];
#pragma omp for private (i,j,k)
for (i = 0; i < limit; ++i)
for (k = 0; k < NUM; ++k)
for (j = limit; j < NUM; ++j)
result[i][j] += a[i][k] * a[k][j];
}
}
void function2( double a[][NUM])
{
int i, j ;
double first = 0.0;
double second = 1.0;
__assume_aligned(a, 64);
a[0][0]=first;
a[0][1]=second;
for (j = 2; j < NUM; j++) {
a[0][j] = first + second;
first = second;
second = a[0][j];
if ( j%20 == 0 && j != 0) {
first = first * (j + 1) / (NUM);
second = second * (j + 1) / (NUM);
}
}
first = NUM / first;
second = NUM / second;
for (i = 1; i < NUM; i++) {
for (j = 0; j < NUM; j++) {
a[i][j] = first + second;
first = second;
second = a[i][j];
if ( j%20 == 0 && j != 0 ) {
first = first * (j + 1) / (NUM);
second = second * (j + 1) / (NUM);
}
}
first = ((i + 1) * NUM) / first;
second = ((i + 1) * NUM) / second;
}
}
__inline void function6(unsigned int second)
{
int imax = XX - 1;
int jmax = YY - 1;
int kmax = ZZ - 1;
int i, j, k;
double kmaxDoubled = kmax * kmax;
for (k = 0; k < kmax; k++)
p[0][0][k] = (k * k) / kmaxDoubled;
for (i = 0; i < imax; i++)
for (j = 0; j < jmax; j++)
for (k = 1; k < kmax; k++)
p[i][j][k] = p[0][0][k];
}
double function4(unsigned int second)
{
int nn = NN;
double gosa = 0.0, val;
int imax = XX - 2;
int jmax = YY - 2;
int kmax = ZZ - 2;
int loop, i, j, k;
function6(second);
srand((unsigned int )second);
val = ( double ) rand() / (10.0 * RAND_MAX);
for (loop = 1; loop < nn; loop++) {
for (k = 1; k < kmax; k++)
for (j = 1; j < jmax; j++)
for (i = 1; i < imax; i++) {
wrk2[i][j][k] = 0.1 * ( p[i + 1][j][k] + p[i][j][k]
+ val * (p[i + 1][j + 1][k] - p[i + 1][j - 1][k]
- p[i - 1][j + 1][k] + p[i - 1][j - 1][k]
+ p[i][j + 1][k + 1] - p[i][j - 1][k + 1]
- p[i][j + 1][k - 1] + p[i][j - 1][k - 1]
+ p[i + 1][j][k + 1] - p[i - 1][j][k + 1]
- p[i + 1][j][k - 1] + p[i - 1][j][k - 1])
+ p[i - 1][j][k] + p[i][j - 1][k] + p[i][j][k - 1]);
}
for (i = 1; i < imax; i++)
for (j = 1; j < jmax; j++)
for (k = 1; k < kmax; k++) p[i][j][k] += wrk2[i][j][k];
}
for (k = 1; k < kmax; k++)
for (j = 1; j < jmax; j++)
for (i = 1; i < imax; i++) {
gosa += p[i + 1][j][k] + p[i][j][k]
+ val * (p[i + 1][j + 1][k] - p[i + 1][j - 1][k]
- p[i - 1][j + 1][k] + p[i - 1][j - 1][k]
+ p[i][j + 1][k + 1] - p[i][j - 1][k + 1]
- p[i][j + 1][k - 1] + p[i][j - 1][k - 1]
+ p[i + 1][j][k + 1] - p[i - 1][j][k + 1]
- p[i + 1][j][k - 1] + p[i - 1][j][k - 1])
+ p[i - 1][j][k] + p[i][j - 1][k] + p[i][j][k - 1];
}
return gosa;
}
double function8( float *a, float p)
{
unsigned int j;
double s = 0.0;
#pragma omp parallel
{
#pragma omp for private (j) reduction (+:s)
for (j = 0; j < ORDER; j++) {
s += (a[j] /(j + 1));
}
}
return (s + p / VAL);
}
* This source code was highlighted with Source Code Highlighter .
これはとても楽しくて非常に役立つ実用的な例です。 コード最適化スキームを要約します。
- アルゴリズム最適化
- アーキテクチャに依存しない最適化：
- パラレルストリームI / O
- 並列化
- 効率的なメモリ処理、効率的なキャッシング
- アーキテクチャ指向の最適化：
- ベクトル命令
- 遷移予測エラーの数を最小限に抑える
- プロセッサメーカーのライブラリ関数
- グラフィックカードで計算を行う
- など
ご清聴ありがとうございました。