找到具有最大值的元素的索引。使用 AVX512 指令的绝对值
Find the INDEX of element having max. absolute value using AVX512 instructions
我是使用 AVX512 指令编码的新手。我的机器是英特尔 KNL 7250。我正在尝试使用 AVX512 指令来查找具有最大绝对值的元素的索引，该元素的双精度和数组大小 % 8 = 0。但它每次都会打印输出索引 = 0。不知道哪里出了问题，求大神指点。
另外，如何将 printf 用于 __m512i 类型？
谢谢。
代码：
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,8,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, maxindices, 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);
}
}
您的函数 returns 0 因为您将 int64 索引视为 double 的 bit-pattern 并将其转换（ tiny) 数字到整数。 double indices_X[8]; 是错误；应该是uint64_t。 还有其他bug，见下文。
如果您在使用变量时声明变量，C99 风格，而不是过时的 C89 风格，则更容易发现此错误。
你 _mm512_storeu_si512 int64_t 的向量索引到 double indices_X[8]，type-punning 它加倍，然后在普通 C 中做 int maxindex = indices_X[0];。这是隐式的 type-conversion，将次正规 double 转换为整数。
(我注意到在 asm https://godbolt.org/z/zsfc36 中有一个神秘的 vcvttsd2si FP->int 转换，同时将代码转换为初始化器旁边的 C99 样式变量声明。这是一个线索：应该没有这个函数中的 FP->int 转换。我注意到大约在同一时间我将 double indices_X[8]; 声明向下移动到使用它的块中，并注意到它的类型为 double。）
实际上可以在 FP bit-patterns
上使用整数运算
但前提是你用对了！ IEEE754 指数偏差意味着编码 / bit-pattern 可以作为 sign/magnitude 整数进行比较。 所以你可以做 abs / min / max 并比较它，但是当然不是整数加法/减法（除非你正在实施nextafter）。
_mm512_abs_epi64是2的补abs，不是sign-magnitude。相反，您必须屏蔽掉符号位。然后您就可以将结果视为无符号整数或有符号 2's-complement。 （两者都有效，因为高位已清除。）
使用整数 max 有一个有趣的 属性，即 NaN 将比任何其他值更高，Inf 低于它，然后是有限值。所以我们基本上免费得到一个NaN-propagatingmax-finder
在 KNL（骑士登陆）上，FP vmaxpd 和 vcmppd 具有与其整数等价物相同的性能：2 个周期延迟，0.5c 吞吐量。（https://agner.org/optimize/）。所以你的方式在 KNL 上的优势为零，但对于主流英特尔来说这是一个巧妙的技巧，比如 Skylake-X 和 IceLake。
错误修复优化版本：
使用 size_t 作为 return 类型和循环计数器/索引来处理潜在的巨大数组，而不是随机混合 int 和 64 位向量元素。 （uint64_t 对于收集 horizontal-max 的临时数组：即使在具有 32 位指针的构建中它也始终是 64 位的/size_t。）
修正：return 0 on N==1, not 1: the only element is 0.
错误修复：return -1 on N%8 != 0，而不是从 non-void 函数的末尾掉落。 （如果调用者在 C 中使用结果，或者在 C++ 中一旦执行结束，则为未定义行为）。
错误修复：FP 值的绝对值 = 清除符号位，而不是 bit-pattern
上的 2 的补码绝对值
某种错误修复：使用无符号整数比较和最大值，因此它适用于具有 _mm512_abs_epi64 的 2 的补码整数（产生无符号结果；记住 -LONG_MIN 溢出到 LONG_MIN 如果您继续将其视为已签名）。
样式改进：if (N%8 != 0) return -1;而不是将主体的大部分放在 if 块中。
风格改进：在首次使用时声明变量，并删除一些未使用的纯噪声。自 C99 以来，这是 C 的惯用语，它在 20 多年前被标准化。
样式改进：对仅包含加载结果的 tmp 向量变量使用更简单的名称。有时您只需要一个 tmp var，因为内部名称太长以至于您不想键入 _mm...load... 作为另一个内部函数的参数。像 v 这样的名称作用于内部循环是一个明确的标志，它只是一个占位符，以后不会使用。 （当你在初始化时声明它时，这种风格效果最好，所以很容易看出它不能在外部作用域中使用。）
优化：SIMD 循环后减少 8 -> 4 个元素：提取高半部分，与现有的低半部分结合。 （与您对更简单的水平缩减（如 sum 或 max）所做的操作相同）。当我们需要的指令只有AVX512有，而KNL没有AVX512VL时很不方便，只好用512位版本，忽略高位垃圾。但是 KNL 确实有 AVX1 / AVX2 所以我们仍然可以存储 256 位向量并做一些事情。
使用 merge-masking _mm512_mask_extracti64x4_epi64 提取直接混合同一向量的高半部分和低半部分是一个很酷的技巧，如果您使用 512 位 mask-blend。 :P
某种错误修复：在 C 中，main 在托管实现中具有 return 类型的 int（运行ning 在 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;
}
On Godbolt：GCC10.2 -O3 -march=knl 的主循环是 8 条指令。因此，即使（最好的情况）KNL 可以在 2/clock 解码并 运行 它，每个向量仍然需要 4 个周期。您可以 运行 Godbolt 上的程序；一世运行s 在 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
循环的另一个选项是屏蔽 vpbroadcastq zmm3{k1}, rax，在循环后添加 [0..7] per-element 偏移量。这实际上会在循环中保存 vpaddq，如果 GCC 无论如何要使用索引的 addressing-mode，我们在寄存器中有正确的 i。 （这对 Skylake-X 不利；击败 memory-source vpandd 的 micro-fusion。）
Agner Fog 没有列出 GP->vector 广播的性能，但希望它至少在 KNL 上只有 single-uop。 （并且 https://uops.info/ 没有 KNL 或 KNM 结果）。
分支策略：当新的最大值非常罕见时
如果您希望找到一个新的最大值非常罕见（例如，数组很大且分布均匀，或者至少没有向上趋势），广播当前最大值和分支找到任何更大的向量可能会更快元素.
找到一个新的最大值意味着跳出循环（这可能会预测错误，所以速度很慢）并广播该元素（可能使用 tzcnt 来查找元素索引，然后使用 broadcast-load ，并更新索引）。
特别是使用 KNL 的 4 路 SMT 来隐藏分支未命中成本，这可能是大型阵列的整体吞吐量胜利；平均每个元素的指令更少。
但对于 做 趋势向上的输入可能更糟，因此平均找到新的最大值 O(n) 次，而不是 sqrt(n) 或 log(n)或者均匀分布给我们的任何频率。
PS：打印向量，存储到数组并重新加载元素。 print a __m128i variable
或者使用调试器向您展示它们的元素。
我是使用 AVX512 指令编码的新手。我的机器是英特尔 KNL 7250。我正在尝试使用 AVX512 指令来查找具有最大绝对值的元素的索引，该元素的双精度和数组大小 % 8 = 0。但它每次都会打印输出索引 = 0。不知道哪里出了问题，求大神指点。 另外，如何将 printf 用于 __m512i 类型？
谢谢。
代码：
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,8,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, maxindices, 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);
}
}
您的函数 returns 0 因为您将 int64 索引视为 double 的 bit-pattern 并将其转换（ tiny) 数字到整数。 double indices_X[8]; 是错误；应该是uint64_t。 还有其他bug，见下文。
如果您在使用变量时声明变量，C99 风格，而不是过时的 C89 风格，则更容易发现此错误。
你 _mm512_storeu_si512 int64_t 的向量索引到 double indices_X[8]，type-punning 它加倍，然后在普通 C 中做 int maxindex = indices_X[0];。这是隐式的 type-conversion，将次正规 double 转换为整数。
(我注意到在 asm https://godbolt.org/z/zsfc36 中有一个神秘的 vcvttsd2si FP->int 转换，同时将代码转换为初始化器旁边的 C99 样式变量声明。这是一个线索：应该没有这个函数中的 FP->int 转换。我注意到大约在同一时间我将 double indices_X[8]; 声明向下移动到使用它的块中，并注意到它的类型为 double。）
实际上可以在 FP bit-patterns
上使用整数运算但前提是你用对了！ IEEE754 指数偏差意味着编码 / bit-pattern 可以作为 sign/magnitude 整数进行比较。 所以你可以做 abs / min / max 并比较它，但是当然不是整数加法/减法（除非你正在实施nextafter）。
_mm512_abs_epi64是2的补abs，不是sign-magnitude。相反，您必须屏蔽掉符号位。然后您就可以将结果视为无符号整数或有符号 2's-complement。 （两者都有效，因为高位已清除。）
使用整数 max 有一个有趣的 属性，即 NaN 将比任何其他值更高，Inf 低于它，然后是有限值。所以我们基本上免费得到一个NaN-propagatingmax-finder
在 KNL（骑士登陆）上，FP vmaxpd 和 vcmppd 具有与其整数等价物相同的性能：2 个周期延迟，0.5c 吞吐量。（https://agner.org/optimize/）。所以你的方式在 KNL 上的优势为零，但对于主流英特尔来说这是一个巧妙的技巧，比如 Skylake-X 和 IceLake。
错误修复优化版本：
使用
size_t作为 return 类型和循环计数器/索引来处理潜在的巨大数组，而不是随机混合int和 64 位向量元素。 （uint64_t对于收集 horizontal-max 的临时数组：即使在具有 32 位指针的构建中它也始终是 64 位的/size_t。）修正：return
0on N==1, not 1: the only element is 0.错误修复：return
-1onN%8 != 0，而不是从 non-void 函数的末尾掉落。 （如果调用者在 C 中使用结果，或者在 C++ 中一旦执行结束，则为未定义行为）。错误修复：FP 值的绝对值 = 清除符号位，而不是 bit-pattern
上的 2 的补码绝对值某种错误修复：使用无符号整数比较和最大值，因此它适用于具有
_mm512_abs_epi64的 2 的补码整数（产生无符号结果；记住-LONG_MIN溢出到LONG_MIN如果您继续将其视为已签名）。样式改进：
if (N%8 != 0) return -1;而不是将主体的大部分放在 if 块中。风格改进：在首次使用时声明变量，并删除一些未使用的纯噪声。自 C99 以来，这是 C 的惯用语，它在 20 多年前被标准化。
样式改进：对仅包含加载结果的 tmp 向量变量使用更简单的名称。有时您只需要一个 tmp var，因为内部名称太长以至于您不想键入
_mm...load...作为另一个内部函数的参数。像v这样的名称作用于内部循环是一个明确的标志，它只是一个占位符，以后不会使用。 （当你在初始化时声明它时，这种风格效果最好，所以很容易看出它不能在外部作用域中使用。）优化：SIMD 循环后减少 8 -> 4 个元素：提取高半部分，与现有的低半部分结合。 （与您对更简单的水平缩减（如 sum 或 max）所做的操作相同）。当我们需要的指令只有AVX512有，而KNL没有AVX512VL时很不方便，只好用512位版本，忽略高位垃圾。但是 KNL 确实有 AVX1 / AVX2 所以我们仍然可以存储 256 位向量并做一些事情。
使用 merge-masking
_mm512_mask_extracti64x4_epi64提取直接混合同一向量的高半部分和低半部分是一个很酷的技巧，如果您使用 512 位 mask-blend。 :P某种错误修复：在 C 中，
main在托管实现中具有 return 类型的int（运行ning 在 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;
}
On Godbolt：GCC10.2 -O3 -march=knl 的主循环是 8 条指令。因此，即使（最好的情况）KNL 可以在 2/clock 解码并 运行 它，每个向量仍然需要 4 个周期。您可以 运行 Godbolt 上的程序；一世运行s 在 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
循环的另一个选项是屏蔽 vpbroadcastq zmm3{k1}, rax，在循环后添加 [0..7] per-element 偏移量。这实际上会在循环中保存 vpaddq，如果 GCC 无论如何要使用索引的 addressing-mode，我们在寄存器中有正确的 i。 （这对 Skylake-X 不利；击败 memory-source vpandd 的 micro-fusion。）
Agner Fog 没有列出 GP->vector 广播的性能，但希望它至少在 KNL 上只有 single-uop。 （并且 https://uops.info/ 没有 KNL 或 KNM 结果）。
分支策略：当新的最大值非常罕见时
如果您希望找到一个新的最大值非常罕见（例如，数组很大且分布均匀，或者至少没有向上趋势），广播当前最大值和分支找到任何更大的向量可能会更快元素.
找到一个新的最大值意味着跳出循环（这可能会预测错误，所以速度很慢）并广播该元素（可能使用 tzcnt 来查找元素索引，然后使用 broadcast-load ，并更新索引）。
特别是使用 KNL 的 4 路 SMT 来隐藏分支未命中成本，这可能是大型阵列的整体吞吐量胜利；平均每个元素的指令更少。
但对于 做 趋势向上的输入可能更糟，因此平均找到新的最大值 O(n) 次，而不是 sqrt(n) 或 log(n)或者均匀分布给我们的任何频率。
PS：打印向量，存储到数组并重新加载元素。 print a __m128i variable
或者使用调试器向您展示它们的元素。