Множество Мандельброта - множество комплексных чисел для которых не расходится (начиная с ), т.е. множество комплексных чисел, для которых остается ограниченным.

Полезный факт: если , то не принадлежит множеству Мондельброта.

Если и , то .

Хочется визуализировать множество таких точек на -плоскости (псевдокод):

mandelbrotContains(x0, y0):
    x, y = x0, y0
    for iteration in 0...MAX_ITERATIONS:
        x, y = x*x - y*y + x0, 2*x*y + y0
        if sqrt(x*x+y*y) > 2:
            return false
    return true

drawMandelbrot(image):
    for x in image:
        for y in image:
            if mandelbrotContains(x, y):
                image[x, y] = true
            else:
                image[x, y] = false

Множество Мандельброта

Давайте реализуем вычисление количества итераций, которое успевает произойти в функции mandelbrotContains(x, y) до ее завершения. Если число итераций достигает MAX_ITERATIONS - эта точка почти наверняка принадлежит нашему множеству (т.к. ряд не расходится достаточно долго).

Реализация на C++ ниже:

void MandelbrotProcessorCPU::process(Vector2f from, Vector2f to,
                                     Image<unsigned short>& iterations) {
    size_t width = iterations.width;
    size_t height = iterations.height;

    float x_step = (to.x() - from.x()) / width;
    float y_step = (to.y() - from.y()) / height;

    for (size_t py = 0; py < height; py++) {
        float y0 = from.y() + y_step * py;

        for (size_t px = 0; px < width; px++) {
            float x0 = from.x() + x_step * px;

            unsigned short iteration;
            float x = 0.0f;
            float y = 0.0f;
            for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
                float xn = x * x - y * y + x0;
                y = 2 * x * y + y0;
                x = xn;
                if (x * x + y * y > INFINITY) {
                    break;
                }
            }
            iterations(py, px) = iteration;
        }
    }
}

Замечания: from и to - координаты углов изображения на комлексной плоскости, INFINITY - порог расхождения ряда (равен 2.0f).

Все замеры производительности ниже реализованы в этом исходнике, но запускается он с разрешением 2048x2048 (MAX_ITERATIONS = 10000 как на github).

На моем Intel i7 6700 данная реализация выполняется 5545 ms.

OpenMP праллелизация

Но в моем CPU четыре ядра, и восемь потоков из-за hyper-threading, поэтому давайте распараллелим вычисления с помощью OpenMP. Код:

void MandelbrotProcessorCPU::process(Vector2f from, Vector2f to,
                                     Image<unsigned short>& iterations) {
    size_t width = iterations.width;
    size_t height = iterations.height;

    float x_step = (to.x() - from.x()) / width;
    float y_step = (to.y() - from.y()) / height;
    
    // The only line needed to add:
    #pragma omp parallel for
    for (size_t py = 0; py < height; py++) {
        float y0 = from.y() + y_step * py;

        for (size_t px = 0; px < width; px++) {
            float x0 = from.x() + x_step * px;

            unsigned short iteration;
            float x = 0.0f;
            float y = 0.0f;
            for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
                float xn = x * x - y * y + x0;
                y = 2 * x * y + y0;
                x = xn;
                if (x * x + y * y > INFINITY) {
                    break;
                }
            }
            iterations(py, px) = iteration;
        }
    }
}

Единственная важная строка кода - #pragma omp parallel for. OpenMP позволяет распараллелить вычисления с помощью таких простых прагм. Данный фреймворк поддерживается почти всеми компиляторами (разве что возможно на Mac OS X он не поддерживается из коробки).

В данном случае OpenMP распределяет вычисления между всем ядрами. Прагма была указана для внешнего for-loop - поэтому блок кода внешнего цикла будет исполнятся с некоторым значением py в одном из потоков. Например для каждого восьмого значения py блок кода может исполняться в Потоке#8, каждый седьмой из восьми - в Потоке#7 и так далее. Поэтому логично ожидать ускорение вплоть до x8 раз.

Теперь время исполнения равно 1274 ms. Ускорение: 5545 ms / 1274 ms = x4.3. Это не те x8 которые ожидались, и вызвано это тем что на самом деле в процессоре только 4 реальных ядра. Hyper-threading одного ядра может привести к ускорение только если два паралелльных потока используют разные ресурсы этого ядра. Узкое место в данном примере в вычислениях - и эти ресурсы ядер не могут быть магическим образом удвоены за счет hyper-threading.

Итого на данный момент:

Speedup  Time    Implementation   Device

 x4.3   1274 ms	 Naive CPU        8 threads		
 x1.0   5545 ms	 Naive CPU        1 thread	

Cpu is Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz (4 cores, 8 threads by hyper-threading)

SIMD, SSE

Уже очень давно все процессоры поддерживают SIMD инструкции. SIMD (Single Instruction, Multiple Data) позволяют произвести одну операцию сразу над несколькими операндами. Например можно сложить с и получить , где . Ориентировочно можно считать, что если инструкция обрабатывает N операндов - она даст прирост вплоть до xN раз.

Наиболее известные и популярные расширения - SSE, SSE2 и вплоть до SSE4.2. Они поддерживаются практически во всех x86 процессорах и оперируют 128-битными данными. В случае работы с 32-битными числами с плавающей точкой это означает ускорение вплоть до x4 раз.

Все современные компиляторы поддерживают intrinsics для таких расширений. Пример SSE инструкции: _mm_add_ps - суммирует четыре числа с плавающей точкой (числа должны располагаться в специальном SSE регистре) с другими четырьмя числами (из другого регистра) и сохраняет результат в третий регистр. Более детальное описание всех инструкций удобно изложено в Intel intrinsics guide. Суффикс инструкции показывает над каким типом данных выполняется работа. Например ps - над числами с плавающей точкой одинарной точности, pi16 - 16-битные целые числа и т.п..

Т.к. эти intrinsics работают с четырьмя 32-битными операндами одновременно - нам надо выполнять вычисления над данными расположенными в специальных больших регистрах, а затем выгрузить результат из регистра в память.

C++ компилятор сам распределит наши вычисления по регистрам, но для этого требуется явным образом работать с четырьмя значениями за раз - с помощью использования специального типа: __m128 - 128-битного регистра.

Итоговый код:

void MandelbrotProcessorCPU_SSE::process(Vector2f from, Vector2f to,
                                         Image<unsigned short>& iterations) {
    float width = iterations.width;
    float height = iterations.height;

    float x_step = (to.x() - from.x()) / width;
    float y_step = (to.y() - from.y()) / height;

    #pragma omp parallel for
    for (size_t py = 0; py < iterations.height; py++) {
        float y0 = from.y() + y_step * py;

        for (size_t px = 0; px < iterations.width / 4 * 4; px += 4) {
            float pxf = (float) px;
                                                                                                                    // Four values in register (__m128):
            __m128 pxs_deltas128 = _mm_mul_ps(_mm_set_ps(0.0f, 1.0f, 2.0f, 3.0f), _mm_set1_ps(x_step));             // 0.0f*x_step, 1.0f*x_step, 2.0f*x_step, 3.0f*x_step 
            __m128 xs0 = _mm_add_ps(_mm_set1_ps(from.x()), _mm_add_ps(_mm_set1_ps(x_step * pxf), pxs_deltas128));   // from.x()+px*x_step, where px takes 4 conseqent values starting from current loop index

            unsigned short iteration;
            __m128i maskAll = _mm_setzero_si128();                                                                  // maskAll is the mask, that stores flag for each value from our four:
            __m128i iters = _mm_setzero_si128();                                                                    // 'is this series divergent?' - sqrt(x*x+y*y) > 2
            __m128 xs = _mm_set1_ps(0.0f);                                                                          // because in that case we should not increment 'iters' for that value
            __m128 ys = _mm_set1_ps(0.0f);
            for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
                __m128 xsn = _mm_add_ps(_mm_sub_ps(_mm_mul_ps(xs, xs), _mm_mul_ps(ys, ys)), xs0);                   // xn = x * x - y * y + x0;
                __m128 ysn = _mm_add_ps(_mm_mul_ps(_mm_mul_ps(_mm_set1_ps(2.0f), xs), ys), _mm_set1_ps(y0));        // yn = 2 * x * y + y0;
                xs = _mm_add_ps(_mm_andnot_ps((__m128) maskAll, xsn), _mm_and_ps((__m128) maskAll, xs));            // updating previous values to new ones
                ys = _mm_add_ps(_mm_andnot_ps((__m128) maskAll, ysn), _mm_and_ps((__m128) maskAll, ys));            // but with taking into account maskAll
                                                                                                                    // (if value divergent once - we should fix its iterations number)

                maskAll = _mm_or_si128(_mm_castps_si128(_mm_cmpge_ps(_mm_add_ps(_mm_mul_ps(xs, xs), _mm_mul_ps(ys, ys)), _mm_set1_ps(INFINITY))), maskAll);  // calculating mask, by checking x*x+y*y > INFINITY for each value
                iters = _mm_add_epi16(iters, _mm_andnot_si128(maskAll, _mm_set1_epi16(1)));                         // increments 'iters' counter for those numbers,
                int mask = _mm_movemask_epi8(maskAll);                                                              // that were taking part in this iteration
                if (mask == 0xffff) {
                    break;                                                                                          // if all values divergent - stop
                }
            }
            iters = _mm_shuffle_epi8(iters, _mm_setr_epi8(12, 13, 8, 9, 4, 5, 0, 1, -1, -1, -1, -1, -1, -1, -1, -1));
            _mm_storel_epi64((__m128i *) &iterations(py, px), iters);                                               // unloading data to memory
        }

        for (size_t px = iterations.width / 4 * 4; px < iterations.width; px++) {
            float x0 = from.x() + (to.x() - from.x()) * px / width;

            unsigned short iteration;
            float x = 0.0f;
            float y = 0.0f;
            for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
                float xn = x * x - y * y + x0;
                y = 2 * x * y + y0;
                x = xn;
                if (x * x + y * y > INFINITY) {
                    break;
                }
            }
            iterations(py, px) = iteration;
        }
    }
}

Да, довольно нечитаемый код! Написание такого кода довольно муторное занятие (если делать это лишь изредко). Intel intrinsics guide помогает очень сильно - в нем есть очень удобные фильтры по наборам расширений (например вам могут быть не интересны инструкции из AVX-набора, т.к. не все процессоры ваших пользователей поддерживают данный набор), там же есть исчерпывающие описания и т.п..

Обратите внимание, что последний цикл является откатом к невекторизованной версии, т.к. объем наших данных может не делиться на 4 (это количество элементов обрабатываемых за раз SSE инструкциями) и это самый простой и очевидный способ обработать такую ситуацию аккуратно.

Результаты производительности (SSE, SSE2 и SSSE3 были использованы):

Speedup   Time   Implementation   Device

x10.1    547 ms	 CPU SSE          8 threads		
 x2.3   2370 ms	 CPU SSE          1 thread	
 x4.3   1274 ms	 Naive CPU        8 threads		
 x1.0   5545 ms	 Naive CPU        1 thread	
 
Cpu is Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz (4 cores, 8 threads by hyper-threading)

SIMD, AVX

Но можно пойти дальше! После SSE было выпущено много расширений, в т.ч. AVX (заявлен в 2008, поддерживается с 2011 - впервые в Sandy Bridge). Главным отличием являются более широкие (в два раза) регистры - 256-битные, т.о. можно ускорить еще в два раза! Обратите внимание, что переход с SSE на AVX довольно тривиальный - достаточно заменить __m128 на __m256 и обновить код с учетом мелких изменений в инструкциях (например иногда порядок аргументов различается).

Код практически идентичен коду с использованием SSE, но на всякий случай:

void MandelbrotProcessorCPU_AVX::process(Vector2f from, Vector2f to,
                                         Image<unsigned short>& iterations) {
    assert (MAX_ITERATIONS < std::numeric_limits<unsigned short>::max());
    assert (iterations.cn == 1);

    assert (((float) iterations.width) + 1 == (float) (iterations.width + 1));
    assert (((float) iterations.height) + 1 == (float) (iterations.height + 1));
    float width = iterations.width;
    float height = iterations.height;

    float x_step = (to.x() - from.x()) / width;
    float y_step = (to.y() - from.y()) / height;

    #pragma omp parallel for
    for (size_t py = 0; py < iterations.height; py++) {
        float y0 = from.y() + y_step * py;

        for (size_t px = 0; px < iterations.width / 8 * 8; px += 8) {
            float pxf = (float) px;
            __m256 pxs_deltas128 = _mm256_mul_ps(_mm256_set_ps(0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f), _mm256_set1_ps(x_step));
            __m256 xs0 = _mm256_add_ps(_mm256_set1_ps(from.x()), _mm256_add_ps(_mm256_set1_ps(x_step * pxf), pxs_deltas128));    // from.x() + x_step * px

            unsigned short iteration;
            __m256i maskAll = _mm256_setzero_si256();
            __m256i iters = _mm256_setzero_si256();
            __m256 xs = _mm256_set1_ps(0.0f);
            __m256 ys = _mm256_set1_ps(0.0f);
            for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
                __m256 xsn = _mm256_add_ps(_mm256_sub_ps(_mm256_mul_ps(xs, xs), _mm256_mul_ps(ys, ys)), xs0);               // xn = x * x - y * y + x0;
                __m256 ysn = _mm256_add_ps(_mm256_mul_ps(_mm256_mul_ps(_mm256_set1_ps(2.0f), xs), ys), _mm256_set1_ps(y0)); // yn = 2 * x * y + y0;
                xs = _mm256_add_ps(_mm256_andnot_ps((__m256) maskAll, xsn), _mm256_and_ps((__m256) maskAll, xs));
                ys = _mm256_add_ps(_mm256_andnot_ps((__m256) maskAll, ysn), _mm256_and_ps((__m256) maskAll, ys));

                maskAll = (__m256i) _mm256_or_ps(_mm256_cmp_ps(_mm256_add_ps(_mm256_mul_ps(xs, xs), _mm256_mul_ps(ys, ys)), _mm256_set1_ps(INFINITY), _CMP_GT_OS), (__m256) maskAll);
                iters = _mm256_add_epi16(iters, _mm256_andnot_si256(maskAll, _mm256_set1_epi16(1)));
                int mask = _mm256_movemask_epi8(maskAll);
                if (mask == (int) 0xffffffff) {
                    break;
                }
            }
            iters = _mm256_shuffle_epi8(iters, _mm256_setr_epi8(0, 1, -1, -1, 4, 5, -1, -1, 8, 9, -1, -1, 12, 13, -1, -1, 16, 17, -1, -1, 20, 21, -1, -1, 24, 25, -1, -1, 28, 29, -1, -1));
            unsigned int tmp[8];
            _mm256_storeu_si256((__m256i *) tmp, iters);
            for (int i = 0; i < 8; i++) {
                iterations(py, px + 7 - i) = (unsigned short) tmp[i];
            }
        }

        for (size_t px = iterations.width / 8 * 8; px < iterations.width; px++) {
            float x0 = from.x() + (to.x() - from.x()) * px / width;

            unsigned short iteration;
            float x = 0.0f;
            float y = 0.0f;
            for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
                float xn = x * x - y * y + x0;
                y = 2 * x * y + y0;
                x = xn;
                if (x * x + y * y > INFINITY) {
                    break;
                }
            }
            iterations(py, px) = iteration;
        }
    }
}

Производительность:

Speedup   Time   Implementation   Device

x20.8    266 ms	 CPU AVX2         8 threads		
 x4.8   1160 ms	 CPU AVX2         1 thread	
x10.1    547 ms	 CPU SSE          8 threads		
 x2.3   2370 ms	 CPU SSE          1 thread	
 x4.3   1274 ms	 Naive CPU        8 threads		
 x1.0   5545 ms	 Naive CPU        1 thread	
 
Cpu is Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz (4 cores, 8 threads by hyper-threading)

AVX действительно дал x2 ускорение относительно SSE.

OpenCL

OpenCL предоставляет единообразный способ реализации паралелльных вычислений для исполнения на множестве многоядерного железа - GPUs и CPUs. Ключевой код содержится в так называемом kernel (ядре), синтаксис в нем основан на C99 и кернелы компилируются прямо на лету (как OpenGL шейдеры). Приложение шлет исходный код OpenCL-драйверу (например драйверу видеокарты, или специальному драйверу процессора) чтобы скомпилировать и затем исполнить kernel.

Код очень похож на изначальную версию:

__kernel void mandelbrotProcess(__global unsigned short* iterations,
                                 int width, int height,
                                 float fromX, float fromY, float toX, float toY) {
    int px = get_global_id(0);
    int py = get_global_id(1);
    if (px >= width || py >= height) return;

    float x0 = fromX + px * (toX - fromX) / width;
    float y0 = fromY + py * (toY - fromY) / height;

    unsigned short iteration;
    float x = 0.0f;
    float y = 0.0f;
    for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
        float xn = x * x - y * y + x0;
        y = 2 * x * y + y0;
        x = xn;
        if (x * x + y * y > INFINITY_THRESHOLD) {
            break;
        }
    }
    iterations[width * py + px] = iteration;
}

Но производительность на AMD R9 390X впечатляет - 6 ms! x924 ускорение! Причина в том, что рассматриваемый пример идеален для архитектуры видеокарты - он идеально параллелизуем и единственное узкое место - вычисления, нет никакого упора в доступе к памяти, нет условного ветвления, только чистые вычисления которые идеально раскладываются по регистрам и вычисления выполняются с полной вычислительной мощностью устройства. В случае AMD R9 390X вычислениями занимаются 2816 ядер работающих на частоте 1 Ghz!

OpenCL kernels так же можно исполнять на процессорах - тест содержит выполнение на i7-6700 через два OpenCL-драйвера: версия от Intel и от AMD (драйвер от AMD так же поддерживает процессоры от Intel):

Speedup   Time   Implementation   Device

 x924      6 ms	 OpenCL (AMD)     Hawaii R9 390X
  x22    252 ms	 OpenCL (Intel)   Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz
 x5.8    942 ms	 OpenCL (AMD)     Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz
x20.8    266 ms	 CPU AVX2         8 threads		
 x4.8   1160 ms	 CPU AVX2         1 thread	
x10.1    547 ms	 CPU SSE          8 threads		
 x2.3   2370 ms	 CPU SSE          1 thread	
 x4.3   1274 ms	 Naive CPU        8 threads		
 x1.0   5545 ms	 Naive CPU        1 thread	
 
Cpu is Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz (4 cores, 8 threads by hyper-threading)

И внезапно Intel драйвер работает быстрее чем рукописная AVX-версия! Это означает, что драйвер успешно векторизовал код с учетом доступных на данном процессоре наборов инструкций: AVX и AVX2. AMD драйвер вероятно исполняет невекторизованную версию на всех доступных ядрах, и это приводит к такой же производительности как и у наивного распараллеливания.

OpenCL позволяет использовать вычислительные мощности центральных процессоров и видеокарт. И это гораздо более простой способ, по сравнению с использованием процессорных intrinsics. Даже на CPU хорошая производительность часто может быть достигнута с прямолинейным и простым кодом благодаря OpenCL-драйверу (т.к. Intel-драйвер часто способен эффективно векторизовать код). Но иногда intrinsics необходимы: если не у всех пользователей есть GPU и слишком сложно заставить их установить OpenCL-драйвер от Intel, или если у вас есть много простаивающих компьютеров с мощными процессорами (например простаивающие сервера с Xeon-процессорами) а Intel-драйвер не сумел векторизовать ваш слишком сложный код, и т.д..