Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 11 additions & 17 deletions src/scf/eigen_solver/eigen_generalized.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -38,14 +38,13 @@ MODULE_RUN(EigenGeneralized) {
auto&& [A, B] = pt::unwrap_inputs(inputs);

// Convert to Eigen buffers
using matrix_alloc_t = tensorwrapper::allocator::Eigen<double, 2>;
using vector_alloc_t = tensorwrapper::allocator::Eigen<double, 1>;
const auto& eigen_A = matrix_alloc_t::rebind(A.buffer());
const auto& eigen_B = matrix_alloc_t::rebind(B.buffer());
tensorwrapper::allocator::Eigen<double> 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.value().data();
const auto* pB = eigen_B.value().data();
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);
Expand All @@ -63,23 +62,18 @@ MODULE_RUN(EigenGeneralized) {
tensorwrapper::layout::Physical vector_layout(vector_shape);
tensorwrapper::layout::Physical matrix_layout(matrix_shape);

using matrix_buffer = typename matrix_alloc_t::eigen_buffer_type;
using vector_buffer = typename vector_alloc_t::eigen_buffer_type;
auto pvalues_buffer = allocator.allocate(vector_layout);
auto pvectors_buffer = allocator.allocate(matrix_layout);

typename vector_buffer::data_type vector_tensor(rows);
typename matrix_buffer::data_type matrix_tensor(rows, cols);
for(auto i = 0; i < rows; ++i) {
vector_tensor(i) = eigen_values(i);
pvalues_buffer->at(i) = eigen_values(i);
for(auto j = 0; j < cols; ++j) {
matrix_tensor(i, j) = eigen_vectors(i, j);
pvectors_buffer->at(i, j) = eigen_vectors(i, j);
}
}

vector_buffer values_buffer(vector_tensor, vector_layout);
matrix_buffer vectors_buffer(matrix_tensor, matrix_layout);

simde::type::tensor values(vector_shape, values_buffer);
simde::type::tensor vectors(matrix_shape, vectors_buffer);
simde::type::tensor values(vector_shape, std::move(pvalues_buffer));
simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer));

auto rv = results();
return pt::wrap_results(rv, values, vectors);
Expand Down
32 changes: 20 additions & 12 deletions src/scf/matrix_builder/density_matrix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
*/

#include "matrix_builder.hpp"
#include <unsupported/Eigen/CXX11/Tensor>

namespace scf::matrix_builder {
namespace {
Expand Down Expand Up @@ -62,27 +63,34 @@ MODULE_RUN(DensityMatrix) {
"Please shuffle your orbitals so that the ensemble is a "
"contiguous slice of the orbitals.");

// Step 2: Grab the orbitals in the ensemble
using allocator_type = tensorwrapper::allocator::Eigen<double, 2>;
using buffer_type = typename allocator_type::eigen_buffer_type;
using eigen_tensor_type = typename buffer_type::data_type;
using allocator_type = tensorwrapper::allocator::Eigen<double>;
using tensor_type = Eigen::Tensor<double, 2, Eigen::RowMajor>;
using const_tensor_type = Eigen::Tensor<const double, 2, Eigen::RowMajor>;
using map_type = Eigen::TensorMap<tensor_type>;
using const_map_type = Eigen::TensorMap<const_tensor_type>;

const auto& c_buffer = allocator_type::rebind(c.buffer());
const auto& c_eigen = c_buffer.value();
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<Eigen::Index, 2>;
slice_t offsets{0, 0};
slice_t extents{n_aos, Eigen::Index(participants.size())};
eigen_tensor_type slice = c_eigen.slice(offsets, extents);
auto slice = c_eigen.slice(offsets, extents);

// Step 3: CC_dagger
using index_pair_t = Eigen::IndexPair<int>;
Eigen::array<index_pair_t, 1> modes{index_pair_t(1, 1)};
eigen_tensor_type p_eigen = slice.contract(slice, modes);
tensorwrapper::shape::Smooth p_shape(n_aos, n_aos);
tensorwrapper::layout::Physical l(p_shape);
buffer_type p_buffer(p_eigen, l);
simde::type::tensor p(p_shape, p_buffer);
p_eigen = slice.contract(slice, modes);

simde::type::tensor p(p_shape, std::move(pp_buffer));

auto rv = results();
return pt::wrap_results(rv, p);
Expand Down
15 changes: 6 additions & 9 deletions src/scf/matrix_builder/determinant_driver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,7 @@ MODULE_CTOR(DeterminantDriver) {
}

MODULE_RUN(DeterminantDriver) {
using float_type = double;
using wf_type = simde::type::rscf_wf;
const auto&& [braket] = pt<wf_type>::unwrap_inputs(inputs);
const auto& bra = braket.bra();
Expand All @@ -126,18 +127,14 @@ MODULE_RUN(DeterminantDriver) {
const auto& t = visitor.m_pt;

// Step 3: Contract
using allocator_type = tensorwrapper::allocator::Eigen<double, 2>;
using scalar_type = tensorwrapper::buffer::Eigen<double, 0>::data_type;
tensorwrapper::Tensor x;
x("") = rho("i,j") * t("i,j");

const auto& t_eigen = allocator_type::rebind(t.buffer()).value();
const auto& p_eigen = allocator_type::rebind(rho.buffer()).value();

using index_pair_t = Eigen::IndexPair<int>;
Eigen::array<index_pair_t, 2> modes{index_pair_t(0, 0), index_pair_t(1, 1)};
scalar_type x = p_eigen.contract(t_eigen, modes);
using allocator_type = tensorwrapper::allocator::Eigen<float_type>;
auto& x_buffer = allocator_type::rebind(x.buffer());

auto rv = results();
return pt<wf_type>::wrap_results(rv, x());
return pt<wf_type>::wrap_results(rv, x_buffer.at());
}

} // namespace scf::matrix_builder
8 changes: 3 additions & 5 deletions src/scf/matrix_builder/fock.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,18 +49,16 @@ MODULE_RUN(Fock) {
chemist::braket::BraKet term0(bra_aos, op_0, ket_aos);
F = ao_dispatcher.run_as<ao_pt>(term0);

using allocator_type = tensorwrapper::allocator::Eigen<double, 2>;
auto& f_buffer = allocator_type::rebind(F.buffer());
f_buffer.value() = f_buffer.value() * c0;
F("i,j") = F("i,j") * c0;

for(decltype(n_terms) i = 1; i < n_terms; ++i) {
const auto ci = f.coefficient(i);
const auto& op_i = f.get_operator(i);
chemist::braket::BraKet termi(bra_aos, op_i, ket_aos);
auto matrix = ao_dispatcher.run_as<ao_pt>(termi);

auto& matrix_buffer = allocator_type::rebind(matrix.buffer());
f_buffer.value() += ci * matrix_buffer.value();
matrix("i,j") = matrix("i,j") * ci;
F("i,j") = F("i,j") + matrix("i,j");
}
}

Expand Down
6 changes: 3 additions & 3 deletions tests/cxx/integration_tests/guess/core.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,11 @@ TEST_CASE("Core") {
REQUIRE(psi.orbital_indices() == occs);
REQUIRE(psi.orbitals().from_space() == aos);
const auto& evals = psi.orbitals().diagonalized_matrix();
using allocator_type = tensorwrapper::allocator::Eigen<double, 1>;
using allocator_type = tensorwrapper::allocator::Eigen<double>;
const auto& eval_buffer = allocator_type::rebind(evals.buffer());

const auto tol = 1E-6;
using Catch::Matchers::WithinAbs;
REQUIRE_THAT(eval_buffer.value()(0), WithinAbs(-1.25330893, tol));
REQUIRE_THAT(eval_buffer.value()(1), WithinAbs(-0.47506974, tol));
REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol));
REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol));
}
20 changes: 10 additions & 10 deletions tests/cxx/integration_tests/matrix_builder/fock.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,27 +35,27 @@ TEST_CASE("Fock Matrix Builder") {
f_e.emplace_back(1.0, std::make_unique<v_en_type>(e, h2));
const auto& F = mod.run_as<pt>(chemist::braket::BraKet(aos, f_e, aos));

using alloc_type = tensorwrapper::allocator::Eigen<double, 2>;
using alloc_type = tensorwrapper::allocator::Eigen<double>;
const auto& F_eigen = alloc_type::rebind(F.buffer());
using Catch::Matchers::WithinAbs;
REQUIRE_THAT(F_eigen.value()(0, 0), WithinAbs(-1.120958, 1E-6));
REQUIRE_THAT(F_eigen.value()(0, 1), WithinAbs(-0.959374, 1E-6));
REQUIRE_THAT(F_eigen.value()(1, 0), WithinAbs(-0.959374, 1E-6));
REQUIRE_THAT(F_eigen.value()(1, 1), WithinAbs(-1.120958, 1E-6));
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));
}

SECTION("With J and K") {
auto f_e = test_scf::h2_fock<simde::type::electron>();

const auto& F = mod.run_as<pt>(chemist::braket::BraKet(aos, f_e, aos));

using alloc_type = tensorwrapper::allocator::Eigen<double, 2>;
using alloc_type = tensorwrapper::allocator::Eigen<double>;
const auto& F_eigen = alloc_type::rebind(F.buffer());

using Catch::Matchers::WithinAbs;
REQUIRE_THAT(F_eigen.value()(0, 0), WithinAbs(-0.319459, 1E-6));
REQUIRE_THAT(F_eigen.value()(0, 1), WithinAbs(-0.571781, 1E-6));
REQUIRE_THAT(F_eigen.value()(1, 0), WithinAbs(-0.571781, 1E-6));
REQUIRE_THAT(F_eigen.value()(1, 1), WithinAbs(-0.319459, 1E-6));
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));
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -23,18 +23,16 @@ namespace {

void compare_matrices(const tensor& A, const tensor& A_corr) {
using Catch::Matchers::WithinAbs;
using alloc_type = tensorwrapper::allocator::Eigen<double, 2>;
using alloc_type = tensorwrapper::allocator::Eigen<double>;
const auto& A_buffer = alloc_type::rebind(A.buffer());
const auto& A_corr_buffer = alloc_type::rebind(A_corr.buffer());
const auto& A_eigen = A_buffer.value();
const auto& A_corr_eigen = A_corr_buffer.value();

const auto tol = 1E-6;

REQUIRE_THAT(A_eigen(0, 0), WithinAbs(A_corr_eigen(0, 0), 1E-6));
REQUIRE_THAT(A_eigen(0, 1), WithinAbs(A_corr_eigen(0, 1), 1E-6));
REQUIRE_THAT(A_eigen(1, 0), WithinAbs(A_corr_eigen(1, 0), 1E-6));
REQUIRE_THAT(A_eigen(1, 1), WithinAbs(A_corr_eigen(1, 1), 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
Expand Down
6 changes: 3 additions & 3 deletions tests/cxx/integration_tests/update/diagonalization.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,12 @@ TEST_CASE("Diagaonalization") {

// Check orbital energies
const auto& evals = psi.orbitals().diagonalized_matrix();
using vector_allocator = tensorwrapper::allocator::Eigen<double, 1>;
using vector_allocator = tensorwrapper::allocator::Eigen<double>;
const auto& eval_buffer = vector_allocator::rebind(evals.buffer());

const auto tol = 1E-6;
using Catch::Matchers::WithinAbs;
REQUIRE_THAT(eval_buffer.value()(0), WithinAbs(-1.25330893, tol));
REQUIRE_THAT(eval_buffer.value()(1), WithinAbs(-0.47506974, tol));
REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol));
REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol));
}
}
6 changes: 3 additions & 3 deletions tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,10 @@ TEST_CASE("EigenGeneralized") {

auto&& [values, vector] = mod.run_as<pt>(A, B);

using value_alloc_t = tensorwrapper::allocator::Eigen<double, 1>;
using value_alloc_t = tensorwrapper::allocator::Eigen<double>;
const auto& eigen_values = value_alloc_t::rebind(values.buffer());

using Catch::Matchers::WithinAbs;
REQUIRE_THAT(eigen_values.value()(0), WithinAbs(-0.236068, 1E-6));
REQUIRE_THAT(eigen_values.value()(1), WithinAbs(4.236068, 1E-6));
REQUIRE_THAT(eigen_values.at(0), WithinAbs(-0.236068, 1E-6));
REQUIRE_THAT(eigen_values.at(1), WithinAbs(4.236068, 1E-6));
}
12 changes: 6 additions & 6 deletions tests/cxx/unit_tests/matrix_builder/density_matrix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,12 +29,12 @@ TEST_CASE("Density Matrix Builder") {

chemist::braket::BraKet p_mn(aos, rho_hat, aos);
const auto& P = mod.run_as<pt>(p_mn);
using allocator_type = tensorwrapper::allocator::Eigen<double, 2>;
const auto& P_eigen = allocator_type::rebind(P.buffer()).value();
using allocator_type = tensorwrapper::allocator::Eigen<double>;
const auto& P_eigen = allocator_type::rebind(P.buffer());

using Catch::Matchers::WithinAbs;
REQUIRE_THAT(P_eigen(0, 0), WithinAbs(0.31980835, 1E-6));
REQUIRE_THAT(P_eigen(0, 1), WithinAbs(0.31980835, 1E-6));
REQUIRE_THAT(P_eigen(1, 0), WithinAbs(0.31980835, 1E-6));
REQUIRE_THAT(P_eigen(1, 1), WithinAbs(0.31980835, 1E-6));
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));
}