Фракталы, Fortran и OpenMP
Когда-то давно я решил «потрогать» Fortran. Единственную задачу которую я придумал — генерация фракталов (заодно и OpenMP в Fortran’е можно было бы попробовать). В процессе написания я часто сталкивался с проблемами, решение которых приходилось додумывать самому (например в интернете не так много примеров использования чисел двойной точности или бинарной записи в файл). Но рано или поздно все проблемы решились, и я хочу написать этот текст, который возможно кому-нибудь поможет.Писать я буду на диалекте Fortran 90, но с GNU расширениями (те же числа двойной точности).
Немного теории о фракталахФрактал (fractus — дробленый, сломанный, разбитый) — в узком смысле, сложная геометрическая фигура, обладающая свойством самоподобия, то есть составленная из нескольких частей, каждая из которых подобна всей фигуре целиком.Существует большая группа фракталов — алгебраические фракталы. Это название вытекает из принципа построения таких фракталов. Для их построения используют рекурсивные функции. Например алгебраическими фракталами являются: Множество Жюлиа, Множество Мандельброта, Бассейн (фрактал) Ньютона.
Построение алгебраических фракталов
Один из методов построения представляет собой итерационный расчет , где
, а
некая функция. Расчет данной функции продолжается до выполнения определенного условия. Когда это условие выполнится на экран выводится точка.Поведение функции для разных точек комплексной плоскости может иметь разное поведение:
Стремится к бесконечности
Стремится к 0
Принимает несколько фиксированных значений и не выходит за их пределы
Хаотичное поведение. Отсутствие каких либо тенденций
Один из простейших (как для понимания, так и для реализации) алгоритмов построения алгебраических фракталов известен как «escape time algorithm». Если кратко, то производится итеративный расчет числа для каждой точки комплексной плоскости.Было доказано, что если окажется больше bailout, то последовательность
. Сравнение
с bailout позволяет находить точки лежащие вне множества. Для точек, лежащих вне множества, последовательность
не будет иметь тенденции к бесконечности, поэтому никогда не достигнет bailout.
Я рассмотрю два простых фрактала — Множество Жюлиа и Множество Мандельброта. Расчитываются они по одной и той-же функции, а отличаются лишь константой, где у Жюлиа это постоянная константа, а у Мандельброба константа зависит от точки комплексной плоскости.
Построение Множества Жюлиа Начало и конец программы простой и тривиальный: program frac end Дальше нам нужно ввести несколько констант: INTEGER, PARAMETER: iterations = 2048! Количество итераций. Чем больше, тем глубже будет идти просчет
REAL (8), PARAMETER: mag_ratio = 1.0_8! Приближение
REAL (8), PARAMETER: x0 = 0.0_8! Смещение комплексной плоскости REAL (8), PARAMETER: y0 = 0.0_8!
INTEGER, PARAMETER: resox = 1024! Разрешение получившегося изображения INTEGER, PARAMETER: resoy = resox
REAL (8), PARAMETER: xshift = (-1.5_8 / mag_ratio) + x0! Смещение центра комплексной оси в REAL (8), PARAMETER: yshift = (-1.5_8 / mag_ratio) + y0! центр изображения
REAL (8), PARAMETER: CXmin = -1.5_8 / mag_ratio! Следующие константы растягивают REAL (8), PARAMETER: CXmax = 1.5_8 / mag_ratio! комплексную плоскость до размеров resox x resoy REAL (8), PARAMETER: CXshag = (DABS (CXmin) + DABS (CXmax)) / resox! Связываем значение одного пикселя и точки на комплексной плоскости REAL (8), PARAMETER: CYmin = -1.5_8 / mag_ratio REAL (8), PARAMETER: CYmax = 1.5_8 / mag_ratio REAL (8), PARAMETER: CYshag = (DABS (CYmin) + DABS (CYmax)) / resoy
COMPLEX (8), PARAMETER: c = DCMPLX (0.285_8, 0.01_8) ! Наша основная константа И вот тут нужно кое-что объяснить. REAL и COMPLEX это число с плавающей точкой одинарной точности и комплексное число, состоящее из двух REAL. REAL (8) это число двойной точности, а COMPLEX (8), соотвественно, комплексное число, состоящее из двух REAL (8). Так-же к стандартным функциям добавляется литера D (как в случае с ABS), что указывает на использование чисел двойной точности.Дальше нам нужно ввести несколько переменных:
CHARACTER, DIMENSION (:), ALLOCATABLE: matrix! Массив точек
COMPLEX (8) :: z INTEGER: x, y, iter REAL (8) :: iter2! Понадобится нам для сглаживания
ALLOCATE (matrix (resox * resoy * 3)) Использование ALLOCATABLE обязательно! Т.к. в ином случае: CHARACTER, DIMENSION (0: resox * resoy * 3) :: matrix! Массив точек Память у нас выделится на стеке, а это не приемлемо при использовании в нескольких потоках. Поэтому мы выделяем память на куче.Так-же я не использую двумерные массивы, т.к. при выделении на куче массив будет занимать больше места, да и записать его в файл будет сложнее.
Дальше нам нужно указать количество потоков, порождаемых OpenMP:
call omp_set_num_threads (16) А тут начинаются непосредственно вычисления. В Fortran’е директивы для OpenMP указываются через коментарии (в C/C++ для это есть специльный макрос #pragma). !$OMP PARALLEL DO do x = 0, resox-1 do y = 0, resoy-1 iter = 0! Обнуляем количество итераций z = DCMPLX (x * CXshag + xshift, y * CYshag + yshift) ! Задаем Z начальное значение do while ((iter .lt. iterations) .and. (CDABS (z) .le. 4.0_8)) ! Обычно за bailout берут 2, но я взял 4, т.к. результат лучше сглаживается z = z**2 + c iter = iter + 1 end do
iter2 = REAL (iter) — DLOG (DLOG (CDABS (z))) / DLOG (2.0_8) ! Формула для сглаживания iter-ln (log2(|Z|)) matrix ((x + y * resox) * 3 + 0) = CHAR (NINT (DMOD (iter2×7.0_8, 256.0_8))) ! Один из способов расскрашивания matrix ((x + y * resox) * 3 + 1) = CHAR (NINT (DMOD (iter2×14.0_8, 256.0_8))) matrix ((x + y * resox) * 3 + 2) = CHAR (NINT (DMOD (iter2×2.0_8, 256.0_8))) end do end do !$OMP END PARALLEL DO Самая трудоемкая часть закончена. Дальше нам остается вывести информацию в удобном для восприятия виде — изображение. Я буду использовать формат PNM, т.к. это самый простой формат изображений.PNM состоит из нескольких секций:
P6 #Комментарий 1024 1024 255 Первая строчка это указание формата записи информации о пикселях: P1/P4 — черно-белое изображение P2/P5 — серое изображение P3/P6 — цветное изображение Первая P это вариант, где цвет пикселя записывается ASCII символом, а вторая P дает возможность записывать цвет пикселя в бинарном виде (что значительно экономит место на диске).Следующей строкой идет комментарий, дальше разрешение изображения и количество цветов на пиксель. После количества цветов идет непосредственно информация о пикселях.
Начинаем записывать изображение в файл:
open (unit=8, file = trim («test.pnm»)) write (8, '(a)') «P6» !(a) это текстовая строка write (8, '(a)') » write (8, '(I0, a, I0)') resox, ' ', resoy write (8, '(I0)') 255 write (8, *) matrix close (8)
DEALLOCATE (matrix) DEALLOCATE мы выполняем для приличия, ибо при выходе из программы, ОС все равно вернет занятую память системе.Для сборки программы можно использовать компилятор от GNU — gfortran:
gfortran -std=gnu frac.f90 -o frac.run -fopenmp
Должно получится следующее изображение:1024×1024
Как сгенерировать Множество Мандельброта
Это сделать несложно, достаточно заменить c и изменить несколько констант. В данном случаем убираем константу c и добавляем переменную c:
COMPLEX (8) :: z, c
z = DCMPLX (x * CXshag + xshift, y * CYshag + yshift) ! Задаем Z начальное значение
c = z
А так-же меняем старые константы:
INTEGER, PARAMETER: MAIN_RES = 1024
INTEGER, PARAMETER: resox = (MAIN_RES / 3) * 3
INTEGER, PARAMETER: resoy = (resox * 2) / 3
REAL (8), PARAMETER: xshift = (-2.0_8 / mag_ratio) + x0 REAL (8), PARAMETER: yshift = (-1.0_8 / mag_ratio) + y0
REAL (8), PARAMETER: CXmin = -2.0_8 / mag_ratio
REAL (8), PARAMETER: CXmax = 1.0_8 / mag_ratio
REAL (8), PARAMETER: CXshag = (DABS (CXmin) + DABS (CXmax)) / resox
REAL (8), PARAMETER: CYmin = -1.0_8 / mag_ratio
REAL (8), PARAMETER: CYmax = 1.0_8 / mag_ratio
REAL (8), PARAMETER: CYshag = (DABS (CYmin) + DABS (CYmax)) / resoy
Должно получится следующее изображение:1024×1024
Используемая литература
en.wikipedia.org/wiki/Fractalen.wikipedia.org/wiki/Julia_seten.wikipedia.org/wiki/Mandelbrot_seten.wikipedia.org/wiki/OpenMPopenmp.org/wp/openmp-specifications/gcc.gnu.org/onlinedocs/gfortran/
