Menggunakan utas OpenMP dan std: (percobaan: :) simd untuk menghitung set Mandelbrot
Saya ingin menerapkan plotter set Mandelbrot sederhana menggunakan berbagai jenis paradigma HPC, menunjukkan kekuatan dan kelemahan mereka dan betapa mudah atau sulitnya implementasinya. Pikirkan GPGPU (CUDA / OpenACC / OpenMP4.5), threading / OpenMP dan MPI. Dan gunakan contoh-contoh ini untuk memberikan pemrogram yang baru mengenal HPC dan untuk melihat apa saja kemungkinannya. Kejelasan kode lebih penting daripada mendapatkan kinerja teratas absolut dari perangkat keras, itulah langkah kedua;)
Karena masalahnya sepele untuk diparalelkan dan CPU modern dapat memperoleh kinerja yang sangat besar menggunakan instruksi vektor, saya juga ingin menggabungkan OpenMP dan SIMD. Sayangnya, hanya menambahkan #pragma omp simd
tidak menghasilkan hasil yang memuaskan dan menggunakan intrinsik sangat tidak ramah pengguna atau bukti masa depan. Atau cantik .
Untungnya, pekerjaan sedang dilakukan untuk standar C ++ sehingga akan lebih mudah untuk mengimplementasikan instruksi vektor secara umum, seperti yang disebutkan di TS: "Ekstensi untuk paralelisme, versi 2" , khususnya bagian 9 tentang tipe data-paralel. Implementasi WIP dapat ditemukan di sini , yang didasarkan pada VC yang dapat ditemukan di sini .
Asumsikan bahwa saya memiliki kelas berikut (yang telah diubah agar sedikit lebih sederhana)
#include <stddef.h>
using Range = std::pair<double, double>;
using Resolution = std::pair<std::size_t, std::size_t>;
class Mandelbrot
{
double* d_iters;
Range d_xrange;
Range d_yrange;
Resolution d_res;
std::size_t d_maxIter;
public:
Mandelbrot(Range xrange, Range yrange, Resolution res, std::size_t maxIter);
~Mandelbrot();
void writeImage(std::string const& fileName);
void computeMandelbrot();
private:
void calculateColors();
};
Dan implementasi berikut computeMandelbrot()
menggunakan OpenMP
void Mandelbrot::computeMandelbrot()
{
double dx = (d_xrange.second - d_xrange.first) / d_res.first;
double dy = (d_yrange.second - d_yrange.first) / d_res.second;
#pragma omp parallel for schedule(dynamic)
for (std::size_t row = 0; row != d_res.second; ++row)
{
double c_imag = d_yrange.first + row * dy;
for (std::size_t col = 0; col != d_res.first; ++col)
{
double real = 0.0;
double imag = 0.0;
double realSquared = 0.0;
double imagSquared = 0.0;
double c_real = d_xrange.first + col * dx;
std::size_t iter = 0;
while (iter < d_maxIter && realSquared + imagSquared < 4.0)
{
realSquared = real * real;
imagSquared = imag * imag;
imag = 2 * real * imag + c_imag;
real = realSquared - imagSquared + c_real;
++iter;
}
d_iters[row * d_res.first + col] = iter;
}
}
}
Kita dapat mengasumsikan bahwa resolusi kedua arah x dan y adalah kelipatan 2/4/8 / .., bergantung pada instruksi SIMD yang kita gunakan.
Sayangnya, hanya ada sedikit informasi yang tersedia secara online std::experimental::simd
. Juga tidak ada contoh non-sepele sejauh yang saya bisa temukan.
Di repositori Vc git, terdapat implementasi kalkulator set Mandelbrot, tetapi cukup berbelit-belit dan karena kurangnya komentar agak sulit diikuti.
Jelas bahwa saya harus mengubah tipe data ganda dalam fungsi computeMandelbrot()
, tapi saya tidak yakin apa. TS menyebutkan dua tipe data baru utama untuk beberapa tipe T,
native_simd = std::experimental::simd<T, std::experimental::simd_abi::native>;
dan
fixed_size_simd = std::experimental::simd<T, std::experimental::simd_abi::fixed_size<N>>;
Menggunakan native_simd
paling masuk akal, karena saya tidak tahu batasan saya pada waktu kompilasi. Tetapi kemudian tidak jelas bagi saya apa yang diwakili oleh tipe-tipe ini, apakah native_simd<double>
double tunggal atau kumpulan ganda di mana instruksi vektor dijalankan? Lalu berapa banyak ganda dalam koleksi ini?
Jika seseorang dapat menunjukkan kepada saya contoh di mana konsep ini digunakan, atau memberi saya beberapa petunjuk tentang cara mengimplementasikan instruksi vektor menggunakan std :: eksperimental :: simd, saya akan sangat berterima kasih.
Jawaban
Berikut adalah implementasi yang sangat mendasar, yang berfungsi (sejauh yang saya tahu). Pengujian elemen vektor mana yang memiliki nilai absolut lebih besar dari 2 dilakukan dengan cara yang sangat rumit dan tidak efisien. Pasti ada cara yang lebih baik untuk melakukan ini, tapi saya belum menemukannya.
Saya mendapatkan peningkatan kinerja sekitar 72% pada AMD Ryzen 5 3600 dan memberikan g ++ opsi -march=znver2
, yang kurang dari yang diharapkan.
template <class T>
void mandelbrot(T xstart, T xend,
T ystart, T yend)
{
namespace stdx = std::experimental;
constexpr auto simdSize = stdx::native_simd<T>().size();
constexpr unsigned size = 4096;
constexpr unsigned maxIter = 250;
assert(size % simdSize == 0);
unsigned* res = new unsigned[size * size];
T dx = (xend - xstart) / size;
T dy = (yend - ystart) / size;
for (std::size_t row = 0; row != size; ++row)
{
T c_imag = ystart + row * dy;
for (std::size_t col = 0; col != size; col += simdSize)
{
stdx::native_simd<T> real{0};
stdx::native_simd<T> imag{0};
stdx::native_simd<T> realSquared{0};
stdx::native_simd<T> imagSquared{0};
stdx::fixed_size_simd<unsigned, simdSize> iters{0};
stdx::native_simd<T> c_real;
for (int idx = 0; idx != simdSize; ++idx)
{
c_real[idx] = xstart + (col + idx) * dx;
}
for (unsigned iter = 0; iter != maxIter; ++iter)
{
realSquared = real * real;
imagSquared = imag * imag;
auto isInside = realSquared + imagSquared > stdx::native_simd<T>{4};
for (int idx = 0; idx != simdSize; ++idx)
{
// if not bigger than 4, increase iters
if (!isInside[idx])
{
iters[idx] += 1;
}
else
{
// prevent that they become inf/nan
real[idx] = static_cast<T>(4);
imag[idx] = static_cast<T>(4);
}
}
if (stdx::all_of(isInside) )
{
break;
}
imag = static_cast<T>(2.0) * real * imag + c_imag;
real = realSquared - imagSquared + c_real;
}
iters.copy_to(res + row * size + col, stdx::element_aligned);
}
}
delete[] res;
}
Seluruh kode pengujian (mulai dari auto test = (...)
) dikompilasi hingga
.L9:
vmulps ymm1, ymm1, ymm1
vmulps ymm13, ymm2, ymm2
xor eax, eax
vaddps ymm2, ymm13, ymm1
vcmpltps ymm2, ymm5, ymm2
vmovaps YMMWORD PTR [rsp+160], ymm2
jmp .L6
.L3:
vmovss DWORD PTR [rsp+32+rax], xmm0
vmovss DWORD PTR [rsp+64+rax], xmm0
add rax, 4
cmp rax, 32
je .L22
.L6:
vucomiss xmm3, DWORD PTR [rsp+160+rax]
jp .L3
jne .L3
inc DWORD PTR [rsp+96+rax]
add rax, 4
cmp rax, 32
jne .L6