|
@@ -0,0 +1,204 @@
|
|
|
|
+#include <immintrin.h>
|
|
|
|
+#include <omp.h>
|
|
|
|
+
|
|
|
|
+#include <cmath>
|
|
|
|
+#include <vector>
|
|
|
|
+#include "CpuGenerators.h"
|
|
|
|
+
|
|
|
|
+#include "LightDoubleDouble.h"
|
|
|
|
+#include "QuadDouble.h"
|
|
|
|
+#include "HexDouble.h"
|
|
|
|
+
|
|
|
|
+class CpuGeneratorFloatAVXFMA : public mnd::MandelGenerator
|
|
|
|
+{
|
|
|
|
+public:
|
|
|
|
+ CpuGeneratorFloatAVXFMA(void) :
|
|
|
|
+ MandelGenerator{ mnd::Precision::FLOAT, mnd::X86_AVX_FMA }
|
|
|
|
+ {
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ virtual void generate(const mnd::MandelInfo& info, float* data) override;
|
|
|
|
+};
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+extern "C" const std::vector<mnd::MandelGenerator*>& mandel_get_generators(void)
|
|
|
|
+{
|
|
|
|
+ static CpuGeneratorFloatAVXFMA instance;
|
|
|
|
+ static std::vector<mnd::MandelGenerator*> vec { &instance };
|
|
|
|
+ return vec;
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+void CpuGeneratorFloatAVXFMA::generate(const mnd::MandelInfo& info, float* data)
|
|
|
|
+{
|
|
|
|
+ const bool parallel = true;
|
|
|
|
+
|
|
|
|
+ using T = float;
|
|
|
|
+
|
|
|
|
+ const auto& view = info.view;
|
|
|
|
+ const T vx = mnd::convert<T>(view.x);
|
|
|
|
+ const T vy = mnd::convert<T>(view.y);
|
|
|
|
+ const T vw = mnd::convert<T>(view.width);
|
|
|
|
+ const T vh = mnd::convert<T>(view.height);
|
|
|
|
+
|
|
|
|
+ const T jX = mnd::convert<T>(info.juliaX);
|
|
|
|
+ const T jY = mnd::convert<T>(info.juliaY);
|
|
|
|
+
|
|
|
|
+ const float dppf = float(vw / info.bWidth);
|
|
|
|
+ const float viewxf = vx;
|
|
|
|
+ __m256 viewx = { viewxf, viewxf, viewxf, viewxf, viewxf, viewxf, viewxf, viewxf };
|
|
|
|
+ __m256 dpp = { dppf, dppf, dppf, dppf, dppf, dppf, dppf, dppf };
|
|
|
|
+
|
|
|
|
+ __m256 juliaX = { jX, jX, jX, jX, jX, jX, jX, jX };
|
|
|
|
+ __m256 juliaY = { jY, jY, jY, jY, jY, jY, jY, jY };
|
|
|
|
+
|
|
|
|
+#if defined(_OPENMP)
|
|
|
|
+ if (parallel)
|
|
|
|
+ omp_set_num_threads(omp_get_num_procs());
|
|
|
|
+# pragma omp parallel for schedule(static, 1) if (parallel)
|
|
|
|
+#endif
|
|
|
|
+ for (long j = 0; j < info.bHeight; j++) {
|
|
|
|
+ T y = vy + T(j) * vw / info.bHeight;
|
|
|
|
+ __m256 ys = {y, y, y, y, y, y, y, y};
|
|
|
|
+ for (long i = 0; i < info.bWidth; i += 24) {
|
|
|
|
+ __m256 pixc = { float(i), float(i + 1), float(i + 2), float(i + 3), float(i + 4), float(i + 5), float(i + 6), float(i + 7) };
|
|
|
|
+ __m256 pixc2 = { float(i + 8), float(i + 9), float(i + 10), float(i + 11), float(i + 12), float(i + 13), float(i + 14), float(i + 15) };
|
|
|
|
+ __m256 pixc3 = { float(i + 16), float(i + 17), float(i + 18), float(i + 19), float(i + 20), float(i + 21), float(i + 22), float(i + 23) };
|
|
|
|
+
|
|
|
|
+ __m256 xs = _mm256_add_ps(_mm256_mul_ps(dpp, pixc), viewx);
|
|
|
|
+ __m256 xs2 = _mm256_add_ps(_mm256_mul_ps(dpp, pixc2), viewx);
|
|
|
|
+ __m256 xs3 = _mm256_add_ps(_mm256_mul_ps(dpp, pixc3), viewx);
|
|
|
|
+
|
|
|
|
+ __m256 counter = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+ __m256 adder = { 1, 1, 1, 1, 1, 1, 1, 1 };
|
|
|
|
+ __m256 resultsa = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+ __m256 resultsb = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+
|
|
|
|
+ __m256 counter2 = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+ __m256 adder2 = { 1, 1, 1, 1, 1, 1, 1, 1 };
|
|
|
|
+ __m256 resultsa2 = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+ __m256 resultsb2 = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+
|
|
|
|
+ __m256 counter3 = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+ __m256 adder3 = { 1, 1, 1, 1, 1, 1, 1, 1 };
|
|
|
|
+ __m256 resultsa3 = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+ __m256 resultsb3 = { 0, 0, 0, 0, 0, 0, 0, 0 };
|
|
|
|
+
|
|
|
|
+ __m256 threshold = { 16.0f, 16.0f, 16.0f, 16.0f, 16.0f, 16.0f, 16.0f, 16.0f };
|
|
|
|
+ __m256 two = { 2, 2, 2, 2, 2, 2, 2, 2 };
|
|
|
|
+
|
|
|
|
+ __m256 a = xs;
|
|
|
|
+ __m256 a2 = xs2;
|
|
|
|
+ __m256 a3 = xs3;
|
|
|
|
+ __m256 b = ys;
|
|
|
|
+ __m256 b2 = ys;
|
|
|
|
+ __m256 b3 = ys;
|
|
|
|
+
|
|
|
|
+ __m256 cx = info.julia ? juliaX : xs;
|
|
|
|
+ __m256 cx2 = info.julia ? juliaX : xs2;
|
|
|
|
+ __m256 cx3 = info.julia ? juliaX : xs3;
|
|
|
|
+ __m256 cy = info.julia ? juliaY : ys;
|
|
|
|
+
|
|
|
|
+ if (info.smooth) {
|
|
|
|
+ __m256 cmp = _mm256_cmp_ps(threshold, threshold, _CMP_LE_OQ);
|
|
|
|
+ __m256 cmp2 = _mm256_cmp_ps(threshold, threshold, _CMP_LE_OQ);
|
|
|
|
+ __m256 cmp3 = _mm256_cmp_ps(threshold, threshold, _CMP_LE_OQ);
|
|
|
|
+ for (int k = 0; k < info.maxIter; k++) {
|
|
|
|
+ __m256 bb = _mm256_mul_ps(b, b);
|
|
|
|
+ __m256 bb2 = _mm256_mul_ps(b2, b2);
|
|
|
|
+ __m256 bb3 = _mm256_mul_ps(b3, b3);
|
|
|
|
+ __m256 ab = _mm256_mul_ps(a, b);
|
|
|
|
+ __m256 ab2 = _mm256_mul_ps(a2, b2);
|
|
|
|
+ __m256 ab3 = _mm256_mul_ps(a3, b3);
|
|
|
|
+ __m256 olda = a;
|
|
|
|
+ __m256 olda2 = a2;
|
|
|
|
+ __m256 olda3 = a3;
|
|
|
|
+ a = _mm256_add_ps(_mm256_fmsub_ps(a, a, bb), cx);
|
|
|
|
+ a2 = _mm256_add_ps(_mm256_fmsub_ps(a2, a2, bb2), cx2);
|
|
|
|
+ a3 = _mm256_add_ps(_mm256_fmsub_ps(a3, a3, bb3), cx3);
|
|
|
|
+ b = _mm256_fmadd_ps(two, ab, cy);
|
|
|
|
+ b2 = _mm256_fmadd_ps(two, ab2, cy);
|
|
|
|
+ b3 = _mm256_fmadd_ps(two, ab3, cy);
|
|
|
|
+ /*resultsa = _mm256_or_ps(_mm256_andnot_ps(cmp, resultsa), _mm256_and_ps(cmp, a));
|
|
|
|
+ resultsb = _mm256_or_ps(_mm256_andnot_ps(cmp, resultsb), _mm256_and_ps(cmp, b));
|
|
|
|
+ resultsa2 = _mm256_or_ps(_mm256_andnot_ps(cmp2, resultsa2), _mm256_and_ps(cmp2, a2));
|
|
|
|
+ resultsb2 = _mm256_or_ps(_mm256_andnot_ps(cmp2, resultsb2), _mm256_and_ps(cmp2, b2));
|
|
|
|
+ resultsa3 = _mm256_or_ps(_mm256_andnot_ps(cmp3, resultsa3), _mm256_and_ps(cmp3, a3));
|
|
|
|
+ resultsb3 = _mm256_or_ps(_mm256_andnot_ps(cmp3, resultsb3), _mm256_and_ps(cmp3, b3));*/
|
|
|
|
+ resultsa = _mm256_blendv_ps(resultsa, a, cmp);
|
|
|
|
+ resultsb = _mm256_blendv_ps(resultsb, b, cmp);
|
|
|
|
+ resultsa2 = _mm256_blendv_ps(resultsa2, a2, cmp2);
|
|
|
|
+ resultsb2 = _mm256_blendv_ps(resultsb2, b2, cmp2);
|
|
|
|
+ resultsa3 = _mm256_blendv_ps(resultsa3, a3, cmp3);
|
|
|
|
+ resultsb3 = _mm256_blendv_ps(resultsb3, b3, cmp3);
|
|
|
|
+ cmp = _mm256_cmp_ps(_mm256_fmadd_ps(olda, olda, bb), threshold, _CMP_LE_OQ);
|
|
|
|
+ cmp2 = _mm256_cmp_ps(_mm256_fmadd_ps(olda2, olda2, bb2), threshold, _CMP_LE_OQ);
|
|
|
|
+ cmp3 = _mm256_cmp_ps(_mm256_fmadd_ps(olda3, olda3, bb3), threshold, _CMP_LE_OQ);
|
|
|
|
+ adder = _mm256_and_ps(adder, cmp);
|
|
|
|
+ counter = _mm256_add_ps(counter, adder);
|
|
|
|
+ adder2 = _mm256_and_ps(adder2, cmp2);
|
|
|
|
+ counter2 = _mm256_add_ps(counter2, adder2);
|
|
|
|
+ adder3 = _mm256_and_ps(adder3, cmp3);
|
|
|
|
+ counter3 = _mm256_add_ps(counter3, adder3);
|
|
|
|
+ if ((k & 0x7) == 0 && _mm256_testz_ps(cmp, cmp) != 0 && _mm256_testz_ps(cmp2, cmp2) != 0 && _mm256_testz_ps(cmp3, cmp3) != 0) {
|
|
|
|
+ break;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ else {
|
|
|
|
+ for (int k = 0; k < info.maxIter; k++) {
|
|
|
|
+ __m256 bb = _mm256_mul_ps(b, b);
|
|
|
|
+ __m256 bb2 = _mm256_mul_ps(b2, b2);
|
|
|
|
+ __m256 bb3 = _mm256_mul_ps(b3, b3);
|
|
|
|
+ __m256 ab = _mm256_mul_ps(a, b);
|
|
|
|
+ __m256 ab2 = _mm256_mul_ps(a2, b2);
|
|
|
|
+ __m256 ab3 = _mm256_mul_ps(a3, b3);
|
|
|
|
+ __m256 cmp = _mm256_cmp_ps(_mm256_fmadd_ps(a, a, bb), threshold, _CMP_LE_OQ);
|
|
|
|
+ __m256 cmp2 = _mm256_cmp_ps(_mm256_fmadd_ps(a2, a2, bb2), threshold, _CMP_LE_OQ);
|
|
|
|
+ __m256 cmp3 = _mm256_cmp_ps(_mm256_fmadd_ps(a3, a3, bb3), threshold, _CMP_LE_OQ);
|
|
|
|
+ a = _mm256_add_ps(_mm256_fmsub_ps(a, a, bb), cx);
|
|
|
|
+ a2 = _mm256_add_ps(_mm256_fmsub_ps(a2, a2, bb2), cx2);
|
|
|
|
+ a3 = _mm256_add_ps(_mm256_fmsub_ps(a3, a3, bb3), cx3);
|
|
|
|
+ b = _mm256_fmadd_ps(two, ab, cy);
|
|
|
|
+ b2 = _mm256_fmadd_ps(two, ab2, cy);
|
|
|
|
+ b3 = _mm256_fmadd_ps(two, ab3, cy);
|
|
|
|
+ adder = _mm256_and_ps(adder, cmp);
|
|
|
|
+ counter = _mm256_add_ps(counter, adder);
|
|
|
|
+ adder2 = _mm256_and_ps(adder2, cmp2);
|
|
|
|
+ counter2 = _mm256_add_ps(counter2, adder2);
|
|
|
|
+ adder3 = _mm256_and_ps(adder3, cmp3);
|
|
|
|
+ counter3 = _mm256_add_ps(counter3, adder3);
|
|
|
|
+ if ((k & 0x7) == 0 && _mm256_testz_ps(cmp, cmp) != 0 && _mm256_testz_ps(cmp2, cmp2) != 0 && _mm256_testz_ps(cmp3, cmp3) != 0) {
|
|
|
|
+ break;
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+
|
|
|
|
+ float resData[96];
|
|
|
|
+ float* ftRes = resData;
|
|
|
|
+ float* resa = ftRes + 24;
|
|
|
|
+ float* resb = resa + 24;
|
|
|
|
+
|
|
|
|
+ _mm256_storeu_ps(ftRes, counter);
|
|
|
|
+ _mm256_storeu_ps(ftRes + 8, counter2);
|
|
|
|
+ _mm256_storeu_ps(ftRes + 16, counter3);
|
|
|
|
+ _mm256_storeu_ps(resa, resultsa);
|
|
|
|
+ _mm256_storeu_ps(resa + 8, resultsa2);
|
|
|
|
+ _mm256_storeu_ps(resa + 16, resultsa3);
|
|
|
|
+ _mm256_storeu_ps(resb, resultsb);
|
|
|
|
+ _mm256_storeu_ps(resb + 8, resultsb2);
|
|
|
|
+ _mm256_storeu_ps(resb + 16, resultsb3);
|
|
|
|
+ for (int k = 0; k < 24 && i + k < info.bWidth; k++) {
|
|
|
|
+ if (info.smooth) {
|
|
|
|
+ data[i + k + j * info.bWidth] = ftRes[k] < 0 ? info.maxIter :
|
|
|
|
+ ftRes[k] >= info.maxIter ? info.maxIter :
|
|
|
|
+ ((float)ftRes[k]) + 1 - ::log2(::log(resa[k] * resa[k] + resb[k] * resb[k]) / 2);
|
|
|
|
+ }
|
|
|
|
+ else {
|
|
|
|
+ data[i + k + j * info.bWidth] = ftRes[k] < 0 ? info.maxIter : ftRes[k];
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+}
|
|
|
|
+
|