fft1d.f

网络带宽测试工具

      Program fft1d
! ----------------------------------------------------------------------
! **********************************************************************
! *** This program is part of the EuroBen Efficiency Benchmark       ***
! ***                                                                ***
! *** Copyright: European Benchmark Group p/o                        ***
! ***            Utrecht University, High Perf. Computing Group      ***
! ***            P.O. Box 80195                                      ***
! ***            3508 TD Utrecht                                     ***
! ***            The Netherlands                                     ***
! ***                                                                ***
! *** Author of this program: Aad J. van der Steen                   ***
! *** Date : Spring 1999                                             ***
! **********************************************************************
! ----------------------------------------------------------------------
! --- Program to measure speed of Radix-4 Complex-to-Complex 1-D FFT.
!     NCF efficiency benchmark.
! ----------------------------------------------------------------------
      Use                      numerics
      Implicit                 None
      
      Real(l_), Allocatable :: a(:), b(:), c(:), ur(:), ui(:), wr(:),
     &                         wi(:)
      Integer, Parameter    :: mxcase = 15
      Integer               :: ml(mxcase), nl(mxcase), ok(mxcase)
      Real(l_)              :: time1(mxcase), time2(mxcase), 
     &                         time3(mxcase), corr, er
      Real(l_)              :: rmflint(mxcase), rmfltrn(mxcase),
     &                         rmfltot(mxcase), wclock
      Real(l_), Parameter   :: tpi = 6.283185307179586_l_
      Integer               :: icase, irep, m, n, ncase, nrep,
     &                         mflint, mfltot, mfltrn

      Integer, Parameter    :: nmax = 50, k = 4
      Integer               :: ic
      Real(l_)              :: size(0:nmax-1), speedcase(0:nmax-1)
      Real(l_)              :: q(0:6,0:nmax+1), aw(0:3,0:3), d(0:3,0:3),
     &                         extra(0:1), result
      Real(l_), Parameter   :: eps = 5.0e-4, fraction = 0.2
      Real(l_)              :: low, up, tpp
      Real(l_)              :: tkn, bcoef
      Common /coeffs/          tkn(0:nmax+5), bcoef(0:nmax+1) 
!
! --- External Functions
!     ------------------
      External                 autint, spoint, wclock
      Real(l_)              :: autint, spoint
! ----------------------------------------------------------------------
! --- Call status routine.

      Call state( 'fft1d   ' )
     
! ----------------------------------------------------------------------
! --- Read problem sizes and initialise arrays.

      ncase = 0
      Open ( 1, File = 'fft1d.in' )
      Read( 1, * ) tpp
      Read( 1, * ) low, up
      Print 900,   low, up, tpp
      ic  = 0
  100 Read ( 1, *, End = 500 ) m, nrep
      ncase     = ncase + 1
      n         = 2**m
      ml(ncase) = m
      nl(ncase) = n
      If ( ncase > mxcase ) Then
         Print *, '*** No. of cases larger than mxcase(=15): ',
     &            'Increase this parameter please.'
         ncase = mxcase
         Go To 500
      End If
      Allocate( a(n), b(n), c(n), ui(n), ur(n), wi(n), wr(n) )
      Call gendat( c, n )
! ----------------------------------------------------------------------
! --- Repeat FFT 'nrep' times for this problem size.

      time1(ncase) = wclock()
      Do irep = 1, nrep
         Call cparrs( c, a, b, n )
         Call cfft4( 0, m, ur, ui, a, b, wr, wi ) !<-- You can insert
         Call cfft4( 1, m, ur, ui, a, b, wr, wi ) !    here your own FFT
      End Do
      time1(ncase) = wclock() - time1(ncase)
! ----------------------------------------------------------------------
! --- Check for errors and correct timing for filling of arrays.

      Call check( a, b, m, n, er )
      ok(ncase) = er
      corr = wclock()
      Do irep = 1, nrep
         Call cparrs( c, a, b, n )
      End Do
      time1(ncase) = ( time1(ncase) - wclock() + corr )/nrep
! ----------------------------------------------------------------------
! --- Measure initialisation time.

      time2(ncase) = wclock()     
      Do irep = 1, nrep
         Call cfft4( 0, m, ur, ui, a, b, wr, wi )
      End Do     
      time2(ncase) = ( wclock() - time2(ncase) )/nrep
      time3(ncase) = time1(ncase) - time2(ncase)

      Deallocate( a, b, c, ui, ur, wi, wr )
! ----------------------------------------------------------------------
! --- Get new problem.

      Go To 100
! ----------------------------------------------------------------------
! --- Calculate no. of flops and Mflop-rate in initialisation and
!     transformation. Print result.

  500 Print 1000
      Do icase = 1, ncase
         size(icase-1) = Log( Real( nl(icase), l_ ) )
         Call nflops ( ml(icase), mflint, mfltrn )
         rmflint(icase) = 1.0e-6_l_*( mflint/time2(icase) )
         rmfltrn(icase) = 1.0e-6_l_*( mfltrn/time3(icase) )
         rmfltot(icase) = 1.0e-6_l_*( ( mflint + mfltrn )
     &                    / time1(icase) )
         speedcase(icase-1) = rmfltot(icase)
         rmflint(icase) = ( 1.0e-6_l_*mflint )/time2(icase)
         Print 1010, nl(icase), time1(icase), time3(icase), 
     &               time2(icase), ok(icase)
         End Do       
         Print 1020
         Print 1030
         Do icase = 1, ncase 
            Print 1040, nl(icase), rmfltot(icase), rmfltrn(icase), 
     &                  rmflint(icase)
      End Do
      ic = ncase
      Call bsplcf( size, speedcase, ic, k, 'FREE', extra, q, aw, d, tkn,
     &             bcoef )
      result = autint( spoint, Log( low ), Log( up ), eps )
      result = result/(Log( (up - low) )*tpp)
      Print 1060, result, fraction

      Print 1050
! ----------------------------------------------------------------------
  900 Format( '1-D FFT test:'/
     &        'Lower bound       =', g13.5/
     &        'Upper bound       =', g13.5/
     &        'Theor. Peak Perf. =', g13.5, ' Mflop/s'/
     &        '-----------------------------------------------',
     &        '-------------------' )
 1000 Format( //'                  FFT results, Radix-4 algorithm',/
     &       1x, 65('-')/
     &       '  Length |  Total Time  | Transform Time |  Init. Time  |',
     &       '  Check |'/
     &       '    N =  |   (sec)      |     (sec)      |     (sec)    |',
     &       ' errors:|'/
     &       1x, 65('-') )
 1010 Format ( i8, ' |', 1pg13.5, ' | ', 1pg13.5, '  |', 1pg13.5, ' | ',
     &       i6, 1x, '|' )
 1020 Format ( 1x, 65('-') )
 1030 Format ( ///'                 FTT results: Mflop-rates',/
     &        1x, 59('-')/
     &        '  Length | Total Transf. | Transf. only  | ',
     &        ' Initialisation |'/
     &        '    N =  |   (Mflop/s)   |   (Mflop/s)   | ',
     &        '    (Mflop/s)   |'/
     &        1x, 59('-') )
 1040 Format ( i8, ' | ', 1pg13.5, ' | ', 1pg13.5, ' |   ', 1pg13.5, 1x,
     &        '|' )
 1050 Format ( 1x, 59('-') )
 1060 Format( '-----------------------------------------------',
     &        '-------------------',
     &        /'Fraction found = ', g11.4, '  Fraction required = ',
     &        g11.4 )
! ----------------------------------------------------------------------
      End Program fft1d