diff --git a/mandelbrot_cpu.png b/mandelbrot_cpu.png new file mode 100644 index 00000000..c714657c Binary files /dev/null and b/mandelbrot_cpu.png differ diff --git a/mandelbrot_gpu.png b/mandelbrot_gpu.png new file mode 100644 index 00000000..bb3c5ec9 Binary files /dev/null and b/mandelbrot_gpu.png differ diff --git a/src/cl/mandelbrot.cl b/src/cl/mandelbrot.cl index 8c2ff5ae..4c8d2aef 100644 --- a/src/cl/mandelbrot.cl +++ b/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; } diff --git a/src/cl/sum.cl b/src/cl/sum.cl index 0ffdd02f..3c2504b8 100644 --- a/src/cl/sum.cl +++ b/src/cl/sum.cl @@ -1 +1,94 @@ -// TODO \ No newline at end of file +#ifdef __CLION_IDE__ +#include +#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]); + } +} diff --git a/src/main_mandelbrot.cpp b/src/main_mandelbrot.cpp index 5ebbbe3f..dc90e385 100644 --- a/src/main_mandelbrot.cpp +++ b/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 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 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 перестанет хватать // Кликами мышки можно смещать ракурс diff --git a/src/main_sum.cpp b/src/main_sum.cpp index 735fd72f..bdf8af7d 100644 --- a/src/main_sum.cpp +++ b/src/main_sum.cpp @@ -1,6 +1,14 @@ #include #include #include +#include +#include + +#include "cl/sum_cl.h" + +#include +#include +#include template @@ -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 &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); } }