diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/f1.f90 | 81 | ||||
| -rw-r--r-- | src/fourier.f90 | 116 |
2 files changed, 116 insertions, 81 deletions
diff --git a/src/f1.f90 b/src/f1.f90 deleted file mode 100644 index 7f07afa..0000000 --- a/src/f1.f90 +++ /dev/null @@ -1,81 +0,0 @@ -PROGRAM test - - use, intrinsic :: iso_c_binding - IMPLICIT none - - INCLUDE 'fftw3.f03' - - real(kind=8), parameter :: pi = acos(-1.0) ! also used from f1() - - ! rest of variables only used from "main" program code - integer, parameter :: npoints = 65536, ncomplex = npoints / 2 + 1 - real(kind=8), parameter :: start = 0, end = 1 - integer :: i - real(kind=8) :: t, & - delta = (end - start) / real(npoints - 1, C_DOUBLE) - - type(C_PTR) :: p_real, p_complex, plan - - real(C_DOUBLE), pointer :: arr_real(:) - complex(C_DOUBLE_COMPLEX), pointer :: arr_complex(:) - - character(100) :: arg - - IF (command_argument_count() < 1) THEN - - DO i = 1, npoints - - t = (i - 1) * delta + start - - write(*,*) t, " ", f1(t) - - END DO - - STOP - END IF - - call get_command_argument(1, arg) - - IF (trim(arg) /= "dft") THEN - - STOP "wrong argument" - END IF - - p_real = fftw_alloc_real(int(npoints, C_SIZE_T)) - - p_complex = fftw_alloc_complex(int(ncomplex, C_SIZE_T)) - - call c_f_pointer(p_real, arr_real, [npoints]) - call c_f_pointer(p_complex, arr_complex, [ncomplex]) - - plan = fftw_plan_dft_r2c_1d(int(npoints, C_INT), arr_real, & - arr_complex, FFTW_ESTIMATE) - - DO i = 1, npoints - - t = (i - 1) * delta + start - arr_real(i) = real(f1(t), C_DOUBLE) - - END DO - - call fftw_execute_dft_r2c(plan, arr_real, arr_complex) - - DO i = 1, ncomplex - - write(*,*) i - 1, " ", abs(arr_complex(i)) - - END DO - - call fftw_free(p_real) - call fftw_free(p_complex) - -CONTAINS - - PURE real(kind=8) FUNCTION f1(t) - real(kind=8), intent(in) :: t - - f1 = sin(2 * pi * t * 200) + 2 * sin(2 * pi * t * 400) - - END FUNCTION f1 - -END PROGRAM test diff --git a/src/fourier.f90 b/src/fourier.f90 new file mode 100644 index 0000000..be0a16f --- /dev/null +++ b/src/fourier.f90 @@ -0,0 +1,116 @@ +PROGRAM test + + use, intrinsic :: iso_c_binding + IMPLICIT none + + INCLUDE 'fftw3.f03' + + real(kind=8), parameter :: pi = acos(-1.0) ! also used from f1() + + ! rest of variables only used from "main" program code + integer, parameter :: npoints = 65536, ncomplex = npoints / 2 + 1 + real(kind=8), parameter :: start = 0, end = 1 + integer :: i + real(kind=8) :: t, & + delta = (end - start) / real(npoints - 1, C_DOUBLE) + + real(kind=8) :: val, randvals(npoints) + + type(C_PTR) :: p_real, p_complex, plan_f, plan_b ! forward,backward + + real(C_DOUBLE), pointer :: arr_real(:) + complex(C_DOUBLE_COMPLEX), pointer :: arr_complex(:) + + character(100) :: arg + + logical :: which_fun + + IF (command_argument_count() < 1) STOP "argument required" + + call get_command_argument(1, arg) + + IF (arg == "f1") THEN + which_fun = .true. + ELSE IF (arg == "f2") THEN + which_fun = .true. + call random_seed() + call random_number(randvals) + ELSE + STOP "bad first argument (should be 'f1' or 'f2')" + END IF + + IF (command_argument_count() < 2) THEN + DO i = 1, npoints + + t = (i - 1) * delta + start + + IF (which_fun) val = f1(t) ! f1 + + IF (.not. which_fun) THEN ! f2 + val = f2(t) + randvals(i) + END IF + + write(*,*) t, " ", val + + END DO + STOP + END IF + + call get_command_argument(2, arg) + + IF ((trim(arg) /= "dft") .and. (trim(arg) /= "transformed")) & + STOP "bad second argument (should be 'dft' or 'transformed')" + + p_real = fftw_alloc_real(int(npoints, C_SIZE_T)) + + p_complex = fftw_alloc_complex(int(ncomplex, C_SIZE_T)) + + call c_f_pointer(p_real, arr_real, [npoints]) + call c_f_pointer(p_complex, arr_complex, [ncomplex]) + + plan_f = fftw_plan_dft_r2c_1d(int(npoints, C_INT), arr_real, & + arr_complex, FFTW_ESTIMATE) + plan_b = fftw_plan_dft_c2r_1d(int(npoints, C_INT), arr_complex,& + arr_real, FFTW_ESTIMATE) + + IF (which_fun) THEN ! f1 + arr_real = real(f1((/(i, i = 1, npoints, 1)/) * delta + start),& + C_DOUBLE) + ELSE ! f2 + arr_real = real(f2((/(i, i = 1, npoints, 1)/) * delta + start & + + randvals), C_DOUBLE) + END IF + + call fftw_execute_dft_r2c(plan_f, arr_real, arr_complex) + + IF (trim(arg) == "dft") THEN + DO i = 1, ncomplex + + write(*,*) i - 1, " ", abs(arr_complex(i)) + + END DO + ELSE + STOP "not implemented" + END IF + + call fftw_free(p_real) + call fftw_free(p_complex) + +CONTAINS + + PURE ELEMENTAL real(kind=8) FUNCTION f1(t) + real(kind=8), intent(in) :: t + + f1 = sin(2 * pi * t * 200) + 2 * sin(2 * pi * t * 400) + + END FUNCTION f1 + + + PURE ELEMENTAL real(kind=8) FUNCTION f2(t) + real(kind=8), intent(in) :: t + + f2 = cos(2 * pi * t) + + END FUNCTION f2 + +END PROGRAM test |