buddhabrot in fortran :)

[ker2x]

I decided to learn fortran yesterday
it's fun and faaaaaaaaaaaaaaaaaaaaaaaaast !

Code:

USE omp_lib
use ifport
IMPLICIT NONE

character ( len = 255 ) :: filename = 'buddhabrot.ppm'


! mix/max grid coordinate
REAL, PARAMETER :: xmin = -2., xmax = 2., ymin = -2., ymax = 2.

! Arbitrary maximum # of iterations
INTEGER, PARAMETER :: n_max=100
INTEGER, PARAMETER :: grid_resolution = 1000
INTEGER*8, PARAMETER :: batchSize = 10000000

INTEGER :: exposureRMap(grid_resolution, grid_resolution)
INTEGER :: exposureGMap(grid_resolution, grid_resolution)
INTEGER :: exposureBMap(grid_resolution, grid_resolution)

REAL :: x,y
INTEGER*8 :: i
COMPLEX :: z, c
INTEGER :: iter
INTEGER :: tempX, tempY

integer ppm_i
integer ppm_j
integer ppm_jhi
integer ppm_jlo
integer file_out_unit

exposureRMap = 1
exposureGMap = 1
exposureBMap = 1

x = RANDOM(1)

DO i=1, batchSize
  x = RANDOM(0) * (xmax-xmin) + xmin
  y = RANDOM(0) * (ymax-ymin) + ymin
  z = CMPLX(x,y)
  IF (notInMSet(z, n_max) .EQV. .TRUE.) THEN
    c = z
    DO iter=1, n_max
      z = z**2 + c
      IF(ABS(z) < 2) THEN
        TempX = INT(grid_resolution * (REAL(z) + xmax) / (xmax - xmin))
        TempY = INT(grid_resolution * (AIMAG(z) + ymax) / (ymax - ymin))
        IF(iter < 10) THEN
            exposureRMap(TempX, TempY) = exposureRMap(TempX, TempY) + 1
          ELSE IF(iter < 50) THEN
            exposureBMap(TempX, TempY) = exposureBMap(TempX, TempY) + 1
          ELSE IF(iter < 100) THEN
            exposureGMap(TempX, TempY) = exposureGMap(TempX, TempY) + 1
        ENDIF
      END IF
    END DO
  END IF
END DO


exposureRMap = exposureRMap - (MINVAL(exposureRMap) )
exposureGMap = exposureGMap - (MINVAL(exposureGMap) )
exposureBMap = exposureBMap - (MINVAL(exposureBMap) )
exposureRMap = exposureRMap / (MAXVAL(exposureRMap)/255. )
exposureGMap = exposureGMap / (MAXVAL(exposureGMap)/128. )
exposureBMap = exposureBMap / (MAXVAL(exposureBMap)/255. )
!write(*,*) MAXVAL(exposureRMap) , MAXVAL(exposureGMap) , MAXVAL(exposureBMap)
!write(*,*) MINVAL(exposureRMap) , MINVAL(exposureGMap) , MINVAL(exposureBMap)

open ( unit = file_out_unit, file = filename, status = 'replace', &
       form = 'formatted', access = 'sequential')
write ( file_out_unit, '(a2)' ) 'P3'
write ( file_out_unit, '(i5,2x,i5)' ) grid_resolution, grid_resolution
write ( file_out_unit, '(i5)' ) 255

do ppm_i = 1, grid_resolution
do ppm_jlo = 1, grid_resolution, 4
    ppm_jhi = min ( ppm_jlo + 3, grid_resolution )
    write ( file_out_unit, '(12i5)' ) ( exposureRMap(ppm_i,ppm_j), exposureGMap(ppm_i,ppm_j),exposureBMap(ppm_i,ppm_j), ppm_j = ppm_jlo,ppm_jhi )
  end do
end do


close ( unit = file_out_unit )

CONTAINS

        PURE FUNCTION notInMset(c, n_max)
          COMPLEX, INTENT(IN) :: c
          INTEGER, INTENT(IN) :: n_max
          INTEGER             :: n
          COMPLEX             :: z
          LOGICAL             :: notInMSet
          z = c
          n = 0

          DO WHILE (ABS(z) < 4.0 .AND. (n < n_max))
            z = z* z + c
            n = n + 1
          END DO

          IF (n >= n_max) THEN
            notInMset = .FALSE.
          ELSE
            notInMset = .TRUE.
          END IF
        END FUNCTION notInMset


END

[ker2x]

actually, it compute 2*100 millions mandelbrot with bailout = 100 in 1mn47s on a single-core in a virtual machine (ubuntu running on win7 using virtualbox)

I need to refactor the code so the compiler can do auto-parallelisation. (yup, that's a nice feature of intel fortran compiler, you can let the compiler do all the multithreading stuff if you want   )

[Sockratease]

I learned to program punching cards in Fortran III - glad to see it's still a toy for some folks!

Nice image too...

[ker2x]

I learned to program punching cards in Fortran III - glad to see it's still a toy for some folks!

and a tools for some others.
Some (most?) Top500 supercomputer are running fortran code.
One of the guy that helped me to learn fortran use fortran on a cluster that is 13th on top500 Supercomputer (~26.000 cores and 80TB of memory ...)
http://www.top500.org/system/10147

[ker2x]

The latest dev version is available here : http://github.com/ker2x/Fortran_buddhabrot

[Timeroot]

Looks cool - I didn't know Fortran handled complex numbers so natively. I think you could get a pretty big speed increase if, instead of checking inside/outside and calculating the orbit separately. Just keep a different counter that you update while finding the orbit the first time, and if it stays inside, then add it to the running total array.

Also, you can save a load of time but immediately removing big chunks, like the main cardioid and period-2 bulb, which can be found very simply. There are some formulas on the Wikipedia page about the MSet that tell how to eliminate them immediately. Good luck! 

[ker2x]

Thx for the tips ! I'll try that
(i'm bad at math, but it's a good way to learn)

fortran have many awesome feature. complex is one of them.
The way it handle array is really awesome too.

a good exemple is :
exposureRMap = (exposureRMap / REAL(maxRExposure))*intensityR

in pseudo-code, it could be something like :
foreach element  in exposureRMap
  exposureRMap[element] = (exposureRMap[element] / REAL(maxRExposure))*intensityR
end

And best of all :
remark: LOOP WAS AUTO-PARALLELIZED.

yup, the compiler transform this single line into some kind of "forearch loop" and handle the multithreading automagically.
Thoses niceties explain why fortran is still used on Top500 supercomputer.
With the right parameters, it can launch the threads on remote computers of a cluster, using openmp.

Sadly, most of my code isn't parallized because of "ANTI" and "FLOW" dependencies.
I'm planning to work on it to parallelize the main loop (the problem is in the exposureMap arrays)

[Nahee_Enterprises]

I decided to learn Fortran yesterday 
it's fun and faaaaaaaaaaaaaaaaaaaaaaaaast !

Way to go!!  It has been at least two decades (probably more) since I last worked with Fortran, and that was on an IBM mainframe.

Are you using the Evaluation copy, or did you fork over the $899 purchase price, or "acquire it" in some other manner??

[Nahee_Enterprises]

actually, it on a single-core in a virtual machine (ubuntu running on win7 using virtualbox)
I need to refactor the code so the compiler can do auto-parallelisation.
(yup, that's a nice feature of Intel Fortran compiler, you can let the compiler
do all the multithreading stuff if you want   )

I noticed that the Intel® Visual Fortran Compiler for Windows Professional Edition also has the option of including the Visual Numerics IMSL Fortran Library.  I am guessing that this might not be part of the Linux version??
    http://software.intel.com/en-us/articles/visual-numerics-imsl-fortran-library/

If you already have a Windows environment, why not just use the Visual version??
 

[ker2x]

actually, it on a single-core in a virtual machine (ubuntu running on win7 using virtualbox)
I need to refactor the code so the compiler can do auto-parallelisation.
(yup, that's a nice feature of Intel Fortran compiler, you can let the compiler
do all the multithreading stuff if you want   )

I noticed that the Intel® Visual Fortran Compiler for Windows Professional Edition also has the option of including the Visual Numerics IMSL Fortran Library.  I am guessing that this might not be part of the Linux version??
    http://software.intel.com/en-us/articles/visual-numerics-imsl-fortran-library/

If you already have a Windows environment, why not just use the Visual version??
 

I use the "free for non commercial use" intel software, that are available only on linux

http://software.intel.com/en-us/articles/non-commercial-software-download/

Intel® C++ Compiler Professional Edition for Linux*
Intel® Fortran Compiler Professional Edition for Linux
Intel® Compiler Suite Professional Edition for Linux
Intel® Math Kernel Library (Intel® MKL) for Linux
Intel® Integrated Performance Primitives (Intel® IPP) for Linux
Intel® VTune™ Performance Analyzer for Linux
Intel® Thread Checker for Linux

[ker2x]

Looks cool - I didn't know Fortran handled complex numbers so natively. I think you could get a pretty big speed increase if, instead of checking inside/outside and calculating the orbit separately. Just keep a different counter that you update while finding the orbit the first time, and if it stays inside, then add it to the running total array.

Also, you can save a load of time but immediately removing big chunks, like the main cardioid and period-2 bulb, which can be found very simply. There are some formulas on the Wikipedia page about the MSet that tell how to eliminate them immediately. Good luck! 

I found an itneresting website about the mandelbrot discs : http://classes.yale.edu/fractals/MandelSet/MandelCombinatorics/N2Scaling/N2Scaling.html

[ker2x]

i added the main cardiod and 2nd circle optimisation.
code pushed to github as usual.

the main change is :

Code:

          IF(((ABS(z - CMPLX(-1,0) )) < 0.25) .OR. (( ABS( 1.0 - SQRT(1-(4*z))  ) < 1.0 ) )  ) THEN
            notInMset = .FALSE.
          ELSE

[ker2x]

updated to run in multithread.
I still have a weird bug...
it works with a grid_size of 500x500 and segfault with 1000x1000

[lkmitch]

I found an itneresting website about the mandelbrot discs : http://classes.yale.edu/fractals/MandelSet/MandelCombinatorics/N2Scaling/N2Scaling.html

Shameless plug--I'm the Kerry Mitchell mentioned on that page.    Michael Frame (who created those class notes) told me about the scaling rule some 20 years ago.  While investigating it numerically (in FORTRAN), I discovered those deviation features.

[ker2x]

i'm wondering if it worth it to pre-compute the 3rd and 4th order disc.

i know they are not real disc, and it certainly doesn't worth it to compute a very good estimation, but i can estimate the size of a rough interal disc contained in the 3rd and 4th order disc/shape. It's very cheap in cpu time. (much cheaper than computing the main cardioid, but considering its size : it worth it)