Subject: Re: tgcm fft on prospect

From: Ben Foster <[email-masked]>
Date: Fri, 14 Apr 2000 12:39:54 -0600 (MDT)

Rich:

No, it would be real to real. We are not spectral, but
we use fft's for time smoothing and filtering in wave
space. We would need a 2d transform (lon x press) with 
8-byte reals (4 or 8-byte ints), forward and inverse, e.g.:

      real :: fx,wfft,trigs
      COMMON/FFT/FX(ZIMX+2,2*ZKMXP),WFFT((ZIMX+1)*2*ZKMXP),
     |  TRIGS(3*ZIMX/2+1)
      integer :: ifax(13)
      common/fft_int/ ifax
!$OMP THREADPRIVATE (/fft/,/fft_int/)

!
! Call fft for forward transform (nlev==zkmxp for most fields):
! (imaxp2 == imax+2 == zimx+2)
      call fftrans(fx,wfft,trigs,ifax,1,imaxp2,imax,nlev,-1)
!
! Multiply by smoothing function
      do i = 1,imax
        do n=1,nlev
!
! babyblue apparently underflows here with cos(pi/2) (?)
          fx(i+2,n) = fx(i+2,n)*cos(((i+1)/2)*c(1)/2.)**(2*nn(kutj))
        enddo
      enddo
!
! Do reverse transform
      call fftrans(fx,wfft,trigs,ifax,1,imaxp2,imax,nlev,+1)
!

These calls are made from inside the main parallel latitude
loop, whether OpenMP or MPI, where nlat=36. So we do not need a
routine that is itself parallel, rather it must be thread-safe,
i.e., can be called simultaneously (re-entrant?). The in/out array 
fx (and others where necessary) are in thread private common. 
The fft is done over maybe 15-20 or more fields. 

There's also a (serial) set99 call before the main time loop to set 
up coefficients (actually, now that I think of it I'd better make 
sure the coeffs from the set99 call are available to the slaves in 
the lat loop -- I assume they are, since it works under SGI and UNICOS)

In this case imax=72 and nlev=29. Imax (nlon) is constant across 
model versions, but nlev may be different (up to 45 in time-gcm).

Sub fftrans is copied below -- it calls the old ecmfft from the lib 
if UNICOS, or uses source if non-UNICOS. I'll attach the fft9.f 
source which is used on SGI/OpenMP and IBM/MPI. It's not changed 
much at all from the original mid-70's code. Its interesting that it
works w/ OpenMP on the Origin, but on the Compaq only the master
proc gets it right (the results are good if run with a single proc).

I'm not sure what was happening w/ the babyblue comment above, but 
its not a big problem on b-forest. (I'm looking forward to trying 
out the new WH II's on b-blue -- might even try OpenMP or hybrid if 
their new compiler will let me)

In fftrans below, we can obviously use something specific for each
platform if necessary, with pre-proc macros, as long as they all get
roughly the same answers. I looked at the dmxl and dmxlp libs on prospect, 
but they were getting very different answers than this in a small test 
program, and was seg faulting with OpenMP anyway, altho maybe I wasnt 
calling it right.

Thanks a lot for looking at this..

--Ben

!-----------------------------------------------------------------------
      subroutine fftrans(a,work,trigs,ifax,inc,jump,n,lot,isign)
!
! Do fft. If unicos, the original library version, FFT991  is used.
! If non-unicos, call FFT999, which is in fft9.f (the cray lib source,
!   with FFT991 changed to FFT999).
! See also setfft below.
!
      implicit none
!
! Args:
      integer,intent(in) :: inc,jump,n,lot,isign
      real,intent(inout) :: a(n)
      real,intent(in) :: work(n),trigs(n)
      integer,intent(in) :: ifax(13) 
!
#ifdef UNICOS
      call fft991(a,work,trigs,ifax,inc,jump,n,lot,isign)
#else
      call fft999(a,work,trigs,ifax,inc,jump,n,lot,isign)
#endif
      end subroutine fftrans
!-----------------------------------------------------------------------

> From [email-masked] Fri Apr 14 08:40 MDT 2000
> Date: Fri, 14 Apr 2000 08:40:00 -0600 (MDT)
> From: Rich Loft <[email-masked]>
> To: [email-masked]
> Subject: Re: tgcm on prospect
> 
> 
>    Ben:
> 
>    This is great news!
> 
>    Yes, I believe that we have a modified fft from fftpack which will 
>    work for you, I assume you want real to complex, complex to real?
>    Any other details will be appreciated!
> 
> 					      
> 								Rich
> 
>    Date: Thu, 13 Apr 2000 16:47:00 -0600 (MDT)
>    From: Ben Foster <[email-masked]>
>    Content-Length: 598
> 
>    Rich:
> 
>    Ok, tgcm is running with OpenMP on Crays and SGI Origin,
>    and MPI on blackforest. Now for Compaq. It runs with OpenMP
>    on prospect, but only the master thread gets it right. I
>    have determined that the problem is with an old fft from
>    ecmfft lib. Altho the old source (fft991, etc) works fine
>    on the other 3 platforms, only the master is correct on
>    prospect (its ok with one proc). I have confirmed this with 
>    a small test program.
> 
>    I talked w/ Dick V. about it, and he mentioned that your
>    group has done some work on reentrant ffts. He has also
>    worked on fftpack5. What do you suggest?
> 
>    --Ben
> 
>