From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 27 Apr 2001 22:01:17 +0000 (+0000)
Subject: Add filtered matrix, test case, doc updates, etc.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9572da587ebd15d1d32d6c153c5c7a01cbd48b6;p=dealii-svn.git

Add filtered matrix, test case, doc updates, etc.


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

diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h
index a7a669649c..ef34ba6759 100644
--- a/deal.II/deal.II/include/numerics/matrices.h
+++ b/deal.II/deal.II/include/numerics/matrices.h
@@ -502,24 +502,27 @@ class MatrixCreator
  *
  * @sect3{Boundary conditions}
  *
- * The @p{apply_boundary_values} function inserts boundary conditions of
- * into a system of equations.  To actually do this you have to specify
- * a list of degree of freedom indices along with the values these degrees of
- * freedom shall assume. To see how to get such a list, see the discussion
- * of the @ref{VectorTools}@p{::interpolate_boundary_values} function.
- *
- * The inclusion into the assemblage process is as follows: when the matrix and
- * vectors are set up, a list of nodes subject to dirichlet bc is made and
- * matrix and vectors are changed accordingly. This is done by deleting all
- * entries in the matrix in the line of this degree of freedom, setting the
- * main diagonal entry to one and the right hand side element to the
- * boundary value at this node. This forces this node's value to be as specified.
- * To decouple the remaining linear system of equations and to make the system
+ * The @p{apply_boundary_values} function inserts boundary conditions
+ * of into a system of equations.  To actually do this you have to
+ * specify a list of degree of freedom indices along with the values
+ * these degrees of freedom shall assume. To see how to get such a
+ * list, see the discussion of the
+ * @ref{VectorTools}@p{::interpolate_boundary_values} function.
+ *
+ * The inclusion into the assemblage process is as follows: when the
+ * matrix and vectors are set up, a list of nodes subject to dirichlet
+ * bc is made and matrix and vectors are changed accordingly. This is
+ * done by deleting all entries in the matrix in the line of this
+ * degree of freedom, setting the main diagonal entry to one and the
+ * right hand side element to the boundary value at this node. This
+ * forces this node's value to be as specified.  To decouple the
+ * remaining linear system of equations and to make the system
  * symmetric again (at least if it was before), one Gauss elimination
- * step is performed with this line, by adding this (now almost empty) line to
- * all other lines which couple with the given degree of freedom and thus
- * eliminating all coupling between this degree of freedom and others. Now
- * also the column consists only of zeroes, apart from the main diagonal entry.
+ * step is performed with this line, by adding this (now almost empty)
+ * line to all other lines which couple with the given degree of
+ * freedom and thus eliminating all coupling between this degree of
+ * freedom and others. Now also the column consists only of zeroes,
+ * apart from the main diagonal entry.
  *
  * Finding which rows contain an entry in the column for which we are
  * presently performing a Gauss elimination step is either difficult
@@ -558,14 +561,15 @@ class MatrixCreator
  * would be @p{O(N*sqrt(N)*log(m))} for the general case; the latter
  * is too expensive to be performed.
  *
- * It seems as if we had to make clear not to overwrite the lines of other
- * boundary nodes when doing the Gauss elimination step. However, since we
- * reset the right hand side when passing such a node, it is not a problem
- * to change the right hand side values of other boundary nodes not yet
- * processed. It would be a problem to change those entries of nodes already
- * processed, but since the matrix entry of the present column on the row
- * of an already processed node is zero, the Gauss step does not change
- * the right hand side. We need therefore not take special care of other
+ * It seems as if we had to make clear not to overwrite the lines of
+ * other boundary nodes when doing the Gauss elimination
+ * step. However, since we reset the right hand side when passing such
+ * a node, it is not a problem to change the right hand side values of
+ * other boundary nodes not yet processed. It would be a problem to
+ * change those entries of nodes already processed, but since the
+ * matrix entry of the present column on the row of an already
+ * processed node is zero, the Gauss step does not change the right
+ * hand side. We need therefore not take special care of other
  * boundary nodes.
  * 
  * To make solving faster, we preset the solution vector with the
@@ -585,6 +589,17 @@ class MatrixCreator
  * set the entry to the mean of the other diagonal entries, but this
  * seems to be too expensive.
  *
+ * In some cases, it might be interesting to solve several times with
+ * the same matrix, but for different right hand sides or boundary
+ * values. However, since the modification for boundary values of the
+ * right hand side vector depends on the original matrix, this is not
+ * possible without storing the original matrix somewhere and applying
+ * the @p{apply_boundary_conditions} function to a copy of it each
+ * time we want to solve. In that case, you can use the
+ * @ref{FilteredMatrix} class in the @p{LAC} sublibrary. There you can
+ * also find a formal (mathematical) description of the process of
+ * modifying the matrix and right hand side vectors for boundary
+ * values.
  * 
  * @author Wolfgang Bangerth, 1998, 2000
  */
@@ -597,6 +612,11 @@ class MatrixTools : public MatrixCreator<dim>
 				      * to the system matrix and vectors
 				      * as described in the general
 				      * documentation.
+				      *
+				      * For a replacement function,
+				      * see the documentation of the
+				      * @ref{FilteredMatrix} class in
+				      * the @p{LAC} sublibrary.
 				      */
     template <typename number>
     static void
@@ -614,6 +634,11 @@ class MatrixTools : public MatrixCreator<dim>
 				      * documentation. This function
 				      * works for block sparse
 				      * matrices and block vectors
+				      *
+				      * For a replacement function,
+				      * see the documentation of the
+				      * @ref{FilteredMatrix} class in
+				      * the @p{LAC} sublibrary.
 				      */
     static void
     apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
diff --git a/deal.II/doc/future.html b/deal.II/doc/future.html
index cf74def2ff..807c0a9151 100644
--- a/deal.II/doc/future.html
+++ b/deal.II/doc/future.html
@@ -130,8 +130,9 @@
          and <code class="class">SolutionTransfer</code> classes, the
          <code class="class">AssertThrow</code> macro, 1d programs,
 	 curved boundaries and different mappings, using the 
-	 <code class="member">MatrixTools::create_*_matrix</code> functions,
-         and output of more than one variable). This is certainly
+	 <code class="member">MatrixTools::create_*_matrix</code>
+         functions, using the <code class="class">FilteredMatrix</code>
+	 class, and output of more than one variable). This is certainly
          something that should be improved, but is rather time
          consuming.
 	 </p>
diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index 30f2e8e513..6e104a6158 100644
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -249,6 +249,17 @@ documentation, etc</a>.
        (WB 2001/04/25)
        </p>
 
+  <li> <p>
+       New: The <code class="class">FilteredMatrix</code> class is a
+       replacement for the <code
+       class="class">MatrixTools::apply_boundary_values</code>
+       function for cases where you would like to solve several times
+       with the same matrix, either for different right hand sides, or
+       for different boundary values.
+       <br>
+       (WB 2001/04/27)
+       </p>
+
   <li> <p>
        New: There is now a function <code
        class="member">Vector::scale(Vector)</code>
diff --git a/deal.II/lac/include/lac/filtered_matrix.h b/deal.II/lac/include/lac/filtered_matrix.h
new file mode 100644
index 0000000000..03dc7f33a4
--- /dev/null
+++ b/deal.II/lac/include/lac/filtered_matrix.h
@@ -0,0 +1,629 @@
+//----------------------------  filtered_matrix.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 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.
+//
+//----------------------------  filtered_matrix.h  ---------------------------
+#ifndef __deal2__filtered_matrix_h
+#define __deal2__filtered_matrix_h
+
+
+
+#include <base/smartpointer.h>
+#include <base/thread_management.h>
+#include <vector>
+#include <algorithm>
+
+
+template <typename number> class Vector;
+
+
+
+/**
+ * This class is a wrapper for linear systems of equations with simple
+ * equality constraints fixing individual degrees of freedom to a
+ * certain value such as when using Dirichlet boundary
+ * values. Mathematically speaking, it is used to represent a system
+ * of linear equations $Ax=b$ with the constraint that $B_D x = g_D$,
+ * where $B_D$ is a rectangular matrix with exactly one $1$ in each
+ * row, and these $1$s in those columns representing constrained
+ * degrees of freedom (e.g. for Dirichlet boundary nodes, thus the
+ * index $D$) and zeroes for all other diagonal entries, and $g_D$
+ * having the requested nodal values for these constrained
+ * nodes. Thus, the underdetermined equation $B_D x = g_D$ fixes only
+ * the constrained nodes and does not impose any condition on the
+ * others. We note that $B_D B_D^T = 1_D$, where $1_D$ is the identity
+ * matrix with dimension as large as the number of constrained degrees
+ * of freedom. Likewise, $B_D^T B_D$ is the diagonal matrix with
+ * diagonal entries $0$ or $1$ that, when applied to a vector, leaves
+ * all constrained nodes untouched and deletes all unconstrained ones.
+ *
+ * For solving such a system of equations, we first write down the
+ * Lagrangian $L=1/2 x^T A x - x^T b + l^T B_D x$, where $l$
+ * is a Lagrange multiplier for the constraints. The stationarity
+ * condition then reads
+ * \begin{verbatim}
+ * [ A   B_D^T ] [x] = [b  ]
+ * [ B_D 0     ] [l] = [g_D]
+ * \begin{verbatim}
+ *
+ * The first equation then reads $B_D^T l = b-Ax$. On the other hand,
+ * if we left-multiply the first equation by $B_D^T B_D$, we obtain
+ * $B_D^T B_D A x + B_D^T l = B_D^T B_D b$ after equating $B_D B_D^T$
+ * to the identity matrix. Inserting the previous equality, this
+ * yields $(A - B_D^T B_D A) x = (1 - B_D^T B_D)b$. Since
+ * $x=(1 - B_D^T B_D) x + B_D^T B_D x = (1 - B_D^T B_D) x + B_D^T g_D$,
+ * we can restate the linear system:
+ * $A_D x = (1 - B_D^T B_D)b - (1 - B_D^T B_D) A B^T g_D$, where
+ * $A_D = (1 - B_D^T B_D) A (1 - B_D^T B_D)$ is the matrix where all
+ * rows and columns corresponding to constrained nodes have been deleted.
+ *
+ * The last system of equation only defines the value of the
+ * unconstrained nodes, while the constrained ones are determined by
+ * the equation $B_D x = g_D$. We can combine these two linear systems
+ * by using the zeroed out rows of $A_D$: if we set the diagonal to
+ * $1$ and the corresponding zeroed out element of the right hand side
+ * to that of $g_D$, then this fixes the constrained elements as
+ * well. We can write this as follows:
+ * $A_X x = (1 - B_D^T B_D)b - (1 - B_D^T B_D) A B^T g_D + B_D^T g_D$,
+ * where $A_X = A_D + B_D^T B_D$. Note that the two parts of the
+ * latter matrix operate on disjoint subspaces (the first on the
+ * unconstrained nodes, the latter on the constrained ones).
+ *
+ * In iterative solvers, it is not actually necessary to compute $A_X$
+ * explicitely, since only matrix-vector operations need to be
+ * performed. This can be done in a three-step procedure that first
+ * clears all elements in the incoming vector that belong to
+ * constrained nodes, then performs the product with the matrix $A$,
+ * then clears again. This class is a wrapper to this procedure, it
+ * takes a pointer to a matrix with which to perform matrix-vector
+ * products, and does the cleaning of constrained elements itself.
+ * This class therefore implements an overloaded @p{vmult} function
+ * that does the matrix-vector product, as well as @p{Tvmult} for
+ * transpose matrix-vector multiplication and @p{residual} for
+ * residual computation, and can thus be used as a matrix replacement
+ * in lineaer solvers.
+ *
+ * It also has the ability to generate the modification of the right
+ * hand side, through the @ref{apply_constraints} function.
+ *
+ *
+ * @sect3{Connection to other classes}
+ *
+ * The function @p{MatrixTools::apply_boundary_values} does exactly
+ * the same that this class does, except for the fact that that
+ * function actually modifies the matrix. Due to this, it is only
+ * possible to solve with a matrix onto which
+ * @p{MatrixTools::apply_boundary_values} was applied for one right
+ * hand side and one set of boundary values since the modification of
+ * the right hand side depends on the original matrix.
+ *
+ * While this is fine (and the recommended way) in cases where only
+ * one solution of the linear system is required, for example in
+ * solving linear stationary systems, one would often like to have the
+ * ability to solve multiply with the same matrix in nonlinear
+ * problems (where one often does not want to update the Hessian
+ * between Newton steps, despite having different right hand sides in
+ * subsequent steps) or time dependent problems, without having to
+ * re-assemble the matrix or copy it to temporary matrices with which
+ * one then can work. For these cases, this class is meant.
+ *
+ *
+ * @sect3{Usage}
+ *
+ * Usage is simple: create an object of this type, point it to a
+ * matrix that shall be used for $A$ above (either through the
+ * constructor, the copy constructor, or the
+ * @ref{set_referenced_matrix} function), specify the list of boundary
+ * values or other constraints (through the @ref{add_constraints}
+ * function), and then for each required solution modify the right
+ * hand side vector (through @ref{apply_constraints}) and use this
+ * object as matrix object in a linear solver. As linear solvers
+ * should only use @ref{vmult} and @ref{residual} functions of a
+ * matrix class, this class should be as a good a matrix as any other
+ * for that purpose.
+ *
+ * Furthermore, also the @ref{precondition_Jacobi} function is
+ * provided (since the computation of diagonal elements of the
+ * filtered matrix $A_X$ is simple), so you can use this as a
+ * preconditioner. Some other function useful for matrices are also
+ * available.
+ *
+ * A typical code snippet showing the above steps is as follows:
+ * @begin{verbatim}
+ *   ... // set up sparse matrix A and right hand side b somehow
+ *
+ *                     // initialize filtered matrix with
+ *                     // matrix and boundary value constraints
+ *   FilteredMatrix<SparseMatrix<double> > filtered_A (A);
+ *   filtered_A.add_constraints (boundary_values);
+ *
+ *                     // set up a linear solver
+ *   SolverControl control (1000, 1.e-10, false, false);
+ *   PrimitiveVectorMemory<Vector<double> > mem;
+ *   SolverCG<Vector<double> > solver (control, mem);
+ *
+ *                     // set up a preconditioner object
+ *   PreconditionJacobi<FilteredMatrix<SparseMatrix<double> > > prec;
+ *   prec.initialize (filtered_A, 1.2);
+ *
+ *                     // compute modification of right hand side
+ *   filtered_A.apply_constraints (b, true);
+ *
+ *                     // solve for solution vector x
+ *   solver.solve (filtered_A, x, b, prec);
+ * @end{verbatim}
+ *
+ *
+ * @sect3{Template arguments}
+ *
+ * This class takes as template arguments a matrix and a vector
+ * class. The former must provide @p{vmult}, @p{Tvmult}, and
+ * @p{residual} member function that operate on the vector type (the
+ * second template argument). The latter template parameter must
+ * provide access to indivual elements through @p{operator()},
+ * assignment through @p{operator=}.
+ *
+ *
+ * @sect3{Thread-safety}
+ *
+ * The functions that operate as a matrix and do not change the
+ * internal state of this object are synchronised and thus
+ * threadsafe. You need not serialize calls to @p{vmult} or
+ * @p{residual} therefore. Because these functions require the use of
+ * a temporary, they block mutual execution, however. It is necessary
+ * to allocate this temporary vector in class space since otherwise we
+ * would have to allocate such a vector each time one of the member
+ * functions is called (which may be very often for slowly converging
+ * linear systems), which would be a serious performance
+ * bottleneck. If you don't want this serialization of operations, you
+ * have to use several objects of this type.
+ *
+ * @author Wolfgang Bangerth 2001
+ */
+template <class Matrix, class Vector=::Vector<typename Matrix::value_type> >
+class FilteredMatrix : public Subscriptor
+{
+  public:
+				     /**
+				      * Type of matrix entries. In
+				      * analogy to the STL container
+				      * classes.
+				      */
+    typedef typename Matrix::value_type value_type;
+
+				     /**
+				      * Typedef defining a type that
+				      * represents a pair of degree of
+				      * freedom index and the value it
+				      * shall have.
+				      */
+    typedef typename std::pair<unsigned int,value_type> IndexValuePair;
+
+				     /**
+				      * Default constructor. You will
+				      * have to set the matrix to be
+				      * used later using the
+				      * @{set_referenced_matrix}
+				      * function.
+				      */
+    FilteredMatrix ();
+
+				     /**
+				      * Copy constructor. Use the
+				      * matrix and the constraints set
+				      * in the given object for the
+				      * present one as well.
+				      */
+    FilteredMatrix (const FilteredMatrix &fm);
+    
+				     /**
+				      * Constructor. Use the given
+				      * matrix for future operations.
+				      */
+    FilteredMatrix (const Matrix &matrix);
+
+				     /**
+				      * Copy operator. Take over
+				      * matrix and constraints from
+				      * the other object.
+				      */
+    FilteredMatrix & operator = (const FilteredMatrix &fm);
+    
+				     /**
+				      * Set the matrix to be used
+				      * further on. You will probably
+				      * also want to call the
+				      * @ref{clear_constraints}
+				      * function if constraits were
+				      * previously added.
+				      */
+    void set_referenced_matrix (const Matrix &m);
+    
+				     /**
+				      * Return a reference to the
+				      * matrix that is used by this
+				      * object.
+				      */
+    const Matrix & get_referenced_matrix () const;
+
+				     /**
+				      * Add a list of constraints to
+				      * the ones already managed by
+				      * this object. The actual data
+				      * type of this list must be so
+				      * that dereferenced iterators
+				      * are pairs of indices and the
+				      * corresponding values to be
+				      * enforced on the respective
+				      * solution vector's entry. Thus,
+				      * the data type might be, for
+				      * example, a @{std::list} or
+				      * @p{std::vector} of
+				      * @ref{IndexValuePair} objects,
+				      * but also a
+				      * @p{std::map<unsigned,value_type>}.
+				      *
+				      * It is an error if the argument
+				      * contains an entry for a degree
+				      * of freedom that has already
+				      * been constrained
+				      * previously. Furthermore, it is
+				      * assumed that the list of
+				      * constraints is sorted. If the
+				      * input argument is a
+				      * @p{std::map<unsigned,value_type>},
+				      * this is automatically the
+				      * case.
+				      */
+    template <class ConstraintList>
+    void add_constraints (const ConstraintList &new_constraints);
+
+				     /**
+				      * Delete the list of constraints
+				      * presently in use.
+				      */
+    void clear_constraints ();
+
+				     /**
+				      * Apply the constraints to a
+				      * right hand side vector. This
+				      * needs to be done before
+				      * starting to solve with the
+				      * filtered matrix. If the matrix
+				      * is symmetric, set the second
+				      * parameter to @p{true} to use a
+				      * faster algorithm.
+				      */
+    void apply_constraints (Vector     &v,
+			    const bool  matrix_is_symmetric) const;
+
+				     /**
+				      * Return the dimension of the
+				      * image space.  To remember: the
+				      * matrix is of dimension
+				      * $m \times n$.
+				      */
+    unsigned int m () const;
+    
+				     /**
+				      * Return the dimension of the
+				      * range space.  To remember: the
+				      * matrix is of dimension
+				      * $m \times n$.
+				      */
+    unsigned int n () const;
+
+				     /**
+				      * Matrix-vector multiplication:
+				      * let $dst = M*src$ with $M$
+				      * being this matrix. (This
+				      * matrix is the filtered one to
+				      * which we store a reference.)
+				      */
+    void vmult (Vector       &dst,
+		const Vector &src) const;
+    
+				     /**
+				      * Matrix-vector multiplication:
+				      * let $dst = M^T*src$ with $M$
+				      * being this matrix. This
+				      * function does the same as
+				      * @p{vmult} but takes the
+				      * transposed matrix. (This
+				      * matrix is the filtered one to
+				      * which we store a reference.)
+				      *
+				      * Because we need to use a
+				      * temporary variable and since
+				      * we only allocate that each
+				      * time the matrix changed, this
+				      * function only works for square
+				      * matrices.
+				      */
+    void Tvmult (Vector       &dst,
+		 const Vector &src) const;
+  
+				     /**
+				      * Return the square of the norm
+				      * of the vector $v$ with respect
+				      * to the norm induced by this
+				      * matrix,
+				      * i.e. $\left(v,Mv\right)$. This
+				      * is useful, e.g. in the finite
+				      * element context, where the
+				      * $L_2$ norm of a function
+				      * equals the matrix norm with
+				      * respect to the mass matrix of
+				      * the vector representing the
+				      * nodal values of the finite
+				      * element function.
+				      *
+				      * Obviously, the matrix needs to
+				      * be square for this operation.
+				      *
+				      * Note that in many cases, you
+				      * will not want to compute the
+				      * norm with respect to the
+				      * filtered matrix, but with
+				      * respect to the original
+				      * one. For example, if you want
+				      * to compute the $L^2$ norm of a
+				      * vector by forming the matrix
+				      * norm with the mass matrix,
+				      * then you want to use the
+				      * original mass matrix, not the
+				      * filtered one where you might
+				      * have eliminated Dirichlet
+				      * boundary values.
+				      */
+    value_type matrix_norm_square (const Vector &v) const;
+
+				     /**
+				      * Compute the residual of an
+				      * equation @p{Mx=b}, where the
+				      * residual is defined to be
+				      * @p{r=b-Mx} with @p{x}
+				      * typically being an approximate
+				      * of the true solution of the
+				      * equation. Write the residual
+				      * into @p{dst}. The l2 norm of
+				      * the residual vector is
+				      * returned.
+				      *
+				      * Note that it is assumed that
+				      * @{b} is a vector that has been
+				      * treated by the
+				      * @ref{modify_rhs} function,
+				      * since we can then assume that
+				      * the components of the residual
+				      * which correspond to
+				      * constrained degrees of freedom
+				      * do not contribute to the
+				      * residual at all.
+				      */
+    value_type residual (Vector       &dst,
+			 const Vector &x,
+			 const Vector &b) const;
+    
+				     /**
+				      * Apply the Jacobi
+				      * preconditioner, which
+				      * multiplies every element of
+				      * the @p{src} vector by the
+				      * inverse of the respective
+				      * diagonal element and
+				      * multiplies the result with the
+				      * damping factor @p{omega}.
+				      */
+    void precondition_Jacobi (Vector           &dst,
+			      const Vector     &src,
+			      const value_type  omega = 1.) const;
+
+				     /**
+				      * Determine an estimate for the
+				      * memory consumption (in bytes)
+				      * of this object. Since we are
+				      * not the owner of the matrix
+				      * referenced, its memory
+				      * consumption is not included.
+				      */
+    unsigned int memory_consumption () const;
+    
+  private:
+				     /**
+				      * Declare an abbreviation for an
+				      * iterator into the array
+				      * constraint pairs, since that
+				      * data type is so often used and
+				      * is rather awkward to write out
+				      * each time.
+				      */
+    typedef typename std::vector<IndexValuePair>::const_iterator const_index_value_iterator;
+    
+				     /**
+				      * Helper class used to sort
+				      * pairs of indices and
+				      * values. Only the index is
+				      * considered as sort key.
+				      */
+    struct PairComparison 
+    {
+					 /**
+					  * Function comparing the
+					  * pairs @p{i1} and @p{i2}
+					  * for their keys.
+					  */
+	bool operator () (const IndexValuePair &i1,
+			  const IndexValuePair &i2) const;
+    };
+    
+				     /**
+				      * Pointer to the sparsity
+				      * pattern used for this
+				      * matrix. In order to guarantee
+				      * that it is not deleted while
+				      * still in use, we subscribe to
+				      * it using the @p{SmartPointer}
+				      * class.
+				      */
+    SmartPointer<const Matrix> matrix;
+
+				     /**
+				      * Sorted list of pairs denoting
+				      * the index of the variable and
+				      * the value to which it shall be
+				      * fixed.
+				      */
+    std::vector<IndexValuePair> constraints;
+
+				     /**
+				      * Vector to be used as temporary
+				      * storage. Since memory
+				      * allocation is expensive, we do
+				      * not want to allocate temporary
+				      * vectors in each call to
+				      * matrix-vector function, so we
+				      * rather allocate it only once
+				      * and then reuse it over and
+				      * over again. Note that in a
+				      * multithreaded environment, we
+				      * have to synchronise access to
+				      * this vector.
+				      */
+    mutable Vector tmp_vector;
+
+				     /**
+				      * Mutex used to synchronise use
+				      * of the temporary vector.
+				      */
+    mutable Threads::ThreadMutex tmp_mutex;
+    
+				     /**
+				      * Do the pre-filtering step,
+				      * i.e. zero out those components
+				      * that belong to constrained
+				      * degrees of freedom.
+				      */
+    void pre_filter (Vector &v) const;
+
+				     /**
+				      * Do the postfiltering step,
+				      * i.e. set constrained degrees
+				      * of freedom to the value of the
+				      * input vector, as the matrix
+				      * contains only ones on the
+				      * diagonal for these degrees of
+				      * freedom.
+				      */
+    void post_filter (const Vector &in,
+		      Vector       &out) const;
+
+				     /**
+				      * Based on the size of the
+				      * matrix and type of the matrix
+				      * and vector, allocate a
+				      * temporary vector. This
+				      * function has to be overloaded
+				      * for the various template
+				      * parameter choices.
+				      */
+    void allocate_tmp_vector ();
+
+				     /**
+				      * Determine all entries in the
+				      * given column of the matrix
+				      * except for the diagonal entry
+				      * and return their index/value
+				      * pairs. If the matrix is
+				      * symmetric, use a faster
+				      * algorithm.
+				      *
+				      * This function needs to be
+				      * specialised for the different
+				      * matrix types.
+				      */
+    void get_column_entries (const unsigned int           index,
+			     std::vector<IndexValuePair> &column_entries,
+			     const bool                   matrix_is_symmetric) const;
+};
+
+
+/*---------------------- Inline functions -----------------------------------*/
+
+
+template <class Matrix, class Vector>
+inline
+bool
+FilteredMatrix<Matrix,Vector>::PairComparison::
+operator () (const IndexValuePair &i1,
+	     const IndexValuePair &i2) const
+{
+  return (i1.first < i2.first);
+};
+
+
+
+template <class Matrix, class Vector>
+template <class ConstraintList>
+void
+FilteredMatrix<Matrix,Vector>::
+add_constraints (const ConstraintList &new_constraints)
+{
+				   // add new constraints to end
+  const unsigned int old_size = constraints.size();
+  constraints.reserve (old_size + new_constraints.size());
+  constraints.insert (constraints.end(),
+		      new_constraints.begin(),
+		      new_constraints.end());
+				   // then merge the two arrays to
+				   // form one sorted one
+  std::inplace_merge (constraints.begin(),
+		      constraints.begin()+old_size,
+		      constraints.end(),
+		      PairComparison());
+//TODO:[WB] Use equal_range etc to assert that the array is indeed sorted  
+};
+
+
+
+template <class Matrix, class Vector>
+inline
+const Matrix &
+FilteredMatrix<Matrix,Vector>::get_referenced_matrix () const
+{
+  return *matrix;
+};
+
+
+
+template <class Matrix, class Vector>
+inline
+unsigned int FilteredMatrix<Matrix,Vector>::m () const
+{
+  return matrix->m();
+};
+
+
+
+template <class Matrix, class Vector>
+inline
+unsigned int FilteredMatrix<Matrix,Vector>::n () const
+{
+  return matrix->n();
+};
+
+
+
+/*----------------------------   filtered_matrix.h     ---------------------------*/
+
+#endif
+/*----------------------------   filtered_matrix.h     ---------------------------*/
+
+
diff --git a/deal.II/lac/include/lac/filtered_matrix.templates.h b/deal.II/lac/include/lac/filtered_matrix.templates.h
new file mode 100644
index 0000000000..f24a79bda1
--- /dev/null
+++ b/deal.II/lac/include/lac/filtered_matrix.templates.h
@@ -0,0 +1,424 @@
+//----------------------------  filtered_matrix.templates.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 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.
+//
+//----------------------------  filtered_matrix.templates.h  ---------------------------
+#ifndef __deal2__filtered_matrix_templates_h
+#define __deal2__filtered_matrix_templates_h
+
+
+#include <base/memory_consumption.h>
+#include <lac/filtered_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector>::
+FilteredMatrix () 
+{};
+
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector>::
+FilteredMatrix (const FilteredMatrix &fm)
+		:
+		Subscriptor (),
+		constraints (fm.constraints)
+{
+  set_referenced_matrix (*fm.matrix);
+};
+
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector>::
+FilteredMatrix (const Matrix &m)
+{
+  set_referenced_matrix (m);
+};
+
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector> &
+FilteredMatrix<Matrix,Vector>::operator = (const FilteredMatrix &fm)
+{
+  set_referenced_matrix (*fm.matrix);
+  constraints = fm.constraints;
+  return *this;
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::
+set_referenced_matrix (const Matrix &m)
+{
+  matrix = &m;
+  allocate_tmp_vector ();
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::clear_constraints () 
+{
+				   // swap vectors to release memory
+  std::vector<IndexValuePair> empty;
+  constraints.swap (empty);
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::
+apply_constraints (Vector     &v,
+		   const bool  matrix_is_symmetric) const
+{
+				   // array that will hold the pairs
+				   // of index/value of all nonzero
+				   // entries in a given column
+  std::vector<IndexValuePair> column_entries;
+  
+				   // iterate over all constraints and
+				   // treat them one after the other
+  const_index_value_iterator       i = constraints.begin();
+  const const_index_value_iterator e = constraints.end();
+  for (; i!=e; ++i)
+    {
+				       // define abbreviations
+      const unsigned   index = i->first;
+      const value_type value = i->second;
+      
+				       // check whether the value is
+				       // zero, since in that case we do
+				       // not have to modify other nodes
+      if (value != 0)
+	{
+					   // first clear array of
+					   // previous content
+	  column_entries.clear ();
+	  
+					   // then get all entries in
+					   // the present column
+	  get_column_entries (index, column_entries, matrix_is_symmetric);
+	  
+					   // modify rhs for each entry
+	  const_index_value_iterator       col     = column_entries.begin();
+	  const const_index_value_iterator col_end = column_entries.end();
+	  for (; col!=col_end; ++col)
+	    v(col->first) -= col->second * value;
+	};
+    };
+
+  
+				       // finally set constrained
+				       // entries themselves. we can't
+				       // do it in the above loop
+				       // since we might end up
+				       // modifying an entry that we
+				       // have already set if
+				       // constrained dofs couple to
+				       // each other
+  for (i=constraints.begin(); i!=e; ++i)
+    v(i->first) = i->second;
+};
+
+
+
+template <>
+void
+FilteredMatrix<SparseMatrix<double>,Vector<double> >::
+get_column_entries (const unsigned int           index,
+		    std::vector<IndexValuePair> &column_entries,
+		    const bool                   matrix_is_symmetric) const
+{
+				   // depending on whether the matrix
+				   // can be assumed symmetric or not,
+				   // either use a fast or a slow
+				   // algorithm
+  if (matrix_is_symmetric == true)
+				     // ok, matrix is symmetric. we
+				     // may determine the matrix
+				     // entries in this column by
+				     // looking at the matrix entries
+				     // in this row which is
+				     // significantly faster since we
+				     // can traverse them linearly and
+				     // do not have to check each row
+				     // for the possible existence of
+				     // a matrix entry
+    {
+      const unsigned int *
+	col_nums   = &(matrix->get_sparsity_pattern().get_column_numbers()
+		       [matrix->get_sparsity_pattern().get_rowstart_indices()[index]]);
+      const unsigned int
+	row_length = matrix->get_sparsity_pattern().row_length(index);
+
+      for (unsigned int i=0; i<row_length; ++i)
+	{
+	  const unsigned int c = *(col_nums+i);
+
+					   // if not diagonal entry,
+					   // add to list
+	  if (c != index)
+	    column_entries.push_back (std::make_pair(c, (*matrix)(c,index)));
+	};
+    }
+  else
+    {
+				       // otherwise check each row for
+				       // occurrence of an entry in
+				       // this column
+      for (unsigned int row=0; row<n(); ++row)
+	if (row != index)
+	  {
+	    const unsigned int
+	      global_index = matrix->get_sparsity_pattern()(row,index);
+	    if (global_index != SparsityPattern::invalid_entry)
+	      column_entries.push_back (std::make_pair(row,
+						       (*matrix)(row,index)));
+	  };
+    };
+};
+
+
+
+template <>
+void
+FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >::
+get_column_entries (const unsigned int           /*index*/,
+		    std::vector<IndexValuePair> &/*column_entries*/,
+		    const bool                   /*matrix_is_symmetric*/) const
+{
+				   // presently not implemented, but
+				   // should be fairly simple to do
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::pre_filter (Vector &v) const
+{
+				   // iterate over all constraints and
+				   // zero out value
+  const_index_value_iterator       i = constraints.begin();
+  const const_index_value_iterator e = constraints.end();
+  for (; i!=e; ++i)
+    v(i->first) = 0;
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::post_filter (const Vector &in,
+					    Vector       &out) const
+{
+				   // iterate over all constraints and
+				   // set value correctly
+  const_index_value_iterator       i = constraints.begin();
+  const const_index_value_iterator e = constraints.end();
+  for (; i!=e; ++i)
+    out(i->first) = in(i->first);
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::vmult (Vector       &dst,
+				      const Vector &src) const
+{
+  tmp_mutex.acquire ();
+				   // first copy over src vector and
+				   // pre-filter
+  tmp_vector = src;
+  pre_filter (tmp_vector);
+				   // then let matrix do its work
+  matrix->vmult (dst, tmp_vector);
+				   // tmp_vector now no more needed
+  tmp_mutex.release ();
+				   // finally do post-filtering
+  post_filter (src, dst);
+};
+
+
+
+template <class Matrix, class Vector>
+typename FilteredMatrix<Matrix,Vector>::value_type
+FilteredMatrix<Matrix,Vector>::residual (Vector       &dst,
+					 const Vector &x,
+					 const Vector &b) const
+{
+  tmp_mutex.acquire ();
+				   // first copy over x vector and
+				   // pre-filter
+  tmp_vector = x;
+  pre_filter (tmp_vector);
+				   // then let matrix do its work
+  value_type res  = matrix->residual (dst, tmp_vector, b);
+  value_type res2 = res*res;
+				   // tmp_vector now no more needed
+  tmp_mutex.release ();
+				   // finally do post-filtering. here,
+				   // we set constrained indices to
+				   // zero, but have to subtract their
+				   // contributions to the residual
+  const_index_value_iterator       i = constraints.begin();
+  const const_index_value_iterator e = constraints.end();
+  for (; i!=e; ++i)
+    {
+      const value_type v = dst(i->first);
+      res2 -= v*v;
+      dst(i->first) = 0;
+    };
+  
+  Assert (res2>=0, ExcInternalError());
+  return std::sqrt (res2);
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::Tvmult (Vector       &dst,
+				       const Vector &src) const
+{
+  tmp_mutex.acquire ();
+				   // first copy over src vector and
+				   // pre-filter
+  tmp_vector = src;
+  pre_filter (tmp_vector);
+				   // then let matrix do its work
+  matrix->Tvmult (dst, tmp_vector);
+				   // tmp_vector now no more needed
+  tmp_mutex.release ();
+				   // finally do post-filtering
+  post_filter (src, dst);
+};
+
+  
+
+template <class Matrix, class Vector>
+typename FilteredMatrix<Matrix,Vector>::value_type
+FilteredMatrix<Matrix,Vector>::matrix_norm_square (const Vector &v) const
+{
+  tmp_mutex.acquire ();
+  tmp_vector = v;
+
+				   // zero out constrained entries and
+				   // form matrix norm with original
+				   // matrix. this is equivalent to
+				   // forming the matrix norm of the
+				   // original vector with the matrix
+				   // where we have zeroed out rows
+				   // and columns
+  pre_filter (tmp_vector);
+  const value_type ret = matrix->matrix_norm_square (tmp_vector);
+  tmp_mutex.release ();
+  return ret;
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::
+precondition_Jacobi (Vector           &dst,
+		     const Vector     &src,
+		     const value_type  omega) const
+{
+				   // first precondition as usual,
+				   // using the fast algorithms of the
+				   // matrix class
+  matrix->precondition_Jacobi (dst, src, omega);
+
+				   // then modify the constrained
+				   // degree of freedom. as the
+				   // diagonal entries of the filtered
+				   // matrix would be 1.0, simply copy
+				   // over old and new values
+  const_index_value_iterator       i = constraints.begin();
+  const const_index_value_iterator e = constraints.end();
+  for (; i!=e; ++i)
+    dst(i->first) = src(i->first);
+};
+
+
+
+template <class Matrix, class Vector>
+unsigned int
+FilteredMatrix<Matrix,Vector>::memory_consumption () const
+{
+  return (MemoryConsumption::memory_consumption (matrix) +
+	  MemoryConsumption::memory_consumption (constraints) +
+	  MemoryConsumption::memory_consumption (tmp_vector));
+};
+
+
+
+template <>
+void
+FilteredMatrix<SparseMatrix<double>,Vector<double> >::
+allocate_tmp_vector () 
+{
+  tmp_mutex.acquire ();
+  tmp_vector.reinit (matrix->n());
+  tmp_mutex.release ();
+};
+
+
+
+template <>
+void
+FilteredMatrix<SparseMatrix<float>,Vector<float> >::
+allocate_tmp_vector () 
+{
+  tmp_mutex.acquire ();
+  tmp_vector.reinit (matrix->n());
+  tmp_mutex.release ();
+};
+
+
+
+template <>
+void
+FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >::
+allocate_tmp_vector () 
+{
+  tmp_mutex.acquire ();
+  tmp_vector.reinit (matrix->n());
+  tmp_mutex.release ();
+};
+
+
+
+template <>
+void
+FilteredMatrix<BlockSparseMatrix<float>,BlockVector<float> >::
+allocate_tmp_vector () 
+{
+  tmp_mutex.acquire ();
+  tmp_vector.reinit (matrix->n());
+  tmp_mutex.release ();
+};
+
+
+#endif
diff --git a/deal.II/lac/source/filtered_matrix.cc b/deal.II/lac/source/filtered_matrix.cc
new file mode 100644
index 0000000000..40dc8c4242
--- /dev/null
+++ b/deal.II/lac/source/filtered_matrix.cc
@@ -0,0 +1,19 @@
+//----------------------------  filtered_matrix.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 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.
+//
+//----------------------------  filtered_matrix.cc  ---------------------------
+
+
+#include <lac/filtered_matrix.templates.h>
+
+
+template class FilteredMatrix<SparseMatrix<double>,Vector<double> >;
+template class FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >;
diff --git a/tests/deal.II/Makefile.in b/tests/deal.II/Makefile.in
index 2f76a715b2..cbd9661400 100644
--- a/tests/deal.II/Makefile.in
+++ b/tests/deal.II/Makefile.in
@@ -40,13 +40,13 @@ mglocal.exe                  : mglocal.go                            $(lib-2d)
 second_derivatives.exe       : second_derivatives.go                 $(lib-2d)           $(libraries)
 wave-test-3.exe              : wave-test-3.go                        $(lib-2d)           $(libraries)
 support_point_map.exe        : support_point_map.go        $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
-
+filtered_matrix.exe          : filtered_matrix.go          $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 
 
 tests = grid_test dof_test data_out derivatives gradients constraints mg \
         mglocal block_matrices second_derivatives derivative_approximation \
         matrices error_estimator intergrid_constraints intergrid_map \
-        wave-test-3 dof_renumbering support_point_map
+        wave-test-3 dof_renumbering support_point_map filtered_matrix
 
 ############################################################
 
diff --git a/tests/deal.II/filtered_matrix.cc b/tests/deal.II/filtered_matrix.cc
new file mode 100644
index 0000000000..7864d44025
--- /dev/null
+++ b/tests/deal.II/filtered_matrix.cc
@@ -0,0 +1,214 @@
+//----------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// (c) 2001 the deal.II authors
+//
+// purpose of this test:
+//
+// compare results with boundary values eliminated from matrix and
+// vector, and with boundary values treated by filtering
+//
+//
+// Method:
+//
+// solve  (u,v) = (f,v)
+//
+//----------------------------------------------------------------------
+
+
+
+#include <base/quadrature_lib.h>
+#include <base/function_lib.h>
+#include <lac/vector.h>
+#include <lac/vector_memory.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/filtered_matrix.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <vector>
+#include <fstream>
+#include <string>
+
+
+
+  
+void
+solve_filtered (std::map<unsigned int,double> &bv,
+		SparseMatrix<double>          &A,
+		Vector<double>                &u,
+		Vector<double>                &f)
+{
+  FilteredMatrix<SparseMatrix<double> > A1 (A);
+  A1.add_constraints (bv);
+  
+  SolverControl control (1000, 1.e-10, false, false);
+  PrimitiveVectorMemory<Vector<double> > mem;
+  SolverCG<Vector<double> > solver (control, mem);
+  PreconditionJacobi<FilteredMatrix<SparseMatrix<double> > > prec;
+  prec.initialize (A1, 1.2);
+  
+  Vector<double> f1 (f.size());
+  f1 = f;
+  A1.apply_constraints (f1, true);
+  
+  solver.solve (A1, u, f1, prec);
+
+  for (typename std::map<unsigned int,double>::const_iterator i=bv.begin();
+       i!=bv.end(); ++i)
+    Assert (std::fabs(u(i->first) - i->second) < 1e-8,
+	    ExcInternalError());
+};
+
+
+
+template <int dim>
+void
+solve_eliminated (std::map<unsigned int,double> &bv,
+		  SparseMatrix<double>          &A,
+		  Vector<double>                &u,
+		  Vector<double>                &f)
+{
+  MatrixTools<dim>::apply_boundary_values (bv, A, u, f);
+  
+  SolverControl control (1000, 1.e-10, false, false);
+  PrimitiveVectorMemory<Vector<double> > mem;
+  SolverCG<Vector<double> > solver (control, mem);
+  PreconditionJacobi<> prec;
+  prec.initialize (A, 1.2);
+  
+  solver.solve (A, u, f, prec);
+};
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+
+  CosineFunction<dim> cosine;
+  
+  if (dim==2)
+    GridGenerator::hyper_ball(tr);
+  else
+  GridGenerator::hyper_cube(tr, -1,1);
+  
+  tr.refine_global (5-dim);
+  
+  MappingQ<dim> mapping(2);
+  FE_Q<dim> element(1);
+  QGauss4<dim> quadrature;
+  
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+  
+  FEValues<dim> fe (mapping, element, quadrature,
+		    update_values | update_gradients
+		    | update_q_points | update_JxW_values);
+
+  vector <unsigned int> global_dofs (element.dofs_per_cell);
+  vector <double> function (quadrature.n_quadrature_points);
+
+  Vector<double> f (dof.n_dofs ());
+
+  SparsityPattern A_pattern (dof.n_dofs (), dof.n_dofs (),
+			     dof.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern(dof, A_pattern);
+  A_pattern.compress ();
+  
+  SparseMatrix<double> A(A_pattern);
+  
+  DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+  const DoFHandler<dim>::cell_iterator end = dof.end();
+
+  for (; cell != end;++cell)
+    {
+      fe.reinit(cell);
+      cell->get_dof_indices (global_dofs);
+      cosine.value_list (fe.get_quadrature_points(), function);
+
+      for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+	{
+	  double dx = fe.JxW (k);
+	  
+	  for (unsigned int i=0;i<element.dofs_per_cell;++i)
+	    {
+	      const double v = fe.shape_value (i,k);
+	      const Tensor<1,dim> grad_v = fe.shape_grad(i,k);
+	      
+	      double rhs = dx * v * (function[k]);
+
+	      f(global_dofs[i]) += rhs;
+	      for (unsigned int j=0;j<element.dofs_per_cell;++j)
+		{
+		  const Tensor<1,dim> grad_u = fe.shape_grad (j,k);
+		  double el = dx * (grad_u*grad_v);
+		  A.add (global_dofs[i], global_dofs[j], el);
+		}
+	    }
+	}
+    }
+
+				   // interpolate boundary values
+  std::map<unsigned int,double> bv;
+  VectorTools::interpolate_boundary_values (mapping, dof, 0, cosine, bv);
+				   // the cosine has too many zero
+				   // values on the boundary of the
+				   // domain, so reset the elements to
+				   // some other value
+  for (typename std::map<unsigned int,double>::iterator i=bv.begin();
+       i!=bv.end(); ++i)
+    i->second = std::sin(i->second+0.5)+1.0;
+
+				   // first solve filtered. this does
+				   // not change the matrix
+  Vector<double> u_filtered (dof.n_dofs ());
+  solve_filtered (bv, A, u_filtered, f);
+  
+				   // then solve by eliminating in the
+				   // matrix. since this changes the
+				   // matrix, this call must come
+				   // second
+  Vector<double> u_eliminated (dof.n_dofs ());
+  solve_eliminated<dim> (bv, A, u_eliminated, f);
+
+				   // output and check
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      deallog << u_filtered(i) << std::endl;
+      Assert (std::fabs(u_filtered(i) - u_eliminated(i)) < 1e-8,
+	      ExcInternalError());
+    };
+}
+
+
+int main ()
+{
+  ofstream logfile ("filtered_matrix.output");
+  logfile.precision (2);
+  logfile.setf(ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/filtered_matrix.checked b/tests/deal.II/filtered_matrix.checked
new file mode 100644
index 0000000000..68f04927b2
--- /dev/null
+++ b/tests/deal.II/filtered_matrix.checked
@@ -0,0 +1,480 @@
+
+DEAL:1d::1.48
+DEAL:1d::1.64
+DEAL:1d::1.79
+DEAL:1d::1.94
+DEAL:1d::2.08
+DEAL:1d::2.21
+DEAL:1d::2.33
+DEAL:1d::2.43
+DEAL:1d::2.52
+DEAL:1d::2.59
+DEAL:1d::2.65
+DEAL:1d::2.69
+DEAL:1d::2.72
+DEAL:1d::2.74
+DEAL:1d::2.75
+DEAL:1d::2.75
+DEAL:1d::2.75
+DEAL:2d::1.64
+DEAL:2d::1.72
+DEAL:2d::1.74
+DEAL:2d::1.68
+DEAL:2d::1.77
+DEAL:2d::1.78
+DEAL:2d::1.80
+DEAL:2d::1.77
+DEAL:2d::1.72
+DEAL:2d::1.80
+DEAL:2d::1.81
+DEAL:2d::1.81
+DEAL:2d::1.82
+DEAL:2d::1.83
+DEAL:2d::1.82
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.79
+DEAL:2d::1.75
+DEAL:2d::1.81
+DEAL:2d::1.78
+DEAL:2d::1.80
+DEAL:2d::1.81
+DEAL:2d::1.77
+DEAL:2d::1.78
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.72
+DEAL:2d::1.74
+DEAL:2d::1.64
+DEAL:2d::1.68
+DEAL:2d::1.72
+DEAL:2d::1.77
+DEAL:2d::1.79
+DEAL:2d::1.81
+DEAL:2d::1.75
+DEAL:2d::1.78
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.83
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.83
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.83
+DEAL:2d::1.80
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.82
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.84
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.74
+DEAL:2d::1.72
+DEAL:2d::1.77
+DEAL:2d::1.80
+DEAL:2d::1.78
+DEAL:2d::1.77
+DEAL:2d::1.79
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.82
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.87
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.84
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.84
+DEAL:2d::1.82
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.80
+DEAL:2d::1.78
+DEAL:2d::1.81
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.80
+DEAL:2d::1.78
+DEAL:2d::1.77
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.79
+DEAL:2d::1.77
+DEAL:2d::1.75
+DEAL:2d::1.72
+DEAL:2d::1.74
+DEAL:2d::1.72
+DEAL:2d::1.68
+DEAL:2d::1.64
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.90
+DEAL:2d::1.89
+DEAL:2d::1.90
+DEAL:2d::1.89
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.88
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.89
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.90
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.90
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.90
+DEAL:2d::1.89
+DEAL:2d::1.89
+DEAL:2d::1.88
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.88
+DEAL:2d::1.87
+DEAL:2d::1.72
+DEAL:2d::1.74
+DEAL:2d::1.77
+DEAL:2d::1.78
+DEAL:2d::1.80
+DEAL:2d::1.77
+DEAL:2d::1.80
+DEAL:2d::1.81
+DEAL:2d::1.81
+DEAL:2d::1.82
+DEAL:2d::1.83
+DEAL:2d::1.82
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.79
+DEAL:2d::1.81
+DEAL:2d::1.80
+DEAL:2d::1.81
+DEAL:2d::1.77
+DEAL:2d::1.78
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.72
+DEAL:2d::1.74
+DEAL:2d::1.64
+DEAL:2d::1.68
+DEAL:2d::1.72
+DEAL:2d::1.77
+DEAL:2d::1.79
+DEAL:2d::1.81
+DEAL:2d::1.75
+DEAL:2d::1.78
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.83
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.83
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.84
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.83
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.74
+DEAL:2d::1.72
+DEAL:2d::1.77
+DEAL:2d::1.80
+DEAL:2d::1.78
+DEAL:2d::1.77
+DEAL:2d::1.79
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.82
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.81
+DEAL:2d::1.83
+DEAL:2d::1.84
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.85
+DEAL:2d::1.86
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.86
+DEAL:2d::1.85
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.87
+DEAL:2d::1.87
+DEAL:2d::1.86
+DEAL:2d::1.84
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.81
+DEAL:2d::1.80
+DEAL:2d::1.82
+DEAL:2d::1.80
+DEAL:2d::1.78
+DEAL:2d::1.77
+DEAL:2d::1.83
+DEAL:2d::1.81
+DEAL:2d::1.79
+DEAL:2d::1.77
+DEAL:2d::1.74
+DEAL:2d::1.72
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.59
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.59
+DEAL:3d::1.63
+DEAL:3d::1.59
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.59
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.59
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.59
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.55
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.53
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48
+DEAL:3d::1.48