using namespace amrex;

using namespace amrex;

#include <reduction.hxx>

        for (int n = 0; n < groupdata.numvars; ++n) {

        for (int vi = 0; vi < groupdata.numvars; ++vi) {

          buf << CCTK_VarName(groupdata.firstvarindex + n);

          buf << CCTK_VarName(groupdata.firstvarindex + vi);

          varnames.at(n) = buf.str();

          varnames.at(vi) = buf.str();

    }
    count_vars += ghext->leveldata.at(0).groupdata.at(gi).numvars;
  }

      count_vars += groupdata.numvars;

  for (int gi = 0; gi < numgroups; ++gi) {
    const int level = 0;
    const GHExt::LevelData &restrict leveldata = ghext->leveldata.at(level);
    const GHExt::LevelData::GroupData &restrict groupdata =
        leveldata.groupdata.at(gi);
    const int numvars = groupdata.numvars;
    for (int vi = 0; vi < numvars; ++vi) {
      const int tl = 0;
      reduction<CCTK_REAL> red = reduce(gi, vi, tl);
      CCTK_VINFO("maxabs(%s)=%g", CCTK_VarName(groupdata.firstvarindex + vi),
                 red.maxabs);

SRCS = driver.cxx io.cxx schedule.cxx

SRCS = driver.cxx io.cxx reduction.cxx schedule.cxx

#include <driver.hxx>
#include <reduction.hxx>
#include <AMReX_Orientation.H>
namespace AMReX {
using namespace amrex;
using namespace std;
namespace {
// Convert a (direction, face) pair to an AMReX Orientation
Orientation orient(int d, int f) {
  return Orientation(d, Orientation::Side(f));
}
} // namespace
template <typename T>
reduction<T> reduce_array(const T *restrict var, const int *restrict ash,
                          const int *restrict lsh) {
  reduction<T> red;
  constexpr int di = 1;
  const int dj = di * ash[0];
  const int dk = dj * ash[1];
  for (int k = 0; k < lsh[2]; ++k) {
    for (int j = 0; j < lsh[1]; ++j) {
#pragma omp simd
      for (int i = 0; i < lsh[0]; ++i) {
        int idx = di * i + dj * j + dk * k;
        red = reduction<T>(red, reduction<T>(var[idx]));
      }
    }
  }
  return red;
}
reduction<CCTK_REAL> reduce(int gi, int vi, int tl) {
  reduction<CCTK_REAL> red;
#pragma omp parallel reduction(reduction : red)
  {
    for (auto &restrict leveldata : ghext->leveldata) {
      auto &restrict groupdata = leveldata.groupdata.at(gi);
      MultiFab &mfab = *groupdata.mfab.at(tl);
      auto mfitinfo =
          MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32});
      for (MFIter mfi(mfab, mfitinfo); mfi.isValid(); ++mfi) {
        const Box &fbx = mfi.fabbox();
        const Box &bx = mfi.tilebox();
        const Dim3 imin = lbound(bx);
        int ash[3], lsh[3];
        for (int d = 0; d < dim; ++d)
          ash[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
        for (int d = 0; d < dim; ++d)
          lsh[d] = bx[orient(d, 1)] - bx[orient(d, 0)] + 1;
        const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
        const CCTK_REAL *restrict var = vars.ptr(imin.x, imin.y, imin.z, vi);
        red += reduce_array(var, ash, lsh);
      }
    }
  }
  return red;
}
} // namespace AMReX

#ifndef REDUCTION_HXX
#define REDUCTION_HXX
#include <cctk.h>
#include <cmath>
namespace AMReX {
using namespace std;
template <typename T> struct reduction {
  T min, max, sum, sum2;
  T count, maxabs, sumabs, sum2abs;
  reduction();
  reduction(const T &x);
  reduction(const reduction &x, const reduction &y);
  reduction operator+(const reduction &y) const { return reduction(*this, y); }
  reduction &operator+=(const reduction &y) { return *this = *this + y; }
};
template <typename T>
reduction<T>::reduction()
    : min(1.0 / 0.0), max(-1.0 / 0.0), sum(0.0), sum2(0.0), count(0.0),
      maxabs(0.0), sumabs(0.0), sum2abs(0.0) {}
template <typename T>
reduction<T>::reduction(const T &x)
    : min(x), max(x), sum(x), sum2(x * x), count(1.0), maxabs(fabs(x)),
      sumabs(fabs(x)), sum2abs(fabs(x) * fabs(x)) {}
template <typename T>
reduction<T>::reduction(const reduction &x, const reduction &y)
    : min(fmin(x.min, y.min)), max(fmax(x.max, y.max)), sum(x.sum + y.sum),
      sum2(x.sum2 + y.sum2), count(x.count + y.count),
      maxabs(fmax(x.maxabs, y.maxabs)), sumabs(x.sumabs + y.sumabs),
      sum2abs(x.sum2abs + y.sum2abs) {}
typedef reduction<CCTK_REAL> reduction_CCTK_REAL;
#pragma omp declare reduction(reduction:reduction_CCTK_REAL : omp_out += omp_in)
reduction<CCTK_REAL> reduce(int gi, int vi, int tl);
} // namespace AMReX
#endif // #ifndef REDUCTION_HXX

////////////////////////////////////////////////////////////////////////////////

vector<cGH> thread_local_cctkGH;

namespace {

} // namespace

vector<cGH> thread_local_cctkGH;

////////////////////////////////////////////////////////////////////////////////

      for (int n = 0; n < groupdata.numvars; ++n) {
        cctkGH->data[groupdata.firstvarindex + n][tl] =
            vars.ptr(imin.x, imin.y, imin.z, n);

      for (int vi = 0; vi < groupdata.numvars; ++vi) {
        cctkGH->data[groupdata.firstvarindex + vi][tl] =
            vars.ptr(imin.x, imin.y, imin.z, vi);

      for (int n = 0; n < groupdata.numvars; ++n)
        cctkGH->data[groupdata.firstvarindex + n][tl] = nullptr;

      for (int vi = 0; vi < groupdata.numvars; ++vi)
        cctkGH->data[groupdata.firstvarindex + vi][tl] = nullptr;