Safer distance table implementation

src/OhmmsPETE/TinyVector.h

@@ -187,6 +187,14 @@ struct TinyVector
   inline Type_t* end() { return X + D; }
   inline const Type_t* end() const { return X + D; }
 
+  TinyVector operator-() const
+  {
+    TinyVector opposite;
+    for (size_t d = 0; d < D; ++d)
+      opposite[d] = -X[d];
+    return opposite;
+  }
+
   /** Elementwise comparison
    *
    *  not optimized but useful for testing

src/Particle/AsymmetricDistanceTableData.h

@@ -46,9 +46,9 @@ struct AsymmetricDTD : public DTD_BConds<T, D, SC>, public DistanceTableData
   {
     if (n1 != N_sources || n2 != N_targets || nactive != N_walkers)
     {
-      N_sources  = n1;
+      N_sources = n1;
       N_targets = n2;
-      int m           = n1 * n2;
+      int m     = n1 * n2;
       if (m)
       {
         M.resize(n1 + 1);
@@ -178,16 +178,8 @@ struct AsymmetricDTD : public DTD_BConds<T, D, SC>, public DistanceTableData
     DTD_BConds<T, D, SC>::apply_bc(dr_m, r_m, rinv_m);
   }
 
-  inline void evaluate(ParticleSet& P, int jat)
-  {
-    APP_ABORT("  No need to call AsymmetricDTD::evaluate(ParticleSet& P, int jat)");
-    //based on full evaluation. Only compute it if jat==0
-    if (jat == 0)
-      evaluate(P);
-  }
-
   ///evaluate the temporary pair relations
-  inline void move(const ParticleSet& P, const PosType& rnew)
+  inline void move(const ParticleSet& P, const PosType& rnew, const IndexType iat, bool prepare_old)
   {
     for (int iat = 0; iat < N_sources; iat++)
     {
@@ -200,7 +192,7 @@ struct AsymmetricDTD : public DTD_BConds<T, D, SC>, public DistanceTableData
     }
   }
 
-  inline void update(IndexType jat)
+  inline void update(IndexType jat, bool forward)
   {
     for (int iat = 0, loc = jat; iat < N_sources; iat++, loc += N_targets)
     {

src/Particle/DistanceTableData.h

@@ -70,17 +70,19 @@ enum DistTableType
  * @brief Abstract class to manage pair data between two ParticleSets.
  *
  * Each DistanceTableData object is fined by Source and Target of ParticleSet types.
+ *
  */
 struct DistanceTableData
 {
   static constexpr unsigned DIM = OHMMS_DIM;
 
-  using IndexType       = QMCTraits::IndexType;
-  using RealType        = QMCTraits::RealType;
-  using PosType         = QMCTraits::PosType;
-  using IndexVectorType = aligned_vector<IndexType>;
-  using RowContainer    = VectorSoaContainer<RealType, DIM>;
+  using IndexType    = QMCTraits::IndexType;
+  using RealType     = QMCTraits::RealType;
+  using PosType      = QMCTraits::PosType;
+  using DistRowType  = Vector<RealType, aligned_allocator<RealType>>;
+  using DisplRowType = VectorSoaContainer<RealType, DIM>;
 #ifndef ENABLE_SOA
+  using IndexVectorType = aligned_vector<IndexType>;
   using TempDistType    = TempDisplacement<RealType, DIM>;
   using ripair          = std::pair<RealType, IndexType>;
 #endif
@@ -95,6 +97,9 @@ struct DistanceTableData
   int N_walkers;
 
 #ifndef ENABLE_SOA
+  ///number of pairs
+  int npairs_m;
+
   /** @brief M.size() = N_sources+1
    *
    * M[i+i] - M[i] = the number of connected points to the i-th source
@@ -128,57 +133,67 @@ struct DistanceTableData
   std::vector<TempDistType> Temp;
 #endif
 
+protected:
   /**defgroup SoA data */
   /*@{*/
   /** Distances[i][j] , [N_targets][N_sources]
-   *  Note: For derived AA, only the lower triangle (j<i) is up-to-date after pbyp move
+   *  Note: Derived classes decide if it is a memory view or the actaully storage
+   *        For derived AA, only the lower triangle (j<i) is up-to-date after pbyp move
    *          The upper triangle is symmetric to the lower one only when the full table is evaluated from scratch.
    *          Avoid using the upper triangle because we may change the code to only allocate the lower triangle part.
    *        For derived BA, the full table is up-to-date after pbyp move
    */
-  Matrix<RealType, aligned_allocator<RealType>> Distances;
+  std::vector<DistRowType> Distances;
 
   /** Displacements[N_targets]x[3][N_sources]
-   *  Note: This is a memory view using the memory space allocated in memoryPool
+   *  Note: Derived classes decide if it is a memory view or the actaully storage
    *        Displacements[i][j] = r_A2[j] - r_A1[i], the opposite sign of AoS dr
-   *        For derived AA, A1=A2=A, only the lower triangle (j<i) is allocated in memoryPool
+   *        For derived AA, A1=A2=A, only the lower triangle (j<i) is allocated in memoryPool_displs_
    *          For this reason, Displacements[i] and Displacements[i+1] overlap in memory
    *          and they must be updated in order during PbyP move.
    *        For derived BA, A1=A, A2=B, the full table is allocated.
    */
-  std::vector<RowContainer> Displacements;
-
-  ///actual memory for Displacements
-  aligned_vector<RealType> memoryPool;
+  std::vector<DisplRowType> Displacements;
 
   /** temp_r */
-  aligned_vector<RealType> Temp_r;
+  DistRowType Temp_r;
 
   /** temp_dr */
-  RowContainer Temp_dr;
-
-  /** true, if full table is needed at loadWalker */
-  bool Need_full_table_loadWalker;
+  DisplRowType Temp_dr;
   /*@}*/
 
+  /** whether full table needs to be ready at anytime or not
+   * Optimization can be implemented during forward PbyP move when the full table is not needed all the time.
+   * DT consumers should know if full table is needed or not and request via addTable.
+   */
+  bool need_full_table_;
+
   ///name of the table
   std::string Name;
+
+public:
   ///constructor using source and target ParticleSet
   DistanceTableData(const ParticleSet& source, const ParticleSet& target)
-      : Origin(&source), N_sources(0), N_targets(0), N_walkers(0), Need_full_table_loadWalker(false)
+      : Origin(&source), N_sources(0), N_targets(0), N_walkers(0), need_full_table_(false)
   {}
 
   ///virutal destructor
   virtual ~DistanceTableData() {}
 
+  ///get need_full_table_
+  inline bool getFullTableNeeds() const { return need_full_table_; }
+
+  ///set need_full_table_
+  inline void setFullTableNeeds(bool is_needed) { need_full_table_ = is_needed; }
+
   ///return the name of table
-  inline std::string getName() const { return Name; }
+  inline const std::string& getName() const { return Name; }
+
   ///set the name of table
   inline void setName(const std::string& tname) { Name = tname; }
 
   ///returns the reference the origin particleset
   const ParticleSet& origin() const { return *Origin; }
-  inline void reset(const ParticleSet* newcenter) { Origin = newcenter; }
 
   inline bool is_same_type(int dt_type) const { return DTType == dt_type; }
 
@@ -199,9 +214,9 @@ struct DistanceTableData
   ///returns the number of source particles
   inline IndexType sources() const { return N_sources; }
 
+#ifndef ENABLE_SOA
   inline IndexType getTotNadj() const { return npairs_m; }
 
-#ifndef ENABLE_SOA
   /// return the distance |R[iadj(i,nj)]-R[i]|
   inline RealType distance(int i, int nj) const { return r_m[M[i] + nj]; }
 
@@ -255,17 +270,67 @@ struct DistanceTableData
   }
 #endif
 
-  ///evaluate the full Distance Table
-  virtual void evaluate(ParticleSet& P) = 0;
+  /** return full table distances
+   */
+  const std::vector<DistRowType>& getDistances() const { return Distances; }
 
-  /// evaluate the Distance Table for a given electron
-  virtual void evaluate(ParticleSet& P, int jat) = 0;
+  /** return full table displacements
+   */
+  const std::vector<DisplRowType>& getDisplacements() const { return Displacements; }
 
-  ///evaluate the temporary pair relations
-  virtual void move(const ParticleSet& P, const PosType& rnew) = 0;
+  /** return a row of distances for a given target particle
+   */
+  const DistRowType& getDistRow(int iel) const { return Distances[iel]; }
 
-  ///update the distance table by the pair relations
-  virtual void update(IndexType jat) = 0;
+  /** return a row of displacements for a given target particle
+   */
+  const DisplRowType& getDisplRow(int iel) const { return Displacements[iel]; }
+
+  /** return old distances set up by move() for optimized distance table consumers
+   */
+  virtual const DistRowType& getOldDists() const
+  {
+    APP_ABORT("DistanceTableData::getOldDists is used incorrectly! Contact developers on github.");
+    return Temp_r; // dummy return to avoid compiler warning.
+  }
+
+  /** return old displacements set up by move() for optimized distance table consumers
+   */
+  virtual const DisplRowType& getOldDispls() const
+  {
+    APP_ABORT("DistanceTableData::getOldDispls is used incorrectly! Contact developers on github.");
+    return Temp_dr; // dummy return to avoid compiler warning.
+  }
+
+  /** return the temporary distances when a move is proposed
+   */
+  const DistRowType& getTemporalDists() const { return Temp_r; }
+
+  /** return the temporary displacements when a move is proposed
+   */
+  const DisplRowType& getTemporalDispls() const { return Temp_dr; }
+
+  /** evaluate the full Distance Table
+   * @param P the target particle set
+   */
+  virtual void evaluate(ParticleSet& P) = 0;
+
+  /** evaluate the temporary pair relations when a move is proposed
+   * @param P the target particle set
+   * @param rnew proposed new position
+   * @param iat the particle to be moved
+   * @param prepare_old if true, prepare old distances and displacements for using getOldDists and getOldDispls functions later.
+   *
+   * Note: some distance table consumers (WaveFunctionComponent) have optimized code paths which require prepare_old = true for accepting a move.
+   * Drivers/Hamiltonians know whether moves will be accepted or not and manager this flag when calling ParticleSet::makeMoveXXX functions.
+   */
+  virtual void move(const ParticleSet& P, const PosType& rnew, const IndexType iat = 0, bool prepare_old = true) = 0;
+
+  /** update the distance table by the pair relations if a move is accepted
+   * @param iat the particle with an accepted move
+   * @param forward If true, rows after iat will not be updated. If false, upon accept a move, the full table should be up-to-date
+   */
+  virtual void update(IndexType jat, bool forward = false) = 0;
 
   /** build a compact list of a neighbor for the iat source
    * @param iat source particle id
@@ -352,9 +417,6 @@ struct DistanceTableData
 #endif
   }
 
-  ///number of pairs
-  int npairs_m;
-
 #ifndef ENABLE_SOA
   /**defgroup storage data for nearest-neighbor relations
    */
@@ -395,7 +457,6 @@ struct DistanceTableData
 #endif
     }
   }
-
 };
 } // namespace qmcplusplus
 #endif

src/Particle/HDFWalkerInput_0_4.cpp

@@ -274,7 +274,7 @@ bool HDFWalkerInput_0_4::read_phdf5(std::string h5name)
   h5name.append(hdf::config_ext);
   std::vector<int> woffsets;
   int woffsets_size = 0;
-  bool success = false;
+  bool success      = false;
 
   { // handle small dataset with master rank
     hdf_archive hin(myComm, false);

src/Particle/MCWalkerConfiguration.h

@@ -87,7 +87,7 @@ class MCWalkerConfiguration : public ParticleSet
   std::vector<int> WalkerOffsets;
 
   MCDataType<FullPrecRealType> EnsembleProperty;
-  
+
   // Data for GPU-acceleration via CUDA
   // These hold a list of pointers to the positions, gradients, and
   // laplacians for each walker.  These vectors .data() is often

src/Particle/ParticleSet.BC.cpp

@@ -211,7 +211,7 @@ void ParticleSet::applyBC(const ParticlePos_t& pin, ParticlePos_t& pout, int fir
     else if (pout.getUnit() == PosUnit::Lattice)
       ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Cart2Unit(pin, Lattice.G, pout, first, last);
     else
-       throw std::runtime_error("Unknown unit conversion");
+      throw std::runtime_error("Unknown unit conversion");
   }
   else if (pin.getUnit() == PosUnit::Lattice)
   {
@@ -220,7 +220,7 @@ void ParticleSet::applyBC(const ParticlePos_t& pin, ParticlePos_t& pout, int fir
     else if (pout.getUnit() == PosUnit::Lattice)
       ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Unit2Unit(pin, pout, first, last);
     else
-       throw std::runtime_error("Unknown unit conversion");
+      throw std::runtime_error("Unknown unit conversion");
   }
   else
     throw std::runtime_error("Unknown unit conversion");