eschnett/cactusamrex - Change Q4BI3KHPY7G275UHURXZ23GQTMCI7OQAVQOFK7XR62KFH7SFHRQQC

CarpetX: Use thorn Arith for small vectors

Created by Erik Schnetter on April 2, 2020

Q4BI3KHPY7G275UHURXZ23GQTMCI7OQAVQOFK7XR62KFH7SFHRQQC

Dependencies

In channels

master

main

Change contents

Insertion in CarpetX/configuration.ccl at line 4 [4.19462]
[3.24]
[2.3571]
```
REQUIRES Arith
```
Insertion in CarpetX/interface.ccl at line 10 [4.17976]
[5.2]
[5.2]
```
USES INCLUDE HEADER: vect.hxx
```

File deletion: vect.hxx

BF:BFD[6.2677] → [2.31772:31782]

BF:BFD[2.31782] → [2.31783:31783]

B:BD[2.31783] → [2.31784:44826]

#ifndef VECT_HXX
#define VECT_HXX
#include <cctk.h>
#undef copysign
#undef fpclassify
#undef isfinite
#undef isinf
#undef isnan
#undef isnormal
#undef signbit
#include <array>
#include <initializer_list>
#include <type_traits>
#include <utility>
#include <vector>
namespace CarpetX {
using namespace std;
constexpr int dim = 3;
////////////////////////////////////////////////////////////////////////////////
// TODO: This is broken (at least on MacOS, with gcc 9.2, with -O0).
// namespace {
// template <typename T, size_t N, typename... Ts>
// CCTK_ATTRIBUTE_ALWAYS_INLINE
// constexpr enable_if_t<(sizeof...(Ts) == N), array<T, N> >
// array_from_initializer_list(const T *const beg, const T *const end,
//                             const Ts &... xs) {
//   return array<T, N>{xs...};
// }
//
// template <typename T, size_t N, typename... Ts>
// CCTK_ATTRIBUTE_ALWAYS_INLINE
// constexpr enable_if_t<(sizeof...(Ts) < N), array<T, N> >
// array_from_initializer_list(const T *const beg, const T *const end,
//                             const Ts &... xs) {
//   return array_from_initializer_list<T, N>(beg + 1, end, xs...,
//                                            beg < end ? *beg : T{});
// }
// } // namespace
//
// template <typename T, size_t N>
//  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE array<T, N>
// array_from_initializer_list(initializer_list<T> l) {
//   return array_from_initializer_list<T, N>(l.begin(), l.end());
// }
// namespace {
// template <typename T, size_t N, typename... Ts>
// CCTK_ATTRIBUTE_ALWAYS_INLINE
// constexpr enable_if_t<(sizeof...(Ts) == N), array<T, N> >
// array_from_initializer_list(const T *const beg, const T *const end,
//                             const Ts &... xs) {
//   return array<T, N>{xs...};
// }
// template <typename T, size_t N, typename... Ts>
// CCTK_ATTRIBUTE_ALWAYS_INLINE
// constexpr enable_if_t<(sizeof...(Ts) < N), array<T, N> >
// array_from_initializer_list(const T *const beg, const T *const end,
//                             const Ts &... xs) {
//   return array_from_initializer_list<T, N>(beg, end - 1,
//                                            beg < end ? *(end - 1) : T{},
//                                            xs...);
// }
// } // namespace
//
// template <typename T, size_t N>
//  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE array<T, N>
// array_from_initializer_list(initializer_list<T> l) {
//   return array_from_initializer_list<T, N>(l.begin(), l.end());
// }
template <typename T, size_t N>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 3), array<T, N> >
array_from_initializer_list(initializer_list<T> l) {
  assert(l.size() >= N);
  const T *restrict const p = l.begin();
  return {p[0], p[1], p[2]};
}
template <typename T, size_t N>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 6), array<T, N> >
array_from_initializer_list(initializer_list<T> l) {
  assert(l.size() >= N);
  const T *restrict const p = l.begin();
  return {p[0], p[1], p[2], p[3], p[4], p[5]};
}
template <typename T, size_t N>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 3), array<T, N> >
array_from_vector(vector<T> &&v) {
  assert(v.size() >= N);
  typename vector<T>::iterator const p = v.begin();
  return {move(p[0]), move(p[1]), move(p[2])};
}
template <typename T, size_t N>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 6), array<T, N> >
array_from_vector(vector<T> &&v) {
  assert(v.size() >= N);
  typename vector<T>::iterator const p = v.begin();
  return {move(p[0]), move(p[1]), move(p[2]),
          move(p[3]), move(p[4]), move(p[5])};
}
static_assert(static_cast<const array<int, 3> &>(
                  array_from_initializer_list<int, 3>({0, 1, 2}))[0] == 0,
              "");
static_assert(static_cast<const array<int, 3> &>(
                  array_from_initializer_list<int, 3>({0, 1, 2}))[1] == 1,
              "");
static_assert(static_cast<const array<int, 3> &>(
                  array_from_initializer_list<int, 3>({0, 1, 2}))[2] == 2,
              "");
////////////////////////////////////////////////////////////////////////////////
template <typename T, int D> struct vect {
  array<T, D> elts;
  // initializes all elts to zero
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect() : elts() {}
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(const array<T, D> &arr)
      : elts(arr) {}
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(initializer_list<T> lst)
      : elts(array_from_initializer_list<T, D>(lst)) {
#ifdef CCTK_DEBUG
    assert(lst.size() == D);
#endif
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(const vector<T> &vec)
      : elts(array_from_vector<T, D>(vec)) {
#ifdef CCTK_DEBUG
    assert(vec.size() == D);
#endif
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(vector<T> &&vec)
      : elts(array_from_vector<T, D>(move(vec))) {}
  static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect pure(T a) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = a;
    return r;
  }
  static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect unit(int dir) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = d == dir;
    return r;
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE const T &operator[](int d) const {
#ifdef CCTK_DEBUG
    assert(d >= 0 && d < D);
#endif
    return elts[d];
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T &operator[](int d) {
#ifdef CCTK_DEBUG
    assert(d >= 0 && d < D);
#endif
    return elts[d];
  }
  template <typename F,
            typename R = remove_cv_t<remove_reference_t<result_of_t<F(T)> > > >
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R map(const F &f) const {
    vect<R, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = f(elts[d]);
    return r;
  }
  template <
      typename F, typename U,
      typename R = remove_cv_t<remove_reference_t<result_of_t<F(T, U)> > > >
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R map2(const F &f,
                                                const vect<U, D> &y) const {
    vect<R, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = f(elts[d], y.elts[d]);
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect &x) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = +x.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect &x) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = -x.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect &x,
                                                               const vect &y) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] + y.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect &x,
                                                               const vect &y) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] - y.elts[d];
    return r;
  }
  // friend  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect
  // &x, const array<T, D> &y) {
  //   return x + vect(y);
  // }
  // friend  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect
  // &x, const array<T, D> &y) {
  //   return x - vect(y);
  // }
  // friend  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect
  // &x, T a) {
  //   vect r;
  //   for (int d = 0; d < D; ++d)
  //     r.elts[d] = x.elts[d] + a;
  //   return r;
  // }
  // friend  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect
  // &x, T a) {
  //   vect r;
  //   for (int d = 0; d < D; ++d)
  //     r.elts[d] = x.elts[d] - a;
  //   return r;
  // }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator*(T a,
                                                               const vect &x) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = a * x.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator*(const vect &x,
                                                               T a) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] * a;
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D> operator!(
      const vect &x) {
    vect<bool, D> r;
    for (int d = 0; d < dim; ++d)
      r.elts[d] = !x.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator&&(const vect &x, const vect &y) {
    vect<bool, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] && y.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator||(const vect &x, const vect &y) {
    vect<bool, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] || y.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator==(const vect &x, const vect &y) {
    vect<bool, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] == y.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator!=(const vect &x, const vect &y) {
    return !(x == y);
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator<(const vect &x, const vect &y) {
    vect<bool, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] < y.elts[d];
    return r;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator>(const vect &x, const vect &y) {
    return y < x;
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator<=(const vect &x, const vect &y) {
    return !(x > y);
  }
  friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>
  operator>=(const vect &x, const vect &y) {
    return !(x < y);
  }
  template <typename U>
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<U, D>
  ifelse(const CarpetX::vect<U, D> &x, const CarpetX::vect<U, D> &y) const {
    vect<U, D> r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = elts[d] ? x.elts[d] : y.elts[d];
    return r;
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE bool all() const {
    bool r = true;
    for (int d = 0; d < D; ++d)
      r &= elts[d];
    return r;
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE bool any() const {
    bool r = false;
    for (int d = 0; d < D; ++d)
      r |= elts[d];
    return r;
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T maxabs() const {
    T r = 0;
    for (int d = 0; d < D; ++d)
      r = fmax(r, fabs(elts[d]));
    return r;
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T prod() const {
    T r = 1;
    for (int d = 0; d < D; ++d)
      r *= elts[d];
    return r;
  }
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T sum() const {
    T r = 0;
    for (int d = 0; d < D; ++d)
      r += elts[d];
    return r;
  }
  friend ostream &operator<<(ostream &os, const vect &x) {
    os << "[";
    for (int d = 0; d < D; ++d) {
      if (d > 0)
        os << ",";
      os << x.elts[d];
    }
    os << "]";
    return os;
  }
};
} // namespace CarpetX
namespace std {
template <typename T, int D> struct equal_to<CarpetX::vect<T, D> > {
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE bool
  operator()(const CarpetX::vect<T, D> &lhs,
             const CarpetX::vect<T, D> &rhs) const {
    return equal_to<array<T, D> >()(lhs.elts, rhs.elts);
  }
};
template <typename T, int D> struct less<CarpetX::vect<T, D> > {
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE bool
  operator()(const CarpetX::vect<T, D> &lhs,
             const CarpetX::vect<T, D> &rhs) const {
    return less<array<T, D> >(lhs.elts, rhs.elts);
  }
};
template <typename T, int D>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CarpetX::vect<T, D>
abs(const CarpetX::vect<T, D> &x) {
  CarpetX::vect<T, D> r;
  for (int d = 0; d < D; ++d)
    r.elts[d] = abs(x.elts[d]);
  return r;
}
template <typename T, int D>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CarpetX::vect<T, D>
max(const CarpetX::vect<T, D> &x, const CarpetX::vect<T, D> &y) {
  CarpetX::vect<T, D> r;
  for (int d = 0; d < D; ++d)
    r.elts[d] = max(x.elts[d], y.elts[d]);
  return r;
}
template <typename T, int D>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CarpetX::vect<T, D>
min(const CarpetX::vect<T, D> &x, const CarpetX::vect<T, D> &y) {
  CarpetX::vect<T, D> r;
  for (int d = 0; d < D; ++d)
    r.elts[d] = min(x.elts[d], y.elts[d]);
  return r;
}
} // namespace std
namespace CarpetX {
namespace static_test {
constexpr vect<int, 3> vect1{0, 1, 2};
static_assert(vect1[0] == 0, "");
static_assert(vect1[1] == 1, "");
static_assert(vect1[2] == 2, "");
constexpr vect<vect<int, 3>, 3> vect2{
    {100, 101, 102},
    {110, 111, 112},
    {120, 121, 122},
};
static_assert(vect2[0][0] == 100, "");
static_assert(vect2[0][1] == 101, "");
static_assert(vect2[0][2] == 102, "");
static_assert(vect2[1][0] == 110, "");
static_assert(vect2[1][1] == 111, "");
static_assert(vect2[1][2] == 112, "");
static_assert(vect2[2][0] == 120, "");
static_assert(vect2[2][1] == 121, "");
static_assert(vect2[2][2] == 122, "");
} // namespace static_test
} // namespace CarpetX
#endif // #ifndef VECT_HXX

Replacement in CarpetX/src/driver.hxx at line 4 [4.15275]
B:BD[2.4006] → [2.4006:4026]
```
#include "vect.hxx"
```
[2.4006]
[5.4443]
```
#include "loop.hxx"
```
Insertion in CarpetX/src/driver.hxx at line 34 [4.15275]
[4.15556]
[7.0]
```
  using Loop::dim;
```
Insertion in CarpetX/src/interpolate.cxx at line 252 [8.373]
[8.4786]
[8.4786]
```
  // CCTK_VINFO("interpolating");
```
Insertion in CarpetX/src/interpolate.cxx at line 255 [8.373]
[8.4876]
[8.4876]
```
    // CCTK_VINFO("interpolating level %d", level);
```

Insertion in CarpetX/src/interpolate.cxx at line 268 [8.373]

[8.5341]

        // CCTK_VINFO("interpolating level %d, variable %d", level, v);

Insertion in CarpetX/src/interpolate.cxx at line 392 [8.373]

[8.8141]


  // CCTK_VINFO("interpolation results");
  // for (const auto &proc_result : results) {
  //   const int p = proc_result.first;
  //   const auto &result = proc_result.second;
  //   CCTK_VINFO("[%d] count=%zd", p, result.size());
  // }

Insertion in CarpetX/src/interpolate.cxx at line 401 [8.373]
[8.8170]
[8.8170]
```
  // CCTK_VINFO("collecting results");
```

Insertion in CarpetX/src/interpolate.cxx at line 418 [8.373]

[8.8716]

  // for (int p = 0; p < nprocs; ++p)
  //   CCTK_VINFO("[%d] senddispl=%d sendcount=%d", p, senddispls.at(p),
  //              sendcounts.at(p));
  // CCTK_VINFO("sendcount=%d", sendcount);

Insertion in CarpetX/src/interpolate.cxx at line 431 [8.373]

[8.9014]

  // for (int p = 0; p < nprocs; ++p)
  //   CCTK_VINFO("[%d] recvdispl=%d recvcount=%d", p, recvdispls.at(p),
  //              recvcounts.at(p));
  // CCTK_VINFO("recvcount=%d", recvcount);

Insertion in CarpetX/src/loop.hxx at line 28 [5.5481]
[5.5630]
[9.398]
```
using Arith::vect;
```
Replacement in CarpetX/src/loop.hxx at line 31 [5.5481]
B:BD[9.399] → [2.18237:18278]
```
using CarpetX::dim;
using CarpetX::vect;
```
[9.399]
[9.665]
```
constexpr int dim = 3;
```
Insertion in CarpetX/src/reduction.hxx at line 22 [10.2777]
[10.2899]
[11.3395]
```
using namespace Arith;
```

CarpetX: Use thorn Arith for small vectors

Dependencies

In channels

Change contents

Insertion in CarpetX/configuration.ccl at line 4 [4.19462]

Insertion in CarpetX/interface.ccl at line 10 [4.17976]

File deletion: vect.hxx

Replacement in CarpetX/src/driver.hxx at line 4 [4.15275]

Insertion in CarpetX/src/driver.hxx at line 34 [4.15275]

Insertion in CarpetX/src/interpolate.cxx at line 252 [8.373]

Insertion in CarpetX/src/interpolate.cxx at line 255 [8.373]

Insertion in CarpetX/src/interpolate.cxx at line 268 [8.373]

Insertion in CarpetX/src/interpolate.cxx at line 392 [8.373]

Insertion in CarpetX/src/interpolate.cxx at line 401 [8.373]

Insertion in CarpetX/src/interpolate.cxx at line 418 [8.373]

Insertion in CarpetX/src/interpolate.cxx at line 431 [8.373]

Insertion in CarpetX/src/loop.hxx at line 28 [5.5481]

Replacement in CarpetX/src/loop.hxx at line 31 [5.5481]

Insertion in CarpetX/src/reduction.hxx at line 22 [10.2777]