Logo by kr0mat1k - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Did you know ? you can use LaTex inside Postings on fractalforums.com!
 
*
Welcome, Guest. Please login or register. September 26, 2022, 07:37:00 AM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: [1] 2   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: buddhabrot in fortran :)  (Read 5876 times)
0 Members and 1 Guest are viewing this topic.
ker2x
Fractal Molossus
**
Posts: 795


WWW
« on: March 07, 2010, 07:33:07 AM »

I decided to learn fortran yesterday smiley
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

Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #1 on: March 07, 2010, 08:02:47 AM »

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  grin )

« Last Edit: March 07, 2010, 08:06:34 AM by ker2x, Reason: 100 millions, not 10 millions » Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
Sockratease
Global Moderator
Fractal Senior
******
Posts: 3181



« Reply #2 on: March 07, 2010, 12:49:33 PM »

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

Nice image too...
Logged

Life is complex - It has real and imaginary components.

The All New Fractal Forums is now in Public Beta Testing! Visit FractalForums.org and check it out!
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #3 on: March 07, 2010, 06:33:52 PM »

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
Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #4 on: March 08, 2010, 03:27:55 AM »

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

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
Timeroot
Fractal Fertilizer
*****
Posts: 362


The pwnge.


WWW
« Reply #5 on: March 08, 2010, 06:15:18 AM »

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!  cheesy
Logged

Someday, man will understand primary theory; how every aspect of our universe has come about. Then we will describe all of physics, build a complete understanding of genetic engineering, catalog all planets, and find intelligent life. And then we'll just puzzle over fractals for eternity.
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #6 on: March 08, 2010, 07:38:40 AM »

Thx for the tips ! I'll try that smiley
(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)
Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
Nahee_Enterprises
World Renowned
Fractal Senior
******
Posts: 2250


use email to contact


nahee_enterprises Nahee.Enterprises NaheeEnterprise
WWW
« Reply #7 on: March 08, 2010, 07:42:24 PM »

I decided to learn Fortran yesterday  smiley
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??
Logged

Nahee_Enterprises
World Renowned
Fractal Senior
******
Posts: 2250


use email to contact


nahee_enterprises Nahee.Enterprises NaheeEnterprise
WWW
« Reply #8 on: March 08, 2010, 07:56:15 PM »

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  grin )

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??
 
Logged

ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #9 on: March 09, 2010, 07:59:11 AM »

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  grin )

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 smiley

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

Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #10 on: March 09, 2010, 06:54:29 PM »

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!  cheesy

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

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #11 on: March 09, 2010, 09:43:34 PM »

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
Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #12 on: March 10, 2010, 01:23:12 AM »

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

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
lkmitch
Fractal Lover
**
Posts: 238



« Reply #13 on: March 10, 2010, 04:32:46 PM »

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.  smiley  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.
Logged
ker2x
Fractal Molossus
**
Posts: 795


WWW
« Reply #14 on: March 10, 2010, 09:01:31 PM »

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)

Logged

often times... there are other approaches which are kinda crappy until you put them in the context of parallel machines
(en) http://www.blog-gpgpu.com/ , (fr) http://www.keru.org/ ,
Sysadmin & DBA @ http://www.over-blog.com/
Pages: [1] 2   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
3D Buddhabrot 3D Fractal Generation David Makin 1 2911 Last post January 12, 2008, 02:36:31 PM
by twinbee
Buddhabrot fractals Images Showcase (Rate My Fractal) « 1 2 ... 5 6 » Buddhi 75 9009 Last post May 26, 2010, 09:30:20 PM
by johandebock
Low iteration anti-buddhabrot Mandelbrot & Julia Set ker2x 1 1522 Last post December 30, 2009, 03:51:32 AM
by ker2x
A generic name for buddhabrot-like fractal ? General Discussion ker2x 3 1755 Last post January 02, 2010, 06:24:55 PM
by Melancholyman
a creative buddhabrot rendering Fractal Programs ker2x 14 4900 Last post May 19, 2010, 08:00:13 PM
by kram1032

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.172 seconds with 24 queries. (Pretty URLs adds 0.01s, 2q)