// Forward declarations:
+namespace internal
+{
+ namespace BlockLinearOperator
+ {
+ template<typename PayloadBlockType = internal::LinearOperator::EmptyPayload>
+ class EmptyBlockPayload;
+ }
+}
+
template <typename Number> class BlockVector;
template <typename Range = BlockVector<double>,
- typename Domain = Range>
+ typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
class BlockLinearOperator;
template <typename Range = BlockVector<double>,
typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>,
typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(const BlockMatrixType &matrix);
template <size_t m, size_t n,
typename Range = BlockVector<double>,
- typename Domain = Range>
-BlockLinearOperator<Range, Domain>
-block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &);
+ typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, n>, m> &);
template <size_t m,
typename Range = BlockVector<double>,
- typename Domain = Range>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &);
+ typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, m> &);
template <size_t m,
typename Range = BlockVector<double>,
- typename Domain = Range>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op);
+ typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> &op);
// This is a workaround for a bug in <=gcc-4.7 that does not like partial
// template default values in combination with local lambda expressions [1]
template <typename Range = BlockVector<double>,
typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>,
typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(const BlockMatrixType &block_matrix);
template <typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Domain, Range>
-block_forward_substitution(const BlockLinearOperator<Range, Domain> &,
- const BlockLinearOperator<Domain, Range> &);
+ typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+LinearOperator<Domain, Range, typename BlockPayload::BlockType>
+block_forward_substitution(const BlockLinearOperator<Range, Domain, BlockPayload> &,
+ const BlockLinearOperator<Domain, Range, BlockPayload> &);
template <typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Domain, Range>
-block_back_substitution(const BlockLinearOperator<Range, Domain> &,
- const BlockLinearOperator<Domain, Range> &);
+ typename Domain = Range,
+ typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+LinearOperator<Domain, Range, typename BlockPayload::BlockType>
+block_back_substitution(const BlockLinearOperator<Range, Domain, BlockPayload> &,
+ const BlockLinearOperator<Domain, Range, BlockPayload> &);
// end of workaround
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-class BlockLinearOperator : public LinearOperator<Range, Domain>
+template <typename Range, typename Domain, typename BlockPayload>
+class BlockLinearOperator : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
{
public:
- typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType> BlockType;
+ typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> BlockType;
/**
* Create an empty BlockLinearOperator object.
* class LinearOperator are initialized with default variants that throw an
* exception upon invocation.
*/
- BlockLinearOperator()
- : LinearOperator<Range, Domain>()
+ BlockLinearOperator(const BlockPayload &payload)
+ : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(typename BlockPayload::BlockType(payload,payload))
{
-
n_block_rows = []() -> unsigned int
{
Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
/**
* Default copy constructor.
*/
- BlockLinearOperator(const BlockLinearOperator<Range, Domain> &) =
+ BlockLinearOperator(const BlockLinearOperator<Range, Domain, BlockPayload> &) =
default;
/**
template<typename Op>
BlockLinearOperator(const Op &op)
{
- *this = block_operator<Range, Domain, Op>(op);
+ *this = block_operator<Range, Domain, BlockPayload, Op>(op);
}
/**
template<size_t m, size_t n>
BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
{
- *this = block_operator<m, n, Range, Domain>(ops);
+ *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
}
/**
template<size_t m>
BlockLinearOperator(const std::array<BlockType, m> &ops)
{
- *this = block_diagonal_operator<m, Range, Domain>(ops);
+ *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
}
/**
* Default copy assignment operator.
*/
- BlockLinearOperator<Range, Domain> &
- operator=(const BlockLinearOperator<Range, Domain> &) = default;
+ BlockLinearOperator<Range, Domain, BlockPayload> &
+ operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
/**
* Templated copy assignment operator for an object @p op for which the
* conversion function <code>block_operator</code> is defined.
*/
template <typename Op>
- BlockLinearOperator<Range, Domain> &operator=(const Op &op)
+ BlockLinearOperator<Range, Domain, BlockPayload> &operator=(const Op &op)
{
- *this = block_operator<Range, Domain, Op>(op);
+ *this = block_operator<Range, Domain, BlockPayload, Op>(op);
return *this;
}
* specialization.
*/
template <size_t m, size_t n>
- BlockLinearOperator<Range, Domain> &
+ BlockLinearOperator<Range, Domain, BlockPayload> &
operator=(const std::array<std::array<BlockType, n>, m> &ops)
{
- *this = block_operator<m, n, Range, Domain>(ops);
+ *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
return *this;
}
* operator calls the corresponding block_operator() specialization.
*/
template <size_t m>
- BlockLinearOperator<Range, Domain> &
+ BlockLinearOperator<Range, Domain, BlockPayload> &
operator=(const std::array<BlockType, m> &ops)
{
- *this = block_diagonal_operator<m, Range, Domain>(ops);
+ *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
return *this;
}
{
// Populate the LinearOperator interfaces with the help of the
// BlockLinearOperator functions
- template <typename Range, typename Domain>
+ template <typename Range, typename Domain, typename BlockPayload>
inline void
populate_linear_operator_functions(
- dealii::BlockLinearOperator<Range, Domain> &op)
+ dealii::BlockLinearOperator<Range, Domain, BlockPayload> &op)
{
op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries)
{
op.block(j, i).Tvmult_add(v.block(i), u.block(j));
};
}
+
+
+ /**
+ * A dummy class for BlockLinearOperators that do not require any extensions
+ * to facilitate the operations of the block matrix or its subblocks.
+ *
+ * This is the Payload class typically associated with deal.II's native
+ * BlockSparseMatrix. To use Trilinos and PETSc BlockSparseMatrices it is
+ * necessary to initialize a BlockLinearOperator with their associated
+ * BlockPayload.
+ *
+ * @author Jean-Paul Pelteret, Matthias Maier, 2016
+ *
+ * @ingroup LAOperators
+ */
+ template<typename PayloadBlockType>
+ class EmptyBlockPayload
+ {
+ public:
+ /**
+ * Type of payload held by each subblock
+ */
+ typedef PayloadBlockType BlockType;
+
+ /**
+ * Default constructor
+ *
+ * Since this class does not do anything in particular and needs no special
+ * configuration, we have only one generic constructor that can be called
+ * under any conditions.
+ */
+ template <typename... Args>
+ EmptyBlockPayload (const Args &...)
+ { }
+ };
+
} /*namespace BlockLinearOperator*/
} /*namespace internal*/
*/
template <typename Range,
typename Domain,
+ typename BlockPayload,
typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(const BlockMatrixType &block_matrix)
{
- typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+ typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
- BlockLinearOperator<Range, Domain> return_op;
+ BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
return_op.n_block_rows = [&block_matrix]() -> unsigned int
{
*
* @ingroup LAOperators
*/
-template <size_t m, size_t n, typename Range, typename Domain>
-BlockLinearOperator<Range, Domain>
-block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
+template <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)
{
static_assert(m > 0 && n > 0,
"a blocked LinearOperator must consist of at least one block");
- typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+ typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
- BlockLinearOperator<Range, Domain> return_op;
+ BlockLinearOperator<Range, Domain, BlockPayload> return_op ((BlockPayload())); // TODO: Create block payload so that this can be initialised correctly
return_op.n_block_rows = []() -> unsigned int
{
*/
template <typename Range,
typename Domain,
+ typename BlockPayload,
typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(const BlockMatrixType &block_matrix)
{
- typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+ typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
- BlockLinearOperator<Range, Domain> return_op;
+ BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
return_op.n_block_rows = [&block_matrix]() -> unsigned int
{
*
* @ingroup LAOperators
*/
-template <size_t m, typename Range, typename Domain>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
+template <size_t m, typename Range, typename Domain, typename BlockPayload>
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, m> &ops)
{
static_assert(m > 0,
"a blockdiagonal LinearOperator must consist of at least one block");
- typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+ typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
std::array<std::array<BlockType, m>, m> new_ops;
*
* @ingroup LAOperators
*/
-template <size_t m, typename Range, typename Domain>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
+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)
{
static_assert(m > 0,
"a blockdiagonal LinearOperator must consist of at least "
"one block");
- typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+ typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
std::array<BlockType, m> new_ops;
new_ops.fill(op);
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-block_forward_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
- const BlockLinearOperator<Domain, Range> &diagonal_inverse)
+template <typename Range, typename Domain, 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)
{
- LinearOperator<Range, Range> return_op;
+ LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-block_back_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
- const BlockLinearOperator<Domain, Range> &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)
{
- LinearOperator<Range, Range> return_op;
+ LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
// Forward declarations:
+namespace internal
+{
+ namespace LinearOperator
+ {
+ class EmptyPayload;
+ }
+}
+
template <typename Number> class Vector;
-template <typename Range = Vector<double>, typename Domain = Range>
+template <typename Range = Vector<double>,
+ typename Domain = Range,
+ typename Payload = internal::LinearOperator::EmptyPayload>
class LinearOperator;
template <typename Range = Vector<double>,
typename Domain = Range,
+ typename Payload = internal::LinearOperator::EmptyPayload,
typename OperatorExemplar,
typename Matrix>
-LinearOperator<Range, Domain> linear_operator (const OperatorExemplar &,
- const Matrix &);
+LinearOperator<Range, Domain, Payload> linear_operator (const OperatorExemplar &,
+ const Matrix &);
template <typename Range = Vector<double>,
typename Domain = Range,
+ typename Payload = internal::LinearOperator::EmptyPayload,
typename Matrix>
-LinearOperator<Range, Domain> linear_operator (const Matrix &);
+LinearOperator<Range, Domain, Payload> linear_operator (const Matrix &);
template <typename Range = Vector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
-null_operator(const LinearOperator<Range, Domain> &);
+ typename Domain = Range,
+ typename Payload = internal::LinearOperator::EmptyPayload>
+LinearOperator<Range, Domain, Payload>
+null_operator(const LinearOperator<Range, Domain, Payload> &);
/**
* this object to encapsulate matrix object of medium to large size (as a rule
* of thumb, sparse matrices with a size $1000\times1000$, or larger).
*
+ * @note In order to use Trilinos or PETSc sparse matrices and preconditioners
+ * in conjunction with the LinearOperator class, it is necessary to extend the
+ * functionality of the LinearOperator class by means of an additional Payload.
+ *
+ * For example, in order to construct an InverseOperator a call to a solver is
+ * required. Naturally these solvers don't have an interface to the
+ * LinearOperator (which, for example, may represent a composite operation).
+ * The TrilinosWrappers::internal::LinearOperator::TrilinosPayload therefore
+ * provides an interface extension to the LinearOperator so that it can be
+ * passed to the solver and used by the solver as if it were a Trilinos
+ * operator. This implies that all of the necessary functionality of the
+ * specific Trilinos operator has been overloaded within the Payload class.
+ * This includes operator-vector multiplication and inverse operator-vector
+ * multiplication, where the operator can be either a
+ * TrilinosWrappers::SparseMatrix or a TrilinosWrappers::PreconditionBase
+ * and the vector is a native Trilinos vector.
+ *
+ * Another case where payloads provide a crucial supplement to the
+ * LinearOperator class are when composite operations are constructed (via
+ * operator overloading). In this instance, it is again necessary to provide
+ * an interface that produces the result of this composite operation that is
+ * compatible with Trilinos operator used by Trilinos solvers.
+ *
+ * @note To ensure that the correct payload is provided, wrapper functions
+ * for linear operators have been provided within the respective
+ * TrilinosWrappers (and, in the future, PetscWrappers) namespaces.
+ *
* @note This class is only available if deal.II was configured with C++11
* support, i.e., if <code>DEAL_II_WITH_CXX11</code> is enabled during cmake
* configure.
*
- * @author Luca Heltai, Matthias Maier, 2015
+ * @author Luca Heltai, Matthias Maier, 2015; Jean-Paul Pelteret, 2016
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain> class LinearOperator
+template <typename Range, typename Domain, typename Payload>
+class LinearOperator
+ : public Payload
{
public:
/**
- * Create an empty LinearOperator object. All <code>std::function</code>
- * member objects are initialized with default variants that throw an
- * exception upon invocation.
+ * Create an empty LinearOperator object.
+ * When a payload is passed to this constructor, the resulting operator is
+ * constructed with a functional payload.
+ * In either case, this constructor yields an object that can not actually
+ * be used for any linear operator operations, and will throw an exception
+ * upon invocation.
*/
- LinearOperator()
+ LinearOperator(const Payload &payload = Payload())
:
+ Payload (payload),
is_null_operator(false)
{
/**
* Default copy constructor.
*/
- LinearOperator (const LinearOperator<Range, Domain> &) = default;
+ LinearOperator (const LinearOperator<Range, Domain, Payload> &) = default;
/**
* Templated copy constructor that creates a LinearOperator object from an
* <code>linear_operator</code> is defined.
*/
template<typename Op,
- typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain>, Op>::value>::type>
+ typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain, Payload>, Op>::value>::type>
LinearOperator (const Op &op)
{
- *this = linear_operator<Range, Domain, Op>(op);
+ *this = linear_operator<Range, Domain, Payload, Op>(op);
}
/**
* Default copy assignment operator.
*/
- LinearOperator<Range, Domain> &
- operator=(const LinearOperator<Range, Domain> &) = default;
+ LinearOperator<Range, Domain, Payload> &
+ operator=(const LinearOperator<Range, Domain, Payload> &) = default;
/**
* Templated copy assignment operator for an object @p op for which the
* conversion function <code>linear_operator</code> is defined.
*/
template <typename Op,
- typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain>, Op>::value>::type>
- LinearOperator<Range, Domain> &operator=(const Op &op)
+ typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain, Payload>, Op>::value>::type>
+ LinearOperator<Range, Domain, Payload> &operator=(const Op &op)
{
- *this = linear_operator<Range, Domain, Op>(op);
+ *this = linear_operator<Range, Domain, Payload, Op>(op);
return *this;
}
* Addition with a LinearOperator @p second_op with the same @p Domain and
* @p Range.
*/
- LinearOperator<Range, Domain> &
- operator+=(const LinearOperator<Range, Domain> &second_op)
+ LinearOperator<Range, Domain, Payload> &
+ operator+=(const LinearOperator<Range, Domain, Payload> &second_op)
{
*this = *this + second_op;
return *this;
* Subtraction with a LinearOperator @p second_op with the same @p Domain
* and @p Range.
*/
- LinearOperator<Range, Domain> &
- operator-=(const LinearOperator<Range, Domain> &second_op)
+ LinearOperator<Range, Domain, Payload> &
+ operator-=(const LinearOperator<Range, Domain, Payload> &second_op)
{
*this = *this - second_op;
return *this;
* Composition of the LinearOperator with an endomorphism @p second_op of
* the @p Domain space.
*/
- LinearOperator<Range, Domain> &
- operator*=(const LinearOperator<Domain, Domain> &second_op)
+ LinearOperator<Range, Domain, Payload> &
+ operator*=(const LinearOperator<Domain, Domain, Payload> &second_op)
{
*this = *this * second_op;
return *this;
* Scalar multiplication of the LinearOperator with @p number from the
* right.
*/
- LinearOperator<Range, Domain>
+ LinearOperator<Range, Domain, Payload>
operator*=(typename Domain::value_type number)
{
*this = *this * number;
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-operator+(const LinearOperator<Range, Domain> &first_op,
- const LinearOperator<Range, Domain> &second_op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator+(const LinearOperator<Range, Domain, Payload> &first_op,
+ const LinearOperator<Range, Domain, Payload> &second_op)
{
if (first_op.is_null_operator)
{
}
else
{
- LinearOperator<Range, Domain> return_op;
+ LinearOperator<Range, Domain, Payload> return_op (
+ static_cast<const Payload &>(first_op) + static_cast<const Payload &>(second_op)
+ );
return_op.reinit_range_vector = first_op.reinit_range_vector;
return_op.reinit_domain_vector = first_op.reinit_domain_vector;
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-operator-(const LinearOperator<Range, Domain> &first_op,
- const LinearOperator<Range, Domain> &second_op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator-(const LinearOperator<Range, Domain, Payload> &first_op,
+ const LinearOperator<Range, Domain, Payload> &second_op)
{
if (first_op.is_null_operator)
{
- return -1 * second_op;
+ return -1. * second_op;
}
else if (second_op.is_null_operator)
{
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
operator*(typename Range::value_type number,
- const LinearOperator<Range, Domain> &op)
+ const LinearOperator<Range, Domain, Payload> &op)
{
static_assert(
std::is_convertible<typename Range::value_type, typename Domain::value_type>::value,
}
else
{
- LinearOperator<Range, Domain> return_op = op;
+ LinearOperator<Range, Domain, Payload> return_op = op;
// ensure to have valid computation objects by catching number and op by
// value
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-operator*(const LinearOperator<Range, Domain> &op,
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator*(const LinearOperator<Range, Domain, Payload> &op,
typename Domain::value_type number)
{
static_assert(
*
* @ingroup LAOperators
*/
-template <typename Range, typename Intermediate, typename Domain>
-LinearOperator<Range, Domain>
-operator*(const LinearOperator<Range, Intermediate> &first_op,
- const LinearOperator<Intermediate, Domain> &second_op)
+template <typename Range, typename Intermediate, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator*(const LinearOperator<Range, Intermediate, Payload> &first_op,
+ const LinearOperator<Intermediate, Domain, Payload> &second_op)
{
if (first_op.is_null_operator || second_op.is_null_operator)
{
- LinearOperator<Range, Domain> return_op;
+ LinearOperator<Range, Domain, Payload> return_op;
return_op.reinit_domain_vector = second_op.reinit_domain_vector;
return_op.reinit_range_vector = first_op.reinit_range_vector;
return null_operator(return_op);
}
else
{
- LinearOperator<Range, Domain> return_op;
+ LinearOperator<Range, Domain, Payload> return_op (
+ static_cast<const Payload &>(first_op) * static_cast<const Payload &>(second_op)
+ );
return_op.reinit_domain_vector = second_op.reinit_domain_vector;
return_op.reinit_range_vector = first_op.reinit_range_vector;
*
* Return the transpose linear operations of @p op.
*
+ * @author Matthias Maier, 2015
+ *
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-transpose_operator(const LinearOperator<Range, Domain> &op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Domain, Range, Payload>
+transpose_operator(const LinearOperator<Range, Domain, Payload> &op)
{
- LinearOperator<Domain, Range> return_op;
+ LinearOperator<Domain, Range, Payload> return_op (
+ op.transpose_payload()
+ );
return_op.reinit_range_vector = op.reinit_domain_vector;
return_op.reinit_domain_vector = op.reinit_range_vector;
return return_op;
}
+
/**
* @relates LinearOperator
*
* of the @p solver object will be modified upon invocation of
* <code>vmult</code> or <code>Tvmult</code>.
*
+ * @author Luca Heltai, Matthias Maier, Jean-Paul Pelteret, 2015
+ *
* @ingroup LAOperators
*/
-template <typename Solver, typename Preconditioner>
-LinearOperator<typename Solver::vector_type, typename Solver::vector_type>
-inverse_operator(const LinearOperator<typename Solver::vector_type, typename Solver::vector_type> &op,
+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)
{
- typedef typename Solver::vector_type Vector;
-
- LinearOperator<Vector, Vector> return_op;
+ LinearOperator<Domain, Range, Payload> return_op (
+ op.inverse_payload(solver, preconditioner)
+ );
return_op.reinit_range_vector = op.reinit_domain_vector;
return_op.reinit_domain_vector = op.reinit_range_vector;
- return_op.vmult = [op, &solver, &preconditioner](Vector &v, const Vector &u)
+ return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u)
{
op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/ false);
solver.solve(op, v, u, preconditioner);
};
return_op.vmult_add =
- [op, &solver, &preconditioner](Vector &v, const Vector &u)
+ [op, &solver, &preconditioner](Range &v, const Domain &u)
{
- static GrowingVectorMemory<typename Solver::vector_type> vector_memory;
+ static GrowingVectorMemory<Range> vector_memory;
- Vector *v2 = vector_memory.alloc();
+ Range *v2 = vector_memory.alloc();
op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/ false);
solver.solve(op, *v2, u, preconditioner);
v += *v2;
vector_memory.free(v2);
};
- return_op.Tvmult = [op, &solver, &preconditioner]( Vector &v, const
- Vector &u)
+ return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u)
{
op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/ false);
solver.solve(transpose_operator(op), v, u, preconditioner);
};
return_op.Tvmult_add =
- [op, &solver, &preconditioner](Vector &v, const Vector &u)
+ [op, &solver, &preconditioner](Range &v, const Domain &u)
{
- static GrowingVectorMemory<typename Solver::vector_type> vector_memory;
+ static GrowingVectorMemory<Range> vector_memory;
- Vector *v2 = vector_memory.alloc();
+ Range *v2 = vector_memory.alloc();
op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/ false);
solver.solve(transpose_operator(op), *v2, u, preconditioner);
v += *v2;
*
* @ingroup LAOperators
*/
-template <typename Range>
-LinearOperator<Range, Range>
+template <typename Range,
+ typename Payload = internal::LinearOperator::EmptyPayload>
+LinearOperator<Range, Range, Payload>
identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
{
- LinearOperator<Range, Range> return_op;
+ LinearOperator<Range, Range, Payload> return_op ((Payload()));
return_op.reinit_range_vector = reinit_vector;
return_op.reinit_domain_vector = reinit_vector;
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-null_operator(const LinearOperator<Range, Domain> &op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+null_operator(const LinearOperator<Range, Domain, Payload> &op)
{
auto return_op = op;
v.reinit(matrix.n(), omit_zeroing_entries);
}
};
+
+
+ /**
+ * A dummy class for LinearOperators that do not require any extensions
+ * to facilitate the operations of the matrix.
+ *
+ * This is the Payload class typically associated with deal.II's native
+ * SparseMatrix. To use Trilinos and PETSc SparseMatrices it is necessary
+ * to initialize a LinearOperator with their associated Payload.
+ *
+ * @author Jean-Paul Pelteret, Matthias Maier, 2016
+ *
+ * @ingroup LAOperators
+ */
+ class EmptyPayload
+ {
+ public:
+
+ /**
+ * Default constructor
+ *
+ * Since this class does not do anything in particular and needs no special
+ * configuration, we have only one generic constructor that can be called
+ * under any conditions.
+ */
+ template <typename... Args>
+ EmptyPayload (const Args &...)
+ { }
+
+
+ /**
+ * Returns a payload configured for transpose operations
+ */
+ EmptyPayload
+ transpose_payload () const
+ {
+ return *this;
+ }
+
+
+ /**
+ * Returns a payload configured for inverse operations
+ */
+ template <typename Solver, typename Preconditioner>
+ EmptyPayload
+ inverse_payload (Solver &, const Preconditioner &) const
+ {
+ return *this;
+ }
+ };
+
+ /**
+ * Operator that returns a payload configured to support the
+ * addition of two LinearOperators
+ */
+ EmptyPayload operator+(const EmptyPayload &,
+ const EmptyPayload &)
+ {
+ return EmptyPayload();
+ }
+
+ /**
+ * Operator that returns a payload configured to support the
+ * multiplication of two LinearOperators
+ */
+ EmptyPayload operator*(const EmptyPayload &,
+ const EmptyPayload &)
+ {
+ return EmptyPayload();
+ }
+
} /* namespace LinearOperator */
} /* namespace internal */
// A helper class to add a reduced matrix interface to a LinearOperator
// (typically provided by Preconditioner classes)
- template <typename Range, typename Domain>
+ template <typename Range, typename Domain, typename Payload>
class MatrixInterfaceWithoutVmultAdd
{
public:
template <typename Matrix>
- void operator()(LinearOperator<Range, Domain> &op, const Matrix &matrix)
+ void operator()(LinearOperator<Range, Domain, Payload> &op, const Matrix &matrix)
{
op.vmult = [&matrix](Range &v, const Domain &u)
{
// A helper class to add the full matrix interface to a LinearOperator
- template <typename Range, typename Domain>
+ template <typename Range, typename Domain, typename Payload>
class MatrixInterfaceWithVmultAdd
{
public:
template <typename Matrix>
- void operator()(LinearOperator<Range, Domain> &op, const Matrix &matrix)
+ void operator()(LinearOperator<Range, Domain, Payload> &op, const Matrix &matrix)
{
// As above ...
- MatrixInterfaceWithoutVmultAdd<Range, Domain>().operator()(op, matrix);
+ MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>().operator()(op, matrix);
// ... but add native vmult_add and Tvmult_add variants:
*
* @ingroup LAOperators
*/
-template <typename Range, typename Domain, typename Matrix>
-LinearOperator<Range, Domain> linear_operator(const Matrix &matrix)
+template <typename Range, typename Domain,
+ typename Payload, typename Matrix>
+LinearOperator<Range, Domain, Payload>
+linear_operator(const Matrix &matrix)
{
// implement with the more generic variant below...
- return linear_operator<Range, Domain, Matrix, Matrix>(matrix, matrix);
+ return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix, matrix);
}
*/
template <typename Range,
typename Domain,
+ typename Payload,
typename OperatorExemplar,
typename Matrix>
-LinearOperator<Range, Domain>
+LinearOperator<Range, Domain, Payload>
linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
{
- LinearOperator<Range, Domain> return_op;
+ // Initialise the payload based on the input exemplar matrix
+ LinearOperator<Range, Domain, Payload> return_op (Payload(operator_exemplar,matrix));
// Always store a reference to matrix and operator_exemplar in the lambda
// functions. This ensures that a modification of the matrix after the
typename std::conditional<
has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
- MatrixInterfaceWithVmultAdd<Range, Domain>,
- MatrixInterfaceWithoutVmultAdd<Range, Domain>>::type().
+ MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
+ MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type().
operator()(return_op, matrix);
return return_op;
}
+
//@}
+
DEAL_II_NAMESPACE_CLOSE
#endif // DEAL_II_WITH_CXX11
projects / dealii.git / commitdiff