cube.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio,
!! G. Bertsch, M. Oliveira, J. Alberdi-Rodriguez
!!
!! 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 cube_oct_m
  use accel_oct_m
  use affine_coordinates_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use fft_oct_m
  use global_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use mesh_cube_map_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use pfft_oct_m
  use profiling_oct_m
  use space_oct_m
 
  implicit none
  private
  public ::             &
    cube_t,             &
    dimensions_t,       &
    cube_init,          &
    cube_init_cube_map, &
    cube_point_to_process, &
    cube_partition,     &
    cube_global2local,  &
    cube_getfftlibrary, &
    cube_end
 
  type cube_t
    ! Components are public by default
    logical :: parallel_in_domains
    type(mpi_grp_t) :: mpi_grp
 
    integer :: rs_n_global(1:3)
    integer :: fs_n_global(1:3)
    integer :: rs_n(1:3)
    integer :: fs_n(1:3)
    integer :: rs_istart(1:3)
    integer :: fs_istart(1:3)
    integer :: center(1:3)
 
    integer, allocatable :: fs_ifx(:)
    integer, allocatable :: fs_ify(:)
    integer, allocatable :: fs_ifz(:)
 
    real(real64), allocatable :: Lrs(:,:)
    real(real64), allocatable :: Lfs(:,:)
 
    integer, allocatable :: np_local(:)
    integer, allocatable :: xlocal(:)
    integer, allocatable :: local(:,:)
    integer, allocatable :: np_local_fs(:)
    integer, allocatable :: xlocal_fs(:)
    integer, allocatable :: local_fs(:,:)
 
 
    type(fft_t), allocatable :: fft
    logical, private :: has_cube_mapping = .false. 
 
    real(real64)            :: spacing(3)
    ! latt is declared as allocatable as a work-around for a bug in gfortran when invoking the finalizer of latt.
    type(lattice_vectors_t), allocatable :: latt
 
    type(mesh_cube_map_t) :: cube_map
    logical :: cube_map_present = .false.
    integer :: batch_capacity
  contains
    procedure :: length => cube_length
  end type cube_t
 
  type dimensions_t
    integer :: start_xyz(1:3)
    integer :: end_xyz(1:3)
  end type dimensions_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine cube_init(cube, nn, namespace, space, spacing, coord_system, fft_type, fft_library, dont_optimize, nn_out, &
    mpi_grp, need_partition, tp_enlarge, blocksize, batch_size, nthreads, single)
    type(cube_t), intent(out) :: cube
    integer, contiguous, intent(in) :: nn(:)
    type(namespace_t), intent(in) :: namespace
    class(space_t), intent(in) :: space
    real(real64), contiguous, intent(in) :: spacing(:)
    class(coordinate_system_t), intent(in) :: coord_system
    integer, optional, intent(in) :: fft_type
    integer, optional, intent(in) :: fft_library
    logical, optional, intent(in) :: dont_optimize
    integer, optional, intent(out) :: nn_out(3)
    type(mpi_grp_t), optional, intent(in) :: mpi_grp
    logical, optional, intent(in) :: need_partition
    real(real64), optional, intent(in) :: tp_enlarge(3)
    integer, optional, intent(in) :: blocksize
    integer, optional, intent(in) :: batch_size
    ! the batch size for the FFTW and cuFFT solvers
    integer, optional, intent(in) :: nthreads
    logical, optional, intent(in) :: single
 
    type(MPI_Comm) :: comm
    integer :: tmp_n(3), fft_type_, optimize_parity(3), fft_library_, nn3d(3)
    integer :: effdim_fft, my_n(3), idir, idir2
    logical :: optimize(3)
    type(mpi_grp_t) :: mpi_grp_
    real(real64)   :: tp_enlarge_(3), lattice_vectors(3, 3)
    type(space_t) :: cube_space
 
    push_sub(cube_init)
 
    assert(all(nn > 0))
    assert(space%dim <= 3)
 
    nn3d(1:space%dim) = nn(1:space%dim)
    nn3d(space%dim+1:3) = 1
 
    cube%spacing(1:space%dim) = spacing(1:space%dim)
    cube%spacing(space%dim+1:3) = -m_one
    cube%batch_capacity = optional_default(batch_size, 1)
 
    fft_type_ = optional_default(fft_type, fft_none)
    tp_enlarge_(:) = (/m_one, m_one, m_one/)
    if (present(tp_enlarge)) tp_enlarge_(:)=tp_enlarge(:)
 
    effdim_fft = min(3, space%dim)
 
    mpi_grp_ = mpi_world
    if (present(mpi_grp)) mpi_grp_ = mpi_grp
 
    if (fft_type_ /= fft_none) then
 
      if (present(fft_library)) then
        fft_library_ = fft_library
      else
        fft_library_ = fft_default_lib
      end if
 
#ifndef HAVE_PFFT
      if (fft_library_ == fftlib_pfft) then
        write(message(1),'(a)')'You have selected the PFFT for FFT, but it was not linked.'
        call messages_fatal(1, namespace=namespace)
      end if
#endif
 
    else
      fft_library_ = fftlib_none
    end if
 
    ! Only FFTW and cuFFT can batch the FFT (howmany > 1); PFFT/NFFT/PNFFT run one transform at a
    ! time, so a batched cube on those backends is not possible. Callers that requested a batch
    ! (batch_size > 1) fall back to a per-function loop (see X(poisson_solve_batch)).
    if (fft_library_ /= fftlib_fftw .and. fft_library_ /= fftlib_accel) then
      cube%batch_capacity = 1
    end if
 
    ! Note: later we set parallel_in_domains if blocksize is given, too
    cube%parallel_in_domains = (fft_library_ == fftlib_pfft .or. fft_library_ == fftlib_pnfft)
    if (present(blocksize)) then
      assert(present(need_partition).and.need_partition)
      assert(fft_library_ == fftlib_none)
      ! For all the different FFT libraries there are strange (?)
      ! rules about how the decomposition is chosen.  What we want
      ! (for libvdwxc) is a cube parallelized according to the simple
      ! but contrary rule "just do what I say".  Hence the blocksize
      ! parameter.  (Later to be expanded to allow 2D distributions.)
      cube%rs_n_global = nn3d
      cube%fs_n_global = nn3d ! not to be used
      cube%fs_n = cube%fs_n_global ! not to be used
      cube%fs_istart = 1 ! not to be used
 
      comm = mpi_grp_%comm
      cube%parallel_in_domains = (mpi_grp_%size > 1) ! XXX whether comm size > 1
      call cube_set_blocksize(cube%rs_n_global, blocksize, mpi_grp_%rank, cube%rs_n, cube%rs_istart)
    else if (fft_library_ == fftlib_none) then
      cube%rs_n_global = nn3d
      cube%fs_n_global = nn3d
      cube%rs_n = cube%rs_n_global
      cube%fs_n = cube%fs_n_global
      cube%rs_istart = 1
      cube%fs_istart = 1
      comm = mpi_comm_undefined
      if (present(nn_out)) nn_out(1:3) = nn3d(1:3)
    else
      safe_allocate(cube%fft)
      tmp_n = nn3d
 
      optimize(1:3) = .false.
      optimize_parity(1:3) = 0
      optimize(space%periodic_dim + 1:effdim_fft) = .true.
      optimize_parity(space%periodic_dim + 1:effdim_fft) = 1
 
      if (present(dont_optimize)) then
        if (dont_optimize) optimize = .false.
      end if
 
      if (present(tp_enlarge)) call cube_tp_fft_defaults(cube, fft_library_)
 
      call fft_init(cube%fft, tmp_n, space%dim, fft_type_, fft_library_, optimize, optimize_parity, &
        comm=comm, mpi_grp = mpi_grp_, use_aligned=.true., howmany=cube%batch_capacity, &
        nthreads=nthreads, single=optional_default(single, .false.))
      if (present(nn_out)) nn_out(1:3) = tmp_n(1:3)
 
      call fft_get_dims(cube%fft, cube%rs_n_global, cube%fs_n_global, cube%rs_n, cube%fs_n, &
        cube%rs_istart, cube%fs_istart)
 
      if (present(tp_enlarge) .or. fft_library_ == fftlib_nfft .or. fft_library_ == fftlib_pnfft) then
        ! the (p)nfft nodes are the cube coordinates, so those are needed before stage1
        call cube_init_coords(cube, tp_enlarge_, cube%spacing, fft_library_)
      end if
 
      if (fft_library_ == fftlib_nfft .or. fft_library_ == fftlib_pnfft) then
        call fft_init_stage1(cube%fft, namespace, cube%Lrs, cube%rs_n_global, &
          dft_nodes = .not. present(tp_enlarge))
        !set local dimensions after stage1 - needed for PNFFT
        call fft_get_dims(cube%fft, cube%rs_n_global, cube%fs_n_global, cube%rs_n, cube%fs_n, &
          cube%rs_istart, cube%fs_istart)
      end if
 
    end if
 
    if (.not. allocated(cube%Lrs)) then
      call cube_init_coords(cube, tp_enlarge_, cube%spacing, fft_library_)
    end if
 
    cube%center(1:3) = cube%rs_n_global(1:3)/2 + 1
 
    call cube_init_fourier_mode_numbers_mapping(cube)
 
    call mpi_grp_init(cube%mpi_grp, comm)
 
    ! Initialize mapping only if needed
    if (present(need_partition) .and. cube%parallel_in_domains) then
      cube%has_cube_mapping = need_partition
    else
      cube%has_cube_mapping = .false.
    end if
    if (cube%has_cube_mapping) then
      call cube_do_mapping(cube, fs = fft_library_ == fftlib_pnfft)
    end if
 
    if (cube%parallel_in_domains) call cube_partition_messages_debug(cube, namespace)
 
    select type (coord_system)
    class is (affine_coordinates_t)
      ! We are constructing a lattice vector for the cube
      ! This differs from the actual lattice vectors if it is not an integer multiple of the spacing
      ! mesh%idx%ll is "general" in aperiodic directions,
      ! but "periodic" in periodic directions.
      my_n(1:space%periodic_dim) = cube%rs_n_global(1:space%periodic_dim) + 1
      my_n(space%periodic_dim + 1:space%dim) = cube%rs_n_global(space%periodic_dim + 1:space%dim)
 
      lattice_vectors = m_zero
      do idir = 1, space%dim
        do idir2 = 1, space%dim
          lattice_vectors(idir2, idir) = cube%spacing(idir) * (my_n(idir) - 1) * coord_system%basis%vectors(idir2, idir)
        end do
      end do
      do idir = space%dim + 1, 3
        lattice_vectors(idir, idir) = m_one
      end do
 
      cube_space%dim = 3
      cube_space%periodic_dim = space%periodic_dim
      safe_allocate(cube%latt)
      cube%latt = lattice_vectors_t(namespace, cube_space, lattice_vectors)
    class default
      message(1) = "The cube only support affine coordinate systems."
      call messages_fatal(1, namespace=namespace)
    end select
 
    pop_sub(cube_init)
  end subroutine cube_init
 
  ! ---------------------------------------------------------
  subroutine cube_end(cube)
    type(cube_t), intent(inout) :: cube
 
    push_sub(cube_end)
 
    if (allocated(cube%fft)) then
      call fft_end(cube%fft)
      safe_deallocate_a(cube%fft)
    end if
 
    if (cube%has_cube_mapping) then
      safe_deallocate_a(cube%np_local)
      safe_deallocate_a(cube%xlocal)
      safe_deallocate_a(cube%local)
 
      safe_deallocate_a(cube%np_local_fs)
      safe_deallocate_a(cube%xlocal_fs)
      safe_deallocate_a(cube%local_fs)
    end if
 
    if (cube%cube_map_present) then
      call mesh_cube_map_end(cube%cube_map)
    end if
 
    safe_deallocate_a(cube%Lrs)
    safe_deallocate_a(cube%Lfs)
 
    safe_deallocate_a(cube%latt)
    safe_deallocate_a(cube%fs_ifx)
    safe_deallocate_a(cube%fs_ify)
    safe_deallocate_a(cube%fs_ifz)
 
    pop_sub(cube_end)
  end subroutine cube_end
 
  pure function cube_length(cube) result(l)
    class(cube_t), intent(in) :: cube
    real(real64) :: l(3)
    l = cube%rs_n_global(1:3) * cube%spacing(1:3)
  end function cube_length
 
  subroutine cube_init_fourier_mode_numbers_mapping(cube)
    class(cube_t), intent(inout) :: cube
 
    integer :: lx, ly, lz
 
    push_sub(cube_init_fourier_mode_numbers_mapping)
 
    safe_allocate(cube%fs_ifx(1:max(1, cube%fs_n(1))))
    safe_allocate(cube%fs_ify(1:max(1, cube%fs_n(2))))
    safe_allocate(cube%fs_ifz(1:max(1, cube%fs_n(3))))
 
    do lx = 1, cube%fs_n(1)
      cube%fs_ifx(lx) = pad_feq(cube%fs_istart(1) + lx - 1, cube%rs_n_global(1), .true.)
    end do
    do ly = 1, cube%fs_n(2)
      cube%fs_ify(ly) = pad_feq(cube%fs_istart(2) + ly - 1, cube%rs_n_global(2), .true.)
    end do
    do lz = 1, cube%fs_n(3)
      cube%fs_ifz(lz) = pad_feq(cube%fs_istart(3) + lz - 1, cube%rs_n_global(3), .true.)
    end do
 
    pop_sub(cube_init_fourier_mode_numbers_mapping)
 
  end subroutine cube_init_fourier_mode_numbers_mapping
 
 
  ! ---------------------------------------------------------
  subroutine cube_tp_fft_defaults(cube, fft_library)
    type(cube_t), intent(inout) :: cube
    integer,      intent(in)    :: fft_library
 
    push_sub(cube_tp_fft_defaults)
    select case (fft_library)
    case (fftlib_pnfft)
      cube%fft%pnfft%set_defaults = .true.
      cube%fft%pnfft%m = 2
      cube%fft%pnfft%sigma = 1.1_real64
 
    case default
      !do nothing
    end select
 
    pop_sub(cube_tp_fft_defaults)
  end subroutine cube_tp_fft_defaults
 
 
  ! ---------------------------------------------------------
  subroutine cube_init_coords(cube, tp_enlarge, spacing, fft_library)
    type(cube_t), intent(inout) :: cube
    real(real64), intent(in) :: tp_enlarge(3)
    real(real64), intent(in) :: spacing(3)
    integer,      intent(in) :: fft_library
 
    real(real64)   :: temp
    integer :: ii, nn(3), maxn, idim
 
    push_sub(cube_init_coords)
 
 
    nn(1:3) = cube%fs_n_global(1:3)
 
    maxn = maxval(nn)
    safe_allocate(cube%Lrs(1:maxn, 1:3))
    cube%Lrs(:,:) = m_zero
 
    !! Real space coordinates
    do idim = 1,3
      if (tp_enlarge(idim) > m_one) then
        do ii = 2, nn(idim) - 1
          cube%Lrs(ii, idim) = (ii - int(nn(idim)/2) -1) * spacing(idim)
        end do
        cube%Lrs(1,        idim) = (-int(nn(idim)/2)) * spacing(idim) * tp_enlarge(idim)
        cube%Lrs(nn(idim), idim) = (int(nn(idim)/2))  * spacing(idim) * tp_enlarge(idim)
      else
        do ii = 1, nn(idim)
          cube%Lrs(ii, idim) = (ii - int(nn(idim)/2) -1) * spacing(idim)
        end do
      end if
    end do
 
 
    !! Fourier space coordinates
    if (fft_library /= fftlib_none) then
 
      safe_allocate(cube%Lfs(1:maxn, 1:3))
      cube%Lfs(:,:) = m_zero
 
      do idim = 1,3
        temp = m_two * m_pi / (nn(idim) * spacing(idim))
!temp = M_PI / (nn * spacing(1))
        do ii = 1, nn(idim)
          if (fft_library == fftlib_nfft .or. fft_library == fftlib_pnfft) then
            !The Fourier space is shrunk by the tp_enlarge factor
            !cube%Lfs(ii, 1:3) = (ii - nn/2 - 1)*temp/tp_enlarge
!HH NOTE:
!not sure this is the right general factor
            cube%Lfs(ii, idim) = (ii - nn(idim)/2 - 1)*temp/tp_enlarge(idim)
          else
            cube%Lfs(ii, idim) = pad_feq(ii,nn(idim), .true.) * temp
          end if
        end do
      end do
    end if
 
    pop_sub(cube_init_coords)
  end subroutine cube_init_coords
 
 
  ! ---------------------------------------------------------
  logical function cube_global2local(cube, ixyz, lxyz) result(is_here)
    type(cube_t), intent(in)  :: cube
    integer,      intent(in)  :: ixyz(3)
    integer,      intent(out) :: lxyz(3)
 
    lxyz(1) = ixyz(1) - cube%rs_istart(1) + 1
    lxyz(2) = ixyz(2) - cube%rs_istart(2) + 1
    lxyz(3) = ixyz(3) - cube%rs_istart(3) + 1
    is_here = lxyz(1) >= 1 .and. lxyz(1) <= cube%rs_n(1) .and. &
      lxyz(2) >= 1 .and. lxyz(2) <= cube%rs_n(2) .and. &
      lxyz(3) >= 1 .and. lxyz(3) <= cube%rs_n(3)
 
  end function cube_global2local
 
 
  ! ---------------------------------------------------------
  integer function cube_getfftlibrary(cube) result(fft_library)
    type(cube_t), intent(in)  :: cube
 
    if (allocated(cube%fft)) then
      fft_library = cube%fft%library
    else
      fft_library = fftlib_none
    end if
  end function cube_getfftlibrary
 
  ! ---------------------------------------------------------
  subroutine cube_do_mapping(cube, fs)
    type(cube_t), intent(inout) :: cube
    logical,      intent(in)    :: fs
 
    integer :: tmp_local(6), position, process, ix, iy, iz, index
    integer, allocatable :: local_sizes(:)
    integer(int64) :: number_points
 
    push_sub(cube_do_mapping)
 
    !!BEGIN:gather the local information into a unique vector.
    !!do a gather in 3d of all the box, into a loop
    tmp_local(1) = cube%rs_istart(1)
    tmp_local(2) = cube%rs_istart(2)
    tmp_local(3) = cube%rs_istart(3)
    tmp_local(4) = cube%rs_n(1)
    tmp_local(5) = cube%rs_n(2)
    tmp_local(6) = cube%rs_n(3)
 
    if (cube%parallel_in_domains) then
      safe_allocate(local_sizes(1:6*cube%mpi_grp%size))
      call profiling_in("CUBE_GAT")
      call cube%mpi_grp%allgather(tmp_local, 6, mpi_integer, local_sizes, 6, mpi_integer)
      call profiling_out("CUBE_GAT")
    else
      safe_allocate(local_sizes(1:6))
      local_sizes = tmp_local
    end if
 
    call profiling_in("CUBE_MAP")
 
    safe_allocate(cube%xlocal(1:cube%mpi_grp%size))
    safe_allocate(cube%np_local(1:cube%mpi_grp%size))
    ! make sure we do not run into integer overflow here
    number_points = cube%rs_n_global(1) * cube%rs_n_global(2)
    number_points = number_points * cube%rs_n_global(3)
    if (number_points >= huge(0)) then
      message(1) = "Error: too many points for the normal cube. Please try to use a distributed FFT."
      call messages_fatal(1)
    end if
    safe_allocate(cube%local(1:cube%rs_n_global(1)*cube%rs_n_global(2)*cube%rs_n_global(3), 1:3))
 
    index = 1
    do process = 1, cube%mpi_grp%size
      position = ((process-1)*6)+1
      if (position == 1) then
        cube%xlocal(1) = 1
        cube%np_local(1) = local_sizes(4)*local_sizes(5)*local_sizes(6)
      else
        ! calculate the begin index and size of each process
        cube%xlocal(process) = cube%xlocal(process-1) + cube%np_local(process-1)
        cube%np_local(process) = local_sizes(position+3)*local_sizes(position+4)*local_sizes(position+5)
      end if
 
      ! save the mapping between the global x,y,z and the global index
      ! and determine which partition the point belongs to
      do iz = local_sizes(position+2), local_sizes(position+2)+local_sizes(position+5)-1
        do iy = local_sizes(position+1), local_sizes(position+1)+local_sizes(position+4)-1
          do ix = local_sizes(position), local_sizes(position)+local_sizes(position+3)-1
            cube%local(index, 1) = ix
            cube%local(index, 2) = iy
            cube%local(index, 3) = iz
            index = index + 1
          end do
        end do
      end do
    end do
 
    call profiling_out("CUBE_MAP")
 
    if (optional_default(fs,.false.)) then
 
      tmp_local(1) = cube%fs_istart(1)
      tmp_local(2) = cube%fs_istart(2)
      tmp_local(3) = cube%fs_istart(3)
      tmp_local(4) = cube%fs_n(1)
      tmp_local(5) = cube%fs_n(2)
      tmp_local(6) = cube%fs_n(3)
 
      local_sizes = 0
      if (cube%parallel_in_domains) then
        call profiling_in("CUBE_GAT_FS")
        call cube%mpi_grp%allgather(tmp_local, 6, mpi_integer, local_sizes, 6, mpi_integer)
        call profiling_out("CUBE_GAT_FS")
      else
        local_sizes = tmp_local
      end if
 
      call profiling_in("CUBE_MAP_FS")
 
      safe_allocate(cube%xlocal_fs(1:cube%mpi_grp%size))
      safe_allocate(cube%np_local_fs(1:cube%mpi_grp%size))
      ! make sure we do not run into integer overflow here
      number_points = cube%fs_n_global(1) * cube%fs_n_global(2)
      number_points = number_points * cube%fs_n_global(3)
      if (number_points >= huge(0)) then
        message(1) = "Error: too many points for the normal cube. Please try to use a distributed FFT."
        call messages_fatal(1)
      end if
      safe_allocate(cube%local_fs(1:cube%fs_n_global(1)*cube%fs_n_global(2)*cube%fs_n_global(3), 1:3))
 
      index = 1
      do process = 1, cube%mpi_grp%size
        position = ((process-1)*6)+1
        if (position == 1) then
          cube%xlocal_fs(1) = 1
          cube%np_local_fs(1) = local_sizes(4)*local_sizes(5)*local_sizes(6)
        else
          ! calculate the begin index and size of each process
          cube%xlocal_fs(process) = cube%xlocal_fs(process-1) + cube%np_local_fs(process-1)
          cube%np_local_fs(process) = local_sizes(position+3)*local_sizes(position+4)*local_sizes(position+5)
        end if
 
        ! save the mapping between the global x,y,z and the global index
        ! and determine which partition the point belongs to
        do iz = local_sizes(position+2), local_sizes(position+2)+local_sizes(position+5)-1
          do iy = local_sizes(position+1), local_sizes(position+1)+local_sizes(position+4)-1
            do ix = local_sizes(position), local_sizes(position)+local_sizes(position+3)-1
              cube%local_fs(index, 1) = ix
              cube%local_fs(index, 2) = iy
              cube%local_fs(index, 3) = iz
              index = index + 1
            end do
          end do
        end do
      end do
 
      call profiling_out("CUBE_MAP_FS")
 
    end if
 
 
 
    safe_deallocate_a(local_sizes)
 
    pop_sub(cube_do_mapping)
  end subroutine cube_do_mapping
 
  !!> Given a x, y, z point of the cube, it returns the corresponding process
  !!
  !! last_found is used to speed-up the search
  integer pure function cube_point_to_process(mpi_grp, xyz, part) result(process)
    type(mpi_grp_t),    intent(in) :: mpi_grp
    integer,            intent(in) :: xyz(1:3)
    type(dimensions_t), contiguous, intent(in) :: part(:)
 
    integer :: proc
    logical :: found
 
    ! No PUSH/POP because it is a PURE function
 
    found = .false.
    do proc = 1, mpi_grp%size
      !Compare XYZ index
      if (all(xyz >= part(proc)%start_xyz) .and. all(xyz <= part(proc)%end_xyz)) then
        process = proc
        found = .true.
        exit
      end if
    end do
 
    ! An error message should be raised, if this point is reached
    if (.not. found) then
      process = -1
    end if
 
  end function cube_point_to_process
 
  ! Sets a 1D decomposition with fixed-size blocks over the last (least-contiguous) axis.
  ! Each core will have <blocksize> slices except the last one which will typically have
  ! less.  (In some cases, there can be multiple trailing cores without any slices.)
  subroutine cube_set_blocksize(rs_n_global, blocksize, rank, rs_n, rs_istart)
    integer,  intent(in) :: rs_n_global(1:3)
    integer,  intent(in) :: blocksize
    integer,  intent(in) :: rank
    integer, intent(out) :: rs_n(1:3)
    integer, intent(out) :: rs_istart(1:3)
 
    integer :: imin, imax
 
    rs_n = rs_n_global
    rs_istart = 1
 
    imin = min(blocksize * rank, rs_n_global(3))
    imax = min(imin + blocksize, rs_n_global(3))
    rs_istart(3) = 1 + imin
    rs_n(3) = imax - imin
  end subroutine cube_set_blocksize
 
  ! ---------------------------------------------------------
  subroutine cube_partition(cube, part)
    type(cube_t),       intent(in)  :: cube
    type(dimensions_t), intent(out) :: part(:)
 
    integer :: tmp_local(6), position, process
    integer, allocatable :: local_sizes(:)
 
    push_sub(cube_partition)
 
    !!gather the local information into a unique vector.
    tmp_local(1) = cube%rs_istart(1)
    tmp_local(2) = cube%rs_istart(2)
    tmp_local(3) = cube%rs_istart(3)
    tmp_local(4) = cube%rs_n(1)
    tmp_local(5) = cube%rs_n(2)
    tmp_local(6) = cube%rs_n(3)
 
    if (cube%parallel_in_domains) then
      safe_allocate(local_sizes(1:6*cube%mpi_grp%size))
      call cube%mpi_grp%allgather(tmp_local, 6, mpi_integer, local_sizes, 6, mpi_integer)
    else
      safe_allocate(local_sizes(1:6))
      local_sizes(:) = tmp_local(:)
    end if
 
    do process = 1, cube%mpi_grp%size
      position = ((process-1)*6)+1
 
      part(process)%start_xyz(1) = local_sizes(position)
      part(process)%start_xyz(2) = local_sizes(position+1)
      part(process)%start_xyz(3) = local_sizes(position+2)
      part(process)%end_xyz(1)   = local_sizes(position)+local_sizes(position+3)-1
      part(process)%end_xyz(2)   = local_sizes(position+1)+local_sizes(position+4)-1
      part(process)%end_xyz(3)   = local_sizes(position+2)+local_sizes(position+5)-1
 
    end do
 
    pop_sub(cube_partition)
  end subroutine cube_partition
 
  ! ---------------------------------------------------------
  subroutine cube_partition_messages_debug(cube, namespace)
    type(cube_t),      intent(in) :: cube
    type(namespace_t), intent(in) :: namespace
 
    integer          :: nn, ii, jj, kk ! Counters.
    integer          :: ixyz(3)        ! Current value of xyz
    integer          :: npart
    integer          :: iunit          ! For debug output to files.
    character(len=3) :: filenum
    type(dimensions_t), allocatable :: part(:)
 
    push_sub(cube_partition_messages_debug)
 
    if (debug%info) then
      safe_allocate(part(1:cube%mpi_grp%size))
      call cube_partition(cube, part)
 
      if (mpi_world%is_root()) then
        call io_mkdir('debug/cube_partition', namespace)
        npart = cube%mpi_grp%size
 
        ! Debug output. Write points of each partition in a different file.
        do nn = 1, npart
 
          write(filenum, '(i3.3)') nn
 
          iunit = io_open('debug/cube_partition/cube_partition.'//filenum, &
            namespace, action='write')
          do kk = 1, cube%rs_n_global(3)
            do jj = 1, cube%rs_n_global(2)
              do ii = 1, cube%rs_n_global(1)
                ixyz(1) = ii
                ixyz(2) = jj
                ixyz(3) = kk
                if (cube_point_to_process(cube%mpi_grp, ixyz, part) == nn) then
                  write(iunit, '(3i8)') ii, jj, kk
                end if
              end do
            end do
          end do
          call io_close(iunit)
        end do
 
 
      end if
 
      safe_deallocate_a(part)
    end if
 
    call cube%mpi_grp%barrier()
 
    pop_sub(cube_partition_messages_debug)
  end subroutine cube_partition_messages_debug
 
  ! ---------------------------------------------------------
  subroutine cube_init_cube_map(cube, mesh)
    type(cube_t),  intent(inout) :: cube
    class(mesh_t), intent(in)    :: mesh
 
    push_sub(cube_init_cube_map)
 
    call mesh_cube_map_init(cube%cube_map, mesh, mesh%np)
    cube%cube_map_present = .true.
 
    pop_sub(cube_init_cube_map)
  end subroutine cube_init_cube_map
end module cube_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: