diff --git a/src/scf/driver/scf_loop.cpp b/src/scf/driver/scf_loop.cpp index c005c96..f4e97d8 100644 --- a/src/scf/driver/scf_loop.cpp +++ b/src/scf/driver/scf_loop.cpp @@ -107,8 +107,8 @@ MODULE_RUN(SCFLoop) { } auto e_nuclear = V_nn_mod.run_as(qs_lhs, qs_rhs); - wf_type psi_old = psi0; - double e_old = 0; + wf_type psi_old = psi0; + simde::type::tensor e_old; const unsigned int max_iter = 3; unsigned int iter = 0; @@ -154,8 +154,23 @@ MODULE_RUN(SCFLoop) { psi_old = new_psi; ++iter; } - auto rv = results(); - return pt::wrap_results(rv, e_old + e_nuclear, psi_old); + simde::type::tensor e_total; + + // e_nuclear is a double. This hack converts it to udouble (if needed) + tensorwrapper::allocator::Eigen dalloc(get_runtime()); + using tensorwrapper::types::udouble; + tensorwrapper::allocator::Eigen ualloc(get_runtime()); + + if(ualloc.can_rebind(e_old.buffer())) { + simde::type::tensor temp(e_old); + auto val = dalloc.rebind(e_nuclear.buffer()).at(); + ualloc.rebind(temp.buffer()).at() = val; + e_nuclear = temp; + } + + e_total("") = e_old("") + e_nuclear(""); + auto rv = results(); + return pt::wrap_results(rv, e_total, psi_old); } } // namespace scf::driver \ No newline at end of file diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp index f84cf14..31014dc 100644 --- a/src/scf/eigen_solver/eigen_generalized.cpp +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -20,6 +20,60 @@ namespace scf::eigen_solver { +namespace { +struct Kernel { + template + auto run(const tensorwrapper::buffer::BufferBase& A, + const tensorwrapper::buffer::BufferBase& B, + parallelzone::runtime::RuntimeView& rv) { + // Convert to Eigen buffers + tensorwrapper::allocator::Eigen allocator(rv); + const auto& eigen_A = allocator.rebind(A); + const auto& eigen_B = allocator.rebind(B); + + // Wrap the tensors in Eigen::Map objects to avoid copy + const auto* pA = eigen_A.data(); + const auto* pB = eigen_B.data(); + const auto& shape_A = eigen_A.layout().shape().as_smooth(); + auto rows = shape_A.extent(0); + auto cols = shape_A.extent(1); + + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using matrix_type = Eigen::Matrix; + using map_type = Eigen::Map; + + map_type A_map(pA, rows, cols); + map_type B_map(pB, rows, cols); + + // Compute + Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); + auto eigen_values = ges.eigenvalues(); + auto eigen_vectors = ges.eigenvectors(); + + // Wrap in TensorWrapper Tensor + tensorwrapper::shape::Smooth vector_shape{rows}; + tensorwrapper::shape::Smooth matrix_shape{rows, cols}; + tensorwrapper::layout::Physical vector_layout(vector_shape); + tensorwrapper::layout::Physical matrix_layout(matrix_shape); + + auto pvalues_buffer = allocator.allocate(vector_layout); + auto pvectors_buffer = allocator.allocate(matrix_layout); + + for(auto i = 0; i < rows; ++i) { + pvalues_buffer->at(i) = eigen_values(i); + for(auto j = 0; j < cols; ++j) { + pvectors_buffer->at(i, j) = eigen_vectors(i, j); + } + } + + simde::type::tensor values(vector_shape, std::move(pvalues_buffer)); + simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer)); + return std::make_pair(values, vectors); + } +}; +} // namespace + using pt = simde::GeneralizedEigenSolve; const auto desc = R"( @@ -37,43 +91,13 @@ MODULE_CTOR(EigenGeneralized) { MODULE_RUN(EigenGeneralized) { auto&& [A, B] = pt::unwrap_inputs(inputs); - // Convert to Eigen buffers - tensorwrapper::allocator::Eigen allocator(get_runtime()); - const auto& eigen_A = allocator.rebind(A.buffer()); - const auto& eigen_B = allocator.rebind(B.buffer()); - - // Wrap the tensors in Eigen::Map objects to avoid copy - const auto* pA = eigen_A.data(); - const auto* pB = eigen_B.data(); - const auto& shape_A = eigen_A.layout().shape().as_smooth(); - auto rows = shape_A.extent(0); - auto cols = shape_A.extent(1); - Eigen::Map A_map(pA, rows, cols); - Eigen::Map B_map(pB, rows, cols); - - // Compute - Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); - auto eigen_values = ges.eigenvalues(); - auto eigen_vectors = ges.eigenvectors(); - - // Wrap in TensorWrapper Tensor - tensorwrapper::shape::Smooth vector_shape{rows}; - tensorwrapper::shape::Smooth matrix_shape{rows, cols}; - tensorwrapper::layout::Physical vector_layout(vector_shape); - tensorwrapper::layout::Physical matrix_layout(matrix_shape); - - auto pvalues_buffer = allocator.allocate(vector_layout); - auto pvectors_buffer = allocator.allocate(matrix_layout); - - for(auto i = 0; i < rows; ++i) { - pvalues_buffer->at(i) = eigen_values(i); - for(auto j = 0; j < cols; ++j) { - pvectors_buffer->at(i, j) = eigen_vectors(i, j); - } - } + using tensorwrapper::utilities::floating_point_dispatch; - simde::type::tensor values(vector_shape, std::move(pvalues_buffer)); - simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer)); + auto r = get_runtime(); + Kernel k; + const auto& A_buffer = A.buffer(); + const auto& B_buffer = B.buffer(); + auto [values, vectors] = floating_point_dispatch(k, A_buffer, B_buffer, r); auto rv = results(); return pt::wrap_results(rv, values, vectors); diff --git a/src/scf/matrix_builder/density_matrix.cpp b/src/scf/matrix_builder/density_matrix.cpp index faa33f6..ecf9fab 100644 --- a/src/scf/matrix_builder/density_matrix.cpp +++ b/src/scf/matrix_builder/density_matrix.cpp @@ -23,7 +23,40 @@ namespace { const auto desc = R"( )"; -} +struct Kernel { + template + auto run(const tensorwrapper::buffer::BufferBase& c, std::size_t n_aos, + std::size_t n_occ) { + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using allocator_type = tensorwrapper::allocator::Eigen; + using tensor_type = Eigen::Matrix; + using map_type = Eigen::Map; + using const_map_type = Eigen::Map; + auto rv = c.allocator().runtime(); + allocator_type alloc(rv); + + tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); + tensorwrapper::layout::Physical l(p_shape); + auto pp_buffer = alloc.allocate(l); + + // Step 2: Grab the orbitals in the ensemble + auto& c_buffer = alloc.rebind(c); + + const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); + map_type p_eigen(pp_buffer->data(), n_aos, n_aos); + auto slice = c_eigen.block(0, 0, n_aos, n_occ); + + // Step 3: CC_dagger + using index_pair_t = Eigen::IndexPair; + Eigen::array modes{index_pair_t(1, 1)}; + p_eigen = slice * slice.transpose(); + + return simde::type::tensor(p_shape, std::move(pp_buffer)); + } +}; + +} // namespace using pt = simde::aos_rho_e_aos; @@ -63,35 +96,10 @@ MODULE_RUN(DensityMatrix) { "Please shuffle your orbitals so that the ensemble is a " "contiguous slice of the orbitals."); - using allocator_type = tensorwrapper::allocator::Eigen; - using tensor_type = Eigen::Tensor; - using const_tensor_type = Eigen::Tensor; - using map_type = Eigen::TensorMap; - using const_map_type = Eigen::TensorMap; - - allocator_type alloc(get_runtime()); - - tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); - tensorwrapper::layout::Physical l(p_shape); - auto pp_buffer = alloc.allocate(l); - - // Step 2: Grab the orbitals in the ensemble - const auto& c_buffer = alloc.rebind(c.buffer()); - const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); - map_type p_eigen(pp_buffer->data(), n_aos, n_aos); - - using slice_t = Eigen::array; - slice_t offsets{0, 0}; - slice_t extents{n_aos, Eigen::Index(participants.size())}; - auto slice = c_eigen.slice(offsets, extents); - - // Step 3: CC_dagger - using index_pair_t = Eigen::IndexPair; - Eigen::array modes{index_pair_t(1, 1)}; - p_eigen = slice.contract(slice, modes); - - simde::type::tensor p(p_shape, std::move(pp_buffer)); - + // TODO: The need to dispatch like this goes away when TW supports slicing + using tensorwrapper::utilities::floating_point_dispatch; + Kernel k; + auto p = floating_point_dispatch(k, c.buffer(), n_aos, participants.size()); auto rv = results(); return pt::wrap_results(rv, p); } diff --git a/src/scf/matrix_builder/determinant_driver.cpp b/src/scf/matrix_builder/determinant_driver.cpp index e0516d5..fd7c276 100644 --- a/src/scf/matrix_builder/determinant_driver.cpp +++ b/src/scf/matrix_builder/determinant_driver.cpp @@ -103,7 +103,6 @@ MODULE_CTOR(DeterminantDriver) { } MODULE_RUN(DeterminantDriver) { - using float_type = double; using wf_type = simde::type::rscf_wf; const auto&& [braket] = pt::unwrap_inputs(inputs); const auto& bra = braket.bra(); @@ -130,11 +129,8 @@ MODULE_RUN(DeterminantDriver) { tensorwrapper::Tensor x; x("") = rho("i,j") * t("i,j"); - using allocator_type = tensorwrapper::allocator::Eigen; - auto& x_buffer = allocator_type::rebind(x.buffer()); - auto rv = results(); - return pt::wrap_results(rv, x_buffer.at()); + return pt::wrap_results(rv, x); } } // namespace scf::matrix_builder diff --git a/src/scf/matrix_builder/electronic_energy.cpp b/src/scf/matrix_builder/electronic_energy.cpp index 66b0fc6..b456549 100644 --- a/src/scf/matrix_builder/electronic_energy.cpp +++ b/src/scf/matrix_builder/electronic_energy.cpp @@ -48,7 +48,7 @@ MODULE_RUN(ElectronicEnergy) { const auto& ket = H_ij.ket(); auto& mod = submods.at("determinant driver"); - double e = 0.0; + simde::type::tensor e; auto n_ops = H.size(); for(decltype(n_ops) i = 0; i < n_ops; ++i) { @@ -57,7 +57,11 @@ MODULE_RUN(ElectronicEnergy) { chemist::braket::BraKet O_ij(bra, O_i, ket); const auto o = mod.run_as>(O_ij); - e += ci * o; + simde::type::tensor temp; + temp("") = o("") * ci; + if(i == 0) + e = temp; + else { e("") = e("") + temp(""); } } auto rv = results(); diff --git a/tests/cxx/integration_tests/coulombs_law.cpp b/tests/cxx/integration_tests/coulombs_law.cpp index bf7525c..d8c73fd 100644 --- a/tests/cxx/integration_tests/coulombs_law.cpp +++ b/tests/cxx/integration_tests/coulombs_law.cpp @@ -19,9 +19,15 @@ using pt = simde::charge_charge_interaction; using Catch::Matchers::WithinAbs; -TEST_CASE("CoulombsLaw") { - auto mm = test_scf::load_modules(); - auto& mod = mm.at("Coulomb's Law"); +TEMPLATE_LIST_TEST_CASE("CoulombsLaw", "", test_scf::float_types) { + using float_type = double; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Coulomb's Law"); + + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; auto h2_nuclei = test_scf::make_h2(); // TODO: Conversions are missing in Chemist. Use those when they're in place @@ -30,5 +36,8 @@ TEST_CASE("CoulombsLaw") { qs.push_back(nucleus.as_point_charge()); auto e_nuclear = mod.run_as(qs, qs); - REQUIRE_THAT(e_nuclear, WithinAbs(0.71510297482837526, 1E-6)); + + pcorr->at() = 0.71510297482837526; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e_nuclear, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/driver/scf_driver.cpp b/tests/cxx/integration_tests/driver/scf_driver.cpp index 417f187..ae80976 100644 --- a/tests/cxx/integration_tests/driver/scf_driver.cpp +++ b/tests/cxx/integration_tests/driver/scf_driver.cpp @@ -16,15 +16,21 @@ #include "../integration_tests.hpp" -using Catch::Matchers::WithinAbs; - using pt = simde::AOEnergy; -TEST_CASE("SCFDriver") { - auto mm = test_scf::load_modules(); - auto h2 = test_scf::make_h2(); - auto aos = test_scf::h2_aos().ao_basis_set(); +TEMPLATE_LIST_TEST_CASE("SCFDriver", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto h2 = test_scf::make_h2(); + auto aos = test_scf::h2_aos().ao_basis_set(); + + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + pcorr->at() = -1.1167592336; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - const auto e = mm.run_as("SCF Driver", aos, h2); - REQUIRE_THAT(e, WithinAbs(-1.1167592336, 1E-6)); + const auto e = mm.template run_as("SCF Driver", aos, h2); + REQUIRE(approximately_equal(corr, e, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/driver/scf_loop.cpp b/tests/cxx/integration_tests/driver/scf_loop.cpp index 5c9cbfd..51870ce 100644 --- a/tests/cxx/integration_tests/driver/scf_loop.cpp +++ b/tests/cxx/integration_tests/driver/scf_loop.cpp @@ -24,17 +24,26 @@ using egy_pt = simde::eval_braket; template using pt = simde::Optimize, WFType>; -TEST_CASE("SCFLoop") { - using wf_type = simde::type::rscf_wf; - auto mm = test_scf::load_modules(); - auto& mod = mm.at("Loop"); +TEMPLATE_LIST_TEST_CASE("SCFLoop", "", test_scf::float_types) { + using float_type = TestType; + using wf_type = simde::type::rscf_wf; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Loop"); + + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; using index_set = typename wf_type::orbital_index_set_type; - wf_type psi0(index_set{0}, test_scf::h2_cmos()); + wf_type psi0(index_set{0}, test_scf::h2_cmos()); auto H = test_scf::h2_hamiltonian(); chemist::braket::BraKet H_00(psi0, H, psi0); - const auto& [e, psi] = mod.run_as>(H_00, psi0); - REQUIRE_THAT(e, WithinAbs(-1.1167592336, 1E-6)); + const auto& [e, psi] = mod.template run_as>(H_00, psi0); + + pcorr->at() = -1.1167592336; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/guess/core.cpp b/tests/cxx/integration_tests/guess/core.cpp index 8ebb6f2..7c12cf6 100644 --- a/tests/cxx/integration_tests/guess/core.cpp +++ b/tests/cxx/integration_tests/guess/core.cpp @@ -20,23 +20,25 @@ using rscf_wf = simde::type::rscf_wf; using pt = simde::InitialGuess; using simde::type::tensor; -TEST_CASE("Core") { - auto mm = test_scf::load_modules(); +TEMPLATE_LIST_TEST_CASE("Core", "", test_scf::float_types) { + using float_type = TestType; + + auto mm = test_scf::load_modules(); auto aos = test_scf::h2_aos(); auto H = test_scf::h2_hamiltonian(); auto mod = mm.at("Core guess"); - auto psi = mod.run_as(H, aos); + auto psi = mod.template run_as(H, aos); typename rscf_wf::orbital_index_set_type occs{0}; REQUIRE(psi.orbital_indices() == occs); REQUIRE(psi.orbitals().from_space() == aos); - const auto& evals = psi.orbitals().diagonalized_matrix(); - using allocator_type = tensorwrapper::allocator::Eigen; - const auto& eval_buffer = allocator_type::rebind(evals.buffer()); + const auto& evals = psi.orbitals().diagonalized_matrix(); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{2}; + auto pbuffer = alloc.construct({-1.25330893, -0.47506974}); + tensorwrapper::Tensor corr(shape_corr, std::move(pbuffer)); - const auto tol = 1E-6; - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(corr, evals, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/integration_tests.hpp b/tests/cxx/integration_tests/integration_tests.hpp index 36ed094..1f879b2 100644 --- a/tests/cxx/integration_tests/integration_tests.hpp +++ b/tests/cxx/integration_tests/integration_tests.hpp @@ -23,8 +23,12 @@ namespace test_scf { +/// Floating point types to test +using float_types = std::tuple; + /// Factors out setting submodules for SCF plugin from other plugins -inline auto load_modules() { +template +pluginplay::ModuleManager load_modules() { auto rv = std::make_shared(); pluginplay::ModuleManager mm(rv, nullptr); scf::load_modules(mm); @@ -44,6 +48,15 @@ inline auto load_modules() { mm.change_submod("Diagonalization Fock update", "Overlap matrix builder", "Overlap"); + if constexpr(!std::is_same_v) { + mm.change_input("Evaluate 2-Index BraKet", "With UQ?", true); + mm.change_input("Evaluate 4-Index BraKet", "With UQ?", true); + mm.change_input("Overlap", "With UQ?", true); + mm.change_input("ERI4", "With UQ?", true); + mm.change_input("Kinetic", "With UQ?", true); + mm.change_input("Nuclear", "With UQ?", true); + } + return mm; } diff --git a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp index 3655ed8..cc9a1cc 100644 --- a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp @@ -26,22 +26,31 @@ template using erased_type = chemist::braket::BraKet; -TEST_CASE("DeterminantDriver") { - auto mm = test_scf::load_modules(); - auto mod = mm.at("Determinant driver"); +TEMPLATE_LIST_TEST_CASE("DeterminantDriver", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto mod = mm.at("Determinant driver"); using wf_type = simde::type::rscf_wf; using index_set = typename wf_type::orbital_index_set_type; - wf_type psi(index_set{0}, test_scf::h2_cmos()); + wf_type psi(index_set{0}, test_scf::h2_cmos()); simde::type::many_electrons es{2}; + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + SECTION("Calling Kinetic") { simde::type::T_e_type T_e(es); chemist::braket::BraKet braket(psi, T_e, psi); erased_type copy_braket(braket); - const auto& T = mod.run_as>(copy_braket); - REQUIRE_THAT(T, WithinAbs(0.637692, 1E-6)); + const auto& T = mod.template run_as>(copy_braket); + + pcorr->at() = 0.637692; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, T, 1E-6)); } SECTION("Calling Electron-Nuclear Attraction") { @@ -49,25 +58,34 @@ TEST_CASE("DeterminantDriver") { simde::type::V_en_type V_en(es, h2_nuclei); chemist::braket::BraKet braket(psi, V_en, psi); erased_type copy_braket(braket); - const auto& V = mod.run_as>(copy_braket); - REQUIRE_THAT(V, WithinAbs(-1.96830777255516853, 1E-6)); + const auto& V = mod.template run_as>(copy_braket); + + pcorr->at() = -1.96830777255516853; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, V, 1E-6)); } SECTION("Calling J") { - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); simde::type::J_e_type J_e(es, rho); chemist::braket::BraKet braket(psi, J_e, psi); erased_type copy_braket(braket); - const auto& J = mod.run_as>(copy_braket); - REQUIRE_THAT(J, WithinAbs(0.76056339681664897, 1E-6)); + const auto& J = mod.template run_as>(copy_braket); + + pcorr->at() = 0.76056339681664897; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, J, 1E-6)); } SECTION("Calling K") { - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); simde::type::K_e_type K_e(es, rho); chemist::braket::BraKet braket(psi, K_e, psi); erased_type copy_braket(braket); - const auto& K = mod.run_as>(copy_braket); - REQUIRE_THAT(K, WithinAbs(0.76056339681664897, 1E-6)); + const auto& K = mod.template run_as>(copy_braket); + + pcorr->at() = 0.76056339681664897; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, K, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp index 28406f2..88b81ac 100644 --- a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp +++ b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp @@ -21,14 +21,20 @@ template using pt = simde::eval_braket; -TEST_CASE("ElectronicEnergy") { - auto mm = test_scf::load_modules(); - auto mod = mm.at("Electronic energy"); +TEMPLATE_LIST_TEST_CASE("ElectronicEnergy", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto mod = mm.at("Electronic energy"); using wf_type = simde::type::rscf_wf; using index_set = typename wf_type::orbital_index_set_type; - wf_type psi(index_set{0}, test_scf::h2_cmos()); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + + wf_type psi(index_set{0}, test_scf::h2_cmos()); simde::type::many_electrons es{2}; simde::type::T_e_type T_e(es); @@ -36,13 +42,16 @@ TEST_CASE("ElectronicEnergy") { auto h2_nuclei = test_scf::make_h2(); simde::type::V_en_type V_en(es, h2_nuclei); - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); simde::type::J_e_type J_e(es, rho); simde::type::K_e_type K_e(es, rho); simde::type::electronic_hamiltonian H_e(T_e * 2.0 + V_en * 2.0 + J_e * 2.0 - K_e); chemist::braket::BraKet braket(psi, H_e, psi); - const auto& E_elec = mod.run_as>(braket); - REQUIRE_THAT(E_elec, WithinAbs(-1.90066758625308307, 1E-6)); + const auto& E_elec = mod.template run_as>(braket); + + pcorr->at() = -1.90066758625308307; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, E_elec, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/fock.cpp b/tests/cxx/integration_tests/matrix_builder/fock.cpp index 1e2b968..b228c7d 100644 --- a/tests/cxx/integration_tests/matrix_builder/fock.cpp +++ b/tests/cxx/integration_tests/matrix_builder/fock.cpp @@ -21,11 +21,19 @@ using pt = simde::aos_f_e_aos; using simde::type::t_e_type; using simde::type::v_en_type; -TEST_CASE("Fock Matrix Builder") { - auto mm = test_scf::load_modules(); +TEMPLATE_LIST_TEST_CASE("Fock Matrix Builder", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Fock Matrix Builder"); auto aos = test_scf::h2_aos(); + using tensorwrapper::operations::approximately_equal; + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{2, 2}; + tensorwrapper::layout::Physical l(shape_corr); + auto pcorr = alloc.allocate(l); + SECTION("No J or K") { auto h2 = test_scf::make_h2(); simde::type::electron e; @@ -33,29 +41,31 @@ TEST_CASE("Fock Matrix Builder") { simde::type::fock f_e; f_e.emplace_back(1.0, std::make_unique(e)); f_e.emplace_back(1.0, std::make_unique(e, h2)); - const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); - - using alloc_type = tensorwrapper::allocator::Eigen; - const auto& F_eigen = alloc_type::rebind(F.buffer()); - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-1.120958, 1E-6)); - REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-1.120958, 1E-6)); + chemist::braket::BraKet f_mn(aos, f_e, aos); + const auto& F = mod.template run_as(f_mn); + + pcorr->at(0, 0) = -1.120958; + pcorr->at(0, 1) = -0.959374; + pcorr->at(1, 0) = -0.959374; + pcorr->at(1, 1) = -1.120958; + + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + + REQUIRE(approximately_equal(F, corr, 1E-6)); } SECTION("With J and K") { - auto f_e = test_scf::h2_fock(); + auto f_e = test_scf::h2_fock(); + chemist::braket::BraKet f_mn(aos, f_e, aos); + const auto& F = mod.template run_as(f_mn); - const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); + pcorr->at(0, 0) = -0.319459; + pcorr->at(0, 1) = -0.571781; + pcorr->at(1, 0) = -0.571781; + pcorr->at(1, 1) = -0.319459; - using alloc_type = tensorwrapper::allocator::Eigen; - const auto& F_eigen = alloc_type::rebind(F.buffer()); + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-0.319459, 1E-6)); - REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-0.319459, 1E-6)); + REQUIRE(approximately_equal(F, corr, 1E-6)); } -} \ No newline at end of file +} diff --git a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp index b2a61d5..d3a0b46 100644 --- a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp @@ -19,43 +19,28 @@ using pt = simde::aos_op_base_aos; using simde::type::tensor; -namespace { - -void compare_matrices(const tensor& A, const tensor& A_corr) { - using Catch::Matchers::WithinAbs; - using alloc_type = tensorwrapper::allocator::Eigen; - const auto& A_buffer = alloc_type::rebind(A.buffer()); - const auto& A_corr_buffer = alloc_type::rebind(A_corr.buffer()); - - const auto tol = 1E-6; - - REQUIRE_THAT(A_buffer.at(0, 0), WithinAbs(A_corr_buffer.at(0, 0), 1E-6)); - REQUIRE_THAT(A_buffer.at(0, 1), WithinAbs(A_corr_buffer.at(0, 1), 1E-6)); - REQUIRE_THAT(A_buffer.at(1, 0), WithinAbs(A_corr_buffer.at(1, 0), 1E-6)); - REQUIRE_THAT(A_buffer.at(1, 1), WithinAbs(A_corr_buffer.at(1, 1), 1E-6)); -} - -} // namespace - using erased_type = chemist::braket::BraKet; -TEST_CASE("SCFIntegralsDriver") { - auto mm = test_scf::load_modules(); - auto aos = test_scf::h2_aos(); - auto mod = mm.at("SCF integral driver"); +TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto aos = test_scf::h2_aos(); + auto mod = mm.at("SCF integral driver"); simde::type::electron e; - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); + + using tensorwrapper::operations::approximately_equal; SECTION("Calling Kinetic") { auto& tmod = mm.at("Kinetic"); simde::type::t_e_type t_e(e); chemist::braket::BraKet braket(aos, t_e, aos); erased_type copy_braket(braket); - const auto& T = mod.run_as(copy_braket); - const auto& T_corr = tmod.run_as(braket); - compare_matrices(T, T_corr); + const auto& T = mod.template run_as(copy_braket); + const auto& T_corr = tmod.template run_as(braket); + REQUIRE(approximately_equal(T, T_corr, 1E-6)); } SECTION("Calling Electron-Nuclear Attraction") { @@ -64,9 +49,9 @@ TEST_CASE("SCFIntegralsDriver") { simde::type::v_en_type v_en(e, h2_nuclei); chemist::braket::BraKet braket(aos, v_en, aos); erased_type copy_braket(braket); - const auto& V = mod.run_as(copy_braket); - const auto& V_corr = tmod.run_as(braket); - compare_matrices(V, V_corr); + const auto& V = mod.template run_as(copy_braket); + const auto& V_corr = tmod.template run_as(braket); + REQUIRE(approximately_equal(V, V_corr, 1E-6)); } SECTION("Calling J Matrix") { @@ -74,9 +59,9 @@ TEST_CASE("SCFIntegralsDriver") { simde::type::j_e_type j_e(e, rho); chemist::braket::BraKet braket(aos, j_e, aos); erased_type copy_braket(braket); - const auto& J = mod.run_as(copy_braket); - const auto& J_corr = jmod.run_as(braket); - compare_matrices(J, J_corr); + const auto& J = mod.template run_as(copy_braket); + const auto& J_corr = jmod.template run_as(braket); + REQUIRE(approximately_equal(J, J_corr, 1E-6)); } SECTION("Calling K Matrix") { @@ -84,22 +69,22 @@ TEST_CASE("SCFIntegralsDriver") { simde::type::k_e_type k_e(e, rho); chemist::braket::BraKet braket(aos, k_e, aos); erased_type copy_braket(braket); - const auto& K = mod.run_as(copy_braket); - const auto& K_corr = kmod.run_as(braket); - compare_matrices(K, K_corr); + const auto& K = mod.template run_as(copy_braket); + const auto& K_corr = kmod.template run_as(braket); + REQUIRE(approximately_equal(K, K_corr, 1E-6)); } SECTION("Calling density matrix") { auto& pmod = mm.at("Density matrix builder"); - auto cmos = test_scf::h2_cmos(); + auto cmos = test_scf::h2_cmos(); std::vector occs{1, 0}; simde::type::rho_e rho_hat(cmos, occs); chemist::braket::BraKet braket(aos, rho_hat, aos); erased_type copy_braket(braket); - const auto& P = mod.run_as(copy_braket); + const auto& P = mod.template run_as(copy_braket); using op_pt = simde::aos_rho_e_aos; - const auto& P_corr = pmod.run_as(braket); - compare_matrices(P, P_corr); + const auto& P_corr = pmod.template run_as(braket); + REQUIRE(approximately_equal(P, P_corr, 1E-6)); } // SECTION("Calling Fock Matrix") { diff --git a/tests/cxx/integration_tests/update/diagonalization.cpp b/tests/cxx/integration_tests/update/diagonalization.cpp index c174044..7a10d3f 100644 --- a/tests/cxx/integration_tests/update/diagonalization.cpp +++ b/tests/cxx/integration_tests/update/diagonalization.cpp @@ -23,9 +23,10 @@ using orbital_index_set = typename wf_t::orbital_index_set_type; using simde::type::t_e_type; using simde::type::v_en_type; -TEST_CASE("Diagaonalization") { - auto mm = test_scf::load_modules(); - auto& mod = mm.at("Diagonalization Fock update"); +TEMPLATE_LIST_TEST_CASE("Diagaonalization", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Diagonalization Fock update"); auto aos = test_scf::h2_aos(); auto h2 = test_scf::make_h2(); @@ -41,18 +42,18 @@ TEST_CASE("Diagaonalization") { simde::type::cmos cmos(empty, aos, empty); simde::type::rscf_wf core_guess(occs, cmos); - const auto& psi = mod.run_as(f_e, core_guess); + const auto& psi = mod.template run_as(f_e, core_guess); REQUIRE(psi.orbital_indices() == occs); REQUIRE(psi.orbitals().from_space() == aos); // Check orbital energies - const auto& evals = psi.orbitals().diagonalized_matrix(); - using vector_allocator = tensorwrapper::allocator::Eigen; - const auto& eval_buffer = vector_allocator::rebind(evals.buffer()); - - const auto tol = 1E-6; - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); + const auto& evals = psi.orbitals().diagonalized_matrix(); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + auto corr_buffer = alloc.construct({-1.25330893, -0.47506974}); + tensorwrapper::shape::Smooth corr_shape{2}; + tensorwrapper::Tensor corr(corr_shape, std::move(corr_buffer)); + + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(evals, corr, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index 62f8c31..08719cb 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -23,6 +23,8 @@ namespace test_scf { +using float_types = std::tuple; + /// Makes a H nucleus at the point @p x, @p y, @p z inline auto h_nucleus(double x, double y, double z) { return simde::type::nucleus("H", 1ul, 1836.15, x, y, z); @@ -91,29 +93,53 @@ inline auto h2_aos() { return simde::type::aos(h_basis(h2)); } +template inline auto h2_mos() { - using mos_type = simde::type::mos; - using tensor_type = typename mos_type::transform_type; - tensor_type c({{-0.565516, -1.07019}, {-0.565516, 1.07019}}); - return mos_type(h2_aos(), std::move(c)); + using mos_type = simde::type::mos; + using tensor_type = typename mos_type::transform_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{2, 2}; + tensorwrapper::layout::Physical l(shape); + auto c_buffer = alloc.allocate(l); + c_buffer->at(0, 0) = -0.565516; + c_buffer->at(0, 1) = -1.07019; + c_buffer->at(1, 0) = -0.565516; + c_buffer->at(1, 1) = 1.07019; + tensor_type t(shape, std::move(c_buffer)); + return mos_type(h2_aos(), std::move(t)); } +template inline auto h2_cmos() { - using cmos_type = simde::type::cmos; - using tensor_type = typename cmos_type::transform_type; - tensor_type e({-1.25330893, -0.47506974}); - return cmos_type(e, h2_aos(), h2_mos().transform()); + using cmos_type = simde::type::cmos; + using tensor_type = typename cmos_type::transform_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{2}; + tensorwrapper::layout::Physical l(shape); + auto e_buffer = alloc.allocate(l); + e_buffer->at(0) = -1.25330893; + e_buffer->at(1) = -0.47506974; + tensor_type e(shape, std::move(e_buffer)); + return cmos_type(std::move(e), h2_aos(), h2_mos().transform()); } +template inline auto h2_density() { - using density_type = simde::type::decomposable_e_density; - typename density_type::value_type rho( - {{0.31980835, 0.31980835}, {0.31980835, 0.31980835}}); - return density_type(rho, h2_mos()); + using density_type = simde::type::decomposable_e_density; + using tensor_type = typename density_type::value_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{2, 2}; + tensorwrapper::layout::Physical l(shape); + auto pbuffer = alloc.construct(l, 0.31980835); + tensor_type t(shape, std::move(pbuffer)); + return density_type(std::move(t), h2_mos()); } /// The Fock matrix consistent with h2_hamiltonian and h2_density -template +template inline auto h2_fock() { ElectronType es; if constexpr(std::is_same_v) { @@ -121,7 +147,7 @@ inline auto h2_fock() { } auto h2 = make_h2(); - auto rho = h2_density(); + auto rho = h2_density(); simde::type::fock F; using namespace chemist::qm_operator; using t_type = Kinetic; diff --git a/tests/cxx/unit_tests/coulombs_law.cpp b/tests/cxx/unit_tests/coulombs_law.cpp index 80038bf..9300d26 100644 --- a/tests/cxx/unit_tests/coulombs_law.cpp +++ b/tests/cxx/unit_tests/coulombs_law.cpp @@ -23,10 +23,16 @@ using pt = simde::charge_charge_interaction; using Catch::Matchers::WithinAbs; TEST_CASE("CoulombsLaw") { + using float_type = double; pluginplay::ModuleManager mm; scf::load_modules(mm); auto& mod = mm.at("Coulomb's Law"); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + using charge_type = typename simde::type::charges::value_type; charge_type q0(1.0, 2.0, 3.0, 4.0); charge_type q1(-1.0, 3.0, 4.0, 5.0); @@ -36,24 +42,32 @@ TEST_CASE("CoulombsLaw") { simde::type::charges qs{q0, q1, q2}; SECTION("empty points") { - auto e = mod.run_as(empty, empty); - REQUIRE(e == 0.0); + auto e = mod.run_as(empty, empty); + pcorr->at() = 0.0; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } SECTION("charges w/ itself") { - auto e = mod.run_as(qs, qs); - REQUIRE_THAT(e, WithinAbs(-1.0103629710818451, 1E-6)); + auto e = mod.run_as(qs, qs); + pcorr->at() = -1.0103629710818451; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } SECTION("charges w/ empty") { - auto e = mod.run_as(qs, empty); - REQUIRE_THAT(e, WithinAbs(0.0, 1E-6)); + auto e = mod.run_as(qs, empty); + pcorr->at() = 0.0; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } SECTION("charges w/ different charges") { simde::type::charges qs0{q0}; simde::type::charges qs12{q1, q2}; - auto e = mod.run_as(qs0, qs12); - REQUIRE_THAT(e, WithinAbs(-0.1443375672974065, 1E-6)); + auto e = mod.run_as(qs0, qs12); + pcorr->at() = -0.1443375672974065; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp index 577cdff..fcbcbad 100644 --- a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp +++ b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp @@ -18,7 +18,8 @@ #include #include -TEST_CASE("EigenGeneralized") { +TEMPLATE_LIST_TEST_CASE("EigenGeneralized", "", test_scf::float_types) { + using float_type = TestType; pluginplay::ModuleManager mm; scf::load_modules(mm); @@ -26,15 +27,31 @@ TEST_CASE("EigenGeneralized") { auto& mod = mm.at("Generalized eigensolve via Eigen"); - simde::type::tensor A({{1.0, 2.0}, {2.0, 3.0}}); - simde::type::tensor B({{1.0, 0.0}, {0.0, 1.0}}); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape{2, 2}; + tensorwrapper::layout::Physical l(shape); + auto A_buffer = alloc.allocate(l); + A_buffer->at(0, 0) = 1.0; + A_buffer->at(0, 1) = 2.0; + A_buffer->at(1, 0) = 2.0; + A_buffer->at(1, 1) = 3.0; + + auto B_buffer = alloc.allocate(l); + B_buffer->at(0, 0) = 1.0; + B_buffer->at(0, 1) = 0.0; + B_buffer->at(1, 0) = 0.0; + B_buffer->at(1, 1) = 1.0; + + simde::type::tensor A(shape, std::move(A_buffer)); + simde::type::tensor B(shape, std::move(B_buffer)); auto&& [values, vector] = mod.run_as(A, B); - using value_alloc_t = tensorwrapper::allocator::Eigen; - const auto& eigen_values = value_alloc_t::rebind(values.buffer()); + const auto& eigen_values = alloc.rebind(values.buffer()); + + auto corr_buffer = alloc.construct({-0.236068, 4.236068}); + tensorwrapper::shape::Smooth corr_shape{2}; + simde::type::tensor corr(corr_shape, std::move(corr_buffer)); - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eigen_values.at(0), WithinAbs(-0.236068, 1E-6)); - REQUIRE_THAT(eigen_values.at(1), WithinAbs(4.236068, 1E-6)); + REQUIRE(tensorwrapper::operations::approximately_equal(corr, values, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/unit_tests/fock_operator/restricted.cpp b/tests/cxx/unit_tests/fock_operator/restricted.cpp index 26f7a17..45cdcf1 100644 --- a/tests/cxx/unit_tests/fock_operator/restricted.cpp +++ b/tests/cxx/unit_tests/fock_operator/restricted.cpp @@ -18,22 +18,22 @@ #include #include -TEST_CASE("Restricted") { +TEMPLATE_LIST_TEST_CASE("Restricted", "", test_scf::float_types) { pluginplay::ModuleManager mm; scf::load_modules(mm); - + using float_type = double; using density_type = simde::type::decomposable_e_density; using pt = simde::FockOperator; auto H = test_scf::h2_hamiltonian(); auto h2 = test_scf::make_h2(); - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); SECTION("Many-electron version") { auto& mod = mm.at("Restricted Fock Op"); SECTION("No density") { - density_type rho; - auto F = mod.run_as(H, rho); + density_type rho_empty; + auto F = mod.run_as(H, rho_empty); simde::type::many_electrons es(2); simde::type::fock F_corr; @@ -45,8 +45,6 @@ TEST_CASE("Restricted") { } SECTION("Density") { - auto rho = test_scf::h2_density(); - auto F = mod.run_as(H, rho); auto F_corr = test_scf::h2_fock(); REQUIRE(F == F_corr); @@ -56,8 +54,8 @@ TEST_CASE("Restricted") { auto& mod = mm.at("Restricted One-Electron Fock Op"); SECTION(" No density") { - density_type rho; - auto f = mod.run_as(H, rho); + density_type rho_empty; + auto f = mod.run_as(H, rho_empty); simde::type::fock f_corr; simde::type::electron e; using simde::type::t_e_type; @@ -68,8 +66,6 @@ TEST_CASE("Restricted") { } SECTION("Density") { - auto rho = test_scf::h2_density(); - auto f = mod.run_as(H, rho); auto f_corr = test_scf::h2_fock(); REQUIRE(f == f_corr); diff --git a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp index b050194..1b5b558 100644 --- a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp +++ b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp @@ -18,23 +18,24 @@ using pt = simde::aos_rho_e_aos; -TEST_CASE("Density Matrix Builder") { +TEMPLATE_LIST_TEST_CASE("Density Matrix Builder", "", test_scf::float_types) { + using float_type = TestType; pluginplay::ModuleManager mm; scf::load_modules(mm); auto& mod = mm.at("Density matrix builder"); auto aos = test_scf::h2_aos(); - auto cmos = test_scf::h2_cmos(); + auto cmos = test_scf::h2_cmos(); std::vector occs{1, 0}; simde::type::rho_e rho_hat(cmos, occs); chemist::braket::BraKet p_mn(aos, rho_hat, aos); - const auto& P = mod.run_as(p_mn); - using allocator_type = tensorwrapper::allocator::Eigen; - const auto& P_eigen = allocator_type::rebind(P.buffer()); + const auto& P = mod.run_as(p_mn); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth corr_shape{2, 2}; + tensorwrapper::layout::Physical l(corr_shape); + auto corr_buffer = alloc.construct(l, 0.31980835); + tensorwrapper::Tensor corr(corr_shape, std::move(corr_buffer)); - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(P_eigen.at(0, 0), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen.at(0, 1), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen.at(1, 0), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen.at(1, 1), WithinAbs(0.31980835, 1E-6)); + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(P, corr, 1E-6)); } \ No newline at end of file