eschnett/cactusamrex - Change HORXVU5QMFSZAMX6HO4O42J73IFN7SH6IWU5PGYE347QIX3TBALQC

Arith: New classes vec, mat; new function sum

Created by Erik Schnetter on June 9, 2020

HORXVU5QMFSZAMX6HO4O42J73IFN7SH6IWU5PGYE347QIX3TBALQC

Dependencies

In channels

main

master

Change contents

Insertion in Arith/interface.ccl at line 7 [3.172]

[2.44]

[3.270]

INCLUDES HEADER: mat.hxx IN mat.hxx
INCLUDES HEADER: sum.hxx IN sum.hxx
INCLUDES HEADER: vec.hxx IN vec.hxx

Replacement in Arith/src/make.code.defn at line 4 [4.6031]
B:BD[4.6108] → [4.6108:6124]
```
SRCS = dual.cxx
```
[4.6108]
[4.6124]
```
SRCS = dual.cxx mat.cxx sum.cxx vec.cxx
```

File addition: mat.cxx (-rw-r--r--)

[3.425]

#include "mat.hxx"
#include <cctk.h>
#include <functional>
namespace Arith {
using namespace std;
// This function is compiled, but not executed. The tests are "run" at
// compile time. If this function compiles, the tests pass.
void TestMat() {
  constexpr equal_to<CCTK_REAL> eq;
  using M3D = mat<CCTK_REAL, 3, DN, DN>;
  constexpr equal_to<M3D> eqm;
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>()(0, 0), 0));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>()(0, 1), 0));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>()(0, 2), 0));
  static_assert(eq(mat<CCTK_REAL, 3, UP, UP>()(1, 1), 0));
  static_assert(eq(mat<CCTK_REAL, 3, UP, UP>()(1, 2), 0));
  static_assert(eq(mat<CCTK_REAL, 3, UP, UP>()(2, 2), 0));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(0, 0), 1));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(0, 1), 2));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(0, 2), 3));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(1, 0), 2));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(1, 1), 4));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(1, 2), 5));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(2, 0), 3));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(2, 1), 5));
  static_assert(eq(mat<CCTK_REAL, 3, DN, DN>{1, 2, 3, 4, 5, 6}(2, 2), 6));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::iota1(), {0, 0, 0, 1, 1, 2}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::iota2(), {0, 1, 2, 1, 2, 2}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(0, 0), {1, 0, 0, 0, 0, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(0, 1), {0, 1, 0, 0, 0, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(0, 2), {0, 0, 1, 0, 0, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(1, 0), {0, 1, 0, 0, 0, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(1, 1), {0, 0, 0, 1, 0, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(1, 2), {0, 0, 0, 0, 1, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(2, 0), {0, 0, 1, 0, 0, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(2, 1), {0, 0, 0, 0, 1, 0}));
  static_assert(eqm(mat<CCTK_REAL, 3, DN, DN>::unit(2, 2), {0, 0, 0, 0, 0, 1}));
  constexpr M3D x = {71, 90, 14, 50, 41, 65};
  static_assert(eqm(+x, {71, 90, 14, 50, 41, 65}));
  static_assert(eqm(-x, {-71, -90, -14, -50, -41, -65}));
  constexpr M3D y = {77, 32, 67, 5, 81, 10};
  static_assert(eqm(x + y, {148, 122, 81, 55, 122, 75}));
  static_assert(eqm(x - y, {-6, 58, -53, 45, -40, 55}));
  static_assert(eqm(3 * x, {213, 270, 42, 150, 123, 195}));
  static_assert(eqm(x * 3, {213, 270, 42, 150, 123, 195}));
}
} // namespace Arith

File addition: mat.hxx (-rw-r--r--)

[3.425]

#ifndef MAT_HXX
#define MAT_HXX
#include "vect.hxx"
// Only for dnup_t
#include "vec.hxx"
#include <cctk.h>
#include <algorithm>
#include <array>
#include <cassert>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <tuple>
#include <utility>
#include <vector>
#ifdef CCTK_DEBUG
#define ARITH_INLINE
#else
#define ARITH_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE
#endif
namespace Arith {
// Symmetric 3-matrix
template <typename T, int D, dnup_t dnup1, dnup_t dnup2> class mat {
  static_assert(dnup1 == dnup2, "");
  template <typename, int, dnup_t, dnup_t> friend class mat;
  constexpr static int N = D * (D - 1);
  vect<T, N> elts;
  static constexpr ARITH_INLINE int symind(const int i, const int j) {
#ifdef CCTK_DEBUG
    assert(i >= 0 && i <= j && j < 3);
#endif
    const int n = 3 * i - i * (i + 1) / 2 + j;
    // i j n
    // 0 0 0
    // 0 1 1
    // 0 2 2
    // 1 1 3
    // 1 2 4
    // 2 2 5
    // const int n = 2 * i + j - (unsigned)i / 2;
#ifdef CCTK_DEBUG
    assert(n >= 0 && n < N);
#endif
    return n;
  }
  static constexpr ARITH_INLINE int ind(const int i, const int j) {
    using std::max, std::min;
    return symind(min(i, j), max(i, j));
  }
  static_assert(symind(0, 0) == 0, "");
  static_assert(symind(0, 1) == 1, "");
  static_assert(symind(0, 2) == 2, "");
  static_assert(symind(1, 1) == 3, "");
  static_assert(symind(1, 2) == 4, "");
  static_assert(symind(2, 2) == 5, "");
  static_assert(ind(1, 0) == ind(0, 1), "");
  static_assert(ind(2, 0) == ind(0, 2), "");
  static_assert(ind(2, 1) == ind(1, 2), "");
public:
  explicit constexpr ARITH_INLINE mat() : elts() {}
  constexpr ARITH_INLINE mat(const mat &) = default;
  constexpr ARITH_INLINE mat(mat &&) = default;
  constexpr ARITH_INLINE mat &operator=(const mat &) = default;
  constexpr ARITH_INLINE mat &operator=(mat &&) = default;
  template <typename U>
  constexpr ARITH_INLINE mat(const mat<U, D, dnup1, dnup2> &x) : elts(x.elts) {}
  template <typename U>
  constexpr ARITH_INLINE mat(mat<U, D, dnup1, dnup2> &&x)
      : elts(move(x.elts)) {}
  constexpr ARITH_INLINE mat(const vect<T, N> &elts) : elts(elts) {}
  constexpr ARITH_INLINE mat(vect<T, N> &&elts) : elts(move(elts)) {}
  constexpr ARITH_INLINE mat(initializer_list<T> A) : elts(A) {}
  constexpr ARITH_INLINE mat(const vector<T> &A) : elts(A) {}
  constexpr ARITH_INLINE mat(vector<T> &&A) : elts(move(A)) {}
  template <typename F, typename = result_of_t<F(int, int)> >
  constexpr ARITH_INLINE mat(F f) : mat(iota1().map(f, iota2())) {
    // #ifdef CCTK_DEBUG
    //     // Check symmetry
    //     const auto is_symmetric{[](const T &fgood, const T &fother) {
    //       return norm1<T>()(fother - fgood) <=
    //              1.0e-12 * (1 + norm1<T>()(fgood) + norm1<T>()(fother));
    //     }};
    //     for (int i = 0; i < D; ++i) {
    //       for (j = i + 1; j < D; ++j) {
    //         const T fij = f(i, j);
    //         const T fji = f(j, i);
    //         if (!is_symmetric(fij, fji)) {
    //           ostringstream buf;
    //           buf << "f(0,1)=" << f(0, 1) << "\n"
    //               << "f(1,0)=" << f10 << "\n"
    //               << "f(0,2)=" << f(0, 2) << "\n"
    //               << "f(2,0)=" << f20 << "\n"
    //               << "f(1,2)=" << f(1, 2) << "\n"
    //               << "f(2,1)=" << f21 << "\n";
    //           CCTK_VERROR("symmetric matrix is not symmetric:\n%s",
    //                       buf.str().c_str());
    //         }
    //         assert(is_symmetric(fij));
    //       }
    //     }
    // #endif
  }
  static constexpr ARITH_INLINE mat unit(int i, int j) {
    mat r;
    r(i, j) = 1;
    return r;
  }
  static constexpr ARITH_INLINE mat<array<int, 2>, D, dnup1, dnup2> iota() {
    mat<array<int, 2>, D, dnup1, dnup2> r;
    for (int i = 0; i < D; ++i)
      for (int j = i; j < D; ++j)
        r(i, j) = {i, j};
    return r;
  }
  static constexpr ARITH_INLINE mat<int, D, dnup1, dnup2> iota1() {
    mat<int, D, dnup1, dnup2> r;
    for (int i = 0; i < D; ++i)
      for (int j = i; j < D; ++j)
        r(i, j) = i;
    return r;
  }
  static constexpr ARITH_INLINE mat<int, D, dnup1, dnup2> iota2() {
    mat<int, D, dnup1, dnup2> r;
    for (int i = 0; i < D; ++i)
      for (int j = i; j < D; ++j)
        r(i, j) = j;
    return r;
  }
  template <typename F,
            typename R = remove_cv_t<remove_reference_t<result_of_t<F(T)> > > >
  constexpr ARITH_INLINE mat<R, D, dnup1, dnup2> map(F f) const {
    mat<R, D, dnup1, dnup2> r;
    for (int i = 0; i < N; ++i)
      r.elts[i] = f(elts[i]);
    return r;
  }
  template <
      typename F, typename U,
      typename R = remove_cv_t<remove_reference_t<result_of_t<F(T, U)> > > >
  constexpr ARITH_INLINE mat<R, D, dnup1, dnup2>
  map(F f, const mat<U, D, dnup1, dnup2> &x) const {
    mat<R, D, dnup1, dnup2> r;
    for (int i = 0; i < N; ++i)
      r.elts[i] = f(elts[i], x.elts[i]);
    return r;
  }
  constexpr ARITH_INLINE const T &operator()(int i, int j) const {
    return elts[ind(i, j)];
  }
  constexpr ARITH_INLINE T &operator()(int i, int j) { return elts[ind(i, j)]; }
  // template <typename U = T>
  // ARITH_INLINE
  //     mat3<remove_cv_t<remove_reference_t<result_of_t<U(vect<int, 3>)> > >,
  //          dnup1, dnup2>
  //     operator()(const vect<int, 3> &I) const {
  //   return {elts[0](I), elts[1](I), elts[2](I),
  //           elts[3](I), elts[4](I), elts[5](I)};
  // }
  friend constexpr ARITH_INLINE mat<T, D, dnup1, dnup2>
  operator+(const mat<T, D, dnup1, dnup2> &x) {
    return {+x.elts};
  }
  friend constexpr ARITH_INLINE mat<T, D, dnup1, dnup2>
  operator-(const mat<T, D, dnup1, dnup2> &x) {
    return {-x.elts};
  }
  friend constexpr ARITH_INLINE mat<T, D, dnup1, dnup2>
  operator+(const mat<T, D, dnup1, dnup2> &x,
            const mat<T, D, dnup1, dnup2> &y) {
    return {x.elts + y.elts};
  }
  friend constexpr ARITH_INLINE mat<T, D, dnup1, dnup2>
  operator-(const mat<T, D, dnup1, dnup2> &x,
            const mat<T, D, dnup1, dnup2> &y) {
    return {x.elts - y.elts};
  }
  friend constexpr ARITH_INLINE mat<T, D, dnup1, dnup2>
  operator*(const T &a, const mat<T, D, dnup1, dnup2> &x) {
    return {a * x.elts};
  }
  friend constexpr ARITH_INLINE mat<T, D, dnup1, dnup2>
  operator*(const mat<T, D, dnup1, dnup2> &x, const T &a) {
    return {x.elts * a};
  }
  constexpr ARITH_INLINE mat operator+=(const mat &x) {
    return *this = *this + x;
  }
  constexpr ARITH_INLINE mat operator-=(const mat &x) {
    return *this = *this - x;
  }
  constexpr ARITH_INLINE mat operator*=(const T &a) {
    return *this = *this * a;
  }
  constexpr ARITH_INLINE mat operator/=(const T &a) {
    return *this = *this / a;
  }
  friend constexpr ARITH_INLINE bool
  operator==(const mat<T, D, dnup1, dnup2> &x,
             const mat<T, D, dnup1, dnup2> &y) {
    return equal_to<vect<T, N> >()(x.elts, y.elts);
  }
  friend constexpr ARITH_INLINE bool
  operator!=(const mat<T, D, dnup1, dnup2> &x,
             const mat<T, D, dnup1, dnup2> &y) {
    return !(x == y);
  }
  constexpr ARITH_INLINE T maxabs() const { return elts.maxabs(); }
  // constexpr ARITH_INLINE T det() const {
  //   const auto &A = *this;
  //   return A(0, 0) * (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)) -
  //          A(1, 0) * (A(0, 1) * A(2, 2) - A(0, 2) * A(2, 1)) +
  //          A(2, 0) * (A(0, 1) * A(1, 2) - A(0, 2) * A(1, 1));
  // }
  // constexpr ARITH_INLINE mat3<T, !dnup1, !dnup2> inv(const T detA) const {
  //   const auto &A = *this;
  //   const T detA1 = 1 / detA;
  //   return mat3<T, !dnup1, !dnup2>{
  //       detA1 * (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)),
  //       detA1 * (A(1, 2) * A(2, 0) - A(1, 0) * A(2, 2)),
  //       detA1 * (A(1, 0) * A(2, 1) - A(1, 1) * A(2, 0)),
  //       detA1 * (A(2, 2) * A(0, 0) - A(2, 0) * A(0, 2)),
  //       detA1 * (A(2, 0) * A(0, 1) - A(2, 1) * A(0, 0)),
  //       detA1 * (A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0))};
  // }
  // constexpr ARITH_INLINE T trace(const mat3<T, !dnup1, !dnup2> &gu) const {
  //   const auto &A = *this;
  //   return sum2([&](int x, int y) ARITH_INLINE { return gu(x, y) * A(x, y);
  //   });
  // }
  // constexpr ARITH_INLINE mat3 trace_free(
  //     const mat<T, D,dnup1, dnup2> &g, const mat3<T, !dnup1, !dnup2> &gu)
  //     const {
  //   const auto &A = *this;
  //   const T trA = A.trace(gu);
  //   return mat3([&](int a, int b)
  //                   ARITH_INLINE { return A(a, b) - trA / 3 * g(a, b); });
  // }
  friend ostream &operator<<(ostream &os, const mat<T, D, dnup1, dnup2> &A) {
    os << "(" << dnup1 << dnup2 << ")[";
    for (int j = 0; j < D; ++j) {
      if (j > 0)
        os << ",";
      os << "[";
      for (int i = 0; i < D; ++i) {
        if (i > 0)
          os << ",";
        os << A(i, j);
      }
      os << "]";
    }
    os << "]";
    return os;
  }
}; // namespace Arith
} // namespace Arith
#endif // #ifndef MAT_HXX

File addition: sum.cxx (-rw-r--r--)

[3.425]

#include "sum.hxx"
#include <cctk.h>
#include <functional>
namespace Arith {
using namespace std;
namespace {
constexpr CCTK_REAL f0(int i) { return 0; }
constexpr CCTK_REAL f1(int i) { return 1; }
constexpr CCTK_REAL fi(int i) { return i; }
constexpr CCTK_REAL g0(int i, int j) { return 0; }
constexpr CCTK_REAL g1(int i, int j) { return 1; }
constexpr CCTK_REAL gi(int i, int j) { return i; }
constexpr CCTK_REAL gj(int i, int j) { return j; }
constexpr CCTK_REAL gij(int i, int j) { return i + 10 * j; }
constexpr CCTK_REAL h0(int i, int j, int k) { return 0; }
constexpr CCTK_REAL h1(int i, int j, int k) { return 1; }
constexpr CCTK_REAL hi(int i, int j, int k) { return i; }
constexpr CCTK_REAL hj(int i, int j, int k) { return j; }
constexpr CCTK_REAL hk(int i, int j, int k) { return j; }
constexpr CCTK_REAL hijk(int i, int j, int k) { return i + 10 * j + 100 * k; }
} // namespace
// This function is compiled, but not executed. The tests are "run" at
// compile time. If this function compiles, the tests pass.
void TestSum() {
  constexpr equal_to<CCTK_REAL> eq;
  static_assert(eq(sum<3>(f0), 0.0));
  static_assert(eq(sum<3>(f1), 3.0));
  static_assert(eq(sum<3>(fi), 3.0));
  static_assert(eq(sum<4>(f0), 0.0));
  static_assert(eq(sum<4>(f1), 4.0));
  static_assert(eq(sum<4>(fi), 6.0));
  static_assert(eq(sum<3>(g0), 0.0));
  static_assert(eq(sum<3>(g1), 9.0));
  static_assert(eq(sum<3>(gi), 9.0));
  static_assert(eq(sum<3>(gj), 9.0));
  static_assert(eq(sum<3>(gij), 99.0));
  static_assert(eq(sum<3>(h0), 0.0));
  static_assert(eq(sum<3>(h1), 27.0));
  static_assert(eq(sum<3>(hi), 27.0));
  static_assert(eq(sum<3>(hj), 27.0));
  static_assert(eq(sum<3>(hk), 27.0));
  static_assert(eq(sum<3>(hijk), 2997.0));
}
} // namespace Arith

File addition: sum.hxx (-rw-r--r--)

[3.425]

#ifndef SUM_HXX
#define SUM_HXX
#include <fixmath.hxx> // include this before <cctk.h>
#include <cctk.h>
#include <functional>
#ifdef CCTK_DEBUG
#define ARITH_INLINE
#else
#define ARITH_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE
#endif
namespace Arith {
using namespace std;
template <int D, typename F>
constexpr ARITH_INLINE remove_cv_t<remove_reference_t<result_of_t<F(int)> > >
sum(F f) {
  using R = remove_cv_t<remove_reference_t<result_of_t<F(int)> > >;
  R s{0};
  for (int x = 0; x < D; ++x)
    s += f(x);
  return s;
}
template <int D, typename F>
constexpr
    ARITH_INLINE remove_cv_t<remove_reference_t<result_of_t<F(int, int)> > >
    sum(F f) {
  using R = remove_cv_t<remove_reference_t<result_of_t<F(int, int)> > >;
  R s{0};
  for (int x = 0; x < D; ++x)
    for (int y = 0; y < D; ++y)
      s += f(x, y);
  return s;
}
template <int D, typename F>
constexpr ARITH_INLINE
    remove_cv_t<remove_reference_t<result_of_t<F(int, int, int)> > >
    sum(F f) {
  using R = remove_cv_t<remove_reference_t<result_of_t<F(int, int, int)> > >;
  R s{0};
  for (int x = 0; x < D; ++x)
    for (int y = 0; y < D; ++y)
      for (int z = 0; z < D; ++z)
        s += f(x, y, z);
  return s;
}
} // namespace Arith
#undef ARITH_INLINE
#endif

File addition: vec.cxx (-rw-r--r--)

[3.425]

#include "vec.hxx"
#include <cctk.h>
#include <functional>
namespace Arith {
using namespace std;
// This function is compiled, but not executed. The tests are "run" at
// compile time. If this function compiles, the tests pass.
void TestVec() {
  constexpr equal_to<CCTK_REAL> eq;
  using V3D = vec<CCTK_REAL, 3, DN>;
  constexpr equal_to<V3D> eqv;
  static_assert(eq(vec<CCTK_REAL, 3, DN>()(0), 0));
  static_assert(eq(vec<CCTK_REAL, 3, DN>()(1), 0));
  static_assert(eq(vec<CCTK_REAL, 3, DN>()(2), 0));
  static_assert(eq(vec<CCTK_REAL, 3, UP>()(0), 0));
  static_assert(eq(vec<CCTK_REAL, 3, UP>()(1), 0));
  static_assert(eq(vec<CCTK_REAL, 3, UP>()(2), 0));
  static_assert(eq(vec<CCTK_REAL, 3, DN>{1, 2, 3}(0), 1));
  static_assert(eq(vec<CCTK_REAL, 3, DN>{1, 2, 3}(1), 2));
  static_assert(eq(vec<CCTK_REAL, 3, DN>{1, 2, 3}(2), 3));
  static_assert(eqv(vec<CCTK_REAL, 3, DN>::iota(), {0, 1, 2}));
  static_assert(eqv(vec<CCTK_REAL, 3, DN>::unit(0), {1, 0, 0}));
  static_assert(eqv(vec<CCTK_REAL, 3, DN>::unit(1), {0, 1, 0}));
  static_assert(eqv(vec<CCTK_REAL, 3, DN>::unit(2), {0, 0, 1}));
  constexpr V3D x = {5, 6, 8};
  static_assert(eqv(+x, {5, 6, 8}));
  static_assert(eqv(-x, {-5, -6, -8}));
  constexpr V3D y = {10, 4, 6};
  static_assert(eqv(x + y, {15, 10, 14}));
  static_assert(eqv(x - y, {-5, 2, 2}));
  static_assert(eqv(3 * x, {15, 18, 24}));
  static_assert(eqv(x * 3, {15, 18, 24}));
}
} // namespace Arith

File addition: vec.hxx (-rw-r--r--)

[3.425]

#ifndef VEC_HXX
#define VEC_HXX
#include "vect.hxx"
#include <cctk.h>
#include <array>
#include <cassert>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <utility>
#include <vector>
#ifdef CCTK_DEBUG
#define ARITH_INLINE
#else
#define ARITH_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE
#endif
namespace Arith {
using namespace std;
enum class dnup_t : bool { dn, up };
constexpr dnup_t DN = dnup_t::dn;
constexpr dnup_t UP = dnup_t::up;
constexpr dnup_t operator!(const dnup_t dnup) { return dnup_t(!bool(dnup)); }
inline ostream &operator<<(ostream &os, const dnup_t dnup) {
  return os << (dnup == DN ? "d" : "u");
}
// vector
template <typename T, int D, dnup_t dnup> class vec {
  template <typename, int, dnup_t> friend class vec;
  vect<T, D> elts;
  static constexpr ARITH_INLINE int ind(const int n) {
#ifdef CCTK_DEBUG
    assert(n >= 0 && n < D);
#endif
    return n;
  }
public:
  // initializes all elements to zero
  constexpr ARITH_INLINE vec() : elts() {}
  constexpr ARITH_INLINE vec(const vec &) = default;
  constexpr ARITH_INLINE vec(vec &&) = default;
  constexpr ARITH_INLINE vec &operator=(const vec &) = default;
  constexpr ARITH_INLINE vec &operator=(vec &&) = default;
  template <typename U>
  constexpr ARITH_INLINE vec(const vec<U, D, dnup> &x) : elts(x.elts) {}
  template <typename U>
  constexpr ARITH_INLINE vec(vec<U, D, dnup> &&x) : elts(move(x.elts)) {}
  constexpr ARITH_INLINE vec(const vect<T, D> &elts) : elts(elts) {}
  constexpr ARITH_INLINE vec(vect<T, D> &&elts) : elts(move(elts)) {}
  constexpr ARITH_INLINE vec(initializer_list<T> v) : elts(v) {}
  constexpr ARITH_INLINE vec(const array<T, D> &v) : elts(v) {}
  constexpr ARITH_INLINE vec(array<T, D> &&v) : elts(move(v)) {}
  constexpr ARITH_INLINE vec(const vector<T> &v) : elts(v) {}
  constexpr ARITH_INLINE vec(vector<T> &&v) : elts(move(v)) {}
  template <typename F, typename = result_of_t<F(int)> >
  constexpr ARITH_INLINE vec(F f) : vec(iota().map(f)) {}
  static constexpr ARITH_INLINE vec unit(int i) {
    vec r;
    r(i) = 1;
    return r;
  }
  static constexpr ARITH_INLINE vec<int, D, dnup> iota() {
    vec<int, D, dnup> r;
    for (int i = 0; i < D; ++i)
      r(i) = i;
    return r;
  }
  template <typename F,
            typename R = remove_cv_t<remove_reference_t<result_of_t<F(T)> > > >
  constexpr ARITH_INLINE vec<R, D, dnup> map(F f) const {
    return vec<R, D, dnup>(elts.map(f));
  }
  template <
      typename F, typename U,
      typename R = remove_cv_t<remove_reference_t<result_of_t<F(T, U)> > > >
  constexpr ARITH_INLINE vec<R, D, dnup> map(F f,
                                             const vec<U, D, dnup> &x) const {
    return vec<R, D, dnup>(elts.map(f, x.elt));
  }
  constexpr ARITH_INLINE const T &operator()(int i) const {
    return elts[ind(i)];
  }
  constexpr ARITH_INLINE T &operator()(int i) { return elts[ind(i)]; }
  // template <typename U = T>
  // ARITH_INLINE
  //     vec3<remove_cv_t<remove_reference_t<result_of_t<U(vect<int, 3>)> > >,
  //          dnup>
  //     operator()(const vect<int, 3> &I) const {
  //   return {elts[0](I), elts[1](I), elts[2](I)};
  // }
  friend constexpr ARITH_INLINE vec<T, D, dnup>
  operator+(const vec<T, D, dnup> &x) {
    return {+x.elts};
  }
  friend constexpr ARITH_INLINE vec<T, D, dnup>
  operator-(const vec<T, D, dnup> &x) {
    return {-x.elts};
  }
  friend constexpr ARITH_INLINE vec<T, D, dnup>
  operator+(const vec<T, D, dnup> &x, const vec<T, D, dnup> &y) {
    return {x.elts + y.elts};
  }
  friend constexpr ARITH_INLINE vec<T, D, dnup>
  operator-(const vec<T, D, dnup> &x, const vec<T, D, dnup> &y) {
    return {x.elts - y.elts};
  }
  friend constexpr ARITH_INLINE vec<T, D, dnup>
  operator*(const T &a, const vec<T, D, dnup> &x) {
    return {a * x.elts};
  }
  friend constexpr ARITH_INLINE vec<T, D, dnup>
  operator*(const vec<T, D, dnup> &x, const T &a) {
    return {x.elts * a};
  }
  constexpr ARITH_INLINE vec operator+=(const vec &x) {
    return *this = *this + x;
  }
  constexpr ARITH_INLINE vec operator-=(const vec &x) {
    return *this = *this - x;
  }
  constexpr ARITH_INLINE vec operator*=(const T &a) {
    return *this = *this * a;
  }
  constexpr ARITH_INLINE vec operator/=(const T &a) {
    return *this = *this / a;
  }
  friend constexpr ARITH_INLINE bool operator==(const vec<T, D, dnup> &x,
                                                const vec<T, D, dnup> &y) {
    return equal_to<vect<T, D> >()(x.elts, y.elts);
  }
  friend constexpr ARITH_INLINE bool operator!=(const vec<T, D, dnup> &x,
                                                const vec<T, D, dnup> &y) {
    return !(x == y);
  }
  constexpr ARITH_INLINE T maxabs() const { return elts.maxabs(); }
  friend ostream &operator<<(ostream &os, const vec<T, D, dnup> &v) {
    os << "(" << dnup << ")[";
    for (int i = 0; i < D; ++i) {
      if (i > 0)
        os << ",";
      os << v(i);
    }
    os << "]";
    return os;
  }
};
} // namespace Arith
#undef ARITH_INLINE
#endif // #ifndef VEC_HXX

Insertion in Arith/src/vect.hxx at line 189 [3.2989]

[3.7228]


  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(const vect &) = default;
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(vect &&) = default;
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect &
  operator=(const vect &) = default;
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect &operator=(vect &&) = default;

Insertion in Arith/src/vect.hxx at line 196 [3.2989]

[3.7229]

  template <typename U>
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(const vect<U, D> &x) : elts() {
    for (int d = 0; d < D; ++d)
      elts[d] = x.elts[d];
  }
  template <typename U>
  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(vect &&x) : elts() {
    for (int d = 0; d < D; ++d)
      elts[d] = move(x.elts[d]);
  }

Replacement in Arith/src/vect.hxx at line 265 [3.2989]

B:BD[3.8501] → [3.8501:8568]

  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R map(const F &f) const {

[3.8501]

[3.8568]

  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<R, D> map(const F &f) const {

Replacement in Arith/src/vect.hxx at line 274 [3.2989]

B:BD[3.8786] → [3.8786:8923]

  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R map2(const F &f,
                                                const vect<U, D> &y) const {

[3.8786]

[3.8923]

  constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<R, D>
  map2(const F &f, const vect<U, D> &y) const {

Arith: New classes vec, mat; new function sum

Dependencies

In channels

Change contents

Insertion in Arith/interface.ccl at line 7 [3.172]

Replacement in Arith/src/make.code.defn at line 4 [4.6031]

File addition: mat.cxx (-rw-r--r--)

File addition: mat.hxx (-rw-r--r--)

File addition: sum.cxx (-rw-r--r--)

File addition: sum.hxx (-rw-r--r--)

File addition: vec.cxx (-rw-r--r--)

File addition: vec.hxx (-rw-r--r--)

Insertion in Arith/src/vect.hxx at line 189 [3.2989]

Insertion in Arith/src/vect.hxx at line 196 [3.2989]

Replacement in Arith/src/vect.hxx at line 265 [3.2989]

Replacement in Arith/src/vect.hxx at line 274 [3.2989]