Ein wenig über die Beschleunigung des Programms: Parallelisierung (manuell oder automatisch) basierend auf superoptimistischen Berechnungen

Hallo liebe Leser. In dieser Veröffentlichung werden wir über solche (bereits bekannten) Dinge sprechen, wie das Programm durch parallele Berechnungen zu beschleunigen. Die Technologien zum Organisieren solcher Berechnungen sind bekannt - dies ist sowohl gewöhnliche Multithread-Programmierung als auch die Verwendung spezieller Schnittstellen: OpenMP, OpenAcc, MPI, DVM und viele andere (in diesem Fall werden Schleifen parallelisiert, Vektorisierung oder Pipelining verwendet, verzögerte Berechnungen werden organisiert, unabhängige Programmblöcke werden zugewiesen, die ausgeführt werden können parallel usw.).



In diesem Fall gehen sie normalerweise davon aus, dass die Parallelisierung die Ergebnisse der Programmausführung nicht beeinflussen sollte. Dies ist eine schwierige Anforderung, aber in vielen Fällen fair. Wenn wir jedoch versuchen, ein Programm zu parallelisieren, das Berechnungen mit numerischen Methoden durchführt (wir trainieren ein neuronales Netzwerk, simulieren die Dynamik eines Fluids oder eines molekularen Systems, lösen gewöhnliche Differentialgleichungen oder Optimierungsprobleme), weist das Ergebnis (in jedem Fall) einen Fehler auf. Wenden Sie daher „riskante“ Parallelisierungstechnologien an, die einen kleinen zusätzlichen Fehler in der mathematischen Lösung verursachen können.aber erlauben Sie Ihnen, etwas mehr Beschleunigung zu bekommen? Eine dieser Technologien - die Aufteilung von Schleifenkörpern mit der Vorhersage von Zwischenergebnissen und das Zurücksetzen bei erfolgloser Vorhersage (tatsächlich handelt es sich um "überoptimistische" Berechnungen im teilweise Transaktionsspeicher) wird diskutiert.



Parallelisierungsidee



Angenommen, wir haben einen Zyklus, dessen Körper aus zwei aufeinanderfolgenden Teilen besteht und dessen zweiter Teil vom ersten abhängt. Lassen Sie einzelne Schleifen des Zyklus voneinander abhängen. Zum Beispiel:



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


Auf den ersten Blick kann eine solche Schleife nicht parallelisiert werden. Wir werden es jedoch versuchen. Versuchen wir, den ersten und den zweiten Operator des Schleifenkörpers parallel auszuführen. Das Problem ist, dass zum Zeitpunkt der Berechnung g (x) x bekannt sein muss, aber erst am Ende des ersten Teils berechnet wird. Nun, lassen Sie uns ein Schema einführen, das zu Beginn des zweiten Teils versucht, den neuen Wert von x vorherzusagen. Dies kann zum Beispiel mit Hilfe der linearen Vorhersage erfolgen, die "lernt", einen neuen Wert von x basierend auf der "Geschichte" seiner Änderung vorherzusagen. Dann kann der zweite Teil parallel zum ersten gelesen werden (dies ist "Überoptimismus"), und wenn beide berechnet sind, vergleichen Sie den vorhergesagten Wert von x mit dem am Ende des ersten Teils erhaltenen realen Wert. Wenn sie ungefähr gleich sind, kann das Ergebnis der Berechnungen beider Teile akzeptiert werden (und mit der nächsten Iteration der Schleife fortfahren).Und wenn sie sehr unterschiedlich sind, muss nur der zweite Teil nachgezählt werden. Mit diesem Schema erhalten wir in einigen Fällen eine reine Parallelisierung, im Übrigen die tatsächliche sequentielle Zählung. Der Schleifenausführungsalgorithmus lautet wie folgt:



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


Der grundlegende Algorithmus ist klar. Die theoretische Beschleunigung ist zweimal, aber in der Praxis wird sie natürlich geringer sein, weil: a) ein Teil der Zeit für Vorhersagen und Koordination aufgewendet wird; b) nicht alle Iterationen werden parallel ausgeführt; c) Der erste und der zweite Teil des Zykluskörpers können unterschiedliche Arbeitsintensitäten aufweisen (idealerweise ist Gleichheit erforderlich). Fahren wir mit der Implementierung fort.



Parallelisierungsimplementierung - Überoptimistisches Computing



Da der Parallelisierungsalgorithmus einige der Berechnungen (im Fehlerfall) abbricht und erneut ausführt, liegt eindeutig etwas in der Idee, im Transaktionsspeicher zu arbeiten. Besser - in einer teilweise transaktionalen, in der bestimmte Variablen gemäß dem Transaktionsspeicherschema arbeiten und der Rest der Variablen wie gewohnt funktioniert. Die Übertragung von Daten vom ersten Teil zum zweiten Teil kann über einen speziellen Kanal organisiert werden. Lassen Sie diesen Kanal vorhersagend sein: a) Wenn zum Zeitpunkt des Empfangs die Daten bereits an den Kanal gesendet wurden, werden sie von diesem gelesen. B) Wenn zum Zeitpunkt des Empfangs die Daten noch nicht im Kanal angekommen sind, versucht er, diese Daten vorherzusagen und gibt das Vorhersageergebnis zurück. Dieser Kanal arbeitet nach einem Schema, das für herkömmliche Transaktionsspeicher etwas ungewöhnlich ist:Wenn am Ende der Transaktion des zweiten Teils des Zyklus eine Diskrepanz zwischen den in den Kanal empfangenen Daten und den von ihm vorhergesagten Daten besteht, wird die Transaktion des zweiten Teils des Zyklus abgebrochen und erneut ausgeführt, und es werden keine "Vorhersagen" aus dem Kanal gelesen, sondern die tatsächlich empfangenen Daten. Der Zyklus sieht folgendermaßen aus:



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


Erledigt. Der Kanal kümmerte sich um die Datenvorhersage, während der Transaktionsspeicher sich teilweise darum kümmerte, Berechnungen im Falle einer zu optimistischen Vorhersage abzubrechen.



Einige nützliche Anwendungen: Neuronale Netze, Partikel-in-Zelle-Methode



Ein solches Schema zur Parallelisierung des Schleifenkörpers kann in einer Reihe von mathematischen Algorithmen verwendet werden, beispielsweise beim Modellieren einer elektrostatischen Linse unter Verwendung der Partikel-in-Zelle-Methode sowie beim Trainieren eines vorwärtsgerichteten neuronalen Netzwerks unter Verwendung der Rückausbreitungsmethode. Die erste Aufgabe ist sehr speziell, daher werde ich hier nicht darauf eingehen. Ich werde nur sagen, dass der beschriebene Ansatz zur Parallelisierung eine Beschleunigung von 10-15% ergab. Aber die zweite Aufgabe ist bereits populärer, so dass es einfach notwendig ist, ein paar Worte darüber zu sagen.



Der Trainingszyklus des neuronalen Netzwerks umfasst einen sequentiellen Durchlauf durch die Trainingspaare, und für jedes Paar werden eine Vorwärtsbewegung (Berechnung der Netzwerkleistung) und eine Rückwärtsbewegung (Korrektur von Gewichten und Vorspannungen) durchgeführt. Dies sind die beiden Teile des Körpers der Schleife für Trainingspaare. Um sie zu parallelisieren, können Sie den obigen Ansatz anwenden (übrigens kann er auch mit einem parallelen Durchlauf durch Trainingspaare angewendet werden, mit geringfügigen Änderungen). Infolgedessen erhielt ich bei einer typischen Trainingsaufgabe für neuronale Netze einen Leistungszuwachs von 50%.



Automatisierung der Parallelisierung von C-Programmen



Die Idee superoptimistischer Berechnungen ist nicht sehr schwierig, daher wurde ein spezielles Übersetzerprogramm geschrieben, das sich mit automatischer Parallelisierung befasst. Es findet Schleifen im ursprünglichen C-Programm, für die eine solche Parallelisierung ein positives Ergebnis liefern kann, und teilt ihre Körper in zwei Teile auf, wobei die erforderlichen OpenMP-Direktiven eingefügt werden Mögliche Variablen für Kanäle, Verbinden einer Bibliothek zum Arbeiten mit teilweise Transaktionsspeicher und Vorhersagekanälen und letztendlich Erzeugen eines parallelisierten Ausgabeprogramms.



Insbesondere wurde ein solcher Übersetzer auf ein Simulationsprogramm für elektrostatische Linsen angewendet. Ich werde beide Programme zitieren - das ursprüngliche (das die Richtlinie enthält, die die Parallelisierung von Schleifen angibt) und das nach der Übersetzung erhaltene.



Originalprogramm:



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


Automatisch paralleles Programm



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


Ergebnis



Manchmal können Sie also versuchen, das Programm zu parallelisieren, selbst wenn es aus streng aufeinanderfolgenden Fragmenten besteht, und sogar positive Ergebnisse bei der Beschleunigung erzielen (in meinen Experimenten wird die Beschleunigung von 15 auf 50% erhöht). Ich hoffe, dieser kurze Artikel wird jemandem nützlich sein.



All Articles