Fractals, Fortran and OpenMP
Once upon a time I decided to “touch” Fortran. The only task I came up with was generating fractals (at the same time, OpenMP in Fortran could be tried). In the process of writing, I often encountered problems that had to be thought out myself (for example, on the Internet there are not many examples of using double-precision numbers or binary writing to a file). But sooner or later, all the problems were solved, and I want to write this text, which perhaps will help someone.
I will write in the Fortran 90 dialect, but with GNU extensions (the same double precision numbers).
Fractal ( fractus - crushed, broken, broken) - in the narrow sense, a complex geometric shape that has the property of self-similarity, that is, composed of several parts, each of which is similar to the entire figure as a whole.
')
There is a large group of fractals - algebraic fractals. This name derives from the principle of constructing such fractals. For their construction use recursive functions. For example, algebraic fractals are: Julia set, Mandelbrot set, Newton's pool (fractal).
One construction method is an iterative calculation.
where 
, but 
some function. The calculation of this function continues until a certain condition is met. When this condition is fulfilled, a dot is displayed on the screen.
The behavior of the function for different points of the complex plane may have different behavior:
One of the simplest (both for understanding and implementation) algorithms for constructing algebraic fractals is known as the “ escape time algorithm ”. In short, an iterative calculation of the number is made for each point of the complex plane.
It has been proven that if
it turns out more bailout , then the sequence 
. Comparison with bailout allows you to find points lying outside the set. For points outside the set, the sequence 
will not tend to infinity, so never reach bailout .
I will look at two simple fractals - Julia set and Mandelbrot set. They are calculated by the same function, and differ only in a constant, where Julia has a constant, and in Mandelbrob the constant depends on the point of the complex plane.
The beginning and end of the program is simple and trivial:
Next we need to enter several constants:
And here you need to explain something. REAL and COMPLEX is a single-precision floating-point number and a complex number consisting of two REALs . REAL (8) is a double precision number, and COMPLEX (8) , respectively, is a complex number consisting of two REAL (8) . Also, the letter D is added to the standard functions (as is the case with ABS ), which indicates the use of double precision numbers.
Next we need to enter a few variables:
Using ALLOCATABLE is a must! Since otherwise:
The memory will be allocated on the stack, and this is not acceptable when used in multiple threads. Therefore, we allocate memory on the heap.
So I do not use two-dimensional arrays, because when selected on the heap, the array will take up more space, and it will be more difficult to write it to the file.
Next we need to specify the number of threads generated by OpenMP:
And here the calculations begin. In Fortran, OpenMP directives are specified via comments (in C / C ++, there is a special #pragma macro for this).
The most time consuming part is over. Next, it remains for us to display information in a convenient form for perception - the image. I will use the PNM format, since This is the easiest image format.
PNM consists of several sections:
The first line is an indication of the format of recording information about pixels:
The first P is an option where the color of a pixel is recorded with an ASCII character, and the second P allows you to record the color of a pixel in binary form (which saves disk space).
The next line is a comment, then the image resolution and the number of colors per pixel. After the number of colors is directly information about the pixels.
Start writing the image to the file:
DEALLOCATE we perform for decency, because when you exit the program, the OS will still return the occupied memory to the system.
You can use the GNU compiler compiler to build the program:
You should get the following image:
It is easy to do, it is enough to replace c and change several constants. In this case, remove the constant c and add the variable c :
And we also change the old constants:
You should get the following image:
upd
After I worked in Ada, I found one mistake - where I used Mod . This function played my role as a crutch, since we have one color channel takes values from 0 to 255, it turned out that I did Mod 'so that the value was not higher than the maximum. As you can see in the images, this leads to sharp transitions in colors.
I will fix this with this formula:
. It turns out a variation of linear interpolation in a linear space.
To make it easier to use, we will write the following function:
And we replace the color assignment to dots:
It turns out the following image:
upd2
Since in Fortran, as in Ada, EOL is automatically selected, then under non-Unix systems, PNM is incorrect. To avoid this, you will have to use the following bike:
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
I will write in the Fortran 90 dialect, but with GNU extensions (the same double precision numbers).
A bit of fractal theory
Fractal ( fractus - crushed, broken, broken) - in the narrow sense, a complex geometric shape that has the property of self-similarity, that is, composed of several parts, each of which is similar to the entire figure as a whole.
')
There is a large group of fractals - algebraic fractals. This name derives from the principle of constructing such fractals. For their construction use recursive functions. For example, algebraic fractals are: Julia set, Mandelbrot set, Newton's pool (fractal).
Construction of algebraic fractals
One construction method is an iterative calculation.
where , but some function. The calculation of this function continues until a certain condition is met. When this condition is fulfilled, a dot is displayed on the screen.The behavior of the function for different points of the complex plane may have different behavior:
- Striving for infinity
- Aims to 0
- Accepts several fixed values and does not go beyond them.
- Chaotic behavior. No trends
One of the simplest (both for understanding and implementation) algorithms for constructing algebraic fractals is known as the “ escape time algorithm ”. In short, an iterative calculation of the number is made for each point of the complex plane.
It has been proven that if
it turns out more bailout , then the sequence . Comparison with bailout allows you to find points lying outside the set. For points outside the set, the sequence will not tend to infinity, so never reach bailout .I will look at two simple fractals - Julia set and Mandelbrot set. They are calculated by the same function, and differ only in a constant, where Julia has a constant, and in Mandelbrob the constant depends on the point of the complex plane.
Building Julia Set
The beginning and end of the program is simple and trivial:
program frac end Next we need to enter several constants:
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) ! And here you need to explain something. REAL and COMPLEX is a single-precision floating-point number and a complex number consisting of two REALs . REAL (8) is a double precision number, and COMPLEX (8) , respectively, is a complex number consisting of two REAL (8) . Also, the letter D is added to the standard functions (as is the case with ABS ), which indicates the use of double precision numbers.
Next we need to enter a few variables:
CHARACTER, DIMENSION(:), ALLOCATABLE :: matrix ! COMPLEX(8) :: z INTEGER :: x, y, iter REAL(8) :: iter2 ! ALLOCATE(matrix(0 : resox * resoy * 3)) Using ALLOCATABLE is a must! Since otherwise:
CHARACTER, DIMENSION(0:resox * resoy * 3) :: matrix ! The memory will be allocated on the stack, and this is not acceptable when used in multiple threads. Therefore, we allocate memory on the heap.
So I do not use two-dimensional arrays, because when selected on the heap, the array will take up more space, and it will be more difficult to write it to the file.
Next we need to specify the number of threads generated by OpenMP:
call omp_set_num_threads(16) And here the calculations begin. In Fortran, OpenMP directives are specified via comments (in C / C ++, there is a special #pragma macro for this).
!$OMP PARALLEL DO PRIVATE(iter, iter2, z) do x = 0, resox-1 do y = 0, resoy-1 iter = 0 ! z = DCMPLX(x * CXshag + xshift, y * CYshag + yshift) ! Z do while ((iter .le. iterations) .and. (CDABS(z) .le. 4.0_8)) ! bailout 2, 4, .. z = z**2 + c iter = iter + 1 end do iter2 = DFLOAT(iter) - DLOG(DLOG(CDABS(z))) / DLOG(2.0_8) ! iter-log2(ln(|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 The most time consuming part is over. Next, it remains for us to display information in a convenient form for perception - the image. I will use the PNM format, since This is the easiest image format.
PNM consists of several sections:
P6 # 1024 1024 255 The first line is an indication of the format of recording information about pixels:
- P1 / P4 - black and white image
- P2 / P5 - gray image
- P3 / P6 - color image
The first P is an option where the color of a pixel is recorded with an ASCII character, and the second P allows you to record the color of a pixel in binary form (which saves disk space).
The next line is a comment, then the image resolution and the number of colors per pixel. After the number of colors is directly information about the pixels.
Start writing the image to the file:
open(unit=8, status = 'replace', 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 we perform for decency, because when you exit the program, the OS will still return the occupied memory to the system.
You can use the GNU compiler compiler to build the program:
gfortran -std=gnu frac.f90 -o frac.run -fopenmp You should get the following image:
1024x1024
How to generate a Mandelbrot Set
It is easy to do, it is enough to replace c and change several constants. In this case, remove the constant c and add the variable c :
COMPLEX(8) :: z, c z = DCMPLX(x * CXshag + xshift, y * CYshag + yshift) ! Z c = z And we also change the old constants:
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 You should get the following image:
1023x682
upd
After I worked in Ada, I found one mistake - where I used Mod . This function played my role as a crutch, since we have one color channel takes values from 0 to 255, it turned out that I did Mod 'so that the value was not higher than the maximum. As you can see in the images, this leads to sharp transitions in colors.
I will fix this with this formula:
. It turns out a variation of linear interpolation in a linear space.To make it easier to use, we will write the following function:
function myround(a) result(ret) REAL(8), intent(in) :: a INTEGER :: ret ! b 255, INTEGER . ret = ABS(INT(a - 255.0_8 * NINT(a / 255.0_8))) !NINT() round() C, .. end function myround program frac INTEGER :: myround And we replace the color assignment to dots:
matrix((x + y * resox) * 3 + 0) = CHAR(myround(iter2 * 7.0_8)) matrix((x + y * resox) * 3 + 1) = CHAR(myround(iter2 * 14.0_8)) matrix((x + y * resox) * 3 + 2) = CHAR(myround(iter2 * 2.0_8)) It turns out the following image:
1024x1024
upd2
Since in Fortran, as in Ada, EOL is automatically selected, then under non-Unix systems, PNM is incorrect. To avoid this, you will have to use the following bike:
write(8, '(a, a, $)') "P6", char(10) write(8, '(a,$)') char(10) write(8, '(I0, a, I0, a, $)') resox, ' ', resoy, char(10) write(8, '(I0, a, $)') 255, char(10) write(8, *) matrix Used Books
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
Source: https://habr.com/ru/post/242221/
All Articles