Let me introduce myself, my name is Dr Keith Evans and I’m a Physicist and Research Software Engineer. I recently started in the BEAR Research Software group as HPC software support for BlueBEAR and HPC Midlands+. On my first day, my boss gave me an objective, to run some code on the HPCs. This meant having to set up my accounts, request access; all the things a new starter needs to do. But there was another problem, I didn’t have any code to run. I needed to write some and since I like maths, pretty pictures and non-trivial parallelisation, the answer was obviously the Mandelbrot set.

For the non-mathematicians, the Mandelbrot set is the set of complex numbers (** c**) for which the function is not divergent.

This is determined iteratively, so let’s take the example where c = 1 + 0i (c = 1):

So, the complex number ** 1 + 0i** is divergent and therefore not a member of the Mandelbrot set. To test for divergence we set a maximum magnitude of the complex number

**; if**

*z***exceeds 2 then it is considered to have diverged. To complete the algorithm for determining the Mandelbrot set, we set a maximum number of iterations at which the function is considered bound. Therefore the algorithm for a determining whether a complex number belongs to the Mandlebrot set is as follows:**

*z**do while ( abs( z_complex ) .le. 2.0_rkind .and. iterations .lt. maximum_iterations ) *

* z_complex = z_complex**2 + c_complex*

* iterations = iterations + 1_ikind*

*end do*

Obviously, my algorithm is written in glorious Fortran as is 70% of all scientific software.

This algorithm returns either the number of iterations it took for the function to exceed a magnitude of 2 **[ abs( z_complex ) ]** or the

**for a given complex number**

*maximum_iterations***. Now we can map out complex number space with orthogonal axes of the real and imaginary parts of the complex number, which we will denote as**

*c***and**

*a***, therefore the complex number**

*b***. This will generate a 2-dimensional dataset of constant**

*c = a + bi***or**

*a***from which we can plot an image using the above defined**

*b***variable.**

*iterations*We do this simply in Fortran, as multidimensional arrays are standard, but we need to define limits in the real and imaginary dimension which we will call ** real_lower**,

**to represent the lower limits and**

*imaginary_lower***,**

*real_upper***to represent the upper limits as well as image dimensions as**

*imaginary_upper***image_width**and

**image_height**.

*do ind = 1, image_width*

* do jnd = 1,image_height*

* c_complex = cmplx( ( ind / image_width ) * ( real_upper – real_lower ) , &*

* ( jnd / image_height ) * ( imaginary_upper – imaginary_lower ) )*

* do while ( abs( z_complex ) .le. 2.0_rkind .and. iterations .lt. maximum_iterations ) *

* z_complex = z_complex**2 + c_complex*

* iterations = iterations + 1_ikind*

* end do*

* mandelbrot_image( jnd , ind ) = iterations*

* end do*

*end do*

With that we have an algorithm and in fact, code to calculate the Mandelbrot set for a given range in complex space. The Mandelbrot set manifests in an attractive multi-coloured fractal pattern, this is because the algorithm above is not a yes/no, it instead returns the number of iterations it took to determine if the complex was bound or not. Strongly divergent complex values will exit quickly and return a low number (red), weakly divergent complex numbers will return a high value (blue). This results in a brightly and varied coloured image and what’s more interesting – infinitely complex – zoom in further and more detail appears, further again and more detail still and so on. So, this presents a computational problem, we want to see as much detail as possible, but double the resolution and quadruple the number of pixels. So, what happens when we want 10^{5} pixels along each dimension? That’s 10^{10 }pixels, if each pixel required a single calculation then 10^{10} operations, but they don’t, they require between 1 and ** maximum_iterations** operations, which I set at 255.

That’s a lot of operations for a single core, in my next blog post I’ll explain how we can we parallelise such a problem to speed it up.