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

Author: Thoemi09

c++/cppdlr/dlr_dyson.hpp

@@ -20,7 +20,7 @@
 #include <cppdlr/dlr_imtime.hpp>
 #include <cppdlr/dlr_kernels.hpp>
 
-#include <nda/linalg/eigenelements.hpp>
+#include <nda/linalg/eigh.hpp>
 #include <type_traits>
 
 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::get_value_t<Tg>>(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<Ht>) {                                             // 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<nda::get_value_t<Ht>, 3>(r, norb, norb);

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<S>) {                         // 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<S>) { // 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<S>) {                            // 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<S>) { // 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 -------------------
 

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<double>(nda::rand<double>(std::array<int, 2>({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<double>(std::array<int, 1>({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<dcomplex>(nda::rand<double>(std::array<int, 2>({n, n})) + 1i * nda::rand<double>(std::array<int, 2>({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<double>(std::array<int, 1>({n}));    // Random eigenvalues in [-1,1]
   auto h         = matmul(matmul(u, nda::diag(eval)), transpose(conj(u))); // Random Hermitian matrix