diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index f88f61e0bcb6dd9a71ae5bb552bfb6bf991334d5..844d03711f0600acaf04af08251a6b75de516156 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1677,6 +1677,19 @@ int main(int argc, char *argv[])
   log << "gamma: " << gamma << std::endl; 
   
 
+  
+  
+  
+  ///////////////////////////////////
+  // --- Choose Solver ---
+  // 1 : CG-Solver
+  // 2 : GMRES
+  // 3 : QR 
+  ///////////////////////////////////
+  unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+  
+  
+  
   ///////////////////////////////////
   // Generate the grid
   ///////////////////////////////////
@@ -1993,82 +2006,95 @@ int main(int argc, char *argv[])
   VectorCT x_2 = load_alpha2;
   VectorCT x_3 = load_alpha3;
   
+  
+  auto load_alpha1BS = load_alpha1;
+  
+//   printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
 //   printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
   
   ////////////////////////////////////////////////////////////////////////////////////
+  
+  if (Solvertype == 1)  // CG - SOLVER
+  {
 
-  // CG - SOLVER
   
-//     MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
-// 
-// 
-// 
-//     // Sequential incomplete LU decomposition as the preconditioner
-//     SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,1.0);
-//     int iter = parameterSet.get<double>("iterations_CG", 999);
-//     // Preconditioned conjugate-gradient solver
-//     CGSolver<VectorCT> solver(op,
-//                             ilu0, //NULL,
-//                             1e-8, // desired residual reduction factorlack
-//                             iter,   // maximum number of iterations
-//                             2);   // verbosity of the solver
-//     InverseOperatorResult statistics;
-//     std::cout << "solve linear system for x_1.\n";
-//     solver.apply(x_1, load_alpha1, statistics);
-//     std::cout << "solve linear system for x_2.\n";
-//     solver.apply(x_2, load_alpha2, statistics);
-//     std::cout << "solve linear system for x_3.\n";
-//     solver.apply(x_3, load_alpha3, statistics);
+    std::cout << "------------ CG - Solver ------------" << std::endl;
+    MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
 
+
+
+    // Sequential incomplete LU decomposition as the preconditioner
+    SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,1.0);
+    int iter = parameterSet.get<double>("iterations_CG", 999);
+    // Preconditioned conjugate-gradient solver
+    CGSolver<VectorCT> solver(op,
+                            ilu0, //NULL,
+                            1e-8, // desired residual reduction factorlack
+                            iter,   // maximum number of iterations
+                            2);   // verbosity of the solver
+    InverseOperatorResult statistics;
+    std::cout << "solve linear system for x_1.\n";
+    solver.apply(x_1, load_alpha1, statistics);
+    std::cout << "solve linear system for x_2.\n";
+    solver.apply(x_2, load_alpha2, statistics);
+    std::cout << "solve linear system for x_3.\n";
+    solver.apply(x_3, load_alpha3, statistics);
+  }
   
   
   ////////////////////////////////////////////////////////////////////////////////////
+ 
   
-  // GMRES - SOLVER
-  
+  else if (Solvertype ==2)  // GMRES - SOLVER
+  {
+      
+    std::cout << "------------ GMRES - Solver ------------" << std::endl;
+    // Turn the matrix into a linear operator
+    MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
+
+    // Fancy (but only) way to not have a preconditioner at all
+    Richardson<VectorCT,VectorCT> preconditioner(1.0);
+
+    // Construct the iterative solver
+    RestartedGMResSolver<VectorCT> solver(
+        stiffnessOperator, // Operator to invert
+        preconditioner,    // Preconditioner
+        1e-10,             // Desired residual reduction factor
+        500,               // Number of iterations between restarts,
+                        // here: no restarting
+        500,               // Maximum number of iterations
+        2);                // Verbosity of the solver
+
+    // Object storing some statistics about the solving process
+    InverseOperatorResult statistics;
+
+    // solve for different Correctors (alpha = 1,2,3)
+    solver.apply(x_1, load_alpha1, statistics);
+    solver.apply(x_2, load_alpha2, statistics);
+    solver.apply(x_3, load_alpha3, statistics);
+  }
+  ////////////////////////////////////////////////////////////////////////////////////
   
-  // Turn the matrix into a linear operator
-  MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
-
-  // Fancy (but only) way to not have a preconditioner at all
-  Richardson<VectorCT,VectorCT> preconditioner(1.0);
+  else if ( Solvertype == 3)// QR - SOLVER
+  {
 
-  // Construct the iterative solver
-  RestartedGMResSolver<VectorCT> solver(
-    stiffnessOperator, // Operator to invert
-    preconditioner,    // Preconditioner
-    1e-10,             // Desired residual reduction factor
-    500,               // Number of iterations between restarts,
-                       // here: no restarting
-    500,               // Maximum number of iterations
-    2);                // Verbosity of the solver
+    std::cout << "------------ QR - Solver ------------" << std::endl;
+    log << "solveLinearSystems: We use QR solver.\n";
+    //TODO install suitesparse
+    SPQR<MatrixCT> sPQR(stiffnessMatrix_CE);
+    sPQR.setVerbosity(1);
+    InverseOperatorResult statisticsQR;
 
-  // Object storing some statistics about the solving process
-  InverseOperatorResult statistics;
+    sPQR.apply(x_1, load_alpha1, statisticsQR);
+    sPQR.apply(x_2, load_alpha2, statisticsQR);
+    sPQR.apply(x_3, load_alpha3, statisticsQR);
 
-  // solve for different Correctors (alpha = 1,2,3)
-  solver.apply(x_1, load_alpha1, statistics);
-  solver.apply(x_2, load_alpha2, statistics);
-  solver.apply(x_3, load_alpha3, statistics);
-  
+  }
+//      printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
+//     printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
   
-    
   
-  ////////////////////////////////////////////////////////////////////////////////////
   
-  // QR - SOLVER
-    
-
-//     std::cout << "We use QR solver.\n";
-//     log << "solveLinearSystems: We use QR solver.\n";
-//     //TODO install suitesparse
-//     SPQR<MatrixCT> sPQR(stiffnessMatrix_CE);
-//     sPQR.setVerbosity(1);
-//     InverseOperatorResult statisticsQR;
-// 
-//     sPQR.apply(x_1, load_alpha1, statisticsQR);
-//     sPQR.apply(x_2, load_alpha2, statisticsQR);
-//     sPQR.apply(x_3, load_alpha3, statisticsQR);
 
     std::cout << "---------load_alpha2 AFTER SOLVER: -------------" << std::endl;
 //     printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
@@ -2312,14 +2338,18 @@ int main(int argc, char *argv[])
   coeffT_3[Basis_CE.size()+1] = m_3[1];
   coeffT_3[Basis_CE.size()+2] = m_3[2];
   
+  
+  
+  
+  
   // TEST
 //   for(int i = 0; i<Basis_CE.size()+3 ; i++)
 //   {
-//         if(x_1[i] < 1e-14)
+//         if(x_1[i] < 1e-10)
 //             x_1[i] = 0.0;
-//         if(x_2[i] < 1e-14)
+//         if(x_2[i] < 1e-10)
 //             x_2[i] = 0.0;
-//         if(x_3[i] < 1e-14)
+//         if(x_3[i] < 1e-10)
 //             x_3[i] = 0.0;
 // 
 //   }