modify the error analysis of big functions to try to calculate as much points...

modify the error analysis of big functions to try to calculate as much points as memory is available.

common/Error_Analysis.h

@@ -33,6 +33,8 @@
 #include <complex>
 #include <valarray>
 #include <algorithm>
+
+#include <unistd.h>
 /*
 //forward declare Observable
 template <class ObsT, class SignType>
@@ -323,23 +325,60 @@ static inline void writeVectorFunctiontoFile(Cont& container, unsigned int len,
     std::vector<Point<FPType> > points;
     points.reserve(4*slices);
     unsigned int nrbins = container.size();//get nr. of bins
-    std::valarray<ScalarType> temp(nrbins);//allocate a valarray for one function    
+    uint availmem = sysconf (_SC_AVPHYS_PAGES) * sysconf(_SC_PAGESIZE);
+    uint memoryperpoint = sizeof(ScalarType)*nrbins;
+    int possibleparallelism = std::min(availmem / memoryperpoint, slices);
+    std::valarray<ScalarType>* temp = NULL;
+    bool memcritical = (possibleparallelism == 0);
+    std::cout<<"Possible parallelism: "<<possibleparallelism<<std::endl;
+    if(!memcritical)
+    {
+      try
+      {
+      temp = new std::valarray<ScalarType>[possibleparallelism];
+      for(int j = 0; j < possibleparallelism; ++j)
+	temp[j].resize(nrbins);
+      }
+      catch(std::bad_alloc& e)
+      {
+	delete [] temp;
+	memcritical = true;
+      }
+    }
+    if(memcritical)
+    {//memory is critical. let's try it anyway and see if the kernel frees memory. 
+      possibleparallelism = 1;
+      temp = new std::valarray<ScalarType>(nrbins);
+    }
+//    std::valarray<ScalarType> temp(nrbins);//allocate a valarray for one function    
     for (unsigned int k = 0; k < len; ++k)//for every index of the function that we stored
     {
-        for (unsigned int j = 0; j < slices; ++j)//for all functionpoints
+        for (unsigned int j = 0; j < slices; j += possibleparallelism)//for all functionpoints
         {
+	  int upperlimit = possibleparallelism;
+	  if(j + possibleparallelism >= slices)
+	    upperlimit = slices - j;//determine the remaining indices
             for (unsigned int l = 0; l < nrbins; ++l)//copy function points
-	        temp[l] = container(l,k)[j];
-            mc_analysis::errordata<ScalarType> edp(handlesign(temp,avsign));//do sign - analysis
+            {
+	      auto tempdatafromfile = container(l, k);
+	      for(int j1 = 0; j1 < upperlimit; ++j1)
+		temp[j1][l] = tempdatafromfile[j+j1];
+//	        temp[l] = container(l, k)[j];
+            }
+            for(int j1 = 0; j1 < upperlimit; ++j1)
+	    {
+            mc_analysis::errordata<ScalarType> edp(handlesign(temp[j1], avsign));//do sign - analysis
             //output data
-            points.push_back(Point<FPType>(j * delta_s, k, real(edp.get_mean()), real(edp.get_error())));
-            fileRe<<j * delta_s<<" "<<real(edp.get_mean())<<" "<<real(edp.get_error())<<std::endl;
-            fileIm<<j * delta_s<<" "<<imag(edp.get_mean())<<" "<<imag(edp.get_error())<<std::endl;
+            points.push_back(Point<FPType>((j+j1) * delta_s, k, real(edp.get_mean()), real(edp.get_error())));
+            fileRe<<(j+j1) * delta_s<<" "<<real(edp.get_mean())<<" "<<real(edp.get_error())<<std::endl;
+            fileIm<<(j+j1) * delta_s<<" "<<imag(edp.get_mean())<<" "<<imag(edp.get_error())<<std::endl;
+	    }
         }
         //output xmgrace separators
         fileRe<<"&"<<std::endl;
         fileIm<<"&"<<std::endl;
     }
+    delete [] temp;
     //now let's generate the transformed data
     sort(points.begin(), points.end());
 #pragma omp parallel for schedule(dynamic)