查找具有最大值的元素的INDEX。使用AVX512指令的绝对值

函数返回0是因为您将int64索引视为的位模式double，并将该（小）数字转换为整数。 double indices_X[8];是错误；应该是uint64_t。还有其他错误，请参见下文。

如果在使用变量时声明了C99样式而不是过时的C89样式，则更容易发现此错误。

您_mm512_storeu_si512将int64_t索引向量转换为double indices_X[8]，将其类型化为两倍，然后以纯C形式进行int maxindex = indices_X[0];。这是隐式类型转换，可将该子范式double转换为整数。

（我注意到，在将代码转换为初始化程序旁边的C99样式变量声明时vcvttsd2si，在asm https://godbolt.org/z/zsfc36中进行了神秘的FP-> int转换。这是一个线索：不应有FP-> int我注意到大约在同一时间，我将double indices_X[8];声明下移到使用它的块中，并注意到它具有类型double。

它是实际上可以使用整数运算上FP位模式

但前提是您使用正确的！IEEE754指数偏差表示可以将编码/位模式作为符号/幅度整数进行比较。因此，您可以进行abs / min / max进行比较，但是当然不能进行整数add / sub（除非您正在实现nextafter）。

_mm512_abs_epi64是2的补码绝对值，不是符号幅度。相反，您必须仅屏蔽符号位。然后，您都准备将结果视为无符号整数或有符号2的补码。（之所以起作用，是因为高位清除了。）

使用整数最大值具有有趣的性质，即NaNs的比较将高于任何其他值，Inf低于该值，然后是有限值。因此，我们基本上免费获得NaN传播的寻星器。

在KNL（Knight's Landing）上，FPvmaxpd和FPvcmppd具有与整数等效值相同的性能：2个周期的延迟，0.5c的吞吐量。（https://agner.org/optimize/）。因此，您的方式在KNL上的优势为零，但这对于Skylake-X和IceLake等主流英特尔来说是一个绝妙的窍门。

错误修正的优化版本：

• 使用size_t用于返回类型和循环计数器/索引来处理潜在的巨大阵列，而不是随机的混合int和64位的向量元素。（uint64_t对于收集horizo​​ntal-max的临时数组：即使在具有32位指针/ size_t的构建中，它也始终为64位。）

• 错误修正：返回0N == 1，而不是1：唯一元素的索引为0。

• 返程：bug修复-1上N%8 != 0的，而不是脱落的非void函数结束。（如果调用者在C中使用结果，或者在C ++中执行一旦结束，则为未定义的行为）。

• 错误修正：FP值的绝对值=清除符号位，而不是位模式上的2的补码绝对值

• 错误修复：使用无符号整数compare和max，因此它适用于2的补码整数_mm512_abs_epi64（产生无符号结果；请记住，如果继续将其视为有符号-LONG_MIN，LONG_MIN则会溢出到）。

• 改善风格：if (N%8 != 0) return -1;不要将大部分身体放在if块中。

• 样式改进：首次使用时声明var，并删除了一些纯噪声未使用的变量。自20年前标准化C99以来，这对于C来说是惯用的。

• 样式改进：为tmp向量var使用更简单的名称，它们仅保存加载结果。有时您只需要一个tmp var，因为内在名称很长，以至于您不想_mm...load...为另一个内在输入arg。像v作用于内部循环的名称是一个明显的标志，它只是一个占位符，以后不再使用。（此样式在初始化时声明时效果最好，因此很容易看出它不能在外部范围中使用。）

• 优化：使用SIMD在循环后减少8-> 4个元素：提取上半部分，并与现有的下半部分合并。（与您进行简单的水平缩减（例如sum或max）相同）。当我们需要仅AVX512具有但KNL没有AVX512VL的指令时，这很不方便，因此我们必须使用512位版本，而忽略高垃圾率。但是KNL确实有AVX1 / AVX2，因此我们仍然可以存储256位向量并做一些事情。

  使用合并屏蔽_mm512_mask_extracti64x4_epi64提取来直接混合同一矢量的下半部分是一个很酷的技巧，如果您使用512位掩码混合，编译器将找不到这些技巧。：P

• 某种错误修复：在C中，托管实现main的返回类型为int（在OS下运行）。

#include <immintrin.h>
#include <stdio.h> 
#include <stdint.h>
#include <stdlib.h>

// bugfix: indices can be larger than an int
size_t avxmax_(size_t N,  double *X )
{
    // return the index of element having maximum absolute value.
    if( N == 1)
        return 0;    // bugfix: 0 is the only valid element in this case, not 1
    if( N % 8 != 0)      // [[unlikely]] // C++20
        return -1;   // bugfix: don't fall off the end of the function in this case

    const __m512i fp_absmask = _mm512_set1_epi64(0x7FFFFFFFFFFFFFFF);
    __m512i indices = _mm512_setr_epi64(0, 1, 2, 3, 4, 5, 6, 7);
    __m512i maxindices = indices;
    __m512i v =  _mm512_loadu_si512(&X[0]);
    __m512i absmax = _mm512_and_si512(v, fp_absmask);
    for(size_t i = 8; i < N; i += 8) // [[likely]]  // C++20
    {
        // advance indices by 8 each.
        indices = _mm512_add_epi64(indices, _mm512_set1_epi64(8));
        // compare
                v    = _mm512_loadu_si512(&X[i]);
        __m512i vabs = _mm512_and_si512(v, fp_absmask);
             // vabs = _mm512_abs_epi64(max_values_tmp);  // for actual integers, not FP bit patterns
        __mmask8 gt  = _mm512_cmpgt_epu64_mask(vabs, absmax);
        // update
        maxindices = _mm512_mask_blend_epi64(gt, maxindices, indices);
        absmax     = _mm512_max_epu64(absmax, vabs);
    }

    // reduce 512->256; KNL doesn't have AVX512VL so some ops require 512-bit vectors
    __m256i absmax_hi = _mm512_extracti64x4_epi64(absmax, 1);
    __m512i absmax_hi512 = _mm512_castsi256_si512(absmax_hi); // free
    __mmask8 gt = _mm512_cmpgt_epu64_mask(absmax_hi512, absmax);
    __m256i abs256 = _mm512_castsi512_si256(_mm512_max_epu64(absmax_hi512, absmax));  // reduced to low 4 elements

    // extract with merge-masking = blend
    __m256i maxindices256 = _mm512_mask_extracti64x4_epi64(
                           _mm512_castsi512_si256(maxindices), gt, maxindices, 1);

    // scalar part
    double values_X[4];
    uint64_t indices_X[4];
    _mm256_storeu_si256((__m256i*)values_X, abs256);
    _mm256_storeu_si256((__m256i*)indices_X, maxindices256);

    size_t maxindex = indices_X[0];
    double  maxvalue = values_X[0];
    for(int i = 1; i < 4; i++)
    {
        if(values_X[i] > maxvalue)
        {
            maxvalue = values_X[i];
            maxindex = indices_X[i];
        }
        
    }
    return maxindex;
}

在Godbolt上：GCC10.2 -O3 -march = knl的主循环是8条指令。因此，即使（最好的情况）KNL可以以2 / clock的频率解码并运行它，每个向量仍然要花费4个周期。您可以在Godbolt上运行该程序；它在Skylake-X服务器上运行，因此可以运行AVX512代码。您可以看到它打印出来10。

.L4:
        vpandd  zmm2, zmm5, ZMMWORD PTR [rsi+rax*8]   # load, folded into AND
        add     rax, 8
        vpcmpuq k1, zmm2, zmm0, 6
        vpaddq  zmm1, zmm1, zmm4                    # increment current indices
        cmp     rdi, rax
        vmovdqa64       zmm3{k1}, zmm1              # blend maxidx using merge-masking
        vpmaxuq zmm0, zmm0, zmm2
        ja      .L4

        vmovapd zmm1, zmm3                        # silly missed optimization related to the case where the loop runs 0 times.
.L3:
        vextracti64x4   ymm2, zmm0, 0x1           # high half of absmax
        vpcmpuq k1, zmm2, zmm0, 6                 # compare high and low
        vpmaxuq zmm0, zmm0, zmm2
      #  vunpckhpd       xmm2, xmm0, xmm0  # setting up for unrolled scalar loop
        vextracti64x4   ymm1{k1}, zmm3, 0x1       # masked extract of indices

循环的另一个选项是masked vpbroadcastq zmm3{k1}, rax，[0..7]在循环之后添加每个元素的偏移量。这实际上将节省vpaddq循环，并且i如果GCC仍将使用索引寻址模式，则我们拥有寄存器的权利。（在Skylake-X上不好；击败了内存源的微融合vpandd。）

Agner Fog没有列出GP-> vector广播的性能，但希望它至少在KNL上只有单播。（并且https://uops.info/没有KNL或KNM结果）。

分支策略：很少有新的最大值

如果您期望找到一个新的最大值非常罕见（例如，数组很大且分布均匀，或者至少没有向上趋势），则广播当前最大值并在找到任何更大的向量元素时分支可能会更快。

找到一个新的max意味着跳出循环（可能会预测错误，因此很慢）并广播该元素（可能使用atzcnt查找元素索引，然后广播加载并更新索引）。

特别是使用KNL的4路SMT来隐藏分支未命中成本时，对于大型阵列而言，这可能是整体吞吐量的胜利；平均每个元素更少的指令。

但对于确实呈上升趋势的输入，可能会变得更糟，因此，平均发现新的最大值是O（n）倍，而不是sqrt（n）或log（n）或均匀分布给我们的频率。

PS：打印矢量，存储到数组中并重新加载元素。打印__m128i变量

或使用调试器向您展示其元素。