bug fix for too large arrays (up from level 9)

mlmc/conf/mlmc.conf

@@ -12,18 +12,18 @@ Problem = StochasticLaplace2D
 StochasticField = LogNormal
 #StochasticField = Gaussian
 
-#Experiment = ConvergenceTest
+Experiment = ConvergenceTest
 #Experiment = MLMCExperiment
-Experiment = MLMCOverEpsilon
+#Experiment = MLMCOverEpsilon
 
-size_init = 3
-l_init = 3 4 5
-M_init = 12 6 3
+size_init = 7
+l_init = 3 4 5 6 7 8 9
+M_init = 100 100 100 100 100 100 100
 Lmax = 9
 
 epsilon = 0.01
-size_lst = 4
-epsilon_lst = 0.01 0.005 0.003 0.001
+size_lst = 3
+epsilon_lst = 0.01 0.005 0.001 #0.005 0.001 0.0005 0.0001
 
 plevel = 2
 functional = L2

mlmc/src/Utils.C

@@ -83,7 +83,8 @@ void fft2(int N2, int N1, fftw_complex *in, fftw_complex *out) {
 void fft2(const vector<vector<complex<double>>> &in, vector<vector<complex<double>>> &out) {
     int N2 = in.size();
     int N1 = in[0].size();
-    fftw_complex in_arr[N2 * N1], out_arr[N2 * N1]; //Todo: Strange crash for N1 = N2 =1023
+    auto in_arr = new fftw_complex[N2 * N1];
+    auto out_arr= new fftw_complex[N2 * N1];
     for (int i = 0; i < N2; i++) {
         for (int j = 0; j < N1; j++) {
             in_arr[i * N2 + j][REAL] = in[i][j].real();
@@ -98,6 +99,8 @@ void fft2(const vector<vector<complex<double>>> &in, vector<vector<complex<doubl
             out[i][j] = {out_arr[i * N2 + j][REAL], out_arr[i * N2 + j][IMAG]};
         }
     }
+    delete[] in_arr;
+    delete[] out_arr;
 }
 
 void fft2(vector<vector<double>> &in, vector<vector<double>> &out) {

mlmc/src/stochfields/CirculantEmbedding.C

@@ -54,7 +54,8 @@ vector<double> CirculantEmbedding1D::generateInterpolField(const vector<double>
     return fieldc;
 }
 
-void CirculantEmbedding1D::createCovarianceMatrix(vector<double> &ar, vector<double> &ac) {
+void
+CirculantEmbedding1D::createCovarianceMatrix(vector<double> &ar, vector<double> &ac) {
     const vector<double> &coordinate1 = linspace(domain_start, domain_end, N);
     ar.resize(N);
     ac.resize(N);
@@ -74,7 +75,8 @@ void CirculantEmbedding1D::createCovarianceMatrix(vector<double> &ar, vector<dou
     }
 }
 
-void CirculantEmbedding1D::embedCovarianceMatrix(vector<double> &ar, vector<double> &ac, vector<double> &br) {
+void CirculantEmbedding1D::embedCovarianceMatrix(vector<double> &ar, vector<double> &ac,
+                                                 vector<double> &br) {
     br.reserve(2 * ar.size() - 1);
     br.insert(br.end(), ar.begin(), ar.end());
     br.insert(br.end(), ac.rbegin(), ac.rend() - 1);
@@ -185,9 +187,11 @@ vector<vector<complex<double>>> CirculantEmbedding2D::generateField() {
     return Xf_;
 }
 
-vector<vector<double>> CirculantEmbedding2D::generateInterpolField(vector<vector<double>> &fieldf) {
+vector<vector<double>>
+CirculantEmbedding2D::generateInterpolField(vector<vector<double>> &fieldf) {
     vector<vector<double>> fieldc;
-    if (N[0] % 2 != 0 || N[1] % 2 != 0) Exit("Interpolation of field does only work for evenly spaced grids.")
+    if (N[0] % 2 != 0 || N[1] % 2 != 0) Exit(
+            "Interpolation of field does only work for evenly spaced grids.")
     int N_interpol[2] = {N[0] / 2, N[1] / 2};
     fieldc.resize(N_interpol[1]);
     for (int i = 0; i < N_interpol[1]; i++) {
@@ -201,7 +205,8 @@ vector<vector<double>> CirculantEmbedding2D::generateInterpolField(vector<vector
     return fieldc;
 }
 
-void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar, vector<vector<double>> &ac) {
+void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar,
+                                                  vector<vector<double>> &ac) {
     const vector<double> &coordinate1 = linspace(domain1_start, domain1_end, N[0]);
     const vector<double> &coordinate2 = linspace(domain2_start, domain2_end, N[1]);
 
@@ -230,7 +235,8 @@ void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar, ve
     }
 }
 
-void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar, vector<vector<double>> &ac,
+void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar,
+                                                 vector<vector<double>> &ac,
                                                  vector<vector<double>> &br) {
     int br_size = 2 * ar.size() - 1;
     br.resize(br_size);
@@ -245,11 +251,13 @@ void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar, vec
     }
 }
 
-void CirculantEmbedding2D::computeEV(vector<vector<double>> &br, vector<vector<double>> &ev) {
+void CirculantEmbedding2D::computeEV(vector<vector<double>> &br,
+                                     vector<vector<double>> &ev) {
     fft2(br, ev);
 }
 
-void CirculantEmbedding2D::padding(vector<vector<double>> &ar, vector<vector<double>> &ac) {
+void
+CirculantEmbedding2D::padding(vector<vector<double>> &ar, vector<vector<double>> &ac) {
     double delta1 = (domain1_end - domain1_start) / (N[0] - 1);
     double delta2 = (domain2_end - domain2_start) / (N[1] - 1);
 

mlmc/src/stochfields/CirculantEmbedding.h

@@ -17,9 +17,6 @@ using namespace std;
  *
  */
 
-
-
-
 class CirculantEmbedding {
 
     /* Circulant Embedding