diff --git a/CMakeLists.txt b/CMakeLists.txt index 4909c16..1421c5b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,11 +49,19 @@ foreach(dependency_i ${DEPENDENCIES}) include(get_${dependency_i}) endforeach() +cmaize_find_or_build_dependency( + eigen + URL https://www.gitlab.com/libeigen/eigen + VERSION 3.4.0 + BUILD_TARGET eigen + FIND_TARGET Eigen3::Eigen +) + cmaize_add_library( scf SOURCE_DIR "${CMAKE_CURRENT_LIST_DIR}/src/scf" INCLUDE_DIRS "${CMAKE_CURRENT_LIST_DIR}/include/scf" - DEPENDS "${DEPENDENCIES}" + DEPENDS "${DEPENDENCIES}" eigen ) if("${BUILD_TAMM_SCF}") diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp new file mode 100644 index 0000000..ec47bed --- /dev/null +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -0,0 +1,88 @@ +/* + * Copyright 2024 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "eigen_solver.hpp" +#include +#include + +namespace scf::eigen_solver { + +using pt = simde::GeneralizedEigenSolve; + +const auto desc = R"( +Generalized Eigen Solve via Eigen +--------------------------------- + +TODO: Write me!!! +)"; + +MODULE_CTOR(EigenGeneralized) { + description(desc); + satisfies_property_type(); +} + +MODULE_RUN(EigenGeneralized) { + auto&& [A, B] = pt::unwrap_inputs(inputs); + + // Convert to Eigen buffers + using matrix_alloc_t = tensorwrapper::allocator::Eigen; + using vector_alloc_t = tensorwrapper::allocator::Eigen; + const auto& eigen_A = matrix_alloc_t::rebind(A.buffer()); + const auto& eigen_B = matrix_alloc_t::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& 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::GeneralizedEigenSolver 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); + + using matrix_buffer = typename matrix_alloc_t::eigen_buffer_type; + using vector_buffer = typename vector_alloc_t::eigen_buffer_type; + + 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).real(); + for(auto j = 0; j < cols; ++j) { + matrix_tensor(i, j) = eigen_vectors(i, j).real(); + } + } + + 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); + + auto rv = results(); + return pt::wrap_results(rv, values, vectors); +} + +} // namespace scf::eigen_solver \ No newline at end of file diff --git a/src/scf/eigen_solver/eigen_solver.hpp b/src/scf/eigen_solver/eigen_solver.hpp new file mode 100644 index 0000000..41d99db --- /dev/null +++ b/src/scf/eigen_solver/eigen_solver.hpp @@ -0,0 +1,28 @@ +/* + * Copyright 2024 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include + +namespace scf::eigen_solver { + +DECLARE_MODULE(EigenGeneralized); + +inline void load_modules(pluginplay::ModuleManager& mm) { + mm.add_module("Generalized eigensolve via Eigen"); +} + +} // namespace scf::eigen_solver \ No newline at end of file diff --git a/src/scf/scf_mm.cpp b/src/scf/scf_mm.cpp index 2dab428..1cc4260 100644 --- a/src/scf/scf_mm.cpp +++ b/src/scf/scf_mm.cpp @@ -14,6 +14,7 @@ * limitations under the License. */ +#include "eigen_solver/eigen_solver.hpp" #include "fock_operator/fock_operator.hpp" #include "scf_modules.hpp" #include @@ -21,6 +22,7 @@ namespace scf { void load_modules(pluginplay::ModuleManager& mm) { + eigen_solver::load_modules(mm); fock_operator::load_modules(mm); mm.add_module("Coulomb's Law"); #ifdef BUILD_TAMM_SCF diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index 45d0e21..40f6391 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include namespace test_scf { diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp new file mode 100644 index 0000000..32e1ad7 --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp @@ -0,0 +1,40 @@ +/* + * Copyright 2024 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "../../test_scf.hpp" +#include +#include + +TEST_CASE("EigenGeneralized") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + + using pt = simde::GeneralizedEigenSolve; + + 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}}); + + 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()); + + 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)); +} \ No newline at end of file