Task03 Ilya Bolkisev ITMO

src/cl/mandelbrot.cl

@@ -4,10 +4,42 @@
 
 #line 6
 
-__kernel void mandelbrot(...)
+__kernel void mandelbrot(__global float* results,
+                         int width, int height,
+                         float fromX, float fromY,
+                         float sizeX, float sizeY,
+                         uint iters,
+                         int smoothing)
 {
-    // TODO если хочется избавиться от зернистости и дрожания при интерактивном погружении, добавьте anti-aliasing:
-    // грубо говоря, при anti-aliasing уровня N вам нужно рассчитать не одно значение в центре пикселя, а N*N значений
-    // в узлах регулярной решетки внутри пикселя, а затем посчитав среднее значение результатов - взять его за результат для всего пикселя
-    // это увеличит число операций в N*N раз, поэтому при рассчетах гигаплопс антиальясинг должен быть выключен
+    int gx = get_global_id(0);
+    int gy = get_global_id(1);
+
+    if (gx >= width || gy >= height)
+        return;
+
+    const float threshold = 256.0f;
+    const float threshold2 = threshold * threshold;
+
+    float x0 = fromX + (gx + 0.5f) * (sizeX / width);
+    float y0 = fromY + (gy + 0.5f) * (sizeY / height);
+
+    float x = x0;
+    float y = y0;
+    int iter = 0;
+
+    for (; iter < iters; ++iter) {
+        float xPrev = x;
+        x = x * x - y * y + x0;
+        y = 2.0f * xPrev * y + y0;
+        if ((x * x + y * y) > threshold2)
+            break;
+
+    }
+
+    float result = (float)iter;
+    if (smoothing == 1 && iter != iters)
+        result -= log(log(sqrt(x * x + y * y)) / log(threshold)) / log(2.0f);
+
+    result /= iters;
+    results[gy * width + gx] = result;
 }

src/cl/sum.cl

@@ -1 +1,94 @@
-// TODO
+#ifdef __CLION_IDE__
+#include <libgpu/opencl/cl/clion_defines.cl>
+#endif
+#line 6
+
+__kernel void sum_atomic(__global const unsigned int *arr, __global unsigned int *sum, unsigned int n) {
+    const unsigned int gid = get_global_id(0);
+    if (gid >= n)
+        return;
+    atomic_add(sum, arr[gid]);
+}
+
+
+__kernel void sum_loop(__global const unsigned int* arr, __global unsigned int* sum, unsigned int n)
+{
+    unsigned int gid = get_global_id(0);
+    unsigned int gsize = get_global_size(0);
+
+    unsigned int partial = 0;
+    for (unsigned int i = gid; i < n; i += gsize)
+        partial += arr[i];
+
+    atomic_add(sum, partial);
+}
+
+
+__kernel void sum_loop_coalesced(__global const unsigned int* arr, __global unsigned int* sum, unsigned int n)
+{
+    unsigned int gid = get_global_id(0);
+    unsigned int gsize = get_global_size(0);
+
+    unsigned int partial = 0;
+    for (unsigned int i = gid; i < n; i += gsize)
+        partial += arr[i];
+
+    atomic_add(sum, partial);
+}
+
+
+__kernel void sum_local(__global const unsigned int* arr, __global unsigned int* sum, unsigned int n)
+{
+    __local unsigned int shared_data[256];
+
+    unsigned int lid = get_local_id(0);
+    unsigned int lsize = get_local_size(0);
+
+    unsigned int gid = get_global_id(0);
+    unsigned int gsize = get_global_size(0);
+
+    unsigned int partial = 0;
+    for (unsigned int i = gid; i < n; i += gsize)
+        partial += arr[i];
+
+    shared_data[lid] = partial;
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    for (unsigned int offset = lsize / 2; offset > 0; offset >>= 1) {
+        if (lid < offset)
+            shared_data[lid] += shared_data[lid + offset];
+        barrier(CLK_LOCAL_MEM_FENCE);
+    }
+
+    if (lid == 0) {
+        atomic_add(sum, shared_data[0]);
+    }
+}
+
+
+__kernel void sum_tree(__global const unsigned int* arr, __global unsigned int* sum, unsigned int n)
+{
+    __local unsigned int local_data[256];
+
+    unsigned int gid = get_global_id(0);
+    unsigned int lid = get_local_id(0);
+    unsigned int group_size = get_local_size(0);
+    unsigned int global_size = get_global_size(0);
+
+    unsigned int partial = 0;
+    for (unsigned int i = gid; i < n; i += global_size)
+        partial += arr[i];
+
+    local_data[lid] = partial;
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    for (unsigned int offset = group_size >> 1; offset > 0; offset >>= 1) {
+        if (lid < offset)
+            local_data[lid] += local_data[lid + offset];
+        barrier(CLK_LOCAL_MEM_FENCE);
+    }
+
+    if (lid == 0) {
+        atomic_add(sum, local_data[0]);
+    }
+}

src/main_mandelbrot.cpp

@@ -102,44 +102,74 @@ int main(int argc, char **argv)
         std::cout << "    Real iterations fraction: " << 100.0 * realIterationsFraction / (width * height) << "%" << std::endl;
 
         renderToColor(cpu_results.ptr(), image.ptr(), width, height);
-        image.savePNG("mandelbrot_cpu.png");
+        image.savePNG("../mandelbrot_cpu.png");
     }
 
 
-//    // Раскомментируйте это:
-//
-//    gpu::Context context;
-//    context.init(device.device_id_opencl);
-//    context.activate();
-//    {
-//        ocl::Kernel kernel(mandelbrot_kernel, mandelbrot_kernel_length, "mandelbrot");
-//        // Если у вас есть интеловский драйвер для запуска на процессоре - попробуйте запустить на нем и взгляните на лог,
-//        // передав printLog=true - скорее всего, в логе будет строчка вроде
-//        // Kernel <mandelbrot> was successfully vectorized (8)
-//        // это означает, что драйвер смог векторизовать вычисления с помощью интринсик, и если множитель векторизации 8, то
-//        // это означает, что одно ядро процессит сразу 8 workItems, а т.к. все вычисления в float, то
-//        // это означает, что используются 8 x float регистры (т.е. 256-битные, т.е. AVX)
-//        // обратите внимание, что и произвдительность относительно референсной ЦПУ реализации выросла почти в восемь раз
-//        bool printLog = false;
-//        kernel.compile(printLog);
-//        // TODO близко к ЦПУ-версии, включая рассчет таймингов, гигафлопс, Real iterations fraction и сохранение в файл
-//        // результат должен оказаться в gpu_results
-//    }
-//
-//    {
-//        double errorAvg = 0.0;
-//        for (int j = 0; j < height; ++j) {
-//            for (int i = 0; i < width; ++i) {
-//                errorAvg += fabs(gpu_results.ptr()[j * width + i] - cpu_results.ptr()[j * width + i]);
-//            }
-//        }
-//        errorAvg /= width * height;
-//        std::cout << "GPU vs CPU average results difference: " << 100.0 * errorAvg << "%" << std::endl;
-//
-//        if (errorAvg > 0.03) {
-//            throw std::runtime_error("Too high difference between CPU and GPU results!");
-//        }
-//    }
+    gpu::Context context;
+    context.init(device.device_id_opencl);
+    context.activate();
+    {
+        ocl::Kernel kernel(mandelbrot_kernel, mandelbrot_kernel_length, "mandelbrot");
+        // Если у вас есть интеловский драйвер для запуска на процессоре - попробуйте запустить на нем и взгляните на лог,
+        // передав printLog=true - скорее всего, в логе будет строчка вроде
+        // Kernel <mandelbrot> was successfully vectorized (8)
+        // это означает, что драйвер смог векторизовать вычисления с помощью интринсик, и если множитель векторизации 8, то
+        // это означает, что одно ядро процессит сразу 8 workItems, а т.к. все вычисления в float, то
+        // это означает, что используются 8 x float регистры (т.е. 256-битные, т.е. AVX)
+        // обратите внимание, что и произвдительность относительно референсной ЦПУ реализации выросла почти в восемь раз
+        bool printLog = false;
+        kernel.compile(printLog);
+
+        gpu::gpu_mem_32f results_gpu;
+        results_gpu.resizeN(width * height);
+
+        timer t;
+        for (int i = 0; i < benchmarkingIters; ++i) {
+            kernel.exec(
+                    gpu::WorkSize(16, 16, width, height),
+                    results_gpu,(int)width, (int)height, centralX - sizeX / 2.0f,
+                    centralY - sizeY / 2.0f, sizeX, sizeY, iterationsLimit, 0);
+            t.nextLap();
+        }
+
+        std::cout << "GPU: " << t.lapAvg() << "+-" << t.lapStd() << " s" << std::endl;
+
+        size_t flopsInLoop = 10;
+        size_t maxApproximateFlops = (size_t)width * height * iterationsLimit * flopsInLoop;
+        size_t gflops = 1000*1000*1000;
+        std::cout << "GPU: " << maxApproximateFlops / gflops / t.lapAvg() << " GFlops" << std::endl;
+
+        results_gpu.readN(gpu_results.ptr(), width * height);
+
+        double realIterationsFraction = 0.0;
+        for (int j = 0; j < (int)height; ++j) {
+            for (int i = 0; i < (int)width; ++i) {
+                realIterationsFraction += gpu_results.ptr()[j * width + i];
+            }
+        }
+        std::cout << "    Real iterations fraction: "
+                  << 100.0 * realIterationsFraction / (width * height)
+                  << "%" << std::endl;
+
+        renderToColor(gpu_results.ptr(), image.ptr(), width, height);
+        image.savePNG("../mandelbrot_gpu.png");
+    }
+
+    {
+        double errorAvg = 0.0;
+        for (int j = 0; j < height; ++j) {
+            for (int i = 0; i < width; ++i) {
+                errorAvg += fabs(gpu_results.ptr()[j * width + i] - cpu_results.ptr()[j * width + i]);
+            }
+        }
+        errorAvg /= width * height;
+        std::cout << "GPU vs CPU average results difference: " << 100.0 * errorAvg << "%" << std::endl;
+
+        if (errorAvg > 0.03) {
+            throw std::runtime_error("Too high difference between CPU and GPU results!");
+        }
+    }
 
     // Это бонус в виде интерактивной отрисовки, не забудьте запустить на ГПУ, чтобы посмотреть, в какой момент числа итераций/точности single float перестанет хватать
     // Кликами мышки можно смещать ракурс

src/main_sum.cpp

@@ -1,6 +1,14 @@
 #include <libutils/misc.h>
 #include <libutils/timer.h>
 #include <libutils/fast_random.h>
+#include <libgpu/context.h>
+#include <libgpu/shared_device_buffer.h>
+
+#include "cl/sum_cl.h"
+
+#include <libutils/misc.h>
+#include <libutils/timer.h>
+#include <libutils/fast_random.h>
 
 
 template<typename T>
@@ -14,9 +22,43 @@ void raiseFail(const T &a, const T &b, std::string message, std::string filename
 
 #define EXPECT_THE_SAME(a, b, message) raiseFail(a, b, message, __FILE__, __LINE__)
 
+void runTest(const std::string &kernel_name, std::vector<unsigned int> &as, unsigned int ref_sum)
+{
+    gpu::gpu_mem_32u as_gpu;
+    unsigned int n = as.size();
+    unsigned int work_group_size = 256;
+    unsigned int global_work_size = n + work_group_size - 1;
+    const int benchmarkingIters = 100;
+
+    as_gpu.resizeN(n);
+    as_gpu.writeN(as.data(), n);
+
+    ocl::Kernel sum_gpu_kernel(sum_kernel, sum_kernel_length, kernel_name);
+    sum_gpu_kernel.compile();
+
+    timer t;
+    for (int iter = 0; iter < benchmarkingIters; ++iter) {
+        unsigned int sum = 0;
+
+        gpu::gpu_mem_32u sum_gpu;
+        sum_gpu.resizeN(1);
+        sum_gpu.writeN(&sum, 1);
+        sum_gpu_kernel.exec(gpu::WorkSize(work_group_size, global_work_size), as_gpu, sum_gpu, n);
+        sum_gpu.readN(&sum, 1);
+
+        EXPECT_THE_SAME(ref_sum, sum, "GPU " + kernel_name + " result should be consistent!");
+
+        t.nextLap();
+    }
+
+    std::cout << "[" << kernel_name << "]" << std::endl;
+    std::cout << "    GPU: " << t.lapAvg() << "+-" << t.lapStd() << " s" << std::endl;
+    std::cout << "    GPU: " << (n / 1000.0 / 1000.0) / t.lapAvg() << " millions/s" << std::endl;
+}
 
 int main(int argc, char **argv)
 {
+
     int benchmarkingIters = 10;
 
     unsigned int reference_sum = 0;
@@ -58,7 +100,15 @@ int main(int argc, char **argv)
     }
 
     {
-        // TODO: implement on OpenCL
-        // gpu::Device device = gpu::chooseGPUDevice(argc, argv);
+        gpu::Device device = gpu::chooseGPUDevice(argc, argv);
+        gpu::Context context;
+        context.init(device.device_id_opencl);
+        context.activate();
+
+        runTest("sum_atomic", as, reference_sum);          // 3.2.2
+        runTest("sum_loop", as, reference_sum);            // 3.2.2
+        runTest("sum_loop_coalesced", as, reference_sum);  // 3.2.3
+        runTest("sum_local", as, reference_sum);           // 3.2.4
+        runTest("sum_tree", as, reference_sum);
     }
 }