From: steigemann <steigemann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 24 Sep 2012 20:43:44 +0000 (+0000)
Subject: //---------------------------------------------------------------------------
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=80cd59fc170a549ce9eb77d93963047a4f6203f6;p=dealii-svn.git

//---------------------------------------------------------------------------
//    $Id: petsc_matrix_free.cc 26043 2012-09-24 19:25:57Z steigemann $
//    Version: $Name$
//
//    Copyright (C) 2012 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//---------------------------------------------------------------------------

#include <deal.II/lac/petsc_matrix_free.h>

#ifdef DEAL_II_USE_PETSC

DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
{
  MatrixFree::MatrixFree ()
            : communicator (PETSC_COMM_SELF)
  {
    const int m=0;
    do_reinit (m, m, m, m);
  }



  MatrixFree::MatrixFree (const MPI_Comm     &communicator,
                          const unsigned int  m,
                          const unsigned int  n,
                          const unsigned int  local_rows,
                          const unsigned int  local_columns)
            : communicator (communicator)
  {
    do_reinit (m, n, local_rows, local_columns);
  }



  MatrixFree::MatrixFree (const MPI_Comm     &communicator,
                          const unsigned int  m,
                          const unsigned int  n,
                          const std::vector<unsigned int> &local_rows_per_process,
                          const std::vector<unsigned int> &local_columns_per_process,
                          const unsigned int  this_process)
            : communicator (communicator)
  {
    Assert (local_rows_per_process.size() == local_columns_per_process.size(),
            ExcDimensionMismatch (local_rows_per_process.size(),
                                  local_columns_per_process.size()));
    Assert (this_process < local_rows_per_process.size(),
            ExcInternalError());

    do_reinit (m, n,
               local_rows_per_process[this_process],
               local_columns_per_process[this_process]);
  }



  MatrixFree::MatrixFree (const unsigned int  m,
                          const unsigned int  n,
                          const unsigned int  local_rows,
                          const unsigned int  local_columns)
            : communicator (MPI_COMM_WORLD)
  {
    do_reinit (m, n, local_rows, local_columns);
  }



  MatrixFree::MatrixFree (const unsigned int  m,
                          const unsigned int  n,
                          const std::vector<unsigned int> &local_rows_per_process,
                          const std::vector<unsigned int> &local_columns_per_process,
                          const unsigned int  this_process)
            : communicator (MPI_COMM_WORLD)
  {
    Assert (local_rows_per_process.size() == local_columns_per_process.size(),
            ExcDimensionMismatch (local_rows_per_process.size(),
                                  local_columns_per_process.size()));
    Assert (this_process < local_rows_per_process.size(),
            ExcInternalError());

    do_reinit (m, n,
               local_rows_per_process[this_process],
               local_columns_per_process[this_process]);
  }



  void MatrixFree::reinit (const MPI_Comm     &communicator,
                           const unsigned int  m,
                           const unsigned int  n,
                           const unsigned int  local_rows,
                           const unsigned int  local_columns)
  {
    this->communicator = communicator;

                                     // destroy the matrix and
                                     // generate a new one
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    int ierr = MatDestroy (matrix);
#else
    int ierr = MatDestroy (&matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    do_reinit (m, n, local_rows, local_columns);
  }



  void MatrixFree::reinit (const MPI_Comm     &communicator,
                           const unsigned int  m,
                           const unsigned int  n,
                           const std::vector<unsigned int> &local_rows_per_process,
                           const std::vector<unsigned int> &local_columns_per_process,
                           const unsigned int  this_process)
  {
    Assert (local_rows_per_process.size() == local_columns_per_process.size(),
            ExcDimensionMismatch (local_rows_per_process.size(),
                                  local_columns_per_process.size()));
    Assert (this_process < local_rows_per_process.size(),
            ExcInternalError());

    this->communicator = communicator;

#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    int ierr = MatDestroy (matrix);
#else
    int ierr = MatDestroy (&matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    do_reinit (m, n,
               local_rows_per_process[this_process],
               local_columns_per_process[this_process]);
  }



  void MatrixFree::reinit (const unsigned int  m,
                           const unsigned int  n,
                           const unsigned int  local_rows,
                           const unsigned int  local_columns)
  {
    reinit (MPI_COMM_WORLD, m, n, local_rows, local_columns);
  }



  void MatrixFree::reinit (const unsigned int  m,
                           const unsigned int  n,
                           const std::vector<unsigned int> &local_rows_per_process,
                           const std::vector<unsigned int> &local_columns_per_process,
                           const unsigned int  this_process)
  {
    reinit (MPI_COMM_WORLD, m, n, local_rows_per_process, local_columns_per_process, this_process);
  }



  void MatrixFree::clear ()
  {
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    int ierr = MatDestroy (matrix);
#else
    int ierr = MatDestroy (&matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    const int m=0;
    do_reinit (m, m, m, m);
  }



  void MatrixFree::vmult (Vec  &dst, const Vec  &src) const
  {

//TODO: Translate the given PETSc Vec* vector into a deal.II
// vector so we can call the vmult function with the usual
// interface; then convert back. This could be much more
// efficient, if the PETScWrappers::*::Vector classes
// had a way to simply generate such a vector object from
// a given PETSc Vec* object without allocating new memory
// and without taking ownership of the Vec*

    VectorBase  *x = 0;
    VectorBase  *y = 0;
                                     // because we do not know,
                                     // if dst and src are sequential
                                     // or distributed vectors,
                                     // we ask for the vector-type
                                     // and reinit x and y with
                                     // dealii::PETScWrappers::*::Vector:
    const char  *vec_type;
    int ierr = VecGetType (src, &vec_type);

    PetscInt  local_size;
    ierr = VecGetLocalSize (src, &local_size);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    if (strcmp(vec_type,"mpi") == 0)
      {
        PetscInt  size;
        ierr = VecGetSize (src, &size);
        AssertThrow (ierr == 0, ExcPETScError(ierr));

        x = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size);
        y = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size);
      }
    else if (strcmp(vec_type,"seq") == 0)
      {
        x = new PETScWrappers::Vector (local_size);
        y = new PETScWrappers::Vector (local_size);
      }
    else
      AssertThrow (false, ExcMessage("PETScWrappers::MPI::MatrixFree::do_matrix_vector_action: "
                                     "This only works for Petsc Vec Type = VECMPI | VECSEQ"));

                                     // copy src to x
    x->equ(1., PETScWrappers::VectorBase(src));
                                     // and call vmult(x,y) which must
                                     // be reimplemented in derived classes
    vmult (*y, *x);

    y->compress();
                                     // copy the result back to dst
    ierr = VecCopy (&(*(*y)), dst);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    delete (x);
    delete (y);
  }



  int MatrixFree::matrix_free_mult (Mat  A, Vec  src, Vec  dst)
  {
                                     // create a pointer to this MatrixFree
                                     // object and link the given matrix A
                                     // to the matrix-vector multiplication
                                     // of this MatrixFree object,
    MatrixFree  *this_object;
    int ierr = MatShellGetContext (A, &this_object);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

                                     // call vmult of this object:
    this_object->vmult (dst, src);

    return (0);
  }



  void MatrixFree::do_reinit (const unsigned int  m,
                              const unsigned int  n,
                              const unsigned int  local_rows,
                              const unsigned int  local_columns)
  {
    Assert (local_rows <= m, ExcDimensionMismatch (local_rows, m));
    Assert (local_columns <= n, ExcDimensionMismatch (local_columns, n));

    int ierr;
                                     // create a PETSc MatShell matrix-type
                                     // object of dimension m x n and local size
                                     // local_rows x local_columns
    ierr = MatCreateShell(communicator, local_rows, local_columns, m, n, (void*)this, &matrix);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
                                     // register the MatrixFree::matrix_free_mult function
                                     // as the matrix multiplication used by this matrix
    ierr = MatShellSetOperation (matrix, MATOP_MULT,
               (void(*)(void))&dealii::PETScWrappers::MatrixFree::matrix_free_mult);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = MatSetFromOptions (matrix);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }
}


DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_USE_PETSC


git-svn-id: https://svn.dealii.org/trunk@26688 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/base/quadrature_lib.h b/deal.II/include/deal.II/base/quadrature_lib.h
index a73a6d2f26..678bb73933 100644
--- a/deal.II/include/deal.II/base/quadrature_lib.h
+++ b/deal.II/include/deal.II/base/quadrature_lib.h
@@ -172,7 +172,7 @@ class QTrapez : public Quadrature<dim>
 
 /**
  * Milne-rule. Closed Newton-Cotes formula, exact for polynomials of degree 5.
- * See Stoer: Einführung in die Numerische Mathematik I, p. 102
+ * See Stoer: Einf�hrung in die Numerische Mathematik I, p. 102
  */
 template <int dim>
 class QMilne : public Quadrature<dim>
@@ -184,7 +184,7 @@ class QMilne : public Quadrature<dim>
 
 /**
  * Weddle-rule. Closed Newton-Cotes formula, exact for polynomials of degree 7.
- * See Stoer: Einführung in die Numerische Mathematik I, p. 102
+ * See Stoer: Einf�hrung in die Numerische Mathematik I, p. 102
  */
 template <int dim>
 class QWeddle : public Quadrature<dim>
@@ -430,6 +430,75 @@ class QGaussOneOverR : public Quadrature<dim>
 
 
 
+/**
+ * Gauss Quadrature Formula with $1/R^{3/2}$ weighting function. This formula
+ * can be used to to integrate $1/R^{3/2} \ f(x)$ on the reference
+ * element $[0,1]^2$, where $f$ is a smooth function without
+ * singularities, and $R$ is the distance from the point $x$ to the vertex
+ * $\xi$, given at construction time by specifying its index. Notice that
+ * this distance is evaluated in the reference element.
+ *
+ * This quadrature formula is a specialization of QGaussOneOverR.
+ * We apply a second transformation $R = t^2$ to cancel a singularity
+ * of order $R^{1/2}$:
+ * \f[
+ * \int_0^1 \int_0^1 \frac{1}{R^{1/2}} f(x,y) \biggr|_{x < y} dxdy =
+ * \int_0^1 \int_0^{r(\pi/4 v)} \frac{1}{R^{1/2}} f(R \cos(\pi/4 v), R \sin(\pi/4 v)) \frac{\pi}{4} R dr dv =
+ * 2 \int_0^1 \int_0^{\sqrt{r(\pi/4 v)}} f(t^2 \cos(\pi/4 v), t^2 \sin(\pi/4 v)) \frac{\pi}{4} t^2 dt dv =
+ * 2 \int_0^1 \int_0^1 f(t^2, t^2 \tan(\pi/4 v)) \frac{\pi}{4} \frac{t^2}{(\cos(\pi/4 v))^{3/2}} dt dv
+ * \f]
+ *
+ * Upon construction it is possible to specify wether we want the
+ * singularity removed, or not. In other words, this quadrature can be
+ * used to integrate $g(x) = 1/R^{3/2}\ f(x)$, or simply $f(x)$, with the $1/R^{3/2}$
+ * factor already included in the quadrature weights.
+ */
+template<int dim>
+class QGaussOneOverRThreeHalfs : public Quadrature<dim>
+{
+  public:
+    /**
+     * The constructor takes three arguments: the order of the Gauss
+     * formula, the index of the vertex where the singularity is
+     * located, and whether we include the weighting singular function
+     * inside the quadrature, or we leave it in the user function to
+     * be integrated. Notice that this constructor only works for the
+     * vertices of the quadrilateral.
+     *
+     * Traditionally, quadrature formulas include their weighting
+     * function, and the last argument is set to false by
+     * default. There are cases, however, where this is undesirable
+     * (for example when you only know that your singularity has the
+     * same order of 1/R, but cannot be written exactly in this
+     * way).
+     *
+     * In other words, you can use this function in either of
+     * the following way, obtaining the same result:
+     *
+     * @code
+     * QGaussOneOverRThreeHalfs singular_quad(order, vertex_id, false);
+     * // This will produce the integral of f(x)/R^{3/2}
+     * for(unsigned int i=0; i<singular_quad.size(); ++i)
+     * 	 integral += f(singular_quad.point(i))*singular_quad.weight(i);
+     *
+     * // And the same here
+     * QGaussOneOverRThreeHalfs singular_quad_noR(order, vertex_id, true);
+     *
+     * // This also will produce the integral of f(x)/R^{3/2}, but 1/R^{3/2} has to
+     * // be specified.
+     * for(unsigned int i=0; i<singular_quad.size(); ++i) {
+     *   double R = (singular_quad_noR.point(i)-cell->vertex(vertex_id)).norm();
+     *   integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R^{3/2};
+     * }
+     * @endcode
+     */
+    QGaussOneOverRThreeHalfs(const unsigned int n,
+                             const unsigned int vertex_index,
+                             const bool factor_out_singular_weight=false);
+};
+
+
+
 /*@}*/
 
 /* -------------- declaration of explicit specializations ------------- */
@@ -460,6 +529,7 @@ template <> QWeddle<1>::QWeddle ();
 template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert);
 template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag);
 template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag);
+template <> QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs (const unsigned int n, const unsigned int index, const bool flag);
 
 
 
diff --git a/deal.II/include/deal.II/base/table_handler.h b/deal.II/include/deal.II/base/table_handler.h
index 0c6366717e..8bcceffebf 100644
--- a/deal.II/include/deal.II/base/table_handler.h
+++ b/deal.II/include/deal.II/base/table_handler.h
@@ -505,6 +505,23 @@ class TableHandler
                                       */
     void write_tex (std::ostream &file, const bool with_header=true) const;
 
+                                     /**
+                                      * Write table as a tex file
+                                      * and use the booktabs package.
+                                      * Values are written in
+                                      * math mode: $value$
+                                      * If with_header is set to false
+                                      * (it is true by default), then
+                                      * no "\documentclass{...}",
+                                      * "\begin{document}" and
+                                      * "\end{document}" are used. In
+                                      * this way the file can be
+                                      * included into an existing tex
+                                      * file using a command like
+                                      * "\input{table_file}".
+                                      */
+    void write_results (std::ostream &file, const bool with_header=true) const;
+
                                      /**
                                       * Read or write the data of this
                                       * object to or from a stream for
@@ -871,6 +888,8 @@ namespace internal
   }
 }
 
+
+
 template <typename T>
 void TableHandler::add_value (const std::string &key,
                               const T            value)
diff --git a/deal.II/include/deal.II/lac/petsc_precondition.h b/deal.II/include/deal.II/lac/petsc_precondition.h
index af3e05acab..43730abcf1 100644
--- a/deal.II/include/deal.II/lac/petsc_precondition.h
+++ b/deal.II/include/deal.II/lac/petsc_precondition.h
@@ -804,8 +804,6 @@ namespace PETScWrappers
           bool output_details;
       };
 
-
-
                                        /**
                                         * Empty Constructor. You need to call
                                         * initialize() before using this
@@ -843,6 +841,263 @@ namespace PETScWrappers
                                         */
       AdditionalData additional_data;
   };
+
+
+
+/**
+ * A class that implements the interface to use the ParaSails sparse
+ * approximate inverse preconditioner from the HYPRE suite. Note that
+ * PETSc has to be configured with HYPRE (e.g. with --download-hypre=1).
+ *
+ * @ingroup PETScWrappers
+ * @author Martin Steigemann, 2011
+ */
+  class PreconditionParaSails : public PreconditionerBase
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional flags to the
+                                        * preconditioner.
+                                        */
+      struct AdditionalData
+      {
+                                           /**
+                                            * Constructor.
+                                            */
+          AdditionalData (
+            const unsigned int symmetric = 0,
+            const unsigned int n_levels = 1,
+            const double threshold = 0.1,
+            const double filter = 0.05,
+            const double load_bal = 0.,
+            const bool output_details = false
+          );
+
+                                           /**
+                                            * This parameter has the following meanings,
+                                            * to indicate the symmetry and definiteness
+                                            * of the problem, and to specify the type
+                                            * of the preconditioner to construct:
+                                            * <ul>
+                                            * <li> @p 0: nonsymmetric and/or indefinite problem, and nonsymmetric preconditioner
+                                            * <li> @p 1: SPD problem, and SPD (factored) preconditioner
+                                            * <li> @p 2: nonsymmetric, definite problem, and SPD (factored) preconditioner
+                                            * </ul>
+                                            *
+                                            */
+          unsigned int symmetric;
+
+          unsigned int n_levels;
+
+          double threshold;
+
+          double filter;
+
+          double load_bal;
+
+                                           /**
+                                            * Setting this flag to true
+                                            * produces output from HYPRE,
+                                            * when the preconditioner
+                                            * is constructed.
+                                            */
+	  bool output_details;
+      };
+
+
+
+                                       /**
+                                        * Empty Constructor. You need to call
+                                        * initialize() before using this
+                                        * object.
+                                        */
+      PreconditionParaSails ();
+
+                                       /**
+                                        * Constructor. Take the matrix which
+                                        * is used to form the preconditioner,
+                                        * and additional flags if there are
+                                        * any.
+                                        */
+      PreconditionParaSails (const MatrixBase     &matrix,
+                             const AdditionalData &additional_data = AdditionalData());
+
+                                       /**
+                                        * Initializes the preconditioner
+                                        * object and calculate all data that
+                                        * is necessary for applying it in a
+                                        * solver. This function is
+                                        * automatically called when calling
+                                        * the constructor with the same
+                                        * arguments and is only used if you
+                                        * create the preconditioner without
+                                        * arguments.
+                                        */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
+    protected:
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular preconditioner.
+                                        */
+      AdditionalData additional_data;
+  };
+
+
+
+/**
+ * A class that implements the interface to use the scalable implementation
+ * of the Parallel ILU algorithm in the Euclid library from the HYPRE suite.
+ * Note that PETSc has to be configured with HYPRE (e.g. with --download-hypre=1).
+ *
+ * @ingroup PETScWrappers
+ * @author Martin Steigemann, 2011
+ */
+  class PreconditionEuclid : public PreconditionerBase
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional flags to the
+                                        * preconditioner.
+                                        */
+      struct AdditionalData
+      {
+                                           /**
+                                            * Constructor.
+                                            */
+          AdditionalData (
+            const unsigned int level = 1,
+            const bool use_block_jacobi = false,
+            const bool output_details = false
+          );
+
+                                           /**
+                                            * Factorization level for ILU(k). Default is 1.
+                                            * For 2D convection-diffusion and similar problems,
+                                            * fastest solution time is typically obtained with
+                                            * levels 4 through 8. For 3D problems, fastest solution
+                                            * time is typically obtained with level 1.
+                                            */
+          unsigned int level;
+
+                                           /**
+                                            * Use Block Jacobi ILU preconditioning
+                                            * instead of PILU. Default is false.
+                                            * If subdomains contain relatively few nodes
+                                            * (less than 1000), or the problem is not
+                                            * well partitioned, Block Jacobi ILU
+                                            * may give faster solution time then PILU.
+                                            */
+          unsigned int use_block_jacobi;
+
+                                           /**
+                                            * Setting this flag to true
+                                            * produces debug output from
+                                            * HYPRE, when the preconditioner
+                                            * is constructed.
+                                            */
+          bool output_details;
+      };
+
+
+
+                                       /**
+                                        * Empty Constructor. You need to call
+                                        * initialize() before using this
+                                        * object.
+                                        */
+      PreconditionEuclid ();
+
+                                       /**
+                                        * Constructor. Take the matrix which
+                                        * is used to form the preconditioner,
+                                        * and additional flags if there are
+                                        * any.
+                                        */
+      PreconditionEuclid (const MatrixBase     &matrix,
+			  const AdditionalData &additional_data = AdditionalData());
+
+                                       /**
+                                        * Initializes the preconditioner
+                                        * object and calculate all data that
+                                        * is necessary for applying it in a
+                                        * solver. This function is
+                                        * automatically called when calling
+                                        * the constructor with the same
+                                        * arguments and is only used if you
+                                        * create the preconditioner without
+                                        * arguments.
+                                        */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
+    protected:
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular preconditioner.
+                                        */
+      AdditionalData additional_data;
+  };
+
+
+
+/**
+ * A class that implements a non-preconditioned Krylov method.
+ *
+ * @ingroup PETScWrappers
+ * @author Martin Steigemann, 2011
+ */
+  class PreconditionNone : public PreconditionerBase
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional flags to the
+                                        * preconditioner.
+                                        */
+      struct AdditionalData
+      {};
+
+                                       /**
+                                        * Empty Constructor. You need to call
+                                        * initialize() before using this
+                                        * object.
+                                        */
+      PreconditionNone ();
+
+                                       /**
+                                        * Constructor. Take the matrix which
+                                        * is used to form the preconditioner,
+                                        * and additional flags if there are
+                                        * any.
+                                        */
+      PreconditionNone (const MatrixBase     &matrix,
+			const AdditionalData &additional_data = AdditionalData());
+
+                                       /**
+                                        * Initializes the preconditioner
+                                        * object and calculate all data that
+                                        * is necessary for applying it in a
+                                        * solver. This function is
+                                        * automatically called when calling
+                                        * the constructor with the same
+                                        * arguments and is only used if you
+                                        * create the preconditioner without
+                                        * arguments.
+                                        */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
+    protected:
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular preconditioner.
+                                        */
+      AdditionalData additional_data;
+  };
 }
 
 
diff --git a/deal.II/include/deal.II/lac/petsc_solver.h b/deal.II/include/deal.II/lac/petsc_solver.h
index 984ea7b6b5..0f3b3b5ab4 100644
--- a/deal.II/include/deal.II/lac/petsc_solver.h
+++ b/deal.II/include/deal.II/lac/petsc_solver.h
@@ -138,6 +138,14 @@ namespace PETScWrappers
              const PreconditionerBase &preconditioner);
 
 
+
+      void
+      solve (const MatrixBase               &A,
+             VectorBase                     &x,
+             const VectorBase               &b,
+             const PreconditionerBase       &preconditioner,
+             const std::vector<VectorBase>  &nullspace);
+
                                        /**
                                         * Resets the contained preconditioner
                                         * and solver object. See class
@@ -169,6 +177,10 @@ namespace PETScWrappers
                       << "An error with error number " << arg1
                       << " occurred while calling a PETSc function");
 
+      DeclException1 (ExcPETScSolverError,
+                      char*,
+                      << "PETSc solver failed: " << arg1);
+
     protected:
 
                                        /**
diff --git a/deal.II/source/base/quadrature_lib.cc b/deal.II/source/base/quadrature_lib.cc
index a38e566433..498f45f66b 100644
--- a/deal.II/source/base/quadrature_lib.cc
+++ b/deal.II/source/base/quadrature_lib.cc
@@ -925,6 +925,93 @@ QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
 }
 
 
+
+template<>
+QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs(const unsigned int n,
+                                                      const unsigned int vertex_index,
+                                                      const bool factor_out_singularity) :
+    Quadrature<2>(2*n*n)
+{
+                                   // This version of the constructor
+                                   // works only for the 4
+                                   // vertices. If you need a more
+                                   // general one, you should use the
+                                   // one with the Point<2> in the
+                                   // constructor.
+    Assert(vertex_index <4, ExcIndexRange(vertex_index, 0, 4));
+
+                                   // Start with the gauss quadrature
+                                   // formula on the (u,v) reference
+                                   // element.
+    QGauss<2> gauss(n);
+
+    Assert(gauss.size() == n*n, ExcInternalError());
+    Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
+
+    // We create only the first one. All other pieces are rotation of
+    // this one.
+    // In this case the transformation is
+    //
+    // (x,y) = (u*u, u*u tan(pi/4 v))
+    //
+    // with Jacobian
+    //
+    // J = 2 pi/4 R*R sqrt(cos(pi/4 v))
+    //
+    // And we get rid of R to take into account the singularity,
+    // unless specified differently in the constructor.
+    std::vector<Point<2> >  &ps = this->quadrature_points;
+    std::vector<double>     &ws = this->weights;
+    double pi4 = numbers::PI/4;
+
+    for(unsigned int q=0; q<gauss.size(); ++q) {
+        const Point<2> &gp = gauss.point(q);
+        ps[q][0] = gp[0]*gp[0];
+        ps[q][1] = gp[0]*gp[0]*std::tan(pi4 *gp[1]);
+        ws[q]    = 2.*gauss.weight(q)*pi4/std::sqrt(std::cos(pi4 *gp[1]));
+        if(factor_out_singularity) {
+            const double  abs_value = (ps[q]-GeometryInfo<2>::unit_cell_vertex(0)).norm();
+            ws[q] *= abs_value;
+            ws[q] *= std::sqrt(abs_value);
+        }
+        // The other half of the quadrilateral is symmetric with
+        // respect to xy plane.
+        ws[gauss.size()+q]    = ws[q];
+        ps[gauss.size()+q][0] = ps[q][1];
+        ps[gauss.size()+q][1] = ps[q][0];
+    }
+
+    // Now we distribute these vertices in the correct manner
+    double theta = 0;
+    switch(vertex_index) {
+    case 0:
+        theta = 0;
+        break;
+    case 1:
+        theta = numbers::PI/2;
+        break;
+    case 2:
+        theta = -numbers::PI/2;
+        break;
+    case 3:
+        theta = numbers::PI;
+        break;
+    }
+
+    double R00 =  std::cos(theta), R01 = -std::sin(theta);
+    double R10 =  std::sin(theta), R11 =  std::cos(theta);
+
+    if(vertex_index != 0)
+        for(unsigned int q=0; q<size(); ++q) {
+            double x = ps[q][0]-.5,  y = ps[q][1]-.5;
+
+            ps[q][0] = R00*x + R01*y + .5;
+            ps[q][1] = R10*x + R11*y + .5;
+        }
+}
+
+
+
 // construct the quadrature formulae in higher dimensions by
 // tensor product of lower dimensions
 
diff --git a/deal.II/source/base/table_handler.cc b/deal.II/source/base/table_handler.cc
index 674651b356..5197347011 100644
--- a/deal.II/source/base/table_handler.cc
+++ b/deal.II/source/base/table_handler.cc
@@ -635,6 +635,136 @@ void TableHandler::write_tex (std::ostream &out, const bool with_header) const
 }
 
 
+void TableHandler::write_results (std::ostream &out, const bool with_header) const
+{
+  bool  math_mode = true;
+
+  AssertThrow (out, ExcIO());
+  if (with_header)
+    out << "\\documentclass[10pt]{report}" << std::endl
+        << "\\usepackage{float}" << std::endl
+        << "\\usepackage{booktabs}" << std::endl
+        << "\\usepackage{amsmath}" << std::endl << std::endl
+        << "\\begin{document}" << std::endl;
+
+  out << "\\begin{table}" << std::endl
+      << "\\begin{center}" << std::endl
+      << "\\begin{tabular}{@{}";
+
+  std::vector<std::string> sel_columns;
+  get_selected_columns(sel_columns);
+
+                                   // write the column formats
+  for (unsigned int j=0; j<column_order.size(); ++j)
+    {
+      std::string key=column_order[j];
+                                       // avoid `supercolumns[key]'
+      const std::map<std::string, std::vector<std::string> >::const_iterator 
+        super_iter=supercolumns.find(key);
+
+      if (super_iter!=supercolumns.end())
+        {
+          const unsigned int n_subcolumns=super_iter->second.size();
+          for (unsigned int k=0; k<n_subcolumns; ++k)
+            {
+                                       // avoid `columns[supercolumns[key]]'
+              const std::map<std::string, Column>::const_iterator
+                col_iter=columns.find(super_iter->second[k]);
+              Assert(col_iter!=columns.end(), ExcInternalError());
+
+              out << col_iter->second.tex_format;
+            }
+        }
+      else
+        {
+                                     // avoid `columns[key]';
+          const std::map<std::string, Column>::const_iterator
+            col_iter=columns.find(key);
+          Assert(col_iter!=columns.end(), ExcInternalError());
+          out << col_iter->second.tex_format;
+        }
+    }
+  out << "@{}} \\toprule" << std::endl;
+
+                                   // write the caption line of the table
+  for (unsigned int j=0; j<column_order.size(); ++j)
+    {
+      std::string key=column_order[j];
+      const std::map<std::string, std::vector<std::string> >::const_iterator 
+        super_iter=supercolumns.find(key);
+
+      if (super_iter!=supercolumns.end())
+        {
+          const unsigned int n_subcolumns=super_iter->second.size();
+                                       // avoid use of `tex_supercaptions[key]'
+          std::map<std::string,std::string>::const_iterator 
+            tex_super_cap_iter=tex_supercaptions.find(key);
+          out << std::endl << "\\multicolumn{" << n_subcolumns << "}{c}{" 
+              << tex_super_cap_iter->second << "}";
+        }
+      else
+        {
+                                       // col_iter->second=columns[col];
+          const std::map<std::string, Column>::const_iterator
+            col_iter=columns.find(key);
+          Assert(col_iter!=columns.end(), ExcInternalError());
+          out << col_iter->second.tex_caption;
+        }
+      if (j<column_order.size()-1)
+        out << " & ";
+    }
+  out << "\\\\ \\midrule" << std::endl;
+
+                                   // write the n rows
+  const unsigned int nrows=n_rows();
+  for (unsigned int i=0; i<nrows; ++i)
+    {
+      const unsigned int n_cols=sel_columns.size();
+
+      for (unsigned int j=0; j<n_cols; ++j)
+        {
+          std::string key=sel_columns[j];
+                                       // avoid `column[key]'
+          const std::map<std::string, Column>::const_iterator
+            col_iter=columns.find(key);
+          Assert(col_iter!=columns.end(), ExcInternalError());
+
+          const Column &column=col_iter->second;
+
+          out << std::setprecision(column.precision);
+
+          if (col_iter->second.scientific)
+            out.setf(std::ios::scientific, std::ios::floatfield);
+          else
+            out.setf(std::ios::fixed, std::ios::floatfield);
+
+          if (math_mode)
+            out << "$";
+
+          out << column.entries[i].value;
+
+          if (math_mode)
+            out << "$";
+
+          if (j<n_cols-1)
+            out << " & ";
+        }
+
+      out << (i == nrows-1 ? "\\\\ \\bottomrule" : "\\\\") << std::endl;
+    }
+
+  out << "\\end{tabular}" << std::endl
+      << "\\end{center}" << std::endl;
+  if(tex_table_caption!="")
+    out << "\\caption{"  << tex_table_caption << "}" << std::endl;
+  if(tex_table_label!="")
+    out << "\\caption{"   << tex_table_label << "}" << std::endl;
+    out << "\\end{table}" << std::endl;
+  if (with_header)
+    out << "\\end{document}" << std::endl;
+}
+
+
 unsigned int TableHandler::n_rows() const
 {
   if (columns.size() == 0)
diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc
index e025665818..8d794272cc 100644
--- a/deal.II/source/distributed/tria.cc
+++ b/deal.II/source/distributed/tria.cc
@@ -3005,9 +3005,79 @@ namespace parallel
     template <int dim, int spacedim>
     void
     Triangulation<dim,spacedim>::
-    copy_triangulation (const dealii::Triangulation<dim, spacedim> &)
+    copy_triangulation (const dealii::Triangulation<dim, spacedim> &old_tria)
     {
-      Assert (false, ExcNotImplemented());
+      clear();
+
+      try
+        {
+          dealii::Triangulation<dim,spacedim>::
+            copy_triangulation (old_tria);
+    }
+      catch (const typename dealii::Triangulation<dim,spacedim>::DistortedCellList &)
+        {
+                                       // the underlying
+                                       // triangulation should not
+                                       // be checking for
+                                       // distorted cells
+          AssertThrow (false, ExcInternalError());
+        }
+
+                                       // note that now we have some content in
+                                       // the p4est objects and call the
+                                       // functions that do the actual work
+                                       // (which are dimension dependent, so
+                                       // separate)
+      triangulation_has_content = true;
+
+      Assert (old_tria.n_levels() == 1,
+              ExcMessage ("Parallel distributed triangulations can only be copied, "
+                          "if they are not refined!"));
+
+      if (dynamic_cast<const dealii::parallel::distributed::Triangulation<dim,spacedim> *>(&old_tria) != 0)
+        {
+          Assert (!(dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+                   (old_tria).refinement_in_progress),
+                  ExcMessage ("Parallel distributed triangulations can only "
+                              "be copied, if no refinement is in progress!"));
+
+          coarse_cell_to_p4est_tree_permutation =
+            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+            (old_tria).coarse_cell_to_p4est_tree_permutation;
+
+          p4est_tree_to_coarse_cell_permutation =
+            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+            (old_tria).p4est_tree_to_coarse_cell_permutation;
+
+          attached_data_size =
+            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+            (old_tria).attached_data_size;
+
+          n_attached_datas   =
+            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+            (old_tria).n_attached_datas;
+        }
+      else
+        {
+          setup_coarse_cell_to_p4est_tree_permutation ();
+        };
+
+      copy_new_triangulation_to_p4est (dealii::internal::int2type<dim>());
+
+      try
+        {
+          copy_local_forest_to_triangulation ();
+        }
+      catch (const typename Triangulation<dim>::DistortedCellList &)
+        {
+                                       // the underlying
+                                       // triangulation should not
+                                       // be checking for
+                                       // distorted cells
+          AssertThrow (false, ExcInternalError());
+        }
+
+      update_number_cache ();
     }
 
 
diff --git a/deal.II/source/lac/petsc_precondition.cc b/deal.II/source/lac/petsc_precondition.cc
index 2df82399d2..370d5c99b5 100644
--- a/deal.II/source/lac/petsc_precondition.cc
+++ b/deal.II/source/lac/petsc_precondition.cc
@@ -483,6 +483,223 @@ namespace PETScWrappers
   }
 
 
+/* ----------------- PreconditionParaSails -------------------- */
+
+  PreconditionParaSails::AdditionalData::
+  AdditionalData(const unsigned int symmetric,
+                 const unsigned int n_levels,
+                 const double threshold,
+                 const double filter,
+                 const double load_bal,
+                 const bool output_details)
+                  :
+                  symmetric(symmetric),
+                  n_levels(n_levels),
+                  threshold(threshold),
+                  filter(filter),
+                  load_bal(load_bal),
+                  output_details(output_details)
+  {}
+
+
+  PreconditionParaSails::PreconditionParaSails ()
+  {}
+
+
+  PreconditionParaSails::PreconditionParaSails (const MatrixBase     &matrix,
+                                                const AdditionalData &additional_data)
+  {
+    initialize(matrix, additional_data);
+  }
+
+
+  void
+  PreconditionParaSails::initialize (const MatrixBase     &matrix_,
+                                     const AdditionalData &additional_data_)
+  {
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+
+#ifdef PETSC_HAVE_HYPRE
+    create_pc();
+
+    int ierr;
+    ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCHYPRESetType(pc, "parasails");
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    if (additional_data.output_details)
+      PetscOptionsSetValue("-pc_hypre_parasails_logging","1");
+
+    Assert ((additional_data.symmetric == 0 ||
+             additional_data.symmetric == 1 ||
+             additional_data.symmetric == 2),
+            ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
+
+    std::stringstream ssStream;
+
+    switch (additional_data.symmetric)
+      {
+        case 0:
+        {
+          ssStream << "nonsymmetric";
+          break;
+        }
+
+        case 1:
+        {
+          ssStream << "SPD";
+          break;
+        }
+
+        case 2:
+        {
+          ssStream << "nonsymmetric,SPD";
+          break;
+        }
+
+        default:
+          Assert (false,
+                  ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
+      };
+
+    PetscOptionsSetValue("-pc_hypre_parasails_sym",ssStream.str().c_str());
+
+    PetscOptionsSetValue("-pc_hypre_parasails_nlevels",
+                         Utilities::int_to_string(
+                           additional_data.n_levels
+                         ).c_str());
+
+    ssStream.str(""); // empty the stringstream
+    ssStream << additional_data.threshold;
+    PetscOptionsSetValue("-pc_hypre_parasails_thresh", ssStream.str().c_str());
+
+    ssStream.str(""); // empty the stringstream
+    ssStream << additional_data.filter;
+    PetscOptionsSetValue("-pc_hypre_parasails_filter", ssStream.str().c_str());
+
+    ssStream.str(""); // empty the stringstream
+    ssStream << additional_data.load_bal;
+    PetscOptionsSetValue("-pc_hypre_parasails_loadbal", ssStream.str().c_str());
+
+    ierr = PCSetFromOptions (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+#else // PETSC_HAVE_HYPRE
+    (void)pc;
+    Assert (false,
+            ExcMessage ("Your PETSc installation does not include a copy of "
+                        "the hypre package necessary for this preconditioner."));
+#endif
+  }
+
+
+/* ----------------- PreconditionEuclid ----------------------- */
+
+  PreconditionEuclid::AdditionalData::
+  AdditionalData(const unsigned int level,
+                 const bool use_block_jacobi,
+                 const bool output_details)
+                  :
+                  level(level),
+                  use_block_jacobi(use_block_jacobi),
+                  output_details(output_details)
+  {}
+
+
+  PreconditionEuclid::PreconditionEuclid ()
+  {}
+
+
+  PreconditionEuclid::PreconditionEuclid (const MatrixBase     &matrix,
+                                                const AdditionalData &additional_data)
+  {
+    initialize(matrix, additional_data);
+  }
+
+
+  void
+  PreconditionEuclid::initialize (const MatrixBase     &matrix_,
+                                  const AdditionalData &additional_data_)
+  {
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+
+#ifdef PETSC_HAVE_HYPRE
+    create_pc();
+
+    int ierr;
+    ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCHYPRESetType(pc, "euclid");
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    PetscOptionsSetValue("-pc_hypre_euclid_levels",
+                         Utilities::int_to_string(
+                           additional_data.level
+                         ).c_str());
+
+    if (additional_data.use_block_jacobi)
+      PetscOptionsSetValue("-pc_hypre_euclid_bj","1");
+
+    if (additional_data.output_details)
+      PetscOptionsSetValue("-pc_hypre_euclid_print_statistics","1");
+
+    ierr = PCSetFromOptions (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+#else // PETSC_HAVE_HYPRE
+    (void)pc;
+    Assert (false,
+            ExcMessage ("Your PETSc installation does not include a copy of "
+                        "the hypre package necessary for this preconditioner."));
+#endif
+  }
+
+
+/* ----------------- PreconditionNone ------------------------- */
+
+  PreconditionNone::PreconditionNone ()
+  {}
+
+
+  PreconditionNone::PreconditionNone (const MatrixBase     &matrix,
+                                      const AdditionalData &additional_data)
+  {
+    initialize(matrix, additional_data);
+  }
+
+
+  void
+  PreconditionNone::initialize (const MatrixBase     &matrix_,
+                                const AdditionalData &additional_data_)
+  {
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+
+    create_pc();
+
+    int ierr;
+    ierr = PCSetType (pc, const_cast<char *>(PCNONE));
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetFromOptions (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+  }
+
+
 /* ----------------- PreconditionLU -------------------- */
 
   PreconditionLU::AdditionalData::
diff --git a/deal.II/source/lac/petsc_solver.cc b/deal.II/source/lac/petsc_solver.cc
index d01dae47f0..0cbb09d2bc 100644
--- a/deal.II/source/lac/petsc_solver.cc
+++ b/deal.II/source/lac/petsc_solver.cc
@@ -129,6 +129,138 @@ namespace PETScWrappers
   }
 
 
+  void
+  SolverBase::solve (const MatrixBase               &A,
+                     VectorBase                     &x,
+                     const VectorBase               &b,
+                     const PreconditionerBase       &preconditioner,
+                     const std::vector<VectorBase>  &nullspace)
+  {
+    int ierr;
+                                     // first create a solver object if this
+                                     // is necessary
+    if (solver_data.get() == 0)
+      {
+        solver_data.reset (new SolverData());
+
+        ierr = KSPCreate (mpi_communicator, &solver_data->ksp);
+        AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+                                         // set the matrices involved. the
+                                         // last argument is irrelevant here,
+                                         // since we use the solver only once
+                                         // anyway
+        ierr = KSPSetOperators (solver_data->ksp, A, preconditioner,
+                                SAME_PRECONDITIONER);
+        AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+                                         // let derived classes set the solver
+                                         // type, and the preconditioning
+                                         // object set the type of
+                                         // preconditioner
+        set_solver_type (solver_data->ksp);
+
+	ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc());
+	AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+        PC  prec;
+        ierr = KSPGetPC (solver_data->ksp, &prec);
+        ierr = PCFactorSetShiftType (prec, MAT_SHIFT_POSITIVE_DEFINITE);
+
+                                         // then a convergence monitor
+                                         // function. that function simply
+                                         // checks with the solver_control
+                                         // object we have in this object for
+                                         // convergence
+#if DEAL_II_PETSC_VERSION_LT(3,0,0)
+        KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
+                               reinterpret_cast<void *>(&solver_control));
+#else
+        KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
+                               reinterpret_cast<void *>(&solver_control),
+                               PETSC_NULL);
+#endif
+
+      }
+
+                                     // then do the real work: set up solver
+                                     // internal data and solve the
+                                     // system. 
+    ierr = KSPSetUp (solver_data->ksp);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    unsigned int dim_nullsp = nullspace.size();
+
+    if (dim_nullsp > 0)
+      {
+        Vec * nullsp_basis;
+
+        PetscMalloc (dim_nullsp*sizeof(Vec), &nullsp_basis);
+
+        for (unsigned int i=0; i<dim_nullsp; ++i)
+          nullsp_basis[i] = &(*nullspace[i]);
+
+        MatNullSpace  nullsp;
+
+        ierr = MatNullSpaceCreate (mpi_communicator, PETSC_FALSE, dim_nullsp, nullsp_basis, &nullsp);
+        AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+        PetscBool  is_nullspace;
+
+        ierr = MatNullSpaceTest (nullsp, A,  &is_nullspace);
+        AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+        ierr = KSPSetNullSpace (solver_data->ksp, nullsp);
+        AssertThrow (ierr == 0, ExcPETScError(ierr));
+      };
+
+    ierr = KSPSolve (solver_data->ksp, b, x);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    KSPConvergedReason  reason;
+    ierr = KSPGetConvergedReason (solver_data->ksp, &reason);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    if (reason < 0)
+      {
+        std::ostringstream error_code;
+
+        if (reason == KSP_DIVERGED_NULL)
+          error_code << "KSP_DIVERGED_NULL";
+        else if (reason == KSP_DIVERGED_ITS)
+          error_code << "KSP_DIVERGED_ITS";
+        else if (reason == KSP_DIVERGED_ITS)
+          error_code << "KSP_DIVERGED_DTOL";
+        else if (reason == KSP_DIVERGED_BREAKDOWN)
+          error_code << "KSP_DIVERGED_BREAKDOWN";
+        else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
+          error_code << "KSP_DIVERGED_BREAKDOWN_BICG";
+        else if (reason == KSP_DIVERGED_NONSYMMETRIC)
+          error_code << "KSP_DIVERGED_NONSYMMETRIC";
+        else if (reason == KSP_DIVERGED_INDEFINITE_PC)
+          error_code << "KSP_DIVERGED_INDEFINITE_PC";
+        else if (reason == KSP_DIVERGED_NAN)
+          error_code << "KSP_DIVERGED_NAN";
+        else if (reason == KSP_DIVERGED_INDEFINITE_MAT)
+          error_code << "KSP_DIVERGED_INDEFINITE_MAT";
+        else
+          error_code << "Unknown Error";
+
+        AssertThrow (false, ExcPETScSolverError(error_code.str().c_str()));
+      };
+
+                                     // do not destroy solver object
+//    solver_data.reset ();
+
+                                     // in case of failure: throw
+                                     // exception
+    if (solver_control.last_check() != SolverControl::success)
+      throw SolverControl::NoConvergence (solver_control.last_step(),
+                                          solver_control.last_value());
+                                     // otherwise exit as normal
+  }
+
+
   void
   SolverBase::set_prefix(const std::string &prefix)
   {
@@ -302,6 +434,9 @@ namespace PETScWrappers
                                      // honor the initial guess in the
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
+
+    KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
+                     PETSC_DEFAULT, this->solver_control.max_steps()+1);
   }
 
 
@@ -404,6 +539,9 @@ namespace PETScWrappers
                                      // honor the initial guess in the
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
+
+    KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
+                     PETSC_DEFAULT, this->solver_control.max_steps()+1);
   }