pfft.F90

Go to the documentation of this file.

!! Copyright (C) 2011 J. Alberdi-Rodriguez, P. Garcia Risueño, 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 pfft_params_oct_m
  use, intrinsic :: iso_c_binding
  use fftw_params_oct_m
  implicit none
 
  private
 
#ifdef HAVE_PFFT
 
  ! make public some symbols from the pfft library
  public ::                     &
    pfft_cleanup,               &
    pfft_create_procmesh_2d,    &
    pfft_destroy_plan,          &
    pfft_execute,               &
    pfft_init,                  &
    pfft_local_size_many_dft_r2c, &
    pfft_local_size_many_dft,   &
    pfft_plan_many_dft_c2r,     &
    pfft_plan_many_dft_r2c,     &
    pfft_plan_many_dft,         &
    pfft_int,                   &
    pfft_ptrdiff_t,             &
    pfft_float,                 &
    pfft_double,                &
    pfft_ldouble,               &
    pfft_unsigned,              &
    pfft_transposed_in,         &
    pfft_transposed_out,        &
    pfft_estimate,              &
    pfft_tune,                  &
    pfft_destroy_input,         &
    pfft_measure,               &
    pfft_forward,               &
    pfft_backward
 
#ifndef HAVE_FFTW3_MPI
  integer(C_INTPTR_T), parameter :: FFTW_MPI_DEFAULT_BLOCK = 0
#endif
  include "pfft.f03"
#endif
 
end module pfft_params_oct_m
 
module pfft_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_c_binding
  use math_oct_m
  use pfft_params_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
 
  public ::                    &
    pfft_decompose,            &
    pfft_prepare_plan_r2c,     &
    pfft_prepare_plan_c2r,     &
    pfft_prepare_plan_c2c,     &
    pfft_get_dims
 
contains
 
  ! ---------------------------------------------------------
  subroutine pfft_decompose(n_proc, dim1, dim2)
    integer, intent(in)  :: n_proc
    integer, intent(out) :: dim1
    integer, intent(out) :: dim2
 
    integer :: np, i
 
    push_sub(pfft_decompose)
 
    assert(n_proc > 0)
 
    dim1 = 1
    dim2 = 1
    np = n_proc
    i = n_proc-1
 
    if (is_prime(n_proc)) then
      dim1 = n_proc
    else
      do
        if (i <= 1) exit
        if (mod(np,i) /= 0.or.(.not. is_prime(i))) then
          i=i-1
          cycle
        end if
        np=np/i
        if (dim1 <= dim2) then
          dim1 = dim1*i
        else
          dim2 = dim2*i
        end if
      end do
    end if
 
    assert(dim1*dim2 == n_proc)
 
    pop_sub(pfft_decompose)
  end subroutine pfft_decompose
 
  ! ---------------------------------------------------------
  subroutine pfft_prepare_plan_r2c(plan, n, in, out, fft_sign, flags, mpi_comm, howmany)
    type(C_PTR),                          intent(out)   :: plan
    integer,                              intent(in)    :: n(:)
    real(real64),    contiguous, pointer, intent(inout) :: in(:,:,:,:)
    complex(real64), contiguous, pointer, intent(inout) :: out(:,:,:,:)
    integer,                              intent(in)    :: fft_sign
    integer,                              intent(in)    :: flags
    integer,                              intent(in)    :: mpi_comm
    integer,                              intent(in)    :: howmany
 
    integer(C_INTPTR_T) :: tmp_n(3)
#ifdef HAVE_PFFT
    integer(C_INTPTR_T) :: default_blocks(3)
#endif
 
    push_sub(pfft_prepare_plan_r2c)
    call profiling_in("PFFT_PLAN_R2C")
 
#ifdef HAVE_PFFT
    assert(fft_sign == pfft_forward)
#endif
 
    tmp_n(1:3) = n(3:1:-1)
 
#ifdef HAVE_PFFT
    ! A per-dimension block of FFTW_MPI_DEFAULT_BLOCK (0) requests the PFFT default block
    ! distribution (the C NULL/PFFT_DEFAULT_BLOCKS pointer cannot be passed through the
    ! dimension(*) Fortran interface).
    default_blocks(1:3) = 0_c_intptr_t
    ! the three tmp_n arguments are the transform size n == n_in (embed) == n_out (embed)
    plan = pfft_plan_many_dft_r2c(3, tmp_n, tmp_n, tmp_n, int(howmany, c_intptr_t), &
      default_blocks, default_blocks, in, out, mpi_comm, fft_sign, &
      pfft_transposed_out + pfft_measure + pfft_destroy_input + pfft_tune)
#else
    plan = c_null_ptr
#endif
 
    call profiling_out("PFFT_PLAN_R2C")
    pop_sub(pfft_prepare_plan_r2c)
  end subroutine pfft_prepare_plan_r2c
 
  ! ---------------------------------------------------------
  subroutine pfft_prepare_plan_c2r(plan, n, in, out, fft_sign, flags, mpi_comm, howmany)
    type(c_ptr),                          intent(out)   :: plan
    integer,                              intent(in)    :: n(:)
    complex(real64), contiguous, pointer, intent(inout) :: in(:,:,:,:)
    real(real64),    contiguous, pointer, intent(inout) :: out(:,:,:,:)
    integer,                              intent(in)    :: fft_sign
    integer,                              intent(in)    :: flags
    integer,                              intent(in)    :: mpi_comm
    integer,                              intent(in)    :: howmany
 
    integer(C_INTPTR_T) :: tmp_n(3)
#ifdef HAVE_PFFT
    integer(C_INTPTR_T) :: default_blocks(3)
#endif
    push_sub(pfft_prepare_plan_c2r)
    call profiling_in("PFFT_PLAN_C2R")
 
#ifdef HAVE_PFFT
    assert(fft_sign == pfft_backward)
#endif
 
    tmp_n(1:3) = n(3:1:-1)
 
#ifdef HAVE_PFFT
    default_blocks(1:3) = 0_c_intptr_t
    ! the three tmp_n arguments are the transform size n == n_in (embed) == n_out (embed)
    plan = pfft_plan_many_dft_c2r(3, tmp_n, tmp_n, tmp_n, int(howmany, c_intptr_t), &
      default_blocks, default_blocks, in, out, mpi_comm, fft_sign, &
      pfft_transposed_in + pfft_measure + pfft_destroy_input + pfft_tune)
#else
    plan = c_null_ptr
#endif
 
    call profiling_out("PFFT_PLAN_C2R")
    pop_sub(pfft_prepare_plan_c2r)
  end subroutine pfft_prepare_plan_c2r
 
  ! ---------------------------------------------------------
  subroutine pfft_prepare_plan_c2c(plan, n, in, out, fft_sign, flags, mpi_comm, howmany)
    type(c_ptr),                          intent(out)   :: plan
    integer,                              intent(in)    :: n(:)
    complex(real64), contiguous, pointer, intent(inout) :: in(:,:,:,:)
    complex(real64), contiguous, pointer, intent(inout) :: out(:,:,:,:)
    integer,                              intent(in)    :: fft_sign
    integer,                              intent(in)    :: flags
    integer,                              intent(in)    :: mpi_comm
    integer,                              intent(in)    :: howmany
 
#ifdef HAVE_PFFT
    integer(C_INT) :: pfft_flags
    integer(C_INTPTR_T) :: default_blocks(3)
#endif
    integer(C_INTPTR_T) :: tmp_n(3)
 
    push_sub(pfft_prepare_plan_c2c)
 
    tmp_n(1:3) = n(3:1:-1)
 
#ifdef HAVE_PFFT
    if (fft_sign == pfft_forward) then
      pfft_flags = pfft_transposed_out
    else
      pfft_flags = pfft_transposed_in
    end if
 
    default_blocks(1:3) = 0_c_intptr_t
    ! the three tmp_n arguments are the transform size n == n_in (embed) == n_out (embed)
    plan = pfft_plan_many_dft(3, tmp_n, tmp_n, tmp_n, int(howmany, c_intptr_t), &
      default_blocks, default_blocks, in, out, mpi_comm, fft_sign, pfft_flags)
#else
    plan = c_null_ptr
#endif
 
    pop_sub(pfft_prepare_plan_c2c)
  end subroutine pfft_prepare_plan_c2c
 
  ! ---------------------------------------------------------
  subroutine pfft_get_dims(rs_n_global, mpi_comm, is_real, alloc_size, fs_n_global, rs_n, fs_n, rs_istart, fs_istart, howmany)
    integer,     intent(in)  :: rs_n_global(1:3)
    integer,     intent(in)  :: mpi_comm
    logical,     intent(in)  :: is_real
    integer(int64), intent(out) :: alloc_size
    integer,     intent(out) :: fs_n_global(1:3)
    integer,     intent(out) :: rs_n(1:3)
    integer,     intent(out) :: fs_n(1:3)
    integer,     intent(out) :: rs_istart(1:3)
    integer,     intent(out) :: fs_istart(1:3)
    integer,     intent(in)  :: howmany
 
    integer(C_INTPTR_T) :: tmp_alloc_size, tmp_n(3)
    integer(C_INTPTR_T) :: tmp_rs_n(3), tmp_rs_istart(3)
    integer(C_INTPTR_T) :: tmp_fs_n(3), tmp_fs_istart(3)
#ifdef HAVE_PFFT
    integer(C_INTPTR_T) :: default_blocks(3)
#endif
 
    push_sub(pfft_get_dims)
 
    fs_n_global(1:3) = rs_n_global(1:3)
    ! if calling the ISO_C_BINDING routines from PFFT, one must take into account that
    ! the array order is the C array order
    tmp_n(1:3) = rs_n_global(3:1:-1)
 
    tmp_alloc_size = 0
#ifdef HAVE_PFFT
    ! A per-dimension block of FFTW_MPI_DEFAULT_BLOCK (0) requests the PFFT default block
    ! distribution (equivalent to passing the C NULL/PFFT_DEFAULT_BLOCKS pointer, which cannot be
    ! expressed through the dimension(*) Fortran interface). The local grid dims
    ! (tmp_rs_n/tmp_fs_n) are independent of howmany; tmp_alloc_size is the total local element
    ! count including the howmany batch.
    default_blocks(1:3) = 0_c_intptr_t
    if (is_real) then
      tmp_alloc_size = pfft_local_size_many_dft_r2c(3, tmp_n, tmp_n, tmp_n, int(howmany, c_intptr_t), &
        default_blocks, default_blocks, mpi_comm, pfft_transposed_out, &
        tmp_rs_n, tmp_rs_istart, tmp_fs_n, tmp_fs_istart)
      fs_n_global(1) = rs_n_global(1)/2 + 1
    else
      tmp_alloc_size = pfft_local_size_many_dft(3, tmp_n, tmp_n, tmp_n, int(howmany, c_intptr_t), &
        default_blocks, default_blocks, mpi_comm, pfft_transposed_out, &
        tmp_rs_n, tmp_rs_istart, tmp_fs_n, tmp_fs_istart)
    end if
#else
    tmp_rs_n = 0
    tmp_fs_n = 0
    tmp_rs_istart = 0
    tmp_fs_istart = 0
#endif
    alloc_size       = int(tmp_alloc_size, int64)
    rs_n(1:3)        = int(tmp_rs_n(3:1:-1))
    fs_n(1:3)        = int(tmp_fs_n(3:1:-1))
    rs_istart(1:3)   = int(tmp_rs_istart(3:1:-1)+1)
    fs_istart(1:3)   = int(tmp_fs_istart(3:1:-1)+1)
 
    pop_sub(pfft_get_dims)
  end subroutine pfft_get_dims
 
end module pfft_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: