diff --git a/.github/workflows/fpm-deployment.yml b/.github/workflows/fpm-deployment.yml index 2e7dc50d..d8368080 100644 --- a/.github/workflows/fpm-deployment.yml +++ b/.github/workflows/fpm-deployment.yml @@ -56,6 +56,11 @@ jobs: - run: | # Use fpm gnu ci to check xdp and qp python config/fypp_deployment.py --with_xdp --with_qp fpm test --profile release --flag '-DWITH_XDP -DWITH_QP -coverage' + fpm run --example example_bench01_zfft --profile release --flag '-DWITH_XDP -DWITH_QP -coverage' + fpm run --example example_bench02_zfft --profile release --flag '-DWITH_XDP -DWITH_QP -coverage' + fpm run --example example_bench03_dfft --profile release --flag '-DWITH_XDP -DWITH_QP -coverage' + fpm run --example example_real_transforms --profile release --flag '-DWITH_XDP -DWITH_QP -coverage' + fpm run --example example_complex_transforms --profile release --flag '-DWITH_XDP -DWITH_QP -coverage' - name: Create coverage report run: | diff --git a/example/example_bench01_zfft.f90 b/example/example_bench01_zfft.f90 index 87ca85e5..1b0793fb 100644 --- a/example/example_bench01_zfft.f90 +++ b/example/example_bench01_zfft.f90 @@ -1,43 +1,43 @@ program bench1 -use fftpack, only: zffti, zfftf, zfftb -use fftpack_kind, only: rk -implicit none -complex(rk), allocatable :: z(:) -real(rk), allocatable :: w(:), x(:) -real(rk) :: err, time_init, time_forward, time_backward, t1, t2 -integer :: N + use fftpack, only: zffti, zfftf, zfftb + use fftpack_kinds, only: dp + implicit none + complex(dp), allocatable :: z(:) + real(dp), allocatable :: w(:), x(:) + real(dp) :: err, time_init, time_forward, time_backward, t1, t2 + integer :: N -N = 1024*1014*16 + N = 1024*1014*16 -allocate(x(N), z(N), w(4*N+15)) -call random_number(x) -z = x + allocate (x(N), z(N), w(4*N + 15)) + call random_number(x) + z = x -print *, "01: Benchmarking zfft" -print *, "Initializing" -call cpu_time(t1) -call zffti(N, w) -call cpu_time(t2) -time_init = t2-t1 + print *, "01: Benchmarking zfft" + print *, "Initializing" + call cpu_time(t1) + call zffti(N, w) + call cpu_time(t2) + time_init = t2 - t1 -print *, "Forward" -call cpu_time(t1) -call zfftf(N, z, w) -call cpu_time(t2) -time_forward = t2-t1 + print *, "Forward" + call cpu_time(t1) + call zfftf(N, z, w) + call cpu_time(t2) + time_forward = t2 - t1 -print *, "Backward" -call cpu_time(t1) -call zfftb(N, z, w) -call cpu_time(t2) -time_backward = t2-t1 -print *, "Done" + print *, "Backward" + call cpu_time(t1) + call zfftb(N, z, w) + call cpu_time(t2) + time_backward = t2 - t1 + print *, "Done" -err = maxval(abs(x-real(z/N,rk))) -print * -print *, "Error: ", err -print *, "Init time: ", time_init -print *, "Forward time: ", time_forward -print *, "Backward time: ", time_backward -print *, "" + err = maxval(abs(x - real(z/N, dp))) + print * + print *, "Error: ", err + print *, "Init time: ", time_init + print *, "Forward time: ", time_forward + print *, "Backward time: ", time_backward + print *, "" end program diff --git a/example/example_bench02_zfft.f90 b/example/example_bench02_zfft.f90 index cb435201..35322fb6 100644 --- a/example/example_bench02_zfft.f90 +++ b/example/example_bench02_zfft.f90 @@ -1,64 +1,64 @@ program bench2 -use fftpack, only: zffti, zfftf, zfftb, fft, ifft -use fftpack_kind, only: rk + use fftpack, only: zffti, zfftf, zfftb, fft, ifft + use fftpack_kinds, only: dp -implicit none -integer, parameter :: N = 1024*135*77 ! (2**10)*(3**3)*5*7*11 + implicit none + integer, parameter :: N = 1024*135*77 ! (2**10)*(3**3)*5*7*11 -complex(rk), dimension(N) :: x, z -real(rk), dimension(4*N+15) :: w -real(rk) :: err, time_i, time_f, time_b, t1, t2 + complex(dp), dimension(N) :: x, z + real(dp), dimension(4*N + 15) :: w + real(dp) :: err, time_i, time_f, time_b, t1, t2 -call random_number(x%re) -z = x + call random_number(x%re) + z = x -print *, "02: Benchmarking zfft & fft" + print *, "02: Benchmarking zfft & fft" -call cpu_time(t1) -call zffti(N, w) -call cpu_time(t2) -time_i = t2-t1 -print *, "Initializing: done" + call cpu_time(t1) + call zffti(N, w) + call cpu_time(t2) + time_i = t2 - t1 + print *, "Initializing: done" -call cpu_time(t1) -call zfftf(N, z, w) -call cpu_time(t2) -time_f = t2-t1 -print *, "Forward: done" + call cpu_time(t1) + call zfftf(N, z, w) + call cpu_time(t2) + time_f = t2 - t1 + print *, "Forward: done" -call cpu_time(t1) -call zfftb(N, z, w) -call cpu_time(t2) -time_b = t2-t1 -print *, "Backward: done" -print *, "" + call cpu_time(t1) + call zfftb(N, z, w) + call cpu_time(t2) + time_b = t2 - t1 + print *, "Backward: done" + print *, "" -err = maxval(abs(x-real(z/N,rk))) -print *, "--zfft" -print *, "Error: ", err -print *, "Init. time: ", time_i -print *, "Forward time: ", time_f -print *, "Backward time: ", time_b -print *, "" + err = maxval(abs(x - real(z/N, dp))) + print *, "--zfft" + print *, "Error: ", err + print *, "Init. time: ", time_i + print *, "Forward time: ", time_f + print *, "Backward time: ", time_b + print *, "" -print *, "Comparing to calls through fft" -call cpu_time(t1) -z = fft(x) -call cpu_time(t2) -time_f = t2-t1 -print *, "Init. & forward: done" + print *, "Comparing to calls through fft" + call cpu_time(t1) + z = fft(x) + call cpu_time(t2) + time_f = t2 - t1 + print *, "Init. & forward: done" -call cpu_time(t1) -z = ifft(z) -call cpu_time(t2) -time_b = t2-t1 -print *, "Backward: done" -print *, "" + call cpu_time(t1) + z = ifft(z) + call cpu_time(t2) + time_b = t2 - t1 + print *, "Backward: done" + print *, "" -err = maxval(abs(x-real(z/N,rk))) -print *, "--fft" -print *, "Error: ", err -print *, "Init. & forward time: ", time_f -print *, "Backward time: ", time_b -print *, "" + err = maxval(abs(x - real(z/N, dp))) + print *, "--fft" + print *, "Error: ", err + print *, "Init. & forward time: ", time_f + print *, "Backward time: ", time_b + print *, "" end program diff --git a/example/example_bench03_dfft.f90 b/example/example_bench03_dfft.f90 index e55e417e..61d53aef 100644 --- a/example/example_bench03_dfft.f90 +++ b/example/example_bench03_dfft.f90 @@ -1,64 +1,64 @@ program bench3 -use fftpack, only: dffti, dfftf, dfftb, rfft, irfft -use fftpack_kind, only: rk + use fftpack, only: dffti, dfftf, dfftb, rfft, irfft + use fftpack_kinds, only: dp -implicit none -integer, parameter :: N = 1024*135*77 ! (2**10)*(3**3)*5*7*11 + implicit none + integer, parameter :: N = 1024*135*77 ! (2**10)*(3**3)*5*7*11 -real(rk), dimension(N) :: x, y -real(rk), dimension(2*N+15) :: w -real(rk) :: err, time_i, time_f, time_b, t1, t2 + real(dp), dimension(N) :: x, y + real(dp), dimension(2*N + 15) :: w + real(dp) :: err, time_i, time_f, time_b, t1, t2 -call random_number(x) -y = x + call random_number(x) + y = x -print *, "03: Benchmarking dfft & rfft" + print *, "03: Benchmarking dfft & rfft" -call cpu_time(t1) -call dffti(N, w) -call cpu_time(t2) -time_i = t2-t1 -print *, "Initializing: done" + call cpu_time(t1) + call dffti(N, w) + call cpu_time(t2) + time_i = t2 - t1 + print *, "Initializing: done" -call cpu_time(t1) -call dfftf(N, y, w) -call cpu_time(t2) -time_f = t2-t1 -print *, "Forward: done" + call cpu_time(t1) + call dfftf(N, y, w) + call cpu_time(t2) + time_f = t2 - t1 + print *, "Forward: done" -call cpu_time(t1) -call dfftb(N, y, w) -call cpu_time(t2) -time_b = t2-t1 -print *, "Backward: done" -print *, "" + call cpu_time(t1) + call dfftb(N, y, w) + call cpu_time(t2) + time_b = t2 - t1 + print *, "Backward: done" + print *, "" -err = maxval(abs(x-y/N)) -print *, "--dfft" -print *, "Error: ", err -print *, "Init. time: ", time_i -print *, "Forward time: ", time_f -print *, "Backward time: ", time_b -print *, "" + err = maxval(abs(x - y/N)) + print *, "--dfft" + print *, "Error: ", err + print *, "Init. time: ", time_i + print *, "Forward time: ", time_f + print *, "Backward time: ", time_b + print *, "" -print *, "Comparing to calls through rfft" -call cpu_time(t1) -y = rfft(x) -call cpu_time(t2) -time_f = t2-t1 -print *, "Init. & forward: done" + print *, "Comparing to calls through rfft" + call cpu_time(t1) + y = rfft(x) + call cpu_time(t2) + time_f = t2 - t1 + print *, "Init. & forward: done" -call cpu_time(t1) -y = irfft(y) -call cpu_time(t2) -time_b = t2-t1 -print *, "Backward: done" -print *, "" + call cpu_time(t1) + y = irfft(y) + call cpu_time(t2) + time_b = t2 - t1 + print *, "Backward: done" + print *, "" -err = maxval(abs(x-y/N)) -print *, "--rfft" -print *, "Error: ", err -print *, "Init. & forward time: ", time_f -print *, "Backward time: ", time_b -print *, "" + err = maxval(abs(x - y/N)) + print *, "--rfft" + print *, "Error: ", err + print *, "Init. & forward time: ", time_f + print *, "Backward time: ", time_b + print *, "" end program diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cd88c3c6..ee57240d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,106 +1,21 @@ #### Pre-process: .fpp -> .f90 via Fypp # Create a list of the files to be preprocessed -# set(fppFiles -# stdlib_ascii.fypp -# stdlib_bitsets.fypp -# stdlib_bitsets_64.fypp -# stdlib_bitsets_large.fypp -# stdlib_codata_type.fypp -# stdlib_constants.fypp -# stdlib_error.fypp -# stdlib_hash_32bit.fypp -# stdlib_hash_32bit_fnv.fypp -# stdlib_hash_32bit_nm.fypp -# stdlib_hash_32bit_water.fypp -# stdlib_hash_64bit.fypp -# stdlib_hash_64bit_fnv.fypp -# stdlib_hash_64bit_pengy.fypp -# stdlib_hash_64bit_spookyv2.fypp -# stdlib_intrinsics_dot_product.fypp -# stdlib_intrinsics_sum.fypp -# stdlib_intrinsics.fypp -# stdlib_io.fypp -# stdlib_io_npy.fypp -# stdlib_io_npy_load.fypp -# stdlib_io_npy_save.fypp -# stdlib_kinds.fypp -# stdlib_linalg.fypp -# stdlib_linalg_diag.fypp -# stdlib_linalg_least_squares.fypp -# stdlib_linalg_outer_product.fypp -# stdlib_linalg_kronecker.fypp -# stdlib_linalg_cross_product.fypp -# stdlib_linalg_eigenvalues.fypp -# stdlib_linalg_solve.fypp -# stdlib_linalg_determinant.fypp -# stdlib_linalg_qr.fypp -# stdlib_linalg_inverse.fypp -# stdlib_linalg_pinv.fypp -# stdlib_linalg_norms.fypp -# stdlib_linalg_state.fypp -# stdlib_linalg_svd.fypp -# stdlib_linalg_cholesky.fypp -# stdlib_linalg_schur.fypp -# stdlib_optval.fypp -# stdlib_selection.fypp -# stdlib_sorting.fypp -# stdlib_sorting_ord_sort.fypp -# stdlib_sorting_sort.fypp -# stdlib_sorting_sort_index.fypp -# stdlib_sparse_constants.fypp -# stdlib_sparse_conversion.fypp -# stdlib_sparse_kinds.fypp -# stdlib_sparse_spmv.fypp -# stdlib_specialfunctions_activations.fypp -# stdlib_specialfunctions_gamma.fypp -# stdlib_specialfunctions.fypp -# stdlib_specialmatrices.fypp -# stdlib_specialmatrices_tridiagonal.fypp -# stdlib_stats.fypp -# stdlib_stats_corr.fypp -# stdlib_stats_cov.fypp -# stdlib_stats_mean.fypp -# stdlib_stats_median.fypp -# stdlib_stats_moment.fypp -# stdlib_stats_moment_all.fypp -# stdlib_stats_moment_mask.fypp -# stdlib_stats_moment_scalar.fypp -# stdlib_stats_distribution_uniform.fypp -# stdlib_stats_distribution_normal.fypp -# stdlib_stats_distribution_exponential.fypp -# stdlib_stats_var.fypp -# stdlib_quadrature.fypp -# stdlib_quadrature_trapz.fypp -# stdlib_quadrature_simps.fypp -# stdlib_random.fypp -# stdlib_math.fypp -# stdlib_math_linspace.fypp -# stdlib_math_logspace.fypp -# stdlib_math_arange.fypp -# stdlib_math_is_close.fypp -# stdlib_math_all_close.fypp -# stdlib_math_diff.fypp -# stdlib_math_meshgrid.fypp -# stdlib_str2num.fypp -# stdlib_string_type.fypp -# stdlib_string_type_constructor.fypp -# stdlib_strings_to_string.fypp -# stdlib_strings.fypp -# stdlib_version.fypp -# ) +set(fppFiles + fftpack_kinds.fypp +) # Preprocessed files to contain preprocessor directives -> .F90 # set(cppFiles # stdlib_linalg_constants.fypp -# stdlib_linalg_blas.fypp -# stdlib_linalg_lapack.fypp # ) -# fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) # fypp_f90pp("${fyppFlags}" "${cppFiles}" outPreprocFiles) set(SRC + fftpack_legacy_drivers_pass.f90 + fftpack_legacy_drivers_rad.f90 cfftb1.f90 cfftf1.f90 cffti1.f90 @@ -132,30 +47,9 @@ set(SRC fftpack_irfft.f90 fftpack_rfft.f90 fftpack_utils.f90 - passb2.f90 - passb3.f90 - passb4.f90 - passb5.f90 - passb.f90 - passf2.f90 - passf3.f90 - passf4.f90 - passf5.f90 - passf.f90 - radb2.f90 - radb3.f90 - radb4.f90 - radb5.f90 - radbg.f90 - radf2.f90 - radf3.f90 - radf4.f90 - radf5.f90 - radfg.f90 rfftb1.f90 rfftf1.f90 rffti1.f90 - rk.f90 sint1.f90 zfftb.f90 zfftf.f90 @@ -166,12 +60,6 @@ set(SRC add_library(${PROJECT_NAME} ${SRC}) -# # Link to BLAS and LAPACK -# if(BLAS_FOUND AND LAPACK_FOUND) -# target_link_libraries(${PROJECT_NAME} "BLAS::BLAS") -# target_link_libraries(${PROJECT_NAME} "LAPACK::LAPACK") -# endif() - set_target_properties( ${PROJECT_NAME} PROPERTIES diff --git a/src/cfftb1.f90 b/src/cfftb1.f90 index 7c0bf79d..21d3d6ff 100644 --- a/src/cfftb1.f90 +++ b/src/cfftb1.f90 @@ -1,5 +1,6 @@ subroutine cfftb1(n, c, ch, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp + use fftpack_legacy_drivers_pass, only: passb, passb2, passb3, passb4, passb5 implicit none integer, intent(in) :: n, ifac(*) real(dp), intent(in) :: wa(*) diff --git a/src/cfftf1.f90 b/src/cfftf1.f90 index f692e3af..81574138 100644 --- a/src/cfftf1.f90 +++ b/src/cfftf1.f90 @@ -1,5 +1,6 @@ subroutine cfftf1(n, c, ch, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp + use fftpack_legacy_drivers_pass, only: passf, passf2, passf3, passf4, passf5 implicit none integer, intent(in) :: n, ifac(*) real(dp), intent(inout) :: c(*), ch(*) diff --git a/src/cffti1.f90 b/src/cffti1.f90 index 6647e9e2..5486d40a 100644 --- a/src/cffti1.f90 +++ b/src/cffti1.f90 @@ -1,5 +1,5 @@ subroutine cffti1(n, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n integer, intent(out) :: ifac(*) diff --git a/src/cosqb1.f90 b/src/cosqb1.f90 index 3d2b699d..db093158 100644 --- a/src/cosqb1.f90 +++ b/src/cosqb1.f90 @@ -1,5 +1,5 @@ subroutine cosqb1(n, x, w, xh) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: x(*) diff --git a/src/cosqf1.f90 b/src/cosqf1.f90 index 21ff74ed..f5faacce 100644 --- a/src/cosqf1.f90 +++ b/src/cosqf1.f90 @@ -1,5 +1,5 @@ subroutine cosqf1(n, x, w, xh) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: x(*) diff --git a/src/dcosqb.f90 b/src/dcosqb.f90 index 9b3c34cc..ddc5b4fd 100644 --- a/src/dcosqb.f90 +++ b/src/dcosqb.f90 @@ -1,5 +1,5 @@ subroutine dcosqb(n, x, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(in) :: wsave(*) diff --git a/src/dcosqf.f90 b/src/dcosqf.f90 index e8233876..3b761e50 100644 --- a/src/dcosqf.f90 +++ b/src/dcosqf.f90 @@ -1,5 +1,5 @@ subroutine dcosqf(n, x, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(in) :: wsave(*) diff --git a/src/dcosqi.f90 b/src/dcosqi.f90 index 184b2145..76be4d14 100644 --- a/src/dcosqi.f90 +++ b/src/dcosqi.f90 @@ -1,5 +1,5 @@ subroutine dcosqi(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/src/dcost.f90 b/src/dcost.f90 index 4273fbb8..199ecb68 100644 --- a/src/dcost.f90 +++ b/src/dcost.f90 @@ -1,10 +1,10 @@ subroutine dcost(n, x, wsave) - use fftpack_kind, only: rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n - real(rk), intent(inout) :: wsave(*) - real(rk), intent(inout) :: x(*) - real(rk) :: c1, t1, t2, tx2, x1h, x1p3, & + real(dp), intent(inout) :: wsave(*) + real(dp), intent(inout) :: x(*) + real(dp) :: c1, t1, t2, tx2, x1h, x1p3, & xi, xim2 integer :: i, k, kc, modn, nm1, np1, ns2 nm1 = n - 1 diff --git a/src/dcosti.f90 b/src/dcosti.f90 index ad68e859..3a52913e 100644 --- a/src/dcosti.f90 +++ b/src/dcosti.f90 @@ -1,5 +1,5 @@ subroutine dcosti(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/src/dfftb.f90 b/src/dfftb.f90 index 8620e63f..b43bb1cd 100644 --- a/src/dfftb.f90 +++ b/src/dfftb.f90 @@ -1,5 +1,5 @@ subroutine dfftb(n, r, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: r(*) diff --git a/src/dfftf.f90 b/src/dfftf.f90 index 39064ef8..501ca3af 100644 --- a/src/dfftf.f90 +++ b/src/dfftf.f90 @@ -1,5 +1,5 @@ subroutine dfftf(n, r, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: r(*) diff --git a/src/dffti.f90 b/src/dffti.f90 index ec7a341e..553e4010 100644 --- a/src/dffti.f90 +++ b/src/dffti.f90 @@ -1,5 +1,5 @@ subroutine dffti(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/src/dsinqb.f90 b/src/dsinqb.f90 index dfc85049..89a36f52 100644 --- a/src/dsinqb.f90 +++ b/src/dsinqb.f90 @@ -1,5 +1,5 @@ subroutine dsinqb(n, x, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: x(*) diff --git a/src/dsinqf.f90 b/src/dsinqf.f90 index 388bb17f..23cf18a3 100644 --- a/src/dsinqf.f90 +++ b/src/dsinqf.f90 @@ -1,5 +1,5 @@ subroutine dsinqf(n, x, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: x(*) diff --git a/src/dsinqi.f90 b/src/dsinqi.f90 index 642b7954..a9fd00d3 100644 --- a/src/dsinqi.f90 +++ b/src/dsinqi.f90 @@ -1,5 +1,5 @@ subroutine dsinqi(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/src/dsint.f90 b/src/dsint.f90 index 3515772d..ee529374 100644 --- a/src/dsint.f90 +++ b/src/dsint.f90 @@ -1,5 +1,5 @@ subroutine dsint(n, x, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: x(*) diff --git a/src/dsinti.f90 b/src/dsinti.f90 index a2b2cfbb..ea425c36 100644 --- a/src/dsinti.f90 +++ b/src/dsinti.f90 @@ -1,5 +1,5 @@ subroutine dsinti(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/src/dzfftb.f90 b/src/dzfftb.f90 index 1acb8c33..5fa42331 100644 --- a/src/dzfftb.f90 +++ b/src/dzfftb.f90 @@ -1,5 +1,5 @@ subroutine dzfftb(n, r, azero, a, b, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: r(*) diff --git a/src/dzfftf.f90 b/src/dzfftf.f90 index 3ce5d8ab..96f5e37c 100644 --- a/src/dzfftf.f90 +++ b/src/dzfftf.f90 @@ -1,6 +1,5 @@ subroutine dzfftf(n, r, azero, a, b, wsave) -! version 3 june 1979 - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(in) :: r(*) diff --git a/src/dzffti.f90 b/src/dzffti.f90 index 46876397..5f5bbef1 100644 --- a/src/dzffti.f90 +++ b/src/dzffti.f90 @@ -1,5 +1,5 @@ subroutine dzffti(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/src/ezfft1.f90 b/src/ezfft1.f90 index 38a89152..cf7507da 100644 --- a/src/ezfft1.f90 +++ b/src/ezfft1.f90 @@ -1,5 +1,5 @@ subroutine ezfft1(n, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wa(*) diff --git a/src/fftpack.f90 b/src/fftpack.f90 index ed14799b..a338861d 100644 --- a/src/fftpack.f90 +++ b/src/fftpack.f90 @@ -1,362 +1,387 @@ module fftpack - use fftpack_kind - - implicit none - private - - public :: zffti, zfftf, zfftb - public :: fft, ifft - public :: fftshift, ifftshift - public :: fftfreq, rfftfreq - - public :: dffti, dfftf, dfftb - public :: rfft, irfft - - public :: dzffti, dzfftf, dzfftb - - public :: dcosqi, dcosqf, dcosqb - public :: dcosti, dcost - public :: dct, idct - public :: dct_t1i, dct_t1 - public :: dct_t23i, dct_t2, dct_t3 - - public :: dsinti, dsint - - public :: rk - - interface - - !> Version: experimental - !> - !> Initialize `zfftf` and `zfftb`. - !> ([Specification](../page/specs/fftpack.html#zffti)) - pure subroutine zffti(n, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: wsave(*) - end subroutine zffti - - !> Version: experimental - !> - !> Forward transform of a complex periodic sequence. - !> ([Specification](../page/specs/fftpack.html#zfftf)) - pure subroutine zfftf(n, c, wsave) - import rk - integer, intent(in) :: n - complex(kind=rk), intent(inout) :: c(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine zfftf - - !> Version: experimental - !> - !> Unnormalized inverse of `zfftf`. - !> ([Specification](../page/specs/fftpack.html#zfftb)) - pure subroutine zfftb(n, c, wsave) - import rk - integer, intent(in) :: n - complex(kind=rk), intent(inout) :: c(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine zfftb - - !> Version: experimental - !> - !> Initialize `dfftf` and `dfftb`. - !> ([Specification](../page/specs/fftpack.html#dffti)) - pure subroutine dffti(n, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: wsave(*) - end subroutine dffti - - !> Version: experimental - !> - !> Forward transform of a real periodic sequence. - !> ([Specification](../page/specs/fftpack.html#dfftf)) - pure subroutine dfftf(n, r, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(inout) :: r(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dfftf - - !> Version: experimental - !> - !> Unnormalized inverse of `dfftf`. - !> ([Specification](../page/specs/fftpack.html#dfftb)) - pure subroutine dfftb(n, r, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(inout) :: r(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dfftb - - !> Version: experimental - !> - !> Initialize `dzfftf` and `dzfftb`. - !> ([Specification](../page/specs/fftpack.html#dzffti)) - pure subroutine dzffti(n, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: wsave(*) - end subroutine dzffti - - !> Version: experimental - !> - !> Simplified forward transform of a real periodic sequence. - !> ([Specification](../page/specs/fftpack.html#dzfftf)) - pure subroutine dzfftf(n, r, azero, a, b, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(in) :: r(*) - real(kind=rk), intent(out) :: azero - real(kind=rk), intent(out) :: a(*), b(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dzfftf - - !> Version: experimental - !> - !> Unnormalized inverse of `dzfftf`. - !> ([Specification](../page/specs/fftpack.html#dzfftb)) - pure subroutine dzfftb(n, r, azero, a, b, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: r(*) - real(kind=rk), intent(in) :: azero - real(kind=rk), intent(in) :: a(*), b(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dzfftb - - !> Version: experimental - !> - !> Initialize `dcosqf` and `dcosqb`. - !> ([Specification](../page/specs/fftpack.html#initialize-dct-2-3-dcosqi-or-dct_t23i)) - pure subroutine dcosqi(n, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: wsave(*) - end subroutine dcosqi - - !> Version: experimental - !> - !> Forward transform of quarter wave data. - !> ([Specification](../page/specs/fftpack.html#compute-dct-3-dcosqf-or-dct_t3)) - pure subroutine dcosqf(n, x, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(inout) :: x(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dcosqf - - !> Version: experimental - !> - !> Unnormalized inverse of `dcosqf`. - !> ([Specification](../page/specs/fftpack.html#compute-dct-2-dcosqb-or-dct_t2)) - pure subroutine dcosqb(n, x, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(inout) :: x(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dcosqb - - !> Version: experimental - !> - !> Initialize `dcost`. - !> ([Specification](../page/specs/fftpack.html#initialize-dct-1-dcosti-or-dct_t1i)) - pure subroutine dcosti(n, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: wsave(*) - end subroutine dcosti - - !> Version: experimental - !> - !> Discrete fourier cosine transform of an even sequence. - !> ([Specification](../page/specs/fftpack.html#compute-dct-1-dcost-or-dct_t1)) - pure subroutine dcost(n, x, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(inout) :: x(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dcost - - - pure subroutine dsinti(n, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(out) :: wsave(*) - end subroutine dsinti - - - pure subroutine dsint(n, x, wsave) - import rk - integer, intent(in) :: n - real(kind=rk), intent(inout) :: x(*) - real(kind=rk), intent(in) :: wsave(*) - end subroutine dsint - - - !> Version: experimental - !> - !> Integer frequency values involved in complex FFT. - !> ([Specifiction](../page/specs/fftpack.html#fftfreq)) - pure module function fftfreq(n) result(out) - integer, intent(in) :: n - integer, dimension(n) :: out - end function fftfreq - - !> Version: experimental - !> - !> Integer frequency values involved in real FFT. - !> ([Specifiction](../page/specs/fftpack.html#rfftfreq)) - pure module function rfftfreq(n) result(out) - integer, intent(in) :: n - integer, dimension(n) :: out - end function rfftfreq - - end interface - - !> Version: experimental - !> - !> Forward transform of a complex periodic sequence. - !> ([Specifiction](../page/specs/fftpack.html#fft)) - interface fft - pure module function fft_rk(x, n) result(result) - complex(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - complex(kind=rk), allocatable :: result(:) - end function fft_rk - end interface fft - - !> Version: experimental - !> - !> Backward transform of a complex periodic sequence. - !> ([Specifiction](../page/specs/fftpack.html#ifft)) - interface ifft - pure module function ifft_rk(x, n) result(result) - complex(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - complex(kind=rk), allocatable :: result(:) - end function ifft_rk - end interface ifft - - !> Version: experimental - !> - !> Forward transform of a real periodic sequence. - !> ([Specifiction](../page/specs/fftpack.html#rfft)) - interface rfft - pure module function rfft_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - end function rfft_rk - end interface rfft - - !> Version: experimental - !> - !> Backward transform of a real periodic sequence. - !> ([Specifiction](../page/specs/fftpack.html#irfft)) - interface irfft - pure module function irfft_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - end function irfft_rk - end interface irfft - - !> Version: experimental - !> - !> Dsicrete cosine transforms. - !> ([Specification](../page/specs/fftpack.html#simplified-dct-of-types-1-2-3-dct)) - interface dct - pure module function dct_rk(x, n, type) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - integer, intent(in), optional :: type - real(kind=rk), allocatable :: result(:) - end function dct_rk - end interface dct - - !> Version: experimental - !> - !> Inverse discrete cosine transforms. - !> ([Specification](../page/specs/fftpack.html#simplified-inverse-dct-of-types-1-2-3-idct)) - interface idct - pure module function idct_rk(x, n, type) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - integer, intent(in), optional :: type - real(kind=rk), allocatable :: result(:) - end function idct_rk - end interface idct - - !> Version: experimental - !> - !> Initialize DCT type-1 - !> ([Specification](../page/specs/fftpack.html#initialize-dct-1-dcosti-or-dct_t1i)) - interface dct_t1i - procedure :: dcosti - end interface dct_t1i - - !> Version: experimental - !> - !> Perform DCT type-1 - !> ([Specification](../page/specs/fftpack.html#compute-dct-1-dcost-or-dct_t1)) - interface dct_t1 - procedure :: dcost - end interface dct_t1 - - !> Version: experimental - !> - !> Initialize DCT types 2, 3 - !> ([Specification](../page/specs/fftpack.html#initialize-dct-2-3-dcosqi-or-dct_t23i)) - interface dct_t23i - procedure :: dcosqi - end interface dct_t23i - - !> Version: experimental - !> - !> Perform DCT type-2 - !> ([Specification](../page/specs/fftpack.html#compute-dct-2-dcosqb-or-dct_t2)) - interface dct_t2 - procedure :: dcosqb - end interface dct_t2 - - !> Version: experimental - !> - !> Perform DCT type-3 - !> ([Specification](../page/specs/fftpack.html#compute-dct-3-dcosqf-or-dct_t3)) - interface dct_t3 - procedure :: dcosqf - end interface dct_t3 - - !> Version: experimental - !> - !> Shifts zero-frequency component to center of spectrum. - !> ([Specifiction](../page/specs/fftpack.html#fftshift)) - interface fftshift - pure module function fftshift_crk(x) result(result) - complex(kind=rk), intent(in) :: x(:) - complex(kind=rk), dimension(size(x)) :: result - end function fftshift_crk - pure module function fftshift_rrk(x) result(result) - real(kind=rk), intent(in) :: x(:) - real(kind=rk), dimension(size(x)) :: result - end function fftshift_rrk - end interface fftshift - - !> Version: experimental - !> - !> Shifts zero-frequency component to beginning of spectrum. - !> ([Specifiction](../page/specs/fftpack.html#ifftshift)) - interface ifftshift - pure module function ifftshift_crk(x) result(result) - complex(kind=rk), intent(in) :: x(:) - complex(kind=rk), dimension(size(x)) :: result - end function ifftshift_crk - pure module function ifftshift_rrk(x) result(result) - real(kind=rk), intent(in) :: x(:) - real(kind=rk), dimension(size(x)) :: result - end function ifftshift_rrk - end interface ifftshift + use fftpack_kinds, only: dp + + implicit none + private + + public :: zffti, zfftf, zfftb + public :: fft, ifft + public :: fftshift, ifftshift + public :: fftfreq, rfftfreq + + public :: dffti, dfftf, dfftb + public :: rfft, irfft + + public :: dzffti, dzfftf, dzfftb + + public :: dcosqi, dcosqf, dcosqb + public :: dcosti, dcost + public :: dct, idct + public :: dct_t1i, dct_t1 + public :: dct_t23i, dct_t2, dct_t3 + + public :: dsinti, dsint + + public :: dp + + interface + + !> Version: experimental + !> + !> Initialize `zfftf` and `zfftb`. + !> ([Specification](../page/specs/fftpack.html#zffti)) + pure subroutine zffti(n, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + end subroutine zffti + + !> Version: experimental + !> + !> Forward transform of a complex periodic sequence. + !> ([Specification](../page/specs/fftpack.html#zfftf)) + pure subroutine zfftf(n, c, wsave) + import dp + implicit none + integer, intent(in) :: n + complex(dp), intent(inout) :: c(*) + real(dp), intent(in) :: wsave(*) + end subroutine zfftf + + !> Version: experimental + !> + !> Unnormalized inverse of `zfftf`. + !> ([Specification](../page/specs/fftpack.html#zfftb)) + pure subroutine zfftb(n, c, wsave) + import dp + implicit none + integer, intent(in) :: n + complex(dp), intent(inout) :: c(*) + real(dp), intent(in) :: wsave(*) + end subroutine zfftb + + !> Version: experimental + !> + !> Initialize `dfftf` and `dfftb`. + !> ([Specification](../page/specs/fftpack.html#dffti)) + pure subroutine dffti(n, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + end subroutine dffti + + !> Version: experimental + !> + !> Forward transform of a real periodic sequence. + !> ([Specification](../page/specs/fftpack.html#dfftf)) + pure subroutine dfftf(n, r, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: r(*) + real(dp), intent(in) :: wsave(*) + end subroutine dfftf + + !> Version: experimental + !> + !> Unnormalized inverse of `dfftf`. + !> ([Specification](../page/specs/fftpack.html#dfftb)) + pure subroutine dfftb(n, r, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: r(*) + real(dp), intent(in) :: wsave(*) + end subroutine dfftb + + !> Version: experimental + !> + !> Initialize `dzfftf` and `dzfftb`. + !> ([Specification](../page/specs/fftpack.html#dzffti)) + pure subroutine dzffti(n, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + end subroutine dzffti + + !> Version: experimental + !> + !> Simplified forward transform of a real periodic sequence. + !> ([Specification](../page/specs/fftpack.html#dzfftf)) + pure subroutine dzfftf(n, r, azero, a, b, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(in) :: r(*) + real(dp), intent(out) :: azero + real(dp), intent(out) :: a(*), b(*) + real(dp), intent(in) :: wsave(*) + end subroutine dzfftf + + !> Version: experimental + !> + !> Unnormalized inverse of `dzfftf`. + !> ([Specification](../page/specs/fftpack.html#dzfftb)) + pure subroutine dzfftb(n, r, azero, a, b, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: r(*) + real(dp), intent(in) :: azero + real(dp), intent(in) :: a(*), b(*) + real(dp), intent(in) :: wsave(*) + end subroutine dzfftb + + !> Version: experimental + !> + !> Initialize `dcosqf` and `dcosqb`. + !> ([Specification](../page/specs/fftpack.html#initialize-dct-2-3-dcosqi-or-dct_t23i)) + pure subroutine dcosqi(n, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + end subroutine dcosqi + + !> Version: experimental + !> + !> Forward transform of quarter wave data. + !> ([Specification](../page/specs/fftpack.html#compute-dct-3-dcosqf-or-dct_t3)) + pure subroutine dcosqf(n, x, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + end subroutine dcosqf + + !> Version: experimental + !> + !> Unnormalized inverse of `dcosqf`. + !> ([Specification](../page/specs/fftpack.html#compute-dct-2-dcosqb-or-dct_t2)) + pure subroutine dcosqb(n, x, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + end subroutine dcosqb + + !> Version: experimental + !> + !> Initialize `dcost`. + !> ([Specification](../page/specs/fftpack.html#initialize-dct-1-dcosti-or-dct_t1i)) + pure subroutine dcosti(n, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + end subroutine dcosti + + !> Version: experimental + !> + !> Discrete fourier cosine transform of an even sequence. + !> ([Specification](../page/specs/fftpack.html#compute-dct-1-dcost-or-dct_t1)) + pure subroutine dcost(n, x, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + end subroutine dcost + + pure subroutine dsinti(n, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(out) :: wsave(*) + end subroutine dsinti + + pure subroutine dsint(n, x, wsave) + import dp + implicit none + integer, intent(in) :: n + real(dp), intent(inout) :: x(*) + real(dp), intent(in) :: wsave(*) + end subroutine dsint + + !> Version: experimental + !> + !> Integer frequency values involved in complex FFT. + !> ([Specifiction](../page/specs/fftpack.html#fftfreq)) + pure module function fftfreq(n) result(out) + implicit none + integer, intent(in) :: n + integer, dimension(n) :: out + end function fftfreq + + !> Version: experimental + !> + !> Integer frequency values involved in real FFT. + !> ([Specifiction](../page/specs/fftpack.html#rfftfreq)) + pure module function rfftfreq(n) result(out) + implicit none + integer, intent(in) :: n + integer, dimension(n) :: out + end function rfftfreq + + end interface + + !> Version: experimental + !> + !> Forward transform of a complex periodic sequence. + !> ([Specifiction](../page/specs/fftpack.html#fft)) + interface fft + pure module function fft_dp(x, n) result(result) + implicit none + complex(dp), intent(in) :: x(:) + integer, intent(in), optional :: n + complex(dp), allocatable :: result(:) + end function fft_dp + end interface fft + + !> Version: experimental + !> + !> Backward transform of a complex periodic sequence. + !> ([Specifiction](../page/specs/fftpack.html#ifft)) + interface ifft + pure module function ifft_dp(x, n) result(result) + implicit none + complex(dp), intent(in) :: x(:) + integer, intent(in), optional :: n + complex(dp), allocatable :: result(:) + end function ifft_dp + end interface ifft + + !> Version: experimental + !> + !> Forward transform of a real periodic sequence. + !> ([Specifiction](../page/specs/fftpack.html#rfft)) + interface rfft + pure module function rfft_dp(x, n) result(result) + implicit none + real(dp), intent(in) :: x(:) + integer, intent(in), optional :: n + real(dp), allocatable :: result(:) + end function rfft_dp + end interface rfft + + !> Version: experimental + !> + !> Backward transform of a real periodic sequence. + !> ([Specifiction](../page/specs/fftpack.html#irfft)) + interface irfft + pure module function irfft_dp(x, n) result(result) + implicit none + real(dp), intent(in) :: x(:) + integer, intent(in), optional :: n + real(dp), allocatable :: result(:) + end function irfft_dp + end interface irfft + + !> Version: experimental + !> + !> Dsicrete cosine transforms. + !> ([Specification](../page/specs/fftpack.html#simplified-dct-of-types-1-2-3-dct)) + interface dct + pure module function dct_dp(x, n, type) result(result) + implicit none + real(dp), intent(in) :: x(:) + integer, intent(in), optional :: n + integer, intent(in), optional :: type + real(dp), allocatable :: result(:) + end function dct_dp + end interface dct + + !> Version: experimental + !> + !> Inverse discrete cosine transforms. + !> ([Specification](../page/specs/fftpack.html#simplified-inverse-dct-of-types-1-2-3-idct)) + interface idct + pure module function idct_dp(x, n, type) result(result) + implicit none + real(dp), intent(in) :: x(:) + integer, intent(in), optional :: n + integer, intent(in), optional :: type + real(dp), allocatable :: result(:) + end function idct_dp + end interface idct + + !> Version: experimental + !> + !> Initialize DCT type-1 + !> ([Specification](../page/specs/fftpack.html#initialize-dct-1-dcosti-or-dct_t1i)) + interface dct_t1i + procedure :: dcosti + end interface dct_t1i + + !> Version: experimental + !> + !> Perform DCT type-1 + !> ([Specification](../page/specs/fftpack.html#compute-dct-1-dcost-or-dct_t1)) + interface dct_t1 + procedure :: dcost + end interface dct_t1 + + !> Version: experimental + !> + !> Initialize DCT types 2, 3 + !> ([Specification](../page/specs/fftpack.html#initialize-dct-2-3-dcosqi-or-dct_t23i)) + interface dct_t23i + procedure :: dcosqi + end interface dct_t23i + + !> Version: experimental + !> + !> Perform DCT type-2 + !> ([Specification](../page/specs/fftpack.html#compute-dct-2-dcosqb-or-dct_t2)) + interface dct_t2 + procedure :: dcosqb + end interface dct_t2 + + !> Version: experimental + !> + !> Perform DCT type-3 + !> ([Specification](../page/specs/fftpack.html#compute-dct-3-dcosqf-or-dct_t3)) + interface dct_t3 + procedure :: dcosqf + end interface dct_t3 + + !> Version: experimental + !> + !> Shifts zero-frequency component to center of spectrum. + !> ([Specifiction](../page/specs/fftpack.html#fftshift)) + interface fftshift + pure module function fftshift_cdp(x) result(result) + implicit none + complex(dp), intent(in) :: x(:) + complex(dp), dimension(size(x)) :: result + end function fftshift_cdp + pure module function fftshift_rdp(x) result(result) + implicit none + real(dp), intent(in) :: x(:) + real(dp), dimension(size(x)) :: result + end function fftshift_rdp + end interface fftshift + + !> Version: experimental + !> + !> Shifts zero-frequency component to beginning of spectrum. + !> ([Specifiction](../page/specs/fftpack.html#ifftshift)) + interface ifftshift + pure module function ifftshift_cdp(x) result(result) + implicit none + complex(dp), intent(in) :: x(:) + complex(dp), dimension(size(x)) :: result + end function ifftshift_cdp + pure module function ifftshift_rdp(x) result(result) + implicit none + real(dp), intent(in) :: x(:) + real(dp), dimension(size(x)) :: result + end function ifftshift_rdp + end interface ifftshift end module fftpack diff --git a/src/fftpack_dct.f90 b/src/fftpack_dct.f90 index e7652157..ad53d69b 100644 --- a/src/fftpack_dct.f90 +++ b/src/fftpack_dct.f90 @@ -1,103 +1,103 @@ submodule(fftpack) fftpack_dct - + implicit none contains - !> Discrete cosine transforms of types 1, 2, 3. - pure module function dct_rk(x, n, type) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - integer, intent(in), optional :: type - real(kind=rk), allocatable :: result(:) + !> Discrete cosine transforms of types 1, 2, 3. + pure module function dct_dp(x, n, type) result(result) + real(kind=dp), intent(in) :: x(:) + integer, intent(in), optional :: n + integer, intent(in), optional :: type + real(kind=dp), allocatable :: result(:) - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) + integer :: lenseq, lensav, i + real(kind=dp), allocatable :: wsave(:) - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, (0.0_rk, i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, (0.0_dp, i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if - ! Default to DCT-2 - if (.not.present(type)) then - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - call dcosqb(lenseq, result, wsave) - return - end if + ! Default to DCT-2 + if (.not. present(type)) then + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqb(lenseq, result, wsave) + return + end if - if (type == 1) then ! DCT-1 - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosti(lenseq, wsave) - call dcost(lenseq, result, wsave) - else if (type == 2) then ! DCT-2 - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - call dcosqb(lenseq, result, wsave) - else if (type == 3) then ! DCT-3 - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - call dcosqf(lenseq, result, wsave) - end if - end function dct_rk + if (type == 1) then ! DCT-1 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosti(lenseq, wsave) + call dcost(lenseq, result, wsave) + else if (type == 2) then ! DCT-2 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqb(lenseq, result, wsave) + else if (type == 3) then ! DCT-3 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqf(lenseq, result, wsave) + end if + end function dct_dp - !> Inverse discrete cosine transforms of types 1, 2, 3. - pure module function idct_rk(x, n, type) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - integer, intent(in), optional :: type - real(kind=rk), allocatable :: result(:) + !> Inverse discrete cosine transforms of types 1, 2, 3. + pure module function idct_dp(x, n, type) result(result) + real(kind=dp), intent(in) :: x(:) + integer, intent(in), optional :: n + integer, intent(in), optional :: type + real(kind=dp), allocatable :: result(:) - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) + integer :: lenseq, lensav, i + real(kind=dp), allocatable :: wsave(:) - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, (0.0_rk, i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, (0.0_dp, i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if - ! Default to t=2; inverse DCT-2 is DCT-3 - if (.not.present(type)) then - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - call dcosqf(lenseq, result, wsave) - return - end if + ! Default to t=2; inverse DCT-2 is DCT-3 + if (.not. present(type)) then + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqf(lenseq, result, wsave) + return + end if - if (type == 1) then ! inverse DCT-1 is DCT-1 - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosti(lenseq, wsave) - call dcost(lenseq, result, wsave) - else if (type == 2) then ! inverse DCT-2 is DCT-3 - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - call dcosqf(lenseq, result, wsave) - else if (type == 3) then ! inverse DCT-3 is DCT-2 - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - call dcosqb(lenseq, result, wsave) - end if - end function idct_rk + if (type == 1) then ! inverse DCT-1 is DCT-1 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosti(lenseq, wsave) + call dcost(lenseq, result, wsave) + else if (type == 2) then ! inverse DCT-2 is DCT-3 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqf(lenseq, result, wsave) + else if (type == 3) then ! inverse DCT-3 is DCT-2 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqb(lenseq, result, wsave) + end if + end function idct_dp end submodule fftpack_dct diff --git a/src/fftpack_fft.f90 b/src/fftpack_fft.f90 index 8f02ee01..9a7d66ab 100644 --- a/src/fftpack_fft.f90 +++ b/src/fftpack_fft.f90 @@ -1,36 +1,36 @@ submodule(fftpack) fftpack_fft - + implicit none contains - !> Forward transform of a complex periodic sequence. - pure module function fft_rk(x, n) result(result) - complex(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - complex(kind=rk), allocatable :: result(:) - - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) - - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, ((0.0_rk, 0.0_rk), i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if - - !> Initialize FFT - lensav = 4*lenseq + 15 - allocate (wsave(lensav)) - call zffti(lenseq, wsave) - - !> Forward transformation - call zfftf(lenseq, result, wsave) - - end function fft_rk + !> Forward transform of a complex periodic sequence. + pure module function fft_dp(x, n) result(result) + complex(kind=dp), intent(in) :: x(:) + integer, intent(in), optional :: n + complex(kind=dp), allocatable :: result(:) + + integer :: lenseq, lensav, i + real(kind=dp), allocatable :: wsave(:) + + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, ((0.0_dp, 0.0_dp), i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if + + !> Initialize FFT + lensav = 4*lenseq + 15 + allocate (wsave(lensav)) + call zffti(lenseq, wsave) + + !> Forward transformation + call zfftf(lenseq, result, wsave) + + end function fft_dp end submodule fftpack_fft diff --git a/src/fftpack_fftshift.f90 b/src/fftpack_fftshift.f90 index 69d0f414..343ad66a 100644 --- a/src/fftpack_fftshift.f90 +++ b/src/fftpack_fftshift.f90 @@ -1,23 +1,23 @@ submodule(fftpack) fftpack_fftshift - + implicit none contains - !> Shifts zero-frequency component to center of spectrum for `complex` type. - pure module function fftshift_crk(x) result(result) - complex(kind=rk), intent(in) :: x(:) - complex(kind=rk), dimension(size(x)) :: result + !> Shifts zero-frequency component to center of spectrum for `complex` type. + pure module function fftshift_cdp(x) result(result) + complex(kind=dp), intent(in) :: x(:) + complex(kind=dp), dimension(size(x)) :: result - result = cshift(x, shift=-floor(0.5_rk*size(x))) + result = cshift(x, shift=-floor(0.5_dp*size(x))) - end function fftshift_crk + end function fftshift_cdp - !> Shifts zero-frequency component to center of spectrum for `real` type. - pure module function fftshift_rrk(x) result(result) - real(kind=rk), intent(in) :: x(:) - real(kind=rk), dimension(size(x)) :: result + !> Shifts zero-frequency component to center of spectrum for `real` type. + pure module function fftshift_rdp(x) result(result) + real(kind=dp), intent(in) :: x(:) + real(kind=dp), dimension(size(x)) :: result - result = cshift(x, shift=-floor(0.5_rk*size(x))) + result = cshift(x, shift=-floor(0.5_dp*size(x))) - end function fftshift_rrk + end function fftshift_rdp end submodule fftpack_fftshift diff --git a/src/fftpack_ifft.f90 b/src/fftpack_ifft.f90 index 680e64b4..2afee0f1 100644 --- a/src/fftpack_ifft.f90 +++ b/src/fftpack_ifft.f90 @@ -1,36 +1,36 @@ submodule(fftpack) fftpack_ifft - + implicit none contains - !> Backward transform of a complex periodic sequence. - pure module function ifft_rk(x, n) result(result) - complex(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - complex(kind=rk), allocatable :: result(:) - - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) - - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, ((0.0_rk, 0.0_rk), i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if - - !> Initialize FFT - lensav = 4*lenseq + 15 - allocate (wsave(lensav)) - call zffti(lenseq, wsave) - - !> Backward transformation - call zfftb(lenseq, result, wsave) - - end function ifft_rk + !> Backward transform of a complex periodic sequence. + pure module function ifft_dp(x, n) result(result) + complex(kind=dp), intent(in) :: x(:) + integer, intent(in), optional :: n + complex(kind=dp), allocatable :: result(:) + + integer :: lenseq, lensav, i + real(kind=dp), allocatable :: wsave(:) + + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, ((0.0_dp, 0.0_dp), i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if + + !> Initialize FFT + lensav = 4*lenseq + 15 + allocate (wsave(lensav)) + call zffti(lenseq, wsave) + + !> Backward transformation + call zfftb(lenseq, result, wsave) + + end function ifft_dp end submodule fftpack_ifft diff --git a/src/fftpack_ifftshift.f90 b/src/fftpack_ifftshift.f90 index b36d2c83..7d20011a 100644 --- a/src/fftpack_ifftshift.f90 +++ b/src/fftpack_ifftshift.f90 @@ -1,23 +1,23 @@ submodule(fftpack) fftpack_ifftshift - + implicit none contains - !> Shifts zero-frequency component to beginning of spectrum for `complex` type. - pure module function ifftshift_crk(x) result(result) - complex(kind=rk), intent(in) :: x(:) - complex(kind=rk), dimension(size(x)) :: result + !> Shifts zero-frequency component to beginning of spectrum for `complex` type. + pure module function ifftshift_cdp(x) result(result) + complex(kind=dp), intent(in) :: x(:) + complex(kind=dp), dimension(size(x)) :: result - result = cshift(x, shift=-ceiling(0.5_rk*size(x))) + result = cshift(x, shift=-ceiling(0.5_dp*size(x))) - end function ifftshift_crk + end function ifftshift_cdp - !> Shifts zero-frequency component to beginning of spectrum for `real` type. - pure module function ifftshift_rrk(x) result(result) - real(kind=rk), intent(in) :: x(:) - real(kind=rk), dimension(size(x)) :: result + !> Shifts zero-frequency component to beginning of spectrum for `real` type. + pure module function ifftshift_rdp(x) result(result) + real(kind=dp), intent(in) :: x(:) + real(kind=dp), dimension(size(x)) :: result - result = cshift(x, shift=-ceiling(0.5_rk*size(x))) + result = cshift(x, shift=-ceiling(0.5_dp*size(x))) - end function ifftshift_rrk + end function ifftshift_rdp end submodule fftpack_ifftshift diff --git a/src/fftpack_irfft.f90 b/src/fftpack_irfft.f90 index 13cdb053..dd496cfa 100644 --- a/src/fftpack_irfft.f90 +++ b/src/fftpack_irfft.f90 @@ -1,36 +1,36 @@ submodule(fftpack) fftpack_irfft - + implicit none contains - !> Backward transform of a real periodic sequence. - pure module function irfft_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) - - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, (0.0_rk, i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if - - !> Initialize FFT - lensav = 2*lenseq + 15 - allocate (wsave(lensav)) - call dffti(lenseq, wsave) - - !> Backward transformation - call dfftb(lenseq, result, wsave) - - end function irfft_rk + !> Backward transform of a real periodic sequence. + pure module function irfft_dp(x, n) result(result) + real(kind=dp), intent(in) :: x(:) + integer, intent(in), optional :: n + real(kind=dp), allocatable :: result(:) + + integer :: lenseq, lensav, i + real(kind=dp), allocatable :: wsave(:) + + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, (0.0_dp, i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if + + !> Initialize FFT + lensav = 2*lenseq + 15 + allocate (wsave(lensav)) + call dffti(lenseq, wsave) + + !> Backward transformation + call dfftb(lenseq, result, wsave) + + end function irfft_dp end submodule fftpack_irfft diff --git a/src/fftpack_kinds.fypp b/src/fftpack_kinds.fypp new file mode 100644 index 00000000..2e39035a --- /dev/null +++ b/src/fftpack_kinds.fypp @@ -0,0 +1,21 @@ +#:include "common.fypp" +module fftpack_kinds + implicit none(type, external) + private + public :: sp, dp, xdp, qp, lk + + !> Single precision real numbers. + integer, parameter :: sp = selected_real_kind(6) + + !> Double precision real numbers. + integer, parameter :: dp = selected_real_kind(8) + + !> Extended double precision real numbers + integer, parameter :: xdp = #{if WITH_XDP}#selected_real_kind(18)#{else}#-1#{endif}# + + !> Quadruple precision real numbers + integer, parameter :: qp = #{if WITH_QP}#selected_real_kind(33)#{else}#-1#{endif}# + + !> Default logical kind parameter + integer, parameter :: lk = kind(.true.) +end module fftpack_kinds diff --git a/src/fftpack_legacy_drivers_pass.f90 b/src/fftpack_legacy_drivers_pass.f90 new file mode 100644 index 00000000..d4b7a37c --- /dev/null +++ b/src/fftpack_legacy_drivers_pass.f90 @@ -0,0 +1,748 @@ +module fftpack_legacy_drivers_pass + use fftpack_kinds, only: dp + implicit none(type, external) + private + + public :: passb, passb2, passb3, passb4, passb5 + public :: passf, passf2, passf3, passf4, passf5 + + interface + pure module subroutine passb(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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) + end subroutine passb + end interface + + interface + pure module subroutine passb2(ido, l1, cc, ch, wa1) + implicit none + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 2, l1), wa1(*) + real(dp), intent(out) :: ch(ido, l1, 2) + end subroutine passb2 + end interface + + interface + pure module subroutine passb3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 3, l1), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, l1, 3) + end subroutine passb3 + end interface + + interface + pure module subroutine passb4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 4, l1), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, l1, 4) + end subroutine passb4 + end interface + + interface + pure module subroutine passb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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) + end subroutine passb5 + end interface + + interface + pure module subroutine passf(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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) + end subroutine passf + end interface + + interface + pure module subroutine passf2(ido, l1, cc, ch, wa1) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 2, l1), wa1(*) + real(dp), intent(out) :: ch(ido, l1, 2) + end subroutine passf2 + end interface + + interface + pure module subroutine passf3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 3, l1), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, l1, 3) + end subroutine passf3 + end interface + + interface + pure module subroutine passf4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 4, l1), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, l1, 4) + end subroutine passf4 + end interface + + interface + pure module subroutine passf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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) + end subroutine passf5 + end interface + +contains + + pure module subroutine passb(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + 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 < l1) then + do concurrent(k=1:l1, j=2:ipph, i=1:ido) + jc = ipp2 - j + ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) + ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) + end do + do concurrent(k=1:l1, i=1:ido) + ch(i, k, 1) = cc(i, 1, k) + end do + else + do concurrent(k=1:l1, j=2:ipph, i=1:ido) + jc = ipp2 - j + ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) + ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) + end do + do concurrent(i=1:ido, k=1:l1) + ch(i, k, 1) = cc(i, 1, k) + end do + end if + idl = 2 - ido + inc = 0 + do l = 2, ipph + lc = ipp2 - l + idl = idl + ido + do concurrent(ik=1:idl1) + c2(ik, l) = ch2(ik, 1) + wa(idl - 1)*ch2(ik, 2) + c2(ik, lc) = wa(idl)*ch2(ik, ip) + end do + 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 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(j=2:ipph, ik=2:idl1:2) + jc = ipp2 - j + 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 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 subroutine passb + + pure module subroutine passb2(ido, l1, cc, ch, wa1) + implicit none(type, external) + 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(k=1:l1, 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 + 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 + + pure module subroutine passb3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + 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(k=1:l1, 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)) + 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 + 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 + + pure module subroutine passb4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + 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(k=1:l1, 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 + cr3 = tr2 - tr3 + 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 + 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 + + pure module subroutine passb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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(k=1:l1, 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 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + dr3 = cr3 - ci4 + dr4 = cr3 + ci4 + di3 = ci3 + cr4 + di4 = ci3 - cr4 + dr5 = cr2 + ci5 + 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 + 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 + + pure module subroutine passf(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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 < l1) then + do concurrent(k=1:l1, j=2:ipph, i=1:ido) + jc = ipp2 - j + ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) + ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) + end do + do concurrent(k=1:l1, i=1:ido) + ch(i, k, 1) = cc(i, 1, k) + end do + else + do concurrent(k=1:l1, j=2:ipph, i=1:ido) + jc = ipp2 - j + ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) + ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) + end do + do concurrent(i=1:ido, k=1:l1) + ch(i, k, 1) = cc(i, 1, k) + end do + end if + idl = 2 - ido + inc = 0 + do l = 2, ipph + lc = ipp2 - l + idl = idl + ido + do concurrent(ik=1:idl1) + c2(ik, l) = ch2(ik, 1) + wa(idl - 1)*ch2(ik, 2) + c2(ik, lc) = -wa(idl)*ch2(ik, ip) + end do + 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 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(j=2:ipph, ik=2:idl1:2) + jc = ipp2 - j + 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 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 + + pure module subroutine passf2(ido, l1, cc, ch, wa1) + implicit none(type, external) + 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(k=1:l1, 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 + 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 + + pure module subroutine passf3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + 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(k=1:l1, 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)) + 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 + 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 + + pure module subroutine passf4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + 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(k=1:l1, 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 + cr3 = tr2 - tr3 + 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 + 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 + + pure module subroutine passf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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(k=1:l1, 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 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + dr3 = cr3 - ci4 + dr4 = cr3 + ci4 + di3 = ci3 + cr4 + di4 = ci3 - cr4 + dr5 = cr2 + ci5 + 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 + 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 +end module fftpack_legacy_drivers_pass diff --git a/src/fftpack_legacy_drivers_rad.f90 b/src/fftpack_legacy_drivers_rad.f90 new file mode 100644 index 00000000..4ea91691 --- /dev/null +++ b/src/fftpack_legacy_drivers_rad.f90 @@ -0,0 +1,811 @@ +module fftpack_legacy_drivers_rad + use fftpack_kinds, only: dp + implicit none(type, external) + private + + public :: radbg, radb2, radb3, radb4, radb5 + public :: radfg, radf2, radf3, radf4, radf5 + + interface + pure module subroutine radbg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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) + end subroutine radbg + end interface + + interface + pure module subroutine radb2(ido, l1, cc, ch, wa1) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 2, l1), wa1(*) + real(dp), intent(out) :: ch(ido, l1, 2) + end subroutine radb2 + end interface + + interface + pure module subroutine radb3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 3, l1), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, l1, 3) + end subroutine radb3 + end interface + + interface + pure module subroutine radb4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, 4, l1), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, l1, 4) + end subroutine radb4 + end interface + + interface + pure module subroutine radb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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) + end subroutine radb5 + end interface + + interface + pure module subroutine radfg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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) + end subroutine radfg + end interface + + interface + pure module subroutine radf2(ido, l1, cc, ch, wa1) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 2), wa1(*) + real(dp), intent(out) :: ch(ido, 2, l1) + end subroutine radf2 + end interface + + interface + pure module subroutine radf3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 3), wa1(*), wa2(*) + real(dp), intent(out) :: ch(ido, 3, l1) + end subroutine radf3 + end interface + + interface + pure module subroutine radf4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + integer, intent(in) :: ido, l1 + real(dp), intent(in) :: cc(ido, l1, 4), wa1(*), wa2(*), wa3(*) + real(dp), intent(out) :: ch(ido, 4, l1) + end subroutine radf4 + end interface + + interface + pure module subroutine radf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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) + end subroutine radf5 + end interface + +contains + + pure module subroutine radbg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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(k=1:l1, i=1:ido) + ch(i, k, 1) = cc(i, 1, k) + end do + else + do concurrent(k=1:l1, i=1:ido) + ch(i, k, 1) = cc(i, 1, k) + end do + end if + do concurrent(k=1:l1, j=2:ipph) + jc = ipp2 - j + j2 = j + j + ch(1, k, j) = cc(ido, j2 - 2, k) + cc(ido, j2 - 2, k) + ch(1, k, jc) = cc(1, j2 - 1, k) + cc(1, j2 - 1, k) + end do + if (ido /= 1) then + if (nbd < l1) then + do concurrent(k=1:l1, j=2:ipph, i=3:ido:2) + jc = ipp2 - j + 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) + end do + else + do concurrent(k=1:l1, j=2:ipph, i=3:ido:2) + jc = ipp2 - j + 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) + end do + end if + end if + ar1 = 1.0_dp + ai1 = 0.0_dp + do l = 2, ipph + lc = ipp2 - l + ar1h = dcp*ar1 - dsp*ai1 + ai1 = dcp*ai1 + dsp*ar1 + ar1 = ar1h + do concurrent(ik=1:idl1) + c2(ik, l) = ch2(ik, 1) + ar1*ch2(ik, 2) + c2(ik, lc) = ai1*ch2(ik, ip) + end do + dc2 = ar1 + ds2 = ai1 + ar2 = ar1 + ai2 = ai1 + do j = 3, ipph + jc = ipp2 - j + ar2h = dc2*ar2 - ds2*ai2 + ai2 = dc2*ai2 + ds2*ar2 + ar2 = ar2h + do concurrent(ik=1:idl1) + c2(ik, l) = c2(ik, l) + ar2*ch2(ik, j) + c2(ik, lc) = c2(ik, lc) + ai2*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(j=2:ipph, k=1:l1) + jc = ipp2 - j + ch(1, k, j) = c1(1, k, j) - c1(1, k, jc) + ch(1, k, jc) = c1(1, k, j) + c1(1, k, jc) + end do + if (ido /= 1) then + if (nbd < l1) then + do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) + jc = ipp2 - j + 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) + end do + else + do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) + jc = ipp2 - j + 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) + end do + end if + end if + if (ido == 1) return + 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) + 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 + 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 + 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) + end do + end do + end do + end if + end subroutine radbg + + pure module subroutine radb2(ido, l1, cc, ch, wa1) + implicit none(type, external) + 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 + 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 + + pure module subroutine radb3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + 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)) + 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 + end do + end subroutine radb3 + + pure module subroutine radb4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + 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 + cr3 = tr2 - tr3 + 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 + 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 + + pure module subroutine radb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + dr3 = cr3 - ci4 + dr4 = cr3 + ci4 + di3 = ci3 + cr4 + di4 = ci3 - cr4 + dr5 = cr2 + ci5 + 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 + end do + end subroutine radb5 + + pure module subroutine radfg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa) + implicit none(type, external) + 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 + 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) + end do + end do + end do + end if + if (nbd < l1) then + do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) + jc = ipp2 - j + c1(i - 1, k, j) = ch(i - 1, k, j) + ch(i - 1, k, jc) + c1(i - 1, k, jc) = ch(i, k, j) - ch(i, k, jc) + c1(i, k, j) = ch(i, k, j) + ch(i, k, jc) + c1(i, k, jc) = ch(i - 1, k, jc) - ch(i - 1, k, j) + end do + else + do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) + jc = ipp2 - j + c1(i - 1, k, j) = ch(i - 1, k, j) + ch(i - 1, k, jc) + c1(i - 1, k, jc) = ch(i, k, j) - ch(i, k, jc) + c1(i, k, j) = ch(i, k, j) + ch(i, k, jc) + c1(i, k, jc) = ch(i - 1, k, jc) - ch(i - 1, k, j) + end do + end if + end if + do concurrent(j=2:ipph, k=1:l1) + jc = ipp2 - j + c1(1, k, j) = ch(1, k, j) + ch(1, k, jc) + c1(1, k, jc) = ch(1, k, jc) - ch(1, k, j) + end do +! + ar1 = 1.0_dp + ai1 = 0.0_dp + do l = 2, ipph + lc = ipp2 - l + ar1h = dcp*ar1 - dsp*ai1 + ai1 = dcp*ai1 + dsp*ar1 + ar1 = ar1h + do concurrent(ik=1:idl1) + ch2(ik, l) = c2(ik, 1) + ar1*c2(ik, 2) + ch2(ik, lc) = ai1*c2(ik, ip) + end do + dc2 = ar1 + ds2 = ai1 + ar2 = ar1 + ai2 = ai1 + do j = 3, ipph + jc = ipp2 - j + ar2h = dc2*ar2 - ds2*ai2 + ai2 = dc2*ai2 + ds2*ar2 + ar2 = ar2h + do concurrent(ik=1:idl1) + ch2(ik, l) = ch2(ik, l) + ar2*c2(ik, j) + ch2(ik, lc) = ch2(ik, lc) + ai2*c2(ik, jc) + end do + end do + end do + do concurrent(j=2:ipph, ik=1:idl1) + ch2(ik, 1) = ch2(ik, 1) + c2(ik, j) + end do +! + if (ido < l1) then + do concurrent(k=1:l1, i=1:ido) + cc(i, 1, k) = ch(i, k, 1) + end do + else + do concurrent(i=1:ido, k=1:l1) + cc(i, 1, k) = ch(i, k, 1) + end do + end if + do concurrent(j=2:ipph, k=1:l1) + jc = ipp2 - j + j2 = j + j + cc(ido, j2 - 2, k) = ch(1, k, j) + cc(1, j2 - 1, k) = ch(1, k, jc) + end do + if (ido == 1) return + if (nbd < l1) then + do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) + jc = ipp2 - j + j2 = j + j + ic = idp2 - i + cc(i - 1, j2 - 1, k) = ch(i - 1, k, j) + ch(i - 1, k, jc) + cc(ic - 1, j2 - 2, k) = ch(i - 1, k, j) - ch(i - 1, k, jc) + cc(i, j2 - 1, k) = ch(i, k, j) + ch(i, k, jc) + cc(ic, j2 - 2, k) = ch(i, k, jc) - ch(i, k, j) + end do + else + do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) + jc = ipp2 - j + j2 = j + j + ic = idp2 - i + cc(i - 1, j2 - 1, k) = ch(i - 1, k, j) + ch(i - 1, k, jc) + cc(ic - 1, j2 - 2, k) = ch(i - 1, k, j) - ch(i - 1, k, jc) + cc(i, j2 - 1, k) = ch(i, k, j) + ch(i, k, jc) + cc(ic, j2 - 2, k) = ch(i, k, jc) - ch(i, k, j) + end do + end if + end subroutine radfg + + pure module subroutine radf2(ido, l1, cc, ch, wa1) + implicit none(type, external) + 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(k=1:l1, i=3:ido:2) + 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 + 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 + + pure module subroutine radf3(ido, l1, cc, ch, wa1, wa2) + implicit none(type, external) + 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(k=1:l1, i=3:ido:2) + 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) + 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 + end do + end subroutine radf3 + + pure module subroutine radf4(ido, l1, cc, ch, wa1, wa2, wa3) + implicit none(type, external) + 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(k=1:l1, i=3:ido:2) + 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) + 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 + 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 + + pure module subroutine radf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4) + implicit none(type, external) + 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(k=1:l1, i=3:ido:2) + 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) + cr2 = dr2 + dr5 + ci5 = dr5 - dr2 + cr5 = di2 - di5 + ci2 = di2 + di5 + cr3 = dr3 + dr4 + 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 + 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 + end do + end subroutine radf5 +end module fftpack_legacy_drivers_rad diff --git a/src/fftpack_rfft.f90 b/src/fftpack_rfft.f90 index 2d107670..f3a6b413 100644 --- a/src/fftpack_rfft.f90 +++ b/src/fftpack_rfft.f90 @@ -1,36 +1,36 @@ submodule(fftpack) fftpack_rfft - + implicit none contains - !> Forward transform of a real periodic sequence. - pure module function rfft_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) - - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, (0.0_rk, i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if - - !> Initialize FFT - lensav = 2*lenseq + 15 - allocate (wsave(lensav)) - call dffti(lenseq, wsave) - - !> Forward transformation - call dfftf(lenseq, result, wsave) - - end function rfft_rk + !> Forward transform of a real periodic sequence. + pure module function rfft_dp(x, n) result(result) + real(kind=dp), intent(in) :: x(:) + integer, intent(in), optional :: n + real(kind=dp), allocatable :: result(:) + + integer :: lenseq, lensav, i + real(kind=dp), allocatable :: wsave(:) + + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, (0.0_dp, i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if + + !> Initialize FFT + lensav = 2*lenseq + 15 + allocate (wsave(lensav)) + call dffti(lenseq, wsave) + + !> Forward transformation + call dfftf(lenseq, result, wsave) + + end function rfft_dp end submodule fftpack_rfft diff --git a/src/fftpack_utils.f90 b/src/fftpack_utils.f90 index 3ea491d1..932f8972 100644 --- a/src/fftpack_utils.f90 +++ b/src/fftpack_utils.f90 @@ -1,60 +1,60 @@ submodule(fftpack) fftpack_utils - + implicit none contains - !> Returns an integer array with the frequency values involved in the - !> performed FFT, ordered in the standard way (zero first, then positive - !> frequencies, then the negative ones). - pure module function fftfreq(n) result(out) - integer, intent(in) :: n - integer, dimension(n) :: out - integer :: i - - out(1) = 0 - if (n == 1) return - - if (mod(n, 2) == 0) then ! n even, smallest n = 2 - do i = 2, n/2 - out(i) = i-1 - end do - out(n/2+1) = -n/2 - do i = n/2+2, n ! only enters if n/2+2 <= n - out(i) = out(i-1) + 1 - end do - else ! n odd, smallest n = 3 - do i = 2, n/2+1 - out(i) = i-1 - end do - out(n/2+2) = -out(n/2+1) - do i = n/2+3, n ! only enters if n/2+3 <= n - out(i) = out(i-1) + 1 - end do - end if - end function fftfreq - - !> Returns an integer array with the frequency values involved in the - !> performed real FFT, ordered in the standard way (zero first, then - !> positive frequencies, then, if applicable, the negative one). - pure module function rfftfreq(n) result(out) - integer, intent(in) :: n - integer, dimension(n) :: out - integer :: i - - out(1) = 0 - if (n == 1) return - - if (mod(n,2) == 0) then ! n even, smallest n = 2 - do i = 2, n-2, 2 - out(i) = out(i-1) + 1 - out(i+1) = out(i) - end do - out(n) = -n/2 - else ! n odd, smallest n = 3 - do i = 2, n-1, 2 - out(i) = out(i-1) + 1 - out(i+1) = out(i) - end do - end if - end function rfftfreq + !> Returns an integer array with the frequency values involved in the + !> performed FFT, ordered in the standard way (zero first, then positive + !> frequencies, then the negative ones). + pure module function fftfreq(n) result(out) + integer, intent(in) :: n + integer, dimension(n) :: out + integer :: i + + out(1) = 0 + if (n == 1) return + + if (mod(n, 2) == 0) then ! n even, smallest n = 2 + do i = 2, n/2 + out(i) = i - 1 + end do + out(n/2 + 1) = -n/2 + do i = n/2 + 2, n ! only enters if n/2+2 <= n + out(i) = out(i - 1) + 1 + end do + else ! n odd, smallest n = 3 + do i = 2, n/2 + 1 + out(i) = i - 1 + end do + out(n/2 + 2) = -out(n/2 + 1) + do i = n/2 + 3, n ! only enters if n/2+3 <= n + out(i) = out(i - 1) + 1 + end do + end if + end function fftfreq + + !> Returns an integer array with the frequency values involved in the + !> performed real FFT, ordered in the standard way (zero first, then + !> positive frequencies, then, if applicable, the negative one). + pure module function rfftfreq(n) result(out) + integer, intent(in) :: n + integer, dimension(n) :: out + integer :: i + + out(1) = 0 + if (n == 1) return + + if (mod(n, 2) == 0) then ! n even, smallest n = 2 + do i = 2, n - 2, 2 + out(i) = out(i - 1) + 1 + out(i + 1) = out(i) + end do + out(n) = -n/2 + else ! n odd, smallest n = 3 + do i = 2, n - 1, 2 + out(i) = out(i - 1) + 1 + out(i + 1) = out(i) + end do + end if + end function rfftfreq end submodule fftpack_utils diff --git a/src/passb.f90 b/src/passb.f90 deleted file mode 100644 index 3a58494c..00000000 --- a/src/passb.f90 +++ /dev/null @@ -1,110 +0,0 @@ - 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 < l1) then - do concurrent(k=1:l1, j=2:ipph, i=1:ido) - jc = ipp2 - j - ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) - ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) - end do - do concurrent(k=1:l1, i=1:ido) - ch(i, k, 1) = cc(i, 1, k) - end do - else - do concurrent(k=1:l1, j=2:ipph, i=1:ido) - jc = ipp2 - j - ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) - ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) - end do - do concurrent(i=1:ido, k=1:l1) - ch(i, k, 1) = cc(i, 1, k) - end do - end if - idl = 2 - ido - inc = 0 - do l = 2, ipph - lc = ipp2 - l - idl = idl + ido - do concurrent(ik=1:idl1) - c2(ik, l) = ch2(ik, 1) + wa(idl - 1)*ch2(ik, 2) - c2(ik, lc) = wa(idl)*ch2(ik, ip) - end do - 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 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(j=2:ipph, ik=2:idl1:2) - jc = ipp2 - j - 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 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 deleted file mode 100644 index 1fd83a87..00000000 --- a/src/passb2.f90 +++ /dev/null @@ -1,26 +0,0 @@ - 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(k=1:l1, 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 - 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 deleted file mode 100644 index b4fd8c4d..00000000 --- a/src/passb3.f90 +++ /dev/null @@ -1,47 +0,0 @@ - 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(k=1:l1, 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)) - 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 - 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 deleted file mode 100644 index 294715a6..00000000 --- a/src/passb4.f90 +++ /dev/null @@ -1,55 +0,0 @@ - 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(k=1:l1, 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 - cr3 = tr2 - tr3 - 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 - 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 deleted file mode 100644 index a39621f5..00000000 --- a/src/passb5.f90 +++ /dev/null @@ -1,85 +0,0 @@ - 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(k=1:l1, 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 - cr5 = ti11*tr5 + ti12*tr4 - ci5 = ti11*ti5 + ti12*ti4 - cr4 = ti12*tr5 - ti11*tr4 - ci4 = ti12*ti5 - ti11*ti4 - dr3 = cr3 - ci4 - dr4 = cr3 + ci4 - di3 = ci3 + cr4 - di4 = ci3 - cr4 - dr5 = cr2 + ci5 - 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 - 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 deleted file mode 100644 index d78af30f..00000000 --- a/src/passf.f90 +++ /dev/null @@ -1,109 +0,0 @@ - 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 < l1) then - do concurrent(k=1:l1, j=2:ipph, i=1:ido) - jc = ipp2 - j - ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) - ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) - end do - do concurrent(k=1:l1, i=1:ido) - ch(i, k, 1) = cc(i, 1, k) - end do - else - do concurrent(k=1:l1, j=2:ipph, i=1:ido) - jc = ipp2 - j - ch(i, k, j) = cc(i, j, k) + cc(i, jc, k) - ch(i, k, jc) = cc(i, j, k) - cc(i, jc, k) - end do - do concurrent(i=1:ido, k=1:l1) - ch(i, k, 1) = cc(i, 1, k) - end do - end if - idl = 2 - ido - inc = 0 - do l = 2, ipph - lc = ipp2 - l - idl = idl + ido - do concurrent(ik=1:idl1) - c2(ik, l) = ch2(ik, 1) + wa(idl - 1)*ch2(ik, 2) - c2(ik, lc) = -wa(idl)*ch2(ik, ip) - end do - 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 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(j=2:ipph, ik=2:idl1:2) - jc = ipp2 - j - 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 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 deleted file mode 100644 index 550d954a..00000000 --- a/src/passf2.f90 +++ /dev/null @@ -1,26 +0,0 @@ - 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(k=1:l1, 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 - 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 deleted file mode 100644 index e987f36b..00000000 --- a/src/passf3.f90 +++ /dev/null @@ -1,47 +0,0 @@ - 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(k=1:l1, 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)) - 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 - 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 deleted file mode 100644 index 483dd117..00000000 --- a/src/passf4.f90 +++ /dev/null @@ -1,55 +0,0 @@ - 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(k=1:l1, 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 - cr3 = tr2 - tr3 - 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 - 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 deleted file mode 100644 index 49fc67db..00000000 --- a/src/passf5.f90 +++ /dev/null @@ -1,85 +0,0 @@ - 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(k=1:l1, 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 - cr5 = ti11*tr5 + ti12*tr4 - ci5 = ti11*ti5 + ti12*ti4 - cr4 = ti12*tr5 - ti11*tr4 - ci4 = ti12*ti5 - ti11*ti4 - dr3 = cr3 - ci4 - dr4 = cr3 + ci4 - di3 = ci3 + cr4 - di4 = ci3 - cr4 - dr5 = cr2 + ci5 - 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 - 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 deleted file mode 100644 index 7fb4c554..00000000 --- a/src/radb2.f90 +++ /dev/null @@ -1,31 +0,0 @@ - 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 - 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 deleted file mode 100644 index 6b9cba68..00000000 --- a/src/radb3.f90 +++ /dev/null @@ -1,41 +0,0 @@ - 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)) - 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 - end do - end subroutine radb3 diff --git a/src/radb4.f90 b/src/radb4.f90 deleted file mode 100644 index 4291e758..00000000 --- a/src/radb4.f90 +++ /dev/null @@ -1,62 +0,0 @@ - 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 - cr3 = tr2 - tr3 - 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 - 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 deleted file mode 100644 index c4de1f6b..00000000 --- a/src/radb5.f90 +++ /dev/null @@ -1,72 +0,0 @@ - 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 - cr5 = ti11*tr5 + ti12*tr4 - ci5 = ti11*ti5 + ti12*ti4 - cr4 = ti12*tr5 - ti11*tr4 - ci4 = ti12*ti5 - ti11*ti4 - dr3 = cr3 - ci4 - dr4 = cr3 + ci4 - di3 = ci3 + cr4 - di4 = ci3 - cr4 - dr5 = cr2 + ci5 - 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 - end do - end subroutine radb5 diff --git a/src/radbg.f90 b/src/radbg.f90 deleted file mode 100644 index 28f28c75..00000000 --- a/src/radbg.f90 +++ /dev/null @@ -1,147 +0,0 @@ - subroutine radbg(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) :: 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(k=1:l1, i=1:ido) - ch(i, k, 1) = cc(i, 1, k) - end do - else - do concurrent(k=1:l1, i=1:ido) - ch(i, k, 1) = cc(i, 1, k) - end do - end if - do concurrent(k=1:l1, j=2:ipph) - jc = ipp2 - j - j2 = j + j - ch(1, k, j) = cc(ido, j2 - 2, k) + cc(ido, j2 - 2, k) - ch(1, k, jc) = cc(1, j2 - 1, k) + cc(1, j2 - 1, k) - end do - if (ido /= 1) then - if (nbd < l1) then - do concurrent(k=1:l1, j=2:ipph, i=3:ido:2) - jc = ipp2 - j - 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) - end do - else - do concurrent(k=1:l1, j=2:ipph, i=3:ido:2) - jc = ipp2 - j - 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) - end do - end if - end if - ar1 = 1.0_dp - ai1 = 0.0_dp - do l = 2, ipph - lc = ipp2 - l - ar1h = dcp*ar1 - dsp*ai1 - ai1 = dcp*ai1 + dsp*ar1 - ar1 = ar1h - do concurrent(ik=1:idl1) - c2(ik, l) = ch2(ik, 1) + ar1*ch2(ik, 2) - c2(ik, lc) = ai1*ch2(ik, ip) - end do - dc2 = ar1 - ds2 = ai1 - ar2 = ar1 - ai2 = ai1 - do j = 3, ipph - jc = ipp2 - j - ar2h = dc2*ar2 - ds2*ai2 - ai2 = dc2*ai2 + ds2*ar2 - ar2 = ar2h - do concurrent(ik=1:idl1) - c2(ik, l) = c2(ik, l) + ar2*ch2(ik, j) - c2(ik, lc) = c2(ik, lc) + ai2*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(j=2:ipph, k=1:l1) - jc = ipp2 - j - ch(1, k, j) = c1(1, k, j) - c1(1, k, jc) - ch(1, k, jc) = c1(1, k, j) + c1(1, k, jc) - end do - if (ido /= 1) then - if (nbd < l1) then - do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) - jc = ipp2 - j - 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) - end do - else - do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) - jc = ipp2 - j - 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) - end do - end if - end if - if (ido == 1) return - 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) - 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 - 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 - 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) - end do - end do - end do - end if - end subroutine radbg diff --git a/src/radf2.f90 b/src/radf2.f90 deleted file mode 100644 index 10d6343e..00000000 --- a/src/radf2.f90 +++ /dev/null @@ -1,31 +0,0 @@ - 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(k=1:l1, i=3:ido:2) - 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 - 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 deleted file mode 100644 index cc27348e..00000000 --- a/src/radf3.f90 +++ /dev/null @@ -1,40 +0,0 @@ - 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(k=1:l1, i=3:ido:2) - 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) - 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 - end do - end subroutine radf3 diff --git a/src/radf4.f90 b/src/radf4.f90 deleted file mode 100644 index 63fe3610..00000000 --- a/src/radf4.f90 +++ /dev/null @@ -1,58 +0,0 @@ - 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(k=1:l1, i=3:ido:2) - 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) - 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 - 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 deleted file mode 100644 index 5ae40d4d..00000000 --- a/src/radf5.f90 +++ /dev/null @@ -1,68 +0,0 @@ - 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(k=1:l1, i=3:ido:2) - 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) - cr2 = dr2 + dr5 - ci5 = dr5 - dr2 - cr5 = di2 - di5 - ci2 = di2 + di5 - cr3 = dr3 + dr4 - 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 - 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 - end do - end subroutine radf5 diff --git a/src/radfg.f90 b/src/radfg.f90 deleted file mode 100644 index 2c9091cc..00000000 --- a/src/radfg.f90 +++ /dev/null @@ -1,156 +0,0 @@ - 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 - 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) - end do - end do - end do - end if - if (nbd < l1) then - do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) - jc = ipp2 - j - c1(i - 1, k, j) = ch(i - 1, k, j) + ch(i - 1, k, jc) - c1(i - 1, k, jc) = ch(i, k, j) - ch(i, k, jc) - c1(i, k, j) = ch(i, k, j) + ch(i, k, jc) - c1(i, k, jc) = ch(i - 1, k, jc) - ch(i - 1, k, j) - end do - else - do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) - jc = ipp2 - j - c1(i - 1, k, j) = ch(i - 1, k, j) + ch(i - 1, k, jc) - c1(i - 1, k, jc) = ch(i, k, j) - ch(i, k, jc) - c1(i, k, j) = ch(i, k, j) + ch(i, k, jc) - c1(i, k, jc) = ch(i - 1, k, jc) - ch(i - 1, k, j) - end do - end if - end if - do concurrent(j=2:ipph, k=1:l1) - jc = ipp2 - j - c1(1, k, j) = ch(1, k, j) + ch(1, k, jc) - c1(1, k, jc) = ch(1, k, jc) - ch(1, k, j) - end do -! - ar1 = 1.0_dp - ai1 = 0.0_dp - do l = 2, ipph - lc = ipp2 - l - ar1h = dcp*ar1 - dsp*ai1 - ai1 = dcp*ai1 + dsp*ar1 - ar1 = ar1h - do concurrent(ik=1:idl1) - ch2(ik, l) = c2(ik, 1) + ar1*c2(ik, 2) - ch2(ik, lc) = ai1*c2(ik, ip) - end do - dc2 = ar1 - ds2 = ai1 - ar2 = ar1 - ai2 = ai1 - do j = 3, ipph - jc = ipp2 - j - ar2h = dc2*ar2 - ds2*ai2 - ai2 = dc2*ai2 + ds2*ar2 - ar2 = ar2h - do concurrent(ik=1:idl1) - ch2(ik, l) = ch2(ik, l) + ar2*c2(ik, j) - ch2(ik, lc) = ch2(ik, lc) + ai2*c2(ik, jc) - end do - end do - end do - do concurrent(j=2:ipph, ik=1:idl1) - ch2(ik, 1) = ch2(ik, 1) + c2(ik, j) - end do -! - if (ido < l1) then - do concurrent(k=1:l1, i=1:ido) - cc(i, 1, k) = ch(i, k, 1) - end do - else - do concurrent(i=1:ido, k=1:l1) - cc(i, 1, k) = ch(i, k, 1) - end do - end if - do concurrent(j=2:ipph, k=1:l1) - jc = ipp2 - j - j2 = j + j - cc(ido, j2 - 2, k) = ch(1, k, j) - cc(1, j2 - 1, k) = ch(1, k, jc) - end do - if (ido == 1) return - if (nbd < l1) then - do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) - jc = ipp2 - j - j2 = j + j - ic = idp2 - i - cc(i - 1, j2 - 1, k) = ch(i - 1, k, j) + ch(i - 1, k, jc) - cc(ic - 1, j2 - 2, k) = ch(i - 1, k, j) - ch(i - 1, k, jc) - cc(i, j2 - 1, k) = ch(i, k, j) + ch(i, k, jc) - cc(ic, j2 - 2, k) = ch(i, k, jc) - ch(i, k, j) - end do - else - do concurrent(j=2:ipph, k=1:l1, i=3:ido:2) - jc = ipp2 - j - j2 = j + j - ic = idp2 - i - cc(i - 1, j2 - 1, k) = ch(i - 1, k, j) + ch(i - 1, k, jc) - cc(ic - 1, j2 - 2, k) = ch(i - 1, k, j) - ch(i - 1, k, jc) - cc(i, j2 - 1, k) = ch(i, k, j) + ch(i, k, jc) - cc(ic, j2 - 2, k) = ch(i, k, jc) - ch(i, k, j) - end do - end if - end subroutine radfg diff --git a/src/rfftb1.f90 b/src/rfftb1.f90 index 9f5ac533..7a8e36dc 100644 --- a/src/rfftb1.f90 +++ b/src/rfftb1.f90 @@ -1,5 +1,6 @@ subroutine rfftb1(n, c, ch, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp + use fftpack_legacy_drivers_rad, only: radbg, radb2, radb3, radb4, radb5 implicit none integer, intent(in) :: n real(dp), intent(inout) :: c(*) diff --git a/src/rfftf1.f90 b/src/rfftf1.f90 index 107e5598..78e144a1 100644 --- a/src/rfftf1.f90 +++ b/src/rfftf1.f90 @@ -1,5 +1,6 @@ subroutine rfftf1(n, c, ch, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp + use fftpack_legacy_drivers_rad, only: radfg, radf2, radf3, radf4, radf5 implicit none integer, intent(in) :: n real(dp), intent(inout) :: c(*) diff --git a/src/rffti1.f90 b/src/rffti1.f90 index 5ae8f418..702dfc04 100644 --- a/src/rffti1.f90 +++ b/src/rffti1.f90 @@ -1,5 +1,5 @@ subroutine rffti1(n, wa, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wa(*) diff --git a/src/rk.f90 b/src/rk.f90 deleted file mode 100644 index 3ab114da..00000000 --- a/src/rk.f90 +++ /dev/null @@ -1,4 +0,0 @@ -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 dddecd11..5252929c 100644 --- a/src/sint1.f90 +++ b/src/sint1.f90 @@ -1,5 +1,5 @@ subroutine sint1(n, war, was, xh, x, ifac) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n, ifac(*) real(dp), intent(in) :: was(*) diff --git a/src/zfftb.f90 b/src/zfftb.f90 index aae8cd5e..7123a963 100644 --- a/src/zfftb.f90 +++ b/src/zfftb.f90 @@ -1,5 +1,5 @@ subroutine zfftb(n, c, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: c(*) diff --git a/src/zfftf.f90 b/src/zfftf.f90 index 675b57e8..a647a420 100644 --- a/src/zfftf.f90 +++ b/src/zfftf.f90 @@ -1,5 +1,5 @@ subroutine zfftf(n, c, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(inout) :: c(*) diff --git a/src/zffti.f90 b/src/zffti.f90 index d6da6d8e..4cd1f849 100644 --- a/src/zffti.f90 +++ b/src/zffti.f90 @@ -1,5 +1,5 @@ subroutine zffti(n, wsave) - use fftpack_kind, only: dp => rk + use fftpack_kinds, only: dp implicit none integer, intent(in) :: n real(dp), intent(out) :: wsave(*) diff --git a/test/test_fftpack_dct.f90 b/test/test_fftpack_dct.f90 index c1f18614..f535949e 100644 --- a/test/test_fftpack_dct.f90 +++ b/test/test_fftpack_dct.f90 @@ -24,15 +24,15 @@ end subroutine collect_dct subroutine test_classic_dct(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: w(3*4 + 15), w2(3*4 + 15) - real(kind=rk) :: x(4) = [1, 2, 3, 4] - real(kind=rk) :: x2(4) - real(kind=rk) :: eps = 1.0e-10_rk + real(kind=dp) :: w(3*4 + 15), w2(3*4 + 15) + real(kind=dp) :: x(4) = [1, 2, 3, 4] + real(kind=dp) :: x2(4) + real(kind=dp) :: eps = 1.0e-10_dp x2 = x call dcosti(4, w) call dcost(4, x, w) - call check(error, sum(abs(x - [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk])) < eps, & + call check(error, sum(abs(x - [real(kind=dp) :: 15, -4, 0, -1.0000000000000009_dp])) < eps, & "`dcost` failed.") if (allocated(error)) return @@ -42,7 +42,7 @@ subroutine test_classic_dct(error) if (allocated(error)) return call dcost(4, x, w) - call check(error, sum(abs(x/(2*3) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & + call check(error, sum(abs(x/(2*3) - [real(kind=dp) :: 1, 2, 3, 4])) < eps, & "2nd `dcost` failed.") if (allocated(error)) return @@ -54,8 +54,8 @@ subroutine test_classic_dct(error) x2 = x call dcosqi(4, w) call dcosqf(4, x, w) - call check(error, sum(abs(x - [11.999626276085150_rk, -9.1029432177492193_rk, & - 2.6176618435106480_rk, -1.5143449018465791_rk])) < eps, & + call check(error, sum(abs(x - [11.999626276085150_dp, -9.1029432177492193_dp, & + 2.6176618435106480_dp, -1.5143449018465791_dp])) < eps, & "`dcosqf` failed.") if (allocated(error)) return @@ -65,7 +65,7 @@ subroutine test_classic_dct(error) if (allocated(error)) return call dcosqb(4, x, w) - call check(error, sum(abs(x/(4*4) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & + call check(error, sum(abs(x/(4*4) - [real(kind=dp) :: 1, 2, 3, 4])) < eps, & "`dcosqb` failed.") if (allocated(error)) return @@ -76,53 +76,53 @@ end subroutine test_classic_dct subroutine test_modernized_dct(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - real(kind=rk) :: x(3) = [9, -9, 3] + real(kind=dp) :: eps = 1.0e-10_dp + real(kind=dp) :: x(3) = [9, -9, 3] ! DCT-1 - call check(error, sum(abs(dct(x, 2, 1) - [0.0_rk, 18.0_rk])) < eps, "`dct(x,2,1)` failed.") + call check(error, sum(abs(dct(x, 2, 1) - [0.0_dp, 18.0_dp])) < eps, "`dct(x,2,1)` failed.") if (allocated(error)) return call check(error, sum(abs(dct(x, 3, 1) - dct(x, type=1))) < eps, "`dct(x,3,1)` failed.") if (allocated(error)) return - call check(error, sum(abs(dct(x, 4, 1) - [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33])) < eps, & + call check(error, sum(abs(dct(x, 4, 1) - [real(kind=dp) :: -3, -3.0000000000000036_dp, 15, 33])) < eps, & "`dct(x,4,1)` failed.") !DCT-2 x = [9, -9, 3] - call check(error, sum(abs(dct(x, 3, 2) - [12.0_rk, 20.784609690826525_rk, 60.0_rk])) < eps, & + call check(error, sum(abs(dct(x, 3, 2) - [12.0_dp, 20.784609690826525_dp, 60.0_dp])) < eps, & "`dct(x,3,2)` failed.") - call check(error, sum(abs(dct(x, 3) - [12.0_rk, 20.784609690826525_rk, 60.0_rk])) < eps, & + call check(error, sum(abs(dct(x, 3) - [12.0_dp, 20.784609690826525_dp, 60.0_dp])) < eps, & "`dct(x,3)` failed.") - call check(error, sum(abs(dct(x, 4, 2) - [12.0_rk, 14.890858416882008_rk, 42.426406871192853_rk, & - 58.122821125684993_rk])) < eps, "`dct(x,4,2)` failed.") + call check(error, sum(abs(dct(x, 4, 2) - [12.0_dp, 14.890858416882008_dp, 42.426406871192853_dp, & + 58.122821125684993_dp])) < eps, "`dct(x,4,2)` failed.") ! DCT-3 x = [9, -9, 3] - call check(error, sum(abs(dct(x, 2, 3) - [-3.7279220613578570_rk, 21.727922061357859_rk])) < eps, & + call check(error, sum(abs(dct(x, 2, 3) - [-3.7279220613578570_dp, 21.727922061357859_dp])) < eps, & "`dct(x,2,3)` failed.") if (allocated(error)) return call check(error, sum(abs(dct(x, 3, 3) - dct(x, type=3))) < eps, & "`dct(x,3,3)` failed.") if (allocated(error)) return - call check(error, sum(abs(dct(x, 4, 3) - [-3.3871908980838743_rk, -2.1309424696909023_rk, & - 11.645661095452331_rk, 29.872472272322447_rk])) < eps, & + call check(error, sum(abs(dct(x, 4, 3) - [-3.3871908980838743_dp, -2.1309424696909023_dp, & + 11.645661095452331_dp, 29.872472272322447_dp])) < eps, & "`dct(x,4,3)` failed.") end subroutine test_modernized_dct subroutine test_modernized_idct(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - real(kind=rk) :: x(4) = [1, 2, 3, 4] + real(kind=dp) :: eps = 1.0e-10_dp + real(kind=dp) :: x(4) = [1, 2, 3, 4] ! inverse DCT-1 call check(error, sum(abs(idct(dct(x, type=1), type=1)/(2*3) - x)) < eps, & "`idct(dct(x,type=1),type=1)/(2*3)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x, type=1), 2, 1)/(2*1) - [5.5_rk, 9.5_rk])) < eps, & + call check(error, sum(abs(idct(dct(x, type=1), 2, 1)/(2*1) - [5.5_dp, 9.5_dp])) < eps, & "`idct(dct(x,type=1), 2, 1)/(2*1)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x, 2, 1), 4, 1)/(2*3) - [0.16666666666666666_rk, 0.33333333333333331_rk, & - 0.66666666666666663_rk, 0.83333333333333315_rk])) < eps, & + call check(error, sum(abs(idct(dct(x, 2, 1), 4, 1)/(2*3) - [0.16666666666666666_dp, 0.33333333333333331_dp, & + 0.66666666666666663_dp, 0.83333333333333315_dp])) < eps, & "`idct(dct(x,2,1), 4, 1)/(2*3)` failed.") ! inverse DCT-2 @@ -130,11 +130,11 @@ subroutine test_modernized_idct(error) call check(error, sum(abs(idct(dct(x, type=2))/(4*4) - x)) < eps, & "`idct(dct(x, type=2))/(4*4)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x, type=2), n=2) - [22.156460020898692_rk, 57.843539979101308_rk])) < eps, & + call check(error, sum(abs(idct(dct(x, type=2), n=2) - [22.156460020898692_dp, 57.843539979101308_dp])) < eps, & "`idct(dct(x, type=2),n=2)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x, n=2, type=2), n=4) - [6.7737481404944937_rk, 9.8352155994152106_rk, & - 14.164784400584789_rk, 17.226251859505506_rk])) < eps, & + call check(error, sum(abs(idct(dct(x, n=2, type=2), n=4) - [6.7737481404944937_dp, 9.8352155994152106_dp, & + 14.164784400584789_dp, 17.226251859505506_dp])) < eps, & "`idct(dct(x, n=2, type=2), n=4)` failed.") ! inverse DCT-3 @@ -143,11 +143,11 @@ subroutine test_modernized_idct(error) "`idct(dct(x, type=3), type=3)/(4*4)` failed.") if (allocated(error)) return call check(error, sum(abs(idct(dct(x, type=3), n=2, type=3)/(4*2) - & - [1.4483415291679655_rk, 7.4608849947753271_rk])) < eps, & + [1.4483415291679655_dp, 7.4608849947753271_dp])) < eps, & "`idct(dct(x, type=3), n=2, type=3)/(4*2)` failed.") if (allocated(error)) return call check(error, sum(abs(idct(dct(x, n=2, type=3), n=4, type=3)/(4*4) - & - [0.5_rk, 0.70932417358418376_rk, 1.0_rk, 0.78858050747473762_rk])) < eps, & + [0.5_dp, 0.70932417358418376_dp, 1.0_dp, 0.78858050747473762_dp])) < eps, & "`idct(dct(x, n=2, type=3), n=4, type=3)/(4*4)` failed.") end subroutine test_modernized_idct diff --git a/test/test_fftpack_fft.f90 b/test/test_fftpack_fft.f90 index 14634f56..6d65498f 100644 --- a/test/test_fftpack_fft.f90 +++ b/test/test_fftpack_fft.f90 @@ -1,6 +1,6 @@ module test_fftpack_fft - use fftpack, only: rk, zffti, zfftf, zfftb, fft, ifft + use fftpack, only: dp, zffti, zfftf, zfftb, fft, ifft use testdrive, only: new_unittest, unittest_type, error_type, check implicit none private @@ -24,48 +24,48 @@ end subroutine collect_fft subroutine test_classic_fft(error) type(error_type), allocatable, intent(out) :: error - complex(kind=rk) :: x(4) = [1, 2, 3, 4] - real(kind=rk) :: w(31) + complex(kind=dp) :: x(4) = [1, 2, 3, 4] + real(kind=dp) :: w(31) call zffti(4, w) call zfftf(4, x, w) - call check(error, all(x == [complex(kind=rk) ::(10, 0), (-2, 2), (-2, 0), (-2, -2)]), & + call check(error, all(x == [complex(kind=dp) ::(10, 0), (-2, 2), (-2, 0), (-2, -2)]), & "`zfftf` failed.") if (allocated(error)) return call zfftb(4, x, w) - call check(error, all(x/4.0_rk == [complex(kind=rk) ::(1, 0), (2, 0), (3, 0), (4, 0)]), & + call check(error, all(x/4.0_dp == [complex(kind=dp) ::(1, 0), (2, 0), (3, 0), (4, 0)]), & "`zfftb` failed.") end subroutine test_classic_fft subroutine test_modernized_fft(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - complex(kind=rk) :: x(3) = [1.0_rk, 2.0_rk, 3.0_rk] + real(kind=dp) :: eps = 1.0e-10_dp + complex(kind=dp) :: x(3) = [1.0_dp, 2.0_dp, 3.0_dp] - call check(error, sum(abs(fft(x, 2) - [(3.0_rk, 0.0_rk), (-1.0_rk, 0.0_rk)])) < eps, & + call check(error, sum(abs(fft(x, 2) - [(3.0_dp, 0.0_dp), (-1.0_dp, 0.0_dp)])) < eps, & "`fft(x, 2)` failed.") if (allocated(error)) return call check(error, sum(abs(fft(x, 3) - fft(x))) < eps, & "`fft(x, 3)` failed.") if (allocated(error)) return - call check(error, sum(abs(fft(x, 4) - [(6.0_rk, 0.0_rk), (-2.0_rk, -2.0_rk), (2.0_rk, 0.0_rk), (-2.0_rk, 2.0_rk)])) < eps, & + call check(error, sum(abs(fft(x, 4) - [(6.0_dp, 0.0_dp), (-2.0_dp, -2.0_dp), (2.0_dp, 0.0_dp), (-2.0_dp, 2.0_dp)])) < eps, & "`fft(x, 4)` failed.") end subroutine test_modernized_fft subroutine test_modernized_ifft(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - complex(kind=rk) :: x(4) = [1, 2, 3, 4] + real(kind=dp) :: eps = 1.0e-10_dp + complex(kind=dp) :: x(4) = [1, 2, 3, 4] - call check(error, sum(abs(ifft(fft(x))/4.0_rk - [complex(kind=rk) :: 1, 2, 3, 4])) < eps, & - "`ifft(fft(x))/4.0_rk` failed.") + call check(error, sum(abs(ifft(fft(x))/4.0_dp - [complex(kind=dp) :: 1, 2, 3, 4])) < eps, & + "`ifft(fft(x))/4.0_dp` failed.") if (allocated(error)) return - call check(error, sum(abs(ifft(fft(x), 2) - [complex(kind=rk) ::(8, 2), (12, -2)])) < eps, & + call check(error, sum(abs(ifft(fft(x), 2) - [complex(kind=dp) ::(8, 2), (12, -2)])) < eps, & "`ifft(fft(x), 2)` failed.") if (allocated(error)) return - call check(error, sum(abs(ifft(fft(x, 2), 4) - [complex(kind=rk) ::(2, 0), (3, -1), (4, 0), (3, 1)])) < eps, & + call check(error, sum(abs(ifft(fft(x, 2), 4) - [complex(kind=dp) ::(2, 0), (3, -1), (4, 0), (3, 1)])) < eps, & "`ifft(fft(x, 2), 4)` failed.") end subroutine test_modernized_ifft diff --git a/test/test_fftpack_original.f90 b/test/test_fftpack_original.f90 index 1cafa1e1..e439426e 100644 --- a/test/test_fftpack_original.f90 +++ b/test/test_fftpack_original.f90 @@ -6,9 +6,9 @@ module test_fftpack_original public :: collect_original - real(rk), parameter :: pi = 4.0_rk*atan(1.0_rk) - real(rk), parameter :: atol = epsilon(1.0_rk) - real(rk), parameter :: rtol = sqrt(atol) + real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp) + real(dp), parameter :: atol = epsilon(1.0_dp) + real(dp), parameter :: rtol = sqrt(atol) integer, parameter :: nd(1:7) = [120, 54, 49, 32, 4, 3, 2] contains @@ -25,10 +25,10 @@ end subroutine collect_original subroutine test_dfft(error) type(error_type), allocatable, intent(out) :: error - real(rk) :: x(200), y(200), xh(200), w(2000) + real(dp) :: x(200), y(200), xh(200), w(2000) integer :: i, j, k, n, np1, nm1, ns2, nz, modn - real(rk) :: dt, sum1, sum2, arg, arg1 - real(rk) :: mismatch, cf + real(dp) :: dt, sum1, sum2, arg, arg1 + real(dp) :: mismatch, cf do nz = 1, size(nd) !> Create multisine signal. @@ -36,7 +36,7 @@ subroutine test_dfft(error) modn = mod(n, 2) np1 = n + 1; nm1 = n - 1 do concurrent(j=1:np1) - x(j) = sin(j*sqrt(2.0_rk)) + x(j) = sin(j*sqrt(2.0_dp)) end do y = x; xh = x @@ -45,7 +45,7 @@ subroutine test_dfft(error) ns2 = (n + 1)/2 if (ns2 >= 2) then do k = 2, ns2 - sum1 = 0.0_rk; sum2 = 0.0_rk + sum1 = 0.0_dp; sum2 = 0.0_dp arg = (k - 1)*dt do i = 1, n arg1 = (i - 1)*arg @@ -74,7 +74,7 @@ subroutine test_dfft(error) !> Inverse Discrete Fourier Transform. x(:n) = xh(:n) ! Restore signal. do i = 1, n - sum1 = 0.5_rk*x(1) + sum1 = 0.5_dp*x(1) arg = (i - 1)*dt if (ns2 >= 2) then do k = 2, ns2 @@ -82,7 +82,7 @@ subroutine test_dfft(error) sum1 = sum1 + x(2*k - 2)*cos(arg1) - x(2*k - 1)*sin(arg1) end do end if - if (modn == 0) sum1 = sum1 + 0.5_rk*(-1)**(i - 1)*x(n) + if (modn == 0) sum1 = sum1 + 0.5_dp*(-1)**(i - 1)*x(n) y(i) = 2*sum1 end do @@ -100,7 +100,7 @@ subroutine test_dfft(error) call dfftf(n, y, w) !> Check error. - cf = 1.0_rk/n + cf = 1.0_dp/n mismatch = maxval(abs(cf*y(:n) - x(:n))) call check(error, mismatch < rtol) if (allocated(error)) return @@ -110,24 +110,24 @@ end subroutine test_dfft subroutine test_zfft(error) type(error_type), allocatable, intent(out) :: error integer :: i, j, k, n, nz - complex(rk) :: cx(200), cy(200) - real(rk) :: w(2000), dt, arg1, arg2, mismatch, cf + complex(dp) :: cx(200), cy(200) + real(dp) :: w(2000), dt, arg1, arg2, mismatch, cf do nz = 1, size(nd) !> Create signal. n = nd(nz) do concurrent(i=1:n) - cx(i) = cmplx(cos(sqrt(2.0_rk)*i), sin(sqrt(2.0_rk)*i**2), kind=rk) + cx(i) = cmplx(cos(sqrt(2.0_dp)*i), sin(sqrt(2.0_dp)*i**2), kind=dp) end do !> Discrete Fourier Transform. dt = 2*pi/n do i = 1, n arg1 = -(i - 1)*dt - cy(i) = cmplx(0.0_rk, 0.0_rk, kind=rk) + cy(i) = cmplx(0.0_dp, 0.0_dp, kind=dp) do k = 1, n arg2 = (k - 1)*arg1 - cy(i) = cy(i) + cmplx(cos(arg2), sin(arg2), kind=rk)*cx(k) + cy(i) = cy(i) + cmplx(cos(arg2), sin(arg2), kind=dp)*cx(k) end do end do @@ -144,10 +144,10 @@ subroutine test_zfft(error) cx(:n) = cx(:n)/n ! Scale signal. do i = 1, n arg1 = (i - 1)*dt - cy(i) = cmplx(0.0_rk, 0.0_rk, kind=rk) + cy(i) = cmplx(0.0_dp, 0.0_dp, kind=dp) do k = 1, n arg2 = (k - 1)*arg1 - cy(i) = cy(i) + cmplx(cos(arg2), sin(arg2), kind=rk)*cx(k) + cy(i) = cy(i) + cmplx(cos(arg2), sin(arg2), kind=dp)*cx(k) end do end do @@ -165,7 +165,7 @@ subroutine test_zfft(error) call zfftb(n, cx, w) !> Check error. - cf = 1.0_rk/n + cf = 1.0_dp/n mismatch = maxval(abs(cf*cx(:n) - cy(:n))) call check(error, mismatch < rtol) if (allocated(error)) return @@ -174,17 +174,17 @@ end subroutine test_zfft subroutine test_sint(error) type(error_type), allocatable, intent(out) :: error - real(rk) :: x(200), y(200), xh(200), w(2000) + real(dp) :: x(200), y(200), xh(200), w(2000) integer :: i, j, k, n, np1, nm1, ns2, nz - real(rk) :: dt, sum1, sum2, arg, arg1 - real(rk) :: mismatch, cf + real(dp) :: dt, sum1, sum2, arg, arg1 + real(dp) :: mismatch, cf do nz = 1, size(nd) !> Create multisine signal. n = nd(nz) np1 = n + 1; nm1 = n - 1 do concurrent(j=1:np1) - x(j) = sin(j*sqrt(2.0_rk)) + x(j) = sin(j*sqrt(2.0_dp)) end do y = x; xh = x @@ -192,7 +192,7 @@ subroutine test_sint(error) dt = pi/n do i = 1, nm1 - y(i) = 0.0_rk + y(i) = 0.0_dp arg1 = i*dt do k = 1, nm1 y(i) = y(i) + x(k)*sin(k*arg1) @@ -205,7 +205,7 @@ subroutine test_sint(error) call dsint(nm1, x, w) !> Check error. - cf = 0.5_rk/n + cf = 0.5_dp/n mismatch = maxval(abs(x(:nm1) - y(:nm1)))*cf call check(error, mismatch < rtol) if (allocated(error)) return @@ -224,17 +224,17 @@ end subroutine test_sint subroutine test_cost(error) type(error_type), allocatable, intent(out) :: error - real(rk) :: x(200), y(200), xh(200), w(2000) + real(dp) :: x(200), y(200), xh(200), w(2000) integer :: i, j, k, n, np1, nm1, ns2, nz - real(rk) :: dt, sum1, sum2, arg, arg1 - real(rk) :: mismatch, cf + real(dp) :: dt, sum1, sum2, arg, arg1 + real(dp) :: mismatch, cf do nz = 1, size(nd) !> Create multisine signal. n = nd(nz) np1 = n + 1; nm1 = n - 1 do concurrent(j=1:np1) - x(j) = sin(j*sqrt(2.0_rk)) + x(j) = sin(j*sqrt(2.0_dp)) end do y = x; xh = x @@ -242,7 +242,7 @@ subroutine test_cost(error) dt = pi/n do i = 1, np1 - y(i) = 0.5_rk*(x(1) + (-1)**(i + 1)*x(n + 1)) + y(i) = 0.5_dp*(x(1) + (-1)**(i + 1)*x(n + 1)) arg = (i - 1)*dt do k = 2, n y(i) = y(i) + cos((k - 1)*arg)*x(k) @@ -255,7 +255,7 @@ subroutine test_cost(error) call dcost(np1, x, w) !> Check error. - cf = 0.5_rk/n + cf = 0.5_dp/n mismatch = maxval(abs(x(:np1) - y(:np1)))*cf call check(error, mismatch < rtol) if (allocated(error)) return @@ -274,17 +274,17 @@ end subroutine test_cost subroutine test_cosqt(error) type(error_type), allocatable, intent(out) :: error - real(rk) :: x(200), y(200), xh(200), w(2000) + real(dp) :: x(200), y(200), xh(200), w(2000) integer :: i, j, k, n, np1, nm1, ns2, nz - real(rk) :: dt, sum1, sum2, arg, arg1 - real(rk) :: mismatch, cf + real(dp) :: dt, sum1, sum2, arg, arg1 + real(dp) :: mismatch, cf do nz = 1, size(nd) !> Create multisine signal. n = nd(nz) np1 = n + 1; nm1 = n - 1 do concurrent(j=1:np1) - x(j) = sin(j*sqrt(2.0_rk)) + x(j) = sin(j*sqrt(2.0_dp)) end do y = x; xh = x @@ -292,7 +292,7 @@ subroutine test_cosqt(error) dt = pi/(2*n) do i = 1, n - x(i) = 0.0_rk + x(i) = 0.0_dp arg = (i - 1)*dt do k = 1, n x(i) = x(i) + y(k)*cos((2*k - 1)*arg) @@ -305,7 +305,7 @@ subroutine test_cosqt(error) call dcosqb(n, y, w) !> Check error. - cf = 0.25_rk/n + cf = 0.25_dp/n mismatch = maxval(abs(x(:n) - y(:n)))*cf call check(error, mismatch < rtol) if (allocated(error)) return @@ -313,7 +313,7 @@ subroutine test_cosqt(error) !> Discrete inverse quarter-cos transform. x(:n) = xh(:n) do i = 1, n - y(i) = 0.5_rk*x(1) + y(i) = 0.5_dp*x(1) arg = (2*i - 1)*dt do k = 2, n y(i) = y(i) + x(k)*cos((k - 1)*arg) @@ -344,12 +344,12 @@ end subroutine test_cosqt subroutine test_dzfft(error) type(error_type), allocatable, intent(out) :: error - real(rk) :: x(200), y(200), xh(200), w(2000) - real(rk) :: a(100), b(100), ah(100), bh(100) - real(rk) :: azero, azeroh + real(dp) :: x(200), y(200), xh(200), w(2000) + real(dp) :: a(100), b(100), ah(100), bh(100) + real(dp) :: azero, azeroh integer :: i, j, k, n, np1, nm1, ns2, ns2m, nz, modn - real(rk) :: dt, sum1, sum2, arg, arg1, arg2 - real(rk) :: mismatch, cf + real(dp) :: dt, sum1, sum2, arg, arg1, arg2 + real(dp) :: mismatch, cf do nz = 1, size(nd) !> Create multisine signal. @@ -357,7 +357,7 @@ subroutine test_dzfft(error) modn = mod(n, 2) np1 = n + 1; nm1 = n - 1 do concurrent(j=1:np1) - x(j) = sin(j*sqrt(2.0_rk)) + x(j) = sin(j*sqrt(2.0_dp)) end do y = x; xh = x @@ -365,10 +365,10 @@ subroutine test_dzfft(error) dt = 2*pi/n ns2 = (n + 1)/2 ns2m = ns2 - 1 - cf = 2.0_rk/n + cf = 2.0_dp/n if (ns2m > 0) then do k = 1, ns2m - sum1 = 0.0_rk; sum2 = 0.0_rk + sum1 = 0.0_dp; sum2 = 0.0_dp arg = k*dt do i = 1, n arg1 = (i - 1)*arg @@ -383,8 +383,8 @@ subroutine test_dzfft(error) sum1 = sum(x(1:nm1:2)) sum2 = sum(x(2:nm1 + 1:2)) if (modn == 1) sum1 = sum1 + x(n) - azero = 0.5_rk*cf*(sum1 + sum2) - if (modn == 0) a(ns2) = 0.5_rk*cf*(sum1 - sum2) + azero = 0.5_dp*cf*(sum1 + sum2) + if (modn == 0) a(ns2) = 0.5_dp*cf*(sum1 - sum2) !> Fast Fourier Transform. call dzffti(n, w) @@ -403,7 +403,7 @@ subroutine test_dzfft(error) !> Inverse Discrete Fourier Transform ns2 = n/2 - if (modn == 0) b(ns2) = 0.0_rk + if (modn == 0) b(ns2) = 0.0_dp do i = 1, n sum1 = azero arg1 = (i - 1)*dt diff --git a/test/test_fftpack_parseval.f90 b/test/test_fftpack_parseval.f90 index 417bac17..f76f6ec4 100644 --- a/test/test_fftpack_parseval.f90 +++ b/test/test_fftpack_parseval.f90 @@ -1,14 +1,14 @@ module test_fftpack_parseval use fftpack, only: fft - use fftpack_kind, only: rk + use fftpack_kinds, only: dp use testdrive, only: new_unittest, unittest_type, error_type, check implicit none private public :: collect_parseval integer, parameter :: n = 1021 - real(rk), parameter :: atol = n*epsilon(1.0_rk) - real(rk), parameter :: rtol = sqrt(atol) + real(dp), parameter :: atol = n*epsilon(1.0_dp) + real(dp), parameter :: rtol = sqrt(atol) contains @@ -25,14 +25,14 @@ subroutine test_plancherel(error) !----- COMPLEX(DP) TRANSFORM ----- !-------------------------------------- block - complex(rk) :: x(n) - real(rk) :: xre(n), xim(n) - real(rk) :: norms(2) - call random_number(xre); call random_number(xim); x = cmplx(xre, xim, kind=rk) + complex(dp) :: x(n) + real(dp) :: xre(n), xim(n) + real(dp) :: norms(2) + call random_number(xre); call random_number(xim); x = cmplx(xre, xim, kind=dp) !> Norm of the signal in the time domain. norms(1) = sum(abs(x)**2) !> Normalize signal. - x = x/sqrt(norms(1)); norms(1) = 1.0_rk + x = x/sqrt(norms(1)); norms(1) = 1.0_dp !> Fourier transform. x = fft(x) !> Norm of the signal in frequency domain. @@ -49,10 +49,10 @@ subroutine test_parseval(error) !----- COMPLEX(DP) TRANSFORM ----- !-------------------------------------- block - complex(rk) :: x(n), y(n), dotprods(2) - real(rk) :: xre(n), xim(n) - call random_number(xre); call random_number(xim); x = cmplx(xre, xim, kind=rk) - call random_number(xre); call random_number(xim); y = cmplx(xre, xim, kind=rk) + complex(dp) :: x(n), y(n), dotprods(2) + real(dp) :: xre(n), xim(n) + call random_number(xre); call random_number(xim); x = cmplx(xre, xim, kind=dp) + call random_number(xre); call random_number(xim); y = cmplx(xre, xim, kind=dp) !> Normalize signals. x = x/sqrt(sum(abs(x)**2)); y = y/sqrt(sum(abs(y)**2)) !> Dot product in time domain. diff --git a/test/test_fftpack_rfft.f90 b/test/test_fftpack_rfft.f90 index 51b81d96..aae949ff 100644 --- a/test/test_fftpack_rfft.f90 +++ b/test/test_fftpack_rfft.f90 @@ -1,6 +1,6 @@ module test_fftpack_rfft - use fftpack, only: rk, dffti, dfftf, dfftb, rfft, irfft + use fftpack, only: dp, dffti, dfftf, dfftb, rfft, irfft use fftpack, only: dzffti, dzfftf, dzfftb use testdrive, only: new_unittest, unittest_type, error_type, check implicit none @@ -25,64 +25,64 @@ end subroutine collect_rfft subroutine test_classic_rfft(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: x(4) = [1, 2, 3, 4] - real(kind=rk) :: w(31) - real(kind=rk) :: azero, a(4/2), b(4/2) + real(kind=dp) :: x(4) = [1, 2, 3, 4] + real(kind=dp) :: w(31) + real(kind=dp) :: azero, a(4/2), b(4/2) call dffti(4, w) call dfftf(4, x, w) - call check(error, all(x == [real(kind=rk) :: 10, -2, 2, -2]), & + call check(error, all(x == [real(kind=dp) :: 10, -2, 2, -2]), & "`dfftf` failed.") if (allocated(error)) return call dfftb(4, x, w) - call check(error, all(x/4.0_rk == [real(kind=rk) :: 1, 2, 3, 4]), & + call check(error, all(x/4.0_dp == [real(kind=dp) :: 1, 2, 3, 4]), & "`dfftb` failed.") if (allocated(error)) return x = [1, 2, 3, 4] call dzffti(4, w) call dzfftf(4, x, azero, a, b, w) - call check(error, azero == 2.5_rk, "dzfftf: azero == 2.5_rk failed.") + call check(error, azero == 2.5_dp, "dzfftf: azero == 2.5_dp failed.") if (allocated(error)) return - call check(error, all(a == [-1.0_rk, -0.5_rk]), "dzfftf: all(a == [-1.0, -0.5]) failed.") + call check(error, all(a == [-1.0_dp, -0.5_dp]), "dzfftf: all(a == [-1.0, -0.5]) failed.") if (allocated(error)) return - call check(error, all(b == [-1.0_rk, 0.0_rk]), "dzfftf: all(b == [-1.0, 0.0]) failed.") + call check(error, all(b == [-1.0_dp, 0.0_dp]), "dzfftf: all(b == [-1.0, 0.0]) failed.") if (allocated(error)) return call dzfftb(4, x, azero, a, b, w) - call check(error, all(x == [real(kind=rk) :: 1, 2, 3, 4]), & - "dzfftb: all(x = [real(kind=rk) :: 1, 2, 3, 4]) failed.") + call check(error, all(x == [real(kind=dp) :: 1, 2, 3, 4]), & + "dzfftb: all(x = [real(kind=dp) :: 1, 2, 3, 4]) failed.") end subroutine test_classic_rfft subroutine test_modernized_rfft(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - real(kind=rk) :: x(3) = [9, -9, 3] + real(kind=dp) :: eps = 1.0e-10_dp + real(kind=dp) :: x(3) = [9, -9, 3] - call check(error, sum(abs(rfft(x, 2) - [real(kind=rk) :: 0, 18])) < eps, & + call check(error, sum(abs(rfft(x, 2) - [real(kind=dp) :: 0, 18])) < eps, & "`rfft(x, 2)` failed.") if (allocated(error)) return call check(error, sum(abs(rfft(x, 3) - rfft(x))) < eps, & "`rfft(x, 3)` failed.") if (allocated(error)) return - call check(error, sum(abs(rfft(x, 4) - [real(kind=rk) :: 3, 6, 9, 21])) < eps, & + call check(error, sum(abs(rfft(x, 4) - [real(kind=dp) :: 3, 6, 9, 21])) < eps, & "`rfft(x, 4)` failed.") end subroutine test_modernized_rfft subroutine test_modernized_irfft(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - real(kind=rk) :: x(4) = [1, 2, 3, 4] + real(kind=dp) :: eps = 1.0e-10_dp + real(kind=dp) :: x(4) = [1, 2, 3, 4] - call check(error, sum(abs(irfft(rfft(x))/4.0_rk - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & - "`irfft(rfft(x))/4.0_rk` failed.") + call check(error, sum(abs(irfft(rfft(x))/4.0_dp - [real(kind=dp) :: 1, 2, 3, 4])) < eps, & + "`irfft(rfft(x))/4.0_dp` failed.") if (allocated(error)) return - call check(error, sum(abs(irfft(rfft(x), 2) - [real(kind=rk) :: 8, 12])) < eps, & + call check(error, sum(abs(irfft(rfft(x), 2) - [real(kind=dp) :: 8, 12])) < eps, & "`irfft(rfft(x), 2)` failed.") if (allocated(error)) return - call check(error, sum(abs(irfft(rfft(x, 2), 4) - [real(kind=rk) :: 1, 3, 5, 3])) < eps, & + call check(error, sum(abs(irfft(rfft(x, 2), 4) - [real(kind=dp) :: 1, 3, 5, 3])) < eps, & "`irfft(rfft(x, 2), 4)` failed.") end subroutine test_modernized_irfft diff --git a/test/test_fftpack_utils.f90 b/test/test_fftpack_utils.f90 index ed7438c2..df21ccd2 100644 --- a/test/test_fftpack_utils.f90 +++ b/test/test_fftpack_utils.f90 @@ -1,6 +1,6 @@ module test_fftpack_utils - use fftpack, only: rk, fft, ifft, fftshift, ifftshift, fftfreq, rfftfreq + use fftpack, only: dp, fft, ifft, fftshift, ifftshift, fftfreq, rfftfreq use testdrive, only: new_unittest, unittest_type, error_type, check implicit none private @@ -29,27 +29,27 @@ end subroutine collect_utils subroutine test_fftshift_complex(error) type(error_type), allocatable, intent(out) :: error - complex(kind=rk) :: xeven(4) = [1, 2, 3, 4] - complex(kind=rk) :: xodd(5) = [1, 2, 3, 4, 5] + complex(kind=dp) :: xeven(4) = [1, 2, 3, 4] + complex(kind=dp) :: xodd(5) = [1, 2, 3, 4, 5] - call check(error, all(fftshift(xeven) == [complex(kind=rk) :: 3, 4, 1, 2]), & - "all(fftshift(xeven) == [complex(kind=rk) :: 3, 4, 1, 2]) failed.") + call check(error, all(fftshift(xeven) == [complex(kind=dp) :: 3, 4, 1, 2]), & + "all(fftshift(xeven) == [complex(kind=dp) :: 3, 4, 1, 2]) failed.") if (allocated(error)) return - call check(error, all(fftshift(xodd) == [complex(kind=rk) :: 4, 5, 1, 2, 3]), & - "all(fftshift(xodd) == [complex(kind=rk) :: 4, 5, 1, 2, 3]) failed.") + call check(error, all(fftshift(xodd) == [complex(kind=dp) :: 4, 5, 1, 2, 3]), & + "all(fftshift(xodd) == [complex(kind=dp) :: 4, 5, 1, 2, 3]) failed.") end subroutine test_fftshift_complex subroutine test_fftshift_real(error) type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: xeven(4) = [1, 2, 3, 4] - real(kind=rk) :: xodd(5) = [1, 2, 3, 4, 5] + real(kind=dp) :: xeven(4) = [1, 2, 3, 4] + real(kind=dp) :: xodd(5) = [1, 2, 3, 4, 5] - call check(error, all(fftshift(xeven) == [real(kind=rk) :: 3, 4, 1, 2]), & - "all(fftshift(xeven) == [real(kind=rk) :: 3, 4, 1, 2]) failed.") + call check(error, all(fftshift(xeven) == [real(kind=dp) :: 3, 4, 1, 2]), & + "all(fftshift(xeven) == [real(kind=dp) :: 3, 4, 1, 2]) failed.") if (allocated(error)) return - call check(error, all(fftshift(xodd) == [real(kind=rk) :: 4, 5, 1, 2, 3]), & - "all(fftshift(xodd) == [real(kind=rk) :: 4, 5, 1, 2, 3]) failed.") + call check(error, all(fftshift(xodd) == [real(kind=dp) :: 4, 5, 1, 2, 3]), & + "all(fftshift(xodd) == [real(kind=dp) :: 4, 5, 1, 2, 3]) failed.") end subroutine test_fftshift_real @@ -57,14 +57,14 @@ subroutine test_ifftshift_complex(error) type(error_type), allocatable, intent(out) :: error integer :: i - complex(kind=rk) :: xeven(4) = [3, 4, 1, 2] - complex(kind=rk) :: xodd(5) = [4, 5, 1, 2, 3] + complex(kind=dp) :: xeven(4) = [3, 4, 1, 2] + complex(kind=dp) :: xodd(5) = [4, 5, 1, 2, 3] - call check(error, all(ifftshift(xeven) == [complex(kind=rk) ::(i, i=1, 4)]), & - "all(ifftshift(xeven) == [complex(kind=rk) ::(i, i=1, 4)]) failed.") + call check(error, all(ifftshift(xeven) == [complex(kind=dp) ::(i, i=1, 4)]), & + "all(ifftshift(xeven) == [complex(kind=dp) ::(i, i=1, 4)]) failed.") if (allocated(error)) return - call check(error, all(ifftshift(xodd) == [complex(kind=rk) ::(i, i=1, 5)]), & - "all(ifftshift(xodd) == [complex(kind=rk) ::(i, i=1, 5)]) failed.") + call check(error, all(ifftshift(xodd) == [complex(kind=dp) ::(i, i=1, 5)]), & + "all(ifftshift(xodd) == [complex(kind=dp) ::(i, i=1, 5)]) failed.") end subroutine test_ifftshift_complex @@ -72,14 +72,14 @@ subroutine test_ifftshift_real(error) type(error_type), allocatable, intent(out) :: error integer :: i - real(kind=rk) :: xeven(4) = [3, 4, 1, 2] - real(kind=rk) :: xodd(5) = [4, 5, 1, 2, 3] + real(kind=dp) :: xeven(4) = [3, 4, 1, 2] + real(kind=dp) :: xodd(5) = [4, 5, 1, 2, 3] - call check(error, all(ifftshift(xeven) == [real(kind=rk) ::(i, i=1, 4)]), & - "all(ifftshift(xeven) == [real(kind=rk) ::(i, i=1, 4)]) failed.") + call check(error, all(ifftshift(xeven) == [real(kind=dp) ::(i, i=1, 4)]), & + "all(ifftshift(xeven) == [real(kind=dp) ::(i, i=1, 4)]) failed.") if (allocated(error)) return - call check(error, all(ifftshift(xodd) == [real(kind=rk) ::(i, i=1, 5)]), & - "all(ifftshift(xodd) == [real(kind=rk) ::(i, i=1, 5)]) failed.") + call check(error, all(ifftshift(xodd) == [real(kind=dp) ::(i, i=1, 5)]), & + "all(ifftshift(xodd) == [real(kind=dp) ::(i, i=1, 5)]) failed.") end subroutine test_ifftshift_real @@ -99,14 +99,14 @@ subroutine test_fftfreq_2(error) implicit none type(error_type), allocatable, intent(out) :: error - real(rk), parameter :: tol = 1.0e-12_rk - real(rk), parameter :: twopi = 8*atan(1.0_rk) ! 2*pi - complex(rk), parameter :: imu = (0, 1) ! imaginary unit + real(dp), parameter :: tol = 1.0e-12_dp + real(dp), parameter :: twopi = 8*atan(1.0_dp) ! 2*pi + complex(dp), parameter :: imu = (0, 1) ! imaginary unit integer, parameter :: n = 128 integer :: i - complex(rk), dimension(n) :: xvec, xfou - real(rk), dimension(n) :: xtrue + complex(dp), dimension(n) :: xvec, xfou + real(dp), dimension(n) :: xtrue do i = 1, n xvec(i) = cos(twopi*(i - 1)/n) @@ -124,14 +124,14 @@ subroutine test_fftfreq_3(error) implicit none type(error_type), allocatable, intent(out) :: error - real(rk), parameter :: tol = 1.0e-12_rk - real(rk), parameter :: twopi = 8*atan(1.0_rk) ! 2*pi - complex(rk), parameter :: imu = (0, 1) ! imaginary unit + real(dp), parameter :: tol = 1.0e-12_dp + real(dp), parameter :: twopi = 8*atan(1.0_dp) ! 2*pi + complex(dp), parameter :: imu = (0, 1) ! imaginary unit integer, parameter :: n = 135 integer :: i - complex(rk), dimension(n) :: xvec, xfou - real(rk), dimension(n) :: xtrue + complex(dp), dimension(n) :: xvec, xfou + real(dp), dimension(n) :: xtrue do i = 1, n xvec(i) = cos(twopi*(i - 1)/n)