source code
! The Computer Language Benchmarks Game
! http://benchmarksgame.alioth.debian.org/
!
! contributed by Pascal Parois
program mandelbrot
implicit none
integer, parameter :: iter = 50, vsize=48
integer, parameter :: int8 = selected_int_kind(2)
integer, parameter :: dp = selected_real_kind(15, 307)
integer :: i, j,k,l, pos
real(dp),dimension(vsize) :: cr,ci,zr,zi,tr,ti
real(dp) cte
integer :: bytepos, grid
integer(int8) :: byte
character(len=:), dimension(:),allocatable :: buf ! character buffer to avoid integer formatting on output
integer w
character(len=8) :: argv
character(len=100) :: sbuffer, formatstr
real(dp), parameter :: threshold = 4.0_dp
call get_command_argument(1, argv)
read(argv, *) w
grid=ceiling(w/real(vsize))*vsize
allocate(character(len=grid/8) :: buf(w))
buf=''
cte=2.0_dp/real(w,dp)
!$omp parallel do default(none) shared(grid, buf, cte, w)&
!$omp& private(i, j, k, bytepos, pos, byte)&
!$omp& private(zr, zi, cr, ci, tr, ti)&
!$omp& schedule(guided)
do i=0, w-1
pos=0
byte=0_int8
bytepos=8
do j = 0, grid-1,vsize
ci = cte*real(i,dp)-1.0_dp
do k=0, vsize-1
cr(k+1) = cte*real(j+k,dp)-1.5_dp
end do
zr=0.0_dp
zi=0.0_dp
tr=0.0_dp
ti=0.0_dp
do k=1, iter
zi=2.0_dp*zr*zi+ci
zr=tr-ti+cr
ti=zi*zi
tr=zr*zr
if (all(tr+ti>threshold)) then
exit
end if
end do
do k=1, vsize
bytepos=bytepos-1
if (.not. isnan(tr(k)+ti(k)) .and. tr(k)+ti(k)<threshold .and. j<=w) then
byte = ibset(byte, bytepos)
end if
if(bytepos==0) then
bytepos=8
pos = pos + 1
buf(i+1)(pos:pos) = transfer(byte, 'a')
byte=0_int8
end if
end do
end do
end do
!$omp end parallel do
! pbm header
write(*,'("P4",/,i0," ",i0)') w, w
! print output
write(formatstr, '(A,I0,A)') '(A',ceiling(w/8.0),')'
do i = 1, w
write(*, formatstr, advance='no') buf(i)(1:ceiling(w/8.0))
end do
end program mandelbrot