En este caso, generalmente parten de la idea de que la paralelización no debería afectar de alguna manera los resultados de la ejecución del programa. Este es un requisito difícil, pero justo en muchos casos. Sin embargo, si intentamos paralelizar un programa que realiza cualquier cálculo por métodos numéricos (entrenamos una red neuronal, simulamos la dinámica de un fluido o un sistema molecular, resolvemos ecuaciones diferenciales ordinarias o problemas de optimización), entonces el resultado (en cualquier caso) tendrá algún error. Por tanto, ¿por qué no aplicar tecnologías de paralelización "riesgosas", que pueden introducir un pequeño error adicional en la solución matemática?pero ¿le permite obtener más aceleración adicional? Se discutirá una de esas tecnologías: la división de cuerpos de bucle con predicción de resultados intermedios y retroceso en caso de predicción fallida (de hecho, se trata de cálculos "demasiado optimistas" en la memoria transaccional parcial).
Idea de paralelización
Supongamos que tenemos un ciclo, cuyo cuerpo consta de dos partes consecutivas y la segunda parte depende de la primera. Deje que los bucles individuales del ciclo dependan unos de otros. Por ejemplo:
for (int i = 0; i < N; i++) {
x = f(y);
y = g(x);
}
A primera vista, dicho bucle no se puede paralelizar. Sin embargo, lo intentaremos. Intentemos ejecutar el primer y segundo operador del cuerpo del bucle en paralelo. El problema es que a la hora de calcular g (x) x se debe conocer, pero solo se calculará al final de la primera parte. Bien, introduzcamos algún esquema, que al comienzo de la segunda parte intentará predecir el nuevo valor de x. Esto se puede hacer, por ejemplo, utilizando predictivo lineal, que "aprende" a predecir un nuevo valor de x, basándose en el "historial" de su cambio. Entonces la segunda parte se puede leer en paralelo con la primera (esto es "exceso de optimismo"), y cuando se calculan ambos, comparar el valor predicho de x con el real obtenido al final de la primera parte. Si son aproximadamente iguales, entonces se puede aceptar el resultado de los cálculos de ambas partes (y pasar a la siguiente iteración del ciclo).Y si son muy diferentes, solo será necesario relatar la segunda parte. Con este esquema, en alguna parte de los casos, obtendremos una paralelización pura, en el resto, el conteo secuencial real. El algoritmo de ejecución del bucle es el siguiente:
for (int i = 0; i < N; i++) {
{
1 – x = f(y). x;
2 – x* y* = g(x*). x x*. , y = y* . , : y = g(x).
}
}
El algoritmo básico es claro. La aceleración teórica es dos veces, pero en la práctica será, por supuesto, menor, porque: a) parte del tiempo se dedica a predicciones y coordinación; b) no todas las iteraciones se ejecutarán en paralelo; c) la primera y segunda partes del cuerpo del ciclo pueden tener diferente intensidad de trabajo (idealmente, se requiere la misma). Pasemos a la implementación.
Implementación de la paralelización: Computación demasiado optimista
Dado que el algoritmo de paralelización se ocupa de cancelar algunos de los cálculos (en caso de falla) y volver a ejecutarlos, claramente hay algo en la idea de trabajar en la memoria transaccional. Mejor: en uno parcialmente transaccional, donde ciertas variables funcionan de acuerdo con el esquema de memoria transaccional y el resto de las variables funcionan como de costumbre. La transferencia de datos de la primera parte a la segunda se puede organizar utilizando algún canal especial. Deje que este canal sea predictivo: a) si en el momento de la recepción los datos ya se han transmitido al canal, entonces se leen de él, b) si en el momento de la recepción los datos aún no han llegado al canal, entonces intenta predecir estos datos y devuelve el resultado de la predicción. Este canal funcionará según un esquema algo inusual para la memoria transaccional convencional:Si al final de la transacción de la segunda parte del ciclo hay una discrepancia entre los datos recibidos en el canal y los datos predichos por él, entonces la transacción de la segunda parte del ciclo se cancela y se vuelve a ejecutar, y no se leerán "predicciones" del canal, sino datos realmente recibidos. El ciclo se verá así:
for (int i = 0; i < N; i++) {
, {
1 ( 1):
x = f(y);
_.put(x);
2 ( 2):
_.get(x);
y = g(x);
}
}
Hecho. El canal se encargó de la predicción de datos, mientras que la memoria transaccional se encargó parcialmente de cancelar los cálculos en caso de una predicción demasiado optimista.
Algunas aplicaciones útiles: redes neuronales, método de partículas en la celda
Un esquema de este tipo para paralelizar el cuerpo del bucle se puede utilizar en una serie de algoritmos matemáticos, por ejemplo, al modelar una lente electrostática utilizando el método de partículas en las células, así como al entrenar una red neuronal de retroalimentación utilizando el método de retropropagación. La primera tarea es muy especial, por lo que no la discutiré aquí, solo diré que el enfoque descrito para la paralelización dio una aceleración del 10-15%. Pero la segunda tarea ya es más popular, por lo que simplemente es necesario decir algunas palabras al respecto.
El ciclo de entrenamiento de la red neuronal incluye un pase secuencial a través de los pares de entrenamiento, y para cada par, se realizan un movimiento hacia adelante (cálculo de la salida de la red) y un movimiento inverso (corrección de pesos y sesgos). Estas son las dos partes del cuerpo del bucle para los pares de entrenamiento, y para paralelizarlos, puede aplicar el enfoque anterior (por cierto, también se puede aplicar con un pase paralelo a través de pares de entrenamiento, con cambios menores). Como resultado, en una tarea típica de entrenamiento de redes neuronales, obtuve una ganancia del 50% en el rendimiento.
Automatización de la paralelización de programas C
La idea de cálculos súper optimistas no es muy difícil, por lo que se escribió un programa traductor especial que se ocupa de la paralelización automática: encuentra bucles en el programa C original para los cuales dicha paralelización puede dar un resultado positivo y divide sus cuerpos en dos partes, insertando las directivas OpenMP necesarias, encontrando posibles variables para canales, conectar una librería para trabajar con memoria transaccional parcial y canales predictivos y, en última instancia, generar un programa de salida en paralelo.
En particular, dicho traductor se aplicó a un programa de simulación de lentes electrostáticos. Daré ambos programas: el original (que incluye una directiva que indica la paralelización de los bucles) y el obtenido después de la traducción.
Programa 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;
}
Programa automáticamente paralelo
#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;
}
Salir
Entonces, a veces puede intentar paralelizar el programa incluso en los casos en que consta de fragmentos estrictamente secuenciales, e incluso obtener resultados positivos en la aceleración (en mis experimentos, la aceleración aumenta del 15 al 50%). Espero que este breve artículo sea de utilidad para alguien.