deal.II version GIT relicensing-2306-ga6fff1538f 2024-12-29 17:30:00+00:00
|
#include <deal.II/lac/block_linear_operator.h>
Public Types | |
using | BlockType = LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > |
Public Member Functions | |
BlockLinearOperator (const BlockPayload &payload) | |
BlockLinearOperator (const BlockLinearOperator< Range, Domain, BlockPayload > &)=default | |
template<typename Op > | |
BlockLinearOperator (const Op &op) | |
template<std::size_t m, std::size_t n> | |
BlockLinearOperator (const std::array< std::array< BlockType, n >, m > &ops) | |
template<std::size_t m> | |
BlockLinearOperator (const std::array< BlockType, m > &ops) | |
BlockLinearOperator< Range, Domain, BlockPayload > & | operator= (const BlockLinearOperator< Range, Domain, BlockPayload > &)=default |
template<typename Op > | |
BlockLinearOperator< Range, Domain, BlockPayload > & | operator= (const Op &op) |
template<std::size_t m, std::size_t n> | |
BlockLinearOperator< Range, Domain, BlockPayload > & | operator= (const std::array< std::array< BlockType, n >, m > &ops) |
template<std::size_t m> | |
BlockLinearOperator< Range, Domain, BlockPayload > & | operator= (const std::array< BlockType, m > &ops) |
Public Attributes | |
std::function< unsigned int()> | n_block_rows |
std::function< unsigned int()> | n_block_cols |
std::function< BlockType(unsigned int, unsigned int)> | block |
std::function< void(Range &v, const Domain &u)> | vmult |
std::function< void(Range &v, const Domain &u)> | vmult_add |
std::function< void(Domain &v, const Range &u)> | Tvmult |
std::function< void(Domain &v, const Range &u)> | Tvmult_add |
std::function< void(Range &v, bool omit_zeroing_entries)> | reinit_range_vector |
std::function< void(Domain &v, bool omit_zeroing_entries)> | reinit_domain_vector |
Related Symbols | |
(Note that these are not member symbols.) | |
Creation of a BlockLinearOperator | |
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 = 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 , 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_back_substitution (const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse) |
LinearOperator< Domain, Range, typename BlockPayload::BlockType > | block_forward_substitution (const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse) |
Creation of a BlockLinearOperator | |
template<typename Range , typename Domain = Range> | |
BlockLinearOperator< Range, Domain, TrilinosWrappers::internal::BlockLinearOperatorImplementation::TrilinosBlockPayload< TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > > | block_operator (const TrilinosWrappers::BlockSparseMatrix &block_matrix) |
template<std::size_t m, std::size_t n, typename Range , typename Domain = Range> | |
BlockLinearOperator< Range, Domain, TrilinosWrappers::internal::BlockLinearOperatorImplementation::TrilinosBlockPayload< TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > > | block_operator (const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload >, n >, m > &ops) |
template<typename Range , typename Domain = Range> | |
BlockLinearOperator< Range, Domain, TrilinosWrappers::internal::BlockLinearOperatorImplementation::TrilinosBlockPayload< TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > > | block_diagonal_operator (const TrilinosWrappers::BlockSparseMatrix &block_matrix) |
template<std::size_t m, typename Range , typename Domain = Range> | |
BlockLinearOperator< Range, Domain, TrilinosWrappers::internal::BlockLinearOperatorImplementation::TrilinosBlockPayload< TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > > | block_diagonal_operator (const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload >, m > &ops) |
Indirectly applying constraints to LinearOperator | |
LinearOperator< Range, Domain, BlockPayload::BlockType > | distribute_constraints_linear_operator (const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, BlockPayload::BlockType > &exemplar) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | project_to_constrained_linear_operator (const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, BlockPayload::BlockType > &exemplar) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | constrained_linear_operator (const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, BlockPayload::BlockType > &linop) |
PackagedOperation< Range > | constrained_right_hand_side (const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, BlockPayload::BlockType > &linop, const Range &right_hand_side) |
Vector space operations | |
LinearOperator< Range, Domain, BlockPayload::BlockType > | operator+ (const LinearOperator< Range, Domain, BlockPayload::BlockType > &first_op, const LinearOperator< Range, Domain, BlockPayload::BlockType > &second_op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | operator- (const LinearOperator< Range, Domain, BlockPayload::BlockType > &first_op, const LinearOperator< Range, Domain, BlockPayload::BlockType > &second_op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | operator* (typename Range::value_type number, const LinearOperator< Range, Domain, BlockPayload::BlockType > &op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | operator* (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op, typename Domain::value_type number) |
Composition and manipulation of a LinearOperator | |
LinearOperator< Range, Domain, BlockPayload::BlockType > | operator* (const LinearOperator< Range, Intermediate, BlockPayload::BlockType > &first_op, const LinearOperator< Intermediate, Domain, BlockPayload::BlockType > &second_op) |
LinearOperator< Domain, Range, BlockPayload::BlockType > | transpose_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op) |
LinearOperator< Domain, Range, BlockPayload::BlockType > | inverse_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op, Solver &solver, const Preconditioner &preconditioner) |
LinearOperator< Domain, Range, BlockPayload::BlockType > | inverse_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op, Solver &solver, const LinearOperator< Range, Domain, BlockPayload::BlockType > &preconditioner) |
LinearOperator< Domain, Range, BlockPayload::BlockType > | inverse_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op, Solver &solver) |
LinearOperator< Domain, Range, BlockPayload::BlockType > | inverse_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op, Solver &solver, const PreconditionIdentity &) |
Creation of a LinearOperator | |
LinearOperator< Range, Range, BlockPayload::BlockType > | identity_operator (const std::function< void(Range &, bool)> &reinit_vector) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | identity_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | null_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op) |
LinearOperator< Range, Range, BlockPayload::BlockType > | mean_value_filter (const std::function< void(Range &, bool)> &reinit_vector) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | mean_value_filter (const LinearOperator< Range, Domain, BlockPayload::BlockType > &op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | linear_operator (const Matrix &matrix) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | linear_operator (const OperatorExemplar &operator_exemplar, const Matrix &matrix) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | linear_operator (const LinearOperator< Range, Domain, BlockPayload::BlockType > &operator_exemplar, const Matrix &matrix) |
Creation of a LinearOperator | |
LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > | linear_operator (const TrilinosWrappers::SparseMatrix &operator_exemplar, const Matrix &matrix) |
LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > | linear_operator (const TrilinosWrappers::SparseMatrix &matrix) |
LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > | linear_operator (const LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > &operator_exemplar, const Matrix &matrix) |
Creation of a LinearOperator related to the Schur Complement | |
LinearOperator< Range_2, Domain_2, BlockPayload::BlockType > | schur_complement (const LinearOperator< Domain_1, Range_1, BlockPayload::BlockType > &A_inv, const LinearOperator< Range_1, Domain_2, BlockPayload::BlockType > &B, const LinearOperator< Range_2, Domain_1, BlockPayload::BlockType > &C, const LinearOperator< Range_2, Domain_2, BlockPayload::BlockType > &D) |
In-place vector space operations | |
bool | is_null_operator |
LinearOperator< Range, Domain, BlockPayload::BlockType > & | operator+= (const LinearOperator< Range, Domain, BlockPayload::BlockType > &second_op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > & | operator-= (const LinearOperator< Range, Domain, BlockPayload::BlockType > &second_op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > & | operator*= (const LinearOperator< Domain, Domain, BlockPayload::BlockType > &second_op) |
LinearOperator< Range, Domain, BlockPayload::BlockType > | operator*= (typename Domain::value_type number) |
A class to store the concept of a block linear operator.
This class increases the interface of LinearOperator (which encapsulates the Matrix
interface) by three additional functions:
that describe the underlying block structure (of an otherwise opaque) linear operator.
Objects of type BlockLinearOperator can be created similarly to LinearOperator with a wrapper function:
Alternatively, there are several helper functions available for creating instances from multiple independent matrices of possibly different types. Here is an example of a block diagonal matrix created from a FullMatrix and a SparseMatrixEZ:
A BlockLinearOperator can be sliced to a LinearOperator at any time. This removes all information about the underlying block structure (because above std::function
objects are no longer available) - the linear operator interface, however, remains intact.
std::function
objects and lambda functions. This flexibility comes with a run-time penalty. Only use this object to encapsulate object with medium to large individual block sizes, and small block structure (as a rule of thumb, matrix blocks greater than \(1000\times1000\)). Definition at line 165 of file block_linear_operator.h.
using BlockLinearOperator< Range, Domain, BlockPayload >::BlockType = LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> |
Definition at line 169 of file block_linear_operator.h.
|
inline |
Create an empty BlockLinearOperator object.
Allstd::function
member objects of this class and its base class LinearOperator are initialized with default variants that throw an exception upon invocation.
Definition at line 180 of file block_linear_operator.h.
|
default |
Default copy constructor.
|
inline |
Templated copy constructor that creates a BlockLinearOperator object from an object op
for which the conversion function block_operator
is defined.
Definition at line 221 of file block_linear_operator.h.
|
inline |
Create a BlockLinearOperator from a two-dimensional array ops
of LinearOperator. This constructor calls the corresponding block_operator() specialization.
Definition at line 232 of file block_linear_operator.h.
|
inline |
Create a block-diagonal BlockLinearOperator from a one-dimensional array ops
of LinearOperator. This constructor calls the corresponding block_operator() specialization.
Definition at line 243 of file block_linear_operator.h.
|
default |
Default copy assignment operator.
|
inline |
Templated copy assignment operator for an object op
for which the conversion function block_operator
is defined.
Definition at line 260 of file block_linear_operator.h.
|
inline |
Copy assignment from a two-dimensional array ops
of LinearOperator. This assignment operator calls the corresponding block_operator() specialization.
Definition at line 273 of file block_linear_operator.h.
|
inline |
Copy assignment from a one-dimensional array ops
of LinearOperator that creates a block-diagonal BlockLinearOperator. This assignment operator calls the corresponding block_operator() specialization.
Definition at line 286 of file block_linear_operator.h.
|
inlineinherited |
Addition with a LinearOperator second_op
with the same Domain
and Range
.
Definition at line 343 of file linear_operator.h.
|
inlineinherited |
Subtraction with a LinearOperator second_op
with the same Domain
and Range
.
Definition at line 354 of file linear_operator.h.
|
inlineinherited |
Composition of the LinearOperator with an endomorphism second_op
of the Domain
space.
Definition at line 365 of file linear_operator.h.
|
inlineinherited |
Scalar multiplication of the LinearOperator with number
from the right.
Definition at line 376 of file linear_operator.h.
|
related |
This function takes an AffineConstraints object constraints
and an operator exemplar exemplar
(this exemplar is usually a linear operator that describes the system matrix - it is only used to create domain and range vectors of appropriate sizes, its action vmult
is never used). A LinearOperator object associated with the "homogeneous
action" of the underlying AffineConstraints object is returned:
Applying the LinearOperator object on a vector u
results in a vector v
that stores the result of calling AffineConstraints::distribute() on u
- with one important difference: inhomogeneities are not applied, but always treated as 0 instead.
The LinearOperator object created by this function is primarily used internally in constrained_linear_operator() to build up a modified system of linear equations. How to solve a linear system of equations with this approach is explained in detail in the Constraints on degrees of freedom topic.
Definition at line 64 of file constrained_linear_operator.h.
|
related |
Given a AffineConstraints constraints
and an operator exemplar exemplar
, return a LinearOperator that is the projection to the subspace of constrained degrees of freedom, i.e. all entries of the result vector that correspond to unconstrained degrees of freedom are set to zero.
Definition at line 157 of file constrained_linear_operator.h.
|
related |
Given a AffineConstraints object constraints
and a LinearOperator linop
, this function creates a LinearOperator object consisting of the composition of three operations and a regularization:
with
and Id_c
is the projection to the subspace consisting of all vector entries associated with constrained degrees of freedom.
This LinearOperator object is used together with constrained_right_hand_side() to build up the following modified system of linear equations:
\[ (C^T A C + Id_c) x = C^T (b - A\,k) \]
with a given (unconstrained) system matrix \(A\), right hand side \(b\), and linear constraints \(C\) with inhomogeneities \(k\).
A detailed explanation of this approach is given in the Constraints on degrees of freedom topic.
Definition at line 245 of file constrained_linear_operator.h.
|
related |
Given a AffineConstraints object constraints
, a LinearOperator linop
and a right-hand side right_hand_side
, this function creates a PackagedOperation that stores the following computation:
with
This LinearOperator object is used together with constrained_right_hand_side() to build up the following modified system of linear equations:
\[ (C^T A C + Id_c) x = C^T (b - A\,k) \]
with a given (unconstrained) system matrix \(A\), right hand side \(b\), and linear constraints \(C\) with inhomogeneities \(k\).
A detailed explanation of this approach is given in the Constraints on degrees of freedom topic.
Definition at line 291 of file constrained_linear_operator.h.
|
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 409 of file linear_operator.h.
|
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 468 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 506 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 573 of file linear_operator.h.
|
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 606 of file linear_operator.h.
|
related |
Return the transpose linear operations of op
.
Definition at line 679 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 720 of file linear_operator.h.
|
related |
Variant of above function that takes a LinearOperator preconditioner
as preconditioner argument.
Definition at line 777 of file linear_operator.h.
|
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 835 of file linear_operator.h.
|
related |
Special overload of above function that takes a PreconditionIdentity argument.
Definition at line 855 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 885 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 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 918 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 938 of file linear_operator.h.
|
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 975 of file linear_operator.h.
|
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 1017 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 1399 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 1427 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 1481 of file linear_operator.h.
|
related |
A function that encapsulates generic matrix
objects, based on an operator_exemplar
, that act on a compatible Vector type into a LinearOperator.
This function is the equivalent of the linear_operator, but ensures full compatibility with Trilinos operations by preselecting the appropriate template parameters.
Definition at line 75 of file trilinos_linear_operator.h.
|
related |
A function that encapsulates generic matrix
objects that act on a compatible Vector type into a LinearOperator.
This function is the equivalent of the linear_operator, but ensures full compatibility with Trilinos operations by preselecting the appropriate template parameters.
Definition at line 105 of file trilinos_linear_operator.h.
|
related |
A function that encapsulates generic matrix
objects, based on an operator_exemplar
, that act on a compatible Vector type into a LinearOperator.
This function is the equivalent of the linear_operator, but ensures full compatibility with Trilinos operations by preselecting the appropriate template parameters.
Definition at line 134 of file trilinos_linear_operator.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 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:
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 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} \).
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 behavior 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 behavior 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
.
Definition at line 247 of file schur_complement.h.
std::function<unsigned int()> BlockLinearOperator< Range, Domain, BlockPayload >::n_block_rows |
Return the number of blocks in a column (i.e, the number of "block rows", or the number \(m\), if interpreted as a \(m\times n\) block system).
Definition at line 296 of file block_linear_operator.h.
std::function<unsigned int()> BlockLinearOperator< Range, Domain, BlockPayload >::n_block_cols |
Return the number of blocks in a row (i.e, the number of "block columns", or the number \(n\), if interpreted as a \(m\times n\) block system).
Definition at line 302 of file block_linear_operator.h.
std::function<BlockType(unsigned int, unsigned int)> BlockLinearOperator< Range, Domain, BlockPayload >::block |
Access the block with the given coordinates. This std::function
object returns a LinearOperator representing the \((i,j)\)-th block of the BlockLinearOperator.
Definition at line 309 of file block_linear_operator.h.
|
inherited |
Application of the LinearOperator object to a vector u of the Domain
space giving a vector v of the Range
space.
Definition at line 294 of file linear_operator.h.
|
inherited |
Application of the LinearOperator object to a vector u of the Domain
space. The result is added to the vector v.
Definition at line 300 of file linear_operator.h.
|
inherited |
Application of the transpose LinearOperator object to a vector u of the Range
space giving a vector v of the Domain
space.
Definition at line 306 of file linear_operator.h.
|
inherited |
Application of the transpose LinearOperator object to a vector u
of the Range
space.The result is added to the vector v
.
Definition at line 312 of file linear_operator.h.
|
inherited |
Initializes a vector v of the Range space to be directly usable as the destination parameter in an application of vmult. Similar to the reinit functions of the vector classes, the boolean determines whether a fast initialization is done, i.e., if it is set to false the content of the vector is set to 0.
Definition at line 321 of file linear_operator.h.
|
inherited |
Initializes a vector of the Domain space to be directly usable as the source parameter in an application of vmult. Similar to the reinit functions of the vector classes, the boolean determines whether a fast initialization is done, i.e., if it is set to false the content of the vector is set to 0.
Definition at line 331 of file linear_operator.h.
|
inherited |
This bool is used to determine whether a linear operator is a null operator. In this case the class is able to optimize some operations like multiplication or addition.
Definition at line 387 of file linear_operator.h.