この場合、彼らは通常、並列化がプログラム実行の結果に何らかの影響を与えるべきではないという考えから進んでいます。これは厳しい要件ですが、多くの場合公平です。ただし、数値的な方法で計算を実行するプログラムを並列化しようとすると（ニューラルネットワークをトレーニングし、流体または分子システムのダイナミクスをシミュレートし、通常の微分方程式または最適化の問題を解決します）、結果には（いずれの場合も）エラーが発生します。したがって、数学的なソリューションに小さな追加のエラーを導入する可能性のある「危険な」並列化テクノロジを適用してみませんか。しかし、あなたはもう少し加速を得ることができますか？そのようなテクノロジーの1つ-中間結果を予測するループ本体の分割と、予測が失敗した場合のロールバック（実際、これは部分的にトランザクションメモリでの「過度に楽観的な」計算です）について説明します。
並列化のアイデア
本体が2つの連続した部分で構成され、2番目の部分が最初の部分に依存するサイクルがあるとします。サイクルの個々のループを相互に依存させます。例えば：
for (int i = 0; i < N; i++) {
x = f(y);
y = g(x);
}
一見、このようなループは並列化できません。ただし、試してみます。ループ本体の1番目と2番目の演算子を並行して実行してみましょう。問題は、g（x）xを計算するときに既知である必要がありますが、最初の部分の終わりでのみ計算されるということです。さて、2番目の部分の始めにxの新しい値を予測しようとするいくつかのスキームを紹介しましょう。これは、たとえば、変化の「履歴」に基づいてxの新しい値を予測することを「学習」する線形予測の助けを借りて行うことができます。次に、2番目の部分を最初の部分と並行して読み取ることができ（これは「過度の楽観主義」です）、両方を計算するときに、xの予測値を最初の部分の最後で得られた実際の値と比較します。それらがほぼ等しい場合、両方の部分の計算結果を受け入れることができます（そしてループの次の反復に進みます）。そして、それらが非常に異なる場合は、2番目の部分だけを再カウントする必要があります。このスキームでは、場合によっては純粋な並列化が得られ、残りの場合は実際のシーケンシャルカウントが得られます。ループ実行アルゴリズムは次のとおりです。
for (int i = 0; i < N; i++) {
{
1 – x = f(y). x;
2 – x* y* = g(x*). x x*. , y = y* . , : y = g(x).
}
}
基本的なアルゴリズムは明確です。理論上の加速は2倍ですが、実際には、もちろん、次の理由により小さくなります。a）時間の一部が予測と調整に費やされている。b）すべての反復が並行して実行されるわけではありません。c）サイクル本体の最初の部分と2番目の部分は、異なる労働強度を持つことができます（理想的には、等しい必要があります）。実装に移りましょう。
並列化の実装-過度に楽観的なコンピューティング
並列化アルゴリズムは、計算の一部をキャンセルして（失敗した場合）それらを再実行することを処理するため、トランザクションメモリで作業するという考えから明らかに何かがあります。より良い-部分的にトランザクションの場合、特定の変数はトランザクションメモリスキームに従って機能し、残りの変数は通常どおりに機能します。最初の部分から2番目の部分へのデータの転送は、特別なチャネルを使用して整理できます。このチャネルを予測可能にします。a）受信時にデータがすでにチャネルに送信されている場合はチャネルから読み取られます。b）受信時にデータがまだチャネルに到着していない場合は、このデータを予測しようとし、予測結果を返します。このチャネルは、従来のトランザクションメモリではやや珍しいスキームに従って機能します。サイクルの2番目の部分のトランザクションの終了時に、チャネルに受信されたデータとそれによって予測されたデータの間に不一致がある場合、サイクルの2番目の部分のトランザクションはキャンセルされて再実行され、「予測」ではなくチャネルから読み取られますが、データは実際に受信されます。サイクルは次のようになります。
for (int i = 0; i < N; i++) {
, {
1 ( 1):
x = f(y);
_.put(x);
2 ( 2):
_.get(x);
y = g(x);
}
}
完了。チャネルはデータ予測を処理し、トランザクションメモリは、過度に楽観的な予測の場合に計算をキャンセルする処理を部分的に処理しました。
いくつかの有用な用途：ニューラルネットワーク、パーティクルインセル法
ループ本体を並列化するためのこのようなスキームは、たとえば、セル内粒子法を使用して静電レンズをモデル化するとき、および逆伝播法を使用してフィードフォワードニューラルネットワークをトレーニングするときに、多くの数学的アルゴリズムで使用できます。最初のタスクは非常に特殊なので、ここでは説明しません。説明した並列化のアプローチでは、10〜15％の加速が得られたとだけ言います。しかし、2番目のタスクはすでに人気があるので、それについていくつかの言葉を言う必要があります。
ニューラルネットワークのトレーニングサイクルには、トレーニングペアのシーケンシャルパスが含まれ、各ペアに対して、順方向移動（ネットワーク出力の計算）と逆方向移動（重みとバイアスの修正）が実行されます。これらはトレーニングペアのループ本体の2つの部分であり、それらを並列化するために、上記のアプローチを適用できます（ちなみに、マイナーな変更を加えて、トレーニングペアの並列パスで適用することもできます）。その結果、一般的なニューラルネットワークトレーニングタスクで、パフォーマンスが50％向上しました。
Cプログラムの並列化の自動化
超楽観的な計算のアイデアはそれほど難しいことではないので、自動並列化を扱う特別なトランスレータプログラムが作成されました-そのような並列化が肯定的な結果をもたらすことができる元のCプログラム内のループを見つけ、それらの本体を2つの部分に分割し、必要なOpenMPディレクティブを挿入し、チャネルの潜在的な変数。部分的なトランザクションメモリと予測チャネルを操作するためのライブラリを接続し、最終的には出力並列化プログラムを生成します。
特に、このようなトランスレータは、静電レンズシミュレーションプログラムに適用されました。元のプログラム（ループの並列化を示すディレクティブを含む）と翻訳後に取得したプログラムの両方を引用します。
オリジナルプログラム：
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#pragma auto parallelize
#pragma auto pure(malloc,fabs,free,sizeof,omp_get_wtime)
#define theta 1.83
#define NX 40
#define NY 40
#define h 0.1
#define NP 15000
//
#define U1 200
#define U2 5000
#define e -1.5E-13
#define m 1E-11
#define e0 8.85E-12
#define V (h*h)
#define tau 0.000015
#define T 0.09
#define POISSON_EPS 0.01
#define TOL_EPS 0.25
int main() {
double * U = (double *)malloc(NY*NX*sizeof(double));
double * UU = (double *)malloc(NY*NX*sizeof(double));
double * EX = (double *)malloc(NY*NX*sizeof(double));
double * EY = (double *)malloc(NY*NX*sizeof(double));
double * PX = (double *)malloc(NP*sizeof(double));
double * PY = (double *)malloc(NP*sizeof(double));
int * X = (int *)malloc(NP*sizeof(int));
int * Y = (int *)malloc(NP*sizeof(int));
double ro[NY][NX];
split_private double t;
split_private double tm;
split_private int i, j;
for (i = 0; i < NY; i++)
for (j = 0; j < NX; j++) {
UU[i*NX+j] = j == NX-1 ? U2 : j == NX/2 && (i < NY/4 || i > 3*NY/4) ? U1 : 0.0;
EX[i*NX+j] = 0.0;
EY[i*NX+j] = 0.0;
}
for (i = 0; i < NP; i++) {
int x, y;
PX[i] = 0.5*NX*h*rand()/RAND_MAX;
PY[i] = NY*h*rand()/RAND_MAX;
x = PX[i]/h;
y = PY[i]/h;
if (x < 0) x = 0;
else if (x > NX-1) x = NX-1;
if (y < 0) y = 0;
else if (y > NY-1) y = NY-1;
X[i] = x;
Y[i] = y;
}
tm = omp_get_wtime();
for (t = 0.0; t < T; t += tau) {
unsigned int n[NY][NX] = { 0 };
double err;
int ptr = 0;
for (i = 0; i < NY; i++)
for (j = 0; j < NX; j++, ptr++)
U[ptr] = UU[ptr];
for (i = 1; i < NY - 1; i++)
for (j = 1; j < NX - 1; j++) {
EX[i*NX+j] = -(U[i*NX+j+1]-U[i*NX+j-1])/2.0/h;
EY[i*NX+j] = -(U[(i+1)*NX+j]-U[(i-1)*NX+j])/2.0/h;
}
for (i = 0; i < NP; i++) {
PX[i] += tau*e*EX[Y[i]*NX+X[i]]/m;
PY[i] += tau*e*EY[Y[i]*NX+X[i]]/m;
}
for (i = 0; i < NP; i++) {
int x = PX[i]/h;
int y = PY[i]/h;
if (x < 0) x = 0;
else if (x > NX-1) x = NX-1;
if (y < 0) y = 0;
else if (y > NY-1) y = NY-1;
Y[i] = y;
X[i] = x;
n[y][x]++;
}
for (i = 0; i < NY; i++)
for (j = 0; j < NX; j++)
ro[i][j] = n[i][j]*e/V;
do {
err = 0.0;
for (i = 1; i < NY - 1; i++)
for (j = 1+(i-1)%2; j < NX - 1; j+=2) {
int ptr = i*NX + j;
if (!(j == NX/2 && (i < NY/4 || i > 3*NY/4))) {
double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);
double loc_err = fabs(UU[ptr] - _new);
if (loc_err > err) err = loc_err;
UU[ptr] = _new;
}
}
for (i = 1; i < NY - 1; i++)
for (j = 1+i%2; j < NX - 1; j+=2) {
int ptr = i*NX + j;
if (!(j == NX/2 && (i < NY/4 || i > 3*NY/4))) {
double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);
double loc_err = fabs(UU[ptr] - _new);
if (loc_err > err) err = loc_err;
UU[ptr] = _new;
}
}
for (j = 0; j < NX; j++) {
UU[j] = UU[NX + j];
UU[(NY-1)*NX + j] = UU[(NY-2)*NX + j];
}
} while (err > POISSON_EPS);
}
for (i = 0; i < NY; i++) {
for (j = 0; j < NX; j++)
printf("%lf\t", UU[i*NX+j]);
printf("\n");
}
return 0;
}
自動的に並列プログラム
#include "transact.h"
#define split_private /* split-private */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define theta 1.83
#define NX 40
#define NY 40
#define h 0.1
#define NP 15000
#define U1 200
#define U2 5000
#define e -1.5E-13
#define m 1E-11
#define e0 8.85E-12
#define V (h*h)
#define tau 0.000015
#define T 0.09
#define POISSON_EPS 0.01
#define TOL_EPS 0.25
int main( ){
double * U = (double *)malloc(NY*NX*sizeof(double));
double * UU = (double *)malloc(NY*NX*sizeof(double));
double * EX = (double *)malloc(NY*NX*sizeof(double));
double * EY = (double *)malloc(NY*NX*sizeof(double));
double * PX = (double *)malloc(NP*sizeof(double));
double * PY = (double *)malloc(NP*sizeof(double));
int * X = (int *)malloc(NP*sizeof(int));
int * Y = (int *)malloc(NP*sizeof(int));
double ro[NY][NX];
split_private double t;
split_private double tm;
split_private int i, j;
for ( i = 0; i < NY; i++ )
for ( j = 0; j < NX; j++ )
{
UU[i*NX+j] = j == NX-1 ? U2 : j == NX/2 && (i < NY/4 || i > 3*NY/4) ? U1 : 0.0;
EX[i*NX+j] = 0.0;
EY[i*NX+j] = 0.0;
}
for ( i = 0; i < NP; i++ )
{
int x, y;
PX[i] = 0.5*NX*h*rand()/RAND_MAX;
PY[i] = NY*h*rand()/RAND_MAX;
x = PX[i]/h;
y = PY[i]/h;
if ( x < 0 )
x = 0;
else
if ( x > NX-1 )
x = NX-1;
if ( y < 0 )
y = 0;
else
if ( y > NY-1 )
y = NY-1;
X[i] = x;
Y[i] = y;
}
tm = omp_get_wtime();
#pragma omp parallel num_threads(2) private(t,tm,i,j)
{
int __id__ = omp_get_thread_num();
TOut<double > * out_ro = __id__ == 0 ? new TOut<double >("ro63", (NY)*(NX), 2, 0.01, -1, "63") : NULL;
TIn<double > * in_ro = __id__ == 1 ? new TIn<double >("ro63", (NY)*(NX), 2, 0.01, -1, "63") : NULL;
for ( t = 0.0; t < T; t += tau )
{
unsigned int n[NY][NX] = { 0 };
double err;
int ptr = 0;
if ( __id__ == 0 )
{
for ( i = 0; i < NY; i++ )
for ( j = 0; j < NX; j++, ptr++ )
U[ptr] = UU[ptr];
}
transaction_atomic("63")
{
if ( __id__ == 0 )
{
for ( i = 1; i < NY - 1; i++ )
for ( j = 1; j < NX - 1; j++ )
{
EX[i*NX+j] = -(U[i*NX+j+1]-U[i*NX+j-1])/2.0/h;
EY[i*NX+j] = -(U[(i+1)*NX+j]-U[(i-1)*NX+j])/2.0/h;
}
for ( i = 0; i < NP; i++ )
{
PX[i] += tau*e*EX[Y[i]*NX+X[i]]/m;
PY[i] += tau*e*EY[Y[i]*NX+X[i]]/m;
}
for ( i = 0; i < NP; i++ )
{
int x = PX[i]/h;
int y = PY[i]/h;
if ( x < 0 )
x = 0;
else
if ( x > NX-1 )
x = NX-1;
if ( y < 0 )
y = 0;
else
if ( y > NY-1 )
y = NY-1;
Y[i] = y;
X[i] = x;
n[y][x]++;
}
for ( i = 0; i < NY; i++ )
for ( j = 0; j < NX; j++ )
ro[i][j] = n[i][j]*e/V;
out_ro->put((double *)ro);
}
else
{
double ro[NY][NX];
in_ro->get((double *)ro, 0);
do
{
err = 0.0;
for ( i = 1; i < NY - 1; i++ )
for ( j = 1+(i-1)%2; j < NX - 1; j+=2 )
{
int ptr = i*NX + j;
if ( !(j == NX/2 && (i < NY/4 || i > 3*NY/4)) )
{
double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);
double loc_err = fabs(UU[ptr] - _new);
if ( loc_err > err )
err = loc_err;
UU[ptr] = _new;
}
}
for ( i = 1; i < NY - 1; i++ )
for ( j = 1+i%2; j < NX - 1; j+=2 )
{
int ptr = i*NX + j;
if ( !(j == NX/2 && (i < NY/4 || i > 3*NY/4)) )
{
double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);
double loc_err = fabs(UU[ptr] - _new);
if ( loc_err > err )
err = loc_err;
UU[ptr] = _new;
}
}
for ( j = 0; j < NX; j++ )
{
UU[j] = UU[NX + j];
UU[(NY-1)*NX + j] = UU[(NY-2)*NX + j];
}
}
while ( err > POISSON_EPS )
;
}
}
}
delete in_ro;
delete out_ro;
}
for ( i = 0; i < NY; i++ )
{
for ( j = 0; j < NX; j++ )
printf("%lf\t", UU[i*NX+j]);
printf("\n");
}
return 0;
}
結果
そのため、厳密にシーケンシャルなフラグメントで構成されている場合でも、プログラムの並列化を試みることができ、スピードアップで肯定的な結果を得ることができます（私の実験では、スピードアップは15％から50％に増加します）。この短い記事が誰かに役立つことを願っています。