UCL WIKI

UCL Logo
Child pages
  • 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