Use slices when possible rather than element-by-element assignment

Author: Martin D. Weinberg

expui/BiorthBasis.H

@@ -273,6 +273,10 @@ namespace BasisClasses
       return varray;
     }
 
+    //!! Sample counts, masses, coefficients, and covariance
+    SubsampleCovariance::CovarData getCoefCovariance(double time)
+    { return covarStore->getCoefCovariance(time); }
+
     //! Write coefficient covariance data to an HDF5 file
     virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag,
 				     double time=0.0)

exputil/EmpCylSL.cc

@@ -4965,6 +4965,7 @@ EmpCylSL::getCoefCovariance()
     std::get<0>(ret) = Eigen::Tensor<std::complex<double>, 3>(sampT, MMAX+1, NMAX);
     std::get<1>(ret) = Eigen::Tensor<std::complex<double>, 4>(sampT, MMAX+1, NMAX, NMAX);
 
+    /*
     for (unsigned T=0; T<sampT; T++) {
       for (int M=0; M<=MMAX; M++)  {
 	for (int n1=0; n1<NMAX; n1++) {
@@ -4974,6 +4975,15 @@ EmpCylSL::getCoefCovariance()
 	}
       }
     }
+    */
+
+    for (unsigned T=0; T<sampT; T++) {
+      for (int M=0; M<=MMAX; M++)  {
+	std::get<0>(ret).chip(T, 0).chip(M, 0) = Eigen::TensorMap<Eigen::Tensor<std::complex<double>, 1>>(VC[0][T][M].data(), NMAX);
+	std::get<1>(ret).chip(T, 0).chip(M, 0) = Eigen::TensorMap<Eigen::Tensor<std::complex<double>, 2>>(MV[0][T][M].data(), NMAX, NMAX);
+      }
+    }
+
   }
 
   return ret;

src/Cube.cc

@@ -1040,6 +1040,7 @@ PotAccel::CovarData Cube::getSubsample()
 
   // Fill the covariance structure with subsamples
   for (int T=0; T<sampT; T++) {
+    /*
     for (int n1=0; n1<osize; n1++) {
       std::get<2>(elem)(T, 0, n1) = meanV[T](n1);
       for (int n2=0; n2<osize; n2++)
@@ -1048,6 +1049,16 @@ PotAccel::CovarData Cube::getSubsample()
 	else
 	  std::get<3>(elem)(T, 0, n1, n2) = 0.0;
     } 
+    */
+
+    std::get<2>(elem).chip(T, 0).chip(0, 0) =
+      Eigen::TensorMap<Eigen::Tensor<std::complex<double>, 1>>
+      (meanV[T].data(), osize);
+	
+    std::get<3>(elem).chip(T, 0).chip(0, 0) =
+      Eigen::TensorMap<Eigen::Tensor<std::complex<double>, 2>>
+      (covrV[T].data(), osize, osize);
+
   }
 
   return elem;