关于加速程序的一些知识:基于超乐观计算的并行化(手动或自动)

尊敬的读者您好。在本出版物中,我们将谈论(已经熟悉)通过使用并行计算来加速程序的事情。组织这种计算的技术是众所周知的-这既是普通的多线程编程,又是特殊接口的使用:OpenMP,OpenAcc,MPI,DVM和许多其他接口(在这种情况下,对循环进行并行化,使用矢量化或流水线,组织惰性计算,分配独立的程序块并可以运行并行等)。



在这种情况下,它们通常是从这样的思想出发的:并行化不应以某种方式影响程序执行的结果。这是一个艰巨的要求,但在许多情况下是公平的。但是,如果我们尝试并行化通过数值方法执行任何计算的程序(我们训练神经网络,模拟流体或分子系统的动力学,解决常微分方程或优化问题),那么结果(无论如何)都会有一些误差。因此,为什么不应用“危险的”并行化技术,该技术会在数学解决方案中引入小的附加误差,但是可以让您获得更多的额外加速吗?将讨论这样的技术之一-在预测不成功的情况下使用循环中间体进行中间结果的预测和回滚(实际上,这是部分事务性存储中的“过分乐观”的计算)。



并行化思想



假设我们有一个循环,循环的主体由两个连续的部分组成,第二个部分取决于第一个部分。让循环的各个循环相互依赖。例如:



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


乍一看,这样的循环无法并行化。但是,我们会尝试。让我们尝试并行执行循环主体的第一和第二运算符。问题在于,在计算g(x)时必须知道x,但是只能在第一部分的末尾进行计算。好吧,让我们介绍一些方案,该方案在第二部分开始时将尝试预测x的新值。例如,这可以使用线性预测来完成,线性预测基于其变化的“历史”来“学习”以预测x的新值。然后,可以与第一部分并行读取第二部分(这是“过分乐观”),并且在计算第二部分时,将x的预测值与在第一部分末尾获得的实数进行比较。如果它们近似相等,则可以接受两个部分的计算结果(并转到循环的下一个迭代)。而且,如果它们之间有很大差异,则只需重新叙述第二部分。使用这种方案,在某些情况下,我们将获得纯并行化,其余的是实际的顺序计数。循环执行算法如下:



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


基本算法很清楚。理论上的加速度是原来的两倍,但实际上,它实际上会更少,因为:a)部分时间用于预测和协调;b)并非所有迭代都将并行执行;c)自行车车体的第一部分和第二部分可以具有不同的劳动强度(理想情况下,要求相等)。让我们继续执行。



并行化实现-超乐观计算



由于并行化算法负责取消一些计算(如果失败)并重新执行它们,因此显然在事务性内存中有一些想法。更好-在部分事务的事务中,某些变量根据事务内存方案工作,而其余变量照常工作。可以使用一些特殊的通道来组织从第一部分到第二部分的数据传输。假设此通道具有预测性:a)如果在接收时已经将数据传输到该通道,则从该通道读取它们; b)如果在接收时尚未将数据到达通道,则它将尝试预测此数据并返回预测结果。该通道将按照一种与常规事务存储不同寻常的方案工作:如果在周期的第二部分的事务结束时,通道中接收到的数据与通道所预测的数据之间存在差异,则取消并重新执行周期的第二部分的事务,并且不会从通道中读取“预测”,而是实际接收到的数据。该周期将如下所示:



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


做完了 通道负责数据预测,而事务存储器则部分负责取消过于乐观的预测的计算。



一些有用的用途:神经网络,单元中的粒子方法



这种用于使环体平行的方案可以应用在许多数学算法中,例如,当通过单元内粒子方法对静电透镜建模时,以及在使用反向传播方法训练前馈神经网络时。第一项任务非常特殊,因此在这里我将不讨论它,我只会说所描述的并行化方法可以加快10-15%。但是第二个任务已经很流行了,因此只需说几句话即可。



神经网络的训练周期包括对训练对的顺序遍历,并且对于每对训练,将进行正向移动(网络输出的计算)和反向移动(权重和偏差的校正)。这是训练对循环主体的两个部分,要使它们并行化,您可以应用上述方法(顺便说一下,也可以将训练对以并行方式应用,但要稍作改动)。结果,在典型的神经网络训练任务中,我的成绩提高了50%。



自动化C程序的并行化



超乐观计算的想法不是很困难,因此编写了一个特殊的转换器程序来处理自动并行化-它在原始C程序中找到循环,这样的并行化可以给出肯定的结果并将其主体分成两部分,插入必要的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%)。我希望这篇简短的文章对某人有用。



All Articles