MODULE globali IMPLICIT NONE ! Global objects: INTEGER, DIMENSION(:,:), ALLOCATABLE :: iters, pixel ! results matrices INTEGER :: XYdots = 0 ! Number of pixels in X, Y directions INTEGER :: NITER = 1000 ! Maximum number of iterations INTEGER :: Procs = 1 ! Activated threads; actually not used CONTAINS subroutine IntData2pixel(s1,s2,idata,pixel) ! Simple subroutine to normalize integer data sets to [0-255] values ! Can be executed multi-threaded implicit none INTEGER, INTENT(IN) :: s1, s2 INTEGER, DIMENSION(s1,s2), intent(in) :: idata INTEGER, DIMENSION(s1,s2), intent(out) :: pixel ! integer :: i, j, dj, j0, j1 real :: rmin, rmax ! Values from 0 to 255 rmin = MINVAL(idata); rmax = MAXVAL(idata) !$omp parallel do default(none) private(i,j) shared(s1,s2,rmin,rmax,idata,pixel) do j = 1, s2 do i = 1, s1 pixel(i,j) = INT ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) ) end do end do !$omp end parallel do return end subroutine IntData2pixel subroutine IntPixel2ppm(s1,s2,pixel,name) ! Simple subroutine to dump integer data in range [0-255] to PPM image format implicit none INTEGER, INTENT(IN) :: s1, s2 INTEGER, DIMENSION(s1,s2), intent(in) :: pixel character(*), intent(in) :: name ! integer :: i, j, ouni integer :: vp, vs, r, g, b real :: rmin, rmax, rvp character(80) :: fname ! ! Write on unit 700 with PGM format ouni = 700 fname = TRIM(name)//".ppm" open(ouni,file=trim(fname),status="replace",iostat=vs) if ( vs /= 0 ) then print*,"!!!! Error write access to file ",trim(fname) end if ! Magic code write(ouni,"(a)") "P3" ! Dimensions write(ouni,*) s1, s2 ! Maximum value write(ouni,"(a)") "255" vs = 0 do i = 1, s1 do j = 1, s2 rvp = pixel(i,j) ! Palette: r,g,b values vp = MOD(INT(rvp),256) SELECT CASE (vp) CASE(1:39) vp = INT ( (rvp - 1.0) * 255.0 / 40.0 ) r = vp; g = 0; b = 0 CASE(40:79) vp = INT ( (rvp - 40.0) * 255.0 / 40.0 ) r = vp; g = vp; b = 0 CASE(80:119) vp = INT ( (rvp - 80.0) * 255.0 / 40.0 ) r = 0; g = vp; b = 0 CASE(120:159) vp = INT ( (rvp - 120.0) * 255.0 / 40.0 ) r = 0; g = vp; b = vp CASE(160:199) vp = INT ( (rvp - 160.0) * 255.0 / 40.0 ) r = 0; g = 0; b = vp CASE(200:239) vp = INT ( (rvp - 200.0) * 255.0 / 40.0 ) r = vp; g = 0; b = vp CASE(240:255) r = 0; g = 0; b = (rvp/10.0) CASE DEFAULT ! 0 value r = 0; g = 0; b = 0 END SELECT vs = vs + 1 write(ouni,"(3(i4,1x))",advance="no") r, g, b if (vs >= 10 ) then write(ouni,"(a)") " " vs = 0 end if end do write(ouni,"(a)") " " vs = 0 end do close(ouni) return end subroutine IntPixel2ppm SUBROUTINE calcola(Sreal,Simag,Range) IMPLICIT NONE REAL(8), INTENT(IN) :: Sreal,Simag,Range INTEGER :: ix, iy, iz REAL(8) :: ca, cb, za, zb REAL(8) ::zsize, znewa, znewb REAL(8) :: Xinc, Yinc INTEGER :: my_id, Procs, omp_get_thread_num, omp_get_num_threads Xinc = Range / real(XYdots) Yinc = Range / real(XYdots) ! Every thread executes ! Insert parallel region, pay attention to shared or private variables ! Do it in parallel ay attention to private variables DO iy = 0, XYdots-1 DO ix = 0, XYdots-1 ca = Xinc * ix + Sreal za = 0.0 zb = 0.0 cb = Yinc * iy + Simag zsize = (ca*ca + cb*cb ) ! iteration loop for each ( z**2 + c ) DO iz = 0, NITER IF ( zsize >= 4 ) EXIT znewa = ( za*za ) - ( zb*zb ) + ca znewb = ( 2*za*zb ) + cb zsize = ( znewa*znewa + znewb*znewb ) za = znewa zb = znewb END DO ! results iters(ix,iy) = iz END DO END DO RETURN END SUBROUTINE calcola END MODULE globali ! ! MAND - Compute the points outside the Mandelbrot set ! ! Computes the iteration number for each pixel in a given squared area ! ! Syntax: ! mandel Sreal Simag XYdots Range procs NITER ! ! INPUT: ! REAL(8) Sreal = Real coord. of lower left vertex ! REAL(8) Simag = Imaginary coord. of lower left vertex ! INTEGER XYdots = Horizontal and vertical pixels ! INTEGER Range = Squared area side length ! INTEGER Procs = Independent computing units ! INTEGER NITER = Maximum iterations ! ! ************************************************************************** PROGRAM Mandel USE globali IMPLICIT NONE REAL(8) :: Sreal = 0.0 ! Beginning coordinate. Real. REAL(8) :: Simag = 0.0 ! Beginning coordinate. Imaginary. REAL(8) :: Range = 0.0 ! Zoom factor REAL :: tempo2, tempo1 REAL(8) :: wall1, wall2 CHARACTER(40) :: fname INTEGER :: dX, dY, t, X0, X1, Y0, Y1 INTEGER :: st REAL(8) :: OMP_GET_WTIME fname = "mandel.inp" OPEN(11,FILE=fname,STATUS="OLD", IOSTAT=st) IF ( st /= 0 ) THEN WRITE(*,*) "Error access file ",TRIM(fname) STOP END IF READ(11,*,IOSTAT=st) Sreal, Simag, XYdots, Range, procs, NITER WRITE(*,*) "Input data are: ", Sreal, Simag, XYdots, Range, procs, NITER IF ( st /= 0 ) THEN WRITE(*,*) "Error reading input data. Examples are:" WRITE(*,*) "-2.000, -2.000, 1000, 4.000, 1, 1000" WRITE(*,*) "-1.500, -0.500, 1000, 0.500, 1, 1000" WRITE(*,*) "-1.250, -0.200, 1000, 0.100, 1, 1000" WRITE(*,*) "-1.250, -0.180, 1000, 0.030, 1, 1000" WRITE(*,*) "-1.235, -0.171, 1000, 0.006, 1, 1000" WRITE(*,*) "-1.233, -0.169, 1000, 0.001, 1, 1000" WRITE(*,*) "-1.233, -0.169, 1000, 0.0001, 1, 1000" STOP END IF ALLOCATE(iters(0:XYdots-1,0:XYdots-1),pixel(0:XYdots-1,0:XYdots-1),STAT=st) IF ( st /= 0 ) THEN WRITE(*,*) "Error allocating ITERS(",XYdots,",",XYdots,")" STOP ELSE WRITE(*,*) "Allocated ITERS(",XYdots,",",XYdots,")" END IF CALL CPU_TIME(tempo1) !This statement is executed only by master thread if compiler openmp flag is used !wall1 = OMP_GET_WTIME() CALL calcola(Sreal, Simag, Range) CALL CPU_TIME(tempo2) !This statement is executed only by master thread if compiler openmp flag is used !wall2 = OMP_GET_WTIME() CALL IntData2pixel(XYdots,XYdots,iters,pixel) CALL IntPixel2ppm(XYdots,XYdots,pixel,"mandel") PRINT*,"Global cpu, elapsed times = ",(tempo2 - tempo1), (wall2-wall1) STOP END PROGRAM MANDEL