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,     &
    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, howmany, batch_axis, is_real, sign, flags, din_, cin_, cout_)
    type(c_ptr), intent(inout) :: plan
    integer,     intent(in)    :: dim
    integer,     intent(in)    :: n(:)
    integer(c_int), intent(in) :: howmany
    integer,     intent(in)    :: batch_axis
    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, 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
 
    push_sub(fftw_prepare_plan)
 
    assert(sign == fftw_forward .or. sign == fftw_backward)
    assert(batch_axis == 1 .or. batch_axis == 4)
 
    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
 
    n_r(1:dim) = n(dim:1:-1)
    inembed(1:dim) = n_r(1:dim)
    onembed(1:dim) = n_r(1:dim)
    ! batch-first (batch_axis == 1): the batch is the fastest-varying axis, so consecutive
    !   transforms are interleaved one element apart (idist = 1) and the elements within a
    !   transform are howmany apart (istride = howmany).
    ! batch-last (batch_axis == 4): 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.
    if (batch_axis == 1) then
      istride = howmany
      ostride = howmany
      idist = 1_c_int
      odist = 1_c_int
    else if (batch_axis == 4) then
      istride = 1_c_int
      ostride = 1_c_int
      idist = 1_c_int   ! set below
      odist = 1_c_int   ! set below
    else
      assert(.false.)
    end if
 
    if (is_real) then
      if (sign == fftw_forward) then
        if (.not. aligned_memory) then
          safe_allocate(rin(1:howmany, 1:n(1), 1:n(2), 1:n(3)))
          safe_allocate(cout(1:howmany, 1:n(1)/2+1, 1:n(2), 1:n(3)))
        else
          rin => din_
          cout => cout_
        end if
 
        onembed(dim) = onembed(dim)/2 + 1
        if (batch_axis == 4) then
          idist = product(inembed(1:dim))
          odist = product(onembed(1:dim))
        end if
        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:howmany, 1:n(1)/2+1, 1:n(2), 1:n(3)))
          safe_allocate(rout(1:howmany, 1:n(1), 1:n(2), 1:n(3)))
        else
          cin => cout_
          rout => din_
        end if
 
        inembed(dim) = inembed(dim)/2 + 1
        if (batch_axis == 4) then
          idist = product(inembed(1:dim))
          odist = product(onembed(1:dim))
        end if
        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:howmany, 1:n(1), 1:n(2), 1:n(3)))
        safe_allocate(cout(1:howmany, 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
 
      if (batch_axis == 4) then
        idist = product(inembed(1:dim))
        odist = product(onembed(1:dim))
      end if
      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
 
    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_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
 
end module fftw_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: