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

Inheritance diagram for FilteredMatrix< VectorType >:

Classes
class	Accessor

class	const_iterator

struct	PairComparison

Public Types
using	size_type = types::global_dof_index

using	IndexValuePair = std::pair< size_type, double >

Public Member Functions
Constructors and initialization
	FilteredMatrix ()

	FilteredMatrix (const FilteredMatrix &fm)

template<typename MatrixType >
	FilteredMatrix (const MatrixType &matrix, const bool expect_constrained_source=false)

FilteredMatrix &	operator= (const FilteredMatrix &fm)

template<typename MatrixType >
void	initialize (const MatrixType &m, const bool expect_constrained_source=false)

void	clear ()

Managing constraints
void	add_constraint (const size_type i, const double v)

template<class ConstraintList >
void	add_constraints (const ConstraintList &new_constraints)

void	clear_constraints ()

void	apply_constraints (VectorType &v, const bool matrix_is_symmetric) const

void	apply_constraints (VectorType &v) const

void	vmult (VectorType &dst, const VectorType &src) const

void	Tvmult (VectorType &dst, const VectorType &src) const

void	vmult_add (VectorType &dst, const VectorType &src) const

void	Tvmult_add (VectorType &dst, const VectorType &src) const

Public Member Functions inherited from Subscriptor
	Subscriptor ()

	Subscriptor (const Subscriptor &)

	Subscriptor (Subscriptor &&) noexcept

virtual	~Subscriptor ()

Subscriptor &	operator= (const Subscriptor &)

Subscriptor &	operator= (Subscriptor &&) noexcept

void	subscribe (std::atomic< bool > *const validity, const std::string &identifier="") const

void	unsubscribe (std::atomic< bool > *const validity, const std::string &identifier="") const

unsigned int	n_subscriptions () const

template<typename StreamType >
void	list_subscribers (StreamType &stream) const

void	list_subscribers () const

template<class Archive >
void	serialize (Archive &ar, const unsigned int version)

Iterators
using	const_index_value_iterator = typename std::vector< IndexValuePair >::const_iterator

bool	expect_constrained_source

std::shared_ptr< PointerMatrixBase< VectorType > >	matrix

std::vector< IndexValuePair >	constraints

class	Accessor

class	FilteredMatrixBlock< VectorType >

const_iterator	begin () const

const_iterator	end () const

std::size_t	memory_consumption () const

void	pre_filter (VectorType &v) const

void	post_filter (const VectorType &in, VectorType &out) const

Additional Inherited Members
Static Public Member Functions inherited from Subscriptor
static ::ExceptionBase &	ExcInUse (int arg1, std::string arg2, std::string arg3)

static ::ExceptionBase &	ExcNoSubscriber (std::string arg1, std::string arg2)

Detailed Description

template<typename VectorType>
class FilteredMatrix< VectorType >

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 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 functions 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. Consequently, it is only possible to solve with a matrix to 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 multiple times 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 individual 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 synchronized and thus threadsafe. Consequently, you do not need to serialize calls to vmult or residual .

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

Deprecated:: Use LinearOperator instead. See the documentation of constrained_linear_operator().

Definition at line 199 of file filtered_matrix.h.

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

deal.II/lac/filtered_matrix.h

Classes

Public Types

Public Member Functions

Iterators

Additional Inherited Members