mengoptimalkan maksimum secara kolom-bijaksana dengan SIMD
Saya memiliki fungsi ini di mana saya menghabiskan banyak waktu dalam kode saya, dan saya ingin mengoptimalkannya dengan vectorization-SIMD-compiler intrinsics, jika memungkinkan.
Ini pada dasarnya menemukan nilai dan lokasi maksimum di atas matriks di atas kolom, dan menyimpannya:
- val_ptr: matriks masukan: kolom-utama (gaya Fortran) n_rows-by-n_cols (biasanya n_rows >> n_cols)
- opt_pos_ptr: vektor int dengan panjang n_rows tempat menyimpan posisi maksimum. Pada entri diisi dengan angka nol.
- max_ptr: vektor float dengan panjang n_rows tempat menyimpan maksimum. Pada entri diisi dengan salinan kolom pertama val_ptr
- Fungsi tersebut akan dipanggil dalam loop paralel
- Wilayah memori dijamin tidak tumpang tindih
- Saya tidak terlalu membutuhkan max_ptr untuk diisi, saat ini hanya digunakan untuk pembukuan dan untuk menghindari alokasi memori
- Saya menggunakan MSVC, C ++ 17 pada Windows 10. Dimaksudkan untuk menjalankan CPU Intel modern
Kode, di mana jenis templat dimaksudkan sebagai float atau 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;
}
}
}
}
Apa yang saya coba sejauh ini:
- Saya mencoba menggunakan OpenMP paralel untuk di loop dalam, tetapi membawa sesuatu hanya pada baris yang sangat besar, sedikit lebih besar dari penggunaan saya saat ini.
- If in inner loop mencegah #pragma omp simd bekerja, dan saya tidak dapat menulis ulang tanpanya.
Jawaban
Berdasarkan contoh kode yang Anda posting, sepertinya Anda ingin menghitung nilai maksimum vertikal, artinya dalam kasus Anda "kolom" berbentuk horizontal. Dalam C / C ++ urutan horizontal elemen (yaitu di mana dua elemen yang berdekatan memiliki jarak satu elemen dalam memori) biasanya disebut baris dan vertikal (di mana dua elemen yang berdekatan memiliki jarak ukuran baris dalam memori) - kolom. Dalam jawaban saya di bawah ini, saya akan menggunakan terminologi tradisional, di mana baris horizontal dan kolom vertikal.
Juga, untuk singkatnya saya akan fokus pada satu kemungkinan jenis elemen matriks - float. Ide dasarnya sama double, dengan perbedaan utama adalah jumlah elemen per vektor dan pilihan _ps/ _pdintrinsik. Saya akan memberikan versi doubledi bagian akhir.
Idenya adalah Anda dapat menghitung maksimum vertikal untuk beberapa kolom secara paralel menggunakan _mm_max_ps/ _mm_max_pd. Untuk juga mencatat posisi maksimum yang ditemukan, Anda dapat membandingkan maksimum sebelumnya dengan elemen saat ini. Hasil perbandingannya adalah topeng, di mana elemen semuanya adalah yang maksimum diperbarui. Topeng itu dapat digunakan untuk memilih posisi mana yang perlu diperbarui juga.
Saya harus mencatat bahwa algoritma di bawah ini mengasumsikan bahwa tidak penting posisi elemen maks mana yang dicatat, jika ada beberapa elemen maks yang sama dalam sebuah kolom. Juga, saya menganggap matriks tidak mengandung nilai NaN, yang akan mempengaruhi perbandingan. Lebih lanjut tentang ini nanti.
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;
}
}
Kode di atas membutuhkan SSE4.1 karena intrinsik pencampuran. Anda dapat menggantinya dengan kombinasi _mm_and_si128/ _ps, _mm_andnot_si128/ _psdan _mm_or_si128/ _ps, dalam hal ini persyaratan akan diturunkan ke SSE2. Lihat Panduan Intrinsik Intel untuk detail lebih lanjut tentang intrinsik tertentu, termasuk ekstensi set instruksi mana yang mereka butuhkan.
Catatan tentang nilai NaN. Jika matriks Anda dapat memiliki NaN, _mm_cmplt_pspengujian akan selalu menghasilkan nilai salah. Adapun _mm_max_ps, umumnya tidak diketahui apa yang akan dikembalikannya. The maxpsinstruksi yang diterjemahkan intrinsik untuk pengembalian yang kedua (sumber) operan jika salah satu operan adalah NaN, sehingga dengan mengatur operan dari instruksi yang Anda dapat mencapai perilaku baik. Namun, tidak didokumentasikan argumen _mm_max_psintrinsik mana yang merepresentasikan operand mana dari instruksi tersebut, dan bahkan mungkin saja kompilator dapat menggunakan asosiasi yang berbeda dalam kasus yang berbeda. Lihat jawaban ini untuk lebih jelasnya.
Untuk memastikan perilaku yang benar wrt. NaN Anda dapat menggunakan assembler sebaris untuk memaksa urutan maxpsoperan yang benar. Sayangnya, itu bukan opsi dengan MSVC untuk target x86-64, yang Anda katakan sedang Anda gunakan, jadi Anda dapat menggunakan kembali _mm_cmplt_pshasilnya untuk campuran kedua seperti ini:
// Compute the maximum
mm_max = _mm_blendv_ps(mm_max, mm_value, mm_mask);
Ini akan menekan NaN dalam nilai maks yang dihasilkan. Jika Anda ingin menyimpan NaN, Anda dapat menggunakan perbandingan kedua untuk mendeteksi 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));
Anda mungkin dapat lebih meningkatkan kinerja algoritme di atas jika Anda menggunakan vektor yang lebih luas ( __m256atau __m512) dan membuka gulungan loop luar dengan faktor kecil, sehingga setidaknya data baris senilai baris cache dimuat pada setiap iterasi loop dalam.
Berikut adalah contoh implementasi untuk double. Hal penting yang perlu diperhatikan di sini adalah karena hanya ada dua doubleelemen per vektor dan masih ada empat posisi per vektor, kita harus membuka gulungan loop luar untuk memproses dua vektor doublesekaligus dan kemudian memampatkan dua topeng dari perbandingan dengan maksimum sebelumnya untuk memadukan posisi 32-bit.
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;
}
}