diff --git a/cxx/src/scf/eigen_solver/eigenvector_uncertainty.hpp b/cxx/src/scf/eigen_solver/eigenvector_uncertainty.hpp index e443dbc..0522983 100644 --- a/cxx/src/scf/eigen_solver/eigenvector_uncertainty.hpp +++ b/cxx/src/scf/eigen_solver/eigenvector_uncertainty.hpp @@ -61,47 +61,57 @@ struct EigenvectorUncertaintyKernel { using clean_t = std::decay_t; // The dispatch instantiates const-qualified spans too; the writing body // is only valid (and only selected at runtime) for the mutable buffer. - if constexpr(!std::is_const_v && - tensorwrapper::types::is_uq_type_v) { + if constexpr(!std::is_const_v) { using tensorwrapper::buffer::get_raw_data; - using tensorwrapper::types::uq_center; - using tensorwrapper::types::uq_upper; - using value_t = decltype(uq_center(std::declval())); - const auto n = m_n; - auto radius = [](const clean_t& x) { - return uq_upper(x) - uq_center(x); - }; - - auto C = get_raw_data(m_C); - auto G = get_raw_data(m_G); - auto eps = get_raw_data(m_eps); - - // Degeneracy threshold: max abs row-sum of off-diagonal coupling - // radii. - value_t sigma(0); - for(std::size_t i = 0; i < n; ++i) { - value_t row(0); - for(std::size_t j = 0; j < n; ++j) { - if(i != j) { row += radius(G[j * n + i]); } + const auto n = m_n; + auto C = get_raw_data(m_C); + + // out starts as a plain copy of C: for non-UQ types this copy IS + // the whole (correctly no-op) result, since out_buf is otherwise + // a freshly-allocated, default-initialized (i.e. all-zero) + // buffer -- without this copy, plain-float callers would silently + // get back an all-zero eigenvector matrix instead of C unchanged. + for(std::size_t k = 0; k < n * n; ++k) { out[k] = C[k]; } + + if constexpr(tensorwrapper::types::is_uq_type_v) { + using tensorwrapper::types::uq_center; + using tensorwrapper::types::uq_upper; + using value_t = decltype(uq_center(std::declval())); + auto radius = [](const clean_t& x) { + return uq_upper(x) - uq_center(x); + }; + + auto G = get_raw_data(m_G); + auto eps = get_raw_data(m_eps); + + // Degeneracy threshold: max abs row-sum of off-diagonal + // coupling radii. + value_t sigma(0); + for(std::size_t i = 0; i < n; ++i) { + value_t row(0); + for(std::size_t j = 0; j < n; ++j) { + if(i != j) { row += radius(G[j * n + i]); } + } + sigma = std::max(sigma, row); } - sigma = std::max(sigma, row); - } - // Eigenvalue centers. - std::vector eps_c(n); - for(std::size_t i = 0; i < n; ++i) { eps_c[i] = uq_center(eps[i]); } + // Eigenvalue centers. + std::vector eps_c(n); + for(std::size_t i = 0; i < n; ++i) { + eps_c[i] = uq_center(eps[i]); + } - // out = C + correction, accumulated from the ORIGINAL columns so - // one column's correction never feeds into another. - for(std::size_t k = 0; k < n * n; ++k) { out[k] = C[k]; } - for(std::size_t i = 0; i < n; ++i) { - for(std::size_t j = 0; j < n; ++j) { - if(i == j) { continue; } - const value_t gap = eps_c[i] - eps_c[j]; - if(std::abs(gap) <= value_t(2) * sigma) { continue; } - const clean_t coeff = G[j * n + i] / clean_t(gap); - for(std::size_t r = 0; r < n; ++r) { - out[r * n + i] += coeff * C[r * n + j]; + // out += correction, accumulated from the ORIGINAL columns so + // one column's correction never feeds into another. + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + if(i == j) { continue; } + const value_t gap = eps_c[i] - eps_c[j]; + if(std::abs(gap) <= value_t(2) * sigma) { continue; } + const clean_t coeff = G[j * n + i] / clean_t(gap); + for(std::size_t r = 0; r < n; ++r) { + out[r * n + i] += coeff * C[r * n + j]; + } } } } diff --git a/tests/cxx/unit_tests/eigen_solver/eigenvector_uncertainty.cpp b/tests/cxx/unit_tests/eigen_solver/eigenvector_uncertainty.cpp index bb461bf..3d908a1 100644 --- a/tests/cxx/unit_tests/eigen_solver/eigenvector_uncertainty.cpp +++ b/tests/cxx/unit_tests/eigen_solver/eigenvector_uncertainty.cpp @@ -29,7 +29,7 @@ using namespace tensorwrapper::generate; // correction adds real uncertainty without blowing up on a single application). using uq_types = std::tuple; + tensorwrapper::types::tadouble, tensorwrapper::types::tmdouble>; TEMPLATE_LIST_TEST_CASE("attach_eigenvector_uncertainty", "", uq_types) { using tensorwrapper::buffer::get_raw_data; diff --git a/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp b/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp index 4bef32f..595a7a0 100644 --- a/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp +++ b/tests/cxx/unit_tests/eigen_solver/jacobi_normal.cpp @@ -18,7 +18,8 @@ using types = std::tuple; + tensorwrapper::types::adouble, tensorwrapper::types::tadouble, + tensorwrapper::types::tmdouble>; using namespace test_eigen_solver; using namespace tensorwrapper::generate;