program chkfft ! Tests and times one-dimensional FFTs computed by FFTW3 use, intrinsic :: iso_c_binding use FFTW3 parameter (NMAX=8*1024*1024) !Maximum FFT length complex a(NMAX),b(NMAX),c(NMAX) real ar(NMAX),br(NMAX),cr(NMAX) real mflops ! integer*8 plan1,plan2 !Pointers to stored plans type(C_PTR) :: plan1,plan2 !Pointers to FFTW plans character infile*12,arg*8 logical list common/patience/npatience equivalence (a,ar),(b,br),(c,cr) ! include 'fftw3.f90' !FFTW definitions nargs=iargc() if(nargs.ne.5) then print*,'Usage: chkfft3 nr nw nc np' print*,' nfft: length of FFT' print*,' nfft=0: do lengths 2^n, n=2^4 to 2^23' print*,' infile: name of file with nfft values, one per line' print*,' nr: 0/1 to not read (or read) wisdom' print*,' nw: 0/1 to not write (or write) wisdom' print*,' nc: 0/1 for real or complex data' print*,' np: 0-4 patience for finding best algorithm' go to 999 endif list=.false. nfft=-1 call getarg(1,infile) open(10,file=infile,status='old',err=1) list=.true. !A valid file name was provided go to 2 1 read(infile,*) nfft !Take first argument to be nfft 2 call getarg(2,arg) read(arg,*) nr call getarg(3,arg) read(arg,*) nw call getarg(4,arg) read(arg,*) ncomplex call getarg(5,arg) read(arg,*) npatience if(list) write(*,1000) infile,nr,nw,ncomplex,npatience 1000 format(/'infile: ',a12,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/) if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience 1002 format(/'nfft: ',i10,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/) nflags=FFTW_ESTIMATE if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT if(npatience.eq.2) nflags=FFTW_MEASURE if(npatience.eq.3) nflags=FFTW_PATIENT if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE open(12,file='chkfft.out',status='unknown') if(nr.ne.0) then isuccess=fftwf_import_wisdom_from_filename('fftwf_wisdom.dat'//char(0)) if(isuccess.eq.1) then write(*,1010) 1010 format('Imported FFTW wisdom.') else write(*,1012) 1012 format('Failed to import FFTW wisdom.') go to 999 endif endif idum=-1 !Set random seed ndim=1 !One-dimensional transforms do i=1,NMAX !Set random data x=gran() y=gran() b(i)=cmplx(x,y) !Generate random data enddo iters=1000000 if(list .or. (nfft.gt.0)) then n1=1 n2=1 if(nfft.eq.-1) n2=999999 write(*,1020) 1020 format(' NFFT Time rms MHz MFlops iters', & ' tplan'/61('-')) else n1=4 n2=23 write(*,1030) 1030 format(' n N=2^n Time rms MHz MFlops iters', & ' tplan'/63('-')) endif do ii=n1,n2 !Test one or more FFT lengths if(list) then read(10,*,end=900) nfft !Read nfft from file else if(n2.gt.n1) then nfft=2**ii !Do powers of 2 endif iformf=1 iformb=1 if(ncomplex.eq.0) then iformf=0 !Real-to-complex transform iformb=-1 !Complex-to-real (inverse) transform endif if(nfft.gt.NMAX) go to 900 a(1:nfft)=b(1:nfft) !Copy test data into a() t0=second() if(ncomplex.ne.0) then plan1=fftwf_plan_dft_1d(nfft,a,c,-1,nflags) plan2=fftwf_plan_dft_1d(nfft,a,c,+1,nflags) else plan1=fftwf_plan_dft_r2c_1d(nfft,ar,c,nflags) plan2=fftwf_plan_dft_c2r_1d(nfft,c,ar,nflags) endif t2=second() tplan=t2-t0 !Total planning time for this length total=0. do iter=1,iters !Now do many iterations a(1:nfft)=b(1:nfft) !Copy test data into a() t0=second() call fftwf_execute_dft(plan1,a,c) call fftwf_execute_dft(plan2,c,a) t1=second() total=total+t1-t0 if(total.ge.1.0) go to 40 !Cut iterations short if t>1 s enddo iter=iters 40 time=0.5*total/iter !Time for one FFT of current length tplan=0.5*tplan-time !Planning time for one FFT if(tplan.lt.0) tplan=0. a(1:nfft)=a(1:nfft)/nfft ! Compute RMS difference between original array and back-transformed array. sq=0. if(ncomplex.eq.1) then do i=1,nfft sq=sq + real(a(i)-b(i))**2 + imag(a(i)-b(i))**2 enddo else do i=1,nfft sq=sq + (ar(i)-br(i))**2 enddo endif rms=sqrt(sq/nfft) freq=1.e-6*nfft/time mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0))) if(n2.eq.1 .or. n2.eq.999999) then write(*,1050) nfft,time,rms,freq,mflops,iter,tplan write(12,1050) nfft,time,rms,freq,mflops,iter,tplan 1050 format(i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1) else write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan write(12,1060) ii,nfft,time,rms,freq,mflops,iter,tplan 1060 format(i2,i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1) endif enddo 900 continue if(nw.eq.1) then ierr=fftwf_export_wisdom_to_filename('fftwf_wisdom.dat'//char(0)) write(*,1070) 1070 format(/'Exported FFTW wisdom') endif call fftwf_destroy_plan(plan1) call fftwf_destroy_plan(plan2) 999 end program chkfft