diff --git a/src/libintegratorSparse/centralDifferencesSparse.cpp b/src/libintegratorSparse/centralDifferencesSparse.cpp
index 84668559e8b4c3b16983eba1bd00bf5bbfcaa78e..150c4e6ce58f706a9cf275614552a76d778dea2e 100644
--- a/src/libintegratorSparse/centralDifferencesSparse.cpp
+++ b/src/libintegratorSparse/centralDifferencesSparse.cpp
@@ -56,39 +56,11 @@ CentralDifferencesSparse::CentralDifferencesSparse(int numDOFs, double timestep,
   systemMatrix->RemoveRowsColumns(numConstrainedDOFs, constrainedDOFs);
   systemMatrix->BuildSuperMatrixIndices(numConstrainedDOFs, constrainedDOFs, tangentStiffnessMatrix);
 
-//   #ifdef PARDISO
-//     printf("Creating Pardiso solver for central differences.\n");
-//     int positiveDefiniteSolver = 0;
-//     pardisoSolver = new PardisoSolver(systemMatrix, numSolverThreads, positiveDefiniteSolver);
-//   #endif
-//
-//   #ifdef SPOOLES
-//     spoolesSolver = NULL;
-//   #endif
-//
-//   #ifdef PCG
-//     printf("Creating Jacobi solver for central differences.\n");
-//     jacobiPreconditionedCGSolver = new CGSolver(systemMatrix);
-//   #endif
-
   DecomposeSystemMatrix();
 }
 
 CentralDifferencesSparse::~CentralDifferencesSparse()
 {
-//   #ifdef PARDISO
-//     delete(pardisoSolver);
-//   #endif
-//
-//   #ifdef SPOOLES
-//     if (spoolesSolver != NULL)
-//       delete(spoolesSolver);
-//   #endif
-//
-//   #ifdef PCG
-//     delete(jacobiPreconditionedCGSolver);
-//   #endif
-
   delete(systemMatrix);
   delete(tangentStiffnessMatrix);
   delete(rayleighDampingMatrix);
@@ -113,35 +85,6 @@ void CentralDifferencesSparse::DecomposeSystemMatrix()
   systemMatrix->AssignSuperMatrix(tangentStiffnessMatrix);
 
   //systemMatrix->SaveToMatlabFormat("system.mat");
-
-//   if(this->linearSolver->getSolverType() == "PARADISO")
-//   {
-//       int info = this->linearSolver->ComputeCholeskyDecomposition(systemMatrix);
-//       if (info != 0)
-//       {
-//           printf("Error: PARDISO solver returned non-zero exit code %d.\n", info);
-//           exit(1);
-//       }
-//   }
-
-//   #ifdef PARDISO
-//     int info = pardisoSolver->ComputeCholeskyDecomposition(systemMatrix);
-//     if (info != 0)
-//     {
-//       printf("Error: PARDISO solver returned non-zero exit code %d.\n", info);
-//       exit(1);
-//     }
-//   #endif
-//
-//   #ifdef SPOOLES
-//     printf("Creating SPOOLES solver for central differences.\n");
-//     if (spoolesSolver != NULL)
-//       delete(spoolesSolver);
-//     if (numSolverThreads > 1)
-//       spoolesSolver = new SPOOLESSolverMT(systemMatrix, numSolverThreads);
-//     else
-//       spoolesSolver = new SPOOLESSolver(systemMatrix);
-//   #endif
 }
 
 int CentralDifferencesSparse::DoTimestep()
@@ -192,23 +135,6 @@ int CentralDifferencesSparse::DoTimestep()
 
   memset(buffer, 0, sizeof(double) * r);
 
-//   #ifdef SPOOLES
-//     int info = spoolesSolver->SolveLinearSystem(buffer, rhsConstrained);
-//     char solverString[16] = "SPOOLES";
-//   #endif
-//
-//   #ifdef PARDISO
-//     int info = pardisoSolver->SolveLinearSystem(buffer, rhsConstrained);
-//     char solverString[16] = "PARDISO";
-//   #endif
-//
-//   #ifdef PCG
-//     int info = jacobiPreconditionedCGSolver->SolveLinearSystemWithJacobiPreconditioner(buffer, rhsConstrained, 1e-6, 10000);
-//     if (info > 0)
-//       info = 0;
-//     char solverString[16] = "PCG";
-//   #endif
-
   if(!this->linearSolver)
   {
       std::cerr << "Error: The linear solver is not set. Attach a solver before calling this method." << std::endl;
@@ -273,4 +199,8 @@ void CentralDifferencesSparse::ResetToRest()
   IntegratorBaseSparse::ResetToRest();
   timestepIndex = 0;
 }
+SparseMatrix* CentralDifferencesSparse::GetSystemMatrix() const
+{
+  return this->systemMatrix;
+}
 
diff --git a/src/libintegratorSparse/centralDifferencesSparse.h b/src/libintegratorSparse/centralDifferencesSparse.h
index d9ef2c107b9f1b859289a4fed63865ab50f47b2e..b2552403e9be07c999b0dcaf1198de639b197b83 100644
--- a/src/libintegratorSparse/centralDifferencesSparse.h
+++ b/src/libintegratorSparse/centralDifferencesSparse.h
@@ -102,11 +102,7 @@ public:
 
   virtual void ResetToRest();
 
-  SparseMatrix *GetSystemMatrix() const
-  {
-   return this->systemMatrix;
-  }
-
+  SparseMatrix *GetSystemMatrix() const;
 
 protected:
   double * rhs;
diff --git a/src/libintegratorSparse/eulerSparse.cpp b/src/libintegratorSparse/eulerSparse.cpp
index feb1d67e6d6dbe6d9b4b2eb72b44999c355932b0..bf3ec585aa0ae25c83d503d7ee5b400d6e34a424 100644
--- a/src/libintegratorSparse/eulerSparse.cpp
+++ b/src/libintegratorSparse/eulerSparse.cpp
@@ -37,45 +37,12 @@
 
 EulerSparse::EulerSparse(int r, double timestep, SparseMatrix * massMatrix_, ForceModel * forceModel_, int symplectic_, int numConstrainedDOFs_, int * constrainedDOFs_, double dampingMassCoef): IntegratorBaseSparse(r, timestep, massMatrix_, forceModel_, numConstrainedDOFs_, constrainedDOFs_, dampingMassCoef, 0.0), symplectic(symplectic_)
 {
-//   #ifdef PARDISO
-//     printf("Creating Pardiso solver for M.\n");
-//     int positiveDefiniteSolver = 1;
-//     pardisoSolver = new PardisoSolver(massMatrix, 1, positiveDefiniteSolver);
-//     int info = pardisoSolver->ComputeCholeskyDecomposition(massMatrix);
-//     if (info != 0)
-//     {
-//       printf("Error: PARDISO solver returned non-zero exit code %d.\n", info);
-//       exit(1);
-//     }
-//     printf("Solver created.\n");
-//   #endif
-//
-//   #ifdef SPOOLES
-//     printf("Creating SPOOLES solver for M.\n");
-//     spoolesSolver = new SPOOLESSolver(massMatrix);
-//     printf("Solver created.\n");
-//   #endif
-//
-//   #ifdef PCG
-//     printf("Creating Jacobi solver for M.\n");
-//     jacobiPreconditionedCGSolver = new CGSolver(massMatrix);
-//     printf("Solver created.\n");
-//   #endif
+
 }
 
 EulerSparse::~EulerSparse()
 {
-//   #ifdef PARDISO
-//     delete(pardisoSolver);
-//   #endif
-//
-//   #ifdef SPOOLES
-//     delete(spoolesSolver);
-//   #endif
-//
-//   #ifdef PCG
-//     delete(jacobiPreconditionedCGSolver);
-//   #endif
+
 }
 
 // sets the state based on given q, qvel
@@ -135,22 +102,6 @@ int EulerSparse::DoTimestep()
 
   memset(qdelta, 0.0, sizeof(double)*r);
 
-//   #ifdef PARDISO
-//     int info = pardisoSolver->SolveLinearSystem(qdelta, qresidual);
-//     char solverString[16] = "PARDISO";
-//   #endif
-//
-//   #ifdef SPOOLES
-//     int info = spoolesSolver->SolveLinearSystem(qdelta, qresidual);
-//     char solverString[16] = "SPOOLES";
-//   #endif
-//
-//   #ifdef PCG
-//     int info = jacobiPreconditionedCGSolver->SolveLinearSystemWithJacobiPreconditioner(qdelta, qresidual, 1e-6, 10000);
-//     if (info > 0)
-//       info = 0;
-//     char solverString[16] = "PCG";
-//   #endif
   if(!this->linearSolver)
   {
       std::cerr << "Error: The linear solver is not set. Attach a solver before calling this method." << std::endl;
diff --git a/src/libintegratorSparse/eulerSparse.h b/src/libintegratorSparse/eulerSparse.h
index e553f59b9920824fbfa9a3d271dffc2d98d333f3..7c2855151a95ee8623aea03cf1d3dce1506d8f1c 100644
--- a/src/libintegratorSparse/eulerSparse.h
+++ b/src/libintegratorSparse/eulerSparse.h
@@ -75,17 +75,6 @@ public:
 protected:
   int symplectic;
 
-//   #ifdef PARDISO
-//     PardisoSolver * pardisoSolver;
-//   #endif
-//
-//   #ifdef SPOOLES
-//     SPOOLESSolver * spoolesSolver;
-//   #endif
-//
-//   #ifdef PCG
-//     CGSolver * jacobiPreconditionedCGSolver;
-//   #endif
 };
 
 #endif
diff --git a/src/libintegratorSparse/implicitBackwardEulerSparse.cpp b/src/libintegratorSparse/implicitBackwardEulerSparse.cpp
index 8f8698a37166b1b4a339a11f03108b850954c98d..d3fae3fbd2b0fa782200b0595ac6c5aa02b282a5 100644
--- a/src/libintegratorSparse/implicitBackwardEulerSparse.cpp
+++ b/src/libintegratorSparse/implicitBackwardEulerSparse.cpp
@@ -215,37 +215,6 @@ int ImplicitBackwardEulerSparse::DoTimestep()
     PerformanceCounter counterSystemSolveTime;
     memset(buffer, 0, sizeof(double) * r);
 
-//     #ifdef SPOOLES
-//       int info;
-//       if (numSolverThreads > 1)
-//       {
-//         SPOOLESSolverMT * solver = new SPOOLESSolverMT(systemMatrix, numSolverThreads);
-//         info = solver->SolveLinearSystem(buffer, bufferConstrained);
-//         delete(solver);
-//       }
-//       else
-//       {
-//         SPOOLESSolver * solver = new SPOOLESSolver(systemMatrix);
-//         info = solver->SolveLinearSystem(buffer, bufferConstrained);
-//         delete(solver);
-//       }
-//       char solverString[16] = "SPOOLES";
-//     #endif
-//
-//     #ifdef PARDISO
-//       int info = pardisoSolver->ComputeCholeskyDecomposition(systemMatrix);
-//       if (info == 0)
-//         info = pardisoSolver->SolveLinearSystem(buffer, bufferConstrained);
-//       char solverString[16] = "PARDISO";
-//     #endif
-//
-//     #ifdef PCG
-//       int info = jacobiPreconditionedCGSolver->SolveLinearSystemWithJacobiPreconditioner(buffer, bufferConstrained, 1e-6, 10000);
-//       if (info > 0)
-//         info = 0;
-//       char solverString[16] = "PCG";
-//     #endif
-
     if(!this->linearSolver)
     {
         std::cerr << "Error: The linear solver is not set. Attach a solver before calling this method." << std::endl;
diff --git a/src/libintegratorSparse/implicitNewmarkSparse.cpp b/src/libintegratorSparse/implicitNewmarkSparse.cpp
index 07c545c73f1f00e1b26e22ea3a8a360620fa3690..fd5382c0d8e786e3f3a9956156e04e2fa88a6fe7 100644
--- a/src/libintegratorSparse/implicitNewmarkSparse.cpp
+++ b/src/libintegratorSparse/implicitNewmarkSparse.cpp
@@ -72,14 +72,6 @@ ImplicitNewmarkSparse::ImplicitNewmarkSparse(int r, double timestep, SparseMatri
   systemMatrix->RemoveRowsColumns(numConstrainedDOFs, constrainedDOFs);
   systemMatrix->BuildSuperMatrixIndices(numConstrainedDOFs, constrainedDOFs, tangentStiffnessMatrix);
 
-//   #ifdef PARDISO
-//     printf("Creating Pardiso solver. Positive-definite solver: %d. Num threads: %d\n", positiveDefiniteSolver, numSolverThreads);
-//     pardisoSolver = new PardisoSolver(systemMatrix, numSolverThreads, positiveDefiniteSolver);
-//   #endif
-//
-//   #ifdef PCG
-//     jacobiPreconditionedCGSolver = new CGSolver(systemMatrix);
-//   #endif
 }
 
 ImplicitNewmarkSparse::~ImplicitNewmarkSparse()
@@ -146,28 +138,6 @@ int ImplicitNewmarkSparse::SetState(double * q_, double * qvel_)
 
   memset(buffer, 0, sizeof(double) * r);
 
-//   #ifdef SPOOLES
-//     SPOOLESSolver solver(systemMatrix);
-//     int info = solver.SolveLinearSystem(buffer, bufferConstrained);
-//     char solverString[16] = "SPOOLES";
-//   #endif
-//
-//   //massMatrix->Save("M");
-//   //systemMatrix->Save("A");
-//
-//   #ifdef PARDISO
-//     pardisoSolver->ComputeCholeskyDecomposition(systemMatrix);
-//     int info = pardisoSolver->SolveLinearSystem(buffer, bufferConstrained);
-//     char solverString[16] = "PARDISO";
-//   #endif
-//
-//   #ifdef PCG
-//     int info = jacobiPreconditionedCGSolver->SolveLinearSystemWithJacobiPreconditioner(buffer, bufferConstrained, 1e-6, 10000);
-//     if (info > 0)
-//       info = 0;
-//     char solverString[16] = "PCG";
-//   #endif
-
   if(!this->linearSolver)
   {
     std::cerr << "Error: The linear solver is not set. Attach a solver before calling this method." << std::endl;
@@ -314,26 +284,6 @@ int ImplicitNewmarkSparse::DoTimestep()
     PerformanceCounter counterSystemSolveTime;
     memset(buffer, 0, sizeof(double) * r);
 
-//     #ifdef SPOOLES
-//       SPOOLESSolver solver(systemMatrix);
-//       int info = solver.SolveLinearSystem(buffer, bufferConstrained);
-//       char solverString[16] = "SPOOLES";
-//     #endif
-//
-//     #ifdef PARDISO
-//       int info = pardisoSolver->ComputeCholeskyDecomposition(systemMatrix);
-//       if (info == 0)
-//         info = pardisoSolver->SolveLinearSystem(buffer, bufferConstrained);
-//       char solverString[16] = "PARDISO";
-//     #endif
-//
-//     #ifdef PCG
-//       int info = jacobiPreconditionedCGSolver->SolveLinearSystemWithJacobiPreconditioner(buffer, bufferConstrained, 1e-6, 10000);
-//       if (info > 0)
-//         info = 0;
-//       char solverString[16] = "PCG";
-//     #endif
-
     if(!this->linearSolver)
     {
         std::cerr << "Error: The linear solver is not set. Attach a solver before calling this method." << std::endl;
@@ -410,3 +360,8 @@ void ImplicitNewmarkSparse::UseStaticSolver(bool useStaticSolver_)
   }
 }
 
+SparseMatrix* ImplicitNewmarkSparse::GetSystemMatrix() const
+{
+  return this->systemMatrix;
+}
+
diff --git a/src/libintegratorSparse/implicitNewmarkSparse.h b/src/libintegratorSparse/implicitNewmarkSparse.h
index c8a32ab32a580280ad13365fb059f9624edf60d1..fbbd828a635b91f1b24df93e5db86a63ba1d444c 100644
--- a/src/libintegratorSparse/implicitNewmarkSparse.h
+++ b/src/libintegratorSparse/implicitNewmarkSparse.h
@@ -100,10 +100,7 @@ public:
   // dynamic solver is default (i.e. useStaticSolver=false)
   virtual void UseStaticSolver(bool useStaticSolver);
 
-  SparseMatrix *GetSystemMatrix() const
-  {
-   return this->systemMatrix;
-  }
+  SparseMatrix *GetSystemMatrix() const;
 
 protected:
   SparseMatrix * rayleighDampingMatrix;
@@ -123,13 +120,7 @@ protected:
 
   int positiveDefiniteSolver;
   int numSolverThreads;
-//   #ifdef PARDISO
-//     PardisoSolver * pardisoSolver;
-//   #endif
-//
-//   #ifdef PCG
-//     CGSolver * jacobiPreconditionedCGSolver;
-//   #endif
+
 };
 
 #endif
diff --git a/src/libintegratorSparse/integratorBaseSparse.cpp b/src/libintegratorSparse/integratorBaseSparse.cpp
index 216224c0605ec30a5df5e1af76314b3c2365be74..59d7c2c7d6254bfb1598e1b5fdc427ae70bd331e 100644
--- a/src/libintegratorSparse/integratorBaseSparse.cpp
+++ b/src/libintegratorSparse/integratorBaseSparse.cpp
@@ -79,3 +79,8 @@ void IntegratorBaseSparse::setLinearSolver ( LinearSolver* solver )
     linearSolver = solver;
 }
 
+SparseMatrix* IntegratorBaseSparse::GetSystemMatrix() const
+{
+  return this->massMatrix;
+}
+
diff --git a/src/libintegratorSparse/integratorBaseSparse.h b/src/libintegratorSparse/integratorBaseSparse.h
index 848a2482d6381914d189a77910b9548686889b88..b9184bc39f5a1d1f52954a9b056cacc16fba9081 100644
--- a/src/libintegratorSparse/integratorBaseSparse.h
+++ b/src/libintegratorSparse/integratorBaseSparse.h
@@ -68,10 +68,7 @@ public:
 
   virtual double GetKineticEnergy();
   virtual double GetTotalMass();
-  virtual SparseMatrix *GetSystemMatrix() const
-  {
-   return this->massMatrix;
-  }
+  virtual SparseMatrix *GetSystemMatrix() const;
 
   const LinearSolver *getLinearSolver() const;
 
diff --git a/src/libsparseSolver/CGSolver.cpp b/src/libsparseSolver/CGSolver.cpp
index 0ec2976ec604c25ac6ffb6e87d7e8fde46f09198..62ccb8963c7aee1c3d37efb8e7f1edfca9914d11 100644
--- a/src/libsparseSolver/CGSolver.cpp
+++ b/src/libsparseSolver/CGSolver.cpp
@@ -31,7 +31,7 @@
 #include <stdio.h>
 #include "CGSolver.h"
 
-CGSolver::CGSolver(SparseMatrix * A_): A(A_)
+CGSolver::CGSolver(SparseMatrix * A_): A(A_), numIterations(10000), epsilon(1e-6)
 {
   numRows = A->GetNumRows();
   InitBuffers();
@@ -41,7 +41,13 @@ CGSolver::CGSolver(SparseMatrix * A_): A(A_)
   this->solverType = "PCG";
 }
 
-CGSolver::CGSolver(int numRows_, blackBoxProductType callBackFunction_, void * data_, double * diagonal): numRows(numRows_), multiplicator(callBackFunction_), multiplicatorData(data_), A(NULL)
+CGSolver::CGSolver(int numRows_, blackBoxProductType callBackFunction_, void * data_, double * diagonal):
+  numRows(numRows_),
+  numIterations(10000),
+  epsilon(1e-6),
+  multiplicator(callBackFunction_),
+  multiplicatorData(data_),
+  A(NULL)
 {
   InitBuffers();
   invDiagonal = (double*) malloc (sizeof(double) * numRows);
@@ -82,10 +88,10 @@ void CGSolver::InitBuffers()
 // implements the virtual method from LinearSolver by calling "SolveLinearSystem" with default parameters
 int CGSolver::SolveLinearSystem(double * x, const double * b)
 {
-  return SolveLinearSystemWithJacobiPreconditioner(x, b, 1E-6, 10000, 0);
+  return SolveLinearSystemWithJacobiPreconditioner(x, b, 0);
 }
 
-int CGSolver::SolveLinearSystemWithoutPreconditioner(double * x, const double * b, double eps, int maxIterations, int verbose)
+int CGSolver::SolveLinearSystemWithoutPreconditioner(double * x, const double * b, int verbose)
 {
   int iteration=1;
   multiplicator(multiplicatorData, x, r); //A->MultiplyVector(x,r);
@@ -97,8 +103,8 @@ int CGSolver::SolveLinearSystemWithoutPreconditioner(double * x, const double *
 
   double residualNorm2 = ComputeDotProduct(r, r);
   double initialResidualNorm2 = residualNorm2;
-
-  while ((residualNorm2 > eps * eps * initialResidualNorm2) && (iteration <= maxIterations))
+  double eps2 = this->epsilon*this->epsilon;
+  while ((residualNorm2 > eps2 * initialResidualNorm2) && (iteration <= this->numIterations))
   {
     if (verbose)
       printf("CG iteration %d: current L2 error vs initial error=%G\n", iteration, sqrt(residualNorm2 / initialResidualNorm2));
@@ -134,10 +140,10 @@ int CGSolver::SolveLinearSystemWithoutPreconditioner(double * x, const double *
     iteration++;
   }
 
-  return (iteration-1) * ((residualNorm2 > eps * eps * initialResidualNorm2) ? 1 : 0);
+  return (iteration-1) * ((residualNorm2 > eps2 * initialResidualNorm2) ? 1 : 0);
 }
 
-int CGSolver::SolveLinearSystemWithJacobiPreconditioner(double * x, const double * b, double eps, int maxIterations, int verbose)
+int CGSolver::SolveLinearSystemWithJacobiPreconditioner(double * x, const double * b, int verbose)
 {
   if (invDiagonal == NULL)
   {
@@ -164,7 +170,8 @@ int CGSolver::SolveLinearSystemWithJacobiPreconditioner(double * x, const double
   double residualNorm2 = ComputeTriDotProduct(r, r, invDiagonal);
   double initialResidualNorm2 = residualNorm2;
 
-  while ((residualNorm2 > eps * eps * initialResidualNorm2) && (iteration <= maxIterations))
+  double eps2 = this->epsilon*this->epsilon;
+  while ((residualNorm2 > eps2 * initialResidualNorm2) && (iteration <= this->numIterations))
   {
     if (verbose)
       printf("CG iteration %d: current M^{-1}-L2 error vs initial error=%G\n", iteration, sqrt(residualNorm2 / initialResidualNorm2));
@@ -204,7 +211,7 @@ int CGSolver::SolveLinearSystemWithJacobiPreconditioner(double * x, const double
     printf("Warning: residualNorm2=%G is negative. Input matrix might not be SPD. Solution could be incorrect.\n", residualNorm2);
   }
 
-  return (iteration-1) * ((residualNorm2 > eps * eps * initialResidualNorm2) ? 1 : 0);
+  return (iteration-1) * ((residualNorm2 > eps2 * initialResidualNorm2) ? 1 : 0);
 }
 
 double CGSolver::ComputeDotProduct(double * v1, double * v2)
@@ -225,3 +232,23 @@ double CGSolver::ComputeTriDotProduct(double * x, double * y, double * z)
   return result;
 }
 
+void CGSolver::setNumberOfIterations ( const size_t iterations )
+{
+  this->numIterations = iterations;
+}
+
+void CGSolver::setEpsilon ( const double eps )
+{
+  this->epsilon = eps;
+}
+
+size_t CGSolver::getNumberOfIterations() const
+{
+  return this->numIterations;
+}
+
+double CGSolver::getEpsilon() const
+{
+  return this->epsilon;
+}
+
diff --git a/src/libsparseSolver/CGSolver.h b/src/libsparseSolver/CGSolver.h
index 70f4ee6591803109d28c105c5a48f7d1583ac6e2..20f9ff6daefed6ea6cdf90408aae34233508f61c 100755
--- a/src/libsparseSolver/CGSolver.h
+++ b/src/libsparseSolver/CGSolver.h
@@ -76,17 +76,25 @@ public:
   // maximum number of conjugate-gradient iterations is set by "maxIterations"
   // return value is the number of iterations performed
   // if solver did not converge, the return value will have a negative sign
-  int SolveLinearSystemWithoutPreconditioner(double * x, const double * b, double eps=1e-6, int maxIterations=1000, int verbose=0);
+  int SolveLinearSystemWithoutPreconditioner(double * x, const double * b, int verbose=0);
 
   // same as above, except it uses Jacobi preconditioning
   // the employed error metric is M^{-1}-weighted L2 residual error (see Shewchuk)
-  int SolveLinearSystemWithJacobiPreconditioner(double * x, const double * b, double eps=1e-6, int maxIterations=1000, int verbose=0);
+  int SolveLinearSystemWithJacobiPreconditioner(double * x, const double * b,  int verbose=0);
 
   virtual int SolveLinearSystem(double * x, const double * b); // implements the virtual method from LinearSolver by calling "SolveLinearSystemWithJacobiPreconditioner" with default parameters
 
   // computes the dot product of two vectors
   double ComputeDotProduct(double * v1, double * v2); // length of vectors v1, v2 equals numRows (dimension of A)
 
+  void setNumberOfIterations(const size_t iterations);
+
+  void setEpsilon(const double eps);
+
+  size_t getNumberOfIterations() const;
+
+  double getEpsilon() const;
+
 protected:
   int numRows;
   blackBoxProductType multiplicator;
@@ -94,6 +102,8 @@ protected:
   SparseMatrix * A;
   double * r, * d, * q; // terminology from Shewchuk's work
   double * invDiagonal;
+  size_t numIterations;
+  double epsilon;
 
   double ComputeTriDotProduct(double * x, double * y, double * z); // sum_i x[i] * y[i] * z[i]
   static void DefaultMultiplicator(const void * data, const double * x, double * Ax);