mirror of https://github.com/saitohirga/WSJT-X.git
A few more tweaks, and add the file nfft.dat of efficient FFT lengths.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@4837 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor
parent
1889b24c78
commit
c628013582
324
lib/chkfft.f90
324
lib/chkfft.f90
|
|
@ -1,162 +1,162 @@
|
||||||
program chkfft
|
program chkfft
|
||||||
|
|
||||||
! Tests and times one-dimensional FFTs computed by four2a().
|
! Tests and times one-dimensional FFTs computed by four2a().
|
||||||
! An all-Fortran version of four2a() is available, but the preferred
|
! An all-Fortran version of four2a() is available, but the preferred
|
||||||
! version uses calls to the FFTW library.
|
! version uses calls to the FFTW library.
|
||||||
|
|
||||||
parameter (NMAX=8*1024*1024) !Maximum FFT length
|
parameter (NMAX=8*1024*1024) !Maximum FFT length
|
||||||
complex a(NMAX),b(NMAX)
|
complex a(NMAX),b(NMAX)
|
||||||
real ar(NMAX),br(NMAX)
|
real ar(NMAX),br(NMAX)
|
||||||
real mflops
|
real mflops
|
||||||
character infile*12,arg*8
|
character infile*12,arg*8
|
||||||
logical list
|
logical list
|
||||||
common/patience/npatience
|
common/patience/npatience
|
||||||
equivalence (a,ar),(b,br)
|
equivalence (a,ar),(b,br)
|
||||||
|
|
||||||
nargs=iargc()
|
nargs=iargc()
|
||||||
if(nargs.ne.5) then
|
if(nargs.ne.5) then
|
||||||
print*,'Usage: chkfft <nfft | infile> nr nw nc np'
|
print*,'Usage: chkfft <nfft | infile> nr nw nc np'
|
||||||
print*,' nfft: length of FFT'
|
print*,' nfft: length of FFT'
|
||||||
print*,' nfft=0: do lengths 2^n, n=2^4 to 2^23'
|
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*,' infile: name of file with nfft values, one per line'
|
||||||
print*,' nr: 0/1 to not read (or read) wisdom'
|
print*,' nr: 0/1 to not read (or read) wisdom'
|
||||||
print*,' nw: 0/1 to not write (or write) wisdom'
|
print*,' nw: 0/1 to not write (or write) wisdom'
|
||||||
print*,' nc: 0/1 for real or complex data'
|
print*,' nc: 0/1 for real or complex data'
|
||||||
print*,' np: 0-4 patience for finding best algorithm'
|
print*,' np: 0-4 patience for finding best algorithm'
|
||||||
go to 999
|
go to 999
|
||||||
endif
|
endif
|
||||||
|
|
||||||
list=.false.
|
list=.false.
|
||||||
nfft=-1
|
nfft=-1
|
||||||
call getarg(1,infile)
|
call getarg(1,infile)
|
||||||
open(10,file=infile,status='old',err=1)
|
open(10,file=infile,status='old',err=1)
|
||||||
list=.true. !A valid file name was provided
|
list=.true. !A valid file name was provided
|
||||||
go to 2
|
go to 2
|
||||||
1 read(infile,*) nfft !Takje first argument to be nfft
|
1 read(infile,*) nfft !Takje first argument to be nfft
|
||||||
2 call getarg(2,arg)
|
2 call getarg(2,arg)
|
||||||
read(arg,*) nr
|
read(arg,*) nr
|
||||||
call getarg(3,arg)
|
call getarg(3,arg)
|
||||||
read(arg,*) nw
|
read(arg,*) nw
|
||||||
call getarg(4,arg)
|
call getarg(4,arg)
|
||||||
read(arg,*) ncomplex
|
read(arg,*) ncomplex
|
||||||
call getarg(5,arg)
|
call getarg(5,arg)
|
||||||
read(arg,*) npatience
|
read(arg,*) npatience
|
||||||
|
|
||||||
if(list) write(*,1000) infile,nr,nw,ncomplex,npatience
|
if(list) write(*,1000) infile,nr,nw,ncomplex,npatience
|
||||||
1000 format(/'infile: ',a12,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/)
|
1000 format(/'infile: ',a12,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/)
|
||||||
if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience
|
if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience
|
||||||
1002 format(/'nfft: ',i10,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/)
|
1002 format(/'nfft: ',i10,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/)
|
||||||
|
|
||||||
open(12,file='chkfft.out',status='unknown')
|
open(12,file='chkfft.out',status='unknown')
|
||||||
open(13,file='fftwf_wisdom.dat',status='unknown')
|
open(13,file='fftwf_wisdom.dat',status='unknown')
|
||||||
|
|
||||||
if(nr.ne.0) then
|
if(nr.ne.0) then
|
||||||
call import_wisdom_from_file(isuccess,13)
|
call import_wisdom_from_file(isuccess,13)
|
||||||
if(isuccess.eq.0) then
|
if(isuccess.eq.0) then
|
||||||
write(*,1010)
|
write(*,1010)
|
||||||
1010 format('Failed to import FFTW wisdom.')
|
1010 format('Failed to import FFTW wisdom.')
|
||||||
go to 999
|
go to 999
|
||||||
endif
|
endif
|
||||||
endif
|
endif
|
||||||
|
|
||||||
idum=-1 !Set random seed
|
idum=-1 !Set random seed
|
||||||
ndim=1 !One-dimensional transforms
|
ndim=1 !One-dimensional transforms
|
||||||
do i=1,NMAX !Set random data
|
do i=1,NMAX !Set random data
|
||||||
x=gran()
|
x=gran()
|
||||||
y=gran()
|
y=gran()
|
||||||
b(i)=cmplx(x,y) !Generate random data
|
b(i)=cmplx(x,y) !Generate random data
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
iters=1000000
|
iters=1000000
|
||||||
if(list .or. (nfft.gt.0)) then
|
if(list .or. (nfft.gt.0)) then
|
||||||
n1=1
|
n1=1
|
||||||
n2=1
|
n2=1
|
||||||
if(nfft.eq.0) n2=999999
|
if(nfft.eq.0) n2=999999
|
||||||
write(*,1020)
|
write(*,1020)
|
||||||
1020 format(' NFFT Time rms MHz MFlops iters', &
|
1020 format(' NFFT Time rms MHz MFlops iters', &
|
||||||
' tplan'/61('-'))
|
' tplan'/61('-'))
|
||||||
else
|
else
|
||||||
n1=4
|
n1=4
|
||||||
n2=23
|
n2=23
|
||||||
write(*,1030)
|
write(*,1030)
|
||||||
1030 format(' n N=2^n Time rms MHz MFlops iters', &
|
1030 format(' n N=2^n Time rms MHz MFlops iters', &
|
||||||
' tplan'/63('-'))
|
' tplan'/63('-'))
|
||||||
endif
|
endif
|
||||||
|
|
||||||
do ii=n1,n2 !Test one or more FFT lengths
|
do ii=n1,n2 !Test one or more FFT lengths
|
||||||
if(list) then
|
if(list) then
|
||||||
read(10,*,end=900) nfft !Read nfft from file
|
read(10,*,end=900) nfft !Read nfft from file
|
||||||
else if(n2.gt.n1) then
|
else if(n2.gt.n1) then
|
||||||
nfft=2**ii !Do powers of 2
|
nfft=2**ii !Do powers of 2
|
||||||
endif
|
endif
|
||||||
|
|
||||||
iformf=1
|
iformf=1
|
||||||
iformb=1
|
iformb=1
|
||||||
if(ncomplex.eq.0) then
|
if(ncomplex.eq.0) then
|
||||||
iformf=0 !Real-to-complex transform
|
iformf=0 !Real-to-complex transform
|
||||||
iformb=-1 !Complex-to-real (inverse) transform
|
iformb=-1 !Complex-to-real (inverse) transform
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if(nfft.gt.NMAX) go to 900
|
if(nfft.gt.NMAX) go to 900
|
||||||
a(1:nfft)=b(1:nfft) !Copy test data into a()
|
a(1:nfft)=b(1:nfft) !Copy test data into a()
|
||||||
t0=second()
|
t0=second()
|
||||||
call four2a(a,nfft,ndim,-1,iformf) !Get planning time for forward FFT
|
call four2a(a,nfft,ndim,-1,iformf) !Get planning time for forward FFT
|
||||||
call four2a(a,nfft,ndim,+1,iformb) !Get planning time for backward FFT
|
call four2a(a,nfft,ndim,+1,iformb) !Get planning time for backward FFT
|
||||||
t2=second()
|
t2=second()
|
||||||
tplan=t2-t0 !Total planning time for this length
|
tplan=t2-t0 !Total planning time for this length
|
||||||
|
|
||||||
total=0.
|
total=0.
|
||||||
do iter=1,iters !Now do many iterations
|
do iter=1,iters !Now do many iterations
|
||||||
a(1:nfft)=b(1:nfft) !Copy test data into a()
|
a(1:nfft)=b(1:nfft) !Copy test data into a()
|
||||||
|
|
||||||
t0=second()
|
t0=second()
|
||||||
call four2a(a,nfft,ndim,-1,iformf) !Forward FFT
|
call four2a(a,nfft,ndim,-1,iformf) !Forward FFT
|
||||||
call four2a(a,nfft,ndim,+1,iformb) !Backward FFT on same data
|
call four2a(a,nfft,ndim,+1,iformb) !Backward FFT on same data
|
||||||
t1=second()
|
t1=second()
|
||||||
total=total+t1-t0
|
total=total+t1-t0
|
||||||
if(total.ge.1.0) go to 40 !Cut iterations short if t>1 s
|
if(total.ge.1.0) go to 40 !Cut iterations short if t>1 s
|
||||||
enddo
|
enddo
|
||||||
iter=iters
|
iter=iters
|
||||||
|
|
||||||
40 time=0.5*total/iter !Time for one FFT of current length
|
40 time=0.5*total/iter !Time for one FFT of current length
|
||||||
tplan=0.5*tplan-time !Planning time for one FFT
|
tplan=0.5*tplan-time !Planning time for one FFT
|
||||||
if(tplan.lt.0) tplan=0.
|
if(tplan.lt.0) tplan=0.
|
||||||
a(1:nfft)=a(1:nfft)/nfft
|
a(1:nfft)=a(1:nfft)/nfft
|
||||||
|
|
||||||
! Compute RMS difference between original array and back-transformed array.
|
! Compute RMS difference between original array and back-transformed array.
|
||||||
sq=0.
|
sq=0.
|
||||||
if(ncomplex.eq.1) then
|
if(ncomplex.eq.1) then
|
||||||
do i=1,nfft
|
do i=1,nfft
|
||||||
sq=sq + real(a(i)-b(i))**2 + imag(a(i)-b(i))**2
|
sq=sq + real(a(i)-b(i))**2 + imag(a(i)-b(i))**2
|
||||||
enddo
|
enddo
|
||||||
else
|
else
|
||||||
do i=1,nfft
|
do i=1,nfft
|
||||||
sq=sq + (ar(i)-br(i))**2
|
sq=sq + (ar(i)-br(i))**2
|
||||||
enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
rms=sqrt(sq/nfft)
|
rms=sqrt(sq/nfft)
|
||||||
|
|
||||||
freq=1.e-6*nfft/time
|
freq=1.e-6*nfft/time
|
||||||
mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0)))
|
mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0)))
|
||||||
if(n2.eq.1 .or. n2.eq.999999) then
|
if(n2.eq.1 .or. n2.eq.999999) then
|
||||||
write(*,1050) nfft,time,rms,freq,mflops,iter,tplan
|
write(*,1050) nfft,time,rms,freq,mflops,iter,tplan
|
||||||
write(12,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)
|
1050 format(i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
|
||||||
else
|
else
|
||||||
write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
|
write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
|
||||||
write(12,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)
|
1060 format(i2,i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
|
||||||
endif
|
endif
|
||||||
if(mod(ii,50).eq.0) call four2a(0,-1,0,0,0)
|
if(mod(ii,50).eq.0) call four2a(0,-1,0,0,0)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
900 continue
|
900 continue
|
||||||
if(nw.eq.1) then
|
if(nw.eq.1) then
|
||||||
rewind 13
|
rewind 13
|
||||||
call export_wisdom_to_file(13)
|
call export_wisdom_to_file(13)
|
||||||
! write(*,1070)
|
! write(*,1070)
|
||||||
!1070 format(/'Exported FFTW wisdom')
|
!1070 format(/'Exported FFTW wisdom')
|
||||||
endif
|
endif
|
||||||
|
|
||||||
999 call four2a(0,-1,0,0,0)
|
999 call four2a(0,-1,0,0,0)
|
||||||
end program chkfft
|
end program chkfft
|
||||||
|
|
|
||||||
|
|
@ -113,3 +113,8 @@ Test data for all transforms is gaussian random noise of zero mean and
|
||||||
standard deviation 1. Tabulated values of "rms" are the
|
standard deviation 1. Tabulated values of "rms" are the
|
||||||
root-mean-square differences between the original data and the
|
root-mean-square differences between the original data and the
|
||||||
back-transfmred data.
|
back-transfmred data.
|
||||||
|
|
||||||
|
File nfft.dat contains all numbers between 2^3 and 2^23 that have
|
||||||
|
no factor greater than 7, followed by their factors. These numbers
|
||||||
|
are good choices for FFT lengths.
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -19,7 +19,7 @@ subroutine four2a(a,nfft,ndim,isign,iform)
|
||||||
! This version of four2a makes calls to the FFTW library to do the
|
! This version of four2a makes calls to the FFTW library to do the
|
||||||
! actual computations.
|
! actual computations.
|
||||||
|
|
||||||
parameter (NPMAX=100) !Max numberf of stored plans
|
parameter (NPMAX=2100) !Max numberf of stored plans
|
||||||
parameter (NSMALL=16384) !Max size of "small" FFTs
|
parameter (NSMALL=16384) !Max size of "small" FFTs
|
||||||
complex a(nfft) !Array to be transformed
|
complex a(nfft) !Array to be transformed
|
||||||
complex aa(NSMALL) !Local copy of "small" a()
|
complex aa(NSMALL) !Local copy of "small" a()
|
||||||
|
|
|
||||||
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue