お馬の写真 別館

「お馬の写真」管理者による徒然なるブログ

マンデルブロ集合: SIMDで更に高速化

      2015/07/16

mandelbrot

前の記事で行なったtbb::parallel_forに加えて、SIMDによる高速化にもチャレンジ。

まず先にハマった点を。
今の開発環境はOSX上でVirtualBox 4.3.28を動かし、その上でWindows 7を動かす、という形。
ホストマシン側ではAVX命令が使えるCPUなのですが、VirtualBox上のWindows 7からはAVX命令使えず。ということでSSSE3までの命令を使う形にしています。(コード上に出てくる命令的にはSSE2まで…のつもりだったけど、AVX命令の_mm_cmp_pdが動作してますね。よくわからん)

参考とさせてもらったのはこちらのページ。あとSIMD命令に関してはIntel Intrinsics Guideが参考になります。

ということでコードをどうぞ。

#include "stdafx.h"
#include "MandelbrotCppLogicSimd.h"
#include <tbb/parallel_for.h>
#include <tmmintrin.h>
using namespace std;
using namespace tbb;
class MandelbrotCppLogicSimd::Impl {
public:
static void CalcPoint(int threshold, const complex<double> *c, int *result)
{
auto n = _mm_setzero_pd();
auto dn = _mm_setr_pd(-1, -1);
auto divergence = _mm_setr_pd(4, 4);
auto cr = _mm_setr_pd(c[0].real(), c[1].real());
auto ci = _mm_setr_pd(c[0].imag(), c[1].imag());
auto zr = _mm_setzero_pd();
auto zi = _mm_setzero_pd();
for (auto i = 0; i < threshold; ++i)
{
auto zr2 = _mm_mul_pd(zr, zr);
auto zi2 = _mm_mul_pd(zi, zi);
auto zrzi = _mm_mul_pd(zr, zi);
zr = _mm_add_pd(_mm_sub_pd(zr2, zi2), cr);
zi = _mm_add_pd(_mm_add_pd(zrzi, zrzi), ci);
auto z2 = _mm_add_pd(zr2, zi2);
auto mask = _mm_cmp_pd(z2, divergence, _CMP_LE_OS);
auto flag = _mm_movemask_pd(mask);
if (!flag) {
break;
}
zr = _mm_and_pd(zr, mask);
zi = _mm_and_pd(zi, mask);
cr = _mm_and_pd(cr, mask);
ci = _mm_and_pd(ci, mask);
dn = _mm_and_pd(dn, mask);
n = _mm_sub_pd(n, dn);
}
__declspec(align(16)) double resultDouble[SimdWidth];
_mm_store_pd(resultDouble, n);
for (auto i = 0; i < SimdWidth; ++i) {
auto item = static_cast<int>(resultDouble[i]);
if (item == threshold) {
item = -1;
}
result[i] = item;
}
return;
}
static const int SimdWidth = 2;
private:
Impl();
~Impl();
};
void MandelbrotCppLogicSimd::Calc(
vector<int> *result,
const complex<double> &center, double step,
int width, int height, int threshold)
{
if (result == nullptr) {
return;
}
auto vectorSize = width * height;
auto padding = Impl::SimdWidth - vectorSize % Impl::SimdWidth;
if (padding == Impl::SimdWidth) {
padding = 0;
}
result->resize(vectorSize + padding);
vector<complex<double>> input(vectorSize + padding);
complex<double> c(0, center.imag() - step * height / 2.0);
int index = 0;
for (auto y = 0; y < height; ++y)
{
c.real(center.real() - step * width / 2.0);
for (auto x = 0; x < width; ++x)
{
input[index++] = c;
c.real(c.real() + step);
}
c.imag(c.imag() + step);
}
if (padding) {
c.real(2);
c.imag(2);
for (auto i = 0; i < padding; ++i)
{
input[vectorSize + i] = c;
}
}
parallel_for(0, index, Impl::SimdWidth,
[&](int i) { Impl::CalcPoint(threshold, input.data() + i, result->data() + i); });
}

 - プログラミング