kssingvo@immd4.informatik.uni-erlangen.de (Klaus Singvogel) (05/15/91)
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 Thank's in advance, Klaus. --- Name : Klaus Singvogel Nickname: Pacman E-Mail : kssingvo@medusa.informatik.uni-erlangen.de
benseb@grumpy.sdsc.edu (Booker Bense) (05/16/91)
In article <1991May15.125842.6521@informatik.uni-erlangen.de> kssingvo@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@grumpy.sdsc.edu "I think it's GOOD that everyone NeXT Mail: benseb@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