From b15d062d02a32ec086f91fffc4db3c5a7872e54b Mon Sep 17 00:00:00 2001 From: "Jonathan M. Waldrop" Date: Thu, 6 Feb 2025 12:27:38 -0600 Subject: [PATCH] Sigma Uncertainty Tracking Option (#196) * add sigma dependency * Compiles, but tests need updating * correct sigma conditionals * optional eigen not working, so making it required * backup: mostly works, issues with eigen tensor contraction * comment out contraction cases for now * comment out failing contraction cases for now * cmake update and uncomment tests * Committing clang-format changes * add sigma tests to CI --------- Co-authored-by: github-actions[bot] --- .github/enable_sigma.cmake | 1 + .github/workflows/pull_request.yaml | 10 + CMakeLists.txt | 46 +- cmake/config.hpp.in | 13 - include/tensorwrapper/allocator/eigen.hpp | 39 +- include/tensorwrapper/backends/eigen.hpp | 12 - include/tensorwrapper/buffer/eigen.hpp | 45 +- src/tensorwrapper/allocator/eigen.cpp | 35 +- src/tensorwrapper/buffer/eigen.cpp | 35 +- .../tensorwrapper/allocator/eigen.cpp | 8 +- .../tensorwrapper/buffer/buffer_base.cpp | 104 +- .../unit_tests/tensorwrapper/buffer/eigen.cpp | 1073 +++++++++-------- .../buffer/eigen_contraction.cpp | 8 +- 13 files changed, 729 insertions(+), 700 deletions(-) create mode 100644 .github/enable_sigma.cmake delete mode 100644 cmake/config.hpp.in diff --git a/.github/enable_sigma.cmake b/.github/enable_sigma.cmake new file mode 100644 index 00000000..2fbe8f6a --- /dev/null +++ b/.github/enable_sigma.cmake @@ -0,0 +1 @@ +set(ENABLE_SIGMA ON) \ No newline at end of file diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index 31fd3f9d..d7585def 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -29,3 +29,13 @@ jobs: compilers: '["gcc-11", "clang-14"]' doc_target: 'tensorwrapper_cxx_api' secrets: inherit + Test-Enable-Sigma: + uses: NWChemEx/.github/.github/workflows/common_pull_request.yaml@master + with: + config_file: '' + source_dir: '' + format_python: false + compilers: '["gcc-11", "clang-14"]' + repo_toolchain: ".github/enable_sigma.cmake" + doc_target: '' + secrets: inherit diff --git a/CMakeLists.txt b/CMakeLists.txt index 4a0a5cb7..28555c61 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,43 +36,40 @@ nwx_cxx_api_docs("${project_inc_dir}" "${project_src_dir}") cmaize_option_list( BUILD_TESTING OFF "Should we build the tests?" BUILD_PYBIND11_PYBINDINGS OFF "Should we build Python3 bindings?" - ENABLE_EIGEN_SUPPORT ON "Should we enable Eigen support?" + ENABLE_SIGMA OFF "Should we enable Sigma for uncertainty tracking?" ) -configure_file(cmake/config.hpp.in config.hpp @ONLY) +### Dependendencies ### +include(get_utilities) -cmaize_find_or_build_dependency( - utilities - URL github.com/NWChemEx/utilities - BUILD_TARGET utilities - FIND_TARGET nwx::utilities - CMAKE_ARGS BUILD_TESTING=OFF -) - -cmaize_find_or_build_dependency( - parallelzone - URL github.com/NWChemEx/ParallelZone - BUILD_TARGET parallelzone - FIND_TARGET nwx::parallelzone - CMAKE_ARGS BUILD_TESTING=OFF -) +include(get_parallelzone) find_package(Boost REQUIRED) -cmaize_find_or_build_optional_dependency( +cmaize_find_or_build_dependency( eigen - ENABLE_EIGEN_SUPPORT URL https://www.gitlab.com/libeigen/eigen VERSION 2e76277bd049f7bec36b0f908c69734a42c5234f BUILD_TARGET eigen FIND_TARGET Eigen3::Eigen ) +cmaize_find_or_build_optional_dependency( + sigma + ENABLE_SIGMA + URL github.com/QCUncertainty/sigma + VERSION ad74a8305c615a825919b3ad067fdcdf4e578ed0 + BUILD_TARGET sigma + FIND_TARGET sigma::sigma + CMAKE_ARGS BUILD_TESTING=OFF + ENABLE_EIGEN_SUPPORT=ON +) + cmaize_add_library( ${PROJECT_NAME} SOURCE_DIR "${project_src_dir}" INCLUDE_DIRS "${project_inc_dir}" - DEPENDS utilities parallelzone Boost::boost eigen + DEPENDS utilities parallelzone Boost::boost eigen sigma ) target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_BINARY_DIR}") @@ -80,13 +77,8 @@ if("${BUILD_TESTING}") set(cxx_test_dir "${CMAKE_CURRENT_LIST_DIR}/tests/cxx") set(cxx_unit_test_dir "${cxx_test_dir}/unit_tests/${PROJECT_NAME}") - cmaize_find_or_build_dependency( - Catch2 - URL github.com/catchorg/Catch2 - BUILD_TARGET Catch2 - FIND_TARGET Catch2::Catch2 - VERSION v3.6.0 - ) + include(get_catch2) + cmaize_add_tests( test_unit_${PROJECT_NAME} SOURCE_DIR "${cxx_unit_test_dir}" diff --git a/cmake/config.hpp.in b/cmake/config.hpp.in deleted file mode 100644 index 39bdd22a..00000000 --- a/cmake/config.hpp.in +++ /dev/null @@ -1,13 +0,0 @@ -#pragma once - -#ifdef ENABLE_EIGEN_SUPPORT -#define TENSORWRAPPER_HAS_EIGEN 1 -#else -#define TENSORWRAPPER_HAS_EIGEN 0 -#endif - -namespace tensorwrapper { - -constexpr bool have_eigen() { return TENSORWRAPPER_HAS_EIGEN; } - -} // namespace tensorwrapper diff --git a/include/tensorwrapper/allocator/eigen.hpp b/include/tensorwrapper/allocator/eigen.hpp index 051046e6..63bf053e 100644 --- a/include/tensorwrapper/allocator/eigen.hpp +++ b/include/tensorwrapper/allocator/eigen.hpp @@ -18,6 +18,10 @@ #include #include +#ifdef ENABLE_SIGMA +#include +#endif + namespace tensorwrapper::allocator { /** @brief Used to allocate buffers which rely on Eigen tensors. @@ -264,21 +268,26 @@ class Eigen : public Replicated { // -- Explicit class template declarations // ----------------------------------------------------------------------------- -#define DECLARE_EIGEN_ALLOCATOR(RANK) \ - extern template class Eigen; \ - extern template class Eigen - -DECLARE_EIGEN_ALLOCATOR(0); -DECLARE_EIGEN_ALLOCATOR(1); -DECLARE_EIGEN_ALLOCATOR(2); -DECLARE_EIGEN_ALLOCATOR(3); -DECLARE_EIGEN_ALLOCATOR(4); -DECLARE_EIGEN_ALLOCATOR(5); -DECLARE_EIGEN_ALLOCATOR(6); -DECLARE_EIGEN_ALLOCATOR(7); -DECLARE_EIGEN_ALLOCATOR(8); -DECLARE_EIGEN_ALLOCATOR(9); -DECLARE_EIGEN_ALLOCATOR(10); +#define DECLARE_EIGEN_ALLOCATOR(TYPE) \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen + +DECLARE_EIGEN_ALLOCATOR(float); +DECLARE_EIGEN_ALLOCATOR(double); + +#ifdef ENABLE_SIGMA +DECLARE_EIGEN_ALLOCATOR(sigma::UFloat); +DECLARE_EIGEN_ALLOCATOR(sigma::UDouble); +#endif #undef DECLARE_EIGEN_ALLOCATOR diff --git a/include/tensorwrapper/backends/eigen.hpp b/include/tensorwrapper/backends/eigen.hpp index e3d95395..fa564773 100644 --- a/include/tensorwrapper/backends/eigen.hpp +++ b/include/tensorwrapper/backends/eigen.hpp @@ -15,23 +15,11 @@ */ #pragma once -#include "config.hpp" -#ifdef TENSORWRAPPER_HAS_EIGEN #include -#endif namespace tensorwrapper::eigen { -#ifdef TENSORWRAPPER_HAS_EIGEN - template using data_type = Eigen::Tensor; -#else - -template -using data_type = std::nullptr_t; - -#endif - } // namespace tensorwrapper::eigen diff --git a/include/tensorwrapper/buffer/eigen.hpp b/include/tensorwrapper/buffer/eigen.hpp index d9aca697..94f68dd2 100644 --- a/include/tensorwrapper/buffer/eigen.hpp +++ b/include/tensorwrapper/buffer/eigen.hpp @@ -18,6 +18,10 @@ #include #include +#ifdef ENABLE_SIGMA +#include +#endif + namespace tensorwrapper::buffer { /** @brief A buffer which wraps an Eigen tensor. @@ -144,6 +148,10 @@ class Eigen : public Replicated { * Two Eigen objects are value equal if they both have the same layout and * they both have the same values. * + * @note For tensors where the @p FloatType is an uncertain floating point + * number, the tensors are required to have the same sources of + * uncertainty. + * * @param[in] rhs The object to compare against. * * @return True if *this is value equal to @p rhs and false otherwise. @@ -153,7 +161,7 @@ class Eigen : public Replicated { bool operator==(const Eigen& rhs) const noexcept { if(my_base_type::operator!=(rhs)) return false; eigen::data_type r = (m_tensor_ - rhs.m_tensor_).sum(); - return r() == 0.0; + return r() == FloatType(0.0); } /** @brief Is *this different from @p rhs? @@ -217,21 +225,26 @@ class Eigen : public Replicated { data_type m_tensor_; }; -#define DECLARE_EIGEN_BUFFER(RANK) \ - extern template class Eigen; \ - extern template class Eigen - -DECLARE_EIGEN_BUFFER(0); -DECLARE_EIGEN_BUFFER(1); -DECLARE_EIGEN_BUFFER(2); -DECLARE_EIGEN_BUFFER(3); -DECLARE_EIGEN_BUFFER(4); -DECLARE_EIGEN_BUFFER(5); -DECLARE_EIGEN_BUFFER(6); -DECLARE_EIGEN_BUFFER(7); -DECLARE_EIGEN_BUFFER(8); -DECLARE_EIGEN_BUFFER(9); -DECLARE_EIGEN_BUFFER(10); +#define DECLARE_EIGEN_BUFFER(TYPE) \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen; \ + extern template class Eigen + +DECLARE_EIGEN_BUFFER(float); +DECLARE_EIGEN_BUFFER(double); + +#ifdef ENABLE_SIGMA +DECLARE_EIGEN_BUFFER(sigma::UFloat); +DECLARE_EIGEN_BUFFER(sigma::UDouble); +#endif #undef DECLARE_EIGEN_BUFFER diff --git a/src/tensorwrapper/allocator/eigen.cpp b/src/tensorwrapper/allocator/eigen.cpp index 67503db3..60905765 100644 --- a/src/tensorwrapper/allocator/eigen.cpp +++ b/src/tensorwrapper/allocator/eigen.cpp @@ -107,21 +107,26 @@ typename EIGEN::buffer_base_pointer EIGEN::allocate_(layout_pointer playout) { // -- Explicit class template instantiation -#define DEFINE_EIGEN_ALLOCATOR(RANK) \ - template class Eigen; \ - template class Eigen - -DEFINE_EIGEN_ALLOCATOR(0); -DEFINE_EIGEN_ALLOCATOR(1); -DEFINE_EIGEN_ALLOCATOR(2); -DEFINE_EIGEN_ALLOCATOR(3); -DEFINE_EIGEN_ALLOCATOR(4); -DEFINE_EIGEN_ALLOCATOR(5); -DEFINE_EIGEN_ALLOCATOR(6); -DEFINE_EIGEN_ALLOCATOR(7); -DEFINE_EIGEN_ALLOCATOR(8); -DEFINE_EIGEN_ALLOCATOR(9); -DEFINE_EIGEN_ALLOCATOR(10); +#define DEFINE_EIGEN_ALLOCATOR(TYPE) \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen + +DEFINE_EIGEN_ALLOCATOR(float); +DEFINE_EIGEN_ALLOCATOR(double); + +#ifdef ENABLE_SIGMA +DEFINE_EIGEN_ALLOCATOR(sigma::UFloat); +DEFINE_EIGEN_ALLOCATOR(sigma::UDouble); +#endif #undef DEFINE_EIGEN_ALLOCATOR diff --git a/src/tensorwrapper/buffer/eigen.cpp b/src/tensorwrapper/buffer/eigen.cpp index 9f778512..a0a35761 100644 --- a/src/tensorwrapper/buffer/eigen.cpp +++ b/src/tensorwrapper/buffer/eigen.cpp @@ -234,21 +234,26 @@ TPARAMS typename EIGEN::dsl_reference EIGEN::contraction_( #undef EIGEN #undef TPARAMS -#define DEFINE_EIGEN_BUFFER(RANK) \ - template class Eigen; \ - template class Eigen - -DEFINE_EIGEN_BUFFER(0); -DEFINE_EIGEN_BUFFER(1); -DEFINE_EIGEN_BUFFER(2); -DEFINE_EIGEN_BUFFER(3); -DEFINE_EIGEN_BUFFER(4); -DEFINE_EIGEN_BUFFER(5); -DEFINE_EIGEN_BUFFER(6); -DEFINE_EIGEN_BUFFER(7); -DEFINE_EIGEN_BUFFER(8); -DEFINE_EIGEN_BUFFER(9); -DEFINE_EIGEN_BUFFER(10); +#define DEFINE_EIGEN_BUFFER(TYPE) \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen; \ + template class Eigen + +DEFINE_EIGEN_BUFFER(float); +DEFINE_EIGEN_BUFFER(double); + +#ifdef ENABLE_SIGMA +DEFINE_EIGEN_BUFFER(sigma::UFloat); +DEFINE_EIGEN_BUFFER(sigma::UDouble); +#endif #undef DEFINE_EIGEN_BUFFER diff --git a/tests/cxx/unit_tests/tensorwrapper/allocator/eigen.cpp b/tests/cxx/unit_tests/tensorwrapper/allocator/eigen.cpp index dedd8e79..a0eab4c7 100644 --- a/tests/cxx/unit_tests/tensorwrapper/allocator/eigen.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/allocator/eigen.cpp @@ -22,7 +22,13 @@ using namespace tensorwrapper; -TEMPLATE_TEST_CASE("EigenAllocator", "", float, double) { +#ifdef ENABLE_SIGMA +using types2test = std::tuple; +#else +using types2test = std::tuple; +#endif + +TEMPLATE_LIST_TEST_CASE("EigenAllocator", "", types2test) { using scalar_alloc_type = allocator::Eigen; using vector_alloc_type = allocator::Eigen; using matrix_alloc_type = allocator::Eigen; diff --git a/tests/cxx/unit_tests/tensorwrapper/buffer/buffer_base.cpp b/tests/cxx/unit_tests/tensorwrapper/buffer/buffer_base.cpp index a27a24c3..e5928b61 100644 --- a/tests/cxx/unit_tests/tensorwrapper/buffer/buffer_base.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/buffer/buffer_base.cpp @@ -34,58 +34,56 @@ using namespace buffer; */ TEST_CASE("BufferBase") { - if constexpr(have_eigen()) { - using scalar_buffer = buffer::Eigen; - using vector_buffer = buffer::Eigen; - - typename scalar_buffer::data_type eigen_scalar; - eigen_scalar() = 1.0; - - typename vector_buffer::data_type eigen_vector(2); - eigen_vector(0) = 1.0; - eigen_vector(1) = 2.0; - - auto scalar_layout = testing::scalar_physical(); - auto vector_layout = testing::vector_physical(2); - - vector_buffer defaulted; - scalar_buffer scalar(eigen_scalar, scalar_layout); - vector_buffer vector(eigen_vector, vector_layout); - - BufferBase& defaulted_base = defaulted; - BufferBase& scalar_base = scalar; - BufferBase& vector_base = vector; - - SECTION("has_layout") { - REQUIRE_FALSE(defaulted_base.has_layout()); - REQUIRE(scalar_base.has_layout()); - REQUIRE(vector_base.has_layout()); - } - - SECTION("layout") { - REQUIRE_THROWS_AS(defaulted_base.layout(), std::runtime_error); - REQUIRE(scalar_base.layout().are_equal(scalar_layout)); - REQUIRE(vector_base.layout().are_equal(vector_layout)); - } - - SECTION("operator==") { - // Defaulted layout == defaulted layout - REQUIRE(defaulted_base == scalar_buffer()); - - // Defaulted layout != non-defaulted layout - REQUIRE_FALSE(defaulted_base == scalar_base); - - // Non-defaulted layout same value - REQUIRE(scalar_base == scalar_buffer(eigen_scalar, scalar_layout)); - - // Non-defaulted layout different value - REQUIRE_FALSE(scalar_base == vector_base); - } - - SECTION("operator!=") { - // Just spot check because it negates operator==, which was tested - REQUIRE(defaulted_base != scalar_base); - REQUIRE_FALSE(defaulted_base != scalar_buffer()); - } + using scalar_buffer = buffer::Eigen; + using vector_buffer = buffer::Eigen; + + typename scalar_buffer::data_type eigen_scalar; + eigen_scalar() = 1.0; + + typename vector_buffer::data_type eigen_vector(2); + eigen_vector(0) = 1.0; + eigen_vector(1) = 2.0; + + auto scalar_layout = testing::scalar_physical(); + auto vector_layout = testing::vector_physical(2); + + vector_buffer defaulted; + scalar_buffer scalar(eigen_scalar, scalar_layout); + vector_buffer vector(eigen_vector, vector_layout); + + BufferBase& defaulted_base = defaulted; + BufferBase& scalar_base = scalar; + BufferBase& vector_base = vector; + + SECTION("has_layout") { + REQUIRE_FALSE(defaulted_base.has_layout()); + REQUIRE(scalar_base.has_layout()); + REQUIRE(vector_base.has_layout()); + } + + SECTION("layout") { + REQUIRE_THROWS_AS(defaulted_base.layout(), std::runtime_error); + REQUIRE(scalar_base.layout().are_equal(scalar_layout)); + REQUIRE(vector_base.layout().are_equal(vector_layout)); + } + + SECTION("operator==") { + // Defaulted layout == defaulted layout + REQUIRE(defaulted_base == scalar_buffer()); + + // Defaulted layout != non-defaulted layout + REQUIRE_FALSE(defaulted_base == scalar_base); + + // Non-defaulted layout same value + REQUIRE(scalar_base == scalar_buffer(eigen_scalar, scalar_layout)); + + // Non-defaulted layout different value + REQUIRE_FALSE(scalar_base == vector_base); + } + + SECTION("operator!=") { + // Just spot check because it negates operator==, which was tested + REQUIRE(defaulted_base != scalar_base); + REQUIRE_FALSE(defaulted_base != scalar_buffer()); } } diff --git a/tests/cxx/unit_tests/tensorwrapper/buffer/eigen.cpp b/tests/cxx/unit_tests/tensorwrapper/buffer/eigen.cpp index abfd5c8b..f3cc81db 100644 --- a/tests/cxx/unit_tests/tensorwrapper/buffer/eigen.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/buffer/eigen.cpp @@ -22,6 +22,12 @@ using namespace tensorwrapper; using namespace testing; +#ifdef ENABLE_SIGMA +using types2test = std::tuple; +#else +using types2test = std::tuple; +#endif + namespace { template @@ -29,668 +35,671 @@ void compare_eigen(const LHSType& lhs, const RHSType& rhs) { using r_type = Eigen::Tensor; auto d = lhs - rhs; r_type r = d.sum(); - REQUIRE_THAT(r() + 1.0, Catch::Matchers::WithinAbs(1.0, 1E-6)); + // sigma's types are the only none fundamental options here currently + // This may need to be adjusted if other custom numeric types are added + if constexpr(std::is_fundamental_v) { + REQUIRE_THAT(r() + 1.0, Catch::Matchers::WithinAbs(1.0, 1E-6)); + } else { + REQUIRE_THAT(r().mean() + 1.0, Catch::Matchers::WithinAbs(1.0, 1E-6)); + REQUIRE_THAT(r().sd() + 1.0, Catch::Matchers::WithinAbs(1.0, 1E-6)); + } } } // namespace -TEMPLATE_TEST_CASE("Eigen", "", float, double) { - if constexpr(have_eigen()) { - using scalar_buffer = buffer::Eigen; - using vector_buffer = buffer::Eigen; - using matrix_buffer = buffer::Eigen; - using tensor_buffer = buffer::Eigen; - - typename scalar_buffer::data_type eigen_scalar; - eigen_scalar() = 10.0; - - typename vector_buffer::data_type eigen_vector(2); - eigen_vector(0) = 10.0; - eigen_vector(1) = 20.0; - - typename matrix_buffer::data_type eigen_matrix(2, 3); - eigen_matrix(0, 0) = 10.0; - eigen_matrix(0, 1) = 20.0; - eigen_matrix(0, 2) = 30.0; - eigen_matrix(1, 0) = 40.0; - eigen_matrix(1, 1) = 50.0; - eigen_matrix(1, 2) = 60.0; - - typename tensor_buffer::data_type eigen_tensor(1, 2, 3); - eigen_tensor(0, 0, 0) = 10.0; - eigen_tensor(0, 0, 1) = 20.0; - eigen_tensor(0, 0, 2) = 30.0; - eigen_tensor(0, 1, 0) = 40.0; - eigen_tensor(0, 1, 1) = 50.0; - eigen_tensor(0, 1, 2) = 60.0; - - auto scalar_layout = scalar_physical(); - auto vector_layout = vector_physical(2); - auto matrix_layout = matrix_physical(2, 3); - auto tensor_layout = tensor_physical(1, 2, 3); - - scalar_buffer scalar(eigen_scalar, scalar_layout); - vector_buffer vector(eigen_vector, vector_layout); - matrix_buffer matrix(eigen_matrix, matrix_layout); - tensor_buffer tensor(eigen_tensor, tensor_layout); - - SECTION("ctors, assignment") { - SECTION("value ctor") { - compare_eigen(scalar.value(), eigen_scalar); - REQUIRE(scalar.layout().are_equal(scalar_layout)); - - compare_eigen(vector.value(), eigen_vector); - REQUIRE(vector.layout().are_equal(vector_layout)); - - compare_eigen(matrix.value(), eigen_matrix); - REQUIRE(matrix.layout().are_equal(matrix_layout)); - } - - test_copy_move_ctor_and_assignment(scalar, vector, matrix); - } - - SECTION("value()") { +TEMPLATE_LIST_TEST_CASE("Eigen", "", types2test) { + using scalar_buffer = buffer::Eigen; + using vector_buffer = buffer::Eigen; + using matrix_buffer = buffer::Eigen; + using tensor_buffer = buffer::Eigen; + + typename scalar_buffer::data_type eigen_scalar; + eigen_scalar() = 10.0; + + typename vector_buffer::data_type eigen_vector(2); + eigen_vector(0) = 10.0; + eigen_vector(1) = 20.0; + + typename matrix_buffer::data_type eigen_matrix(2, 3); + eigen_matrix(0, 0) = 10.0; + eigen_matrix(0, 1) = 20.0; + eigen_matrix(0, 2) = 30.0; + eigen_matrix(1, 0) = 40.0; + eigen_matrix(1, 1) = 50.0; + eigen_matrix(1, 2) = 60.0; + + typename tensor_buffer::data_type eigen_tensor(1, 2, 3); + eigen_tensor(0, 0, 0) = 10.0; + eigen_tensor(0, 0, 1) = 20.0; + eigen_tensor(0, 0, 2) = 30.0; + eigen_tensor(0, 1, 0) = 40.0; + eigen_tensor(0, 1, 1) = 50.0; + eigen_tensor(0, 1, 2) = 60.0; + + auto scalar_layout = scalar_physical(); + auto vector_layout = vector_physical(2); + auto matrix_layout = matrix_physical(2, 3); + auto tensor_layout = tensor_physical(1, 2, 3); + + scalar_buffer scalar(eigen_scalar, scalar_layout); + vector_buffer vector(eigen_vector, vector_layout); + matrix_buffer matrix(eigen_matrix, matrix_layout); + tensor_buffer tensor(eigen_tensor, tensor_layout); + + SECTION("ctors, assignment") { + SECTION("value ctor") { compare_eigen(scalar.value(), eigen_scalar); + REQUIRE(scalar.layout().are_equal(scalar_layout)); + compare_eigen(vector.value(), eigen_vector); + REQUIRE(vector.layout().are_equal(vector_layout)); + compare_eigen(matrix.value(), eigen_matrix); + REQUIRE(matrix.layout().are_equal(matrix_layout)); } - SECTION("value() const") { - const auto& cscalar = scalar; - const auto& cvector = vector; - const auto& cmatrix = matrix; - compare_eigen(cscalar.value(), eigen_scalar); - compare_eigen(cvector.value(), eigen_vector); - compare_eigen(cmatrix.value(), eigen_matrix); + test_copy_move_ctor_and_assignment(scalar, vector, matrix); + } + + SECTION("value()") { + compare_eigen(scalar.value(), eigen_scalar); + compare_eigen(vector.value(), eigen_vector); + compare_eigen(matrix.value(), eigen_matrix); + } + + SECTION("value() const") { + const auto& cscalar = scalar; + const auto& cvector = vector; + const auto& cmatrix = matrix; + compare_eigen(cscalar.value(), eigen_scalar); + compare_eigen(cvector.value(), eigen_vector); + compare_eigen(cmatrix.value(), eigen_matrix); + } + + SECTION("operator==") { + // We assume the eigen tensor and the layout objects work. So when + // comparing Eigen objects we have four states: same everything, + // different everything, same tensor different layout, and different + // tensor same layout. + + typename scalar_buffer::data_type eigen_scalar2; + eigen_scalar2() = 10.0; + + // Everything the same + REQUIRE(scalar == scalar_buffer(eigen_scalar2, scalar_layout)); + + SECTION("Different scalar") { + eigen_scalar2() = 2.0; + scalar_buffer scalar2(eigen_scalar2, scalar_layout); + REQUIRE_FALSE(scalar == scalar2); } - SECTION("operator==") { - // We assume the eigen tensor and the layout objects work. So when - // comparing Eigen objects we have four states: same everything, - // different everything, same tensor different layout, and different - // tensor same layout. + SECTION("Different layout") { + scalar_buffer scalar2(eigen_scalar, vector_layout); + REQUIRE_FALSE(scalar == scalar2); + } - typename scalar_buffer::data_type eigen_scalar2; - eigen_scalar2() = 10.0; + SECTION("Different tensor and layout") { + eigen_scalar2() = 2.0; + scalar_buffer scalar2(eigen_scalar2, vector_layout); + REQUIRE_FALSE(scalar == scalar2); + } + } - // Everything the same - REQUIRE(scalar == scalar_buffer(eigen_scalar2, scalar_layout)); + SECTION("operator!=") { + // This just negates operator== so spot-checking is okay - SECTION("Different scalar") { - eigen_scalar2() = 2.0; - scalar_buffer scalar2(eigen_scalar2, scalar_layout); - REQUIRE_FALSE(scalar == scalar2); - } + typename scalar_buffer::data_type eigen_scalar2; + eigen_scalar2() = 10.0; - SECTION("Different layout") { - scalar_buffer scalar2(eigen_scalar, vector_layout); - REQUIRE_FALSE(scalar == scalar2); - } + // Everything the same + scalar_buffer scalar2(eigen_scalar2, scalar_layout); + REQUIRE_FALSE(scalar != scalar2); - SECTION("Different tensor and layout") { - eigen_scalar2() = 2.0; - scalar_buffer scalar2(eigen_scalar2, vector_layout); - REQUIRE_FALSE(scalar == scalar2); - } + eigen_scalar2() = 2.0; + scalar_buffer scalar3(eigen_scalar2, scalar_layout); + REQUIRE(scalar3 != scalar); + } + + SECTION("virtual method overrides") { + using const_reference = + typename scalar_buffer::const_buffer_base_reference; + const_reference pscalar = scalar; + const_reference pvector = vector; + const_reference pmatrix = matrix; + + SECTION("clone") { + REQUIRE(pscalar.clone()->are_equal(pscalar)); + REQUIRE(pvector.clone()->are_equal(pvector)); + REQUIRE(pmatrix.clone()->are_equal(pmatrix)); } - SECTION("operator!=") { - // This just negates operator== so spot-checking is okay + SECTION("are_equal") { + scalar_buffer scalar2(eigen_scalar, scalar_layout); + REQUIRE(pscalar.are_equal(scalar2)); + REQUIRE_FALSE(pmatrix.are_equal(scalar2)); + } + } - typename scalar_buffer::data_type eigen_scalar2; - eigen_scalar2() = 10.0; + SECTION("addition_assignment_") { + SECTION("scalar") { + scalar_buffer scalar2(eigen_scalar, scalar_layout); + scalar2.value()() = 42.0; - // Everything the same - scalar_buffer scalar2(eigen_scalar2, scalar_layout); - REQUIRE_FALSE(scalar != scalar2); + auto s = scalar(""); + auto pscalar2 = &(scalar2.addition_assignment("", s, s)); - eigen_scalar2() = 2.0; - scalar_buffer scalar3(eigen_scalar2, scalar_layout); - REQUIRE(scalar3 != scalar); + scalar_buffer scalar_corr(eigen_scalar, scalar_layout); + scalar_corr.value()() = 20.0; + REQUIRE(pscalar2 == &scalar2); + REQUIRE(scalar2 == scalar_corr); } - SECTION("virtual method overrides") { - using const_reference = - typename scalar_buffer::const_buffer_base_reference; - const_reference pscalar = scalar; - const_reference pvector = vector; - const_reference pmatrix = matrix; + SECTION("vector") { + auto vector2 = testing::eigen_vector(); - SECTION("clone") { - REQUIRE(pscalar.clone()->are_equal(pscalar)); - REQUIRE(pvector.clone()->are_equal(pvector)); - REQUIRE(pmatrix.clone()->are_equal(pmatrix)); - } + auto vi = vector("i"); + auto pvector2 = &(vector2.addition_assignment("i", vi, vi)); - SECTION("are_equal") { - scalar_buffer scalar2(eigen_scalar, scalar_layout); - REQUIRE(pscalar.are_equal(scalar2)); - REQUIRE_FALSE(pmatrix.are_equal(scalar2)); - } + vector_buffer vector_corr(eigen_vector, vector_layout); + vector_corr.value()(0) = 20.0; + vector_corr.value()(1) = 40.0; + + REQUIRE(pvector2 == &vector2); + REQUIRE(vector2 == vector_corr); } - SECTION("addition_assignment_") { - SECTION("scalar") { - scalar_buffer scalar2(eigen_scalar, scalar_layout); - scalar2.value()() = 42.0; + SECTION("matrix : no permutation") { + auto matrix2 = testing::eigen_matrix(); - auto s = scalar(""); - auto pscalar2 = &(scalar2.addition_assignment("", s, s)); + auto mij = matrix("i,j"); + auto pmatrix2 = &(matrix2.addition_assignment("i,j", mij, mij)); - scalar_buffer scalar_corr(eigen_scalar, scalar_layout); - scalar_corr.value()() = 20.0; - REQUIRE(pscalar2 == &scalar2); - REQUIRE(scalar2 == scalar_corr); - } + matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - SECTION("vector") { - auto vector2 = testing::eigen_vector(); + matrix_corr.value()(0, 0) = 20.0; + matrix_corr.value()(0, 1) = 40.0; + matrix_corr.value()(0, 2) = 60.0; + matrix_corr.value()(1, 0) = 80.0; + matrix_corr.value()(1, 1) = 100.0; + matrix_corr.value()(1, 2) = 120.0; - auto vi = vector("i"); - auto pvector2 = &(vector2.addition_assignment("i", vi, vi)); + REQUIRE(pmatrix2 == &matrix2); + REQUIRE(matrix2 == matrix_corr); + } - vector_buffer vector_corr(eigen_vector, vector_layout); - vector_corr.value()(0) = 20.0; - vector_corr.value()(1) = 40.0; + SECTION("matrix: permutations") { + auto matrix2 = testing::eigen_matrix(); + auto l = testing::matrix_physical(3, 2); + std::array p10{1, 0}; + auto eigen_matrix_t = eigen_matrix.shuffle(p10); + matrix_buffer matrix1(eigen_matrix_t, l); - REQUIRE(pvector2 == &vector2); - REQUIRE(vector2 == vector_corr); - } + auto mij = matrix("i,j"); + auto mji = matrix1("j,i"); - SECTION("matrix : no permutation") { - auto matrix2 = testing::eigen_matrix(); + matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - auto mij = matrix("i,j"); - auto pmatrix2 = &(matrix2.addition_assignment("i,j", mij, mij)); + matrix_corr.value()(0, 0) = 20.0; + matrix_corr.value()(0, 1) = 40.0; + matrix_corr.value()(0, 2) = 60.0; + matrix_corr.value()(1, 0) = 80.0; + matrix_corr.value()(1, 1) = 100.0; + matrix_corr.value()(1, 2) = 120.0; - matrix_buffer matrix_corr(eigen_matrix, matrix_layout); + SECTION("permute this") { + matrix2.addition_assignment("j,i", mij, mij); - matrix_corr.value()(0, 0) = 20.0; - matrix_corr.value()(0, 1) = 40.0; - matrix_corr.value()(0, 2) = 60.0; - matrix_corr.value()(1, 0) = 80.0; - matrix_corr.value()(1, 1) = 100.0; - matrix_corr.value()(1, 2) = 120.0; + matrix_buffer corr(eigen_matrix_t, l); + corr.value()(0, 0) = 20.0; + corr.value()(0, 1) = 80.0; + corr.value()(1, 0) = 40.0; + corr.value()(1, 1) = 100.0; + corr.value()(2, 0) = 60.0; + corr.value()(2, 1) = 120.0; + + REQUIRE(matrix2 == corr); + } - REQUIRE(pmatrix2 == &matrix2); + SECTION("permute LHS") { + matrix2.addition_assignment("i,j", mji, mij); REQUIRE(matrix2 == matrix_corr); } - SECTION("matrix: permutations") { - auto matrix2 = testing::eigen_matrix(); - auto l = testing::matrix_physical(3, 2); - std::array p10{1, 0}; - auto eigen_matrix_t = eigen_matrix.shuffle(p10); - matrix_buffer matrix1(eigen_matrix_t, l); - - auto mij = matrix("i,j"); - auto mji = matrix1("j,i"); - - matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - - matrix_corr.value()(0, 0) = 20.0; - matrix_corr.value()(0, 1) = 40.0; - matrix_corr.value()(0, 2) = 60.0; - matrix_corr.value()(1, 0) = 80.0; - matrix_corr.value()(1, 1) = 100.0; - matrix_corr.value()(1, 2) = 120.0; - - SECTION("permute this") { - matrix2.addition_assignment("j,i", mij, mij); - - matrix_buffer corr(eigen_matrix_t, l); - corr.value()(0, 0) = 20.0; - corr.value()(0, 1) = 80.0; - corr.value()(1, 0) = 40.0; - corr.value()(1, 1) = 100.0; - corr.value()(2, 0) = 60.0; - corr.value()(2, 1) = 120.0; - - REQUIRE(matrix2 == corr); - } - - SECTION("permute LHS") { - matrix2.addition_assignment("i,j", mji, mij); - REQUIRE(matrix2 == matrix_corr); - } - - SECTION("permute RHS") { - matrix2.addition_assignment("i,j", mij, mji); - REQUIRE(matrix2 == matrix_corr); - } + SECTION("permute RHS") { + matrix2.addition_assignment("i,j", mij, mji); + REQUIRE(matrix2 == matrix_corr); } + } - SECTION("tensor (must permute all)") { - auto tensor2 = testing::eigen_tensor3(); + SECTION("tensor (must permute all)") { + auto tensor2 = testing::eigen_tensor3(); - std::array p102{1, 0, 2}; - auto l102 = testing::tensor_physical(2, 1, 3); - tensor_buffer tensor102(eigen_tensor.shuffle(p102), l102); + std::array p102{1, 0, 2}; + auto l102 = testing::tensor_physical(2, 1, 3); + tensor_buffer tensor102(eigen_tensor.shuffle(p102), l102); - auto tijk = tensor("i,j,k"); - auto tjik = tensor102("j,i,k"); + auto tijk = tensor("i,j,k"); + auto tjik = tensor102("j,i,k"); - tensor2.addition_assignment("k,j,i", tijk, tjik); + tensor2.addition_assignment("k,j,i", tijk, tjik); - std::array p210{2, 1, 0}; - auto l210 = testing::tensor_physical(3, 2, 1); - tensor_buffer corr(eigen_tensor.shuffle(p210), l210); - corr.value()(0, 0, 0) = 20.0; - corr.value()(0, 1, 0) = 80.0; - corr.value()(1, 0, 0) = 40.0; - corr.value()(1, 1, 0) = 100.0; - corr.value()(2, 0, 0) = 60.0; - corr.value()(2, 1, 0) = 120.0; - REQUIRE(tensor2 == corr); - } + std::array p210{2, 1, 0}; + auto l210 = testing::tensor_physical(3, 2, 1); + tensor_buffer corr(eigen_tensor.shuffle(p210), l210); + corr.value()(0, 0, 0) = 20.0; + corr.value()(0, 1, 0) = 80.0; + corr.value()(1, 0, 0) = 40.0; + corr.value()(1, 1, 0) = 100.0; + corr.value()(2, 0, 0) = 60.0; + corr.value()(2, 1, 0) = 120.0; + REQUIRE(tensor2 == corr); } + } - SECTION("subtraction_assignment_") { - SECTION("scalar") { - scalar_buffer scalar2(eigen_scalar, scalar_layout); - scalar2.value()() = 42.0; + SECTION("subtraction_assignment_") { + SECTION("scalar") { + scalar_buffer scalar2(eigen_scalar, scalar_layout); + scalar2.value()() = 42.0; - auto s = scalar(""); - auto pscalar2 = &(scalar2.subtraction_assignment("", s, s)); + auto s = scalar(""); + auto pscalar2 = &(scalar2.subtraction_assignment("", s, s)); - scalar_buffer scalar_corr(eigen_scalar, scalar_layout); - scalar_corr.value()() = 0.0; - REQUIRE(pscalar2 == &scalar2); - REQUIRE(scalar2 == scalar_corr); - } + scalar_buffer scalar_corr(eigen_scalar, scalar_layout); + scalar_corr.value()() = 0.0; + REQUIRE(pscalar2 == &scalar2); + REQUIRE(scalar2 == scalar_corr); + } - SECTION("vector") { - auto vector2 = testing::eigen_vector(); + SECTION("vector") { + auto vector2 = testing::eigen_vector(); - auto vi = vector("i"); - auto pvector2 = &(vector2.subtraction_assignment("i", vi, vi)); + auto vi = vector("i"); + auto pvector2 = &(vector2.subtraction_assignment("i", vi, vi)); - vector_buffer vector_corr(eigen_vector, vector_layout); - vector_corr.value()(0) = 0.0; - vector_corr.value()(1) = 0.0; + vector_buffer vector_corr(eigen_vector, vector_layout); + vector_corr.value()(0) = 0.0; + vector_corr.value()(1) = 0.0; - REQUIRE(pvector2 == &vector2); - REQUIRE(vector2 == vector_corr); - } + REQUIRE(pvector2 == &vector2); + REQUIRE(vector2 == vector_corr); + } + + SECTION("matrix : no permutation") { + auto matrix2 = testing::eigen_matrix(); + + auto mij = matrix("i,j"); + auto pmatrix2 = &(matrix2.subtraction_assignment("i,j", mij, mij)); + + matrix_buffer matrix_corr(eigen_matrix, matrix_layout); + + matrix_corr.value()(0, 0) = 0.0; + matrix_corr.value()(0, 1) = 0.0; + matrix_corr.value()(0, 2) = 0.0; + matrix_corr.value()(1, 0) = 0.0; + matrix_corr.value()(1, 1) = 0.0; + matrix_corr.value()(1, 2) = 0.0; + + REQUIRE(pmatrix2 == &matrix2); + REQUIRE(matrix2 == matrix_corr); + } + + SECTION("matrix: permutations") { + auto matrix2 = testing::eigen_matrix(); + auto l = testing::matrix_physical(3, 2); + std::array p10{1, 0}; + auto eigen_matrix_t = eigen_matrix.shuffle(p10); + matrix_buffer matrix1(eigen_matrix_t, l); + + auto mij = matrix("i,j"); + auto mji = matrix1("j,i"); - SECTION("matrix : no permutation") { - auto matrix2 = testing::eigen_matrix(); + matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - auto mij = matrix("i,j"); - auto pmatrix2 = - &(matrix2.subtraction_assignment("i,j", mij, mij)); + matrix_corr.value()(0, 0) = 0.0; + matrix_corr.value()(0, 1) = 0.0; + matrix_corr.value()(0, 2) = 0.0; + matrix_corr.value()(1, 0) = 0.0; + matrix_corr.value()(1, 1) = 0.0; + matrix_corr.value()(1, 2) = 0.0; - matrix_buffer matrix_corr(eigen_matrix, matrix_layout); + SECTION("permute this") { + matrix2.subtraction_assignment("j,i", mij, mij); - matrix_corr.value()(0, 0) = 0.0; - matrix_corr.value()(0, 1) = 0.0; - matrix_corr.value()(0, 2) = 0.0; - matrix_corr.value()(1, 0) = 0.0; - matrix_corr.value()(1, 1) = 0.0; - matrix_corr.value()(1, 2) = 0.0; + matrix_buffer corr(eigen_matrix_t, l); + corr.value()(0, 0) = 0.0; + corr.value()(0, 1) = 0.0; + corr.value()(1, 0) = 0.0; + corr.value()(1, 1) = 0.0; + corr.value()(2, 0) = 0.0; + corr.value()(2, 1) = 0.0; - REQUIRE(pmatrix2 == &matrix2); + REQUIRE(matrix2 == corr); + } + + SECTION("permute LHS") { + matrix2.subtraction_assignment("i,j", mji, mij); REQUIRE(matrix2 == matrix_corr); } - SECTION("matrix: permutations") { - auto matrix2 = testing::eigen_matrix(); - auto l = testing::matrix_physical(3, 2); - std::array p10{1, 0}; - auto eigen_matrix_t = eigen_matrix.shuffle(p10); - matrix_buffer matrix1(eigen_matrix_t, l); - - auto mij = matrix("i,j"); - auto mji = matrix1("j,i"); - - matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - - matrix_corr.value()(0, 0) = 0.0; - matrix_corr.value()(0, 1) = 0.0; - matrix_corr.value()(0, 2) = 0.0; - matrix_corr.value()(1, 0) = 0.0; - matrix_corr.value()(1, 1) = 0.0; - matrix_corr.value()(1, 2) = 0.0; - - SECTION("permute this") { - matrix2.subtraction_assignment("j,i", mij, mij); - - matrix_buffer corr(eigen_matrix_t, l); - corr.value()(0, 0) = 0.0; - corr.value()(0, 1) = 0.0; - corr.value()(1, 0) = 0.0; - corr.value()(1, 1) = 0.0; - corr.value()(2, 0) = 0.0; - corr.value()(2, 1) = 0.0; - - REQUIRE(matrix2 == corr); - } - - SECTION("permute LHS") { - matrix2.subtraction_assignment("i,j", mji, mij); - REQUIRE(matrix2 == matrix_corr); - } - - SECTION("permute RHS") { - matrix2.subtraction_assignment("i,j", mij, mji); - REQUIRE(matrix2 == matrix_corr); - } + SECTION("permute RHS") { + matrix2.subtraction_assignment("i,j", mij, mji); + REQUIRE(matrix2 == matrix_corr); } + } - SECTION("tensor (must permute all)") { - auto tensor2 = testing::eigen_tensor3(); + SECTION("tensor (must permute all)") { + auto tensor2 = testing::eigen_tensor3(); - std::array p102{1, 0, 2}; - auto l102 = testing::tensor_physical(2, 1, 3); - tensor_buffer tensor102(eigen_tensor.shuffle(p102), l102); + std::array p102{1, 0, 2}; + auto l102 = testing::tensor_physical(2, 1, 3); + tensor_buffer tensor102(eigen_tensor.shuffle(p102), l102); - auto tijk = tensor("i,j,k"); - auto tjik = tensor102("j,i,k"); + auto tijk = tensor("i,j,k"); + auto tjik = tensor102("j,i,k"); - tensor2.subtraction_assignment("k,j,i", tijk, tjik); + tensor2.subtraction_assignment("k,j,i", tijk, tjik); - std::array p210{2, 1, 0}; - auto l210 = testing::tensor_physical(3, 2, 1); - tensor_buffer corr(eigen_tensor.shuffle(p210), l210); - corr.value()(0, 0, 0) = 0.0; - corr.value()(0, 1, 0) = 0.0; - corr.value()(1, 0, 0) = 0.0; - corr.value()(1, 1, 0) = 0.0; - corr.value()(2, 0, 0) = 0.0; - corr.value()(2, 1, 0) = 0.0; - REQUIRE(tensor2 == corr); - } + std::array p210{2, 1, 0}; + auto l210 = testing::tensor_physical(3, 2, 1); + tensor_buffer corr(eigen_tensor.shuffle(p210), l210); + corr.value()(0, 0, 0) = 0.0; + corr.value()(0, 1, 0) = 0.0; + corr.value()(1, 0, 0) = 0.0; + corr.value()(1, 1, 0) = 0.0; + corr.value()(2, 0, 0) = 0.0; + corr.value()(2, 1, 0) = 0.0; + REQUIRE(tensor2 == corr); } + } - SECTION("multiplication_assignment_") { - // Multiplication just dispatches to hadamard_ or contraction_ - // Here we test the error-handling + SECTION("multiplication_assignment_") { + // Multiplication just dispatches to hadamard_ or contraction_ + // Here we test the error-handling - // Must be either a pure hadamard or a pure contraction - auto matrix2 = testing::eigen_matrix(); - auto mij = matrix("i,j"); + // Must be either a pure hadamard or a pure contraction + auto matrix2 = testing::eigen_matrix(); + auto mij = matrix("i,j"); - REQUIRE_THROWS_AS(matrix2.subtraction_assignment("i", mij, mij), - std::runtime_error); + REQUIRE_THROWS_AS(matrix2.subtraction_assignment("i", mij, mij), + std::runtime_error); + } + + SECTION("permute_assignment_") { + SECTION("scalar") { + auto scalar2 = testing::eigen_scalar(); + scalar2.value()() = 42.0; + + auto s = scalar(""); + auto pscalar2 = &(scalar2.permute_assignment("", s)); + REQUIRE(pscalar2 == &scalar2); + REQUIRE(scalar2 == scalar); } - SECTION("permute_assignment_") { - SECTION("scalar") { - auto scalar2 = testing::eigen_scalar(); - scalar2.value()() = 42.0; + SECTION("vector") { + auto vector2 = testing::eigen_vector(); - auto s = scalar(""); - auto pscalar2 = &(scalar2.permute_assignment("", s)); - REQUIRE(pscalar2 == &scalar2); - REQUIRE(scalar2 == scalar); - } + auto vi = vector("i"); + auto pvector2 = &(vector2.permute_assignment("i", vi)); - SECTION("vector") { - auto vector2 = testing::eigen_vector(); + REQUIRE(pvector2 == &vector2); + REQUIRE(vector2 == vector); + } + + SECTION("matrix : no permutation") { + auto matrix2 = testing::eigen_matrix(); - auto vi = vector("i"); - auto pvector2 = &(vector2.permute_assignment("i", vi)); + auto mij = matrix("i,j"); + auto pmatrix2 = &(matrix2.permute_assignment("i,j", mij)); - REQUIRE(pvector2 == &vector2); - REQUIRE(vector2 == vector); - } + REQUIRE(pmatrix2 == &matrix2); + REQUIRE(matrix2 == matrix); + } - SECTION("matrix : no permutation") { - auto matrix2 = testing::eigen_matrix(); + SECTION("matrix: permutation") { + auto matrix2 = testing::eigen_matrix(); + auto p = &(matrix2.permute_assignment("j,i", matrix("i,j"))); + + auto corr = testing::eigen_matrix(3, 2); + corr.value()(0, 0) = 10.0; + corr.value()(1, 0) = 20.0; + corr.value()(2, 0) = 30.0; + corr.value()(0, 1) = 40.0; + corr.value()(1, 1) = 50.0; + corr.value()(2, 1) = 60.0; + REQUIRE(p == &matrix2); + compare_eigen(corr.value(), matrix2.value()); + } + } - auto mij = matrix("i,j"); - auto pmatrix2 = &(matrix2.permute_assignment("i,j", mij)); + SECTION("scalar_multiplication_") { + SECTION("scalar") { + auto scalar2 = testing::eigen_scalar(); + scalar2.value()() = 42.0; - REQUIRE(pmatrix2 == &matrix2); - REQUIRE(matrix2 == matrix); - } + auto s = scalar(""); + auto pscalar2 = &(scalar2.scalar_multiplication("", 2.0, s)); - SECTION("matrix: permutation") { - auto matrix2 = testing::eigen_matrix(); - auto p = &(matrix2.permute_assignment("j,i", matrix("i,j"))); - - auto corr = testing::eigen_matrix(3, 2); - corr.value()(0, 0) = 10.0; - corr.value()(1, 0) = 20.0; - corr.value()(2, 0) = 30.0; - corr.value()(0, 1) = 40.0; - corr.value()(1, 1) = 50.0; - corr.value()(2, 1) = 60.0; - REQUIRE(p == &matrix2); - compare_eigen(corr.value(), matrix2.value()); - } + auto corr = testing::eigen_scalar(); + corr.value()() = 20.0; + REQUIRE(pscalar2 == &scalar2); + REQUIRE(scalar2 == corr); } - SECTION("scalar_multiplication_") { - SECTION("scalar") { - auto scalar2 = testing::eigen_scalar(); - scalar2.value()() = 42.0; + SECTION("vector") { + auto vector2 = testing::eigen_vector(); - auto s = scalar(""); - auto pscalar2 = &(scalar2.scalar_multiplication("", 2.0, s)); + auto vi = vector("i"); + auto pvector2 = &(vector2.scalar_multiplication("i", 2.0, vi)); - auto corr = testing::eigen_scalar(); - corr.value()() = 20.0; - REQUIRE(pscalar2 == &scalar2); - REQUIRE(scalar2 == corr); - } + auto corr = testing::eigen_vector(2); + corr.value()(0) = 20.0; + corr.value()(1) = 40.0; - SECTION("vector") { - auto vector2 = testing::eigen_vector(); + REQUIRE(pvector2 == &vector2); + REQUIRE(vector2 == corr); + } - auto vi = vector("i"); - auto pvector2 = &(vector2.scalar_multiplication("i", 2.0, vi)); + SECTION("matrix : no permutation") { + auto matrix2 = testing::eigen_matrix(); - auto corr = testing::eigen_vector(2); - corr.value()(0) = 20.0; - corr.value()(1) = 40.0; + auto mij = matrix("i,j"); + auto p = &(matrix2.scalar_multiplication("i,j", 2.0, mij)); + + auto corr = testing::eigen_matrix(2, 3); + corr.value()(0, 0) = 20.0; + corr.value()(0, 1) = 40.0; + corr.value()(0, 2) = 60.0; + corr.value()(1, 0) = 80.0; + corr.value()(1, 1) = 100.0; + corr.value()(1, 2) = 120.0; + + REQUIRE(p == &matrix2); + REQUIRE(matrix2 == corr); + } - REQUIRE(pvector2 == &vector2); - REQUIRE(vector2 == corr); - } + SECTION("matrix: permutation") { + auto matrix2 = testing::eigen_matrix(); + auto mij = matrix("i,j"); + auto p = &(matrix2.scalar_multiplication("j,i", 2.0, mij)); + + auto corr = testing::eigen_matrix(3, 2); + corr.value()(0, 0) = 20.0; + corr.value()(1, 0) = 40.0; + corr.value()(2, 0) = 60.0; + corr.value()(0, 1) = 80.0; + corr.value()(1, 1) = 100.0; + corr.value()(2, 1) = 120.0; + REQUIRE(p == &matrix2); + compare_eigen(corr.value(), matrix2.value()); + } + } - SECTION("matrix : no permutation") { - auto matrix2 = testing::eigen_matrix(); + SECTION("hadamard_") { + SECTION("scalar") { + scalar_buffer scalar2(eigen_scalar, scalar_layout); + scalar2.value()() = 42.0; - auto mij = matrix("i,j"); - auto p = &(matrix2.scalar_multiplication("i,j", 2.0, mij)); + auto s = scalar(""); + auto pscalar2 = &(scalar2.multiplication_assignment("", s, s)); - auto corr = testing::eigen_matrix(2, 3); - corr.value()(0, 0) = 20.0; - corr.value()(0, 1) = 40.0; - corr.value()(0, 2) = 60.0; - corr.value()(1, 0) = 80.0; - corr.value()(1, 1) = 100.0; - corr.value()(1, 2) = 120.0; + scalar_buffer scalar_corr(eigen_scalar, scalar_layout); + scalar_corr.value()() = 100.0; + REQUIRE(pscalar2 == &scalar2); + REQUIRE(scalar2 == scalar_corr); + } - REQUIRE(p == &matrix2); - REQUIRE(matrix2 == corr); - } + SECTION("vector") { + auto vector2 = testing::eigen_vector(); - SECTION("matrix: permutation") { - auto matrix2 = testing::eigen_matrix(); - auto mij = matrix("i,j"); - auto p = &(matrix2.scalar_multiplication("j,i", 2.0, mij)); + auto vi = vector("i"); + auto pvector2 = &(vector2.multiplication_assignment("i", vi, vi)); - auto corr = testing::eigen_matrix(3, 2); - corr.value()(0, 0) = 20.0; - corr.value()(1, 0) = 40.0; - corr.value()(2, 0) = 60.0; - corr.value()(0, 1) = 80.0; - corr.value()(1, 1) = 100.0; - corr.value()(2, 1) = 120.0; - REQUIRE(p == &matrix2); - compare_eigen(corr.value(), matrix2.value()); - } + vector_buffer vector_corr(eigen_vector, vector_layout); + vector_corr.value()(0) = 100.0; + vector_corr.value()(1) = 400.0; + + REQUIRE(pvector2 == &vector2); + REQUIRE(vector2 == vector_corr); } - SECTION("hadamard_") { - SECTION("scalar") { - scalar_buffer scalar2(eigen_scalar, scalar_layout); - scalar2.value()() = 42.0; + SECTION("matrix : no permutation") { + auto matrix2 = testing::eigen_matrix(); + + auto mij = matrix("i,j"); + auto pmatrix2 = + &(matrix2.multiplication_assignment("i,j", mij, mij)); - auto s = scalar(""); - auto pscalar2 = &(scalar2.multiplication_assignment("", s, s)); + matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - scalar_buffer scalar_corr(eigen_scalar, scalar_layout); - scalar_corr.value()() = 100.0; - REQUIRE(pscalar2 == &scalar2); - REQUIRE(scalar2 == scalar_corr); - } + matrix_corr.value()(0, 0) = 100.0; + matrix_corr.value()(0, 1) = 400.0; + matrix_corr.value()(0, 2) = 900.0; + matrix_corr.value()(1, 0) = 1600.0; + matrix_corr.value()(1, 1) = 2500.0; + matrix_corr.value()(1, 2) = 3600.0; - SECTION("vector") { - auto vector2 = testing::eigen_vector(); + REQUIRE(pmatrix2 == &matrix2); + REQUIRE(matrix2 == matrix_corr); + } - auto vi = vector("i"); - auto pvector2 = - &(vector2.multiplication_assignment("i", vi, vi)); + SECTION("matrix: permutations") { + auto matrix2 = testing::eigen_matrix(); + auto l = testing::matrix_physical(3, 2); + std::array p10{1, 0}; + auto eigen_matrix_t = eigen_matrix.shuffle(p10); + matrix_buffer matrix1(eigen_matrix_t, l); - vector_buffer vector_corr(eigen_vector, vector_layout); - vector_corr.value()(0) = 100.0; - vector_corr.value()(1) = 400.0; + auto mij = matrix("i,j"); + auto mji = matrix1("j,i"); - REQUIRE(pvector2 == &vector2); - REQUIRE(vector2 == vector_corr); - } + matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - SECTION("matrix : no permutation") { - auto matrix2 = testing::eigen_matrix(); + matrix_corr.value()(0, 0) = 100.0; + matrix_corr.value()(0, 1) = 400.0; + matrix_corr.value()(0, 2) = 900.0; + matrix_corr.value()(1, 0) = 1600.0; + matrix_corr.value()(1, 1) = 2500.0; + matrix_corr.value()(1, 2) = 3600.0; - auto mij = matrix("i,j"); - auto pmatrix2 = - &(matrix2.multiplication_assignment("i,j", mij, mij)); + SECTION("permute this") { + matrix2.multiplication_assignment("j,i", mij, mij); - matrix_buffer matrix_corr(eigen_matrix, matrix_layout); + matrix_buffer corr(eigen_matrix_t, l); + corr.value()(0, 0) = 100.0; + corr.value()(0, 1) = 1600.0; + corr.value()(1, 0) = 400.0; + corr.value()(1, 1) = 2500.0; + corr.value()(2, 0) = 900.0; + corr.value()(2, 1) = 3600.0; - matrix_corr.value()(0, 0) = 100.0; - matrix_corr.value()(0, 1) = 400.0; - matrix_corr.value()(0, 2) = 900.0; - matrix_corr.value()(1, 0) = 1600.0; - matrix_corr.value()(1, 1) = 2500.0; - matrix_corr.value()(1, 2) = 3600.0; + REQUIRE(matrix2 == corr); + } - REQUIRE(pmatrix2 == &matrix2); + SECTION("permute LHS") { + matrix2.multiplication_assignment("i,j", mji, mij); REQUIRE(matrix2 == matrix_corr); } - SECTION("matrix: permutations") { - auto matrix2 = testing::eigen_matrix(); - auto l = testing::matrix_physical(3, 2); - std::array p10{1, 0}; - auto eigen_matrix_t = eigen_matrix.shuffle(p10); - matrix_buffer matrix1(eigen_matrix_t, l); - - auto mij = matrix("i,j"); - auto mji = matrix1("j,i"); - - matrix_buffer matrix_corr(eigen_matrix, matrix_layout); - - matrix_corr.value()(0, 0) = 100.0; - matrix_corr.value()(0, 1) = 400.0; - matrix_corr.value()(0, 2) = 900.0; - matrix_corr.value()(1, 0) = 1600.0; - matrix_corr.value()(1, 1) = 2500.0; - matrix_corr.value()(1, 2) = 3600.0; - - SECTION("permute this") { - matrix2.multiplication_assignment("j,i", mij, mij); - - matrix_buffer corr(eigen_matrix_t, l); - corr.value()(0, 0) = 100.0; - corr.value()(0, 1) = 1600.0; - corr.value()(1, 0) = 400.0; - corr.value()(1, 1) = 2500.0; - corr.value()(2, 0) = 900.0; - corr.value()(2, 1) = 3600.0; - - REQUIRE(matrix2 == corr); - } - - SECTION("permute LHS") { - matrix2.multiplication_assignment("i,j", mji, mij); - REQUIRE(matrix2 == matrix_corr); - } - - SECTION("permute RHS") { - matrix2.multiplication_assignment("i,j", mij, mji); - REQUIRE(matrix2 == matrix_corr); - } + SECTION("permute RHS") { + matrix2.multiplication_assignment("i,j", mij, mji); + REQUIRE(matrix2 == matrix_corr); } + } - SECTION("tensor (must permute all)") { - auto tensor2 = testing::eigen_tensor3(); + SECTION("tensor (must permute all)") { + auto tensor2 = testing::eigen_tensor3(); - std::array p102{1, 0, 2}; - auto l102 = testing::tensor_physical(2, 1, 3); - tensor_buffer tensor102(eigen_tensor.shuffle(p102), l102); + std::array p102{1, 0, 2}; + auto l102 = testing::tensor_physical(2, 1, 3); + tensor_buffer tensor102(eigen_tensor.shuffle(p102), l102); - auto tijk = tensor("i,j,k"); - auto tjik = tensor102("j,i,k"); + auto tijk = tensor("i,j,k"); + auto tjik = tensor102("j,i,k"); - tensor2.multiplication_assignment("k,j,i", tijk, tjik); + tensor2.multiplication_assignment("k,j,i", tijk, tjik); - std::array p210{2, 1, 0}; - auto l210 = testing::tensor_physical(3, 2, 1); - tensor_buffer corr(eigen_tensor.shuffle(p210), l210); - corr.value()(0, 0, 0) = 100.0; - corr.value()(0, 1, 0) = 1600.0; - corr.value()(1, 0, 0) = 400.0; - corr.value()(1, 1, 0) = 2500.0; - corr.value()(2, 0, 0) = 900.0; - corr.value()(2, 1, 0) = 3600.0; - REQUIRE(tensor2 == corr); - } + std::array p210{2, 1, 0}; + auto l210 = testing::tensor_physical(3, 2, 1); + tensor_buffer corr(eigen_tensor.shuffle(p210), l210); + corr.value()(0, 0, 0) = 100.0; + corr.value()(0, 1, 0) = 1600.0; + corr.value()(1, 0, 0) = 400.0; + corr.value()(1, 1, 0) = 2500.0; + corr.value()(2, 0, 0) = 900.0; + corr.value()(2, 1, 0) = 3600.0; + REQUIRE(tensor2 == corr); } + } - SECTION("contraction_") { - auto vi = vector("i"); - auto mij = matrix("i,j"); - auto mik = matrix("i,k"); - auto mjk = matrix("j,k"); + SECTION("contraction_") { + auto vi = vector("i"); + auto mij = matrix("i,j"); + auto mik = matrix("i,k"); + auto mjk = matrix("j,k"); - SECTION("vector with vector") { - auto p = &(scalar.multiplication_assignment("", vi, vi)); + SECTION("vector with vector") { + auto p = &(scalar.multiplication_assignment("", vi, vi)); - auto scalar_corr = testing::eigen_scalar(); - scalar_corr.value()() = 500.0; // 10*10 + 20*20 - REQUIRE(p == &scalar); - compare_eigen(scalar_corr.value(), scalar.value()); - } + auto scalar_corr = testing::eigen_scalar(); + scalar_corr.value()() = 500.0; // 10*10 + 20*20 + REQUIRE(p == &scalar); + compare_eigen(scalar_corr.value(), scalar.value()); + } - SECTION("ij,ij->") { - auto p = &(scalar.multiplication_assignment("", mij, mij)); + SECTION("ij,ij->") { + auto p = &(scalar.multiplication_assignment("", mij, mij)); - auto scalar_corr = testing::eigen_scalar(); - scalar_corr.value()() = 9100.0; // 1400 + 7700 - REQUIRE(p == &scalar); - compare_eigen(scalar_corr.value(), scalar.value()); - } + auto scalar_corr = testing::eigen_scalar(); + scalar_corr.value()() = 9100.0; // 1400 + 7700 + REQUIRE(p == &scalar); + compare_eigen(scalar_corr.value(), scalar.value()); + } - SECTION("ki,kj->ij") { - auto buffer2 = testing::eigen_matrix(); - auto p = &(buffer2.multiplication_assignment("i,j", mik, mjk)); + SECTION("ki,kj->ij") { + auto buffer2 = testing::eigen_matrix(); + auto p = &(buffer2.multiplication_assignment("i,j", mik, mjk)); - auto matrix_corr = testing::eigen_matrix(2, 2); - matrix_corr.value()(0, 0) = 1400.0; // 100 + 400 + 900 - matrix_corr.value()(0, 1) = 3200.0; // 400 + 1000 + 1800 - matrix_corr.value()(1, 0) = 3200.0; // 400 + 1000 + 1800 - matrix_corr.value()(1, 1) = 7700.0; // 1600 + 2500 + 3600 + auto matrix_corr = testing::eigen_matrix(2, 2); + matrix_corr.value()(0, 0) = 1400.0; // 100 + 400 + 900 + matrix_corr.value()(0, 1) = 3200.0; // 400 + 1000 + 1800 + matrix_corr.value()(1, 0) = 3200.0; // 400 + 1000 + 1800 + matrix_corr.value()(1, 1) = 7700.0; // 1600 + 2500 + 3600 - REQUIRE(p == &buffer2); - compare_eigen(matrix_corr.value(), buffer2.value()); - } + REQUIRE(p == &buffer2); + compare_eigen(matrix_corr.value(), buffer2.value()); + } - SECTION("ij,i->j") { - auto buffer1 = testing::eigen_vector(); - auto p = &(buffer1.multiplication_assignment("j", mij, vi)); + SECTION("ij,i->j") { + auto buffer1 = testing::eigen_vector(); + auto p = &(buffer1.multiplication_assignment("j", mij, vi)); - auto vector_corr = testing::eigen_vector(3); - vector_corr.value()(0) = 900.0; // 10(10) + 20(40) - vector_corr.value()(1) = 1200.0; // 10(20) + 20(50) - vector_corr.value()(2) = 1500.0; // 10(30) + 20(60) - REQUIRE(p == &buffer1); - compare_eigen(vector_corr.value(), buffer1.value()); - } + auto vector_corr = testing::eigen_vector(3); + vector_corr.value()(0) = 900.0; // 10(10) + 20(40) + vector_corr.value()(1) = 1200.0; // 10(20) + 20(50) + vector_corr.value()(2) = 1500.0; // 10(30) + 20(60) + REQUIRE(p == &buffer1); + compare_eigen(vector_corr.value(), buffer1.value()); } } } diff --git a/tests/cxx/unit_tests/tensorwrapper/buffer/eigen_contraction.cpp b/tests/cxx/unit_tests/tensorwrapper/buffer/eigen_contraction.cpp index c75fe2be..eef3f23a 100644 --- a/tests/cxx/unit_tests/tensorwrapper/buffer/eigen_contraction.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/buffer/eigen_contraction.cpp @@ -21,7 +21,13 @@ using namespace tensorwrapper; using namespace buffer; -TEMPLATE_TEST_CASE("eigen_contraction", "", float, double) { +#ifdef ENABLE_SIGMA +using types2test = std::tuple; +#else +using types2test = std::tuple; +#endif + +TEMPLATE_LIST_TEST_CASE("eigen_contraction", "", types2test) { using float_t = TestType; using label_type = typename BufferBase::label_type;