Index: ThirdParty/vmg/test/unit_test/Makefile.am
===================================================================
--- ThirdParty/vmg/test/unit_test/Makefile.am	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/Makefile.am	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,163 @@
+if COND_HAVE_BOOST_UNIT_TEST
+
+AM_CPPFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/test/unit_test -I$(top_srcdir)/test/interfaces $(VTK_CXXFLAGS) $(BOOST_CPPFLAGS)
+AM_LDFLAGS = $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS)
+TESTLIBS = \
+	$(BOOST_UNIT_TEST_FRAMEWORK_LIBS) \
+	$(LIBVMG) \
+	libinterfaces.a \
+	$(VTK_LDFLAGS) \
+	$(VTK_LIBS)
+
+if HAVE_LAPACK
+TESTLIBS += \
+	$(LAPACK_LIBS) \
+	$(BLAS_LIBS) \
+	$(FCLIBS)
+endif
+
+noinst_LTLIBRARIES = libvmg_testhelper.la
+libvmg_testhelper_la_SOURCES = \
+	unit_test/unit_test/smoother_fixture.cpp \
+	unit_test/unit_test/smoother_fixture.hpp
+libvmg_testhelper_la_CPPFLAGS = $(AM_CPPFLAGS)
+libvmg_testhelper_la_LDFLAGS = $(AM_LDFLAGS)
+
+BOOSTTESTS = \
+	InterpolateSuite \
+	InterpolateFieldSuite \
+	LibraryDirichletFASLRTest \
+	LibraryDirichletCSTest \
+	LibraryPeriodicCSTest \
+	LibraryPeriodicFASTest \
+	LibraryForceSinusTest \
+	LibraryDirichletCSMPITest \
+	LibraryPeriodicFASMPITest \
+	LibraryPeriodicCSMPITest \
+	LibraryDirichletFASTest \
+	LibraryDirichletFASLRMPITest \
+	LibraryOpenFASMPITest \
+	LibraryDirichletFASMPITest \
+	SmootherGSLexTest \
+	SmootherGSRBTest
+
+NONFUNCTIONINGTESTS = \
+	LibraryForceSinusUnitTest
+
+noinst_PROGRAMS += $(BOOSTTESTS) \
+	LibraryForceSinusUnitTest
+
+check_PROGRAMS += $(BOOSTTESTS) \
+	LibraryForceSinusUnitTest
+
+TESTS += $(BOOSTTESTS)
+
+LibraryDirichletFASLRTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/dirichlet_fas_lr.cpp
+LibraryDirichletFASLRTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryDirichletFASLRTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryDirichletFASLRTest_LDADD = $(TESTLIBS)
+LibraryDirichletCSTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/dirichlet_cs.cpp
+LibraryDirichletCSTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryDirichletCSTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryDirichletCSTest_LDADD = $(TESTLIBS)
+LibraryPeriodicCSTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/periodic_cs.cpp
+LibraryPeriodicCSTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryPeriodicCSTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryPeriodicCSTest_LDADD = $(TESTLIBS)
+LibraryPeriodicFASTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/periodic_fas.cpp
+LibraryPeriodicFASTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryPeriodicFASTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryPeriodicFASTest_LDADD = $(TESTLIBS)
+LibraryForceSinusTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/force_sinus.cpp
+LibraryForceSinusTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryForceSinusTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryForceSinusTest_LDADD = $(TESTLIBS)
+LibraryDirichletCSMPITest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/dirichlet_cs_mpi.cpp
+LibraryDirichletCSMPITest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryDirichletCSMPITest_LDFLAGS = $(AM_LDFLAGS)
+LibraryDirichletCSMPITest_LDADD = $(TESTLIBS)
+LibraryPeriodicFASMPITest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/periodic_fas_mpi.cpp
+LibraryPeriodicFASMPITest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryPeriodicFASMPITest_LDFLAGS = $(AM_LDFLAGS)
+LibraryPeriodicFASMPITest_LDADD = $(TESTLIBS)
+LibraryPeriodicCSMPITest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/periodic_cs_mpi.cpp
+LibraryPeriodicCSMPITest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryPeriodicCSMPITest_LDFLAGS = $(AM_LDFLAGS)
+LibraryPeriodicCSMPITest_LDADD = $(TESTLIBS)
+LibraryDirichletFASTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/dirichlet_fas.cpp
+LibraryDirichletFASTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryDirichletFASTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryDirichletFASTest_LDADD = $(TESTLIBS)
+LibraryDirichletFASLRMPITest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/dirichlet_fas_lr_mpi.cpp
+LibraryDirichletFASLRMPITest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryDirichletFASLRMPITest_LDFLAGS = $(AM_LDFLAGS)
+LibraryDirichletFASLRMPITest_LDADD = $(TESTLIBS)
+LibraryOpenFASMPITest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/open_fas_mpi.cpp
+LibraryOpenFASMPITest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryOpenFASMPITest_LDFLAGS = $(AM_LDFLAGS)
+LibraryOpenFASMPITest_LDADD = $(TESTLIBS)
+LibraryDirichletFASMPITest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/library/dirichlet_fas_mpi.cpp
+LibraryDirichletFASMPITest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryDirichletFASMPITest_LDFLAGS = $(AM_LDFLAGS)
+LibraryDirichletFASMPITest_LDADD = $(TESTLIBS)
+InterpolateSuite_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/unit_test/interpolate_test.cpp
+InterpolateSuite_CPPFLAGS=$(AM_CPPFLAGS)
+InterpolateSuite_LDFLAGS = $(AM_LDFLAGS)
+InterpolateSuite_LDADD = $(TESTLIBS)
+InterpolateFieldSuite_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/unit_test/interpolate_field_test.cpp
+InterpolateFieldSuite_CPPFLAGS=$(AM_CPPFLAGS)
+InterpolateFieldSuite_LDFLAGS = $(AM_LDFLAGS)
+InterpolateFieldSuite_LDADD = $(TESTLIBS)
+LibraryForceSinusUnitTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/unit_test/force_sinus_unit_test.cpp
+LibraryForceSinusUnitTest_CPPFLAGS=$(AM_CPPFLAGS)
+LibraryForceSinusUnitTest_LDFLAGS = $(AM_LDFLAGS)
+LibraryForceSinusUnitTest_LDADD = $(TESTLIBS)
+SmootherGSLexTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/unit_test/smoother_lex_test.cpp
+SmootherGSLexTest_CPPFLAGS=$(AM_CPPFLAGS)
+SmootherGSLexTest_LDFLAGS = $(AM_LDFLAGS)
+SmootherGSLexTest_LDADD = \
+	libvmg_testhelper.la \
+	$(TESTLIBS)
+SmootherGSRBTest_SOURCES = \
+	unit_test/boost_test.cpp \
+	unit_test/unit_test/smoother_rb_test.cpp
+SmootherGSRBTest_CPPFLAGS=$(AM_CPPFLAGS)
+SmootherGSRBTest_LDFLAGS = $(AM_LDFLAGS)
+SmootherGSRBTest_LDADD = \
+	libvmg_testhelper.la \
+	$(TESTLIBS)
+
+endif
+
Index: ThirdParty/vmg/test/unit_test/boost_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/boost_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/boost_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,77 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#endif
+
+#ifndef BOOST_TEST_NO_MAIN
+#define BOOST_TEST_NO_MAIN
+#endif
+
+#ifdef HAVE_BOOST_TEST_UNIT_TEST_HPP
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#else
+#include <boost/test/included/unit_test.hpp>
+#endif
+
+#include <boost/test/framework.hpp>
+#include <boost/test/detail/unit_test_parameters.hpp>
+#include <boost/test/impl/unit_test_main.ipp>
+
+#ifdef BOOST_TEST_ALTERNATIVE_INIT_API
+bool init_unit_test()
+#else
+::boost::unit_test::test_suite* init_unit_test_suite(int argc, char* argv[])
+#endif
+{
+  boost::unit_test::framework::master_test_suite().p_name.value = "VMG Test Suite";
+#ifdef BOOST_TEST_ALTERNATIVE_INIT_API
+  return true;
+#else
+  return NULL;
+#endif
+}
+
+int BOOST_TEST_CALL_DECL
+main(int argc, char* argv[])
+{
+#ifdef HAVE_MPI
+  MPI_Init(&argc, &argv);
+#endif
+
+  // prototype for user's unit test init function
+#ifdef BOOST_TEST_ALTERNATIVE_INIT_API
+  boost::unit_test::init_unit_test_func init_func = &init_unit_test;
+#else
+  boost::unit_test::init_unit_test_func init_func = &init_unit_test_suite;
+#endif
+
+  int rval = ::boost::unit_test::unit_test_main(init_func, argc, argv);
+
+#ifdef HAVE_MPI
+  MPI_Finalize();
+#endif
+
+  return rval;
+}
Index: ThirdParty/vmg/test/unit_test/library/dirichlet_cs.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/dirichlet_cs.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/dirichlet_cs.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,94 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * dirichlet_cs.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "cycles/cycle_cs_dirichlet.hpp"
+#include "comm/comm_serial.hpp"
+#include "level/level_operator_cs.hpp"
+#include "level/level_operator.hpp"
+#include "level/stencils.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#ifdef HAVE_LAPACK
+#include "solver/dgesv.hpp"
+#endif
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryDirichletCSFixture
+{
+  LibraryDirichletCSFixture()
+  {
+    const Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+    new CommSerial(boundary);
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleCSDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryDirichletCSFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryDirichletCSTest, LibraryDirichletCSFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
Index: ThirdParty/vmg/test/unit_test/library/dirichlet_cs_mpi.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/dirichlet_cs_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/dirichlet_cs_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,100 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * dirichlet_cs_mpi.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_cs_dirichlet.hpp"
+#include "level/level_operator_cs.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryDirichletCSMPIFixture
+{
+  LibraryDirichletCSMPIFixture()
+  {
+    const Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleCSDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryDirichletCSMPIFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryDirichletCSMPITest, LibraryDirichletCSMPIFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/library/dirichlet_fas.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/dirichlet_fas.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/dirichlet_fas.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,90 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * dirichlet_fas.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_serial.hpp"
+#include "cycles/cycle_fas_dirichlet.hpp"
+#include "level/level_operator_fas.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryDirichletFASFixture
+{
+  LibraryDirichletFASFixture()
+  {
+    const Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+    new CommSerial(boundary);
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleFASDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryDirichletFASFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryDirichletFASTest, LibraryDirichletFASFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
Index: ThirdParty/vmg/test/unit_test/library/dirichlet_fas_lr.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/dirichlet_fas_lr.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/dirichlet_fas_lr.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,90 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * dirichlet_fas_lr.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_serial.hpp"
+#include "cycles/cycle_fas_dirichlet.hpp"
+#include "level/level_operator_fas.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fv.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryDirichletFASLRFixture
+{
+  LibraryDirichletFASLRFixture()
+  {
+    Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+    new CommSerial(boundary);
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0, 2, 1.6);
+    new DiscretizationPoissonFV(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleFASDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 8);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryDirichletFASLRFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryDirichletFASLRTest, LibraryDirichletFASLRFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
Index: ThirdParty/vmg/test/unit_test/library/dirichlet_fas_lr_mpi.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/dirichlet_fas_lr_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/dirichlet_fas_lr_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,101 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * dirichlet_fas_lr_mpi.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_fas_dirichlet.hpp"
+#include "level/level_operator_fas.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fv.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryDirichletFASLRMPIFixture
+{
+  LibraryDirichletFASLRMPIFixture()
+  {
+    const Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0, 2, 1.6);
+    new DiscretizationPoissonFV(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleFASDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 8);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryDirichletFASLRMPIFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryDirichletFASLRMPITest, LibraryDirichletFASLRMPIFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/library/dirichlet_fas_mpi.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/dirichlet_fas_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/dirichlet_fas_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,101 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * dirichlet_fas_mpi.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_fas_dirichlet.hpp"
+#include "level/level_operator_fas.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryDirichletFASMPIFixture
+{
+  LibraryDirichletFASMPIFixture()
+  {
+    const Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleFASDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryDirichletFASMPIFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryDirichletFASMPITest, LibraryDirichletFASMPIFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/library/force_sinus.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/force_sinus.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/force_sinus.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,126 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * force_sinus.cpp
+ *
+ *  Created on: 13.06.2012
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include <cmath>
+#include <iostream>
+
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_cs_periodic.hpp"
+#include "level/level_operator_cs.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_singular.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/particle.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryForceSinusFixture
+{
+  LibraryForceSinusFixture()
+  {
+    const Boundary boundary(Periodic, Periodic, Periodic);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverSingular>();
+    new CycleCSPeriodic(2);
+
+    new ObjectStorage<vmg_int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<vmg_int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<vmg_int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryForceSinusFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryForceSinusTest, LibraryForceSinusFixture)
+{
+  Particle::Interpolation ip(5);
+  Particle::Particle p;
+  Vector pos;
+
+  MG::Solve();
+
+  const Grid& sol = (*MG::GetSol())(MG::GetSol()->MaxLevel());
+
+  for (pos.X()=0.25; pos.X()<=0.75; pos.X()+=0.01)
+    for (pos.Y()=0.25; pos.Y()<=0.75; pos.Y()+=0.01)
+      for (pos.Z()=0.25; pos.Z()<=0.75; pos.Z()+=0.01) {
+
+	p.Pos() = pos;
+
+	const Index index = (p.Pos()-sol.Extent().Begin())/sol.Extent().MeshWidth()-sol.Global().LocalBegin()+sol.Local().Begin();
+
+	ip.ComputeCoefficients(sol, index);
+	ip.Evaluate(p);
+
+	Vector ref = Vector(std::cos(sine_factor*p.Pos().X())*std::sin(sine_factor*p.Pos().Y())*std::sin(sine_factor*p.Pos().Z()),
+			    std::sin(sine_factor*p.Pos().X())*std::cos(sine_factor*p.Pos().Y())*std::sin(sine_factor*p.Pos().Z()),
+			    std::sin(sine_factor*p.Pos().X())*std::sin(sine_factor*p.Pos().Y())*std::cos(sine_factor*p.Pos().Z()));
+	ref *= -1.0 * sine_factor;
+
+	//for (int i=0; i<3; ++i)
+	  //BOOST_CHECK_SMALL(p.Field()[i] - ref[i], vmg_float(0.01));
+
+      }
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/library/open_fas_mpi.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/open_fas_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/open_fas_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,99 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * open_fas_mpi.cpp
+ *
+ *  Created on: 25.01.2011
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_fas_dirichlet.hpp"
+#include "level/level_operator_fas.hpp"
+#include "level/level_operator.hpp"
+#include "discretization/discretization_poisson_fv.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "mg.hpp"
+
+#include "interface_gaussian.hpp"
+
+using namespace VMG;
+
+struct LibraryOpenFASMPIFixture
+{
+  LibraryOpenFASMPIFixture()
+  {
+    const Boundary boundary(Open, Open, Open);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceGaussian(1.0, boundary, 2, 3, 0.0, 20.0, 9, 1.6);
+    new DiscretizationPoissonFV(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverRegular>();
+    new CycleFASDirichlet(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 12);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-12);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryOpenFASMPIFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryOpenFASMPITest, LibraryOpenFASMPIFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/library/periodic_cs.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/periodic_cs.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/periodic_cs.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,90 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * periodic_cs.cpp
+ *
+ *  Created on: Oct 20, 2010
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_serial.hpp"
+#include "cycles/cycle_cs_periodic.hpp"
+#include "level/level_operator_cs.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gs.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_singular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryPeriodicCSFixture
+{
+  LibraryPeriodicCSFixture()
+  {
+    const Boundary boundary(Periodic, Periodic, Periodic);
+
+    new CommSerial(boundary);
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverSingular>();
+    new CycleCSPeriodic(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryPeriodicCSFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryPeriodicCSTest, LibraryPeriodicCSFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
Index: ThirdParty/vmg/test/unit_test/library/periodic_cs_mpi.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/periodic_cs_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/periodic_cs_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,101 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * periodic_cs_mpi.cpp
+ *
+ *  Created on: Oct 20, 2010
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_cs_periodic.hpp"
+#include "cycles/cycle_cs_periodic_debug.hpp"
+#include "level/level_operator_cs.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gs.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_singular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryPeriodicCSMPIFixture
+{
+  LibraryPeriodicCSMPIFixture()
+  {
+    const Boundary boundary(Periodic, Periodic, Periodic);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverSingular>();
+    new CycleCSPeriodic(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryPeriodicCSMPIFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryPeriodicCSMPITest, LibraryPeriodicCSMPIFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 1e-10);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/library/periodic_fas.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/periodic_fas.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/periodic_fas.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,88 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * periodic_cs.cpp
+ *
+ *  Created on: Oct 20, 2010
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "comm/comm_serial.hpp"
+#include "cycles/cycle_fas_periodic.hpp"
+#include "level/level_operator_fas.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_singular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryPeriodicFASFixture
+{
+  LibraryPeriodicFASFixture()
+  {
+    const Boundary boundary(Periodic, Periodic, Periodic);
+
+    new CommSerial(boundary);
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverSingular>();
+    new CycleFASPeriodic(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryPeriodicFASFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryPeriodicFASTest, LibraryPeriodicFASFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 2e-7);
+}
Index: ThirdParty/vmg/test/unit_test/library/periodic_fas_mpi.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/library/periodic_fas_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/library/periodic_fas_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,99 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * periodic_fas_mpi.cpp
+ *
+ *  Created on: Oct 20, 2010
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/factory.hpp"
+#include "base/math.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "cycles/cycle_fas_periodic.hpp"
+#include "level/level_operator_fas.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gsrb.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_singular.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryPeriodicFASMPIFixture
+{
+  LibraryPeriodicFASMPIFixture()
+  {
+    const Boundary boundary(Periodic, Periodic, Periodic);
+
+    new CommMPI(boundary, new DomainDecompositionMPI());
+    new VMGInterfaces::InterfaceSinus(sine_factor, boundary, 2, 5, 0.0, 1.0);
+    new DiscretizationPoissonFD(2);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new GaussSeidelRB();
+    new Givens<SolverSingular>();
+    new CycleFASPeriodic(2);
+
+    new ObjectStorage<int>("PRESMOOTHSTEPS", 3);
+    new ObjectStorage<int>("POSTSMOOTHSTEPS", 3);
+    new ObjectStorage<int>("MAX_ITERATION", 7);
+    new ObjectStorage<vmg_float>("PRECISION", 1.0e-10);
+
+    MG::PostInit();
+
+    MG::IsInitialized();
+  }
+
+  ~LibraryPeriodicFASMPIFixture()
+  {
+    MG::Destroy();
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryPeriodicFASMPITest, LibraryPeriodicFASMPIFixture)
+{
+  MG::Solve();
+
+  double res_init = MG::GetFactory().Get("INITIAL_RESIDUAL")->Cast< ObjectStorage<double> >()->Val();
+  double res = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+
+  BOOST_CHECK_SMALL(res/res_init, 2e-7);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/unit_test/force_sinus_unit_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/force_sinus_unit_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/force_sinus_unit_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,123 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * force_sinus.cpp
+ *
+ *  Created on: 19.06.2012
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/math.hpp"
+#include "grid/tempgrid.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/particle.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const vmg_float sine_factor = static_cast<vmg_float>(2.0 * Math::pi);
+
+struct LibraryForceSinusFixture
+{
+  LibraryForceSinusFixture()
+  {
+    grid = new TempGrid(129, 0, 0.0, 1.0);
+  }
+
+  ~LibraryForceSinusFixture()
+  {
+    delete grid;
+  }
+
+  TempGrid* grid;
+};
+
+BOOST_FIXTURE_TEST_CASE(LibraryForceSinusUnitTest, LibraryForceSinusFixture)
+{
+  Particle::Interpolation ip(5);
+  Particle::Particle p;
+  Vector pos;
+
+  MG::Solve();
+
+  Index i;
+  for (i.X()=grid->Local().Begin().X(); i.X()<grid->Local().End().X(); ++i.X())
+      for (i.Y()=grid->Local().Begin().Y(); i.Y()<grid->Local().End().Y(); ++i.Y())
+        for (i.Z()=grid->Local().Begin().Z(); i.Z()<grid->Local().End().Z(); ++i.Z()) {
+          (*grid)(i) = std::sin(sine_factor * (grid->Extent().Begin().X() + (i.X() - grid->Local().Begin().X() + grid->Global().LocalBegin().X()) * grid->Extent().MeshWidth().X())) *
+                       std::sin(sine_factor * (grid->Extent().Begin().Y() + (i.Y() - grid->Local().Begin().Y() + grid->Global().LocalBegin().Y()) * grid->Extent().MeshWidth().Y())) *
+                       std::sin(sine_factor * (grid->Extent().Begin().Z() + (i.Z() - grid->Local().Begin().Z() + grid->Global().LocalBegin().Z()) * grid->Extent().MeshWidth().Z()));
+        }
+
+  for (pos.X()=0.25; pos.X()<=0.75; pos.X()+=0.01)
+    for (pos.Y()=0.25; pos.Y()<=0.75; pos.Y()+=0.01)
+      for (pos.Z()=0.25; pos.Z()<=0.75; pos.Z()+=0.01) {
+
+	p.Pos() = pos;
+
+	const Index index = (p.Pos()-grid->Extent().Begin())/grid->Extent().MeshWidth()-grid->Global().LocalBegin()+grid->Local().Begin();
+
+	ip.ComputeCoefficients(*grid, index);
+	ip.Evaluate(p);
+
+	Vector ref = Vector(std::cos(sine_factor*p.Pos().X())*std::sin(sine_factor*p.Pos().Y())*std::sin(sine_factor*p.Pos().Z()),
+			    std::sin(sine_factor*p.Pos().X())*std::cos(sine_factor*p.Pos().Y())*std::sin(sine_factor*p.Pos().Z()),
+			    std::sin(sine_factor*p.Pos().X())*std::sin(sine_factor*p.Pos().Y())*std::cos(sine_factor*p.Pos().Z()));
+	ref *= -1.0 * sine_factor;
+
+	for (int i=0; i<3; ++i)
+	  BOOST_CHECK_SMALL(p.Field()[i] - ref[i], 0.01);
+
+      }
+
+  pos.X() = 0.25;
+  pos.Y() = 0.25;
+  for (pos.Z() = 0.25; pos.Z()<=0.75; pos.Z() += 0.001) {
+	p.Pos() = pos;
+
+	const Index index = (p.Pos()-grid->Extent().Begin())/grid->Extent().MeshWidth()-grid->Global().LocalBegin()+grid->Local().Begin();
+
+	ip.ComputeCoefficients(*grid, index);
+	ip.Evaluate(p);
+
+	Vector ref = Vector(std::cos(sine_factor*p.Pos().X())*std::sin(sine_factor*p.Pos().Y())*std::sin(sine_factor*p.Pos().Z()),
+			    std::sin(sine_factor*p.Pos().X())*std::cos(sine_factor*p.Pos().Y())*std::sin(sine_factor*p.Pos().Z()),
+			    std::sin(sine_factor*p.Pos().X())*std::sin(sine_factor*p.Pos().Y())*std::cos(sine_factor*p.Pos().Z()));
+	ref *= -1.0 * sine_factor;
+        std::cout << pos.Z() << " " << std::abs((ref.Z() - p.Field().Z()) / ref.Z()) << std::endl;
+  }
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/test/unit_test/unit_test/interpolate_field_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/interpolate_field_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/interpolate_field_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,134 @@
+/*
+ * interpolate_test.cpp
+ *
+ *  Created on: 17.04.2012
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/index.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_serial.hpp"
+#include "grid/grid.hpp"
+#include "grid/multigrid.hpp"
+#include "grid/tempgrid.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/particle.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const Vector sine_factor(1.5 * Math::pi, 3.7 * Math::pi, 2.213 * Math::pi);
+
+static inline vmg_float f(const Vector& pos)
+{
+  return std::sin(sine_factor[0]*pos.X())
+       * std::sin(sine_factor[1]*pos.Y())
+       * std::sin(sine_factor[2]*pos.Z());
+}
+
+static inline Vector grad_f(const Vector& pos)
+{
+  return -1.0 * sine_factor * Vector(std::cos(sine_factor[0]*pos.X()) * std::sin(sine_factor[1]*pos.Y()) * std::sin(sine_factor[2]*pos.Z()),
+                                     std::sin(sine_factor[0]*pos.X()) * std::cos(sine_factor[1]*pos.Y()) * std::sin(sine_factor[2]*pos.Z()),
+                                     std::sin(sine_factor[0]*pos.X()) * std::sin(sine_factor[1]*pos.Y()) * std::cos(sine_factor[2]*pos.Z()));
+}
+
+struct InterpolateFixture
+{
+  InterpolateFixture()
+  {
+    grid = new TempGrid(32, 3, 0.0, 1.0);
+
+    Index i;
+    for (i.X() = 0; i.X()<grid->Local().SizeTotal().X(); ++i.X()) {
+      for (i.Y() = 0; i.Y()<grid->Local().SizeTotal().Y(); ++i.Y()) {
+        for (i.Z() = 0; i.Z()<grid->Local().SizeTotal().Z(); ++i.Z()) {
+          const Vector pos = grid->Extent().Begin() + (i - grid->Local().HaloSize1()) * grid->Extent().MeshWidth();
+          (*grid)(i) = f(pos);
+        }
+      }
+    }
+  }
+
+  ~InterpolateFixture()
+  {
+    delete grid;
+  }
+
+  Grid* grid;
+};
+
+BOOST_FIXTURE_TEST_SUITE(InterpolateSuite, InterpolateFixture)
+
+BOOST_AUTO_TEST_CASE(InterpolateFieldTest)
+{
+  Index i;
+  Index i_old = -1;
+  Particle::Interpolation ip(5);
+  Particle::Particle p;
+
+  vmg_float max_error = std::numeric_limits<vmg_float>::min();
+
+  for(p.Pos().X()=grid->Extent().Begin().X(); p.Pos().X()<grid->Extent().End().X(); p.Pos().X() += 0.01 * (grid->Extent().End().X() - grid->Extent().Begin().X()))
+    for(p.Pos().Y()=grid->Extent().Begin().Y(); p.Pos().Y()<grid->Extent().End().Y(); p.Pos().Y() += 0.01 * (grid->Extent().End().Y() - grid->Extent().Begin().Y()))
+      for(p.Pos().Z()=grid->Extent().Begin().Z(); p.Pos().Z()<grid->Extent().End().Z(); p.Pos().Z() += 0.01 * (grid->Extent().End().Z() - grid->Extent().Begin().Z())) {
+
+        i = (p.Pos()-grid->Extent().Begin())/grid->Extent().MeshWidth()-grid->Global().LocalBegin()+grid->Local().Begin();
+
+        if (i != i_old) {
+          ip.ComputeCoefficients(*grid, i);
+          i_old = i;
+        }
+
+        ip.Evaluate(p);
+
+        for (int j=0; j<3; ++j) {
+          BOOST_CHECK_SMALL(p.Field()[j] - grad_f(p.Pos())[j], vmg_float(2.0e-3));
+          if (fabs(grad_f(p.Pos())[j]) > 1e-10)
+          	max_error = std::max(max_error, fabs((p.Field()[j] - grad_f(p.Pos())[j])/grad_f(p.Pos())[j]));
+        }
+      }
+
+  std::printf("Max rel. error: %e\n", max_error);
+
+  max_error = std::numeric_limits<vmg_float>::min();
+  p.Pos().X() = 0.25;
+  p.Pos().Y() = 0.25;
+  for(p.Pos().Z()=0.25; p.Pos().Z()<0.75; p.Pos().Z() += 0.001) {
+
+    i = (p.Pos()-grid->Extent().Begin())/grid->Extent().MeshWidth()-grid->Global().LocalBegin()+grid->Local().Begin();
+
+    if (i != i_old) {
+      ip.ComputeCoefficients(*grid, i);
+      i_old = i;
+    }
+
+    ip.Evaluate(p);
+
+    for (int j=0; j<3; ++j) {
+	  BOOST_CHECK_SMALL(p.Field()[j] - grad_f(p.Pos())[j], vmg_float(2.0e-3));
+      if (fabs(grad_f(p.Pos())[j]) > 1e-10)
+        max_error = std::max(max_error, fabs((p.Field()[j] - grad_f(p.Pos())[j])/grad_f(p.Pos())[j]));
+    }
+  }
+  std::printf("Max rel. error moving along z: %e\n", max_error);
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
Index: ThirdParty/vmg/test/unit_test/unit_test/interpolate_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/interpolate_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/interpolate_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,107 @@
+/*
+ * interpolate_test.cpp
+ *
+ *  Created on: 17.04.2012
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/index.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_serial.hpp"
+#include "grid/grid.hpp"
+#include "grid/multigrid.hpp"
+#include "grid/tempgrid.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/particle.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+
+using namespace VMG;
+
+const Vector sine_factor(1.5 * Math::pi, 3.7 * Math::pi, 2.213 * Math::pi);
+
+static inline vmg_float f(const Vector& pos)
+{
+  return std::sin(sine_factor[0]*pos.X())
+       * std::sin(sine_factor[1]*pos.Y())
+       * std::sin(sine_factor[2]*pos.Z());
+}
+
+static inline Vector grad_f(const Vector& pos)
+{
+  return -1.0 * sine_factor * Vector(std::cos(sine_factor[0]*pos.X()) * std::sin(sine_factor[1]*pos.Y()) * std::sin(sine_factor[2]*pos.Z()),
+                                     std::sin(sine_factor[0]*pos.X()) * std::cos(sine_factor[1]*pos.Y()) * std::sin(sine_factor[2]*pos.Z()),
+                                     std::sin(sine_factor[0]*pos.X()) * std::sin(sine_factor[1]*pos.Y()) * std::cos(sine_factor[2]*pos.Z()));
+}
+
+struct InterpolateFixture
+{
+  InterpolateFixture()
+  {
+    grid = new TempGrid(32, 3, 0.0, 1.0);
+
+    Index i;
+    for (i.X() = 0; i.X()<grid->Local().SizeTotal().X(); ++i.X()) {
+      for (i.Y() = 0; i.Y()<grid->Local().SizeTotal().Y(); ++i.Y()) {
+        for (i.Z() = 0; i.Z()<grid->Local().SizeTotal().Z(); ++i.Z()) {
+          const Vector pos = grid->Extent().Begin() + (i - grid->Local().HaloSize1()) * grid->Extent().MeshWidth();
+          (*grid)(i) = f(pos);
+        }
+      }
+    }
+  }
+
+  ~InterpolateFixture()
+  {
+    delete grid;
+  }
+
+  Grid* grid;
+};
+
+BOOST_FIXTURE_TEST_SUITE(InterpolateSuite, InterpolateFixture)
+
+BOOST_AUTO_TEST_CASE(InterpolateTest)
+{
+  Index i;
+  Particle::Interpolation ip(5);
+  Particle::Particle p;
+
+  for (i.X() = grid->Local().Begin().X(); i.X()<grid->Local().End().X(); ++i.X())
+    for (i.Y() = grid->Local().Begin().Y(); i.Y()<grid->Local().End().Y(); ++i.Y())
+      for (i.Z() = grid->Local().Begin().Z(); i.Z()<grid->Local().End().Z(); ++i.Z()) {
+
+        ip.ComputeCoefficients(*grid, i);
+
+        const Vector pos_begin = grid->Extent().Begin()
+          + (i - grid->Local().Begin() + grid->Global().LocalBegin()) * grid->Extent().MeshWidth();
+
+        p.Pos() = pos_begin;
+
+        ip.Evaluate(p);
+        BOOST_CHECK_SMALL(p.Pot() - grid->GetVal(i), vmg_float(1.0e-12));
+
+        p.Pos() += 0.5 * grid->Extent().MeshWidth();
+        ip.Evaluate(p);
+        BOOST_CHECK_SMALL(p.Pot() - f(p.Pos()), vmg_float(1.0e-4));
+
+      }
+
+}
+
+BOOST_AUTO_TEST_SUITE_END()
Index: ThirdParty/vmg/test/unit_test/unit_test/polynomial_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/polynomial_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/polynomial_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,70 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "base/polynomial.hpp"
+
+using namespace VMG;
+
+BOOST_AUTO_TEST_SUITE(PolynomialTest)
+
+BOOST_AUTO_TEST_CASE(PolynomialTest1)
+{
+  Polynomial p(2, 0.0, 0.0, 1.0);
+  BOOST_REQUIRE_CLOSE(p(2.0), 4.0, 0.0000001);
+}
+
+BOOST_AUTO_TEST_CASE(PolynomialTest2)
+{
+  Polynomial p(0, 4.0);
+  BOOST_REQUIRE_CLOSE(p(0.0), 4.0, 1.0e-16);
+  BOOST_REQUIRE_CLOSE(p(1.0), 4.0, 1.0e-16);
+  BOOST_REQUIRE_CLOSE(p(1000.0), 4.0, 1.0e-16);
+}
+
+BOOST_AUTO_TEST_CASE(PolynomialTest3)
+{
+  Polynomial p(4, -5.0, 7.0, -3.0, 4.5, -4.0);
+  BOOST_REQUIRE_CLOSE(p(1.5), -6.3125, 1.0e-16);
+  BOOST_REQUIRE_CLOSE(p(-7.5), -14780.9375, 1.0e-16);
+}
+
+BOOST_AUTO_TEST_CASE(PolynomialTestCopyConstructor)
+{
+  Polynomial p1(3, 2.5, 33.2, 7.0, 3.0);
+  Polynomial p2(p1);
+  BOOST_REQUIRE_CLOSE(p1(0.0), p2(0.0), 1.0e-16);
+  BOOST_REQUIRE_CLOSE(p1(32.3), p2(32.3), 1.0e-16);
+}
+
+BOOST_AUTO_TEST_CASE(PolynomialTestAssignment)
+{
+  Polynomial p1(3, 2.1, 5.3, -2.3, -4.3);
+  Polynomial p2;
+  p2 = p1;
+  BOOST_REQUIRE_CLOSE(p1(0.0), p2(0.0), 1.0e-16);
+  BOOST_REQUIRE_CLOSE(p1(32.3), p2(32.3), 1.0e-16);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
Index: ThirdParty/vmg/test/unit_test/unit_test/smoother_fixture.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/smoother_fixture.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/smoother_fixture.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,66 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * smoother_fixture.cpp
+ *
+ *  Created on: 20.09.2010
+ *      Author: Julian Iseringhausen, Frederik Heber
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include "base/interface.hpp"
+#include "base/math.hpp"
+#include "comm/comm_serial.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gs.hpp"
+#include "smoother/gsrb.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+#include "unit_test/smoother_fixture.hpp"
+
+using namespace VMG;
+
+VMG::SmootherFixture::SmootherFixture()
+{
+  const Boundary boundary(Dirichlet, Dirichlet, Dirichlet);
+
+  new CommSerial(boundary);
+  new DiscretizationPoissonFD(2);
+  new VMGInterfaces::InterfaceSinus(2.0*Math::pi, boundary, 4, 4, 0.0, 1.0);
+
+  MG::PostInit();
+
+  MG::GetInterface()->ImportRightHandSide(*MG::GetRhs());
+
+  gs = new GaussSeidel(false);
+  gsrb = new GaussSeidelRB(false);
+}
+
+VMG::SmootherFixture::~SmootherFixture()
+{
+  MG::Destroy();
+
+  delete gs;
+  delete gsrb;
+}
+
Index: ThirdParty/vmg/test/unit_test/unit_test/smoother_fixture.hpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/smoother_fixture.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/smoother_fixture.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,45 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * smoother_fixture.hpp
+ *
+ *  Created on: 16.02.2017
+ *      Author: Frederik Heber
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+using namespace VMG;
+
+namespace VMG {
+
+class Smoother;
+
+struct SmootherFixture
+{
+  SmootherFixture();
+  ~SmootherFixture();
+
+  Smoother* gs;
+  Smoother* gsrb;
+};
+
+} /* VMG */
Index: ThirdParty/vmg/test/unit_test/unit_test/smoother_lex_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/smoother_lex_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/smoother_lex_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,67 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * smoother_lex_test.cpp
+ *
+ *  Created on: 16.02.2017
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/interface.hpp"
+#include "base/math.hpp"
+#include "comm/comm_serial.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gs.hpp"
+#include "smoother/gsrb.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+#include "unit_test/smoother_fixture.hpp"
+
+using namespace VMG;
+
+
+BOOST_FIXTURE_TEST_CASE(SmootherGSLexTest, SmootherFixture)
+{
+  double norm;
+  Multigrid& sol = *MG::GetSol();
+  Multigrid& rhs = *MG::GetRhs();
+  Comm& comm = *MG::GetComm();
+
+  sol.ClearAll();
+
+  for (int i=0; i<20; ++i) {
+
+    gs->Run(sol, rhs, 50);
+    norm = comm.ComputeResidualNorm(sol, rhs);
+
+    if (norm < 1e-10)
+      break;
+  }
+
+  BOOST_CHECK_SMALL(norm, 1.0e-10);
+}
+
Index: ThirdParty/vmg/test/unit_test/unit_test/smoother_rb_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/smoother_rb_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/smoother_rb_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,61 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * smoother_rb_test.cpp
+ *
+ *  Created on: 16.02.2017
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include "base/interface.hpp"
+#include "base/math.hpp"
+#include "comm/comm_serial.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "smoother/gs.hpp"
+#include "smoother/gsrb.hpp"
+#include "mg.hpp"
+
+#include "interface_sinus.hpp"
+#include "unit_test/smoother_fixture.hpp"
+
+using namespace VMG;
+
+BOOST_FIXTURE_TEST_CASE(SmootherGSRBTest, SmootherFixture)
+{
+  double norm;
+
+  MG::GetSol()->ClearAll();
+
+  for (int i=0; i<20; ++i) {
+    gsrb->Run(*MG::GetSol(), *MG::GetRhs(), 50);
+    norm = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());
+    if (norm < 1e-10)
+      break;
+  }
+
+  BOOST_CHECK_SMALL(norm, 1e-10);
+}
+
Index: ThirdParty/vmg/test/unit_test/unit_test/solver_test.cpp
===================================================================
--- ThirdParty/vmg/test/unit_test/unit_test/solver_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/test/unit_test/unit_test/solver_test.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,322 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg 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 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * solver_test.cpp
+ *
+ *  Created on: 21.07.2010
+ *      Author: Julian Iseringhausen
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#ifdef HAVE_LAPACK
+#include "solver/dgesv.hpp"
+#include "solver/dsysv.hpp"
+#endif
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+
+using namespace VMG;
+
+#ifdef HAVE_LAPACK
+class DGESVTest : public DGESV<SolverRegular>
+{
+public:
+  DGESVTest() :
+    DGESV<SolverRegular>(5)
+  {
+    this->Mat(0, 0) = 1.0;
+    this->Mat(0, 1) = -4.0;
+    this->Mat(0, 2) = 0.0;
+    this->Mat(0, 3) = -5.0;
+    this->Mat(0, 4) = 0.0;
+
+    this->Mat(1, 0) = 0.0;
+    this->Mat(1, 1) = 2.0;
+    this->Mat(1, 2) = 2.0;
+    this->Mat(1, 3) = 5.0;
+    this->Mat(1, 4) = 4.0;
+
+    this->Mat(2, 0) = 4.0;
+    this->Mat(2, 1) = -16.0;
+    this->Mat(2, 2) = 0.0;
+    this->Mat(2, 3) = -25.0;
+    this->Mat(2, 4) = 0.0;
+
+    this->Mat(3, 0) = -1.0;
+    this->Mat(3, 1) = -2.0;
+    this->Mat(3, 2) = -4.0;
+    this->Mat(3, 3) = -5.0;
+    this->Mat(3, 4) = -8.0;
+
+    this->Mat(4, 0) = 0.0;
+    this->Mat(4, 1) = 0.0;
+    this->Mat(4, 2) = 2.0;
+    this->Mat(4, 3) = 0.0;
+    this->Mat(4, 4) = 5.0;
+
+    this->Rhs(0) = 12.5;
+    this->Rhs(1) = 8.0;
+    this->Rhs(2) = -10.0;
+    this->Rhs(3) = 3.5;
+    this->Rhs(4) = -16.0;
+
+    this->Compute();
+
+    for (int i=0; i<5; ++i)
+      solution_array[i] = this->Sol(i);
+  }
+
+  double solution_array[5];
+};
+
+class DSYSVTest : public DSYSV<SolverRegular>
+{
+public:
+  DSYSVTest() :
+    DSYSV<SolverRegular>(10)
+  {
+    this->Mat(0, 0) = 4.0;
+    this->Mat(0, 1) = -1.0;
+    this->Mat(0, 2) = -1.0;
+    this->Mat(0, 3) = -1.0;
+    this->Mat(0, 4) = 0.0;
+    this->Mat(0, 5) = 0.0;
+    this->Mat(0, 6) = -1.0;
+    this->Mat(0, 7) = 0.0;
+    this->Mat(0, 8) = 0.0;
+    this->Mat(0, 9) = 1.0;
+
+    this->Mat(1, 0) = -1.0;
+    this->Mat(1, 1) = 4.0;
+    this->Mat(1, 2) = -1.0;
+    this->Mat(1, 3) = 0.0;
+    this->Mat(1, 4) = -1.0;
+    this->Mat(1, 5) = 0.0;
+    this->Mat(1, 6) = 0.0;
+    this->Mat(1, 7) = -1.0;
+    this->Mat(1, 8) = 0.0;
+    this->Mat(1, 9) = 1.0;
+
+    this->Mat(2, 0) = -1.0;
+    this->Mat(2, 1) = -1.0;
+    this->Mat(2, 2) = 4.0;
+    this->Mat(2, 3) = 0.0;
+    this->Mat(2, 4) = 0.0;
+    this->Mat(2, 5) = -1.0;
+    this->Mat(2, 6) = 0.0;
+    this->Mat(2, 7) = 0.0;
+    this->Mat(2, 8) = -1.0;
+    this->Mat(2, 9) = 1.0;
+
+    this->Mat(3, 0) = -1.0;
+    this->Mat(3, 1) = 0.0;
+    this->Mat(3, 2) = 0.0;
+    this->Mat(3, 3) = 4.0;
+    this->Mat(3, 4) = -1.0;
+    this->Mat(3, 5) = -1.0;
+    this->Mat(3, 6) = -1.0;
+    this->Mat(3, 7) = 0.0;
+    this->Mat(3, 8) = 0.0;
+    this->Mat(3, 9) = 1.0;
+
+    this->Mat(4, 0) = 0.0;
+    this->Mat(4, 1) = -1.0;
+    this->Mat(4, 2) = 0.0;
+    this->Mat(4, 3) = -1.0;
+    this->Mat(4, 4) = 4.0;
+    this->Mat(4, 5) = -1.0;
+    this->Mat(4, 6) = 0.0;
+    this->Mat(4, 7) = -1.0;
+    this->Mat(4, 8) = 0.0;
+    this->Mat(4, 9) = 1.0;
+
+    this->Mat(5, 0) = 0.0;
+    this->Mat(5, 1) = 0.0;
+    this->Mat(5, 2) = -1.0;
+    this->Mat(5, 3) = -1.0;
+    this->Mat(5, 4) = -1.0;
+    this->Mat(5, 5) = 4.0;
+    this->Mat(5, 6) = 0.0;
+    this->Mat(5, 7) = 0.0;
+    this->Mat(5, 8) = -1.0;
+    this->Mat(5, 9) = 1.0;
+
+    this->Mat(6, 0) = -1.0;
+    this->Mat(6, 1) = 0.0;
+    this->Mat(6, 2) = 0.0;
+    this->Mat(6, 3) = -1.0;
+    this->Mat(6, 4) = 0.0;
+    this->Mat(6, 5) = 0.0;
+    this->Mat(6, 6) = 4.0;
+    this->Mat(6, 7) = -1.0;
+    this->Mat(6, 8) = -1.0;
+    this->Mat(6, 9) = 1.0;
+
+    this->Mat(7, 0) = 0.0;
+    this->Mat(7, 1) = -1.0;
+    this->Mat(7, 2) = 0.0;
+    this->Mat(7, 3) = 0.0;
+    this->Mat(7, 4) = -1.0;
+    this->Mat(7, 5) = 0.0;
+    this->Mat(7, 6) = -1.0;
+    this->Mat(7, 7) = 4.0;
+    this->Mat(7, 8) = -1.0;
+    this->Mat(7, 9) = 1.0;
+
+    this->Mat(8, 0) = 0.0;
+    this->Mat(8, 1) = 0.0;
+    this->Mat(8, 2) = -1.0;
+    this->Mat(8, 3) = 0.0;
+    this->Mat(8, 4) = 0.0;
+    this->Mat(8, 5) = -1.0;
+    this->Mat(8, 6) = -1.0;
+    this->Mat(8, 7) = -1.0;
+    this->Mat(8, 8) = 4.0;
+    this->Mat(8, 9) = 1.0;
+
+    this->Mat(9, 0) = 1.0;
+    this->Mat(9, 1) = 1.0;
+    this->Mat(9, 2) = 1.0;
+    this->Mat(9, 3) = 1.0;
+    this->Mat(9, 4) = 1.0;
+    this->Mat(9, 5) = 1.0;
+    this->Mat(9, 6) = 1.0;
+    this->Mat(9, 7) = 1.0;
+    this->Mat(9, 8) = 1.0;
+    this->Mat(9, 9) = 0.0;
+
+    this->Rhs(0) = -0.125;
+    this->Rhs(1) = -0.125;
+    this->Rhs(2) = -0.125;
+    this->Rhs(3) = -0.125;
+    this->Rhs(4) = 1.0;
+    this->Rhs(5) = -0.125;
+    this->Rhs(6) = -0.125;
+    this->Rhs(7) = -0.125;
+    this->Rhs(8) = -0.125;
+    this->Rhs(9) = 0.0;
+
+    this->Compute();
+
+    for (int i=0; i<10; ++i)
+      solution_array[i] = this->Sol(i);
+  }
+
+  double solution_array[10];
+};
+#endif /* HAVE_LAPACK */
+
+class GivensTest : public Givens<SolverRegular>
+{
+public:
+  GivensTest() :
+    Givens<SolverRegular>(5)
+  {
+    this->Mat(0, 0) = 1.0;
+    this->Mat(0, 1) = -4.0;
+    this->Mat(0, 2) = 0.0;
+    this->Mat(0, 3) = -5.0;
+    this->Mat(0, 4) = 0.0;
+
+    this->Mat(1, 0) = 0.0;
+    this->Mat(1, 1) = 2.0;
+    this->Mat(1, 2) = 2.0;
+    this->Mat(1, 3) = 5.0;
+    this->Mat(1, 4) = 4.0;
+
+    this->Mat(2, 0) = 4.0;
+    this->Mat(2, 1) = -16.0;
+    this->Mat(2, 2) = 0.0;
+    this->Mat(2, 3) = -25.0;
+    this->Mat(2, 4) = 0.0;
+
+    this->Mat(3, 0) = -1.0;
+    this->Mat(3, 1) = -2.0;
+    this->Mat(3, 2) = -4.0;
+    this->Mat(3, 3) = -5.0;
+    this->Mat(3, 4) = -8.0;
+
+    this->Mat(4, 0) = 0.0;
+    this->Mat(4, 1) = 0.0;
+    this->Mat(4, 2) = 2.0;
+    this->Mat(4, 3) = 0.0;
+    this->Mat(4, 4) = 5.0;
+
+    this->Rhs(0) = 12.5;
+    this->Rhs(1) = 8.0;
+    this->Rhs(2) = -10.0;
+    this->Rhs(3) = 3.5;
+    this->Rhs(4) = -16.0;
+
+    this->Compute();
+
+    for (int i=0; i<5; ++i)
+      solution_array[i] = this->Sol(i);
+  }
+
+  double solution_array[5];
+};
+
+BOOST_AUTO_TEST_SUITE(SolverTestSuite)
+
+BOOST_AUTO_TEST_CASE(SolverDGESVTest)
+{
+#ifdef HAVE_LAPACK
+  DGESVTest* dgesv = new DGESVTest();
+  BOOST_CHECK_CLOSE(  8.5, dgesv->solution_array[0], 1.0e-12);
+  BOOST_CHECK_CLOSE(-16.0, dgesv->solution_array[1], 1.0e-12);
+  BOOST_CHECK_CLOSE(-18.0, dgesv->solution_array[2], 1.0e-12);
+  BOOST_CHECK_CLOSE( 12.0, dgesv->solution_array[3], 1.0e-12);
+  BOOST_CHECK_CLOSE(  4.0, dgesv->solution_array[4], 1.0e-12);
+#endif /* HAVE_LAPACK */
+}
+
+BOOST_AUTO_TEST_CASE(SolverDSYSVTest)
+{
+#ifdef HAVE_LAPACK
+  DSYSVTest* dsysv = new DSYSVTest();
+  BOOST_CHECK_CLOSE(-0.0625, dsysv->solution_array[0], 1.0e-12);
+  BOOST_CHECK_SMALL(         dsysv->solution_array[1], 1.0e-12);
+  BOOST_CHECK_CLOSE(-0.0625, dsysv->solution_array[2], 1.0e-12);
+  BOOST_CHECK_SMALL(         dsysv->solution_array[3], 1.0e-12);
+  BOOST_CHECK_CLOSE(   0.25, dsysv->solution_array[4], 1.0e-12);
+  BOOST_CHECK_SMALL(         dsysv->solution_array[5], 1.0e-12);
+  BOOST_CHECK_CLOSE(-0.0625, dsysv->solution_array[6], 1.0e-12);
+  BOOST_CHECK_SMALL(         dsysv->solution_array[7], 1.0e-12);
+  BOOST_CHECK_CLOSE(-0.0625, dsysv->solution_array[8], 1.0e-12);
+#endif
+}
+
+BOOST_AUTO_TEST_CASE(SolverGivensTest)
+{
+  GivensTest* givens = new GivensTest();
+  BOOST_CHECK_CLOSE(  8.5, givens->solution_array[0], 1.0e-12);
+  BOOST_CHECK_CLOSE(-16.0, givens->solution_array[1], 1.0e-12);
+  BOOST_CHECK_CLOSE(-18.0, givens->solution_array[2], 1.0e-12);
+  BOOST_CHECK_CLOSE( 12.0, givens->solution_array[3], 1.0e-12);
+  BOOST_CHECK_CLOSE(  4.0, givens->solution_array[4], 1.0e-12);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
