最近，我通过有限差异方法在计算电动力学上的数值求解器上工作。
求解器实现非常简单，但是很难达到现代处理器的理论吞吐量，因为加载数据上只有1个数学操作，例如：
#pragma ivdep
for(int ii=0;ii<Large_Number;ii++)
{ Z[ii] = C1*Z[ii] + C2*D[ii];}
大_number约为1,000,000，但不超过10,000,000
我试图手动展开循环并编写AVX代码，但无法使其更快：
int Vec_Size = 8;
int Unroll_Num = 6;
int remainder = Large_Number%(Vec_Size*Unroll_Num);
int iter = Large_Number/(Vec_Size*Unroll_Num);
int addr_incr = Vec_Size*Unroll_Num;
__m256 AVX_Div1, AVX_Div2, AVX_Div3, AVX_Div4, AVX_Div5, AVX_Div6;
__m256 AVX_Z1, AVX_Z2, AVX_Z3, AVX_Z4, AVX_Z5, AVX_Z6;
__m256 AVX_Zb = _mm256_set1_ps(Zb);
__m256 AVX_Za = _mm256_set1_ps(Za);
for(int it=0;it<iter;it++)
{
int addr = addr + addr_incr;
AVX_Div1 = _mm256_loadu_ps(&Div1[addr]);
AVX_Z1 = _mm256_loadu_ps(&Z[addr]);
AVX_Z1 = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div1),_mm256_mul_ps(AVX_Za,AVX_Z1));
_mm256_storeu_ps(&Z[addr],AVX_Z1);
AVX_Div2 = _mm256_loadu_ps(&Div1[addr+8]);
AVX_Z2 = _mm256_loadu_ps(&Z[addr+8]);
AVX_Z2 = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div2),_mm256_mul_ps(AVX_Za,AVX_Z2));
_mm256_storeu_ps(&Z[addr+8],AVX_Z2);
AVX_Div3 = _mm256_loadu_ps(&Div1[addr+16]);
AVX_Z3 = _mm256_loadu_ps(&Z[addr+16]);
AVX_Z3 = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div3),_mm256_mul_ps(AVX_Za,AVX_Z3));
_mm256_storeu_ps(&Z[addr+16],AVX_Z3);
AVX_Div4 = _mm256_loadu_ps(&Div1[addr+24]);
AVX_Z4 = _mm256_loadu_ps(&Z[addr+24]);
AVX_Z4 = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div4),_mm256_mul_ps(AVX_Za,AVX_Z4));
_mm256_storeu_ps(&Z[addr+24],AVX_Z4);
AVX_Div5 = _mm256_loadu_ps(&Div1[addr+32]);
AVX_Z5 = _mm256_loadu_ps(&Z[addr+32]);
AVX_Z5 = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div5),_mm256_mul_ps(AVX_Za,AVX_Z5));
_mm256_storeu_ps(&Z[addr+32],AVX_Z5);
AVX_Div6 = _mm256_loadu_ps(&Div1[addr+40]);
AVX_Z6 = _mm256_loadu_ps(&Z[addr+40]);
AVX_Z6 = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div6),_mm256_mul_ps(AVX_Za,AVX_Z6));
_mm256_storeu_ps(&Z[addr+40],AVX_Z6);
}
上面的AVX循环实际上比编译器生成的代码慢一点。
编译器生成的代码可以达到约8G的拖鞋/s，约有3GHz Ivybridge处理器的单线理论吞吐量的25％。我想知道是否有可能达到这样的简单循环的吞吐量。
谢谢！
改善像您这样的代码的性能是"探索良好的"且仍然流行的区域。查看DOT-Product（Z Boson已经提供的完美链接）或（D）Axpy优化讨论（https://scicomp.stackexchange.com/questions/questions/1932/are-daxpy-daxpy-dcopy-dcopy-dscal-overkills）
通常，要探索和考虑应用的关键主题是：
- 由于 fma ，AVX2优于avx的优势，以及更好的负载/存储端口U-Architection u-Architection a haswell
- 预取。某些平台的"流式商店"，"非时空存储"。
- 螺纹并行性达到最大带宽
- 展开（您已经完成了；英特尔编译器也能够使用#pragma展开（x））。对于"流"代码而言，这不是很大的区别。
- 最终确定要优化代码 的一组硬件平台是什么
最后一颗子弹特别重要，因为对于"流"和整体内存的代码，重要的是要了解有关目标内存综合系统的更多信息；例如，使用现有的，尤其是未来的高端HPC服务器（第二代Xeon Phi代号为骑士登陆为例），您可能在带宽和计算之间具有不同的"屋顶线模型"平衡，甚至与优化的技术之间的平衡相比平均台式机。
您确定8 Gflops/s约为3 GHz Ivybridge处理器的最大吞吐量的25％吗？让我们进行计算。
每8个元素需要两个单精度的AVX乘法和一个AVX添加。Ivybridge处理器每个周期只能执行一个8宽的AVX添加和一个8宽的AVX乘法。同样，由于添加取决于两个乘法，因此需要3个周期才能处理8个元素。由于添加可以与下一个乘法重叠，因此我们可以将其减少到每个8个元素的2个周期。对于十亿个元素，需要2*10^9/8 = 10^9/4周期。考虑到3 GHz时钟，我们得到10^9/4 * 10^-9/3 = 1/12 = 0.08秒。因此，最大的理论吞吐量为12 GLOPS/S，编译器生成的代码达到66％，这很好。
另外一件事，通过将循环展开8次，可以有效地将其进行矢量化。我怀疑，如果您将此特定循环展开的速度超过16次，您会获得任何显着的速度。
我认为真正的瓶颈是每2个乘法有2个负载和1个存储说明，而1个添加。也许计算是内存带宽有限的。每个元素都需要传输12个数据的数据，如果每秒处理2G元素（6G FLOPS），即24GB/s数据传输，则达到IVY桥的理论带宽。我想知道这个论点是否存在，确实没有解决这个问题的方法。
我要回答自己的问题的原因是希望有人可以纠正我，然后才能轻松放弃优化。这个简单的循环对于许多科学求解器来说非常重要，它是有限元和有限差异方法的骨干。如果一个人甚至因为计算是内存bandwith Limited而无法喂养一个处理器，为什么要使用多项核心？高频段gpu或Xeon Phi应该是更好的解决方案。