MODULE global IMPLICIT NONE INTEGER, DIMENSION(:,:), ALLOCATABLE :: iters ! global objects: REAL(8), DIMENSION(2) :: XYINC = 0.0 ! Distance between pixels REAL(8), DIMENSION(2) :: SWV = 0.0 ! Beginning coordinate: real and imaginery. INTEGER :: XYdots = 0 ! Number of display pixels in X, Y directions INTEGER :: Niter = 1000! Maximum number of iterations INTEGER :: Procs = 1 CONTAINS subroutine IntData2pgm(s1,s2,idata,name) ! Simple subroutine to dump integer data in a PGM format implicit none INTEGER, INTENT(IN) :: s1, s2 INTEGER, DIMENSION(s1,s2), intent(in) :: idata character(*), intent(in) :: name ! integer :: i, j, ouni integer :: vp, vs real :: rmin, rmax character(80) :: fname ! ! Write on unit 700 with PGM format ouni = 700 fname = TRIM(name)//".pgm" 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)") "P2" ! Dimensions write(ouni,*) s1, s2 ! Maximum value write(ouni,"(a)") "255" ! Values from 0 to 255 rmin = MINVAL(idata); rmax = MAXVAL(idata) vs = 0 do j = 1, s2 do i = 1, s1 vp = INT ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) ) vs = vs + 1 write(ouni,"(i4)",advance="no") vp 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 IntData2pgm subroutine IntData2ppm(s1,s2,idata,name) ! Simple subroutine to dump integer data in a PPM format implicit none INTEGER, INTENT(IN) :: s1, s2 INTEGER, DIMENSION(s1,s2), intent(in) :: idata 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 PPM 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" ! Values from 0 to 255 rmin = MINVAL(idata); rmax = MAXVAL(idata) vs = 0 do j = 1, s2 do i = 1, s1 rvp = ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) ) ! 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 IntData2ppm END MODULE global ! ! 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 = Square area side length ! INTEGER procs = Independent computing units ! INTEGER NITER = Maximum iterations ! ! ************************************************************************** PROGRAM Mandel USE global IMPLICIT NONE !!! MPI Fortran include file ! Mandel variables REAL(8) :: Sreal, Simag, Range = 0.0 ! Beginning coords and zoom factor REAL :: tempo2, tempo1, cpu2, cpu1, tempro, cpupro CHARACTER(40) :: fname INTEGER :: dX, dY, t, X0, X1, Y0, Y1 INTEGER :: st LOGICAL :: LetsGo ! Add MPI variables INTEGER :: my_rank, numprocs ! !!! MPI initialization: define my_rank and numprocs my_rank=0; numprocs=1 ! LetsGo = .TRUE. IF ( my_rank == 0 ) THEN ! Master process looks for input data CALL CPU_TIME(tempo1) fname = "mandel.inp" OPEN(11,FILE=fname,STATUS="OLD", IOSTAT=st) IF ( st /= 0 ) THEN WRITE(*,*) "Error access file ",TRIM(fname) LetsGo = .FALSE. ELSE READ(11,*,IOSTAT=st) 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" LetsGo = .FALSE. END IF END IF SWV(1) = Sreal; SWV(2) = Simag IF ( LetsGo ) THEN ! Master allocates entire matrix ALLOCATE(iters(0:XYdots-1,0:XYdots-1),STAT=st) IF ( st /= 0 ) THEN WRITE(*,*) "Error allocating ITERS(",XYdots,",",XYdots,")" LetsGo = .FALSE. END IF END IF END IF ! Send start flag IF ( .NOT.LetsGo ) STOP "Error exit" ! Abort program !!! Send data to processes: SWV, XYdots, Range, NITER WRITE(*,*) " Input data are: ", & & SWV(1), SWV(2), XYdots, Range, procs, NITER ! Every process executes procs = numprocs ! procs data required for compatibility reasons ! ---> Compute complex spacing between pixels <--- XYinc(1) = Range / real(XYdots) XYinc(2) = Range / real(XYdots) dY = ( XYdots + (Procs - 1) ) / Procs X0 = 0; X1 = XYdots Y0 = dY * my_rank Y1 = dY * ( my_rank + 1) if ( Y1 > XYdots ) Y1 = XYdots ! IF ( my_rank /= 0 ) THEN ALLOCATE(iters(0:XYdots-1,0:XYdots-1),STAT=st) END IF CALL calcola(X0, X1, Y0, Y1) !!!! Root process gathers results CALL CPU_TIME(tempo2) tempro = (tempo2 - tempo1) CALL IntData2pgm(XYdots,XYdots,iters,"mandel") CALL IntData2ppm(XYdots,XYdots,iters,"mandel") WRITE(*,*) "CPU time is ",tempro !!! Terminate MPI processes STOP END PROGRAM MANDEL SUBROUTINE calcola( X0, X1, Y0, Y1) USE global IMPLICIT NONE INTEGER, INTENT(IN) :: x0, x1, y0, y1 INTEGER :: pos INTEGER :: IX = 0 ! General purpose counter. INTEGER :: IY = 0 ! General purpose counter. INTEGER :: IZ = 0 ! General purpose counter. REAL(8) :: CA = 0.0 ! Current coordinate. Real. REAL(8) :: CB = 0.0 ! Current coordinate. Imaginary. REAL(8) :: ZA = 0.0 ! Current itteration. Real. REAL(8) :: ZB = 0.0 ! Current itteration. Imaginary. REAL(8) :: ZSIZE = 0.0 ! Real distance of Z(A,B) from 0,0 REAL(8) :: ZNEWA = 0.0 ! Next itteration. Real. REAL(8) :: ZNEWB = 0.0 ! Next itteration. Immaginary. DO IX = X0, X1-1 ! New column CA = XYinc(1) * IX + SWV(1) DO IY = Y0, Y1-1 ZA = 0.0 ZB = 0.0 CB = XYinc(2) * IY + SWV(2) 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 END SUBROUTINE calcola