--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2010 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__linear_operator_h
+#define __deal2__linear_operator_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <functional>
+
+#ifdef DEAL_II_WITH_CXX11
+
+DEAL_II_NAMESPACE_OPEN
+
+// Forward declarations:
+
+template <typename Number> class Vector;
+
+template <typename Range, typename Domain> class LinearOperator;
+template <typename Range = Vector<double>,
+ typename Domain = Range,
+ typename Matrix>
+LinearOperator<Range, Domain> linop(const Matrix &matrix);
+
+/**
+ * A class to store the abstract concept of a linear operator.
+ *
+ * This is achieved by storing the concept of such operations in an
+ * "abstract" class LinearOperator that only holds knowledge of how to
+ * apply a linear operation via an std::function object vmult (vmult_add,
+ * Tvmult, Tvmult_add), as well as, information on how to create an
+ * appropriate Domain and Range vector.
+ *
+ * The primary purpose of this class is to provide syntactic sugar for
+ * complex matrix-vector operations and free the user from having to
+ * create, set up and handle intermediate storage locations by hand.
+ *
+ * As an example consider the operation $(A+k\,B)\,C$, where $A$, $B$ and
+ * $C$ denote (possible different) matrices. In order to construct a
+ * LinearOperator <code>op</code> that stores the knowledge of this
+ * operation, one can write:
+ *
+ * @code
+ * dealii::SparseMatrix<double> A, B, C;
+ * const double k = ...;
+ *
+ * // Setup and assembly of matrices
+ *
+ * const auto op_a = linop(A);
+ * const auto op_b = linop(B);
+ * const auto op_c = linop(C);
+ *
+ * const auto op = (op_a + k * op_b) * op_c;
+ * @endcode
+ *
+ * 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
+ */
+template <typename Range, typename Domain> class LinearOperator
+{
+public:
+
+ /**
+ * Create an empty LinearOperator object. All std::function objects are
+ * initialized with default variants that throw an exception upon
+ * invocation.
+ */
+ LinearOperator ()
+ {
+ vmult = [](Range &, const Domain &)
+ {
+ Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+ "Domain>::vmult called"));
+ };
+
+ vmult_add = [](Range &, const Domain &)
+ {
+ Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+ "Domain>::vmult_add called"));
+ };
+
+ Tvmult = [](Domain &, const Range &)
+ {
+ Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+ "Domain>::Tvmult called"));
+ };
+
+ Tvmult_add = [](Domain &, const Range &)
+ {
+ Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+ "Domain>::Tvmult_add called"));
+ };
+
+ reinit_range_vector = [](Range &, bool)
+ {
+ Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+ "Domain>::reinit_range_vector method called"));
+ };
+
+ reinit_domain_vector = [](Domain &, bool)
+ {
+ Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+ "Domain>::reinit_domain_vector method called"));
+ };
+ }
+
+ /**
+ * Create a LinearOperator object from an object @p op for which the
+ * conversion function <code>linop</code> is defined.
+ */
+ template<typename Op>
+ LinearOperator (const Op &op)
+ {
+ *this = linop<Range, Domain, Op>(op);
+ }
+
+ /**
+ * Copy assignment operator for an object @p op for which the conversion
+ * function <code>linop</code> is defined.
+ */
+ template <typename Op>
+ LinearOperator<Range, Domain> &operator=(const Op &op)
+ {
+ return *this = linop<Range, Domain, Op>(op);
+ }
+
+ /**
+ * Application of the LinearOperator object to a vector @p u of the
+ * Domain space giving a vector @p v of the Range space.
+ */
+ std::function<const void(Range &v, const Domain &u)> vmult;
+
+ /**
+ * Application of the LinearOperator object to a vector @p u of the
+ * Domain space. The result is added to the vector @p v.
+ */
+ std::function<const void(Range &v, const Domain &u)> vmult_add;
+
+ /**
+ * Application of the transpose LinearOperator object to a vector @p u of
+ * the Range space giving a vector @p v of the Domain space.
+ */
+ std::function<const void(Domain &v, const Range &u)> Tvmult;
+
+ /**
+ * Application of the transpose LinearOperator object to a vector @p u of
+ * the Range space. The result is added to the vector @p v.
+ */
+ std::function<const void(Domain &v, const Range &u)> Tvmult_add;
+
+ /**
+ * Initializes a vector 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 initalization is done, i.e., if it is set to false the content
+ * of the vector is set to 0.
+ */
+ std::function<void(Range &, bool)> reinit_range_vector;
+
+ /**
+ * 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
+ * initalization is done, i.e., if it is set to false the content of the
+ * vector is set to 0.
+ */
+ std::function<void(Domain &, bool)> reinit_domain_vector;
+
+ /**
+ * A memory object for vectors of Range space used for intermediate
+ * storage during computations in vmult, etc.
+ */
+ mutable GrowingVectorMemory<Range> range_vector_memory;
+
+ /**
+ * A memory object for vectors of Domain space used for intermediate
+ * storage during computations in vmult, etc.
+ */
+ mutable GrowingVectorMemory<Domain> domain_vector_memory;
+
+
+ /**
+ * Addition with a LinearOperator @p second_op with the same Domain
+ * and Range.
+ */
+ LinearOperator<Range, Domain> &
+ operator+=(const LinearOperator<Range, Domain> &second_op)
+ {
+ *this = *this + second_op;
+ return *this;
+ }
+
+ /**
+ * Subtraction with a LinearOperator @p second_op with the same Domain
+ * and Range.
+ */
+ LinearOperator<Range, Domain> &
+ operator-=(const LinearOperator<Range, Domain> &second_op)
+ {
+ *this = *this - second_op;
+ return *this;
+ }
+
+ /**
+ * Concatenation of the LinearOperator with an endormorphism @p second_op
+ * on the Domain space.
+ */
+ LinearOperator<Range, Domain> &
+ operator*=(const LinearOperator<Domain, Domain> &second_op)
+ {
+ *this = *this * second_op;
+ return *this;
+ }
+};
+
+
+
+/**
+ * Addition of two linear operators @p first_op and @p second_op given
+ * by $(\text{first_op}+\text{second_op})x:=\text{first_op}(x)+\text{second_op}(x)$
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator+(const LinearOperator<Range, Domain> &first_op,
+ const LinearOperator<Range, Domain> &second_op)
+{
+ LinearOperator<Range, Domain> return_op;
+
+ return_op.reinit_range_vector = first_op.reinit_range_vector;
+ return_op.reinit_domain_vector = first_op.reinit_domain_vector;
+
+ // ensure to have valid computation objects by catching first_op and
+ // second_op by value
+
+ return_op.vmult = [first_op, second_op](Range &v, const Domain &u)
+ {
+ first_op.vmult(v, u);
+ second_op.vmult_add(v, u);
+ };
+
+ return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u)
+ {
+ first_op.vmult_add(v, u);
+ second_op.vmult_add(v, u);
+ };
+
+ return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u)
+ {
+ second_op.Tvmult(v, u);
+ first_op.Tvmult_add(v, u);
+ };
+
+ return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u)
+ {
+ second_op.Tvmult_add(v, u);
+ first_op.Tvmult_add(v, u);
+ };
+
+ return return_op;
+}
+
+
+/**
+ * Subtraction of two linear operators @p first_op and @p second_op given
+ * by $(\text{first_op}-\text{second_op})x:=\text{first_op}(x)-\text{second_op}(x)$
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator-(const LinearOperator<Range, Domain> &first_op,
+ const LinearOperator<Range, Domain> &second_op)
+{
+ // implement with addition and scalar multiplication
+ return first_op + (-1. * second_op);
+}
+
+
+/**
+ * Concatenation of two linear operators @p first_op and @p second_op given
+ * by $(\text{first_op}*\text{second_op})x:=\text{first_op}(\text{second_op}(x))$
+ */
+template <typename Range, typename Intermediate, typename Domain>
+LinearOperator<Range, Domain>
+operator*(const LinearOperator<Range, Intermediate> &first_op,
+ const LinearOperator<Intermediate, Domain> &second_op)
+{
+ LinearOperator<Range, Domain> return_op;
+
+ return_op.reinit_domain_vector = second_op.reinit_domain_vector;
+ return_op.reinit_range_vector = first_op.reinit_range_vector;
+
+ // ensure to have valid computation objects by catching first_op and
+ // second_op by value
+
+ // Note: The range_vector_memory and domain_vector_memory objects inside
+ // the lambda exrepessions are the one belonging to the by-value captured
+ // objects first_op and second_op. Thus, they belong in essence to the
+ // std::function object vmult (vmult_add, ...) - and therefore have the
+ // same lifetime as the function object and not the LinearOperator
+ // parameters first_op and second_op.
+
+ return_op.vmult = [first_op, second_op](Range &v, const Domain &u)
+ {
+ Intermediate *i = second_op.range_vector_memory.alloc();
+ second_op.reinit_range_vector(*i, /*bool fast =*/ true);
+ second_op.vmult(*i, u);
+ first_op.vmult(v, *i);
+ second_op.range_vector_memory.free(i);
+ };
+
+ return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u)
+ {
+ Intermediate *i = second_op.range_vector_memory.alloc();
+ second_op.reinit_range_vector(*i, /*bool fast =*/ true);
+ second_op.vmult(*i, u);
+ first_op.vmult_add(v, *i);
+ second_op.range_vector_memory.free(i);
+ };
+
+ return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u)
+ {
+ Intermediate *i = first_op.domain_vector_memory.alloc();
+ first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
+ first_op.Tvmult(*i, u);
+ second_op.Tvmult(v, *i);
+ first_op.domain_vector_memory.free(i);
+ };
+
+ return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u)
+ {
+ Intermediate *i = first_op.domain_vector_memory.alloc();
+ first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
+ first_op.Tvmult(*i, u);
+ second_op.Tvmult_add(v, *i);
+ first_op.domain_vector_memory.free(i);
+ };
+
+ return return_op;
+}
+
+
+/**
+ * Scalar multiplication of a ScalarOperator object from the left.
+ *
+ * The Domain and Range types must implement the following
+ * <code>operator*=</code> member functions accepting the appropriate
+ * scalar Number type for rescaling:
+ *
+ * @code
+ * Domain & operator *=(Number);
+ * Range & operator *=(Number);
+ * @endcode
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator*(typename Range::value_type number,
+ const LinearOperator<Range, Domain> &op)
+{
+ static_assert(
+ std::is_convertible<typename Range::value_type, typename Domain::value_type>::value,
+ "Range and Domain must have implicitly convertible 'value_type's");
+
+ LinearOperator<Range, Domain> return_op = op;
+
+ // ensure to have valid computation objects by catching number and op by
+ // value
+
+ return_op.vmult = [number, op](Range &v, const Domain &u)
+ {
+ op.vmult(v,u);
+ v *= number;
+ };
+
+ return_op.vmult_add = [number, op](Range &v, const Domain &u)
+ {
+ v /= number;
+ op.vmult_add(v,u);
+ v *= number;
+ };
+
+ return_op.Tvmult = [number, op](Domain &v, const Range &u)
+ {
+ op.Tvmult(v,u);
+ v *= number;
+ };
+
+ return_op.Tvmult_add = [number, op](Domain &v, const Range &u)
+ {
+ v /= number;
+ op.Tvmult_add(v,u);
+ v *= number;
+ };
+
+ return return_op;
+}
+
+
+/**
+ * Scalar multiplication of a ScalarOperator object from the right.
+ *
+ * The Domain and Range types must implement the following
+ * <code>operator*=</code> member functions for rescaling:
+ *
+ * @code
+ * Domain & operator *=(Number);
+ * Range & operator *=(Number);
+ * @endcode
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator*(const LinearOperator<Range, Domain> &op,
+ typename Domain::value_type number)
+{
+ static_assert(
+ std::is_convertible<typename Range::value_type, typename Domain::value_type>::value,
+ "Range and Domain must have implicitly convertible 'value_type's");
+
+ return number * op;
+}
+
+
+/**
+ * Returns the transpose linear operations of @p op.
+ */
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+transpose_linop(const LinearOperator<Range, Domain> &op)
+{
+ LinearOperator<Domain, Range> return_op;
+
+ return_op.reinit_range_vector = op.reinit_domain_vector;
+ return_op.reinit_domain_vector = op.reinit_range_vector;
+
+ return_op.vmult = op.Tvmult;
+ return_op.vmult_add = op.Tvmult_add;
+ return_op.Tvmult = op.vmult;
+ return_op.Tvmult_add = op.vmult_add;
+
+ return return_op;
+}
+
+
+/**
+ * Returns a LinearOperator that is the identity of the vector space
+ * Domain.
+ *
+ * The function takes an <code>std::function</code> object @p
+ * exemplar as an argument to initialize the reinit_range_vector
+ * and reinit_domain_vector objects of the LinearOperator object.
+ */
+template <typename Range>
+LinearOperator<Range, Range>
+identity_linop(const std::function<void(Range &, bool)> &exemplar)
+{
+ LinearOperator<Range, Range> return_op;
+
+ return_op.reinit_range_vector = exemplar;
+ return_op.reinit_domain_vector = exemplar;
+
+ return_op.vmult = [](Range &v, const Range &u)
+ {
+ v = u;
+ };
+
+ return_op.vmult_add = [](Range &v, const Range &u)
+ {
+ v += u;
+ };
+
+ return_op.Tvmult = [](Range &v, const Range &u)
+ {
+ v = u;
+ };
+
+ return_op.Tvmult_add = [](Range &v, const Range &u)
+ {
+ v += u;
+ };
+
+ return return_op;
+}
+
+
+/**
+ * Returns an object representing the inverse of the LinearOperator @p op.
+ *
+ * The function takes references @p solver and @p preconditioner to an
+ * iterative solver and a preconditioner that are used in the
+ * <code>vmult</code> and <code>Tvmult</code> implementations of the
+ * LinearOperator object.
+ *
+ * The LinearOperator object that is created stores a reference to @p
+ * solver and @p preconditioner. Thus, both objects must remain a valid
+ * reference for the whole lifetime of the LinearOperator object. Internal
+ * data structures of the @p solver object will be modified upon invocation
+ * of <code>vmult</code> or <code>Tvmult</code>.
+ */
+template <typename Solver, typename Preconditioner>
+LinearOperator<typename Solver::vector_type, typename Solver::vector_type>
+inverse_linop(Solver &solver,
+ const Preconditioner &preconditioner,
+ const LinearOperator<typename Solver::vector_type, typename Solver::vector_type> &op)
+{
+ typedef typename Solver::vector_type Vector;
+
+ LinearOperator<Vector, Vector> return_op;
+
+ 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)
+ {
+ solver.solve(op, v, u, preconditioner);
+ };
+
+ // Note: The range_vector_memory and domain_vector_memory objects inside
+ // the lambda exrepessions are the one belonging to the by-value captured
+ // objects first_op and second_op. Thus, they belong in essence to the
+ // std::function object vmult (vmult_add, ...) - and therefore have the
+ // same lifetime as the function object and not the LinearOperator
+ // parameter op.
+
+ return_op.vmult_add =
+ [op, &solver, &preconditioner](Vector &v, const Vector &u)
+ {
+ Vector *v2 = op.range_vector_memory.alloc();
+ op.reinit_range_vector(*v2, /*bool fast =*/ true);
+ solver.solve(op, *v2, u, preconditioner);
+ v += *v2;
+ op.range_vector_memory.free(v2);
+ };
+
+ return_op.Tvmult = [op, &solver, &preconditioner]( Vector &v, const
+ Vector &u)
+ {
+ solver.solve(transpose_linop(op), v, u, preconditioner);
+ };
+
+ return_op.Tvmult =
+ [op, &solver, &preconditioner](Vector &v, const Vector &u)
+ {
+ Vector *v2 = op.range_vector_memory.alloc();
+ op.reinit_range_vector(*v2, /*bool fast =*/ true);
+ solver.solve(transpose_linop(op), *v2, u, preconditioner);
+ v += *v2;
+ op.range_vector_memory.free(v2);
+ };
+
+ return return_op;
+}
+
+
+/**
+ * A factory class that is responsible to create a reinit_range_vector
+ * object for a given pair of vector type Range and matrix type Matrix.
+ *
+ * The generic version of this class just calls Range::reinit() with the
+ * result of Matrix::m(). This class is specialized for more complicated
+ * data structures, such as TrilinosWrappers::MPI::Vector, etc.
+ */
+template<typename Range>
+class ReinitRangeFactory
+{
+public:
+
+ template <typename Matrix>
+ std::function<void(Range &, bool)>
+ operator()(const Matrix &matrix)
+ {
+ return [&matrix](Range &v, bool fast)
+ {
+ v.reinit(matrix.m(), fast);
+ };
+ }
+};
+
+
+/**
+ * A factory class that is responsible to create a reinit_domain_vector
+ * object for a given pair of vector type Domain and matrix type Matrix.
+ *
+ * The generic version of this class just calls Domain::reinit() with the
+ * result of Matrix::n(). This class is specialized for more complicated
+ * data structures, such as TrilinosWrappers::MPI::Vector, etc.
+ */
+template<typename Domain>
+class ReinitDomainFactory
+{
+public:
+
+ template <typename Matrix>
+ std::function<void(Domain &, bool)>
+ operator()(const Matrix &matrix)
+ {
+ return [&matrix](Domain &v, bool fast)
+ {
+ v.reinit(matrix.n(), fast);
+ };
+ }
+};
+
+
+/**
+ * A function that encapsulates generic @p 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, @p matrix
+ * must remain a valid reference for the whole lifetime of the
+ * LinearOperator object.
+ *
+ * All changes made on @p 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:
+ *
+ * @code
+ * 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(VECTOR &dst, const VECTOR &src);
+ *
+ * // Application of matrix to vector src, adds the result to dst.
+ * vmult_add(VECTOR &dst, const VECTOR &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(VECTOR &dst, const VECTOR &src);
+ *
+ * // Application of the transpose of matrix to vector src, adds the
+ * // result to dst. (Depending on the usage of the linear operator class
+ * // this can be a dummy implementation throwing an error.)
+ * Tvmult_add(VECTOR &dst, const VECTOR &src);
+ * };
+ * @endcode
+ *
+ * @author: Matthias Maier, 2015
+ */
+template <typename Range = Vector<double>,
+ typename Domain = Range,
+ typename Matrix>
+LinearOperator<Range, Domain> linop(const Matrix &matrix)
+{
+ LinearOperator<Range, Domain> return_op;
+
+ // Always store a reference to matrix in the lambda functions.
+ // This ensures that a modification of the matrix after the creation of a
+ // LinearOperator wrapper is respected - further a matrix cannot usually
+ // be copied...
+
+ return_op.reinit_range_vector =
+ ReinitRangeFactory<Range>().operator()(matrix);
+
+ return_op.reinit_domain_vector =
+ ReinitDomainFactory<Domain>().operator()(matrix);
+
+ return_op.vmult = [&matrix](Range &v, const Domain &u)
+ {
+ matrix.vmult(v,u);
+ };
+
+ return_op.vmult_add = [&matrix](Range &v, const Domain &u)
+ {
+ matrix.vmult_add(v,u);
+ };
+
+ return_op.Tvmult = [&matrix](Domain &v, const Range &u)
+ {
+ matrix.Tvmult(v,u);
+ };
+
+ return_op.Tvmult_add = [&matrix](Domain &v, const Range &u)
+ {
+ matrix.Tvmult_add(v,u);
+ };
+
+ return return_op;
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_CXX11
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test the LinearOperator template on a trivial vector implementation
+// :: RightVector -> LeftVector
+
+#include "../tests.h"
+
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/vector_memory.templates.h>
+
+using namespace dealii;
+
+// Dummy vectors with different, non convertible types:
+
+struct LeftVector
+{
+ typedef double value_type;
+ double value;
+
+ LeftVector & operator *= (double scale)
+ {
+ value *= scale;
+ return *this;
+ }
+ LeftVector & operator /= (double scale)
+ {
+ value /= scale;
+ return *this;
+ }
+ LeftVector & operator += (const LeftVector &u)
+ {
+ value += u.value;
+ return *this;
+ }
+ int size() const
+ {
+ return 1;
+ }
+ std::size_t memory_consumption () const
+ {
+ return 1;
+ }
+};
+
+struct RightVector
+{
+ typedef double value_type;
+ double value;
+
+ RightVector & operator *= (double scale)
+ {
+ value *= scale;
+ return *this;
+ }
+ RightVector & operator /= (double scale)
+ {
+ value /= scale;
+ return *this;
+ }
+ RightVector & operator += (const RightVector &u)
+ {
+ value += u.value;
+ return *this;
+ }
+ int size() const
+ {
+ return 1;
+ }
+ std::size_t memory_consumption () const
+ {
+ return 1;
+ }
+};
+
+int main()
+{
+ initlog();
+
+ // Create to distinct linear operators:
+
+ LinearOperator<LeftVector, RightVector> multiply2;
+ multiply2.vmult = [](LeftVector &v, const RightVector &u)
+ {
+ v.value = 2 * u.value;
+ };
+ multiply2.vmult_add = [](LeftVector &v, const RightVector &u)
+ {
+ v.value += 2 * u.value;
+ };
+ multiply2.Tvmult = [](RightVector &v, const LeftVector &u)
+ {
+ v.value = 2 * u.value;
+ };
+ multiply2.Tvmult_add = [](RightVector &v, const LeftVector &u)
+ {
+ v.value += 2 * u.value;
+ };
+ multiply2.reinit_range_vector = [](LeftVector &, bool)
+ {
+ };
+ multiply2.reinit_domain_vector = [](RightVector &, bool)
+ {
+ };
+
+ auto multiply4 = multiply2;
+ multiply4.vmult = [](LeftVector &v, const RightVector &u)
+ {
+ v.value = 4 * u.value;
+ };
+ multiply4.vmult_add = [](LeftVector &v, const RightVector &u)
+ {
+ v.value += 4 * u.value;
+ };
+ multiply4.Tvmult = [](RightVector &v, const LeftVector &u)
+ {
+ v.value = 4 * u.value;
+ };
+ multiply4.Tvmult_add = [](RightVector &v, const LeftVector &u)
+ {
+ v.value += 4 * u.value;
+ };
+
+
+ // Small unit tests for all functions:
+
+ RightVector u = { 4. };
+ LeftVector v = { 0. };
+
+ // vmult, vmult_add
+
+ multiply2.vmult(v, u);
+ deallog << "2 * " << u.value << " = " << v.value << std::endl;
+
+ multiply4.vmult_add(v, u);
+ deallog << "... + 4 * " << u.value << " = " << v.value << std::endl;
+
+ multiply4.vmult(v, u);
+ deallog << "4 * " << u.value << " = " << v.value << std::endl;
+
+ multiply2.vmult_add(v, u);
+ deallog << "... + 2 * " << u.value << " = " << v.value << std::endl;
+
+ // Tvmult, Tvmult_add
+
+ v.value = 4.;
+
+ multiply2.Tvmult(u, v);
+ deallog << "2 * " << v.value << " = " << u.value << std::endl;
+
+ multiply4.Tvmult_add(u, v);
+ deallog << "... + 4 * " << v.value << " = " << u.value << std::endl;
+
+ multiply4.Tvmult(u, v);
+ deallog << "4 * " << v.value << " = " << u.value << std::endl;
+
+ multiply2.Tvmult_add(u, v);
+ deallog << "... + 2 * " << v.value << " = " << u.value << std::endl;
+
+ // operator+, operator-, operator+=, operator-=
+
+ auto test = multiply2 + multiply4;
+ test.vmult(v, u);
+ deallog << "(2 + 4) * " << u.value << " = " << v.value << std::endl;
+
+ test = multiply2 - multiply4;
+ test.vmult(v, u);
+ deallog << "(2 - 4) * " << u.value << " = " << v.value << std::endl;
+
+ test += multiply2;
+ test.vmult(v, u);
+ deallog << "(2 - 4 + 2) * " << u.value << " = " << v.value << std::endl;
+
+ test -= multiply4;
+ test.vmult(v, u);
+ deallog << "(2 - 4 + 2 - 4) * " << u.value << " = " << v.value << std::endl;
+
+ // operator* with scalar
+
+ test = 4. * multiply4;
+ test.vmult(v, u);
+ deallog << "(4 * 4) * " << u.value << " = " << v.value << std::endl;
+
+ test = multiply4 * 4.;
+ test.vmult(v, u);
+ deallog << "(4 * 4) * " << u.value << " = " << v.value << std::endl;
+
+ // operator* and transpose
+
+ auto test2 = transpose_linop(multiply2) * multiply4;
+ RightVector w = { 0. };
+ test2.vmult(w, u);
+ deallog << "(2 * 4) * " << u.value << " = " << w.value << std::endl;
+
+ // identity
+
+ auto test3 = identity_linop(test2.reinit_range_vector) + test2;
+ test3.vmult(w, u);
+ deallog << "(1 + 2 * 4) * " << u.value << " = " << w.value << std::endl;
+}
--- /dev/null
+
+DEAL::2 * 4.00000 = 8.00000
+DEAL::... + 4 * 4.00000 = 24.0000
+DEAL::4 * 4.00000 = 16.0000
+DEAL::... + 2 * 4.00000 = 24.0000
+DEAL::2 * 4.00000 = 8.00000
+DEAL::... + 4 * 4.00000 = 24.0000
+DEAL::4 * 4.00000 = 16.0000
+DEAL::... + 2 * 4.00000 = 24.0000
+DEAL::(2 + 4) * 24.0000 = 144.000
+DEAL::(2 - 4) * 24.0000 = -48.0000
+DEAL::(2 - 4 + 2) * 24.0000 = 0.00000
+DEAL::(2 - 4 + 2 - 4) * 24.0000 = -96.0000
+DEAL::(4 * 4) * 24.0000 = 384.000
+DEAL::(4 * 4) * 24.0000 = 384.000
+DEAL::(2 * 4) * 24.0000 = 192.000
+DEAL::(1 + 2 * 4) * 24.0000 = 216.000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Tests for the LinearOperator template with
+// dealii::Vector<double>
+// dealii::SparseMatrix<double>
+// dealii::FullMatrix<double>
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+using namespace dealii;
+
+int main()
+{
+ initlog();
+ deallog << std::setprecision(10);
+
+ static const int dim = 2;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global(2);
+
+ MappingQ1<dim> mapping_q1;
+ FE_Q<dim> q1(1);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(q1);
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp);
+ dsp.compress();
+ SparsityPattern sparsity_pattern;
+ sparsity_pattern.copy_from(dsp);
+ sparsity_pattern.compress();
+
+ SparseMatrix<double> a (sparsity_pattern);
+ SparseMatrix<double> b (sparsity_pattern);
+
+ QGauss<dim> quadrature(4);
+ MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, a);
+ MatrixCreator::create_mass_matrix(mapping_q1, dof_handler, quadrature, b);
+
+
+ // Constructors and assignment:
+
+ auto op_a = linop(a);
+ auto op_b = linop(b);
+
+ {
+ LinearOperator<dealii::Vector<double>, dealii::Vector<double>> op_x (a);
+ op_a = a;
+ op_b = b;
+ }
+
+ // vmult:
+
+ Vector<double> u;
+ op_a.reinit_domain_vector(u, true);
+ for (unsigned int i = 0; i < u.size(); ++i) {
+ u[i] = (double)(i+1);
+ }
+
+ deallog << "u: " << u << std::endl;
+
+ Vector<double> v;
+ op_a.reinit_domain_vector(v, false);
+ Vector<double> w;
+ op_a.reinit_domain_vector(w, false);
+ Vector<double> x;
+ op_a.reinit_domain_vector(x, false);
+
+ op_a.vmult(v, u);
+ deallog << "Au: " << v << std::endl;
+
+ op_b.vmult(w, u);
+ deallog << "Bu: " << w << std::endl;
+
+ // operator+, operator-, operator+=, operator-=:
+
+ x = v;
+ x += w;
+ deallog << "Au+Bu: " << x << std::endl;
+
+ (op_a + op_b).vmult(x, u);
+ deallog << "(A+B)u: " << x << std::endl;
+
+ auto op_x = op_a;
+ op_x += op_b;
+ op_x.vmult(x, u);
+ deallog << "(A+=B)u: " << x << std::endl;
+
+ x = v;
+ x -= w;
+ deallog << "Au-Bu: " << x << std::endl;
+
+ (op_a - op_b).vmult(x, u);
+ deallog << "(A-B)u: " << x << std::endl;
+
+ op_x = op_a;
+ op_x -= op_b;
+ op_x.vmult(x, u);
+ deallog << "(A-=B)u: " << x << std::endl;
+
+ // operator*, operator*=
+
+ op_b.vmult(v,u);
+ op_a.vmult(w,v);
+ deallog << "(A(Bu)): " << w << std::endl;
+
+ (op_a * op_b).vmult(x, u);
+ deallog << "(A*B)u: " << x << std::endl;
+
+ op_x = op_a;
+ op_x *= op_b;
+ op_x.vmult(x, u);
+ deallog << "(A*=B)u: " << x << std::endl;
+
+ // solver interface:
+
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> solver (solver_control);
+
+ deallog.depth_file(0);
+ solver.solve(op_b, v, u, PreconditionIdentity());
+ deallog.depth_file(3);
+ deallog << "solve(B, v, u): " << v << std::endl;
+
+ deallog.depth_file(0);
+ inverse_linop(solver, PreconditionIdentity(), op_b).vmult(v, u);
+ deallog.depth_file(3);
+ deallog << "inverse_linop(B)u: " << v << std::endl;
+
+ deallog.depth_file(0);
+ op_b.vmult(w, v);
+ deallog.depth_file(3);
+ deallog << "B(inverse_linop(B)u): " << w << std::endl;
+
+ deallog.depth_file(0);
+ (op_b * inverse_linop(solver, PreconditionIdentity(), op_b)).vmult(w, u);
+ deallog.depth_file(3);
+ deallog << "(B*inverse_linop(B))u: " << w << std::endl;
+
+ deallog.depth_file(0);
+ (inverse_linop(solver, PreconditionIdentity(), op_b) * op_b).vmult(w, u);
+ deallog.depth_file(3);
+ deallog << "(inverse_linop(B)*B)u: " << w << std::endl;
+
+ SolverControl inner_solver_control (1000, 1e-12);
+ SolverCG<> inner_solver (solver_control);
+
+ deallog.depth_file(0);
+ solver.solve(inverse_linop(inner_solver, PreconditionIdentity(), op_b), v, u,
+ PreconditionIdentity());
+ deallog.depth_file(3);
+ deallog << "solve(inverse_linop(B), v, u) == Bu: " << v << std::endl;
+
+ deallog.depth_file(0);
+ solver.solve(op_b + 0.005 * op_a, v, u, PreconditionIdentity());
+ deallog.depth_file(3);
+ deallog << "solve(B+0.5*0.01*A, v, u): " << v << std::endl;
+}
--- /dev/null
+
+DEAL::u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL::
+DEAL::Au: -1.500000000 -2.666666667 -2.000000000 -3.000000000 -2.333333333 -5.000000000 -3.500000000 -5.333333333 -9.333333333 0.5000000000 1.333333333 0.5000000000 1.166666667 -1.666666667 -1.666666667 2.000000000 6.000000000 3.000000000 1.000000000 3.000000000 1.666666667 9.000000000 4.000000000 3.333333333 1.500000000
+DEAL::
+DEAL::Bu: 0.03125000000 0.09201388889 0.1145833333 0.2812500000 0.1788194444 0.4270833333 0.2552083333 0.5538194444 0.6631944444 0.3072916667 0.6753472222 0.1822916667 0.3923611111 0.8888888889 0.4878472222 0.4791666667 1.000000000 1.093750000 0.2864583333 0.5937500000 0.6371527778 1.281250000 0.6770833333 0.7170138889 0.3750000000
+DEAL::
+DEAL::Au+Bu: -1.468750000 -2.574652778 -1.885416667 -2.718750000 -2.154513889 -4.572916667 -3.244791667 -4.779513889 -8.670138889 0.8072916667 2.008680556 0.6822916667 1.559027778 -0.7777777778 -1.178819444 2.479166667 7.000000000 4.093750000 1.286458333 3.593750000 2.303819444 10.28125000 4.677083333 4.050347222 1.875000000
+DEAL::
+DEAL::(A+B)u: -1.468750000 -2.574652778 -1.885416667 -2.718750000 -2.154513889 -4.572916667 -3.244791667 -4.779513889 -8.670138889 0.8072916667 2.008680556 0.6822916667 1.559027778 -0.7777777778 -1.178819444 2.479166667 7.000000000 4.093750000 1.286458333 3.593750000 2.303819444 10.28125000 4.677083333 4.050347222 1.875000000
+DEAL::
+DEAL::(A+=B)u: -1.468750000 -2.574652778 -1.885416667 -2.718750000 -2.154513889 -4.572916667 -3.244791667 -4.779513889 -8.670138889 0.8072916667 2.008680556 0.6822916667 1.559027778 -0.7777777778 -1.178819444 2.479166667 7.000000000 4.093750000 1.286458333 3.593750000 2.303819444 10.28125000 4.677083333 4.050347222 1.875000000
+DEAL::
+DEAL::Au-Bu: -1.531250000 -2.758680556 -2.114583333 -3.281250000 -2.512152778 -5.427083333 -3.755208333 -5.887152778 -9.996527778 0.1927083333 0.6579861111 0.3177083333 0.7743055556 -2.555555556 -2.154513889 1.520833333 5.000000000 1.906250000 0.7135416667 2.406250000 1.029513889 7.718750000 3.322916667 2.616319444 1.125000000
+DEAL::
+DEAL::(A-B)u: -1.531250000 -2.758680556 -2.114583333 -3.281250000 -2.512152778 -5.427083333 -3.755208333 -5.887152778 -9.996527778 0.1927083333 0.6579861111 0.3177083333 0.7743055556 -2.555555556 -2.154513889 1.520833333 5.000000000 1.906250000 0.7135416667 2.406250000 1.029513889 7.718750000 3.322916667 2.616319444 1.125000000
+DEAL::
+DEAL::(A-=B)u: -1.531250000 -2.758680556 -2.114583333 -3.281250000 -2.512152778 -5.427083333 -3.755208333 -5.887152778 -9.996527778 0.1927083333 0.6579861111 0.3177083333 0.7743055556 -2.555555556 -2.154513889 1.520833333 5.000000000 1.906250000 0.7135416667 2.406250000 1.029513889 7.718750000 3.322916667 2.616319444 1.125000000
+DEAL::
+DEAL::(A(Bu)): -0.1073495370 -0.1866319444 -0.2039930556 -0.02199074074 -0.2893518519 -0.07465277778 -0.3703703704 0.03877314815 -0.2986111111 -0.1487268519 0.6250000000 -0.2201967593 -0.2123842593 0.4710648148 -0.4762731481 -0.1672453704 1.145833333 0.8049768519 -0.3211805556 -0.2199074074 -0.4939236111 1.570023148 -0.2034143519 -0.2300347222 -0.4094328704
+DEAL::
+DEAL::(A*B)u: -0.1073495370 -0.1866319444 -0.2039930556 -0.02199074074 -0.2893518519 -0.07465277778 -0.3703703704 0.03877314815 -0.2986111111 -0.1487268519 0.6250000000 -0.2201967593 -0.2123842593 0.4710648148 -0.4762731481 -0.1672453704 1.145833333 0.8049768519 -0.3211805556 -0.2199074074 -0.4939236111 1.570023148 -0.2034143519 -0.2300347222 -0.4094328704
+DEAL::
+DEAL::(A*=B)u: -0.1073495370 -0.1866319444 -0.2039930556 -0.02199074074 -0.2893518519 -0.07465277778 -0.3703703704 0.03877314815 -0.2986111111 -0.1487268519 0.6250000000 -0.2201967593 -0.2123842593 0.4710648148 -0.4762731481 -0.1672453704 1.145833333 0.8049768519 -0.3211805556 -0.2199074074 -0.4939236111 1.570023148 -0.2034143519 -0.2300347222 -0.4094328704
+DEAL::
+DEAL::solve(B, v, u): 80.16326531 28.24489796 86.53061224 25.79591837 259.4285714 57.14285714 328.0000000 88.00000000 184.0000000 226.6122449 40.48979592 1446.693878 315.7551020 160.0000000 760.0000000 387.7551020 81.63265306 190.8571429 2286.693878 470.0408163 1140.571429 108.0816327 569.9591837 539.1020408 3018.448980
+DEAL::
+DEAL::inverse_linop(B)u: 80.16326531 28.24489796 86.53061224 25.79591837 259.4285714 57.14285714 328.0000000 88.00000000 184.0000000 226.6122449 40.48979592 1446.693878 315.7551020 160.0000000 760.0000000 387.7551020 81.63265306 190.8571429 2286.693878 470.0408163 1140.571429 108.0816327 569.9591837 539.1020408 3018.448980
+DEAL::
+DEAL::B(inverse_linop(B)u): 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL::
+DEAL::(B*inverse_linop(B))u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL::
+DEAL::(inverse_linop(B)*B)u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL::
+DEAL::solve(inverse_linop(B), v, u) == Bu: 0.03125000000 0.09201388889 0.1145833333 0.2812500000 0.1788194444 0.4270833333 0.2552083333 0.5538194444 0.6631944444 0.3072916667 0.6753472222 0.1822916667 0.3923611111 0.8888888889 0.4878472222 0.4791666667 1.000000000 1.093750000 0.2864583333 0.5937500000 0.6371527778 1.281250000 0.6770833333 0.7170138889 0.3750000000
+DEAL::
+DEAL::solve(B+0.5*0.01*A, v, u): 59.31059595 52.77843700 93.77361796 53.05371114 176.7862626 80.04635988 236.6925420 110.0293599 135.7027276 299.7370750 139.0997515 987.0719769 399.5852319 197.1836670 535.4087997 487.8703686 221.2510235 251.7282805 1564.703948 612.1971228 783.5848738 287.5180747 708.4784200 725.1386389 2062.398876
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Tests for the LinearOperator template with
+// dealii::BlockVector<double>
+// dealii::BlockSparseMatrix<double>
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#define PRINTME(name, var) \
+ deallog << name << ": [block 0] " << var.block(0) \
+ << " [block 1] " << var.block(1) << std::endl;
+
+
+using namespace dealii;
+
+int main()
+{
+ initlog();
+ deallog << std::setprecision(10);
+
+ static const int dim = 2;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global(2);
+
+ MappingQ1<dim> mapping_q1;
+ FESystem<dim> fe(FE_Q<dim>(1), 1, FE_Q<dim>(1), 1);
+ DoFHandler<dim> dof_handler(triangulation);
+
+ dof_handler.distribute_dofs(fe);
+
+ std::vector<types::global_dof_index> dofs_per_component (2);
+ DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
+ const unsigned int n_u = dofs_per_component[0],
+ n_p = dofs_per_component[1];
+
+ BlockDynamicSparsityPattern dsp(2, 2);
+ dsp.block(0, 0).reinit (n_u, n_u);
+ dsp.block(1, 0).reinit (n_p, n_u);
+ dsp.block(0, 1).reinit (n_u, n_p);
+ dsp.block(1, 1).reinit (n_p, n_p);
+ dsp.collect_sizes ();
+ DoFTools::make_sparsity_pattern (dof_handler, dsp);
+
+ BlockSparsityPattern sparsity_pattern;
+ sparsity_pattern.copy_from(dsp);
+ sparsity_pattern.compress();
+
+ BlockSparseMatrix<double> a (sparsity_pattern);
+ BlockSparseMatrix<double> b (sparsity_pattern);
+
+ for(unsigned int i = 0; i < a.n(); ++i)
+ {
+ a.set(i, i, 1.);
+ b.set(i, i, 5.);
+ }
+
+ // Constructors and assignment:
+
+ auto op_a = linop<BlockVector<double>>(a);
+ auto op_b = linop<BlockVector<double>>(b);
+
+ {
+ decltype(op_a) op_x (a);
+ op_a = a;
+ op_b = b;
+ }
+
+ // vmult:
+
+ BlockVector<double> u;
+ op_a.reinit_domain_vector(u, true);
+ for (unsigned int i = 0; i < u.size(); ++i) {
+ u[i] = (double)(i+1);
+ }
+
+ PRINTME("u", u);
+
+ BlockVector<double> v;
+ op_a.reinit_domain_vector(v, false);
+ BlockVector<double> w;
+ op_a.reinit_domain_vector(w, false);
+ BlockVector<double> x;
+ op_a.reinit_domain_vector(x, false);
+
+ op_a.vmult(v, u);
+ PRINTME("Au", v);
+
+ op_b.vmult(w, u);
+ PRINTME("Bu", w);
+
+ // operator+, operator-, operator+=, operator-=:
+
+ x = v;
+ x += w;
+ PRINTME("Au+Bu", x);
+
+ (op_a + op_b).vmult(x, u);
+ PRINTME("(A+B)u", x);
+
+ auto op_x = op_a;
+ op_x += op_b;
+ op_x.vmult(x, u);
+ PRINTME("(A+=B)u", x);
+
+ x = v;
+ x -= w;
+ PRINTME("Au-Bu", x);
+
+ (op_a - op_b).vmult(x, u);
+ PRINTME("(A-B)u", x);
+
+ op_x = op_a;
+ op_x -= op_b;
+ op_x.vmult(x, u);
+ PRINTME("(A-=B)u", x);
+
+ // operator*, operator*=
+
+ op_b.vmult(v,u);
+ op_a.vmult(w,v);
+ PRINTME("(A(Bu))", x);
+
+ (op_a * op_b).vmult(x, u);
+ PRINTME("(A*B)u", x);
+
+ op_x = op_a;
+ op_x *= op_b;
+ op_x.vmult(x, u);
+ PRINTME("(A*=B)u", x);
+
+ // solver interface:
+
+ SolverControl solver_control (1000, 1e-10);
+ SolverGMRES<BlockVector<double>> solver (solver_control);
+
+ deallog.depth_file(0);
+ solver.solve(op_b, v, u, PreconditionIdentity());
+ deallog.depth_file(3);
+ PRINTME("solve(B, v, u)", v);
+
+ deallog.depth_file(0);
+ inverse_linop(solver, PreconditionIdentity(), op_b).vmult(v, u);
+ deallog.depth_file(3);
+ PRINTME("inverse_linop(B)u", v);
+
+ deallog.depth_file(0);
+ op_b.vmult(w, v);
+ deallog.depth_file(3);
+ PRINTME("B(inverse_linop(B)u)", w);
+
+ deallog.depth_file(0);
+ (op_b * inverse_linop(solver, PreconditionIdentity(), op_b)).vmult(w, u);
+ deallog.depth_file(3);
+ PRINTME("(B*inverse_linop(B))u", w);
+
+ deallog.depth_file(0);
+ (inverse_linop(solver, PreconditionIdentity(), op_b) * op_b).vmult(w, u);
+ deallog.depth_file(3);
+ PRINTME("(inverse_linop(B)*B)u", w);
+
+ SolverControl inner_solver_control (1000, 1e-12);
+ SolverCG<BlockVector<double>> inner_solver (solver_control);
+
+ deallog.depth_file(0);
+ solver.solve(inverse_linop(inner_solver, PreconditionIdentity(), op_b), v, u,
+ PreconditionIdentity());
+ deallog.depth_file(3);
+ PRINTME("solve(inverse_linop(B), v, u) == Bu", v);
+
+ deallog.depth_file(0);
+ solver.solve(op_b + 0.005 * op_a, v, u, PreconditionIdentity());
+ deallog.depth_file(3);
+ PRINTME("solve(B+0.5*0.01*A, v, u)", v);
+}
--- /dev/null
+
+DEAL::u: [block 0] 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL:: [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000
+DEAL::
+DEAL::Au: [block 0] 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL:: [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000
+DEAL::
+DEAL::Bu: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000
+DEAL:: [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000
+DEAL::
+DEAL::Au+Bu: [block 0] 6.000000000 12.00000000 18.00000000 24.00000000 30.00000000 36.00000000 42.00000000 48.00000000 54.00000000 60.00000000 66.00000000 72.00000000 78.00000000 84.00000000 90.00000000 96.00000000 102.0000000 108.0000000 114.0000000 120.0000000 126.0000000 132.0000000 138.0000000 144.0000000 150.0000000
+DEAL:: [block 1] 156.0000000 162.0000000 168.0000000 174.0000000 180.0000000 186.0000000 192.0000000 198.0000000 204.0000000 210.0000000 216.0000000 222.0000000 228.0000000 234.0000000 240.0000000 246.0000000 252.0000000 258.0000000 264.0000000 270.0000000 276.0000000 282.0000000 288.0000000 294.0000000 300.0000000
+DEAL::
+DEAL::(A+B)u: [block 0] 6.000000000 12.00000000 18.00000000 24.00000000 30.00000000 36.00000000 42.00000000 48.00000000 54.00000000 60.00000000 66.00000000 72.00000000 78.00000000 84.00000000 90.00000000 96.00000000 102.0000000 108.0000000 114.0000000 120.0000000 126.0000000 132.0000000 138.0000000 144.0000000 150.0000000
+DEAL:: [block 1] 156.0000000 162.0000000 168.0000000 174.0000000 180.0000000 186.0000000 192.0000000 198.0000000 204.0000000 210.0000000 216.0000000 222.0000000 228.0000000 234.0000000 240.0000000 246.0000000 252.0000000 258.0000000 264.0000000 270.0000000 276.0000000 282.0000000 288.0000000 294.0000000 300.0000000
+DEAL::
+DEAL::(A+=B)u: [block 0] 6.000000000 12.00000000 18.00000000 24.00000000 30.00000000 36.00000000 42.00000000 48.00000000 54.00000000 60.00000000 66.00000000 72.00000000 78.00000000 84.00000000 90.00000000 96.00000000 102.0000000 108.0000000 114.0000000 120.0000000 126.0000000 132.0000000 138.0000000 144.0000000 150.0000000
+DEAL:: [block 1] 156.0000000 162.0000000 168.0000000 174.0000000 180.0000000 186.0000000 192.0000000 198.0000000 204.0000000 210.0000000 216.0000000 222.0000000 228.0000000 234.0000000 240.0000000 246.0000000 252.0000000 258.0000000 264.0000000 270.0000000 276.0000000 282.0000000 288.0000000 294.0000000 300.0000000
+DEAL::
+DEAL::Au-Bu: [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000
+DEAL:: [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000
+DEAL::
+DEAL::(A-B)u: [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000
+DEAL:: [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000
+DEAL::
+DEAL::(A-=B)u: [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000
+DEAL:: [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000
+DEAL::
+DEAL::(A(Bu)): [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000
+DEAL:: [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000
+DEAL::
+DEAL::(A*B)u: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000
+DEAL:: [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000
+DEAL::
+DEAL::(A*=B)u: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000
+DEAL:: [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000
+DEAL::
+DEAL::solve(B, v, u): [block 0] 0.2000000000 0.4000000000 0.6000000000 0.8000000000 1.000000000 1.200000000 1.400000000 1.600000000 1.800000000 2.000000000 2.200000000 2.400000000 2.600000000 2.800000000 3.000000000 3.200000000 3.400000000 3.600000000 3.800000000 4.000000000 4.200000000 4.400000000 4.600000000 4.800000000 5.000000000
+DEAL:: [block 1] 5.200000000 5.400000000 5.600000000 5.800000000 6.000000000 6.200000000 6.400000000 6.600000000 6.800000000 7.000000000 7.200000000 7.400000000 7.600000000 7.800000000 8.000000000 8.200000000 8.400000000 8.600000000 8.800000000 9.000000000 9.200000000 9.400000000 9.600000000 9.800000000 10.00000000
+DEAL::
+DEAL::inverse_linop(B)u: [block 0] 0.2000000000 0.4000000000 0.6000000000 0.8000000000 1.000000000 1.200000000 1.400000000 1.600000000 1.800000000 2.000000000 2.200000000 2.400000000 2.600000000 2.800000000 3.000000000 3.200000000 3.400000000 3.600000000 3.800000000 4.000000000 4.200000000 4.400000000 4.600000000 4.800000000 5.000000000
+DEAL:: [block 1] 5.200000000 5.400000000 5.600000000 5.800000000 6.000000000 6.200000000 6.400000000 6.600000000 6.800000000 7.000000000 7.200000000 7.400000000 7.600000000 7.800000000 8.000000000 8.200000000 8.400000000 8.600000000 8.800000000 9.000000000 9.200000000 9.400000000 9.600000000 9.800000000 10.00000000
+DEAL::
+DEAL::B(inverse_linop(B)u): [block 0] 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL:: [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000
+DEAL::
+DEAL::(B*inverse_linop(B))u: [block 0] 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL:: [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000
+DEAL::
+DEAL::(inverse_linop(B)*B)u: [block 0] 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000
+DEAL:: [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000
+DEAL::
+DEAL::solve(inverse_linop(B), v, u) == Bu: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000
+DEAL:: [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000
+DEAL::
+DEAL::solve(B+0.5*0.01*A, v, u): [block 0] 0.1998001998 0.3996003996 0.5994005994 0.7992007992 0.9990009990 1.198801199 1.398601399 1.598401598 1.798201798 1.998001998 2.197802198 2.397602398 2.597402597 2.797202797 2.997002997 3.196803197 3.396603397 3.596403596 3.796203796 3.996003996 4.195804196 4.395604396 4.595404595 4.795204795 4.995004995
+DEAL:: [block 1] 5.194805195 5.394605395 5.594405594 5.794205794 5.994005994 6.193806194 6.393606394 6.593406593 6.793206793 6.993006993 7.192807193 7.392607393 7.592407592 7.792207792 7.992007992 8.191808192 8.391608392 8.591408591 8.791208791 8.991008991 9.190809191 9.390609391 9.590409590 9.790209790 9.990009990
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test that it is possible to instantiate a LinearOperator object for all
+// different kinds of Trilinos matrices and vectors
+
+#include "../tests.h"
+
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+using namespace dealii;
+
+int main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ initlog();
+ deallog << std::setprecision(10);
+
+ TrilinosWrappers::SparseMatrix a;
+
+ auto op_a = linop<TrilinosWrappers::MPI::Vector>(a);
+ auto op_a2 = linop<TrilinosWrappers::Vector>(a);
+
+ TrilinosWrappers::BlockSparseMatrix b;
+
+ auto op_b = linop<TrilinosWrappers::MPI::BlockVector>(b);
+ auto op_b3 = linop<TrilinosWrappers::BlockVector>(b);
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::OK
projects / dealii.git / commitdiff