Move from std::vector to Eigen::Vector returns for getAccel

Martin D. Weinberg · Martin D. Weinberg · commit 29b6a6fac1df · 2025-08-20T16:55:32.000-04:00
diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H
@@ -85,15 +85,15 @@ namespace BasisClasses
     virtual const std::string harmonic() = 0;
 
     //@{
-    //! Compute acceleration helper
+    //! Compute acceleration helper for each derived class
     virtual void computeAccel(double x, double y, double z,
-			      std::vector<double>& vstore) = 0;
+			      Eigen::Ref<Eigen::Vector3d> vstore) = 0;
 
-    //! C++ local storage for getAccel
-    std::vector<double> vstore {0.0, 0.0, 0.0};
+    //! C++ local storage for computeAccel
+    Eigen::Vector3d vstore;
 
     //! C++ local storage for getAccelArray
-    std::vector<std::vector<double>> vstore_array;
+    RowMatrixXd varray;
     //@}
 
   public:
@@ -184,43 +184,31 @@ namespace BasisClasses
 
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
-    virtual std::vector<double>& getAccel(double x, double y, double z)
+    virtual Eigen::Vector3d& getAccel(double x, double y, double z)
     {
       computeAccel(x, y, z, vstore);
       return vstore;
     }
 
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system for a collections of points
-    virtual RowMatrixXd getAccelArray(Eigen::VectorXd& x,
-				      Eigen::VectorXd& y,
-				      Eigen::VectorXd& z)
+    virtual RowMatrixXd& getAccelArray(Eigen::VectorXd& x,
+				       Eigen::VectorXd& y,
+				       Eigen::VectorXd& z)
     {
       // Sanity check
       if (x.size() != y.size() || x.size() != z.size())
 	throw std::runtime_error("BiorthBasis::getAccelArray: "
 				 "x, y, z vectors must be of the same size");
       
-      int max_threads = 1;
-#if defined(_OPENMP)
-      max_threads = omp_get_max_threads();
-#endif
-      vstore_array.resize(max_threads);
-      for (auto & v : vstore_array) v.resize(3);			
-
-      RowMatrixXd ret(x.size(), 3);
-      std::vector<double> v(3);
-      int tid = 0;		// Default thread ID if no OpenMP
+      // Resize the return array
+      varray.resize(x.size(), 3);
 #pragma omp parallel for
       for (int i=0; i<x.size(); ++i) {
-#if defined(_OPENMP)
-	tid = omp_get_thread_num();
-#endif
-	computeAccel(x(i), y(i), z(i), vstore_array[tid]);
-	for (int k=0; k<3; k++) ret(i, k) = vstore_array[tid][k];
+	computeAccel(x(i), y(i), z(i), varray.row(i));
       }
 
-      return ret;
+      return varray;
     }
   };
 
@@ -315,7 +303,7 @@ namespace BasisClasses
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
     virtual void computeAccel(double x, double y, double z, 
-			  std::vector<double>& vstore);
+			      Eigen::Ref<Eigen::Vector3d> vstore);
 
   public:
     
@@ -558,7 +546,7 @@ namespace BasisClasses
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& vstore);
+		      Eigen::Ref<Eigen::Vector3d> vstore);
 
   public:
     
@@ -714,7 +702,7 @@ namespace BasisClasses
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& vstore);
+		      Eigen::Ref<Eigen::Vector3d> vstore);
 
   public:
     
@@ -865,7 +853,7 @@ namespace BasisClasses
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& vstore);
+		      Eigen::Ref<Eigen::Vector3d> vstore);
 
   public:
     
@@ -1024,7 +1012,7 @@ namespace BasisClasses
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& vstore);
+		      Eigen::Ref<Eigen::Vector3d> vstore);
 
   public:
     
@@ -1159,7 +1147,7 @@ namespace BasisClasses
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system
     void computeAccel(double x, double y, double z, 
-		      std::vector<double>& vstore);
+		      Eigen::Ref<Eigen::Vector3d> vstore);
 
 
   public:
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
@@ -693,7 +693,7 @@ namespace BasisClasses
   }
 
   void Spherical::computeAccel(double x, double y, double z,
-			       std::vector<double>& acc)
+			       Eigen::Ref<Eigen::Vector3d> acc)
   {
     // Get polar coordinates
     double R     = sqrt(x*x + y*y);
@@ -795,7 +795,7 @@ namespace BasisClasses
 
     // Return force not potential gradient
     //
-    acc = {tpotx, tpoty, potz};
+    acc << tpotx, tpoty, potz;
   }
 
 
@@ -1600,7 +1600,7 @@ namespace BasisClasses
   
   // Evaluate in cartesian coordinates
   void Cylindrical::computeAccel(double x, double y, double z,
-				 std::vector<double>& acc)
+				 Eigen::Ref<Eigen::Vector3d> acc)
   {
     double R = sqrt(x*x + y*y);
     double phi = atan2(y, x);
@@ -1614,7 +1614,7 @@ namespace BasisClasses
     double tpotx = tpotR*x/R - tpotp*y/R ;
     double tpoty = tpotR*y/R + tpotp*x/R ;
 
-    acc = {tpotx, tpoty, tpotz};
+    acc << tpotx, tpoty, tpotz;
   }
 
   // Evaluate in cylindrical coordinates
@@ -2182,7 +2182,7 @@ namespace BasisClasses
   }
 
   void FlatDisk::computeAccel(double x, double y, double z, 
-			      std::vector<double>& acc)
+			      Eigen::Ref<Eigen::Vector3d> acc)
   {
     // Get thread id
     int tid = omp_get_thread_num();
@@ -2206,7 +2206,7 @@ namespace BasisClasses
       rpot = -totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
       zpot = -totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
       
-      acc = {rpot, zpot, ppot};
+      acc << rpot, zpot, ppot;
     }
 
     // Get the basis fields
@@ -2270,7 +2270,7 @@ namespace BasisClasses
     double potx = rpot*x/R - ppot*y/R;
     double poty = rpot*y/R + ppot*x/R;
 
-    acc = {potx, poty, zpot};
+    acc << potx, poty, zpot;
   }
 
 
@@ -2958,7 +2958,7 @@ namespace BasisClasses
 
 
   void CBDisk::computeAccel(double x, double y, double z,
-			    std::vector<double>& acc)
+			    Eigen::Ref<Eigen::Vector3d> acc)
   {
     // Get thread id
     int tid = omp_get_thread_num();
@@ -3023,7 +3023,7 @@ namespace BasisClasses
     double potx = rpot*x/R - ppot*y/R;
     double poty = rpot*y/R + ppot*x/R;
 
-    acc = {potx, poty, zpot};
+    acc << potx, poty, zpot;
   }
 
 
@@ -3431,7 +3431,7 @@ namespace BasisClasses
 
 
   void Slab::computeAccel(double x, double y, double z,
-			  std::vector<double>& acc)
+			  Eigen::Ref<Eigen::Vector3d> acc)
   {
     // Loop indices
     //
@@ -3504,7 +3504,7 @@ namespace BasisClasses
       }
     }
 
-    acc = {accx.real(), accy.real(), accz.real()};
+    acc << accx.real(), accy.real(), accz.real();
   }
 
 
@@ -3870,7 +3870,7 @@ namespace BasisClasses
   }
 
   void Cube::computeAccel(double x, double y, double z,
-			  std::vector<double>& acc)
+			  Eigen::Ref<Eigen::Vector3d> acc)
   {
     // Get thread id
     int tid = omp_get_thread_num();
@@ -3881,7 +3881,7 @@ namespace BasisClasses
     // Get the basis fields
     auto frc = ortho->get_force(expcoef, pos);
     
-    acc = {-frc(0).real(), -frc(1).real(), -frc(2).real()};
+    acc << -frc(0).real(), -frc(1).real(), -frc(2).real();
   }
 
   std::vector<double> Cube::cyl_eval(double R, double z, double phi)
diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
@@ -408,7 +408,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE_PURE(void, BiorthBasis, computeAccel, x, y, z, v);
     }
@@ -468,7 +468,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE(void, Spherical, computeAccel, x, y, z, v);
     }
@@ -522,7 +522,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE(void, Cylindrical, computeAccel, x, y, z, v);
     }
@@ -590,7 +590,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE(void, FlatDisk, computeAccel, x, y, z, v);
     }
@@ -662,7 +662,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		  std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE(void, CBDisk, computeAccel, x, y, z, v);
     }
@@ -736,7 +736,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE(void, Slab, computeAccel, x, y, z, v);
     }
@@ -810,7 +810,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void computeAccel(double x, double y, double z,
-		      std::vector<double>& v) override
+		      Eigen::Ref<Eigen::Vector3d> v) override
     {
       PYBIND11_OVERRIDE(void, Cube, computeAccel, x, y, z, v);
     }

Original file line number	Diff line number	Diff line change
`@@ -693,7 +693,7 @@ namespace BasisClasses`
`693`	`693`	`}`
`694`	`694`
`695`	`695`	`void Spherical::computeAccel(double x, double y, double z,`
`696`		`- std::vector<double>& acc)`
	`696`	`+ Eigen::Ref<Eigen::Vector3d> acc)`
`697`	`697`	`{`
`698`	`698`	`// Get polar coordinates`
`699`	`699`	`double R = sqrt(xx + yy);`
`@@ -795,7 +795,7 @@ namespace BasisClasses`
`795`	`795`
`796`	`796`	`// Return force not potential gradient`
`797`	`797`	`//`
`798`		`- acc = {tpotx, tpoty, potz};`
	`798`	`+ acc << tpotx, tpoty, potz;`
`799`	`799`	`}`
`800`	`800`
`801`	`801`
`@@ -1600,7 +1600,7 @@ namespace BasisClasses`
`1600`	`1600`
`1601`	`1601`	`// Evaluate in cartesian coordinates`
`1602`	`1602`	`void Cylindrical::computeAccel(double x, double y, double z,`
`1603`		`- std::vector<double>& acc)`
	`1603`	`+ Eigen::Ref<Eigen::Vector3d> acc)`
`1604`	`1604`	`{`
`1605`	`1605`	`double R = sqrt(xx + yy);`
`1606`	`1606`	`double phi = atan2(y, x);`
`@@ -1614,7 +1614,7 @@ namespace BasisClasses`
`1614`	`1614`	`double tpotx = tpotRx/R - tpotpy/R ;`
`1615`	`1615`	`double tpoty = tpotRy/R + tpotpx/R ;`
`1616`	`1616`
`1617`		`- acc = {tpotx, tpoty, tpotz};`
	`1617`	`+ acc << tpotx, tpoty, tpotz;`
`1618`	`1618`	`}`
`1619`	`1619`
`1620`	`1620`	`// Evaluate in cylindrical coordinates`
`@@ -2182,7 +2182,7 @@ namespace BasisClasses`
`2182`	`2182`	`}`
`2183`	`2183`
`2184`	`2184`	`void FlatDisk::computeAccel(double x, double y, double z,`
`2185`		`- std::vector<double>& acc)`
	`2185`	`+ Eigen::Ref<Eigen::Vector3d> acc)`
`2186`	`2186`	`{`
`2187`	`2187`	`// Get thread id`
`2188`	`2188`	`int tid = omp_get_thread_num();`
`@@ -2206,7 +2206,7 @@ namespace BasisClasses`
`2206`	`2206`	`rpot = -totalMassR/(rr2 + 10.0*std::numeric_limits<double>::min());`
`2207`	`2207`	`zpot = -totalMassz/(rr2 + 10.0*std::numeric_limits<double>::min());`
`2208`	`2208`
`2209`		`- acc = {rpot, zpot, ppot};`
	`2209`	`+ acc << rpot, zpot, ppot;`
`2210`	`2210`	`}`
`2211`	`2211`
`2212`	`2212`	`// Get the basis fields`
`@@ -2270,7 +2270,7 @@ namespace BasisClasses`
`2270`	`2270`	`double potx = rpotx/R - ppoty/R;`
`2271`	`2271`	`double poty = rpoty/R + ppotx/R;`
`2272`	`2272`
`2273`		`- acc = {potx, poty, zpot};`
	`2273`	`+ acc << potx, poty, zpot;`
`2274`	`2274`	`}`
`2275`	`2275`
`2276`	`2276`
`@@ -2958,7 +2958,7 @@ namespace BasisClasses`
`2958`	`2958`
`2959`	`2959`
`2960`	`2960`	`void CBDisk::computeAccel(double x, double y, double z,`
`2961`		`- std::vector<double>& acc)`
	`2961`	`+ Eigen::Ref<Eigen::Vector3d> acc)`
`2962`	`2962`	`{`
`2963`	`2963`	`// Get thread id`
`2964`	`2964`	`int tid = omp_get_thread_num();`
`@@ -3023,7 +3023,7 @@ namespace BasisClasses`
`3023`	`3023`	`double potx = rpotx/R - ppoty/R;`
`3024`	`3024`	`double poty = rpoty/R + ppotx/R;`
`3025`	`3025`
`3026`		`- acc = {potx, poty, zpot};`
	`3026`	`+ acc << potx, poty, zpot;`
`3027`	`3027`	`}`
`3028`	`3028`
`3029`	`3029`
`@@ -3431,7 +3431,7 @@ namespace BasisClasses`
`3431`	`3431`
`3432`	`3432`
`3433`	`3433`	`void Slab::computeAccel(double x, double y, double z,`
`3434`		`- std::vector<double>& acc)`
	`3434`	`+ Eigen::Ref<Eigen::Vector3d> acc)`
`3435`	`3435`	`{`
`3436`	`3436`	`// Loop indices`
`3437`	`3437`	`//`
`@@ -3504,7 +3504,7 @@ namespace BasisClasses`
`3504`	`3504`	`}`
`3505`	`3505`	`}`
`3506`	`3506`
`3507`		`- acc = {accx.real(), accy.real(), accz.real()};`
	`3507`	`+ acc << accx.real(), accy.real(), accz.real();`
`3508`	`3508`	`}`
`3509`	`3509`
`3510`	`3510`
`@@ -3870,7 +3870,7 @@ namespace BasisClasses`
`3870`	`3870`	`}`
`3871`	`3871`
`3872`	`3872`	`void Cube::computeAccel(double x, double y, double z,`
`3873`		`- std::vector<double>& acc)`
	`3873`	`+ Eigen::Ref<Eigen::Vector3d> acc)`
`3874`	`3874`	`{`
`3875`	`3875`	`// Get thread id`
`3876`	`3876`	`int tid = omp_get_thread_num();`
`@@ -3881,7 +3881,7 @@ namespace BasisClasses`
`3881`	`3881`	`// Get the basis fields`
`3882`	`3882`	`auto frc = ortho->get_force(expcoef, pos);`
`3883`	`3883`
`3884`		`- acc = {-frc(0).real(), -frc(1).real(), -frc(2).real()};`
	`3884`	`+ acc << -frc(0).real(), -frc(1).real(), -frc(2).real();`
`3885`	`3885`	`}`
`3886`	`3886`
`3887`	`3887`	`std::vector<double> Cube::cyl_eval(double R, double z, double phi)`
Original file line number	Diff line number	Diff line change
`@@ -408,7 +408,7 @@ void BasisFactoryClasses(py::module &m)`
`408`	`408`	`}`
`409`	`409`
`410`	`410`	`void computeAccel(double x, double y, double z,`
`411`		`- std::vector<double>& v) override`
	`411`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`412`	`412`	`{`
`413`	`413`	`PYBIND11_OVERRIDE_PURE(void, BiorthBasis, computeAccel, x, y, z, v);`
`414`	`414`	`}`
`@@ -468,7 +468,7 @@ void BasisFactoryClasses(py::module &m)`
`468`	`468`	`}`
`469`	`469`
`470`	`470`	`void computeAccel(double x, double y, double z,`
`471`		`- std::vector<double>& v) override`
	`471`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`472`	`472`	`{`
`473`	`473`	`PYBIND11_OVERRIDE(void, Spherical, computeAccel, x, y, z, v);`
`474`	`474`	`}`
`@@ -522,7 +522,7 @@ void BasisFactoryClasses(py::module &m)`
`522`	`522`	`}`
`523`	`523`
`524`	`524`	`void computeAccel(double x, double y, double z,`
`525`		`- std::vector<double>& v) override`
	`525`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`526`	`526`	`{`
`527`	`527`	`PYBIND11_OVERRIDE(void, Cylindrical, computeAccel, x, y, z, v);`
`528`	`528`	`}`
`@@ -590,7 +590,7 @@ void BasisFactoryClasses(py::module &m)`
`590`	`590`	`}`
`591`	`591`
`592`	`592`	`void computeAccel(double x, double y, double z,`
`593`		`- std::vector<double>& v) override`
	`593`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`594`	`594`	`{`
`595`	`595`	`PYBIND11_OVERRIDE(void, FlatDisk, computeAccel, x, y, z, v);`
`596`	`596`	`}`
`@@ -662,7 +662,7 @@ void BasisFactoryClasses(py::module &m)`
`662`	`662`	`}`
`663`	`663`
`664`	`664`	`void computeAccel(double x, double y, double z,`
`665`		`- std::vector<double>& v) override`
	`665`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`666`	`666`	`{`
`667`	`667`	`PYBIND11_OVERRIDE(void, CBDisk, computeAccel, x, y, z, v);`
`668`	`668`	`}`
`@@ -736,7 +736,7 @@ void BasisFactoryClasses(py::module &m)`
`736`	`736`	`}`
`737`	`737`
`738`	`738`	`void computeAccel(double x, double y, double z,`
`739`		`- std::vector<double>& v) override`
	`739`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`740`	`740`	`{`
`741`	`741`	`PYBIND11_OVERRIDE(void, Slab, computeAccel, x, y, z, v);`
`742`	`742`	`}`
`@@ -810,7 +810,7 @@ void BasisFactoryClasses(py::module &m)`
`810`	`810`	`}`
`811`	`811`
`812`	`812`	`void computeAccel(double x, double y, double z,`
`813`		`- std::vector<double>& v) override`
	`813`	`+ Eigen::Ref<Eigen::Vector3d> v) override`
`814`	`814`	`{`
`815`	`815`	`PYBIND11_OVERRIDE(void, Cube, computeAccel, x, y, z, v);`
`816`	`816`	`}`