optimize calculation of covariance matrix.

common/jackknife.h

@@ -152,6 +152,7 @@ errordata<typename Func::res_t > vecjackknife(Func func, const Cont& data, unsig
 	RetType temp = jackmean - jackbin[k];
         fun_error += temp * temp;
     }
+    typedef typename ToScalar<RetType>::RetType ScalarType;
     auto fac = typename ToScalar<RetType>::RetType(data.size() - 1);
     fun_error *= fac/static_cast<double>(data.size());
     fun_error = std::sqrt(fun_error);
@@ -162,17 +163,21 @@ errordata<typename Func::res_t > vecjackknife(Func func, const Cont& data, unsig
 #ifdef _OPENMP
     double start = omp_get_wtime();
 #endif
-    std::valarray<typename ToScalar<RetType>::RetType> cov(numcont*numcont);
+    std::valarray<ScalarType> cov(numcont*numcont);
 #pragma omp parallel for schedule(dynamic)
     for(uint y = 0; y < numcont; ++y)
     {
+      ScalarType fmy = fun_mean[y];
       for(uint x = 0; x < numcont; ++x)
       {
+	ScalarType fmx = fun_mean[x];
+	ScalarType coventry(0);
 	for(uint k = 0; k < data.size(); ++k)
 	{
-	  cov[y * numcont + x] += (jackbin[k][x] - fun_mean[x]) * (jackbin[k][y] - fun_mean[y]);
+	   coventry += (jackbin[k][x] - fmx) * (jackbin[k][y] - fmy);
 	}
-	cov[y * numcont + x] *= (data.size() - 1.0)/data.size();
+	coventry *= (data.size() - 1.0)/data.size();
+	cov[y * numcont + x] = coventry;
       }
     }
 #ifdef _OPENMP