This section describes the development of the Mandelbrot example program, and explains many of the Fortran 90 and HPF features used in this program.
At this point, you may want to compile and run the Mandelbrot exam-
ple program provided with your PSE software. Simple instructions can
be found in the file /usr/examples/hpf/README.mandelbrot .
Figure 4-1 shows the visualization that
is displayed when the example program first starts. A window is
displayed showing the Mandelbrot set in black, surrounded by the
Mandelbrot complement set shown in multiple colors representing
various ranges of potential according to the electrostatic model
explained in Section 4.1.3.
The window and image are sized so that the center of the window
represents the origin of the complex plane. The axes of the plane
intersect the edges of the window at a distance of 2 from the ori-
gin. The point representing -2, the point in the Mandelbrot set most
distant from the origin, is located in the middle of the left side
of the window. The size of the display area is 625 x 625.
The mouse buttons are used to zoom in or zoom out, creating a new
image of different scale with each click, as explained in the file
/usr/examples/hpf/README.mandelbrot .
Figure 4-2 is a zoomed-in view of a
particularly complex area of the border between the Mandelbrot set
and its complement:
Example 4-1 shows the iteration of the function that determines whether a given complex number is in the Mandelbrot set.
COMPLEX :: z, c
INTEGER :: n, esc_time=0
INTEGER, PARAMETER :: n_max=1000 ! Arbitrary maximum # of iterations
INTEGER, PARAMETER :: escape_radius=400 ! Arbitrary criterion for escape
LOGICAL :: in_the_mandel_set
z=0
n = 0
DO WHILE (ABS(z) < escape_radius .AND. (n <= n_max) )
z = z**2 + c
n = n + 1
END DO
esc_time = n
IF (n <= n_max) THEN
in_the_mandel_set = .TRUE.
ELSE
in_the_mandel_set = .FALSE.
END IF
Some of the Fortran 90 features used in Example 4-1 are:
Example 4-1 tests whether any
given value of c is in the Mandelbrot set. The loop
condition uses the ABS intrinsic, because complex numbers can only
be compared by absolute value, which is defined as the distance from
the origin in the complex plane. If n_max iterations
are performed without the absolute value of z exceeding
the escape radius, the given value of c is presumed
(although not proven) to be part of the Mandelbrot set. If the
the loop is exited before n_max iterations have been
completed, the given value of c has been proven to lie
outside of the Mandelbrot set.
For points proven to be outside of the set, the value of n when the loop is exited is the escape time. The escape time can be used to plot equipotential lines and to color in regions of varying potential.
Any value greater than or equal to 2 could have been chosen for the escape radius, because the Mandelbrot set is entirely contained within a circle of radius 2. However, a substantially larger value (400) was chosen because it will cause the equipotential lines in the Mandelbrot complement set to be considerably more accurate.
Even though it makes the complement set more accurate, using a larger escape radius causes a very slight degradation in the accuracy of the shape of the Mandelbrot set itself. However, the effect of this degradation is barely noticable in visual terms, because values tend to escalate very rapidly once their absolute value exceeds 2.
The image of the Mandelbrot set is plotted on a grid, with each pixel of a window of a computer monitor representing one point on the grid. The escape time is calculated for each point proven to lie outside of the Mandelbrot set, with all points having the same escape time assigned the same color in the image. Points not proven to lie outside of the set are left black. Example 4-1, which calculates the escape time for a single point, can be expanded to generate the entire grid simply by putting nested DO loops around the calculation. See Example 4-2:
COMPLEX :: z, c
INTEGER :: n, esc_time=0, target(grid_height, grid_width)
INTEGER, PARAMETER :: n_max=1000 ! Arbitrary maximum # of iterations
INTEGER, PARAMETER :: escape_radius=400 ! Arbitrary criterion for escape
INTEGER, PARAMETER :: grid_height=625, grid_width=625
DO x = 1, grid_width
DO y = 1, grid_height
c = CMPLX(x, y)
z=0
n = 0
DO WHILE (ABS(z) < escape_radius .AND. (n <= n_max) )
z = z**2 + c
n = n + 1
END DO
esc_time = n
target(x, y) = esc_time
END DO
END DO
As a simplification, Example 4-2 assumes the origin of the complex plane is in the lower left-hand corner of the image.
DO loops prescribe that calculations be done in a certain order. Therefore, Example 4-2 prescribes the order in which the grid points are calculated. However, careful examination of Example 4-2 reveals that the computation for each grid point is completely independent and unrelated to the computation for any other point on the grid. Thus, the order of the calculation has no effect on the result of the program. The same result would be produced if the grid points were calculated in the opposite order, or even in random order. This means that this routine is an excellent candidate for parallelizing with HPF. When the routine is converted to HPF, several grid points will be calculated simultaneously, depending upon the number of processors available. Generating Mandelbrot visualizations is a completely (or "embarrassingly") parallel computation.
To allow parallel execution of this routine, the target array must be distributed across processors using the DISTRIBUTE directive, and the two outer DO loops must either be replaced with a FORALL structure, or marked with the INDEPENDENT directive.
Replacing DO loops with a FORALL structure presents a problem, however: FORALL is not a loop, but an assignment statement. An assignment statement cannot contain the assignments to multiple variables and flow control constructs (such as DO WHILE) that occur in Example 4-2. A FORALL structure is limited to assigning values to elements of a single array.
The solution to this problem is to package the bulk of the
calculation into a user-defined function. Function calls inside
assignment statements are permitted, and in this way the entire
routine can be parallelized. Example 4-3 shows the FORALL structure containing a call
to the user-defined function escape_time , and Example 4-4 shows the function, which
contains the calculation for a single grid point.
INTEGER :: target(grid_height, grid_width)
INTEGER, PARAMETER :: n_max=1000 ! Arbitrary maximum # of iterations
INTEGER, PARAMETER :: grid_height=625, grid_width=625
FORALL(x=1:grid_width, y=1:grid_height)
target(x, y) = escape_time( CMPLX(x, y), n_max )
END FORALL
PURE FUNCTION escape_time(c, n_max)
COMPLEX, INTENT(IN) :: c
INTEGER, INTENT(IN) :: nmax
INTEGER :: n
COMPLEX :: z
n = 0
z = c
DO WHILE (ABS(z) < 2.0 .AND. (n < n_max))
z = z * z + c
n = n + 1
END DO
IF (n >= n_max) THEN
escape_time = nmax
ELSE
escape_time = n
END IF
END FUNCTION escape_time
The escape_time function is given the PURE attribute.
The PURE attribute is an assertion that a function has no side
effects and makes no reference to mapped variables other than its
actual argument. A PURE function's only effect on the state of the
program is to return a value.
User-defined functions may be called inside a FORALL structure only if they are PURE functions. The reason for this rule is that iterations of a FORALL structure occur in an indeterminate order. Therefore, allowing functions that have side effects (such as modifying the value of a global variable) to be called from within a FORALL structure could lead to indeterminate program results.
More precise discussion of PURE and side effects can be found in Section 6.4.5 of this manual, and in the High Performance Fortran Language Specification.
Source code of the program on which this chapter is based can be
found in the file /usr/examples/hpf/mandelbrot.f90 .
