diff --git a/src/scf/driver/scf_loop.cpp b/src/scf/driver/scf_loop.cpp index 3b850f8..d455e65 100644 --- a/src/scf/driver/scf_loop.cpp +++ b/src/scf/driver/scf_loop.cpp @@ -28,23 +28,14 @@ struct Kernel { template auto operator()(const std::span& a) { - using tensorwrapper::types::fabs; - using clean_type = std::decay_t; - bool in_noise = false; - bool converged = false; - if constexpr(tensorwrapper::types::is_interval_v) { - auto abs_val = fabs(a[0]); - // If uncertainty is greater than the tolerance, return true - in_noise = abs_val.radius() > m_tol; - converged = abs_val.upper() < m_tol; - } else if constexpr(tensorwrapper::types::is_uncertain_v) { - auto abs_val = fabs(a[0]); - in_noise = abs_val.sd() > m_tol; - converged = (abs_val.mean() + abs_val.sd()) < m_tol; - } else { - converged = fabs(a[0]) < m_tol; - } - return std::make_pair(converged, in_noise); + using tensorwrapper::types::uq_center; + /// For UQ, the center/median of the difference is our best estimate + /// for the uncertainty caused by incomplete convergence. It plus/minus + // the radius/standard deviation is the total range of uncertainty. + // Convergence happens when the center/median is within the tolerance, + // regardless of the radius/standard deviation. + auto abs_val = std::fabs(uq_center(a[0])); + return abs_val < m_tol; } }; @@ -290,11 +281,7 @@ MODULE_RUN(SCFLoop) { auto e_conv = check_tolerance(de.buffer(), e_tol); auto g_conv = check_tolerance(grad_norm.buffer(), g_tol); auto dp_conv = check_tolerance(dp_norm.buffer(), dp_tol); - if(e_conv.first && g_conv.first && dp_conv.first) converged = true; - if(e_conv.second) { - logger.log(" Energy has converged to within uncertainty"); - converged = true; - } + if(e_conv && g_conv && dp_conv) converged = true; // If using DIIS and not converged, extrapolate new Fock matrix if(diis_on && !converged) { F = diis.extrapolate(F, grad); } diff --git a/src/scf/eigen_solver/ball_generalized.cpp b/src/scf/eigen_solver/ball_generalized.cpp deleted file mode 100644 index c95b1ca..0000000 --- a/src/scf/eigen_solver/ball_generalized.cpp +++ /dev/null @@ -1,282 +0,0 @@ -/* - * Copyright 2026 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. - */ - -/* - * Copyright 2026 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 "submodule_request.hpp" -#include "wtf/fp/float_view.hpp" - -#include -#include -#ifdef ENABLE_SIGMA -#include - -using pt = simde::GeneralizedEigenSolve; - -namespace scf::eigen_solver { -namespace { -const auto desc = R"( -Generalized Eigen Solve via Ball Arithmetic -------------------------------------------- -https://www.texmacs.org/joris/ball/ball.html -TODO: Write me!!! -H = TMT - L -n = number of rows (or columns) -k = max[max(1/(L_jj -Lii)), max(1/L_ii)] -Omega_n = n by n matrix 0+/-1 -eta = 6 sqrt(n) k ||H|| -T_ball = T(1 + eta Omega_n) -L_ball = B(1, eta) * L -For the generalized problem, AC=BCL, we'll try: -H = TAC - TBCL -where T = C^T -)"; - -template -auto wrap_subdiagonalization(const tensorwrapper::buffer::Contiguous& A_buf, - const tensorwrapper::buffer::Contiguous& B_buf, - pluginplay::SubmoduleRequest& eigen_solver_mod) { - using namespace tensorwrapper::buffer; - using namespace tensorwrapper::shape; - - auto A_data = get_raw_data(A_buf); - auto B_data = get_raw_data(B_buf); - - auto A_shape = A_buf.shape().make_smooth(); - auto B_shape = B_buf.shape().make_smooth(); - assert(A_shape.rank() == 2); - assert(B_shape.rank() == 2); - auto n_rows = A_shape.extent(0); - auto n_cols = A_shape.extent(1); - assert(n_rows == B_shape.extent(0)); - assert(n_cols == B_shape.extent(1)); - Smooth shape({n_rows, n_cols}); - - std::vector A_median(shape.size()); - std::vector B_median(shape.size()); - - for(std::size_t i = 0; i < shape.size(); ++i) { - A_median[i] = A_data[i].median(); - B_median[i] = B_data[i].median(); - } - - Contiguous A_median_buffer(std::move(A_median), shape); - Contiguous B_median_buffer(std::move(B_median), shape); - tensorwrapper::Tensor A_tensor(shape, std::move(A_median_buffer)); - tensorwrapper::Tensor B_tensor(shape, std::move(B_median_buffer)); - - return eigen_solver_mod.run_as(A_tensor, B_tensor); -}; - -template -auto convert_to_uq(const tensorwrapper::Tensor& values, - const tensorwrapper::Tensor& vectors) { - using float_type = typename UQType::value_t; - auto values_buf = make_contiguous(values.buffer()); - auto vectors_buf = make_contiguous(vectors.buffer()); - auto n = values_buf.shape().extent(0); - tensorwrapper::shape::Smooth vector_shape{n}; - tensorwrapper::shape::Smooth matrix_shape{n, n}; - using tensorwrapper::buffer::make_contiguous; - auto uq_values_buffer = make_contiguous(vector_shape); - auto uq_vectors_buffer = make_contiguous(matrix_shape); - using wtf::fp::float_cast; - for(std::size_t i = 0; i < n; ++i) { - auto vi = float_cast(values_buf.get_elem({i})); - uq_values_buffer.set_elem({i}, UQType(vi)); - for(std::size_t j = 0; j < n; ++j) { - auto vij = float_cast(vectors_buf.get_elem({i, j})); - uq_vectors_buffer.set_elem({i, j}, UQType(vij)); - } - } - tensorwrapper::Tensor uq_values(vector_shape, std::move(uq_values_buffer)); - tensorwrapper::Tensor uq_vectors(matrix_shape, - std::move(uq_vectors_buffer)); - return std::make_pair(uq_values, uq_vectors); -} - -template -auto compute_residual(const tensorwrapper::Tensor& A, - const tensorwrapper::Tensor& B, - const tensorwrapper::Tensor& C, - const tensorwrapper::Tensor& L) { - using label_type = tensorwrapper::Tensor::label_type; - using namespace tensorwrapper::buffer; - label_type ij("i,j"); - label_type ji("j,i"); - label_type jk("j,k"); - label_type ik("i,k"); - label_type j("j"); - - tensorwrapper::Tensor TA, TAC; - TA(ik) = C(ji) * A(jk); - TAC(ik) = TA(ij) * C(jk); - - tensorwrapper::Tensor TB, TBC; - TB(ik) = C(ji) * B(jk); - TBC(ik) = TB(ij) * C(jk); - - // TODO: Replace when batch contraction works... - auto TBC_buf = make_contiguous(TBC.buffer()); - auto L_buf = make_contiguous(L.buffer()); - auto TBCL_shape = TBC_buf.shape().make_smooth(); - auto TBCL_buffer = make_contiguous(TBCL_shape); - - auto n_rows = TBCL_shape.extent(0); - auto n_cols = TBCL_shape.extent(1); - using wtf::fp::float_cast; - for(std::size_t i = 0; i < n_rows; ++i) { - auto Li = float_cast(L_buf.get_elem({i})); - for(std::size_t j = 0; j < n_cols; ++j) { - auto TBCij = float_cast(TBC_buf.get_elem({i, j})); - TBCL_buffer.set_elem({i, j}, TBCij * Li); - } - } - tensorwrapper::Tensor TBCL(TBCL_shape, std::move(TBCL_buffer)); - - tensorwrapper::Tensor H; - H(ij) = TAC(ij) - TBCL(ij); - return H; -} - -template -auto compute_eta(const tensorwrapper::Tensor& H) { - auto H_buffer = make_contiguous(H.buffer()); - auto shape = H_buffer.shape(); - auto n_rows = shape.extent(0); - - using float_type = typename UQType::value_t; - using wtf::fp::float_cast; - float_type norm = 0.0; - for(std::size_t i = 0; i < n_rows; ++i) { - float_type col_norm = 0.0; - for(std::size_t j = 0; j < n_rows; ++j) { - auto Hji = H_buffer.get_elem({j, i}); - auto Hji_uq = float_cast(Hji); - auto Hji_abs = - std::max(std::abs(Hji_uq.lower()), std::abs(Hji_uq.upper())); - col_norm += Hji_abs * Hji_abs; - } - col_norm = std::sqrt(col_norm / n_rows); - norm = std::max(norm, col_norm); - } - return 6.0 * std::sqrt(n_rows) * norm; -} - -template -auto t_ball(const tensorwrapper::Tensor& C, - const typename UQType::value_t& eta) { - using namespace tensorwrapper::buffer; - UQType pm1(-1.0, 1.0); - auto shape = make_contiguous(C.buffer()).shape().make_smooth(); - std::vector eta_omega(shape.size(), eta * pm1); - Contiguous eta_omega_buffer(std::move(eta_omega), shape); - tensorwrapper::Tensor eta_omega_tensor(shape, std::move(eta_omega_buffer)); - tensorwrapper::Tensor T_ball, COmega; - COmega("i,k") = C("i,j") * eta_omega_tensor("j,k"); - T_ball("i,j") = C("i,j") + COmega("i,j"); - return T_ball; -} - -template -auto l_ball(const tensorwrapper::Tensor& L, - const typename UQType::value_t& eta) { - using namespace tensorwrapper::buffer; - UQType one_pm_eta(1.0 - eta, 1.0 + eta); - auto shape = make_contiguous(L.buffer()).shape().make_smooth(); - std::vector one_pm_eta_data(shape.size(), one_pm_eta); - Contiguous one_pm_eta_buffer(std::move(one_pm_eta_data), shape); - tensorwrapper::Tensor one_pm_eta_tensor(shape, - std::move(one_pm_eta_buffer)); - tensorwrapper::Tensor L_ball; - L_ball("i") = L("i") * one_pm_eta_tensor("i"); - return L_ball; -} - -} // namespace - -MODULE_CTOR(BallGeneralized) { - description(desc); - satisfies_property_type(); - - add_submodule("Eigen Solve"); -} - -MODULE_RUN(BallGeneralized) { - auto&& [A, B] = pt::unwrap_inputs(inputs); - - using uq_type = tensorwrapper::types::idouble; - using namespace tensorwrapper::buffer; - using namespace tensorwrapper::shape; - - auto A_buf = make_contiguous(A.buffer()); - auto B_buf = make_contiguous(B.buffer()); - auto shape = A_buf.shape().make_smooth(); - - // N.b., wrap_subdiagonalization will verify shapes match - auto n = shape.extent(0); - - Smooth vector_shape{n}; - Smooth matrix_shape{n, n}; - - auto eigen_solver_mod = submods.at("Eigen Solve"); - auto [values, vectors] = - wrap_subdiagonalization(A_buf, B_buf, eigen_solver_mod); - - // Here we need to convert the Eigen values and vectors to UQ type - auto [uq_values, uq_vectors] = convert_to_uq(values, vectors); - - auto H = compute_residual(A, B, uq_vectors, uq_values); - - auto eta_val = compute_eta(H); - - auto C_ball = t_ball(uq_vectors, eta_val); - auto L_ball = l_ball(uq_values, eta_val); - - auto rv = results(); - return pt::wrap_results(rv, L_ball, C_ball); -} -} // namespace scf::eigen_solver -#else -namespace scf::eigen_solver { -using pt = simde::GeneralizedEigenSolve; -MODULE_CTOR(BallGeneralized) { - description("Sigma was not enabled."); - satisfies_property_type(); - - add_submodule("Eigen Solve"); -} - -MODULE_RUN(BallGeneralized) { - throw std::runtime_error("Sigma was not enabled."); -} -} // namespace scf::eigen_solver -#endif diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp index 03166a1..73af583 100644 --- a/src/scf/eigen_solver/eigen_generalized.cpp +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -44,7 +44,7 @@ struct Kernel { using clean_t = std::decay_t; // Convert to Eigen buffers - if constexpr(tensorwrapper::types::is_interval_v) { + if constexpr(tensorwrapper::types::is_uq_type_v) { throw std::runtime_error( "EigenGeneralized Kernel: Interval types not supported"); } else { diff --git a/src/scf/eigen_solver/eigen_normal.cpp b/src/scf/eigen_solver/eigen_normal.cpp new file mode 100644 index 0000000..ae244ad --- /dev/null +++ b/src/scf/eigen_solver/eigen_normal.cpp @@ -0,0 +1,117 @@ +/* + * 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 { + +namespace { +struct Kernel { + std::size_t m_n_rows; + std::size_t m_n_cols; + + using tensor_t = simde::type::tensor; + using return_t = std::pair; + + Kernel(std::size_t n_rows, std::size_t n_cols) : + m_n_rows(n_rows), m_n_cols(n_cols) {} + + template + return_t operator()(const std::span& A) { + using clean_t = std::decay_t; + // Convert to Eigen buffers + + if constexpr(tensorwrapper::types::is_uq_type_v) { + throw std::runtime_error( + "EigenNormalKernel: Interval types not supported"); + } else { + // Wrap the tensors in Eigen::Map objects to avoid copy + const auto* pA = A.data(); + auto rows = m_n_rows; + auto cols = m_n_cols; + + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using clean_type = std::decay_t; + using matrix_type = + Eigen::Matrix; + using map_type = Eigen::Map; + + map_type A_map(pA, rows, cols); + + // Compute + Eigen::SelfAdjointEigenSolver es(A_map); + auto eigen_values = es.eigenvalues(); + auto eigen_vectors = es.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 tensorwrapper::buffer::make_contiguous; + + auto pvalues_buffer = make_contiguous(vector_shape); + auto pvectors_buffer = make_contiguous(matrix_shape); + + for(decltype(rows) i = 0; i < rows; ++i) { + pvalues_buffer.set_elem({i}, eigen_values(i)); + for(decltype(cols) j = 0; j < cols; ++j) { + pvectors_buffer.set_elem({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::EigenSolve; + +const auto desc = R"( + Eigen Solve via Eigen + --------------------------------- + + TODO: Write me!!! + )"; + +MODULE_CTOR(EigenNormal) { + description(desc); + satisfies_property_type(); +} + +MODULE_RUN(EigenNormal) { + auto&& [A] = pt::unwrap_inputs(inputs); + + using tensorwrapper::buffer::make_contiguous; + const auto& A_buffer = make_contiguous(A.buffer()); + const auto& A_shape = A_buffer.shape(); + Kernel k(A_shape.extent(0), A_shape.extent(1)); + using tensorwrapper::buffer::visit_contiguous_buffer; + auto [values, vectors] = visit_contiguous_buffer(k, A_buffer); + + auto rv = results(); + return pt::wrap_results(rv, values, vectors); +} + +} // namespace scf::eigen_solver diff --git a/src/scf/eigen_solver/eigen_solve_driver.cpp b/src/scf/eigen_solver/eigen_solve_driver.cpp new file mode 100644 index 0000000..2f81962 --- /dev/null +++ b/src/scf/eigen_solver/eigen_solve_driver.cpp @@ -0,0 +1,92 @@ +/* + * Copyright 2026 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 +#include +#include + +namespace scf::eigen_solver { + +namespace { + +const auto desc = R"( +Eigen Solve Driver +------------------ + +Routes normal eigenvalue problems to the appropriate backend based on the +floating-point type stored in the input matrix. Plain ``float``/``double`` +matrices use the deterministic Eigen backend; uncertain and interval types use +ball-arithmetic certification. +)"; + +constexpr const char* kNoneSubmodule = "none"; +constexpr const char* kUncertainSubmodule = "uncertain"; +constexpr const char* kIntervalSubmodule = "interval"; + +using pt = simde::EigenSolve; + +struct Router { + using return_t = std::tuple; + + const simde::type::tensor& m_A; + pluginplay::type::submodule_map& m_submods; + + Router(const simde::type::tensor& A, + pluginplay::type::submodule_map& submods) : + m_A(A), m_submods(submods) {} + + template + return_t operator()(const std::span&) { + using clean_t = std::decay_t; + + if constexpr(tensorwrapper::types::is_uncertain_v) { + return m_submods.at(kUncertainSubmodule).run_as(m_A); + } else if constexpr(tensorwrapper::types::is_uq_type_v) { + return m_submods.at(kIntervalSubmodule).run_as(m_A); + } else { + return m_submods.at(kNoneSubmodule).run_as(m_A); + } + } +}; + +} // namespace + +MODULE_CTOR(EigenSolveDriver) { + description(desc); + satisfies_property_type(); + + add_submodule(kNoneSubmodule); + add_submodule(kUncertainSubmodule); + add_submodule(kIntervalSubmodule); +} + +MODULE_RUN(EigenSolveDriver) { + auto&& [A] = pt::unwrap_inputs(inputs); + + using tensorwrapper::buffer::make_contiguous; + using tensorwrapper::buffer::visit_contiguous_buffer; + const auto& A_buffer = make_contiguous(A.buffer()); + Router router(A, submods); + auto [values, vectors] = visit_contiguous_buffer(router, A_buffer); + + auto rv = results(); + return pt::wrap_results(rv, values, vectors); +} + +} // namespace scf::eigen_solver diff --git a/src/scf/eigen_solver/eigen_solver.hpp b/src/scf/eigen_solver/eigen_solver.hpp index c9b753c..f5dd10e 100644 --- a/src/scf/eigen_solver/eigen_solver.hpp +++ b/src/scf/eigen_solver/eigen_solver.hpp @@ -19,18 +19,25 @@ namespace scf::eigen_solver { -DECLARE_MODULE(BallGeneralized); +DECLARE_MODULE(GeneralizedEigenSolver); +DECLARE_MODULE(EigenSolveDriver); DECLARE_MODULE(EigenGeneralized); +DECLARE_MODULE(EigenNormal); +DECLARE_MODULE(JacobiNormal); inline void set_defaults(pluginplay::ModuleManager& mm) { - mm.change_submod("Generalized eigensolve via Ball arithmetic", - "Eigen Solve", "Generalized eigensolve via Eigen"); + mm.change_submod("Eigen Solve", "none", "Eigen Solve via Eigen"); + mm.change_submod("Eigen Solve", "uncertain", "Eigen Solve via Jacobi"); + mm.change_submod("Eigen Solve", "interval", "Eigen Solve via Jacobi"); + mm.change_submod("Generalized eigensolve", "Eigen Solve", "Eigen Solve"); } inline void load_modules(pluginplay::ModuleManager& mm) { + mm.add_module("Eigen Solve"); + mm.add_module("Eigen Solve via Eigen"); + mm.add_module("Eigen Solve via Jacobi"); mm.add_module("Generalized eigensolve via Eigen"); - mm.add_module( - "Generalized eigensolve via Ball arithmetic"); + mm.add_module("Generalized eigensolve"); set_defaults(mm); } diff --git a/src/scf/eigen_solver/generalized_eigen_solver.cpp b/src/scf/eigen_solver/generalized_eigen_solver.cpp new file mode 100644 index 0000000..fc39964 --- /dev/null +++ b/src/scf/eigen_solver/generalized_eigen_solver.cpp @@ -0,0 +1,75 @@ +/* + * Copyright 2026 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 +#include + +using pt = simde::GeneralizedEigenSolve; +using pt_normal = simde::EigenSolve; +namespace scf::eigen_solver { +namespace { +const auto desc = R"( +Generalized Eigen Solve +----------------------- + +TODO: Write me!!! +)"; +} + +MODULE_CTOR(GeneralizedEigenSolver) { + description(desc); + satisfies_property_type(); + + add_submodule("Eigen Solve"); +} + +MODULE_RUN(GeneralizedEigenSolver) { + auto&& [A, B] = pt::unwrap_inputs(inputs); + + auto eigen_solver_mod = submods.at("Eigen Solve"); + + // Step 1: Diagonalize B to get B_values and B_vectors + auto [B_values, B_vectors] = eigen_solver_mod.run_as(B); + + // Step 2: Compute (B_values)**-1/2 + using tensorwrapper::operations::power; + using tensorwrapper::utilities::diagonal_matrix; + auto B_values_inv_sqrt = power(B_values, -0.5); + auto B_values_matrix = diagonal_matrix(B_values_inv_sqrt); + + // Step 3: Compute C = B_vectors * (B_values)**-1/2 + simde::type::tensor C; + C("i,k") = B_vectors("i,j") * B_values_matrix("j,k"); + + // Step 4: A' = C^T * A * C + simde::type::tensor CA, A_prime; + CA("i,k") = C("j,i") * A("j,k"); + A_prime("i,k") = CA("i,j") * C("j,k"); + + // Step 5: Diagonalize A' to get A_values and A'_vectors + auto [A_values, A_vectors] = eigen_solver_mod.run_as(A_prime); + + // Step 6: A_vectors = C * A'_vectors + simde::type::tensor evectors; + evectors("i,k") = C("i,j") * A_vectors("j,k"); + + auto rv = results(); + return pt::wrap_results(rv, A_values, evectors); +} +} // namespace scf::eigen_solver diff --git a/src/scf/eigen_solver/jacobi_eigen_helpers.hpp b/src/scf/eigen_solver/jacobi_eigen_helpers.hpp new file mode 100644 index 0000000..d9dcffa --- /dev/null +++ b/src/scf/eigen_solver/jacobi_eigen_helpers.hpp @@ -0,0 +1,255 @@ +/* + * Copyright 2026 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +using dynamic_matrix = Eigen::Matrix; + +namespace scf::eigen_solver::detail { + +template> +struct value_type_impl { + using type = T; +}; + +template +struct value_type_impl { + using type = typename T::value_t; +}; + +template +using value_type = typename value_type_impl::type; + +// Computes the Frobenius norm of the off-diagonal elements of n by n matrix S. +template +T off_diagonal_frobenius(const dynamic_matrix& S, std::size_t n) { + using tensorwrapper::types::uq_center; + value_type frob(0); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + if(i != j) { + auto x = uq_center(S(i, j)); + frob += x * x; + } + } + } + return T(std::sqrt(frob)); +} + +/** @brief Computes the cosine and sine of a Jacobi rotation for the p, q plane. + * + * @param[in] s_pp_uq The (p, p) element of the matrix to be diagonalized. + * @param[in] s_pq_uq The (p, q) element of the matrix to be diagonalized. + * @param[in] s_qq_uq The (q, q) element of the matrix to be diagonalized. + * @param[out] c_out Output parameter for the cosine of the rotation. + * @param[out] s_out Output parameter for the sine of the rotation. + * + * @return True if the rotation was successfully computed, false if the off- + * diagonal element is too small. + */ +template +bool make_jacobi(const T& s_pp_uq, const T& s_pq_uq, const T& s_qq_uq, T& c_out, + T& s_out) { + using value_t = value_type; + using tensorwrapper::types::uq_center; + + const value_t s_pp = uq_center(s_pp_uq); + const value_t s_pq = uq_center(s_pq_uq); + const value_t s_qq = uq_center(s_qq_uq); + + const value_t deno = value_t(2) * std::abs(s_pq); + // If the off-diagonal element is already almost zero, the angle will be + // zero, so we can skip the rotation by just returning cos(0) = 1 and + // sin(0)=0. + if(deno < std::numeric_limits::min()) { + c_out = T(1); + s_out = T(0); + return false; + } + + // tau = cot(2*theta) + const value_t tau = (s_pp - s_qq) / deno; + + // w = |csc(2*theta)| = sqrt(tau^2 + 1) + const value_t w = std::sqrt(tau * tau + value_t(1)); + + // t = tan(theta) + value_t t; + if(tau > value_t(0)) { + t = value_t(1) / (tau + w); + } else { + t = value_t(1) / (tau - w); + } + const value_t sign_t = t > value_t(0) ? value_t(1) : value_t(-1); + const value_t n = value_t(1) / std::sqrt(t * t + value_t(1)); + const value_t y_sign = s_pq >= value_t(0) ? value_t(1) : value_t(-1); + const value_t s = -sign_t * y_sign * std::abs(t) * n; + const value_t c = n; + c_out = T(c); + s_out = T(s); + return true; +} + +/** @brief Applies the adjoint of the Jacobi rotation to the left of a matrix. + * + * @param[in,out] mat The matrix to which the rotation is applied. + * @param[in] p The first index of the rotation plane. + * @param[in] q The second index of the rotation plane. + * @param[in] c The cosine of the rotation. + * @param[in] s The sine of the rotation. + */ +template +void apply_on_the_left_adjoint(dynamic_matrix& mat, Eigen::Index p, + Eigen::Index q, const T& c, const T& s) { + const Eigen::Index n_cols = mat.cols(); + for(Eigen::Index j = 0; j < n_cols; ++j) { + const T xp = mat(p, j); + const T xq = mat(q, j); + const T new_p = c * xp - s * xq; + const T new_q = s * xp + c * xq; + mat(p, j) = new_p; + mat(q, j) = new_q; + } +} + +/** @brief Applies the Jacobi rotation to the right of a matrix. + * + * @param[in,out] mat The matrix to which the rotation is applied. + * @param[in] p The first index of the rotation plane. + * @param[in] q The second index of the rotation plane. + * @param[in] c The cosine of the rotation. + * @param[in] s The sine of the rotation. + */ +template +void apply_on_the_right(dynamic_matrix& mat, Eigen::Index p, Eigen::Index q, + const T& c, const T& s) { + const Eigen::Index n_rows = mat.rows(); + for(Eigen::Index i = 0; i < n_rows; ++i) { + const T xp = mat(i, p); + const T xq = mat(i, q); + const T new_p = c * xp - s * xq; + const T new_q = s * xp + c * xq; + mat(i, p) = new_p; + mat(i, q) = new_q; + } +} + +// Diagaonalizes a matrix via V^T A V, extracts diagonal as a vector. +template +std::vector rayleigh_diagonal(const dynamic_matrix& V, + const dynamic_matrix& A_orig, + std::size_t n) { + std::vector diag(n); + const dynamic_matrix rayleigh = V.transpose() * A_orig * V; + for(std::size_t j = 0; j < n; ++j) { diag[j] = rayleigh(j, j); } + return diag; +} + +/** @brief Performs one sweep of Jacobi rotations. + * + * A sweep consists of applying a Jacobi rotation to every pair of indices + * (p, q). + * + * @param[in,out] S The matrix to be diagonalized. Gets modified in-place. + * @param[in,out] V The matrix that accumulates the rotations. Gets modified + * in-place. + * @param[in] n The size of the matrix S (assumed to be n by n). + */ +template +void jacobi_sweep(dynamic_matrix& S, dynamic_matrix& V, std::size_t n) { + for(std::size_t p = 0; p < n; ++p) { + for(std::size_t q = p + 1; q < n; ++q) { + T c; + T s; + if(make_jacobi(S(p, p), S(p, q), S(q, q), c, s)) { + apply_on_the_left_adjoint(S, p, q, c, s); + apply_on_the_right(S, p, q, c, s); + apply_on_the_right(V, p, q, c, s); + S(p, q) = T(0); + S(q, p) = T(0); + } + } + } +} + +template +inline std::pair, std::vector> symmetric_jacobi_eigen( + std::span A, std::size_t n, double tol, std::size_t max_sweeps) { + using matrix_type = dynamic_matrix; + + // Copy of A, used to get eigenvalues via V^T A V at the end. We can't just + // use S because the Jacobi rotations are applied in-place to S, so S + // doesn't equal V^T A V until convergence. + matrix_type A_orig(n, n); + + // Iterative approximation to the eigenvalues of A. Converges to a diagonal + // matrix. + matrix_type S(n, n); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + A_orig(i, j) = A[i * n + j]; + S(i, j) = A[i * n + j]; + } + } + + // Accumulates the rotations. Converges to the eigenvectors of A. + matrix_type V = matrix_type::Identity(n, n); + using tensorwrapper::types::strictly_less; + using tensorwrapper::types::uq_center; + for(std::size_t sweep = 0; sweep < max_sweeps; ++sweep) { + const auto frob_sqrt = off_diagonal_frobenius(S, n); + if(strictly_less(frob_sqrt, tol)) { break; } + if(sweep == max_sweeps - 1) { + throw std::runtime_error("Jacobi algorithm did not converge"); + } + + jacobi_sweep(S, V, n); + // Symmetrize S to prevent numerical issues from destroying symmetry. + S = (S + S.transpose()) / T(2.0); + } + + const auto diag = rayleigh_diagonal(V, A_orig, n); + + // Order eigenvalue in increasing order + std::vector order(n); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](std::size_t i, std::size_t j) { + return uq_center(diag[i]) < uq_center(diag[j]); + }); + + // Fill in return values/vectors + std::vector values(n); + std::vector vectors(n * n); + for(std::size_t k = 0; k < n; ++k) { + const auto idx = order[k]; + values[k] = diag[idx]; + for(std::size_t i = 0; i < n; ++i) { vectors[i * n + k] = V(i, idx); } + } + return {std::move(values), std::move(vectors)}; +} + +} // namespace scf::eigen_solver::detail diff --git a/src/scf/eigen_solver/jacobi_normal.cpp b/src/scf/eigen_solver/jacobi_normal.cpp new file mode 100644 index 0000000..aa51616 --- /dev/null +++ b/src/scf/eigen_solver/jacobi_normal.cpp @@ -0,0 +1,94 @@ +/* + * Copyright 2026 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 "jacobi_eigen_helpers.hpp" +#include +#include + +namespace scf::eigen_solver { + +namespace { + +struct Kernel { + std::size_t m_n_rows; + std::size_t m_n_cols; + + using tensor_t = simde::type::tensor; + using return_t = std::pair; + + Kernel(std::size_t n_rows, std::size_t n_cols) : + m_n_rows(n_rows), m_n_cols(n_cols) {} + + template + return_t operator()(const std::span& A) { + using clean_t = std::decay_t; + const auto n = m_n_rows; + if(m_n_rows != m_n_cols) { + throw std::runtime_error("JacobiNormal: matrix must be square"); + } + double tol = []() { + if constexpr(tensorwrapper::types::is_uq_type_v) { + using value_t = typename clean_t::value_t; + return value_t(10) * std::numeric_limits::epsilon(); + } else { + return 10 * std::numeric_limits::epsilon(); + } + }(); + + const auto max_sweeps = + tensorwrapper::types::is_uq_type_v ? 2000 * n : 50 * n; + auto [evals, evecs] = + detail::symmetric_jacobi_eigen(A, n, tol, max_sweeps); + + using tensorwrapper::utilities::make_tensor; + auto values = make_tensor({n}, evals); + auto vectors = make_tensor({n, n}, evecs); + return std::make_pair(values, vectors); + } +}; +} // namespace + +using pt = simde::EigenSolve; + +const auto desc = R"( + Eigen Solve via Jacobi + --------------------------------- + + Symmetric eigen solve by cyclic Jacobi rotations using Eigen's + JacobiRotation class. + )"; + +MODULE_CTOR(JacobiNormal) { + description(desc); + satisfies_property_type(); +} + +MODULE_RUN(JacobiNormal) { + auto&& [A] = pt::unwrap_inputs(inputs); + + using tensorwrapper::buffer::make_contiguous; + const auto& A_buffer = make_contiguous(A.buffer()); + const auto& A_shape = A_buffer.shape(); + Kernel k(A_shape.extent(0), A_shape.extent(1)); + using tensorwrapper::buffer::visit_contiguous_buffer; + auto [values, vectors] = visit_contiguous_buffer(k, A_buffer); + + auto rv = results(); + return pt::wrap_results(rv, values, vectors); +} + +} // namespace scf::eigen_solver diff --git a/src/scf/xc/libxc/libxc.cpp b/src/scf/xc/libxc/libxc.cpp index 998ff18..9e5d43f 100644 --- a/src/scf/xc/libxc/libxc.cpp +++ b/src/scf/xc/libxc/libxc.cpp @@ -293,7 +293,7 @@ struct BatchedDotKernel { // AOs on rows, grid points on columns if(m_sum_row) { for(std::size_t grid_i = 0; grid_i < m_n_grid; ++grid_i) { - clean_type sum = 0.0; + clean_type sum(0.0); for(std::size_t ao_i = 0; ao_i < m_n_aos; ++ao_i) { const auto idx = ao_i * m_n_grid + grid_i; sum += aos_on_grid[idx] * X[idx]; @@ -302,7 +302,7 @@ struct BatchedDotKernel { } } else { for(std::size_t ao_i = 0; ao_i < m_n_aos; ++ao_i) { - clean_type sum = 0.0; + clean_type sum(0.0); for(std::size_t grid_i = 0; grid_i < m_n_grid; ++grid_i) { const auto idx = ao_i * m_n_grid + grid_i; sum += aos_on_grid[idx] * X[idx]; diff --git a/tests/cxx/integration_tests/driver/scf_driver.cpp b/tests/cxx/integration_tests/driver/scf_driver.cpp index ee70754..94be890 100644 --- a/tests/cxx/integration_tests/driver/scf_driver.cpp +++ b/tests/cxx/integration_tests/driver/scf_driver.cpp @@ -20,7 +20,8 @@ using pt = simde::AOEnergy; using tensorwrapper::operations::approximately_equal; using types = std::tuple>; + tensorwrapper::types::interval_type, + tensorwrapper::types::thresholded_affine_type>; TEMPLATE_LIST_TEST_CASE("SCFDriver", "", types) { using float_type = TestType; @@ -53,8 +54,7 @@ TEMPLATE_LIST_TEST_CASE("SCFDriver", "", types) { } SECTION("DFT") { // GauXC not currently compatible with Uncertain values - if constexpr(!tensorwrapper::types::is_uncertain_v && - !tensorwrapper::types::is_interval_v) { + if constexpr(!tensorwrapper::types::is_uq_type_v) { auto func = chemist::qm_operator::xc_functional::PBE; const auto RKS_op = "Restricted Kohn-Sham Op"; const auto rks_op = "Restricted One-Electron Kohn-Sham Op"; diff --git a/tests/cxx/integration_tests/integration_tests.hpp b/tests/cxx/integration_tests/integration_tests.hpp index 3b0b979..0e5ddc6 100644 --- a/tests/cxx/integration_tests/integration_tests.hpp +++ b/tests/cxx/integration_tests/integration_tests.hpp @@ -33,12 +33,9 @@ namespace { mm.change_input("Kinetic", "UQ Type", uq_type); mm.change_input("Nuclear", "UQ Type", uq_type); mm.change_input("sto-3g atomic density matrix", "With UQ?", true); - if(uq_type == "interval") { - mm.change_submod("Loop", "Diagonalizer", - "Generalized eigensolve via Ball Arithmetic"); - mm.change_submod("Diagonalization Fock Update", "Diagonalizer", - "Generalized eigensolve via Ball Arithmetic"); - } + mm.change_submod("Loop", "Diagonalizer", "Generalized eigensolve"); + mm.change_submod("Diagonalization Fock update", "Diagonalizer", + "Generalized eigensolve"); } } // namespace @@ -71,17 +68,15 @@ pluginplay::ModuleManager load_modules() { configure_uq(mm, "uncertain"); } else if constexpr(tensorwrapper::types::is_interval_v) { std::string key = "interval"; - mm.at("Generalized eigensolve via Eigen").turn_off_memoization(); - mm.at("Density matrix builder").turn_off_memoization(); - mm.at("Electronic energy").turn_off_memoization(); - mm.at("Fock matrix builder").turn_off_memoization(); - mm.at("Restricted One-Electron Fock Op").turn_off_memoization(); - mm.at("Restricted Fock Op").turn_off_memoization(); - mm.at("Four center J builder").turn_off_memoization(); - mm.at("Four center K builder").turn_off_memoization(); - configure_uq(mm, "interval"); + configure_uq(mm, key); + } else if constexpr(tensorwrapper::types::is_affine_v) { + std::string key = "affine"; + configure_uq(mm, key); + } else if constexpr(tensorwrapper::types::is_thresholded_affine_v< + FloatType>) { + std::string key = "thresholded affine"; + configure_uq(mm, key); } - return mm; } diff --git a/tests/cxx/unit_tests/eigen_solver/ball_generalized.cpp b/tests/cxx/unit_tests/eigen_solver/ball_generalized.cpp deleted file mode 100644 index 691371d..0000000 --- a/tests/cxx/unit_tests/eigen_solver/ball_generalized.cpp +++ /dev/null @@ -1,318 +0,0 @@ -/* - * 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 -#ifdef ENABLE_SIGMA -namespace { - -// 7x7 Fock / overlap-style fixtures and expected generalized eigenvalues -// (L_ball) transcribed from console output: each printed "m+/-r" token is -// stored as idouble(lower, upper) with lower = m - r, upper = m + r. - -inline simde::type::tensor uncertain_A() { - using uq_type = tensorwrapper::types::idouble; - tensorwrapper::shape::Smooth shape{7, 7}; - auto buf = tensorwrapper::buffer::make_contiguous(shape); - buf.set_elem({0, 0}, - uq_type{-5.9761941965020215e-01, -5.9761890608307133e-01}); - buf.set_elem({0, 1}, - uq_type{-4.0307029666324207e-01, -4.0306978309613190e-01}); - buf.set_elem({0, 2}, - uq_type{-1.2374285072777236e+00, -1.2374279937106170e+00}); - buf.set_elem({0, 3}, - uq_type{-1.0494685522547358e+00, -1.0494680386876154e+00}); - buf.set_elem({0, 4}, - uq_type{2.1281939270258100e-02, 2.1282452837359501e-02}); - buf.set_elem({0, 5}, - uq_type{-5.0779585063444665e-01, -5.0779533706733271e-01}); - buf.set_elem({0, 6}, - uq_type{1.9650251601842519e-01, 1.9650302958553259e-01}); - buf.set_elem({1, 0}, - uq_type{-4.0307029666324190e-01, -4.0306978309613173e-01}); - buf.set_elem({1, 1}, - uq_type{-5.9047593004366428e-01, -5.9047541647653523e-01}); - buf.set_elem({1, 2}, - uq_type{-1.2227376779992984e+00, -1.2227371644321914e+00}); - buf.set_elem({1, 3}, - uq_type{-1.0408334164996009e+00, -1.0408329029324810e+00}); - buf.set_elem({1, 4}, - uq_type{-1.7104644978368001e-03, -1.7099509307360000e-03}); - buf.set_elem({1, 5}, - uq_type{-1.4859304502705128e-01, -1.4859253145994750e-01}); - buf.set_elem({1, 6}, - uq_type{-5.2200750954699371e-01, -5.2200699597987887e-01}); - buf.set_elem({2, 0}, - uq_type{-1.2374285072777242e+00, -1.2374279937106172e+00}); - buf.set_elem({2, 1}, - uq_type{-1.2227376779992984e+00, -1.2227371644321918e+00}); - buf.set_elem({2, 2}, - uq_type{-2.0229432721235455e+01, -2.0229432207668175e+01}); - buf.set_elem({2, 3}, - uq_type{-5.1628583778067823e+00, -5.1628578642396556e+00}); - buf.set_elem({2, 4}, - uq_type{7.9385631994439998e-04, 7.9436988704520005e-04}); - buf.set_elem({2, 5}, - uq_type{-2.6603689531534101e-02, -2.6603175964432901e-02}); - buf.set_elem({2, 6}, - uq_type{-1.3130169293639199e-02, -1.3129655726538200e-02}); - buf.set_elem({3, 0}, - uq_type{-1.0494685522547353e+00, -1.0494680386876150e+00}); - buf.set_elem({3, 1}, - uq_type{-1.0408334164996000e+00, -1.0408329029324805e+00}); - buf.set_elem({3, 2}, - uq_type{-5.1628583778067814e+00, -5.1628578642396548e+00}); - buf.set_elem({3, 3}, - uq_type{-2.4488877696649971e+00, -2.4488872560978492e+00}); - buf.set_elem({3, 4}, - uq_type{3.5224479263691001e-03, 3.5229614934699002e-03}); - buf.set_elem({3, 5}, - uq_type{-1.1753004101803359e-01, -1.1752952745093080e-01}); - buf.set_elem({3, 6}, - uq_type{-5.6952104344840299e-02, -5.6951590777738301e-02}); - buf.set_elem({4, 0}, - uq_type{2.1281939270258100e-02, 2.1282452837359501e-02}); - buf.set_elem({4, 1}, - uq_type{-1.7104644978368001e-03, -1.7099509307360000e-03}); - buf.set_elem({4, 2}, - uq_type{7.9385631994439998e-04, 7.9436988704520005e-04}); - buf.set_elem({4, 3}, - uq_type{3.5224479263691001e-03, 3.5229614934699002e-03}); - buf.set_elem({4, 4}, - uq_type{-3.9094388879395919e-01, -3.9094337522681283e-01}); - buf.set_elem({4, 5}, - uq_type{-1.5239204150682001e-03, -1.5234068479674000e-03}); - buf.set_elem({4, 6}, - uq_type{1.1589217834762000e-03, 1.1594353505770000e-03}); - buf.set_elem({5, 0}, - uq_type{-5.0779585063444688e-01, -5.0779533706733360e-01}); - buf.set_elem({5, 1}, - uq_type{-1.4859304502705112e-01, -1.4859253145994750e-01}); - buf.set_elem({5, 2}, - uq_type{-2.6603689531534198e-02, -2.6603175964433200e-02}); - buf.set_elem({5, 3}, - uq_type{-1.1753004101803349e-01, -1.1752952745093090e-01}); - buf.set_elem({5, 4}, - uq_type{-1.5239204150682001e-03, -1.5234068479674000e-03}); - buf.set_elem({5, 5}, - uq_type{-3.5326421481592180e-01, -3.5326370124877843e-01}); - buf.set_elem({5, 6}, - uq_type{-1.0335912128416299e-02, -1.0335398561315099e-02}); - buf.set_elem({6, 0}, - uq_type{1.9650251601842419e-01, 1.9650302958553159e-01}); - buf.set_elem({6, 1}, - uq_type{-5.2200750954699304e-01, -5.2200699597987910e-01}); - buf.set_elem({6, 2}, - uq_type{-1.3130169293640000e-02, -1.3129655726539001e-02}); - buf.set_elem({6, 3}, - uq_type{-5.6952104344839098e-02, -5.6951590777737100e-02}); - buf.set_elem({6, 4}, - uq_type{1.1589217834762000e-03, 1.1594353505770000e-03}); - buf.set_elem({6, 5}, - uq_type{-1.0335912128416299e-02, -1.0335398561315099e-02}); - buf.set_elem({6, 6}, - uq_type{-3.3401382547702690e-01, -3.3401331190988387e-01}); - return simde::type::tensor(shape, std::move(buf)); -} - -inline simde::type::tensor uncertain_B() { - using uq_type = tensorwrapper::types::idouble; - tensorwrapper::shape::Smooth shape{7, 7}; - auto buf = tensorwrapper::buffer::make_contiguous(shape); - buf.set_elem({0, 0}, - uq_type{9.9999999999999978e-01, 1.0000000000000002e+00}); - buf.set_elem({0, 1}, - uq_type{2.5358257133619461e-01, 2.5358257133619483e-01}); - buf.set_elem({0, 2}, - uq_type{5.5900804544533403e-02, 5.5900804544533597e-02}); - buf.set_elem({0, 3}, - uq_type{4.8454762713240557e-01, 4.8454762713240579e-01}); - buf.set_elem({0, 4}, - uq_type{-1.5482109990821101e-02, -1.5482109990820900e-02}); - buf.set_elem({0, 5}, - uq_type{3.5644036709221749e-01, 3.5644036709221771e-01}); - buf.set_elem({0, 6}, - uq_type{-1.7760045089221962e-01, -1.7760045089221940e-01}); - buf.set_elem({1, 0}, - uq_type{2.5358257133619461e-01, 2.5358257133619483e-01}); - buf.set_elem({1, 1}, - uq_type{9.9999999999999978e-01, 1.0000000000000002e+00}); - buf.set_elem({1, 2}, - uq_type{5.5230922941791903e-02, 5.5230922941792097e-02}); - buf.set_elem({1, 3}, - uq_type{4.8117972220632926e-01, 4.8117972220632949e-01}); - buf.set_elem({1, 4}, - uq_type{2.6830766568749997e-03, 2.6830766568752000e-03}); - buf.set_elem({1, 5}, - uq_type{7.2658318900858498e-02, 7.2658318900858693e-02}); - buf.set_elem({1, 6}, - uq_type{3.9007755239467901e-01, 3.9007755239467923e-01}); - buf.set_elem({2, 0}, - uq_type{5.5900804544533403e-02, 5.5900804544533597e-02}); - buf.set_elem({2, 1}, - uq_type{5.5230922941791903e-02, 5.5230922941792097e-02}); - buf.set_elem({2, 2}, - uq_type{1.0000000000000002e+00, 1.0000000000000007e+00}); - buf.set_elem({2, 3}, - uq_type{2.3670392057272610e-01, 2.3670392057272632e-01}); - buf.set_elem({2, 4}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({2, 5}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({2, 6}, - uq_type{-2.0000000000000000e-16, 0.0000000000000000e+00}); - buf.set_elem({3, 0}, - uq_type{4.8454762713240557e-01, 4.8454762713240579e-01}); - buf.set_elem({3, 1}, - uq_type{4.8117972220632926e-01, 4.8117972220632949e-01}); - buf.set_elem({3, 2}, - uq_type{2.3670392057272618e-01, 2.3670392057272641e-01}); - buf.set_elem({3, 3}, - uq_type{9.9999999999999967e-01, 9.9999999999999989e-01}); - buf.set_elem({3, 4}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({3, 5}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({3, 6}, - uq_type{-7.9999999999999998e-16, -6.0000000000000009e-16}); - buf.set_elem({4, 0}, - uq_type{-1.5482109990821101e-02, -1.5482109990820900e-02}); - buf.set_elem({4, 1}, - uq_type{2.6830766568749997e-03, 2.6830766568752000e-03}); - buf.set_elem({4, 2}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({4, 3}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({4, 4}, - uq_type{1.0000000000000000e+00, 1.0000000000000004e+00}); - buf.set_elem({4, 5}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({4, 6}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({5, 0}, - uq_type{3.5644036709221738e-01, 3.5644036709221760e-01}); - buf.set_elem({5, 1}, - uq_type{7.2658318900858498e-02, 7.2658318900858693e-02}); - buf.set_elem({5, 2}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({5, 3}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({5, 4}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({5, 5}, - uq_type{9.9999999999999978e-01, 1.0000000000000002e+00}); - buf.set_elem({5, 6}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({6, 0}, - uq_type{-1.7760045089221932e-01, -1.7760045089221910e-01}); - buf.set_elem({6, 1}, - uq_type{3.9007755239467889e-01, 3.9007755239467912e-01}); - buf.set_elem({6, 2}, - uq_type{-2.9999999999999999e-16, -9.9999999999999998e-17}); - buf.set_elem({6, 3}, - uq_type{-9.0000000000000003e-16, -6.9999999999999994e-16}); - buf.set_elem({6, 4}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({6, 5}, - uq_type{-9.9999999999999998e-17, 9.9999999999999998e-17}); - buf.set_elem({6, 6}, - uq_type{1.0000000000000000e+00, 1.0000000000000004e+00}); - return simde::type::tensor(shape, std::move(buf)); -} - -// Reference L_ball from pasted console output (m+/-0 -> idouble(m,m)). -inline simde::type::tensor l_ball_corr() { - using uq_type = tensorwrapper::types::idouble; - tensorwrapper::shape::Smooth shape{7}; - auto buf = tensorwrapper::buffer::make_contiguous(shape); - buf.set_elem({0}, uq_type{-20.2383215358238289, -20.2383215358238289}); - buf.set_elem({1}, uq_type{-1.2730342062134834, -1.2730342062134834}); - buf.set_elem({2}, uq_type{-0.6267577709783282, -0.6267577709783282}); - buf.set_elem({3}, uq_type{-0.4510422475343768, -0.4510422475343768}); - buf.set_elem({4}, uq_type{-0.3910154399987634, -0.3910154399987634}); - buf.set_elem({5}, uq_type{0.6153219563024107, 0.6153219563024107}); - buf.set_elem({6}, uq_type{0.7613504962941617, 0.7613504962941617}); - return simde::type::tensor(shape, std::move(buf)); -} - -} // namespace - -using types = std::tuple; -TEMPLATE_LIST_TEST_CASE("BallGeneralized", "", types) { - using uq_type = TestType; - using float_type = typename uq_type::value_t; - pluginplay::ModuleManager mm; - scf::load_modules(mm); - - using pt = simde::GeneralizedEigenSolve; - - auto& mod = mm.at("Generalized eigensolve via Ball Arithmetic"); - - SECTION("2 by 2 test case") { - tensorwrapper::shape::Smooth shape{2, 2}; - using tensorwrapper::buffer::make_contiguous; - auto A_buffer = make_contiguous(shape); - float_type noise = 0.001; - - A_buffer.set_elem({0, 0}, uq_type{1.0 - noise, 1.0 + noise}); - A_buffer.set_elem({0, 1}, uq_type{2.0 - noise, 2.0 + noise}); - A_buffer.set_elem({1, 0}, uq_type{2.0 - noise, 2.0 + noise}); - A_buffer.set_elem({1, 1}, uq_type{3.0 - noise, 3.0 + noise}); - - auto B_buffer = make_contiguous(shape); - B_buffer.set_elem({0, 0}, uq_type{1.0 - noise, 1 + noise}); - B_buffer.set_elem({0, 1}, uq_type{0.0 - noise, 0.0 + noise}); - B_buffer.set_elem({1, 0}, uq_type{0.0 - noise, 0.0 + noise}); - B_buffer.set_elem({1, 1}, uq_type{1.0 - noise, 1.0 + noise}); - - 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); - - std::vector expected_values{-0.236068, 4.236068}; - tensorwrapper::shape::Smooth corr_shape{2}; - tensorwrapper::buffer::Contiguous corr_buffer(expected_values, - corr_shape); - auto value_buffer = make_contiguous(values.buffer()); - for(std::size_t i = 0; i < corr_buffer.shape().size(); ++i) { - auto corr_value = corr_buffer.get_elem({i}); - auto value = value_buffer.get_elem({i}); - auto corr_uq = wtf::fp::float_cast(corr_value); - auto value_uq = wtf::fp::float_cast(value); - REQUIRE(value_uq.contains(corr_uq.median())); - } - } - - SECTION("7 by 7 test case") { - auto A = uncertain_A(); - auto B = uncertain_B(); - auto&& [values, vector] = mod.run_as(A, B); - auto L_ball = l_ball_corr(); - - auto value_buffer = make_contiguous(values.buffer()); - auto corr_buffer = make_contiguous(L_ball.buffer()); - for(std::size_t i = 0; i < corr_buffer.shape().size(); ++i) { - auto corr_value = corr_buffer.get_elem({i}); - auto value = value_buffer.get_elem({i}); - auto corr_uq = wtf::fp::float_cast(corr_value); - auto value_uq = wtf::fp::float_cast(value); - REQUIRE(value_uq.contains(corr_uq.median())); - } - } -} -#endif diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp index fc0a287..9c8388f 100644 --- a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp +++ b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp @@ -14,43 +14,24 @@ * limitations under the License. */ -#include "../../test_scf.hpp" -#include -#include +#include "h2_dimer_pencil.hpp" +#include "test_eigen_solver.hpp" using types = std::tuple; -TEMPLATE_LIST_TEST_CASE("EigenGeneralized", "", types) { - using float_type = TestType; - pluginplay::ModuleManager mm; - scf::load_modules(mm); +using namespace test_eigen_solver; +TEMPLATE_LIST_TEST_CASE("EigenGeneralized H2 dimer", "", types) { using pt = simde::GeneralizedEigenSolve; + pluginplay::ModuleManager mm; + scf::load_modules(mm); - auto& mod = mm.at("Generalized eigensolve via Eigen"); - - tensorwrapper::shape::Smooth shape{2, 2}; - using tensorwrapper::buffer::make_contiguous; - auto A_buffer = make_contiguous(shape); - A_buffer.set_elem({0, 0}, float_type{1.0}); - A_buffer.set_elem({0, 1}, float_type{2.0}); - A_buffer.set_elem({1, 0}, float_type{2.0}); - A_buffer.set_elem({1, 1}, float_type{3.0}); - - auto B_buffer = make_contiguous(shape); - B_buffer.set_elem({0, 0}, float_type{1.0}); - B_buffer.set_elem({0, 1}, float_type{0.0}); - B_buffer.set_elem({1, 0}, float_type{0.0}); - B_buffer.set_elem({1, 1}, float_type{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); - - std::vector expected_values{-0.236068, 4.236068}; - tensorwrapper::shape::Smooth corr_shape{2}; - tensorwrapper::buffer::Contiguous corr_buffer(expected_values, corr_shape); - simde::type::tensor corr(corr_shape, std::move(corr_buffer)); + auto rtol = std::is_same_v ? 5e-4 : 1e-5; + auto A = h2_dimer_fock_as(); + auto B = h2_dimer_overlap_as(); - REQUIRE(tensorwrapper::operations::approximately_equal(corr, values, 1E-6)); + auto& mod = mm.at("Generalized eigensolve via Eigen"); + auto [values, vectors] = mod.run_as(A, B); + auto eval_corr = h2_dimer_evals(); + require_eigenvalues_approx(values, eval_corr, rtol); + require_eigenpair_residual(A, values, vectors, rtol); } diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_normal.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_normal.cpp new file mode 100644 index 0000000..8c1b88c --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/eigen_normal.cpp @@ -0,0 +1,63 @@ +/* + * 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_eigen_solver.hpp" + +using types = std::tuple; +using namespace test_eigen_solver; +using namespace tensorwrapper::generate; + +TEMPLATE_LIST_TEST_CASE("EigenNormal", "", types) { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + + auto& mod = mm.at("Eigen Solve via Eigen"); + auto rtol = std::is_same_v ? 1e-3 : 1e-5; + using pt = simde::EigenSolve; + SECTION("classic 2 by 2") { + auto system = test_eigen_solver::classic_2x2(); + auto [values, vectors] = mod.run_as(system.matrix); + + require_eigenvalues_approx(values, system.eigenvalues, rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); + } + + SECTION("generated n=4 condition number 1e3") { + test_eigen_solver::SymmetricMatrixSpec spec; + spec.n = 4; + spec.condition_number = 1e3; + spec.spacing = test_eigen_solver::EigenvalueSpacing::Linear; + spec.seed = 11; + auto system = generate_eigen_system(spec); + auto [values, vectors] = mod.run_as(system.matrix); + require_eigenvalues_approx(values, system.eigenvalues, rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); + } + + SECTION("generated clustered n=6") { + test_eigen_solver::SymmetricMatrixSpec spec; + spec.n = 6; + spec.condition_number = 100.0; + spec.spacing = test_eigen_solver::EigenvalueSpacing::Clustered; + spec.n_clusters = 3; + spec.cluster_width = 1e-8; + spec.seed = 23; + auto system = generate_eigen_system(spec); + auto [values, vectors] = mod.run_as(system.matrix); + require_eigenvalues_approx(values, system.eigenvalues, rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); + } +} diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_solve_driver.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_solve_driver.cpp new file mode 100644 index 0000000..5538868 --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/eigen_solve_driver.cpp @@ -0,0 +1,33 @@ +/* + * Copyright 2026 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_eigen_solver.hpp" + +using types = std::tuple; +using namespace test_eigen_solver; + +TEMPLATE_LIST_TEST_CASE("EigenSolveDriver", "", types) { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("Eigen Solve"); + using pt = simde::EigenSolve; + + auto rtol = std::is_same_v ? 5e-4 : 1e-5; + auto system = classic_2x2(); + auto [values, vectors] = mod.run_as(system.matrix); + require_eigenvalues_approx(values, system.eigenvalues, rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); +} diff --git a/tests/cxx/unit_tests/eigen_solver/generalized_eigen_solver.cpp b/tests/cxx/unit_tests/eigen_solver/generalized_eigen_solver.cpp new file mode 100644 index 0000000..89cd2b1 --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/generalized_eigen_solver.cpp @@ -0,0 +1,37 @@ +/* + * 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 "h2_dimer_pencil.hpp" +#include "test_eigen_solver.hpp" + +using types = std::tuple; +using namespace test_eigen_solver; + +TEMPLATE_LIST_TEST_CASE("GeneralizedEigenSolver H2 dimer", "", types) { + using pt = simde::GeneralizedEigenSolve; + pluginplay::ModuleManager mm; + scf::load_modules(mm); + + auto rtol = std::is_same_v ? 5e-4 : 1e-5; + auto A = h2_dimer_fock_as(); + auto B = h2_dimer_overlap_as(); + + auto& mod = mm.at("Generalized eigensolve"); + auto [values, vectors] = mod.run_as(A, B); + auto eval_corr = h2_dimer_evals(); + require_eigenvalues_approx(values, eval_corr, rtol); + require_eigenpair_residual(A, values, vectors, rtol); +} diff --git a/tests/cxx/unit_tests/eigen_solver/h2_dimer_pencil.hpp b/tests/cxx/unit_tests/eigen_solver/h2_dimer_pencil.hpp new file mode 100644 index 0000000..b7b79dc --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/h2_dimer_pencil.hpp @@ -0,0 +1,77 @@ +/* + * Copyright 2026 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 +#include +#include +#include +#include +#include + +namespace test_eigen_solver { + +// H2 dimer (4 H, STO-3G) Fock and overlap from the core-guess +// diagonalization step. Generated by the SCFDriver integration test dump. + +inline constexpr std::array h2_dimer_fock_data{ + -1.5208138831078950e+00, -1.2688563238876407e+00, -1.4288125943995872e-01, + -3.5691640230145288e-02, -1.2688563238876407e+00, -1.6807172836050259e+00, + -4.4451594013857459e-01, -1.4288125943995875e-01, -1.4288125943995875e-01, + -4.4451594013857459e-01, -1.6807172836050259e+00, -1.2688563238876405e+00, + -3.5691640230145288e-02, -1.4288125943995872e-01, -1.2688563238876402e+00, + -1.5208138831078946e+00}; + +inline constexpr std::array h2_dimer_overlap_data{ + 1.0000000000000000e+00, 6.5987566030170408e-01, 6.7949947315035300e-02, + 1.5819246523949011e-02, 6.5987566030170408e-01, 1.0000000000000000e+00, + 2.2618964447525747e-01, 6.7949947315035300e-02, 6.7949947315035300e-02, + 2.2618964447525747e-01, 1.0000000000000000e+00, 6.5987566030170408e-01, + 1.5819246523949011e-02, 6.7949947315035300e-02, 6.5987566030170408e-01, + 1.0000000000000000e+00}; + +// Input Fock matrix +template +inline simde::type::tensor h2_dimer_fock_as() { + using tensorwrapper::utilities::make_tensor; + std::vector elements(h2_dimer_fock_data.begin(), + h2_dimer_fock_data.end()); + for(auto& x : elements) { x = static_cast(x); } + return make_tensor({4, 4}, std::move(elements)); +} + +// Input overlap matrix +template +inline simde::type::tensor h2_dimer_overlap_as() { + using tensorwrapper::utilities::make_tensor; + std::vector elements(h2_dimer_overlap_data.begin(), + h2_dimer_overlap_data.end()); + for(auto& x : elements) { x = static_cast(x); } + return make_tensor({4, 4}, std::move(elements)); +} + +// Correct eigenvalues +template +inline simde::type::tensor h2_dimer_evals() { + using tensorwrapper::utilities::make_tensor; + std::vector corr_data{static_cast(-1.7782489061355591), + static_cast(-1.6984246969223022), + static_cast(-1.0329644680023193), + static_cast(-0.8134391903877258)}; + return make_tensor({4}, corr_data); +} + +} // namespace test_eigen_solver diff --git a/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp b/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp new file mode 100644 index 0000000..4bef32f --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp @@ -0,0 +1,118 @@ +/* + * Copyright 2026 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_eigen_solver.hpp" + +using types = + std::tuple; +using namespace test_eigen_solver; +using namespace tensorwrapper::generate; + +TEMPLATE_LIST_TEST_CASE("JacobiNormal", "", types) { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + + auto& mod = mm.at("Eigen Solve via Jacobi"); + auto rtol = std::is_same_v ? 5e-4 : 1e-5; + using pt = simde::EigenSolve; + + SECTION("classic 2 by 2") { + auto system = classic_2x2(); + auto [values, vectors] = mod.run_as(system.matrix); + + require_eigenvalues_approx(values, system.eigenvalues, rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); + } + + SECTION("generated n=4 condition number 1e3") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.condition_number = 1e3; + spec.spacing = EigenvalueSpacing::Linear; + spec.seed = 11; + auto system = generate_eigen_system(spec); + auto [values, vectors] = mod.run_as(system.matrix); + require_eigenvalues_approx(values, system.eigenvalues, 10 * rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); + } + SECTION("generated clustered n=6") { + SymmetricMatrixSpec spec; + spec.n = 6; + spec.condition_number = 100.0; + spec.spacing = EigenvalueSpacing::Clustered; + spec.n_clusters = 3; + spec.cluster_width = 1e-3; + spec.seed = 23; + auto system = generate_eigen_system(spec); + auto [values, vectors] = mod.run_as(system.matrix); + require_eigenvalues_approx(values, system.eigenvalues, rtol); + require_eigenpair_residual(system.matrix, values, vectors, rtol); + } +} + +#ifdef ENABLE_SIGMA +using types2 = std::tuple; +TEMPLATE_LIST_TEST_CASE("JacobiNormal with noise", "", types2) { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + + auto& mod = mm.at("Eigen Solve via Jacobi"); + auto noise_level = 1e-6; + using pt = simde::EigenSolve; + SECTION("classic 2 by 2") { + auto system = classic_2x2(); + auto noisy_matrix = add_noise(system.matrix, noise_level); + auto [values, vectors] = mod.run_as(noisy_matrix); + require_uq_eigenvalues_contain(values, system.eigenvalues, + noise_level); + require_uq_eigenpair_residual(noisy_matrix, values, vectors, + noise_level); + } + + SECTION("generated n=4 condition number 1e3") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.condition_number = 1e3; + spec.spacing = EigenvalueSpacing::Linear; + spec.seed = 11; + auto system = generate_eigen_system(spec); + auto noisy_matrix = add_noise(system.matrix, noise_level); + auto [values, vectors] = mod.run_as(noisy_matrix); + require_uq_eigenvalues_contain(values, system.eigenvalues, + noise_level); + require_uq_eigenpair_residual(noisy_matrix, values, vectors, + noise_level); + } + + SECTION("generated clustered n=6") { + SymmetricMatrixSpec spec; + spec.n = 6; + spec.condition_number = 100.0; + spec.spacing = EigenvalueSpacing::Clustered; + spec.n_clusters = 3; + spec.cluster_width = 1e-8; + spec.seed = 23; + auto system = generate_eigen_system(spec); + auto noisy_matrix = add_noise(system.matrix, noise_level); + auto [values, vectors] = mod.run_as(noisy_matrix); + require_uq_eigenvalues_contain(values, system.eigenvalues, + noise_level); + require_uq_eigenpair_residual(noisy_matrix, values, vectors, + noise_level); + } +} +#endif diff --git a/tests/cxx/unit_tests/eigen_solver/test_eigen_solver.hpp b/tests/cxx/unit_tests/eigen_solver/test_eigen_solver.hpp new file mode 100644 index 0000000..7b86581 --- /dev/null +++ b/tests/cxx/unit_tests/eigen_solver/test_eigen_solver.hpp @@ -0,0 +1,125 @@ +/* + * Copyright 2026 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 "../../test_scf.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace test_eigen_solver { + +using EigenSystem = tensorwrapper::generate::EigenSystem; +using EigenvalueSpacing = tensorwrapper::generate::EigenvalueSpacing; +using SymmetricMatrixSpec = tensorwrapper::generate::SymmetricMatrixSpec; +constexpr std::size_t kMaxMatrixDim = tensorwrapper::generate::kMaxMatrixDim; + +// EigenSystem consistent with A = [[1,2],[2,1]] +template +inline EigenSystem classic_2x2() { + using tensorwrapper::utilities::make_tensor; + std::vector evals{static_cast(-0.2360680401325226), + static_cast(4.2360687255859375)}; + std::vector evecs{static_cast(-0.8506508469581604), + static_cast(-0.5257311463356018), + static_cast(0.5257311463356018), + static_cast(-0.8506508469581604)}; + std::vector mat{static_cast(1.0), static_cast(2.0), + static_cast(2.0), static_cast(3.0)}; + EigenSystem rv; + rv.n = 2; + rv.eigenvalues = make_tensor({2}, evals); + rv.matrix = make_tensor({2, 2}, mat); + rv.eigenvectors = make_tensor({2, 2}, evecs); + return rv; +} + +// Checks that the eigenvalues are in the same order and approximately equal. +inline void require_eigenvalues_approx(const simde::type::tensor& values, + const simde::type::tensor& expected, + double rtol = 1e-6) { + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(values, expected, rtol)); +} + +// Ensures residual H = C^T A C - diag(L) is small. +inline void require_eigenpair_residual(const simde::type::tensor& A, + const simde::type::tensor& L, + const simde::type::tensor& C, + double rtol = 1e-6) { + using tensorwrapper::operations::approximately_equal; + using tensorwrapper::utilities::diagonal_matrix; + auto L_diag = diagonal_matrix(L); + + tensorwrapper::Tensor VA, VAV; + VA("i,k") = C("j,i") * A("j,k"); + VAV("i,k") = VA("i,j") * C("j,k"); + REQUIRE(approximately_equal(VAV, L_diag, rtol)); +} + +#ifdef ENABLE_SIGMA +template +inline void require_uq_eigenvalues_contain(const simde::type::tensor& values, + const simde::type::tensor& expected, + double noise_level) { + using tensorwrapper::buffer::make_contiguous; + using wtf::fp::float_cast; + + auto values_data = get_raw_data(values.buffer()); + auto expected_data = get_raw_data(expected.buffer()); + REQUIRE(values_data.size() == expected_data.size()); + + for(std::size_t i = 0; i < values_data.size(); ++i) { + REQUIRE(values_data[i].contains(expected_data[i].median())); + // Ensure width did not grow too much + REQUIRE(values_data[i].width() <= 50 * noise_level); + } +} + +template +inline void require_uq_eigenpair_residual(const simde::type::tensor& A, + const simde::type::tensor& L, + const simde::type::tensor& C, + double noise_level) { + using tensorwrapper::operations::approximately_equal; + using tensorwrapper::utilities::diagonal_matrix; + auto L_diag = diagonal_matrix(L); + + tensorwrapper::Tensor VA, VAV, H; + VA("i,k") = C("j,i") * A("j,k"); + VAV("i,k") = VA("i,j") * C("j,k"); + H("i,j") = VAV("i,j") - L_diag("i,j"); + + auto H_data = get_raw_data(H.buffer()); + for(std::size_t i = 0; i < H_data.size(); ++i) { + REQUIRE(H_data[i].contains(0.0)); + REQUIRE(H_data[i].width() <= 50 * noise_level); + } +} + +#endif + +} // namespace test_eigen_solver