Added Intel MPI continuous integration workflow

.github/workflows/intelmpi.yml

@@ -0,0 +1,40 @@
+# this method does not use Github Actions cache--good for infrequent simple runs
+jobs:
+
+  linux-intel-oneapi:
+    runs-on: ubuntu-latest
+
+    steps:
+    - name: Intel Apt repository
+      timeout-minutes: 1
+      run: |
+        wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB
+        sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB
+        rm GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB
+        echo "deb https://apt.repos.intel.com/oneapi all main" | sudo tee /etc/apt/sources.list.d/oneAPI.list
+        sudo apt-get update
+        
+    - name: Install Intel oneAPI compilers
+      timeout-minutes: 5
+      run: sudo apt-get install intel-oneapi-compiler-fortran intel-oneapi-compiler-dpcpp-cpp 
+
+    # optional
+    - name: Install Intel MPI and MKL
+      timeout-minutes: 5
+      run: intel-oneapi-mpi intel-oneapi-mpi-devel intel-oneapi-mkl
+
+    - name: Setup Intel oneAPI environment
+      run: |
+        source /opt/intel/oneapi/setvars.sh
+        printenv >> $GITHUB_ENV
+    - name: checkout project code
+      uses: actions/checkout@v4
+
+    - name: CMake Configure
+      run:  |
+        cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=icpx -D BUILD_TESTING=On
+        cmake --build build
+
+    - name: run tests
+      run: |
+        ./build/tests/tests_distributed

README.md

@@ -3,6 +3,7 @@
 [![docs](https://img.shields.io/badge/Documentation-here-red.svg)](https://awietek.github.io/xdiag)
 [![Linux CI](https://github.com/awietek/xdiag/actions/workflows/linux.yml/badge.svg?style=for-the-badge)](https://github.com/awietek/xdiag/actions/workflows/linux.yml)
 [![Mac OSX CI](https://github.com/awietek/xdiag/actions/workflows/osx.yml/badge.svg?style=for-the-badge)](https://github.com/awietek/xdiag/actions/workflows/osx.yml)
+[![[Intel MPI CI](https://github.com/awietek/xdiag/actions/workflows/intelmpi.yml/badge.svg?style=for-the-badge)](https://github.com/awietek/xdiag/actions/workflows/intelmpi.yml)
 [![DOI](https://zenodo.org/badge/169422780.svg)](https://zenodo.org/badge/latestdoi/169422780)
 
 

tests/algorithms/time_evolution/test_time_evolution_distributed.cpp

@@ -15,7 +15,7 @@ TEST_CASE("time_evolution_distributed", "[time_evolution]") try {
   using namespace xdiag;
 
   Log("Test time_evolution_distributed");
-  
+
   int L = 3;
   int n_sites = L * L;
 
@@ -29,14 +29,13 @@ TEST_CASE("time_evolution_distributed", "[time_evolution]") try {
       int site = y * L + x;
       int right = y * L + nx;
       int top = ny * L + x;
-      ops << Op("HOP", "T", {site, right});
-      ops << Op("TJISING", "JZ", {site, right});
-      ops << Op("EXCHANGE", "JEX", {site, right});
-
-      ops << Op("HOP", "T", {site, top});
-      ops << Op("TJISING", "JZ", {site, top});
-      ops << Op("EXCHANGE", "JEX", {site, top});
+      ops += Op("HOP", "T", {site, right});
+      ops += Op("TJISING", "JZ", {site, right});
+      ops += Op("EXCHANGE", "JEX", {site, right});
 
+      ops += Op("HOP", "T", {site, top});
+      ops += Op("TJISING", "JZ", {site, top});
+      ops += Op("EXCHANGE", "JEX", {site, top});
     }
   }
   ops["T"] = 1.0 + 0.2i;
@@ -63,15 +62,14 @@ TEST_CASE("time_evolution_distributed", "[time_evolution]") try {
   fill(psi_0, pstate);
   fill(psi_0d, pstate);
 
-
   auto H_psi_0 = State(block, false);
   apply(ops, psi_0, H_psi_0);
   auto H_psi_0d = State(blockd, false);
   apply(ops, psi_0d, H_psi_0d);
-  
+
   for (int s = 0; s < n_sites; ++s) {
-    auto n = innerC(Op("NUMBER", s), H_psi_0);
-    auto nd = innerC(Op("NUMBER", s), H_psi_0d);
+    auto n = innerC(Op("NUMBER", 1.0, s), H_psi_0);
+    auto nd = innerC(Op("NUMBER", 1.0, s), H_psi_0d);
     // Log("i {} {} {}", s, n, nd);
     REQUIRE(close(n, nd));
   }
@@ -83,14 +81,14 @@ TEST_CASE("time_evolution_distributed", "[time_evolution]") try {
     auto psi = time_evolve(ops, psi_0, time, tol);
     auto psid = time_evolve(ops, psi_0d, time, tol);
     for (int s = 0; s < n_sites; ++s) {
-      auto n = innerC(Op("NUMBER", s), psi);
-      auto nd = innerC(Op("NUMBER", s), psid);
+      auto n = innerC(Op("NUMBER", 1.0, s), psi);
+      auto nd = innerC(Op("NUMBER", 1.0, s), psid);
       Log("{} {} {} {:.6f}", s, n, nd, time);
       REQUIRE(std::abs(n - nd) < 1e-6);
     }
     Log("\n");
   }
 
- } catch (xdiag::Error e) {
+} catch (xdiag::Error e) {
   error_trace(e);
 }

tests/blocks/tj_distributed/test_tj_distributed_apply.cpp

@@ -14,8 +14,7 @@
 
 using namespace xdiag;
 
-void test_tjdistributed_e0_real(OpSum ops, int nup, int ndn, double e0) {
-  int n_sites = ops.n_sites();
+void test_tjdistributed_e0_real(OpSum ops, int n_sites, int nup, int ndn, double e0) {
   auto block = tJDistributed(n_sites, nup, ndn);
   double e0c = eigval0(ops, block);
   REQUIRE(std::abs(e0 - e0c) < 1e-6);
@@ -33,10 +32,10 @@ TEST_CASE("tj_distributed_apply", "[tj_distributed]") try {
       auto block = tJDistributed(N, nup, ndn);
       OpSum ops;
       for (int i = 0; i < N; ++i) {
-        ops << Op("ISING", "Jz", {i, (i + 1) % N});
-        ops << Op("EXCHANGE", "Jx", {i, (i + 1) % N});
-        ops << Op("HOPDN", "TDN", {i, (i + 1) % N});
-        ops << Op("HOPUP", "TUP", {i, (i + 1) % N});
+        ops += Op("ISING", "Jz", {i, (i + 1) % N});
+        ops += Op("EXCHANGE", "Jx", {i, (i + 1) % N});
+        ops += Op("HOPDN", "TDN", {i, (i + 1) % N});
+        ops += Op("HOPUP", "TUP", {i, (i + 1) % N});
       }
       ops["Jz"] = 1.32;
       ops["Jx"] = complex(.432, .576);
@@ -82,7 +81,7 @@ TEST_CASE("tj_distributed_apply", "[tj_distributed]") try {
         {4, 2, -2.11803398}, {5, 0, -0.99999999}, {5, 1, -0.49999999},
         {6, 0, 1.500000000}};
     for (auto [nup, ndn, e0] : nup_ndn_e0) {
-      test_tjdistributed_e0_real(ops, nup, ndn, e0);
+      test_tjdistributed_e0_real(ops, 6, nup, ndn, e0);
     }
   }
 
@@ -101,7 +100,7 @@ TEST_CASE("tj_distributed_apply", "[tj_distributed]") try {
         {4, 2, 0.000000000}, {5, 0, -2.00000000}, {5, 1, 0.000000000},
         {6, 0, 0.000000000}};
     for (auto [nup, ndn, e0] : nup_ndn_e0) {
-      test_tjdistributed_e0_real(ops, nup, ndn, e0);
+      test_tjdistributed_e0_real(ops, 6, nup, ndn, e0);
     }
   }
 

xdiag/basis/spinhalf_distributed/basis_sz.cpp

@@ -117,7 +117,7 @@ BasisSz<bit_t>::BasisSz(int64_t n_sites, int64_t n_up)
   // Create the transpose communicator
   std::vector<int64_t> n_states_i_send(mpi_size_, 0);
   for (bit_t prefix : prefixes()) {
-    for (bit_t postfix : postfixes(prefix)) {
+    for (bit_t postfix : postfix_states(prefix)) {
       int target_rank = rank(postfix);
       ++n_states_i_send[target_rank];
     }
@@ -127,9 +127,9 @@ BasisSz<bit_t>::BasisSz(int64_t n_sites, int64_t n_up)
   // Create the transpose communicator (reverse)
   std::vector<int64_t> n_states_i_send_reverse(mpi_size_, 0);
   for (bit_t postfix : postfixes()) {
-    for (bit_t prefix : prefixes(postfix)) {
+    for (bit_t prefix : prefix_states(postfix)) {
       int target_rank = rank(prefix);
-        ++n_states_i_send_reverse[target_rank;
+      ++n_states_i_send_reverse[target_rank];
     }
   }
   transpose_communicator_reverse_ = mpi::Communicator(n_states_i_send_reverse);
@@ -204,7 +204,11 @@ std::vector<bit_t> const &BasisSz<bit_t>::prefix_states(bit_t postfix) const {
   return prefix_states_[n_up_prefix];
 }
 
-mpi::CommPattern &BasisSz<bit_t>::comm_pattern() { return comm_pattern_; }
+template <typename bit_t> mpi::CommPattern &BasisSz<bit_t>::comm_pattern() {
+  return comm_pattern_;
+}
+
+template <typename bit_t>
 mpi::Communicator BasisSz<bit_t>::transpose_communicator(bool reverse) const {
   return reverse ? transpose_communicator_reverse_ : transpose_communicator_;
 }

xdiag/basis/spinhalf_distributed/basis_sz.hpp

@@ -6,6 +6,7 @@
 #include <xdiag/bits/bitops.hpp>
 #include <xdiag/combinatorics/lin_table.hpp>
 #include <xdiag/common.hpp>
+#include <xdiag/parallel/mpi/comm_pattern.hpp>
 #include <xdiag/parallel/mpi/communicator.hpp>
 #include <xdiag/random/hash_functions.hpp>
 

xdiag/blocks/tj_distributed/terms/apply_exchange.hpp

@@ -14,6 +14,8 @@ namespace xdiag::tj_distributed {
 template <typename bit_t, typename coeff_t, class Basis>
 void apply_exchange(Op const &op, Basis &&basis, const coeff_t *vec_in,
                     coeff_t *vec_out) {
+  using namespace bits;
+
   assert(op.coupling_defined());
   assert(op.type_defined());
   assert(op.size() == 2);
@@ -24,13 +26,15 @@ void apply_exchange(Op const &op, Basis &&basis, const coeff_t *vec_in,
   int mpi_size;
   MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
 
-  using namespace bits;
+  Coupling cpl = op.coupling();
+  assert(cpl.isexplicit() && !cpl.ismatrix());
+  coeff_t J = cpl.as<coeff_t>();
+  coeff_t Jhalf = J / 2.;
+  coeff_t Jhalf_conj = conj(Jhalf);
 
   int64_t s1 = op[0];
   int64_t s2 = op[1];
-  coeff_t J = op.coupling<coeff_t>();
-  coeff_t Jhalf = J / 2.;
-  coeff_t Jhalf_conj = xdiag::conj(Jhalf);
+
   bit_t flipmask = ((bit_t)1 << s1) | ((bit_t)1 << s2);
 
   int64_t l = std::min(s1, s2);

xdiag/blocks/tj_distributed/terms/apply_hopping.hpp

@@ -20,13 +20,16 @@ void apply_hopping(Op const &op, Basis &&basis, const coeff_t *vec_in,
   std::string type = op.type();
   assert((type == "HOPUP") || (type == "HOPDN"));
 
+  Coupling cpl = op.coupling();
+  assert(cpl.isexplicit() && !cpl.ismatrix());
+  coeff_t t = cpl.as<coeff_t>();
+
   int64_t s1 = op[0];
   int64_t s2 = op[1];
   bit_t flipmask = ((bit_t)1 << s1) | ((bit_t)1 << s2);
   int64_t l = std::min(s1, s2);
   int64_t u = std::max(s1, s2);
   bit_t fermimask = (((bit_t)1 << (u - l - 1)) - 1) << (l + 1);
-  coeff_t t = op.coupling<coeff_t>();
 
   auto term_action = [&](bit_t spins) -> std::pair<bit_t, coeff_t> {
     bool fermi = bits::popcnt(spins & fermimask) & 1;

xdiag/blocks/tj_distributed/terms/apply_ising.hpp

@@ -17,7 +17,10 @@ void apply_ising(Op const &op, Basis &&basis, const coeff_t *vec_in,
   std::string type = op.type();
   assert((type == "ISING") || (type == "TJISING"));
 
-  coeff_t J = op.coupling<coeff_t>();
+  Coupling cpl = op.coupling();
+  assert(cpl.isexplicit() && !cpl.ismatrix());
+  coeff_t J = cpl.as<coeff_t>();
+
   int64_t s1 = op[0];
   int64_t s2 = op[1];
   bit_t s1_mask = (bit_t)1 << s1;

xdiag/blocks/tj_distributed/terms/apply_number.hpp

@@ -16,8 +16,11 @@ void apply_number(Op const &op, Basis &&basis, const coeff_t *vec_in,
   std::string type = op.type();
   assert((type == "NUMBERUP") || (type == "NUMBERDN"));
 
-  coeff_t mu = op.coupling<coeff_t>();
-  int64_t s = op.site(0);
+  Coupling cpl = op.coupling();
+  assert(cpl.isexplicit() && !cpl.ismatrix());
+  coeff_t mu = cpl.as<coeff_t>();
+
+  int64_t s = op[0];
   bit_t mask = (bit_t)1 << s;
 
   if (type == "NUMBERUP") {

xdiag/blocks/tj_distributed/terms/apply_raise_lower.hpp

@@ -10,8 +10,8 @@ namespace xdiag::tj_distributed {
 
 template <typename bit_t, typename coeff_t, class BasisIn, class BasisOut,
           class Fill>
-void apply_raise_lower(Op const &op, BasisIn &&basis_in,
-                       BasisOut &&basis_out, Fill &&fill) {
+void apply_raise_lower(Op const &op, BasisIn &&basis_in, BasisOut &&basis_out,
+                       Fill &&fill) {
   assert(op.coupling_defined());
   assert(op.type_defined());
   assert(op.size() == 1);
@@ -23,7 +23,10 @@ void apply_raise_lower(Op const &op, BasisIn &&basis_in,
   int64_t s = op[0];
   bit_t site_mask = (bit_t)1 << s;
   bit_t fermi_mask = site_mask - 1;
-  coeff_t c = op.coupling<coeff_t>();
+
+  Coupling cpl = op.coupling();
+  assert(cpl.isexplicit() && !cpl.ismatrix());
+  coeff_t c = cpl.as<coeff_t>();
 
   // Raising operators
   if ((type == "CDAGUP") || (type == "CDAGDN")) {

xdiag/blocks/tj_distributed/terms/apply_terms.hpp

@@ -10,7 +10,7 @@
 namespace xdiag::tj_distributed {
 
 template <typename bit_t, typename coeff_t, class BasisIn, class BasisOut>
-void apply_terms(OpList const &ops, BasisIn const &basis_in,
+void apply_terms(OpSum const &ops, BasisIn const &basis_in,
                  arma::Col<coeff_t> const &vec_in, BasisOut const &basis_out,
                  arma::Col<coeff_t> &vec_out) try {
   (void)basis_out;

xdiag/blocks/tj_distributed/tj_distributed_apply.cpp

@@ -17,12 +17,12 @@ void apply(OpSum const &ops, tJDistributed const &block_in,
   XDIAG_RETHROW(e);
 }
 
-template void apply<double>(OpSum const &, tJDistributed const &,
-                            arma::Col<double> const &, tJDistributed const &,
-                            arma::Col<double> &, double);
+template void apply(OpSum const &, tJDistributed const &,
+                    arma::Col<double> const &, tJDistributed const &,
+                    arma::Col<double> &, double);
 
-template void apply<complex>(OpSum const &, tJDistributed const &,
-                             arma::Col<complex> const &, tJDistributed const &,
-                             arma::Col<complex> &, double);
+template void apply(OpSum const &, tJDistributed const &,
+                    arma::Col<complex> const &, tJDistributed const &,
+                    arma::Col<complex> &, double);
 
 } // namespace xdiag

xdiag/blocks/tj_distributed/tj_distributed_apply.hpp

@@ -11,7 +11,7 @@ namespace xdiag {
 template <typename coeff_t>
 void apply(OpSum const &ops, tJDistributed const &block_in,
            arma::Col<coeff_t> const &vec_in, tJDistributed const &block_out,
-           arma::Col<coeff_t> &vec_out);
+           arma::Col<coeff_t> &vec_out, double zero_precision=1e-12);
 
 } // namespace xdiag
 #endif

xdiag/parallel/mpi/comm_pattern.hpp

@@ -8,7 +8,7 @@ namespace xdiag::mpi {
 
 class CommPattern {
 public:
-  CommPattern = default();
+  CommPattern() = default;
   bool contains(Op const &op);
   Communicator const &operator[](Op const &op) const;
   Communicator &operator[](Op const &op);

xdiag/utils/error.cpp

@@ -97,7 +97,7 @@ void check_dimension_works_with_blas_int_size(int64_t dim) try {
                 "backend, for example the ILP64 interface of IntelMKL.");
   }
   // Backend 64 bit Blas implementation
-  else if ((sizeof(arma::blas_int) == 8) && (dim > ((int64_t)1 << 63) - 2)) {
+  else if ((sizeof(arma::blas_int) == 8) && (dim > ((int64_t)1 << 62) - 2)) {
     XDIAG_THROW(
         "Trying to create a block whose dimension is too large for the "
         "backend BLAS routines. The block dimension requested is larger than"