#include <sstream>
#include <string>
#include <type_traits>
#include <vector>

template <typename T> MPI_Datatype reduction_mpi_datatype() {
  static MPI_Datatype datatype = MPI_DATATYPE_NULL;
  if (datatype == MPI_DATATYPE_NULL) {
    MPI_Type_contiguous(sizeof(reduction<T>) / sizeof(T),
                        mpi_datatype<T>::value, &datatype);
    char name[MPI_MAX_OBJECT_NAME];
    int namelen;
    MPI_Type_get_name(mpi_datatype<T>::value, name, &namelen);
    ostringstream buf;
    buf << "reduction<" << name << ">";
    string newname = buf.str();
    MPI_Type_set_name(datatype, newname.c_str());
    MPI_Type_commit(&datatype);
  }
  return datatype;
}

  if (*datatype == MPI_FLOAT)
    return mpi_reduce_typed<float>(x, y, length);
  if (*datatype == MPI_DOUBLE)
    return mpi_reduce_typed<double>(x, y, length);
  if (*datatype == MPI_LONG_DOUBLE)
    return mpi_reduce_typed<long double>(x, y, length);

  // Analyze MPI datatype
  int num_integers, num_addresses, num_datatypes, combiner;
  MPI_Type_get_envelope(*datatype, &num_integers, &num_addresses,
                        &num_datatypes, &combiner);
  assert(combiner == MPI_COMBINER_CONTIGUOUS);
  assert(num_integers == 1);
  assert(num_addresses = 1);
  assert(num_datatypes == 1);
  vector<int> integers(num_integers);
  vector<MPI_Aint> addresses(num_addresses);
  vector<MPI_Datatype> datatypes(num_datatypes);
  MPI_Type_get_contents(*datatype, num_integers, num_addresses, num_datatypes,
                        integers.data(), addresses.data(), datatypes.data());
  MPI_Datatype inner_datatype = datatypes.at(0);
  if (inner_datatype == MPI_FLOAT)
    return mpi_reduce_typed<reduction<float> >(x, y, length);
  if (inner_datatype == MPI_DOUBLE)
    return mpi_reduce_typed<reduction<double> >(x, y, length);
  if (inner_datatype == MPI_LONG_DOUBLE)
    return mpi_reduce_typed<reduction<long double> >(x, y, length);

    MPI_Op_create(mpi_reduce, 1, &op);

    MPI_Op_create(mpi_reduce, 1 /*commutes*/, &op);

namespace {

reduction<T> reduce_array(const Geometry &geom,
                          const GHExt::LevelData::GroupData &groupdata,
                          const GridPtrDesc1 &grid, const GF3D1<const T> ptr_,
                          const Array4<const int> *restrict finemask_array4) {
  const CCTK_REAL *restrict dx = geom.CellSize();
  CCTK_REAL dV = 1.0;
  for (int d = 0; d < dim; ++d)
    dV *= dx[d];

reduction<T> reduce_array(const Array4<const T> &restrict vars, int n,
                          const array<int, dim> &imin,
                          const array<int, dim> &imax,
                          const Array4<const int> *restrict finemask, T dV) {

  grid.loop_idx(
      where_t::interior, groupdata.indextype, [&](const Loop::PointDesc &p) {
        if (!finemask_array4 || !(*finemask_array4)(grid.cactus_offset.x + p.i,
                                                    grid.cactus_offset.y + p.j,
                                                    grid.cactus_offset.z + p.k))
          red = reduction<T>(red, reduction<T>(dV, ptr_(p.I)));
      });

  for (int k = imin[2]; k < imax[2]; ++k) {
    for (int j = imin[1]; j < imax[1]; ++j) {
      // #pragma omp simd reduction(red : reduction_CCTK_REAL)
      for (int i = imin[0]; i < imax[0]; ++i) {
        bool is_masked = finemask && (*finemask)(i, j, k);
        if (!is_masked)
          red += reduction<T>(dV, vars(i, j, k, n));
      }
    }
  }

} // namespace

    const auto &restrict geom = ghext->amrcore->Geom(leveldata.level);

    unique_ptr<iMultiFab> finemask;

    unique_ptr<iMultiFab> finemask_imfab;

    const auto &restrict geom = ghext->amrcore->Geom(leveldata.level);
    const CCTK_REAL *restrict dx = geom.CellSize();
    CCTK_REAL dV = 1.0;
    for (int d = 0; d < dim; ++d)
      dV *= dx[d];

      finemask = make_unique<iMultiFab>(

      finemask_imfab = make_unique<iMultiFab>(

    for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
      GridPtrDesc1 grid(leveldata, groupdata, mfi);
      unique_ptr<Array4<const int> > finemask_array4;
      if (finemask)
        finemask_array4 = make_unique<Array4<const int> >(finemask->array(mfi));

    for (MFIter mfi(mfab, mfitinfo); mfi.isValid(); ++mfi) {
      const Box &bx = mfi.tilebox(); // current region (without ghosts)
      const array<int, dim> imin{bx.smallEnd(0), bx.smallEnd(1),
                                 bx.smallEnd(2)};
      const array<int, dim> imax{bx.bigEnd(0) + 1, bx.bigEnd(1) + 1,
                                 bx.bigEnd(2) + 1};

      const GF3D1<const CCTK_REAL> ptr_ = grid.gf3d(vars, vi);
      red += reduce_array(geom, groupdata, grid, ptr_, finemask_array4.get());

      unique_ptr<Array4<const int> > finemask;
      if (finemask_imfab)
        finemask = make_unique<Array4<const int> >(finemask_imfab->array(mfi));
      red += reduce_array(vars, vi, imin, imax, finemask.get(), dV);

  MPI_Datatype datatype = reduction_mpi_datatype<CCTK_REAL>();

  MPI_Allreduce(MPI_IN_PLACE, &red, 1, MPI_DOUBLE, op, MPI_COMM_WORLD);

  MPI_Allreduce(MPI_IN_PLACE, &red, 1, datatype, op, MPI_COMM_WORLD);

#include <ostream>

template <typename T> struct mpi_datatype;
template <> struct mpi_datatype<float> {
  static constexpr MPI_Datatype value = MPI_FLOAT;
};
template <> struct mpi_datatype<double> {
  static constexpr MPI_Datatype value = MPI_DOUBLE;
};
template <> struct mpi_datatype<long double> {
  static constexpr MPI_Datatype value = MPI_LONG_DOUBLE;
};

  if (CCTK_isnan(x))

  if (CCTK_isnan(CCTK_REAL(x)))

  if (CCTK_isnan(y))

  if (CCTK_isnan(CCTK_REAL(y)))

  if (CCTK_isnan(x))

  if (CCTK_isnan(CCTK_REAL(x)))

  if (CCTK_isnan(y))

  if (CCTK_isnan(CCTK_REAL(y)))

  friend ostream &operator<<(ostream &os, const reduction &red) {
    return os << "reduction{min:" << red.min << ",max:" << red.max
              << ",sum:" << red.sum << ",sum2:" << red.sum2
              << ",vol:" << red.vol << ",maxabs:" << red.maxabs
              << ",sumabs:" << red.sumabs << ",sum2abs:" << red.sum2abs << "}";
  }

MPI_Datatype reduction_mpi_datatype_CCTK_REAL();