查找具有最大值的元素的INDEX使用AVX512指令的绝对值实际上可以在FP位图上使用整数运算修正的优化版本：随机策略：当新的最大值非常罕见时

问题描述

我是使用AVX512指令进行编码的新手。我的机器是Intel KNL7250。我试图使用AVX512指令来查找具有最大绝对值的元素的INDEX，该元素的双精度和数组％8 = 0的大小加倍。但是，每次打印输出索引= 0。我不知道哪里有问题，请帮助我。另外，如何对__m512i类型使用 printf ？

谢谢。

代码：

void main()
{
    int i;
    int N=160;
    double vec[N];

    for(i=0;i<N;i++)
    {
        vec[i]=(double)(-i) ;
        if(i==10)
        {
            vec[i] = -1127;
        }
    }

    int max = avxmax_(N,vec);
    printf("maxindex=%d\n",max);
}

int avxmax_(int   N,double *X )
{
    // return the index of element having maximum absolute value.
    int maxindex,ix,i,M;
    register __m512i increment,indices,maxindices,maxvalues,absmax,max_values_tmp,abs_max_tmp,tmp;
    register __mmask8 gt;
    double values_X[8];
    double indices_X[8];
    double maxvalue;
    maxindex = 1;
    if( N == 1) return(maxindex);


    M = N % 8;
    if( M == 0)
    {
        increment =  _mm512_set1_epi64(8); // [8,8,8]
        indices = _mm512_setr_epi64(0,1,2,3,4,5,6,7);
        maxindices = indices;
        maxvalues =  _mm512_loadu_si512(&X[0]);
        absmax = _mm512_abs_epi64(maxvalues); 

        for( i = 8; i < N; i += 8)
        {
            // advance scalar indices: indices[0] + 8,indices[1] + 8,...,indices[7] + 8
            indices = _mm512_add_epi64(indices,increment);

            // compare
            max_values_tmp  = _mm512_loadu_si512(&X[i]);
            abs_max_tmp = _mm512_abs_epi64(max_values_tmp);
            gt           = _mm512_cmpgt_epi64_mask(abs_max_tmp,absmax);

            // update
            maxindices = _mm512_mask_blend_epi64(gt,indices);
            absmax     = _mm512_max_epi64(absmax,abs_max_tmp);
        }

        // scalar part
        _mm512_storeu_si512((__m512i*)values_X,absmax);
        _mm512_storeu_si512((__m512i*)indices_X,maxindices);
        maxindex = indices_X[0];
        maxvalue = values_X[0];
        for(i = 1; i < 8; i++)
        {
            if(values_X[i] > maxvalue)
            {
                maxvalue = values_X[i];
                maxindex = indices_X[i];
            }
            
        }
        return(maxindex);
    }
}

解决方法

您的函数返回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转换为整数。

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

实际上可以在FP位图上使用整数运算

但仅当您使用正确的选项时！ IEEE754指数偏差表示可以将编码/位模式作为符号/幅度整数进行比较。因此，您可以进行abs / min / max进行比较，但不能进行整数加/减（除非您要实现nextafter）。

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

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

在KNL（骑士降落）上，FP vmaxpd和vcmppd的性能与其整数等效值相同：2个周期的延迟，0.5c的吞吐量。（{{3 }}。因此，您的方式在KNL上的优势为零，但是对于Skylake-X和IceLake等主流英特尔来说，这是一个绝妙的窍门。

修正的优化版本：

对返回类型和循环计数器/索引使用size_t处理潜在的巨大数组，而不是int和64位向量元素的随机混合。（{uint64_t用于收集水平最大值的临时数组：即使在具有32位指针/ size_t的构建中，它也始终为64位。）
修正：在N == 1而不是1上返回0，唯一元素的索引为0。
错误修正：在-1上返回N%8 != 0，而不是掉到非void函数的末尾。（如果调用者在C中使用结果，或者在C ++中执行一旦结束，则为未定义的行为。）
错误修正：FP值的绝对值=清除符号位，而不是位模式上的2的补码绝对值
某种错误修复：使用无符号整数比较和最大值，因此它将与_mm512_abs_epi64用于2的补码整数（产生无符号结果；请记住-LONG_MIN溢出到{{1 }}（如果您继续将其视为已签名）。
样式改进：LONG_MIN而不是将大部分内容放在if块中。
样式改进：首次使用时声明vars ，并删除了一些纯噪声未使用的变量。自20年前标准化C99以来，这对于C来说是惯用的。
样式改进：对仅保存加载结果的tmp矢量变量使用更简单的名称。有时，您只需要一个tmp var，因为内部名称太长了，您不想将if (N%8 != 0) return -1;用作另一个内部参数的arg。像_mm...load...这样的名称作用于内部循环是一个明显的标志，它只是一个占位符，以后不再使用。（此样式在初始化时声明时效果最好，因此很容易看出它不能在外部范围中使用。）
优化：使用SIMD循环后减少8-> 4个元素：提取上半部分，并与现有的下半部分合并。（与用于简单的水平缩减（例如sum或max）的方法相同）。当我们需要仅AVX512拥有但KNL没有AVX512VL的指令时很不方便，因此我们必须使用512位版本，而忽略高垃圾率。但是KNL确实有AVX1 / AVX2，因此我们仍然可以存储256位向量并做一些事情。

使用合并遮罩v提取来直接混合同一向量的下半部分是一个很酷的技巧，如果您使用512位的遮罩混合，编译器将找不到这些技巧。：P
一种错误修复：在C中，_mm512_mask_extracti64x4_epi64在托管实现（在OS中运行）中的返回类型为main。

int

https://agner.org/optimize/ ：来自GCC10.2 -O3 -march = knl的主循环为8条指令。因此，即使（最好的情况）KNL可以以2 / clock的频率解码并运行它，每个向量仍然要花费4个周期。您可以在Godbolt上运行该程序；它在Skylake-X服务器上运行，因此可以运行AVX512代码。您可以看到它打印了#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,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; }。

循环的另一个选项是屏蔽的.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,zmm2 ja .L4 vmovapd zmm1,zmm3 # silly missed optimization related to the case where the loop runs 0 times. .L3: vextracti64x4 ymm2,0x1 # high half of absmax vpcmpuq k1,6 # compare high and low vpmaxuq zmm0,zmm2 # vunpckhpd xmm2,xmm0,xmm0 # setting up for unrolled scalar loop vextracti64x4 ymm1{k1},zmm3,0x1 # masked extract of indices，在循环之后添加每个元素偏移量vpbroadcastq zmm3{k1},rax。这实际上可以将[0..7]保存在循环中，并且如果GCC仍将使用索引寻址模式，则我们在寄存器中拥有正确的vpaddq。（在Skylake-X上不好；击败了内存源i的微融合。）

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

随机策略：当新的最大值非常罕见时

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

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

尤其是使用KNL的4向SMT来隐藏分支未命中成本时，这可能是大型阵列在总体吞吐量方面的优势；平均每个元素的指令更少。

但对于 do 呈上升趋势的输入，可能会严重恶化，因此平均发现新的最大值是O（n）倍，而不是sqrt（n）或log（n）或统一频率分配会给我们。

PS：打印矢量，存储到数组中并重新加载元素。 https://uops.info/

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

avx512 c instructions max