C  Timing tests for fft demonstrations.
C
C  File:  fftm.f | fftm.fcm


C  Initialize with small test polynomials in
C  axis 1.
C
C   xx = 2 + z +   z^2 + z^3
C   yy = 3 - z + 2*z^2 + z^3
C
C   product = 6 + z + 6*z^2 + 6*z^3 + 2*z^4 + 3*z^5 + z^6

      subroutine poly(xx, yy)
      INCLUDE 'fft.hp'
      complex, array(0:R-1, 0:C-1) :: xx, yy
CMF$  LAYOUT xx(,), yy(,)
CMPF  ONDPU xx, yy

      print *, "-- 2 + z +   z^2 + z^3 --"
      xx = (0.0, 0.0)
      yy = (0.0, 0.0)
      xx( 0 , :) =  2.0
      xx(1:3, :) =  1.0
      call prc8(xx, 1, 1)
      print *
      print *, "-- 3 - z + 2*z^2 + z^3 --"
      yy( 0 , :) =  3.0
      yy( 1 , :) = -1.0
      yy( 2 , :) =  2.0
      yy( 3 , :) =  1.0
      call prc8(yy, 1, 1)
      end subroutine poly


C  Print first 8 entries of complex array z along
C  axis ax at position pn .

      subroutine prc8(z, ax, pn)
      INCLUDE 'fft.hp'
      complex z(0:R-1, 0:C-1)
CMF$  LAYOUT z(,)
CMPF  ONDPU z
      integer ax, pn

      complex u(0:7)
CMF$  LAYOUT u(:serial)
CMPF  ONFE u
      integer k

      if ( ax == 1 ) then
	do, k = 0, 7
         u(k) = z(k, pn)
	end do
       else
	do, k = 0, 7
	  u(k) = z(pn, k)
	end do
      end if
      print 21, real(u(0)), aimag(u(0)), real(u(1)), aimag(u(1))
      print 21, real(u(2)), aimag(u(2)), real(u(3)), aimag(u(3))
      print 21, real(u(4)), aimag(u(4)), real(u(5)), aimag(u(5))
      print 21, real(u(6)), aimag(u(6)), real(u(7)), aimag(u(7))
 21   FORMAT (1X, 2("(", G14.6, ", ", G14.6, ")   "))
      end subroutine prc8


C  Basic tests and timing for fft demonstrations.

      program fftest
      INCLUDE 'fft.hc'
      complex, array(0:R-1, 0:C-1) ::
     *               xx, yy, zz, xt, yt, zt, ww
CMF$  LAYOUT xx(,), yy(,), zz(,), xt(,), yt(,), zt(,), ww(,)
CMPF  ONDPU  xx, yy, zz, xt, yt, zt, ww
      integer eltm

C  Compiler crashes if we put in interfaces from fft.hf,
C  but it needs some to generate correct code.

      INTERFACE
        function Ufft1(Z, N, r, c, ud, inv)
        integer N, r, c
        complex, array(0:r-1, 0:c-1) :: Ufft1, Z
CMF$    LAYOUT Z(,)
CMPF    ONDPU Z
        logical ud, inv
      END INTERFACE

      call fftsetup        ! Initialize
      print 21, R-1, C-1
 21   FORMAT (" Arrays are (0:", I4, ",  0:", I4, ")")
      print *, "-- XX, YY --"
      print *, "Timing with ", R, "-point transforms."
      call poly(xx, yy)    ! Test example
C      call CM_timer_clear(1)  ! *** CM5 ***
C      call CM_timer_start(1)  ! *** CM5 ***
      call mpTimerStart()
      xt = Ufft1(xx, R, R, C, .FALSE., .FALSE.)
      yt = Ufft1(yy, R, R, C, .FALSE., .FALSE.)
      zt = xt*yt           ! Product of transforms
      zz = Ufft1(zt, R, R, C, .TRUE., .TRUE.)/R
C      call CM_timer_stop(1)   ! *** CM5 ***
C      call CM_timer_print(1)  ! *** CM5 ***
      eltm = mpTimerElapsed()
      print *, "Elapsed time = ", eltm, " msec."
      print *, " -- XT --"
      call prc8(xt, 1, 0)
      print *, " -- YT --"
      call prc8(yt, 1, 0)
      print *, " -- ZT --"
      call prc8(zt, 1, 0)
      print *, " -- ZZ --"
      call prc8(zz, 1, 0)
      end program fftest
