SymmLTensor.hpp

Go to the documentation of this file.

/**
 * @file SymmLTensor.hpp
 *
 */
 
#pragma once
 
namespace FTensor {
 
/**
 * @brief Mapping from symmetric tensor indices to packed storage index
 *
 * @tparam T
 * @tparam Dim
 */
template <class T = int, int Dim = 3> class SymmLTensor {
public:
  using value_type = typename promote<T, double>::V;
  static constexpr int size_symm = (Dim * (Dim + 1)) / 2;
 
  constexpr value_type operator()(const int N0, const int N1,
                                  const int N2) const {
    const auto hi = N0 > N1 ? N0 : N1;
    const auto lo = N0 > N1 ? N1 : N0;
    const auto packed = hi + (lo * (2 * Dim - lo - 1)) / 2;
    return N2 == packed ? value_type(1) : value_type(0);
  }
 
  template <char i, char j, char k, int Dim01, int Dim2>
  typename std::enable_if<
      (Dim01 == Dim && Dim2 == size_symm),
      Dg_Expr<SymmLTensor<T, Dim>, value_type, Dim01, Dim2, i, j, k>>::type
  operator()(const Index<i, Dim01> &, const Index<j, Dim01> &,
             const Index<k, Dim2> &) const {
    return Dg_Expr<SymmLTensor<T, Dim>, value_type, Dim01, Dim2, i, j, k>(
        *this);
  };
};
 
template <class T = int, int Dim> constexpr auto symm_l_tensor() {
  return SymmLTensor<T, Dim>();
}
 
template <class T = int, int Dim>
constexpr auto symm_l_tensor(const Number<Dim> &) {
  return SymmLTensor<T, Dim>();
}
 
template <class T = int, char i, char j, char k, int Dim01, int Dim2>
typename std::enable_if<
    (Dim2 == (Dim01 * (Dim01 + 1)) / 2),
    Dg_Expr<SymmLTensor<T, Dim01>, typename SymmLTensor<T, Dim01>::value_type,
            Dim01, Dim2, i, j, k>>::type
symm_l_tensor(const Index<i, Dim01> &, const Index<j, Dim01> &,
              const Index<k, Dim2> &) {
  return SymmLTensor<T, Dim01>()(Index<i, Dim01>(), Index<j, Dim01>(),
                                 Index<k, Dim2>());
}
 
} // namespace FTensor