Wanted: Mandelbrot vectorizeable in Fortran
Booker Bense
benseb at grumpy.sdsc.edu
Thu May 16 05:05:21 AEST 1991
In article <1991May15.125842.6521 at informatik.uni-erlangen.de> kssingvo at immd4.informatik.uni-erlangen.de (Klaus Singvogel) writes:
>I'm searching for the Mandelbrot-algorithm in Fortran, which can be vectorized
>on a Cray-Ymp.
>Please don't send me the usual algorithm - 'cause it cannot be vectorized with
>the following depedecy of x: ;-)
>
> ...
> DO 30 i=1, X_PIXEL
> DO 30 j=1, Y_PIXEL
> DO 10 k=1, max_iteration
> x = x * x + c
> IF (abs(x).gt.4) GOTO 20
>10 CONTINUE
>20 field(i,j)=k
> ...
>30 CONTINUE
>
>
- This is an easy one !!! Basically all you have to do to get it
to run at about 60 mflops is promote x to a vector and split it into
real and imaginary parts.
i.e.
do 30 k=1,max_it
do 20 j = 1,y_pixel
do 10 i= 1,x_pixel
if ( abs( x(i,j) ) .lt. 5 ) then
x(i,j) = x(i,j) *x(i,j) + c
field(i,j) = field(i,j) + 1
endif
10 continue
20 continue
30 continue
- this loop will vectorize and run at an average of about 60 mflops
the speed depends on the number of times the if condition is true.
- You can play more games with the if test and get the algorithm to run
at about 90-100 mflops. The code below plays games with pointers and
doing ``safe and unsafe'' versions of the update loops. It also takes
advantage of the ``quadratic '' nature of a single scanline through
the M-set.
- Booker C. Bense
prefered: benseb at grumpy.sdsc.edu "I think it's GOOD that everyone
NeXT Mail: benseb at next.sdsc.edu becomes food " - Hobbes
--------cut-here---------------------------------------------------
c
c Your mission, should you choose to accept it, is to incorparate
c this subroutine into your favorite M-set viewer.
c If asked SDSC will disavow all knowledge of this code.
c This article will self-destruct in 30 seconds....
c
c It might work , it might not. It's your problem if it doesn't.
c WARNING: Playing with M-sets is hazardous to your thesis,
c work projects, cpu allocation, etc.......
c
subroutine mandelsr(p0,q0,pinc,qinc,iter,limit,height,
$ width,size)
C...TRANSLATED BY FPP 3.00Z36 09/28/90 13:01:29 OPTOFF=C
integer height,width,size
integer limit,how_far
integer iter(size)
integer i,j,k,index,outdex,backdex
integer fcol,lcol
real p0, q0
real mandelinf
integer maxsize ,maxside
parameter( mandelinf = 5 )
parameter(maxsize = 16385)
parameter( maxside = 129 )
real realst(maxsize), imagst(maxsize)
real realsq(maxsize) , imagsq(maxsize)
real ptemp(maxsize),qtemp(maxsize)
real qrow,pcol(maxside)
c
c A pointer into iter
c
integer iter_ptr1(maxside),iter_ptr2(maxside)
logical first_fnd(maxside), last_fnd(maxside)
logical done,loop_sw(maxside)
integer count_min(maxside)
c
c Job control parameters
c
integer next,too_far
parameter( next = 9 )
parameter( too_far = 10 - next)
CFPP$ ITERATIONS ( height = 64, width = 64 , size = 4096, limit = 248)
CDIR@ IVDEP
do 25 j = 1,height
iter_ptr1(j) = 1
iter_ptr2(j) = width
25 continue
CDIR@ IVDEP
do 50 k = 1,width
pcol(k) = p0 + (k-1)*pinc
50 continue
do 100 j = 1,height
outdex = (j-1)*width
qrow = q0 + (j-1)* qinc
CDIR@ IVDEP
do 75 k = 1,width
index = outdex + k
ptemp(index) = pcol(k)
qtemp(index) = qrow
75 continue
100 continue
fcol = 1
lcol = width
CDIR@ IVDEP
do 200 k = 1,size
realst(k) = 0
imagst(k) = 0
realsq(k) = 0
imagsq(k) = 0
iter(k) = 0
c
c a bias so you can see which ones the cray draws should be iter = 0
c
200 continue
c
c We can get away with this until i > 10
c
how_far = 1
c
do 400 i = how_far,how_far+9
do 300 j = 1,height
outdex = (j-1) *width
CDIR@ IVDEP
do 250 k = 1,width
index = outdex + k
if ( (realsq(index) + imagsq(index) ) .LT. mandelinf )
$ then
iter(index) = iter(index) + 1
endif
imagst(index) = 2.0 * realst(index) *
$ imagst(index)+qtemp(index)
realst(index) = realsq(index) - imagsq(index)
$ + ptemp(index)
realsq(index) = realst(index) * realst(index)
imagsq(index) = imagst(index) * imagst(index)
250 continue
300 continue
400 continue
how_far = how_far+9
c
c Now we check and set iter_ptr!!!
c
c top of safe limit loop
c
401 continue
CDIR@ IVDEP
do 410 j = 1, height
first_fnd(j) = .FALSE.
last_fnd(j) = .FALSE.
count_min(j) = how_far
410 continue
do 425 k = fcol,lcol
CDIR@ IVDEP
do 450 j = 1,height
index = k + (j-1) *width
if ( first_fnd(j) .AND. last_fnd(j)) then
backdex = width - (k-1) + (j-1)*height
if(iter(index).EQ. how_far.AND.
$ .NOT.(first_fnd(j))) then
iter_ptr1(j) = k
first_fnd(j) = .TRUE.
endif
if(iter(backdex).EQ. how_far.AND.
$ .NOT.(last_fnd(j))) then
iter_ptr2(j) = width - (k-1)
first_fnd(j) = .TRUE.
endif
endif
if( iter(index) .LT. count_min(j) ) then
count_min(j) = iter(index)
endif
450 continue
425 continue
c
c Now check if we're done
c
done = .FALSE.
do 475 j = 1,height
done = done.OR.first_fnd(j).OR.last_fnd(j)
475 continue
how_far = how_far + 1
if (( done) .OR. ( how_far .GT. LIMIT) ) then
return
endif
c
c Get first and last col for next time though.
c
lcol = 1
fcol = width
CDIR@ IVDEP
DO 490 J = 1, HEIGHT
FCOL = MIN0(ITER_PTR1(J),FCOL)
LCOL = MAX0(ITER_PTR2(J),LCOL)
IF (HOW_FAR - COUNT_MIN(J) .GT. 1) THEN
LOOP_SW(J) = .TRUE.
ELSE
LOOP_SW(J) = .FALSE.
ENDIF
490 CONTINUE
c
c
c
do 800 i = how_far,how_far+next
do 700 j = 1,height
outdex = (j-1) * width
if ( loop_sw(j) ) then
c
c Safe Loop
c
CDIR@ IVDEP
do 600 k = iter_ptr1(j),iter_ptr2(j)
index = outdex + k
if ( (realsq(index) + imagsq(index) )
$ .LT. mandelinf ) then
iter(index) = iter(index) + 1
imagst(index) = 2.0 * realst(index) *
$ imagst(index)+qtemp(index)
realst(index) = realsq(index) - imagsq(index)
$ + ptemp(index)
realsq(index) = realst(index) * realst(index)
imagsq(index) = imagst(index) * imagst(index)
endif
600 continue
else
c
c Unsafe Loop
c
CDIR@ IVDEP
do 650 k = iter_ptr1(j),iter_ptr2(j)
index = outdex + k
if ( (realsq(index) + imagsq(index) )
$ .LT. mandelinf ) then
iter(index) = iter(index) + 1
endif
imagst(index) = 2.0 * realst(index) *
$ imagst(index)+qtemp(index)
realst(index) = realsq(index) - imagsq(index)
$ + ptemp(index)
realsq(index) = realst(index) * realst(index)
imagsq(index) = imagst(index) * imagst(index)
650 continue
endif
700 continue
800 continue
how_far = how_far+next
if ( how_far .LT. limit ) go to 401
return
end
More information about the Comp.unix.cray
mailing list