c - 一些从c到sse2的mandelbrot绘制例程

Question

我想将这样简单的例程重写为SSE2代码，（最好
在nasm中），但我不太确定该怎么做，两件事
不清楚（如何表达计算（内部循环以及来自
外循环）以及如何调用c代码函数“ SetPixelInDibInt（i，j，palette [n]）;”
从静态链接的汇编代码下

    void DrawMandelbrotD(double ox, double oy, double lx, int N_ITER)
    {
     double ly = lx * double(CLIENT_Y)/double(CLIENT_X);
     double dx = lx / CLIENT_X;
     double dy = ly / CLIENT_Y;
     double ax = ox - lx * 0.5 + dx * 0.5;
     double ay = oy - ly * 0.5 + dy * 0.5;
    static  double re, im, re_n, im_n, c_re, c_im, rere, imim, int n;

    for(int j=0; j<CLIENT_Y; j+=1)
    {
     for(int i=0; i<CLIENT_X; i+=1)
     {
      c_re = ax + i * dx;
      c_im = ay + j * dy;
      re = c_re;
      im = c_im;
      rere=re*re;
      imim=im*im;
      n=1;

      for(int k=0;k<N_ITER;k++)
      {
        im =  (re+re)*im    + c_im;
        re =   rere - imim  + c_re;
        rere=re*re;
        imim=im*im;
        if ( (rere + imim) > 4.0 ) break;
        n++;
       }
        SetPixelInDibInt(i ,j, palette[n]);
      }
     }
    }

有人可以帮忙，我不想看其他代码
的实现，但是上面的实现只是nasm-sse的翻译
-这对我来说是最有帮助的-有人可以帮忙吗？

Answer 1

英特尔作为AVX示例具有完整的实现。见下文。

使Mandelbrot棘手的是，集合中每个点（即像素）的提前提取条件都不同。您可以保持一对或四像素重复，直到两个像素的大小都超过2.0（或达到最大迭代次数）。否则，将需要跟踪哪个像素点位于哪个矢量元素中。

无论如何，一次简单地实现以2的向量（或带有AVX的4）加倍的向量实现的吞吐量将受到依赖链的等待时间的限制。您需要并行执行多个依赖项链，以保持同时向Haswell的FMA单元供电。因此，您需要复制变量，并在内部循环内对外部循环的两次迭代进行交错操作。

跟踪正在计算哪些像素会有些棘手。我认为使用一组寄存器用于一行像素，而另一组寄存器用于另一行可能会减少开销。 （因此，您始终可以仅向右移动4个像素，而不用检查其他dep链是否已经在处理该矢量。）

我怀疑仅每4次迭代检查一次循环退出条件可能会获胜。根据打包矢量比较来使代码分支，比标量情况要贵一些。所需的额外FP添加也很昂贵。 （Haswell每个周期可以执行两个FMA，（延迟= 5）。单个FP添加单元与FMA单元之一是同一端口。两个FP mul单元位于可以运行FMA的相同端口上。）

可以使用压缩比较检查循环条件以生成零和一的掩码，以及该寄存器与自身的(V)PTEST来查看其是否全部为零。 （编辑：movmskps然后test+jcc的uops较少，但可能会有更高的延迟。）然后，显然je或jne适当，这取决于您是否进行了FP比较，当退出时是零，还是零当你不应该。 NAN应该不可能，但是没有理由不选择比较操作，使得NAN将导致退出条件为true。

const __mm256d const_four = _mm256_set1_pd(4.0);  // outside the loop

__m256i cmp_result = _mm256_cmp_pd(mag_squared, const_four, _CMP_LE_OQ);  // vcmppd.  result is non-zero if at least one element < 4.0
if (_mm256_testz_si256(cmp_result, cmp_result))
    break;

可能有一些方法可以直接在压缩双精度上使用 PTEST，并且可以使用一些位修改AND掩码，以选择FP值> 4.0时要设置的位。像是指数中的一些位？也许值得考虑。我找到了一个 forum post，但是没有尝试。

嗯，糟糕，当循环条件失败时，这不会记录每个矢量元素的单独情况，目的是为Mandelbrot集外的点着色。也许测试是否有任何符合条件的元素（而不是全部），记录结果，然后将该元素（以及该元素的 c）设置为0.0，这样就不会再次触发退出条件。也许将像素调度到矢量元素中毕竟是要走的路。这段代码在超线程CPU上可能会做得很好，因为每个元素都会触发很多分支错误预测，从而分别触发提前条件。

这可能会浪费大量的吞吐量，并且每个周期可以执行4微指令，但是FP mul / add / FMA中只有2微指令可以使用，因此有大量的整数代码可以将点调度到向量元素中。 （在没有FMA的Sandybridge / Ivybrideg上，FP吞吐量较低。但是只有3个端口可以处理整数op，其中2个是FP mul和FP add单元的端口。）

由于您不必读取任何源数据，因此每个dep链只有1个内存访问流，而这是一个写流。 （它的带宽很低，因为大多数点在准备编写单个像素值之前需要进行很多迭代。）因此，硬件预取流的数量并不是并行运行的Dep链数量的限制因素。 。缓存未命中延迟应由写缓冲区隐藏。

如果有人仍然对此感兴趣，我可以编写一些代码（只需发表评论）。我停在高级设计阶段，因为这是一个老问题。

==============

我还发现，英特尔已经使用Mandelbrot集作为其 AVX tutorials之一的示例。他们将mask-off-vector-elements方法用于循环条件。 （使用 vcmpps直接生成的掩码与AND）。他们的结果表明，AVX（单精度）的速度比标量浮点数提高了7倍，因此，相邻像素以不同的迭代次数达到提前条件的情况显然并不常见。 （至少对于他们测试的缩放/平移。）

他们只是让FP结果继续为那些未能通过提前条件的元素进行累积。他们只是停止增加该元素的计数器。希望大多数系统默认将控制字设置为零异常，如果异常仍然需要额外的周期。

但是，他们的代码在一种方式上很愚蠢：他们使用浮点向量跟踪每个向量元素的迭代计数，然后在使用前将其转换为int。为此，使用打包整数会更快，并且不会占用FP执行单元。哦，我知道他们为什么这样做：AVX（不带AVX2）不支持256位整数矢量操作。他们本来可以使用打包的16位int循环计数器，但这可能会溢出。 （而且他们必须将掩码从256b压缩到128b）。

他们还使用 movmskps测试所有大于4.0的元素，然后测试而不是使用 ptest。我猜想 test / jcc可以宏融合，并且可以在与FP vector ops不同的执行单元上运行，所以它甚至可能不会更慢。哦，当然，AVX（没有AVX2）没有256位 PTEST。同样， PTEST为2微秒，因此，实际上 movmskps + test / jcc比 ptest + jcc少微秒。 （ PTEST是SnB上的1个融合域uop，但对于执行端口仍是2个未融合的uops。在IvB / HSW上，即使在融合域中也是2个uops。）因此，看起来 movmskps是最佳方法，除非您可以利用 PTEST中按位与的优势，或者需要测试的不仅仅是每个元素的高位。如果分支是不可预测的，则 ptest可能会降低等待时间，因此可以通过更快地捕获错误的预测来值得。