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_dft_r2c_3d, &
    pfft_local_size_dft_3d,     &
    pfft_plan_dft_c2r_3d,       &
    pfft_plan_dft_r2c_3d,       &
    pfft_plan_dft_3d,           &
    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, sign, flags, mpi_comm)
    type(C_PTR),      intent(out)   :: plan
    integer,          intent(in)    :: n(:)
    real(real64),   pointer, intent(inout) :: in(:,:,:)
    complex(real64),   pointer, intent(inout) :: out(:,:,:)
    integer,          intent(in)    :: sign
    integer,          intent(in)    :: flags
    integer,          intent(in)    :: mpi_comm
 
    integer(C_INTPTR_T) :: tmp_n(3)
 
    push_sub(pfft_prepare_plan_r2c)
    call profiling_in("PFFT_PLAN_R2C")
 
#ifdef HAVE_PFFT
    assert(sign == pfft_forward)
#endif
 
    tmp_n(1:3) = n(3:1:-1)
 
#ifdef HAVE_PFFT
    plan = pfft_plan_dft_r2c_3d(tmp_n,in,out,mpi_comm,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, sign, flags, mpi_comm)
    type(c_ptr),      intent(out)   :: plan
    integer,          intent(in)    :: n(:)
    complex(real64),   pointer, intent(inout) :: in(:,:,:)
    real(real64),   pointer, intent(inout) :: out(:,:,:)
    integer,          intent(in)    :: sign
    integer,          intent(in)    :: flags
    integer,          intent(in)    :: mpi_comm
 
    integer(C_INTPTR_T) :: tmp_n(3)
    push_sub(pfft_prepare_plan_c2r)
    call profiling_in("PFFT_PLAN_C2R")
 
#ifdef HAVE_PFFT
    assert(sign == pfft_backward)
#endif
 
    tmp_n(1:3) = n(3:1:-1)
 
#ifdef HAVE_PFFT
    plan = pfft_plan_dft_c2r_3d(tmp_n,in,out,mpi_comm, sign, &
      pfft_transposed_in + pfft_measure + pfft_destroy_input + pfft_tune)
    !call PDFFT(plan_dft_c2r_3d) (plan, tmp_n(1), in(1,1,1), out(1,1,1), mpi_comm, 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, sign, flags, mpi_comm)
    type(c_ptr),      intent(out)   :: plan
    integer,          intent(in)    :: n(:)
    complex(real64),   pointer, intent(inout) :: in(:,:,:)
    complex(real64),   pointer, intent(inout) :: out(:,:,:)
    integer,          intent(in)    :: sign
    ! 1                                              Has to be FFTW_FORWARD or FFTW_BACKWARD
    integer,          intent(in)    :: flags
    integer,          intent(in)    :: mpi_comm
 
#ifdef HAVE_PFFT
    integer(C_INT) :: pfft_flags
#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 (sign == pfft_forward) then
      pfft_flags = pfft_transposed_out
    else
      pfft_flags = pfft_transposed_in
    end if
 
    plan = pfft_plan_dft_3d(tmp_n,in,out,mpi_comm, sign, pfft_flags)
    !call PDFFT(plan_dft_3d) (plan, tmp_n(1), in(1,1,1), out(1,1,1), mpi_comm, 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)
    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(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)
 
    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
    if (is_real) then
      tmp_alloc_size = pfft_local_size_dft_r2c_3d(tmp_n, mpi_comm, pfft_transposed_out, &
        tmp_rs_n, tmp_rs_istart, tmp_fs_n, tmp_fs_istart)
      !call PDFFT(local_size_dft_r2c_3d)(tmp_alloc_size, tmp_n(1), mpi_comm, PFFT_TRANSPOSED_OUT, &
      !       tmp_rs_n(1), tmp_rs_istart(1), tmp_fs_n(1), tmp_fs_istart(1))
      fs_n_global(1) = rs_n_global(1)/2 + 1
    else
      tmp_alloc_size = pfft_local_size_dft_3d(tmp_n,mpi_comm, pfft_transposed_out, &
        tmp_rs_n,tmp_rs_istart,tmp_fs_n,tmp_fs_istart)
      !call PDFFT(local_size_dft_3d)(tmp_alloc_size, tmp_n(1), mpi_comm, PFFT_TRANSPOSED_OUT, &
      !             tmp_rs_n(1), tmp_rs_istart(1), tmp_fs_n(1), tmp_fs_istart(1))
    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: