File blockgeometry.h

File List > code_source > templat > include > TempLat > lattice > IO > HDF5 > helpers > blockgeometry.h

Go to the documentation of this file

#ifndef TEMPLAT_LATTICE_IO_HDF5_HELPERS_BLOCKGEOMETRY_H
#define TEMPLAT_LATTICE_IO_HDF5_HELPERS_BLOCKGEOMETRY_H

/* This file is part of TempLat, available at https://cosmolattice.github.io/templat .
   Copyright 2021-2026 The TempLat authors, see AUTHORS.md.
   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Adrien Florio, Year: 2026

#include <cstddef>

#include "TempLat/parallel/device.h"

namespace TempLat
{
  struct BlockGeometryDim {
    device::Idx count = 0;  
    device::Idx first = 0;  
    device::Idx offset = 0; 
  };

  inline BlockGeometryDim computeBlockGeometryDim(device::Idx start, device::Idx size, device::Idx ini,
                                                  device::Idx end, device::Idx step, device::Idx datasetExtent)
  {
    const device::Idx stp = step < 1 ? 1 : step;

    const device::Idx lo = start > ini ? start : ini;                          // max(start, ini)
    const device::Idx hi = (start + size) < end ? (start + size) : end;        // min(start+size, end)

    // Ceiling division; both numerators are non-negative where they are used.
    const device::Idx kFirst = (lo <= ini) ? 0 : (lo - ini + stp - 1) / stp;
    const device::Idx kEnd = (hi <= ini) ? 0 : (hi - ini + stp - 1) / stp;
    const device::Idx kStop = kEnd < datasetExtent ? kEnd : datasetExtent;

    BlockGeometryDim result;
    result.count = kStop > kFirst ? kStop - kFirst : 0;
    result.first = ini + kFirst * stp;
    result.offset = kFirst;
    return result;
  }

  template <size_t NDim> struct SaveBlockGeometry {
    device::IdxArray<NDim> count{};  
    device::IdxArray<NDim> first{};  
    device::IdxArray<NDim> offset{}; 
    device::IdxArray<NDim> step{};   

    bool isEmpty() const
    {
      for (size_t d = 0; d < NDim; ++d)
        if (count[d] <= 0) return true;
      return false;
    }

    device::Idx totalCount() const
    {
      if (isEmpty()) return 0;
      device::Idx total = 1;
      for (size_t d = 0; d < NDim; ++d)
        total *= count[d];
      return total;
    }
  };

  template <size_t NDim, typename CStarts, typename CSizes, typename CGrid, typename CLimits>
  SaveBlockGeometry<NDim> computeSaveBlockGeometry(const CStarts &starts, const CSizes &sizes, const CGrid &nGrid,
                                                   bool sparsesave, const CLimits &down, const CLimits &up,
                                                   const CLimits &step, const CLimits &datasetSizes)
  {
    SaveBlockGeometry<NDim> geo{};
    for (size_t d = 0; d < NDim; ++d) {
      const device::Idx ini = sparsesave ? (device::Idx)down[d] : 0;
      const device::Idx end = sparsesave ? (device::Idx)up[d] : (device::Idx)nGrid[d];
      const device::Idx stp = sparsesave ? (device::Idx)step[d] : 1;

      const auto dim = computeBlockGeometryDim((device::Idx)starts[d], (device::Idx)sizes[d], ini, end, stp,
                                               (device::Idx)datasetSizes[d]);
      geo.count[d] = dim.count;
      geo.first[d] = dim.first;
      geo.offset[d] = dim.offset;
      geo.step[d] = stp;
    }
    return geo;
  }

  template <size_t NDim, typename CStarts, typename CSizes, typename CGrid>
  SaveBlockGeometry<NDim> computeDenseBlockGeometry(const CStarts &starts, const CSizes &sizes, const CGrid &nGrid)
  {
    // The limit containers are never read when sparsesave is false; nGrid stands in so that all
    // four share a type.
    return computeSaveBlockGeometry<NDim>(starts, sizes, nGrid, false, nGrid, nGrid, nGrid, nGrid);
  }

} // namespace TempLat

#endif