source code
! The Computer Language Benchmarks Game
! http://benchmarksgame.alioth.debian.org/
!
! Simon Geard, 6/1/05
! Modified by Waldemar Spomer, 10/1/09: openmp, one I/O
program mandelbrot
use omp_lib
implicit none
integer N, x, y, bit_num, i, incr, width_bytes
integer(kind=1) byte_acc, state, mask, res, maskbita, maskbitb
integer(kind=1), parameter :: K0 = 0, K1 = 1
integer, parameter :: iter = 50
real*8, parameter :: limit2 = 4.0d0
double precision :: absZ1, absZ2, invert
character(len=8) argv
complex(kind=8) Z1, Z2, C1, C2
logical in_mandelbrot
! Modified to use pointers
character(len=1), pointer, dimension(:) :: whole_data, pdata
double precision, dimension(1:2) :: m1, m2 ! for mask
nullify(pdata, whole_data)
call getarg(1,argv)
read(argv,*) N
allocate(whole_data(N**2/8),STAT=state)
! Output pbm header
write(*,'(a)') 'P4'
write(*,'(i0,a,i0)') N,' ',N
width_bytes = ishft(N,-3)
invert = 2.0d0/N
! Modified
!$omp parallel default(private) shared(whole_data,N,width_bytes,invert)
!$omp do schedule(dynamic)
do y=1,N-1
bit_num = 0
byte_acc = K0
! Adopted form c++ example
pdata => whole_data((y-1)*width_bytes:(y)*width_bytes)
incr=0
do x=1,N-1,2
C1 = cmplx(x*invert-1.5d0,y*invert-1.0d0)
C2 = cmplx((x+1)*invert-1.5d0,y*invert-1.0d0)
Z1 = C1
Z2 = C2
in_mandelbrot = .true.
res=3
do i=1,iter
! Adopted from C/C++ example
Z1 = Z1**2 + C1
Z2 = Z2**2 + C2
absZ1 = real(Z1*conjg(Z1))
absZ2 = real(Z2*conjg(Z2))
if (absZ2 <= limit2 .AND. absZ1 <= limit2) then
mask = 3
else if (absZ2 > limit2 .AND. absZ1 <= limit2) then
mask = 2
else if (absZ2 <= limit2 .AND. absZ1 > limit2) then
mask = 1
else if (absZ2 > limit2 .AND. absZ1 > limit2) then
mask = 0
end if
res = iand(res,mask)
if (res==0) exit
end do
bit_num = bit_num + 2
byte_acc = ior(ishft(byte_acc,2),res)
if (bit_num == 8) then
! All bits set so output them
incr=incr+1
pdata(incr) = char(byte_acc)
byte_acc = K0
bit_num = 0
end if
end do
! End of a row so left-justify the bits we have and output them
byte_acc = ishft(byte_acc,8-mod(N,8))
pdata(incr) = char(byte_acc)
end do
!$omp end do nowait
!$omp end parallel
write(*,*) whole_data
deallocate(whole_data)
nullify(pdata, whole_data)
end program mandelbrot