fftw.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2011 M. Oliveira
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module fftw_params_oct_m
  use, intrinsic :: iso_c_binding
  implicit none
 
  private
 
  ! make public some needed symbols from the FFTW library
  public ::               &
    fftw_cleanup,         &
    fftw_execute_dft,     &
    fftw_execute_dft_c2r, &
    fftw_execute_dft_r2c, &
    fftw_execute_r2r,     &
    fftw_destroy_plan,    &
    fftw_alloc_real,      &
    fftw_alloc_complex,   &
    fftw_free,            &
    fftw_plan_dft_1d,     &
    fftw_plan_dft_2d,     &
    fftw_plan_dft_3d,     &
    fftw_plan_dft_r2c_1d, &
    fftw_plan_dft_r2c_2d, &
    fftw_plan_dft_r2c_3d, &
    fftw_plan_dft_c2r_1d, &
    fftw_plan_dft_c2r_2d, &
    fftw_plan_dft_c2r_3d, &
    fftw_plan_r2r_1d,     &
    c_fftw_r2r_kind,      &
    fftw_r2hc,            &
    fftw_hc2r,            &
    fftw_dht,             &
    fftw_redft00,         &
    fftw_redft01,         &
    fftw_redft10,         &
    fftw_redft11,         &
    fftw_rodft00,         &
    fftw_rodft01,         &
    fftw_rodft10,         &
    fftw_rodft11,         &
    fftw_forward,         &
    fftw_backward,        &
    fftw_measure,         &
    fftw_estimate,        &
    fftw_destroy_input,   &
    fftw_unaligned
#ifdef HAVE_FFTW3_MPI
  public :: fftw_mpi_default_block
#endif
#if defined(HAVE_OPENMP) && defined(HAVE_FFTW3_THREADS)
  public ::                  &
    fftw_init_threads,       &
    fftw_plan_with_nthreads, &
    fftw_cleanup_threads
#endif
 
 
#ifdef HAVE_FFTW3_MPI
  include "fftw3-mpi.f03"
#else
  include "fftw3.f03"
#endif
end module fftw_params_oct_m
 
module fftw_oct_m
  use fftw_params_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_c_binding
  use messages_oct_m
  use profiling_oct_m
  implicit none
 
  private
 
  public :: &
    fftw_prepare_plan,       &
    fftw_get_dims,           &
    fftw_alloc_memory,       &
    fftw_free_memory
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine fftw_prepare_plan(plan, dim, n, is_real, sign, flags, din_, cin_, cout_)
    type(c_ptr), intent(inout) :: plan
    integer,     intent(in)    :: dim
    integer,     intent(in)    :: n(:)
    logical,     intent(in)    :: is_real
    integer,     intent(in)    :: sign
    integer,     intent(in)    :: flags
    real(real64), optional, target, intent(in) :: din_(:,:,:)
    complex(real64), optional, target, intent(in) :: cin_(:,:,:)
    complex(real64), optional, target, intent(in) :: cout_(:,:,:)
 
    real(real64), pointer :: rin(:,:,:)
    real(real64), pointer :: rout(:,:,:)
    complex(real64), pointer :: cin(:,:,:)
    complex(real64), pointer :: cout(:,:,:)
    logical :: aligned_memory
 
    push_sub(fftw_prepare_plan)
 
    assert(sign == fftw_forward .or. sign == fftw_backward)
 
    aligned_memory = .false.
    if (present(din_) .or. present(cin_)) then
      assert(present(cout_))
      assert(present(din_) .neqv. present(cin_))
      aligned_memory = .true.
      if (is_real) then
        assert(present(din_))
      else
        assert(present(cin_))
      end if
    end if
 
    if (is_real) then
      if (sign == fftw_forward) then
        if (.not. aligned_memory) then
          safe_allocate(rin(1:n(1), 1:n(2), 1:n(3)))
          safe_allocate(cout(1:n(1)/2+1, 1:n(2), 1:n(3)))
        else
          rin => din_
          cout => cout_
        end if
 
        select case (dim)
        case (1)
          plan = fftw_plan_dft_r2c_1d(n(1), rin, cout, flags)
        case (2)
          plan = fftw_plan_dft_r2c_2d(n(2), n(1),rin, cout, flags)
        case (3)
          plan = fftw_plan_dft_r2c_3d(n(3), n(2), n(1), rin, cout, flags)
        end select
 
        if (.not. aligned_memory) then
          safe_deallocate_p(rin)
          safe_deallocate_p(cout)
        else
          nullify(rin, cout)
        end if
      else
        if (.not. aligned_memory) then
          safe_allocate(cin(1:n(1)/2+1, 1:n(2), 1:n(3)))
          safe_allocate(rout(1:n(1), 1:n(2), 1:n(3)))
        else
          cin => cout_
          rout => din_
        end if
 
        select case (dim)
        case (1)
          plan = fftw_plan_dft_c2r_1d(n(1), cin, rout, flags)
        case (2)
          plan = fftw_plan_dft_c2r_2d(n(2), n(1), cin, rout, flags)
        case (3)
          plan = fftw_plan_dft_c2r_3d(n(3), n(2), n(1), cin, rout, flags)
        end select
 
        if (.not. aligned_memory) then
          safe_deallocate_p(cin)
          safe_deallocate_p(rout)
        else
          nullify(cin,rout)
        end if
      end if
    else
      if (.not. aligned_memory) then
        safe_allocate(cin(1:n(1), 1:n(2), 1:n(3)))
        safe_allocate(cout(1:n(1), 1:n(2), 1:n(3)))
      else
        if (sign == fftw_forward) then
          cin => cin_
          cout => cout_
        else
          cout => cin_
          cin => cout_
        end if
      end if
 
      select case (dim)
      case (1)
        plan = fftw_plan_dft_1d(n(1), cin, cout, sign, flags)
      case (2)
        plan = fftw_plan_dft_2d(n(2), n(1), cin, cout, sign, flags)
      case (3)
        plan = fftw_plan_dft_3d(n(3), n(2), n(1), cin, cout, sign, flags)
      end select
 
      if (.not. aligned_memory) then
        safe_deallocate_p(cin)
        safe_deallocate_p(cout)
      else
        nullify(cin, cout)
      end if
    end if
 
    pop_sub(fftw_prepare_plan)
  end subroutine fftw_prepare_plan
 
  ! ---------------------------------------------------------
  subroutine fftw_get_dims(rs_n, is_real, fs_n)
    integer, intent(in)  :: rs_n(:)
    logical, intent(in)  :: is_real
    integer, intent(out) :: fs_n(:)
 
    push_sub(fftw_get_dims)
 
    fs_n = rs_n
    if (is_real) fs_n(1) = rs_n(1)/2 + 1
 
    pop_sub(fftw_get_dims)
  end subroutine fftw_get_dims
 
  ! ---------------------------------------------------------
  subroutine fftw_alloc_memory(rs_n, is_real, fs_n, drs_data, zrs_data, fs_data)
    integer, intent(in)    :: rs_n(:)
    logical, intent(in)    :: is_real
    integer, intent(in)    :: fs_n(:)
    real(real64), pointer, intent(inout) :: drs_data(:,:,:)
    complex(real64), pointer, intent(inout) :: zrs_data(:,:,:)
    complex(real64), pointer, intent(inout) ::  fs_data(:,:,:)
    type(c_ptr) :: address
 
    push_sub(fftw_alloc_memory)
 
    if (is_real) then
      address  = fftw_alloc_real(int(product(rs_n(1:3)), c_size_t))
      call c_f_pointer(address, drs_data, rs_n)
    else
      address = fftw_alloc_complex(int(product(rs_n(1:3)), c_size_t))
      call c_f_pointer(address, zrs_data, rs_n)
    end if
    address = fftw_alloc_complex(int(product(fs_n(1:3)), c_size_t))
    call c_f_pointer(address, fs_data, fs_n)
 
    pop_sub(fftw_alloc_memory)
  end subroutine fftw_alloc_memory
 
  ! ---------------------------------------------------------
  subroutine fftw_free_memory(is_real, drs_data, zrs_data, fs_data)
    logical, intent(in)    :: is_real
    real(real64), pointer, intent(inout) :: drs_data(:,:,:)
    complex(real64), pointer, intent(inout) :: zrs_data(:,:,:)
    complex(real64), pointer, intent(inout) ::  fs_data(:,:,:)
 
    push_sub(fftw_free_memory)
 
    if (is_real) then
      call fftw_free(c_loc(drs_data))
      nullify(drs_data)
    else
      call fftw_free(c_loc(zrs_data))
      nullify(zrs_data)
    end if
    call fftw_free(c_loc(fs_data))
    nullify(fs_data)
 
    pop_sub(fftw_free_memory)
  end subroutine fftw_free_memory
 
end module fftw_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: