diff --git a/bin/demonstrator1.cc b/bin/demonstrator1.cc
index c549752..4ea6053 100644
--- a/bin/demonstrator1.cc
+++ b/bin/demonstrator1.cc
@@ -1,142 +1,143 @@
 /**
  * file   demonstrator1.cc
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   03 Jan 2018
  *
  * @brief  larger problem to show off
  *
  * @section LICENSE
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <iostream>
 #include <memory>
 #include <chrono>
 #include "external/cxxopts.hpp"
 
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "system/system_factory.hh"
 #include "materials/material_hyper_elastic1.hh"
 #include "solver/solvers.hh"
 #include "solver/solver_cg.hh"
 
 
 using opt_ptr = std::unique_ptr<cxxopts::Options>;
 
 
 opt_ptr parse_args(int argc, char **argv) {
   opt_ptr options =
     std::make_unique<cxxopts::Options>(argv[0],  "Tests MPI fft scalability");
 
   try {
     options->add_options()
       ("0,N0", "number of rows", cxxopts::value<int>(), "N0")
       ("h,help", "print help")
       ("positional",
        "Positional arguments: these are the arguments that are entered "
        "without an option", cxxopts::value<std::vector<std::string>>());
     options->parse_positional(std::vector<std::string>{"N0", "positional"});
     options->parse(argc, argv);
     if (options->count("help")) {
       std::cout << options->help({"", "Group"}) << std::endl;
       exit(0);
     }
     if (options->count("N0") != 1 ) {
       throw cxxopts::OptionException("Parameter N0 missing");
     } else if ((*options)["N0"].as<int>()%2 != 1) {
       throw cxxopts::OptionException("N0 must be odd");
     } else if (options->count("positional") > 0) {
       throw cxxopts::OptionException("There are too many positional arguments");
     }
   } catch (const cxxopts::OptionException & e) {
     std::cout << "Error parsing options: " << e.what() << std::endl;
     exit(1);
   }
   return options;
 }
 
 
 using namespace muSpectre;
 
 int main(int argc, char *argv[])
 {
   banner("demonstrator1", 2018, "Till Junge <till.junge@epfl.ch>");
   auto options{parse_args(argc, argv)};
   auto & opt{*options};
   const Dim_t size{opt["N0"].as<int>()};
   constexpr Real fsize{1.};
   constexpr Dim_t dim{3};
   const Dim_t nb_dofs{ipow(size, dim)*ipow(dim, 2)};
   std::cout << "Number of dofs: " <<  nb_dofs << std::endl;
 
   constexpr Formulation form{Formulation::finite_strain};
 
   const Rcoord_t<dim> lengths{CcoordOps::get_cube<dim>(fsize)};
   const Ccoord_t<dim> resolutions{CcoordOps::get_cube<dim>(size)};
 
   auto system{make_system<dim, dim>(resolutions, lengths, form)};
 
   constexpr Real E{1.0030648180242636};
   constexpr Real nu{0.29930675909878679};
 
   using Material_t = MaterialHyperElastic1<dim, dim>;
   auto Material_soft{std::make_unique<Material_t>("soft",    E, nu)};
   auto Material_hard{std::make_unique<Material_t>("hard", 10*E, nu)};
 
   int counter{0};
   for (const auto && pixel:system) {
 
     int sum = 0;
     for (Dim_t i = 0; i < dim; ++i) {
       sum  += pixel[i]*2 / resolutions[i];
     }
 
     if (sum == 0) {
       Material_hard->add_pixel(pixel);
       counter ++;
     } else {
       Material_soft->add_pixel(pixel);
     }
   }
   std::cout << counter << " Pixel out of " << system.size()
             << " are in the hard material" << std::endl;
 
   system.add_material(std::move(Material_soft));
   system.add_material(std::move(Material_hard));
   system.initialise(FFT_PlanFlags::measure);
 
   constexpr Real newton_tol{1e-4};
   constexpr Real cg_tol{1e-7};
-  const size_t maxiter = nb_dofs;
+  const Uint maxiter = nb_dofs;
 
   Grad_t<dim> DeltaF{Grad_t<dim>::Zero()};
   DeltaF(0, 1) = .1;
   Dim_t verbose {1};
 
   auto start = std::chrono::high_resolution_clock::now();
   GradIncrements<dim> grads{DeltaF};
-  de_geus(system, grads, cg_tol, newton_tol, maxiter, verbose);
+  SolverCG<dim> cg{system, cg_tol, maxiter, bool(verbose)};
+  de_geus(system, grads, cg, newton_tol, verbose);
   std::chrono::duration<Real> dur = std::chrono::high_resolution_clock::now() - start;
   std::cout << "Resolution time = " << dur.count() << "s" << std::endl;
 
   return 0;
 }
diff --git a/bin/demonstrator2.cc b/bin/demonstrator2.cc
index ca548ea..c9ee12d 100644
--- a/bin/demonstrator2.cc
+++ b/bin/demonstrator2.cc
@@ -1,91 +1,92 @@
 /**
  * file   demonstrator1.cc
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   03 Jan 2018
  *
  * @brief  larger problem to show off
  *
  * @section LICENSE
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <iostream>
 #include <memory>
 #include <chrono>
 
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "system/system_factory.hh"
 #include "materials/material_hyper_elastic1.hh"
 #include "solver/solvers.hh"
 #include "solver/solver_cg.hh"
 
 using namespace muSpectre;
 
 int main()
 {
   banner("demonstrator1", 2018, "Till Junge <till.junge@epfl.ch>");
   constexpr Dim_t dim{2};
 
   constexpr Formulation form{Formulation::finite_strain};
 
   const Rcoord_t<dim> lengths{5.2, 8.3};
   const Ccoord_t<dim> resolutions{5, 7};
 
   auto system{make_system<dim, dim>(resolutions, lengths, form)};
 
   constexpr Real E{1.0030648180242636};
   constexpr Real nu{0.29930675909878679};
 
   using Material_t = MaterialHyperElastic1<dim, dim>;
   auto & soft{Material_t::make(system, "soft",    E, nu)};
   auto & hard{Material_t::make(system, "hard", 10*E, nu)};
 
   int counter{0};
   for (const auto && pixel:system) {
     if (counter < 3) {
       hard.add_pixel(pixel);
       counter++;
     } else {
       soft.add_pixel(pixel);
     }
   }
   std::cout << counter << " Pixel out of " << system.size()
             << " are in the hard material" << std::endl;
 
   system.initialise();
 
   constexpr Real newton_tol{1e-4};
   constexpr Real cg_tol{1e-7};
   const size_t maxiter = 100;
 
   Grad_t<dim> DeltaF{Grad_t<dim>::Zero()};
   DeltaF(0, 1) = .1;
   Dim_t verbose {1};
 
   auto start = std::chrono::high_resolution_clock::now();
-  auto res = de_geus(system, DeltaF, cg_tol, newton_tol, maxiter, verbose);
+  SolverCG<dim, dim> cg{system, cg_tol, maxiter, bool(verbose)};
+  auto res = de_geus(system, DeltaF, cg, newton_tol, verbose);
   std::chrono::duration<Real> dur = std::chrono::high_resolution_clock::now() - start;
   std::cout << "Resolution time = " << dur.count() << "s" << std::endl;
 
   std::cout << res.grad.transpose() << std::endl;
   return 0;
 }
diff --git a/bin/hyper-elasticity.cc b/bin/hyper-elasticity.cc
index acfa0ad..011cdf2 100644
--- a/bin/hyper-elasticity.cc
+++ b/bin/hyper-elasticity.cc
@@ -1,91 +1,92 @@
 /**
  * file   hyper-elasticity.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   16 Jan 2018
  *
  * @brief  Recreation of GooseFFT's hyper-elasticity.py calculation
  *
  * @section LICENSE
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 
 #include "system/system_factory.hh"
 #include "materials/material_hyper_elastic1.hh"
 #include "solver/solvers.hh"
+#include "solver/solver_cg.hh"
 
 #include <iostream>
 #include <iomanip>
 
 using namespace muSpectre;
 
 
 int main()
 {
   constexpr Dim_t dim{3};
   constexpr Ccoord_t<dim> N{CcoordOps::get_cube<dim>(11)};
   constexpr Rcoord_t<dim> lens{CcoordOps::get_cube<dim>(1.)};
   constexpr Dim_t incl_size{3};
 
   auto cell{make_system(N, lens, Formulation::small_strain)};
 
   // constexpr Real K_hard{8.33}, K_soft{.833};
   // constexpr Real mu_hard{3.86}, mu_soft{.386};
   // auto E = [](Real K, Real G) {return 9*K*G / (3*K+G);}; //G is mu
   // auto nu= [](Real K, Real G) {return (3*K-2*G) / (2*(3*K+G));};
 
   // auto & hard{MaterialHyperElastic1<dim, dim>::make(cell, "hard",
   //                                                   E(K_hard, mu_hard),
   //                                                   nu(K_hard, mu_hard))};
   // auto & soft{MaterialHyperElastic1<dim, dim>::make(cell, "soft",
   //                                                   E(K_soft, mu_soft),
   //                                                   nu(K_soft, mu_soft))};
   Real ex{1e-5};
   using Mat_t = MaterialHyperElastic1<dim, dim>;
   auto & hard{Mat_t::make(cell, "hard",
                           210.*ex, .33)};
   auto & soft{Mat_t::make(cell, "soft",
                           70.*ex, .33)};
 
   for (auto pixel: cell) {
     if ((pixel[0] >= N[0]-incl_size) &&
         (pixel[1] < incl_size) &&
         (pixel[2] >= N[2]-incl_size)) {
       hard.add_pixel(pixel);
     } else {
       soft.add_pixel(pixel);
     }
   }
   std::cout << hard.size() << " pixels in the inclusion" << std::endl;
   cell.initialise();
   constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
   constexpr Dim_t maxiter{200};
   constexpr Dim_t verbose{1};
 
   Grad_t<dim> dF_bar{Grad_t<dim>::Zero()};
   dF_bar(0, 1) = 1.;
-
-  auto optimize_res = de_geus(cell, dF_bar, cg_tol, newton_tol, maxiter, verbose);
+  SolverCG<dim> cg{cell, cg_tol, maxiter, verbose};
+  auto optimize_res = de_geus(cell, dF_bar, cg, newton_tol, verbose);
 
   std::cout << "nb_cg: " << optimize_res.nb_fev << std::endl;
   std::cout << optimize_res.grad.transpose().block(0,0,10,9) << std::endl;
   return 0;
 }
diff --git a/bin/small_case.cc b/bin/small_case.cc
index fcad3d6..a0aad7d 100644
--- a/bin/small_case.cc
+++ b/bin/small_case.cc
@@ -1,83 +1,85 @@
 /**
  * file   small_case.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   12 Jan 2018
  *
  * @brief  small case for debugging
  *
  * @section LICENSE
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "common/common.hh"
 #include "common/iterators.hh"
 #include "system/system_factory.hh"
 #include "materials/material_hyper_elastic1.hh"
 #include "solver/solvers.hh"
+#include "solver/solver_cg.hh"
 
 #include <iostream>
 
 
 using namespace muSpectre;
 
 
 int main()
 {
   constexpr Dim_t dim{twoD};
 
   Ccoord_t<dim> resolution{11, 11};
 
   Rcoord_t<dim> lengths{CcoordOps::get_cube<dim>(11.)};//{5.2e-9, 8.3e-9, 8.3e-9};
   Formulation form{Formulation::finite_strain};
 
   auto rve{make_system(resolution,
                        lengths,
                        form)};
 
   auto & hard{MaterialHyperElastic1<dim, dim>::make
       (rve, "hard", 210., .33)};
   auto & soft{MaterialHyperElastic1<dim, dim>::make
       (rve, "soft",  70., .33)};
 
   for (auto && tup: akantu::enumerate(rve)) {
     auto & i = std::get<0>(tup);
     auto & pixel = std::get<1>(tup);
     if (i < 3) {
       hard.add_pixel(pixel);
     } else {
       soft.add_pixel(pixel);
     }
   }
 
   rve.initialise();
 
   Real tol{1e-6};
   Grad_t<dim> Del0{};
   Del0 <<  0, .1,
            0,  0;
 
-  Dim_t maxiter{31};
+  Uint maxiter{31};
   Dim_t verbose{3};
 
-  auto res = de_geus(rve, Del0, tol, tol, maxiter, verbose);
+  SolverCG<dim> cg{rve, tol, maxiter, bool(verbose)};
+  auto res = de_geus(rve, Del0, cg, tol, verbose);
   std::cout << res.grad.transpose() << std::endl;
   return 0;
 }
diff --git a/bin/small_case.py b/bin/small_case.py
index eaae229..1dfd5f6 100755
--- a/bin/small_case.py
+++ b/bin/small_case.py
@@ -1,92 +1,112 @@
 #!/usr/bin/env python3
 """
 file   small_case.py
 
 @author Till Junge <till.junge@epfl.ch>
 
 @date   12 Jan 2018
 
 @brief  small case for debugging
 
 @section LICENSE
 
 Copyright © 2018 Till Junge
 
 µSpectre is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License as
 published by the Free Software Foundation, either version 3, or (at
 your option) any later version.
 
 µSpectre is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with GNU Emacs; see the file COPYING. If not, write to the
 Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 Boston, MA 02111-1307, USA.
 """
 
 import sys
 import os
 import numpy as np
 
 sys.path.append(os.path.join(os.getcwd(), "language_bindings/python"))
 import pyMuSpectre as µ
 
 
-resolution = [251, 251]
+resolution = [51, 51]
 center = np.array([r//2 for r in resolution])
 incl = resolution[0]//5
 
 lengths = [7., 5.]
 formulation = µ.Formulation.small_strain
 
 rve = µ.SystemFactory(resolution,
                       lengths,
                       formulation)
 hard = µ.material.MaterialHooke2d.make(
-    rve, "hard", 210e9, .33)
+    rve, "hard", 10e9, .33)
 soft = µ.material.MaterialHooke2d.make(
     rve, "soft",  70e9, .33)
 
 
 for i, pixel in enumerate(rve):
     if np.linalg.norm(center - np.array(pixel),2)<incl:
     #if (abs(center - np.array(pixel)).max()<incl or
     #    np.linalg.norm(center/2 - np.array(pixel))<incl):
         hard.add_pixel(pixel)
     else:
         soft.add_pixel(pixel)
 
-tol = 1e-6
+tol = 1e-5
+cg_tol = 1e-8
 
 Del0 = np.array([[.0, .0],
                  [0,  .03]])
 if formulation == µ.Formulation.small_strain:
     Del0 = .5*(Del0 + Del0.T)
-maxiter = 301
-verbose = 3
-
-r = µ.solvers.de_geus(rve, Del0, tol, tol, maxiter, verbose)
-print(r.grad.T[:3])
-print(r.stress.T[:3])
-print(r.nb_fev)
-
-print(r.grad.T.shape)
-import matplotlib.pyplot as plt
-stress = r.stress.T.reshape(*resolution, 2, 2)
-def comp_von_mises(arr):
-    out_arr = np.zeros(resolution)
-    s11 = arr[:,:,0,0]
-    s22 = arr[:,:,1,1]
-    s21_2 = arr[:,:,0,1]*arr[:,:,1,0]
-
-    out_arr[:] = np.sqrt(.5*((s11-s22)**2) + s11**2 + s22**2 + 6*s21_2)
-    return out_arr
-
-von_mises = comp_von_mises(stress)
-plt.pcolormesh(von_mises)
-plt.colorbar()
-plt.show()
+maxiter = 401
+verbose = 2
+
+for solvclass in (µ.solvers.SolverCG,
+                  µ.solvers.SolverCGEigen,
+                  µ.solvers.SolverBiCGSTABEigen,
+                  µ.solvers.SolverGMRESEigen,
+                  µ.solvers.SolverDGMRESEigen,
+                  µ.solvers.SolverMINRESEigen):
+    print()
+    try:
+        solver = solvclass(rve, cg_tol, maxiter, verbose=False)
+        r = µ.solvers.newton_cg(rve, Del0, solver, tol, verbose)
+        print("nb of {} iterations: {}".format(solver.name(), r.nb_fev))
+    except RuntimeError as err:
+        print(err)
+    try:
+        solver = solvclass(rve, cg_tol, maxiter, verbose=False)
+        r = µ.solvers.de_geus(rve, Del0, solver, tol, verbose)
+        print("nb of {} iterations: {}".format(solver.name(), r.nb_fev))
+    except RuntimeError as err:
+        print(err)
+
+
+# print(r.grad.T[:3])
+# print(r.stress.T[:3])
+# 
+# print(r.grad.T.shape)
+# import matplotlib.pyplot as plt
+# stress = r.stress.T.reshape(*resolution, 2, 2)
+# def comp_von_mises(arr):
+#     out_arr = np.zeros(resolution)
+#     s11 = arr[:,:,0,0]
+#     s22 = arr[:,:,1,1]
+#     s21_2 = arr[:,:,0,1]*arr[:,:,1,0]
+# 
+#     out_arr[:] = np.sqrt(.5*((s11-s22)**2) + s11**2 + s22**2 + 6*s21_2)
+#     return out_arr
+# 
+# von_mises = comp_von_mises(stress)
+# plt.pcolormesh(von_mises)
+# plt.colorbar()
+# plt.show()
diff --git a/cmake/pybind11.cmake b/cmake/pybind11.cmake
index dcbcd04..2cc45d2 100644
--- a/cmake/pybind11.cmake
+++ b/cmake/pybind11.cmake
@@ -1,19 +1,20 @@
 download_external_project(pybind11
   URL "https://github.com/pybind/pybind11.git"
   TAG "v${_pybind11_version}"
   BACKEND GIT
   THIRD_PARTY_SRC_DIR ${_pybind11_external_dir}
+  NO_UPDATE
   )
 
 set(PYBIND11_LTO_CXX_FLAGS "" CACHE INTERNAL "")
 set(PYBIND11_LTO_LINKER_FLAGS "" CACHE INTERNAL "")
 
 add_subdirectory(${_pybind11_external_dir}/pybind11)
 
 include_directories(SYSTEM ${PYBIND11_INCLUDE_DIR})
 
 mark_as_advanced_prefix(PYBIND11)
 mark_as_advanced(USE_PYTHON_INCLUDE_DIR)
 
 set(${package}_FOUND TRUE CACHE INTERNAL "To avoid cyclic search" FORCE)
 set(${package}_FOUND_EXTERNAL TRUE CACHE INTERNAL "" FORCE)
diff --git a/language_bindings/python/bind_py_solvers.cc b/language_bindings/python/bind_py_solvers.cc
index 9bf7c39..7f2354e 100644
--- a/language_bindings/python/bind_py_solvers.cc
+++ b/language_bindings/python/bind_py_solvers.cc
@@ -1,136 +1,175 @@
 /**
  * file   bind_py_solver.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   09 Jan 2018
  *
  * @brief  python bindings for the muSpectre solvers
  *
  * @section LICENSE
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "common/common.hh"
 #include "solver/solvers.hh"
+#include "solver/solver_cg.hh"
+#include "solver/solver_cg_eigen.hh"
 
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include <pybind11/eigen.h>
 
 using namespace muSpectre;
 namespace py=pybind11;
 using namespace pybind11::literals;
 
 /**
  * Solvers instanciated for systems with equal spatial and material dimension
  */
 
+template <Dim_t sdim, class Solver>
+void add_iterative_solver_helper(py::module & mod, std::string name_start) {
+  using sys = SystemBase<sdim>;
+  std::stringstream name{};
+  name << name_start << '_' << sdim << 'd';
+  py::class_<Solver, typename Solver::Parent>(mod, name.str().c_str())
+    .def(py::init<sys&, Real, Uint, bool>(),
+         "cell"_a,
+         "tol"_a,
+         "maxiter"_a,
+         "verbose"_a)
+    .def("name", &Solver::name);
+  mod.def(name_start.c_str(),
+          [](sys& cell, Real tol, Uint maxiter, bool verbose) {
+            return std::make_unique<Solver>(cell, tol, maxiter, verbose);
+          },
+          "cell"_a,
+          "tol"_a,
+          "maxiter"_a,
+          "verbose"_a);
+}
+
+template <Dim_t sdim>
+void add_iterative_solver_dispatcher(py::module & mod) {
+  std::stringstream name{};
+  name << "SolverBase_" << sdim << 'd';
+  py::class_<SolverBase<sdim>>(mod, name.str().c_str());
+  add_iterative_solver_helper<sdim, SolverCG<sdim>>(mod, "SolverCG");
+  add_iterative_solver_helper<sdim, SolverCGEigen<sdim>>(mod, "SolverCGEigen");
+  add_iterative_solver_helper<sdim, SolverGMRESEigen<sdim>>(mod, "SolverGMRESEigen");
+  add_iterative_solver_helper<sdim, SolverBiCGSTABEigen<sdim>>(mod, "SolverBiCGSTABEigen");
+  add_iterative_solver_helper<sdim, SolverDGMRESEigen<sdim>>(mod, "SolverDGMRESEigen");
+  add_iterative_solver_helper<sdim, SolverMINRESEigen<sdim>>(mod, "SolverMINRESEigen");
+}
 
+void add_iterative_solver(py::module & mod) {
+  add_iterative_solver_dispatcher<  twoD>(mod);
+  add_iterative_solver_dispatcher<threeD>(mod);
+}
 
 template <Dim_t sdim>
 void add_newton_cg_helper(py::module & mod) {
 
   const char name []{"newton_cg"};
   constexpr Dim_t mdim{sdim};
   using sys = SystemBase<sdim, mdim>;
+  using solver = SolverBase<sdim, mdim>;
   using grad = Grad_t<sdim>;
   using grad_vec = GradIncrements<sdim>;
 
   mod.def(name,
-          [](sys & s, const grad & g, Real ct, Real nt,
-             Uint max, Dim_t verb) -> OptimizeResult {
-            return newton_cg(s, g, ct, nt, max, verb);
+          [](sys & s, const grad & g, solver & so, Real nt,
+             Dim_t verb) -> OptimizeResult {
+            return newton_cg(s, g, so, nt, verb);
 
           },
           "system"_a,
           "ΔF₀"_a,
-          "cg_tol"_a,
+          "solver"_a,
           "newton_tol"_a,
-          "maxiter"_a=0,
           "verbose"_a=0);
   mod.def(name,
-          [](sys & s, const grad_vec & g, Real ct, Real nt,
-             Uint max, Dim_t verb) -> std::vector<OptimizeResult> {
-            return newton_cg(s, g, ct, nt, max, verb);
+          [](sys & s, const grad_vec & g, solver & so, Real nt,
+             Dim_t verb) -> std::vector<OptimizeResult> {
+            return newton_cg(s, g, so, nt, verb);
           },
           "system"_a,
           "ΔF₀"_a,
-          "cg_tol"_a,
+          "solver"_a,
           "newton_tol"_a,
-          "maxiter"_a=0,
           "verbose"_a=0);
 }
 
 template <Dim_t sdim>
 void add_de_geus_helper(py::module & mod) {
   const char name []{"de_geus"};
   constexpr Dim_t mdim{sdim};
   using sys = SystemBase<sdim, mdim>;
+  using solver = SolverBase<sdim, mdim>; 
   using grad = Grad_t<sdim>;
   using grad_vec = GradIncrements<sdim>;
 
   mod.def(name,
-          [](sys & s, const grad & g, Real ct, Real nt,
-             Uint max, Dim_t verb) -> OptimizeResult {
-            return de_geus(s, g, ct, nt, max, verb);
+          [](sys & s, const grad & g, solver & so, Real nt,
+             Dim_t verb) -> OptimizeResult {
+            return de_geus(s, g, so, nt, verb);
+
           },
           "system"_a,
           "ΔF₀"_a,
-          "cg_tol"_a,
+          "solver"_a,
           "newton_tol"_a,
-          "maxiter"_a=0,
           "verbose"_a=0);
   mod.def(name,
-          [](sys & s, const grad_vec & g, Real ct, Real nt,
-             Uint max, Dim_t verb) -> std::vector<OptimizeResult> {
-            return de_geus(s, g, ct, nt, max, verb);
+          [](sys & s, const grad_vec & g, solver & so, Real nt,
+             Dim_t verb) -> std::vector<OptimizeResult> {
+            return de_geus(s, g, so, nt, verb);
           },
           "system"_a,
           "ΔF₀"_a,
-          "cg_tol"_a,
+          "solver"_a,
           "newton_tol"_a,
-          "maxiter"_a=0,
           "verbose"_a=0);
 }
 
 template <Dim_t dim>
 void add_solver_helper(py::module & mod) {
   add_newton_cg_helper<dim>(mod);
   add_de_geus_helper  <dim>(mod);
 }
 
 void add_solvers(py::module & mod) {
   auto solvers{mod.def_submodule("solvers")};
   solvers.doc() = "bindings for solvers";
 
   py::class_<OptimizeResult>(mod, "OptimizeResult")
     .def_readwrite("grad", &OptimizeResult::grad)
     .def_readwrite("stress", &OptimizeResult::stress)
     .def_readwrite("success", &OptimizeResult::success)
     .def_readwrite("status", &OptimizeResult::status)
     .def_readwrite("message", &OptimizeResult::message)
     .def_readwrite("nb_it", &OptimizeResult::nb_it)
     .def_readwrite("nb_fev", &OptimizeResult::nb_fev);
 
-
+  add_iterative_solver(solvers);
 
   add_solver_helper<twoD  >(solvers);
   add_solver_helper<threeD>(solvers);
 }
diff --git a/src/common/field.hh b/src/common/field.hh
index 205c53e..24d33f8 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,807 +1,817 @@
 /**
  * file   field.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   07 Sep 2017
  *
  * @brief  header-only implementation of a field for field collections
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 
 #ifndef FIELD_H
 #define FIELD_H
 
 #include "common/T4_map_proxy.hh"
 
 #include <Eigen/Dense>
 
 #include <string>
 #include <sstream>
 #include <utility>
 #include <typeinfo>
 #include <vector>
 #include <algorithm>
 #include <cmath>
 #include <memory>
 #include <type_traits>
 
 namespace muSpectre {
 
   /* ---------------------------------------------------------------------- */
   class FieldCollectionError: public std::runtime_error {
   public:
     explicit FieldCollectionError(const std::string& what)
       :std::runtime_error(what){}
     explicit FieldCollectionError(const char * what)
       :std::runtime_error(what){}
   };
 
   class FieldError: public FieldCollectionError {
     using Parent = FieldCollectionError;
   public:
     explicit FieldError(const std::string& what)
       :Parent(what){}
     explicit FieldError(const char * what)
       :Parent(what){}
   };
   class FieldInterpretationError: public FieldError
   {
   public:
     explicit FieldInterpretationError(const std::string & what)
       :FieldError(what){}
     explicit FieldInterpretationError(const char * what)
       :FieldError(what){}
   };
 
 
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     class FieldBase
     {
 
     protected:
       //! constructor
       //! unique name (whithin Collection)
       //! number of components
       //! collection to which this field belongs (eg, material, system)
       FieldBase(std::string unique_name,
                 size_t nb_components,
                 FieldCollection & collection);
 
     public:
       using collection_t = FieldCollection;
 
       //! Copy constructor
       FieldBase(const FieldBase &other) = delete;
 
       //! Move constructor
       FieldBase(FieldBase &&other) = delete;
 
       //! Destructor
       virtual ~FieldBase() = default;
 
       //! Copy assignment operator
       FieldBase& operator=(const FieldBase &other) = delete;
 
       //! Move assignment operator
       FieldBase& operator=(FieldBase &&other) = delete;
 
       /* ---------------------------------------------------------------------- */
       //!Identifying accessors
       //! return field name
       inline const std::string & get_name() const;
       //! return field type
       //inline const Field_t & get_type() const;
       //! return my collection (for iterating)
       inline const FieldCollection & get_collection() const;
       //! return my collection (for iterating)
       inline const size_t & get_nb_components() const;
       //! return type_id of stored type
       virtual const std::type_info & get_stored_typeid() const = 0;
 
       virtual size_t size() const = 0;
 
       //! initialise field to zero (do more complicated initialisations through
       //! fully typed maps)
       virtual void set_zero() = 0;
 
       //! give access to collections
       friend FieldCollection;
       friend typename FieldCollection::Parent;
 
     protected:
       /* ---------------------------------------------------------------------- */
       //! allocate memory etc
       virtual void resize(size_t size) = 0;
       const std::string name;
       const size_t nb_components;
       const FieldCollection & collection;
     private:
     };
 
 
     /**
      * dummy intermediate class to provide a run-time polymorphic
      * typed field. Mainly for binding python
      */
     template <class FieldCollection, typename T>
     class TypedFieldBase: public FieldBase<FieldCollection>
     {
     public:
       using Parent = FieldBase<FieldCollection>;
       using collection_t = typename Parent::collection_t;
       using Scalar = T;
       using Base = Parent;
       using EigenRep = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>;
       using EigenMap = Eigen::Map<EigenRep>;
       using EigenVec = Eigen::Map<Eigen::VectorXd>;
+      using EigenVecConst = Eigen::Map<const Eigen::VectorXd>;
       //! Default constructor
       TypedFieldBase() = delete;
 
       TypedFieldBase(std::string unique_name,
                      size_t nb_components,
                      FieldCollection& collection);
 
       //! Copy constructor
       TypedFieldBase(const TypedFieldBase &other) = delete;
 
       //! Move constructor
       TypedFieldBase(TypedFieldBase &&other) = delete;
 
       //! Destructor
       virtual ~TypedFieldBase() = default;
 
       //! Copy assignment operator
       TypedFieldBase& operator=(const TypedFieldBase &other) = delete;
 
       //! Move assignment operator
       TypedFieldBase& operator=(TypedFieldBase &&other) = delete;
 
       //! return type_id of stored type
       virtual const std::type_info & get_stored_typeid() const override final;
 
       virtual size_t size() const override = 0;
 
       //! initialise field to zero (do more complicated initialisations through
       //! fully typed maps)
       virtual void set_zero() override = 0;
       virtual T* data() = 0;
       virtual const T* data() const = 0;
 
       EigenMap eigen();
       EigenVec eigenvec();
+      EigenVecConst eigenvec() const;
 
     protected:
     private:
     };
 
     /* ---------------------------------------------------------------------- */
     //! declaraton for friending
     template <class FieldCollection, typename T, Dim_t NbComponents, bool isConst>
     class FieldMap;
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool ArrayStore=false>
     class TypedSizedFieldBase: public TypedFieldBase<FieldCollection, T>
     {
       friend class FieldMap<FieldCollection, T, NbComponents, true>;
       friend class FieldMap<FieldCollection, T, NbComponents, false>;
     public:
       constexpr static auto nb_components{NbComponents};
       using Parent = TypedFieldBase<FieldCollection, T>;
       using collection_t = typename Parent::collection_t;
       using Scalar = T;
       using Base = typename Parent::Base;
       //using storage_type = Eigen::Array<T, Eigen::Dynamic, NbComponents>;
       using StoredType = Eigen::Array<T, NbComponents, 1>;
       using StorageType = std::conditional_t
         <ArrayStore,
          std::vector<StoredType,
                      Eigen::aligned_allocator<StoredType>>,
          std::vector<T,Eigen::aligned_allocator<T>>>;
 
       using EigenRep = Eigen::Array<T, NbComponents, Eigen::Dynamic>;
       using EigenMap = std::conditional_t<
         ArrayStore,
         Eigen::Map<EigenRep, Eigen::Aligned,
                    Eigen::OuterStride<sizeof(StoredType)/sizeof(T)>>,
         Eigen::Map<EigenRep>>;
 
       using ConstEigenMap = std::conditional_t<
         ArrayStore,
         Eigen::Map<const EigenRep, Eigen::Aligned,
                    Eigen::OuterStride<sizeof(StoredType)/sizeof(T)>>,
         Eigen::Map<const EigenRep>>;
 
       TypedSizedFieldBase(std::string unique_name,
                      FieldCollection& collection);
       virtual ~TypedSizedFieldBase() = default;
 
       //! initialise field to zero (do more complicated initialisations through
       //! fully typed maps)
       inline void set_zero() override final;
       template <bool isArrayStore = ArrayStore>
       inline void push_back(const
                             std::enable_if_t<isArrayStore, StoredType> &
                             value);
       template <bool componentStore = !ArrayStore>
       inline std::enable_if_t<componentStore> push_back(const StoredType & value);
 
       //! Number of stored arrays (i.e. total number of stored
       //! scalars/NbComponents)
       size_t size() const override final;
 
       static TypedSizedFieldBase & check_ref(Base & other);
       static const TypedSizedFieldBase & check_ref(const Base & parent);
 
       inline T* data() override final {return this->get_ptr_to_entry(0);}
       inline const T* data() const override final {return this->get_ptr_to_entry(0);}
 
       inline EigenMap eigen();
       inline ConstEigenMap eigen() const;
       inline typename Parent::EigenMap dyn_eigen() {return Parent::eigen();}
 
       template<typename otherT>
       inline Real inner_product(const TypedSizedFieldBase<FieldCollection, otherT, NbComponents,
                                 ArrayStore> & other) const;
     protected:
 
       template <bool isArray=ArrayStore>
       inline std::enable_if_t<isArray, T*> get_ptr_to_entry(const size_t&& index);
 
       template <bool isArray=ArrayStore>
       inline std::enable_if_t<isArray, const T*>
       get_ptr_to_entry(const size_t&& index) const;
 
       template <bool noArray = !ArrayStore>
       inline T* get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index);
 
       template <bool noArray = !ArrayStore>
       inline const T*
       get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const;
 
       inline virtual void resize(size_t size) override final;
       StorageType values{};
     };
 
   }  // internal
 
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   class TensorField: public internal::TypedSizedFieldBase<FieldCollection,
                                                      T,
                                                      ipow(dim,order)>
   {
   public:
     using Parent = internal::TypedSizedFieldBase<FieldCollection,
                                             T,
                                             ipow(dim,order)>;
     using Base = typename Parent::Base;
     using Field_p = typename FieldCollection::Field_p;
     using component_type = T;
     using Scalar = typename Parent::Scalar;
     //! Copy constructor
     TensorField(const TensorField &other) = delete;
 
     //! Move constructor
     TensorField(TensorField &&other) = delete;
 
     //! Destructor
     virtual ~TensorField() = default;
 
     //! Copy assignment operator
     TensorField& operator=(const TensorField &other) = delete;
 
     //! Move assignment operator
     TensorField& operator=(TensorField &&other) = delete;
 
     //! accessors
     inline Dim_t get_order() const;
     inline Dim_t get_dim() const;
 
 
     //! factory function
     template<class FieldType, class CollectionType, typename... Args>
     friend FieldType& make_field(std::string unique_name,
                              CollectionType & collection,
                              Args&&... args);
 
     static TensorField & check_ref(Base & other) {
       return static_cast<TensorField &>(Parent::check_ref(other));}
     static const TensorField & check_ref(const Base & other) {
       return static_cast<const TensorField &>(Parent::check_ref(other));}
 
     /**
      * Pure convenience functions to get a MatrixFieldMap of
      * appropriate dimensions mapped to this field. You can also
      * create other types of maps, as long as they have the right
      * fundamental type (T) and the correct size (nbComponents).
      */
     decltype(auto) get_map();
     decltype(auto) get_map() const;
 
   protected:
     //! constructor protected!
     TensorField(std::string unique_name, FieldCollection & collection);
 
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol=NbRow>
   class MatrixField: public internal::TypedSizedFieldBase<FieldCollection,
                                                      T,
                                                      NbRow*NbCol>
   {
   public:
     using Parent = internal::TypedSizedFieldBase<FieldCollection,
                                             T,
                                             NbRow*NbCol>;
     using Base = typename Parent::Base;
     using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
     using component_type = T;
     //! Copy constructor
     MatrixField(const MatrixField &other) = delete;
 
     //! Move constructor
     MatrixField(MatrixField &&other) = delete;
 
     //! Destructor
     virtual ~MatrixField() = default;
 
     //! Copy assignment operator
     MatrixField& operator=(const MatrixField &other) = delete;
 
     //! Move assignment operator
     MatrixField& operator=(MatrixField &&other) = delete;
 
     //! accessors
     inline Dim_t get_nb_row() const;
     inline Dim_t get_nb_col() const;
 
 
     //! factory function
     template<class FieldType, class CollectionType, typename... Args>
     friend FieldType&  make_field(std::string unique_name,
                              CollectionType & collection,
                              Args&&... args);
 
     static MatrixField & check_ref(Base & other) {
       return static_cast<MatrixField &>(Parent::check_ref(other));}
     static const MatrixField & check_ref(const Base & other) {
       return static_cast<const MatrixField &>(Parent::check_ref(other));}
 
     decltype(auto) get_map();
     decltype(auto) get_map() const;
 
 
 protected:
     //! constructor protected!
     MatrixField(std::string unique_name, FieldCollection & collection);
 
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   //! convenience alias (
   template <class FieldCollection, typename T>
   using ScalarField = MatrixField<FieldCollection, T, 1, 1>;
   /* ---------------------------------------------------------------------- */
   // Implementations
   namespace internal {
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     FieldBase<FieldCollection>::FieldBase(std::string unique_name,
                                           size_t nb_components_,
                                           FieldCollection & collection_)
       :name(unique_name), nb_components(nb_components_),
        collection(collection_) {}
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     inline const std::string & FieldBase<FieldCollection>::get_name() const {
       return this->name;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     inline const FieldCollection & FieldBase<FieldCollection>::
     get_collection() const {
       return this->collection;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     inline const size_t & FieldBase<FieldCollection>::
     get_nb_components() const {
       return this->nb_components;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T>
     TypedFieldBase<FieldCollection, T>::
     TypedFieldBase(std::string unique_name, size_t nb_components,
                    FieldCollection & collection)
       :Parent(unique_name, nb_components, collection)
     {}
 
     /* ---------------------------------------------------------------------- */
     //! return type_id of stored type
     template <class FieldCollection, typename T>
     const std::type_info & TypedFieldBase<FieldCollection, T>::
     get_stored_typeid() const {
       return typeid(T);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T>
     typename TypedFieldBase<FieldCollection, T>::EigenMap
     TypedFieldBase<FieldCollection, T>::
     eigen() {
       return EigenMap(this->data(), this->get_nb_components(), this->size());
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T>
     typename TypedFieldBase<FieldCollection, T>::EigenVec
     TypedFieldBase<FieldCollection, T>::
     eigenvec() {
       return EigenVec(this->data(), this->get_nb_components() * this->size());
     }
 
+    /* ---------------------------------------------------------------------- */
+    template <class FieldCollection, typename T>
+    typename TypedFieldBase<FieldCollection, T>::EigenVecConst
+    TypedFieldBase<FieldCollection, T>::
+    eigenvec() const{
+      return EigenVecConst(this->data(), this->get_nb_components() * this->size());
+    }
+
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     TypedSizedFieldBase(std::string unique_name, FieldCollection & collection)
       :Parent(unique_name, NbComponents, collection){
       static_assert
         ((std::is_arithmetic<T>::value ||
           std::is_same<T, Complex>::value),
          "Use TypedSizedFieldBase for integer, real or complex scalars for T");
       static_assert(NbComponents > 0,
                     "Only fields with more than 0 components");
     }
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     void
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     set_zero() {
       std::fill(this->values.begin(), this->values.end(), T{});
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     size_t TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     size() const {
       if (ArrayStore) {
         return this->values.size();
       } else  {
         return this->values.size()/NbComponents;
       }
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     check_ref(Base & other) {
       if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
         std::string err ="Cannot create a Reference of requested type " +(
            "for field '" + other.get_name() + "' of type '" +
            other.get_stored_typeid().name() + "'");
         throw std::runtime_error
           (err);
       }
       //check size compatibility
       if (NbComponents != other.get_nb_components()) {
         throw std::runtime_error
           ("Cannot create a Reference to a field with " +
            std::to_string(NbComponents) + " components " +
            "for field '" + other.get_name() + "' with " +
            std::to_string(other.get_nb_components()) + " components");
       }
       return static_cast<TypedSizedFieldBase&>(other);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool ArrayStore>
     const TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     check_ref(const Base & other) {
       if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
         std::stringstream err_str{};
         err_str << "Cannot create a Reference of requested type "
                 << "for field '"  << other.get_name() << "' of type '"
                 << other.get_stored_typeid().name() << "'";
         throw std::runtime_error
           (err_str.str());
       }
       //check size compatibility
       if (NbComponents != other.get_nb_components()) {
         throw std::runtime_error
           ("Cannot create a Reference to a field with " +
            std::to_string(NbComponents) + " components " +
            "for field '" + other.get_name() + "' with " +
            std::to_string(other.get_nb_components()) + " components");
       }
       return static_cast<const TypedSizedFieldBase&>(other);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool ArrayStore>
     typename TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::EigenMap
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     eigen() {
       return EigenMap(this->data(), NbComponents, this->size());
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool ArrayStore>
     typename TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::ConstEigenMap
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     eigen() const {
       return ConstEigenMap(this->data(), NbComponents, this->size());
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool ArrayStore>
     template<typename otherT>
     Real
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     inner_product(const TypedSizedFieldBase<FieldCollection, otherT, NbComponents,
                   ArrayStore> & other) const {
       return (this->eigen() * other.eigen()).sum();
     }
 
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool isArray>
     std::enable_if_t<isArray, T*>
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(const size_t&& index) {
       static_assert (isArray == ArrayStore, "SFINAE");
       return &this->values[std::move(index)](0, 0);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool noArray>
     T* TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) {
       static_assert (noArray != ArrayStore, "SFINAE");
       return &this->values[NbComponents*std::move(index)];
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool isArray>
     std::enable_if_t<isArray, const T*>
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(const size_t&& index) const {
       static_assert (isArray == ArrayStore, "SFINAE");
       return &this->values[std::move(index)](0, 0);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool noArray>
     const T* TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const {
       static_assert (noArray != ArrayStore, "SFINAE");
       return &this->values[NbComponents*std::move(index)];
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     void TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     resize(size_t size) {
       if (ArrayStore) {
         this->values.resize(size);
       } else {
         this->values.resize(size*NbComponents);
       }
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool isArrayStore>
     void
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     push_back(const std::enable_if_t<isArrayStore,StoredType> & value) {
       static_assert(isArrayStore == ArrayStore, "SFINAE");
       this->values.push_back(value);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool componentStore>
     std::enable_if_t<componentStore>
     TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     push_back(const StoredType & value) {
       static_assert(componentStore != ArrayStore, "SFINAE");
       for (Dim_t i = 0; i < NbComponents; ++i) {
         this->values.push_back(value(i));
       }
     }
 
   }  // internal
 
   /* ---------------------------------------------------------------------- */
   //! Factory function, guarantees that only fields get created that are
   //! properly registered and linked to a collection.
   template<class FieldType, class FieldCollection, typename... Args>
   FieldType &
   make_field(std::string unique_name,
              FieldCollection & collection,
              Args&&... args) {
     auto ptr = std::unique_ptr<FieldType>{
       new FieldType(unique_name, collection, args...)};
     auto& retref = *ptr;
     collection.register_field(std::move(ptr));
     return retref;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   TensorField<FieldCollection, T, order, dim>::
   TensorField(std::string unique_name, FieldCollection & collection)
     :Parent(unique_name, collection) {}
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   Dim_t TensorField<FieldCollection, T, order, dim>::
   get_order() const {
     return order;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   Dim_t TensorField<FieldCollection, T, order, dim>::
   get_dim() const {
     return dim;
   }
 
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   MatrixField<FieldCollection, T, NbRow, NbCol>::
   MatrixField(std::string unique_name, FieldCollection & collection)
     :Parent(unique_name, collection) {}
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
   get_nb_col() const {
     return NbCol;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
   get_nb_row() const {
     return NbRow;
   }
 
 }  // muSpectre
 
 #include "common/field_map.hh"
 
 namespace muSpectre {
 
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, size_t order, Dim_t dim,
               bool ConstMap>
     struct tensor_map_type {
     };
 
     template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
     struct tensor_map_type<FieldCollection, T, firstOrder, dim, ConstMap> {
       using type = MatrixFieldMap<FieldCollection, T, dim, 1, ConstMap>;
     };
 
     template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
     struct tensor_map_type<FieldCollection, T, secondOrder, dim, ConstMap> {
       using type = MatrixFieldMap<FieldCollection, T, dim, dim, ConstMap>;
     };
 
     template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
     struct tensor_map_type<FieldCollection, T, fourthOrder, dim, ConstMap> {
       using type = T4MatrixFieldMap<FieldCollection, T, dim, ConstMap>;
     };
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol,
               bool ConstMap>
     struct matrix_map_type {
       using type =
         MatrixFieldMap<FieldCollection, T, NbRow, NbCol, ConstMap>;
     };
 
     template <class FieldCollection, typename T, bool ConstMap>
     struct matrix_map_type<FieldCollection, T, oneD, oneD, ConstMap> {
       using type = ScalarFieldMap<FieldCollection, T, ConstMap>;
     };
 
   }  // internal
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   decltype(auto) TensorField<FieldCollection, T, order, dim>::
   get_map() {
     constexpr bool map_constness{false};
     using RawMap_t =
       typename internal::tensor_map_type<FieldCollection, T, order, dim,
                                          map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   decltype(auto) TensorField<FieldCollection, T, order, dim>::
   get_map() const {
     constexpr bool map_constness{true};
     using RawMap_t =
       typename internal::tensor_map_type<FieldCollection, T, order, dim,
                                          map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
   get_map() {
     constexpr bool map_constness{false};
     using RawMap_t =
       typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
                                          map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
   get_map() const {
     constexpr bool map_constness{true};
     using RawMap_t =
       typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
                                          map_constness>::type;
     return RawMap_t(*this);
   }
 
 
 }  // muSpectre
 
 #endif /* FIELD_H */
diff --git a/src/solver/CMakeLists.txt b/src/solver/CMakeLists.txt
index bd2c1d0..e68a7d7 100644
--- a/src/solver/CMakeLists.txt
+++ b/src/solver/CMakeLists.txt
@@ -1,39 +1,40 @@
 # =============================================================================
 # file   CMakeLists.txt
 #
 # @author Till Junge <till.junge@epfl.ch>
 #
 # @date   08 Jan 2018
 #
 # @brief  configuration for solvers
 #
 # @section LICENSE
 #
 # Copyright © 2018 Till Junge
 #
 # µSpectre is free software; you can redistribute it and/or
 # modify it under the terms of the GNU General Public License as
 # published by the Free Software Foundation, either version 3, or (at
 # your option) any later version.
 #
 # µSpectre is distributed in the hope that it will be useful, but
 # WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 # General Public License for more details.
 #
 # You should have received a copy of the GNU General Public License
 # along with GNU Emacs; see the file COPYING. If not, write to the
 # Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 # Boston, MA 02111-1307, USA.
 # =============================================================================
 
 set (solvers_SRC
   ${CMAKE_CURRENT_SOURCE_DIR}/solver_base.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/solver_cg.cc
+  ${CMAKE_CURRENT_SOURCE_DIR}/solver_cg_eigen.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/solvers.cc
   )
 
 set (muSpectre_SRC
   ${muSpectre_SRC}
   ${solvers_SRC}
   PARENT_SCOPE)
diff --git a/src/solver/solver_base.cc b/src/solver/solver_base.cc
index 2d30f0c..ee2be48 100644
--- a/src/solver/solver_base.cc
+++ b/src/solver/solver_base.cc
@@ -1,40 +1,65 @@
 /**
  * file   solver_base.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   18 Dec 2017
  *
  * @brief  definitions for solvers
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "solver/solver_base.hh"
 #include "solver/solver_cg.hh"
 #include "common/field.hh"
 #include "common/iterators.hh"
 
 #include <iostream>
 #include <memory>
 
 
 namespace muSpectre {
+
+  //----------------------------------------------------------------------------//
+  template <Dim_t DimS, Dim_t DimM>
+  SolverBase<DimS, DimM>::SolverBase(Sys_t & sys, Real tol, Uint maxiter,
+                                     bool verbose )
+    : sys{sys}, tol{tol}, maxiter{maxiter}, verbose{verbose}
+  {}
+
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void SolverBase<DimS, DimM>::reset_counter() {
+    this->counter = 0;
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  Uint SolverBase<DimS, DimM>::get_counter() const {
+    return this->counter;
+  }
+
+  template class SolverBase<twoD, twoD>;
+  //template class SolverBase<twoD, threeD>;
+  template class SolverBase<threeD, threeD>;
+
 }  // muSpectre
diff --git a/src/solver/solver_base.hh b/src/solver/solver_base.hh
index e829c25..f405590 100644
--- a/src/solver/solver_base.hh
+++ b/src/solver/solver_base.hh
@@ -1,98 +1,125 @@
 /**
  * file   solver_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   18 Dec 2017
  *
  * @brief  Base class for solvers
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef SOLVER_BASE_H
 #define SOLVER_BASE_H
 
 #include "solver/solver_error.hh"
 #include "common/common.hh"
 #include "system/system_base.hh"
 #include "common/tensor_algebra.hh"
 
 #include <Eigen/Dense>
 
 #include <vector>
 
 namespace muSpectre {
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM=DimS>
   class SolverBase
   {
   public:
     enum class TangentRequirement{NoNeed, NeedEffect, NeedTangents};
+    using Sys_t = SystemBase<DimS, DimM>;
     using Ccoord = Ccoord_t<DimS>;
     using Collection_t = GlobalFieldCollection<DimS, DimM>;
+    using SolvVectorIn = Eigen::Ref<Eigen::VectorXd>;
+    using SolvVectorInC = Eigen::Ref<const Eigen::VectorXd>;
+    using SolvVectorOut = Eigen::VectorXd;
+
 
     //! Default constructor
     SolverBase() = delete;
 
     //! Constructor with domain resolutions
-    SolverBase(Ccoord resolutions): resolutions{resolutions} {}
+    SolverBase(Sys_t & sys, Real tol, Uint maxiter=0, bool verbose =false);
 
     //! Copy constructor
     SolverBase(const SolverBase &other) = delete;
 
     //! Move constructor
     SolverBase(SolverBase &&other) = default;
 
     //! Destructor
     virtual ~SolverBase() = default;
 
     //! Copy assignment operator
     SolverBase& operator=(const SolverBase &other) = delete;
 
     //! Move assignment operator
     SolverBase& operator=(SolverBase &&other) = default;
 
     //! Allocate fields used during the solution
-    void initialise() {this->collection.initialise(resolutions);}
+    virtual void initialise() {this->collection.initialise(this->sys.get_resolutions());}
 
     bool need_tangents() const {
       return (this->get_tangent_req() == TangentRequirement::NeedTangents);}
 
     bool need_effect() const {
       return (this->get_tangent_req() == TangentRequirement::NeedEffect);}
 
     bool no_need_tangent() const {
       return (this->get_tangent_req() == TangentRequirement::NoNeed);}
 
-    bool has_converged() const {return this->converged;}
+    virtual bool has_converged() const = 0;
+
+    //! reset the iteration counter to zero
+    void reset_counter();
+
+    //! get the count of how many solve steps have been executed since
+    //! construction of most recent counter reset
+    Uint get_counter() const;
+
+    virtual SolvVectorOut solve(const SolvVectorInC rhs, SolvVectorIn x_0) = 0;
+
+    Sys_t & get_system() {return sys;}
+
+    Uint get_maxiter() const {return this->maxiter;}
+    void set_maxiter(Uint val) {this->maxiter = val;}
+
+    Real get_tol() const {return this->tol;}
+    void set_tol(Real val) {this->tol = val;}
+
+    virtual std::string name() const = 0;
 
   protected:
     virtual TangentRequirement get_tangent_req() const = 0;
-    Ccoord resolutions;
+    Sys_t & sys;
+    Real tol;
+    Uint maxiter;
+    bool verbose;
+    Uint counter{0};
     //! storage for internal fields to avoid reallocations between calls
     Collection_t collection{};
-    bool converged{false};
   private:
   };
 
 }  // muSpectre
 
 #endif /* SOLVER_BASE_H */
diff --git a/src/solver/solver_cg.cc b/src/solver/solver_cg.cc
index 3d21e90..84deeb3 100644
--- a/src/solver/solver_cg.cc
+++ b/src/solver/solver_cg.cc
@@ -1,127 +1,131 @@
 /**
  * file   solver_cg.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief  Implementation of cg solver
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "solver/solver_cg.hh"
 #include "solver/solver_error.hh"
 
 #include <iomanip>
 #include <cmath>
+#include <sstream>
 
 namespace muSpectre {
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
-  SolverCG<DimS, DimM>::SolverCG(Ccoord resolutions, Real tol, Uint maxiter,
+  SolverCG<DimS, DimM>::SolverCG(Sys_t& sys, Real tol, Uint maxiter,
                                  bool verbose)
-    :Parent(resolutions),
+    :Parent(sys, tol, maxiter, verbose),
      r_k{make_field<Field_t>("residual r_k", this->collection)},
      p_k{make_field<Field_t>("search direction r_k", this->collection)},
-     Ap_k{make_field<Field_t>("Effect of tangent A*p_k", this->collection)},
-     tol{tol}, maxiter{maxiter}, verbose{verbose}, counter{0}
+     Ap_k{make_field<Field_t>("Effect of tangent A*p_k", this->collection)}
   {}
 
+
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
-  void SolverCG<DimS, DimM>::solve(const Fun_t& tangent_effect,
-                                   const Field_t & rhs,
-                              Field_t & x_f) {
+  void SolverCG<DimS, DimM>::solve(const Field_t & rhs,
+                                   Field_t & x_f) {
+    x_f.eigenvec() = this-> solve(rhs.eigenvec(), x_f.eigenvec());
+  };
+
+  //----------------------------------------------------------------------------//
+  template <Dim_t DimS, Dim_t DimM>
+  typename SolverCG<DimS, DimM>::SolvVectorOut
+  SolverCG<DimS, DimM>::solve(const SolvVectorInC rhs, SolvVectorIn x_0) {
+
     // Following implementation of algorithm 5.2 in Nocedal's Numerical Optimization (p. 112)
 
     auto r = this->r_k.eigen();
     auto p = this->p_k.eigen();
     auto Ap = this->Ap_k.eigen();
-    auto x = x_f.eigen();
+    auto x = typename Field_t::EigenMap(x_0.data(), r.rows(), r.cols());
 
     // initialisation of algo
-    tangent_effect(x_f, this->r_k);
+    r = this->sys.directional_stiffness_with_copy(x);
+
 
-    r -= rhs.eigen();
+    r -= typename Field_t::ConstEigenMap(rhs.data(), r.rows(), r.cols());
     p = -r;
 
     this->converged = false;
     Real rdr = (r*r).sum();
-    Real rhs_norm2 = (rhs.eigen()*rhs.eigen()).sum();
+    Real rhs_norm2 = rhs.squaredNorm();
     Real tol2 = ipow(this->tol,2)*rhs_norm2;
 
 
     size_t count_width{}; // for output formatting in verbose case
     if (this->verbose) {
       count_width = size_t(std::log10(this->maxiter))+1;
     }
 
     for (Uint i = 0; i < this->maxiter && rdr > tol2; ++i, ++this->counter) {
-      tangent_effect(this->p_k, this->Ap_k);
+      Ap = this->sys.directional_stiffness_with_copy(p);
 
       Real alpha = rdr/(p*Ap).sum();
 
       x += alpha * p;
       r += alpha * Ap;
 
       Real new_rdr = (r*r).sum();
       Real beta = new_rdr/rdr;
       rdr = new_rdr;
 
       if (this->verbose) {
         std::cout << "  at CG step " << std::setw(count_width) << i
                   << ": |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2)
                   << ", cg_tol = " << this->tol << std::endl;
       }
 
       p = -r+beta*p;
     }
     if (rdr < tol2) {
       this->converged=true;
     } else {
-      throw ConvergenceError("Conjugate gradient has not converged");
+      std::stringstream err {};
+      err << " After " << this->counter << " steps, the solver "
+          << " FAILED with  |r|/|b| = "
+          << std::setw(15) << std::sqrt(rdr/rhs_norm2)
+          << ", cg_tol = " << this->tol << std::endl;
+      throw ConvergenceError("Conjugate gradient has not converged." + err.str());
     }
+    return x_0;
   }
 
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM>
-  void SolverCG<DimS, DimM>::reset_counter() {
-    this->counter = 0;
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM>
-  Uint SolverCG<DimS, DimM>::get_counter() const {
-    return this->counter;
-  }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   typename SolverCG<DimS, DimM>::Tg_req_t
   SolverCG<DimS, DimM>::get_tangent_req() const {
     return tangent_requirement;
   }
 
   template class SolverCG<twoD, twoD>;
   //template class SolverCG<twoD, threeD>;
   template class SolverCG<threeD, threeD>;
 }  // muSpectre
diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh
index 0920fbb..46ceca5 100644
--- a/src/solver/solver_cg.hh
+++ b/src/solver/solver_cg.hh
@@ -1,104 +1,103 @@
 /**
  * file   solver_cg.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief class for a simple implementation of a conjugate gradient
  *        solver. This follows algorithm 5.2 in Nocedal's Numerical
  *        Optimization (p 112)
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef SOLVER_CG_H
 #define SOLVER_CG_H
 
 #include "solver/solver_base.hh"
 #include "common/field.hh"
 
 #include <functional>
 
 namespace muSpectre {
 
   template <Dim_t DimS, Dim_t DimM=DimS>
   class SolverCG: public SolverBase<DimS, DimM>
   {
   public:
     using Parent = SolverBase<DimS, DimM>;
+    using SolvVectorIn = typename Parent::SolvVectorIn;
+    using SolvVectorInC = typename Parent::SolvVectorInC;
+    using SolvVectorOut = typename Parent::SolvVectorOut;
+    using Sys_t = typename Parent::Sys_t;
     using Ccoord = typename Parent::Ccoord;
     using Tg_req_t = typename Parent::TangentRequirement;
     // cg only needs to handle fields that look like strain and stress
     using Field_t = TensorField<
       typename Parent::Collection_t, Real, secondOrder, DimM>;
     using Fun_t = std::function<void(const Field_t& in, Field_t & out)>;
 
     constexpr static Tg_req_t tangent_requirement{Tg_req_t::NeedEffect};
     //! Default constructor
     SolverCG() = delete;
 
     //! Constructor with domain resolutions, etc,
-    SolverCG(Ccoord resolutions, Real tol, Uint maxiter=0, bool verbose =false);
+    SolverCG(Sys_t& sys, Real tol, Uint maxiter=0, bool verbose =false);
 
     //! Copy constructor
     SolverCG(const SolverCG &other) = delete;
 
     //! Move constructor
     SolverCG(SolverCG &&other) = default;
 
     //! Destructor
     virtual ~SolverCG() = default;
 
     //! Copy assignment operator
     SolverCG& operator=(const SolverCG &other) = delete;
 
     //! Move assignment operator
     SolverCG& operator=(SolverCG &&other) = default;
 
+    bool has_converged() const override final {return this->converged;}
+
     //! actual solver
-    void solve(const Fun_t & tangent_effect,
-               const Field_t & rhs,
+    void solve(const Field_t & rhs,
                Field_t & x);
 
-    //! reset the iteration counter to zero
-    void reset_counter();
-
-    //! get the count of how many solve steps have been executed since
-    //! construction of most recent counter reset
-    Uint get_counter() const;
+    // this simplistic implementation has no initialisation phase so the default is ok
 
+    SolvVectorOut solve(const SolvVectorInC rhs, SolvVectorIn x_0) override final;
 
+    std::string name() const override final {return "CG";}
 
   protected:
     Tg_req_t get_tangent_req() const override final;
     Field_t & r_k; // residual
     Field_t & p_k; // search direction
     Field_t & Ap_k; // effect of tangent on search direction
-    const Real tol;
-    const Uint maxiter;
-    const bool verbose;
-    Uint counter{};
+    bool converged{false};
   private:
   };
 
 }  // muSpectre
 
 #endif /* SOLVER_CG_H */
diff --git a/src/solver/solver_cg_eigen.cc b/src/solver/solver_cg_eigen.cc
new file mode 100644
index 0000000..4c571ce
--- /dev/null
+++ b/src/solver/solver_cg_eigen.cc
@@ -0,0 +1,112 @@
+/**
+ * file   solver_cg_eigen.cc
+ *
+ * @author Till Junge <till.junge@epfl.ch>
+ *
+ * @date   19 Jan 2018
+ *
+ * @brief  implementation for binding to Eigen's conjugate gradient solver
+ *
+ * @section LICENCE
+ *
+ * Copyright (C) 2018 Till Junge
+ *
+ * µSpectre is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3, or (at
+ * your option) any later version.
+ *
+ * µSpectre is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Emacs; see the file COPYING. If not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+#include "solver_cg_eigen.hh"
+
+#include <sstream>
+#include <iomanip>
+
+namespace muSpectre {
+
+  //----------------------------------------------------------------------------//
+  template <class SolverType, Dim_t DimS, Dim_t DimM>
+  SolverEigen<SolverType, DimS, DimM>::SolverEigen(Sys_t& sys, Real tol, Uint maxiter,
+                                                   bool verbose)
+    :Parent(sys, tol, maxiter, verbose),
+     adaptor{sys.get_adaptor()},
+     solver{}
+  {}
+
+  //----------------------------------------------------------------------------//
+  template <class SolverType, Dim_t DimS, Dim_t DimM>
+  void
+  SolverEigen<SolverType, DimS, DimM>::initialise() {
+    this->solver.setTolerance(this->tol);
+    this->solver.setMaxIterations(this->maxiter);
+    this->solver.compute(this->adaptor);
+  }
+
+  //----------------------------------------------------------------------------//
+  template <class SolverType, Dim_t DimS, Dim_t DimM>
+  typename SolverEigen<SolverType, DimS, DimM>::SolvVectorOut
+  SolverEigen<SolverType, DimS, DimM>::solve(const SolvVectorInC rhs, SolvVectorIn x_0) {
+    auto & this_solver = static_cast<SolverType&> (*this);
+    SolvVectorOut retval = this->solver.solveWithGuess(rhs, x_0);
+    this->counter += this->solver.iterations();
+    if (this->solver.info() != Eigen::Success) {
+      std::stringstream err {};
+      err << this_solver.name() << " has not converged,"
+          << " After " << this->solver.iterations() << " steps, the solver "
+          << " FAILED with  |r|/|b| = "
+          << std::setw(15) << this->solver.error()
+          << ", cg_tol = " << this->tol << std::endl;
+      throw ConvergenceError(err.str());
+    }
+    if (this->verbose) {
+      std::cout << " After " << this->solver.iterations() << " "
+                << this_solver.name() << " steps, |r|/|b| = "
+                << std::setw(15) << this->solver.error()
+                << ", cg_tol = " << this->tol << std::endl;
+    }
+    return retval;
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <class SolverType, Dim_t DimS, Dim_t DimM>
+  typename SolverEigen<SolverType, DimS, DimM>::Tg_req_t
+  SolverEigen<SolverType, DimS, DimM>::get_tangent_req() const {
+    return tangent_requirement;
+  }
+
+  template class SolverEigen<SolverCGEigen<twoD, twoD>, twoD, twoD>;
+  template class SolverEigen<SolverCGEigen<threeD, threeD>, threeD, threeD>;
+  template class SolverCGEigen<twoD, twoD>;
+  template class SolverCGEigen<threeD, threeD>;
+
+  template class SolverEigen<SolverGMRESEigen<twoD, twoD>, twoD, twoD>;
+  template class SolverEigen<SolverGMRESEigen<threeD, threeD>, threeD, threeD>;
+  template class SolverGMRESEigen<twoD, twoD>;
+  template class SolverGMRESEigen<threeD, threeD>;
+
+  template class SolverEigen<SolverBiCGSTABEigen<twoD, twoD>, twoD, twoD>;
+  template class SolverEigen<SolverBiCGSTABEigen<threeD, threeD>, threeD, threeD>;
+  template class SolverBiCGSTABEigen<twoD, twoD>;
+  template class SolverBiCGSTABEigen<threeD, threeD>;
+
+  template class SolverEigen<SolverDGMRESEigen<twoD, twoD>, twoD, twoD>;
+  template class SolverEigen<SolverDGMRESEigen<threeD, threeD>, threeD, threeD>;
+  template class SolverDGMRESEigen<twoD, twoD>;
+  template class SolverDGMRESEigen<threeD, threeD>;
+
+  template class SolverEigen<SolverMINRESEigen<twoD, twoD>, twoD, twoD>;
+  template class SolverEigen<SolverMINRESEigen<threeD, threeD>, threeD, threeD>;
+  template class SolverMINRESEigen<twoD, twoD>;
+  template class SolverMINRESEigen<threeD, threeD>;
+
+} // muSpectre
diff --git a/src/solver/solver_cg_eigen.hh b/src/solver/solver_cg_eigen.hh
new file mode 100644
index 0000000..7b8cc73
--- /dev/null
+++ b/src/solver/solver_cg_eigen.hh
@@ -0,0 +1,203 @@
+/**
+ * file   solver_cg_eigen.hh
+ *
+ * @author Till Junge <till.junge@epfl.ch>
+ *
+ * @date   19 Jan 2018
+ *
+ * @brief  binding to Eigen's conjugate gradient solver
+ *
+ * @section LICENCE
+ *
+ * Copyright © 2018 Till Junge
+ *
+ * µSpectre is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3, or (at
+ * your option) any later version.
+ *
+ * µSpectre is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Emacs; see the file COPYING. If not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+#ifndef SOLVER_EIGEN_H
+#define SOLVER_EIGEN_H
+
+#include "solver/solver_base.hh"
+
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/IterativeSolvers>
+
+#include <type_traits>
+
+namespace muSpectre {
+  template <class SolverType, Dim_t DimS, Dim_t DimM>
+  class SolverEigen;
+
+  template <Dim_t DimS, Dim_t DimM=DimS>
+  class SolverCGEigen;
+
+  template <Dim_t DimS, Dim_t DimM=DimS>
+  class SolverGMRESEigen;
+
+  template <Dim_t DimS, Dim_t DimM=DimS>
+  class SolverBiCGSTABEigen;
+
+  template <Dim_t DimS, Dim_t DimM=DimS>
+  class SolverDGMRESEigen;
+
+  template <Dim_t DimS, Dim_t DimM=DimS>
+  class SolverMINRESEigen;
+
+  namespace internal {
+
+    template <class Solver>
+    struct Solver_traits {
+    };
+
+
+    template <Dim_t DimS, Dim_t DimM>
+    struct Solver_traits<SolverCGEigen<DimS, DimM>> {
+      using Solver =
+        Eigen::ConjugateGradient<typename SolverEigen<SolverCGEigen<DimS, DimM>,
+                                                      DimS, DimM>::Adaptor,
+                                 Eigen::Lower|Eigen::Upper,
+                                 Eigen::IdentityPreconditioner>;
+    };
+
+    template <Dim_t DimS, Dim_t DimM>
+    struct Solver_traits<SolverGMRESEigen<DimS, DimM>> {
+      using Solver =
+        Eigen::GMRES<typename SolverEigen<SolverGMRESEigen<DimS, DimM>,
+                                          DimS, DimM>::Adaptor,
+                     Eigen::IdentityPreconditioner>;
+    };
+
+    template <Dim_t DimS, Dim_t DimM>
+    struct Solver_traits<SolverBiCGSTABEigen<DimS, DimM>> {
+      using Solver =
+        Eigen::BiCGSTAB<typename SolverEigen<SolverBiCGSTABEigen<DimS, DimM>,
+                                             DimS, DimM>::Adaptor,
+                        Eigen::IdentityPreconditioner>;
+    };
+
+    template <Dim_t DimS, Dim_t DimM>
+    struct Solver_traits<SolverDGMRESEigen<DimS, DimM>> {
+      using Solver =
+        Eigen::DGMRES<typename SolverEigen<SolverDGMRESEigen<DimS, DimM>,
+                                           DimS, DimM>::Adaptor,
+                      Eigen::IdentityPreconditioner>;
+    };
+
+    template <Dim_t DimS, Dim_t DimM>
+    struct Solver_traits<SolverMINRESEigen<DimS, DimM>> {
+      using Solver =
+        Eigen::MINRES<typename SolverEigen<SolverMINRESEigen<DimS, DimM>,
+                                           DimS, DimM>::Adaptor,
+                      Eigen::Lower|Eigen::Upper,
+                      Eigen::IdentityPreconditioner>;
+    };
+
+  }  // internal
+
+  template <class SolverType, Dim_t DimS, Dim_t DimM=DimS>
+  class SolverEigen: public SolverBase<DimS, DimM>
+  {
+  public:
+    using Parent = SolverBase<DimS, DimM>;
+    using SolvVectorIn = typename Parent::SolvVectorIn;
+    using SolvVectorInC = typename Parent::SolvVectorInC;
+    using SolvVectorOut = typename Parent::SolvVectorOut;
+    using Sys_t = typename Parent::Sys_t;
+    using Ccoord = typename Parent::Ccoord;
+    using Tg_req_t = typename Parent::TangentRequirement;
+    using Adaptor = typename Sys_t::Adaptor;
+    using Solver = typename internal::Solver_traits<SolverType>::Solver;
+
+    constexpr static Tg_req_t tangent_requirement{Tg_req_t::NeedEffect};
+
+        //! Default constructor
+    SolverEigen() = delete;
+
+    //! Constructor with domain resolutions, etc,
+    SolverEigen(Sys_t& sys, Real tol, Uint maxiter=0, bool verbose =false);
+
+    //! Copy constructor
+    SolverEigen(const SolverEigen &other) = delete;
+
+    //! Move constructor
+    SolverEigen(SolverEigen &&other) = default;
+
+    //! Destructor
+    virtual ~SolverEigen() = default;
+
+    //! Copy assignment operator
+    SolverEigen& operator=(const SolverEigen &other) = delete;
+
+    //! Move assignment operator
+    SolverEigen& operator=(SolverEigen &&other) = default;
+
+    bool has_converged() const override final {return this->solver.info() == Eigen::Success;}
+
+    void initialise() override final;
+
+    SolvVectorOut solve(const SolvVectorInC rhs, SolvVectorIn x_0) override final;
+
+
+  protected:
+    Tg_req_t get_tangent_req() const override final;
+    Adaptor adaptor;
+    Solver solver;
+
+  };
+
+  template <Dim_t DimS, Dim_t DimM>
+  class SolverCGEigen:
+    public SolverEigen<SolverCGEigen<DimS, DimM>, DimS, DimM> {
+  public:
+    using SolverEigen<SolverCGEigen<DimS, DimM>, DimS, DimM>::SolverEigen;
+    std::string name() const override final {return "CG";}
+  };
+
+  template <Dim_t DimS, Dim_t DimM>
+  class SolverGMRESEigen:
+    public SolverEigen<SolverGMRESEigen<DimS, DimM>, DimS, DimM> {
+  public:
+    using SolverEigen<SolverGMRESEigen<DimS, DimM>, DimS, DimM>::SolverEigen;
+    std::string name() const override final {return "GMRES";}
+  };
+
+  template <Dim_t DimS, Dim_t DimM>
+  class SolverBiCGSTABEigen:
+    public SolverEigen<SolverBiCGSTABEigen<DimS, DimM>, DimS, DimM> {
+  public:
+    using SolverEigen<SolverBiCGSTABEigen<DimS, DimM>, DimS, DimM>::SolverEigen;
+    std::string name() const override final {return "BiCGSTAB";}
+  };
+
+  template <Dim_t DimS, Dim_t DimM>
+  class SolverDGMRESEigen:
+    public SolverEigen<SolverDGMRESEigen<DimS, DimM>, DimS, DimM> {
+  public:
+    using SolverEigen<SolverDGMRESEigen<DimS, DimM>, DimS, DimM>::SolverEigen;
+    std::string name() const override final {return "DGMRES";}
+  };
+
+  template <Dim_t DimS, Dim_t DimM>
+  class SolverMINRESEigen:
+    public SolverEigen<SolverMINRESEigen<DimS, DimM>, DimS, DimM> {
+  public:
+    using SolverEigen<SolverMINRESEigen<DimS, DimM>, DimS, DimM>::SolverEigen;
+    std::string name() const override final {return "MINRES";}
+  };
+
+} // muSpectre
+
+#endif /* SOLVER_EIGEN_H */
diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc
index 59e9b4f..fd90dcf 100644
--- a/src/solver/solvers.cc
+++ b/src/solver/solvers.cc
@@ -1,397 +1,376 @@
 /**
  * file   solvers.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief  implementation of solver functions
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "solver/solvers.hh"
 #include "solver/solver_cg.hh"
 #include "common/iterators.hh"
 
 #include <Eigen/IterativeLinearSolvers>
 
 #include <iomanip>
 #include <cmath>
 
 namespace muSpectre {
 
   template <Dim_t DimS, Dim_t DimM>
   std::vector<OptimizeResult>
   de_geus (SystemBase<DimS, DimM> & sys, const GradIncrements<DimM> & delFs,
-           const Real cg_tol, const Real newton_tol, Uint maxiter,
+           SolverBase<DimS, DimM> & solver, Real newton_tol,
            Dim_t verbose) {
 
     using Field_t = typename MaterialBase<DimS, DimM>::StrainField_t;
     auto solver_fields{std::make_unique<GlobalFieldCollection<DimS, DimM>>()};
     solver_fields->initialise(sys.get_resolutions());
 
     // Corresponds to symbol δF or δε
     auto & incrF{make_field<Field_t>("δF", *solver_fields)};
 
     // Corresponds to symbol ΔF or Δε
     auto & DeltaF{make_field<Field_t>("ΔF", *solver_fields)};
 
     // field to store the rhs for cg calculations
     auto & rhs{make_field<Field_t>("rhs", *solver_fields)};
 
-    // SolverCG<DimS, DimM> cg(sys.get_resolutions(),
-    //                         cg_tol, maxiter, verbose-1>0);
-    // cg.initialise();
-    Eigen::ConjugateGradient<typename SystemBase<DimS, DimM>::Adaptor, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg{};
-    cg.setTolerance(cg_tol);
-    cg.setMaxIterations(maxiter);
-    auto adaptor = sys.get_adaptor();
-    cg.compute(adaptor);
+    solver.initialise();
 
 
-    if (maxiter == 0) {
-      maxiter = sys.size()*DimM*DimM*10;
+    if (solver.get_maxiter() == 0) {
+      solver.set_maxiter(sys.size()*DimM*DimM*10);
     }
 
     size_t count_width{};
     const auto form{sys.get_formulation()};
     std::string strain_symb{};
     if (verbose > 0) {
       //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
-      std::cout << "de Geus for ";
+      std::cout << "de Geus-" << solver.name() << " for ";
       switch (form) {
       case Formulation::small_strain: {
         strain_symb = "ε";
         std::cout << "small";
         break;
       }
       case Formulation::finite_strain: {
         strain_symb = "F";
         std::cout << "finite";
         break;
       }
       default:
         throw SolverError("unknown formulation");
         break;
       }
       std::cout << " strain with" << std::endl
                 << "newton_tol = " << newton_tol << ", cg_tol = "
-                << cg_tol << " maxiter = " << maxiter << " and Δ"
+                << solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ"
                 << strain_symb << " =" <<std::endl;
       for (auto&& tup: akantu::enumerate(delFs)) {
         auto && counter{std::get<0>(tup)};
         auto && grad{std::get<1>(tup)};
         std::cout << "Step " << counter + 1 << ":" << std::endl
                   << grad << std::endl;
       }
-      count_width = size_t(std::log10(maxiter))+1;
+      count_width = size_t(std::log10(solver.get_maxiter()))+1;
     }
 
     // initialise F = I or ε = 0
     auto & F{sys.get_strain()};
     switch (form) {
     case Formulation::finite_strain: {
       F.get_map() = Matrices::I2<DimM>();
       break;
     }
     case Formulation::small_strain: {
       F.get_map() = Matrices::I2<DimM>().Zero();
       for (const auto & delF: delFs) {
         if (!check_symmetry(delF)) {
           throw SolverError("all Δε must be symmetric!");
         }
       }
       break;
     }
     default:
       throw SolverError("Unknown formulation");
       break;
     }
 
     // initialise return value
     std::vector<OptimizeResult> ret_val{};
 
     // initialise materials
     constexpr bool need_tangent{true};
     sys.initialise_materials(need_tangent);
 
     Grad_t<DimM> previous_grad{Grad_t<DimM>::Zero()};
-    Uint cg_counter{0};
     for (const auto & delF: delFs) { //incremental loop
 
       Real incrNorm{2*newton_tol}, gradNorm{1};
       auto convergence_test = [&incrNorm, &gradNorm, &newton_tol] () {
         return incrNorm/gradNorm <= newton_tol;
       };
       Uint newt_iter{0};
       for (;
-           (newt_iter < maxiter) && (!convergence_test() ||
+           (newt_iter < solver.get_maxiter()) && (!convergence_test() ||
                                      (newt_iter==1));
            ++newt_iter) {
 
         // obtain material response
         auto res_tup{sys.evaluate_stress_tangent(F)};
         auto & P{std::get<0>(res_tup)};
         auto & K{std::get<1>(res_tup)};
 
         auto tangent_effect = [&sys, &K] (const Field_t & dF, Field_t & dP) {
           sys.directional_stiffness(K, dF, dP);
         };
 
 
         if (newt_iter == 0) {
           DeltaF.get_map() = -(delF-previous_grad); // neg sign because rhs
           tangent_effect(DeltaF, rhs);
-          // cg.solve(tangent_effect, rhs, incrF);
-          incrF.eigenvec() = cg.solve(rhs.eigenvec());
+          incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
           F.eigen() -= DeltaF.eigen();
         } else {
           rhs.eigen() = -P.eigen();
           sys.project(rhs);
-          //cg.solve(tangent_effect, rhs, incrF);
-          incrF.eigenvec() = cg.solve(rhs.eigenvec());
-        }
-        if (cg.info() != Eigen::Success) {
-          throw SolverError("cg did not converge");
+          incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
         }
 
-        cg_counter += cg.iterations();
         F.eigen() += incrF.eigen();
 
         incrNorm = incrF.eigen().matrix().norm();
         gradNorm = F.eigen().matrix().norm();
         if (verbose>0) {
           std::cout << "at Newton step " << std::setw(count_width) << newt_iter
                     << ", |δ" << strain_symb << "|/|Δ" << strain_symb
                     << "| = " << std::setw(17) << incrNorm/gradNorm
                     << ", tol = " << newton_tol << std::endl;
           if (verbose-1>1) {
             std::cout << "<" << strain_symb << "> =" << std::endl
                       << F.get_map().mean() << std::endl;
           }
         }
       }
       // update previous gradient
       previous_grad = delF;
 
       ret_val.push_back(OptimizeResult{F.eigen(), sys.get_stress().eigen(),
             convergence_test(), Int(convergence_test()),
             "message not yet implemented",
-            newt_iter, cg_counter});
+            newt_iter, solver.get_counter()});
 
 
       //!store history variables here
 
     }
 
     return ret_val;
 
   }
 
   template std::vector<OptimizeResult>
   de_geus (SystemBase<twoD, twoD> & sys, const GradIncrements<twoD>& delF0,
-           const Real cg_tol, const Real newton_tol, Uint maxiter,
+           SolverBase<twoD, twoD> & solver, Real newton_tol,
            Dim_t verbose);
 
   // template typename SystemBase<twoD, threeD>::StrainField_t &
   // de_geus (SystemBase<twoD, threeD> & sys, const GradIncrements<threeD>& delF0,
   //            const Real cg_tol, const Real newton_tol, Uint maxiter,
   //            Dim_t verbose);
 
   template std::vector<OptimizeResult>
   de_geus (SystemBase<threeD, threeD> & sys, const GradIncrements<threeD>& delF0,
-           const Real cg_tol, const Real newton_tol, Uint maxiter,
+           SolverBase<threeD, threeD> & solver, Real newton_tol,
            Dim_t verbose);
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   std::vector<OptimizeResult>
   newton_cg (SystemBase<DimS, DimM> & sys, const GradIncrements<DimM> & delFs,
-             const Real cg_tol, const Real newton_tol, Uint maxiter,
+             SolverBase<DimS, DimM> & solver, Real newton_tol,
              Dim_t verbose) {
     using Field_t = typename MaterialBase<DimS, DimM>::StrainField_t;
     auto solver_fields{std::make_unique<GlobalFieldCollection<DimS, DimM>>()};
     solver_fields->initialise(sys.get_resolutions());
 
     // Corresponds to symbol δF or δε
     auto & incrF{make_field<Field_t>("δF", *solver_fields)};
 
     // field to store the rhs for cg calculations
     auto & rhs{make_field<Field_t>("rhs", *solver_fields)};
 
-    SolverCG<DimS, DimM> cg(sys.get_resolutions(),
-                            cg_tol, maxiter, verbose-1>0);
-    cg.initialise();
-
+    solver.initialise();
 
-    if (maxiter == 0) {
-      maxiter = sys.size()*DimM*DimM*10;
+    if (solver.get_maxiter() == 0) {
+      solver.set_maxiter(sys.size()*DimM*DimM*10);
     }
 
     size_t count_width{};
     const auto form{sys.get_formulation()};
     std::string strain_symb{};
     if (verbose > 0) {
       //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
-      std::cout << "Newton-CG for ";
+      std::cout << "Newton-" << solver.name() << " for ";
       switch (form) {
       case Formulation::small_strain: {
         strain_symb = "ε";
         std::cout << "small";
         break;
       }
       case Formulation::finite_strain: {
         strain_symb = "Fy";
         std::cout << "finite";
         break;
       }
       default:
         throw SolverError("unknown formulation");
         break;
       }
       std::cout << " strain with" << std::endl
                 << "newton_tol = " << newton_tol << ", cg_tol = "
-                << cg_tol << " maxiter = " << maxiter << " and Δ"
+                << solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ"
                 << strain_symb << " =" <<std::endl;
       for (auto&& tup: akantu::enumerate(delFs)) {
         auto && counter{std::get<0>(tup)};
         auto && grad{std::get<1>(tup)};
         std::cout << "Step " << counter + 1 << ":" << std::endl
                   << grad << std::endl;
       }
-      count_width = size_t(std::log10(maxiter))+1;
+      count_width = size_t(std::log10(solver.get_maxiter()))+1;
     }
 
     // initialise F = I or ε = 0
     auto & F{sys.get_strain()};
     switch (sys.get_formulation()) {
     case Formulation::finite_strain: {
       F.get_map() = Matrices::I2<DimM>();
       break;
     }
     case Formulation::small_strain: {
       F.get_map() = Matrices::I2<DimM>().Zero();
       for (const auto & delF: delFs) {
         if (!check_symmetry(delF)) {
           throw SolverError("all Δε must be symmetric!");
         }
       }
       break;
     }
     default:
       throw SolverError("Unknown formulation");
       break;
     }
 
     // initialise return value
     std::vector<OptimizeResult> ret_val{};
 
     // initialise materials
     constexpr bool need_tangent{true};
     sys.initialise_materials(need_tangent);
 
     Grad_t<DimM> previous_grad{Grad_t<DimM>::Zero()};
     for (const auto & delF: delFs) { //incremental loop
       // apply macroscopic strain increment
       for (auto && grad: F.get_map()) {
         grad += delF - previous_grad;
       }
 
       Real incrNorm{2*newton_tol}, gradNorm{1};
       auto convergence_test = [&incrNorm, &gradNorm, &newton_tol] () {
         return incrNorm/gradNorm <= newton_tol;
       };
       Uint newt_iter{0};
 
       for (;
-           newt_iter < maxiter && !convergence_test();
+           newt_iter < solver.get_maxiter() && !convergence_test();
            ++newt_iter) {
 
         // obtain material response
         auto res_tup{sys.evaluate_stress_tangent(F)};
         auto & P{std::get<0>(res_tup)};
-        auto & K{std::get<1>(res_tup)};
-
-        auto tangent_effect = [&sys, &K] (const Field_t & dF, Field_t & dP) {
-          sys.directional_stiffness(K, dF, dP);
-        };
 
         rhs.eigen() = -P.eigen();
         sys.project(rhs);
-        cg.solve(tangent_effect, rhs, incrF);
+        incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
+
 
         F.eigen() += incrF.eigen();
 
         incrNorm = incrF.eigen().matrix().norm();
         gradNorm = F.eigen().matrix().norm();
         if (verbose > 0) {
           std::cout << "at Newton step " << std::setw(count_width) << newt_iter
                     << ", |δ" << strain_symb << "|/|Δ" << strain_symb
                     << "| = " << std::setw(17) << incrNorm/gradNorm
                     << ", tol = " << newton_tol << std::endl;
 
           if (verbose-1>1) {
             std::cout << "<" << strain_symb << "> =" << std::endl
                       << F.get_map().mean() << std::endl;
           }
         }
 
       }
       // update previous gradient
       previous_grad = delF;
 
       ret_val.push_back(OptimizeResult{F.eigen(), sys.get_stress().eigen(),
             convergence_test(), Int(convergence_test()),
             "message not yet implemented",
-            newt_iter, cg.get_counter()});
+            newt_iter, solver.get_counter()});
 
       //store history variables here
 
     }
 
     return ret_val;
 
   }
 
   template std::vector<OptimizeResult>
   newton_cg (SystemBase<twoD, twoD> & sys, const GradIncrements<twoD>& delF0,
-             const Real cg_tol, const Real newton_tol, Uint maxiter,
+             SolverBase<twoD, twoD> & solver, Real newton_tol,
              Dim_t verbose);
 
   // template typename SystemBase<twoD, threeD>::StrainField_t &
   // newton_cg (SystemBase<twoD, threeD> & sys, const GradIncrements<threeD>& delF0,
   //            const Real cg_tol, const Real newton_tol, Uint maxiter,
   //            Dim_t verbose);
 
   template std::vector<OptimizeResult>
   newton_cg (SystemBase<threeD, threeD> & sys, const GradIncrements<threeD>& delF0,
-             const Real cg_tol, const Real newton_tol, Uint maxiter,
+             SolverBase<threeD, threeD> & solver, Real newton_tol,
              Dim_t verbose);
 
 
   /* ---------------------------------------------------------------------- */
   bool check_symmetry(const Eigen::Ref<const Eigen::ArrayXXd>& eps,
                       Real rel_tol){
     return rel_tol >= (eps-eps.transpose()).matrix().norm()/eps.matrix().norm();
   }
 
 
 }  // muSpectre
diff --git a/src/solver/solvers.hh b/src/solver/solvers.hh
index bb9f7ac..c9a15e4 100644
--- a/src/solver/solvers.hh
+++ b/src/solver/solvers.hh
@@ -1,116 +1,116 @@
 /**
  * file   solvers.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief  Free functions for solving
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef SOLVERS_H
 #define SOLVERS_H
 
-#include "solver/solver_cg.hh"
+#include "solver/solver_base.hh"
 
 #include <Eigen/Dense>
 
 #include <vector>
 #include <string>
 
 namespace muSpectre {
 
   /**
    * emulates scipy.optimize.OptimizeResult
    */
   struct OptimizeResult
   {
     //! Strain ε or Gradient F at solution
     Eigen::ArrayXXd grad;
     //! Cauchy stress σ or first Piola-Kirchhoff stress P at solution
     Eigen::ArrayXXd stress;
     //! whether or not the solver exited successfully
     bool success;
     //! Termination status of the optimizer. Its value depends on the
     //! underlying solver. Refer to message for details.
     Int status;
     //! Description of the cause of the termination.
     std::string message;
     //! number of iterations
     Uint nb_it;
     //! number of system evaluations
     Uint nb_fev;
   };
 
   template <Dim_t Dim>
   using Grad_t = Matrices::Tens2_t<Dim>;
   template <Dim_t Dim>
   using GradIncrements = std::vector<Grad_t<Dim>,
                                      Eigen::aligned_allocator<Grad_t<Dim>>>;
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM=DimS>
   std::vector<OptimizeResult>
   newton_cg (SystemBase<DimS, DimM> & sys,
              const GradIncrements<DimM> & delF0,
-             const Real cg_tol, const Real newton_tol, const Uint maxiter=0,
+             SolverBase<DimS, DimM> & solver, Real newton_tol,
              Dim_t verbose = 0);
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM=DimS>
   inline OptimizeResult
   newton_cg (SystemBase<DimS, DimM> & sys, const Grad_t<DimM> & delF0,
-             const Real cg_tol, const Real newton_tol, const Uint maxiter=0,
+             SolverBase<DimS, DimM> & solver, Real newton_tol,
              Dim_t verbose = 0){
     return newton_cg(sys, GradIncrements<DimM>{delF0},
-                     cg_tol, newton_tol, maxiter, verbose)[0];
+                     solver, newton_tol, verbose)[0];
   }
 
     /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM=DimS>
   std::vector<OptimizeResult>
   de_geus (SystemBase<DimS, DimM> & sys,
            const GradIncrements<DimM> & delF0,
-           const Real cg_tol, const Real newton_tol, const Uint maxiter=0,
+           SolverBase<DimS, DimM> & solver, Real newton_tol,
            Dim_t verbose = 0);
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM=DimS>
   OptimizeResult
   de_geus (SystemBase<DimS, DimM> & sys, const Grad_t<DimM> & delF0,
-           const Real cg_tol, const Real newton_tol, const Uint maxiter=0,
+           SolverBase<DimS, DimM> & solver, Real newton_tol,
            Dim_t verbose = 0){
     return de_geus(sys, GradIncrements<DimM>{delF0},
-                   cg_tol, newton_tol, maxiter, verbose)[0];
+                   solver, newton_tol, verbose)[0];
   }
 
   /* ---------------------------------------------------------------------- */
   /**
    * check whether a strain is symmetric, for the purposes of small
    * strain problems
    */
   bool check_symmetry(const Eigen::Ref<const Eigen::ArrayXXd>& eps,
                       Real rel_tol = 1e-8);
 
 }  // muSpectre
 
 #endif /* SOLVERS_H */
diff --git a/src/system/system_base.hh b/src/system/system_base.hh
index 7451533..f12e59e 100644
--- a/src/system/system_base.hh
+++ b/src/system/system_base.hh
@@ -1,251 +1,251 @@
 /**
  * file   system_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 Nov 2017
  *
  * @brief Base class representing a unit cell system with single
  *        projection operator
  *
  * @section LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef SYSTEM_BASE_H
 #define SYSTEM_BASE_H
 
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "common/field.hh"
 #include "common/utilities.hh"
 #include "materials/material_base.hh"
 #include "fft/projection_base.hh"
 #include "system/system_traits.hh"
 
 #include <vector>
 #include <memory>
 #include <tuple>
 #include <functional>
 
 namespace muSpectre {
   template <class System>
   class SystemAdaptor;
   //! DimS spatial dimension (dimension of problem
   //! DimM material_dimension (dimension of constitutive law)
-  template <Dim_t DimS, Dim_t DimM>
+  template <Dim_t DimS, Dim_t DimM=DimS>
   class SystemBase
   {
   public:
     using Ccoord = Ccoord_t<DimS>;
     using Rcoord = Rcoord_t<DimS>;
     using FieldCollection_t = FieldCollection<DimS, DimM>;
     using Collection_ptr = std::unique_ptr<FieldCollection_t>;
     using Material_t = MaterialBase<DimS, DimM>;
     using Material_ptr = std::unique_ptr<Material_t>;
     using Projection_t = ProjectionBase<DimS, DimM>;
     using Projection_ptr = std::unique_ptr<Projection_t>;
     using StrainField_t =
       TensorField<FieldCollection_t, Real, secondOrder, DimM>;
     using StressField_t =
       TensorField<FieldCollection_t, Real, secondOrder, DimM>;
     using TangentField_t =
       TensorField<FieldCollection_t, Real, fourthOrder, DimM>;
     using FullResponse_t =
       std::tuple<const StressField_t&, const TangentField_t&>;
     using iterator = typename CcoordOps::Pixels<DimS>::iterator;
     using SolvVectorIn = Eigen::Ref<const Eigen::VectorXd>;
     using SolvVectorOut = Eigen::Map<Eigen::VectorXd>;
     using Adaptor = SystemAdaptor<SystemBase>;
 
     //! Default constructor
     SystemBase() = delete;
 
     //! constructor using sizes and resolution
     SystemBase(Projection_ptr projection);
 
     //! Copy constructor
     SystemBase(const SystemBase &other) = delete;
 
     //! Move constructor
     SystemBase(SystemBase &&other) = default;
 
     //! Destructor
     virtual ~SystemBase() = default;
 
     //! Copy assignment operator
     SystemBase& operator=(const SystemBase &other) = delete;
 
     //! Move assignment operator
     SystemBase& operator=(SystemBase &&other) = default;
 
     /**
      * Materials can only be moved. This is to assure exclusive
      * ownership of any material by this system
      */
     Material_t & add_material(Material_ptr mat);
 
     /**
      * evaluates all materials
      */
     FullResponse_t evaluate_stress_tangent(StrainField_t & F);
 
     /**
      * evaluate directional stiffness (i.e. G:K:δF or G:K:δε)
      */
 
     StressField_t & directional_stiffness(const TangentField_t & K,
                                           const StrainField_t & delF,
                                           StressField_t & delP);
     /**
      * vectorized version for eigen solvers, no copy, but only works
      * when fields have ArrayStore=false
      */
      SolvVectorOut directional_stiffness_vec(const SolvVectorIn & delF);
     /**
      * Evaluate directional stiffness into a temporary array and
      * return a copy. This is a costly and wasteful interface to
      * directional_stiffness and should only be used for debugging or
      * in the python interface
      */
     Eigen::ArrayXXd directional_stiffness_with_copy
       (Eigen::Ref<Eigen::ArrayXXd> delF);
 
     /**
      * Convenience function circumventing the neeed to use the
      * underlying projection
      */
     StressField_t & project(StressField_t & field);
 
     StrainField_t & get_strain();
 
     const StressField_t & get_stress() const;
 
     const TangentField_t & get_tangent(bool create = false);
 
     StrainField_t & get_managed_field(std::string unique_name);
 
     /**
      * general initialisation; initialises the projection and
      * fft_engine (i.e. infrastructure) but not the materials. These
      * need to be initialised separately
      */
     void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
     /**
      * initialise materials (including resetting any history variables)
      */
     void initialise_materials(bool stiffness=false);
 
     iterator begin();
     iterator end();
     size_t size() const {return pixels.size();}
 
     const Ccoord & get_resolutions() const {return this->resolutions;}
     const Rcoord & get_lengths() const {return this->lengths;}
 
     /**
      * formulation is hard set by the choice of the projection class
      */
     const Formulation & get_formulation() const {
       return this->projection->get_formulation();}
 
     bool is_initialised() const {return this->initialised;}
 
     Eigen::Map<Eigen::ArrayXXd> get_projection() {
       return this->projection->get_operator();}
 
     constexpr static Dim_t get_sdim() {return DimS;};
 
     Adaptor get_adaptor();
     Dim_t nb_dof() const {return this->size()*ipow(DimS, 2);};
 
 
   protected:
     //! make sure that every pixel is assigned to one and only one material
     void check_material_coverage();
 
     const Ccoord & resolutions;
     CcoordOps::Pixels<DimS> pixels;
     const Rcoord & lengths;
     Collection_ptr fields;
     StrainField_t & F;
     StressField_t & P;
     //! Tangent field might not even be required; so this is an
     //! optional ref_wrapper instead of a ref
     optional<std::reference_wrapper<TangentField_t>> K{};
     std::vector<Material_ptr> materials{};
     Projection_ptr projection;
     bool initialised{false};
     const Formulation form;
   private:
   };
 
 
   template <class System>
   class SystemAdaptor: public Eigen::EigenBase<SystemAdaptor<System>> {
 
   public:
     using Scalar = double;
     using RealScalar = double;
     using StorageIndex = int;
     enum {
       ColsAtCompileTime = Eigen::Dynamic,
       MaxColsAtCompileTime = Eigen::Dynamic,
       RowsAtCompileTime = Eigen::Dynamic,
       MaxRowsAtCompileTime = Eigen::Dynamic,
       IsRowMajor = false
     };
 
     SystemAdaptor(System & sys):sys{sys}{}
     Eigen::Index rows() const {return this->sys.nb_dof();}
     Eigen::Index cols() const {return this->rows();}
 
     template<typename Rhs>
     Eigen::Product<SystemAdaptor,Rhs,Eigen::AliasFreeProduct>
     operator*(const Eigen::MatrixBase<Rhs>& x) const {
       return Eigen::Product<SystemAdaptor,Rhs,Eigen::AliasFreeProduct>
         (*this, x.derived());
     }
     System & sys;
   };
 
 }  // muSpectre
 
 // Implementation of SystemAdaptor * Eigen::DenseVector though a specialization of internal::generic_product_impl:
 namespace Eigen {
   namespace internal {
     template<typename Rhs, class SystemAdaptor>
     struct generic_product_impl<SystemAdaptor, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
       : generic_product_impl_base<SystemAdaptor,Rhs,generic_product_impl<SystemAdaptor,Rhs> >
     {
       typedef typename Product<SystemAdaptor,Rhs>::Scalar Scalar;
       template<typename Dest>
       static void scaleAndAddTo(Dest& dst, const SystemAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/)
       {
         // This method should implement "dst += alpha * lhs * rhs" inplace,
         // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
         // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
         dst.noalias() += const_cast<SystemAdaptor&>(lhs).sys.directional_stiffness_vec(rhs);
       }
     };
   }
 }
 
 
 #endif /* SYSTEM_BASE_H */
diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py
index 4b3a5f5..24c6347 100755
--- a/tests/python_binding_tests.py
+++ b/tests/python_binding_tests.py
@@ -1,106 +1,107 @@
 #!/usr/bin/env python3
 """
 file   python_binding_tests.py
 
 @author Till Junge <till.junge@epfl.ch>
 
 @date   09 Jan 2018
 
 @brief  Unit tests for python bindings
 
 @section LICENCE
 
 Copyright © 2018 Till Junge
 
 µSpectre is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License as
 published by the Free Software Foundation, either version 3, or (at
 your option) any later version.
 
 µSpectre is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with GNU Emacs; see the file COPYING. If not, write to the
 Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 Boston, MA 02111-1307, USA.
 """
 
 import unittest
 import numpy as np
 
 from python_test_imports import µ
 
 from python_fft_tests import FFT_Check
 from python_projection_tests import *
 
 class SystemCheck(unittest.TestCase):
     def test_Construction(self):
         """
         Simple check for system constructors
         """
         resolution = [5,7]
         lengths = [5.2, 8.3]
         formulation = µ.Formulation.small_strain
         try:
             sys = µ.SystemFactory(resolution,
                                   lengths,
                                   formulation)
             mat = µ.material.MaterialHooke2d.make(sys, "material", 210e9, .33)
         except Exception as err:
             print(err)
             raise err
 
 class MaterialHooke2dCheck(unittest.TestCase):
     def setUp(self):
         self.resolution = [5,7]
         self.lengths = [5.2, 8.3]
         self.formulation = µ.Formulation.small_strain
         self.sys = µ.SystemFactory(self.resolution,
                                    self.lengths,
                                    self.formulation)
         self.mat = µ.material.MaterialHooke2d.make(
             self.sys, "material", 210e9, .33)
 
     def test_add_material(self):
         self.mat.add_pixel([2,1])
 
 
 class SolverCheck(unittest.TestCase):
     def setUp(self):
         self.resolution = [3, 3]#[5,7]
         self.lengths = [3., 3.]#[5.2, 8.3]
         self.formulation = µ.Formulation.finite_strain
         self.sys = µ.SystemFactory(self.resolution,
                                    self.lengths,
                                    self.formulation)
         self.hard = µ.material.MaterialHooke2d.make(
             self.sys, "hard", 210e9, .33)
         self.soft = µ.material.MaterialHooke2d.make(
             self.sys, "soft",  70e9, .33)
 
     def test_solve(self):
         for i, pixel in enumerate(self.sys):
             if i < 3:
                 self.hard.add_pixel(pixel)
             else:
                 self.soft.add_pixel(pixel)
 
         self.sys.initialise()
         tol = 1e-6
         Del0 = np.array([[0, .1],
                          [0,  0]])
         maxiter = 100
         verbose = 0
         # the following segfaults:
+        solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose)
         r = µ.solvers.de_geus(self.sys, Del0,
-                              tol, tol, maxiter, verbose)
+                              solver,tol, verbose)
         #print(r)
 
 
 
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc
index 20a6485..c97f751 100644
--- a/tests/test_solver_newton_cg.cc
+++ b/tests/test_solver_newton_cg.cc
@@ -1,175 +1,196 @@
 /**
  * file   test_solver_newton_cg.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief  Tests for the standard Newton-Raphson + Conjugate Gradient solver
  *
  * @section LICENCE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "tests.hh"
 #include "solver/solvers.hh"
+#include "solver/solver_cg.hh"
+#include "solver/solver_cg_eigen.hh"
 #include "fft/fftw_engine.hh"
 #include "fft/projection_finite_strain_fast.hh"
 #include "materials/material_hyper_elastic1.hh"
 #include "common/iterators.hh"
 #include "common/ccoord_operations.hh"
 #include "system/system_factory.hh"
 
 
 namespace muSpectre {
 
   BOOST_AUTO_TEST_SUITE(newton_cg_tests);
 
   BOOST_AUTO_TEST_CASE(manual_construction_test) {
     // constexpr Dim_t dim{twoD};
     constexpr Dim_t dim{threeD};
 
     // constexpr Ccoord_t<dim> resolutions{3, 3};
     // constexpr Rcoord_t<dim> lengths{2.3, 2.7};
     constexpr Ccoord_t<dim> resolutions{5, 5, 5};
     constexpr Rcoord_t<dim> lengths{5, 5, 5};
     auto fft_ptr{std::make_unique<FFTW_Engine<dim, dim>>(resolutions, lengths)};
     auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr))};
     SystemBase<dim, dim> sys(std::move(proj_ptr));
 
     using Mat_t = MaterialHyperElastic1<dim, dim>;
     //const Real Young{210e9}, Poisson{.33};
     const Real Young{1.0030648180242636}, Poisson{0.29930675909878679};
     // const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))};
     // const Real mu{Young/(2*(1+Poisson))};
 
     auto& Material_hard = Mat_t::make(sys, "hard", 10*Young, Poisson);
     auto& Material_soft = Mat_t::make(sys, "soft",    Young, Poisson);
 
     for (auto && tup: akantu::enumerate(sys)) {
       auto && pixel = std::get<1>(tup);
       if (std::get<0>(tup) == 0) {
         Material_hard.add_pixel(pixel);
       } else {
         Material_soft.add_pixel(pixel);
       }
     }
     sys.initialise();
 
     Grad_t<dim> delF0;
     delF0 << 0, 1., 0, 0, 0, 0, 0, 0, 0;
     constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
     constexpr Uint maxiter{CcoordOps::get_size(resolutions)*ipow(dim, secondOrder)*10};
     constexpr bool verbose{false};
 
     GradIncrements<dim> grads; grads.push_back(delF0);
-    Eigen::ArrayXXd res1{de_geus(sys, grads, cg_tol, newton_tol, maxiter, verbose)[0].grad};
+    SolverCG<dim> cg{sys, cg_tol, maxiter, bool(verbose)};
+    Eigen::ArrayXXd res1{de_geus(sys, grads, cg, newton_tol, verbose)[0].grad};
 
-    Eigen::ArrayXXd res2{newton_cg(sys, grads, cg_tol, newton_tol, maxiter, verbose)[0].grad};
+    SolverCG<dim> cg2{sys, cg_tol, maxiter, bool(verbose)};
+    Eigen::ArrayXXd res2{newton_cg(sys, grads, cg2, newton_tol, verbose)[0].grad};
     BOOST_CHECK_LE(abs(res1-res2).mean(), cg_tol);
   }
 
   BOOST_AUTO_TEST_CASE(small_strain_patch_test) {
     constexpr Dim_t dim{twoD};
     using Ccoord = Ccoord_t<dim>;
     using Rcoord = Rcoord_t<dim>;
     constexpr Ccoord resolutions{CcoordOps::get_cube<dim>(3)};
     constexpr Rcoord lengths{CcoordOps::get_cube<dim>(1.)};
     constexpr Formulation form{Formulation::small_strain};
 
     // number of layers in the hard material
     constexpr Uint nb_lays{1};
     constexpr Real contrast{2};
     static_assert(nb_lays < resolutions[0],
                   "the number or layers in the hard material must be smaller "
                   "than the total number of layers in dimension 0");
 
     auto sys{make_system(resolutions, lengths, form)};
 
     using Mat_t = MaterialHyperElastic1<dim, dim>;
     constexpr Real Young{2.}, Poisson{.33};
     auto material_hard{std::make_unique<Mat_t>("hard", contrast*Young, Poisson)};
     auto material_soft{std::make_unique<Mat_t>("soft",          Young, Poisson)};
 
     for (const auto & pixel: sys) {
       if (pixel[0] < Dim_t(nb_lays)) {
         material_hard->add_pixel(pixel);
       } else {
         material_soft->add_pixel(pixel);
       }
     }
 
     sys.add_material(std::move(material_hard));
     sys.add_material(std::move(material_soft));
     sys.initialise();
 
     Grad_t<dim> delEps0{Grad_t<dim>::Zero()};
     constexpr Real eps0 = 1.;
-    delEps0(0, 1) = delEps0(1, 0) = eps0;
+    //delEps0(0, 1) = delEps0(1, 0) = eps0;
+    delEps0(0, 0) = eps0;
 
     constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
     constexpr Uint maxiter{dim*10};
-    constexpr Dim_t verbose{3};
+    constexpr Dim_t verbose{0};
 
-    auto result = de_geus(sys, delEps0, cg_tol, newton_tol, maxiter, verbose);
+    SolverCGEigen<dim> cg{sys, cg_tol, maxiter, bool(verbose)};
+    auto result = de_geus(sys, delEps0, cg, newton_tol, verbose);
     if (verbose) {
       std::cout << "result:" << std::endl << result.grad << std::endl;
       std::cout << "mean strain = " << std::endl
                 << sys.get_strain().get_map().mean() << std::endl;
     }
 
     /**
      *  verification of resultant strains: subscript ₕ for hard and ₛ
      *  for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast
      *
      *     Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ
      *  => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ
      *
      *  σ is constant across all layers
      *        σₕ = σₛ
      *  => Eₕ εₕ = Eₛ εₛ
      *  => εₕ = 1/k εₛ
      *  => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ
      */
     constexpr Real factor{1/contrast * Real(nb_lays)/resolutions[0]
         + 1.-nb_lays/Real(resolutions[0])};
     constexpr Real eps_soft{eps0/factor};
     constexpr Real eps_hard{eps_soft/contrast};
     if (verbose) {
       std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl;
       std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl;
     }
     Grad_t<dim> Eps_hard; Eps_hard << eps_hard, 0, 0, 0;
     Grad_t<dim> Eps_soft; Eps_soft << eps_soft, 0, 0, 0;
 
-    //TODO small strain projection incorrect
-    // for (const auto & pixel: sys) {
-    //   if (pixel[0] < Dim_t(nb_lays)) {
-    //     BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
-    //   } else {
-    //     BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
-    //   }
-    // }
+    // verify uniaxial tension patch test
+    for (const auto & pixel: sys) {
+      if (pixel[0] < Dim_t(nb_lays)) {
+        BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
+      } else {
+        BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
+      }
+    }
 
+    delEps0 = Grad_t<dim>::Zero();
+    delEps0(0, 1) = delEps0(1, 0) = eps0;
 
+    SolverCG<dim> cg2{sys, cg_tol, maxiter, bool(verbose)};
+    result = newton_cg(sys, delEps0, cg2, newton_tol, verbose);
+    Eps_hard << 0, eps_hard, eps_hard, 0;
+    Eps_soft << 0, eps_soft, eps_soft, 0;
+
+    // verify pure shear patch test
+    for (const auto & pixel: sys) {
+      if (pixel[0] < Dim_t(nb_lays)) {
+        BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
+      } else {
+        BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
+      }
+    }
   }
 
   BOOST_AUTO_TEST_SUITE_END();
 
 }  // muSpectre