Subject: Re: coding question

From: Rich Loft <[email-masked]>
Date: Wed, 14 Jun 2000 21:31:04 -0600 (MDT)

Ben:

   I have been out of town at a conference since last Friday...
   Regarding this day night business FIRST of all, completely independently
   of your THAT, just looking at

   do lam=16,lmax ! wavelength 
     do i=1,nlon ! nlon=76   
       if (idn(i)==1) then
	 do k=1,npress ! npress=29   
	   a(i,k) = s(lam)*b(i,k)
	   ...
	 enddo ! pressure
       endif
     enddo ! longitude   
   enddo ! wavelength 

   With lam on the outside then you are streaming M*(lmax-15)*npress*nlon words 
   through the system for the M ik slice variable like a or b that you manipulate 
   inside the triple loop. Additionally you are loading N*(lmax-15) values of the N
   variables like s. So you are loading  

      (lmax-15)*(M*(npress*nlon)+N) words.

   If you put lam as the inside loop, then you are streaming 

       npress*nlon*(M + N*(lmax-15)) words 

   through the system in the case where lam is innermost. 

   It is very difficult to see what would be best without seeing the whole routine,
   i.e. I have no idea what M and N and lmax are. 

   Now in all cases (as you have discovered) you never want to stride through memory, 
   particularly using a value like 29, since this assures using something like 1/16th 
   of the memory bw, at least on the IBM.

   Ideally, I would suggest trying the following loop order on the code, especially if lmax 
   << npress*nlon, since then we would see cache bandwidths (3+ GB/sec) for these loads:

   do i=1,nlon
      if (idn(i)==1) then
         do k=1,npress
            do lam=16,lmax
               a(k,i)=a(k,i)+s(lam)*b(k,i)
            end do
         end do
      end if
   end do

   Although this will mean transposes, as you observe, it may still be worth it, because we can 
   unroll the transpose for optimal performance. Perhaps you can construct a standalone example of this
   code that can be experimented on independently. Eventually the transposes could go away as 
   more physics gets columnized...

										Rich

=================================================================================================

   Date: Wed, 14 Jun 2000 12:35:36 -0600 (MDT)
   From: Ben Foster <[email-masked]>
   Content-Type: text/plain; charset=us-ascii
   Content-MD5: OOE+2M2X56ViDGQewZV7Cg==
   Content-Length: 2080

   Rich:

   I have something like the following:

   !$OMP PARALLEL DO ...
   do j=1,nlat ! nlat=36
     ...
     call qrj
     ...
   enddo

   sub qrj
   lmax = 59
   do lam=1,lmax ! wavelength
   !
   ! Several of these 2d loops in the lamdas loop:
     do k=1,npress ! npress=29
       do i=1,nlon ! nlon=76
   !
   ! Blocks of (often reduction) statements like:
   ! (a and b are thread-private)
	 b(i,k) = blah,blah
	 a(i,k) = a(i,k)+s(lam)*b(i,k) ! sum over wavelength
       enddo ! longitude
     enddo ! pressure
     ...
   enddo ! wavelength

   The qrj routine is one of the highest cpu consumers
   in the code. In an attempt to speed it up I realized
   that most of this code does not need to be executed
   at night, which is dependent on j and i. So I made an
   index array idn(nlon) (a private one for each j), in
   which idn(i)==0 for nighttime, or 1 for daytime, then 
   did the following in qrj:

   do lam=16,lmax ! wavelength
     do k=1,npress
       do i=1,nlon
	 if (idn(i)==1) then ! daytime
	   ...
	   a(i,k) = s(lam)*b(i,k) 
	   ...
	 endif
       enddo
     enddo
   enddo

   This did not speed it up much, if at all, even tho there 
   were many (altho < 1/2) night-time iterations. So I figured 
   to save overhead on at least the k-loop by reversing the
   i and k loops:

   do lam=16,lmax ! wavelength 
     do i=1,nlon ! nlon=76   
       if (idn(i)==1) then
	 do k=1,npress ! npress=29   
	   a(i,k) = s(lam)*b(i,k)
	   ...
	 enddo ! pressure
       endif
     enddo ! longitude   
   enddo ! wavelength 

   This appeared to be even slower, presumably because of
   the non-stride-1 on a and b. When I put the i-loop outside
   of the lamdas loop, it was *definitely* much slower. 
   There appears to be a trade-off between saving work on 
   the diurnal conditional vs stride-1 on the arrays. 

   I've considered introducing local arrays c(k,i) and
   d(k,i) but some (not all) of the a and b arrays are in 
   common, so for those I would have to dump c=trans(a) and 
   d=trans(b) at the beginning of the routine, and a=trans(c) 
   and b=trans(d) at the end.

   Do you have any suggestions?

   --Ben

   PS: did you get my tgcm benchmarks msg last week?