deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40:00+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Public Attributes | List of all members
BlockLinearOperator< Range, Domain, BlockPayload > Class Template Reference

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

Inheritance diagram for BlockLinearOperator< Range, Domain, BlockPayload >:
Inheritance graph
[legend]

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::TrilinosPayloadlinear_operator (const TrilinosWrappers::SparseMatrix &operator_exemplar, const Matrix &matrix)
 
LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayloadlinear_operator (const TrilinosWrappers::SparseMatrix &matrix)
 
LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayloadlinear_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)
 

Detailed Description

template<typename Range, typename Domain, typename BlockPayload>
class BlockLinearOperator< Range, Domain, BlockPayload >

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:

std::function<unsigned int()> n_block_rows;
std::function<unsigned int()> n_block_cols;
std::function<BlockType(unsigned int, unsigned int)> block;
LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > BlockType
std::function< BlockType(unsigned int, unsigned int)> block
std::function< unsigned int()> n_block_cols
std::function< unsigned int()> n_block_rows

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:

const auto block_op_a = block_operator(A);
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &block_matrix)

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:

FullMatrix<double> top_left(2, 2);
top_left(0, 0) = 2.0;
top_left(0, 1) = -1.0;
top_left(1, 0) = -1.0;
top_left(1, 1) = 2.0;
SparseMatrixEZ<double> bottom_right(4, 4, 4);
for (std::size_t row_n = 0; row_n < 4; ++row_n)
{
bottom_right.add(row_n, row_n, 1.0);
if (row_n < 3)
bottom_right.add(row_n, row_n + 1, -1.0);
}
auto top_left_op = linear_operator(top_left);
auto bottom_right_op = linear_operator(bottom_right);
std::array<decltype(top_left_op), 2> operators {{top_left_op,
bottom_right_op}};
auto block_op = block_diagonal_operator (operators);
std::vector<BlockVector<double>::size_type> block_sizes {2, 4};
BlockVector<double> src(block_sizes);
src = 2.0;
BlockVector<double> dst(block_sizes);
block_op.vmult(dst, src); // now equal to 2, 2, 0, 0, 0, 2
LinearOperator< Range, Domain, BlockPayload::BlockType > linear_operator(const Matrix &matrix)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const BlockMatrixType &block_matrix)

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.

Note
This class makes heavy use of 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.

Member Typedef Documentation

◆ BlockType

template<typename Range , typename Domain , typename BlockPayload >
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.

Constructor & Destructor Documentation

◆ BlockLinearOperator() [1/5]

template<typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload >::BlockLinearOperator ( const BlockPayload &  payload)
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.

◆ BlockLinearOperator() [2/5]

template<typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload >::BlockLinearOperator ( const BlockLinearOperator< Range, Domain, BlockPayload > &  )
default

Default copy constructor.

◆ BlockLinearOperator() [3/5]

template<typename Range , typename Domain , typename BlockPayload >
template<typename Op >
BlockLinearOperator< Range, Domain, BlockPayload >::BlockLinearOperator ( const Op &  op)
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.

◆ BlockLinearOperator() [4/5]

template<typename Range , typename Domain , typename BlockPayload >
template<std::size_t m, std::size_t n>
BlockLinearOperator< Range, Domain, BlockPayload >::BlockLinearOperator ( const std::array< std::array< BlockType, n >, m > &  ops)
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.

◆ BlockLinearOperator() [5/5]

template<typename Range , typename Domain , typename BlockPayload >
template<std::size_t m>
BlockLinearOperator< Range, Domain, BlockPayload >::BlockLinearOperator ( const std::array< BlockType, m > &  ops)
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.

Member Function Documentation

◆ operator=() [1/4]

template<typename Range , typename Domain , typename BlockPayload >
BlockLinearOperator< Range, Domain, BlockPayload > & BlockLinearOperator< Range, Domain, BlockPayload >::operator= ( const BlockLinearOperator< Range, Domain, BlockPayload > &  )
default

Default copy assignment operator.

◆ operator=() [2/4]

template<typename Range , typename Domain , typename BlockPayload >
template<typename Op >
BlockLinearOperator< Range, Domain, BlockPayload > & BlockLinearOperator< Range, Domain, BlockPayload >::operator= ( const Op &  op)
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.

◆ operator=() [3/4]

template<typename Range , typename Domain , typename BlockPayload >
template<std::size_t m, std::size_t n>
BlockLinearOperator< Range, Domain, BlockPayload > & BlockLinearOperator< Range, Domain, BlockPayload >::operator= ( const std::array< std::array< BlockType, n >, m > &  ops)
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.

◆ operator=() [4/4]

template<typename Range , typename Domain , typename BlockPayload >
template<std::size_t m>
BlockLinearOperator< Range, Domain, BlockPayload > & BlockLinearOperator< Range, Domain, BlockPayload >::operator= ( const std::array< BlockType, m > &  ops)
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.

◆ operator+=()

LinearOperator< Range, Domain, BlockPayload::BlockType > & LinearOperator< Range, Domain, BlockPayload::BlockType >::operator+= ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  second_op)
inlineinherited

Addition with a LinearOperator second_op with the same Domain and Range.

Definition at line 343 of file linear_operator.h.

◆ operator-=()

LinearOperator< Range, Domain, BlockPayload::BlockType > & LinearOperator< Range, Domain, BlockPayload::BlockType >::operator-= ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  second_op)
inlineinherited

Subtraction with a LinearOperator second_op with the same Domain and Range.

Definition at line 354 of file linear_operator.h.

◆ operator*=() [1/2]

LinearOperator< Range, Domain, BlockPayload::BlockType > & LinearOperator< Range, Domain, BlockPayload::BlockType >::operator*= ( const LinearOperator< Domain, Domain, BlockPayload::BlockType > &  second_op)
inlineinherited

Composition of the LinearOperator with an endomorphism second_op of the Domain space.

Definition at line 365 of file linear_operator.h.

◆ operator*=() [2/2]

LinearOperator< Range, Domain, BlockPayload::BlockType > LinearOperator< Range, Domain, BlockPayload::BlockType >::operator*= ( typename Domain::value_type  number)
inlineinherited

Scalar multiplication of the LinearOperator with number from the right.

Definition at line 376 of file linear_operator.h.

Friends And Related Symbol Documentation

◆ distribute_constraints_linear_operator()

LinearOperator< Range, Domain, BlockPayload::BlockType > distribute_constraints_linear_operator ( const AffineConstraints< typename Range::value_type > &  constraints,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  exemplar 
)
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.

Note
Currently, this function may not work correctly for distributed data structures.

Definition at line 64 of file constrained_linear_operator.h.

◆ project_to_constrained_linear_operator()

LinearOperator< Range, Domain, BlockPayload::BlockType > project_to_constrained_linear_operator ( const AffineConstraints< typename Range::value_type > &  constraints,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  exemplar 
)
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.

◆ constrained_linear_operator()

LinearOperator< Range, Domain, BlockPayload::BlockType > constrained_linear_operator ( const AffineConstraints< typename Range::value_type > &  constraints,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  linop 
)
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:

Ct * linop * C + Id_c;

with

Id_c = project_to_constrained_linear_operator(constraints, linop);
LinearOperator< Domain, Range, BlockPayload::BlockType > transpose_operator(const LinearOperator< Range, Domain, BlockPayload::BlockType > &op)
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 > distribute_constraints_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, BlockPayload::BlockType > &exemplar)

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.

Note
Currently, this function may not work correctly for distributed data structures.

Definition at line 245 of file constrained_linear_operator.h.

◆ constrained_right_hand_side()

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

Ct * (right_hand_side - linop * k)

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.

Note
Currently, this function may not work correctly for distributed data structures.

Definition at line 291 of file constrained_linear_operator.h.

◆ operator+()

LinearOperator< Range, Domain, BlockPayload::BlockType > operator+ ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  first_op,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 409 of file linear_operator.h.

◆ operator-()

LinearOperator< Range, Domain, BlockPayload::BlockType > operator- ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  first_op,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 468 of file linear_operator.h.

◆ operator*() [1/3]

LinearOperator< Range, Domain, BlockPayload::BlockType > operator* ( typename Range::value_type  number,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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);
LinearOperator< Range, Domain, BlockPayload::BlockType > & operator*=(const LinearOperator< Domain, Domain, BlockPayload::BlockType > &second_op)

Definition at line 506 of file linear_operator.h.

◆ operator*() [2/3]

LinearOperator< Range, Domain, BlockPayload::BlockType > operator* ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 573 of file linear_operator.h.

◆ operator*() [3/3]

LinearOperator< Range, Domain, BlockPayload::BlockType > operator* ( const LinearOperator< Range, Intermediate, BlockPayload::BlockType > &  first_op,
const LinearOperator< Intermediate, Domain, BlockPayload::BlockType > &  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 606 of file linear_operator.h.

◆ transpose_operator()

LinearOperator< Domain, Range, BlockPayload::BlockType > transpose_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  op)
related

Return the transpose linear operations of op.

Definition at line 679 of file linear_operator.h.

◆ inverse_operator() [1/4]

LinearOperator< Domain, Range, BlockPayload::BlockType > inverse_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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.

Definition at line 720 of file linear_operator.h.

◆ inverse_operator() [2/4]

LinearOperator< Domain, Range, BlockPayload::BlockType > inverse_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  op,
Solver &  solver,
const LinearOperator< Range, Domain, BlockPayload::BlockType > &  preconditioner 
)
related

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

Definition at line 777 of file linear_operator.h.

◆ inverse_operator() [3/4]

LinearOperator< Domain, Range, BlockPayload::BlockType > inverse_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 835 of file linear_operator.h.

◆ inverse_operator() [4/4]

LinearOperator< Domain, Range, BlockPayload::BlockType > inverse_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  op,
Solver &  solver,
const PreconditionIdentity  
)
related

Special overload of above function that takes a PreconditionIdentity argument.

Definition at line 855 of file linear_operator.h.

◆ identity_operator() [1/2]

LinearOperator< Range, Range, BlockPayload::BlockType > 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 885 of file linear_operator.h.

◆ identity_operator() [2/2]

LinearOperator< Range, Domain, BlockPayload::BlockType > identity_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 918 of file linear_operator.h.

◆ null_operator()

LinearOperator< Range, Domain, BlockPayload::BlockType > null_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 938 of file linear_operator.h.

◆ mean_value_filter() [1/2]

LinearOperator< Range, Range, BlockPayload::BlockType > 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 975 of file linear_operator.h.

◆ mean_value_filter() [2/2]

LinearOperator< Range, Domain, BlockPayload::BlockType > mean_value_filter ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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 1017 of file linear_operator.h.

◆ linear_operator() [1/6]

LinearOperator< Range, Domain, BlockPayload::BlockType > linear_operator ( const Matrix &  matrix)
related

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

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

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

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

The following (optional) interface is used if available:

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

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

Definition at line 1399 of file linear_operator.h.

◆ linear_operator() [2/6]

LinearOperator< Range, Domain, BlockPayload::BlockType > linear_operator ( const OperatorExemplar &  operator_exemplar,
const Matrix &  matrix 
)
related

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

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

Definition at line 1427 of file linear_operator.h.

◆ linear_operator() [3/6]

LinearOperator< Range, Domain, BlockPayload::BlockType > linear_operator ( const LinearOperator< Range, Domain, BlockPayload::BlockType > &  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).

Definition at line 1481 of file linear_operator.h.

◆ linear_operator() [4/6]

LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > linear_operator ( const TrilinosWrappers::SparseMatrix operator_exemplar,
const Matrix &  matrix 
)
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.

◆ linear_operator() [5/6]

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.

◆ linear_operator() [6/6]

LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > linear_operator ( const LinearOperator< Range, Domain, TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload > &  operator_exemplar,
const Matrix &  matrix 
)
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.

◆ 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 
)
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);
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)
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)
LinearOperator< Domain, Range, BlockPayload::BlockType > inverse_operator(const LinearOperator< Range, Domain, BlockPayload::BlockType > &op, Solver &solver, const Preconditioner &preconditioner)
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)

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

See also
Block (linear algebra)

Definition at line 247 of file schur_complement.h.

Member Data Documentation

◆ n_block_rows

template<typename Range , typename Domain , typename BlockPayload >
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.

◆ n_block_cols

template<typename Range , typename Domain , typename BlockPayload >
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.

◆ block

template<typename Range , typename Domain , typename BlockPayload >
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.

◆ vmult

std::function<void(Range &v, const Domain &u)> LinearOperator< Range, Domain, BlockPayload::BlockType >::vmult
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.

◆ vmult_add

std::function<void(Range &v, const Domain &u)> LinearOperator< Range, Domain, BlockPayload::BlockType >::vmult_add
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.

◆ Tvmult

std::function<void(Domain &v, const Range &u)> LinearOperator< Range, Domain, BlockPayload::BlockType >::Tvmult
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.

◆ Tvmult_add

std::function<void(Domain &v, const Range &u)> LinearOperator< Range, Domain, BlockPayload::BlockType >::Tvmult_add
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.

◆ reinit_range_vector

std::function<void(Range &v, bool omit_zeroing_entries)> LinearOperator< Range, Domain, BlockPayload::BlockType >::reinit_range_vector
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.

◆ reinit_domain_vector

std::function<void(Domain &v, bool omit_zeroing_entries)> LinearOperator< Range, Domain, BlockPayload::BlockType >::reinit_domain_vector
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.

◆ is_null_operator

bool LinearOperator< Range, Domain, BlockPayload::BlockType >::is_null_operator
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.


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