From 8c7764777d434f0e9b31de6d86eaccd76913c6d3 Mon Sep 17 00:00:00 2001 From: Thomas Hahn Date: Wed, 19 Nov 2025 09:53:25 -0500 Subject: [PATCH] Updates to account for changes in nda - getrs requries B to be F-layout - eigenelements was rename to eigh and moved to nda/linalg/eigh.hpp --- c++/cppdlr/dlr_dyson.hpp | 8 ++++---- c++/cppdlr/dlr_imtime.hpp | 22 +++++++++++++--------- test/c++/dyson_it.cpp | 4 ++-- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/c++/cppdlr/dlr_dyson.hpp b/c++/cppdlr/dlr_dyson.hpp index 02cdaf5..20042eb 100644 --- a/c++/cppdlr/dlr_dyson.hpp +++ b/c++/cppdlr/dlr_dyson.hpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include namespace cppdlr { @@ -86,7 +86,7 @@ namespace cppdlr { * \note Hamiltonian must either be a symmetric matrix, a Hermitian matrix, * or a real scalar. */ - dyson_it(double beta, imtime_ops itops, Ht const &h, bool time_order) : dyson_it(beta, itops, h, 0, time_order){}; + dyson_it(double beta, imtime_ops itops, Ht const &h, bool time_order) : dyson_it(beta, itops, h, 0, time_order) {}; /** * @brief Solve Dyson equation for given self-energy @@ -121,7 +121,7 @@ namespace cppdlr { auto g = Tg(sig.shape()); // Declare Green's function g = rhs; // Get right hand side of Dyson equation auto g_rs = nda::matrix_view>(nda::reshape(g, norb, r * norb)); // Reshape g to be compatible w/ LAPACK - nda::lapack::getrs(sysmat, g_rs, ipiv); // Back solve + nda::lapack::getrs(sysmat, transpose(g_rs), ipiv); // Back solve if constexpr (std::floating_point) { // If h is scalar, g is scalar-valued return g; } else { // Otherwise, g is matrix-valued, and need to transpose some indices after solve to undo LAPACK formatting @@ -181,7 +181,7 @@ namespace cppdlr { int norb = h.shape(0); // Diagonalize Hamiltonian - auto [eval, evec] = nda::linalg::eigenelements(h); + auto [eval, evec] = nda::linalg::eigh(h); // Get free Green's function auto g = nda::array, 3>(r, norb, norb); diff --git a/c++/cppdlr/dlr_imtime.hpp b/c++/cppdlr/dlr_imtime.hpp index f8f67eb..1aca5b2 100644 --- a/c++/cppdlr/dlr_imtime.hpp +++ b/c++/cppdlr/dlr_imtime.hpp @@ -580,11 +580,12 @@ namespace cppdlr { fconv += matmul(cf2it, tmp2); // Then precompose with DLR grid values to DLR coefficients matrix - if constexpr (nda::is_complex_v) { // getrs requires matrix and rhs to have same value type - nda::lapack::getrs(transpose(it2cf.zlu), fconv, it2cf.piv); // Note: lapack effectively tranposes fconv by fortran reordering here + if constexpr (nda::is_complex_v) { // getrs requires matrix and rhs to have same value type + nda::lapack::getrs(transpose(it2cf.zlu), transpose(fconv), + it2cf.piv); // Note: lapack effectively tranposes fconv by fortran reordering here } else { // getrs requires matrix and rhs to have same value type - nda::lapack::getrs(transpose(it2cf.lu), fconv, it2cf.piv); + nda::lapack::getrs(transpose(it2cf.lu), transpose(fconv), it2cf.piv); } fconv *= beta; @@ -627,11 +628,12 @@ namespace cppdlr { } // Do the solve - if constexpr (nda::is_complex_v) { // getrs requires matrix and rhs to have same value type - nda::lapack::getrs(transpose(it2cf.zlu), fconvtmp, it2cf.piv); // Note: lapack effectively tranposes fconv by fortran reordering here + if constexpr (nda::is_complex_v) { // getrs requires matrix and rhs to have same value type + nda::lapack::getrs(transpose(it2cf.zlu), transpose(fconvtmp), + it2cf.piv); // Note: lapack effectively tranposes fconv by fortran reordering here } else { // getrs requires matrix and rhs to have same value type - nda::lapack::getrs(transpose(it2cf.lu), fconvtmp, it2cf.piv); + nda::lapack::getrs(transpose(it2cf.lu), transpose(fconvtmp), it2cf.piv); } // Transpose back @@ -903,7 +905,7 @@ namespace cppdlr { } // Precompose with DLR values to coefficients matrix - nda::lapack::getrs(transpose(it2cf.lu), refl, it2cf.piv); // Lapack effectively transposes refl by fortran reordering here + nda::lapack::getrs(transpose(it2cf.lu), transpose(refl), it2cf.piv); // Lapack effectively transposes refl by fortran reordering here } private: @@ -946,7 +948,7 @@ namespace cppdlr { * @param[in] ar Archive to serialize into */ void serialize(auto &ar) const { - ar &lambda_ &r &dlr_rf &dlr_it &cf2it &it2cf.lu &it2cf.zlu &it2cf.piv &hilb &tcf2it &thilb &ttcf2it &ipmat &refl; + ar & lambda_ & r & dlr_rf & dlr_it & cf2it & it2cf.lu & it2cf.zlu & it2cf.piv & hilb & tcf2it & thilb & ttcf2it & ipmat & refl; } /** @@ -956,7 +958,9 @@ namespace cppdlr { * * @param[in] ar Archive to deserialize from */ - void deserialize(auto &ar) { ar &lambda_ &r &dlr_rf &dlr_it &cf2it &it2cf.lu &it2cf.zlu &it2cf.piv &hilb &tcf2it &thilb &ttcf2it &ipmat &refl; } + void deserialize(auto &ar) { + ar & lambda_ & r & dlr_rf & dlr_it & cf2it & it2cf.lu & it2cf.zlu & it2cf.piv & hilb & tcf2it & thilb & ttcf2it & ipmat & refl; + } // -------------------- hdf5 ------------------- diff --git a/test/c++/dyson_it.cpp b/test/c++/dyson_it.cpp index 952a2b5..a3f3329 100644 --- a/test/c++/dyson_it.cpp +++ b/test/c++/dyson_it.cpp @@ -62,7 +62,7 @@ TEST(dyson_it, dyson_vs_ed_real) { // Get random nxn Hamiltonian w/ eigenvalues in [-1 1] auto a = nda::matrix(nda::rand(std::array({n, n}))); // Random matrix a = (a + transpose(a)) / 2; // Make symmetric - auto [eval, u] = nda::linalg::eigenelements(a); // Random orthogonal matrix + auto [eval, u] = nda::linalg::eigh(a); // Random orthogonal matrix eval = -1 + 2 * nda::rand(std::array({n})); // Random eigenvalues in [-1,1] auto h = matmul(matmul(u, nda::diag(eval)), transpose(conj(u))); // Random symmetric matrix @@ -151,7 +151,7 @@ TEST(dyson_it, dyson_vs_ed_cmplx) { // Get random nxn Hamiltonian w/ eigenvalues in [-1 1] auto a = nda::matrix(nda::rand(std::array({n, n})) + 1i * nda::rand(std::array({n, n}))); // Random matrix a = (a + transpose(conj(a))) / 2; // Make Hermitian - auto [eval, u] = nda::linalg::eigenelements(a); // Random unitary matrix + auto [eval, u] = nda::linalg::eigh(a); // Random unitary matrix eval = -1 + 2 * nda::rand(std::array({n})); // Random eigenvalues in [-1,1] auto h = matmul(matmul(u, nda::diag(eval)), transpose(conj(u))); // Random Hermitian matrix