aya

libmandel/include/ClGenerators.h

@@ -17,6 +17,7 @@ namespace mnd
     class ClGeneratorFloat;
     class ClGeneratorDouble;
     class ClGeneratorDoubleDouble;
+    class ClGeneratorQuadDouble;
     class ClGenerator128;
 }
 
@@ -73,6 +74,20 @@ protected:
     virtual std::string getKernelCode(bool smooth) const;
 };
 
+
+class mnd::ClGeneratorQuadDouble : public ClGenerator
+{
+    bool smooth;
+public:
+    ClGeneratorQuadDouble(cl::Device device);
+    virtual ~ClGeneratorQuadDouble(void) = default;
+
+    virtual void generate(const MandelInfo& info, float* data);
+protected:
+    virtual std::string getKernelCode(bool smooth) const;
+};
+
+
 class mnd::ClGenerator128 : public ClGenerator
 {
 public:

libmandel/src/ClGenerators.cpp

@@ -13,6 +13,7 @@ using mnd::ClGenerator;
 using mnd::ClGeneratorFloat;
 using mnd::ClGeneratorDouble;
 using mnd::ClGeneratorDoubleDouble;
+using mnd::ClGeneratorQuadDouble;
 using mnd::ClGenerator128;
 
 Platform getPlatform() {
@@ -20,7 +21,7 @@ Platform getPlatform() {
     std::vector<Platform> all_platforms;
     Platform::get(&all_platforms);
 
-    if (all_platforms.size()==0) {
+    if (all_platforms.size() == 0) {
         std::cout << "No platforms found. Check OpenCL installation!\n";
         exit(1);
     }
@@ -300,6 +301,62 @@ std::string ClGeneratorDoubleDouble::getKernelCode(bool smooth) const
 }
 
 
+ClGeneratorQuadDouble::ClGeneratorQuadDouble(cl::Device device) :
+    ClGenerator{ device }
+{
+    context = Context{ device };
+    Program::Sources sources;
+
+    std::string kcode = this->getKernelCode(false);
+
+    sources.push_back({ kcode.c_str(), kcode.length() });
+
+    program = Program{ context, sources };
+    if (program.build({ device }) != CL_SUCCESS) {
+        throw std::string(program.getBuildInfo<CL_PROGRAM_BUILD_LOG>(device));
+    }
+
+    queue = CommandQueue(context, device);
+}
+
+
+void ClGeneratorQuadDouble::generate(const mnd::MandelInfo& info, float* data)
+{
+    ::size_t bufferSize = info.bWidth * info.bHeight * sizeof(float);
+
+    Buffer buffer_A(context, CL_MEM_WRITE_ONLY, bufferSize);
+
+    mnd::DoubleDouble x = mnd::convert<mnd::DoubleDouble>(info.view.x);
+    mnd::DoubleDouble y = mnd::convert<mnd::DoubleDouble>(info.view.y);
+
+    mnd::DoubleDouble psx = mnd::convert<mnd::DoubleDouble>(info.view.width / info.bWidth);
+    mnd::DoubleDouble psy = mnd::convert<mnd::DoubleDouble>(info.view.height / info.bHeight);
+
+    Kernel iterate = Kernel(program, "iterate");
+    iterate.setArg(0, buffer_A);
+    iterate.setArg(1, int(info.bWidth));
+    iterate.setArg(2, x.x[0]);
+    iterate.setArg(3, x.x[1]);
+    iterate.setArg(4, y.x[0]);
+    iterate.setArg(5, y.x[1]);
+    iterate.setArg(6, psx.x[0]);
+    iterate.setArg(7, psx.x[1]);
+    iterate.setArg(8, psy.x[0]);
+    iterate.setArg(9, psy.x[1]);
+    iterate.setArg(10, int(info.maxIter));
+    iterate.setArg(11, int(info.smooth ? 1 : 0));
+
+    cl_int result = queue.enqueueNDRangeKernel(iterate, 0, NDRange(info.bWidth * info.bHeight));
+    queue.enqueueReadBuffer(buffer_A, CL_TRUE, 0, bufferSize, data);
+}
+
+
+std::string ClGeneratorQuadDouble::getKernelCode(bool smooth) const
+{
+    return (char*) doubledouble_cl;
+}
+
+
 ClGenerator128::ClGenerator128(cl::Device device) :
     ClGenerator{ device }
 {

libmandel/src/CpuGeneratorsAVXFMA.cpp

@@ -0,0 +1,267 @@
+#include "CpuGenerators.h"
+
+#include <immintrin.h>
+#include <omp.h>
+#include <cmath>
+
+#include <utility>
+#include <memory>
+
+using mnd::CpuGenerator;
+
+namespace mnd
+{
+    template class CpuGenerator<DoubleDouble, mnd::X86_AVX_FMA, false>;
+    template class CpuGenerator<DoubleDouble, mnd::X86_AVX_FMA, true>;
+}
+
+
+struct VecPair
+{
+    __m256d a;
+    __m256d b;
+};
+
+
+static inline VecPair quickTwoSum(__m256d a, __m256d b)
+{
+    __m256d s = _mm256_add_pd(a, b);
+    __m256d e = _mm256_sub_pd(b, _mm256_sub_pd(s, a));
+    return { s, e };
+}
+
+static inline VecPair quickTwoDiff(__m256d a, __m256d b)
+{
+    __m256d s = _mm256_sub_pd(a, b);
+    __m256d e = _mm256_sub_pd(_mm256_sub_pd(a, s), b);
+    return { s, e };
+}
+
+static inline VecPair twoSum(__m256d a, __m256d b)
+{
+    __m256d s = _mm256_add_pd(a, b);
+    __m256d bb = _mm256_sub_pd(s, a);
+    __m256d e = _mm256_add_pd(_mm256_sub_pd(a, _mm256_sub_pd(s, bb)), _mm256_sub_pd(b, bb));
+    return { s, e };
+}
+
+static inline VecPair twoDiff(__m256d a, __m256d b)
+{
+    __m256d s = _mm256_sub_pd(a, b);
+    __m256d bb = _mm256_sub_pd(s, a);
+    __m256d e = _mm256_sub_pd(_mm256_sub_pd(a, _mm256_sub_pd(s, bb)), _mm256_add_pd(b, bb));
+    return { s, e };
+}
+
+/*
+static inline VecPair split(__m256d a)
+{
+    static const __m256d SPLIT_THRESH = { 6.69692879491417e+299, 6.69692879491417e+299, 6.69692879491417e+299, 6.69692879491417e+299 };
+    static const __m256d MINUS_SPLIT_THRESH = { -6.69692879491417e+299, -6.69692879491417e+299, -6.69692879491417e+299, -6.69692879491417e+299 };
+    static const __m256d SPLITTER = { 134217729.0, 134217729.0, 134217729.0, 134217729.0};
+    __m256d temp;
+    __m256i cmp1 = _mm256_castpd_si256(_mm256_cmp_pd(a, SPLIT_THRESH, _CMP_GT_OQ));
+    __m256i cmp2 = _mm256_castpd_si256(_mm256_cmp_pd(a, MINUS_SPLIT_THRESH, _CMP_LT_OQ));
+    __m256i cmp = _mm256_or_si256
+}*/
+
+static inline VecPair twoProd(__m256d a, __m256d b)
+{
+//#ifdef CPUID_FMA
+    __m256d p = _mm256_mul_pd(a, b);
+    __m256d e = _mm256_fmsub_pd(a, b, p);
+    return { p, e };
+//#else
+/*    double a_hi, a_lo, b_hi, b_lo;
+    __m256d p = _mm256_mul_ps(a, b);
+    split(a, a_hi, a_lo);
+    split(b, b_hi, b_lo);
+    err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
+    return p;*/
+//#endif
+}
+
+struct AvxDoubleDouble
+{
+    __m256d x[2];
+
+    inline AvxDoubleDouble(__m256d a, __m256d b) :
+        x{ a, b }
+    {}
+
+
+    inline AvxDoubleDouble operator + (const AvxDoubleDouble& sm) const
+    {
+        auto[s, e] = twoSum(x[0], sm.x[0]);
+        e = _mm256_add_pd(e, _mm256_add_pd(x[1], sm.x[1]));
+        auto[r1, r2] = quickTwoSum(s, e);
+        return AvxDoubleDouble{ r1, r2 };
+    }
+
+    inline AvxDoubleDouble operator - (const AvxDoubleDouble& sm) const
+    {
+        auto[s, e] = twoDiff(x[0], sm.x[0]);
+        e = _mm256_add_pd(e, x[1]);
+        e = _mm256_sub_pd(e, sm.x[1]);
+        auto[r1, r2] = quickTwoSum(s, e);
+        return AvxDoubleDouble{ r1, r2 };
+    }
+
+    inline AvxDoubleDouble operator * (const AvxDoubleDouble& sm) const
+    {
+        auto[p1, p2] = twoProd(this->x[0], sm.x[0]);
+        p2 = _mm256_add_pd(p2,
+            _mm256_add_pd(_mm256_mul_pd(sm.x[1], x[0]), _mm256_mul_pd(sm.x[0], x[1])) );
+        auto[r1, r2] = quickTwoSum(p1, p2);
+        return AvxDoubleDouble{ r1, r2 };
+    }
+};
+
+template<bool parallel>
+void CpuGenerator<mnd::DoubleDouble, mnd::X86_AVX_FMA, parallel>::generate(const mnd::MandelInfo& info, float* data)
+{
+    const MandelViewport& view = info.view;
+
+    using T = DoubleDouble;
+
+    T viewx = mnd::convert<T>(view.x);
+    T viewy = mnd::convert<T>(view.y);
+    T wpp = mnd::convert<T>(view.width / info.bWidth);
+    T hpp = mnd::convert<T>(view.height / info.bHeight);
+
+//    if constexpr(parallel)
+//        omp_set_num_threads(2 * omp_get_num_procs());
+//#pragma omp parallel for schedule(static, 1) if (parallel)
+
+    for (long j = 0; j < info.bHeight; j++) {
+        T y = viewy + T(double(j)) * hpp;
+        __m256d y0s = { y.x[0], y.x[0], y.x[0], y.x[0] };
+        __m256d y1s = { y.x[1], y.x[1], y.x[1], y.x[1] };
+        AvxDoubleDouble ys{ y0s, y1s };
+        long i = 0;
+        for (i; i < info.bWidth; i += 4) {
+            T x1 = viewx + T(double(i)) * wpp;
+            T x2 = x1 + wpp;
+            T x3 = x2 + wpp;
+            T x4 = x3 + wpp;
+
+            __m256d x0s = {
+                x1.x[0], x2.x[0], x3.x[0], x4.x[0],
+            };
+
+            __m256d x1s = {
+                x1.x[1], x2.x[1], x3.x[1], x4.x[1],
+            };
+
+            AvxDoubleDouble xs{ x0s, x1s };
+
+            int itRes[4] = { 0, 0, 0, 0 };
+
+            __m256d threshold = { 16.0, 16.0, 16.0, 16.0 };
+            __m256d counter = { 0, 0, 0, 0 };
+            __m256d adder = { 1, 1, 1, 1 };
+
+            AvxDoubleDouble a = xs;
+            AvxDoubleDouble b = ys;
+
+            for (int k = 0; k < info.maxIter; k++) {
+                AvxDoubleDouble aa = a * a;
+                AvxDoubleDouble bb = b * b;
+                AvxDoubleDouble abab = a * b; abab = abab + abab;
+                a = aa - bb + xs;
+                b = abab + ys;
+                __m256i cmp = _mm256_castpd_si256(_mm256_cmp_pd(_mm256_add_pd(aa.x[0], bb.x[0]), threshold, _CMP_LE_OQ));
+                /*if (info.smooth) {
+                    resultsa = _mm256_or_pd(_mm256_andnot_ps(cmp, resultsa), _mm256_and_ps(cmp, a));
+                    resultsb = _mm256_or_ps(_mm256_andnot_ps(cmp, resultsb), _mm256_and_ps(cmp, b));
+                }*/
+                adder = _mm256_and_pd(adder, _mm256_castsi256_pd(cmp));
+                counter = _mm256_add_pd(counter, adder);
+                if (_mm256_testz_si256(cmp, cmp) != 0) {
+                    break;
+                }
+            }
+
+            auto alignVec = [](double* data) -> double* {
+                void* aligned = data;
+                ::size_t length = 64;
+                std::align(32, 4 * sizeof(double), aligned, length);
+                return static_cast<double*>(aligned);
+            };
+
+            double resData[8];
+            double* ftRes = alignVec(resData);
+            _mm256_store_pd(ftRes, counter);
+            for (int k = 0; k < 4 && i + k < info.bWidth; k++)
+                data[i + k + j * info.bWidth] = ftRes[k] > 0 ? float(ftRes[k]) : info.maxIter;
+        }
+    }
+    return;
+
+
+
+    for (long j = 0; j < info.bHeight; j++) {
+        T y = viewy + T(double(j)) * hpp;
+        __m256d y0s = { y.x[0], y.x[0], y.x[0], y.x[0] };
+        __m256d y1s = { y.x[1], y.x[1], y.x[1], y.x[1] };
+        AvxDoubleDouble ys{ y0s, y1s };
+        long i = 0;
+        for (i; i < info.bWidth; i += 4) {
+            T x1 = viewx + T(double(i)) * wpp;
+            T x2 = viewx + T(double(i + 1)) * wpp;
+            T x3 = viewx + T(double(i + 2)) * wpp;
+            T x4 = viewx + T(double(i + 3)) * wpp;
+
+            __m256d x0s = {
+                x1.x[0], x2.x[0], x3.x[0], x4.x[0],
+            };
+
+            __m256d x1s = {
+                x1.x[1], x2.x[1], x3.x[1], x4.x[1],
+            };
+
+            AvxDoubleDouble xs{ x0s, x1s };
+
+            int itRes[4] = { 0, 0, 0, 0 };
+
+            __m256d threshold = { 16.0, 16.0, 16.0, 16.0 };
+            __m256d counter = { 0, 0, 0, 0 };
+            __m256d adder = { 1, 1, 1, 1 };
+
+            AvxDoubleDouble a = xs;
+            AvxDoubleDouble b = ys;
+
+            for (int k = 0; k < info.maxIter; k++) {
+                AvxDoubleDouble aa = a * a;
+                AvxDoubleDouble bb = b * b;
+                AvxDoubleDouble abab = a * b; abab = abab + abab;
+                a = aa - bb + xs;
+                b = abab + ys;
+                __m256i cmp = _mm256_castpd_si256(_mm256_cmp_pd(_mm256_add_pd(aa.x[0], bb.x[0]), threshold, _CMP_LE_OQ));
+                /*if (info.smooth) {
+                    resultsa = _mm256_or_pd(_mm256_andnot_ps(cmp, resultsa), _mm256_and_ps(cmp, a));
+                    resultsb = _mm256_or_ps(_mm256_andnot_ps(cmp, resultsb), _mm256_and_ps(cmp, b));
+                }*/
+                adder = _mm256_and_pd(adder, _mm256_castsi256_pd(cmp));
+                counter = _mm256_add_pd(counter, adder);
+                if ((k & 7) == 0 && _mm256_testz_si256(cmp, cmp) != 0) {
+                    break;
+                }
+            }
+
+            auto alignVec = [](double* data) -> double* {
+                void* aligned = data;
+                ::size_t length = 64;
+                std::align(32, 4 * sizeof(double), aligned, length);
+                return static_cast<double*>(aligned);
+            };
+
+            double resData[8];
+            double* ftRes = alignVec(resData);
+            _mm256_store_pd(ftRes, counter);
+            for (int k = 0; k < 4 && i + k < info.bWidth; k++)
+                data[i + k + j * info.bWidth] = ftRes[k] > 0 ? float(ftRes[k]) : info.maxIter;
+        }
+    }
+}
+

libmandel/src/Mandel.cpp

@@ -157,12 +157,14 @@ std::vector<MandelDevice> MandelContext::createDevices(void)
         std::string name = platform.getInfo<CL_PLATFORM_NAME>();
         std::string profile = platform.getInfo<CL_PLATFORM_PROFILE>();
 
+        printf("using opencl platform: %s\n", name.c_str());
+
         //std::string ext = platform.getInfo<CL_PLATFORM_EXTENSIONS>();
         //printf("Platform extensions: %s\n", ext.c_str());
         //printf("Platform: %s, %s\n", name.c_str(), profile.c_str());
 
         std::vector<cl::Device> devices;
-        platform.getDevices(CL_DEVICE_TYPE_GPU, &devices);
+        platform.getDevices(CL_DEVICE_TYPE_GPU | CL_DEVICE_TYPE_CPU, &devices);
         for (auto& device : devices) {
             //printf("Device: %s\n", device.getInfo<CL_DEVICE_NAME>().c_str());
             //printf("preferred float width: %d\n", device.getInfo<CL_DEVICE_PREFERRED_VECTOR_WIDTH_FLOAT>());
@@ -178,6 +180,7 @@ std::vector<MandelDevice> MandelContext::createDevices(void)
 
             md.name = device.getInfo<CL_DEVICE_NAME>();
             md.vendor = device.getInfo<CL_DEVICE_VENDOR>();
+            printf("    using opencl device: %s\n", md.name.c_str());
             try {
                 md.generators.insert({ GeneratorType::FLOAT, std::make_unique<ClGeneratorFloat>(device) });
             }
@@ -189,6 +192,7 @@ std::vector<MandelDevice> MandelContext::createDevices(void)
                 try {
                     md.generators.insert({ GeneratorType::DOUBLE, std::make_unique<ClGeneratorDouble>(device) });
                     md.generators.insert({ GeneratorType::DOUBLE_DOUBLE, std::make_unique<ClGeneratorDoubleDouble>(device) });
+                    md.generators.insert({ GeneratorType::QUAD_DOUBLE, std::make_unique<ClGeneratorQuadDouble>(device) });
                 }
                 catch (const std::string& err) {
                     printf("err: %s", err.c_str());

libmandel/src/quaddouble.cl

@@ -14,20 +14,55 @@ inline double2 quickTwoSum(double a, double b) {
 }
 
 inline double2 twoProd(double a, double b) {
-//#ifdef QD_FMS
     double p = a * b;
     double e = fma(a, b, -p);
     return (double2)(p, e);
-//#else
-//  double a_hi, a_lo, b_hi, b_lo;
-//  double p = a * b;
-//  split(a, a_hi, a_lo);
-//  split(b, b_hi, b_lo);
-//  err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
-//  return p;
-//#endif
 }
 
+
+inline void threeSum(double* a, double* b, double* c) {
+    double2 t = twoSum(*a, *b);
+    double2 at3 = twoSum(*c, t.s0);
+    double2 bc = twoSum(t.s1, at3.s1);
+    *a = at3.s0;
+    *b = bc.s0;
+    *c = bc.s1;
+}
+
+
+inline void threeSum2(double* a, double* b, double* c) {
+    double2 t = twoSum(*a, *b);
+    double2 at3 = twoSum(*c, t.s0);
+    *a = at3.s0;
+    *b = t.s1 + at3.s1;
+}
+
+
+inline double4 add(double4 a, double4 b) {
+    double2 su0 = twoSum(a.s0, b.s0);
+    double2 su1 = twoSum(a.s1, b.s1);
+    double2 su2 = twoSum(a.s2, b.s2);
+    double2 su3 = twoSum(a.s3, b.s3);
+
+    double2 s1t0 = twoSum(su1.s0, su0.s1);
+    threeSum(&su2.s0, &su0.s1, &su1.s1);
+    threeSum2(&su3.s0, &su0.s1, &su2.s1);
+    su0.s1 = su0.s1 + su1.s1 + su3.s1;
+
+    renorm(su0.s0, su1.s0, su2.s0, su3.s0, su0.s1);
+    return (double4)(su0.s0, su1.s0, su2.s0, su3.s0);
+}
+
+
+
+
+
+
+
+
+
+
+
 inline double2 mul(double2 a, double2 b) {
     double2 p = twoProd(a.s0, b.s0);
     p.s1 += (a.s0 * b.s1 + a.s1 * b.s0);