From 0167b6550bc9e39b8fdb6250e57a6ec10fc87d01 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 8 Apr 2015 12:13:36 +0200 Subject: [PATCH] Implement a LinearOperator class --- include/deal.II/lac/linear_operator.h | 706 ++++++++++++++++++ tests/lac/linear_operator_01.cc | 212 ++++++ .../linear_operator_01.with_cxx11=on.output | 17 + tests/lac/linear_operator_02.cc | 188 +++++ .../linear_operator_02.with_cxx11=on.output | 39 + tests/lac/linear_operator_03.cc | 208 ++++++ .../linear_operator_03.with_cxx11=on.output | 58 ++ tests/lac/linear_operator_04.cc | 48 ++ ...r_04.with_cxx11=on.with_trilinos=on.output | 2 + 9 files changed, 1478 insertions(+) create mode 100644 include/deal.II/lac/linear_operator.h create mode 100644 tests/lac/linear_operator_01.cc create mode 100644 tests/lac/linear_operator_01.with_cxx11=on.output create mode 100644 tests/lac/linear_operator_02.cc create mode 100644 tests/lac/linear_operator_02.with_cxx11=on.output create mode 100644 tests/lac/linear_operator_03.cc create mode 100644 tests/lac/linear_operator_03.with_cxx11=on.output create mode 100644 tests/lac/linear_operator_04.cc create mode 100644 tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h new file mode 100644 index 0000000000..e08f06c17f --- /dev/null +++ b/include/deal.II/lac/linear_operator.h @@ -0,0 +1,706 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include + +#include + +#ifdef DEAL_II_WITH_CXX11 + +DEAL_II_NAMESPACE_OPEN + +// Forward declarations: + +template class Vector; + +template class LinearOperator; +template , + typename Domain = Range, + typename Matrix> +LinearOperator 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 op that stores the knowledge of this + * operation, one can write: + * + * @code + * dealii::SparseMatrix 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 DEAL_II_WITH_CXX11 is enabled during + * cmake configure. + * + * @author: Luca Heltai, Matthias Maier, 2015 + */ +template 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::vmult called")); + }; + + vmult_add = [](Range &, const Domain &) + { + Assert(false, ExcMessage("Uninitialized LinearOperator::vmult_add called")); + }; + + Tvmult = [](Domain &, const Range &) + { + Assert(false, ExcMessage("Uninitialized LinearOperator::Tvmult called")); + }; + + Tvmult_add = [](Domain &, const Range &) + { + Assert(false, ExcMessage("Uninitialized LinearOperator::Tvmult_add called")); + }; + + reinit_range_vector = [](Range &, bool) + { + Assert(false, ExcMessage("Uninitialized LinearOperator::reinit_range_vector method called")); + }; + + reinit_domain_vector = [](Domain &, bool) + { + Assert(false, ExcMessage("Uninitialized LinearOperator::reinit_domain_vector method called")); + }; + } + + /** + * Create a LinearOperator object from an object @p op for which the + * conversion function linop is defined. + */ + template + LinearOperator (const Op &op) + { + *this = linop(op); + } + + /** + * Copy assignment operator for an object @p op for which the conversion + * function linop is defined. + */ + template + LinearOperator &operator=(const Op &op) + { + return *this = linop(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 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 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 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 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 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 reinit_domain_vector; + + /** + * A memory object for vectors of Range space used for intermediate + * storage during computations in vmult, etc. + */ + mutable GrowingVectorMemory range_vector_memory; + + /** + * A memory object for vectors of Domain space used for intermediate + * storage during computations in vmult, etc. + */ + mutable GrowingVectorMemory domain_vector_memory; + + + /** + * Addition with a LinearOperator @p second_op with the same Domain + * and Range. + */ + LinearOperator & + operator+=(const LinearOperator &second_op) + { + *this = *this + second_op; + return *this; + } + + /** + * Subtraction with a LinearOperator @p second_op with the same Domain + * and Range. + */ + LinearOperator & + operator-=(const LinearOperator &second_op) + { + *this = *this - second_op; + return *this; + } + + /** + * Concatenation of the LinearOperator with an endormorphism @p second_op + * on the Domain space. + */ + LinearOperator & + operator*=(const LinearOperator &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 +LinearOperator +operator+(const LinearOperator &first_op, + const LinearOperator &second_op) +{ + LinearOperator 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 +LinearOperator +operator-(const LinearOperator &first_op, + const LinearOperator &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 +LinearOperator +operator*(const LinearOperator &first_op, + const LinearOperator &second_op) +{ + LinearOperator 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 + * operator*= member functions accepting the appropriate + * scalar Number type for rescaling: + * + * @code + * Domain & operator *=(Number); + * Range & operator *=(Number); + * @endcode + */ +template +LinearOperator +operator*(typename Range::value_type number, + const LinearOperator &op) +{ + static_assert( + std::is_convertible::value, + "Range and Domain must have implicitly convertible 'value_type's"); + + LinearOperator 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 + * operator*= member functions for rescaling: + * + * @code + * Domain & operator *=(Number); + * Range & operator *=(Number); + * @endcode + */ +template +LinearOperator +operator*(const LinearOperator &op, + typename Domain::value_type number) +{ + static_assert( + std::is_convertible::value, + "Range and Domain must have implicitly convertible 'value_type's"); + + return number * op; +} + + +/** + * Returns the transpose linear operations of @p op. + */ +template +LinearOperator +transpose_linop(const LinearOperator &op) +{ + LinearOperator 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 std::function object @p + * exemplar as an argument to initialize the reinit_range_vector + * and reinit_domain_vector objects of the LinearOperator object. + */ +template +LinearOperator +identity_linop(const std::function &exemplar) +{ + LinearOperator 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 + * vmult and Tvmult 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 vmult or Tvmult. + */ +template +LinearOperator +inverse_linop(Solver &solver, + const Preconditioner &preconditioner, + const LinearOperator &op) +{ + typedef typename Solver::vector_type Vector; + + LinearOperator 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 +class ReinitRangeFactory +{ +public: + + template + std::function + 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 +class ReinitDomainFactory +{ +public: + + template + std::function + 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 Domain = Range, + typename Matrix> +LinearOperator linop(const Matrix &matrix) +{ + LinearOperator 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().operator()(matrix); + + return_op.reinit_domain_vector = + ReinitDomainFactory().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 diff --git a/tests/lac/linear_operator_01.cc b/tests/lac/linear_operator_01.cc new file mode 100644 index 0000000000..4d601e31c3 --- /dev/null +++ b/tests/lac/linear_operator_01.cc @@ -0,0 +1,212 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +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 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; +} diff --git a/tests/lac/linear_operator_01.with_cxx11=on.output b/tests/lac/linear_operator_01.with_cxx11=on.output new file mode 100644 index 0000000000..95b24c5f4c --- /dev/null +++ b/tests/lac/linear_operator_01.with_cxx11=on.output @@ -0,0 +1,17 @@ + +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 diff --git a/tests/lac/linear_operator_02.cc b/tests/lac/linear_operator_02.cc new file mode 100644 index 0000000000..b2c34a19b8 --- /dev/null +++ b/tests/lac/linear_operator_02.cc @@ -0,0 +1,188 @@ +// --------------------------------------------------------------------- +// +// 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 +// dealii::SparseMatrix +// dealii::FullMatrix + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +int main() +{ + initlog(); + deallog << std::setprecision(10); + + static const int dim = 2; + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global(2); + + MappingQ1 mapping_q1; + FE_Q q1(1); + DoFHandler 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 a (sparsity_pattern); + SparseMatrix b (sparsity_pattern); + + QGauss 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> op_x (a); + op_a = a; + op_b = b; + } + + // vmult: + + Vector 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 v; + op_a.reinit_domain_vector(v, false); + Vector w; + op_a.reinit_domain_vector(w, false); + Vector 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; +} diff --git a/tests/lac/linear_operator_02.with_cxx11=on.output b/tests/lac/linear_operator_02.with_cxx11=on.output new file mode 100644 index 0000000000..f9f16141e6 --- /dev/null +++ b/tests/lac/linear_operator_02.with_cxx11=on.output @@ -0,0 +1,39 @@ + +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:: diff --git a/tests/lac/linear_operator_03.cc b/tests/lac/linear_operator_03.cc new file mode 100644 index 0000000000..85416d844a --- /dev/null +++ b/tests/lac/linear_operator_03.cc @@ -0,0 +1,208 @@ +// --------------------------------------------------------------------- +// +// 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 +// dealii::BlockSparseMatrix + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global(2); + + MappingQ1 mapping_q1; + FESystem fe(FE_Q(1), 1, FE_Q(1), 1); + DoFHandler dof_handler(triangulation); + + dof_handler.distribute_dofs(fe); + + std::vector 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 a (sparsity_pattern); + BlockSparseMatrix 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>(a); + auto op_b = linop>(b); + + { + decltype(op_a) op_x (a); + op_a = a; + op_b = b; + } + + // vmult: + + BlockVector 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 v; + op_a.reinit_domain_vector(v, false); + BlockVector w; + op_a.reinit_domain_vector(w, false); + BlockVector 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> 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> 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); +} diff --git a/tests/lac/linear_operator_03.with_cxx11=on.output b/tests/lac/linear_operator_03.with_cxx11=on.output new file mode 100644 index 0000000000..7d006da0c1 --- /dev/null +++ b/tests/lac/linear_operator_03.with_cxx11=on.output @@ -0,0 +1,58 @@ + +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:: diff --git a/tests/lac/linear_operator_04.cc b/tests/lac/linear_operator_04.cc new file mode 100644 index 0000000000..6f7f12d1a2 --- /dev/null +++ b/tests/lac/linear_operator_04.cc @@ -0,0 +1,48 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include + +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(a); + auto op_a2 = linop(a); + + TrilinosWrappers::BlockSparseMatrix b; + + auto op_b = linop(b); + auto op_b3 = linop(b); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output b/tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5