1program timefft 2 3! Tests and times one-dimensional FFTs computed by FFTW3 4 5 use, intrinsic :: iso_c_binding 6 use FFTW3 7 8 complex(C_FLOAT_COMPLEX),pointer :: a(:),b(:),c(:) 9 real(C_FLOAT),pointer :: ar(:),br(:) 10 type(C_PTR) :: plan1,plan2 !Pointers to FFTW plans 11 type(C_PTR) :: pa,pb,pc 12 integer(C_INT) iret 13 integer*8 count0,count1,clkfreq 14 character problem*9 15 logical linplace,lcomplex,lthreading 16 17! Get command-line parameters 18 call timefft_opts(npatience,maxthreads,linplace,lcomplex,nfft,problem,nflags) 19 lthreading=maxthreads.ge.1 20 maxthreads=max(1,maxthreads) 21 22 call sgran() ! see C rand generator (used in gran) 23 24! Allocate data arrays 25 pa=fftwf_alloc_complex(int(nfft,C_SIZE_T)) 26 call c_f_pointer(pa,a,[nfft]) 27 call c_f_pointer(pa,ar,[nfft]) 28 29 pb=fftwf_alloc_complex(int(nfft,C_SIZE_T)) 30 call c_f_pointer(pb,b,[nfft]) 31 call c_f_pointer(pb,br,[nfft]) 32 33 pc=fftwf_alloc_complex(int(nfft,C_SIZE_T)) 34 call c_f_pointer(pc,c,[nfft]) 35 36! Initialize FFTW threading 37 if(lthreading) iret=fftwf_init_threads() 38 39! Import FFTW wisdom, if available 40 iret=fftwf_import_wisdom_from_filename(C_CHAR_'wis.dat' // C_NULL_CHAR) 41 42 do i=1,nfft !Generate random data 43 x=gran() 44 y=gran() 45 b(i)=cmplx(x,y) 46 enddo 47 iters=100 48 49 write(*,1000) 501000 format(/'Problem Threads Plan Time Gflops RMS iters'/ & 51 '--------------------------------------------------------') 52 53! Try nthreads = 1,maxthreads 54 do nthreads=1,maxthreads 55 a(1:nfft)=b(1:nfft) !Copy test data into a() 56 call system_clock(count0,clkfreq) 57! Make the plans 58 if(lthreading) call fftwf_plan_with_nthreads(nthreads) 59 if(lcomplex) then 60 if(linplace) then 61 plan1=fftwf_plan_dft_1d(nfft,a,a,-1,nflags) 62 plan2=fftwf_plan_dft_1d(nfft,a,a,+1,nflags) 63 else 64 plan1=fftwf_plan_dft_1d(nfft,a,c,-1,nflags) 65 plan2=fftwf_plan_dft_1d(nfft,c,a,+1,nflags) 66 endif 67 else 68 if(linplace) then 69 plan1=fftwf_plan_dft_r2c_1d(nfft,ar,a,nflags) 70 plan2=fftwf_plan_dft_c2r_1d(nfft,a,ar,nflags) 71 else 72 plan1=fftwf_plan_dft_r2c_1d(nfft,ar,c,nflags) 73 plan2=fftwf_plan_dft_c2r_1d(nfft,c,ar,nflags) 74 endif 75 endif 76 call system_clock(count1,clkfreq) 77 tplan=0.5*float(count1-count0)/float(clkfreq) !Plan time for one transform 78 79 total=0. 80 do iter=1,iters !Do many iterations 81 a=b !Copy test data into a() 82 call system_clock(count0,clkfreq) 83! Compute the transforms 84 if(lcomplex) then 85 if(linplace) then 86 call fftwf_execute_dft(plan1,a,a) 87 call fftwf_execute_dft(plan2,a,a) 88 else 89 call fftwf_execute_dft(plan1,a,c) 90 call fftwf_execute_dft(plan2,c,a) 91 endif 92 else 93 if(linplace) then 94 call fftwf_execute_dft_r2c(plan1,ar,a) 95 call fftwf_execute_dft_c2r(plan2,a,ar) 96 else 97 call fftwf_execute_dft_r2c(plan1,ar,c) 98 call fftwf_execute_dft_c2r(plan2,c,ar) 99 endif 100 endif 101 call system_clock(count1,clkfreq) 102 total=total + float(count1-count0)/float(clkfreq) 103 if(total>=1.0 .and. iter>=10) go to 40 !Cut iterations short ? 104 enddo 105 iter=iters 106 10740 time=0.5*total/iter !Time for one FFT 108 gflops=5.0/(1.e9*time/(nfft*log(float(nfft))/log(2.0))) 109 a(1:nfft)=a(1:nfft)/nfft !Normalize the back-transformed data 110 111! Compute RMS difference between original data and back-transformed data. 112 sq=0. 113 if(lcomplex) then 114 do i=1,nfft 115 sq=sq + real(a(i)-b(i))**2 + aimag(a(i)-b(i))**2 116 enddo 117 else 118 do i=1,nfft 119 sq=sq + (ar(i)-br(i))**2 120 enddo 121 endif 122 rms=sqrt(sq/nfft) 123 124! Display results 125 write(*,1050) problem,nthreads,tplan,time,gflops,rms,iter 1261050 format(a9,i4,f8.3,f10.6,f7.2,f11.7,i5) 127 enddo 128 129! Export accumulated FFTW wisdom 130 iret=fftwf_export_wisdom_to_filename(C_CHAR_'wis.dat' // C_NULL_CHAR) 131 132! Clean up 133 call fftwf_destroy_plan(plan1) 134 call fftwf_destroy_plan(plan2) 135 call fftwf_free(pa) 136 call fftwf_free(pb) 137 call fftwf_free(pc) 138 call fftwf_cleanup_threads() 139 call fftwf_cleanup() 140 141end program timefft 142