diff --git a/CMakeLists.txt b/CMakeLists.txt index cbea4cef92..b0eec97af9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -129,6 +129,8 @@ option(SEQUANT_BTAS "Enable BTAS eval backend" ${SEQUANT_EVAL_TESTS}) add_feature_info(SEQUANT_EVAL_BTAS SEQUANT_BTAS "Enable BTAS eval backend") option(SEQUANT_TILEDARRAY "Enable TiledArray eval backend" ${SEQUANT_EVAL_TESTS}) add_feature_info(SEQUANT_EVAL_TILEDARRAY SEQUANT_TILEDARRAY "Enable TiledArray eval backend") +option(SEQUANT_TAPP "Enable TAPP eval backend" ${SEQUANT_EVAL_TESTS}) +add_feature_info(SEQUANT_EVAL_TAPP SEQUANT_TAPP "Enable TAPP eval backend") option(SEQUANT_IWYU "Whether to use the include-what-you-use tool (if found)" OFF) option(SEQUANT_WARNINGS_AS_ERRORS "Whether to treat compiler warnings as errors" ${PROJECT_IS_TOP_LEVEL}) @@ -460,6 +462,14 @@ set(SeQuant_eval_btas_src SeQuant/core/eval/backends/btas/result.hpp ) +# TAPP backend sources (header-only) +set(SeQuant_eval_tapp_src + SeQuant/core/eval/backends/tapp/tensor.hpp + SeQuant/core/eval/backends/tapp/ops.hpp + SeQuant/core/eval/backends/tapp/eval_expr.hpp + SeQuant/core/eval/backends/tapp/result.hpp +) + # Core sources (just version info) set(SeQuant_core_src ${PROJECT_BINARY_DIR}/SeQuant/version.hpp @@ -486,6 +496,12 @@ if (SEQUANT_BTAS) endif() set(SEQUANT_HAS_BTAS ON) endif () +if (SEQUANT_TAPP) + if (NOT TARGET tapp::reference) + include(FindOrFetchTAPP) + endif() + set(SEQUANT_HAS_TAPP ON) +endif () if (NOT PROJECT_IS_TOP_LEVEL) set(SEQUANT_SYSTEM "SYSTEM") @@ -578,6 +594,17 @@ if (SEQUANT_HAS_BTAS) target_compile_definitions(SeQuant-eval-btas INTERFACE SEQUANT_HAS_BTAS=1) endif() +########################## +# SeQuant-eval-tapp: TAPP backend (header-only) +########################## +if (SEQUANT_HAS_TAPP) + add_library(SeQuant-eval-tapp INTERFACE) + set_target_properties(SeQuant-eval-tapp PROPERTIES EXPORT_NAME eval::tapp) + add_library(SeQuant::eval::tapp ALIAS SeQuant-eval-tapp) + target_link_libraries(SeQuant-eval-tapp INTERFACE SeQuant-eval SeQuant-symb tapp::reference) + target_compile_definitions(SeQuant-eval-tapp INTERFACE SEQUANT_HAS_TAPP=1) +endif() + ########################## # SeQuant-export: code generation/export (depends on eval) ########################## @@ -616,6 +643,9 @@ endif() if (SEQUANT_HAS_BTAS) list(APPEND SEQUANT_MODULES SeQuant-eval-btas) endif() +if (SEQUANT_HAS_TAPP) + list(APPEND SEQUANT_MODULES SeQuant-eval-tapp) +endif() ########################## # SeQuant: umbrella INTERFACE target (backwards compatibility) diff --git a/SeQuant/core/eval/backends/btas/result.hpp b/SeQuant/core/eval/backends/btas/result.hpp index 5762141236..aa7cff67ed 100644 --- a/SeQuant/core/eval/backends/btas/result.hpp +++ b/SeQuant/core/eval/backends/btas/result.hpp @@ -210,7 +210,6 @@ class ResultTensorBTAS final : public Result { private: [[nodiscard]] std::size_t size_in_bytes() const final { - static_assert(std::is_arithmetic_v); const auto& tensor = get(); // only count data return tensor.range().volume() * sizeof(T); diff --git a/SeQuant/core/eval/backends/tapp/eval_expr.hpp b/SeQuant/core/eval/backends/tapp/eval_expr.hpp new file mode 100644 index 0000000000..b21f2a7218 --- /dev/null +++ b/SeQuant/core/eval/backends/tapp/eval_expr.hpp @@ -0,0 +1,64 @@ +#ifndef SEQUANT_EVAL_BACKENDS_TAPP_EVAL_EXPR_HPP +#define SEQUANT_EVAL_BACKENDS_TAPP_EVAL_EXPR_HPP + +#ifdef SEQUANT_HAS_TAPP + +#include + +#include +#include +#include + +#include + +#include + +namespace sequant { + +/// +/// \brief This class extends the EvalExpr class by adding an annot() method so +/// that it can be used to evaluate using TAPP. +/// +class EvalExprTAPP final : public EvalExpr { + public: + using annot_t = container::svector; + + /// + /// \param bk iterable of Index objects. + /// \return vector of int64_t-type hash values + /// of the labels of indices in \c bk + /// + template + static auto index_hash(Iterable&& bk) { + return ranges::views::transform( + std::forward(bk), [](auto const& idx) { + return static_cast(sequant::hash::value(Index{idx}.label())); + }); + } + + template >> + EvalExprTAPP(Args&&... args) : EvalExpr{std::forward(args)...} { + annot_ = index_hash(canon_indices()) | ranges::to; + } + + /// + /// \return Annotation (container::svector) for TAPP tensors. + /// + [[nodiscard]] inline annot_t const& annot() const noexcept { return annot_; } + + private: + annot_t annot_; +}; + +/// Type alias for TAPP evaluation nodes +using EvalNodeTAPP = EvalNode; + +static_assert(meta::eval_node); +static_assert(meta::can_evaluate); + +} // namespace sequant + +#endif // SEQUANT_HAS_TAPP + +#endif // SEQUANT_EVAL_BACKENDS_TAPP_EVAL_EXPR_HPP diff --git a/SeQuant/core/eval/backends/tapp/ops.hpp b/SeQuant/core/eval/backends/tapp/ops.hpp new file mode 100644 index 0000000000..017b1b7b93 --- /dev/null +++ b/SeQuant/core/eval/backends/tapp/ops.hpp @@ -0,0 +1,267 @@ +#ifndef SEQUANT_EVAL_BACKENDS_TAPP_OPS_HPP +#define SEQUANT_EVAL_BACKENDS_TAPP_OPS_HPP + +#ifdef SEQUANT_HAS_TAPP + +#include + +#include + +#include +#include +#include +#include + +namespace sequant::tapp_ops { + +namespace detail { + +/// RAII wrapper for TAPP_handle. +struct HandleGuard { + TAPP_handle handle = 0; + HandleGuard() { tapp_detail::check_error(TAPP_create_handle(&handle)); } + ~HandleGuard() { + if (handle) TAPP_destroy_handle(handle); + } + HandleGuard(HandleGuard const&) = delete; + HandleGuard& operator=(HandleGuard const&) = delete; +}; + +/// RAII wrapper for TAPP_executor. +struct ExecutorGuard { + TAPP_executor executor = 0; + ExecutorGuard() { tapp_detail::check_error(TAPP_create_executor(&executor)); } + ~ExecutorGuard() { + if (executor) TAPP_destroy_executor(executor); + } + ExecutorGuard(ExecutorGuard const&) = delete; + ExecutorGuard& operator=(ExecutorGuard const&) = delete; +}; + +/// RAII wrapper for TAPP_tensor_product. +struct ProductGuard { + TAPP_tensor_product plan = 0; + ProductGuard() = default; + explicit ProductGuard(TAPP_tensor_product p) : plan{p} {} + ~ProductGuard() { + if (plan) TAPP_destroy_tensor_product(plan); + } + ProductGuard(ProductGuard const&) = delete; + ProductGuard& operator=(ProductGuard const&) = delete; +}; + +/// RAII wrapper for TAPP_tensor_info. +struct TensorInfoGuard { + TAPP_tensor_info info = 0; + TensorInfoGuard() = default; + explicit TensorInfoGuard(TAPP_tensor_info i) : info{i} {} + ~TensorInfoGuard() { + if (info) TAPP_destroy_tensor_info(info); + } + TensorInfoGuard(TensorInfoGuard const&) = delete; + TensorInfoGuard& operator=(TensorInfoGuard const&) = delete; +}; + +/// Create a scalar (0-mode) TAPP_tensor_info. +template +TensorInfoGuard make_scalar_info() { + TAPP_tensor_info info; + tapp_detail::check_error( + TAPP_create_tensor_info(&info, tapp_detail::datatype(), 0, + /*extents=*/nullptr, /*strides=*/nullptr)); + return TensorInfoGuard{info}; +} + +/// Compute the extents of the result tensor from a contraction. +/// Indices appearing in idx_result but not contracted are "free" indices. +/// Their extents come from A or B. +template +container::svector compute_result_extents( + TAPPTensor const& A, container::svector const& idx_A, + TAPPTensor const& B, container::svector const& idx_B, + container::svector const& idx_result) { + container::svector result_extents(idx_result.size()); + for (size_t r = 0; r < idx_result.size(); ++r) { + bool found = false; + for (size_t a = 0; a < idx_A.size(); ++a) { + if (idx_A[a] == idx_result[r]) { + result_extents[r] = A.extents()[a]; + found = true; + break; + } + } + if (!found) { + for (size_t b = 0; b < idx_B.size(); ++b) { + if (idx_B[b] == idx_result[r]) { + result_extents[r] = B.extents()[b]; + found = true; + break; + } + } + } + assert(found && "Result index not found in A or B"); + } + return result_extents; +} + +} // namespace detail + +/// +/// \brief Tensor contraction using the TAPP API. +/// +/// Computes: result[idx_result] = alpha * A[idx_A] * B[idx_B] +/// + beta * result[idx_result] +/// +/// When beta == 0 and the result tensor is empty, it is automatically +/// allocated with the correct shape inferred from the index labels and +/// the operand extents. +/// +template +void contract(T alpha, // + TAPPTensor const& A, // + container::svector const& idx_A, // + TAPPTensor const& B, // + container::svector const& idx_B, // + T beta, // + TAPPTensor& result, // + container::svector const& idx_result) { + // If result is empty and beta==0, allocate it + if (result.volume() == 0) { + auto extents = + detail::compute_result_extents(A, idx_A, B, idx_B, idx_result); + result = TAPPTensor(std::move(extents)); + result.fill(T{0}); + } + + detail::HandleGuard handle; + detail::ExecutorGuard executor; + + // Create the product plan: + // D = alpha * op_A(A) * op_B(B) + beta * op_C(C) + // We set C = D (in-place accumulation when beta != 0) + TAPP_tensor_product plan; + tapp_detail::check_error(TAPP_create_tensor_product( + &plan, handle.handle, // + TAPP_IDENTITY, A.info(), idx_A.data(), // + TAPP_IDENTITY, B.info(), idx_B.data(), // + TAPP_IDENTITY, result.info(), idx_result.data(), // C = result + TAPP_IDENTITY, result.info(), idx_result.data(), // D = result + TAPP_DEFAULT_PREC)); + detail::ProductGuard plan_guard{plan}; + + TAPP_status status = 0; + tapp_detail::check_error(TAPP_execute_product( + plan, executor.executor, &status, // + &alpha, // + A.data(), B.data(), // + &beta, // + beta != T{0} ? result.data() : nullptr, // C data (null if beta==0) + result.data())); // D data +} + +/// +/// \brief Permute a tensor using the TAPP API. +/// +/// Computes: dst[dst_idx] = src[src_idx] +/// where src_idx and dst_idx define the index correspondence. +/// +template +void permute(TAPPTensor const& src, // + container::svector const& src_idx, // + TAPPTensor& dst, // + container::svector const& dst_idx) { + assert(src_idx.size() == static_cast(src.rank())); + assert(dst_idx.size() == static_cast(src.rank())); + + // Determine destination extents from index mapping + container::svector dst_extents(dst_idx.size()); + for (size_t d = 0; d < dst_idx.size(); ++d) { + for (size_t s = 0; s < src_idx.size(); ++s) { + if (src_idx[s] == dst_idx[d]) { + dst_extents[d] = src.extents()[s]; + break; + } + } + } + + dst = TAPPTensor(std::move(dst_extents)); + + // Use TAPP product: D = 1.0 * A * scalar(1.0) + 0.0 * C + // with A's indices reordered to D's layout + detail::HandleGuard handle; + detail::ExecutorGuard executor; + auto scalar_info = detail::make_scalar_info(); + + T one{1}; + T zero{0}; + + TAPP_tensor_product plan; + tapp_detail::check_error(TAPP_create_tensor_product( + &plan, handle.handle, // + TAPP_IDENTITY, src.info(), src_idx.data(), // A + TAPP_IDENTITY, scalar_info.info, nullptr, // B = scalar + TAPP_IDENTITY, dst.info(), dst_idx.data(), // C = D + TAPP_IDENTITY, dst.info(), dst_idx.data(), // D + TAPP_DEFAULT_PREC)); + detail::ProductGuard plan_guard{plan}; + + TAPP_status status = 0; + tapp_detail::check_error( + TAPP_execute_product(plan, executor.executor, &status, // + &one, // alpha + src.data(), &one, // A data, B data (scalar 1) + &zero, // beta + nullptr, // C data (unused, beta=0) + dst.data())); // D data +} + +/// +/// \brief Scale a tensor in place: tensor *= alpha. +/// +template +void scal(T alpha, TAPPTensor& tensor) { + T* d = tensor.data(); + size_t n = tensor.volume(); + for (size_t i = 0; i < n; ++i) d[i] *= alpha; +} + +/// +/// \brief Dot product of two tensors (full contraction). +/// +/// \return The sum of element-wise products of A and B. +/// \pre A and B must have the same shape. +/// +template +T dot(TAPPTensor const& A, TAPPTensor const& B) { + assert(A.extents() == B.extents()); + T const* a = A.data(); + T const* b = B.data(); + size_t n = A.volume(); + T result{0}; + for (size_t i = 0; i < n; ++i) result += a[i] * b[i]; + return result; +} + +/// +/// \brief Compute the dot product with complex conjugation: sum(conj(A)*B). +/// +template +T dotc(TAPPTensor const& A, TAPPTensor const& B) { + if constexpr (std::is_arithmetic_v) { + return dot(A, B); + } else { + assert(A.extents() == B.extents()); + T const* a = A.data(); + T const* b = B.data(); + size_t n = A.volume(); + T result{0}; + for (size_t i = 0; i < n; ++i) result += std::conj(a[i]) * b[i]; + return result; + } +} + +} // namespace sequant::tapp_ops + +#endif // SEQUANT_HAS_TAPP + +#endif // SEQUANT_EVAL_BACKENDS_TAPP_OPS_HPP diff --git a/SeQuant/core/eval/backends/tapp/result.hpp b/SeQuant/core/eval/backends/tapp/result.hpp new file mode 100644 index 0000000000..111710a01a --- /dev/null +++ b/SeQuant/core/eval/backends/tapp/result.hpp @@ -0,0 +1,244 @@ +#ifndef SEQUANT_EVAL_BACKENDS_TAPP_RESULT_HPP +#define SEQUANT_EVAL_BACKENDS_TAPP_RESULT_HPP + +#ifdef SEQUANT_HAS_TAPP + +#include +#include +#include +#include + +namespace sequant { + +namespace { + +/// +/// \brief Symmetrize a TAPPTensor by summing over all symmetric particle +/// permutations. +/// +/// \param arr The tensor to be symmetrized. +/// \pre The rank of the tensor must be even. +/// \return The symmetrized TAPPTensor. +/// +template +auto column_symmetrize_tapp(TAPPTensor const& arr) { + using ranges::views::iota; + + size_t const rank = arr.rank(); + + if (rank % 2 != 0) + throw std::domain_error("This function only supports even-ranked tensors"); + + perm_t perm = iota(size_t{0}, rank) | ranges::to; + + // lannot is identity ordering + auto const lannot = perm | ranges::views::transform([](size_t v) { + return static_cast(v); + }) | + ranges::to>; + + auto result = TAPPTensor(arr.extents()); + result.fill(T{0}); + + auto call_back = [&result, &lannot, &arr, &perm = std::as_const(perm)]() { + auto perm_annot = perm | ranges::views::transform([](size_t v) { + return static_cast(v); + }) | + ranges::to>; + TAPPTensor temp; + tapp_ops::permute(arr, lannot, temp, perm_annot); + result += temp; + }; + + auto const nparticles = rank / 2; + symmetric_permutation(SymmetricParticleRange{perm.begin(), // + perm.begin() + nparticles, // + nparticles}, + call_back); + auto const nf = static_cast(rational{1, factorial(nparticles)}); + tapp_ops::scal(static_cast(nf), result); + + return result; +} + +/// +/// \brief Antisymmetrize a TAPPTensor over bra and ket index groups. +/// +/// \param arr The tensor to be antisymmetrized. +/// \param bra_rank The rank of the bra indices. +/// \return The antisymmetrized TAPPTensor. +/// +template +auto particle_antisymmetrize_tapp(TAPPTensor const& arr, + size_t bra_rank) { + using ranges::views::concat; + using ranges::views::iota; + size_t const rank = arr.rank(); + SEQUANT_ASSERT(bra_rank <= rank); + size_t const ket_rank = rank - bra_rank; + + perm_t bra_perm = iota(size_t{0}, bra_rank) | ranges::to; + perm_t ket_perm = iota(bra_rank, rank) | ranges::to; + const auto lannot_perm = iota(size_t{0}, rank) | ranges::to; + + // Convert perm_t (svector) to annot (svector) + auto to_annot = [](perm_t const& p) { + return p | ranges::views::transform([](size_t v) { + return static_cast(v); + }) | + ranges::to>; + }; + + auto const lannot = to_annot(lannot_perm); + + auto process_permutations = [&lannot, &to_annot]( + const TAPPTensor& input_arr, + size_t range_rank, perm_t range_perm, + const perm_t& other_perm, bool is_bra) { + if (range_rank <= 1) return input_arr; + TAPPTensor result{input_arr.extents()}; + result.fill(T{0}); + + auto callback = [&](int parity) { + const auto annot_perm = + is_bra ? concat(range_perm, other_perm) | ranges::to() + : concat(other_perm, range_perm) | ranges::to(); + + auto const annot = to_annot(annot_perm); + + T p_ = parity == 0 ? T{1} : T{-1}; + TAPPTensor temp; + tapp_ops::permute(input_arr, lannot, temp, annot); + tapp_ops::scal(p_, temp); + result += temp; + }; + + antisymmetric_permutation(ParticleRange{range_perm.begin(), range_rank}, + callback); + return result; + }; + + // Process bra permutations first + const auto ket_annot = ket_rank == 0 ? perm_t{} : ket_perm; + auto result = process_permutations(arr, bra_rank, bra_perm, ket_annot, true); + + // Process ket permutations if needed + const auto bra_annot = bra_rank == 0 ? perm_t{} : bra_perm; + result = process_permutations(result, ket_rank, ket_perm, bra_annot, false); + + auto const nf = static_cast( + rational{1, factorial(bra_rank) * factorial(ket_rank)}); + tapp_ops::scal(static_cast(nf), result); + + return result; +} + +} // namespace + +/// +/// \brief Result for a tensor value of TAPPTensor type. +/// \tparam T TAPPTensor type. Must be a specialization of TAPPTensor. +/// +template +class ResultTensorTAPP final : public Result { + public: + using Result::id_t; + using numeric_type = typename T::numeric_type; + + explicit ResultTensorTAPP(T arr) : Result{std::move(arr)} {} + + private: + using annot_t = container::svector; + using annot_wrap = Annot; + + [[nodiscard]] id_t type_id() const noexcept override { + return id_for_type>(); + } + + [[nodiscard]] ResultPtr sum( + Result const& other, + std::array const& annot) const override { + SEQUANT_ASSERT(other.is>()); + auto const a = annot_wrap{annot}; + + T lres, rres; + tapp_ops::permute(get(), a.lannot, lres, a.this_annot); + tapp_ops::permute(other.get(), a.rannot, rres, a.this_annot); + lres += rres; + return eval_result>(std::move(lres)); + } + + [[nodiscard]] ResultPtr prod(Result const& other, + std::array const& annot, + DeNest /*DeNestFlag*/) const override { + auto const a = annot_wrap{annot}; + + if (other.is>()) { + T result; + tapp_ops::permute(get(), a.lannot, result, a.this_annot); + tapp_ops::scal(other.as>().value(), result); + return eval_result>(std::move(result)); + } + + SEQUANT_ASSERT(other.is>()); + + if (a.this_annot.empty()) { + // Full contraction -> scalar result (dot product) + T rres; + tapp_ops::permute(other.get(), a.rannot, rres, a.lannot); + return eval_result>( + tapp_ops::dot(get(), rres)); + } + + T result; + tapp_ops::contract(numeric_type{1}, // + get(), a.lannot, // + other.get(), a.rannot, // + numeric_type{0}, // + result, a.this_annot); + return eval_result>(std::move(result)); + } + + [[nodiscard]] ResultPtr mult_by_phase(std::int8_t factor) const override { + auto pre = get(); + tapp_ops::scal(numeric_type(factor), pre); + return eval_result>(std::move(pre)); + } + + [[nodiscard]] ResultPtr permute( + std::array const& ann) const override { + auto const pre_annot = std::any_cast(ann[0]); + auto const post_annot = std::any_cast(ann[1]); + T result; + tapp_ops::permute(get(), pre_annot, result, post_annot); + return eval_result>(std::move(result)); + } + + void add_inplace(Result const& other) override { + auto& t = get(); + auto const& o = other.get(); + SEQUANT_ASSERT(t.extents() == o.extents()); + t += o; + } + + [[nodiscard]] ResultPtr symmetrize() const override { + return eval_result>(column_symmetrize_tapp(get())); + } + + [[nodiscard]] ResultPtr antisymmetrize(size_t bra_rank) const override { + return eval_result>( + particle_antisymmetrize_tapp(get(), bra_rank)); + } + + private: + [[nodiscard]] std::size_t size_in_bytes() const final { + const auto& tensor = get(); + return tensor.volume() * sizeof(typename T::value_type); + } +}; + +} // namespace sequant + +#endif // SEQUANT_HAS_TAPP + +#endif // SEQUANT_EVAL_BACKENDS_TAPP_RESULT_HPP diff --git a/SeQuant/core/eval/backends/tapp/tensor.hpp b/SeQuant/core/eval/backends/tapp/tensor.hpp new file mode 100644 index 0000000000..a6a9abd2c2 --- /dev/null +++ b/SeQuant/core/eval/backends/tapp/tensor.hpp @@ -0,0 +1,279 @@ +#ifndef SEQUANT_EVAL_BACKENDS_TAPP_TENSOR_HPP +#define SEQUANT_EVAL_BACKENDS_TAPP_TENSOR_HPP + +#ifdef SEQUANT_HAS_TAPP + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sequant { + +namespace tapp_detail { + +/// Maps C++ numeric types to TAPP_datatype constants. +template +constexpr TAPP_datatype datatype() noexcept { + if constexpr (std::is_same_v) + return TAPP_F32; + else if constexpr (std::is_same_v) + return TAPP_F64; + else if constexpr (std::is_same_v>) + return TAPP_C32; + else if constexpr (std::is_same_v>) + return TAPP_C64; + else + static_assert(!sizeof(T), "Unsupported type for TAPP"); +} + +/// Compute row-major (C-contiguous) strides from extents. +inline container::svector compute_strides( + container::svector const& extents) { + container::svector strides(extents.size()); + if (extents.empty()) return strides; + strides.back() = 1; + for (int i = static_cast(extents.size()) - 2; i >= 0; --i) { + strides[i] = strides[i + 1] * extents[i + 1]; + } + return strides; +} + +/// Compute the volume (total number of elements) from extents. +inline size_t compute_volume(container::svector const& extents) { + if (extents.empty()) return 1; // scalar + return static_cast(std::accumulate(extents.begin(), extents.end(), + int64_t{1}, std::multiplies<>{})); +} + +/// Check a TAPP error code and throw on failure. +inline void check_error(TAPP_error err) { + if (!TAPP_check_success(err)) { + char msg[256]; + TAPP_explain_error(err, sizeof(msg), msg); + throw std::runtime_error(std::string("TAPP error: ") + msg); + } +} + +} // namespace tapp_detail + +/// +/// \brief A lightweight C++ tensor container backed by the TAPP C API. +/// +/// Owns the data buffer and manages a TAPP_tensor_info opaque handle that +/// describes the tensor's shape and layout for use with TAPP operations. +/// +/// \tparam T The numeric element type (float, double, complex, +/// complex). +/// \tparam Allocator The allocator for the data buffer. +/// +template > +class TAPPTensor { + public: + using numeric_type = T; + using value_type = T; + using allocator_type = Allocator; + + /// Default constructor: creates an empty (scalar-like) tensor. + TAPPTensor() = default; + + /// Construct from extents (row-major layout). + explicit TAPPTensor(container::svector extents, + Allocator alloc = Allocator{}) + : extents_{std::move(extents)}, + strides_{tapp_detail::compute_strides(extents_)}, + data_(tapp_detail::compute_volume(extents_), T{}, std::move(alloc)) { + create_info(); + } + + /// Construct from any range of extents. + template < + typename Range, + typename = std::enable_if_t< + !std::is_same_v, container::svector> && + !std::is_same_v, TAPPTensor>>> + explicit TAPPTensor(Range const& range, Allocator alloc = Allocator{}) + : TAPPTensor( + container::svector(std::begin(range), std::end(range)), + std::move(alloc)) {} + + ~TAPPTensor() { destroy_info(); } + + TAPPTensor(TAPPTensor const& other) + : extents_{other.extents_}, strides_{other.strides_}, data_{other.data_} { + create_info(); + } + + TAPPTensor(TAPPTensor&& other) noexcept + : extents_{std::move(other.extents_)}, + strides_{std::move(other.strides_)}, + data_{std::move(other.data_)}, + info_{other.info_} { + other.info_ = 0; + } + + TAPPTensor& operator=(TAPPTensor const& other) { + if (this != &other) { + destroy_info(); + extents_ = other.extents_; + strides_ = other.strides_; + data_ = other.data_; + create_info(); + } + return *this; + } + + TAPPTensor& operator=(TAPPTensor&& other) noexcept { + if (this != &other) { + destroy_info(); + extents_ = std::move(other.extents_); + strides_ = std::move(other.strides_); + data_ = std::move(other.data_); + info_ = other.info_; + other.info_ = 0; + } + return *this; + } + + /// \return The number of modes (dimensions). + [[nodiscard]] int rank() const noexcept { + return static_cast(extents_.size()); + } + + [[nodiscard]] container::svector const& extents() const noexcept { + return extents_; + } + + [[nodiscard]] container::svector const& strides() const noexcept { + return strides_; + } + + [[nodiscard]] T* data() noexcept { return data_.data(); } + + [[nodiscard]] T const* data() const noexcept { return data_.data(); } + + /// \return The total number of elements. + [[nodiscard]] size_t volume() const noexcept { return data_.size(); } + + /// \return The TAPP_tensor_info opaque handle for this tensor. + [[nodiscard]] TAPP_tensor_info info() const noexcept { return info_; } + + // Iterator support (for compatibility with std::begin/std::end) + using iterator = typename std::vector::iterator; + using const_iterator = typename std::vector::const_iterator; + iterator begin() noexcept { return data_.begin(); } + iterator end() noexcept { return data_.end(); } + const_iterator begin() const noexcept { return data_.begin(); } + const_iterator end() const noexcept { return data_.end(); } + const_iterator cbegin() const noexcept { return data_.cbegin(); } + const_iterator cend() const noexcept { return data_.cend(); } + + /// Fill all elements with the given value. + void fill(T val) { std::fill(data_.begin(), data_.end(), val); } + + /// Generate elements using a callable. + template + void generate(Gen gen) { + std::generate(data_.begin(), data_.end(), gen); + } + + /// Element-wise in-place addition. + TAPPTensor& operator+=(TAPPTensor const& other) { + assert(extents_ == other.extents_); + for (size_t i = 0; i < data_.size(); ++i) data_[i] += other.data_[i]; + return *this; + } + + /// Element-wise subtraction returning a new tensor. + TAPPTensor operator-(TAPPTensor const& other) const { + assert(extents_ == other.extents_); + TAPPTensor result(extents_); + for (size_t i = 0; i < data_.size(); ++i) + result.data_[i] = data_[i] - other.data_[i]; + return result; + } + + /// Variadic element access (row-major). + template + [[nodiscard]] T& operator()(Indices... indices) { + return data_[offset(indices...)]; + } + + template + [[nodiscard]] T const& operator()(Indices... indices) const { + return data_[offset(indices...)]; + } + + /// Element access via index vector. + template + requires(!std::is_integral_v) + [[nodiscard]] T& operator()(Vec const& idx_vec) { + return data_[offset_from_vec(idx_vec)]; + } + + template + requires(!std::is_integral_v) + [[nodiscard]] T const& operator()(Vec const& idx_vec) const { + return data_[offset_from_vec(idx_vec)]; + } + + /// Reset to empty state. + void clear() { + data_.clear(); + extents_.clear(); + strides_.clear(); + destroy_info(); + } + + private: + container::svector extents_; + container::svector strides_; + std::vector data_; + TAPP_tensor_info info_ = 0; + + void create_info() { + tapp_detail::check_error( + TAPP_create_tensor_info(&info_, tapp_detail::datatype(), rank(), + extents_.data(), strides_.data())); + } + + void destroy_info() { + if (info_) { + TAPP_destroy_tensor_info(info_); + info_ = 0; + } + } + + template + [[nodiscard]] size_t offset(Indices... indices) const { + static_assert((std::is_integral_v && ...)); + assert(sizeof...(indices) == extents_.size()); + size_t off = 0; + size_t i = 0; + ((off += static_cast(indices) * strides_[i++]), ...); + return off; + } + + template + [[nodiscard]] size_t offset_from_vec(Vec const& idx_vec) const { + size_t off = 0; + for (size_t i = 0; i < extents_.size(); ++i) + off += static_cast(idx_vec[i]) * strides_[i]; + return off; + } +}; + +} // namespace sequant + +#endif // SEQUANT_HAS_TAPP + +#endif // SEQUANT_EVAL_BACKENDS_TAPP_TENSOR_HPP diff --git a/SeQuant/domain/mbpt/biorthogonalization.hpp b/SeQuant/domain/mbpt/biorthogonalization.hpp index 5593b2d5f1..f4fabdb4f7 100644 --- a/SeQuant/domain/mbpt/biorthogonalization.hpp +++ b/SeQuant/domain/mbpt/biorthogonalization.hpp @@ -14,6 +14,10 @@ #include #include #endif +#if defined(SEQUANT_HAS_TAPP) +#include +#include +#endif #include #include @@ -377,6 +381,83 @@ auto biorthogonal_nns_project(btas::Tensor const& arr, #endif // defined(SEQUANT_HAS_BTAS) +#if defined(SEQUANT_HAS_TAPP) + +/// \brief This function is used to implement +/// ResultPtr::biorthogonal_nns_project for TAPPTensor +/// +/// \param arr The tensor to be "cleaned up" +/// \param bra_rank The rank of the bra indices +/// +/// \return The cleaned TAPPTensor. +template +auto biorthogonal_nns_project_tapp(TAPPTensor const& arr, + size_t bra_rank) { + using ranges::views::iota; + size_t const rank = arr.rank(); + SEQUANT_ASSERT(bra_rank <= rank); + size_t const ket_rank = rank - bra_rank; + + // Residuals of rank 4 or less have no redundancy and don't require NNS + // projection + if (rank <= 4) return arr; + + using numeric_type = T; + + const auto& nns_p_coeffs = + detail::nns_projection_weights(ket_rank); + + using perm_type = container::svector; + + TAPPTensor result; + + perm_type perm = iota(size_t{0}, rank) | ranges::to; + perm_type bra_perm = iota(size_t{0}, bra_rank) | ranges::to; + perm_type ket_perm = iota(bra_rank, rank) | ranges::to; + + if (ket_rank > 2 && !nns_p_coeffs.empty()) { + bool result_initialized = false; + + size_t num_perms = nns_p_coeffs.size(); + for (size_t perm_rank = 0; perm_rank < num_perms; ++perm_rank) { + perm_type permuted_ket = + detail::compute_permuted_indices(ket_perm, perm_rank, ket_rank); + + numeric_type coeff = nns_p_coeffs[perm_rank]; + + perm_type annot = bra_perm; + annot.insert(annot.end(), permuted_ket.begin(), permuted_ket.end()); + + container::svector annot_i64(annot.begin(), annot.end()); + container::svector perm_i64(perm.begin(), perm.end()); + + TAPPTensor temp; + tapp_ops::permute(arr, annot_i64, temp, perm_i64); + tapp_ops::scal(coeff, temp); + + if (result_initialized) { + result += temp; + } else { + result = temp; + result_initialized = true; + } + } + + } else { + result = arr; + } + + return result; +} + +template +auto biorthogonal_nns_project(TAPPTensor const& arr, + size_t bra_rank) { + return biorthogonal_nns_project_tapp(arr, bra_rank); +} + +#endif // defined(SEQUANT_HAS_TAPP) + } // namespace sequant::mbpt #endif diff --git a/cmake/modules/FindOrFetchBTAS.cmake b/cmake/modules/FindOrFetchBTAS.cmake index 22253cf3b4..e6dac87e86 100644 --- a/cmake/modules/FindOrFetchBTAS.cmake +++ b/cmake/modules/FindOrFetchBTAS.cmake @@ -21,7 +21,7 @@ if (NOT TARGET BTAS::BTAS) FetchContent_Declare( BTAS GIT_REPOSITORY https://github.com/BTAS/btas.git - GIT_TAG ${TA_TRACKED_BTAS_TAG} + GIT_TAG ${SEQUANT_TRACKED_BTAS_TAG} EXCLUDE_FROM_ALL SYSTEM ) diff --git a/cmake/modules/FindOrFetchTAPP.cmake b/cmake/modules/FindOrFetchTAPP.cmake new file mode 100644 index 0000000000..0c1c6faa88 --- /dev/null +++ b/cmake/modules/FindOrFetchTAPP.cmake @@ -0,0 +1,29 @@ +# try find_package +if (NOT TARGET tapp::reference) + include(FindPackageRegimport) + find_package_regimport(tapp QUIET CONFIG COMPONENTS reference) + if (TARGET tapp::reference) + message(STATUS "Found tapp CONFIG at ${tapp_CONFIG}") + endif() +endif() + +# if not found, build via FetchContent +if (NOT TARGET tapp::reference) + include(FetchContent) + FetchContent_Declare( + tapp + GIT_REPOSITORY https://github.com/TAPPorg/reference-implementation.git + GIT_TAG ${SEQUANT_TRACKED_TAPP_TAG} + SYSTEM + ) + FetchContent_MakeAvailable(tapp) + + # set tapp_CONFIG to the install location so that sequant-config.cmake knows where to find it + set(tapp_CONFIG ${CMAKE_INSTALL_PREFIX}/${TAPP_INSTALL_CMAKEDIR}/tapp-config.cmake) + +endif() + +# postcond check +if (NOT TARGET tapp::reference) + message(FATAL_ERROR "FindOrFetchTAPP could not make tapp::reference target available") +endif() diff --git a/cmake/sequant-config.cmake.in b/cmake/sequant-config.cmake.in index 7597f69fbf..c6da4ec1c3 100644 --- a/cmake/sequant-config.cmake.in +++ b/cmake/sequant-config.cmake.in @@ -44,6 +44,16 @@ if(SEQUANT_HAS_BTAS AND NOT TARGET BTAS::BTAS) find_dependency(BTAS CONFIG QUIET REQUIRED PATHS ${BTAS_DIR} NO_DEFAULT_PATH) endif() +set(SEQUANT_HAS_TAPP @SEQUANT_HAS_TAPP@) +if(SEQUANT_HAS_TAPP AND NOT TARGET tapp::reference) + set(tapp_CONFIG @tapp_CONFIG@) + if (NOT tapp_CONFIG OR NOT EXISTS ${tapp_CONFIG}) + message(FATAL_ERROR "Expected tapp config at ${tapp_CONFIG}; directory moved since SeQuant configuration?") + endif() + get_filename_component(tapp_DIR ${tapp_CONFIG} DIRECTORY) + find_dependency(tapp CONFIG QUIET REQUIRED COMPONENTS reference PATHS ${tapp_DIR} NO_DEFAULT_PATH) +endif() + set(SEQUANT_HAS_EIGEN @SEQUANT_HAS_EIGEN@) if (NOT TARGET Eigen3::Eigen AND SEQUANT_HAS_EIGEN) if (TARGET TiledArray_Eigen) @@ -86,6 +96,9 @@ if(NOT TARGET SeQuant::SeQuant) if (SEQUANT_HAS_BTAS) list(APPEND grandTargetList eval::btas) endif() + if (SEQUANT_HAS_TAPP) + list(APPEND grandTargetList eval::tapp) + endif() foreach(_target ${grandTargetList}) if(NOT TARGET SeQuant::${_target}) message(FATAL_ERROR "expected SeQuant::${_target} among imported SeQuant targets") diff --git a/external/versions.cmake b/external/versions.cmake index 7379bbb1aa..79c59501e5 100644 --- a/external/versions.cmake +++ b/external/versions.cmake @@ -5,7 +5,7 @@ set(SEQUANT_TRACKED_VGCMAKEKIT_TAG a65dfddae0c90bae8e1826b21e3b10cb9a98a9e8) set(SEQUANT_TRACKED_RANGEV3_TAG 0.12.0) -set(SEQUANT_TRACKED_TILEDARRAY_TAG 5477174cd8e857f76d29cf6c0dcc80277393a9eb) +set(SEQUANT_TRACKED_TILEDARRAY_TAG cdbd9bb4a1907843593ca90a5252938a188569b6) set(SEQUANT_TRACKED_LIBPERM_TAG cada3e185549896203cf4d0c7f26ea22c7de428f) @@ -41,3 +41,6 @@ set(SEQUANT_OLDEST_DOXYGEN_VERSION 1.9.6) set(SEQUANT_TRACKED_EIGEN_TAG 5.0.1) set(SEQUANT_OLDEST_EIGEN_VERSION 3.0...5) + +set(SEQUANT_TRACKED_BTAS_TAG 9c8c8f68fee2b82e64755270a8348e4612cf9941) +set(SEQUANT_TRACKED_TAPP_TAG 5d6fb56f7d4cbb4daefe23f408bbdde8cf9fc015) diff --git a/tests/integration/eval/CMakeLists.txt b/tests/integration/eval/CMakeLists.txt index 8b389424f2..04d84a4adf 100644 --- a/tests/integration/eval/CMakeLists.txt +++ b/tests/integration/eval/CMakeLists.txt @@ -49,3 +49,21 @@ if (SEQUANT_HAS_BTAS) build_test_as_needed(eval_btas "${test_name}" test_name) endif (SEQUANT_HAS_BTAS) + +if (SEQUANT_HAS_TAPP) + add_executable(eval_tapp ${BUILD_BY_DEFAULT} + "tapp/data_world_tapp.hpp" + "tapp/scf_tapp.hpp" + "tapp/main.cpp" + ) + target_link_libraries(eval_tapp PRIVATE eval_shared tapp::reference) + + set(test_name "sequant/integration/eval_tapp") + add_test( + NAME "${test_name}" + COMMAND eval_tapp + WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}" + ) + + build_test_as_needed(eval_tapp "${test_name}" test_name) +endif (SEQUANT_HAS_TAPP) diff --git a/tests/integration/eval/tapp/data_world_tapp.hpp b/tests/integration/eval/tapp/data_world_tapp.hpp new file mode 100644 index 0000000000..a0a8dbbae8 --- /dev/null +++ b/tests/integration/eval/tapp/data_world_tapp.hpp @@ -0,0 +1,215 @@ +// +// TAPP backend data world for integration tests. +// Mirrors btas/data_world_btas.hpp +// + +#ifndef SEQUANT_EVAL_DATA_WORLD_TAPP_HPP +#define SEQUANT_EVAL_DATA_WORLD_TAPP_HPP + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sequant::eval::tapp { + +template +class DataWorldTAPP { + private: + size_t const nocc; + + size_t const nvirt; + + std::string const fock_file; + + std::string const eri_file; + + Tensor_t G_pqrs; + + Tensor_t F_pq; + + container::vector Ts; + + mutable container::map cache_; + + public: + DataWorldTAPP(DataInfo const& info, size_t excit) + : nocc{info.nocc()}, + nvirt{info.nvirt()}, + fock_file{info.fock_file()}, + eri_file{info.eri_file()} { + using namespace ranges::views; + + //-------------------------------------------------------- + // read ERI and Fock tensors + //-------------------------------------------------------- + auto const nobs = nocc + nvirt; + F_pq = Tensor_t{repeat_n(static_cast(nobs), DataInfo::fock_rank) | + ranges::to>}; + F_pq.fill(0); + G_pqrs = Tensor_t{repeat_n(static_cast(nobs), DataInfo::eri_rank) | + ranges::to>}; + G_pqrs.fill(0); + read_tensor(fock_file, F_pq); + read_tensor(eri_file, G_pqrs); + //-------- + + //-------------------------------------------------------- + // Initialize amplitude tensors + //-------------------------------------------------------- + Ts.reserve(excit); + for (size_t i = 0; i < excit; ++i) { + size_t const bk_rank = i + 1; + + auto extents = concat(repeat_n(static_cast(nvirt), bk_rank), + repeat_n(static_cast(nocc), bk_rank)) | + ranges::to>; + Ts.emplace_back( + eval_result>(Tensor_t{extents})); + auto& t = Ts[i]->template get(); + t.fill(0); + } + //-------- + } // ctor + + /// + /// Update an amplitude tensor using a corresponding residual tensor. + /// + /// \param T Amplitude tensor with the layout {virt, virt, ..., occ, occ,} + /// \param R Residual tensor with the layout {virt, virt, ..., occ, occ,} + /// \note The ranks and the layouts of T and R must match exactly. + /// \note The rank is even. + /// + void update_single_T(Tensor_t& T, Tensor_t const& R) { + using namespace ranges::views; + SEQUANT_ASSERT(T.rank() == R.rank() && "Incompatible ranks of R and T"); + + SEQUANT_ASSERT(T.rank() % 2 == 0 && "Odd rank not supported"); + + size_t const n = T.rank() / 2; + + auto const& F = F_pq; + auto updater = [&T, &R, &F, n, this](auto const& idx_vec) { + double diag = 0.; + for (auto x : idx_vec | take(n)) diag -= F(x + nocc, x + nocc); + for (auto x : idx_vec | drop(n) | take(n)) diag += F(x, x); + + T(idx_vec) += R(idx_vec) / diag; + }; + + auto virt_occ = + concat(repeat_n(iota(size_t{0}, nvirt) | ranges::to_vector, n), + repeat_n(iota(size_t{0}, nocc) | ranges::to_vector, n)) | + ranges::to_vector; + + cartesian_foreach(virt_occ, updater); + } + + Tensor_t operator()(sequant::Tensor const& tensor) const { + using namespace ranges::views; + + SEQUANT_ASSERT(tensor.label() != L"t"); + + auto r1_limits = range1_limits(tensor, nocc, nvirt); + SEQUANT_ASSERT(r1_limits.size() == DataInfo::fock_rank || + r1_limits.size() == DataInfo::eri_rank); + + auto const iter_limits = + r1_limits | ranges::views::transform([this](auto x) { + return x == nocc ? std::pair{size_t{0}, nocc} + : std::pair{nocc, nocc + nvirt}; + }); + + auto const tlabel = tensor.label(); + SEQUANT_ASSERT(tlabel == L"g" || tlabel == L"f"); + + auto const& big_tensor = tlabel == L"g" ? G_pqrs : F_pq; + + auto r1_extents = r1_limits | ranges::views::transform([](size_t v) { + return static_cast(v); + }) | + ranges::to>; + + auto slice = Tensor_t{r1_extents}; + if (iter_limits.size() == DataInfo::fock_rank) { + auto loop1 = iter_limits[0]; + auto loop2 = iter_limits[1]; + for (auto i = loop1.first; i < loop1.second; ++i) + for (auto j = loop2.first; j < loop2.second; ++j) + slice(i - loop1.first, j - loop2.first) = big_tensor(i, j); + + } else { // DataInfo::eri_rank + auto loop1 = iter_limits[0]; + auto loop2 = iter_limits[1]; + auto loop3 = iter_limits[2]; + auto loop4 = iter_limits[3]; + for (auto i = loop1.first; i < loop1.second; ++i) + for (auto j = loop2.first; j < loop2.second; ++j) + for (auto k = loop3.first; k < loop3.second; ++k) + for (auto l = loop4.first; l < loop4.second; ++l) + slice(i - loop1.first, // + j - loop2.first, // + k - loop3.first, // + l - loop4.first) = // + big_tensor(i, j, k, l); + } + return slice; + } + + ResultPtr operator()(meta::can_evaluate auto const& n) const { + using numeric_type = typename Tensor_t::numeric_type; + if (n->result_type() == ResultType::Scalar) { + SEQUANT_ASSERT(n->expr()->template is()); + auto d = n->as_constant().template value(); + return eval_result>(d); + } + + SEQUANT_ASSERT(n->result_type() == ResultType::Tensor && + n->expr()->template is()); + + if (auto t = n->as_tensor(); t.label() == L"t") { + auto rank = t.rank(); + SEQUANT_ASSERT(rank <= Ts.size()); + return Ts[rank - 1]; + } + auto h = hash::value(*n); + if (auto exists = cache_.find(h); exists != cache_.end()) + return exists->second; + else { + auto tnsr = + eval_result>((*this)(n->as_tensor())); + auto stored = cache_.emplace(h, std::move(tnsr)); + SEQUANT_ASSERT(stored.second && "failed to store tensor"); + return stored.first->second; + } + } + + void update_amplitudes(std::vector const& rs) { + using ranges::views::transform; + using ranges::views::zip; + + SEQUANT_ASSERT(rs.size() == Ts.size() && "Unequal number of Rs and Ts!"); + + for (auto&& [t, r] : zip(Ts | transform([](auto&& res) -> Tensor_t& { + return res->template get(); + }), + rs)) + update_single_T(t, r); + } + + Tensor_t const& amplitude(size_t excitation) const { + return Ts[excitation - 1]->template get(); + } + +}; // class + +} // namespace sequant::eval::tapp + +#endif // SEQUANT_EVAL_DATA_WORLD_TAPP_HPP diff --git a/tests/integration/eval/tapp/main.cpp b/tests/integration/eval/tapp/main.cpp new file mode 100644 index 0000000000..6dbde1d045 --- /dev/null +++ b/tests/integration/eval/tapp/main.cpp @@ -0,0 +1,103 @@ +// +// TAPP backend integration test entry point. +// Mirrors btas/main.cpp +// + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#define runtime_assert(tf) \ + if (!(tf)) { \ + std::ostringstream oss; \ + oss << "failed assert at line " << __LINE__ << " in eval_tapp example"; \ + throw std::runtime_error(oss.str().c_str()); \ + } + +/// +/// excitation (2 = ccsd (default) through 6 supported) +/// fock(/eri) tensor data file name (default "fock.dat"/"eri.dat") +/// eri(/fock) tensor data file name (default "fock.dat"/"eri.dat") +/// \note Format of data files: +/// ---------------------------- +/// size_t size_t size_t # rank, nocc, nvirt +/// double # data ------ +/// ... # data | +/// ... # .... | no. of double entries = +/// (nocc+nvirt)^rank +/// ... # data | +/// double # data ------ +/// ---------------------------- +/// \note The rank of fock tensor is 2 +/// \note The rank of eri tensor is 4 +/// \note The nocc and nvirt in both files must match +/// +int main(int argc, char* argv[]) { + if (argc < 4) { + std::cout << "\nHelp:\n" + << " " + " []" + << "\n\n" + << "Config file format\n" + << "----\n" + << sequant::eval::ParseConfigFile{}.help() << "\n----\n"; + // return 1; + } + + sequant::set_locale(); + + using namespace sequant; + detail::OpIdRegistrar op_id_registrar; + sequant::set_default_context( + {.index_space_registry_shared_ptr = mbpt::make_min_sr_spaces(), + .vacuum = Vacuum::SingleProduct, + .canonicalization_options = + CanonicalizeOptions::default_options().copy_and_set( + CanonicalizationMethod::Complete)}); + TensorCanonicalizer::register_instance( + std::make_shared()); + mbpt::set_default_mbpt_context( + {.op_registry_ptr = mbpt::make_minimal_registry()}); + + // for optimization tests, set occupied and unoccupied index extents + { + auto reg = get_default_context().mutable_index_space_registry(); + auto occ = reg->retrieve_ptr(L"i"); + auto uocc = reg->retrieve_ptr(L"a"); + SEQUANT_ASSERT(occ); + SEQUANT_ASSERT(uocc); + occ->approximate_size(10); + uocc->approximate_size(100); + SEQUANT_ASSERT(uocc->approximate_size() == 100); + } + + std::string calc_config = argc > 1 ? argv[1] : "calc.inp"; + std::string fock_file = argc > 2 ? argv[2] : "fock_so.dat"; + std::string eri_file = argc > 3 ? argv[3] : "eri_so.dat"; + // not yet implemented: + std::string out_file = argc > 4 ? argv[4] : ""; + // + + using TAPPTensorD = sequant::TAPPTensor; + auto const calc_info = + eval::make_calc_info(calc_config, fock_file, eri_file, out_file); + auto scf_tapp = eval::tapp::SequantEvalScfTAPP{calc_info}; + scf_tapp.scf(std::wcout); + + double const expected{-0.07068045196165902}; + double const threshold = calc_info.scf_opts.conv; + double const ediff = std::fabs(expected - scf_tapp.energy()); + + runtime_assert((ediff <= threshold)); + return 0; +} diff --git a/tests/integration/eval/tapp/scf_tapp.hpp b/tests/integration/eval/tapp/scf_tapp.hpp new file mode 100644 index 0000000000..70baff9ee4 --- /dev/null +++ b/tests/integration/eval/tapp/scf_tapp.hpp @@ -0,0 +1,145 @@ +// +// TAPP backend SCF solver for integration tests. +// Mirrors btas/scf_btas.hpp +// + +#ifndef SEQUANT_EVAL_SCF_TAPP_HPP +#define SEQUANT_EVAL_SCF_TAPP_HPP + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +namespace sequant::eval::tapp { + +template +class SequantEvalScfTAPP final : public SequantEvalScf { + public: + using ExprT = EvalExprTAPP; + using EvalNodeTAPP = EvalNode; + + private: + container::vector> nodes_; + CacheManager cman_; + DataWorldTAPP data_world_; + + Tensor_t const& f_vo() const { + static Tensor_t tnsr = + data_world_(parse_expr(L"f{a1;i1}", Symmetry::Nonsymm)->as()); + return tnsr; + } + + Tensor_t const& g_vvoo() const { + static Tensor_t tnsr = data_world_( + parse_expr(L"g{a1,a2;i1,i2}", Symmetry::Nonsymm)->as()); + return tnsr; + } + + double energy_spin_orbital() { + static const std::wstring_view energy_expr = + L"f{i1;a1} * t{a1;i1} + g{i1,i2;a1,a2} * " + L"(1/4 * t{a1,a2;i1,i2} + 1/2 t{a1;i1} * t{a2;i2})"; + static auto const node = + binarize(parse_expr(energy_expr, Symmetry::Antisymm)); + + return evaluate(node, data_world_)->template get(); + } + + double energy_spin_free_orbital() const { + auto const& T1 = data_world_.amplitude(1); + auto const& T2 = data_world_.amplitude(2); + auto const& G_vvoo = g_vvoo(); + auto const& F_vo = f_vo(); + + Tensor_t temp; + tapp_ops::contract(1.0, T1, {'a', 'i'}, T1, {'b', 'j'}, 0.0, temp, + {'a', 'b', 'i', 'j'}); + + return 2.0 * (tapp_ops::dot(F_vo, T1) + tapp_ops::dot(G_vvoo, T2) + + tapp_ops::dot(G_vvoo, temp)) - + tapp_ops::dot(G_vvoo, temp); + } + + void reset_cache() override { cman_.reset(); } + + double norm() const override { + // todo use all Ts instead of only T2 + auto const& T2 = data_world_.amplitude(2); + return std::sqrt(tapp_ops::dot(T2, T2)); + } + + double solve() override { + using ranges::views::concat; + using ranges::views::repeat_n; + using ranges::views::transform; + using ranges::views::zip; + + using Evaluator = ResultPtr (*)( + std::vector const&, EvalExprTAPP::annot_t const&, + DataWorldTAPP const&, CacheManager&); + + constexpr auto funcs = + std::array{std::array{evaluate_antisymm, + evaluate_antisymm}, + std::array{evaluate_symm, + evaluate_symm}}; + + auto sorted_annot = [](Tensor const& tnsr) { + auto b = tnsr.bra() | ranges::to_vector; + auto k = tnsr.ket() | ranges::to_vector; + ranges::sort(b, Index::LabelCompare{}); + ranges::sort(k, Index::LabelCompare{}); + return EvalExprTAPP::index_hash(concat(b, k)) | + ranges::to; + }; + + auto rs = repeat_n(Tensor_t{}, info_.eqn_opts.excit) | ranges::to_vector; + auto st = info_.eqn_opts.spintrace; + auto log = Logger::instance().eval.level > 0; + + for (auto&& [r, n] : zip(rs, nodes_)) { + SEQUANT_ASSERT(!n.empty()); + // update_amplitudes assumes that the residuals are in specific layout + auto const target_indices = sorted_annot(ranges::front(n)->as_tensor()); + r = funcs[st][log](n, target_indices, data_world_, cman_) + ->template get(); + } + + const auto E = info_.eqn_opts.spintrace ? energy_spin_free_orbital() + : energy_spin_orbital(); + data_world_.update_amplitudes(rs); + return E; + } + + public: + SequantEvalScfTAPP(CalcInfo const& calc_info) + : SequantEvalScf{calc_info}, + cman_{{}}, + data_world_{calc_info.fock_eri, calc_info.eqn_opts.excit} { + assert(info_.eqn_opts.excit >= 2 && + "At least double excitation (CCSD) is required!"); + // todo time it + auto const exprs = info_.exprs(); + + nodes_ = info_.nodes(exprs); + + cman_ = info_.cache_manager_scf(nodes_); + } +}; + +} // namespace sequant::eval::tapp + +#endif // SEQUANT_EVAL_SCF_TAPP_HPP diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index bc3b9f67b0..90aac1ef2e 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -132,6 +132,21 @@ if (SEQUANT_HAS_BTAS) PROPERTIES SKIP_UNITY_BUILD_INCLUSION ON) endif() +########################## +# OBJECT library: TAPP backend tests (conditional) +########################## +if (SEQUANT_HAS_TAPP) + set(eval_tapp_test_sources "test_eval_tapp.cpp") + add_library(unit_tests-sequant-eval-tapp-obj OBJECT ${eval_tapp_test_sources}) + target_link_libraries(unit_tests-sequant-eval-tapp-obj + PUBLIC SeQuant::eval::tapp SeQuant::mbpt + PRIVATE Catch2::Catch2 dtl::dtl) + target_compile_definitions(unit_tests-sequant-eval-tapp-obj PRIVATE + SEQUANT_UNIT_TESTS_SOURCE_DIR="${CMAKE_CURRENT_SOURCE_DIR}") + set_source_files_properties(${eval_tapp_test_sources} + PROPERTIES SKIP_UNITY_BUILD_INCLUSION ON) +endif() + ########################## # Single executable combining all OBJECT libraries # NB linking OBJECT libraries (since CMake 3.21) both includes their object @@ -155,6 +170,9 @@ endif() if (SEQUANT_HAS_BTAS) target_link_libraries(unit_tests-sequant PRIVATE unit_tests-sequant-eval-btas-obj) endif() +if (SEQUANT_HAS_TAPP) + target_link_libraries(unit_tests-sequant PRIVATE unit_tests-sequant-eval-tapp-obj) +endif() target_compile_definitions(unit_tests-sequant PRIVATE SEQUANT_UNIT_TESTS_SOURCE_DIR="${CMAKE_CURRENT_SOURCE_DIR}") diff --git a/tests/unit/test_eval_btas.cpp b/tests/unit/test_eval_btas.cpp index a406b50b5e..d65305d285 100644 --- a/tests/unit/test_eval_btas.cpp +++ b/tests/unit/test_eval_btas.cpp @@ -119,7 +119,7 @@ class rand_tensor_yield { /// otherwise throws assertion error. To avoid that use the other /// overload of operator() that takes sequant::Tensor const& sequant::ResultPtr operator()(std::wstring_view label) const { - auto&& found = label_to_tnsr_.find(label.data()); + auto&& found = label_to_tnsr_.find(std::wstring{label}); if (found == label_to_tnsr_.end()) { SEQUANT_ASSERT(false && "attempted access of non-existent tensor!"); } @@ -199,7 +199,6 @@ TEST_CASE("eval_with_btas", "[eval_btas]") { auto eval_symm = [&yield_](sequant::ExprPtr const& expr, container::svector const& target_labels) { - auto cache = sequant::CacheManager::empty(); return evaluate_symm(eval_node(expr), target_labels, yield_) ->get(); }; @@ -440,6 +439,11 @@ TEST_CASE("eval_with_btas", "[eval_btas]") { auto const eval2 = evaluate(nodes1, tidx1, yield_)->get(); - REQUIRE(norm(eval1) == Catch::Approx(norm(eval1))); + REQUIRE(norm(eval1) == Catch::Approx(norm(eval2))); + + BTensorD zero2{eval2.range()}; + zero2 = eval1 - eval2; + REQUIRE(norm(zero2) == Catch::Approx(0).margin( + 100 * std::numeric_limits::epsilon())); } } diff --git a/tests/unit/test_eval_tapp.cpp b/tests/unit/test_eval_tapp.cpp new file mode 100644 index 0000000000..5d655cddf7 --- /dev/null +++ b/tests/unit/test_eval_tapp.cpp @@ -0,0 +1,451 @@ +#include +#include + +#include "catch2_sequant.hpp" + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace { + +auto eval_node(sequant::ExprPtr const& expr) { + using namespace sequant; + return binarize(expr); +} + +static auto const idx_rgx = boost::wregex{L"([ia])([↑↓])?_?(\\d+)"}; +auto tensor_to_key(sequant::Tensor const& tnsr) { + auto formatter = [](boost::wsmatch mo) -> std::wstring { + return (mo[1].str() == L"i" ? L"o" : L"v") + mo[2].str(); + }; + + auto const tnsr_deparsed = sequant::deparse(tnsr.clone(), false); + return boost::regex_replace(tnsr_deparsed, idx_rgx, formatter); +} + +[[maybe_unused]] auto tensor_to_key(std::wstring_view spec) { + return tensor_to_key(sequant::parse_expr(spec, sequant::Symmetry::Nonsymm) + ->as()); +} + +using TAPPTensorD = sequant::TAPPTensor; + +template +class rand_tensor_yield { + private: + size_t const nocc_; + size_t const nvirt_; + mutable std::map label_to_tnsr_; + + public: + rand_tensor_yield(size_t noccupied, size_t nvirtual) + : nocc_{noccupied}, nvirt_{nvirtual} {} + + [[nodiscard]] Tensor_t make_rand_tensor(sequant::Tensor const& tnsr) const { + using ranges::views::transform; + using sequant::IndexSpace; + auto isr = sequant::get_default_context().index_space_registry(); + + SEQUANT_ASSERT( + ranges::all_of(tnsr.const_braket_indices(), + [&isr](auto const& idx) { + return idx.space() == isr->retrieve(L"i") || + idx.space() == isr->retrieve(L"a"); + }) && + "Unsupported IndexSpace type found while generating tensor."); + + auto extents = tnsr.const_braket_indices() | + transform([this, &isr](auto const& idx) { + return static_cast( + idx.space() == isr->retrieve(L"i") ? nocc_ : nvirt_); + }) | + ranges::to>; + + auto result = Tensor_t{extents}; + result.generate( + []() { return static_cast(std::rand()) / RAND_MAX; }); + return result; + } + + sequant::ResultPtr operator()(sequant::Tensor const& tnsr) const { + using namespace sequant; + std::wstring const label = tensor_to_key(tnsr); + if (auto&& found = label_to_tnsr_.find(label); + found != label_to_tnsr_.end()) { + return found->second; + } + auto t = make_rand_tensor(tnsr); + auto&& success = label_to_tnsr_.emplace( + label, eval_result>(std::move(t))); + SEQUANT_ASSERT(success.second && "couldn't store tensor!"); + return success.first->second; + } + + sequant::ResultPtr operator()( + sequant::meta::can_evaluate auto const& node) const { + using namespace sequant; + if (node->result_type() == sequant::ResultType::Tensor) { + SEQUANT_ASSERT(node->expr()->template is()); + return (*this)(node->expr()->template as()); + } + + using result_t = ResultScalar; + + SEQUANT_ASSERT(node->expr()->template is()); + auto d = node->as_constant().template value(); + return eval_result(d); + } + + /// + /// \param label eg. t{v,v;o,o}, f{o;v} + /// \return const ref to Tensor_t type tensor + /// \note The tensor should be already present in the yielder cache + /// otherwise throws assertion error. To avoid that use the other + /// overload of operator() that takes sequant::Tensor const& + sequant::ResultPtr operator()(std::wstring_view label) const { + auto&& found = label_to_tnsr_.find(std::wstring{label}); + if (found == label_to_tnsr_.end()) { + SEQUANT_ASSERT(false && "attempted access of non-existent tensor!"); + } + return found->second; + } +}; + +using namespace sequant; + +template < + typename Iterable, + std::enable_if_t< + std::is_convertible_v, Index> && + !sequant::meta::is_statically_castable_v< + Iterable const&, std::wstring> // prefer the ctor taking the + // std::wstring + , + bool> = true> +container::svector tidxs(Iterable const& indices) noexcept { + return sequant::EvalExprTAPP::index_hash(indices) | + ranges::to>; +} + +container::svector tidxs(Tensor const& tnsr) noexcept { + return sequant::EvalExprTAPP::index_hash(tnsr.const_braket_indices()) | + ranges::to>; +} + +container::svector tidxs( + ExprPtr expr, std::initializer_list tnsr_coords) noexcept { + auto tnsr_p = expr; + for (auto i : tnsr_coords) tnsr_p = tnsr_p->at(i); + SEQUANT_ASSERT(tnsr_p->is()); + return tidxs(tnsr_p->as()); +} + +auto terse_index = [](std::wstring const& spec) { + auto formatter = [](boost::wsmatch mo) -> std::wstring { + return mo[1].str() + mo[2].str() + L"_" + mo[3].str(); + }; + return boost::regex_replace(spec, idx_rgx, formatter); +}; + +container::svector tidxs(std::wstring const& csv) noexcept { + using ranges::views::all; + using ranges::views::split; + using ranges::views::transform; + auto const detersed = terse_index(csv); + return tidxs(detersed | split(L',') | + transform([](auto&& v) { return ranges::to(v); })); +} + +} // namespace + +TEST_CASE("eval_with_tapp", "[eval_tapp]") { + using ranges::views::transform; + using namespace sequant; + + auto norm = [](TAPPTensorD const& tnsr) { + return std::sqrt(tapp_ops::dotc(tnsr, tnsr)); + }; + + std::srand(2023); + const size_t nocc = 2, nvirt = 20; + auto yield_ = rand_tensor_yield{nocc, nvirt}; + auto yield = [&yield_](std::wstring_view lbl) -> TAPPTensorD const& { + return yield_(lbl)->get(); + }; + + auto eval = [&yield_](sequant::ExprPtr const& expr, + container::svector const& target_labels) { + return evaluate(eval_node(expr), target_labels, yield_)->get(); + }; + + auto eval_symm = [&yield_](sequant::ExprPtr const& expr, + container::svector const& target_labels) { + return evaluate_symm(eval_node(expr), target_labels, yield_) + ->get(); + }; + + auto eval_antisymm = [&yield_]( + sequant::ExprPtr const& expr, + container::svector const& target_labels) { + return evaluate_antisymm(eval_node(expr), target_labels, yield_) + ->get(); + }; + + auto eval_biorthogonal_nns_project = + [&yield_](sequant::ExprPtr const& expr, + container::svector const& target_labels) { + auto result = evaluate(eval_node(expr), target_labels, yield_); + return mbpt::biorthogonal_nns_project( + result->get(), + eval_node(expr)->as_tensor().bra_rank()); + }; + + auto parse_antisymm = [](auto const& xpr) { + return parse_expr(xpr, sequant::Symmetry::Antisymm); + }; + + SECTION("Summation") { + auto expr1 = parse_antisymm(L"t_{a1}^{i1} + f_{i1}^{a1}"); + auto const tidx1 = tidxs(expr1, {0}); + auto eval1 = eval(expr1, tidx1); + + auto man1 = yield(L"t{v;o}"); + TAPPTensorD f_permuted; + tapp_ops::permute(yield(L"f{o;v}"), {0, 1}, f_permuted, {1, 0}); + man1 += f_permuted; + + REQUIRE(norm(eval1) == Catch::Approx(norm(man1))); + + auto expr2 = parse_antisymm(L"2 * t_{a1}^{i1} + 3/2 * f_{i1}^{a1}"); + auto const tidx2 = tidxs(expr2, {0, 0}); + auto eval2 = eval(expr2, tidx2); + + auto man2 = yield(L"t{v;o}"); + tapp_ops::scal(2.0, man2); + TAPPTensorD temp; + tapp_ops::permute(yield(L"f{o;v}"), {0, 1}, temp, {1, 0}); + tapp_ops::scal(1.5, temp); + man2 += temp; + + REQUIRE(norm(eval2) == Catch::Approx(norm(man2))); + } + + SECTION("Product") { + auto expr1 = + parse_antisymm(L"1/2 * g_{i2,i4}^{a2,a4} * t_{a1,a2}^{ i1, i2}"); + auto const tidx1 = tidxs(L"i1,i4,a1,a4"); + auto eval1 = eval(expr1, tidx1); + + // mnemonics + // === + // i looks like 1 + // a looks like 7 + + TAPPTensorD man1; + auto const& g = yield(L"g{o,o;v,v}"); + auto const& t2 = yield(L"t{v,v;o,o}"); + tapp_ops::contract(0.5, g, {12, 14, 72, 74}, t2, {71, 72, 11, 12}, 0.0, + man1, {11, 14, 71, 74}); + REQUIRE(norm(eval1) == Catch::Approx(norm(man1))); + + auto expr2 = parse_antisymm( + L"-1/4 * g_{i3,i4}^{a3,a4} * t_{a1,a3}^{i3,i4} * t_{a2,a4}^{ i1, " + L"i2}"); + auto tidx2 = tidxs(L"i1,i2,a1,a2"); + auto eval2 = eval(expr2, tidx2); + + TAPPTensorD man2, temp2; + tapp_ops::contract(1.0, g, {13, 14, 73, 74}, t2, {71, 73, 13, 14}, 0.0, + temp2, {71, 74}); + tapp_ops::contract(-0.25, temp2, {71, 74}, t2, {72, 74, 11, 12}, 0.0, man2, + {11, 12, 71, 72}); + REQUIRE(norm(eval2) == Catch::Approx(norm(man2))); + } + + SECTION("Summation and Product") { + auto expr1 = parse_antisymm( + L"-1/4 * g_{i3,i4}^{a3,a4} * t_{a2,a4}^{i1,i2} * t_{a1,a3}^{i3,i4}" + " + " + " 1/16 * g_{i3,i4}^{a3,a4} * t_{a1,a2}^{i3,i4} * " + "t_{a3,a4}^{i1,i2}"); + auto tidx1 = tidxs(L"i1,i2,a1,a2"); + auto eval1 = eval(expr1, tidx1); + + auto const& g = yield(L"g{o,o;v,v}"); + auto const& t2 = yield(L"t{v,v;o,o}"); + TAPPTensorD temp1, man1; + tapp_ops::contract(1.0, g, {13, 14, 73, 74}, t2, {71, 73, 13, 14}, 0.0, + temp1, {71, 74}); + tapp_ops::contract((-1 / 4.0), temp1, {71, 74}, t2, {72, 74, 11, 12}, 0.0, + man1, {11, 12, 71, 72}); + temp1.clear(); + tapp_ops::contract(1.0, g, {13, 14, 73, 74}, t2, {73, 74, 11, 12}, 0.0, + temp1, {11, 12, 13, 14}); + TAPPTensorD temp2; + tapp_ops::contract((1 / 16.0), temp1, {11, 12, 13, 14}, t2, + {71, 72, 13, 14}, 0.0, temp2, {11, 12, 71, 72}); + man1 += temp2; + temp1.clear(); + temp2.clear(); + REQUIRE(norm(eval1) == Catch::Approx(norm(man1))); + } + + SECTION("Antisymmetrization") { + auto expr1 = parse_antisymm(L"g_{i1, i2}^{a1, a2}"); + auto tidx1 = tidxs(L"i_1,i_2,a_1,a_2"); + auto eval1 = eval_antisymm(expr1, tidx1); + + auto const& g = yield(L"g{o,o;v,v}"); + TAPPTensorD man1{g.extents()}; + man1.fill(0); + TAPPTensorD temp_as{g.extents()}; + temp_as.fill(0); + + TAPPTensorD perm_res; + tapp_ops::permute(g, {0, 1, 2, 3}, perm_res, {0, 1, 2, 3}); + man1 += perm_res; + + tapp_ops::permute(g, {0, 1, 2, 3}, perm_res, {1, 0, 2, 3}); + tapp_ops::scal(-1.0, perm_res); + man1 += perm_res; + + tapp_ops::permute(g, {0, 1, 2, 3}, perm_res, {0, 1, 3, 2}); + tapp_ops::scal(-1.0, perm_res); + man1 += perm_res; + + tapp_ops::permute(g, {0, 1, 2, 3}, perm_res, {1, 0, 3, 2}); + man1 += perm_res; + + tapp_ops::scal(0.25, man1); + + REQUIRE(norm(eval1) == Catch::Approx(norm(man1))); + + auto expr2 = parse_antisymm(L"R_{a1,a2}^{i1}"); + auto tidx2 = tidxs(L"a_1,a_2,i_1"); + auto eval2 = eval_antisymm(expr2, tidx2); + + auto const& r = yield(L"R{v,v;o}"); + TAPPTensorD man2{r.extents()}; + man2.fill(0.0); + TAPPTensorD temp2_as{r.extents()}; + temp2_as.fill(0.0); + + tapp_ops::permute(r, {0, 1, 2}, perm_res, {0, 1, 2}); + man2 += perm_res; + + tapp_ops::permute(r, {0, 1, 2}, perm_res, {1, 0, 2}); + tapp_ops::scal(-1.0, perm_res); + man2 += perm_res; + tapp_ops::scal(0.5, man2); + + REQUIRE(norm(eval2) == Catch::Approx(norm(man2))); + } + + SECTION("Symmetrization") { + auto expr1 = parse_antisymm(L"g_{i1, i2}^{a1, a2}"); + auto tidx1 = tidxs(L"i_1,i_2,a_1,a_2"); + auto eval1 = eval_symm(expr1, tidx1); + + auto const& g = yield(L"g{o,o;v,v}"); + + TAPPTensorD man1{g.extents()}; + man1.fill(0); + + TAPPTensorD perm_res; + tapp_ops::permute(g, {0, 1, 2, 3}, perm_res, {0, 1, 2, 3}); + man1 += perm_res; + tapp_ops::permute(g, {0, 1, 2, 3}, perm_res, {1, 0, 3, 2}); + man1 += perm_res; + tapp_ops::scal(0.5, man1); + + REQUIRE(norm(eval1) == Catch::Approx(norm(man1))); + } + + SECTION("Biorthogonal Cleanup") { + // low-rank residuals: skip cleanup + auto expr1 = parse_antisymm(L"R_{a1, a2}^{i1, i2}"); + auto tidx1 = tidxs(L"a_1,a_2,i_1,i_2"); + auto eval1 = eval_biorthogonal_nns_project(expr1, tidx1); + auto const& r1 = yield(L"R{v,v;o,o}"); + + TAPPTensorD man1{r1.extents()}; + man1.fill(0); + // copy r1 data + for (size_t i = 0; i < r1.volume(); ++i) + const_cast(man1.data())[i] = r1.data()[i]; + + REQUIRE(norm(eval1) == Catch::Approx(norm(man1))); + + TAPPTensorD zero1 = man1 - eval1; + REQUIRE(norm(zero1) == Catch::Approx(0).margin( + 100 * std::numeric_limits::epsilon())); + + // high-rank residuals: cleanup applies: + // result = identity - (1/ket_rank!) * sum_of_ket_permutations + auto expr2 = parse_antisymm(L"R_{a1, a2, a3}^{i1, i2, i3}"); + auto tidx2 = tidxs(L"a_1,a_2,a_3,i_1,i_2,i_3"); + auto eval2 = eval_biorthogonal_nns_project(expr2, tidx2); + auto const& r2 = yield(L"R{v,v,v;o,o,o}"); + + TAPPTensorD man2{r2.extents()}; + for (size_t i = 0; i < r2.volume(); ++i) + const_cast(man2.data())[i] = r2.data()[i]; + + TAPPTensorD perm_sum{r2.extents()}; + perm_sum.fill(0); + + TAPPTensorD perm_res; + tapp_ops::permute(r2, {0, 1, 2, 3, 4, 5}, perm_res, {0, 1, 2, 3, 5, 4}); + perm_sum += perm_res; + tapp_ops::permute(r2, {0, 1, 2, 3, 4, 5}, perm_res, {0, 1, 2, 4, 3, 5}); + perm_sum += perm_res; + tapp_ops::permute(r2, {0, 1, 2, 3, 4, 5}, perm_res, {0, 1, 2, 4, 5, 3}); + perm_sum += perm_res; + tapp_ops::permute(r2, {0, 1, 2, 3, 4, 5}, perm_res, {0, 1, 2, 5, 3, 4}); + perm_sum += perm_res; + tapp_ops::permute(r2, {0, 1, 2, 3, 4, 5}, perm_res, {0, 1, 2, 5, 4, 3}); + perm_sum += perm_res; + + tapp_ops::scal(1.0 / 5.0, perm_sum); + auto man2_result = man2 - perm_sum; + REQUIRE(norm(eval2) == Catch::Approx(norm(man2_result))); + + TAPPTensorD zero2 = man2_result - eval2; + REQUIRE(norm(zero2) == Catch::Approx(0).margin( + 100 * std::numeric_limits::epsilon())); + } + + SECTION("Others") { + auto expr1 = parse_antisymm( + L"-1/4 * g_{i3,i4}^{a3,a4} * t_{a2,a4}^{i1,i2} * t_{a1,a3}^{i3,i4}" + " + " + " 1/16 * g_{i3,i4}^{a3,a4} * t_{a1,a2}^{i3,i4} * " + "t_{a3,a4}^{i1,i2}"); + + auto tidx1 = tidxs(L"i1,i2,a1,a2"); + + auto const eval1 = + evaluate(eval_node(expr1), tidx1, yield_)->get(); + + auto nodes1 = *expr1 | ranges::views::transform([](auto&& x) { + return eval_node(x); + }) | ranges::to_vector; + + auto const eval2 = evaluate(nodes1, tidx1, yield_)->get(); + + REQUIRE(norm(eval1) == Catch::Approx(norm(eval2))); + + TAPPTensorD zero2 = eval1 - eval2; + REQUIRE(norm(zero2) == Catch::Approx(0).margin( + 100 * std::numeric_limits::epsilon())); + } +}