Skip to content

Refactor #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 11 commits into from
Jan 27, 2022
Merged
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 0 additions & 62 deletions src/cfftb1.f

This file was deleted.

68 changes: 68 additions & 0 deletions src/cfftb1.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
subroutine cfftb1(n,c,Ch,Wa,Ifac)
use fftpack_kind
implicit none
real(rk) :: c , Ch , Wa
integer :: i , idl1 , ido , idot , Ifac , ip , iw , ix2 , ix3 , ix4, &
k1 , l1 , l2 , n , n2 , na , nac , nf
dimension Ch(*) , c(*) , Wa(*) , Ifac(*)
nf = Ifac(2)
na = 0
l1 = 1
iw = 1
do k1 = 1 , nf
ip = Ifac(k1+2)
l2 = ip*l1
ido = n/l2
idot = ido + ido
idl1 = idot*l1
if ( ip==4 ) then
ix2 = iw + idot
ix3 = ix2 + idot
if ( na/=0 ) then
call passb4(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3))
else
call passb4(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3))
endif
na = 1 - na
elseif ( ip==2 ) then
if ( na/=0 ) then
call passb2(idot,l1,Ch,c,Wa(iw))
else
call passb2(idot,l1,c,Ch,Wa(iw))
endif
na = 1 - na
elseif ( ip==3 ) then
ix2 = iw + idot
if ( na/=0 ) then
call passb3(idot,l1,Ch,c,Wa(iw),Wa(ix2))
else
call passb3(idot,l1,c,Ch,Wa(iw),Wa(ix2))
endif
na = 1 - na
elseif ( ip/=5 ) then
if ( na/=0 ) then
call passb(nac,idot,ip,l1,idl1,Ch,Ch,Ch,c,c,Wa(iw))
else
call passb(nac,idot,ip,l1,idl1,c,c,c,Ch,Ch,Wa(iw))
endif
if ( nac/=0 ) na = 1 - na
else
ix2 = iw + idot
ix3 = ix2 + idot
ix4 = ix3 + idot
if ( na/=0 ) then
call passb5(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4))
else
call passb5(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4))
endif
na = 1 - na
endif
l1 = l2
iw = iw + (ip-1)*idot
enddo
if ( na==0 ) return
n2 = n + n
do i = 1 , n2
c(i) = Ch(i)
enddo
end subroutine cfftb1
62 changes: 0 additions & 62 deletions src/cfftf1.f

This file was deleted.

68 changes: 68 additions & 0 deletions src/cfftf1.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
subroutine cfftf1(n,c,Ch,Wa,Ifac)
use fftpack_kind
implicit none
real(rk) :: c , Ch , Wa
integer :: i , idl1 , ido , idot , Ifac , ip , iw , ix2 , ix3 , ix4, &
k1 , l1 , l2 , n , n2 , na , nac , nf
dimension Ch(*) , c(*) , Wa(*) , Ifac(*)
nf = Ifac(2)
na = 0
l1 = 1
iw = 1
do k1 = 1 , nf
ip = Ifac(k1+2)
l2 = ip*l1
ido = n/l2
idot = ido + ido
idl1 = idot*l1
if ( ip==4 ) then
ix2 = iw + idot
ix3 = ix2 + idot
if ( na/=0 ) then
call passf4(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3))
else
call passf4(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3))
endif
na = 1 - na
elseif ( ip==2 ) then
if ( na/=0 ) then
call passf2(idot,l1,Ch,c,Wa(iw))
else
call passf2(idot,l1,c,Ch,Wa(iw))
endif
na = 1 - na
elseif ( ip==3 ) then
ix2 = iw + idot
if ( na/=0 ) then
call passf3(idot,l1,Ch,c,Wa(iw),Wa(ix2))
else
call passf3(idot,l1,c,Ch,Wa(iw),Wa(ix2))
endif
na = 1 - na
elseif ( ip/=5 ) then
if ( na/=0 ) then
call passf(nac,idot,ip,l1,idl1,Ch,Ch,Ch,c,c,Wa(iw))
else
call passf(nac,idot,ip,l1,idl1,c,c,c,Ch,Ch,Wa(iw))
endif
if ( nac/=0 ) na = 1 - na
else
ix2 = iw + idot
ix3 = ix2 + idot
ix4 = ix3 + idot
if ( na/=0 ) then
call passf5(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4))
else
call passf5(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4))
endif
na = 1 - na
endif
l1 = l2
iw = iw + (ip-1)*idot
enddo
if ( na==0 ) return
n2 = n + n
do i = 1 , n2
c(i) = Ch(i)
enddo
end subroutine cfftf1
61 changes: 0 additions & 61 deletions src/cffti1.f

This file was deleted.

68 changes: 68 additions & 0 deletions src/cffti1.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
subroutine cffti1(n,Wa,Ifac)
use fftpack_kind
implicit none
real(rk) :: arg , argh , argld , fi , Wa
integer :: i , i1 , ib , ido , idot , Ifac , ii , ip , ipm , j , k1, &
l1 , l2 , ld , n , nf , nl , nq , nr , ntry
dimension Wa(*) , Ifac(*)
integer,dimension(4),parameter :: ntryh = [3 , 4 , 2 , 5]
real(rk),parameter :: tpi = 2.0_rk * acos(-1.0_rk) ! 2 * pi
nl = n
nf = 0
j = 0
100 j = j + 1
if ( j<=4 ) then
ntry = ntryh(j)
else
ntry = ntry + 2
endif
200 nq = nl/ntry
nr = nl - ntry*nq
if ( nr/=0 ) goto 100
nf = nf + 1
Ifac(nf+2) = ntry
nl = nq
if ( ntry==2 ) then
if ( nf/=1 ) then
do i = 2 , nf
ib = nf - i + 2
Ifac(ib+2) = Ifac(ib+1)
enddo
Ifac(3) = 2
endif
endif
if ( nl/=1 ) goto 200
Ifac(1) = n
Ifac(2) = nf
argh = tpi/real(n, rk)
i = 2
l1 = 1
do k1 = 1 , nf
ip = Ifac(k1+2)
ld = 0
l2 = l1*ip
ido = n/l2
idot = ido + ido + 2
ipm = ip - 1
do j = 1 , ipm
i1 = i
Wa(i-1) = 1.0_rk
Wa(i) = 0.0_rk
ld = ld + l1
fi = 0.0_rk
argld = real(ld, rk)*argh
do ii = 4 , idot , 2
i = i + 2
fi = fi + 1.0_rk
arg = fi*argld
Wa(i-1) = cos(arg)
Wa(i) = sin(arg)
enddo
if ( ip>5 ) then
Wa(i1-1) = Wa(i-1)
Wa(i1) = Wa(i)
endif
enddo
l1 = l2
enddo
end subroutine cffti1
Loading