]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add filtered matrix, test case, doc updates, etc.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 27 Apr 2001 22:01:17 +0000 (22:01 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 27 Apr 2001 22:01:17 +0000 (22:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@4498 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/matrices.h
deal.II/doc/future.html
deal.II/doc/news/2001/c-3-1.html
deal.II/lac/include/lac/filtered_matrix.h [new file with mode: 0644]
deal.II/lac/include/lac/filtered_matrix.templates.h [new file with mode: 0644]
deal.II/lac/source/filtered_matrix.cc [new file with mode: 0644]
tests/deal.II/Makefile.in
tests/deal.II/filtered_matrix.cc [new file with mode: 0644]
tests/deal.II/filtered_matrix.checked [new file with mode: 0644]

index a7a669649cb733b52deaec299cdf597b4bad6de7..ef34ba675932a48260ce33799844fa245f53e6b4 100644 (file)
@@ -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,
index cf74def2ff9d3c3a1faf37ae133a162a1aec2c09..807c0a9151c3ec5c82ace47e830f9412018d98a7 100644 (file)
          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>
index 30f2e8e513303acfe2b0e77494048b822090f91d..6e104a6158a351f41941e10e6b156ea643551f4a 100644 (file)
@@ -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 (file)
index 0000000..03dc7f3
--- /dev/null
@@ -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 (file)
index 0000000..f24a79b
--- /dev/null
@@ -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 (file)
index 0000000..40dc8c4
--- /dev/null
@@ -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> >;
index 2f76a715b2b9631e04b669a557d38c773e185c5b..cbd9661400737feb9fcd9fada1dfcb2a984082c1 100644 (file)
@@ -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 (file)
index 0000000..7864d44
--- /dev/null
@@ -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 (file)
index 0000000..68f0492
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.