caop thrust

apps/caop_selected_basis_diagonalization/main.cc

@@ -87,7 +87,7 @@ int main(int argc, char * argv[]) {
   }
 
   double energy;
-  std::vector<std::vector<size_t>> cobasis;
+  sbd::det_vector<size_t> cobasis;
   sbd::caop::diag<ElemType>(comm,sbd_data,hamfile,basisfiles,
 			    loadname,savename,energy,cobasis);
 

include/sbd/caop/basic/basis.h

@@ -204,18 +204,39 @@ namespace sbd {
 			    size_t bit_length,
 			    size_t total_bit_length) {
     if( get_extension(filename) == std::string("txt") ) {
-      std::ifstream ifs(filename);
-      if( !ifs.is_open() ) {
+      // Fast path: text is regular fixed-width records (total_bit_length chars of
+      // '0'/'1' + '\n').  Read the whole file in one shot and parse in-place.
+      std::ifstream ifs(filename, std::ios::binary);
+      if( !ifs.is_open() )
 	throw std::runtime_error("Failed to open basis bit-string file.");
-      }
-      std::string line;
-      std::vector<std::string> lines;
-      while( std::getline(ifs,line) ) {
-	lines.push_back(line);
-      }
-      config.resize(lines.size());
-      for(size_t i=0; i < lines.size(); i++) {
-	config[i] = from_string(lines[i],bit_length,total_bit_length);
+
+      ifs.seekg(0, std::ios::end);
+      const size_t file_size = ifs.tellg();
+      ifs.seekg(0, std::ios::beg);
+
+      const size_t record_size = total_bit_length + 1;  // L chars + '\n'
+      const size_t num_records = (file_size + 1) / record_size;
+      const size_t num_words   = (total_bit_length + bit_length - 1) / bit_length;
+
+      std::vector<char> buf(file_size);
+      ifs.read(buf.data(), file_size);
+      if( !ifs )
+	throw std::runtime_error("Failed to read basis bit-string file.");
+
+      config.resize(num_records);
+      for(size_t i = 0; i < num_records; i++) {
+	const char* rec = buf.data() + i * record_size;
+	config[i].resize(num_words);
+	std::fill(config[i].begin(), config[i].end(), size_t(0));
+	for(size_t j = 0; j < total_bit_length; j++) {
+	  if( rec[j] == '1' ) {
+	    size_t bit_idx = total_bit_length - 1 - j;
+	    config[i][bit_idx / bit_length] |= (size_t(1) << (bit_idx % bit_length));
+	  } else if( rec[j] != '0' ) {
+	    throw std::runtime_error(
+	      "Unexpected character in basis file: expected '0' or '1'");
+	  }
+	}
       }
     } else if ( get_extension(filename) == std::string("bin") ) {
       std::ifstream ifs(filename, std::ios::binary);

include/sbd/caop/basic/carryover.h

@@ -20,10 +20,10 @@ namespace sbd {
 
   template <typename ElemT, typename RealT>
   void CarryOverBasis(const std::vector<ElemT> & w,
-		      const std::vector<std::vector<size_t>> & bs,
+		      const det_vector<size_t> & bs,
 		      MPI_Comm b_comm,
 		      size_t kept,
-		      std::vector<std::vector<size_t>> & rbs,
+		      det_vector<size_t> & rbs,
 		      RealT & discarted_weight) {
     // using RealT = typename GetRealType<ElemT>::RealT;
     std::vector<RealT> r(w.size());

include/sbd/caop/basic/davidson.h

@@ -14,7 +14,7 @@ namespace sbd {
 
   template <typename ElemT>
   void InitVectorCAOP(std::vector<ElemT> & w,
-		      const std::vector<std::vector<size_t>> & basis,
+		      const det_vector<size_t> & basis,
 		      MPI_Comm h_comm,
 		      MPI_Comm b_comm,
 		      MPI_Comm t_comm,
@@ -38,8 +38,8 @@ namespace sbd {
   template <typename ElemT, typename RealT>
   void Davidson(const std::vector<ElemT> & hii,
 		std::vector<ElemT> & W,
-		const std::vector<std::vector<size_t>> & bs,
-		const size_t bit_length,
+		const det_vector<size_t> & bs,
+		const int bit_length,
 		const std::vector<int> & slide,
 		const GeneralOp<ElemT> & Ham,
 		bool sign,

include/sbd/caop/basic/generalop.h

@@ -80,8 +80,8 @@ namespace sbd {
     friend void mult(const std::vector<ElemT_> & hd,
 		     const std::vector<ElemT_> & wk,
 		     std::vector<ElemT_> & wb,
-		     const std::vector<std::vector<size_t>> & bs,
-		     const size_t bit_length,
+		     const det_vector<size_t> & bs,
+		     const int bit_length,
 		     const std::vector<int> & slide,
 		     const GeneralOp<ElemT_> & H,
 		     bool sign,
@@ -90,14 +90,14 @@ namespace sbd {
 		     MPI_Comm t_comm);
 
     template <typename ElemT_>
-    friend void makeCAOpHamDiagTerms(const std::vector<std::vector<size_t>> & bs,
+    friend void makeCAOpHamDiagTerms(const det_vector<size_t> & bs,
 				     const size_t bit_length,
 				     const std::vector<int> & slide,
 				     const GeneralOp<ElemT_> & H,
 				     std::vector<ElemT_> & hii);
 
     template <typename ElemT_>
-    friend void makeCAOpHam(const std::vector<std::vector<size_t>> & bs,
+    friend void makeCAOpHam(const det_vector<size_t> & bs,
 			    const size_t bit_length,
 			    const std::vector<int> & slide,
 			    const GeneralOp<ElemT_> & H,
@@ -330,8 +330,8 @@ namespace sbd {
     friend void mult(const std::vector<ElemT_> & hd,
 		     const std::vector<ElemT_> & wk,
 		     std::vector<ElemT_> & wb,
-		     const std::vector<std::vector<size_t>> & bs,
-		     const size_t bit_length,
+		     const det_vector<size_t> & bs,
+		     const int bit_length,
 		     const std::vector<int> & slide,
 		     const GeneralOp<ElemT_> & H,
 		     bool sign,
@@ -340,14 +340,14 @@ namespace sbd {
 		     MPI_Comm t_comm);
 
     template <typename ElemT_>
-    friend void makeCAOpHamDiagTerms(const std::vector<std::vector<size_t>> & bs,
+    friend void makeCAOpHamDiagTerms(const det_vector<size_t> & bs,
 				     const size_t bit_length,
 				     const std::vector<int> & slide,
 				     const GeneralOp<ElemT_> & H,
 				     std::vector<ElemT_> & hii);
 
     template <typename ElemT_>
-    friend void makeCAOpHam(const std::vector<std::vector<size_t>> & bs,
+    friend void makeCAOpHam(const det_vector<size_t> & bs,
 			    const size_t bit_length,
 			    const std::vector<int> & slide,
 			    const GeneralOp<ElemT_> & H,
@@ -751,8 +751,8 @@ namespace sbd {
     friend void mult(const std::vector<ElemT_> & hd,
 		     const std::vector<ElemT_> & wk,
 		     std::vector<ElemT_> & wb,
-		     const std::vector<std::vector<size_t>> & bs,
-		     const size_t bit_length,
+		     const det_vector<size_t> & bs,
+		     const int bit_length,
 		     const std::vector<int> & slide,
 		     const GeneralOp<ElemT_> & H,
 		     bool sign,
@@ -761,14 +761,14 @@ namespace sbd {
 		     MPI_Comm t_comm);
 
     template <typename ElemT_>
-    friend void makeCAOpHamDiagTerms(const std::vector<std::vector<size_t>> & bs,
+    friend void makeCAOpHamDiagTerms(const det_vector<size_t> & bs,
 				     const size_t bit_length,
 				     const std::vector<int> & slide,
 				     const GeneralOp<ElemT_> & H,
 				     std::vector<ElemT_> & hii);
 
     template <typename ElemT_>
-    friend void makeCAOpHam(const std::vector<std::vector<size_t>> & bs,
+    friend void makeCAOpHam(const det_vector<size_t> & bs,
 			    const size_t bit_length,
 			    const std::vector<int> & slide,
 			    const GeneralOp<ElemT_> & H,
@@ -804,23 +804,55 @@ namespace sbd {
     friend std::ostream & operator << (std::ostream & s,
 				       const GeneralOp<ElemT_> & gop);
 
+    // Precompute per-term bitmasks for fast rejection in mult (method 0).
+    // m1_[n]: bits that must be occupied in bra (creation targets in o_[n]).
+    // m2_[n]: bits that must be unoccupied in bra (annihilation targets in o_[n]).
+    // Mutable so mult() can populate lazily on first call.
+    void PrecomputeMasks(size_t bit_length) const {
+      size_t det_words = (size_t)(max_index() / (int)bit_length) + 1;
+      std::vector<size_t> zero(det_words, 0);
+      m1_ = det_vector<size_t>(o_.size(), zero);
+      m2_ = det_vector<size_t>(o_.size(), zero);
+      for (size_t n = 0; n < o_.size(); n++) {
+        for (int k = 0; k < o_[n].n_dag_; k++) {
+          size_t q = (size_t)o_[n].fops_[k].q_;
+          m1_[n][q / bit_length] |= (size_t(1) << (q % bit_length));
+        }
+        for (int k = o_[n].n_dag_; k < (int)o_[n].fops_.size(); k++) {
+          size_t q = (size_t)o_[n].fops_[k].q_;
+          m2_[n][q / bit_length] |= (size_t(1) << (q % bit_length));
+        }
+        // bits already in m1_ will be cleared before the undo-annihilation pass
+        // checks them; no bra-level check is needed for those bits in m2_.
+        for (size_t w = 0; w < det_words; w++) {
+          m2_[n][w] &= ~m1_[n][w];
+        }
+      }
+    }
+
   private:
     std::vector<ElemT> e_; // coefficient for diagonal part
     std::vector<ElemT> c_; // coefficient for off-diagonal part
     std::vector<ProductOp> d_; // diagonal part
     std::vector<ProductOp> o_; // off-diagonal part
+    mutable det_vector<size_t> m1_; // per-term occupied-required mask
+    mutable det_vector<size_t> m2_; // per-term unoccupied-required mask
     void copy(const GeneralOp & other) {
       e_ = other.e_;
       c_ = other.c_;
       d_ = other.d_;
       o_ = other.o_;
+      m1_ = other.m1_;
+      m2_ = other.m2_;
     }
 
     void init() {
       e_.resize(0);
       c_.resize(0);
       d_.resize(0);
       o_.resize(0);
+      m1_ = det_vector<size_t>();
+      m2_ = det_vector<size_t>();
     }
 
   }; // end GeneralOp

include/sbd/caop/basic/lanczos.h

@@ -14,8 +14,8 @@ namespace sbd {
   template <typename ElemT, typename RealT>
   void Lanczos(const std::vector<ElemT> & hii,
 	       std::vector<ElemT> & W,
-	       const std::vector<std::vector<size_t>> & bs,
-	       const size_t bit_length,
+	       const det_vector<size_t> & bs,
+	       const int bit_length,
 	       const std::vector<int> & slide,
 	       const GeneralOp<ElemT> & Ham,
 	       bool sign,

include/sbd/caop/basic/mkham.h

@@ -10,7 +10,7 @@
 namespace sbd {
 
   template <typename ElemT>
-  void makeCAOpHamDiagTerms(const std::vector<std::vector<size_t>> & bs,
+  void makeCAOpHamDiagTerms(const det_vector<size_t> & bs,
 			       const size_t bit_length,
 			       const std::vector<int> & slide,
 			       const GeneralOp<ElemT> & H,
@@ -44,7 +44,7 @@ namespace sbd {
   }
 
   template <typename ElemT>
-  void makeCAOpHam(const std::vector<std::vector<size_t>> & bs,
+  void makeCAOpHam(const det_vector<size_t> & bs,
 		   const size_t bit_length,
 		   const std::vector<int> & slide,
 		   const GeneralOp<ElemT> & H,
@@ -71,8 +71,8 @@ namespace sbd {
     size_t num_task = slide.size();
     size_t num_terms = H.NumOpTerms();
 
-    std::vector<std::vector<size_t>> tbs;
-    std::vector<std::vector<size_t>> rbs;
+    det_vector<size_t> tbs;
+    det_vector<size_t> rbs;
 
     if( slide.size() != 0 ) {
       if( slide[0] != 0 ) {
@@ -158,8 +158,8 @@ namespace sbd {
 					 });
 	    */
 	    auto itik = std::lower_bound(tbs.begin(),tbs.end(),vk,
-					 [](const std::vector<size_t> & x,
-					    const std::vector<size_t> & y) {
+					 [](const auto & x,
+					    const auto & y) {
 					   return sbd::less_from_back(x,y);
 					 });
 	    if( itik == tbs.end() ) continue;