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

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); });
}
スポンサーリンク

フォローする

スポンサーリンク