SIMD로 컬럼 별 최대 값 최적화
이 함수는 코드에서 상당한 시간을 보냈으며 가능한 경우 벡터화 -SIMD- 컴파일러 내장 함수로 최적화하고 싶습니다.
기본적으로 열에 대한 행렬에서 최대 값과 최대 위치를 찾아 저장합니다.
- val_ptr : 입력 행렬 : column-major (Fortran 스타일) n_rows-by-n_cols (일반적으로 n_rows >> n_cols)
- opt_pos_ptr : 최대 위치를 저장할 n_rows 길이의 int 벡터. 0으로 채워진 항목.
- max_ptr : 최대 값을 저장할 n_rows 길이의 부동 벡터. val_ptr의 첫 번째 열 사본으로 채워진 항목
- 함수는 병렬 루프에서 호출됩니다.
- 메모리 영역은 겹치지 않도록 보장됩니다.
- 나는 실제로 max_ptr을 채울 필요가 없으며 현재는 단지 부기 및 메모리 할당을 피하기 위해 사용됩니다.
- Windows 10에서 MSVC, C ++ 17을 사용합니다. 최신 Intel 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 병렬을 사용하려고 시도했지만 현재 사용보다 약간 큰 매우 큰 행에서만 무언가를 가져옵니다.
- if in inner 루프는 #pragma omp simd가 작동하는 것을 막았고 그것 없이는 다시 쓸 수 없었습니다.
답변
게시 한 코드 샘플에 따르면 수직 최대 값을 계산하려는 것 같습니다. 즉, "열"이 수평임을 의미합니다. C / C ++에서 요소의 수평 시퀀스 (즉, 인접한 두 요소가 메모리에서 한 요소의 거리를 가짐)는 일반적으로 행과 수직 (두 개의 인접한 요소가 메모리에서 행 크기의 거리를 가짐)-열이라고합니다. 아래의 내 대답에서는 행이 수평이고 열이 수직 인 전통적인 용어를 사용합니다.
또한 간결함을 위해 가능한 한 가지 유형의 행렬 요소 인 float. 기본 아이디어는에서도 동일 double하지만 주요 차이점은 벡터 당 요소 수와 _ps/ _pd내장 함수 선택입니다. double마지막에 버전을 제공하겠습니다 .
아이디어는 _mm_max_ps/를 사용하여 병렬로 여러 열에 대한 수직 최대 값을 계산할 수 있다는 것 _mm_max_pd입니다. 발견 된 최대 값의 위치도 기록하기 위해 이전 최대 값을 현재 요소와 비교할 수 있습니다. 비교 결과는 최대 값이 업데이트되는 요소가 모두 1 인 마스크입니다. 이 마스크를 사용하여 업데이트해야하는 위치를 선택할 수도 있습니다.
아래 알고리즘은 열에 동일한 최대 요소가 여러 개있는 경우 최대 요소의 위치가 기록되는 것이 중요하지 않다고 가정합니다. 또한 행렬에 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;
}
}
The code above requires SSE4.1 because of the blending intrinsics. You can replace those with a combination of _mm_and_si128/_ps, _mm_andnot_si128/_ps and _mm_or_si128/_ps, in which case the requirements will be lowered to SSE2. See Intel Intrinsics Guide for more details on the particular intrinsics, including which instruction set extensions they require.
NaN 값에 대한 참고 사항. 행렬에 NaN이있을 수있는 경우 _mm_cmplt_ps테스트는 항상 false를 반환합니다. 의 _mm_max_ps경우 일반적으로 무엇을 반환할지 알 수 없습니다. maxps피연산자 중 하나에 따라서 명령의 오퍼랜드를 배치함으로써, NaN의 경우 복귀에 대한 극한를 번역 번째 (소스)는 피연산자가 어느 동작을 달성 할 수 있음을 지시. 그러나 _mm_max_ps내장 함수 의 어떤 인수가 명령어의 피연산자를 나타내는지는 문서화되어 있지 않으며 컴파일러가 다른 경우에 다른 연관을 사용할 수도 있습니다. 자세한 내용은 이 답변을 참조하십시오.
올바른 동작을 보장하기 위해 wrt. NaN은 인라인 어셈블러를 사용하여 maxps피연산자 의 올바른 순서를 적용 할 수 있습니다 . 불행히도 x86-64 대상 용 MSVC의 옵션이 아니므로 대신 다음 _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));
더 넓은 벡터 ( __m256또는 __m512) 를 사용 하고 외부 루프를 작은 요소로 풀면 위의 알고리즘 성능을 더욱 향상시킬 수 있으므로 내부 루프가 반복 될 때마다 최소한 행 데이터의 캐시 라인 가치가로드됩니다.
다음은 double. 여기서 주목해야 할 중요한 점 double은 벡터 당 두 개의 요소 만 있고 벡터 당 여전히 4 개의 위치 가 있기 때문에 외부 루프를 펼쳐서 double한 번에 두 개의 벡터를 처리 한 다음 두 개의 마스크를 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;
}
}