Merge pull request #12 from mulaga/main

awietek · Aug 26, 2024 · 65e51f9 · 65e51f9
2 parents d9d3fca + a48ebb9
commit 65e51f9
Show file tree

Hide file tree

Showing 18 changed files with 286 additions and 5 deletions.
diff --git a/docs/documentation/algebra/apply.md b/docs/documentation/algebra/apply.md
@@ -2,7 +2,7 @@
 title: apply
 ---
 
-Applies an operator to a state $\vert w \rangle = O \vert v\rangle$.
+Applies an operator to a state $\vert w \rangle = O \vert v\rangle$ or block of states $\left( \vert w_1 \rangle \dots \vert w_M \rangle \right) = O \left( \vert v_1 \rangle \dots \vert v_M \rangle \right)$. 
 
 **Source** [apply.hpp](https://github.com/awietek/xdiag/blob/main/xdiag/states/apply.hpp)
 

diff --git a/docs/documentation/algorithms/eigs_lanczos.md b/docs/documentation/algorithms/eigs_lanczos.md
@@ -31,6 +31,15 @@ Performs an iterative eigenvalue calculation building eigenvectors using the Lan
                 int64_t random_seed = 42);
 	```
 
+=== "C++ (explicit state initialization)"
+    ```c++
+    eigs_lanczos_result_t 
+	eigs_lanczos(OpSum const &ops, Block const &block,
+				State const &state, int64_t neigenvalues = 1,
+				double precision = 1e-12, int64_t max_iterations = 1000,
+				bool force_complex = false);
+	```
+
 ## Parameters
 
 | Name           | Description                                                                       | Default |

diff --git a/docs/documentation/algorithms/eigvals_lanczos.md b/docs/documentation/algorithms/eigvals_lanczos.md
@@ -31,6 +31,15 @@ Performs an iterative eigenvalue calculation using the Lanczos algorithm.
                 int64_t random_seed = 42);
 	```
 
+=== "C++ (explicit state initialization)"
+    ```c++
+    eigvals_lanczos_result_t
+	eigvals_lanczos(OpSum const &ops, Block const &block, 
+				State const &state, int64_t neigvals = 1,
+	            double precision = 1e-12, int64_t max_iterations = 1000,
+			    bool force_complex = false);
+	```
+
 ## Parameters
 
 | Name           | Description                                                                       | Default |

diff --git a/tests/blocks/electron/test_electron_apply.cpp b/tests/blocks/electron/test_electron_apply.cpp
@@ -56,6 +56,11 @@ TEST_CASE("electron_apply", "[electron]") {
         arma::cx_vec w2(block.size(), arma::fill::zeros);
         apply(ops, block, v, block, w2);
         REQUIRE(close(w1, w2));
+        arma::mat m(block.size(), 5, arma::fill::randn);
+        arma::mat n1 = H * m;
+        arma::mat n2(block.size(), 5, arma::fill::zeros);
+        apply(ops, block, m, block, n2);
+        REQUIRE(close(n1, n2));
 
         // Compute eigenvalues and compare
         arma::vec evals_mat;

diff --git a/tests/blocks/electron_symmetric/test_electron_symmetric_apply.cpp b/tests/blocks/electron_symmetric/test_electron_symmetric_apply.cpp
@@ -46,6 +46,11 @@ void test_electron_symmetric_apply(OpSum ops,
           arma::cx_vec w2(block.size(), arma::fill::randn);
           apply(ops, block, v, block, w2);
           REQUIRE(close(w1, w2));
+          arma::cx_mat m(block.size(), 5, arma::fill::randn);
+          arma::cx_mat n1 = H * m;
+          arma::cx_mat n2(block.size(), 5, arma::fill::zeros);
+          apply(ops, block, m, block, n2);
+          REQUIRE(close(n1, n2));
           // toc("apply");
 
           // Compute eigenvalues and compare

diff --git a/tests/blocks/spinhalf/test_spinhalf_apply.cpp b/tests/blocks/spinhalf/test_spinhalf_apply.cpp
@@ -8,6 +8,7 @@
 #include <xdiag/algebra/matrix.hpp>
 #include <xdiag/algorithms/sparse_diag.hpp>
 #include <xdiag/utils/close.hpp>
+#include <xdiag/utils/timing.hpp>
 
 using namespace xdiag;
 
@@ -20,7 +21,9 @@ void test_apply(int N, OpSum ops) {
     arma::vec v(block.size(), arma::fill::randn);
     arma::vec w1 = H * v;
     arma::vec w2(block.size(), arma::fill::zeros);
+    // tic();
     apply(ops, block, v, block, w2);
+    // toc("1 M-V-M");
 
     arma::vec w3 = H * H * v;
     arma::vec w4(block.size(), arma::fill::zeros);
@@ -38,22 +41,53 @@ void test_apply(int N, OpSum ops) {
   }
 }
 
+void test_apply_mat(int N, OpSum ops) {
+  for (int nup = 1; nup <= N; ++nup) {
+    auto block = Spinhalf(N, nup);
+    auto H = matrix(ops, block, block);
+    REQUIRE(H.is_hermitian(1e-8));
+
+    arma::mat v(block.size(), 5, arma::fill::randn);
+    arma::mat w1 = H * v;
+    arma::mat w2(block.size(), 5, arma::fill::zeros);
+    // tic();
+    apply(ops, block, v, block, w2);
+    // toc("5 column M-M-M");
+
+    arma::mat w3 = H * H * v;
+    arma::mat w4(block.size(), 5, arma::fill::zeros);
+    apply(ops, block, w2, block, w4);
+    REQUIRE(close(w3, w4));
+
+    arma::vec evals_mat;
+    arma::eig_sym(evals_mat, H);
+
+    double e0_mat = evals_mat(0);
+    // double e0_app = eigval0(ops, block);
+    auto [e0_app, ev] = eig0(ops, block);
+    // Log("H: {}, nup: {}, mat: {:.5f} app: {:.5f}", N, nup, e0_mat, e0_app);
+    REQUIRE(close(e0_mat, e0_app));
+  }
+}
+
 TEST_CASE("spinhalf_apply", "[spinhalf]") {
   using namespace xdiag::testcases::spinhalf;
 
-  Log.out("spinhalf_apply: Heisenberg chain apply test, J=1.0, N=2,..,6");
+  Log.out("spinhalf_apply: Heisenberg chain apply test, J=1.0, N=2,..,6, matrix-matrix multiplication");
   for (int N = 2; N <= 6; ++N) {
     auto ops = HBchain(N, 1.0);
     test_apply(N, ops);
+    test_apply_mat(N, ops);
   }
 
-  Log.out("spinhalf_apply: Heisenberg alltoall apply test, N=2,..,6");
+  Log.out("spinhalf_apply: Heisenberg alltoall apply test, N=2,..,6, matrix-matrix multiplication");
   for (int N = 2; N <= 6; ++N) {
     auto ops = HB_alltoall(N);
     test_apply(N, ops);
+    test_apply_mat(N, ops); 
   }
 
-  Log.out("spinhalf_apply: Heisenberg all-to-all Sz <-> NoSz comparison");
+  Log.out("spinhalf_apply: Heisenberg all-to-all Sz <-> NoSz comparison, matrix-matrix multiplication");
   for (int n_sites = 2; n_sites <= 6; ++n_sites) {
     auto ops = HB_alltoall(n_sites);
     auto block_no_sz = Spinhalf(n_sites);

diff --git a/tests/blocks/spinhalf_symmetric/test_spinhalf_symmetric_apply.cpp b/tests/blocks/spinhalf_symmetric/test_spinhalf_symmetric_apply.cpp
@@ -30,6 +30,11 @@ void test_spinhalf_symmetric_apply(OpSum ops, PermutationGroup space_group,
         arma::cx_vec w2(block.size(), arma::fill::zeros);
         apply(ops, block, v, block, w2);
         REQUIRE(close(w1, w2));
+        arma::cx_mat m(block.size(), 5, arma::fill::randn);
+        arma::cx_mat n1 = H * m;
+        arma::cx_mat n2(block.size(), 5, arma::fill::zeros);
+        apply(ops, block, m, block, n2);
+        REQUIRE(close(n1, n2));
 
         // Compute eigenvalues and compare
         arma::vec evals_mat;
@@ -74,6 +79,11 @@ void test_spinhalf_symmetric_apply_no_sz(OpSum ops,
       arma::cx_vec w2(block.size(), arma::fill::zeros);
       apply(ops, block, v, block, w2);
       REQUIRE(close(w1, w2));
+      arma::cx_mat m(block.size(), 5, arma::fill::randn);
+      arma::cx_mat n1 = H * m;
+      arma::cx_mat n2(block.size(), 5, arma::fill::zeros);
+      apply(ops, block, m, block, n2);
+      REQUIRE(close(n1, n2));
       // Compute eigenvalues and compare
       arma::vec evals_mat;
       arma::eig_sym(evals_mat, H);

diff --git a/tests/blocks/tj/test_tj_apply.cpp b/tests/blocks/tj/test_tj_apply.cpp
@@ -29,6 +29,11 @@ TEST_CASE("tj_apply", "[tj]") {
         arma::vec w2(block.size(), arma::fill::zeros);
         apply(ops, block, v, block, w2);
         REQUIRE(close(w1, w2));
+        arma::mat m(block.size(), 5, arma::fill::randn);
+        arma::mat n1 = H * m;
+        arma::mat n2(block.size(), 5, arma::fill::zeros);
+        apply(ops, block, m, block, n2);
+        REQUIRE(close(n1, n2));
 
         arma::vec evals_mat;
         arma::eig_sym(evals_mat, H);

diff --git a/tests/blocks/tj_symmetric/test_tj_symmetric_apply.cpp b/tests/blocks/tj_symmetric/test_tj_symmetric_apply.cpp
@@ -35,6 +35,11 @@ void test_apply_tj_symmetric(OpSum ops, PermutationGroup space_group,
           arma::cx_vec w2(block.size(), arma::fill::zeros);
           apply(ops, block, v, block, w2);
           REQUIRE(close(w1, w2));
+          arma::cx_mat m(block.size(), 5, arma::fill::randn);
+          arma::cx_mat n1 = H_sym * m;
+          arma::cx_mat n2(block.size(), 5, arma::fill::zeros);
+          apply(ops, block, m, block, n2);
+          REQUIRE(close(n1, n2));
 
           // Compute eigenvalues and compare
           arma::vec evals_mat;

diff --git a/xdiag/algebra/apply.cpp b/xdiag/algebra/apply.cpp
@@ -61,9 +61,46 @@ void apply(OpSum const &ops, State const &v, State &w, double precision) try {
         apply(ops, v.block(), vvec, w.block(), wvec, precision);
       }
     }
+  } else if (v.n_cols() == w.n_cols()) {
+    if (ops.isreal()) {
+      if (v.isreal() && w.isreal()) {
+        arma::mat vmat = v.matrix(false);
+        arma::mat wmat = w.matrix(false);
+        apply(ops, v.block(), vmat, w.block(), wmat, precision);
+      } } else if (v.isreal() && !w.isreal()) {
+        auto w2 = State(w.block(), true);
+        arma::mat vmat = v.matrix(false);
+        arma::mat wmat = w2.matrix(false);
+        apply(ops, v.block(), vmat, w.block(), wmat, precision);
+        w = w2;
+      } else if (!v.isreal() && w.isreal()) {
+        w.make_complex();
+        arma::cx_mat vmat = v.matrixC(false);
+        arma::cx_mat wmat = w.matrixC(false);
+        apply(ops, v.block(), vmat, w.block(), wmat, precision);
+      } else if (!v.isreal() && !w.isreal()) {
+        arma::cx_mat vmat = v.matrixC(false);
+        arma::cx_mat wmat = w.matrixC(false);
+        apply(ops, v.block(), vmat, w.block(), wmat, precision);
+      } else {
+      if (v.isreal()) {
+        auto v2 = v;
+        v2.make_complex();
+        w.make_complex();
+        arma::cx_mat vmat = v2.matrixC(false);
+        arma::cx_mat wmat = w.matrixC(false);
+        apply(ops, v.block(), vmat, w.block(), wmat, precision);
+      } else {
+        w.make_complex();
+        arma::cx_vec vvec = v.matrixC(false);
+        arma::cx_vec wvec = w.matrixC(false);
+        apply(ops, v.block(), vvec, w.block(), wvec, precision);
+      } 
+    }
   } else {
     XDIAG_THROW("Applying a OpSum to a state with multiple "
-                "columns not yet implemented");
+                "columns generically not yet implemented "
+                "(are the States of the same size?)");
   }
 } catch (Error const &error) {
   XDIAG_RETHROW(error);
@@ -89,6 +126,25 @@ template void apply<complex>(OpSum const &, Spinhalf const &,
                              arma::Col<complex> const &, Spinhalf const &,
                              arma::Col<complex> &, double);
 
+template <typename coeff_t>
+void apply(OpSum const &ops, Spinhalf const &block_in,
+           arma::Mat<coeff_t> const &mat_in, Spinhalf const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision) try {
+  int64_t n_sites = block_in.n_sites();
+  OpSum opsc = operators::compile_spinhalf(ops, n_sites, precision);
+  mat_out.zeros();
+  basis::spinhalf::dispatch_apply(opsc, block_in, mat_in, block_out, mat_out);
+} catch (Error const &e) {
+  XDIAG_RETHROW(e);
+}
+
+template void apply<double>(OpSum const &, Spinhalf const &,
+                            arma::Mat<double> const &, Spinhalf const &,
+                            arma::Mat<double> &, double);
+template void apply<complex>(OpSum const &, Spinhalf const &,
+                            arma::Mat<complex> const &, Spinhalf const &,
+                            arma::Mat<complex> &, double);
+
 template <typename coeff_t>
 void apply(OpSum const &ops, tJ const &block_in,
            arma::Col<coeff_t> const &vec_in, tJ const &block_out,
@@ -109,6 +165,25 @@ template void apply<complex>(OpSum const &, tJ const &,
                              arma::Col<complex> const &, tJ const &,
                              arma::Col<complex> &, double);
 
+template <typename coeff_t>
+void apply(OpSum const &ops, tJ const &block_in,
+           arma::Mat<coeff_t> const &mat_in, tJ const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision) try {
+  int64_t n_sites = block_in.n_sites();
+  OpSum opsc = operators::compile_tj(ops, n_sites, precision);
+  mat_out.zeros();
+  basis::tj::dispatch_apply(opsc, block_in, mat_in, block_out, mat_out);
+} catch (Error const &e) {
+  XDIAG_RETHROW(e);
+}
+
+template void apply<double>(OpSum const &, tJ const &,
+                            arma::Mat<double> const &, tJ const &,
+                            arma::Mat<double> &, double);
+template void apply<complex>(OpSum const &, tJ const &,
+                            arma::Mat<complex> const &, tJ const &,
+                            arma::Mat<complex> &, double);
+
 template <typename coeff_t>
 void apply(OpSum const &ops, Electron const &block_in,
            arma::Col<coeff_t> const &vec_in, Electron const &block_out,
@@ -129,6 +204,25 @@ template void apply<complex>(OpSum const &, Electron const &,
                              arma::Col<complex> const &, Electron const &,
                              arma::Col<complex> &, double);
 
+template <typename coeff_t>
+void apply(OpSum const &ops, Electron const &block_in,
+           arma::Mat<coeff_t> const &mat_in, Electron const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision) try {
+  int64_t n_sites = block_in.n_sites();
+  OpSum opsc = operators::compile_electron(ops, n_sites, precision);
+  mat_out.zeros();
+  basis::electron::dispatch_apply(opsc, block_in, mat_in, block_out, mat_out);
+} catch (Error const &e) {
+  XDIAG_RETHROW(e);
+}
+
+template void apply<double>(OpSum const &, Electron const &,
+                            arma::Mat<double> const &, Electron const &,
+                            arma::Mat<double> &, double);
+template void apply<complex>(OpSum const &, Electron const &,
+                            arma::Mat<complex> const &, Electron const &,
+                            arma::Mat<complex> &, double);
+
 #ifdef XDIAG_USE_MPI
 
 template <typename coeff_t>
@@ -190,6 +284,18 @@ void apply(OpSum const &ops, Block const &block_in,
 } catch (Error const &error) {
   XDIAG_RETHROW(error);
 }
+template <typename coeff_t>
+void apply(OpSum const &ops, Block const &block_in,
+           arma::Mat<coeff_t> const &mat_in, Block const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision) try {
+  std::visit(
+      [&](auto &&block_in, auto &&block_out) {
+        apply(ops, block_in, mat_in, block_out, mat_out, precision);
+      },
+      block_in, block_out);
+} catch (Error const &error) {
+  XDIAG_RETHROW(error);
+}
 
 template void apply(OpSum const &, Block const &, arma::vec const &,
                     Block const &, arma::vec &, double);

diff --git a/xdiag/algebra/apply.hpp b/xdiag/algebra/apply.hpp
@@ -27,16 +27,28 @@ template <typename coeff_t>
 void apply(OpSum const &op, Spinhalf const &block_in,
            arma::Col<coeff_t> const &vec_in, Spinhalf const &block_out,
            arma::Col<coeff_t> &vec_out, double precision = 1e-12);
+template <typename coeff_t>
+void apply(OpSum const &op, Spinhalf const &block_in,
+           arma::Mat<coeff_t> const &mat_in, Spinhalf const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision = 1e-12);
 
 template <typename coeff_t>
 void apply(OpSum const &op, tJ const &block_in,
            arma::Col<coeff_t> const &vec_in, tJ const &block_out,
            arma::Col<coeff_t> &vec_out, double precision = 1e-12);
+template <typename coeff_t>
+void apply(OpSum const &op, tJ const &block_in,
+           arma::Mat<coeff_t> const &mat_in, tJ const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision = 1e-12);
 
 template <typename coeff_t>
 void apply(OpSum const &op, Electron const &block_in,
            arma::Col<coeff_t> const &vec_in, Electron const &block_out,
            arma::Col<coeff_t> &vec_out, double precision = 1e-12);
+           template <typename coeff_t>
+void apply(OpSum const &op, Electron const &block_in,
+           arma::Mat<coeff_t> const &mat_in, Electron const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision = 1e-12);
 
 #ifdef XDIAG_USE_MPI
 template <typename coeff_t>
@@ -56,4 +68,9 @@ void apply(OpSum const &ops, Block const &block_in,
            arma::Col<coeff_t> const &vec_in, Block const &block_out,
            arma::Col<coeff_t> &vec_out, double precision = 1e-12);
 
+template <typename coeff_t>
+void apply(OpSum const &ops, Block const &block_in,
+           arma::Mat<coeff_t> const &mat_in, Block const &block_out,
+           arma::Mat<coeff_t> &mat_out, double precision = 1e-12);
+
 } // namespace xdiag
diff --git a/xdiag/algebra/fill.hpp b/xdiag/algebra/fill.hpp
@@ -45,5 +45,14 @@ inline void fill_apply(arma::Col<coeff_t> const &vec_in,
                        int64_t idx_out, coeff_t val) {
   fill_apply(vec_in.memptr(), vec_out.memptr(), idx_in, idx_out, val);
 }
+template <typename coeff_t>
+inline void fill_apply(arma::Mat<coeff_t> const &mat_in,
+                       arma::Mat<coeff_t> &mat_out, int64_t idx_in,
+                       int64_t idx_out, coeff_t val) {
+  // for each column call the usual fill_apply.
+  for (int i=0; i < mat_in.n_cols; i++) {
+    fill_apply(mat_in.colptr(i), mat_out.colptr(i), idx_in, idx_out, val);
+  }
+}
 
 } // namespace xdiag