FilteredMatrix< VECTOR > Class Template Reference
[Derived matrices]

Inheritance diagram for FilteredMatrix< VECTOR >:
Inheritance graph
[legend]

List of all members.

Classes

class  Accessor
class  const_iterator
struct  PairComparison

Public Types

typedef std::pair< unsigned
int, double
IndexValuePair

Public Member Functions

unsigned int memory_consumption () const
Constructors and initialization



 FilteredMatrix ()
 FilteredMatrix (const FilteredMatrix &fm)
template<class MATRIX >
 FilteredMatrix (const MATRIX &matrix, bool expect_constrained_source=false)
FilteredMatrixoperator= (const FilteredMatrix &fm)
template<class MATRIX >
void initialize (const MATRIX &m, bool expect_constrained_source=false)
void clear ()
Managing constraints



void add_constraint (const unsigned int i, const double v)
template<class ConstraintList >
void add_constraints (const ConstraintList &new_constraints)
void clear_constraints ()



void apply_constraints (VECTOR &v, const bool matrix_is_symmetric) const
void vmult (VECTOR &dst, const VECTOR &src) const
void Tvmult (VECTOR &dst, const VECTOR &src) const
void vmult_add (VECTOR &dst, const VECTOR &src) const
void Tvmult_add (VECTOR &dst, const VECTOR &src) const
Iterators



const_iterator begin () const
const_iterator end () const

Private Types

typedef std::vector
< IndexValuePair >
::const_iterator 
const_index_value_iterator

Private Member Functions

void pre_filter (VECTOR &v) const
void post_filter (const VECTOR &in, VECTOR &out) const

Private Attributes

bool expect_constrained_source
std_cxx1x::shared_ptr
< PointerMatrixBase< VECTOR > > 
matrix
std::vector< IndexValuePairconstraints

Friends

class Accessor
class FilteredMatrixBlock< VECTOR >

Detailed Description

template<class VECTOR>
class FilteredMatrix< 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.

In order to accomplish this, the vmult(), Tvmult(), vmult_add() and Tvmult_add functions modify the same function of the original matrix such as if all constrained entries of the source vector were zero. Additionally, all constrained entries of the destination vector are set to zero.

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 set_referenced_matrix() function), specify the list of boundary values or other constraints (through the add_constraints() function), and then for each required solution modify the right hand side vector (through apply_constraints()) and use this object as matrix object in a linear solver. As linear solvers should only use vmult() and residual() functions of a matrix class, this class should be as a good a matrix as any other for that purpose.

Furthermore, also the 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:

 *   ... // set up sparse matrix A and right hand side b somehow
 *
 *                     // initialize filtered matrix with
 *                     // matrix and boundary value constraints
 *   FilteredMatrix<Vector<double> > filtered_A (A);
 *   filtered_A.add_constraints (boundary_values);
 *
 *                     // set up a linear solver
 *   SolverControl control (1000, 1.e-10, false, false);
 *   GrowingVectorMemory<Vector<double> > mem;
 *   SolverCG<Vector<double> > solver (control, mem);
 *
 *                     // set up a preconditioner object
 *   PreconditionJacobi<SparseMatrix<double> > prec;
 *   prec.initialize (A, 1.2);
 *   FilteredMatrix<Vector<double> > filtered_prec (prec);
 *   filtered_prec.add_constraints (boundary_values);
 *
 *                     // compute modification of right hand side
 *   filtered_A.apply_constraints (b, true);
 *
 *                     // solve for solution vector x
 *   solver.solve (filtered_A, x, b, filtered_prec);
 * 

Connection to other classes

The function 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 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 a feasible method 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.

Some background

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

 [ A   B_D^T ] [x] = [b  ]
 [ B_D 0     ] [l] = [g_D]

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$ explicitly, 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 vmult function that does the matrix-vector product, as well as Tvmult for transpose matrix-vector multiplication and residual for residual computation, and can thus be used as a matrix replacement in linear solvers.

It also has the ability to generate the modification of the right hand side, through the apply_constraints() function.

Template arguments

This class takes as template arguments a matrix and a vector class. The former must provide vmult, vmult_add, Tvmult, and residual member function that operate on the vector type (the second template argument). The latter template parameter must provide access to indivual elements through operator(), assignment through operator=.

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 vmult or residual therefore.

Author:
Wolfgang Bangerth 2001, Luca Heltai 2006, Guido Kanschat 2007, 2008

Member Typedef Documentation

template<class VECTOR>
typedef std::pair<unsigned int, double> FilteredMatrix< VECTOR >::IndexValuePair

Typedef defining a type that represents a pair of degree of freedom index and the value it shall have.

template<class VECTOR>
typedef std::vector<IndexValuePair>::const_iterator FilteredMatrix< VECTOR >::const_index_value_iterator [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.


Constructor & Destructor Documentation

template<class VECTOR >
FilteredMatrix< VECTOR >::FilteredMatrix (  )  [inline]

Default constructor. You will have to set the matrix to be used later using initialize().

template<class VECTOR >
FilteredMatrix< VECTOR >::FilteredMatrix ( const FilteredMatrix< VECTOR > &  fm  )  [inline]

Copy constructor. Use the matrix and the constraints set in the given object for the present one as well.

template<class VECTOR >
template<class MATRIX >
FilteredMatrix< VECTOR >::FilteredMatrix ( const MATRIX &  matrix,
bool  expect_constrained_source = false 
) [inline]

Constructor. Use the given matrix for future operations.

  • m: The matrix being used in multiplications.

References FilteredMatrix< VECTOR >::initialize().


Member Function Documentation

template<class VECTOR >
FilteredMatrix< VECTOR > & FilteredMatrix< VECTOR >::operator= ( const FilteredMatrix< VECTOR > &  fm  )  [inline]

Copy operator. Take over matrix and constraints from the other object.

Reimplemented from Subscriptor.

References FilteredMatrix< VECTOR >::constraints, FilteredMatrix< VECTOR >::expect_constrained_source, and FilteredMatrix< VECTOR >::matrix.

template<class VECTOR >
template<class MATRIX >
void FilteredMatrix< VECTOR >::initialize ( const MATRIX &  m,
bool  expect_constrained_source = false 
) [inline]

Set the matrix to be used further on. You will probably also want to call the clear_constraints() function if constraits were previously added.

  • m: The matrix being used in multiplications.

References FilteredMatrix< VECTOR >::expect_constrained_source, and FilteredMatrix< VECTOR >::matrix.

Referenced by FilteredMatrix< VECTOR >::FilteredMatrix().

template<class VECTOR >
void FilteredMatrix< VECTOR >::clear (  )  [inline]

Delete all constraints and the matrix pointer.

References FilteredMatrix< VECTOR >::clear_constraints(), and FilteredMatrix< VECTOR >::matrix.

template<class VECTOR >
void FilteredMatrix< VECTOR >::add_constraint ( const unsigned int  i,
const double  v 
) [inline]

Add the constraint that the value with index i should have the value v.

References FilteredMatrix< VECTOR >::constraints.

template<class VECTOR >
template<class ConstraintList >
void FilteredMatrix< VECTOR >::add_constraints ( const ConstraintList &  new_constraints  )  [inline]

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 std::vector of IndexValuePair objects, but also a std::map<unsigned, double>.

The second component of these pairs will only be used in apply_constraints(). The first is used to set values to zero in matrix vector multiplications.

It is an error if the argument contains an entry for a degree of freedom that has already been constrained previously.

References FilteredMatrix< VECTOR >::constraints.

template<class VECTOR >
void FilteredMatrix< VECTOR >::clear_constraints (  )  [inline]

Delete the list of constraints presently in use.

References FilteredMatrix< VECTOR >::constraints.

Referenced by FilteredMatrix< VECTOR >::clear().

template<class VECTOR >
void FilteredMatrix< VECTOR >::apply_constraints ( VECTOR &  v,
const bool  matrix_is_symmetric 
) const [inline]

Vector operations 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 (i.e. the matrix itself, not only its sparsity pattern), set the second parameter to true to use a faster algorithm.

References GrowingVectorMemory< VECTOR >::alloc(), FilteredMatrix< VECTOR >::constraints, GrowingVectorMemory< VECTOR >::free(), and FilteredMatrix< VECTOR >::matrix.

template<class VECTOR >
void FilteredMatrix< VECTOR >::vmult ( VECTOR &  dst,
const VECTOR &  src 
) const [inline]
template<class VECTOR >
void FilteredMatrix< VECTOR >::Tvmult ( VECTOR &  dst,
const VECTOR &  src 
) const [inline]
template<class VECTOR >
void FilteredMatrix< VECTOR >::vmult_add ( VECTOR &  dst,
const VECTOR &  src 
) const [inline]

Adding matrix-vector multiplication.

Note:
The result vector of this multiplication will have the constraint entries set to zero, independent of the previous value of dst. We excpect that in most cases this is the required behavior.

References GrowingVectorMemory< VECTOR >::alloc(), FilteredMatrix< VECTOR >::expect_constrained_source, GrowingVectorMemory< VECTOR >::free(), FilteredMatrix< VECTOR >::matrix, FilteredMatrix< VECTOR >::post_filter(), and FilteredMatrix< VECTOR >::pre_filter().

template<class VECTOR >
void FilteredMatrix< VECTOR >::Tvmult_add ( VECTOR &  dst,
const VECTOR &  src 
) const [inline]

Adding transpose matrix-vector multiplication:

Note:
The result vector of this multiplication will have the constraint entries set to zero, independent of the previous value of dst. We excpect that in most cases this is the required behavior.

References GrowingVectorMemory< VECTOR >::alloc(), FilteredMatrix< VECTOR >::expect_constrained_source, GrowingVectorMemory< VECTOR >::free(), FilteredMatrix< VECTOR >::matrix, FilteredMatrix< VECTOR >::post_filter(), and FilteredMatrix< VECTOR >::pre_filter().

template<typename number >
FilteredMatrix< number >::const_iterator FilteredMatrix< number >::begin (  )  const [inline]

Iterator to the first constraint.

template<typename number >
FilteredMatrix< number >::const_iterator FilteredMatrix< number >::end (  )  const [inline]

Final iterator.

References FilteredMatrix< VECTOR >::constraints.

template<class VECTOR >
unsigned int FilteredMatrix< VECTOR >::memory_consumption (  )  const [inline]

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.

References FilteredMatrix< VECTOR >::constraints, FilteredMatrix< VECTOR >::matrix, and MemoryConsumption::memory_consumption().

template<class VECTOR >
void FilteredMatrix< VECTOR >::pre_filter ( VECTOR &  v  )  const [inline, private]

Do the pre-filtering step, i.e. zero out those components that belong to constrained degrees of freedom.

References FilteredMatrix< VECTOR >::constraints.

Referenced by FilteredMatrix< VECTOR >::Tvmult(), FilteredMatrix< VECTOR >::Tvmult_add(), FilteredMatrix< VECTOR >::vmult(), and FilteredMatrix< VECTOR >::vmult_add().

template<class VECTOR >
void FilteredMatrix< VECTOR >::post_filter ( const VECTOR &  in,
VECTOR &  out 
) const [inline, private]

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.

References FilteredMatrix< VECTOR >::constraints.

Referenced by FilteredMatrix< VECTOR >::Tvmult(), FilteredMatrix< VECTOR >::Tvmult_add(), FilteredMatrix< VECTOR >::vmult(), and FilteredMatrix< VECTOR >::vmult_add().


Friends And Related Function Documentation

template<class VECTOR>
friend class Accessor [friend]
template<class VECTOR>
friend class FilteredMatrixBlock< VECTOR > [friend]

FilteredMatrixBlock accesses pre_filter() and post_filter().


Member Data Documentation

template<class VECTOR>
bool FilteredMatrix< VECTOR >::expect_constrained_source [private]

Determine, whether multiplications can expect that the source vector has all constrained entries set to zero.

If so, the auxiliary vector can be avoided and memory as well as time can be saved.

We expect this for instance in Newton's method, where the residual already should be zero on constrained nodes. This is, because there is no testfunction in these nodes.

Referenced by FilteredMatrix< VECTOR >::initialize(), FilteredMatrix< VECTOR >::operator=(), FilteredMatrix< VECTOR >::Tvmult(), FilteredMatrix< VECTOR >::Tvmult_add(), FilteredMatrix< VECTOR >::vmult(), and FilteredMatrix< VECTOR >::vmult_add().

template<class VECTOR>
std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > FilteredMatrix< VECTOR >::matrix [private]
template<class VECTOR>
std::vector<IndexValuePair> FilteredMatrix< VECTOR >::constraints [private]

The documentation for this class was generated from the following file:

deal.II documentation generated on Mon Nov 23 22:57:45 2009 by doxygen 1.6.1