I am writing a container class for particles to be used in a N-body simulation. In the code below (C++11 flavor), particle data (positions, charge, mass etc.) are stored in a slightly expanded Eigen matrix with an interface that makes it easy to replace an already existing particle container. My question is if such a design is expected to give optimal performance when looping over the particles? One could alternatively split the major matrix into matrices for each property, but I suspect this data partitioning could cause cache misses. Any tips are greatly appreciated!

Code: Select all

#ifdef EIGEN_CORE_H
// Expand Eigen::DenseBase
enum id : Index {X=0,Y,Z,ID,CHARGE,MASS};
inline Scalar charge() const { return (*this)[CHARGE]; }
inline  Scalar& charge() { return (*this)[CHARGE]; }
inline typename FixedSegmentReturnType<3>::Type pos() { return head<3>(); }
inline typename ConstFixedSegmentReturnType<3>::Type pos() const { return head<3>(); }
#else

#define EIGEN_DENSEBASE_PLUGIN "pvec.cpp"
#include <Eigen/Core>
#include <iostream>

template<int DIM=6>
class Particles {
  Eigen::Matrix<float,DIM,Eigen::Dynamic> p_;
  public:
  Particles(int N) {
    p_.resize(DIM,N);
    p_.setRandom(); // fill w. arbitrary data
  }
  int size() const { return p_.cols();}
  inline auto operator[](int i) -> decltype(p_.col(i)) { return p_.col(i); }
};

int main() {
  Particles<> p(10000); // 10000 particles
  float u=0;            // coulomb energy
  for (int i=0; i<p.size()-1; ++i)
    for (int j=i+1; j<p.size(); ++j)
      u += p[i].charge() * p[j].charge()
        / (p[i].pos()-p[j].pos()).norm();
  std::cout << u;
}
#endif


it's ok, but it might be better to store the particles into a rowmajor matrix with accessors to the vector of charges and positions in Particles. This way you could remove one loop, and exploit vectorization:

for (int i=0; i<p.size()-1; ++i)
u += p[i].charge() * (p.charges() .cwiseQuotient( (p.positions().rowwise() - p[i]).rowwise().norm() )).sum();

Thanks, but would this solution not evaluate N^2 pairs while I need only N^2/2 ?

oh yes, you can truncate the computations:

(....).tail(p.size()-i-1).sum();

Aha, I see. So if I get it correctly, you suggest the following matrix layout

Code: Select all

Matrix<float, DIM, N, RowMajor> p;


where each row corresponds to a property (x,y,z,q,...), tightly linked in memory?

The key point is to provide accessors to the set of charges, positions, etc, replace for loops with Eigen expressions, and then you can easily change between a SoA and AoS by switching between row-major and column major to see what's best.