Dans ce cas, ils partent généralement de l'idée que la parallélisation ne devrait pas affecter d'une manière ou d'une autre les résultats de l'exécution du programme. C'est une exigence difficile, mais juste dans de nombreux cas. Cependant, si nous essayons de paralléliser un programme qui effectue certains calculs par des méthodes numériques (nous entraînons un réseau de neurones, simulons la dynamique d'un fluide ou d'un système moléculaire, résolvons des équations différentielles ordinaires ou des problèmes d'optimisation), alors le résultat aura (dans tous les cas) une erreur. Par conséquent, pourquoi ne pas appliquer des technologies de parallélisation «risquées», qui peuvent introduire une petite erreur supplémentaire dans la solution mathématique,mais vous permettre d'obtenir une accélération supplémentaire supplémentaire? Une de ces technologies - le fractionnement des corps de boucle avec prédiction des résultats intermédiaires et retour en arrière en cas de prédiction infructueuse (en fait, il s'agit de calculs "trop optimistes" dans une mémoire partiellement transactionnelle) sera discutée.
Idée de parallélisation
Supposons que nous ayons un cycle dont le corps se compose de deux parties consécutives, et la seconde partie dépend de la première. Laissez les différentes boucles du cycle dépendre les unes des autres. Par exemple:
for (int i = 0; i < N; i++) {
x = f(y);
y = g(x);
}
À première vue, une telle boucle ne peut pas être parallélisée. Cependant, nous essaierons. Essayons d'exécuter les premier et deuxième opérateurs du corps de la boucle en parallèle. Le problème est qu'au moment du calcul g (x) x doit être connu, mais il ne sera calculé qu'à la fin de la première partie. Eh bien, introduisons un schéma qui, au début de la deuxième partie, tentera de prédire la nouvelle valeur de x. Cela peut être fait, par exemple, à l'aide du prédictif linéaire, qui "apprend" à prédire une nouvelle valeur de x, en fonction de "l'historique" de son changement. Ensuite, la deuxième partie peut être lue en parallèle avec la première (c'est un "sur-optimisme"), et lorsque les deux sont calculées, comparez la valeur prédite de x avec la valeur réelle obtenue à la fin de la première partie. S'ils sont approximativement égaux, alors le résultat des calculs des deux parties peut être accepté (et passer à l'itération suivante de la boucle).Et s'ils sont très différents, alors seule la deuxième partie devra être racontée. Avec ce schéma, dans certaines parties des cas, nous obtiendrons une parallélisation pure, dans le reste - le décompte séquentiel réel. L'algorithme d'exécution de la boucle est le suivant:
for (int i = 0; i < N; i++) {
{
1 – x = f(y). x;
2 – x* y* = g(x*). x x*. , y = y* . , : y = g(x).
}
}
L'algorithme de base est clair. L'accélération théorique est le double, mais en pratique, elle sera bien entendu moindre, car: a) une partie du temps est consacrée aux prévisions et à la coordination; b) toutes les itérations ne seront pas exécutées en parallèle; c) les première et deuxième parties du corps du cycle peuvent avoir une intensité de travail différente (idéalement, une égalité est requise). Passons à la mise en œuvre.
Mise en œuvre de la parallélisation - Calcul trop optimiste
Puisque l'algorithme de parallélisation traite d'annuler certains des calculs (en cas d'échec) et de les réexécuter, il y a clairement quelque chose dans l'idée de travailler en mémoire transactionnelle. Mieux - dans une situation partiellement transactionnelle, où certaines variables fonctionnent selon le schéma de mémoire transactionnelle, et le reste des variables fonctionnent comme d'habitude. Le transfert de données de la première partie à la seconde peut être organisé en utilisant un canal spécial. Que ce canal soit prédictif: a) si au moment de la réception les données ont déjà été transmises au canal, alors elles sont lues à partir de celui-ci, b) si au moment de la réception les données ne sont pas encore arrivées dans le canal, alors il essaie de prédire ces données et renvoie le résultat de la prédiction. Ce canal fonctionnera selon un schéma quelque peu inhabituel pour la mémoire transactionnelle conventionnelle:Si à la fin de la transaction de la deuxième partie du cycle il y a un écart entre les données reçues dans le canal et les données prédites par celui-ci, alors la transaction de la deuxième partie du cycle est annulée et réexécutée, et non les «prédictions» seront lues à partir du canal, mais les données effectivement reçues. Le cycle ressemblera à:
for (int i = 0; i < N; i++) {
, {
1 ( 1):
x = f(y);
_.put(x);
2 ( 2):
_.get(x);
y = g(x);
}
}
Terminé. Le canal s'est occupé de la prédiction des données, tandis que la mémoire de transaction s'est occupée de l'annulation des calculs en cas de prédiction trop optimiste.
Quelques applications utiles: réseaux de neurones, méthode particule-dans-cellule
Un tel schéma de parallélisation du corps de boucle peut être utilisé dans un certain nombre d'algorithmes mathématiques, par exemple, lors de la modélisation d'une lentille électrostatique en utilisant la méthode des particules dans la cellule, ainsi que lors de la formation d'un réseau de neurones à anticipation en utilisant la méthode de rétropropagation. La première tâche est très spéciale, donc je ne la discuterai pas ici, je dirai seulement que l'approche décrite de la parallélisation a donné une accélération de 10-15%. Mais la deuxième tâche est déjà plus populaire, il faut donc simplement en dire quelques mots.
Le cycle d'apprentissage du réseau neuronal comprend un passage séquentiel à travers les paires d'apprentissage, et pour chaque paire, un mouvement vers l'avant (calcul de la sortie du réseau) et un mouvement inverse (correction des poids et des biais) sont effectués. Ce sont les deux parties du corps de la boucle pour les paires d'entraînement, et pour les paralléliser, vous pouvez appliquer l'approche ci-dessus (en passant, elle peut également être appliquée avec un passage parallèle à travers des paires d'entraînement, avec des changements mineurs). En conséquence, sur une tâche d'entraînement de réseau neuronal typique, j'ai obtenu un gain de performances de 50%.
Automatisation de la parallélisation des programmes C
L'idée de calculs super-optimistes n'est pas très difficile, donc un programme de traduction spécial a été écrit qui traite de la parallélisation automatique - il trouve des boucles dans le programme C original pour lequel une telle parallélisation peut donner un résultat positif et divise leurs corps en deux parties, insérant les directives OpenMP nécessaires, trouvant variables potentielles pour les canaux, connexion d'une bibliothèque pour travailler avec une mémoire transactionnelle partielle et des canaux prédictifs et, finalement, générer un programme parallélisé de sortie.
En particulier, un tel traducteur a été appliqué à un programme de simulation de lentilles électrostatiques. Je citerai les deux programmes - celui d'origine (qui comprend la directive indiquant la parallélisation des boucles) et celui obtenu après traduction.
Programme original:
#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;
}
Programme automatiquement parallèle
#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;
}
Résultat
Ainsi, parfois, vous pouvez essayer de paralléliser le programme même dans les cas où il se compose de fragments strictement séquentiels, et même obtenir des résultats positifs en accélérant (dans mes expériences, l'accélération est augmentée de 15 à 50%). J'espère que ce court article sera utile à quelqu'un.