我在代码中花费了大量时间来实现此功能，如果可能的话，我想通过矢量化-SIMD-编译器内部函数对其进行优化。
它本质上是在列矩阵中找到最大值和最大值的位置，并将它们存储起来:
- val_ptr:输入矩阵:列优先(Fortran 风格)n_rows-by-n_cols(通常是 n_rows>>n_cols)
- opt_pos_ptr:长度为 n_rows 的 int vector ，用于存储最大值的位置。在条目中填入零。
- max_ptr:长度为 n_rows 的浮点 vector ，用于存储最大值。在条目中填充了 val_ptr 第一列的拷贝
- 该函数将在并行循环中调用
- 保证内存区域不重叠
- 我真的不需要填充 max_ptr，目前它只是用于簿记和避免内存分配
- 我在 Windows 10 上使用 MSVC、C++17。旨在运行现代英特尔 CPU
模板类型为 float 或 double 的代码:
template <typename eT>
find_max(const int n_cols,
const int n_rows,
const eT* val_ptr,
int* opt_pos_ptr,
eT* max_ptr){
for (int col = 1; col < n_cols; ++col)
{
//Getting the pointer to the beginning of the column
const auto* value_col = val_ptr + col * n_rows;
//Looping over the rows
for (int row = 0; row < n_rows; ++row)
{
//If the value is larger than the current maximum, we replace and we store its positions
if (value_col[row] > max_ptr[row])
{
max_ptr[row] = value_col[row];
opt_pos_ptr[row] = col;
}
}
}
}
到目前为止我尝试了什么:
- 我尝试在内部循环中使用 OpenMP parallel for，但只在非常大的行上带来一些东西，比我当前的使用量大一点。
- 内部循环中的 if 阻止了#pragma omp simd 的工作，没有它我无法重写它。
最佳答案
根据您发布的代码示例，您似乎想要计算垂直最大值，这意味着在您的情况下“列”是水平的。在 C/C++ 中，元素的水平序列(即两个相邻元素在内存中的距离为一个元素的距离)通常称为行和垂直(其中两个相邻元素在内存中的距离为行大小)-列。在下面的回答中，我将使用传统术语，其中行是水平的，列是垂直的。
此外，为简洁起见，我将重点关注一种可能的矩阵元素类型 - float . double 的基本思想相同，主要区别在于每个 vector 的元素数量和 _ps/_pd内在选择。我将为 double 提供一个版本在最后。
想法是您可以使用 _mm_max_ps 并行计算多列的垂直最大值/_mm_max_pd .为了同时记录找到的最大值的位置，您可以将先前的最大值与当前元素进行比较。比较的结果是一个掩码，其中元素是更新最大值的全一。该掩码也可用于选择需要更新的位置。
我必须注意，如果一列中有多个相等的最大元素，下面的算法假设记录哪个最大元素的位置并不重要。另外，我假设矩阵不包含 NaN 值，这会影响比较。稍后会详细介绍。
void find_max(const int n_cols,
const int n_rows,
const float* val_ptr,
int* opt_pos_ptr,
float* max_ptr){
const __m128i mm_one = _mm_set1_epi32(1);
// Pre-compute the number of rows that can be processed in full vector width.
// In a 128-bit vector there are 4 floats or 2 doubles
int tail_size = n_rows & 3;
int n_rows_aligned = n_rows - tail_size;
int row = 0;
for (; row < n_rows_aligned; row += 4)
{
const auto* col_ptr = val_ptr + row;
__m128 mm_max = _mm_loadu_ps(col_ptr);
__m128i mm_max_pos = _mm_setzero_si128();
__m128i mm_pos = mm_one;
col_ptr += n_rows;
for (int col = 1; col < n_cols; ++col)
{
__m128 mm_value = _mm_loadu_ps(col_ptr);
// See if this value is greater than the old maximum
__m128 mm_mask = _mm_cmplt_ps(mm_max, mm_value);
// If it is, save its position
mm_max_pos = _mm_blendv_epi8(mm_max_pos, mm_pos, _mm_castps_si128(mm_mask));
// Compute the maximum
mm_max = _mm_max_ps(mm_value, mm_max);
mm_pos = _mm_add_epi32(mm_pos, mm_one);
col_ptr += n_rows;
}
// Store the results
_mm_storeu_ps(max_ptr + row, mm_max);
_mm_storeu_si128(reinterpret_cast< __m128i* >(opt_pos_ptr + row), mm_max_pos);
}
// Process tail serially
for (; row < n_rows; ++row)
{
const auto* col_ptr = val_ptr + row;
auto max = *col_ptr;
int max_pos = 0;
col_ptr += n_rows;
for (int col = 1; col < n_cols; ++col)
{
auto value = *col_ptr;
if (value > max)
{
max = value;
max_pos = col;
}
col_ptr += n_rows;
}
max_ptr[row] = max;
opt_pos_ptr[row] = max_pos;
}
}
由于混合内在函数，上面的代码需要 SSE4.1。您可以将它们替换为 _mm_and_si128 的组合/_ps , _mm_andnot_si128/_ps和 _mm_or_si128/_ps ，在这种情况下，要求将降低到 SSE2。参见 Intel Intrinsics Guide有关特定内在函数的更多详细信息，包括它们需要哪些指令集扩展。
关于 NaN 值的注释。如果您的矩阵可以有 NaN，则 _mm_cmplt_ps测试将始终返回 false。至于_mm_max_ps ，通常不知道它会返回什么。 maxps如果其中一个操作数是 NaN，则内在转换为的指令返回其第二个(源)操作数，因此通过排列指令的操作数，您可以实现任一行为。但是，没有记录 _mm_max_ps 的哪个参数intrinsic 表示指令的哪个操作数，编译器甚至可能在不同情况下使用不同的关联。参见 this回答更多详情。
为了确保正确的行为。 NaNs 你可以使用内联汇编器来强制 maxps 的正确顺序操作数。不幸的是，这不是 MSVC for x86-64 target 的一个选项，你说你正在使用它，所以你可以重用 _mm_cmplt_ps第二次混合的结果如下:
// Compute the maximum
mm_max = _mm_blendv_ps(mm_max, mm_value, mm_mask);
这将抑制结果最大值中的 NaN。如果您想保留 NaN，则可以使用第二次比较来检测 NaN:
// Detect NaNs
__m128 mm_nan_mask = _mm_cmpunord_ps(mm_value, mm_value);
// Compute the maximum
mm_max = _mm_blendv_ps(mm_max, mm_value, _mm_or_ps(mm_mask, mm_nan_mask));
如果您使用更宽的 vector (__m256 或 __m512)并以小因子展开外循环，那么您可能会进一步提高上述算法的性能，以便至少加载一个缓存行值的行数据在内循环的每次迭代中。
这是 double 的实现示例.这里要注意的重点是因为只有两个double每个 vector 的元素并且每个 vector 仍然有四个位置，我们必须展开外循环来处理 double 的两个 vector 一次，然后压缩与先前最大值比较的两个掩码以混合 32 位位置。
void find_max(const int n_cols,
const int n_rows,
const double* val_ptr,
int* opt_pos_ptr,
double* max_ptr){
const __m128i mm_one = _mm_set1_epi32(1);
// Pre-compute the number of rows that can be processed in full vector width.
// In a 128-bit vector there are 2 doubles, but we want to process
// two vectors at a time.
int tail_size = n_rows & 3;
int n_rows_aligned = n_rows - tail_size;
int row = 0;
for (; row < n_rows_aligned; row += 4)
{
const auto* col_ptr = val_ptr + row;
__m128d mm_max1 = _mm_loadu_pd(col_ptr);
__m128d mm_max2 = _mm_loadu_pd(col_ptr + 2);
__m128i mm_max_pos = _mm_setzero_si128();
__m128i mm_pos = mm_one;
col_ptr += n_rows;
for (int col = 1; col < n_cols; ++col)
{
__m128d mm_value1 = _mm_loadu_pd(col_ptr);
__m128d mm_value2 = _mm_loadu_pd(col_ptr + 2);
// See if this value is greater than the old maximum
__m128d mm_mask1 = _mm_cmplt_pd(mm_max1, mm_value1);
__m128d mm_mask2 = _mm_cmplt_pd(mm_max2, mm_value2);
// Compress the 2 masks into one
__m128i mm_mask = _mm_packs_epi32(
_mm_castpd_si128(mm_mask1), _mm_castpd_si128(mm_mask2));
// If it is, save its position
mm_max_pos = _mm_blendv_epi8(mm_max_pos, mm_pos, mm_mask);
// Compute the maximum
mm_max1 = _mm_max_pd(mm_value1, mm_max1);
mm_max2 = _mm_max_pd(mm_value2, mm_max2);
mm_pos = _mm_add_epi32(mm_pos, mm_one);
col_ptr += n_rows;
}
// Store the results
_mm_storeu_pd(max_ptr + row, mm_max1);
_mm_storeu_pd(max_ptr + row + 2, mm_max2);
_mm_storeu_si128(reinterpret_cast< __m128i* >(opt_pos_ptr + row), mm_max_pos);
}
// Process 2 doubles at once
if (tail_size >= 2)
{
const auto* col_ptr = val_ptr + row;
__m128d mm_max1 = _mm_loadu_pd(col_ptr);
__m128i mm_max_pos = _mm_setzero_si128();
__m128i mm_pos = mm_one;
col_ptr += n_rows;
for (int col = 1; col < n_cols; ++col)
{
__m128d mm_value1 = _mm_loadu_pd(col_ptr);
// See if this value is greater than the old maximum
__m128d mm_mask1 = _mm_cmplt_pd(mm_max1, mm_value1);
// Compress the mask. The upper half doesn't matter.
__m128i mm_mask = _mm_packs_epi32(
_mm_castpd_si128(mm_mask1), _mm_castpd_si128(mm_mask1));
// If it is, save its position
mm_max_pos = _mm_blendv_epi8(mm_max_pos, mm_pos, mm_mask);
// Compute the maximum
mm_max1 = _mm_max_pd(mm_value1, mm_max1);
mm_pos = _mm_add_epi32(mm_pos, mm_one);
col_ptr += n_rows;
}
// Store the results
_mm_storeu_pd(max_ptr + row, mm_max1);
// Only store the lower two positions
_mm_storel_epi64(reinterpret_cast< __m128i* >(opt_pos_ptr + row), mm_max_pos);
row += 2;
}
// Process tail serially
for (; row < n_rows; ++row)
{
const auto* col_ptr = val_ptr + row;
auto max = *col_ptr;
int max_pos = 0;
col_ptr += n_rows;
for (int col = 1; col < n_cols; ++col)
{
auto value = *col_ptr;
if (value > max)
{
max = value;
max_pos = col;
}
col_ptr += n_rows;
}
max_ptr[row] = max;
opt_pos_ptr[row] = max_pos;
}
}
关于c++ - 使用 SIMD 优化列式最大值，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/63424833/