Reference documentation for deal.II version 8.4.1
Linear Operators
Collaboration diagram for Linear Operators:


This browser is not able to show SVG: try Firefox, Chrome, Safari, or Opera instead.

## Classes

class  BlockLinearOperator< Range, Domain >

class  LinearOperator< Range, Domain >

class  PackagedOperation< Range >

## Creation of a BlockLinearOperator

template<typename Range , typename Domain , typename BlockMatrixType >
BlockLinearOperator< Range, Domain > block_operator (const BlockMatrixType &block_matrix)

template<size_t m, size_t n, typename Range , typename Domain >
BlockLinearOperator< Range, Domain > block_operator (const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType >, n >, m > &ops)

template<typename Range , typename Domain , typename BlockMatrixType >
BlockLinearOperator< Range, Domain > block_diagonal_operator (const BlockMatrixType &block_matrix)

template<size_t m, typename Range , typename Domain >
BlockLinearOperator< Range, Domain > block_diagonal_operator (const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType >, m > &ops)

template<size_t m, typename Range , typename Domain >
BlockLinearOperator< Range, Domain > block_diagonal_operator (const LinearOperator< typename Range::BlockType, typename Domain::BlockType > &op)

## Manipulation of a BlockLinearOperator

template<typename Range , typename Domain >
LinearOperator< Domain, Range > block_forward_substitution (const BlockLinearOperator< Range, Domain > &block_operator, const BlockLinearOperator< Domain, Range > &diagonal_inverse)

template<typename Range , typename Domain >
LinearOperator< Domain, Range > block_back_substitution (const BlockLinearOperator< Range, Domain > &block_operator, const BlockLinearOperator< Domain, Range > &diagonal_inverse)

## Vector space operations

template<typename Range , typename Domain >
LinearOperator< Range, Domain > operator+ (const LinearOperator< Range, Domain > &first_op, const LinearOperator< Range, Domain > &second_op)

template<typename Range , typename Domain >
LinearOperator< Range, Domain > operator- (const LinearOperator< Range, Domain > &first_op, const LinearOperator< Range, Domain > &second_op)

template<typename Range , typename Domain >
LinearOperator< Range, Domain > operator* (typename Range::value_type number, const LinearOperator< Range, Domain > &op)

template<typename Range , typename Domain >
LinearOperator< Range, Domain > operator* (const LinearOperator< Range, Domain > &op, typename Domain::value_type number)

## Composition and manipulation of a LinearOperator

template<typename Range , typename Intermediate , typename Domain >
LinearOperator< Range, Domain > operator* (const LinearOperator< Range, Intermediate > &first_op, const LinearOperator< Intermediate, Domain > &second_op)

template<typename Range , typename Domain >
LinearOperator< Domain, Range > transpose_operator (const LinearOperator< Range, Domain > &op)

template<typename Solver , typename Preconditioner >
LinearOperator< typename Solver::vector_type, typename Solver::vector_typeinverse_operator (const LinearOperator< typename Solver::vector_type, typename Solver::vector_type > &op, Solver &solver, const Preconditioner &preconditioner)

## Creation of a LinearOperator

template<typename Range >
LinearOperator< Range, Range > identity_operator (const std::function< void(Range &, bool)> &reinit_vector)

template<typename Range , typename Domain >
LinearOperator< Range, Domain > null_operator (const LinearOperator< Range, Domain > &op)

template<typename Range , typename Domain , typename Matrix >
LinearOperator< Range, Domain > linear_operator (const Matrix &matrix)

template<typename Range , typename Domain , typename OperatorExemplar , typename Matrix >
LinearOperator< Range, Domain > linear_operator (const OperatorExemplar &operator_exemplar, const Matrix &matrix)

## Vector space operations

template<typename Range >
PackagedOperation< Range > operator+ (const PackagedOperation< Range > &first_comp, const PackagedOperation< Range > &second_comp)

template<typename Range >
PackagedOperation< Range > operator- (const PackagedOperation< Range > &first_comp, const PackagedOperation< Range > &second_comp)

template<typename Range >
PackagedOperation< Range > operator* (const PackagedOperation< Range > &comp, typename Range::value_type number)

template<typename Range >
PackagedOperation< Range > operator* (typename Range::value_type number, const PackagedOperation< Range > &comp)

template<typename Range >
PackagedOperation< Range > operator+ (const PackagedOperation< Range > &comp, const Range &offset)

template<typename Range >
PackagedOperation< Range > operator+ (const Range &offset, const PackagedOperation< Range > &comp)

template<typename Range >
PackagedOperation< Range > operator- (const PackagedOperation< Range > &comp, const Range &offset)

template<typename Range >
PackagedOperation< Range > operator- (const Range &offset, const PackagedOperation< Range > &comp)

## Creation of a PackagedOperation object

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation< Range > operator+ (const Range &u, const Range &v)

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation< Range > operator- (const Range &u, const Range &v)

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation< Range > operator* (const Range &u, typename Range::value_type number)

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation< Range > operator* (typename Range::value_type number, const Range &u)

template<typename Range , typename Domain >
PackagedOperation< Range > operator* (const LinearOperator< Range, Domain > &op, const Domain &u)

template<typename Range , typename Domain >
PackagedOperation< Domain > operator* (const Range &u, const LinearOperator< Range, Domain > &op)

template<typename Range , typename Domain >
PackagedOperation< Range > operator* (const LinearOperator< Range, Domain > &op, const PackagedOperation< Domain > &comp)

template<typename Range , typename Domain >
PackagedOperation< Domain > operator* (const PackagedOperation< Range > &comp, const LinearOperator< Range, Domain > &op)

## Creation of a LinearOperator related to the Schur Complement

template<typename Range_1 , typename Domain_1 , typename Range_2 , typename Domain_2 >
LinearOperator< Range_2, Domain_2 > schur_complement (const LinearOperator< Domain_1, Range_1 > &A_inv, const LinearOperator< Range_1, Domain_2 > &B, const LinearOperator< Range_2, Domain_1 > &C, const LinearOperator< Range_2, Domain_2 > &D)

## Creation of PackagedOperation objects related to the Schur Complement

template<typename Range_1 , typename Domain_1 , typename Range_2 >
PackagedOperation< Range_2 > condense_schur_rhs (const LinearOperator< Range_1, Domain_1 > &A_inv, const LinearOperator< Range_2, Domain_1 > &C, const Range_1 &f, const Range_2 &g)

template<typename Range_1 , typename Domain_1 , typename Domain_2 >
PackagedOperation< Domain_1 > postprocess_schur_solution (const LinearOperator< Range_1, Domain_1 > &A_inv, const LinearOperator< Range_1, Domain_2 > &B, const Domain_2 &y, const Range_1 &f)

## Detailed Description

### Linear Operator

If deal.II is configured with C++11 support (i.e., DEAL_II_WITH_CXX11=ON or DEAL_II_WITH_CXX14=ON during configuration) a versatile mechanism for storing the concept of a linear operator is available. (For questions about C++11, see deal.II and the C++11 standard .)

This is done with a LinearOperator class that, like the MatrixType concept, defines a minimal interface for applying a linear operation on a vector.

std::function<void(Range &, const Domain &)> vmult;
std::function<void(Range &, const Domain &)> vmult_add;
std::function<void(Domain &, const Range &)> Tvmult;
std::function<void(Domain &, const Range &)> Tvmult_add;

Thus, such an object can be used as a matrix object in all iterative solver classes, either as a matrix object, or as preconditioner.

The big advantage of the LinearOperator class is that it provides syntactic sugar for complex matrix-vector operations. As an example consider the operation $$(A+k\,B)\,C$$, where $$A$$, $$B$$ and $$C$$ denote (possibly different) SparseMatrix objects. In order to construct a LinearOperator op that performs above computation when applied on a vector, one can write:

double k;
// Setup and assembly...
const auto op_a = linear_operator(A);
const auto op_b = linear_operator(B);
const auto op_c = linear_operator(C);
const auto op = (op_a + k * op_b) * op_c;

Now, op can be used as a matrix object for further computation.

The linear_operator() function can be used to wrap an ordinary matrix or preconditioner object into a LinearOperator. A linear operator can be transposed with transpose_operator(), or inverted by using the inverse_operator() together with an iterative solver.

For objects of type LinearOperator, all vector space operations, i.e., addition and subtraction, scalar multiplication and composition (of compatible linear operators) are implemented:

::LinearOperator<> op_a, op_b;
double k;
// vector space addition, subtraction and scalar multiplication
op_a + op_b;
op_a - op_b;
k * op_a;
op_a * k;
// in-place variants
op_a += op_b;
op_a -= op_b;
op_a *= k;
// operator composition
op_a * op_b;
op_a *= op_b; // If op_b is an endomorphism of the domain space of op_a

block_operator() and block_diagonal_operator() provide further encapsulation of individual linear operators into blocked linear operator variants.

Note
The LinearOperator facility obsoletes some of the derived matrix classes, such as BlockDiagonalMatrix, IterativeInverse, ProductMatrix, ScaledMatrix, ProductSparseMatrix, InverseMatrixRichardson, SchurMatrix, ShiftedMatrix, ShiftedMatrixGeneralized, TransposeMatrix

### Packaged Operation

An application of a LinearOperator object to a vector via operator* yields a PackagedOperation object that stores this computation.

The PackagedOperation class allows lazy evaluation of expressions involving vectors and linear operators. This is done by storing the computational expression and only performing the computation when either the object is implicitly converted to a vector object, or PackagedOperation::apply() (or PackagedOperation::apply_add()) is invoked by hand. This avoids unnecessary temporary storage of intermediate results.

As an example consider the addition of multiple vectors:

::Vector<double> a, b, c, d;
// ..
::Vector<double> result = a + b - c + d;

Converting the PackagedOperation a + b - c + d to a vector results in code equivalent to the following code

::Vector<double> a, b, c, d;
// ..
::Vector<double> result = a;
result += b;
result -= c;
result += d;

that avoids any intermediate storage. As a second example (involving a LinearOperator object) consider the computation of a residual $$b-Ax$$:

// ..
const auto op_a = linear_operator(A);
::Vector<double> residual = b - op_a * x;

Here, the expression b - op_a * x results again in an object of type PackagedOperation that stores the sequence of operations that should be performed using the two vectors and the linear operator. Converting the expression to a vector (as happens here with the assignment to the vector residual) executes the computation (see the following note).

Note
Lazy evaluation of a computational expression necessarily involves references to the underlying vector and matrix objects. For example, the creation of a residual_expr object
auto residual_expr = b - op_a * x;
stores the computational expression of the residual with references to the vector b and matrix A. It does not perform any computation at this point. In particular, if b or A are changed after the creation of residual_expr every subsequent evaluation of the expression is performed with the new values
auto residual_expr = b - op_a * x;
residual_expr.apply(tmp); // tmp is a Vector<double>
// modify b, or A
residual_expr.apply(tmp2); // tmp2 is a Vector<double>
// tmp and tmp2 are different
Thus, as a safeguard, if you want to compute the result of an expression right away, always explicitly use a vector type on the left side (and not auto):
Vector<double> residual = b - op_a * x; // computes the residual at this point

## Function Documentation

template<typename Range , typename Domain , typename BlockMatrixType >
 BlockLinearOperator< Range, Domain > block_operator ( const BlockMatrixType & block_matrix )
related

A function that encapsulates a block_matrix into a BlockLinearOperator.

All changes made on the block structure and individual blocks of block_matrix after the creation of the BlockLinearOperator object are reflected by the operator object.

Definition at line 387 of file block_linear_operator.h.

template<size_t m, size_t n, typename Range , typename Domain >
 BlockLinearOperator< Range, Domain > block_operator ( const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType >, n >, m > & ops )
related

A variant of above function that encapsulates a given collection ops of LinearOperators into a block structure. Here, it is assumed that Range and Domain are blockvectors, i.e., derived from BlockVectorBase. The individual linear operators in ops must act on the underlying vector type of the block vectors, i.e., on Domain::BlockType yielding a result in Range::BlockType.

The list ops is best passed as an initializer list. Consider for example a linear operator block (acting on Vector<double>)

op_a00 | op_a01
|
---------------
|
op_a10 | op_a11

The corresponding block_operator invocation takes the form

block_operator<2, 2, BlockVector<double>>({op_a00, op_a01, op_a10, op_a11});

Definition at line 450 of file block_linear_operator.h.

template<typename Range , typename Domain , typename BlockMatrixType >
 BlockLinearOperator< Range, Domain > block_diagonal_operator ( const BlockMatrixType & block_matrix )
related

This function extracts the diagonal blocks of block_matrix (either a block matrix type or a BlockLinearOperator) and creates a BlockLinearOperator with the diagonal. Off-diagonal elements are initialized as null_operator (with correct reinit_range_vector and reinit_domain_vector methods).

All changes made on the individual diagonal blocks of block_matrix after the creation of the BlockLinearOperator object are reflected by the operator object.

Definition at line 502 of file block_linear_operator.h.

template<size_t m, typename Range , typename Domain >
 BlockLinearOperator< Range, Domain > block_diagonal_operator ( const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType >, m > & ops )
related

A variant of above function that builds up a block diagonal linear operator from an array ops of diagonal elements (off-diagonal blocks are assumed to be 0).

The list ops is best passed as an initializer list. Consider for example a linear operator block (acting on Vector<double>) diag(op_a0, op_a1, ..., op_am). The corresponding block_operator invocation takes the form

block_diagonal_operator<m, BlockVector<double>>({op_00, op_a1, ..., op_am});

Definition at line 558 of file block_linear_operator.h.

template<size_t m, typename Range , typename Domain >
 BlockLinearOperator< Range, Domain > block_diagonal_operator ( const LinearOperator< typename Range::BlockType, typename Domain::BlockType > & op )
related

A variant of above function that only takes a single LinearOperator argument op and creates a blockdiagonal linear operator with m copies of it.

Definition at line 602 of file block_linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Domain, Range > block_forward_substitution ( const BlockLinearOperator< Range, Domain > & block_operator, const BlockLinearOperator< Domain, Range > & diagonal_inverse )
related

This function implements forward substitution to invert a lower block triangular matrix. As arguments, it takes a BlockLinearOperator block_operator representing a block lower triangular matrix, as well as a BlockLinearOperator diagonal_inverse representing inverses of diagonal blocks of block_operator.

Let us assume we have a linear system with the following block structure:

A00 x0 + ... = y0
A01 x0 + A11 x1 + ... = y1
... ...
A0n x0 + A1n x1 + ... + Ann xn = yn

First of all, x0 = A00^-1 y0. Then, we can use x0 to recover x1:

x1 = A11^-1 ( y1 - A01 x0 )

and therefore:

xn = Ann^-1 ( yn - A0n x0 - ... - A(n-1)n x(n-1) )
Note
We are not using all blocks of the BlockLinearOperator arguments: Just the lower triangular block matrix of block_operator is used as well as the diagonal of diagonal_inverse.

Definition at line 660 of file block_linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Domain, Range > block_back_substitution ( const BlockLinearOperator< Range, Domain > & block_operator, const BlockLinearOperator< Domain, Range > & diagonal_inverse )
related

This function implements back substitution to invert an upper block triangular matrix. As arguments, it takes a BlockLinearOperator block_operator representing an upper block triangular matrix, as well as a BlockLinearOperator diagonal_inverse representing inverses of diagonal blocks of block_operator.

Let us assume we have a linear system with the following block structure:

A00 x0 + A01 x1 + ... + A0n xn = yn
A11 x1 + ... = y1
... ..
Ann xn = yn

First of all, xn = Ann^-1 yn. Then, we can use xn to recover x(n-1):

x(n-1) = A(n-1)(n-1)^-1 ( y(n-1) - A(n-1)n x(n-1) )

and therefore:

x0 = A00^-1 ( y0 - A0n xn - ... - A01 x1 )
Note
We are not using all blocks of the BlockLinearOperator arguments: Just the upper triangular block matrix of block_operator is used as well as the diagonal of diagonal_inverse.

Definition at line 772 of file block_linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Range, Domain > operator+ ( const LinearOperator< Range, Domain > & first_op, const LinearOperator< Range, Domain > & second_op )
related

Addition of two linear operators first_op and second_op given by $$(\text{first\_op}+\text{second\_op})x := \text{first\_op}(x) + \text{second\_op}(x)$$

Definition at line 319 of file linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Range, Domain > operator- ( const LinearOperator< Range, Domain > & first_op, const LinearOperator< Range, Domain > & second_op )
related

Subtraction of two linear operators first_op and second_op given by $$(\text{first\_op}-\text{second\_op})x := \text{first\_op}(x) - \text{second\_op}(x)$$

Definition at line 380 of file linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Range, Domain > operator* ( typename Range::value_type number, const LinearOperator< Range, Domain > & op )
related

Scalar multiplication of a ScalarOperator object op with number from the left.

The Domain and Range types must implement the following operator*= member functions accepting the appropriate scalar Number type for rescaling:

Domain & operator *=(Domain::value_type);
Range & operator *=(Range::value_type);

Definition at line 418 of file linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Range, Domain > operator* ( const LinearOperator< Range, Domain > & op, typename Domain::value_type number )
related

Scalar multiplication of a ScalarOperator object from the right.

The Domain and Range types must implement the following operator*= member functions for rescaling:

Domain & operator *=(Domain::value_type);
Range & operator *=(Range::value_type);

Definition at line 488 of file linear_operator.h.

template<typename Range , typename Intermediate , typename Domain >
 LinearOperator< Range, Domain > operator* ( const LinearOperator< Range, Intermediate > & first_op, const LinearOperator< Intermediate, Domain > & second_op )
related

Composition of two linear operators first_op and second_op given by $$(\text{first\_op}*\text{second\_op})x := \text{first\_op}(\text{second\_op}(x))$$

Definition at line 517 of file linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Domain, Range > transpose_operator ( const LinearOperator< Range, Domain > & op )
related

Returns the transpose linear operations of op.

Definition at line 594 of file linear_operator.h.

template<typename Solver , typename Preconditioner >
 LinearOperator< typename Solver::vector_type, typename Solver::vector_type > inverse_operator ( const LinearOperator< typename Solver::vector_type, typename Solver::vector_type > & op, Solver & solver, const Preconditioner & preconditioner )
related

Returns an object representing the inverse of the LinearOperator op.

The function takes references solver and preconditioner to an iterative solver and a preconditioner that are used in the vmult and Tvmult implementations of the LinearOperator object.

The LinearOperator object that is created stores a reference to solver and preconditioner. Thus, both objects must remain a valid reference for the whole lifetime of the LinearOperator object. Internal data structures of the solver object will be modified upon invocation of vmult or Tvmult.

Definition at line 629 of file linear_operator.h.

template<typename Range >
 LinearOperator< Range, Range > identity_operator ( const std::function< void(Range &, bool)> & reinit_vector )
related

Returns a LinearOperator that is the identity of the vector space Range.

The function takes an std::function object reinit_vector as an argument to initialize the reinit_range_vector and reinit_domain_vector objects of the LinearOperator object.

Definition at line 701 of file linear_operator.h.

template<typename Range , typename Domain >
 LinearOperator< Range, Domain > null_operator ( const LinearOperator< Range, Domain > & op )
related

Returns a nulled variant of the LinearOperator op, i.e. with optimized LinearOperator::vmult, LinearOperator::vmult_add, etc. functions and with LinearOperator::is_null_operator set to true.

Definition at line 743 of file linear_operator.h.

template<typename Range , typename Domain , typename Matrix >
 LinearOperator< Range, Domain > linear_operator ( const Matrix & matrix )
related

A function that encapsulates generic matrix objects that act on a compatible Vector type into a LinearOperator. The LinearOperator object that is created stores a reference to the matrix object. Thus, matrix must remain a valid reference for the whole lifetime of the LinearOperator object.

All changes made on matrix after the creation of the LinearOperator object are reflected by the operator object. For example, it is a valid procedure to first create a LinearOperator and resize, reassemble the matrix later.

The Matrix class in question must provide the following minimal interface:

class Matrix
{
public:
// (type specific) information how to create a Range and Domain vector
// with appropriate size and internal layout
// Application of matrix to vector src, writes the result into dst.
vmult(Range &dst, const Domain &src);
// Application of the transpose of matrix to vector src, writes the
// result into dst. (Depending on the usage of the linear operator
// class this can be a dummy implementation throwing an error.)
Tvmult(Range &dst, const Domain &src);
};

The following (optional) interface is used if available:

class Matrix
{
public:
// Application of matrix to vector src, adds the result to dst.
vmult_add(Range &dst, const Domain &src);
// Application of the transpose of matrix to vector src, adds the
// result to dst.
Tvmult_add(Range &dst, const Domain &src);
};

If the Matrix does not provide vmult_add and Tvmult_add, they are implemented in terms of vmult and Tvmult (requiring intermediate storage).

Definition at line 1061 of file linear_operator.h.

template<typename Range , typename Domain , typename OperatorExemplar , typename Matrix >
 LinearOperator< Range, Domain > linear_operator ( const OperatorExemplar & operator_exemplar, const Matrix & matrix )
related

Variant of above function that takes an operator object operator_exemplar as an additional reference. This object is used to populate the reinit_domain_vector and reinit_range_vector function objects. The reference matrix is used to construct vmult, Tvmult, etc.

This variant can, for example, be used to encapsulate preconditioners (that typically do not expose any information about the underlying matrix).

Definition at line 1088 of file linear_operator.h.

template<typename Range >
 PackagedOperation< Range > operator+ ( const PackagedOperation< Range > & first_comp, const PackagedOperation< Range > & second_comp )
related

Addition of two PackagedOperation objects first_comp and second_comp given by vector space addition of the corresponding results.

Definition at line 289 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator- ( const PackagedOperation< Range > & first_comp, const PackagedOperation< Range > & second_comp )
related

Subtraction of two PackagedOperation objects first_comp and second_comp given by vector space addition of the corresponding results.

Definition at line 324 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator* ( const PackagedOperation< Range > & comp, typename Range::value_type number )
related

Scalar multiplication of a PackagedOperation objects comp with a scalar number given by a scaling PackagedOperation result with number.

Definition at line 362 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator* ( typename Range::value_type number, const PackagedOperation< Range > & comp )
related

Scalar multiplication of a PackagedOperation objects comp with a scalar number given by a scaling PackagedOperation result with number.

Definition at line 410 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator+ ( const PackagedOperation< Range > & comp, const Range & offset )
related

Add a constant offset (of the Range space) to the result of a PackagedOperation.

Definition at line 425 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator+ ( const Range & offset, const PackagedOperation< Range > & comp )
related

Add a constant offset (of the Range space) to the result of a PackagedOperation.

Definition at line 440 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator- ( const PackagedOperation< Range > & comp, const Range & offset )
related

Subtract a constant offset (of the Range space) from the result of a PackagedOperation.

Definition at line 455 of file packaged_operation.h.

template<typename Range >
 PackagedOperation< Range > operator- ( const Range & offset, const PackagedOperation< Range > & comp )
related

Subtract a computational result from a constant offset (of the Range space). The result is a PackagedOperation object that applies this computation.

Definition at line 472 of file packaged_operation.h.

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
 PackagedOperation< Range > operator+ ( const Range & u, const Range & v )
related

Create a PackagedOperation object that stores the addition of two vectors.

The PackagedOperation object that is created stores a reference to u and v. Thus, the vectors must remain valid references for the whole lifetime of the PackagedOperation object. All changes made on u or v after the creation of the PackagedOperation object are reflected by the operator object.

Definition at line 528 of file packaged_operation.h.

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
 PackagedOperation< Range > operator- ( const Range & u, const Range & v )
related

Create a PackagedOperation object that stores the subtraction of two vectors.

The PackagedOperation object that is created stores a reference to u and v. Thus, the vectors must remain valid references for the whole lifetime of the PackagedOperation object. All changes made on u or v after the creation of the PackagedOperation object are reflected by the operator object.

Definition at line 573 of file packaged_operation.h.

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
 PackagedOperation< Range > operator* ( const Range & u, typename Range::value_type number )
related

Create a PackagedOperation object that stores the scaling of a vector with a number.

The PackagedOperation object that is created stores a reference to u. Thus, the vectors must remain valid references for the whole lifetime of the PackagedOperation object. All changes made on u or v after the creation of the PackagedOperation object are reflected by the operator object.

Definition at line 617 of file packaged_operation.h.

template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
 PackagedOperation< Range > operator* ( typename Range::value_type number, const Range & u )
related

Create a PackagedOperation object that stores the scaling of a vector with a number.

The PackagedOperation object that is created stores a reference to u. Thus, the vectors must remain valid references for the whole lifetime of the PackagedOperation object. All changes made on u or v after the creation of the PackagedOperation object are reflected by the operator object.

Definition at line 640 of file packaged_operation.h.

template<typename Range , typename Domain >
 PackagedOperation< Range > operator* ( const LinearOperator< Range, Domain > & op, const Domain & u )
related

Create a PackagedOperation object from a LinearOperator and a reference to a vector u of the Domain space. The object stores the PackagedOperation $$\text{op} \,u$$ (in matrix notation). return (return_add) are implemented with vmult(__1,u) (vmult_add(__1,u)).

The PackagedOperation object that is created stores a reference to u. Thus, the vector must remain a valid reference for the whole lifetime of the PackagedOperation object. All changes made on u after the creation of the PackagedOperation object are reflected by the operator object.

Definition at line 665 of file packaged_operation.h.

template<typename Range , typename Domain >
 PackagedOperation< Domain > operator* ( const Range & u, const LinearOperator< Range, Domain > & op )
related

Create a PackagedOperation object from a LinearOperator and a reference to a vector u of the Range space. The object stores the PackagedOperation $$\text{op}^T \,u$$ (in matrix notation). return (return_add) are implemented with Tvmult(__1,u) (Tvmult_add(__1,u)).

The PackagedOperation object that is created stores a reference to u. Thus, the vector must remain a valid reference for the whole lifetime of the PackagedOperation object. All changes made on u after the creation of the PackagedOperation object are reflected by the operator object.

Definition at line 707 of file packaged_operation.h.

template<typename Range , typename Domain >
 PackagedOperation< Range > operator* ( const LinearOperator< Range, Domain > & op, const PackagedOperation< Domain > & comp )
related

Composition of a PackagedOperation object with a LinearOperator. The object stores the computation $$\text{op} \,comp$$ (in matrix notation).

Definition at line 741 of file packaged_operation.h.

template<typename Range , typename Domain >
 PackagedOperation< Domain > operator* ( const PackagedOperation< Range > & comp, const LinearOperator< Range, Domain > & op )
related

Composition of a PackagedOperation object with a LinearOperator. The object stores the computation $$\text{op}^T \,comp$$ (in matrix notation).

Definition at line 791 of file packaged_operation.h.

template<typename Range_1 , typename Domain_1 , typename Range_2 , typename Domain_2 >
 LinearOperator< Range_2, Domain_2 > schur_complement ( const LinearOperator< Domain_1, Range_1 > & A_inv, const LinearOperator< Range_1, Domain_2 > & B, const LinearOperator< Range_2, Domain_1 > & C, const LinearOperator< Range_2, Domain_2 > & D )
related

Returns a LinearOperator that performs the operations associated with the Schur complement. There are two additional helper functions, condense_schur_rhs() and postprocess_schur_solution(), that are likely necessary to be used in order to perform any useful tasks in linear algebra with this operator.

We construct the definition of the Schur complement in the following way:

Consider a general system of linear equations that can be decomposed into two major sets of equations:

\begin{eqnarray*} \mathbf{K}\mathbf{d} = \mathbf{f} \quad \Rightarrow\quad \left(\begin{array}{cc} A & B \\ C & D \end{array}\right) \left(\begin{array}{cc} x \\ y \end{array}\right) = \left(\begin{array}{cc} f \\ g \end{array}\right), \end{eqnarray*}

where $$A,B,C,D$$ represent general subblocks of the matrix $$\mathbf{K}$$ and, similarly, general subvectors of $$\mathbf{d},\mathbf{f}$$ are given by $$x,y,f,g$$ .

This is equivalent to the following two statements:

\begin{eqnarray*} (1) \quad Ax + By &=& f \\ (2) \quad Cx + Dy &=& g \quad . \end{eqnarray*}

Assuming that $$A,D$$ are both square and invertible, we could then perform one of two possible substitutions,

\begin{eqnarray*} (3) \quad x &=& A^{-1}(f - By) \quad \text{from} \quad (1) \\ (4) \quad y &=& D^{-1}(g - Cx) \quad \text{from} \quad (2) , \end{eqnarray*}

which amount to performing block Gaussian elimination on this system of equations.

For the purpose of the current implementation, we choose to substitute (3) into (2)

\begin{eqnarray*} C \: A^{-1}(f - By) + Dy &=& g \\ -C \: A^{-1} \: By + Dy &=& g - C \: A^{-1} \: f \quad . \end{eqnarray*}

This leads to the result

$(5) \quad (D - C\: A^{-1} \:B)y = g - C \: A^{-1} f \quad \Rightarrow \quad Sy = g'$

with $$S = (D - C\: A^{-1} \:B)$$ being the Schur complement and the modified right-hand side vector $$g' = g - C \: A^{-1} f$$ arising from the condensation step. Note that for this choice of $$S$$, submatrix $$D$$ need not be invertible and may thus be the null matrix. Ideally $$A$$ should be well-conditioned.

So for any arbitrary vector $$a$$, the Schur complement performs the following operation:

$(6) \quad Sa = (D - C \: A^{-1} \: B)a$

A typical set of steps needed the solve a linear system (1),(2) would be:

1. Define the inverse matrix A_inv (using inverse_operator()).
2. Define the Schur complement $$S$$ (using schur_complement()).
3. Define iterative inverse matrix $$S^{-1}$$ such that (6) holds. It is necessary to use a solver with a preconditioner to compute the approximate inverse operation of $$S$$ since we never compute $$S$$ directly, but rather the result of its operation. To achieve this, one may again use the inverse_operator() in conjunction with the Schur complement that we've just constructed. Observe that the both $$S$$ and its preconditioner operate over the same space as $$D$$.
4. Perform pre-processing step on the RHS of (5) using condense_schur_rhs():

$g' = g - C \: A^{-1} \: f$

5. Solve for $$y$$ in (5):

$y = S^{-1} g'$

6. Perform the post-processing step from (3) using postprocess_schur_solution():

$x = A^{-1} (f - By)$

An illustration of typical usage of this operator for a fully coupled system is given below.

#include<deal.II/lac/schur_complement.h>
// Given BlockMatrix K and BlockVectors d,F
// Decomposition of tangent matrix
const auto A = linear_operator(K.block(0,0));
const auto B = linear_operator(K.block(0,1));
const auto C = linear_operator(K.block(1,0));
const auto D = linear_operator(K.block(1,1));
// Decomposition of solution vector
auto x = d.block(0);
auto y = d.block(1);
// Decomposition of RHS vector
auto f = F.block(0);
auto g = F.block(1);
// Construction of inverse of Schur complement
const auto prec_A = PreconditionSelector<...>(A);
const auto A_inv = inverse_operator<...>(A,prec_A);
const auto S = schur_complement(A_inv,B,C,D);
const auto S_prec = PreconditionSelector<...>(D); // D and S operate on same space
const auto S_inv = inverse_operator<...>(S,...,prec_S);
// Solve reduced block system
auto rhs = condense_schur_rhs (A_inv,C,f,g); // PackagedOperation that represents the condensed form of g
y = S_inv * rhs; // Solve for y
x = postprocess_schur_solution (A_inv,B,y,f); // Compute x using resolved solution y

In the above example, the preconditioner for $$S$$ was defined as the preconditioner for $$D$$, which is valid since they operate on the same space. However, if $$D$$ and $$S$$ are too dissimilar, then this may lead to a large number of solver iterations as $$\text{prec}(D)$$ is not a good approximation for $$S^{-1}$$.

A better preconditioner in such a case would be one that provides a more representative approximation for $$S^{-1}$$. One approach is shown in step-22, where $$D$$ is the null matrix and the preconditioner for $$S^{-1}$$ is derived from the mass matrix over this space.

From another viewpoint, a similar result can be achieved by first constructing an object that represents an approximation for $$S$$ wherein expensive operation, namely $$A^{-1}$$, is approximated. Thereafter we construct the approximate inverse operator $$\tilde{S}^{-1}$$ which is then used as the preconditioner for computing $$S^{-1}$$.

// Construction of approximate inverse of Schur complement
const auto A_inv_approx = linear_operator(preconditioner_A);
const auto S_approx = schur_complement(A_inv_approx,B,C,D);
const auto S_approx_prec = PreconditionSelector<...>(D); // D and S_approx operate on same space
const auto S_inv_approx = inverse_operator(S_approx,...,S_approx_prec); // Inner solver: Typically limited to few iterations using IterationNumberControl
// Construction of exact inverse of Schur complement
const auto S = schur_complement(A_inv,B,C,D);
const auto S_inv = inverse_operator(S,...,S_inv_approx); // Outer solver
// Solve reduced block system
auto rhs = condense_schur_rhs (A_inv,C,f,g);
y = S_inv * rhs; // Solve for y
x = postprocess_schur_solution (A_inv,B,y,f);

Note that due to the construction of S_inv_approx and subsequently S_inv, there are a pair of nested iterative solvers which could collectively consume a lot of resources. Therefore care should be taken in the choices leading to the construction of the iterative inverse_operators. One might consider the use of a IterationNumberControl (or a similar mechanism) to limit the number of inner solver iterations. This controls the accuracy of the approximate inverse operation $$\tilde{S}^{-1}$$ which acts only as the preconditioner for $$S^{-1}$$. Furthermore, the preconditioner to $$\tilde{S}^{-1}$$, which in this example is $$\text{prec}(D)$$, should ideally be computationally inexpensive.

However, if an iterative solver based on IterationNumberControl is used as a preconditioner then the preconditioning operation is not a linear operation. Here a flexible solver like SolverFGMRES (flexible GMRES) is best employed as an outer solver in order to deal with the variable behaviour of the preconditioner. Otherwise the iterative solver can stagnate somewhere near the tolerance of the preconditioner or generally behave erratically. Alternatively, using a ReductionControl would ensure that the preconditioner always solves to the same tolerance, thereby rendering its behaviour constant.

Further examples of this functionality can be found in the test-suite, such as tests/lac/schur_complement_01.cc . The solution of a multi- component problem (namely step-22) using the schur_complement can be found in tests/lac/schur_complement_03.cc .

Block (linear algebra)

Definition at line 229 of file schur_complement.h.

template<typename Range_1 , typename Domain_1 , typename Range_2 >
 PackagedOperation< Range_2 > condense_schur_rhs ( const LinearOperator< Range_1, Domain_1 > & A_inv, const LinearOperator< Range_2, Domain_1 > & C, const Range_1 & f, const Range_2 & g )
related

For the system of equations

\begin{eqnarray*} Ax + By &=& f \\ Cx + Dy &=& g \quad , \end{eqnarray*}

this operation performs the pre-processing (condensation) step on the RHS subvector g so that the Schur complement can be used to solve this system of equations. More specifically, it produces an object that represents the condensed form of the subvector g, namely

$g' = g - C \: A^{-1} \: f$

Block (linear algebra)

Definition at line 366 of file schur_complement.h.

template<typename Range_1 , typename Domain_1 , typename Domain_2 >
 PackagedOperation< Domain_1 > postprocess_schur_solution ( const LinearOperator< Range_1, Domain_1 > & A_inv, const LinearOperator< Range_1, Domain_2 > & B, const Domain_2 & y, const Range_1 & f )
related

For the system of equations

\begin{eqnarray*} Ax + By &=& f \\ Cx + Dy &=& g \quad , \end{eqnarray*}

this operation performs the post-processing step of the Schur complement to solve for the second subvector x once subvector y is known, with the result that

$x = A^{-1}(f - By)$