# UCL WIKI

##### Child pages
• Exercises - Fractals
Skip to end of banner
Go to start of banner

# Exercises - Fractals

Skip to end of metadata
Go to start of metadata

## Exercise - “Fractal Mandelbrot”

The following program produces an image of the Mandelbrot fractal.

```program mandelbrot

implicit none
integer  , parameter :: i_max    =  800
integer  , parameter :: j_max    =  600
integer  , parameter :: n_max    =  1000
real (8), parameter :: x_centre = -0.5d0
real (8), parameter :: y_centre =  0.0d0
real (8), parameter :: width    =  4.0d0
real (8), parameter :: height   =  3.0d0
real (8), parameter :: dx_di    =   width / i_max
real (8), parameter :: dy_dj    = -height / j_max
real (8), parameter :: x_offset = x_centre - 0.5d0 * (i_max + 1) * dx_di
real (8), parameter :: y_offset = y_centre - 0.5d0 * (j_max + 1) * dy_dj
integer, dimension (i_max, j_max) :: image
integer   :: i
integer   :: j
integer   :: n
real (8) :: x
real (8) :: y
real (8) :: x_0
real (8) :: y_0
real (8) :: x_sqr
real (8) :: y_sqr

do j = 1, j_max
y_0 = y_offset + dy_dj * j
do i = 1, i_max
x_0 = x_offset + dx_di * i
x = 0.0d0
y = 0.0d0
n = 0
do
x_sqr = x ** 2
y_sqr = y ** 2
if (x_sqr + y_sqr > 4.0) then
image (i, j) = (x_sqr + y_sqr)/4.0*255.0
exit
end if
if (n == n_max) then
image (i, j) = 0
exit
end if
y = y_0 + 2.0 * x * y
x = x_0 + x_sqr - y_sqr
n = n + 1
end do
end do
end do
open  (10, file = 'out.pgm')
write (10, '(a/ i0, 1x, i0/ i0)') 'P2', i_max, j_max, 255
write (10, '(i0)') image
close (10)

end program mandelbrot

```
1. Using a Code::Blockes editor type the above program into a file called mandelbrot.f90, compile and run. You can visualize the image 'out.pgm' using gimp.
1. Parallelize the loop by adding the \$omp .. statement, see Lecture slides and openmp.org. The env. parameter OMP_NUM_THREADS should be set to the number of processors by, e.g.

`OMP_NUM_THREADS=4`
2. During running time check top processors

`top`
3. Record time and compare with the serial code. How does it scale?

```program mandelbrot

implicit none
integer  , parameter :: i_max    =  8000
integer  , parameter :: j_max    =  6000
integer  , parameter :: n_max    =  1000
real (8), parameter :: x_centre = -0.5d0
real (8), parameter :: y_centre =  0.0d0
real (8), parameter :: width    =  4.0d0
real (8), parameter :: height   =  3.0d0
real (8), parameter :: dx_di    =   width / i_max
real (8), parameter :: dy_dj    = -height / j_max
real (8), parameter :: x_offset = x_centre - 0.5d0 * (i_max + 1) * dx_di
real (8), parameter :: y_offset = y_centre - 0.5d0 * (j_max + 1) * dy_dj
integer, dimension (i_max, j_max) :: image
integer   :: i
integer   :: j
integer   :: n
real (8) :: x
real (8) :: y
real (8) :: x_0
real (8) :: y_0
real (8) :: x_sqr
real (8) :: y_sqr

!\$omp parallel do private(j,y_0,i,x_0,x,y,n,x_sqr,y_sqr) shared(image)
do j = 1, j_max
y_0 = y_offset + dy_dj * j
do i = 1, i_max
x_0 = x_offset + dx_di * i
x = 0.0d0
y = 0.0d0
n = 0
do
x_sqr = x ** 2
y_sqr = y ** 2
if (x_sqr + y_sqr > 4.0) then
image (i, j) = (x_sqr + y_sqr)/4.0*255.0
exit
end if
if (n == n_max) then
image (i, j) = 0
exit
end if
y = y_0 + 2.0 * x * y
x = x_0 + x_sqr - y_sqr
n = n + 1
end do
end do
end do
!\$omp end parallel do
!
open  (10, file = 'out.pgm')
write (10, '(a/ i0, 1x, i0/ i0)') 'P2', i_max, j_max, 255
write (10, '(i0)') image
close (10)

end program mandelbrot

```

## Exercise - “Fractal Julia”

the following program to include the openmp parallelization. You will need an input file fractal.input

which can be used in Code::Blocks by adding the command < fractal.input in the Program Arguments section.

Record and compare execution times using the serial and OpenMP versions.

```!this is the FORTRAN 90 version of the program that puts a pretty fractal
!image into julia.ppm
!The contolling parameters are input in a namelist file, run juliaf with:
!juliaf < fractal.input
!Note that the use of argv and system is not part of standard fortran
!The namelist input is a very useful aspect of fortran
program julia
!the following means mrl=8, or double precision.  kind(1.0e0) is single
integer, parameter :: mrl=kind(1.0d0)
integer itermax,magnify,hxres,hyres
real(mrl) :: breakout,cr,ci,x0,y0
namelist/params/itermax,breakout,magnify,cr,ci,x0,y0,&
hxres,hyres,fname !note the use of & to continue this statement
real(mrl) x,xx,y,xl,yl,zsq,zm,rb
integer iter,hx,hy,info
integer red,green,blue
integer,allocatable :: RGB(:,:,:)
character fname*50
fname="fromfort.ppm" !default filename
read (*,nml=params) ! the * means read from standard input
!next write the header to the output file
open (unit=20,file=trim(fname),status="unknown",recl=3*hxres*hyres+10)
!open (unit=20,file="julia.ppm",status="unknown") !for P3
write(20,'(a)') "P6"  !edit to P3 if using P3
write(20,'("# magnify=",i6,"  itermax=",i6)') magnify,itermax !PPM comment line
write(20,*) hxres,hyres ! the * means let the compiler figure out the format
write(20,'(a)') "255"
!
allocate(RGB(hxres,hyres,3),stat=info)
!
if (info/=0) stop 'Allocation problem, RGB'
!
rb=sqrt(breakout)
do hy=1,hyres
do hx=1,hxres
y = (4.*hyres/hxres)*((hyres+1-hy-.5)/hyres-0.5)/magnify+y0
x = 4*((hx-.5)/hxres-0.5)/magnify+x0
zm=0
do iter=1,itermax-1
xl=x
yl=y
xx = x*x-y*y+cr
y = 2.0*x*y+ci
x = xx
zsq=x*x+y*y
if (zsq > zm) zm=zsq
if (zsq>breakout) exit
end do
if (iter>=itermax) then
red=0
green=255.*zm/breakout
blue=255.*zm/breakout
else
red=255*(rb+xl)/(2*rb)
green=0
blue=.5*255*(rb+yl)/(2*rb)
end if
RGB(hx,hy,1) = red
RGB(hx,hy,2) = green
RGB(hx,hy,3) = blue
!
continue
!
end do
end do
!
do hy=1,hyres
do hx=1,hxres
!
write(20,'(3a)',advance='no') char(RGB(hx,hy,1)),char(RGB(hx,hy,2)),char(RGB(hx,hy,3))
!
continue
!
end do
end do
!
deallocate(RGB)
close(20)
!
print *,  "wrote ",hxres,"by",hyres, "ppm file: ",trim(fname)
end program julia

```

Visualize julia.ppm   using any graphical software available.

• No labels