使用SIMD优化1D热方程

我正在使用CFD代码（用于计算流体动态）。 我最近有机会在我的一个循环中看到英特尔编译器使用SSE，在这个循环中为计算性能增加了近2倍的因子。 但是，使用SSE和SIMD指令似乎更像是运气。 大多数情况下，编译器什么都不做。

我正在尝试强制使用SSE，考虑到AVX指令将在不久的将来强化这一方面。

我做了一个简单的1D热转印代码。 它由两个阶段组成，使用另一个阶段的结果（U0 – > U1，然后是U1 – > U0，然后是U0 – > U1等）。 当它迭代时，它会收敛到稳定的解决方案。 主代码中的大部分循环都使用相同的计算方式。 （有限差分）。

但是，我的代码比普通循环慢两倍。 结果是相同的，因此计算是一致的。

我弄错了吗？ 在使用Westmer进行测试之前，我正在使用Core 2来测试循环。

这是代码，带有SSE循环，然后是引用循环：

#include  #include  #include  //#include  #define n1 1004 #define niter 200000 int i,j,t; double U0[n1] __attribute__ ((aligned(16))); double U1[n1] __attribute__ ((aligned(16))); double Dx,Dy,Lx,Ly,InvDxDx,Dt,alpha,totaltime,Stab,DtAlpha,DxDx; __m128d vmmx00; __m128d vmmx01; __m128d vmmx02; __m128d vmmx10; __m128d va; __m128d vb; __m128d vc; __m128d vd; clock_t time0,time1; FILE *f1; int main() { /* ---- GENERAL ---- */ alpha = 0.4; totaltime = 1.0/100.0; Dt = totaltime/((niter-1)*1.0); Lx = 1.0; Dx = Lx/((n1-1)*1.0); InvDxDx = 1.0/(Dx*Dx); DxDx = Dx*Dx; Stab = alpha*Dt*(InvDxDx); DtAlpha = Dt*alpha; /* Stability if result <= 0.5 */ printf("Stability factor : %f \n",Stab); for( i = 0; i < n1; i++){U0[i] = 0.0;} U0[1] = 1.0; U0[2] = 1.0; U0[3] = 1.0; U0[n1-2] = 2.0; // for ( i = 0; i < n1; i++) { // for ( j = i + 1; j < n2; j++) { // std::swap(U0[i][j], U0[j][i]); // } //} va = _mm_set1_pd(-2.0); vb = _mm_set1_pd(InvDxDx); vd = _mm_set1_pd(DtAlpha); time0=clock(); for( t = 0; t < niter; t++) { for( i = 2; i < n1-2; i+=2) { //printf("%d %d \n",i,j); //fflush(stdout); vmmx00 = _mm_load_pd(&U0[i]); vmmx01 = _mm_loadu_pd(&U0[i+1]); vmmx02 = _mm_loadu_pd(&U0[i-1]); vmmx10 = _mm_mul_pd(va,vmmx00); // U1[i][j] = -2.0*U0[i][j]; vmmx10 = _mm_add_pd(vmmx10,vmmx01); // U1[i][j] = U1[i][j] + U0[i+1][j]; vmmx10 = _mm_add_pd(vmmx10,vmmx02); // U1[i][j] = U1[i][j] + U0[i-1][j]; vmmx10 = _mm_mul_pd(vb,vmmx10); // U1[i][j] = U1[i][j] * InvDxDx; vmmx10 = _mm_mul_pd(vd,vmmx10); // U1[i][j] = U1[i][j] * DtAlpha; vmmx10 = _mm_add_pd(vmmx10,vmmx00); // U1[i][j] = U1[i][j] + U0[i][j]; _mm_store_pd(&U1[i],vmmx10); // U1[i][j] = U0[i][j] + DtAlpha*( (U0[i+1][j]-2.0*U0[i][j]+U0[i-1][j])*InvDxDx } for( i = 2; i < n1-2; i+=2) { //printf("%d %d \n",i,j); //fflush(stdout); vmmx00 = _mm_load_pd(&U1[i]); vmmx01 = _mm_loadu_pd(&U1[i+1]); vmmx02 = _mm_loadu_pd(&U1[i-1]); vmmx10 = _mm_mul_pd(va,vmmx00); // U0[i][j] = -2.0*U1[i][j]; vmmx10 = _mm_add_pd(vmmx10,vmmx01); // U0[i][j] = U0[i][j] + U1[i+1][j]; vmmx10 = _mm_add_pd(vmmx10,vmmx02); // U0[i][j] = U0[i][j] + U1[i-1][j]; vmmx10 = _mm_mul_pd(vb,vmmx10); // U0[i][j] = U0[i][j] * InvDxDx; vmmx10 = _mm_mul_pd(vd,vmmx10); // U0[i][j] = U0[i][j] * DtAlpha; vmmx10 = _mm_add_pd(vmmx10,vmmx00); // U0[i][j] = U0[i][j] + U1[i][j]; _mm_store_pd(&U0[i],vmmx10); // U1[i][j] = U0[i][j] + DtAlpha*( (U0[i+1][j]-2.0*U0[i][j]+U0[i-1][j])*InvDxDx } } time1=clock(); printf("Loop 0, total time : %f \n", (double) time1-time0); f1 = fopen ("out0.dat", "wt"); for( i = 1; i < n1-1; i++) { fprintf (f1, "%d\t%f\n", i, U0[i]); } // REF for( i = 0; i < n1; i++){U0[i] = 0.0;} U0[1] = 1.0; U0[2] = 1.0; U0[3] = 1.0; U0[n1-2] = 2.0; time0=clock(); for( t = 0; t < niter; t++) { for( i = 2; i < n1-2; i++) { U1[i] = U0[i] + DtAlpha* (U0[i+1]-2.0*U0[i]+U0[i-1])*InvDxDx; } for( i = 2; i < n1-2; i++) { U0[i] = U1[i] + DtAlpha* (U1[i+1]-2.0*U1[i]+U1[i-1])*InvDxDx; } } time1=clock(); printf("Loop 0, total time : %f \n", (double) time1-time0); f1 = fopen ("outref.dat", "wt"); for( i = 1; i < n1-1; i++) { fprintf (f1, "%d\t%f\n", i, U0[i]); } } 

++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++

编辑：

++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++

考虑到你的答案，我找到了讨论这个问题的正确位置，所以我将扩大主题并解释我的目标。 如果您接受，我们将一个接一个地讨论所有循环。 这可能很长，但对于我的域中的很多人，以及OpenFoam等OpenSource求解器，它可能非常有用。 没有考虑能源消耗的影响（我们都使用大型超级计算器）。

我使用的CFD代码需要在512个Westmer Core上运行超过1个月。 我正在使用MPI（消息传递接口）在procs之间进行通信。 Physic字段可以视为网格，因此可以将1D，2D或3D数组视为模拟类型。 但3D是你能想象到的最好的。

完整代码在Fortran 95中，实际上是C简化版。 它很容易与C接口，C例程可以直接从Fortran调用，类型是相同的（int，double，long等）。 但是，Fortran不允许这样的优化：它的设计很简单。 这就是我正在调查C指令的原因。

在所有CFD代码中，我们都面临着同样的问题：3种类型的循环和MPI内存分布。 我们先来讨论一下循环：

• 空间导数（称为有限差分）：循环由所有情况的1D卷积组成（1D，2D，3D，您一次仅在一个轴上导出）（DF = F [i-1] * A + F [i ] * B + F [i + 1] * C）。 但是，当使用1D以上时，内存访问将变为以下内容：

   // x1 derivative for i 1 -> n1 for j 1 ->n2 DF_x1[i][j] = F[i-1][j]*A + F[i][j]*B + F[i+1][j]*C // x2 derivative for i 1 -> n1 for j 1 ->n2 DF_x2[i][j] = F[i][j-1]*D + F[i][j]*E + F[i][j+1]*G 

  在第一个循环中，内存访问不会继续（Fortran中的反转，内存被反转）。 这是第一个问题。 在使用3Darrays时同意。

• 泊松方程分辨率，即矩阵乘法： 根据模拟 ，循环由1D，2D或3D卷积组成。 这实际上是二阶导数（DDF = D（DF））。

   for i 1 -> n1 for j 1 ->n2 DDF[i][j] = F[i-1][j]*A + F[i][j]*B + F[i+1][j]*C + F[i][j-1]*D + F[i][j]*E + F[i][j+1]*G 

  这个循环与我第一次给你的循环相同，但它是直接计算的，不是偶数和奇数。

• 加权高斯赛德尔分辨率，即与下面相同的循环，但具有依赖性：

   // even for i 1 -> n1 for j 1 ->n2 F1[i][j] = F0[i-1][j]*A + F0[i][j]*B + F0[i+1][j]*C + F0[i][j-1]*D + F0[i][j]*E + F0[i][j+1]*G //odd for i 1 -> n1 for j 1 ->n2 F0[i][j] = F1[i-1][j]*A + F1[i][j]*B + F1[i+1][j]*C + F1[i][j-1]*D + F1[i][j]*E + F1[i][j+1]*G 

  这是您之前调查过的循环。

然后，我们面临另一个问题：内存分配。 每个核心都有自己的内存，需要与其他人共享。 让我们考虑最后一个循环，但简化：

  for t 1 -> niter // even for i 1 -> n1-2 F1[i] = F0[i-1]*A + F0[i]*B + F0[i+1]*C //odd for i 1 -> n1-2 F0[i] = F1[i-1]*A + F1[i]*B + F1[i+1]*C 

考虑到n1 = 512，但由于RAM容量较低，因此无法将其存储在本地存储器中。 内存分布在core0（1-> 255）和core1（256-512）上，这些内核不在同一台计算机上，而是在网络上。 在这种情况下，i = 256中的导数需要知道点i = 255，但该值在另一个proc上。 包含其他处理器值的内存称为GHOST内存。 所以循环是：

  ! update boundary memory : Share to ghost : core0 : F0[255] -> Network -> F0[0] : core1 (don't forget that for core1, the array restart from 0) Share to ghost : core1 : F0[1] -> Network -> F0[256] : core0 (you understand that F0[256] is the ghost for core0, and F0[0] is the ghost for core1) // even, each core do this loop. for i 1 -> n1-2 F1[i] = F0[i-1]*A + F0[i]*B + F0[i+1]*C ! update boundary memory : Share to ghost : core0 : F1[255] -> Network -> F1[0] : core1 Share to ghost : core1 : F1[1] -> Network -> F1[256] : core0 //odd, each core do this loop. for i 1 -> n1-2 F0[i] = F1[i-1]*A + F1[i]*B + F1[i+1]*C 

我们需要处理这个问题。 Mysticial ，你现在看到我要去的地方：循环隔行扫描需要考虑到这一点。 我认为这可以这样做：

  ! update boundary memory : Share to ghost : core0 : F0[255] -> Network -> F0[0] : core1 Share to ghost : core1 : F0[1] -> Network -> F0[256] : core0 for t 1 -> niter ! compute borders in advance : core0 only : F1[255] = F0[254]*A + F0[255]*B + F0[256]*C core1 only : F1[1] = F0[0]*A + F0[1]*B + F0[2]*C Launch Share to ghost asynchronous : core0 : F1[255] -> Network -> F1[0] : core1 Launch Share to ghost asynchronous : core1 : F1[1] -> Network -> F1[256] : core0 During the same time (this can be done at the same time because MPI support asynchronous communications) // even for i 2 -> n1-3 (note the reduced domain) F1[i] = F0[i-1]*A + F0[i]*B + F0[i+1]*C Check that communications are done. ! compute borders in advance : core0 only : F0[255] = F1[254]*A + F1[255]*B + F1[256]*C core1 only : F0[1] = F1[0]*A + F1[1]*B + F1[2]*C Launch Share to ghost asynchronous : core0 : F0[255] -> Network -> F0[0] : core1 Launch Share to ghost asynchronous : core1 : F0[1] -> Network -> F0[256] : core0 //odd, each core do this loop. for i 2 -> n1-3 F0[i] = F1[i-1]*A + F1[i]*B + F1[i+1]*C Check that communications are done. 

希望我某处的指数没有错。 让我们考虑当下的第一种类型的循环，这是最简单的，我们可以通过循环2和3之后，它们是相似的。 目的是这样做（类似于图像处理）：

  // x1 derivative for i 1 -> n1 for j 1 ->n2 DF_x1[i][j] = F[i-1][j]*A + F[i][j]*B + F[i+1][j]*C // x2 derivative for i 1 -> n1 for j 1 ->n2 DF_x2[i][j] = F[i][j-1]*D + F[i][j]*E + F[i][j+1]*G 

我正在努力，我将在几个小时内发布结果代码，同时考虑您的建议。

Answer 1

代码中Mysticial的明智之词：

  xmm7 = InvDxDx * dtAlpha; // precalculate xmm6 = -2*xmm7; xmm0 = *ptr++; // has values d0, d1 xmm1 = *ptr++; // has values d2, d3 while (loop--) { xmm2 = xmm0; // take a copy of d0,d1 xmm0 += xmm1; // d0+d2, d1+d3 xmm2 = shufps (xmm2,xmm1, 0x47); // the middle elements d1,d2 ?? xmm0 *= xmm7; // sum of outer elements * factor xmm2 *= xmm6; // -2 * center element * factor // here's still a nasty dependency xmm2 += xmm0; xmm0 = xmm1; // shift registers *out_ptr ++= xmm2; // flush xmm1 = *ptr++; // read in new values } 

这可以通过将循环分成2个段来改进，每个段分开502个条目或者来自不同的任务并交织结果，这打破了依赖链。

此外，移位寄存器接近xmm0 < - xmm1，xmm1 < - new可以通过展开循环两次并在每种情况下改变xmm0和xmm1的含义来避免。

这指出了另一个大问题：当使用内在函数时，寄存器一直溢出到内存中。