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_many_dft,   &
    fftw_plan_many_dft_r2c, &
    fftw_plan_many_dft_c2r, &
    fftw_plan_r2r_1d,     &
    fftwf_cleanup,         &
    fftwf_execute_dft,     &
    fftwf_execute_dft_c2r, &
    fftwf_execute_dft_r2c, &
    fftwf_destroy_plan,    &
    fftwf_alloc_real,      &
    fftwf_alloc_complex,   &
    fftwf_free,            &
    fftwf_plan_many_dft,   &
    fftwf_plan_many_dft_r2c, &
    fftwf_plan_many_dft_c2r, &
    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
#ifndef NDEBUG
  use, intrinsic :: ieee_exceptions
#endif
  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,         &
    fftw_prepare_plan_single, &
    fftw_alloc_memory_single, &
    fftw_free_memory_single
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine fftw_prepare_plan(plan, dim, n, howmany, is_real, sign, flags_, din_, cin_, cout_)
    type(c_ptr), intent(inout) :: plan
    integer, intent(in) :: dim
    integer, contiguous, intent(in) :: n(:)
    integer(c_int), intent(in) :: howmany
    logical, intent(in) :: is_real
    integer, intent(in) :: sign
    integer, intent(in) :: flags_
    real(real64), contiguous, optional, target, intent(in) :: din_(:,:,:,:)
    complex(real64), contiguous, optional, target, intent(in) :: cin_(:,:,:,:)
    complex(real64), contiguous, optional, target, intent(in) :: cout_(:,:,:,:)
 
    real(real64), pointer, contiguous :: rin(:,:,:,:)
    real(real64), pointer, contiguous :: rout(:,:,:,:)
    complex(real64), pointer, contiguous :: cin(:,:,:,:)
    complex(real64), pointer, contiguous :: cout(:,:,:,:)
    logical :: aligned_memory
    integer(c_int) :: n_r(1:dim)
    integer(c_int) :: inembed(1:dim), onembed(1:dim)
    integer(c_int) :: istride, ostride, idist, odist, flags
#ifndef NDEBUG
    logical :: halting_mode(size(ieee_all))
#endif
 
    push_sub(fftw_prepare_plan)
 
    assert(sign == fftw_forward .or. sign == fftw_backward)
 
#ifndef NDEBUG
    ! The oneMKL FFTW wrappers raise spurious floating-point exceptions while
    ! committing a DFTI descriptor, so suspend trapping for builds that enable it
    call ieee_get_halting_mode(ieee_all, halting_mode)
    call ieee_set_halting_mode(ieee_all, .false.)
#endif
 
    aligned_memory = .false.
    ! at the calling site, one of these is always passed, but present is only true
    ! if the corresponding pointer is associated
    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
      flags = flags_
    else
      flags = flags_ + fftw_unaligned
    end if
 
    n_r(1:dim) = n(dim:1:-1)
    inembed(1:dim) = n_r(1:dim)
    onembed(1:dim) = n_r(1:dim)
    ! batch-last: each transform is contiguous (istride = 1) and the transforms are a whole grid
    ! apart (idist = product of the per-transform extents). Because the real and complex grids
    ! differ in size for r2c/c2r, idist/odist are finalized per branch below, once inembed/onembed
    ! have been halved.
    istride = 1_c_int
    ostride = 1_c_int
    idist = 1_c_int   ! set below
    odist = 1_c_int   ! set below
 
    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), 1:howmany))
          safe_allocate(cout(1:n(1)/2+1, 1:n(2), 1:n(3), 1:howmany))
        else
          rin => din_
          cout => cout_
        end if
 
        onembed(dim) = onembed(dim)/2 + 1
        idist = product(inembed(1:dim))
        odist = product(onembed(1:dim))
        plan = fftw_plan_many_dft_r2c(dim, n_r, howmany, &
          rin, inembed, istride, idist, &
          cout, onembed, ostride, odist, &
          flags)
        assert(c_associated(plan))
 
        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), 1:howmany))
          safe_allocate(rout(1:n(1), 1:n(2), 1:n(3), 1:howmany))
        else
          cin => cout_
          rout => din_
        end if
 
        inembed(dim) = inembed(dim)/2 + 1
        idist = product(inembed(1:dim))
        odist = product(onembed(1:dim))
        plan = fftw_plan_many_dft_c2r(dim, n_r, howmany, &
          cin, inembed, istride, idist, &
          rout, onembed, ostride, odist, &
          flags)
        assert(c_associated(plan))
 
        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), 1:howmany))
        safe_allocate(cout(1:n(1), 1:n(2), 1:n(3), 1:howmany))
      else
        if (sign == fftw_forward) then
          cin => cin_
          cout => cout_
        else
          cout => cin_
          cin => cout_
        end if
      end if
 
      idist = product(inembed(1:dim))
      odist = product(onembed(1:dim))
      plan = fftw_plan_many_dft(dim, n_r, howmany, &
        cin, inembed, istride, idist, &
        cout, onembed, ostride, odist, &
        sign, flags)
      assert(c_associated(plan))
 
      if (.not. aligned_memory) then
        safe_deallocate_p(cin)
        safe_deallocate_p(cout)
      else
        nullify(cin, cout)
      end if
    end if
 
#ifndef NDEBUG
    call ieee_set_halting_mode(ieee_all, halting_mode)
#endif
 
    pop_sub(fftw_prepare_plan)
  end subroutine fftw_prepare_plan
 
  ! ---------------------------------------------------------
  subroutine fftw_get_dims(rs_n, is_real, fs_n)
    integer, contiguous, intent(in) :: rs_n(:)
    logical, intent(in) :: is_real
    integer, contiguous, 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_dims, is_real, fs_dims, drs_data, zrs_data, fs_data)
    integer, intent(in)    :: rs_dims(4)
    logical, intent(in)    :: is_real
    integer, intent(in)    :: fs_dims(4)
    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_dims), c_size_t))
      call c_f_pointer(address, drs_data, rs_dims)
    else
      address = fftw_alloc_complex(int(product(rs_dims), c_size_t))
      call c_f_pointer(address, zrs_data, rs_dims)
    end if
 
    address = fftw_alloc_complex(int(product(fs_dims), c_size_t))
    call c_f_pointer(address, fs_data, fs_dims)
 
    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
      if (associated(drs_data)) then
        call fftw_free(c_loc(drs_data))
        nullify(drs_data)
      end if
    else
      if (associated(zrs_data)) then
        call fftw_free(c_loc(zrs_data))
        nullify(zrs_data)
      end if
    end if
    if (associated(fs_data)) then
      call fftw_free(c_loc(fs_data))
      nullify(fs_data)
    end if
 
    pop_sub(fftw_free_memory)
  end subroutine fftw_free_memory
 
  ! ---------------------------------------------------------
  subroutine fftw_prepare_plan_single(plan, dim, n, howmany, is_real, sign, flags, sin_, cin_, cout_)
    type(c_ptr), intent(inout) :: plan
    integer,     intent(in)    :: dim
    integer, contiguous, intent(in) :: n(:)
    integer(c_int), intent(in) :: howmany
    logical,     intent(in)    :: is_real
    integer,     intent(in)    :: sign
    integer,     intent(in)    :: flags
    real(real32), contiguous, optional, target, intent(in) :: sin_(:,:,:,:)
    complex(real32), contiguous, optional, target, intent(in) :: cin_(:,:,:,:)
    complex(real32), contiguous, optional, target, intent(in) :: cout_(:,:,:,:)
 
    real(real32),    pointer, contiguous :: rin(:,:,:,:), rout(:,:,:,:)
    complex(real32), pointer, contiguous :: cin(:,:,:,:), cout(:,:,:,:)
    integer(c_int) :: n_r(1:dim)
    integer(c_int) :: inembed(1:dim), onembed(1:dim)
    integer(c_int) :: istride, ostride, idist, odist
#ifndef NDEBUG
    logical :: halting_mode(size(ieee_all))
#endif
 
    push_sub(fftw_prepare_plan_single)
 
    assert(sign == fftw_forward .or. sign == fftw_backward)
 
#ifndef NDEBUG
    ! See fftw_prepare_plan
    call ieee_get_halting_mode(ieee_all, halting_mode)
    call ieee_set_halting_mode(ieee_all, .false.)
#endif
 
    n_r(1:dim) = n(dim:1:-1)
    inembed(1:dim) = n_r(1:dim)
    onembed(1:dim) = n_r(1:dim)
    istride = 1_c_int
    ostride = 1_c_int
 
    if (is_real) then
      if (sign == fftw_forward) then
        rin => sin_
        cout => cout_
        onembed(dim) = onembed(dim)/2 + 1
        idist = product(inembed(1:dim))
        odist = product(onembed(1:dim))
        plan = fftwf_plan_many_dft_r2c(dim, n_r, howmany, &
          rin, inembed, istride, idist, cout, onembed, ostride, odist, flags)
      else
        cin => cout_
        rout => sin_
        inembed(dim) = inembed(dim)/2 + 1
        idist = product(inembed(1:dim))
        odist = product(onembed(1:dim))
        plan = fftwf_plan_many_dft_c2r(dim, n_r, howmany, &
          cin, inembed, istride, idist, rout, onembed, ostride, odist, flags)
      end if
    else
      if (sign == fftw_forward) then
        cin => cin_
        cout => cout_
      else
        cout => cin_
        cin => cout_
      end if
      idist = product(inembed(1:dim))
      odist = product(onembed(1:dim))
      plan = fftwf_plan_many_dft(dim, n_r, howmany, &
        cin, inembed, istride, idist, cout, onembed, ostride, odist, sign, flags)
    end if
    assert(c_associated(plan))
 
#ifndef NDEBUG
    call ieee_set_halting_mode(ieee_all, halting_mode)
#endif
 
    pop_sub(fftw_prepare_plan_single)
  end subroutine fftw_prepare_plan_single
 
  ! ---------------------------------------------------------
  subroutine fftw_alloc_memory_single(rs_dims, is_real, fs_dims, srs_data, crs_data, cfs_data)
    integer, intent(in)    :: rs_dims(4)
    logical, intent(in)    :: is_real
    integer, intent(in)    :: fs_dims(4)
    real(real32),    pointer, intent(inout) :: srs_data(:,:,:,:)
    complex(real32), pointer, intent(inout) :: crs_data(:,:,:,:)
    complex(real32), pointer, intent(inout) :: cfs_data(:,:,:,:)
 
    type(c_ptr) :: address
 
    push_sub(fftw_alloc_memory_single)
 
    if (is_real) then
      address = fftwf_alloc_real(int(product(rs_dims), c_size_t))
      call c_f_pointer(address, srs_data, rs_dims)
    else
      address = fftwf_alloc_complex(int(product(rs_dims), c_size_t))
      call c_f_pointer(address, crs_data, rs_dims)
    end if
 
    address = fftwf_alloc_complex(int(product(fs_dims), c_size_t))
    call c_f_pointer(address, cfs_data, fs_dims)
 
    pop_sub(fftw_alloc_memory_single)
  end subroutine fftw_alloc_memory_single
 
  ! ---------------------------------------------------------
  subroutine fftw_free_memory_single(is_real, srs_data, crs_data, cfs_data)
    logical, intent(in)    :: is_real
    real(real32),    pointer, intent(inout) :: srs_data(:,:,:,:)
    complex(real32), pointer, intent(inout) :: crs_data(:,:,:,:)
    complex(real32), pointer, intent(inout) :: cfs_data(:,:,:,:)
 
    push_sub(fftw_free_memory_single)
 
    if (is_real) then
      if (associated(srs_data)) then
        call fftwf_free(c_loc(srs_data))
        nullify(srs_data)
      end if
    else
      if (associated(crs_data)) then
        call fftwf_free(c_loc(crs_data))
        nullify(crs_data)
      end if
    end if
    if (associated(cfs_data)) then
      call fftwf_free(c_loc(cfs_data))
      nullify(cfs_data)
    end if
 
    pop_sub(fftw_free_memory_single)
  end subroutine fftw_free_memory_single
 
end module fftw_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: