Um pouco sobre como acelerar o programa: paralelização (manual ou automática) baseada em cálculos superotimistas

Olá queridos leitores. Nesta publicação, vamos falar sobre algo (já conhecido) como acelerar o programa usando cálculos paralelos. As tecnologias para organizar tais cálculos são conhecidas - isto é tanto programação multithread comum quanto o uso de interfaces especiais: OpenMP, OpenAcc, MPI, DVM e muitos outros (neste caso, loops são paralelizados, vetorização ou pipelining são usados, computações preguiçosas são organizadas, blocos de programas independentes são alocados para serem executados em paralelo, etc.).



Nesse caso, geralmente partem da ideia de que a paralelização não deve afetar de alguma forma os resultados da execução do programa. Este é um requisito difícil, mas justo em muitos casos. No entanto, se tentarmos paralelizar um programa que realiza qualquer cálculo por métodos numéricos (treinamos uma rede neural, simulamos a dinâmica de um fluido ou sistema molecular, resolvemos equações diferenciais ordinárias ou problemas de otimização), então o resultado (em qualquer caso) terá algum erro. Portanto, por que não aplicar tecnologias de paralelização "arriscadas", que podem introduzir um pequeno erro adicional na solução matemática,mas permite que você obtenha mais aceleração adicional? Uma dessas tecnologias - a divisão de corpos de loop com previsão de resultados intermediários e reversão em caso de predição malsucedida (na verdade, isso é cálculos "superotimistas" em memória parcialmente transacional) será discutida.



Ideia de paralelização



Suponha que temos um ciclo cujo corpo consiste em duas partes consecutivas e a segunda parte depende da primeira. Deixe os loops individuais do ciclo dependerem uns dos outros. Por exemplo:



for (int i = 0; i < N; i++) {
	x = f(y);
	y = g(x);
}


À primeira vista, tal loop não pode ser paralelizado. No entanto, vamos tentar. Vamos tentar executar o primeiro e o segundo operadores do corpo do loop em paralelo. O problema é que na hora de calcular g (x) x deve ser conhecido, mas só será calculado no final da primeira parte. Bem, vamos introduzir algum esquema, que no início da segunda parte tentará prever o novo valor de x. Isso pode ser feito, por exemplo, usando preditiva linear, que "aprende" a prever um novo valor de x, com base no "histórico" de sua mudança. Então, a segunda parte pode ser lida em paralelo com a primeira (isso é "superotimismo"), e quando ambos são calculados, compare o valor previsto de x com o real obtido no final da primeira parte. Se eles forem aproximadamente iguais, o resultado dos cálculos de ambas as partes pode ser aceito (e ir para a próxima iteração do loop).E se eles forem muito diferentes, então apenas a segunda parte precisará ser recontada. Com este esquema, em alguma parte dos casos, obteremos paralelização pura, no resto - a contagem sequencial real. O algoritmo de execução do loop é o seguinte:



for (int i = 0; i < N; i++) {
	    {
		  1 –  x = f(y).        x;
		  2 –   x*   y* = g(x*).   x        x*.   ,  y = y*    .   ,     : y = g(x). 
	}
}


O algoritmo básico é claro. A aceleração teórica é dupla, mas na prática será, claro, menor, porque: a) parte do tempo é gasto em previsões e coordenação; b) nem todas as iterações serão executadas em paralelo; c) a primeira e a segunda partes do corpo do ciclo podem ter diferentes intensidades de trabalho (idealmente, igual é necessário). Vamos prosseguir para a implementação.



Implementação de paralelização - computação super otimista



Uma vez que o algoritmo de paralelização trata do cancelamento de alguns cálculos (em caso de falha) e da sua reexecução, há algo claro na ideia de trabalhar na memória transacional. Melhor - em um parcialmente transacional, onde certas variáveis ​​funcionam de acordo com o esquema de memória transacional, e o resto das variáveis ​​funcionam normalmente. A transferência de dados da primeira parte para a segunda pode ser organizada usando algum canal especial. Deixe este canal ser preditivo: a) se no momento da recepção os dados já foram transmitidos ao canal, então eles são lidos a partir dele, b) se no momento da recepção os dados ainda não chegaram ao canal, então ele tenta prever esses dados e retorna o resultado da previsão. Este canal funcionará de acordo com um esquema um tanto incomum para a memória transacional convencional:Se no final da transação da segunda parte do ciclo houver uma discrepância entre os dados recebidos no canal e os dados previstos por ele, então a transação da segunda parte do ciclo é cancelada e reexecutada, e não "predições" serão lidas do canal, mas os dados realmente recebidos. O ciclo será semelhante a:



for (int i = 0; i < N; i++) {
	   ,     {
		 1 ( 1):
			x = f(y);
			_.put(x);
		 2 ( 2):
			_.get(x);
			y = g(x);
	}
}


Feito. O canal cuidava da previsão dos dados, enquanto a memória transacional parcialmente cuidava do cancelamento dos cálculos no caso de uma previsão excessivamente otimista.



Alguns usos úteis: redes neurais, método de partícula na célula



Tal esquema para paralelizar o corpo do loop pode ser usado em uma série de algoritmos matemáticos, por exemplo, ao modelar uma lente eletrostática usando o método de partícula na célula, bem como ao treinar uma rede neural feedforward usando o método de retropropagação. A primeira tarefa é muito especial, portanto não a discutirei aqui, apenas direi que a abordagem descrita para a paralelização deu uma aceleração de 10-15%. Mas a segunda tarefa já é mais popular, portanto, basta dizer algumas palavras a respeito.



O ciclo de treinamento da rede neural inclui uma passagem sequencial pelos pares de treinamento, e para cada par, um movimento para frente (cálculo da saída da rede) e um movimento reverso (correção de pesos e vieses) são executados. Essas são as duas partes do corpo do loop para pares de treinamento e, para paralelizá-las, você pode aplicar a abordagem acima (aliás, também pode ser aplicada com uma passagem paralela por pares de treinamento, com pequenas alterações). Como resultado, em uma tarefa típica de treinamento de rede neural, obtive um ganho de 50% no desempenho.



Automação de paralelização de programas C



A ideia de cálculos superotimistas não é muito difícil, então um programa tradutor especial foi escrito que lida com paralelização automática - ele encontra loops no programa C original para os quais tal paralelização pode dar um resultado positivo e divide seus corpos em duas partes, inserindo as diretivas OpenMP necessárias, encontrando variáveis ​​potenciais para canais, conectando uma biblioteca para trabalhar com memória parcialmente transacional e canais preditivos e, em última instância, gerar um programa paralelizado de saída.



Em particular, tal tradutor foi aplicado a um programa de simulação de lentes eletrostáticas. Darei os dois programas - o original (que inclui uma diretiva que indica a paralelização de loops) e o obtido após a tradução.



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 automaticamente 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;
}


Resultado



Então, às vezes você pode tentar paralelizar o programa mesmo nos casos em que ele consiste em fragmentos estritamente sequenciais, e até obter resultados positivos na aceleração (em meus experimentos, a aceleração é aumentada de 15 para 50%). Espero que este breve artigo seja útil para alguém.



All Articles