Фрактали, Fortran та OpenMP

Коли-то давно я вирішив «помацати» Fortran. Єдину задачу яку я придумав — генерація фракталів (заодно і OpenMP в Fortran'е можна було б спробувати). В процесі написання я часто стикався з проблемами, вирішення яких доводилося додумувати самому (наприклад в інтернеті не так багато прикладів використання чисел подвійної точності або бінарної запису у файл). Але рано чи пізно всі проблеми вирішилися, і я хочу написати цей текст, який можливо комусь допоможе.

Писати я буду на діалекті Fortran 90, але з GNU розширеннями (ті ж числа подвійної точності).


Трохи теорії про фрактали

Фрактал (fractus — подрібнений, зламаний, розбитий) — у вузькому сенсі, складна геометрична фігура, що володіє властивістю самоподібності, тобто складена з кількох частин, кожна з яких подібна всій фігурі

Існує велика група фракталів — алгебраїчні фрактали. Це назва випливає з принципу побудови таких фракталів. Для їх побудови використовують рекурсивні функції. Наприклад алгебраїчними фракталами є: Безліч Жюлиа, Безліч Мандельброта, Басейн(фрактал) Ньютона.

Побудова алгебраїчних фракталів

Один з методів побудови являє собою ітераційний розрахунок image, де imageimageякась функція. Розрахунок даної функції триває до виконання певної умови. Коли це умова виконується на екран виводиться крапка.

Поведінка функції для різних точок комплексної площини може мати різну поведінку:
  • Прагне до нескінченності
  • Прагне до 0
  • Приймає кілька фіксованих значень і не виходить за їх межі
  • Хаотичне поводження. Відсутність будь-яких тенденцій
Один з найпростіших (як для розуміння, так і для реалізації) алгоритмів побудови алгебраїчних фракталів відомий як «escape time algorithm». Якщо коротко, то проводиться ітеративний розрахунок числа для кожної точки комплексної площини.

Було доведено, що якщо imageвиявиться більше bailout, то послідовність image. Порівняння imagebailout дозволяє знаходити точки лежать поза безлічі. Для точок, що лежать поза безлічі, послідовність imageне буде мати тенденції до нескінченності, тому ніколи не досягне 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) пор. (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

Повинно вийде наступне зображення:
1024x1024

Як згенерувати Множину Мандельброта

Це зробити нескладно, достатньо замінити 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

Повинно вийде наступне зображення:
1024x1024

Використана література

en.wikipedia.org/wiki/Fractal
en.wikipedia.org/wiki/Julia_set
en.wikipedia.org/wiki/Mandelbrot_set
en.wikipedia.org/wiki/OpenMP
openmp.org/wp/openmp-specifications/
gcc.gnu.org/onlinedocs/gfortran/

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.