diff --git a/src/cfftb1.f90 b/src/cfftb1.f90 index 43f5b54..7c0bf79 100644 --- a/src/cfftb1.f90 +++ b/src/cfftb1.f90 @@ -1,68 +1,69 @@ - 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)) + subroutine cfftb1(n, c, ch, wa, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n, ifac(*) + real(dp), intent(in) :: wa(*) + real(dp), intent(inout) :: c(*), ch(*) + integer :: i, idl1, ido, idot, ip, iw, ix2, ix3, ix4, & + k1, l1, l2, n2, na, nac, nf + 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)) + end if + 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)) + end if + 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)) + end if + 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)) + end if + if (nac /= 0) na = 1 - na 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 \ No newline at end of file + 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)) + end if + na = 1 - na + end if + l1 = l2 + iw = iw + (ip - 1)*idot + end do + if (na == 0) return + n2 = n + n + do concurrent(i=1:n2) + c(i) = ch(i) + end do + end subroutine cfftb1 diff --git a/src/cfftf1.f90 b/src/cfftf1.f90 index 0139f39..f692e3a 100644 --- a/src/cfftf1.f90 +++ b/src/cfftf1.f90 @@ -1,68 +1,69 @@ - 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)) + subroutine cfftf1(n, c, ch, wa, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n, ifac(*) + real(dp), intent(inout) :: c(*), ch(*) + real(dp), intent(in) :: wa(*) + integer :: i, idl1, ido, idot, ip, iw, ix2, ix3, ix4, & + k1, l1, l2, n2, na, nac, nf + 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)) + end if + 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)) + end if + 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)) + end if + 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)) + end if + if (nac /= 0) na = 1 - na 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 \ No newline at end of file + 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)) + end if + na = 1 - na + end if + l1 = l2 + iw = iw + (ip - 1)*idot + end do + if (na == 0) return + n2 = n + n + do concurrent(i=1:n2) + c(i) = ch(i) + end do + end subroutine cfftf1 diff --git a/src/cffti1.f90 b/src/cffti1.f90 index 1335697..6647e9e 100644 --- a/src/cffti1.f90 +++ b/src/cffti1.f90 @@ -1,68 +1,70 @@ - 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 \ No newline at end of file + subroutine cffti1(n, wa, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + integer, intent(out) :: ifac(*) + real(dp), intent(out) :: wa(*) + real(dp) :: arg, argh, argld, fi + integer :: i, i1, ib, ido, idot, ii, ip, ipm, j, k1, & + l1, l2, ld, nf, nl, nq, nr, ntry + integer, dimension(4), parameter :: ntryh = [3, 4, 2, 5] + real(dp), parameter :: tpi = 2.0_dp*acos(-1.0_dp) ! 2 * pi + nl = n + nf = 0 + j = 0 +100 j = j + 1 + if (j <= 4) then + ntry = ntryh(j) + else + ntry = ntry + 2 + end if +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) + end do + ifac(3) = 2 + end if + end if + if (nl /= 1) goto 200 + ifac(1) = n + ifac(2) = nf + argh = tpi/real(n, kind=dp) + 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_dp + wa(i) = 0.0_dp + ld = ld + l1 + fi = 0.0_dp + argld = real(ld, kind=dp)*argh + do ii = 4, idot, 2 + i = i + 2 + fi = fi + 1.0_dp + arg = fi*argld + wa(i - 1) = cos(arg) + wa(i) = sin(arg) + end do + if (ip > 5) then + wa(i1 - 1) = wa(i - 1) + wa(i1) = wa(i) + end if + end do + l1 = l2 + end do + end subroutine cffti1 diff --git a/src/cosqb1.f90 b/src/cosqb1.f90 index 11d7f4a..3d2b699 100644 --- a/src/cosqb1.f90 +++ b/src/cosqb1.f90 @@ -1,30 +1,33 @@ - subroutine cosqb1(n,x,w,Xh) - use fftpack_kind - implicit none - integer :: i , k , kc , modn , n , np2 , ns2 - real(rk) :: w , x , Xh , xim1 - dimension x(*) , w(*) , Xh(*) - ns2 = (n+1)/2 - np2 = n + 2 - do i = 3 , n , 2 - xim1 = x(i-1) + x(i) - x(i) = x(i) - x(i-1) - x(i-1) = xim1 - enddo - x(1) = x(1) + x(1) - modn = mod(n,2) - if ( modn==0 ) x(n) = x(n) + x(n) - call dfftb(n,x,Xh) - do k = 2 , ns2 - kc = np2 - k - Xh(k) = w(k-1)*x(kc) + w(kc-1)*x(k) - Xh(kc) = w(k-1)*x(k) - w(kc-1)*x(kc) - enddo - if ( modn==0 ) x(ns2+1) = w(ns2)*(x(ns2+1)+x(ns2+1)) - do k = 2 , ns2 - kc = np2 - k - x(k) = Xh(k) + Xh(kc) - x(kc) = Xh(k) - Xh(kc) - enddo - x(1) = x(1) + x(1) - end subroutine cosqb1 \ No newline at end of file + subroutine cosqb1(n, x, w, xh) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: w(*) + real(dp), intent(out) :: xh(*) + integer :: i, k, kc, modn, np2, ns2 + real(dp) :: xim1 + ns2 = (n + 1)/2 + np2 = n + 2 + do i = 3, n, 2 + xim1 = x(i - 1) + x(i) + x(i) = x(i) - x(i - 1) + x(i - 1) = xim1 + end do + x(1) = x(1) + x(1) + modn = mod(n, 2) + if (modn == 0) x(n) = x(n) + x(n) + call dfftb(n, x, xh) + do k = 2, ns2 + kc = np2 - k + xh(k) = w(k - 1)*x(kc) + w(kc - 1)*x(k) + xh(kc) = w(k - 1)*x(k) - w(kc - 1)*x(kc) + end do + if (modn == 0) x(ns2 + 1) = w(ns2)*(x(ns2 + 1) + x(ns2 + 1)) + do k = 2, ns2 + kc = np2 - k + x(k) = xh(k) + xh(kc) + x(kc) = xh(k) - xh(kc) + end do + x(1) = x(1) + x(1) + end subroutine cosqb1 diff --git a/src/cosqf1.f90 b/src/cosqf1.f90 index 8ef39f4..21ff74e 100644 --- a/src/cosqf1.f90 +++ b/src/cosqf1.f90 @@ -1,28 +1,31 @@ - subroutine cosqf1(n,x,w,Xh) - use fftpack_kind - implicit none - integer :: i , k , kc , modn , n , np2 , ns2 - real(rk) :: w , x , Xh , xim1 - dimension x(*) , w(*) , Xh(*) - ns2 = (n+1)/2 - np2 = n + 2 - do k = 2 , ns2 - kc = np2 - k - Xh(k) = x(k) + x(kc) - Xh(kc) = x(k) - x(kc) - enddo - modn = mod(n,2) - if ( modn==0 ) Xh(ns2+1) = x(ns2+1) + x(ns2+1) - do k = 2 , ns2 - kc = np2 - k - x(k) = w(k-1)*Xh(kc) + w(kc-1)*Xh(k) - x(kc) = w(k-1)*Xh(k) - w(kc-1)*Xh(kc) - enddo - if ( modn==0 ) x(ns2+1) = w(ns2)*Xh(ns2+1) - call dfftf(n,x,Xh) - do i = 3 , n , 2 - xim1 = x(i-1) - x(i) - x(i) = x(i-1) + x(i) - x(i-1) = xim1 - enddo - end subroutine cosqf1 \ No newline at end of file + subroutine cosqf1(n, x, w, xh) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: w(*) + real(dp), intent(out) :: xh(*) + integer :: i, k, kc, modn, np2, ns2 + real(dp) :: xim1 + ns2 = (n + 1)/2 + np2 = n + 2 + do k = 2, ns2 + kc = np2 - k + xh(k) = x(k) + x(kc) + xh(kc) = x(k) - x(kc) + end do + modn = mod(n, 2) + if (modn == 0) xh(ns2 + 1) = x(ns2 + 1) + x(ns2 + 1) + do k = 2, ns2 + kc = np2 - k + x(k) = w(k - 1)*xh(kc) + w(kc - 1)*xh(k) + x(kc) = w(k - 1)*xh(k) - w(kc - 1)*xh(kc) + end do + if (modn == 0) x(ns2 + 1) = w(ns2)*xh(ns2 + 1) + call dfftf(n, x, xh) + do i = 3, n, 2 + xim1 = x(i - 1) - x(i) + x(i) = x(i - 1) + x(i) + x(i - 1) = xim1 + end do + end subroutine cosqf1 diff --git a/src/dcosqb.f90 b/src/dcosqb.f90 index 1ab8de7..9b3c34c 100644 --- a/src/dcosqb.f90 +++ b/src/dcosqb.f90 @@ -1,19 +1,20 @@ - subroutine dcosqb(n,x,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: Wsave , x , x1 - dimension x(*) , Wsave(*) - real(rk),parameter :: tsqrt2 = 2.0_rk * sqrt(2.0_rk) - if ( n<2 ) then - x(1) = 4.0_rk*x(1) - return - elseif ( n==2 ) then - x1 = 4.0_rk*(x(1)+x(2)) - x(2) = tsqrt2*(x(1)-x(2)) - x(1) = x1 - return - else - call cosqb1(n,x,Wsave,Wsave(n+1)) - endif - end subroutine dcosqb \ No newline at end of file + subroutine dcosqb(n, x, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(in) :: wsave(*) + real(dp), intent(inout) :: x(*) + real(dp) :: x1 + real(dp), parameter :: tsqrt2 = 2.0_dp*sqrt(2.0_dp) + if (n < 2) then + x(1) = 4.0_dp*x(1) + return + elseif (n == 2) then + x1 = 4.0_dp*(x(1) + x(2)) + x(2) = tsqrt2*(x(1) - x(2)) + x(1) = x1 + return + else + call cosqb1(n, x, wsave, wsave(n + 1)) + end if + end subroutine dcosqb diff --git a/src/dcosqf.f90 b/src/dcosqf.f90 index 27f4ceb..e823387 100644 --- a/src/dcosqf.f90 +++ b/src/dcosqf.f90 @@ -1,17 +1,18 @@ - subroutine dcosqf(n,x,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: tsqx , Wsave , x - dimension x(*) , Wsave(*) - real(rk),parameter :: sqrt2 = sqrt(2.0_rk) - if ( n<2 ) then - return - elseif ( n==2 ) then - tsqx = sqrt2*x(2) - x(2) = x(1) - tsqx - x(1) = x(1) + tsqx - else - call cosqf1(n,x,Wsave,Wsave(n+1)) - endif - end subroutine dcosqf \ No newline at end of file + subroutine dcosqf(n, x, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(in) :: wsave(*) + real(dp), intent(inout) :: x(*) + real(dp) :: tsqx + real(dp), parameter :: sqrt2 = sqrt(2.0_dp) + if (n < 2) then + return + elseif (n == 2) then + tsqx = sqrt2*x(2) + x(2) = x(1) - tsqx + x(1) = x(1) + tsqx + else + call cosqf1(n, x, wsave, wsave(n + 1)) + end if + end subroutine dcosqf diff --git a/src/dcosqi.f90 b/src/dcosqi.f90 index 1faf8d0..184b214 100644 --- a/src/dcosqi.f90 +++ b/src/dcosqi.f90 @@ -1,15 +1,16 @@ - subroutine dcosqi(n,Wsave) - use fftpack_kind - implicit none - real(rk) :: dt , fk , Wsave - integer :: k , n - dimension Wsave(*) - real(rk),parameter :: pih = acos(-1.0_rk) / 2.0_rk ! pi / 2 - dt = pih/real(n, rk) - fk = 0.0_rk - do k = 1 , n - fk = fk + 1.0_rk - Wsave(k) = cos(fk*dt) - enddo - call dffti(n,Wsave(n+1)) - end subroutine dcosqi \ No newline at end of file + subroutine dcosqi(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + real(dp) :: dt, fk + integer :: k + real(dp), parameter :: pih = acos(-1.0_dp)/2.0_dp ! pi / 2 + dt = pih/real(n, kind=dp) + fk = 0.0_dp + do k = 1, n + fk = fk + 1.0_dp + wsave(k) = cos(fk*dt) + end do + call dffti(n, wsave(n + 1)) + end subroutine dcosqi diff --git a/src/dcost.f90 b/src/dcost.f90 index fdc431d..4273fbb 100644 --- a/src/dcost.f90 +++ b/src/dcost.f90 @@ -1,48 +1,50 @@ - subroutine dcost(n,x,Wsave) - use fftpack_kind - implicit none - real(rk) :: c1 , t1 , t2 , tx2 , Wsave , x , x1h , x1p3 , & - xi , xim2 - integer :: i , k , kc , modn , n , nm1 , np1 , ns2 - dimension x(*) , Wsave(*) - nm1 = n - 1 - np1 = n + 1 - ns2 = n/2 - if ( n<2 ) return - if ( n==2 ) then - x1h = x(1) + x(2) - x(2) = x(1) - x(2) - x(1) = x1h - return - elseif ( n>3 ) then - c1 = x(1) - x(n) - x(1) = x(1) + x(n) - do k = 2 , ns2 - kc = np1 - k - t1 = x(k) + x(kc) - t2 = x(k) - x(kc) - c1 = c1 + Wsave(kc)*t2 - t2 = Wsave(k)*t2 - x(k) = t1 - t2 - x(kc) = t1 + t2 - enddo - modn = mod(n,2) - if ( modn/=0 ) x(ns2+1) = x(ns2+1) + x(ns2+1) - call dfftf(nm1,x,Wsave(n+1)) - xim2 = x(2) - x(2) = c1 - do i = 4 , n , 2 - xi = x(i) - x(i) = x(i-2) - x(i-1) - x(i-1) = xim2 - xim2 = xi - enddo - if ( modn/=0 ) x(n) = xim2 - return - endif - x1p3 = x(1) + x(3) - tx2 = x(2) + x(2) - x(2) = x(1) - x(3) - x(1) = x1p3 + tx2 - x(3) = x1p3 - tx2 - end subroutine dcost \ No newline at end of file + subroutine dcost(n, x, wsave) + use fftpack_kind, only: rk + implicit none + integer, intent(in) :: n + real(rk), intent(inout) :: wsave(*) + real(rk), intent(inout) :: x(*) + real(rk) :: c1, t1, t2, tx2, x1h, x1p3, & + xi, xim2 + integer :: i, k, kc, modn, nm1, np1, ns2 + nm1 = n - 1 + np1 = n + 1 + ns2 = n/2 + if (n < 2) return + if (n == 2) then + x1h = x(1) + x(2) + x(2) = x(1) - x(2) + x(1) = x1h + return + elseif (n > 3) then + c1 = x(1) - x(n) + x(1) = x(1) + x(n) + do k = 2, ns2 + kc = np1 - k + t1 = x(k) + x(kc) + t2 = x(k) - x(kc) + c1 = c1 + wsave(kc)*t2 + t2 = wsave(k)*t2 + x(k) = t1 - t2 + x(kc) = t1 + t2 + end do + modn = mod(n, 2) + if (modn /= 0) x(ns2 + 1) = x(ns2 + 1) + x(ns2 + 1) + call dfftf(nm1, x, wsave(n + 1)) + xim2 = x(2) + x(2) = c1 + do i = 4, n, 2 + xi = x(i) + x(i) = x(i - 2) - x(i - 1) + x(i - 1) = xim2 + xim2 = xi + end do + if (modn /= 0) x(n) = xim2 + return + end if + x1p3 = x(1) + x(3) + tx2 = x(2) + x(2) + x(2) = x(1) - x(3) + x(1) = x1p3 + tx2 + x(3) = x1p3 - tx2 + end subroutine dcost diff --git a/src/dcosti.f90 b/src/dcosti.f90 index ed4fc1d..ad68e85 100644 --- a/src/dcosti.f90 +++ b/src/dcosti.f90 @@ -1,21 +1,22 @@ - subroutine dcosti(n,Wsave) - use fftpack_kind - implicit none - real(rk) :: dt , fk , Wsave - integer :: k , kc , n , nm1 , np1 , ns2 - dimension Wsave(*) - real(rk),parameter :: pi = acos(-1.0_rk) - if ( n<=3 ) return - nm1 = n - 1 - np1 = n + 1 - ns2 = n/2 - dt = pi/real(nm1, rk) - fk = 0.0_rk - do k = 2 , ns2 - kc = np1 - k - fk = fk + 1.0_rk - Wsave(k) = 2.0_rk*sin(fk*dt) - Wsave(kc) = 2.0_rk*cos(fk*dt) - enddo - call dffti(nm1,Wsave(n+1)) - end subroutine dcosti \ No newline at end of file + subroutine dcosti(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + real(dp) :: dt, fk + integer :: k, kc, nm1, np1, ns2 + real(dp), parameter :: pi = acos(-1.0_dp) + if (n <= 3) return + nm1 = n - 1 + np1 = n + 1 + ns2 = n/2 + dt = pi/real(nm1, kind=dp) + fk = 0.0_dp + do k = 2, ns2 + kc = np1 - k + fk = fk + 1.0_dp + wsave(k) = 2.0_dp*sin(fk*dt) + wsave(kc) = 2.0_dp*cos(fk*dt) + end do + call dffti(nm1, wsave(n + 1)) + end subroutine dcosti diff --git a/src/dfftb.f90 b/src/dfftb.f90 index 045edc1..8620e63 100644 --- a/src/dfftb.f90 +++ b/src/dfftb.f90 @@ -1,9 +1,9 @@ - subroutine dfftb(n,r,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: r , Wsave - dimension r(*) , Wsave(*) - if ( n==1 ) return - call rfftb1(n,r,Wsave,Wsave(n+1),Wsave(2*n+1)) - end subroutine dfftb \ No newline at end of file + subroutine dfftb(n, r, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: r(*) + real(dp), intent(inout) :: wsave(*) + if (n == 1) return + call rfftb1(n, r, wsave, wsave(n + 1), wsave(2*n + 1)) + end subroutine dfftb diff --git a/src/dfftf.f90 b/src/dfftf.f90 index d23437e..39064ef 100644 --- a/src/dfftf.f90 +++ b/src/dfftf.f90 @@ -1,9 +1,9 @@ - subroutine dfftf(n,r,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: r , Wsave - dimension r(*) , Wsave(*) - if ( n==1 ) return - call rfftf1(n,r,Wsave,Wsave(n+1),Wsave(2*n+1)) - end subroutine dfftf \ No newline at end of file + subroutine dfftf(n, r, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: r(*) + real(dp), intent(inout) :: wsave(*) + if (n == 1) return + call rfftf1(n, r, wsave, wsave(n + 1), wsave(2*n + 1)) + end subroutine dfftf diff --git a/src/dffti.f90 b/src/dffti.f90 index d9f07db..ec7a341 100644 --- a/src/dffti.f90 +++ b/src/dffti.f90 @@ -1,9 +1,8 @@ - subroutine dffti(n,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: Wsave - dimension Wsave(*) - if ( n==1 ) return - call rffti1(n,Wsave(n+1),Wsave(2*n+1)) - end subroutine dffti \ No newline at end of file + subroutine dffti(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + if (n == 1) return + call rffti1(n, wsave(n + 1), wsave(2*n + 1)) + end subroutine dffti diff --git a/src/dsinqb.f90 b/src/dsinqb.f90 index 1832e6b..dfc8504 100644 --- a/src/dsinqb.f90 +++ b/src/dsinqb.f90 @@ -1,23 +1,25 @@ - subroutine dsinqb(n,x,Wsave) - use fftpack_kind - implicit none - integer :: k , kc , n , ns2 - real(rk) :: Wsave , x , xhold - dimension x(*) , Wsave(*) - if ( n>1 ) then - ns2 = n/2 - do k = 2 , n , 2 - x(k) = -x(k) - enddo - call dcosqb(n,x,Wsave) - do k = 1 , ns2 - kc = n - k - xhold = x(k) - x(k) = x(kc+1) - x(kc+1) = xhold - enddo + subroutine dsinqb(n, x, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + integer :: k, kc, ns2 + real(dp) :: xhold + if (n > 1) then + ns2 = n/2 + do k = 2, n, 2 + x(k) = -x(k) + end do + call dcosqb(n, x, wsave) + do k = 1, ns2 + kc = n - k + xhold = x(k) + x(k) = x(kc + 1) + x(kc + 1) = xhold + end do + return + end if + x(1) = 4.0_dp*x(1) return - endif - x(1) = 4.0_rk*x(1) - return - end subroutine dsinqb \ No newline at end of file + end subroutine dsinqb diff --git a/src/dsinqf.f90 b/src/dsinqf.f90 index 66e7312..388bb17 100644 --- a/src/dsinqf.f90 +++ b/src/dsinqf.f90 @@ -1,19 +1,21 @@ - subroutine dsinqf(n,x,Wsave) - use fftpack_kind - implicit none - integer :: k , kc , n , ns2 - real(rk) :: Wsave , x , xhold - dimension x(*) , Wsave(*) - if ( n==1 ) return - ns2 = n/2 - do k = 1 , ns2 - kc = n - k - xhold = x(k) - x(k) = x(kc+1) - x(kc+1) = xhold - enddo - call dcosqf(n,x,Wsave) - do k = 2 , n , 2 - x(k) = -x(k) - enddo - end subroutine dsinqf \ No newline at end of file + subroutine dsinqf(n, x, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + integer :: k, kc, ns2 + real(dp) :: xhold + if (n == 1) return + ns2 = n/2 + do k = 1, ns2 + kc = n - k + xhold = x(k) + x(k) = x(kc + 1) + x(kc + 1) = xhold + end do + call dcosqf(n, x, wsave) + do k = 2, n, 2 + x(k) = -x(k) + end do + end subroutine dsinqf diff --git a/src/dsinqi.f90 b/src/dsinqi.f90 index d95e2db..642b795 100644 --- a/src/dsinqi.f90 +++ b/src/dsinqi.f90 @@ -1,8 +1,7 @@ - subroutine dsinqi(n,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: Wsave - dimension Wsave(*) - call dcosqi(n,Wsave) - end subroutine dsinqi \ No newline at end of file + subroutine dsinqi(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + call dcosqi(n, wsave) + end subroutine dsinqi diff --git a/src/dsint.f90 b/src/dsint.f90 index 8c97f2f..3515772 100644 --- a/src/dsint.f90 +++ b/src/dsint.f90 @@ -1,12 +1,13 @@ - subroutine dsint(n,x,Wsave) - use fftpack_kind - implicit none - integer :: iw1 , iw2 , iw3 , n , np1 - real(rk) :: Wsave , x - dimension x(*) , Wsave(*) - np1 = n + 1 - iw1 = n/2 + 1 - iw2 = iw1 + np1 - iw3 = iw2 + np1 - call sint1(n,x,Wsave,Wsave(iw1),Wsave(iw2),Wsave(iw3)) - end subroutine dsint \ No newline at end of file + subroutine dsint(n, x, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + integer :: iw1, iw2, iw3, np1 + np1 = n + 1 + iw1 = n/2 + 1 + iw2 = iw1 + np1 + iw3 = iw2 + np1 + call sint1(n, x, wsave, wsave(iw1), wsave(iw2), wsave(iw3)) + end subroutine dsint diff --git a/src/dsinti.f90 b/src/dsinti.f90 index 2069076..a2b2cfb 100644 --- a/src/dsinti.f90 +++ b/src/dsinti.f90 @@ -1,16 +1,17 @@ - subroutine dsinti(n,Wsave) - use fftpack_kind - implicit none - real(rk) :: dt , Wsave - integer :: k , n , np1 , ns2 - dimension Wsave(*) - real(rk),parameter :: pi = acos(-1.0_rk) - if ( n<=1 ) return - ns2 = n/2 - np1 = n + 1 - dt = pi/real(np1, rk) - do k = 1 , ns2 - Wsave(k) = 2.0_rk*sin(k*dt) - enddo - call dffti(np1,Wsave(ns2+1)) - end subroutine dsinti \ No newline at end of file + subroutine dsinti(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + real(dp) :: dt + integer :: k, np1, ns2 + real(dp), parameter :: pi = acos(-1.0_dp) + if (n <= 1) return + ns2 = n/2 + np1 = n + 1 + dt = pi/real(np1, dp) + do k = 1, ns2 + wsave(k) = 2.0_dp*sin(k*dt) + end do + call dffti(np1, wsave(ns2 + 1)) + end subroutine dsinti diff --git a/src/dzfftb.f90 b/src/dzfftb.f90 index e02dcc1..1acb8c3 100644 --- a/src/dzfftb.f90 +++ b/src/dzfftb.f90 @@ -1,24 +1,26 @@ - subroutine dzfftb(n,r,Azero,a,b,Wsave) - use fftpack_kind - implicit none - real(rk) :: a , Azero , b , r , Wsave - integer :: i , n , ns2 - dimension r(*) , a(*) , b(*) , Wsave(*) - if ( n<2 ) then - r(1) = Azero - return - elseif ( n==2 ) then - r(1) = Azero + a(1) - r(2) = Azero - a(1) - return - else - ns2 = (n-1)/2 - do i = 1 , ns2 - r(2*i) = 0.5_rk*a(i) - r(2*i+1) = -0.5_rk*b(i) - enddo - r(1) = Azero - if ( mod(n,2)==0 ) r(n) = a(ns2+1) - call dfftb(n,r,Wsave(n+1)) - endif - end subroutine dzfftb \ No newline at end of file + subroutine dzfftb(n, r, azero, a, b, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: r(*) + real(dp), intent(inout) :: wsave(*) + real(dp), intent(in) :: azero, a(*), b(*) + integer :: i, ns2 + if (n < 2) then + r(1) = azero + return + elseif (n == 2) then + r(1) = azero + a(1) + r(2) = azero - a(1) + return + else + ns2 = (n - 1)/2 + do concurrent(i=1:ns2) + r(2*i) = 0.5_dp*a(i) + r(2*i + 1) = -0.5_dp*b(i) + end do + r(1) = azero + if (mod(n, 2) == 0) r(n) = a(ns2 + 1) + call dfftb(n, r, wsave(n + 1)) + end if + end subroutine dzfftb diff --git a/src/dzfftf.f90 b/src/dzfftf.f90 index d122b9b..3ce5d8a 100644 --- a/src/dzfftf.f90 +++ b/src/dzfftf.f90 @@ -1,35 +1,36 @@ - subroutine dzfftf(n,r,Azero,a,b,Wsave) -! -! VERSION 3 JUNE 1979 -! - use fftpack_kind - implicit none - real(rk) :: a , Azero , b , cf , cfm , r , Wsave - integer :: i , n , ns2 , ns2m - dimension r(*) , a(*) , b(*) , Wsave(*) - if ( n<2 ) then - Azero = r(1) - return - elseif ( n==2 ) then - Azero = 0.5_rk*(r(1)+r(2)) - a(1) = 0.5_rk*(r(1)-r(2)) - return - else - do i = 1 , n - Wsave(i) = r(i) - enddo - call dfftf(n,Wsave,Wsave(n+1)) - cf = 2.0_rk/real(n, rk) - cfm = -cf - Azero = 0.5_rk*cf*Wsave(1) - ns2 = (n+1)/2 - ns2m = ns2 - 1 - do i = 1 , ns2m - a(i) = cf*Wsave(2*i) - b(i) = cfm*Wsave(2*i+1) - enddo - if ( mod(n,2)==1 ) return - a(ns2) = 0.5_rk*cf*Wsave(n) - b(ns2) = 0.0_rk - endif - end subroutine dzfftf \ No newline at end of file + subroutine dzfftf(n, r, azero, a, b, wsave) +! version 3 june 1979 + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(in) :: r(*) + real(dp), intent(out) :: azero, a(*), b(*) + real(dp), intent(inout) :: wsave(*) + real(dp) :: cf, cfm + integer :: i, ns2, ns2m + if (n < 2) then + azero = r(1) + return + elseif (n == 2) then + azero = 0.5_dp*(r(1) + r(2)) + a(1) = 0.5_dp*(r(1) - r(2)) + return + else + do concurrent(i=1:n) + wsave(i) = r(i) + end do + call dfftf(n, wsave, wsave(n + 1)) + cf = 2.0_dp/real(n, dp) + cfm = -cf + azero = 0.5_dp*cf*wsave(1) + ns2 = (n + 1)/2 + ns2m = ns2 - 1 + do concurrent(i=1:ns2m) + a(i) = cf*wsave(2*i) + b(i) = cfm*wsave(2*i + 1) + end do + if (mod(n, 2) == 1) return + a(ns2) = 0.5_dp*cf*wsave(n) + b(ns2) = 0.0_dp + end if + end subroutine dzfftf diff --git a/src/dzffti.f90 b/src/dzffti.f90 index 8b51821..4687639 100644 --- a/src/dzffti.f90 +++ b/src/dzffti.f90 @@ -1,9 +1,8 @@ - subroutine dzffti(n,Wsave) - use fftpack_kind - implicit none - integer :: n - real(rk) :: Wsave - dimension Wsave(*) - if ( n==1 ) return - call ezfft1(n,Wsave(2*n+1),Wsave(3*n+1)) - end subroutine dzffti \ No newline at end of file + subroutine dzffti(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + if (n == 1) return + call ezfft1(n, wsave(2*n + 1), wsave(3*n + 1)) + end subroutine dzffti diff --git a/src/ezfft1.f90 b/src/ezfft1.f90 index 45e1c7d..38a8915 100644 --- a/src/ezfft1.f90 +++ b/src/ezfft1.f90 @@ -1,71 +1,72 @@ - subroutine ezfft1(n,Wa,Ifac) - use fftpack_kind - implicit none - real(rk) :: arg1 , argh , ch1 , ch1h , dch1 , dsh1 , sh1 , & - Wa - integer :: i , ib , ido , Ifac , ii , ip , ipm , is , j , k1 , l1 , & - l2 , n , nf , nfm1 , nl , nq , nr , ntry - dimension Wa(*) , Ifac(*) - integer,dimension(4),parameter :: ntryh = [4 , 2 , 3 , 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) - is = 0 - nfm1 = nf - 1 - l1 = 1 - if ( nfm1==0 ) return - do k1 = 1 , nfm1 - ip = Ifac(k1+2) - l2 = l1*ip - ido = n/l2 - ipm = ip - 1 - arg1 = real(l1, rk)*argh - ch1 = 1.0_rk - sh1 = 0.0_rk - dch1 = cos(arg1) - dsh1 = sin(arg1) - do j = 1 , ipm - ch1h = dch1*ch1 - dsh1*sh1 - sh1 = dch1*sh1 + dsh1*ch1 - ch1 = ch1h - i = is + 2 - Wa(i-1) = ch1 - Wa(i) = sh1 - if ( ido>=5 ) then - do ii = 5 , ido , 2 - i = i + 2 - Wa(i-1) = ch1*Wa(i-3) - sh1*Wa(i-2) - Wa(i) = ch1*Wa(i-2) + sh1*Wa(i-3) - enddo - endif - is = is + ido - enddo - l1 = l2 - enddo - end subroutine ezfft1 \ No newline at end of file + subroutine ezfft1(n, wa, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wa(*) + integer, intent(out) :: ifac(*) + real(dp) :: arg1, argh, ch1, ch1h, dch1, dsh1, sh1 + integer :: i, ib, ido, ii, ip, ipm, is, j, k1, l1, & + l2, nf, nfm1, nl, nq, nr, ntry + integer, dimension(4), parameter :: ntryh = [4, 2, 3, 5] + real(dp), parameter :: tpi = 2.0_dp*acos(-1.0_dp) ! 2 * pi + nl = n + nf = 0 + j = 0 +100 j = j + 1 + if (j <= 4) then + ntry = ntryh(j) + else + ntry = ntry + 2 + end if +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) + end do + ifac(3) = 2 + end if + end if + if (nl /= 1) goto 200 + ifac(1) = n + ifac(2) = nf + argh = tpi/real(n, dp) + is = 0 + nfm1 = nf - 1 + l1 = 1 + if (nfm1 == 0) return + do k1 = 1, nfm1 + ip = ifac(k1 + 2) + l2 = l1*ip + ido = n/l2 + ipm = ip - 1 + arg1 = real(l1, dp)*argh + ch1 = 1.0_dp + sh1 = 0.0_dp + dch1 = cos(arg1) + dsh1 = sin(arg1) + do j = 1, ipm + ch1h = dch1*ch1 - dsh1*sh1 + sh1 = dch1*sh1 + dsh1*ch1 + ch1 = ch1h + i = is + 2 + wa(i - 1) = ch1 + wa(i) = sh1 + if (ido >= 5) then + do ii = 5, ido, 2 + i = i + 2 + wa(i - 1) = ch1*wa(i - 3) - sh1*wa(i - 2) + wa(i) = ch1*wa(i - 2) + sh1*wa(i - 3) + end do + end if + is = is + ido + end do + l1 = l2 + end do + end subroutine ezfft1 diff --git a/src/passb.f90 b/src/passb.f90 index 12c3928..0b43577 100644 --- a/src/passb.f90 +++ b/src/passb.f90 @@ -1,125 +1,110 @@ - subroutine passb(Nac,Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) - use fftpack_kind - implicit none - real(rk) :: c1 , c2 , Cc , Ch , Ch2 , Wa , wai , war - integer :: i , idij , idj , idl , Idl1 , idlj , Ido , idot , idp , & - ik , inc , Ip , ipp2 , ipph , j , jc , k , l , l1 , lc - integer :: Nac , nt - dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , Wa(*) , & - c2(Idl1,Ip) , Ch2(Idl1,Ip) - idot = Ido/2 - nt = Ip*Idl1 - ipp2 = Ip + 2 - ipph = (Ip+1)/2 - idp = Ip*Ido + subroutine passb(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(out) :: nac + integer, intent(in) :: ido, ip, l1, idl1 + real(dp), intent(in) :: cc(ido, ip, l1), wa(*) + real(dp), intent(out) :: c1(ido, l1, ip), c2(idl1, ip), ch(ido, l1, ip) + real(dp), intent(inout) :: ch2(idl1, ip) + real(dp) :: wai, war + integer :: i, idij, idj, idl, idlj, idot, idp, & + ik, inc, ipp2, ipph, j, jc, k, l, lc + integer :: nt + idot = ido/2 + nt = ip*idl1 + ipp2 = ip + 2 + ipph = (ip + 1)/2 + idp = ip*ido ! - if ( Ido idp) idlj = idlj - idp + war = wa(idlj - 1) + wai = wa(idlj) + do concurrent(ik=1:idl1) + c2(ik, l) = c2(ik, l) + war*ch2(ik, j) + c2(ik, lc) = c2(ik, lc) + wai*ch2(ik, jc) + end do + end do + end do + do concurrent(ik=1:idl1, j=2:ipph) + ch2(ik, 1) = ch2(ik, 1) + ch2(ik, j) + end do + do concurrent(ik=2:idl1:2, j=2:ipph) jc = ipp2 - j - do i = 1 , Ido - do k = 1 , l1 - Ch(i,k,j) = Cc(i,j,k) + Cc(i,jc,k) - Ch(i,k,jc) = Cc(i,j,k) - Cc(i,jc,k) - enddo - enddo - enddo - do i = 1 , Ido - do k = 1 , l1 - Ch(i,k,1) = Cc(i,1,k) - enddo - enddo - else - do j = 2 , ipph - jc = ipp2 - j - do k = 1 , l1 - do i = 1 , Ido - Ch(i,k,j) = Cc(i,j,k) + Cc(i,jc,k) - Ch(i,k,jc) = Cc(i,j,k) - Cc(i,jc,k) - enddo - enddo - enddo - do k = 1 , l1 - do i = 1 , Ido - Ch(i,k,1) = Cc(i,1,k) - enddo - enddo - endif - idl = 2 - Ido - inc = 0 - do l = 2 , ipph - lc = ipp2 - l - idl = idl + Ido - do ik = 1 , Idl1 - c2(ik,l) = Ch2(ik,1) + Wa(idl-1)*Ch2(ik,2) - c2(ik,lc) = Wa(idl)*Ch2(ik,Ip) - enddo - idlj = idl - inc = inc + Ido - do j = 3 , ipph - jc = ipp2 - j - idlj = idlj + inc - if ( idlj>idp ) idlj = idlj - idp - war = Wa(idlj-1) - wai = Wa(idlj) - do ik = 1 , Idl1 - c2(ik,l) = c2(ik,l) + war*Ch2(ik,j) - c2(ik,lc) = c2(ik,lc) + wai*Ch2(ik,jc) - enddo - enddo - enddo - do j = 2 , ipph - do ik = 1 , Idl1 - Ch2(ik,1) = Ch2(ik,1) + Ch2(ik,j) - enddo - enddo - do j = 2 , ipph - jc = ipp2 - j - do ik = 2 , Idl1 , 2 - Ch2(ik-1,j) = c2(ik-1,j) - c2(ik,jc) - Ch2(ik-1,jc) = c2(ik-1,j) + c2(ik,jc) - Ch2(ik,j) = c2(ik,j) + c2(ik-1,jc) - Ch2(ik,jc) = c2(ik,j) - c2(ik-1,jc) - enddo - enddo - Nac = 1 - if ( Ido==2 ) return - Nac = 0 - do ik = 1 , Idl1 - c2(ik,1) = Ch2(ik,1) - enddo - do j = 2 , Ip - do k = 1 , l1 - c1(1,k,j) = Ch(1,k,j) - c1(2,k,j) = Ch(2,k,j) - enddo - enddo - if ( idot>l1 ) then - idj = 2 - Ido - do j = 2 , Ip - idj = idj + Ido - do k = 1 , l1 - idij = idj - do i = 4 , Ido , 2 - idij = idij + 2 - c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij) & - *Ch(i,k,j) - c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij) & - *Ch(i-1,k,j) - enddo - enddo - enddo - return - endif - idij = 0 - do j = 2 , Ip - idij = idij + 2 - do i = 4 , Ido , 2 + ch2(ik - 1, j) = c2(ik - 1, j) - c2(ik, jc) + ch2(ik - 1, jc) = c2(ik - 1, j) + c2(ik, jc) + ch2(ik, j) = c2(ik, j) + c2(ik - 1, jc) + ch2(ik, jc) = c2(ik, j) - c2(ik - 1, jc) + end do + nac = 1 + if (ido == 2) return + nac = 0 + do concurrent(ik=1:idl1) + c2(ik, 1) = ch2(ik, 1) + end do + do concurrent(j=2:ip, k=1:l1) + c1(1, k, j) = ch(1, k, j) + c1(2, k, j) = ch(2, k, j) + end do + if (idot > l1) then + idj = 2 - ido + do j = 2, ip + idj = idj + ido + do k = 1, l1 + idij = idj + do i = 4, ido, 2 + idij = idij + 2 + c1(i - 1, k, j) = wa(idij - 1)*ch(i - 1, k, j) - wa(idij) & + *ch(i, k, j) + c1(i, k, j) = wa(idij - 1)*ch(i, k, j) + wa(idij) & + *ch(i - 1, k, j) + end do + end do + end do + return + end if + idij = 0 + do j = 2, ip idij = idij + 2 - do k = 1 , l1 - c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij)*Ch(i,k,j) - c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij)*Ch(i-1,k,j) - enddo - enddo - enddo - return - end subroutine passb \ No newline at end of file + do i = 4, ido, 2 + idij = idij + 2 + do concurrent(k=1:l1) + c1(i - 1, k, j) = wa(idij - 1)*ch(i - 1, k, j) - wa(idij)*ch(i, k, j) + c1(i, k, j) = wa(idij - 1)*ch(i, k, j) + wa(idij)*ch(i - 1, k, j) + end do + end do + end do + return + end subroutine passb diff --git a/src/passb2.f90 b/src/passb2.f90 index f75bf41..3d68d5a 100644 --- a/src/passb2.f90 +++ b/src/passb2.f90 @@ -1,26 +1,26 @@ - subroutine passb2(Ido,l1,Cc,Ch,Wa1) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 - integer :: i , Ido , k , l1 - dimension Cc(Ido,2,l1) , Ch(Ido,l1,2) , Wa1(*) - if ( Ido>2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - Ch(i-1,k,1) = Cc(i-1,1,k) + Cc(i-1,2,k) - tr2 = Cc(i-1,1,k) - Cc(i-1,2,k) - Ch(i,k,1) = Cc(i,1,k) + Cc(i,2,k) - ti2 = Cc(i,1,k) - Cc(i,2,k) - Ch(i,k,2) = Wa1(i-1)*ti2 + Wa1(i)*tr2 - Ch(i-1,k,2) = Wa1(i-1)*tr2 - Wa1(i)*ti2 - enddo - enddo - else - do k = 1 , l1 - Ch(1,k,1) = Cc(1,1,k) + Cc(1,2,k) - Ch(1,k,2) = Cc(1,1,k) - Cc(1,2,k) - Ch(2,k,1) = Cc(2,1,k) + Cc(2,2,k) - Ch(2,k,2) = Cc(2,1,k) - Cc(2,2,k) - enddo - end if - end subroutine passb2 \ No newline at end of file + subroutine passb2(ido, l1, cc, ch, wa1) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 2, l1), wa1(*) + real(dp), intent(out) :: ch(ido, l1, 2) + real(dp) :: ti2, tr2 + integer :: i, k + if (ido > 2) then + do concurrent(i=2:ido:2, k=1:l1) + ch(i - 1, k, 1) = cc(i - 1, 1, k) + cc(i - 1, 2, k) + tr2 = cc(i - 1, 1, k) - cc(i - 1, 2, k) + ch(i, k, 1) = cc(i, 1, k) + cc(i, 2, k) + ti2 = cc(i, 1, k) - cc(i, 2, k) + ch(i, k, 2) = wa1(i - 1)*ti2 + wa1(i)*tr2 + ch(i - 1, k, 2) = wa1(i - 1)*tr2 - wa1(i)*ti2 + end do + else + do concurrent(k=1:l1) + ch(1, k, 1) = cc(1, 1, k) + cc(1, 2, k) + ch(1, k, 2) = cc(1, 1, k) - cc(1, 2, k) + ch(2, k, 1) = cc(2, 1, k) + cc(2, 2, k) + ch(2, k, 2) = cc(2, 1, k) - cc(2, 2, k) + end do + end if + end subroutine passb2 diff --git a/src/passb3.f90 b/src/passb3.f90 index 0161f73..f0b3c75 100644 --- a/src/passb3.f90 +++ b/src/passb3.f90 @@ -1,47 +1,47 @@ - subroutine passb3(Ido,l1,Cc,Ch,Wa1,Wa2) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , cr2 , cr3 , di2 , di3 , & - dr2 , dr3 , ti2 , tr2 , Wa1 , Wa2 - integer :: i , Ido , k , l1 - dimension Cc(Ido,3,l1) , Ch(Ido,l1,3) , Wa1(*) , Wa2(*) - real(rk),parameter :: taur = -0.5_rk - real(rk),parameter :: taui = sqrt(3.0_rk) / 2.0_rk - if ( Ido/=2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - tr2 = Cc(i-1,2,k) + Cc(i-1,3,k) - cr2 = Cc(i-1,1,k) + taur*tr2 - Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 - ti2 = Cc(i,2,k) + Cc(i,3,k) - ci2 = Cc(i,1,k) + taur*ti2 - Ch(i,k,1) = Cc(i,1,k) + ti2 - cr3 = taui*(Cc(i-1,2,k)-Cc(i-1,3,k)) - ci3 = taui*(Cc(i,2,k)-Cc(i,3,k)) + subroutine passb3(ido, l1, cc, ch, wa1, wa2) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 3, l1), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, l1, 3) + real(dp) :: ci2, ci3, cr2, cr3, di2, di3, & + dr2, dr3, ti2, tr2 + integer :: i, k + real(dp), parameter :: taur = -0.5_dp + real(dp), parameter :: taui = sqrt(3.0_dp)/2.0_dp + if (ido /= 2) then + do concurrent(i=2:ido:2, k=1:l1) + tr2 = cc(i - 1, 2, k) + cc(i - 1, 3, k) + cr2 = cc(i - 1, 1, k) + taur*tr2 + ch(i - 1, k, 1) = cc(i - 1, 1, k) + tr2 + ti2 = cc(i, 2, k) + cc(i, 3, k) + ci2 = cc(i, 1, k) + taur*ti2 + ch(i, k, 1) = cc(i, 1, k) + ti2 + cr3 = taui*(cc(i - 1, 2, k) - cc(i - 1, 3, k)) + ci3 = taui*(cc(i, 2, k) - cc(i, 3, k)) dr2 = cr2 - ci3 dr3 = cr2 + ci3 di2 = ci2 + cr3 di3 = ci2 - cr3 - Ch(i,k,2) = Wa1(i-1)*di2 + Wa1(i)*dr2 - Ch(i-1,k,2) = Wa1(i-1)*dr2 - Wa1(i)*di2 - Ch(i,k,3) = Wa2(i-1)*di3 + Wa2(i)*dr3 - Ch(i-1,k,3) = Wa2(i-1)*dr3 - Wa2(i)*di3 - enddo - enddo - else - do k = 1 , l1 - tr2 = Cc(1,2,k) + Cc(1,3,k) - cr2 = Cc(1,1,k) + taur*tr2 - Ch(1,k,1) = Cc(1,1,k) + tr2 - ti2 = Cc(2,2,k) + Cc(2,3,k) - ci2 = Cc(2,1,k) + taur*ti2 - Ch(2,k,1) = Cc(2,1,k) + ti2 - cr3 = taui*(Cc(1,2,k)-Cc(1,3,k)) - ci3 = taui*(Cc(2,2,k)-Cc(2,3,k)) - Ch(1,k,2) = cr2 - ci3 - Ch(1,k,3) = cr2 + ci3 - Ch(2,k,2) = ci2 + cr3 - Ch(2,k,3) = ci2 - cr3 - enddo - end if - end subroutine passb3 \ No newline at end of file + ch(i, k, 2) = wa1(i - 1)*di2 + wa1(i)*dr2 + ch(i - 1, k, 2) = wa1(i - 1)*dr2 - wa1(i)*di2 + ch(i, k, 3) = wa2(i - 1)*di3 + wa2(i)*dr3 + ch(i - 1, k, 3) = wa2(i - 1)*dr3 - wa2(i)*di3 + end do + else + do concurrent(k=1:l1) + tr2 = cc(1, 2, k) + cc(1, 3, k) + cr2 = cc(1, 1, k) + taur*tr2 + ch(1, k, 1) = cc(1, 1, k) + tr2 + ti2 = cc(2, 2, k) + cc(2, 3, k) + ci2 = cc(2, 1, k) + taur*ti2 + ch(2, k, 1) = cc(2, 1, k) + ti2 + cr3 = taui*(cc(1, 2, k) - cc(1, 3, k)) + ci3 = taui*(cc(2, 2, k) - cc(2, 3, k)) + ch(1, k, 2) = cr2 - ci3 + ch(1, k, 3) = cr2 + ci3 + ch(2, k, 2) = ci2 + cr3 + ch(2, k, 3) = ci2 - cr3 + end do + end if + end subroutine passb3 diff --git a/src/passb4.f90 b/src/passb4.f90 index 0c78a1d..9db877b 100644 --- a/src/passb4.f90 +++ b/src/passb4.f90 @@ -1,56 +1,55 @@ - subroutine passb4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & - & ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3 , tr4 , & - & Wa1 , Wa2 , Wa3 - integer :: i , Ido , k , l1 - dimension Cc(Ido,4,l1) , Ch(Ido,l1,4) , Wa1(*) , Wa2(*) , Wa3(*) - if ( Ido/=2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - ti1 = Cc(i,1,k) - Cc(i,3,k) - ti2 = Cc(i,1,k) + Cc(i,3,k) - ti3 = Cc(i,2,k) + Cc(i,4,k) - tr4 = Cc(i,4,k) - Cc(i,2,k) - tr1 = Cc(i-1,1,k) - Cc(i-1,3,k) - tr2 = Cc(i-1,1,k) + Cc(i-1,3,k) - ti4 = Cc(i-1,2,k) - Cc(i-1,4,k) - tr3 = Cc(i-1,2,k) + Cc(i-1,4,k) - Ch(i-1,k,1) = tr2 + tr3 + subroutine passb4(ido, l1, cc, ch, wa1, wa2, wa3) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 4, l1), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, l1, 4) + real(dp) :: ci2, ci3, ci4, cr2, cr3, cr4, & + & ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4 + integer :: i, k + if (ido /= 2) then + do concurrent(i=2:ido:2, k=1:l1) + ti1 = cc(i, 1, k) - cc(i, 3, k) + ti2 = cc(i, 1, k) + cc(i, 3, k) + ti3 = cc(i, 2, k) + cc(i, 4, k) + tr4 = cc(i, 4, k) - cc(i, 2, k) + tr1 = cc(i - 1, 1, k) - cc(i - 1, 3, k) + tr2 = cc(i - 1, 1, k) + cc(i - 1, 3, k) + ti4 = cc(i - 1, 2, k) - cc(i - 1, 4, k) + tr3 = cc(i - 1, 2, k) + cc(i - 1, 4, k) + ch(i - 1, k, 1) = tr2 + tr3 cr3 = tr2 - tr3 - Ch(i,k,1) = ti2 + ti3 + ch(i, k, 1) = ti2 + ti3 ci3 = ti2 - ti3 cr2 = tr1 + tr4 cr4 = tr1 - tr4 ci2 = ti1 + ti4 ci4 = ti1 - ti4 - Ch(i-1,k,2) = Wa1(i-1)*cr2 - Wa1(i)*ci2 - Ch(i,k,2) = Wa1(i-1)*ci2 + Wa1(i)*cr2 - Ch(i-1,k,3) = Wa2(i-1)*cr3 - Wa2(i)*ci3 - Ch(i,k,3) = Wa2(i-1)*ci3 + Wa2(i)*cr3 - Ch(i-1,k,4) = Wa3(i-1)*cr4 - Wa3(i)*ci4 - Ch(i,k,4) = Wa3(i-1)*ci4 + Wa3(i)*cr4 - enddo - enddo - else - do k = 1 , l1 - ti1 = Cc(2,1,k) - Cc(2,3,k) - ti2 = Cc(2,1,k) + Cc(2,3,k) - tr4 = Cc(2,4,k) - Cc(2,2,k) - ti3 = Cc(2,2,k) + Cc(2,4,k) - tr1 = Cc(1,1,k) - Cc(1,3,k) - tr2 = Cc(1,1,k) + Cc(1,3,k) - ti4 = Cc(1,2,k) - Cc(1,4,k) - tr3 = Cc(1,2,k) + Cc(1,4,k) - Ch(1,k,1) = tr2 + tr3 - Ch(1,k,3) = tr2 - tr3 - Ch(2,k,1) = ti2 + ti3 - Ch(2,k,3) = ti2 - ti3 - Ch(1,k,2) = tr1 + tr4 - Ch(1,k,4) = tr1 - tr4 - Ch(2,k,2) = ti1 + ti4 - Ch(2,k,4) = ti1 - ti4 - enddo - end if - end subroutine passb4 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 1)*cr2 - wa1(i)*ci2 + ch(i, k, 2) = wa1(i - 1)*ci2 + wa1(i)*cr2 + ch(i - 1, k, 3) = wa2(i - 1)*cr3 - wa2(i)*ci3 + ch(i, k, 3) = wa2(i - 1)*ci3 + wa2(i)*cr3 + ch(i - 1, k, 4) = wa3(i - 1)*cr4 - wa3(i)*ci4 + ch(i, k, 4) = wa3(i - 1)*ci4 + wa3(i)*cr4 + end do + else + do concurrent(k=1:l1) + ti1 = cc(2, 1, k) - cc(2, 3, k) + ti2 = cc(2, 1, k) + cc(2, 3, k) + tr4 = cc(2, 4, k) - cc(2, 2, k) + ti3 = cc(2, 2, k) + cc(2, 4, k) + tr1 = cc(1, 1, k) - cc(1, 3, k) + tr2 = cc(1, 1, k) + cc(1, 3, k) + ti4 = cc(1, 2, k) - cc(1, 4, k) + tr3 = cc(1, 2, k) + cc(1, 4, k) + ch(1, k, 1) = tr2 + tr3 + ch(1, k, 3) = tr2 - tr3 + ch(2, k, 1) = ti2 + ti3 + ch(2, k, 3) = ti2 - ti3 + ch(1, k, 2) = tr1 + tr4 + ch(1, k, 4) = tr1 - tr4 + ch(2, k, 2) = ti1 + ti4 + ch(2, k, 4) = ti1 - ti4 + end do + end if + end subroutine passb4 diff --git a/src/passb5.f90 b/src/passb5.f90 index acdfbfe..cdb4185 100644 --- a/src/passb5.f90 +++ b/src/passb5.f90 @@ -1,36 +1,36 @@ - subroutine passb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & - cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & - dr4 , dr5 - real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & - tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 - integer :: i , Ido , k , l1 - dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), & - Wa4(*) - real(rk),parameter :: pi = acos(-1.0_rk) - real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk) - if ( Ido/=2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - ti5 = Cc(i,2,k) - Cc(i,5,k) - ti2 = Cc(i,2,k) + Cc(i,5,k) - ti4 = Cc(i,3,k) - Cc(i,4,k) - ti3 = Cc(i,3,k) + Cc(i,4,k) - tr5 = Cc(i-1,2,k) - Cc(i-1,5,k) - tr2 = Cc(i-1,2,k) + Cc(i-1,5,k) - tr4 = Cc(i-1,3,k) - Cc(i-1,4,k) - tr3 = Cc(i-1,3,k) + Cc(i-1,4,k) - Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3 - Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3 - cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3 - ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3 - cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3 - ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3 + subroutine passb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 5, l1), wa1(*), wa2(*), wa3(*), wa4(*) + real(dp), intent(out) :: ch(ido, l1, 5) + real(dp) :: ci2, ci3, ci4, ci5, cr2, cr3, & + cr4, cr5, di2, di3, di4, di5, dr2, dr3, & + dr4, dr5 + real(dp) :: ti2, ti3, ti4, ti5, tr2, tr3, & + tr4, tr5 + integer :: i, k + real(dp), parameter :: pi = acos(-1.0_dp) + real(dp), parameter :: tr11 = cos(2.0_dp*pi/5.0_dp) + real(dp), parameter :: ti11 = sin(2.0_dp*pi/5.0_dp) + real(dp), parameter :: tr12 = cos(4.0_dp*pi/5.0_dp) + real(dp), parameter :: ti12 = sin(4.0_dp*pi/5.0_dp) + if (ido /= 2) then + do concurrent(i=2:ido:2, k=1:l1) + ti5 = cc(i, 2, k) - cc(i, 5, k) + ti2 = cc(i, 2, k) + cc(i, 5, k) + ti4 = cc(i, 3, k) - cc(i, 4, k) + ti3 = cc(i, 3, k) + cc(i, 4, k) + tr5 = cc(i - 1, 2, k) - cc(i - 1, 5, k) + tr2 = cc(i - 1, 2, k) + cc(i - 1, 5, k) + tr4 = cc(i - 1, 3, k) - cc(i - 1, 4, k) + tr3 = cc(i - 1, 3, k) + cc(i - 1, 4, k) + ch(i - 1, k, 1) = cc(i - 1, 1, k) + tr2 + tr3 + ch(i, k, 1) = cc(i, 1, k) + ti2 + ti3 + cr2 = cc(i - 1, 1, k) + tr11*tr2 + tr12*tr3 + ci2 = cc(i, 1, k) + tr11*ti2 + tr12*ti3 + cr3 = cc(i - 1, 1, k) + tr12*tr2 + tr11*tr3 + ci3 = cc(i, 1, k) + tr12*ti2 + tr11*ti3 cr5 = ti11*tr5 + ti12*tr4 ci5 = ti11*ti5 + ti12*ti4 cr4 = ti12*tr5 - ti11*tr4 @@ -43,44 +43,43 @@ subroutine passb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) dr2 = cr2 - ci5 di5 = ci2 - cr5 di2 = ci2 + cr5 - Ch(i-1,k,2) = Wa1(i-1)*dr2 - Wa1(i)*di2 - Ch(i,k,2) = Wa1(i-1)*di2 + Wa1(i)*dr2 - Ch(i-1,k,3) = Wa2(i-1)*dr3 - Wa2(i)*di3 - Ch(i,k,3) = Wa2(i-1)*di3 + Wa2(i)*dr3 - Ch(i-1,k,4) = Wa3(i-1)*dr4 - Wa3(i)*di4 - Ch(i,k,4) = Wa3(i-1)*di4 + Wa3(i)*dr4 - Ch(i-1,k,5) = Wa4(i-1)*dr5 - Wa4(i)*di5 - Ch(i,k,5) = Wa4(i-1)*di5 + Wa4(i)*dr5 - enddo - enddo - else - do k = 1 , l1 - ti5 = Cc(2,2,k) - Cc(2,5,k) - ti2 = Cc(2,2,k) + Cc(2,5,k) - ti4 = Cc(2,3,k) - Cc(2,4,k) - ti3 = Cc(2,3,k) + Cc(2,4,k) - tr5 = Cc(1,2,k) - Cc(1,5,k) - tr2 = Cc(1,2,k) + Cc(1,5,k) - tr4 = Cc(1,3,k) - Cc(1,4,k) - tr3 = Cc(1,3,k) + Cc(1,4,k) - Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3 - Ch(2,k,1) = Cc(2,1,k) + ti2 + ti3 - cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3 - ci2 = Cc(2,1,k) + tr11*ti2 + tr12*ti3 - cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3 - ci3 = Cc(2,1,k) + tr12*ti2 + tr11*ti3 - cr5 = ti11*tr5 + ti12*tr4 - ci5 = ti11*ti5 + ti12*ti4 - cr4 = ti12*tr5 - ti11*tr4 - ci4 = ti12*ti5 - ti11*ti4 - Ch(1,k,2) = cr2 - ci5 - Ch(1,k,5) = cr2 + ci5 - Ch(2,k,2) = ci2 + cr5 - Ch(2,k,3) = ci3 + cr4 - Ch(1,k,3) = cr3 - ci4 - Ch(1,k,4) = cr3 + ci4 - Ch(2,k,4) = ci3 - cr4 - Ch(2,k,5) = ci2 - cr5 - enddo - end if - end subroutine passb5 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 1)*dr2 - wa1(i)*di2 + ch(i, k, 2) = wa1(i - 1)*di2 + wa1(i)*dr2 + ch(i - 1, k, 3) = wa2(i - 1)*dr3 - wa2(i)*di3 + ch(i, k, 3) = wa2(i - 1)*di3 + wa2(i)*dr3 + ch(i - 1, k, 4) = wa3(i - 1)*dr4 - wa3(i)*di4 + ch(i, k, 4) = wa3(i - 1)*di4 + wa3(i)*dr4 + ch(i - 1, k, 5) = wa4(i - 1)*dr5 - wa4(i)*di5 + ch(i, k, 5) = wa4(i - 1)*di5 + wa4(i)*dr5 + end do + else + do concurrent(k=1:l1) + ti5 = cc(2, 2, k) - cc(2, 5, k) + ti2 = cc(2, 2, k) + cc(2, 5, k) + ti4 = cc(2, 3, k) - cc(2, 4, k) + ti3 = cc(2, 3, k) + cc(2, 4, k) + tr5 = cc(1, 2, k) - cc(1, 5, k) + tr2 = cc(1, 2, k) + cc(1, 5, k) + tr4 = cc(1, 3, k) - cc(1, 4, k) + tr3 = cc(1, 3, k) + cc(1, 4, k) + ch(1, k, 1) = cc(1, 1, k) + tr2 + tr3 + ch(2, k, 1) = cc(2, 1, k) + ti2 + ti3 + cr2 = cc(1, 1, k) + tr11*tr2 + tr12*tr3 + ci2 = cc(2, 1, k) + tr11*ti2 + tr12*ti3 + cr3 = cc(1, 1, k) + tr12*tr2 + tr11*tr3 + ci3 = cc(2, 1, k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + ch(1, k, 2) = cr2 - ci5 + ch(1, k, 5) = cr2 + ci5 + ch(2, k, 2) = ci2 + cr5 + ch(2, k, 3) = ci3 + cr4 + ch(1, k, 3) = cr3 - ci4 + ch(1, k, 4) = cr3 + ci4 + ch(2, k, 4) = ci3 - cr4 + ch(2, k, 5) = ci2 - cr5 + end do + end if + end subroutine passb5 diff --git a/src/passf.f90 b/src/passf.f90 index ebf9278..21bc102 100644 --- a/src/passf.f90 +++ b/src/passf.f90 @@ -1,124 +1,109 @@ - subroutine passf(Nac,Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) - use fftpack_kind - implicit none - real(rk) :: c1 , c2 , Cc , Ch , Ch2 , Wa , wai , war - integer :: i , idij , idj , idl , Idl1 , idlj , Ido , idot , idp , & - ik , inc , Ip , ipp2 , ipph , j , jc , k , l , l1 , lc - integer :: Nac , nt - dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , Wa(*) , & - & c2(Idl1,Ip) , Ch2(Idl1,Ip) - idot = Ido/2 - nt = Ip*Idl1 - ipp2 = Ip + 2 - ipph = (Ip+1)/2 - idp = Ip*Ido + subroutine passf(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(out) :: nac + integer, intent(in) :: ido, ip, l1, idl1 + real(dp), intent(in) :: cc(ido, ip, l1), wa(*) + real(dp), intent(out) :: c1(ido, l1, ip), c2(idl1, ip), ch(ido, l1, ip) + real(dp), intent(inout) :: ch2(idl1, ip) + real(dp) :: wai, war + integer :: i, idij, idj, idl, idlj, idot, idp, & + ik, inc, ipp2, ipph, j, jc, k, l, lc + integer :: nt + idot = ido/2 + nt = ip*idl1 + ipp2 = ip + 2 + ipph = (ip + 1)/2 + idp = ip*ido ! - if ( Ido idp) idlj = idlj - idp + war = wa(idlj - 1) + wai = wa(idlj) + do concurrent(ik=1:idl1) + c2(ik, l) = c2(ik, l) + war*ch2(ik, j) + c2(ik, lc) = c2(ik, lc) - wai*ch2(ik, jc) + end do + end do + end do + do concurrent(ik=1:idl1, j=2:ipph) + ch2(ik, 1) = ch2(ik, 1) + ch2(ik, j) + end do + do concurrent(ik=2:idl1:2, j=2:ipph) jc = ipp2 - j - do i = 1 , Ido - do k = 1 , l1 - Ch(i,k,j) = Cc(i,j,k) + Cc(i,jc,k) - Ch(i,k,jc) = Cc(i,j,k) - Cc(i,jc,k) - enddo - enddo - enddo - do i = 1 , Ido - do k = 1 , l1 - Ch(i,k,1) = Cc(i,1,k) - enddo - enddo - else - do j = 2 , ipph - jc = ipp2 - j - do k = 1 , l1 - do i = 1 , Ido - Ch(i,k,j) = Cc(i,j,k) + Cc(i,jc,k) - Ch(i,k,jc) = Cc(i,j,k) - Cc(i,jc,k) - enddo - enddo - enddo - do k = 1 , l1 - do i = 1 , Ido - Ch(i,k,1) = Cc(i,1,k) - enddo - enddo - endif - idl = 2 - Ido - inc = 0 - do l = 2 , ipph - lc = ipp2 - l - idl = idl + Ido - do ik = 1 , Idl1 - c2(ik,l) = Ch2(ik,1) + Wa(idl-1)*Ch2(ik,2) - c2(ik,lc) = -Wa(idl)*Ch2(ik,Ip) - enddo - idlj = idl - inc = inc + Ido - do j = 3 , ipph - jc = ipp2 - j - idlj = idlj + inc - if ( idlj>idp ) idlj = idlj - idp - war = Wa(idlj-1) - wai = Wa(idlj) - do ik = 1 , Idl1 - c2(ik,l) = c2(ik,l) + war*Ch2(ik,j) - c2(ik,lc) = c2(ik,lc) - wai*Ch2(ik,jc) - enddo - enddo - enddo - do j = 2 , ipph - do ik = 1 , Idl1 - Ch2(ik,1) = Ch2(ik,1) + Ch2(ik,j) - enddo - enddo - do j = 2 , ipph - jc = ipp2 - j - do ik = 2 , Idl1 , 2 - Ch2(ik-1,j) = c2(ik-1,j) - c2(ik,jc) - Ch2(ik-1,jc) = c2(ik-1,j) + c2(ik,jc) - Ch2(ik,j) = c2(ik,j) + c2(ik-1,jc) - Ch2(ik,jc) = c2(ik,j) - c2(ik-1,jc) - enddo - enddo - Nac = 1 - if ( Ido==2 ) return - Nac = 0 - do ik = 1 , Idl1 - c2(ik,1) = Ch2(ik,1) - enddo - do j = 2 , Ip - do k = 1 , l1 - c1(1,k,j) = Ch(1,k,j) - c1(2,k,j) = Ch(2,k,j) - enddo - enddo - if ( idot>l1 ) then - idj = 2 - Ido - do j = 2 , Ip - idj = idj + Ido - do k = 1 , l1 - idij = idj - do i = 4 , Ido , 2 - idij = idij + 2 - c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) + Wa(idij) & - *Ch(i,k,j) - c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) - Wa(idij) & - *Ch(i-1,k,j) - enddo - enddo - enddo - else - idij = 0 - do j = 2 , Ip - idij = idij + 2 - do i = 4 , Ido , 2 + ch2(ik - 1, j) = c2(ik - 1, j) - c2(ik, jc) + ch2(ik - 1, jc) = c2(ik - 1, j) + c2(ik, jc) + ch2(ik, j) = c2(ik, j) + c2(ik - 1, jc) + ch2(ik, jc) = c2(ik, j) - c2(ik - 1, jc) + end do + nac = 1 + if (ido == 2) return + nac = 0 + do concurrent(ik=1:idl1) + c2(ik, 1) = ch2(ik, 1) + end do + do concurrent(j=2:ip, k=1:l1) + c1(1, k, j) = ch(1, k, j) + c1(2, k, j) = ch(2, k, j) + end do + if (idot > l1) then + idj = 2 - ido + do j = 2, ip + idj = idj + ido + do k = 1, l1 + idij = idj + do i = 4, ido, 2 + idij = idij + 2 + c1(i - 1, k, j) = wa(idij - 1)*ch(i - 1, k, j) + wa(idij) & + *ch(i, k, j) + c1(i, k, j) = wa(idij - 1)*ch(i, k, j) - wa(idij) & + *ch(i - 1, k, j) + end do + end do + end do + else + idij = 0 + do j = 2, ip idij = idij + 2 - do k = 1 , l1 - c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) + Wa(idij)*Ch(i,k,j) - c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) - Wa(idij)*Ch(i-1,k,j) - enddo - enddo - enddo - end if - end subroutine passf \ No newline at end of file + do i = 4, ido, 2 + idij = idij + 2 + do concurrent(k=1:l1) + c1(i - 1, k, j) = wa(idij - 1)*ch(i - 1, k, j) + wa(idij)*ch(i, k, j) + c1(i, k, j) = wa(idij - 1)*ch(i, k, j) - wa(idij)*ch(i - 1, k, j) + end do + end do + end do + end if + end subroutine passf diff --git a/src/passf2.f90 b/src/passf2.f90 index 4c9f17c..eb96540 100644 --- a/src/passf2.f90 +++ b/src/passf2.f90 @@ -1,26 +1,26 @@ - subroutine passf2(Ido,l1,Cc,Ch,Wa1) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 - integer :: i , Ido , k , l1 - dimension Cc(Ido,2,l1) , Ch(Ido,l1,2) , Wa1(*) - if ( Ido>2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - Ch(i-1,k,1) = Cc(i-1,1,k) + Cc(i-1,2,k) - tr2 = Cc(i-1,1,k) - Cc(i-1,2,k) - Ch(i,k,1) = Cc(i,1,k) + Cc(i,2,k) - ti2 = Cc(i,1,k) - Cc(i,2,k) - Ch(i,k,2) = Wa1(i-1)*ti2 - Wa1(i)*tr2 - Ch(i-1,k,2) = Wa1(i-1)*tr2 + Wa1(i)*ti2 - enddo - enddo - else - do k = 1 , l1 - Ch(1,k,1) = Cc(1,1,k) + Cc(1,2,k) - Ch(1,k,2) = Cc(1,1,k) - Cc(1,2,k) - Ch(2,k,1) = Cc(2,1,k) + Cc(2,2,k) - Ch(2,k,2) = Cc(2,1,k) - Cc(2,2,k) - enddo - end if - end subroutine passf2 \ No newline at end of file + subroutine passf2(ido, l1, cc, ch, wa1) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 2, l1), wa1(*) + real(dp), intent(out) :: ch(ido, l1, 2) + real(dp) :: ti2, tr2 + integer :: i, k + if (ido > 2) then + do concurrent(i=2:ido:2, k=1:l1) + ch(i - 1, k, 1) = cc(i - 1, 1, k) + cc(i - 1, 2, k) + tr2 = cc(i - 1, 1, k) - cc(i - 1, 2, k) + ch(i, k, 1) = cc(i, 1, k) + cc(i, 2, k) + ti2 = cc(i, 1, k) - cc(i, 2, k) + ch(i, k, 2) = wa1(i - 1)*ti2 - wa1(i)*tr2 + ch(i - 1, k, 2) = wa1(i - 1)*tr2 + wa1(i)*ti2 + end do + else + do concurrent(k=1:l1) + ch(1, k, 1) = cc(1, 1, k) + cc(1, 2, k) + ch(1, k, 2) = cc(1, 1, k) - cc(1, 2, k) + ch(2, k, 1) = cc(2, 1, k) + cc(2, 2, k) + ch(2, k, 2) = cc(2, 1, k) - cc(2, 2, k) + end do + end if + end subroutine passf2 diff --git a/src/passf3.f90 b/src/passf3.f90 index 165b286..19f78ff 100644 --- a/src/passf3.f90 +++ b/src/passf3.f90 @@ -1,47 +1,47 @@ - subroutine passf3(Ido,l1,Cc,Ch,Wa1,Wa2) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , cr2 , cr3 , di2 , di3 , & - & dr2 , dr3 , ti2 , tr2 , Wa1 , Wa2 - integer :: i , Ido , k , l1 - dimension Cc(Ido,3,l1) , Ch(Ido,l1,3) , Wa1(*) , Wa2(*) - real(rk),parameter :: taur = -0.5_rk - real(rk),parameter :: taui = -sqrt(3.0_rk) / 2.0_rk - if ( Ido/=2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - tr2 = Cc(i-1,2,k) + Cc(i-1,3,k) - cr2 = Cc(i-1,1,k) + taur*tr2 - Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 - ti2 = Cc(i,2,k) + Cc(i,3,k) - ci2 = Cc(i,1,k) + taur*ti2 - Ch(i,k,1) = Cc(i,1,k) + ti2 - cr3 = taui*(Cc(i-1,2,k)-Cc(i-1,3,k)) - ci3 = taui*(Cc(i,2,k)-Cc(i,3,k)) + subroutine passf3(ido, l1, cc, ch, wa1, wa2) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 3, l1), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, l1, 3) + real(dp) :: ci2, ci3, cr2, cr3, di2, di3, & + & dr2, dr3, ti2, tr2 + integer :: i, k + real(dp), parameter :: taur = -0.5_dp + real(dp), parameter :: taui = -sqrt(3.0_dp)/2.0_dp + if (ido /= 2) then + do concurrent(i=2:ido:2, k=1:l1) + tr2 = cc(i - 1, 2, k) + cc(i - 1, 3, k) + cr2 = cc(i - 1, 1, k) + taur*tr2 + ch(i - 1, k, 1) = cc(i - 1, 1, k) + tr2 + ti2 = cc(i, 2, k) + cc(i, 3, k) + ci2 = cc(i, 1, k) + taur*ti2 + ch(i, k, 1) = cc(i, 1, k) + ti2 + cr3 = taui*(cc(i - 1, 2, k) - cc(i - 1, 3, k)) + ci3 = taui*(cc(i, 2, k) - cc(i, 3, k)) dr2 = cr2 - ci3 dr3 = cr2 + ci3 di2 = ci2 + cr3 di3 = ci2 - cr3 - Ch(i,k,2) = Wa1(i-1)*di2 - Wa1(i)*dr2 - Ch(i-1,k,2) = Wa1(i-1)*dr2 + Wa1(i)*di2 - Ch(i,k,3) = Wa2(i-1)*di3 - Wa2(i)*dr3 - Ch(i-1,k,3) = Wa2(i-1)*dr3 + Wa2(i)*di3 - enddo - enddo - else - do k = 1 , l1 - tr2 = Cc(1,2,k) + Cc(1,3,k) - cr2 = Cc(1,1,k) + taur*tr2 - Ch(1,k,1) = Cc(1,1,k) + tr2 - ti2 = Cc(2,2,k) + Cc(2,3,k) - ci2 = Cc(2,1,k) + taur*ti2 - Ch(2,k,1) = Cc(2,1,k) + ti2 - cr3 = taui*(Cc(1,2,k)-Cc(1,3,k)) - ci3 = taui*(Cc(2,2,k)-Cc(2,3,k)) - Ch(1,k,2) = cr2 - ci3 - Ch(1,k,3) = cr2 + ci3 - Ch(2,k,2) = ci2 + cr3 - Ch(2,k,3) = ci2 - cr3 - enddo - end if - end subroutine passf3 \ No newline at end of file + ch(i, k, 2) = wa1(i - 1)*di2 - wa1(i)*dr2 + ch(i - 1, k, 2) = wa1(i - 1)*dr2 + wa1(i)*di2 + ch(i, k, 3) = wa2(i - 1)*di3 - wa2(i)*dr3 + ch(i - 1, k, 3) = wa2(i - 1)*dr3 + wa2(i)*di3 + end do + else + do concurrent(k=1:l1) + tr2 = cc(1, 2, k) + cc(1, 3, k) + cr2 = cc(1, 1, k) + taur*tr2 + ch(1, k, 1) = cc(1, 1, k) + tr2 + ti2 = cc(2, 2, k) + cc(2, 3, k) + ci2 = cc(2, 1, k) + taur*ti2 + ch(2, k, 1) = cc(2, 1, k) + ti2 + cr3 = taui*(cc(1, 2, k) - cc(1, 3, k)) + ci3 = taui*(cc(2, 2, k) - cc(2, 3, k)) + ch(1, k, 2) = cr2 - ci3 + ch(1, k, 3) = cr2 + ci3 + ch(2, k, 2) = ci2 + cr3 + ch(2, k, 3) = ci2 - cr3 + end do + end if + end subroutine passf3 diff --git a/src/passf4.f90 b/src/passf4.f90 index 110dea8..4deaa6f 100644 --- a/src/passf4.f90 +++ b/src/passf4.f90 @@ -1,56 +1,55 @@ - subroutine passf4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & - & ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3 , tr4 , & - & Wa1 , Wa2 , Wa3 - integer :: i , Ido , k , l1 - dimension Cc(Ido,4,l1) , Ch(Ido,l1,4) , Wa1(*) , Wa2(*) , Wa3(*) - if ( Ido/=2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - ti1 = Cc(i,1,k) - Cc(i,3,k) - ti2 = Cc(i,1,k) + Cc(i,3,k) - ti3 = Cc(i,2,k) + Cc(i,4,k) - tr4 = Cc(i,2,k) - Cc(i,4,k) - tr1 = Cc(i-1,1,k) - Cc(i-1,3,k) - tr2 = Cc(i-1,1,k) + Cc(i-1,3,k) - ti4 = Cc(i-1,4,k) - Cc(i-1,2,k) - tr3 = Cc(i-1,2,k) + Cc(i-1,4,k) - Ch(i-1,k,1) = tr2 + tr3 + subroutine passf4(ido, l1, cc, ch, wa1, wa2, wa3) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 4, l1), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, l1, 4) + real(dp) :: ci2, ci3, ci4, cr2, cr3, cr4, & + & ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4 + integer :: i, k + if (ido /= 2) then + do concurrent(i=2:ido:2, k=1:l1) + ti1 = cc(i, 1, k) - cc(i, 3, k) + ti2 = cc(i, 1, k) + cc(i, 3, k) + ti3 = cc(i, 2, k) + cc(i, 4, k) + tr4 = cc(i, 2, k) - cc(i, 4, k) + tr1 = cc(i - 1, 1, k) - cc(i - 1, 3, k) + tr2 = cc(i - 1, 1, k) + cc(i - 1, 3, k) + ti4 = cc(i - 1, 4, k) - cc(i - 1, 2, k) + tr3 = cc(i - 1, 2, k) + cc(i - 1, 4, k) + ch(i - 1, k, 1) = tr2 + tr3 cr3 = tr2 - tr3 - Ch(i,k,1) = ti2 + ti3 + ch(i, k, 1) = ti2 + ti3 ci3 = ti2 - ti3 cr2 = tr1 + tr4 cr4 = tr1 - tr4 ci2 = ti1 + ti4 ci4 = ti1 - ti4 - Ch(i-1,k,2) = Wa1(i-1)*cr2 + Wa1(i)*ci2 - Ch(i,k,2) = Wa1(i-1)*ci2 - Wa1(i)*cr2 - Ch(i-1,k,3) = Wa2(i-1)*cr3 + Wa2(i)*ci3 - Ch(i,k,3) = Wa2(i-1)*ci3 - Wa2(i)*cr3 - Ch(i-1,k,4) = Wa3(i-1)*cr4 + Wa3(i)*ci4 - Ch(i,k,4) = Wa3(i-1)*ci4 - Wa3(i)*cr4 - enddo - enddo - else - do k = 1 , l1 - ti1 = Cc(2,1,k) - Cc(2,3,k) - ti2 = Cc(2,1,k) + Cc(2,3,k) - tr4 = Cc(2,2,k) - Cc(2,4,k) - ti3 = Cc(2,2,k) + Cc(2,4,k) - tr1 = Cc(1,1,k) - Cc(1,3,k) - tr2 = Cc(1,1,k) + Cc(1,3,k) - ti4 = Cc(1,4,k) - Cc(1,2,k) - tr3 = Cc(1,2,k) + Cc(1,4,k) - Ch(1,k,1) = tr2 + tr3 - Ch(1,k,3) = tr2 - tr3 - Ch(2,k,1) = ti2 + ti3 - Ch(2,k,3) = ti2 - ti3 - Ch(1,k,2) = tr1 + tr4 - Ch(1,k,4) = tr1 - tr4 - Ch(2,k,2) = ti1 + ti4 - Ch(2,k,4) = ti1 - ti4 - enddo - end if - end subroutine passf4 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 1)*cr2 + wa1(i)*ci2 + ch(i, k, 2) = wa1(i - 1)*ci2 - wa1(i)*cr2 + ch(i - 1, k, 3) = wa2(i - 1)*cr3 + wa2(i)*ci3 + ch(i, k, 3) = wa2(i - 1)*ci3 - wa2(i)*cr3 + ch(i - 1, k, 4) = wa3(i - 1)*cr4 + wa3(i)*ci4 + ch(i, k, 4) = wa3(i - 1)*ci4 - wa3(i)*cr4 + end do + else + do concurrent(k=1:l1) + ti1 = cc(2, 1, k) - cc(2, 3, k) + ti2 = cc(2, 1, k) + cc(2, 3, k) + tr4 = cc(2, 2, k) - cc(2, 4, k) + ti3 = cc(2, 2, k) + cc(2, 4, k) + tr1 = cc(1, 1, k) - cc(1, 3, k) + tr2 = cc(1, 1, k) + cc(1, 3, k) + ti4 = cc(1, 4, k) - cc(1, 2, k) + tr3 = cc(1, 2, k) + cc(1, 4, k) + ch(1, k, 1) = tr2 + tr3 + ch(1, k, 3) = tr2 - tr3 + ch(2, k, 1) = ti2 + ti3 + ch(2, k, 3) = ti2 - ti3 + ch(1, k, 2) = tr1 + tr4 + ch(1, k, 4) = tr1 - tr4 + ch(2, k, 2) = ti1 + ti4 + ch(2, k, 4) = ti1 - ti4 + end do + end if + end subroutine passf4 diff --git a/src/passf5.f90 b/src/passf5.f90 index 6c29b44..a209649 100644 --- a/src/passf5.f90 +++ b/src/passf5.f90 @@ -1,36 +1,36 @@ - subroutine passf5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & - cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & - dr4 , dr5 - real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & - tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 - integer :: i , Ido , k , l1 - dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), & - Wa4(*) - real(rk),parameter :: pi = acos(-1.0_rk) - real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti11 = -sin(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti12 = -sin(4.0_rk * pi / 5.0_rk) - if ( Ido/=2 ) then - do k = 1 , l1 - do i = 2 , Ido , 2 - ti5 = Cc(i,2,k) - Cc(i,5,k) - ti2 = Cc(i,2,k) + Cc(i,5,k) - ti4 = Cc(i,3,k) - Cc(i,4,k) - ti3 = Cc(i,3,k) + Cc(i,4,k) - tr5 = Cc(i-1,2,k) - Cc(i-1,5,k) - tr2 = Cc(i-1,2,k) + Cc(i-1,5,k) - tr4 = Cc(i-1,3,k) - Cc(i-1,4,k) - tr3 = Cc(i-1,3,k) + Cc(i-1,4,k) - Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3 - Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3 - cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3 - ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3 - cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3 - ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3 + subroutine passf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 5, l1), wa1(*), wa2(*), wa3(*), wa4(*) + real(dp), intent(out) :: ch(ido, l1, 5) + real(dp) :: ci2, ci3, ci4, ci5, cr2, cr3, & + cr4, cr5, di2, di3, di4, di5, dr2, dr3, & + dr4, dr5 + real(dp) :: ti2, ti3, ti4, ti5, tr2, tr3, & + tr4, tr5 + integer :: i, k + real(dp), parameter :: pi = acos(-1.0_dp) + real(dp), parameter :: tr11 = cos(2.0_dp*pi/5.0_dp) + real(dp), parameter :: ti11 = -sin(2.0_dp*pi/5.0_dp) + real(dp), parameter :: tr12 = cos(4.0_dp*pi/5.0_dp) + real(dp), parameter :: ti12 = -sin(4.0_dp*pi/5.0_dp) + if (ido /= 2) then + do concurrent(i=2:ido:2, k=1:l1) + ti5 = cc(i, 2, k) - cc(i, 5, k) + ti2 = cc(i, 2, k) + cc(i, 5, k) + ti4 = cc(i, 3, k) - cc(i, 4, k) + ti3 = cc(i, 3, k) + cc(i, 4, k) + tr5 = cc(i - 1, 2, k) - cc(i - 1, 5, k) + tr2 = cc(i - 1, 2, k) + cc(i - 1, 5, k) + tr4 = cc(i - 1, 3, k) - cc(i - 1, 4, k) + tr3 = cc(i - 1, 3, k) + cc(i - 1, 4, k) + ch(i - 1, k, 1) = cc(i - 1, 1, k) + tr2 + tr3 + ch(i, k, 1) = cc(i, 1, k) + ti2 + ti3 + cr2 = cc(i - 1, 1, k) + tr11*tr2 + tr12*tr3 + ci2 = cc(i, 1, k) + tr11*ti2 + tr12*ti3 + cr3 = cc(i - 1, 1, k) + tr12*tr2 + tr11*tr3 + ci3 = cc(i, 1, k) + tr12*ti2 + tr11*ti3 cr5 = ti11*tr5 + ti12*tr4 ci5 = ti11*ti5 + ti12*ti4 cr4 = ti12*tr5 - ti11*tr4 @@ -43,44 +43,43 @@ subroutine passf5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) dr2 = cr2 - ci5 di5 = ci2 - cr5 di2 = ci2 + cr5 - Ch(i-1,k,2) = Wa1(i-1)*dr2 + Wa1(i)*di2 - Ch(i,k,2) = Wa1(i-1)*di2 - Wa1(i)*dr2 - Ch(i-1,k,3) = Wa2(i-1)*dr3 + Wa2(i)*di3 - Ch(i,k,3) = Wa2(i-1)*di3 - Wa2(i)*dr3 - Ch(i-1,k,4) = Wa3(i-1)*dr4 + Wa3(i)*di4 - Ch(i,k,4) = Wa3(i-1)*di4 - Wa3(i)*dr4 - Ch(i-1,k,5) = Wa4(i-1)*dr5 + Wa4(i)*di5 - Ch(i,k,5) = Wa4(i-1)*di5 - Wa4(i)*dr5 - enddo - enddo - else - do k = 1 , l1 - ti5 = Cc(2,2,k) - Cc(2,5,k) - ti2 = Cc(2,2,k) + Cc(2,5,k) - ti4 = Cc(2,3,k) - Cc(2,4,k) - ti3 = Cc(2,3,k) + Cc(2,4,k) - tr5 = Cc(1,2,k) - Cc(1,5,k) - tr2 = Cc(1,2,k) + Cc(1,5,k) - tr4 = Cc(1,3,k) - Cc(1,4,k) - tr3 = Cc(1,3,k) + Cc(1,4,k) - Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3 - Ch(2,k,1) = Cc(2,1,k) + ti2 + ti3 - cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3 - ci2 = Cc(2,1,k) + tr11*ti2 + tr12*ti3 - cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3 - ci3 = Cc(2,1,k) + tr12*ti2 + tr11*ti3 - cr5 = ti11*tr5 + ti12*tr4 - ci5 = ti11*ti5 + ti12*ti4 - cr4 = ti12*tr5 - ti11*tr4 - ci4 = ti12*ti5 - ti11*ti4 - Ch(1,k,2) = cr2 - ci5 - Ch(1,k,5) = cr2 + ci5 - Ch(2,k,2) = ci2 + cr5 - Ch(2,k,3) = ci3 + cr4 - Ch(1,k,3) = cr3 - ci4 - Ch(1,k,4) = cr3 + ci4 - Ch(2,k,4) = ci3 - cr4 - Ch(2,k,5) = ci2 - cr5 - enddo - end if - end subroutine passf5 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 1)*dr2 + wa1(i)*di2 + ch(i, k, 2) = wa1(i - 1)*di2 - wa1(i)*dr2 + ch(i - 1, k, 3) = wa2(i - 1)*dr3 + wa2(i)*di3 + ch(i, k, 3) = wa2(i - 1)*di3 - wa2(i)*dr3 + ch(i - 1, k, 4) = wa3(i - 1)*dr4 + wa3(i)*di4 + ch(i, k, 4) = wa3(i - 1)*di4 - wa3(i)*dr4 + ch(i - 1, k, 5) = wa4(i - 1)*dr5 + wa4(i)*di5 + ch(i, k, 5) = wa4(i - 1)*di5 - wa4(i)*dr5 + end do + else + do concurrent(k=1:l1) + ti5 = cc(2, 2, k) - cc(2, 5, k) + ti2 = cc(2, 2, k) + cc(2, 5, k) + ti4 = cc(2, 3, k) - cc(2, 4, k) + ti3 = cc(2, 3, k) + cc(2, 4, k) + tr5 = cc(1, 2, k) - cc(1, 5, k) + tr2 = cc(1, 2, k) + cc(1, 5, k) + tr4 = cc(1, 3, k) - cc(1, 4, k) + tr3 = cc(1, 3, k) + cc(1, 4, k) + ch(1, k, 1) = cc(1, 1, k) + tr2 + tr3 + ch(2, k, 1) = cc(2, 1, k) + ti2 + ti3 + cr2 = cc(1, 1, k) + tr11*tr2 + tr12*tr3 + ci2 = cc(2, 1, k) + tr11*ti2 + tr12*ti3 + cr3 = cc(1, 1, k) + tr12*tr2 + tr11*tr3 + ci3 = cc(2, 1, k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + ch(1, k, 2) = cr2 - ci5 + ch(1, k, 5) = cr2 + ci5 + ch(2, k, 2) = ci2 + cr5 + ch(2, k, 3) = ci3 + cr4 + ch(1, k, 3) = cr3 - ci4 + ch(1, k, 4) = cr3 + ci4 + ch(2, k, 4) = ci3 - cr4 + ch(2, k, 5) = ci2 - cr5 + end do + end if + end subroutine passf5 diff --git a/src/radb2.f90 b/src/radb2.f90 index b37708a..7fb4c55 100644 --- a/src/radb2.f90 +++ b/src/radb2.f90 @@ -1,31 +1,31 @@ - subroutine radb2(Ido,l1,Cc,Ch,Wa1) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Cc(Ido,2,l1) , Ch(Ido,l1,2) , Wa1(*) - do k = 1 , l1 - Ch(1,k,1) = Cc(1,1,k) + Cc(Ido,2,k) - Ch(1,k,2) = Cc(1,1,k) - Cc(Ido,2,k) - enddo - if ( Ido<2 ) return - if ( Ido/=2 ) then - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radb2(ido, l1, cc, ch, wa1) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 2, l1), wa1(*) + real(dp), intent(out) :: ch(ido, l1, 2) + real(dp) :: ti2, tr2 + integer :: i, ic, idp2, k + do concurrent(k=1:l1) + ch(1, k, 1) = cc(1, 1, k) + cc(ido, 2, k) + ch(1, k, 2) = cc(1, 1, k) - cc(ido, 2, k) + end do + if (ido < 2) return + if (ido /= 2) then + idp2 = ido + 2 + do concurrent(k=1:l1, i=3:ido:2) ic = idp2 - i - Ch(i-1,k,1) = Cc(i-1,1,k) + Cc(ic-1,2,k) - tr2 = Cc(i-1,1,k) - Cc(ic-1,2,k) - Ch(i,k,1) = Cc(i,1,k) - Cc(ic,2,k) - ti2 = Cc(i,1,k) + Cc(ic,2,k) - Ch(i-1,k,2) = Wa1(i-2)*tr2 - Wa1(i-1)*ti2 - Ch(i,k,2) = Wa1(i-2)*ti2 + Wa1(i-1)*tr2 - enddo - enddo - if ( mod(Ido,2)==1 ) return - endif - do k = 1 , l1 - Ch(Ido,k,1) = Cc(Ido,1,k) + Cc(Ido,1,k) - Ch(Ido,k,2) = -(Cc(1,2,k)+Cc(1,2,k)) - enddo - end subroutine radb2 \ No newline at end of file + ch(i - 1, k, 1) = cc(i - 1, 1, k) + cc(ic - 1, 2, k) + tr2 = cc(i - 1, 1, k) - cc(ic - 1, 2, k) + ch(i, k, 1) = cc(i, 1, k) - cc(ic, 2, k) + ti2 = cc(i, 1, k) + cc(ic, 2, k) + ch(i - 1, k, 2) = wa1(i - 2)*tr2 - wa1(i - 1)*ti2 + ch(i, k, 2) = wa1(i - 2)*ti2 + wa1(i - 1)*tr2 + end do + if (mod(ido, 2) == 1) return + end if + do concurrent(k=1:l1) + ch(ido, k, 1) = cc(ido, 1, k) + cc(ido, 1, k) + ch(ido, k, 2) = -(cc(1, 2, k) + cc(1, 2, k)) + end do + end subroutine radb2 diff --git a/src/radb3.f90 b/src/radb3.f90 index f2c3e23..6b9cba6 100644 --- a/src/radb3.f90 +++ b/src/radb3.f90 @@ -1,41 +1,41 @@ - subroutine radb3(Ido,l1,Cc,Ch,Wa1,Wa2) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , cr2 , cr3 , di2 , di3 , & - dr2 , dr3 , ti2 , tr2 , Wa1 , Wa2 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Cc(Ido,3,l1) , Ch(Ido,l1,3) , Wa1(*) , Wa2(*) - real(rk),parameter :: taur = - 0.5_rk - real(rk),parameter :: taui = sqrt(3.0_rk) / 2.0_rk - do k = 1 , l1 - tr2 = Cc(Ido,2,k) + Cc(Ido,2,k) - cr2 = Cc(1,1,k) + taur*tr2 - Ch(1,k,1) = Cc(1,1,k) + tr2 - ci3 = taui*(Cc(1,3,k)+Cc(1,3,k)) - Ch(1,k,2) = cr2 - ci3 - Ch(1,k,3) = cr2 + ci3 - enddo - if ( Ido==1 ) return - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radb3(ido, l1, cc, ch, wa1, wa2) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 3, l1), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, l1, 3) + real(dp) :: ci2, ci3, cr2, cr3, di2, di3, & + dr2, dr3, ti2, tr2 + integer :: i, ic, idp2, k + real(dp), parameter :: taur = -0.5_dp + real(dp), parameter :: taui = sqrt(3.0_dp)/2.0_dp + do concurrent(k=1:l1) + tr2 = cc(ido, 2, k) + cc(ido, 2, k) + cr2 = cc(1, 1, k) + taur*tr2 + ch(1, k, 1) = cc(1, 1, k) + tr2 + ci3 = taui*(cc(1, 3, k) + cc(1, 3, k)) + ch(1, k, 2) = cr2 - ci3 + ch(1, k, 3) = cr2 + ci3 + end do + if (ido == 1) return + idp2 = ido + 2 + do concurrent(k=1:l1, i=3:ido:2) ic = idp2 - i - tr2 = Cc(i-1,3,k) + Cc(ic-1,2,k) - cr2 = Cc(i-1,1,k) + taur*tr2 - Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 - ti2 = Cc(i,3,k) - Cc(ic,2,k) - ci2 = Cc(i,1,k) + taur*ti2 - Ch(i,k,1) = Cc(i,1,k) + ti2 - cr3 = taui*(Cc(i-1,3,k)-Cc(ic-1,2,k)) - ci3 = taui*(Cc(i,3,k)+Cc(ic,2,k)) + tr2 = cc(i - 1, 3, k) + cc(ic - 1, 2, k) + cr2 = cc(i - 1, 1, k) + taur*tr2 + ch(i - 1, k, 1) = cc(i - 1, 1, k) + tr2 + ti2 = cc(i, 3, k) - cc(ic, 2, k) + ci2 = cc(i, 1, k) + taur*ti2 + ch(i, k, 1) = cc(i, 1, k) + ti2 + cr3 = taui*(cc(i - 1, 3, k) - cc(ic - 1, 2, k)) + ci3 = taui*(cc(i, 3, k) + cc(ic, 2, k)) dr2 = cr2 - ci3 dr3 = cr2 + ci3 di2 = ci2 + cr3 di3 = ci2 - cr3 - Ch(i-1,k,2) = Wa1(i-2)*dr2 - Wa1(i-1)*di2 - Ch(i,k,2) = Wa1(i-2)*di2 + Wa1(i-1)*dr2 - Ch(i-1,k,3) = Wa2(i-2)*dr3 - Wa2(i-1)*di3 - Ch(i,k,3) = Wa2(i-2)*di3 + Wa2(i-1)*dr3 - enddo - enddo - end subroutine radb3 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 2)*dr2 - wa1(i - 1)*di2 + ch(i, k, 2) = wa1(i - 2)*di2 + wa1(i - 1)*dr2 + ch(i - 1, k, 3) = wa2(i - 2)*dr3 - wa2(i - 1)*di3 + ch(i, k, 3) = wa2(i - 2)*di3 + wa2(i - 1)*dr3 + end do + end subroutine radb3 diff --git a/src/radb4.f90 b/src/radb4.f90 index e6fadf5..4291e75 100644 --- a/src/radb4.f90 +++ b/src/radb4.f90 @@ -1,62 +1,62 @@ - subroutine radb4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & - ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3, & - tr4 , Wa1 , Wa2 , Wa3 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Cc(Ido,4,l1) , Ch(Ido,l1,4) , Wa1(*) , Wa2(*) , Wa3(*) - real(rk),parameter :: sqrt2 = sqrt(2.0_rk) - do k = 1 , l1 - tr1 = Cc(1,1,k) - Cc(Ido,4,k) - tr2 = Cc(1,1,k) + Cc(Ido,4,k) - tr3 = Cc(Ido,2,k) + Cc(Ido,2,k) - tr4 = Cc(1,3,k) + Cc(1,3,k) - Ch(1,k,1) = tr2 + tr3 - Ch(1,k,2) = tr1 - tr4 - Ch(1,k,3) = tr2 - tr3 - Ch(1,k,4) = tr1 + tr4 - enddo - if ( Ido<2 ) return - if ( Ido/=2 ) then - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radb4(ido, l1, cc, ch, wa1, wa2, wa3) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 4, l1), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, l1, 4) + real(dp) :: ci2, ci3, ci4, cr2, cr3, cr4, & + ti1, ti2, ti3, ti4, tr1, tr2, tr3, & + tr4 + integer :: i, ic, idp2, k + real(dp), parameter :: sqrt2 = sqrt(2.0_dp) + do concurrent(k=1:l1) + tr1 = cc(1, 1, k) - cc(ido, 4, k) + tr2 = cc(1, 1, k) + cc(ido, 4, k) + tr3 = cc(ido, 2, k) + cc(ido, 2, k) + tr4 = cc(1, 3, k) + cc(1, 3, k) + ch(1, k, 1) = tr2 + tr3 + ch(1, k, 2) = tr1 - tr4 + ch(1, k, 3) = tr2 - tr3 + ch(1, k, 4) = tr1 + tr4 + end do + if (ido < 2) return + if (ido /= 2) then + idp2 = ido + 2 + do concurrent(k=1:l1, i=3:ido:2) ic = idp2 - i - ti1 = Cc(i,1,k) + Cc(ic,4,k) - ti2 = Cc(i,1,k) - Cc(ic,4,k) - ti3 = Cc(i,3,k) - Cc(ic,2,k) - tr4 = Cc(i,3,k) + Cc(ic,2,k) - tr1 = Cc(i-1,1,k) - Cc(ic-1,4,k) - tr2 = Cc(i-1,1,k) + Cc(ic-1,4,k) - ti4 = Cc(i-1,3,k) - Cc(ic-1,2,k) - tr3 = Cc(i-1,3,k) + Cc(ic-1,2,k) - Ch(i-1,k,1) = tr2 + tr3 + ti1 = cc(i, 1, k) + cc(ic, 4, k) + ti2 = cc(i, 1, k) - cc(ic, 4, k) + ti3 = cc(i, 3, k) - cc(ic, 2, k) + tr4 = cc(i, 3, k) + cc(ic, 2, k) + tr1 = cc(i - 1, 1, k) - cc(ic - 1, 4, k) + tr2 = cc(i - 1, 1, k) + cc(ic - 1, 4, k) + ti4 = cc(i - 1, 3, k) - cc(ic - 1, 2, k) + tr3 = cc(i - 1, 3, k) + cc(ic - 1, 2, k) + ch(i - 1, k, 1) = tr2 + tr3 cr3 = tr2 - tr3 - Ch(i,k,1) = ti2 + ti3 + ch(i, k, 1) = ti2 + ti3 ci3 = ti2 - ti3 cr2 = tr1 - tr4 cr4 = tr1 + tr4 ci2 = ti1 + ti4 ci4 = ti1 - ti4 - Ch(i-1,k,2) = Wa1(i-2)*cr2 - Wa1(i-1)*ci2 - Ch(i,k,2) = Wa1(i-2)*ci2 + Wa1(i-1)*cr2 - Ch(i-1,k,3) = Wa2(i-2)*cr3 - Wa2(i-1)*ci3 - Ch(i,k,3) = Wa2(i-2)*ci3 + Wa2(i-1)*cr3 - Ch(i-1,k,4) = Wa3(i-2)*cr4 - Wa3(i-1)*ci4 - Ch(i,k,4) = Wa3(i-2)*ci4 + Wa3(i-1)*cr4 - enddo - enddo - if ( mod(Ido,2)==1 ) return - endif - do k = 1 , l1 - ti1 = Cc(1,2,k) + Cc(1,4,k) - ti2 = Cc(1,4,k) - Cc(1,2,k) - tr1 = Cc(Ido,1,k) - Cc(Ido,3,k) - tr2 = Cc(Ido,1,k) + Cc(Ido,3,k) - Ch(Ido,k,1) = tr2 + tr2 - Ch(Ido,k,2) = sqrt2*(tr1-ti1) - Ch(Ido,k,3) = ti2 + ti2 - Ch(Ido,k,4) = -sqrt2*(tr1+ti1) - enddo - end subroutine radb4 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 2)*cr2 - wa1(i - 1)*ci2 + ch(i, k, 2) = wa1(i - 2)*ci2 + wa1(i - 1)*cr2 + ch(i - 1, k, 3) = wa2(i - 2)*cr3 - wa2(i - 1)*ci3 + ch(i, k, 3) = wa2(i - 2)*ci3 + wa2(i - 1)*cr3 + ch(i - 1, k, 4) = wa3(i - 2)*cr4 - wa3(i - 1)*ci4 + ch(i, k, 4) = wa3(i - 2)*ci4 + wa3(i - 1)*cr4 + end do + if (mod(ido, 2) == 1) return + end if + do concurrent(k=1:l1) + ti1 = cc(1, 2, k) + cc(1, 4, k) + ti2 = cc(1, 4, k) - cc(1, 2, k) + tr1 = cc(ido, 1, k) - cc(ido, 3, k) + tr2 = cc(ido, 1, k) + cc(ido, 3, k) + ch(ido, k, 1) = tr2 + tr2 + ch(ido, k, 2) = sqrt2*(tr1 - ti1) + ch(ido, k, 3) = ti2 + ti2 + ch(ido, k, 4) = -sqrt2*(tr1 + ti1) + end do + end subroutine radb4 diff --git a/src/radb5.f90 b/src/radb5.f90 index e90a4d8..c4de1f6 100644 --- a/src/radb5.f90 +++ b/src/radb5.f90 @@ -1,53 +1,53 @@ - subroutine radb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & - cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & - dr4 , dr5 - real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & - tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), & - Wa4(*) - real(rk),parameter :: pi = acos(-1.0_rk) - real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk) - do k = 1 , l1 - ti5 = Cc(1,3,k) + Cc(1,3,k) - ti4 = Cc(1,5,k) + Cc(1,5,k) - tr2 = Cc(Ido,2,k) + Cc(Ido,2,k) - tr3 = Cc(Ido,4,k) + Cc(Ido,4,k) - Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3 - cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3 - cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3 - ci5 = ti11*ti5 + ti12*ti4 - ci4 = ti12*ti5 - ti11*ti4 - Ch(1,k,2) = cr2 - ci5 - Ch(1,k,3) = cr3 - ci4 - Ch(1,k,4) = cr3 + ci4 - Ch(1,k,5) = cr2 + ci5 - enddo - if ( Ido==1 ) return - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 5, l1), wa1(*), wa2(*), wa3(*), wa4(*) + real(dp), intent(out) :: ch(ido, l1, 5) + real(dp) :: ci2, ci3, ci4, ci5, cr2, cr3, & + cr4, cr5, di2, di3, di4, di5, dr2, dr3, & + dr4, dr5 + real(dp) :: ti2, ti3, ti4, ti5, tr2, tr3, & + tr4, tr5 + integer :: i, ic, idp2, k + real(dp), parameter :: pi = acos(-1.0_dp) + real(dp), parameter :: tr11 = cos(2.0_dp*pi/5.0_dp) + real(dp), parameter :: ti11 = sin(2.0_dp*pi/5.0_dp) + real(dp), parameter :: tr12 = cos(4.0_dp*pi/5.0_dp) + real(dp), parameter :: ti12 = sin(4.0_dp*pi/5.0_dp) + do concurrent(k=1:l1) + ti5 = cc(1, 3, k) + cc(1, 3, k) + ti4 = cc(1, 5, k) + cc(1, 5, k) + tr2 = cc(ido, 2, k) + cc(ido, 2, k) + tr3 = cc(ido, 4, k) + cc(ido, 4, k) + ch(1, k, 1) = cc(1, 1, k) + tr2 + tr3 + cr2 = cc(1, 1, k) + tr11*tr2 + tr12*tr3 + cr3 = cc(1, 1, k) + tr12*tr2 + tr11*tr3 + ci5 = ti11*ti5 + ti12*ti4 + ci4 = ti12*ti5 - ti11*ti4 + ch(1, k, 2) = cr2 - ci5 + ch(1, k, 3) = cr3 - ci4 + ch(1, k, 4) = cr3 + ci4 + ch(1, k, 5) = cr2 + ci5 + end do + if (ido == 1) return + idp2 = ido + 2 + do concurrent(k=1:l1, i=3:ido:2) ic = idp2 - i - ti5 = Cc(i,3,k) + Cc(ic,2,k) - ti2 = Cc(i,3,k) - Cc(ic,2,k) - ti4 = Cc(i,5,k) + Cc(ic,4,k) - ti3 = Cc(i,5,k) - Cc(ic,4,k) - tr5 = Cc(i-1,3,k) - Cc(ic-1,2,k) - tr2 = Cc(i-1,3,k) + Cc(ic-1,2,k) - tr4 = Cc(i-1,5,k) - Cc(ic-1,4,k) - tr3 = Cc(i-1,5,k) + Cc(ic-1,4,k) - Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3 - Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3 - cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3 - ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3 - cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3 - ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3 + ti5 = cc(i, 3, k) + cc(ic, 2, k) + ti2 = cc(i, 3, k) - cc(ic, 2, k) + ti4 = cc(i, 5, k) + cc(ic, 4, k) + ti3 = cc(i, 5, k) - cc(ic, 4, k) + tr5 = cc(i - 1, 3, k) - cc(ic - 1, 2, k) + tr2 = cc(i - 1, 3, k) + cc(ic - 1, 2, k) + tr4 = cc(i - 1, 5, k) - cc(ic - 1, 4, k) + tr3 = cc(i - 1, 5, k) + cc(ic - 1, 4, k) + ch(i - 1, k, 1) = cc(i - 1, 1, k) + tr2 + tr3 + ch(i, k, 1) = cc(i, 1, k) + ti2 + ti3 + cr2 = cc(i - 1, 1, k) + tr11*tr2 + tr12*tr3 + ci2 = cc(i, 1, k) + tr11*ti2 + tr12*ti3 + cr3 = cc(i - 1, 1, k) + tr12*tr2 + tr11*tr3 + ci3 = cc(i, 1, k) + tr12*ti2 + tr11*ti3 cr5 = ti11*tr5 + ti12*tr4 ci5 = ti11*ti5 + ti12*ti4 cr4 = ti12*tr5 - ti11*tr4 @@ -60,14 +60,13 @@ subroutine radb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) dr2 = cr2 - ci5 di5 = ci2 - cr5 di2 = ci2 + cr5 - Ch(i-1,k,2) = Wa1(i-2)*dr2 - Wa1(i-1)*di2 - Ch(i,k,2) = Wa1(i-2)*di2 + Wa1(i-1)*dr2 - Ch(i-1,k,3) = Wa2(i-2)*dr3 - Wa2(i-1)*di3 - Ch(i,k,3) = Wa2(i-2)*di3 + Wa2(i-1)*dr3 - Ch(i-1,k,4) = Wa3(i-2)*dr4 - Wa3(i-1)*di4 - Ch(i,k,4) = Wa3(i-2)*di4 + Wa3(i-1)*dr4 - Ch(i-1,k,5) = Wa4(i-2)*dr5 - Wa4(i-1)*di5 - Ch(i,k,5) = Wa4(i-2)*di5 + Wa4(i-1)*dr5 - enddo - enddo - end subroutine radb5 \ No newline at end of file + ch(i - 1, k, 2) = wa1(i - 2)*dr2 - wa1(i - 1)*di2 + ch(i, k, 2) = wa1(i - 2)*di2 + wa1(i - 1)*dr2 + ch(i - 1, k, 3) = wa2(i - 2)*dr3 - wa2(i - 1)*di3 + ch(i, k, 3) = wa2(i - 2)*di3 + wa2(i - 1)*dr3 + ch(i - 1, k, 4) = wa3(i - 2)*dr4 - wa3(i - 1)*di4 + ch(i, k, 4) = wa3(i - 2)*di4 + wa3(i - 1)*dr4 + ch(i - 1, k, 5) = wa4(i - 2)*dr5 - wa4(i - 1)*di5 + ch(i, k, 5) = wa4(i - 2)*di5 + wa4(i - 1)*dr5 + end do + end subroutine radb5 diff --git a/src/radbg.f90 b/src/radbg.f90 index 55a6af4..128d363 100644 --- a/src/radbg.f90 +++ b/src/radbg.f90 @@ -1,174 +1,147 @@ - subroutine radbg(Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) - use fftpack_kind - implicit none - real(rk) :: ai1 , ai2 , ar1 , ar1h , ar2 , ar2h , arg , c1 , & - c2 , Cc , Ch , Ch2 , dc2 , dcp , ds2 , dsp , & - Wa - integer :: i , ic , idij , Idl1 , Ido , idp2 , ik , Ip , ipp2 , & - ipph , is , j , j2 , jc , k , l , l1 , lc , nbd - dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , & - c2(Idl1,Ip) , Ch2(Idl1,Ip) , Wa(*) - real(rk),parameter :: tpi = 2*acos(-1.0_rk) ! 2 * pi - arg = tpi/real(Ip, rk) - dcp = cos(arg) - dsp = sin(arg) - idp2 = Ido + 2 - nbd = (Ido-1)/2 - ipp2 = Ip + 2 - ipph = (Ip+1)/2 - if ( Ido rk + implicit none + integer, intent(in) :: ido, ip, l1, idl1 + real(dp), intent(in) :: cc(ido, ip, l1), wa(*) + real(dp), intent(inout) :: c1(ido, l1, ip), ch2(idl1, ip) + real(dp), intent(out) :: c2(idl1, ip), ch(ido, l1, ip) + real(dp) :: ai1, ai2, ar1, ar1h, ar2, ar2h, arg, & + dc2, dcp, ds2, dsp + integer :: i, ic, idij, idp2, ik, ipp2, & + ipph, is, j, j2, jc, k, l, lc, nbd + real(dp), parameter :: tpi = 2*acos(-1.0_dp) ! 2 * pi + arg = tpi/real(ip, dp) + dcp = cos(arg) + dsp = sin(arg) + idp2 = ido + 2 + nbd = (ido - 1)/2 + ipp2 = ip + 2 + ipph = (ip + 1)/2 + if (ido < l1) then + do concurrent(i=1:ido, k=1:l1) + ch(i, k, 1) = cc(i, 1, k) + end do else - do j = 2 , ipph - jc = ipp2 - j - do k = 1 , l1 - do i = 3 , Ido , 2 - ic = idp2 - i - Ch(i-1,k,j) = Cc(i-1,2*j-1,k) + Cc(ic-1,2*j-2,k) - Ch(i-1,k,jc) = Cc(i-1,2*j-1,k) - Cc(ic-1,2*j-2,k) - Ch(i,k,j) = Cc(i,2*j-1,k) - Cc(ic,2*j-2,k) - Ch(i,k,jc) = Cc(i,2*j-1,k) + Cc(ic,2*j-2,k) - enddo - enddo - enddo - endif - endif - ar1 = 1.0_rk - ai1 = 0.0_rk - do l = 2 , ipph - lc = ipp2 - l - ar1h = dcp*ar1 - dsp*ai1 - ai1 = dcp*ai1 + dsp*ar1 - ar1 = ar1h - do ik = 1 , Idl1 - c2(ik,l) = Ch2(ik,1) + ar1*Ch2(ik,2) - c2(ik,lc) = ai1*Ch2(ik,Ip) - enddo - dc2 = ar1 - ds2 = ai1 - ar2 = ar1 - ai2 = ai1 - do j = 3 , ipph + do concurrent(k=1:l1, i=1:ido) + ch(i, k, 1) = cc(i, 1, k) + end do + end if + do concurrent(j=2:ipph, k=1:l1) jc = ipp2 - j - ar2h = dc2*ar2 - ds2*ai2 - ai2 = dc2*ai2 + ds2*ar2 - ar2 = ar2h - do ik = 1 , Idl1 - c2(ik,l) = c2(ik,l) + ar2*Ch2(ik,j) - c2(ik,lc) = c2(ik,lc) + ai2*Ch2(ik,jc) - enddo - enddo - enddo - do j = 2 , ipph - do ik = 1 , Idl1 - Ch2(ik,1) = Ch2(ik,1) + Ch2(ik,j) - enddo - enddo - do j = 2 , ipph - jc = ipp2 - j - do k = 1 , l1 - Ch(1,k,j) = c1(1,k,j) - c1(1,k,jc) - Ch(1,k,jc) = c1(1,k,j) + c1(1,k,jc) - enddo - enddo - if ( Ido/=1 ) then - if ( nbd l1) then + is = -ido + do j = 2, ip + is = is + ido + do k = 1, l1 + idij = is + do i = 3, ido, 2 + idij = idij + 2 + c1(i - 1, k, j) = wa(idij - 1)*ch(i - 1, k, j) - wa(idij) & + *ch(i, k, j) + c1(i, k, j) = wa(idij - 1)*ch(i, k, j) + wa(idij) & + *ch(i - 1, k, j) + end do + end do + end do else - do j = 2 , ipph - jc = ipp2 - j - do k = 1 , l1 - do i = 3 , Ido , 2 - Ch(i-1,k,j) = c1(i-1,k,j) - c1(i,k,jc) - Ch(i-1,k,jc) = c1(i-1,k,j) + c1(i,k,jc) - Ch(i,k,j) = c1(i,k,j) + c1(i-1,k,jc) - Ch(i,k,jc) = c1(i,k,j) - c1(i-1,k,jc) - enddo - enddo - enddo - endif - endif - if ( Ido==1 ) return - do ik = 1 , Idl1 - c2(ik,1) = Ch2(ik,1) - enddo - do j = 2 , Ip - do k = 1 , l1 - c1(1,k,j) = Ch(1,k,j) - enddo - enddo - if ( nbd>l1 ) then - is = -Ido - do j = 2 , Ip - is = is + Ido - do k = 1 , l1 + is = -ido + do j = 2, ip + is = is + ido idij = is - do i = 3 , Ido , 2 + do i = 3, ido, 2 idij = idij + 2 - c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij) & - *Ch(i,k,j) - c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij) & - *Ch(i-1,k,j) - enddo - enddo - enddo - else - is = -Ido - do j = 2 , Ip - is = is + Ido - idij = is - do i = 3 , Ido , 2 - idij = idij + 2 - do k = 1 , l1 - c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij) & - *Ch(i,k,j) - c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij) & - *Ch(i-1,k,j) - enddo - enddo - enddo - endif - end subroutine radbg \ No newline at end of file + do k = 1, l1 + c1(i - 1, k, j) = wa(idij - 1)*ch(i - 1, k, j) - wa(idij) & + *ch(i, k, j) + c1(i, k, j) = wa(idij - 1)*ch(i, k, j) + wa(idij) & + *ch(i - 1, k, j) + end do + end do + end do + end if + end subroutine radbg diff --git a/src/radf2.f90 b/src/radf2.f90 index 51a6fc4..d556f9b 100644 --- a/src/radf2.f90 +++ b/src/radf2.f90 @@ -1,31 +1,31 @@ - subroutine radf2(Ido,l1,Cc,Ch,Wa1) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Ch(Ido,2,l1) , Cc(Ido,l1,2) , Wa1(*) - do k = 1 , l1 - Ch(1,1,k) = Cc(1,k,1) + Cc(1,k,2) - Ch(Ido,2,k) = Cc(1,k,1) - Cc(1,k,2) - enddo - if ( Ido<2 ) return - if ( Ido/=2 ) then - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radf2(ido, l1, cc, ch, wa1) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 2), wa1(*) + real(dp), intent(out) :: ch(ido, 2, l1) + real(dp) :: ti2, tr2 + integer :: i, ic, idp2, k + do concurrent(k=1:l1) + ch(1, 1, k) = cc(1, k, 1) + cc(1, k, 2) + ch(ido, 2, k) = cc(1, k, 1) - cc(1, k, 2) + end do + if (ido < 2) return + if (ido /= 2) then + idp2 = ido + 2 + do concurrent(i=3:ido:2, k=1:l1) ic = idp2 - i - tr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) - ti2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) - Ch(i,1,k) = Cc(i,k,1) + ti2 - Ch(ic,2,k) = ti2 - Cc(i,k,1) - Ch(i-1,1,k) = Cc(i-1,k,1) + tr2 - Ch(ic-1,2,k) = Cc(i-1,k,1) - tr2 - enddo - enddo - if ( mod(Ido,2)==1 ) return - endif - do k = 1 , l1 - Ch(1,2,k) = -Cc(Ido,k,2) - Ch(Ido,1,k) = Cc(Ido,k,1) - enddo - end subroutine radf2 \ No newline at end of file + tr2 = wa1(i - 2)*cc(i - 1, k, 2) + wa1(i - 1)*cc(i, k, 2) + ti2 = wa1(i - 2)*cc(i, k, 2) - wa1(i - 1)*cc(i - 1, k, 2) + ch(i, 1, k) = cc(i, k, 1) + ti2 + ch(ic, 2, k) = ti2 - cc(i, k, 1) + ch(i - 1, 1, k) = cc(i - 1, k, 1) + tr2 + ch(ic - 1, 2, k) = cc(i - 1, k, 1) - tr2 + end do + if (mod(ido, 2) == 1) return + end if + do concurrent(k=1:l1) + ch(1, 2, k) = -cc(ido, k, 2) + ch(ido, 1, k) = cc(ido, k, 1) + end do + end subroutine radf2 diff --git a/src/radf3.f90 b/src/radf3.f90 index 472489a..6c44f29 100644 --- a/src/radf3.f90 +++ b/src/radf3.f90 @@ -1,40 +1,40 @@ - subroutine radf3(Ido,l1,Cc,Ch,Wa1,Wa2) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , cr2 , di2 , di3 , dr2 , dr3 , & - ti2 , ti3 , tr2 , tr3 , Wa1 , Wa2 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Ch(Ido,3,l1) , Cc(Ido,l1,3) , Wa1(*) , Wa2(*) - real(rk),parameter :: taur = -0.5_rk - ! note: original comment said this was -SQRT(3)/2 but value was 0.86602540378443864676d0 - real(rk),parameter :: taui = sqrt(3.0_rk) / 2.0_rk - do k = 1 , l1 - cr2 = Cc(1,k,2) + Cc(1,k,3) - Ch(1,1,k) = Cc(1,k,1) + cr2 - Ch(1,3,k) = taui*(Cc(1,k,3)-Cc(1,k,2)) - Ch(Ido,2,k) = Cc(1,k,1) + taur*cr2 - enddo - if ( Ido==1 ) return - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radf3(ido, l1, cc, ch, wa1, wa2) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 3), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, 3, l1) + real(dp) :: ci2, cr2, di2, di3, dr2, dr3, & + ti2, ti3, tr2, tr3 + integer :: i, ic, idp2, k + real(dp), parameter :: taur = -0.5_dp + ! note: original comment said this was -sqrt(3)/2 but value was 0.86602540378443864676d0 + real(dp), parameter :: taui = sqrt(3.0_dp)/2.0_dp + do concurrent(k=1:l1) + cr2 = cc(1, k, 2) + cc(1, k, 3) + ch(1, 1, k) = cc(1, k, 1) + cr2 + ch(1, 3, k) = taui*(cc(1, k, 3) - cc(1, k, 2)) + ch(ido, 2, k) = cc(1, k, 1) + taur*cr2 + end do + if (ido == 1) return + idp2 = ido + 2 + do concurrent(i=3:ido:2, k=1:l1) ic = idp2 - i - dr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) - di2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) - dr3 = Wa2(i-2)*Cc(i-1,k,3) + Wa2(i-1)*Cc(i,k,3) - di3 = Wa2(i-2)*Cc(i,k,3) - Wa2(i-1)*Cc(i-1,k,3) + dr2 = wa1(i - 2)*cc(i - 1, k, 2) + wa1(i - 1)*cc(i, k, 2) + di2 = wa1(i - 2)*cc(i, k, 2) - wa1(i - 1)*cc(i - 1, k, 2) + dr3 = wa2(i - 2)*cc(i - 1, k, 3) + wa2(i - 1)*cc(i, k, 3) + di3 = wa2(i - 2)*cc(i, k, 3) - wa2(i - 1)*cc(i - 1, k, 3) cr2 = dr2 + dr3 ci2 = di2 + di3 - Ch(i-1,1,k) = Cc(i-1,k,1) + cr2 - Ch(i,1,k) = Cc(i,k,1) + ci2 - tr2 = Cc(i-1,k,1) + taur*cr2 - ti2 = Cc(i,k,1) + taur*ci2 - tr3 = taui*(di2-di3) - ti3 = taui*(dr3-dr2) - Ch(i-1,3,k) = tr2 + tr3 - Ch(ic-1,2,k) = tr2 - tr3 - Ch(i,3,k) = ti2 + ti3 - Ch(ic,2,k) = ti3 - ti2 - enddo - enddo - end subroutine radf3 \ No newline at end of file + ch(i - 1, 1, k) = cc(i - 1, k, 1) + cr2 + ch(i, 1, k) = cc(i, k, 1) + ci2 + tr2 = cc(i - 1, k, 1) + taur*cr2 + ti2 = cc(i, k, 1) + taur*ci2 + tr3 = taui*(di2 - di3) + ti3 = taui*(dr3 - dr2) + ch(i - 1, 3, k) = tr2 + tr3 + ch(ic - 1, 2, k) = tr2 - tr3 + ch(i, 3, k) = ti2 + ti3 + ch(ic, 2, k) = ti3 - ti2 + end do + end subroutine radf3 diff --git a/src/radf4.f90 b/src/radf4.f90 index c0b62bb..e908cfd 100644 --- a/src/radf4.f90 +++ b/src/radf4.f90 @@ -1,58 +1,58 @@ - subroutine radf4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & - ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3, & - tr4 , Wa1 , Wa2 , Wa3 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Cc(Ido,l1,4) , Ch(Ido,4,l1) , Wa1(*) , Wa2(*) , Wa3(*) - real(rk),parameter :: hsqt2 = sqrt(2.0_rk) / 2.0_rk - do k = 1 , l1 - tr1 = Cc(1,k,2) + Cc(1,k,4) - tr2 = Cc(1,k,1) + Cc(1,k,3) - Ch(1,1,k) = tr1 + tr2 - Ch(Ido,4,k) = tr2 - tr1 - Ch(Ido,2,k) = Cc(1,k,1) - Cc(1,k,3) - Ch(1,3,k) = Cc(1,k,4) - Cc(1,k,2) - enddo - if ( Ido<2 ) return - if ( Ido/=2 ) then - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radf4(ido, l1, cc, ch, wa1, wa2, wa3) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 4), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, 4, l1) + real(dp) :: ci2, ci3, ci4, cr2, cr3, cr4, & + ti1, ti2, ti3, ti4, tr1, tr2, tr3, & + tr4 + integer :: i, ic, idp2, k + real(dp), parameter :: hsqt2 = sqrt(2.0_dp)/2.0_dp + do concurrent(k=1:l1) + tr1 = cc(1, k, 2) + cc(1, k, 4) + tr2 = cc(1, k, 1) + cc(1, k, 3) + ch(1, 1, k) = tr1 + tr2 + ch(ido, 4, k) = tr2 - tr1 + ch(ido, 2, k) = cc(1, k, 1) - cc(1, k, 3) + ch(1, 3, k) = cc(1, k, 4) - cc(1, k, 2) + end do + if (ido < 2) return + if (ido /= 2) then + idp2 = ido + 2 + do concurrent(i=3:ido:2, k=1:l1) ic = idp2 - i - cr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) - ci2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) - cr3 = Wa2(i-2)*Cc(i-1,k,3) + Wa2(i-1)*Cc(i,k,3) - ci3 = Wa2(i-2)*Cc(i,k,3) - Wa2(i-1)*Cc(i-1,k,3) - cr4 = Wa3(i-2)*Cc(i-1,k,4) + Wa3(i-1)*Cc(i,k,4) - ci4 = Wa3(i-2)*Cc(i,k,4) - Wa3(i-1)*Cc(i-1,k,4) + cr2 = wa1(i - 2)*cc(i - 1, k, 2) + wa1(i - 1)*cc(i, k, 2) + ci2 = wa1(i - 2)*cc(i, k, 2) - wa1(i - 1)*cc(i - 1, k, 2) + cr3 = wa2(i - 2)*cc(i - 1, k, 3) + wa2(i - 1)*cc(i, k, 3) + ci3 = wa2(i - 2)*cc(i, k, 3) - wa2(i - 1)*cc(i - 1, k, 3) + cr4 = wa3(i - 2)*cc(i - 1, k, 4) + wa3(i - 1)*cc(i, k, 4) + ci4 = wa3(i - 2)*cc(i, k, 4) - wa3(i - 1)*cc(i - 1, k, 4) tr1 = cr2 + cr4 tr4 = cr4 - cr2 ti1 = ci2 + ci4 ti4 = ci2 - ci4 - ti2 = Cc(i,k,1) + ci3 - ti3 = Cc(i,k,1) - ci3 - tr2 = Cc(i-1,k,1) + cr3 - tr3 = Cc(i-1,k,1) - cr3 - Ch(i-1,1,k) = tr1 + tr2 - Ch(ic-1,4,k) = tr2 - tr1 - Ch(i,1,k) = ti1 + ti2 - Ch(ic,4,k) = ti1 - ti2 - Ch(i-1,3,k) = ti4 + tr3 - Ch(ic-1,2,k) = tr3 - ti4 - Ch(i,3,k) = tr4 + ti3 - Ch(ic,2,k) = tr4 - ti3 - enddo - enddo - if ( mod(Ido,2)==1 ) return - endif - do k = 1 , l1 - ti1 = -hsqt2*(Cc(Ido,k,2)+Cc(Ido,k,4)) - tr1 = hsqt2*(Cc(Ido,k,2)-Cc(Ido,k,4)) - Ch(Ido,1,k) = tr1 + Cc(Ido,k,1) - Ch(Ido,3,k) = Cc(Ido,k,1) - tr1 - Ch(1,2,k) = ti1 - Cc(Ido,k,3) - Ch(1,4,k) = ti1 + Cc(Ido,k,3) - enddo - end subroutine radf4 \ No newline at end of file + ti2 = cc(i, k, 1) + ci3 + ti3 = cc(i, k, 1) - ci3 + tr2 = cc(i - 1, k, 1) + cr3 + tr3 = cc(i - 1, k, 1) - cr3 + ch(i - 1, 1, k) = tr1 + tr2 + ch(ic - 1, 4, k) = tr2 - tr1 + ch(i, 1, k) = ti1 + ti2 + ch(ic, 4, k) = ti1 - ti2 + ch(i - 1, 3, k) = ti4 + tr3 + ch(ic - 1, 2, k) = tr3 - ti4 + ch(i, 3, k) = tr4 + ti3 + ch(ic, 2, k) = tr4 - ti3 + end do + if (mod(ido, 2) == 1) return + end if + do concurrent(k=1:l1) + ti1 = -hsqt2*(cc(ido, k, 2) + cc(ido, k, 4)) + tr1 = hsqt2*(cc(ido, k, 2) - cc(ido, k, 4)) + ch(ido, 1, k) = tr1 + cc(ido, k, 1) + ch(ido, 3, k) = cc(ido, k, 1) - tr1 + ch(1, 2, k) = ti1 - cc(ido, k, 3) + ch(1, 4, k) = ti1 + cc(ido, k, 3) + end do + end subroutine radf4 diff --git a/src/radf5.f90 b/src/radf5.f90 index 5fbcb11..399ebf7 100644 --- a/src/radf5.f90 +++ b/src/radf5.f90 @@ -1,43 +1,43 @@ - subroutine radf5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) - use fftpack_kind - implicit none - real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & - cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & - dr4 , dr5 - real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & - tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 - integer :: i , ic , Ido , idp2 , k , l1 - dimension Cc(Ido,l1,5) , Ch(Ido,5,l1) , Wa1(*) , Wa2(*) , Wa3(*), & - Wa4(*) - real(rk),parameter :: pi = acos(-1.0_rk) - real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk) - real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) - real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk) - do k = 1 , l1 - cr2 = Cc(1,k,5) + Cc(1,k,2) - ci5 = Cc(1,k,5) - Cc(1,k,2) - cr3 = Cc(1,k,4) + Cc(1,k,3) - ci4 = Cc(1,k,4) - Cc(1,k,3) - Ch(1,1,k) = Cc(1,k,1) + cr2 + cr3 - Ch(Ido,2,k) = Cc(1,k,1) + tr11*cr2 + tr12*cr3 - Ch(1,3,k) = ti11*ci5 + ti12*ci4 - Ch(Ido,4,k) = Cc(1,k,1) + tr12*cr2 + tr11*cr3 - Ch(1,5,k) = ti12*ci5 - ti11*ci4 - enddo - if ( Ido==1 ) return - idp2 = Ido + 2 - do k = 1 , l1 - do i = 3 , Ido , 2 + subroutine radf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 5), wa1(*), wa2(*), wa3(*), wa4(*) + real(dp), intent(out) :: ch(ido, 5, l1) + real(dp) :: ci2, ci3, ci4, ci5, cr2, cr3, & + cr4, cr5, di2, di3, di4, di5, dr2, dr3, & + dr4, dr5 + real(dp) :: ti2, ti3, ti4, ti5, tr2, tr3, & + tr4, tr5 + integer :: i, ic, idp2, k + real(dp), parameter :: pi = acos(-1.0_dp) + real(dp), parameter :: tr11 = cos(2.0_dp*pi/5.0_dp) + real(dp), parameter :: ti11 = sin(2.0_dp*pi/5.0_dp) + real(dp), parameter :: tr12 = cos(4.0_dp*pi/5.0_dp) + real(dp), parameter :: ti12 = sin(4.0_dp*pi/5.0_dp) + do concurrent(k=1:l1) + cr2 = cc(1, k, 5) + cc(1, k, 2) + ci5 = cc(1, k, 5) - cc(1, k, 2) + cr3 = cc(1, k, 4) + cc(1, k, 3) + ci4 = cc(1, k, 4) - cc(1, k, 3) + ch(1, 1, k) = cc(1, k, 1) + cr2 + cr3 + ch(ido, 2, k) = cc(1, k, 1) + tr11*cr2 + tr12*cr3 + ch(1, 3, k) = ti11*ci5 + ti12*ci4 + ch(ido, 4, k) = cc(1, k, 1) + tr12*cr2 + tr11*cr3 + ch(1, 5, k) = ti12*ci5 - ti11*ci4 + end do + if (ido == 1) return + idp2 = ido + 2 + do concurrent(i=3:ido:2, k=1:l1) ic = idp2 - i - dr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) - di2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) - dr3 = Wa2(i-2)*Cc(i-1,k,3) + Wa2(i-1)*Cc(i,k,3) - di3 = Wa2(i-2)*Cc(i,k,3) - Wa2(i-1)*Cc(i-1,k,3) - dr4 = Wa3(i-2)*Cc(i-1,k,4) + Wa3(i-1)*Cc(i,k,4) - di4 = Wa3(i-2)*Cc(i,k,4) - Wa3(i-1)*Cc(i-1,k,4) - dr5 = Wa4(i-2)*Cc(i-1,k,5) + Wa4(i-1)*Cc(i,k,5) - di5 = Wa4(i-2)*Cc(i,k,5) - Wa4(i-1)*Cc(i-1,k,5) + dr2 = wa1(i - 2)*cc(i - 1, k, 2) + wa1(i - 1)*cc(i, k, 2) + di2 = wa1(i - 2)*cc(i, k, 2) - wa1(i - 1)*cc(i - 1, k, 2) + dr3 = wa2(i - 2)*cc(i - 1, k, 3) + wa2(i - 1)*cc(i, k, 3) + di3 = wa2(i - 2)*cc(i, k, 3) - wa2(i - 1)*cc(i - 1, k, 3) + dr4 = wa3(i - 2)*cc(i - 1, k, 4) + wa3(i - 1)*cc(i, k, 4) + di4 = wa3(i - 2)*cc(i, k, 4) - wa3(i - 1)*cc(i - 1, k, 4) + dr5 = wa4(i - 2)*cc(i - 1, k, 5) + wa4(i - 1)*cc(i, k, 5) + di5 = wa4(i - 2)*cc(i, k, 5) - wa4(i - 1)*cc(i - 1, k, 5) cr2 = dr2 + dr5 ci5 = dr5 - dr2 cr5 = di2 - di5 @@ -46,24 +46,23 @@ subroutine radf5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) ci4 = dr4 - dr3 cr4 = di3 - di4 ci3 = di3 + di4 - Ch(i-1,1,k) = Cc(i-1,k,1) + cr2 + cr3 - Ch(i,1,k) = Cc(i,k,1) + ci2 + ci3 - tr2 = Cc(i-1,k,1) + tr11*cr2 + tr12*cr3 - ti2 = Cc(i,k,1) + tr11*ci2 + tr12*ci3 - tr3 = Cc(i-1,k,1) + tr12*cr2 + tr11*cr3 - ti3 = Cc(i,k,1) + tr12*ci2 + tr11*ci3 + ch(i - 1, 1, k) = cc(i - 1, k, 1) + cr2 + cr3 + ch(i, 1, k) = cc(i, k, 1) + ci2 + ci3 + tr2 = cc(i - 1, k, 1) + tr11*cr2 + tr12*cr3 + ti2 = cc(i, k, 1) + tr11*ci2 + tr12*ci3 + tr3 = cc(i - 1, k, 1) + tr12*cr2 + tr11*cr3 + ti3 = cc(i, k, 1) + tr12*ci2 + tr11*ci3 tr5 = ti11*cr5 + ti12*cr4 ti5 = ti11*ci5 + ti12*ci4 tr4 = ti12*cr5 - ti11*cr4 ti4 = ti12*ci5 - ti11*ci4 - Ch(i-1,3,k) = tr2 + tr5 - Ch(ic-1,2,k) = tr2 - tr5 - Ch(i,3,k) = ti2 + ti5 - Ch(ic,2,k) = ti5 - ti2 - Ch(i-1,5,k) = tr3 + tr4 - Ch(ic-1,4,k) = tr3 - tr4 - Ch(i,5,k) = ti3 + ti4 - Ch(ic,4,k) = ti4 - ti3 - enddo - enddo - end subroutine radf5 \ No newline at end of file + ch(i - 1, 3, k) = tr2 + tr5 + ch(ic - 1, 2, k) = tr2 - tr5 + ch(i, 3, k) = ti2 + ti5 + ch(ic, 2, k) = ti5 - ti2 + ch(i - 1, 5, k) = tr3 + tr4 + ch(ic - 1, 4, k) = tr3 - tr4 + ch(i, 5, k) = ti3 + ti4 + ch(ic, 4, k) = ti4 - ti3 + end do + end subroutine radf5 diff --git a/src/radfg.f90 b/src/radfg.f90 index 589a466..3fbf3c6 100644 --- a/src/radfg.f90 +++ b/src/radfg.f90 @@ -1,180 +1,156 @@ - subroutine radfg(Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) - use fftpack_kind - implicit none - real(rk) :: ai1 , ai2 , ar1 , ar1h , ar2 , ar2h , arg , c1 , & - c2 , Cc , Ch , Ch2 , dc2 , dcp , ds2 , dsp , & - Wa - integer :: i , ic , idij , Idl1 , Ido , idp2 , ik , Ip , ipp2 , & - ipph , is , j , j2 , jc , k , l , l1 , lc , nbd - dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , & - c2(Idl1,Ip) , Ch2(Idl1,Ip) , Wa(*) - real(rk),parameter :: tpi = 2.0_rk * acos(-1.0_rk) ! 2 * pi - arg = tpi/real(Ip, rk) - dcp = cos(arg) - dsp = sin(arg) - ipph = (Ip+1)/2 - ipp2 = Ip + 2 - idp2 = Ido + 2 - nbd = (Ido-1)/2 - if ( Ido==1 ) then - do ik = 1 , Idl1 - c2(ik,1) = Ch2(ik,1) - enddo - else - do ik = 1 , Idl1 - Ch2(ik,1) = c2(ik,1) - enddo - do j = 2 , Ip - do k = 1 , l1 - Ch(1,k,j) = c1(1,k,j) - enddo - enddo - if ( nbd>l1 ) then - is = -Ido - do j = 2 , Ip - is = is + Ido - do k = 1 , l1 + subroutine radfg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: ido, ip, l1, idl1 + real(dp), intent(in) :: wa(*) + real(dp), intent(inout) :: cc(ido, ip, l1) + real(dp), intent(inout) :: c1(ido, l1, ip) + real(dp), intent(inout) :: c2(idl1, ip) + real(dp), intent(out) :: ch(ido, l1, ip) + real(dp), intent(inout) :: ch2(idl1, ip) + real(dp) :: ai1, ai2, ar1, ar1h, ar2, ar2h, arg, & + dc2, dcp, ds2, dsp + integer :: i, ic, idij, idp2, ik, ipp2, & + ipph, is, j, j2, jc, k, l, lc, nbd + real(dp), parameter :: tpi = 2.0_dp*acos(-1.0_dp) ! 2 * pi + arg = tpi/real(ip, kind=dp) + dcp = cos(arg) + dsp = sin(arg) + ipph = (ip + 1)/2 + ipp2 = ip + 2 + idp2 = ido + 2 + nbd = (ido - 1)/2 + if (ido == 1) then + do concurrent(ik=1:idl1) + c2(ik, 1) = ch2(ik, 1) + end do + else + do concurrent(ik=1:idl1) + ch2(ik, 1) = c2(ik, 1) + end do + do concurrent(j=2:ip, k=1:l1) + ch(1, k, j) = c1(1, k, j) + end do + if (nbd > l1) then + is = -ido + do j = 2, ip + is = is + ido + do k = 1, l1 + idij = is + do i = 3, ido, 2 + idij = idij + 2 + ch(i - 1, k, j) = wa(idij - 1)*c1(i - 1, k, j) + wa(idij) & + *c1(i, k, j) + ch(i, k, j) = wa(idij - 1)*c1(i, k, j) - wa(idij) & + *c1(i - 1, k, j) + end do + end do + end do + else + is = -ido + do j = 2, ip + is = is + ido idij = is - do i = 3 , Ido , 2 + do i = 3, ido, 2 idij = idij + 2 - Ch(i-1,k,j) = Wa(idij-1)*c1(i-1,k,j) + Wa(idij) & - *c1(i,k,j) - Ch(i,k,j) = Wa(idij-1)*c1(i,k,j) - Wa(idij) & - *c1(i-1,k,j) - enddo - enddo - enddo - else - is = -Ido - do j = 2 , Ip - is = is + Ido - idij = is - do i = 3 , Ido , 2 - idij = idij + 2 - do k = 1 , l1 - Ch(i-1,k,j) = Wa(idij-1)*c1(i-1,k,j) + Wa(idij) & - *c1(i,k,j) - Ch(i,k,j) = Wa(idij-1)*c1(i,k,j) - Wa(idij) & - *c1(i-1,k,j) - enddo - enddo - enddo - endif - if ( nbd rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: c(*) + real(dp), intent(in) :: wa(*) + real(dp), intent(inout) :: ch(*) + integer, intent(in) :: ifac(*) + integer :: i, idl1, ido, ip, iw, ix2, ix3, ix4, k1, & + l1, l2, na, nf + nf = ifac(2) + na = 0 + l1 = 1 + iw = 1 + do k1 = 1, nf + ip = ifac(k1 + 2) + l2 = ip*l1 + ido = n/l2 + idl1 = ido*l1 + if (ip == 4) then + ix2 = iw + ido + ix3 = ix2 + ido + if (na /= 0) then + call radb4(ido, l1, ch, c, wa(iw), wa(ix2), wa(ix3)) + else + call radb4(ido, l1, c, ch, wa(iw), wa(ix2), wa(ix3)) + end if + na = 1 - na + elseif (ip == 2) then + if (na /= 0) then + call radb2(ido, l1, ch, c, wa(iw)) + else + call radb2(ido, l1, c, ch, wa(iw)) + end if + na = 1 - na + elseif (ip == 3) then + ix2 = iw + ido + if (na /= 0) then + call radb3(ido, l1, ch, c, wa(iw), wa(ix2)) + else + call radb3(ido, l1, c, ch, wa(iw), wa(ix2)) + end if + na = 1 - na + elseif (ip /= 5) then + if (na /= 0) then + call radbg(ido, ip, l1, idl1, ch, ch, ch, c, c, wa(iw)) + else + call radbg(ido, ip, l1, idl1, c, c, c, ch, ch, wa(iw)) + end if + if (ido == 1) na = 1 - na else - call radb4(ido,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3)) - endif - na = 1 - na - elseif ( ip==2 ) then - if ( na/=0 ) then - call radb2(ido,l1,Ch,c,Wa(iw)) - else - call radb2(ido,l1,c,Ch,Wa(iw)) - endif - na = 1 - na - elseif ( ip==3 ) then - ix2 = iw + ido - if ( na/=0 ) then - call radb3(ido,l1,Ch,c,Wa(iw),Wa(ix2)) - else - call radb3(ido,l1,c,Ch,Wa(iw),Wa(ix2)) - endif - na = 1 - na - elseif ( ip/=5 ) then - if ( na/=0 ) then - call radbg(ido,ip,l1,idl1,Ch,Ch,Ch,c,c,Wa(iw)) - else - call radbg(ido,ip,l1,idl1,c,c,c,Ch,Ch,Wa(iw)) - endif - if ( ido==1 ) na = 1 - na - else - ix2 = iw + ido - ix3 = ix2 + ido - ix4 = ix3 + ido - if ( na/=0 ) then - call radb5(ido,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) - else - call radb5(ido,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) - endif - na = 1 - na - endif - l1 = l2 - iw = iw + (ip-1)*ido - enddo - if ( na==0 ) return - do i = 1 , n - c(i) = Ch(i) - enddo - end subroutine rfftb1 \ No newline at end of file + ix2 = iw + ido + ix3 = ix2 + ido + ix4 = ix3 + ido + if (na /= 0) then + call radb5(ido, l1, ch, c, wa(iw), wa(ix2), wa(ix3), wa(ix4)) + else + call radb5(ido, l1, c, ch, wa(iw), wa(ix2), wa(ix3), wa(ix4)) + end if + na = 1 - na + end if + l1 = l2 + iw = iw + (ip - 1)*ido + end do + if (na == 0) return + do concurrent(i=1:n) + c(i) = ch(i) + end do + end subroutine rfftb1 diff --git a/src/rfftf1.f90 b/src/rfftf1.f90 index acbde50..107e559 100644 --- a/src/rfftf1.f90 +++ b/src/rfftf1.f90 @@ -1,66 +1,69 @@ - subroutine rfftf1(n,c,Ch,Wa,Ifac) - use fftpack_kind - implicit none - real(rk) :: c , Ch , Wa - integer :: i , idl1 , ido , Ifac , ip , iw , ix2 , ix3 , ix4 , k1 , & - kh , l1 , l2 , n , na , nf - dimension Ch(*) , c(*) , Wa(*) , Ifac(*) - nf = Ifac(2) - na = 1 - l2 = n - iw = n - do k1 = 1 , nf - kh = nf - k1 - ip = Ifac(kh+3) - l1 = l2/ip - ido = n/l2 - idl1 = ido*l1 - iw = iw - (ip-1)*ido - na = 1 - na - if ( ip==4 ) then - ix2 = iw + ido - ix3 = ix2 + ido - if ( na/=0 ) then - call radf4(ido,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3)) - else - call radf4(ido,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3)) - endif - elseif ( ip/=2 ) then - if ( ip==3 ) then + subroutine rfftf1(n, c, ch, wa, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: c(*) + real(dp), intent(in) :: wa(*) + real(dp), intent(inout) :: ch(*) + integer, intent(in) :: ifac(*) + integer :: i, idl1, ido, ip, iw, ix2, ix3, ix4, k1, & + kh, l1, l2, na, nf + nf = ifac(2) + na = 1 + l2 = n + iw = n + do k1 = 1, nf + kh = nf - k1 + ip = ifac(kh + 3) + l1 = l2/ip + ido = n/l2 + idl1 = ido*l1 + iw = iw - (ip - 1)*ido + na = 1 - na + if (ip == 4) then ix2 = iw + ido - if ( na/=0 ) then - call radf3(ido,l1,Ch,c,Wa(iw),Wa(ix2)) + ix3 = ix2 + ido + if (na /= 0) then + call radf4(ido, l1, ch, c, wa(iw), wa(ix2), wa(ix3)) else - call radf3(ido,l1,c,Ch,Wa(iw),Wa(ix2)) - endif - elseif ( ip/=5 ) then - if ( ido==1 ) na = 1 - na - if ( na/=0 ) then - call radfg(ido,ip,l1,idl1,Ch,Ch,Ch,c,c,Wa(iw)) - na = 0 + call radf4(ido, l1, c, ch, wa(iw), wa(ix2), wa(ix3)) + end if + elseif (ip /= 2) then + if (ip == 3) then + ix2 = iw + ido + if (na /= 0) then + call radf3(ido, l1, ch, c, wa(iw), wa(ix2)) + else + call radf3(ido, l1, c, ch, wa(iw), wa(ix2)) + end if + elseif (ip /= 5) then + if (ido == 1) na = 1 - na + if (na /= 0) then + call radfg(ido, ip, l1, idl1, ch, ch, ch, c, c, wa(iw)) + na = 0 + else + call radfg(ido, ip, l1, idl1, c, c, c, ch, ch, wa(iw)) + na = 1 + end if else - call radfg(ido,ip,l1,idl1,c,c,c,Ch,Ch,Wa(iw)) - na = 1 - endif + ix2 = iw + ido + ix3 = ix2 + ido + ix4 = ix3 + ido + if (na /= 0) then + call radf5(ido, l1, ch, c, wa(iw), wa(ix2), wa(ix3), wa(ix4)) + else + call radf5(ido, l1, c, ch, wa(iw), wa(ix2), wa(ix3), wa(ix4)) + end if + end if + elseif (na /= 0) then + call radf2(ido, l1, ch, c, wa(iw)) else - ix2 = iw + ido - ix3 = ix2 + ido - ix4 = ix3 + ido - if ( na/=0 ) then - call radf5(ido,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) - else - call radf5(ido,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) - endif - endif - elseif ( na/=0 ) then - call radf2(ido,l1,Ch,c,Wa(iw)) - else - call radf2(ido,l1,c,Ch,Wa(iw)) - endif - l2 = l1 - enddo - if ( na==1 ) return - do i = 1 , n - c(i) = Ch(i) - enddo - end subroutine rfftf1 \ No newline at end of file + call radf2(ido, l1, c, ch, wa(iw)) + end if + l2 = l1 + end do + if (na == 1) return + do concurrent(i=1:n) + c(i) = ch(i) + end do + end subroutine rfftf1 diff --git a/src/rffti1.f90 b/src/rffti1.f90 index 167eb5d..5ae8f41 100644 --- a/src/rffti1.f90 +++ b/src/rffti1.f90 @@ -1,64 +1,66 @@ - subroutine rffti1(n,Wa,Ifac) - use fftpack_kind - implicit none - real(rk) :: arg , argh , argld , fi , Wa - integer :: i , ib , ido , Ifac , ii , ip , ipm , is , j , k1 , l1 , & - l2 , ld , n , nf , nfm1 , nl , nq , nr , ntry - dimension Wa(*) , Ifac(*) - integer,dimension(4),parameter :: ntryh = [4 , 2 , 3 , 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) - is = 0 - nfm1 = nf - 1 - l1 = 1 - if ( nfm1==0 ) return - do k1 = 1 , nfm1 - ip = Ifac(k1+2) - ld = 0 - l2 = l1*ip - ido = n/l2 - ipm = ip - 1 - do j = 1 , ipm - ld = ld + l1 - i = is - argld = real(ld, rk)*argh - fi = 0.0_rk - do ii = 3 , ido , 2 - i = i + 2 - fi = fi + 1.0_rk - arg = fi*argld - Wa(i-1) = cos(arg) - Wa(i) = sin(arg) - enddo - is = is + ido - enddo - l1 = l2 - enddo - end subroutine rffti1 \ No newline at end of file + subroutine rffti1(n, wa, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wa(*) + integer, intent(out) :: ifac(*) + real(dp) :: arg, argh, argld, fi + integer :: i, ib, ido, ii, ip, ipm, is, j, k1, l1, & + l2, ld, nf, nfm1, nl, nq, nr, ntry + integer, dimension(4), parameter :: ntryh = [4, 2, 3, 5] + real(dp), parameter :: tpi = 2.0_dp*acos(-1.0_dp) ! 2 * pi + nl = n + nf = 0 + j = 0 +100 j = j + 1 + if (j <= 4) then + ntry = ntryh(j) + else + ntry = ntry + 2 + end if +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) + end do + ifac(3) = 2 + end if + end if + if (nl /= 1) goto 200 + ifac(1) = n + ifac(2) = nf + argh = tpi/real(n, dp) + is = 0 + nfm1 = nf - 1 + l1 = 1 + if (nfm1 == 0) return + do k1 = 1, nfm1 + ip = ifac(k1 + 2) + ld = 0 + l2 = l1*ip + ido = n/l2 + ipm = ip - 1 + do j = 1, ipm + ld = ld + l1 + i = is + argld = real(ld, dp)*argh + fi = 0.0_dp + do ii = 3, ido, 2 + i = i + 2 + fi = fi + 1.0_dp + arg = fi*argld + wa(i - 1) = cos(arg) + wa(i) = sin(arg) + end do + is = is + ido + end do + l1 = l2 + end do + end subroutine rffti1 diff --git a/src/rk.f90 b/src/rk.f90 index 663df9c..3ab114d 100644 --- a/src/rk.f90 +++ b/src/rk.f90 @@ -1,4 +1,4 @@ - module fftpack_kind - implicit none - integer,parameter :: rk = kind(1.0d0) - end module fftpack_kind +module fftpack_kind + use, intrinsic :: iso_fortran_env, only: rk => real64 + implicit none(type, external) +end module fftpack_kind diff --git a/src/sint1.f90 b/src/sint1.f90 index 09cad5d..dddecd1 100644 --- a/src/sint1.f90 +++ b/src/sint1.f90 @@ -1,43 +1,46 @@ - subroutine sint1(n,War,Was,Xh,x,Ifac) - use fftpack_kind - implicit none - integer :: i , Ifac , k , kc , modn , n , np1 , ns2 - real(rk) :: t1 , t2 , War , Was , x , Xh , xhold - dimension War(*) , Was(*) , x(*) , Xh(*) , Ifac(*) - real(rk),parameter :: sqrt3 = sqrt(3.0_rk) - do i = 1 , n - Xh(i) = War(i) - War(i) = x(i) - enddo - if ( n<2 ) then - Xh(1) = Xh(1) + Xh(1) - elseif ( n==2 ) then - xhold = sqrt3*(Xh(1)+Xh(2)) - Xh(2) = sqrt3*(Xh(1)-Xh(2)) - Xh(1) = xhold - else - np1 = n + 1 - ns2 = n/2 - x(1) = 0.0_rk - do k = 1 , ns2 - kc = np1 - k - t1 = Xh(k) - Xh(kc) - t2 = Was(k)*(Xh(k)+Xh(kc)) - x(k+1) = t1 + t2 - x(kc+1) = t2 - t1 - enddo - modn = mod(n,2) - if ( modn/=0 ) x(ns2+2) = 4.0_rk*Xh(ns2+1) - call rfftf1(np1,x,Xh,War,Ifac) - Xh(1) = 0.5_rk*x(1) - do i = 3 , n , 2 - Xh(i-1) = -x(i) - Xh(i) = Xh(i-2) + x(i-1) - enddo - if ( modn==0 ) Xh(n) = -x(n+1) - endif - do i = 1 , n - x(i) = War(i) - War(i) = Xh(i) - enddo - end subroutine sint1 \ No newline at end of file + subroutine sint1(n, war, was, xh, x, ifac) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n, ifac(*) + real(dp), intent(in) :: was(*) + real(dp), intent(inout) :: war(*), x(*) + real(dp), intent(out) :: xh(*) + integer :: i, k, kc, modn, np1, ns2 + real(dp) :: t1, t2, xhold + real(dp), parameter :: sqrt3 = sqrt(3.0_dp) + do i = 1, n + xh(i) = war(i) + war(i) = x(i) + end do + if (n < 2) then + xh(1) = xh(1) + xh(1) + elseif (n == 2) then + xhold = sqrt3*(xh(1) + xh(2)) + xh(2) = sqrt3*(xh(1) - xh(2)) + xh(1) = xhold + else + np1 = n + 1 + ns2 = n/2 + x(1) = 0.0_dp + do k = 1, ns2 + kc = np1 - k + t1 = xh(k) - xh(kc) + t2 = was(k)*(xh(k) + xh(kc)) + x(k + 1) = t1 + t2 + x(kc + 1) = t2 - t1 + end do + modn = mod(n, 2) + if (modn /= 0) x(ns2 + 2) = 4.0_dp*xh(ns2 + 1) + call rfftf1(np1, x, xh, war, ifac) + xh(1) = 0.5_dp*x(1) + do i = 3, n, 2 + xh(i - 1) = -x(i) + xh(i) = xh(i - 2) + x(i - 1) + end do + if (modn == 0) xh(n) = -x(n + 1) + end if + do i = 1, n + x(i) = war(i) + war(i) = xh(i) + end do + end subroutine sint1 diff --git a/src/zfftb.f90 b/src/zfftb.f90 index 792526e..aae8cd5 100644 --- a/src/zfftb.f90 +++ b/src/zfftb.f90 @@ -1,11 +1,12 @@ - subroutine zfftb(n,c,Wsave) - use fftpack_kind - implicit none - real(rk) :: c , Wsave - integer :: iw1 , iw2 , n - dimension c(*) , Wsave(*) - if ( n==1 ) return - iw1 = n + n + 1 - iw2 = iw1 + n + n - call cfftb1(n,c,Wsave,Wsave(iw1),Wsave(iw2)) - end subroutine zfftb \ No newline at end of file + subroutine zfftb(n, c, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: c(*) + real(dp), intent(inout) :: wsave(*) + integer :: iw1, iw2 + if (n == 1) return + iw1 = n + n + 1 + iw2 = iw1 + n + n + call cfftb1(n, c, wsave, wsave(iw1), wsave(iw2)) + end subroutine zfftb diff --git a/src/zfftf.f90 b/src/zfftf.f90 index 1403702..675b57e 100644 --- a/src/zfftf.f90 +++ b/src/zfftf.f90 @@ -1,11 +1,12 @@ - subroutine zfftf(n,c,Wsave) - use fftpack_kind - implicit none - real(rk) :: c , Wsave - integer :: iw1 , iw2 , n - dimension c(*) , Wsave(*) - if ( n==1 ) return - iw1 = n + n + 1 - iw2 = iw1 + n + n - call cfftf1(n,c,Wsave,Wsave(iw1),Wsave(iw2)) - end subroutine zfftf \ No newline at end of file + subroutine zfftf(n, c, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: c(*) + real(dp), intent(inout) :: wsave(*) + integer :: iw1, iw2 + if (n == 1) return + iw1 = n + n + 1 + iw2 = iw1 + n + n + call cfftf1(n, c, wsave, wsave(iw1), wsave(iw2)) + end subroutine zfftf diff --git a/src/zffti.f90 b/src/zffti.f90 index 5cc7cca..d6da6d8 100644 --- a/src/zffti.f90 +++ b/src/zffti.f90 @@ -1,11 +1,11 @@ - subroutine zffti(n,Wsave) - use fftpack_kind - implicit none - integer :: iw1 , iw2 , n - real(rk) :: Wsave - dimension Wsave(*) - if ( n==1 ) return - iw1 = n + n + 1 - iw2 = iw1 + n + n - call cffti1(n,Wsave(iw1),Wsave(iw2)) - end subroutine zffti \ No newline at end of file + subroutine zffti(n, wsave) + use fftpack_kind, only: dp => rk + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + integer :: iw1, iw2 + if (n == 1) return + iw1 = n + n + 1 + iw2 = iw1 + n + n + call cffti1(n, wsave(iw1), wsave(iw2)) + end subroutine zffti