Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Classes
Linear Operators
Collaboration diagram for Linear Operators:

Classes

class  BlockLinearOperator< Range, Domain, BlockPayload >
 
class  internal::BlockLinearOperatorImplementation::EmptyBlockPayload< PayloadBlockType >
 
class  LinearOperator< Range, Domain, Payload >
 
class  internal::LinearOperatorImplementation::EmptyPayload
 
class  PackagedOperation< Range >
 

Creation of a BlockLinearOperator

template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>, typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > block_operator (const BlockMatrixType &matrix)
 
template<std::size_t m, std::size_t n, typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator< Range, Domain, BlockPayload > block_operator (const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &)
 
template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>, typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator (const BlockMatrixType &block_matrix)
 
template<std::size_t m, typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator (const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &)
 
template<std::size_t m, typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator (const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &op)
 
template<typename Range , typename Domain , typename BlockPayload , typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > block_operator (const BlockMatrixType &block_matrix)
 
template<std::size_t m, std::size_t n, typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > block_operator (const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &ops)
 
template<typename Range , typename Domain , typename BlockPayload , typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator (const BlockMatrixType &block_matrix)
 
template<std::size_t m, typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator (const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &ops)
 
template<std::size_t m, typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator (const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &op)
 

Manipulation of a BlockLinearOperator

template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution (const BlockLinearOperator< Range, Domain, BlockPayload > &, const BlockLinearOperator< Domain, Range, BlockPayload > &)
 
template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution (const BlockLinearOperator< Range, Domain, BlockPayload > &, const BlockLinearOperator< Domain, Range, BlockPayload > &)
 
template<typename Range , typename Domain , typename BlockPayload >
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution (const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
 
template<typename Range , typename Domain , typename BlockPayload >
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution (const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
 

Vector space operations

template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator+ (const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator- (const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* (typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* (const LinearOperator< Range, Domain, Payload > &op, typename Domain::value_type number)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator+ (const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator- (const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* (typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* (const LinearOperator< Range, Domain, Payload > &op, typename Domain::value_type number)
 

Composition and manipulation of a LinearOperator

template<typename Range , typename Intermediate , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* (const LinearOperator< Range, Intermediate, Payload > &first_op, const LinearOperator< Intermediate, Domain, Payload > &second_op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Domain, Range, Payload > transpose_operator (const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Payload , typename Solver , typename Preconditioner , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
 
template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const LinearOperator< Range, Domain, Payload > &preconditioner)
 
template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver)
 
template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const PreconditionIdentity &)
 
template<typename Range , typename Intermediate , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* (const LinearOperator< Range, Intermediate, Payload > &first_op, const LinearOperator< Intermediate, Domain, Payload > &second_op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Domain, Range, Payload > transpose_operator (const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Payload , typename Solver , typename Preconditioner , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
 
template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const LinearOperator< Range, Domain, Payload > &preconditioner)
 
template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver)
 
template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator (const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const PreconditionIdentity &)
 

Creation of a LinearOperator

template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > identity_operator (const std::function< void(Range &, bool)> &reinit_vector)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > identity_operator (const LinearOperator< Range, Domain, Payload > &)
 
template<typename Range = Vector<double>, typename Domain = Range, typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Domain, Payload > null_operator (const LinearOperator< Range, Domain, Payload > &)
 
template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > mean_value_filter (const std::function< void(Range &, bool)> &reinit_vector)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > mean_value_filter (const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range = Vector<double>, typename Domain = Range, typename Payload = internal::LinearOperatorImplementation::EmptyPayload, typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator (const Matrix &)
 
template<typename Range = Vector<double>, typename Domain = Range, typename Payload = internal::LinearOperatorImplementation::EmptyPayload, typename OperatorExemplar , typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator (const OperatorExemplar &, const Matrix &)
 
template<typename Range , typename Domain , typename Payload , typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator (const LinearOperator< Range, Domain, Payload > &operator_exemplar, const Matrix &matrix)
 
template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > identity_operator (const std::function< void(Range &, bool)> &reinit_vector)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > identity_operator (const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > null_operator (const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > mean_value_filter (const std::function< void(Range &, bool)> &reinit_vector)
 
template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > mean_value_filter (const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Domain , typename Payload , typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator (const Matrix &matrix)
 
template<typename Range , typename Domain , typename Payload , typename OperatorExemplar , typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator (const OperatorExemplar &operator_exemplar, const Matrix &matrix)
 
template<typename Range , typename Domain , typename Payload , typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator (const LinearOperator< Range, Domain, Payload > &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)
 
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< internal::PackagedOperationImplementation::has_vector_interface< Range>::type::value>::type>
PackagedOperation< Range > operator+ (const Range &u, const Range &v)
 
template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::has_vector_interface< Range>::type::value>::type>
PackagedOperation< Range > operator- (const Range &u, const Range &v)
 
template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::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< internal::PackagedOperationImplementation::has_vector_interface< Range>::type::value>::type>
PackagedOperation< Range > operator* (typename Range::value_type number, const Range &u)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Range > operator* (const LinearOperator< Range, Domain, Payload > &op, const Domain &u)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Domain > operator* (const Range &u, const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Range > operator* (const LinearOperator< Range, Domain, Payload > &op, const PackagedOperation< Domain > &comp)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Domain > operator* (const PackagedOperation< Range > &comp, const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::has_vector_interface< Range>::type::value>::type>
PackagedOperation< Range > operator+ (const Range &u, const Range &v)
 
template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::has_vector_interface< Range>::type::value>::type>
PackagedOperation< Range > operator- (const Range &u, const Range &v)
 
template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::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< internal::PackagedOperationImplementation::has_vector_interface< Range>::type::value>::type>
PackagedOperation< Range > operator* (typename Range::value_type number, const Range &u)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Range > operator* (const LinearOperator< Range, Domain, Payload > &op, const Domain &u)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Domain > operator* (const Range &u, const LinearOperator< Range, Domain, Payload > &op)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Range > operator* (const LinearOperator< Range, Domain, Payload > &op, const PackagedOperation< Domain > &comp)
 
template<typename Range , typename Domain , typename Payload >
PackagedOperation< Domain > operator* (const PackagedOperation< Range > &comp, const LinearOperator< Range, Domain, Payload > &op)
 

Creation of a LinearOperator related to the Schur Complement

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

Creation of PackagedOperation objects related to the Schur Complement

template<typename Range_1 , typename Domain_1 , typename Range_2 , typename Payload >
PackagedOperation< Range_2 > condense_schur_rhs (const LinearOperator< Range_1, Domain_1, Payload > &A_inv, const LinearOperator< Range_2, Domain_1, Payload > &C, const Range_1 &f, const Range_2 &g)
 
template<typename Range_1 , typename Domain_1 , typename Domain_2 , typename Payload >
PackagedOperation< Domain_1 > postprocess_schur_solution (const LinearOperator< Range_1, Domain_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const Domain_2 &y, const Range_1 &f)
 
template<typename Range_1 , typename Domain_1 , typename Range_2 , typename Payload >
PackagedOperation< Range_2 > condense_schur_rhs (const LinearOperator< Range_1, Domain_1, Payload > &A_inv, const LinearOperator< Range_2, Domain_1, Payload > &C, const Range_1 &f, const Range_2 &g)
 
template<typename Range_1 , typename Domain_1 , typename Domain_2 , typename Payload >
PackagedOperation< Domain_1 > postprocess_schur_solution (const LinearOperator< Range_1, Domain_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const Domain_2 &y, const Range_1 &f)
 

Detailed Description

Linear Operator

deal.II includes support for describing linear transformations in a very general way. 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.

The step-20 tutorial program has a detailed usage example of the LinearOperator class.

Note
As explained below, when using LinearOperator as res = op_a*x a PackagedOperation class instance is generated behind-the-curtains. Consequently, the user program has to include header files for both classes for compilation to be successful. In an attempt to make easier the decision of which headers to include in what circumstances and to prevent hidden templates-related compiler errors, all headers relevant to LinearOperator are grouped in the <deal.II/lac/linear_operator_tools.h> header file.

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> 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> 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

The step-20 tutorial program has a detailed usage example of the PackagedOperation class.

Function Documentation

◆ block_operator() [1/4]

template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>, typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > block_operator ( const BlockMatrixType &  block_matrix)

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 622 of file block_linear_operator.h.

◆ block_operator() [2/4]

template<std::size_t m, std::size_t n, typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator< Range, Domain, BlockPayload > block_operator ( const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &  ops)

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 block vectors, 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 689 of file block_linear_operator.h.

◆ block_diagonal_operator() [1/6]

template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>, typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator ( const BlockMatrixType &  block_matrix)

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 742 of file block_linear_operator.h.

◆ block_diagonal_operator() [2/6]

template<std::size_t m, typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator ( const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &  ops)

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 798 of file block_linear_operator.h.

◆ block_diagonal_operator() [3/6]

template<std::size_t m, typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator ( const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &  op)

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 847 of file block_linear_operator.h.

◆ block_forward_substitution() [1/2]

template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution ( const BlockLinearOperator< Range, Domain, BlockPayload > &  block_operator,
const BlockLinearOperator< Domain, Range, BlockPayload > &  diagonal_inverse 
)

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 909 of file block_linear_operator.h.

◆ block_back_substitution() [1/2]

template<typename Range = BlockVector<double>, typename Domain = Range, typename BlockPayload = internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution ( const BlockLinearOperator< Range, Domain, BlockPayload > &  block_operator,
const BlockLinearOperator< Domain, Range, BlockPayload > &  diagonal_inverse 
)

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 1024 of file block_linear_operator.h.

◆ operator+() [1/10]

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

Addition of two linear operators first_op and second_op given by \((\mathrm{first\_op}+\mathrm{second\_op})x \dealcoloneq \mathrm{first\_op}(x) + \mathrm{second\_op}(x)\)

Definition at line 380 of file linear_operator.h.

◆ operator-() [1/10]

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

Subtraction of two linear operators first_op and second_op given by \((\mathrm{first\_op}-\mathrm{second\_op})x \dealcoloneq \mathrm{first\_op}(x) - \mathrm{second\_op}(x)\)

Definition at line 439 of file linear_operator.h.

◆ operator*() [1/22]

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

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 477 of file linear_operator.h.

◆ operator*() [2/22]

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

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 544 of file linear_operator.h.

◆ operator*() [3/22]

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

Composition of two linear operators first_op and second_op given by \((\mathrm{first\_op}*\mathrm{second\_op})x \dealcoloneq \mathrm{first\_op}(\mathrm{second\_op}(x))\)

Definition at line 577 of file linear_operator.h.

◆ transpose_operator() [1/2]

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

Return the transpose linear operations of op.

Author
Matthias Maier, 2015

Definition at line 651 of file linear_operator.h.

◆ inverse_operator() [1/8]

template<typename Payload , typename Solver , typename Preconditioner , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver,
const Preconditioner &  preconditioner 
)

Return 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.

Author
Luca Heltai, Matthias Maier, Jean-Paul Pelteret, 2015

Definition at line 693 of file linear_operator.h.

◆ inverse_operator() [2/8]

template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver,
const LinearOperator< Range, Domain, Payload > &  preconditioner 
)

Variant of above function that takes a LinearOperator preconditioner as preconditioner argument.

Definition at line 750 of file linear_operator.h.

◆ inverse_operator() [3/8]

template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver 
)

Variant of above function without a preconditioner argument. In this case the identity_operator() of the op argument is used as a preconditioner. This is equivalent to using PreconditionIdentity.

Definition at line 808 of file linear_operator.h.

◆ inverse_operator() [4/8]

template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver,
const PreconditionIdentity  
)

Special overload of above function that takes a PreconditionIdentity argument.

Definition at line 828 of file linear_operator.h.

◆ identity_operator() [1/4]

template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > identity_operator ( const std::function< void(Range &, bool)> &  reinit_vector)

Return 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 858 of file linear_operator.h.

◆ identity_operator() [2/4]

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

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

The function takes a LinearOperator op and uses its range initializer to create an identity operator. In contrast to the function above, this function also ensures that the underlying Payload matches that of the input op.

Definition at line 891 of file linear_operator.h.

◆ null_operator() [1/2]

template<typename Range = Vector<double>, typename Domain = Range, typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Domain, Payload > null_operator ( const LinearOperator< Range, Domain, Payload > &  op)

Return 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 911 of file linear_operator.h.

◆ mean_value_filter() [1/4]

template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > mean_value_filter ( const std::function< void(Range &, bool)> &  reinit_vector)

Return a LinearOperator that acts as a mean value filter. The vmult() functions of this matrix subtract the mean values of the vector.

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 948 of file linear_operator.h.

◆ mean_value_filter() [2/4]

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

Return a LinearOperator that acts as a mean value filter. The vmult() functions of this matrix subtract the mean values of the vector.

The function takes a LinearOperator op and uses its range initializer to create an mean value filter operator. The function also ensures that the underlying Payload matches that of the input op.

Definition at line 990 of file linear_operator.h.

◆ linear_operator() [1/6]

template<typename Range = Vector<double>, typename Domain = Range, typename Payload = internal::LinearOperatorImplementation::EmptyPayload, typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator ( const Matrix &  matrix)

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).

Author
Matthias Maier, 2015

Definition at line 1373 of file linear_operator.h.

◆ linear_operator() [2/6]

template<typename Range = Vector<double>, typename Domain = Range, typename Payload = internal::LinearOperatorImplementation::EmptyPayload, typename OperatorExemplar , typename Matrix >
LinearOperator< Range, Domain, Payload > linear_operator ( const OperatorExemplar &  operator_exemplar,
const Matrix &  matrix 
)

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).

Author
Matthias Maier, 2015

Definition at line 1402 of file linear_operator.h.

◆ linear_operator() [3/6]

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

Variant of above function that takes a LinearOperator operator_exemplar as an additional reference. The reinit_domain_vector and reinit_range_vector function are copied from the operator_exemplar object.

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).

Author
Matthias Maier, 2017

Definition at line 1457 of file linear_operator.h.

◆ operator+() [2/10]

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

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

Definition at line 296 of file packaged_operation.h.

◆ operator-() [2/10]

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

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

Definition at line 329 of file packaged_operation.h.

◆ operator*() [4/22]

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

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

Definition at line 364 of file packaged_operation.h.

◆ operator*() [5/22]

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

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

Definition at line 404 of file packaged_operation.h.

◆ operator+() [3/10]

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

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

Definition at line 420 of file packaged_operation.h.

◆ operator+() [4/10]

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

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

Definition at line 435 of file packaged_operation.h.

◆ operator-() [3/10]

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

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

Definition at line 450 of file packaged_operation.h.

◆ operator-() [4/10]

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

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 467 of file packaged_operation.h.

◆ operator+() [5/10]

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

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 530 of file packaged_operation.h.

◆ operator-() [5/10]

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

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 575 of file packaged_operation.h.

◆ operator*() [6/22]

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

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 618 of file packaged_operation.h.

◆ operator*() [7/22]

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

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 643 of file packaged_operation.h.

◆ operator*() [8/22]

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

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 668 of file packaged_operation.h.

◆ operator*() [9/22]

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

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 703 of file packaged_operation.h.

◆ operator*() [10/22]

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

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

Definition at line 730 of file packaged_operation.h.

◆ operator*() [11/22]

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

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

Definition at line 774 of file packaged_operation.h.

◆ schur_complement() [1/2]

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

Return 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.

// 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);
// D and S operate on same space
const auto S_prec = PreconditionSelector<...>(D);
const auto S_inv = inverse_operator<...>(S,...,prec_S);
// Solve reduced block system
// PackagedOperation that represents the condensed form of g
auto rhs = condense_schur_rhs (A_inv,C,f,g);
// Solve for y
y = S_inv * rhs;
// Compute x using resolved solution y
x = postprocess_schur_solution (A_inv,B,y,f);

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

See also
Block (linear algebra)
Author
Jean-Paul Pelteret, Matthias Maier, Martin Kronbichler, 2015, 2017

Definition at line 249 of file schur_complement.h.

◆ condense_schur_rhs() [1/2]

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

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 \]

See also
Block (linear algebra)
Author
Jean-Paul Pelteret, Matthias Maier, 2015, 2017

Definition at line 300 of file schur_complement.h.

◆ postprocess_schur_solution() [1/2]

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

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) \]

See also
Block (linear algebra)
Author
Jean-Paul Pelteret, Matthias Maier, 2015, 2017

Definition at line 339 of file schur_complement.h.

◆ block_operator() [3/4]

template<typename Range , typename Domain , typename BlockPayload , typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > 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 622 of file block_linear_operator.h.

◆ block_operator() [4/4]

template<std::size_t m, std::size_t n, typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > block_operator ( const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::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 block vectors, 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 689 of file block_linear_operator.h.

◆ block_diagonal_operator() [4/6]

template<typename Range , typename Domain , typename BlockPayload , typename BlockMatrixType >
BlockLinearOperator< Range, Domain, BlockPayload > 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 742 of file block_linear_operator.h.

◆ block_diagonal_operator() [5/6]

template<std::size_t m, typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator ( const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::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 798 of file block_linear_operator.h.

◆ block_diagonal_operator() [6/6]

template<std::size_t m, typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator ( const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::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 847 of file block_linear_operator.h.

◆ block_forward_substitution() [2/2]

template<typename Range , typename Domain , typename BlockPayload >
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution ( const BlockLinearOperator< Range, Domain, BlockPayload > &  block_operator,
const BlockLinearOperator< Domain, Range, BlockPayload > &  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 909 of file block_linear_operator.h.

◆ block_back_substitution() [2/2]

template<typename Range , typename Domain , typename BlockPayload >
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution ( const BlockLinearOperator< Range, Domain, BlockPayload > &  block_operator,
const BlockLinearOperator< Domain, Range, BlockPayload > &  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 1024 of file block_linear_operator.h.

◆ operator+() [6/10]

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

Addition of two linear operators first_op and second_op given by \((\mathrm{first\_op}+\mathrm{second\_op})x \dealcoloneq \mathrm{first\_op}(x) + \mathrm{second\_op}(x)\)

Definition at line 380 of file linear_operator.h.

◆ operator-() [6/10]

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

Subtraction of two linear operators first_op and second_op given by \((\mathrm{first\_op}-\mathrm{second\_op})x \dealcoloneq \mathrm{first\_op}(x) - \mathrm{second\_op}(x)\)

Definition at line 439 of file linear_operator.h.

◆ operator*() [12/22]

template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* ( typename Range::value_type  number,
const LinearOperator< Range, Domain, Payload > &  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 477 of file linear_operator.h.

◆ operator*() [13/22]

template<typename Range , typename Domain , typename Payload >
LinearOperator< Range, Domain, Payload > operator* ( const LinearOperator< Range, Domain, Payload > &  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 544 of file linear_operator.h.

◆ operator*() [14/22]

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

Composition of two linear operators first_op and second_op given by \((\mathrm{first\_op}*\mathrm{second\_op})x \dealcoloneq \mathrm{first\_op}(\mathrm{second\_op}(x))\)

Definition at line 577 of file linear_operator.h.

◆ transpose_operator() [2/2]

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

Return the transpose linear operations of op.

Author
Matthias Maier, 2015

Definition at line 651 of file linear_operator.h.

◆ inverse_operator() [5/8]

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

Return 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.

Author
Luca Heltai, Matthias Maier, Jean-Paul Pelteret, 2015

Definition at line 693 of file linear_operator.h.

◆ inverse_operator() [6/8]

template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver,
const LinearOperator< Range, Domain, Payload > &  preconditioner 
)
related

Variant of above function that takes a LinearOperator preconditioner as preconditioner argument.

Definition at line 750 of file linear_operator.h.

◆ inverse_operator() [7/8]

template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver 
)
related

Variant of above function without a preconditioner argument. In this case the identity_operator() of the op argument is used as a preconditioner. This is equivalent to using PreconditionIdentity.

Definition at line 808 of file linear_operator.h.

◆ inverse_operator() [8/8]

template<typename Payload , typename Solver , typename Range = typename Solver::vector_type, typename Domain = Range>
LinearOperator< Domain, Range, Payload > inverse_operator ( const LinearOperator< Range, Domain, Payload > &  op,
Solver solver,
const PreconditionIdentity  
)
related

Special overload of above function that takes a PreconditionIdentity argument.

Definition at line 828 of file linear_operator.h.

◆ identity_operator() [3/4]

template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > identity_operator ( const std::function< void(Range &, bool)> &  reinit_vector)
related

Return 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 858 of file linear_operator.h.

◆ identity_operator() [4/4]

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

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

The function takes a LinearOperator op and uses its range initializer to create an identity operator. In contrast to the function above, this function also ensures that the underlying Payload matches that of the input op.

Definition at line 891 of file linear_operator.h.

◆ null_operator() [2/2]

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

Return 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 911 of file linear_operator.h.

◆ mean_value_filter() [3/4]

template<typename Range , typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator< Range, Range, Payload > mean_value_filter ( const std::function< void(Range &, bool)> &  reinit_vector)
related

Return a LinearOperator that acts as a mean value filter. The vmult() functions of this matrix subtract the mean values of the vector.

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 948 of file linear_operator.h.

◆ mean_value_filter() [4/4]

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

Return a LinearOperator that acts as a mean value filter. The vmult() functions of this matrix subtract the mean values of the vector.

The function takes a LinearOperator op and uses its range initializer to create an mean value filter operator. The function also ensures that the underlying Payload matches that of the input op.

Definition at line 990 of file linear_operator.h.

◆ linear_operator() [4/6]

template<typename Range , typename Domain , typename Payload , typename Matrix >
LinearOperator< Range, Domain, Payload > 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).

Author
Matthias Maier, 2015

Definition at line 1373 of file linear_operator.h.

◆ linear_operator() [5/6]

template<typename Range , typename Domain , typename Payload , typename OperatorExemplar , typename Matrix >
LinearOperator< Range, Domain, Payload > 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).

Author
Matthias Maier, 2015

Definition at line 1402 of file linear_operator.h.

◆ linear_operator() [6/6]

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

Variant of above function that takes a LinearOperator operator_exemplar as an additional reference. The reinit_domain_vector and reinit_range_vector function are copied from the operator_exemplar object.

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).

Author
Matthias Maier, 2017

Definition at line 1457 of file linear_operator.h.

◆ operator+() [7/10]

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 296 of file packaged_operation.h.

◆ operator-() [7/10]

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 329 of file packaged_operation.h.

◆ operator*() [15/22]

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 364 of file packaged_operation.h.

◆ operator*() [16/22]

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 404 of file packaged_operation.h.

◆ operator+() [8/10]

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 420 of file packaged_operation.h.

◆ operator+() [9/10]

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 435 of file packaged_operation.h.

◆ operator-() [8/10]

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 450 of file packaged_operation.h.

◆ operator-() [9/10]

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 467 of file packaged_operation.h.

◆ operator+() [10/10]

template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::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 530 of file packaged_operation.h.

◆ operator-() [10/10]

template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::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 575 of file packaged_operation.h.

◆ operator*() [17/22]

template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::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 618 of file packaged_operation.h.

◆ operator*() [18/22]

template<typename Range , typename = typename std::enable_if< internal::PackagedOperationImplementation::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 643 of file packaged_operation.h.

◆ operator*() [19/22]

template<typename Range , typename Domain , typename Payload >
PackagedOperation< Range > operator* ( const LinearOperator< Range, Domain, Payload > &  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 668 of file packaged_operation.h.

◆ operator*() [20/22]

template<typename Range , typename Domain , typename Payload >
PackagedOperation< Domain > operator* ( const Range &  u,
const LinearOperator< Range, Domain, Payload > &  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 703 of file packaged_operation.h.

◆ operator*() [21/22]

template<typename Range , typename Domain , typename Payload >
PackagedOperation< Range > operator* ( const LinearOperator< Range, Domain, Payload > &  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 730 of file packaged_operation.h.

◆ operator*() [22/22]

template<typename Range , typename Domain , typename Payload >
PackagedOperation< Domain > operator* ( const PackagedOperation< Range > &  comp,
const LinearOperator< Range, Domain, Payload > &  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 774 of file packaged_operation.h.

◆ schur_complement() [2/2]

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

Return 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.

// 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);
// D and S operate on same space
const auto S_prec = PreconditionSelector<...>(D);
const auto S_inv = inverse_operator<...>(S,...,prec_S);
// Solve reduced block system
// PackagedOperation that represents the condensed form of g
auto rhs = condense_schur_rhs (A_inv,C,f,g);
// Solve for y
y = S_inv * rhs;
// Compute x using resolved solution y
x = postprocess_schur_solution (A_inv,B,y,f);

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

See also
Block (linear algebra)
Author
Jean-Paul Pelteret, Matthias Maier, Martin Kronbichler, 2015, 2017

Definition at line 249 of file schur_complement.h.

◆ condense_schur_rhs() [2/2]

template<typename Range_1 , typename Domain_1 , typename Range_2 , typename Payload >
PackagedOperation< Range_2 > condense_schur_rhs ( const LinearOperator< Range_1, Domain_1, Payload > &  A_inv,
const LinearOperator< Range_2, Domain_1, Payload > &  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 \]

See also
Block (linear algebra)
Author
Jean-Paul Pelteret, Matthias Maier, 2015, 2017

Definition at line 300 of file schur_complement.h.

◆ postprocess_schur_solution() [2/2]

template<typename Range_1 , typename Domain_1 , typename Domain_2 , typename Payload >
PackagedOperation< Domain_1 > postprocess_schur_solution ( const LinearOperator< Range_1, Domain_1, Payload > &  A_inv,
const LinearOperator< Range_1, Domain_2, Payload > &  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) \]

See also
Block (linear algebra)
Author
Jean-Paul Pelteret, Matthias Maier, 2015, 2017

Definition at line 339 of file schur_complement.h.

Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
LinearOperator::operator*=
LinearOperator< Range, Domain, Payload > & operator*=(const LinearOperator< Domain, Domain, Payload > &second_op)
Definition: linear_operator.h:336
schur_complement.h
SparseMatrix< double >
PreconditionSelector
Definition: precondition_selector.h:103
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
LinearOperator::vmult_add
std::function< void(Range &v, const Domain &u)> vmult_add
Definition: linear_operator.h:271
LinearOperator::vmult
std::function< void(Range &v, const Domain &u)> vmult
Definition: linear_operator.h:265
LinearOperator::inverse_operator
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
Definition: linear_operator.h:693
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
LinearOperator::schur_complement
LinearOperator< Range_2, Domain_2, Payload > schur_complement(const LinearOperator< Domain_1, Range_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const LinearOperator< Range_2, Domain_1, Payload > &C, const LinearOperator< Range_2, Domain_2, Payload > &D)
Definition: schur_complement.h:249
LinearOperator::linear_operator
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
Definition: linear_operator.h:1373
LinearOperator::Tvmult_add
std::function< void(Domain &v, const Range &u)> Tvmult_add
Definition: linear_operator.h:283
PackagedOperation::condense_schur_rhs
PackagedOperation< Range_2 > condense_schur_rhs(const LinearOperator< Range_1, Domain_1, Payload > &A_inv, const LinearOperator< Range_2, Domain_1, Payload > &C, const Range_1 &f, const Range_2 &g)
Definition: schur_complement.h:300
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
LinearOperator::Tvmult
std::function< void(Domain &v, const Range &u)> Tvmult
Definition: linear_operator.h:277
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
PackagedOperation::postprocess_schur_solution
PackagedOperation< Domain_1 > postprocess_schur_solution(const LinearOperator< Range_1, Domain_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const Domain_2 &y, const Range_1 &f)
Definition: schur_complement.h:339
Vector< double >
LinearOperator
Definition: linear_operator.h:169