医学图像重建中的局部化和减少缓存污染

Question

我正在为我的大学做一项与医学使用的图像重建算法相关的研究。

我被困在长达3周的时间里，我需要改进以下代码的性能：

for (lor=lor0[mypid]; lor <= lor1[mypid]; lor++)
{
  LOR_X = P.symmLOR[lor].x;
  LOR_Y = P.symmLOR[lor].y;
  LOR_XY = P.symmLOR[lor].xy;
  lor_z = P.symmLOR[lor].z;
  LOR_Z_X = P.symmLOR[lor_z].x;
  LOR_Z_Y = P.symmLOR[lor_z].y;
  LOR_Z_XY = P.symmLOR[lor_z].xy;  

  s0 = P.a2r[lor];
  s1 = P.a2r[lor+1];

  for (s=s0; s < s1; s++)
  {
    pixel     = P.a2b[s];
    v         = P.a2p[s]; 

    b[lor]    += v * x[pixel];

    p          = P.symm_Xpixel[pixel];
    b[LOR_X]  += v * x[p];

    p          = P.symm_Ypixel[pixel];
    b[LOR_Y]  += v * x[p];

    p          = P.symm_XYpixel[pixel];
    b[LOR_XY] += v * x[p];


    // do Z symmetry.
    pixel_z    = P.symm_Zpixel[pixel];
    b[lor_z]  += v * x[pixel_z];


    p          = P.symm_Xpixel[pixel_z];
    b[LOR_Z_X]  += v * x[p];


    p          = P.symm_Ypixel[pixel_z];
    b[LOR_Z_Y]  += v * x[p];

    p          = P.symm_XYpixel[pixel_z];
    b[LOR_Z_XY] += v * x[p];

   }

}

对于任何想知道的人，代码实现了MLEM前向函数，并且所有变量都是浮点。

经过几次测试后，我注意到代码的这一部分出现了很大的延迟。(你知道，90-10规则)。

后来，我使用Papi (http://cl.cs.utk.edu/papi/)来测量L1D缓存丢失。正如我所想的，Papi确认性能下降是由于更多的错误，特别是对b矢量的随机访问(体积很大)。

在互联网上阅读信息，我只知道到目前为止提高性能的两个选项:提高数据局部性或减少数据污染。

为了进行第一个改进，我将尝试更改代码，使其具有缓存感知能力，就像Ulrich Drepper提出的关于每个程序员应该了解内存的内容一样(www.akkadia.org/drepper/cpumymy.pdf) A.1矩阵乘法。

我相信阻断SpMV (稀疏矩阵-向量乘法)将提高性能。

另一方面，每当程序试图访问b向量时，我们就会受到所谓的缓存污染。

有没有办法在不使用缓存的情况下，用SIMD指令从b向量加载一个值？

此外，还可以使用像void ( float *p，__m128 a)这样的函数在向量b上存储一个浮点数，而不会污染缓存吗？

我不能使用_mm_stream_ps，因为总是存储4个浮点数，但是对b向量的访问显然是随机的。

我希望能在我的困境中说清楚。

更多信息:v是CRS格式稀疏矩阵存储的列值。我意识到，如果我尝试将CRS格式更改为其他格式，还可以进行其他优化，但是，正如我之前说过的，我已经做了几个月的测试，并且我知道性能下降与向量b上的随机访问有关。当我不存储在向量b中时，我可以将400.000.000个L1D丢失到100~漏出。

谢谢。

Answer 1

我想说，一开始试着帮助你的编译器一点。

在循环之前，size_t s0 = P.a2r[lor];

• This

• 将外部循环的边界声明为const。

• 将所有可能的变量(所有LOR_..)声明为局部变量，类似于：float LOR_X = P.symmLOR[lor].x;或循环变量，特别是对于循环变量，如果您碰巧有符合C99的现代编译器：for (size_t s=s0; s < s1; s++)

• Separate load和存储您的b矢量。您在那里访问的项目的位置不依赖于s。因此，创建一个局部变量来保存内环之前处理的所有不同情况的实际值，更新内环中的这些局部变量，并在内部loop.

• Perhaps将内部循环分隔成几个后存储结果。索引计算相对便宜，然后您的系统可能会更好地识别编译器生成的汇编程序上对某些vectors.

• Look的流访问，并识别内部循环的代码。尝试一下，移动尽可能多的“远”负载，并将存储从循环中移出。

编辑:在重读引力子的答案和你的评论之后，这里最重要的事情是尽可能地声明变量，并检查汇编程序编译器是否成功地加载了缺少缓存的负载，并将其存储在内部循环之外。

Answer 2

要减少向量b的随机访问，一个简单的优化就是不要在内部for循环中写入向量b。

相反，将向量B所需的所有值加载到临时变量中，在更新这些临时变量时执行整个内部for循环，然后将临时变量写回向量B。

在最坏的情况下，临时变量将位于相同的缓存行上，取决于编译器和环境，您可能还提示编译器对这些变量使用寄存器。

Answer 3

我甚至不会假装我知道代码在做什么:)，但是一些额外内存访问的一个可能原因是别名:如果编译器不能确定b、x和各种P.symm数组不重叠，那么写到b将影响从x和P.symm的读取如何被调度。如果编译器特别悲观，它甚至可能会强制从内存中重取P字段。所有这些都会导致你所看到的缓存丢失。改善这一状况的两个简单方法是：

1. 在b上使用__restrict。这保证了arrays.
2. Reorder不会与其他数组重叠，因此对它的写入不会影响对其他b的读取(或写入)，从而使来自P.symm的所有读取首先是从P.symm的读取，然后是从x的读取，最后是对b的所有写入。这将打破读取和编译器调度中的一些依赖关系--对P.symm的并行读取，然后并行地从x读取，并希望能够明智地完成对b的写入。

另一件风格化的事情(对#2有帮助)是不要重用变量，所以p就这么多了。您没有理由不能拥有例如p_x、p_y、p_xy等，这将使重新排序代码变得更容易。

一旦一切就绪，您就可以在已知的缓存错误之前开始喷洒预取指令(即gcc上的__builtin_prefetch )。

希望这能有所帮助。