Reference documentation for deal.II version Git d228448 20171017 21:17:08 +0200

Creation of a BlockLinearOperator  
template<typename Range , typename Domain , typename BlockPayload , typename BlockMatrixType >  
BlockLinearOperator< Range, Domain, BlockPayload >  block_operator (const BlockMatrixType &block_matrix) 
template<size_t m, 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<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<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 , 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) 
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) 
Creation of a LinearOperator  
template<typename Range , typename Payload = internal::LinearOperator::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 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) 
Creation of a PackagedOperation object  
template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>  
PackagedOperation< Range >  operator+ (const Range &u, const Range &v) 
template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>  
PackagedOperation< Range >  operator (const Range &u, const Range &v) 
template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>  
PackagedOperation< Range >  operator* (const Range &u, typename Range::value_type number) 
template<typename Range , typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>  
PackagedOperation< Range >  operator* (typename Range::value_type number, const Range &u) 
template<typename Range , typename Domain , 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) 
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) 
If deal.II is configured with C++11 support (i.e., DEAL_II_WITH_CXX11=ON
or DEAL_II_WITH_CXX14=ON
during configuration) a versatile mechanism for storing the concept of a linear operator is available. (For questions about C++11, see deal.II and the C++11 standard .)
This is done with a LinearOperator class that, like the MatrixType concept, defines a minimal interface for applying a linear operation on a vector.
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 matrixvector 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:
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:
block_operator() and block_diagonal_operator() provide further encapsulation of individual linear operators into blocked linear operator variants.
res = op_a*x
a PackagedOperation class instance is generated behindthecurtains. 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 templatesrelated compiler errors, all headers relevant to LinearOperator are grouped in <deal.ii/lac/linear_operator_tools.h>.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:
Converting the PackagedOperation a + b  c + d
to a vector results in code equivalent to the following code
that avoids any intermediate storage. As a second example (involving a LinearOperator object) consider the computation of a residual \(bAx\):
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).
residual_expr
object 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
):

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

related 
A variant of above function that encapsulates a given collection ops
of LinearOperators into a block structure. Here, it is assumed that Range and Domain are blockvectors, i.e., derived from BlockVectorBase. The individual linear operators in ops
must act on the underlying vector type of the block vectors, i.e., on Domain::BlockType yielding a result in Range::BlockType.
The list ops
is best passed as an initializer list. Consider for example a linear operator block (acting on Vector<double>)
The corresponding block_operator invocation takes the form
Definition at line 499 of file block_linear_operator.h.

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

related 
A variant of above function that builds up a block diagonal linear operator from an array ops
of diagonal elements (offdiagonal 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
Definition at line 608 of file block_linear_operator.h.

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

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:
First of all, x0 = A00^1 y0
. Then, we can use x0 to recover x1:
and therefore:
block_operator
is used as well as the diagonal of diagonal_inverse
. Definition at line 710 of file block_linear_operator.h.

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:
First of all, xn = Ann^1 yn
. Then, we can use xn to recover x(n1):
and therefore:
block_operator
is used as well as the diagonal of diagonal_inverse
. Definition at line 820 of file block_linear_operator.h.

related 
Addition of two linear operators first_op
and second_op
given by \((\text{first\_op}+\text{second\_op})x := \text{first\_op}(x) + \text{second\_op}(x)\)
Definition at line 359 of file linear_operator.h.

related 
Subtraction of two linear operators first_op
and second_op
given by \((\text{first\_op}\text{second\_op})x := \text{first\_op}(x)  \text{second\_op}(x)\)
Definition at line 422 of file linear_operator.h.

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:
Definition at line 460 of file linear_operator.h.

related 
Scalar multiplication of a ScalarOperator object from the right.
The Domain
and Range
types must implement the following operator*=
member functions for rescaling:
Definition at line 530 of file linear_operator.h.

related 
Composition of two linear operators first_op
and second_op
given by \((\text{first\_op}*\text{second\_op})x := \text{first\_op}(\text{second\_op}(x))\)
Definition at line 559 of file linear_operator.h.

related 
Return the transpose linear operations of op
.
Definition at line 636 of file linear_operator.h.

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
.
Definition at line 678 of file linear_operator.h.

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

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 identiy operator. In contrast to the function above, this function also ensures that the underlying Payload matches that of the input op
.
Definition at line 793 of file linear_operator.h.

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

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:
The following (optional) interface is used if available:
If the Matrix does not provide vmult_add
and Tvmult_add
, they are implemented in terms of vmult
and Tvmult
(requiring intermediate storage).
Definition at line 1228 of file linear_operator.h.

related 
Variant of above function that takes an operator object operator_exemplar
as an additional reference. This object is used to populate the reinit_domain_vector and reinit_range_vector function objects. The reference matrix
is used to construct vmult, Tvmult, etc.
This variant can, for example, be used to encapsulate preconditioners (that typically do not expose any information about the underlying matrix).
Definition at line 1256 of file linear_operator.h.

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).
Definition at line 1309 of file linear_operator.h.

related 
Addition of two PackagedOperation objects first_comp
and second_comp
given by vector space addition of the corresponding results.
Definition at line 287 of file packaged_operation.h.

related 
Subtraction of two PackagedOperation objects first_comp
and second_comp
given by vector space addition of the corresponding results.
Definition at line 322 of file packaged_operation.h.

related 
Scalar multiplication of a PackagedOperation objects comp
with a scalar number
given by a scaling PackagedOperation result with number
.
Definition at line 360 of file packaged_operation.h.

related 
Scalar multiplication of a PackagedOperation objects comp
with a scalar number
given by a scaling PackagedOperation result with number
.
Definition at line 408 of file packaged_operation.h.

related 
Add a constant offset
(of the Range
space) to the result of a PackagedOperation.
Definition at line 423 of file packaged_operation.h.

related 
Add a constant offset
(of the Range
space) to the result of a PackagedOperation.
Definition at line 438 of file packaged_operation.h.

related 
Subtract a constant offset
(of the Range
space) from the result of a PackagedOperation.
Definition at line 453 of file packaged_operation.h.

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

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

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

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

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

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

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

related 
Composition of a PackagedOperation object with a LinearOperator. The object stores the computation \(\text{op} \,comp\) (in matrix notation).
Definition at line 739 of file packaged_operation.h.

related 
Composition of a PackagedOperation object with a LinearOperator. The object stores the computation \(\text{op}^T \,comp\) (in matrix notation).
Definition at line 785 of file packaged_operation.h.

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 righthand 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 wellconditioned.
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:
A_inv
(using inverse_operator()).\[ g' = g  C \: A^{1} \: f \]
\[ y = S^{1} g' \]
\[ x = A^{1} (f  By) \]
An illustration of typical usage of this operator for a fully coupled system is given below.
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 step22, 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} \).
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 testsuite, such as tests/lac/schur_complement_01.cc
. The solution of a multi component problem (namely step22) using the schur_complement can be found in tests/lac/schur_complement_03.cc
.
Definition at line 229 of file schur_complement.h.

related 
For the system of equations
\begin{eqnarray*} Ax + By &=& f \\ Cx + Dy &=& g \quad , \end{eqnarray*}
this operation performs the preprocessing (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 \]
Definition at line 278 of file schur_complement.h.

related 
For the system of equations
\begin{eqnarray*} Ax + By &=& f \\ Cx + Dy &=& g \quad , \end{eqnarray*}
this operation performs the postprocessing 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) \]
Definition at line 315 of file schur_complement.h.