task03 done

src/cl/mandelbrot.cl

@@ -1,13 +1,41 @@
-#ifdef __CLION_IDE__
-#include <libgpu/opencl/cl/clion_defines.cl>
-#endif
+__kernel void mandelbrot(__global float* results,
+                         unsigned int width, unsigned int height,
+                         float fromX, float fromY,
+                         float sizeX, float sizeY,
+                         unsigned int iters)
+{
+    const float threshold = 256.0f;
+    const float threshold2 = threshold * threshold;
 
-#line 6
+    int i = get_global_id(0);
+    int j = get_global_id(1);
 
-__kernel void mandelbrot(...)
-{
-    // TODO если хочется избавиться от зернистости и дрожания при интерактивном погружении, добавьте anti-aliasing:
-    // грубо говоря, при anti-aliasing уровня N вам нужно рассчитать не одно значение в центре пикселя, а N*N значений
-    // в узлах регулярной решетки внутри пикселя, а затем посчитав среднее значение результатов - взять его за результат для всего пикселя
-    // это увеличит число операций в N*N раз, поэтому при рассчетах гигаплопс антиальясинг должен быть выключен
+    if(i >= width)
+        return;
+    if(j >= height)
+        return;
+
+    float x0 = fromX + (i + 0.5f) * sizeX / width;
+    float y0 = fromY + (j + 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 = iter;
+
+    result = 1.0f * result / iters;
+    results[j * width + i] = result;
+//     TODO если хочется избавиться от зернистости и дрожания при интерактивном погружении, добавьте anti-aliasing:
+//     грубо говоря, при anti-aliasing уровня N вам нужно рассчитать не одно значение в центре пикселя, а N*N значений
+//     в узлах регулярной решетки внутри пикселя, а затем посчитав среднее значение результатов - взять его за результат для всего пикселя
+//     это увеличит число операций в N*N раз, поэтому при рассчетах гигаплопс антиальясинг должен быть выключен
 }

src/cl/sum.cl

@@ -1 +1,106 @@
-// TODO
+#define VALUES_PER_WORK_ITEM 32
+#define WORKGROUP_SIZE 128
+
+__kernel void atomic_sum(__global const int *arr,
+                         __global unsigned int *sum,
+                         unsigned int n)
+{
+    unsigned int id = get_global_id(0);
+    if (id < n)
+    {
+        atomic_add(sum, arr[id]);
+    }
+}
+
+__kernel void loop_sum(__global const int *arr,
+                      __global unsigned int *sum,
+                      unsigned int n)
+{
+    const unsigned int idx = get_global_id(0);
+    unsigned int res = 0;
+    for (int i = idx * VALUES_PER_WORK_ITEM; i < (idx + 1) * VALUES_PER_WORK_ITEM; ++i)
+    {
+        if (i < n)
+        {
+            res += arr[i];
+        }
+    }
+
+    atomic_add(sum, res);
+}
+
+__kernel void loop_coalesced_sum(__global const int *arr,
+                                 __global unsigned int *sum,
+                                 unsigned int n)
+{
+    const unsigned int lid = get_local_id(0);
+    const unsigned int wid = get_group_id(0);
+    const unsigned int grs = get_local_size(0);
+
+    unsigned int res = 0;
+    for (int i = 0; i < VALUES_PER_WORK_ITEM; ++i)
+    {
+        int idx = wid * grs * VALUES_PER_WORK_ITEM + i * grs + lid;
+        if (idx < n)
+        {
+            res += arr[idx];
+        }
+    }
+
+    atomic_add(sum, res);
+}
+
+__kernel void sum_local_mem(__global const int *arr,
+                            __global unsigned int *sum,
+                            unsigned int n)
+{
+    const unsigned int gid = get_global_id(0);
+    const unsigned int lid = get_local_id(0);
+
+    __local unsigned int buf[WORKGROUP_SIZE];
+
+    buf[lid] = gid < n ? arr[gid] : 0;
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    if (lid == 0)
+    {
+        unsigned int group_res = 0;
+        for (unsigned int i = 0; i < WORKGROUP_SIZE; ++i)
+        {
+            group_res += buf[i];
+        }
+
+        atomic_add(sum, group_res);
+    }
+}
+
+__kernel void tree_sum(__global const int *arr,
+                       __global unsigned int *sum,
+                       const unsigned int n)
+{
+    const unsigned int lid = get_local_id(0);
+    const unsigned int gid = get_global_id(0);
+
+    __local unsigned int buf [WORKGROUP_SIZE];
+    buf[lid] = gid < n ? arr[gid] : 0;
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    for (int nValues = WORKGROUP_SIZE; nValues > 1; nValues /= 2)
+    {
+        if (2 * lid < nValues)
+        {
+            unsigned int a = buf[lid];
+            unsigned int b = buf[lid + nValues / 2];
+            buf[lid] = a + b;
+        }
+        barrier(CLK_LOCAL_MEM_FENCE);
+    }
+
+    if (lid == 0)
+    {
+        atomic_add(sum, buf[0]);
+    }
+}
+
+

src/main_mandelbrot.cpp

@@ -106,44 +106,71 @@ int main(int argc, char **argv)
     }
 
 
-//    // Раскомментируйте это:
-//
-//    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);
+        // TODO близко к ЦПУ-версии, включая рассчет таймингов, гигафлопс, Real iterations fraction и сохранение в файл
+        // результат должен оказаться в gpu_results
+
+        gpu::gpu_mem_32f gpu_img;
+        gpu_img.resizeN(width * height);
+
+        unsigned int workGroupSize = 16;
+        unsigned int global_work_size_w = (width + workGroupSize - 1) / workGroupSize * workGroupSize;
+        unsigned int global_work_size_h = (height + workGroupSize - 1) / workGroupSize * workGroupSize;
+        timer t;
+        for (int i = 0; i < benchmarkingIters; ++i) {
+            kernel.exec(gpu::WorkSize(workGroupSize, workGroupSize, global_work_size_w, global_work_size_h),
+                        gpu_img,
+                        width, height,
+                        centralX - sizeX / 2.0f, centralY - sizeY / 2.0f,
+                        sizeX, sizeY,
+                        iterationsLimit);
+            t.nextLap();
+        }
+
+        size_t flopsInLoop = 10;
+        size_t maxApproximateFlops = width * height * iterationsLimit * flopsInLoop;
+        size_t gflops = 1000*1000*1000;
+        std::cout << "GPU: " << t.lapAvg() << "+-" << t.lapStd() << " s" << std::endl;
+        std::cout << "GPU: " << maxApproximateFlops / gflops / t.lapAvg() << " GFlops" << std::endl;
+
+        gpu_img.readN(gpu_results.ptr(), width * height);
+        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 перестанет хватать
     // Кликами мышки можно смещать ракурс
-    // Но в Pull-request эти две строки должны быть закомментированы, т.к. на автоматическом тестировании нет оконной подсистемы 
+    // Но в Pull-request эти две строки должны быть закомментированы, т.к. на автоматическом тестировании нет оконной подсистемы
 //    bool useGPU = false;
 //    renderInWindow(centralX, centralY, iterationsLimit, useGPU);
 

src/main_sum.cpp

@@ -1,6 +1,12 @@
 #include <libutils/misc.h>
 #include <libutils/timer.h>
 #include <libutils/fast_random.h>
+#include <libgpu/context.h>
+#include <libgpu/shared_device_buffer.h>
+#include <list>
+#include <string>
+
+#include "cl/sum_cl.h"
 
 
 template<typename T>
@@ -59,6 +65,47 @@ 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();
+
+        gpu::gpu_mem_32u as_buffer;
+        gpu::gpu_mem_32u sum_buffer;
+
+        as_buffer.resizeN(n);
+        sum_buffer.resizeN(1);
+
+        unsigned int workGroupSize = 128;
+        unsigned int global_work_size = (n + workGroupSize - 1) / workGroupSize * workGroupSize;
+        unsigned int n_work_groups = global_work_size / workGroupSize;
+
+        as_buffer.writeN(as.data(), n);
+        const unsigned int init = 0;
+        std::list<std::string> kernel_names = {"atomic_sum",
+                                               "loop_sum",
+                                               "loop_coalesced_sum",
+                                               "sum_local_mem",
+                                               "tree_sum"};
+
+        std::cout << std::endl;
+        for(auto it_kernel_name = kernel_names.begin(); it_kernel_name != kernel_names.end(); ++it_kernel_name)
+        {
+            ocl::Kernel kernel(sum_kernel, sum_kernel_length, it_kernel_name->data());
+            kernel.compile();
+
+            timer t;
+            unsigned int sum = 0;
+            for (int iter = 0; iter < benchmarkingIters; ++iter) {
+                sum_buffer.writeN(&init, 1);
+                kernel.exec(gpu::WorkSize(workGroupSize, global_work_size),
+                            as_buffer, sum_buffer, n);
+                t.nextLap();
+            }
+            sum_buffer.readN(&sum, 1);
+            EXPECT_THE_SAME(reference_sum, sum, "GPU result should be consistent!");
+            std::cout << "GPU " << *it_kernel_name << " : " <<  t.lapAvg() << "+-" << t.lapStd() << " s" << std::endl;
+            std::cout << "GPU " << *it_kernel_name << " : " << (n/1000.0/1000.0) / t.lapAvg() << " millions/s" << std::endl << std::endl;
+        }
     }
 }