! ! ****************** ! * pfft.f90 * ! * R. J. Purser * ! * NOAA/NCEP/EMC * ! * 1994 * ! ****************** ! Latest revisions: April 2003 ! April 2011 (include radix-7, 11, 13) ! Sept 2012 (replace "entry" statements by subroutines) ! ! NOAA/NCEP/Environmental Modeling Center. jim.purser@noaa.gov ! Suite of Fast Fourier Transform routines based on factors 2, 3, 5, 7, ! 11 ("B") and 13 ("D").. ! Data assumed from 0 to n-1 and transforms performed in place. ! For complex transforms of real data, real components of wavenumber k ! are found at positions k [0,..,n/2], imag components of wavenumber -k ! are found at positions n-k [n/2+1,..n-1]. ! ! DEPENDENCIES ! Modules: pkind, pietc !============================================================================= !module fft23457bd !============================================================================= !implicit none !integer ln2,ln3,ln4,ln5,ln7,lnb,lnd,ln,n,nm,nh !end module fft23457bd !============================================================================= module pfft !============================================================================= use pkind, only: sp,dp implicit none private public :: fftco,fftcop,cfft,dfft,csft,dsft,dsfe,hsfe, & rfft,hfft,fftcnv,fconv interface snfftln; module procedure infftln; end interface interface sget2357bd; module procedure iget2357bd; end interface interface get23457bd; module procedure iget23457bd; end interface interface rumble; module procedure srumble,drumble; end interface interface fftco; module procedure sfftco, dfftco; end interface interface fftcop; module procedure sfftcop, dfftcop; end interface interface cfft module procedure scfft, dcfft, scfftp, dcfftp, zcfft; end interface interface dfft module procedure sdfft, ddfft, sdfftp, ddfftp, zdfft; end interface interface gfft; module procedure sgfft, dgfft; end interface interface csft; module procedure scsft, dcsft; end interface interface dsft; module procedure sdsft, ddsft; end interface interface dsfe; module procedure sdsfe, ddsfe; end interface interface hsfe; module procedure shsfe, dhsfe; end interface interface rfft module procedure srfft, drfft, srfftp, drfftp; end interface interface hfft module procedure shfft, dhfft, shfftp, dhfftp; end interface interface fftcnv module procedure sfftcnv,dfftcnv,sfftcnvp, dfftcnvp; end interface interface fconv; module procedure sfconv, dfconv; end interface contains !============================================================================= subroutine infftln(n1,m,nfftln)! [snfftln] !============================================================================= ! Divide out as many factors, m, of given n1 as is possible, returning ! this number of factors as nfftln, and returning n1 as the quotient. !============================================================================= integer,intent(INOUT):: n1 integer,intent(IN ):: m integer,intent(OUT ):: nfftln !----------------------------------------------------------------------------- integer :: i,n2 !============================================================================= nfftln=0 do i=1,30 n2=n1/m if(n2*m /= n1)return n1=n2 nfftln=i enddo end subroutine infftln !============================================================================= subroutine iget2357bd(nin,get2357bd)! [sget2357bd] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! ! Factorize NIN in terms of 2's, 3's, 4's, 5's, 7's, 11's and 13's and ! verify that these include the only prime factors. Set LN2, LN3, LN4, LN5, ! LN7, LNB, LND respectively to the number of powers of: 2, 3, 4, 5, 7, ! 11, 13; set LN to be the sum, ! LN2+LN3+LN4*2+LN5+LN7+LNB+LND, and initialize n, nm, nh, ! all in fft23457bd. ! ! --> NIN Number of data along the line of Fourier transformation !============================================================================= use fft23457bd logical,intent(OUT):: get2357bd integer,intent(IN ):: nin !----------------------------------------------------------------------------- integer :: k !============================================================================= k=nin call snfftln(k,13,lnd) call snfftln(k,11,lnb) call snfftln(k,7,ln7) call snfftln(k,5,ln5) call snfftln(k,4,ln4) call snfftln(k,3,ln3) call snfftln(k,2,ln2) ln=ln2+ln4*2+ln3+ln5+ln7+lnb+lnd n=(2**ln2)*(3**ln3)*(4**ln4)*(5**ln5)*(7**ln7)*(11**lnb)*(13**lnd) nm=n-1; nh=n/2 if(n /= nin)stop 'prime factors of fft period are not only 2, 3, 5, 7, 11, 13' get2357bd=(k==1) end subroutine iget2357bd !============================================================================= subroutine iget23457bd(pl2,pl3,pl4,pl5,pl7,plb,pld)! [get23457bd] !============================================================================= use fft23457bd integer, intent(OUT):: pl2,pl3,pl4,pl5,pl7,plb,pld pl2=ln2; pl3=ln3; pl4=ln4; pl5=ln5; pl7=ln7; plb=lnb; pld=lnd end subroutine iget23457bd !============================================================================= subroutine srumble(jumble,tumble)! [rumble] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! ! Initialize coefficient arrays JUMBLE and TUMBLE for use with the ! fast-Fourier-transform routines when the number N of data has the prime- ! factorization 2**LN2*3**LN3*5**LN5*7**LN7*11**LNB*13**LND ! ! <-- JUMBLE: Permutation of N data encoded as a sequence of transpositions ! <-- TUMBLE: Trigonometric coefficients for use in FFT. The first half are ! the cosines, the second half the sines, of uniformly increasing ! relevant angles. !============================================================================= use pietc, only: pi2 use fft23457bd integer, dimension(0:*),intent(OUT):: jumble real(sp),dimension(0:*),intent(OUT):: tumble !----------------------------------------------------------------------------- integer, parameter :: ml=30 integer :: i,j,l,id,is,ir,kd integer, dimension(ml) :: nd,md real(sp) :: ang,pi2on !============================================================================= pi2on=pi2/n do i=0,nh-1 ang=pi2on*i tumble(i) =cos(ang); tumble(i+nh)=sin(ang) enddo id=1; is=0 do i=1,lnd; is=is+1; md(is)=id; id=id*13; enddo do i=1,lnb; is=is+1; md(is)=id; id=id*11; enddo do i=1,ln7; is=is+1; md(is)=id; id=id*7; enddo do i=1,ln5; is=is+1; md(is)=id; id=id*5; enddo do i=1,ln3; is=is+1; md(is)=id; id=id*3; enddo do i=1,ln2+ln4*2; is=is+1; md(is)=id; id=id*2; enddo id=1 do i=1,ln2+ln4*2; nd(is)=id; id=id*2; is=is-1; enddo do i=1,ln3; nd(is)=id; id=id*3; is=is-1; enddo do i=1,ln5; nd(is)=id; id=id*5; is=is-1; enddo do i=1,ln7; nd(is)=id; id=id*7; is=is-1; enddo do i=1,lnb; nd(is)=id; id=id*11; is=is-1; enddo do i=1,lnd; nd(is)=id; id=id*13; is=is-1; enddo jumble(0)=n do i=1,nm ir=i; j=0 do l=1,ln; kd=ir/nd(l); ir=ir-kd*nd(l); j=j+kd*md(l); enddo jumble(i)=j enddo do i=1,nm j=jumble(i) if(j < i)then 400 j=jumble(j); if(j < i)goto 400 jumble(i)=j endif enddo end subroutine srumble !============================================================================= subroutine drumble(jumble,tumble)! [rumble] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 !============================================================================= use pietc, only: pi2 use fft23457bd integer, dimension(0:*),intent(OUT):: jumble real(dp),dimension(0:*),intent(OUT):: tumble !----------------------------------------------------------------------------- integer, parameter :: ml=30 integer :: i,j,l,id,is,ir,kd integer, dimension(ml) :: nd,md real(dp) :: ang,pi2on !============================================================================= pi2on=pi2/n do i=0,nh-1 ang=pi2on*i tumble(i) =cos(ang); tumble(i+nh)=sin(ang) enddo id=1; is=0 do i=1,lnd; is=is+1; md(is)=id; id=id*13; enddo do i=1,lnb; is=is+1; md(is)=id; id=id*11; enddo do i=1,ln7; is=is+1; md(is)=id; id=id*7; enddo do i=1,ln5; is=is+1; md(is)=id; id=id*5; enddo do i=1,ln3; is=is+1; md(is)=id; id=id*3; enddo do i=1,ln2+ln4*2; is=is+1; md(is)=id; id=id*2; enddo id=1 do i=1,ln2+ln4*2; nd(is)=id; id=id*2; is=is-1; enddo do i=1,ln3; nd(is)=id; id=id*3; is=is-1; enddo do i=1,ln5; nd(is)=id; id=id*5; is=is-1; enddo do i=1,ln7; nd(is)=id; id=id*7; is=is-1; enddo do i=1,lnb; nd(is)=id; id=id*11; is=is-1; enddo do i=1,lnd; nd(is)=id; id=id*13; is=is-1; enddo jumble(0)=n do i=1,nm ir=i; j=0 do l=1,ln; kd=ir/nd(l); ir=ir-kd*nd(l); j=j+kd*md(l); enddo jumble(i)=j enddo do i=1,nm j=jumble(i) if(j < i)then 400 j=jumble(j); if(j < i)goto 400 jumble(i)=j endif enddo end subroutine drumble !============================================================================= subroutine sfftco(n,j,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd)! [fftco] !============================================================================= ! R. J. Purser, NCEP, Washington D.C. 1999. jim.purser@noaa.gov ! SUBROUTINE SFFTCO ! Provide the FFT integer coefficients, j, and real coefficients, w, ! corresponding to a period, n, which factors into powers of 2, 3, 5, 7, 11, 13 ! only. For short period transforms, retain (in js, ws) up to two previously ! used sets of these coefficients to avoid the cost of recalculating them ! in repeated applications. Then the first element of the relevant ! column of js is always the period, n, making it easy to recognize ! whether or not the new requirements for j and w are met in previously ! recorded pairs, js and ws. ! For long period transforms, it is worth supplying an integer array, js, ! and a real array, ws, both of size n, to store the transform parameters ! between applications so that they do not have to be recomputed each time. ! This is found to save significant time in large scale repeated applications. ! In that case, the relevant version of this routine is SFFTCOP. ! ! --> n: period of data ! --> j: array of n permutation indices to unscramble the fft output. ! --> w: array of n real trigonometric coefficients for the fft. !============================================================================= integer, parameter :: bsize=2048,bsizem=bsize-1,bsize2=bsize*2 integer, intent(IN) :: n integer, dimension(0:n-1),intent(OUT):: j real(sp),dimension(0:n-1),intent(OUT):: w integer, intent(OUT):: ln2,ln3,ln4,ln5,ln7,lnb,lnd !----------------------------------------------------------------------------- real(sp),parameter :: zero=0 logical :: get2357bd integer :: ult,nm integer, dimension(0:bsizem,2):: js real(sp), dimension(0:bsizem,2):: ws data ult/1/ ! Column index of js and ws for latest fft coefficients used. data js/bsize2*0/,ws/bsize2*zero/ !============================================================================= nm=n-1 call sget2357bd(n,get2357bd) if(.not. get2357bd)stop 'prime factors are not only 2, 3, 5, 7, 11, 13' call get23457bd(ln2,ln3,ln4,ln5,ln7,lnb,lnd) if(n <= bsize)then if(n == js(0,ult))then j(0:nm)=js(0:nm,ult); w(0:nm)=ws(0:nm,ult) ! copy existing values return endif ult=3-ult ! Try the alternative location if(n == js(0,ult))then j(0:nm)=js(0:nm,ult); w(0:nm)=ws(0:nm,ult) ! copy existing values return endif call rumble(j,w) js(0:nm,ult)=j(0:nm); ws(0:nm,ult)=w(0:nm) else call rumble(j,w) ! Too big to record in js and ws endif end subroutine sfftco !============================================================================= subroutine dfftco(n,j,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd)! [fftco] !============================================================================= ! R. J. Purser, NCEP, Washington D.C. 1999. jim.purser@noaa.gov ! SUBROUTINE DFFTCO: double precision version of fftco !============================================================================= integer, parameter:: bsize=2048,bsizem=bsize-1,bsize2=bsize*2 integer, intent(IN) :: n integer, dimension(0:n-1),intent(OUT):: j real(dp),dimension(0:n-1),intent(OUT):: w integer, intent(OUT):: ln2,ln3,ln4,ln5,ln7,lnb,lnd !----------------------------------------------------------------------------- real(dp),parameter :: zero=0 logical :: get2357bd integer :: ult,nm integer, dimension(0:bsizem,2) :: js real(dp),dimension(0:bsizem,2) :: ws data ult/1/ ! Column index of js and ws for latest fft coefficients used. data js/bsize2*0/,ws/bsize2*zero/ !============================================================================= nm=n-1 call sget2357bd(n,get2357bd) if(.not. get2357bd)stop 'prime factors are not only 2, 3, 5, 7, 11, 13' call get23457bd(ln2,ln3,ln4,ln5,ln7,lnb,lnd) if(n <= bsize)then if(n == js(0,ult))then j(0:nm)=js(0:nm,ult); w(0:nm)=ws(0:nm,ult) ! copy existing values return endif ult=3-ult ! Try the alternative location if(n == js(0,ult))then j(0:nm)=js(0:nm,ult); w(0:nm)=ws(0:nm,ult) ! copy existing values return endif call rumble(j,w) js(0:nm,ult)=j(0:nm); ws(0:nm,ult)=w(0:nm) else call rumble(j,w) ! Too big to record in js and ws endif end subroutine dfftco !============================================================================= subroutine sfftcop(n,js,ws,ln2,ln3,ln4,ln5,ln7,lnb,lnd)! [fftcop] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1), intent(INOUT):: js real(sp),dimension(0:n-1), intent(INOUT):: ws integer, intent( OUT):: ln2,ln3,ln4,ln5,ln7,lnb,lnd !----------------------------------------------------------------------------- logical :: get2357bd !============================================================================= call sget2357bd(n,get2357bd) if(.not. get2357bd)stop 'prime factors are not only 2, 3, 5, 7, 11, 13' call get23457bd(ln2,ln3,ln4,ln5,ln7,lnb,lnd) if(n /= js(0))call rumble(js,ws) end subroutine sfftcop !============================================================================= subroutine dfftcop(n,js,ws,ln2,ln3,ln4,ln5,ln7,lnb,lnd)! [fftcop] !============================================================================= ! R. J. Purser, NCEP, Washington D.C. 1999. jim.purser@noaa.gov ! double precision version of fftcop with the saved parameter arrays ! provided through the parameter list. This allows the user to provide ! these parameter arrays of the appropriate size for the duration of the ! fft tasks and does not require permanent storage to be set aside for them. ! This is advantageous when the FFT problem size might be very large, or ! where many different-sized FFT tasks are alternated repetitively (which ! would necessitate redundant recalculations of these parameters is only ! a fixed size pair of work arrays were reserved). !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1), intent(INOUT):: js real(dp),dimension(0:n-1), intent(INOUT):: ws integer, intent( OUT):: ln2,ln3,ln4,ln5,ln7,lnb,lnd !----------------------------------------------------------------------------- logical :: get2357bd !============================================================================= call sget2357bd(n,get2357bd) if(.not. get2357bd)stop 'prime factors are not only 2, 3, 5, 7, 11, 13' call get23457bd(ln2,ln3,ln4,ln5,ln7,lnb,lnd) if(n /= js(0))call rumble(js,ws) end subroutine dfftcop !============================================================================= subroutine zcfft(n,z)! [cfft] !============================================================================= ! Like cfft, but with r and q combined into complex z. !============================================================================= use pkind,only: dp,dpc integer, intent(IN ):: n complex(dpc),dimension(0:n-1),intent(INOUT):: z !----------------------------------------------------------------------------- real(dp), dimension(0:n-1) :: r,q !============================================================================= r=real(z); q=aimag(z) call cfft(n,r,q) z=cmplx(r,q,dpc) end subroutine zcfft !============================================================================= subroutine zdfft(n,z)! [dfft] !============================================================================= ! Like dfft, but with r and q combined into complex z. !============================================================================= use pkind,only: dp,dpc integer, intent(IN ):: n complex(dpc),dimension(0:n-1),intent(INOUT):: z !----------------------------------------------------------------------------- real(dp), dimension(0:n-1) :: r,q !============================================================================= r=real(z); q=aimag(z) call dfft(n,r,q) z=cmplx(r,q,dpc) end subroutine zdfft !============================================================================= subroutine sdfft(n,r,q)! [dfft] !============================================================================= use pietc, only: one integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer :: nm real(sp) :: rfac !============================================================================= rfac=one/n; nm=n-1 ! for fourier synthesis, scale, and reverse the order of wavenumbers: r(0)=r(0)*rfac; r(1:nm)=r(nm:1:-1)*rfac q(0)=q(0)*rfac; q(1:nm)=q(nm:1:-1)*rfac call scfft(n,r,q) end subroutine sdfft !============================================================================= subroutine scfft(n,r,q)! [cfft] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! (Revised for f90, 1999; radix-7, 11, 13 included, 2011) !============================================================================= integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer :: ln2,ln3,ln4,ln5,ln7,lnb,lnd integer, dimension(0:n-1) :: jumble real(sp),dimension(0:n-1) :: w !============================================================================= call fftco(n,jumble,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd) call gfft(n,jumble,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd,r,q) end subroutine scfft !============================================================================= subroutine ddfft(n,r,q)! [dfft] !============================================================================= use pietc, only: one integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer :: nm real(dp) :: rfac !============================================================================= rfac=one/n; nm=n-1 ! for fourier synthesis, scale, and reverse the order of wavenumbers: r(0)=r(0)*rfac; r(1:nm)=r(nm:1:-1)*rfac q(0)=q(0)*rfac; q(1:nm)=q(nm:1:-1)*rfac call dcfft(n,r,q) end subroutine ddfft !============================================================================= subroutine dcfft(n,r,q)! [cfft] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! (Revised for f90, 1999; radix-7, 11, 13 included, 2011) ! SUBROUTINES DCFFT, DDFFT ! Double precision versions of cfft, dfft !============================================================================= integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer :: ln2,ln3,ln4,ln5,ln7,lnb,lnd integer, dimension(0:n-1) :: jumble real(dp),dimension(0:n-1) :: w !============================================================================= call fftco(n,jumble,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd) call gfft(n,jumble,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd,r,q) end subroutine dcfft !============================================================================= subroutine sdfftp(n,js,ws,r,q)! [dfft] !============================================================================= use pietc, only: one integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(sp),dimension(0:n-1),intent(INOUT):: ws real(sp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer :: nm real(sp) :: rfac !============================================================================= rfac=one/n; nm=n-1 ! for fourier synthesis, scale, and reverse the order of wavenumbers: r(0)=r(0)*rfac; r(1:nm)=r(nm:1:-1)*rfac q(0)=q(0)*rfac; q(1:nm)=q(nm:1:-1)*rfac call scfftp(n,js,ws,r,q) end subroutine sdfftp !============================================================================= subroutine ddfftp(n,js,ws,r,q)! [dfft] !============================================================================= use pietc, only: one integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(dp),dimension(0:n-1),intent(INOUT):: ws real(dp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer :: nm real(dp) :: rfac !============================================================================= rfac=one/n; nm=n-1 ! for fourier synthesis, scale, and reverse the order of wavenumbers: r(0)=r(0)*rfac; r(1:nm)=r(nm:1:-1)*rfac q(0)=q(0)*rfac; q(1:nm)=q(nm:1:-1)*rfac call dcfftp(n,js,ws,r,q) end subroutine ddfftp !============================================================================= subroutine scfftp(n,js,ws,r,q)! [cfft] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! (Revised for f90, 1999; radix-7, 11, 13, included and user-provided parameter ! arrays an option, 2011) !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(sp),dimension(0:n-1),intent(INOUT):: ws real(sp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer:: ln2,ln3,ln4,ln5,ln7,lnb,lnd !============================================================================= call fftcop(n,js,ws,ln2,ln3,ln4,ln5,ln7,lnb,lnd) call gfft(n,js,ws,ln2,ln3,ln4,ln5,ln7,lnb,lnd,r,q) end subroutine scfftp !============================================================================= subroutine dcfftp(n,js,ws,r,q)! [cfft] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! (Revised for f90, 1999; radix-7, 11, 13, included and user-provided parameter ! arrays an option, 2011) ! Double precision versions of cfft, dfft !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(dp),dimension(0:n-1),intent(INOUT):: ws real(dp),dimension(0:n-1),intent(INOUT):: r,q !------------------------------------------------------------------------------ integer:: ln2,ln3,ln4,ln5,ln7,lnb,lnd !============================================================================= call fftcop(n,js,ws,ln2,ln3,ln4,ln5,ln7,lnb,lnd) call gfft(n,js,ws,ln2,ln3,ln4,ln5,ln7,lnb,lnd,r,q) end subroutine dcfftp !============================================================================= subroutine sgfft(n,jumble,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd,r,q)! [gfft] !============================================================================= use pietc, only: one,half integer, intent(IN ):: n integer, dimension(0:n-1),intent(IN ):: jumble real(sp),dimension(0:n-1),intent(IN ):: w integer, intent(IN ):: ln2,ln3,ln4,ln5,ln7,lnb,lnd real(sp),dimension(0:n-1),intent(INOUT):: r,q !----------------------------------------------------------------------------- integer :: nm,nh,i,j,l,k0,k1,k2,k3,k4,k5,k6,k7,k8,k9,ka,kb,kc, & ma,ma2,ma3,ma4,ma5,ma7,mab,mad,mb,mb2, & jmb,jmb2,jmb3,jmb4,jmb5,jmb6, & nze real(sp):: r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,ra,rb,rc, & q0,q1,q2,q3,q4,q5,q6,q7,q8,q9,qa,qb,qc, & rf1,rf2,rf3,rf4,rf5,rf6,rf7,rf8,rf9,rfa,rfb,rfc,& qf1,qf2,qf3,qf4,qf5,qf6,qf7,qf8,qf9,qfa,qfb,qfc, & rep,rze,rzc,ret,rec,rth,rtc,rio,ric,rka,rkc,rla,rlc, & qep,qze,qet,qth,qio,qka,qla, & t1,t2,t3,t4,t5,t6 !============================================================================= nm=n-1 nh=n/2 ! PERMUTE THE DATA: do i=1,nm j=jumble(i) if(j > i)then t1=r(i); r(i)=r(j); r(j)=t1 t1=q(i); q(i)=q(j); q(j)=t1 endif enddo ! TRANSFORM THE DATA: ma=1; mb=n ! RADIX 4 do l=1,ln4 mb=mb/4; ma4=ma*4; mb2=mb*2 do j=0,ma-1 jmb=j*mb; jmb2=j*mb2 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=rf1*rf2-qf1*qf2; qf3=rf1*qf2+qf1*rf2 ! f**3 do i=0,nm,ma4 k0=i+j; k1=k0+ma; k2=k1+ma; k3=k2+ma r0=r(k0); r1=r(k1); r2=r(k2); r3=r(k3) q0=q(k0); q1=q(k1); q2=q(k2); q3=q(k3) t1 =r3*rf3-q3*qf3 ! q13 q3 =r3*qf3+q3*rf3 ! r13 r3 =r2*rf1-q2*qf1 ! r12 r2 =r2*qf1+q2*rf1 ! q12 q2 =q3 -r2 ! r23 r2 =q3 +r2 ! q22 q3 =r3 +t1 ! r22 t1 =r3 -t1 ! q23 r3 =r1*rf2-q1*qf2 ! r11 q1 =r1*qf2+q1*rf2 ! q11 r1 =r0 -r3 ! r21 r0 =r0 +r3 ! r20 r(k3) =r1 -q2 ! r43 r(k1) =r1 +q2 ! r41 q2 =q0 +q1 ! q20 q1 =q0 -q1 ! q21 q(k0) =q2 +r2 ! q40 q(k2) =q2 -r2 ! q42 r(k2) =r0 -q3 ! r42 r(k0) =r0 +q3 ! r40 q(k3) =q1 -t1 ! q43 q(k1) =q1 +t1 ! q41 enddo enddo ma=ma4 enddo if(ln2==1)then ! RADIX 2 mb=mb/2; ma2=ma*2 do j=0,ma-1 jmb=j*mb do i=0,nm,ma2 k0=j+i; k1=k0+ma rf1=w(jmb); qf1=w(nh+jmb) ! f**1 r0=r(k0); q0=q(k0) r1=r(k1); q1=q(k1) t1 =r1*qf1+q1*rf1 ! q11 q1 =r1*rf1-q1*qf1 ! r11 r(k1)=r0 -q1 ! r41 r(k0)=r0 +q1 ! r40 q(k1)=q0 -t1 ! q41 q(k0)=q0 +t1 ! q40 enddo enddo ma=ma2 endif ! RADIX 3 rep=-half; rec=3*half; qep=half*sqrt(3*one) ! <- epsilon do l=1,ln3 mb=mb/3; ma3=ma*3 do j=0,ma-1 jmb=j*mb rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=rf1*rf1-qf1*qf1; qf2=2*rf1*qf1 ! f**2 do i=0,nm,ma3 k0=i+j; k1=k0+ma; k2=k1+ma r1=r(k1); q1=q(k1) r2=r(k2); q2=q(k2) t1 = r2*qf2+q2 *rf2 ! r12 q2 = r2*rf2-q2 *qf2 ! q12 r2 = r1*qf1+q1 *rf1 ! q11 r1 = r1*rf1-q1 *qf1 ! r11 q1 = r2 +t1 ! q21 r2 =(r2 -t1)*qep ! r22 t1 = r1 +q2 ! r21 r1 =(r1 -q2)*qep ! q22 r(k0) = r(k0)+t1 ! r40 q(k0) = q(k0)+q1 ! q40 t1 = r(k0)-t1 *rec ! r21 q1 = q(k0)-q1 *rec ! q21 q(k2) = q1 -r1 ! q42 q(k1) = q1 +r1 ! q41 r(k1) = t1 -r2 ! r41 r(k2) = t1 +r2 ! r42 enddo enddo ma=ma3 enddo if(ln5 > 0)then ! RADIX 5 nze=n/5; rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta do l=1,ln5 mb=mb/5; ma5=ma*5 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=rf1*rf2-qf1*qf2; qf3=rf1*qf2+qf1*rf2 ! f**3 rf4=rf2*rf2-qf2*qf2; qf4=2*rf2*qf2 ! f**4 do i=0,nm,ma5 k0=i+j; k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma r1=r(k1); r2=r(k2); r3=r(k3); r4=r(k4) q1=q(k1); q2=q(k2); q3=q(k3); q4=q(k4) t1 =r1*qf1+q1*rf1 ! q11 r1 =r1*rf1-q1*qf1 ! r11 q1 =r4*rf4-q4*qf4 ! r14 r4 =r4*qf4+q4*rf4 ! q14 q4 =r1 -q1 ! q24 r1 =r1 +q1 ! r21 q1 =t1 +r4 ! q21 r4 =t1 -r4 ! q24 t1 =r3*rf3-q3*qf3 ! r13 r3 =r3*qf3+q3*rf3 ! q13 q3 =r2*qf2+q2*rf2 ! q12 r2 =r2*rf2-q2*qf2 ! r12 q2 =q3 +r3 ! q22 r3 =q3 -r3 ! q23 q3 =r2 -t1 ! r23 r2 =r2 +t1 ! r22 r(k0) =r(k0)+r1 +r2 ! r40 q(k0) =q(k0)+q1 +q2 ! q40 t1 =r4*qze+r3*qet ! q34 r3 =r4*qet-r3*qze ! q33 r4 =r(k0)-r1*rec-r2*rzc ! r32 r1 =r(k0)-r1*rzc-r2*rec ! r31 r(k2) =r4 -r3 ! r42 r(k3) =r4 +r3 ! r43 r(k4) =r1 +t1 ! r44 r(k1) =r1 -t1 ! r41 t1 =q(k0)-q1*rzc-q2*rec ! q31 q2 =q(k0)-q1*rec-q2*rzc ! q32 q1 =q4*qet-q3*qze ! r33 q4 =q4*qze+q3*qet ! r34 q(k3) =q2 -q1 ! q43 q(k2) =q2 +q1 ! q42 q(k1) =t1 +q4 ! q41 q(k4) =t1 -q4 ! q44 enddo enddo ma=ma5 enddo endif if(ln7 > 0)then ! RADIX 7 nze=n/7 rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta rth=rze*ret-qze*qet; qth=rze*qet+qze*ret; rtc=one-rth ! <- theta do L=1,Ln7 mb=mb/7; ma7=ma*7 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2; jmb3=jmb*3 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=w(jmb3); qf3=w(nh+jmb3) ! f**3 rf4=rf2*rf2-qf2*qf2; qf4=2*rf2*qf2 ! f**4 rf5=rf2*rf3-qf2*qf3; qf5=rf2*qf3+qf2*rf3 ! f**5 rf6=rf3*rf3-qf3*qf3; qf6=2*rf3*qf3 ! f**6 do i=0,nm,ma7 k0=i+j k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma; k5=k4+ma; k6=k5+ma r0=r(k0) r1=r(k1);r2=r(k2);r3=r(k3);r4=r(k4);r5=r(k5);r6=r(k6) q0=q(k0) q1=q(k1);q2=q(k2);q3=q(k3);q4=q(k4);q5=q(k5);q6=q(k6) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=r6*rf6-q6*qf6 ! r16 r6=r6*qf6+q6*rf6 ! q16 q6=r1-q1 ! r26 r1=r1+q1 ! r21 q1=t1+r6 ! q21 r6=t1-r6 ! q26 t2=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=r5*rf5-q5*qf5 ! r15 r5=r5*qf5+q5*rf5 ! q15 q5=r2-q2 ! r25 r2=r2+q2 ! r22 q2=t2+r5 ! q22 r5=t2-r5 ! q25 t3=r3*qf3+q3*rf3 ! q13 r3=r3*rf3-q3*qf3 ! r13 q3=r4*rf4-q4*qf4 ! r14 r4=r4*qf4+q4*rf4 ! q14 q4=r3-q3 ! r24 r3=r3+q3 ! r23 q3=t3+r4 ! q23 r4=t3-r4 ! q24 r0=r0+r1+r2+r3 ! r40 q0=q0+q1+q2+q3 ! q40 t1= r6*qze+r5*qet+r4*qth ! q36 t2= r6*qet-r5*qth-r4*qze ! q35 t3= r6*qth-r5*qze+r4*qet ! q34 r6=r0-r1*rzc-r2*rec-r3*rtc ! r31 r5=r0-r1*rec-r2*rtc-r3*rzc ! r32 r4=r0-r1*rtc-r2*rzc-r3*rec ! r33 r(k0)=r0 ! r40 r(k1)=r6-t1 ! r41 r(k6)=r6+t1 ! r46 r(k2)=r5-t2 ! r42 r(k5)=r5+t2 ! r45 r(k3)=r4-t3 ! r43 r(k4)=r4+t3 ! r44 t1= q6*qze+q5*qet+q4*qth ! r36 t2= q6*qet-q5*qth-q4*qze ! r35 t3= q6*qth-q5*qze+q4*qet ! r34 q6=q0-q1*rzc-q2*rec-q3*rtc ! q31 q5=q0-q1*rec-q2*rtc-q3*rzc ! q32 q4=q0-q1*rtc-q2*rzc-q3*rec ! q33 q(k0)=q0 ! q40 q(k1)=q6+t1 ! q41 q(k6)=q6-t1 ! q46 q(k2)=q5+t2 ! q42 q(k5)=q5-t2 ! q45 q(k3)=q4+t3 ! q43 q(k4)=q4-t3 ! q44 enddo enddo ma=ma7 enddo endif if(lnb > 0)then ! RADIX 11 nze=n/11 rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta rth=rze*ret-qze*qet; qth=rze*qet+qze*ret; rtc=one-rth ! <- theta rio=ret*ret-qet*qet; qio=2*ret*qet; ric=one-rio ! <- iota rka=ret*rth-qet*qth; qka=ret*qth+qet*rth; rkc=one-rka ! <- kappa do L=1,Lnb mb=mb/11; mab=ma*11 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2; jmb3=jmb*3; jmb4=jmb*4; jmb5=jmb*5 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=w(jmb3); qf3=w(nh+jmb3) ! f**3 rf4=w(jmb4); qf4=w(nh+jmb4) ! f**4 rf5=w(jmb5); qf5=w(nh+jmb5) ! f**5 rf6=rf3*rf3-qf3*qf3; qf6=2*rf3*qf3 ! f**6 rf7=rf3*rf4-qf3*qf4; qf7=rf3*qf4+qf3*rf4 ! f**7 rf8=rf4*rf4-qf4*qf4; qf8=2*rf4*qf4 ! f**8 rf9=rf4*rf5-qf4*qf5; qf9=rf4*qf5+qf4*rf5 ! f**9 rfa=rf5*rf5-qf5*qf5; qfa=2*rf5*qf5 ! f**10 do i=0,nm,mab k0=i+j k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma; k5=k4+ma; k6=k5+ma k7=k6+ma; k8=k7+ma; k9=k8+ma; ka=k9+ma r0=r(k0) r1=r(k1);r2=r(k2);r3=r(k3);r4=r(k4);r5=r(k5);r6=r(k6) r7=r(k7);r8=r(k8);r9=r(k9);ra=r(ka) q0=q(k0) q1=q(k1);q2=q(k2);q3=q(k3);q4=q(k4);q5=q(k5);q6=q(k6) q7=q(k7);q8=q(k8);q9=q(k9);qa=q(ka) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=ra*rfa-qa*qfa ! r1a ra=ra*qfa+qa*rfa ! q1a qa=r1-q1 ! r2a r1=r1+q1 ! r21 q1=t1+ra ! q21 ra=t1-ra ! q2a t2=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=r9*rf9-q9*qf9 ! r19 r9=r9*qf9+q9*rf9 ! q19 q9=r2-q2 ! r29 r2=r2+q2 ! r22 q2=t2+r9 ! q22 r9=t2-r9 ! q29 t3=r3*qf3+q3*rf3 ! q13 r3=r3*rf3-q3*qf3 ! r13 q3=r8*rf8-q8*qf8 ! r18 r8=r8*qf8+q8*rf8 ! q18 q8=r3-q3 ! r28 r3=r3+q3 ! r23 q3=t3+r8 ! q23 r8=t3-r8 ! q28 t4=r4*qf4+q4*rf4 ! q14 r4=r4*rf4-q4*qf4 ! r14 q4=r7*rf7-q7*qf7 ! r17 r7=r7*qf7+q7*rf7 ! q17 q7=r4-q4 ! r27 r4=r4+q4 ! r24 q4=t4+r7 ! q24 r7=t4-r7 ! q27 t5=r5*qf5+q5*rf5 ! q15 r5=r5*rf5-q5*qf5 ! r15 q5=r6*rf6-q6*qf6 ! r16 r6=r6*qf6+q6*rf6 ! q16 q6=r5-q5 ! r26 r5=r5+q5 ! r25 q5=t5+r6 ! q25 r6=t5-r6 ! q26 r0=r0+r1+r2+r3+r4+r5 ! r40 q0=q0+q1+q2+q3+q4+q5 ! q40 t1= ra*qze+r9*qet+r8*qth+r7*qio+r6*qka ! q3a t2= ra*qet+r9*qio-r8*qka-r7*qth-r6*qze ! q39 t3= ra*qth-r9*qka-r8*qet+r7*qze+r6*qio ! q38 t4= ra*qio-r9*qth+r8*qze+r7*qka-r6*qet ! q37 t5= ra*qka-r9*qze+r8*qio-r7*qet+r6*qth ! q36 ra=r0-r1*rzc-r2*rec-r3*rtc-r4*ric-r5*rkc ! r31 r9=r0-r1*rec-r2*ric-r3*rkc-r4*rtc-r5*rzc ! r32 r8=r0-r1*rtc-r2*rkc-r3*rec-r4*rzc-r5*ric ! r33 r7=r0-r1*ric-r2*rtc-r3*rzc-r4*rkc-r5*rec ! r34 r6=r0-r1*rkc-r2*rzc-r3*ric-r4*rec-r5*rtc ! r35 r(k0)=r0 ! r40 r(k1)=ra-t1 ! r41 r(ka)=ra+t1 ! r4a r(k2)=r9-t2 ! r42 r(k9)=r9+t2 ! r49 r(k3)=r8-t3 ! r43 r(k8)=r8+t3 ! r48 r(k4)=r7-t4 ! r44 r(k7)=r7+t4 ! r47 r(k5)=r6-t5 ! r45 r(k6)=r6+t5 ! r46 t1= qa*qze+q9*qet+q8*qth+q7*qio+q6*qka ! r3a t2= qa*qet+q9*qio-q8*qka-q7*qth-q6*qze ! r39 t3= qa*qth-q9*qka-q8*qet+q7*qze+q6*qio ! r38 t4= qa*qio-q9*qth+q8*qze+q7*qka-q6*qet ! r37 t5= qa*qka-q9*qze+q8*qio-q7*qet+q6*qth ! r36 qa=q0-q1*rzc-q2*rec-q3*rtc-q4*ric-q5*rkc ! q31 q9=q0-q1*rec-q2*ric-q3*rkc-q4*rtc-q5*rzc ! q32 q8=q0-q1*rtc-q2*rkc-q3*rec-q4*rzc-q5*ric ! q33 q7=q0-q1*ric-q2*rtc-q3*rzc-q4*rkc-q5*rec ! q34 q6=q0-q1*rkc-q2*rzc-q3*ric-q4*rec-q5*rtc ! q35 q(k0)=q0 ! q40 q(k1)=qa+t1 ! q41 q(ka)=qa-t1 ! q4a q(k2)=q9+t2 ! q42 q(k9)=q9-t2 ! q49 q(k3)=q8+t3 ! q43 q(k8)=q8-t3 ! q48 q(k4)=q7+t4 ! q44 q(k7)=q7-t4 ! q47 q(k5)=q6+t5 ! q45 q(k6)=q6-t5 ! q46 enddo enddo ma=mab enddo endif if(lnd > 0)then ! RADIX 13 nze=n/13 rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta rth=rze*ret-qze*qet; qth=rze*qet+qze*ret; rtc=one-rth ! <- theta rio=ret*ret-qet*qet; qio=2*ret*qet; ric=one-rio ! <- iota rka=ret*rth-qet*qth; qka=ret*qth+qet*rth; rkc=one-rka ! <- kappa rla=rth*rth-qth*qth; qla=2*rth*qth; rlc=one-rla ! <- lambda do L=1,Lnd mb=mb/13; mad=ma*13 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2; jmb3=jmb*3; jmb4=jmb*4; jmb5=jmb*5; jmb6=jmb*6 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=w(jmb3); qf3=w(nh+jmb3) ! f**3 rf4=w(jmb4); qf4=w(nh+jmb4) ! f**4 rf5=w(jmb5); qf5=w(nh+jmb5) ! f**5 rf6=w(jmb6); qf6=w(nh+jmb6) ! f**6 rf7=rf3*rf4-qf3*qf4; qf7=rf3*qf4+qf3*rf4 ! f**7 rf8=rf4*rf4-qf4*qf4; qf8=2*rf4*qf4 ! f**8 rf9=rf4*rf5-qf4*qf5; qf9=rf4*qf5+qf4*rf5 ! f**9 rfa=rf5*rf5-qf5*qf5; qfa=2*rf5*qf5 ! f**10 rfb=rf5*rf6-qf5*qf6; qfb=rf5*qf6+qf5*rf6 ! f**11 rfc=rf6*rf6-qf6*qf6; qfc=2*rf6*qf6 ! f**12 do i=0,nm,mad k0=i+j k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma; k5=k4+ma; k6=k5+ma k7=k6+ma; k8=k7+ma; k9=k8+ma; ka=k9+ma; kb=ka+ma; kc=kb+ma r0=r(k0) r1=r(k1);r2=r(k2);r3=r(k3);r4=r(k4);r5=r(k5);r6=r(k6) r7=r(k7);r8=r(k8);r9=r(k9);ra=r(ka);rb=r(kb);rc=r(kc) q0=q(k0) q1=q(k1);q2=q(k2);q3=q(k3);q4=q(k4);q5=q(k5);q6=q(k6) q7=q(k7);q8=q(k8);q9=q(k9);qa=q(ka);qb=q(kb);qc=q(kc) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=rc*rfc-qc*qfc ! r1c rc=rc*qfc+qc*rfc ! q1c qc=r1-q1 ! r2c r1=r1+q1 ! r21 q1=t1+rc ! q21 rc=t1-rc ! q2c t2=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=rb*rfb-qb*qfb ! r1b rb=rb*qfb+qb*rfb ! q1b qb=r2-q2 ! r2b r2=r2+q2 ! r22 q2=t2+rb ! q22 rb=t2-rb ! q2b t3=r3*qf3+q3*rf3 ! q13 r3=r3*rf3-q3*qf3 ! r13 q3=ra*rfa-qa*qfa ! r1a ra=ra*qfa+qa*rfa ! q1a qa=r3-q3 ! r2a r3=r3+q3 ! r23 q3=t3+ra ! q23 ra=t3-ra ! q2a t4=r4*qf4+q4*rf4 ! q14 r4=r4*rf4-q4*qf4 ! r14 q4=r9*rf9-q9*qf9 ! r19 r9=r9*qf9+q9*rf9 ! q19 q9=r4-q4 ! r29 r4=r4+q4 ! r24 q4=t4+r9 ! q24 r9=t4-r9 ! q29 t5=r5*qf5+q5*rf5 ! q15 r5=r5*rf5-q5*qf5 ! r15 q5=r8*rf8-q8*qf8 ! r18 r8=r8*qf8+q8*rf8 ! q18 q8=r5-q5 ! r28 r5=r5+q5 ! r25 q5=t5+r8 ! q25 r8=t5-r8 ! q28 t6=r6*qf6+q6*rf6 ! q16 r6=r6*rf6-q6*qf6 ! r16 q6=r7*rf7-q7*qf7 ! r17 r7=r7*qf7+q7*rf7 ! q17 q7=r6-q6 ! r27 r6=r6+q6 ! r26 q6=t6+r7 ! q26 r7=t6-r7 ! q27 r0=r0+r1+r2+r3+r4+r5+r6 ! r40 q0=q0+q1+q2+q3+q4+q5+q6 ! q40 t1= rc*qze+rb*qet+ra*qth+r9*qio+r8*qka+r7*qla ! q3c t2= rc*qet+rb*qio+ra*qla-r9*qka-r8*qth-r7*qze ! q3b t3= rc*qth+rb*qla-ra*qio-r9*qze+r8*qet+r7*qka ! q3a t4= rc*qio-rb*qka-ra*qze+r9*qth-r8*qla-r7*qet ! q39 t5= rc*qka-rb*qth+ra*qet-r9*qla-r8*qze+r7*qio ! q38 t6= rc*qla-rb*qze+ra*qka-r9*qet+r8*qio-r7*qth ! q37 rc=r0-r1*rzc-r2*rec-r3*rtc-r4*ric-r5*rkc-r6*rlc ! r31 rb=r0-r1*rec-r2*ric-r3*rlc-r4*rkc-r5*rtc-r6*rzc ! r32 ra=r0-r1*rtc-r2*rlc-r3*ric-r4*rzc-r5*rec-r6*rkc ! r33 r9=r0-r1*ric-r2*rkc-r3*rzc-r4*rtc-r5*rlc-r6*rec ! r34 r8=r0-r1*rkc-r2*rtc-r3*rec-r4*rlc-r5*rzc-r6*ric ! r35 r7=r0-r1*rlc-r2*rzc-r3*rkc-r4*rec-r5*ric-r6*rtc ! r36 r(k0)=r0 ! r40 r(k1)=rc-t1 ! r41 r(kc)=rc+t1 ! r4c r(k2)=rb-t2 ! r42 r(kb)=rb+t2 ! r4b r(k3)=ra-t3 ! r43 r(ka)=ra+t3 ! r4a r(k4)=r9-t4 ! r44 r(k9)=r9+t4 ! r49 r(k5)=r8-t5 ! r45 r(k8)=r8+t5 ! r48 r(k6)=r7-t6 ! r46 r(k7)=r7+t6 ! r47 t1= qc*qze+qb*qet+qa*qth+q9*qio+q8*qka+q7*qla ! r3c t2= qc*qet+qb*qio+qa*qla-q9*qka-q8*qth-q7*qze ! r3b t3= qc*qth+qb*qla-qa*qio-q9*qze+q8*qet+q7*qka ! r3a t4= qc*qio-qb*qka-qa*qze+q9*qth-q8*qla-q7*qet ! r39 t5= qc*qka-qb*qth+qa*qet-q9*qla-q8*qze+q7*qio ! r38 t6= qc*qla-qb*qze+qa*qka-q9*qet+q8*qio-q7*qth ! r37 qc=q0-q1*rzc-q2*rec-q3*rtc-q4*ric-q5*rkc-q6*rlc ! q31 qb=q0-q1*rec-q2*ric-q3*rlc-q4*rkc-q5*rtc-q6*rzc ! q32 qa=q0-q1*rtc-q2*rlc-q3*ric-q4*rzc-q5*rec-q6*rkc ! q33 q9=q0-q1*ric-q2*rkc-q3*rzc-q4*rtc-q5*rlc-q6*rec ! q34 q8=q0-q1*rkc-q2*rtc-q3*rec-q4*rlc-q5*rzc-q6*ric ! q35 q7=q0-q1*rlc-q2*rzc-q3*rkc-q4*rec-q5*ric-q6*rtc ! q36 q(k0)=q0 ! q40 q(k1)=qc+t1 ! q41 q(kc)=qc-t1 ! q4c q(k2)=qb+t2 ! q42 q(kb)=qb-t2 ! q4b q(k3)=qa+t3 ! q43 q(ka)=qa-t3 ! q4a q(k4)=q9+t4 ! q44 q(k9)=q9-t4 ! q49 q(k5)=q8+t5 ! q45 q(k8)=q8-t5 ! q48 q(k6)=q7+t6 ! q46 q(k7)=q7-t6 ! q47 enddo enddo ma=mad enddo endif end subroutine sgfft !============================================================================= subroutine dgfft(n,jumble,w,ln2,ln3,ln4,ln5,ln7,lnb,lnd,r,q)! [gfft] !============================================================================= use pietc, only: one,half integer, intent(IN ):: n integer, dimension(0:n-1),intent(IN ):: jumble real(dp),dimension(0:n-1),intent(IN ):: w integer, intent(IN ):: ln2,ln3,ln4,ln5,ln7,lnb,lnd real(dp),dimension(0:n-1),intent(INOUT):: r,q !----------------------------------------------------------------------------- integer :: nm,nh,i,j,l,k0,k1,k2,k3,k4,k5,k6,k7,k8,k9,ka,kb,kc, & ma,ma2,ma3,ma4,ma5,ma7,mab,mad,mb,mb2, & jmb,jmb2,jmb3,jmb4,jmb5,jmb6, & nze real(dp):: r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,ra,rb,rc, & q0,q1,q2,q3,q4,q5,q6,q7,q8,q9,qa,qb,qc, & rf1,rf2,rf3,rf4,rf5,rf6,rf7,rf8,rf9,rfa,rfb,rfc,& qf1,qf2,qf3,qf4,qf5,qf6,qf7,qf8,qf9,qfa,qfb,qfc, & rep,rze,rzc,ret,rec,rth,rtc,rio,ric,rka,rkc,rla,rlc, & qep,qze,qet,qth,qio,qka,qla, & t1,t2,t3,t4,t5,t6 !============================================================================= nm=n-1 nh=n/2 ! PERMUTE THE DATA: do i=1,nm j=jumble(i) if(j > i)then t1=r(i); r(i)=r(j); r(j)=t1 t1=q(i); q(i)=q(j); q(j)=t1 endif enddo ! TRANSFORM THE DATA: ma=1; mb=n ! RADIX 4 do l=1,ln4 mb=mb/4; ma4=ma*4; mb2=mb*2 do j=0,ma-1 jmb=j*mb; jmb2=j*mb2 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=rf1*rf2-qf1*qf2; qf3=rf1*qf2+qf1*rf2 ! f**3 do i=0,nm,ma4 k0=i+j; k1=k0+ma; k2=k1+ma; k3=k2+ma r0=r(k0); r1=r(k1); r2=r(k2); r3=r(k3) q0=q(k0); q1=q(k1); q2=q(k2); q3=q(k3) t1=r3*rf3-q3*qf3 ! q13 q3=r3*qf3+q3*rf3 ! r13 r3=r2*rf1-q2*qf1 ! r12 r2=r2*qf1+q2*rf1 ! q12 q2=q3 -r2 ! r23 r2=q3 +r2 ! q22 q3=r3 +t1 ! r22 t1=r3 -t1 ! q23 r3=r1*rf2-q1*qf2 ! r11 q1=r1*qf2+q1*rf2 ! q11 r1=r0 -r3 ! r21 r0=r0 +r3 ! r20 r(k3)=r1 -q2 ! r43 r(k1)=r1 +q2 ! r41 q2=q0 +q1 ! q20 q1=q0 -q1 ! q21 q(k0)=q2 +r2 ! q40 q(k2)=q2 -r2 ! q42 r(k2)=r0 -q3 ! r42 r(k0)=r0 +q3 ! r40 q(k3)=q1 -t1 ! q43 q(k1)=q1 +t1 ! q41 enddo enddo ma=ma4 enddo if(ln2==1)then ! RADIX 2 mb=mb/2; ma2=ma*2 do j=0,ma-1 jmb=j*mb do i=0,nm,ma2 k0=j+i; k1=k0+ma rf1=w(jmb); qf1=w(nh+jmb) ! f**1 r0=r(k0); q0=q(k0) r1=r(k1); q1=q(k1) t1=r1*qf1+q1*rf1 ! q11 q1=r1*rf1-q1*qf1 ! r11 r(k1)=r0 -q1 ! r41 r(k0)=r0 +q1 ! r40 q(k1)=q0 -t1 ! q41 q(k0)=q0 +t1 ! q40 enddo enddo ma=ma2 endif ! RADIX 3 rep=-half; rec=3*half; qep=half*sqrt(3*one) ! <- epsilon do l=1,ln3 mb=mb/3; ma3=ma*3 do j=0,ma-1 jmb=j*mb rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=rf1*rf1-qf1*qf1; qf2=2*rf1*qf1 ! f**2 do i=0,nm,ma3 k0=i+j; k1=k0+ma; k2=k1+ma r1=r(k1); q1=q(k1) r2=r(k2); q2=q(k2) t1= r2*qf2+q2 *rf2 ! r12 q2= r2*rf2-q2 *qf2 ! q12 r2= r1*qf1+q1 *rf1 ! q11 r1= r1*rf1-q1 *qf1 ! r11 q1= r2 +t1 ! q21 r2=(r2 -t1)*qep ! r22 t1= r1 +q2 ! r21 r1=(r1 -q2)*qep ! q22 r(k0)= r(k0)+t1 ! r40 q(k0)= q(k0)+q1 ! q40 t1= r(k0)-t1 *rec ! r21 q1= q(k0)-q1 *rec ! q21 q(k2)= q1 -r1 ! q42 q(k1)= q1 +r1 ! q41 r(k1)= t1 -r2 ! r41 r(k2)= t1 +r2 ! r42 enddo enddo ma=ma3 enddo if(ln5 > 0)then ! RADIX 5 nze=n/5; rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta do l=1,ln5 mb=mb/5; ma5=ma*5 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=rf1*rf2-qf1*qf2; qf3=rf1*qf2+qf1*rf2 ! f**3 rf4=rf2*rf2-qf2*qf2; qf4=2*rf2*qf2 ! f**4 do i=0,nm,ma5 k0=i+j; k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma r1=r(k1); r2=r(k2); r3=r(k3); r4=r(k4) q1=q(k1); q2=q(k2); q3=q(k3); q4=q(k4) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=r4*rf4-q4*qf4 ! r14 r4=r4*qf4+q4*rf4 ! q14 q4=r1 -q1 ! q24 r1=r1 +q1 ! r21 q1=t1 +r4 ! q21 r4=t1 -r4 ! q24 t1=r3*rf3-q3*qf3 ! r13 r3=r3*qf3+q3*rf3 ! q13 q3=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=q3 +r3 ! q22 r3=q3 -r3 ! q23 q3=r2 -t1 ! r23 r2=r2 +t1 ! r22 r(k0)=r(k0)+r1+r2 ! r40 q(k0)=q(k0)+q1+q2 ! q40 t1= r4*qze+r3*qet ! q34 r3= r4*qet-r3*qze ! q33 r4=r(k0)-r1*rec-r2*rzc ! r32 r1=r(k0)-r1*rzc-r2*rec ! r31 r(k2)=r4 -r3 ! r42 r(k3)=r4 +r3 ! r43 r(k4)=r1 +t1 ! r44 r(k1)=r1 -t1 ! r41 t1=q(k0)-q1*rzc-q2*rec ! q31 q2=q(k0)-q1*rec-q2*rzc ! q32 q1= q4*qet-q3*qze ! r33 q4= q4*qze+q3*qet ! r34 q(k3)=q2 -q1 ! q43 q(k2)=q2 +q1 ! q42 q(k1)=t1 +q4 ! q41 q(k4)=t1 -q4 ! q44 enddo enddo ma=ma5 enddo endif if(ln7 > 0)then ! RADIX 7 nze=n/7 rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta rth=rze*ret-qze*qet; qth=rze*qet+qze*ret; rtc=one-rth ! <- theta do L=1,Ln7 mb=mb/7; ma7=ma*7 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2; jmb3=jmb*3 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=w(jmb3); qf3=w(nh+jmb3) ! f**3 rf4=rf2*rf2-qf2*qf2; qf4=2*rf2*qf2 ! f**4 rf5=rf2*rf3-qf2*qf3; qf5=rf2*qf3+qf2*rf3 ! f**5 rf6=rf3*rf3-qf3*qf3; qf6=2*rf3*qf3 ! f**6 do i=0,nm,ma7 k0=i+j k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma; k5=k4+ma; k6=k5+ma r0=r(k0) r1=r(k1);r2=r(k2);r3=r(k3);r4=r(k4);r5=r(k5);r6=r(k6) q0=q(k0) q1=q(k1);q2=q(k2);q3=q(k3);q4=q(k4);q5=q(k5);q6=q(k6) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=r6*rf6-q6*qf6 ! r16 r6=r6*qf6+q6*rf6 ! q16 q6=r1-q1 ! r26 r1=r1+q1 ! r21 q1=t1+r6 ! q21 r6=t1-r6 ! q26 t2=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=r5*rf5-q5*qf5 ! r15 r5=r5*qf5+q5*rf5 ! q15 q5=r2-q2 ! r25 r2=r2+q2 ! r22 q2=t2+r5 ! q22 r5=t2-r5 ! q25 t3=r3*qf3+q3*rf3 ! q13 r3=r3*rf3-q3*qf3 ! r13 q3=r4*rf4-q4*qf4 ! r14 r4=r4*qf4+q4*rf4 ! q14 q4=r3-q3 ! r24 r3=r3+q3 ! r23 q3=t3+r4 ! q23 r4=t3-r4 ! q24 r0=r0+r1+r2+r3 ! r40 q0=q0+q1+q2+q3 ! q40 t1= r6*qze+r5*qet+r4*qth ! q36 t2= r6*qet-r5*qth-r4*qze ! q35 t3= r6*qth-r5*qze+r4*qet ! q34 r6=r0-r1*rzc-r2*rec-r3*rtc ! r31 r5=r0-r1*rec-r2*rtc-r3*rzc ! r32 r4=r0-r1*rtc-r2*rzc-r3*rec ! r33 r(k0)=r0 ! r40 r(k1)=r6-t1 ! r41 r(k6)=r6+t1 ! r46 r(k2)=r5-t2 ! r42 r(k5)=r5+t2 ! r45 r(k3)=r4-t3 ! r43 r(k4)=r4+t3 ! r44 t1= q6*qze+q5*qet+q4*qth ! r36 t2= q6*qet-q5*qth-q4*qze ! r35 t3= q6*qth-q5*qze+q4*qet ! r34 q6=q0-q1*rzc-q2*rec-q3*rtc ! q31 q5=q0-q1*rec-q2*rtc-q3*rzc ! q32 q4=q0-q1*rtc-q2*rzc-q3*rec ! q33 q(k0)=q0 ! q40 q(k1)=q6+t1 ! q41 q(k6)=q6-t1 ! q46 q(k2)=q5+t2 ! q42 q(k5)=q5-t2 ! q45 q(k3)=q4+t3 ! q43 q(k4)=q4-t3 ! q44 enddo enddo ma=ma7 enddo endif if(lnb > 0)then ! RADIX 11 nze=n/11 rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta rth=rze*ret-qze*qet; qth=rze*qet+qze*ret; rtc=one-rth ! <- theta rio=ret*ret-qet*qet; qio=2*ret*qet; ric=one-rio ! <- iota rka=ret*rth-qet*qth; qka=ret*qth+qet*rth; rkc=one-rka ! <- kappa do L=1,Lnb mb=mb/11; mab=ma*11 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2; jmb3=jmb*3; jmb4=jmb*4; jmb5=jmb*5 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=w(jmb3); qf3=w(nh+jmb3) ! f**3 rf4=w(jmb4); qf4=w(nh+jmb4) ! f**4 rf5=w(jmb5); qf5=w(nh+jmb5) ! f**5 rf6=rf3*rf3-qf3*qf3; qf6=2*rf3*qf3 ! f**6 rf7=rf3*rf4-qf3*qf4; qf7=rf3*qf4+qf3*rf4 ! f**7 rf8=rf4*rf4-qf4*qf4; qf8=2*rf4*qf4 ! f**8 rf9=rf4*rf5-qf4*qf5; qf9=rf4*qf5+qf4*rf5 ! f**9 rfa=rf5*rf5-qf5*qf5; qfa=2*rf5*qf5 ! f**10 do i=0,nm,mab k0=i+j k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma; k5=k4+ma; k6=k5+ma k7=k6+ma; k8=k7+ma; k9=k8+ma; ka=k9+ma r0=r(k0) r1=r(k1);r2=r(k2);r3=r(k3);r4=r(k4);r5=r(k5);r6=r(k6) r7=r(k7);r8=r(k8);r9=r(k9);ra=r(ka) q0=q(k0) q1=q(k1);q2=q(k2);q3=q(k3);q4=q(k4);q5=q(k5);q6=q(k6) q7=q(k7);q8=q(k8);q9=q(k9);qa=q(ka) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=ra*rfa-qa*qfa ! r1a ra=ra*qfa+qa*rfa ! q1a qa=r1-q1 ! r2a r1=r1+q1 ! r21 q1=t1+ra ! q21 ra=t1-ra ! q2a t2=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=r9*rf9-q9*qf9 ! r19 r9=r9*qf9+q9*rf9 ! q19 q9=r2-q2 ! r29 r2=r2+q2 ! r22 q2=t2+r9 ! q22 r9=t2-r9 ! q29 t3=r3*qf3+q3*rf3 ! q13 r3=r3*rf3-q3*qf3 ! r13 q3=r8*rf8-q8*qf8 ! r18 r8=r8*qf8+q8*rf8 ! q18 q8=r3-q3 ! r28 r3=r3+q3 ! r23 q3=t3+r8 ! q23 r8=t3-r8 ! q28 t4=r4*qf4+q4*rf4 ! q14 r4=r4*rf4-q4*qf4 ! r14 q4=r7*rf7-q7*qf7 ! r17 r7=r7*qf7+q7*rf7 ! q17 q7=r4-q4 ! r27 r4=r4+q4 ! r24 q4=t4+r7 ! q24 r7=t4-r7 ! q27 t5=r5*qf5+q5*rf5 ! q15 r5=r5*rf5-q5*qf5 ! r15 q5=r6*rf6-q6*qf6 ! r16 r6=r6*qf6+q6*rf6 ! q16 q6=r5-q5 ! r26 r5=r5+q5 ! r25 q5=t5+r6 ! q25 r6=t5-r6 ! q26 r0=r0+r1+r2+r3+r4+r5 ! r40 q0=q0+q1+q2+q3+q4+q5 ! q40 t1= ra*qze+r9*qet+r8*qth+r7*qio+r6*qka ! q3a t2= ra*qet+r9*qio-r8*qka-r7*qth-r6*qze ! q39 t3= ra*qth-r9*qka-r8*qet+r7*qze+r6*qio ! q38 t4= ra*qio-r9*qth+r8*qze+r7*qka-r6*qet ! q37 t5= ra*qka-r9*qze+r8*qio-r7*qet+r6*qth ! q36 ra=r0-r1*rzc-r2*rec-r3*rtc-r4*ric-r5*rkc ! r31 r9=r0-r1*rec-r2*ric-r3*rkc-r4*rtc-r5*rzc ! r32 r8=r0-r1*rtc-r2*rkc-r3*rec-r4*rzc-r5*ric ! r33 r7=r0-r1*ric-r2*rtc-r3*rzc-r4*rkc-r5*rec ! r34 r6=r0-r1*rkc-r2*rzc-r3*ric-r4*rec-r5*rtc ! r35 r(k0)=r0 ! r40 r(k1)=ra-t1 ! r41 r(ka)=ra+t1 ! r4a r(k2)=r9-t2 ! r42 r(k9)=r9+t2 ! r49 r(k3)=r8-t3 ! r43 r(k8)=r8+t3 ! r48 r(k4)=r7-t4 ! r44 r(k7)=r7+t4 ! r47 r(k5)=r6-t5 ! r45 r(k6)=r6+t5 ! r46 t1= qa*qze+q9*qet+q8*qth+q7*qio+q6*qka ! r3a t2= qa*qet+q9*qio-q8*qka-q7*qth-q6*qze ! r39 t3= qa*qth-q9*qka-q8*qet+q7*qze+q6*qio ! r38 t4= qa*qio-q9*qth+q8*qze+q7*qka-q6*qet ! r37 t5= qa*qka-q9*qze+q8*qio-q7*qet+q6*qth ! r36 qa=q0-q1*rzc-q2*rec-q3*rtc-q4*ric-q5*rkc ! q31 q9=q0-q1*rec-q2*ric-q3*rkc-q4*rtc-q5*rzc ! q32 q8=q0-q1*rtc-q2*rkc-q3*rec-q4*rzc-q5*ric ! q33 q7=q0-q1*ric-q2*rtc-q3*rzc-q4*rkc-q5*rec ! q34 q6=q0-q1*rkc-q2*rzc-q3*ric-q4*rec-q5*rtc ! q35 q(k0)=q0 ! q40 q(k1)=qa+t1 ! q41 q(ka)=qa-t1 ! q4a q(k2)=q9+t2 ! q42 q(k9)=q9-t2 ! q49 q(k3)=q8+t3 ! q43 q(k8)=q8-t3 ! q48 q(k4)=q7+t4 ! q44 q(k7)=q7-t4 ! q47 q(k5)=q6+t5 ! q45 q(k6)=q6-t5 ! q46 enddo enddo ma=mab enddo endif if(lnd > 0)then ! RADIX 13 nze=n/13 rze=w(nze); qze=w(nh+nze); rzc=one-rze ! <- zeta ret=rze*rze-qze*qze; qet=2*rze*qze; rec=one-ret ! <- eta rth=rze*ret-qze*qet; qth=rze*qet+qze*ret; rtc=one-rth ! <- theta rio=ret*ret-qet*qet; qio=2*ret*qet; ric=one-rio ! <- iota rka=ret*rth-qet*qth; qka=ret*qth+qet*rth; rkc=one-rka ! <- kappa rla=rth*rth-qth*qth; qla=2*rth*qth; rlc=one-rla ! <- lambda do L=1,Lnd mb=mb/13; mad=ma*13 do j=0,ma-1 jmb=j*mb; jmb2=jmb*2; jmb3=jmb*3; jmb4=jmb*4; jmb5=jmb*5; jmb6=jmb*6 rf1=w(jmb); qf1=w(nh+jmb) ! f**1 rf2=w(jmb2); qf2=w(nh+jmb2) ! f**2 rf3=w(jmb3); qf3=w(nh+jmb3) ! f**3 rf4=w(jmb4); qf4=w(nh+jmb4) ! f**4 rf5=w(jmb5); qf5=w(nh+jmb5) ! f**5 rf6=w(jmb6); qf6=w(nh+jmb6) ! f**6 rf7=rf3*rf4-qf3*qf4; qf7=rf3*qf4+qf3*rf4 ! f**7 rf8=rf4*rf4-qf4*qf4; qf8=2*rf4*qf4 ! f**8 rf9=rf4*rf5-qf4*qf5; qf9=rf4*qf5+qf4*rf5 ! f**9 rfa=rf5*rf5-qf5*qf5; qfa=2*rf5*qf5 ! f**10 rfb=rf5*rf6-qf5*qf6; qfb=rf5*qf6+qf5*rf6 ! f**11 rfc=rf6*rf6-qf6*qf6; qfc=2*rf6*qf6 ! f**12 do i=0,nm,mad k0=i+j k1=k0+ma; k2=k1+ma; k3=k2+ma; k4=k3+ma; k5=k4+ma; k6=k5+ma k7=k6+ma; k8=k7+ma; k9=k8+ma; ka=k9+ma; kb=ka+ma; kc=kb+ma r0=r(k0) r1=r(k1);r2=r(k2);r3=r(k3);r4=r(k4);r5=r(k5);r6=r(k6) r7=r(k7);r8=r(k8);r9=r(k9);ra=r(ka);rb=r(kb);rc=r(kc) q0=q(k0) q1=q(k1);q2=q(k2);q3=q(k3);q4=q(k4);q5=q(k5);q6=q(k6) q7=q(k7);q8=q(k8);q9=q(k9);qa=q(ka);qb=q(kb);qc=q(kc) t1=r1*qf1+q1*rf1 ! q11 r1=r1*rf1-q1*qf1 ! r11 q1=rc*rfc-qc*qfc ! r1c rc=rc*qfc+qc*rfc ! q1c qc=r1-q1 ! r2c r1=r1+q1 ! r21 q1=t1+rc ! q21 rc=t1-rc ! q2c t2=r2*qf2+q2*rf2 ! q12 r2=r2*rf2-q2*qf2 ! r12 q2=rb*rfb-qb*qfb ! r1b rb=rb*qfb+qb*rfb ! q1b qb=r2-q2 ! r2b r2=r2+q2 ! r22 q2=t2+rb ! q22 rb=t2-rb ! q2b t3=r3*qf3+q3*rf3 ! q13 r3=r3*rf3-q3*qf3 ! r13 q3=ra*rfa-qa*qfa ! r1a ra=ra*qfa+qa*rfa ! q1a qa=r3-q3 ! r2a r3=r3+q3 ! r23 q3=t3+ra ! q23 ra=t3-ra ! q2a t4=r4*qf4+q4*rf4 ! q14 r4=r4*rf4-q4*qf4 ! r14 q4=r9*rf9-q9*qf9 ! r19 r9=r9*qf9+q9*rf9 ! q19 q9=r4-q4 ! r29 r4=r4+q4 ! r24 q4=t4+r9 ! q24 r9=t4-r9 ! q29 t5=r5*qf5+q5*rf5 ! q15 r5=r5*rf5-q5*qf5 ! r15 q5=r8*rf8-q8*qf8 ! r18 r8=r8*qf8+q8*rf8 ! q18 q8=r5-q5 ! r28 r5=r5+q5 ! r25 q5=t5+r8 ! q25 r8=t5-r8 ! q28 t6=r6*qf6+q6*rf6 ! q16 r6=r6*rf6-q6*qf6 ! r16 q6=r7*rf7-q7*qf7 ! r17 r7=r7*qf7+q7*rf7 ! q17 q7=r6-q6 ! r27 r6=r6+q6 ! r26 q6=t6+r7 ! q26 r7=t6-r7 ! q27 r0=r0+r1+r2+r3+r4+r5+r6 ! r40 q0=q0+q1+q2+q3+q4+q5+q6 ! q40 t1= rc*qze+rb*qet+ra*qth+r9*qio+r8*qka+r7*qla ! q3c t2= rc*qet+rb*qio+ra*qla-r9*qka-r8*qth-r7*qze ! q3b t3= rc*qth+rb*qla-ra*qio-r9*qze+r8*qet+r7*qka ! q3a t4= rc*qio-rb*qka-ra*qze+r9*qth-r8*qla-r7*qet ! q39 t5= rc*qka-rb*qth+ra*qet-r9*qla-r8*qze+r7*qio ! q38 t6= rc*qla-rb*qze+ra*qka-r9*qet+r8*qio-r7*qth ! q37 rc=r0-r1*rzc-r2*rec-r3*rtc-r4*ric-r5*rkc-r6*rlc ! r31 rb=r0-r1*rec-r2*ric-r3*rlc-r4*rkc-r5*rtc-r6*rzc ! r32 ra=r0-r1*rtc-r2*rlc-r3*ric-r4*rzc-r5*rec-r6*rkc ! r33 r9=r0-r1*ric-r2*rkc-r3*rzc-r4*rtc-r5*rlc-r6*rec ! r34 r8=r0-r1*rkc-r2*rtc-r3*rec-r4*rlc-r5*rzc-r6*ric ! r35 r7=r0-r1*rlc-r2*rzc-r3*rkc-r4*rec-r5*ric-r6*rtc ! r36 r(k0)=r0 ! r40 r(k1)=rc-t1 ! r41 r(kc)=rc+t1 ! r4c r(k2)=rb-t2 ! r42 r(kb)=rb+t2 ! r4b r(k3)=ra-t3 ! r43 r(ka)=ra+t3 ! r4a r(k4)=r9-t4 ! r44 r(k9)=r9+t4 ! r49 r(k5)=r8-t5 ! r45 r(k8)=r8+t5 ! r48 r(k6)=r7-t6 ! r46 r(k7)=r7+t6 ! r47 t1= qc*qze+qb*qet+qa*qth+q9*qio+q8*qka+q7*qla ! r3c t2= qc*qet+qb*qio+qa*qla-q9*qka-q8*qth-q7*qze ! r3b t3= qc*qth+qb*qla-qa*qio-q9*qze+q8*qet+q7*qka ! r3a t4= qc*qio-qb*qka-qa*qze+q9*qth-q8*qla-q7*qet ! r39 t5= qc*qka-qb*qth+qa*qet-q9*qla-q8*qze+q7*qio ! r38 t6= qc*qla-qb*qze+qa*qka-q9*qet+q8*qio-q7*qth ! r37 qc=q0-q1*rzc-q2*rec-q3*rtc-q4*ric-q5*rkc-q6*rlc ! q31 qb=q0-q1*rec-q2*ric-q3*rlc-q4*rkc-q5*rtc-q6*rzc ! q32 qa=q0-q1*rtc-q2*rlc-q3*ric-q4*rzc-q5*rec-q6*rkc ! q33 q9=q0-q1*ric-q2*rkc-q3*rzc-q4*rtc-q5*rlc-q6*rec ! q34 q8=q0-q1*rkc-q2*rtc-q3*rec-q4*rlc-q5*rzc-q6*ric ! q35 q7=q0-q1*rlc-q2*rzc-q3*rkc-q4*rec-q5*ric-q6*rtc ! q36 q(k0)=q0 ! q40 q(k1)=qc+t1 ! q41 q(kc)=qc-t1 ! q4c q(k2)=qb+t2 ! q42 q(kb)=qb-t2 ! q4b q(k3)=qa+t3 ! q43 q(ka)=qa-t3 ! q4a q(k4)=q9+t4 ! q44 q(k9)=q9-t4 ! q49 q(k5)=q8+t5 ! q45 q(k8)=q8-t5 ! q48 q(k6)=q7+t6 ! q46 q(k7)=q7-t6 ! q47 enddo enddo ma=mad enddo endif end subroutine dgfft !============================================================================= subroutine scsft(n,rb,qb)! [csft] !============================================================================= ! SUBROUTINE SCSFT ! Complex slow Fourier transform. For inverse, see SDSFT. ! --> n: Period ! <-> rb,qb: Real and imaginary parts of data and transform. !============================================================================= use pietc, only: one,pi2 integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(INOUT):: rb,qb !----------------------------------------------------------------------------- integer :: nm,i,k,ik real(sp) :: dang,ang real(sp),dimension(0:n-1):: ra,qa,wc,ws !============================================================================= nm=n-1 dang=pi2/n do k=0,nm; ang=k*dang; wc(k)=cos(ang); ws(k)=sin(ang); enddo ra=0; qa=0 do k=0,nm do i=0,nm ik=mod(i*k,n) ra(k)=ra(k)+wc(ik)*rb(i)-ws(ik)*qb(i) qa(k)=qa(k)+wc(ik)*qb(i)+ws(ik)*rb(i) enddo enddo rb=ra; qb=qa end subroutine scsft !============================================================================= subroutine dcsft(n,rb,qb)! [csft] !============================================================================= ! SUBROUTINE DCSFT: Double precision version of csft !============================================================================= use pietc, only: pi2 integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(INOUT):: rb,qb !----------------------------------------------------------------------------- integer :: nm,i,k,ik real(dp) :: dang,ang real(dp),dimension(0:n-1):: ra,qa,wc,ws !============================================================================= nm=n-1 dang=pi2/n do k=0,nm; ang=k*dang; wc(k)=cos(ang); ws(k)=sin(ang); enddo ra=0; qa=0 do k=0,nm do i=0,nm ik=mod(i*k,n) ra(k)=ra(k)+wc(ik)*rb(i)-ws(ik)*qb(i) qa(k)=qa(k)+wc(ik)*qb(i)+ws(ik)*rb(i) enddo enddo rb=ra; qb=qa end subroutine dcsft !============================================================================= subroutine sdsft(n,rb,qb)! [dsft] !============================================================================= ! SUBROUTINE SDSFT ! Complex slow inverse Fourier transform. For direct transform, see CSFT. ! --> n: Period ! <-> rb,qb: Real and imaginary parts of transform and synthesis. !============================================================================= use pietc, only: one,pi2 integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(INOUT):: rb,qb !----------------------------------------------------------------------------- integer :: nm,i,k,ik real(sp) :: dang,ang,rfac real(sp),dimension(0:n-1):: ra,qa,wc,ws !============================================================================= nm=n-1 rfac=one/n dang=pi2/n do k=0,nm; ang=k*dang; wc(k)=cos(ang); ws(k)=sin(ang); enddo ra=0; qa=0 do k=0,nm do i=0,nm ik=mod(i*k,n) ra(k)=ra(k)+wc(ik)*rb(i)+ws(ik)*qb(i) qa(k)=qa(k)+wc(ik)*qb(i)-ws(ik)*rb(i) enddo enddo rb=ra*rfac; qb=qa*rfac end subroutine sdsft !============================================================================= subroutine ddsft(n,rb,qb)! [dsft] !============================================================================= ! SUBROUTINE DDSFT: Double precision version of dsft !============================================================================= use pietc, only: one,pi2 integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(INOUT):: rb,qb !----------------------------------------------------------------------------- integer :: nm,i,k,ik real(dp) :: dang,ang,rfac real(dp),dimension(0:n-1):: ra,qa,wc,ws !============================================================================= nm=n-1 rfac=one/n dang=pi2/n do k=0,nm; ang=k*dang; wc(k)=cos(ang); ws(k)=sin(ang); enddo ra=0; qa=0 do k=0,nm do i=0,nm ik=mod(i*k,n) ra(k)=ra(k)+wc(ik)*rb(i)+ws(ik)*qb(i) qa(k)=qa(k)+wc(ik)*qb(i)-ws(ik)*rb(i) enddo enddo rb=ra*rfac; qb=qa*rfac end subroutine ddsft !============================================================================= subroutine sdsfe(n,rb,qb,x,y,r,q,dr,dq)! [dsfe] !============================================================================= ! Single precision complex "slow" Fourier evaluation. ! Evaluate the (complex) Fourier transform at a single point, together with ! its complex derivative there. Assume the period of the series is 2*pi. ! ! --> n : number of imaginary Fourier coefficients of a complex function ! --> rb,qb : real and imaginary Fourier coefficients ! --> x,y : real and imaginary components of complex representation of point. ! <-- r,q : real and imaginary components of the function at the point. ! <-- dr,dq : real and imaginary components of derivative of function there. !============================================================================= use pietc, only: one integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(IN ):: rb,qb real(sp), intent(IN ):: x,y real(sp), intent(OUT):: r,q,dr,dq !----------------------------------------------------------------------------- real(sp) :: rfac,angx,angy,cangx,sangx,eangy,& rbk,qbk,rbl,qbl,rterm,qterm integer :: k,nm,nh,nhm,L !============================================================================ nm=n-1; nh=(n+1)/2; nhm=nh-1 r=0; q=0; dr=0; dq=0 rfac=one/n do k=0,nhm L=mod(n-k,n); rbk=rb(k); qbk=qb(k); rbl=rb(L); qbl=qb(L) angx=k*x; cangx=cos(angx); sangx=sin(angx); angy=k*y if(rbk==0 .and. qbk==0)goto 3 eangy=exp(angy) rterm=(rbk*cangx+qbk*sangx)*eangy; qterm=(-rbk*sangx+qbk*cangx)*eangy r =r +rterm; q =q +qterm; if(k==0)cycle dr=dr+qterm*k; dq=dq-rterm*k 3 if(rbl==0 .and. qbl==0)cycle eangy=exp(-angy) rterm=(rbl*cangx-qbl*sangx)*eangy; qterm=( rbl*sangx+qbl*cangx)*eangy r=r +rterm; q =q +qterm dr=dr-qterm*k; dq=dq+rterm*k enddo r=r*rfac; q=q*rfac; dr=dr*rfac; dq=dq*rfac end subroutine sdsfe !============================================================================= subroutine ddsfe(n,rb,qb,x,y,r,q,dr,dq)! [dsfe] !============================================================================= ! Double precision complex "slow" Fourier evaluation. ! Evaluate the (complex) Fourier transform at a single point, together with ! its complex derivative there. Assume the period of the series is 2*pi. ! ! --> n : number of imaginary Fourier coefficients of a complex function ! --> rb,qb : real and imaginary Fourier coefficients ! --> x,y : real and imaginary components of complex representation of point. ! <-- r,q : real and imaginary components of the function at the point. ! <-- dr,dq : real and imaginary components of derivative of function there. !============================================================================= use pietc, only: one integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(IN ):: rb,qb real(dp), intent(IN ):: x,y real(dp), intent(OUT):: r,q,dr,dq !----------------------------------------------------------------------------- real(dp) :: rfac,angx,angy,cangx,sangx,eangy,& rbk,qbk,rbl,qbl,rterm,qterm integer :: k,nm,nh,nhm,L !============================================================================ nm=n-1; nh=(n+1)/2; nhm=nh-1 r=0; q=0; dr=0; dq=0 rfac=one/n do k=0,nhm L=mod(n-k,n); rbk=rb(k); qbk=qb(k); rbl=rb(L); qbl=qb(L) angx=k*x; cangx=cos(angx); sangx=sin(angx); angy=k*y if(rbk==0 .and. qbk==0)goto 3 eangy=exp(angy) rterm=(rbk*cangx+qbk*sangx)*eangy; qterm=(-rbk*sangx+qbk*cangx)*eangy r =r +rterm; q =q +qterm; if(k==0)cycle dr=dr+qterm*k; dq=dq-rterm*k 3 if(rbl==0 .and. qbl==0)cycle eangy=exp(-angy) rterm=(rbl*cangx-qbl*sangx)*eangy; qterm=( rbl*sangx+qbl*cangx)*eangy r=r +rterm; q =q +qterm dr=dr-qterm*k; dq=dq+rterm*k enddo r=r*rfac; q=q*rfac; dr=dr*rfac; dq=dq*rfac end subroutine ddsfe !============================================================================ subroutine shsfe(n,b,x,r,dr)! [hsfe] !============================================================================ ! Slow Hermitian Fourier synthesis evaluation and derivative !============================================================================ use pietc, only: two integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(IN ):: b real(sp), intent(IN ):: x real(sp), intent(OUT):: r,dr !---------------------------------------------------------------------------- integer :: nm,nh,nhm,k,l real(sp) :: ang,cang,sang,bk,bl,rfac !============================================================================ nm=n-1; nh=(n+1)/2; nhm=nh-1 rfac=two/n r=(b(0)+b(nh)*cos(x*nh))/2 dr=0 do k=1,nhm l=n-k bk=b(k); bl=b(l) ang=x*k; cang=cos(ang); sang=sin(ang) r=r+bk*cang-bl*sang dr=dr-(bk*sang+bl*cang)*k enddo r=r*rfac; dr=dr*rfac end subroutine shsfe ! !============================================================================ subroutine dhsfe(n,b,x,r,dr)! [hsfe] !============================================================================ use pietc, only: two integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(IN ):: b real(dp), intent(IN ):: x real(dp), intent(OUT):: r,dr !---------------------------------------------------------------------------- integer :: nm,nh,nhm,k,l real(dp) :: ang,cang,sang,bk,bl,rfac !============================================================================ nm=n-1; nh=(n+1)/2; nhm=nh-1 rfac=two/n r=(b(0)+b(nh)*cos(x*nh))/2 dr=0 do k=1,nhm l=n-k bk=b(k); bl=b(l) ang=x*k; cang=cos(ang); sang=sin(ang) r=r+bk*cang-bl*sang dr=dr-(bk*sang+bl*cang)*k enddo r=r*rfac; dr=dr*rfac end subroutine dhsfe !============================================================================= subroutine srfft(n,r)! [rfft] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! (Revised for f90, 1999) ! SUBROUTINES SRFFT ! Fourier analyze (RFFT) a line of real data (for inverse, see HFFT) ! The real coefficients of the transforms are stored in direct order ! of wavenumber in elements 0 through n/2. The imaginary coefficients ! of the transform are stored in reverse order of wavenumber in ! the remaining elements. In this way, the structure of the transform ! of real data is made to conform more immediately to that of the ! complex data transform (provided in cfft). ! ! --> N Number of data of series (product of 2s,3s,5s,7s,11s,13s) ! <-> R Data and transform !============================================================================= integer, intent(IN ):: n real(sp), dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nhp real(sp), dimension(0:n-1) :: t !============================================================================= nm=n-1; nhp=n/2+1 t=0 call cfft(n,r,t) r(nhp:nm)=t(nhp:nm) end subroutine srfft !============================================================================= subroutine drfft(n,r)! [rfft] !============================================================================= integer, intent(IN ):: n real(dp), dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nhp real(dp), dimension(0:n-1) :: t !============================================================================= nm=n-1; nhp=n/2+1 t=0 call cfft(n,r,t) r(nhp:nm)=t(nhp:nm) end subroutine drfft !============================================================================= subroutine srfftp(n,js,ws,r)! [rfft] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(sp),dimension(0:n-1),intent(INOUT):: ws real(sp),dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nhp real(sp), dimension(0:n-1) :: t !============================================================================= nm=n-1; nhp=n/2+1 t=0 call cfft(n,js,ws,r,t) r(nhp:nm)=t(nhp:nm) end subroutine srfftp !============================================================================= subroutine drfftp(n,js,ws,r)! [rfft] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(dp),dimension(0:n-1),intent(INOUT):: ws real(dp),dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nhp real(dp), dimension(0:n-1) :: t !============================================================================= nm=n-1; nhp=n/2+1 t=0 call cfft(n,js,ws,r,t) r(nhp:nm)=t(nhp:nm) end subroutine drfftp !============================================================================= subroutine shfft(n,r)! [hfft] !============================================================================= ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 ! (Revised for f90, 1999) ! SUBROUTINES SHFFT ! Fourier synthesize (HFFT) a line of real data ! ! --> N Number of data of series (product of 2s,3s,5s,7s,11s,13s) ! <-> R Data and transform !============================================================================= integer, intent(IN ):: n real(sp), dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nh,nmh,nhp real(sp), dimension(0:n-1) :: t !============================================================================= nm=n-1; nh=n/2; nmh=nm/2; nhp=nh+1 t(0)=0; t(nh)=0 t(nhp:nm)=r(nhp:nm); t(1:nmh)=-t(nm:nhp:-1); r(nhp:nm)=r(nmh:1:-1) call dfft(n,r,t) end subroutine shfft !============================================================================= subroutine dhfft(n,r)! [hfft] !============================================================================= integer, intent(IN ):: n real(dp), dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nh,nmh,nhp real(dp),dimension(0:n-1) :: t !============================================================================= nm=n-1; nh=n/2; nmh=nm/2; nhp=nh+1 t(0)=0; t(nh)=0 t(nhp:nm)=r(nhp:nm); t(1:nmh)=-t(nm:nhp:-1); r(nhp:nm)=r(nmh:1:-1) call dfft(n,r,t) end subroutine dhfft !============================================================================= subroutine shfftp(n,js,ws,r)! [hfft] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(sp),dimension(0:n-1),intent(INOUT):: ws real(sp),dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nh,nmh,nhp real(sp),dimension(0:n-1) :: t !============================================================================= nm=n-1; nh=n/2; nmh=nm/2; nhp=nh+1 t(0)=0; t(nh)=0 t(nhp:nm)=r(nhp:nm); t(1:nmh)=-t(nm:nhp:-1); r(nhp:nm)=r(nmh:1:-1) call dfft(n,js,ws,r,t) end subroutine shfftp !============================================================================= subroutine dhfftp(n,js,ws,r)! [hfft] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(dp),dimension(0:n-1),intent(INOUT):: ws real(dp),dimension(0:n-1),intent(INOUT):: r !----------------------------------------------------------------------------- integer :: nm,nh,nmh,nhp real(dp),dimension(0:n-1) :: t !============================================================================= nm=n-1; nh=n/2; nmh=nm/2; nhp=nh+1 t(0)=0; t(nh)=0 t(nhp:nm)=r(nhp:nm); t(1:nmh)=-t(nm:nhp:-1); r(nhp:nm)=r(nmh:1:-1) call dfft(n,js,ws,r,t) end subroutine dhfftp !============================================================================= subroutine sfftcnv(n,a,b,d)! [fftcnv] !============================================================================= integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(IN ):: a,b real(sp),dimension(0:n-1),intent(OUT):: d real(sp),dimension(0:n-1) :: ws integer, dimension(0:n-1) :: js !============================================================================= js(0)=0 call fftcnv(n,js,ws,a,b,d) end subroutine sfftcnv !============================================================================= subroutine dfftcnv(n,a,b,d)! [fftcnv] !============================================================================= integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(IN ):: a,b real(dp),dimension(0:n-1),intent(OUT):: d real(dp),dimension(0:n-1) :: ws integer, dimension(0:n-1) :: js !============================================================================= js(0)=0 call fftcnv(n,js,ws,a,b,d) end subroutine dfftcnv !============================================================================= subroutine sfftcnvp(n,js,ws,a,b,d)! [fftcnv] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(sp),dimension(0:n-1),intent(INOUT):: ws real(sp),dimension(0:n-1),intent(IN ):: a,b real(sp),dimension(0:n-1),intent( OUT):: d real(sp),dimension(0:n-1) :: at,bt !============================================================================= at=a bt=b call rfft(n,js,ws,at) call rfft(n,js,ws,bt) call fconv(n,at,bt,d) call hfft(n,js,ws,d) end subroutine sfftcnvp !============================================================================= subroutine dfftcnvp(n,js,ws,a,b,d)! [fftcnv] !============================================================================= integer, intent(IN ):: n integer, dimension(0:n-1),intent(INOUT):: js real(dp),dimension(0:n-1),intent(INOUT):: ws real(dp),dimension(0:n-1),intent(IN ):: a,b real(dp),dimension(0:n-1),intent( OUT):: d real(dp),dimension(0:n-1) :: at,bt !============================================================================= at=a bt=b call rfft(n,js,ws,at) call rfft(n,js,ws,bt) call fconv(n,at,bt,d) call hfft(n,js,ws,d) end subroutine dfftcnvp !============================================================================= subroutine sfconv(n,a,b,c)! [fconv] !============================================================================= ! Assuming Hermetian representations, A and B, of period n, form their ! product, C, that represents the transform of the convolution of the a and b ! whose transforms are aforementioned A and B. !============================================================================= integer, intent(IN ):: n real(sp),dimension(0:n-1),intent(IN ):: a,b real(sp),dimension(0:n-1),intent(OUT):: c !----------------------------------------------------------------------------- integer :: i,ic,nmnh,nh !============================================================================= nh=n/2 nmnh=n-nh ! =nh for n even; =nh+1 for n odd do i=0,nh,nmnh ! <- trick index stride and limit allows for no 2-wave at n odd c(i)=a(i)*b(i) enddo do i=1,nmnh-1; ic=n-i c(i )=a(i)*b(i )-a(ic)*b(ic) c(ic)=a(i)*b(ic)+a(ic)*b(i ) enddo end subroutine sfconv !============================================================================= subroutine dfconv(n,a,b,c)! [fconv] !============================================================================= integer, intent(IN ):: n real(dp),dimension(0:n-1),intent(IN ):: a,b real(dp),dimension(0:n-1),intent(OUT):: c !----------------------------------------------------------------------------- integer :: i,ic,nmnh,nh !============================================================================= nh=n/2 nmnh=n-nh ! =nh for n even; =nh+1 for n odd do i=0,nh,nmnh ! <- trick index stride and limit allows for no 2-wave at n odd c(i)=a(i)*b(i) enddo do i=1,nmnh-1; ic=n-i c(i )=a(i)*b(i )-a(ic)*b(ic) c(ic)=a(i)*b(ic)+a(ic)*b(i ) enddo end subroutine dfconv end module pfft