From: Matthias Maier Date: Wed, 27 May 2015 12:23:37 +0000 (+0200) Subject: Refactoring II: Reintroduce PackagedOperation and related tests X-Git-Tag: v8.3.0-rc1~150^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8a0558ecc71e266756c70477fd8f623bf2251734;p=dealii.git Refactoring II: Reintroduce PackagedOperation and related tests This reintroduces PackagedOperation and related tests under a new header file "packaged_operation.h". Tests renamed accordingly. --- diff --git a/include/deal.II/lac/packaged_operation.h b/include/deal.II/lac/packaged_operation.h new file mode 100644 index 0000000000..13192cce59 --- /dev/null +++ b/include/deal.II/lac/packaged_operation.h @@ -0,0 +1,830 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 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__packaged_operation_h +#define __deal2__packaged_operation_h + +#include +#include +#include + +#ifdef DEAL_II_WITH_CXX11 + +#include + +DEAL_II_NAMESPACE_OPEN + +// Forward declarations: +template class Vector; +template class LinearOperator; +template > class PackagedOperation; + + +/** + * A class to store a computation. + * + * The PackagedOperation class allows lazy evaluation of expressions + * involving vectors and linear operators. This is done by storing the + * computational expression and only performing the computation when either + * the object is implicitly converted to a vector object, or + * apply (or apply_add) is invoked by hand. This + * avoids unnecessary temporary storage of intermediate results. + * + * The class essentially consists of std::function objects + * that store the knowledge of how to generate the result of a computation + * and store it in a vector: + * @code + * std::function apply; + * std::function apply_add; + * @endcode + * + * Similar to the LinearOperator class it also has knowledge about how to + * initialize a vector of the @p Range space: + * @code + * std::function reinit_vector; + * @endcode + * + * As an example consider the addition of multiple vectors + * @code + * dealii::Vector a, b, c, d; + * // .. + * dealii::Vector result = a + b - c + d; + * @endcode + * or the computation of a residual $b-Ax$: + * @code + * dealii::SparseMatrix A; + * dealii::Vector b, x; + * // .. + * const auto op_a = linear_operator(A); + * + * dealii::Vector residual = b - op_a * x; + * @endcode + * The expression residual is of type + * PackagedOperation>. It stores + * references to A, b and x and + * defers the actual computation until apply, or + * apply_add are explicitly invoked, + * @code + * dealii::Vector y; + * residual.reinit_vector(y); + * residual.apply(y); + * residual.apply_add(y); + * @endcode + * or until the @p PackagedOperation object is implictly converted: + * @code + * dealii::Vector y; + * y = residual; + * y += residual; + * y -= residual; + * @endcode + * + * @author Matthias Maier, 2015 + * + * @ingroup LAOperators + */ +template class PackagedOperation +{ +public: + + /** + * Create an empty PackagedOperation object. All std::function + * member objects are initialized with default variants that throw an + * exception upon invocation. + */ + PackagedOperation () + { + apply = [](Range &) + { + Assert(false, + ExcMessage("Uninitialized PackagedOperation::apply called")); + }; + + apply_add = [](Range &) + { + Assert(false, + ExcMessage("Uninitialized PackagedOperation::apply_add called")); + }; + + reinit_vector = [](Range &, bool) + { + Assert(false, + ExcMessage("Uninitialized PackagedOperation::reinit_vector " + "method called")); + }; + } + + /** + * Default copy constructor. + */ + PackagedOperation (const PackagedOperation &) = default; + + /** + * Constructor that creates a PackagedOperation object from a reference vector + * @p u. The PackagedOperation returns @p u. + * + * The PackagedOperation object that is created stores a reference to @p + * u. Thus, the vector must remain a valid reference for the whole + * lifetime of the PackagedOperation object. All changes made on @p u + * after the creation of the PackagedOperation object are reflected by + * the operator object. + */ + PackagedOperation (const Range &u) + { + *this = u; + } + + /** + * Default copy assignment operator. + */ + PackagedOperation &operator=(const PackagedOperation &) = default; + + /** Copy assignment operator that creates a PackagedOperation object from a + * reference vector @p u. The PackagedOperation returns @p u. + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vector must remain a valid reference for the whole lifetime + * of the PackagedOperation object. All changes made on @p u after the creation + * of the PackagedOperation object are reflected by the operator object. + */ + PackagedOperation &operator=(const Range &u) + { + apply = [&u](Range &v) + { + v = u; + }; + + apply_add = [&u](Range &v) + { + v += u; + }; + + reinit_vector = [&u](Range &v, bool fast) + { + v.reinit(u, fast); + }; + + return *this; + } + + /** + * Convert a PackagedOperation to its result. + * + * This conversion operator creates a vector of the Range space and calls + * apply() on it. + */ + operator Range() const + { + Range result_vector; + + reinit_vector(result_vector, /*bool fast=*/ true); + apply(result_vector); + + return result_vector; + } + + /** + * @name In-place vector space operations + */ + //@{ + + /** + * Addition with a PackagedOperation @p second_comp with the same @p Range. + */ + PackagedOperation &operator+=(const PackagedOperation &second_comp) + { + return *this = *this + second_comp; + } + + /** + * Subtraction with a PackagedOperation @p second_comp with the same @p Range. + */ + PackagedOperation &operator-=(const PackagedOperation &second_comp) + { + return *this = *this - second_comp; + } + + /** + * Add a constant @p offset (of the @p Range space) to the result of a + * PackagedOperation. + */ + PackagedOperation &operator+=(const Range &offset) + { + return *this = *this + PackagedOperation(offset); + } + + /** + * Subract a constant @p offset (of the @p Range space) from the result + * of a PackagedOperation. + */ + PackagedOperation &operator-=(const Range &offset) + { + return *this = *this - PackagedOperation(offset); + } + + /** + * Scalar multiplication of the PackagedOperation with a @p number. + */ + PackagedOperation &operator*=(typename Range::value_type number) + { + return *this = *this * number; + } + //@} + + /** + * Store the result of the PackagedOperation in a vector v of the @p Range + * space. + */ + std::function apply; + + /** + * Add the result of the PackagedOperation to a vector v of the @p Range space. + */ + std::function apply_add; + + /** + * Initializes a vector v of the Range space to be directly usable as the + * destination parameter in an application of apply, or apply_add. + * Similar to the reinit functions of the vector classes, the boolean + * determines whether a fast initialization is done, i.e., if it is set + * to false the content of the vector is set to 0. + */ + std::function reinit_vector; +}; + + +/** + * @name Vector space operations + */ +//@{ + +/** + * @relates PackagedOperation + * + * Addition of two PackagedOperation objects @p first_comp and + * @p second_comp given by vector space addition of the corresponding results. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator+(const PackagedOperation &first_comp, + const PackagedOperation &second_comp) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = first_comp.reinit_vector; + + // ensure to have valid PackagedOperation objects by catching first_comp and + // second_comp by value + + return_comp.apply = [first_comp, second_comp](Range &v) + { + first_comp.apply(v); + second_comp.apply_add(v); + }; + + return_comp.apply_add = [first_comp, second_comp](Range &v) + { + first_comp.apply_add(v); + second_comp.apply_add(v); + }; + + return return_comp; +} + +/** + * @relates PackagedOperation + * + * Subtraction of two PackagedOperation objects @p first_comp and @p second_comp + * given by vector space addition of the corresponding results. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator-(const PackagedOperation &first_comp, + const PackagedOperation &second_comp) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = first_comp.reinit_vector; + + // ensure to have valid PackagedOperation objects by catching first_comp and + // second_comp by value + + return_comp.apply = [first_comp, second_comp](Range &v) + { + second_comp.apply(v); + v *= -1.; + first_comp.apply_add(v); + }; + + return_comp.apply_add = [first_comp, second_comp](Range &v) + { + first_comp.apply_add(v); + v *= -1.; + second_comp.apply_add(v); + v *= -1.; + }; + + return return_comp; +} + +/** + * @relates PackagedOperation + * + * Scalar multiplication of a PackagedOperation objects @p comp with a + * scalar @p number given by a scaling PackagedOperation result with + * @p number. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const PackagedOperation &comp, + typename Range::value_type number) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = comp.reinit_vector; + + // the trivial case: number is zero + if (number == 0.) + { + return_comp.apply = [](Range &v) + { + v = 0.; + }; + + return_comp.apply_add = [](Range &) + { + }; + } + else + { + return_comp.apply = [comp, number](Range &v) + { + comp.apply(v); + v *= number; + }; + + return_comp.apply_add = [comp, number](Range &v) + { + v /= number; + comp.apply_add(v); + v *= number; + }; + } + + return return_comp; +} + +/** + * @relates PackagedOperation + * + * Scalar multiplication of a PackagedOperation objects @p comp with a + * scalar @p number given by a scaling PackagedOperation result with + * @p number. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(typename Range::value_type number, + const PackagedOperation &comp) +{ + return comp * number; +} + +/** + * @relates PackagedOperation + * + * Add a constant @p offset (of the @p Range space) to the result of a + * PackagedOperation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator+(const PackagedOperation &comp, + const Range &offset) +{ + return comp + PackagedOperation(offset); +} + +/** + * @relates PackagedOperation + * + * Add a constant @p offset (of the @p Range space) to the result of a + * PackagedOperation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator+(const Range &offset, + const PackagedOperation &comp) +{ + return PackagedOperation(offset) + comp; +} + +/** + * @relates PackagedOperation + * + * Subtract a constant @p offset (of the @p Range space) from the result of a + * PackagedOperation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator-(const PackagedOperation &comp, + const Range &offset) +{ + return comp - PackagedOperation(offset); +} + + +/** + * @relates PackagedOperation + * + * Subtract a computational result from a constant @p offset (of the @p + * Range space). The result is a PackagedOperation object that applies this + * computation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator-(const Range &offset, + const PackagedOperation &comp) +{ + return PackagedOperation(offset) - comp; +} + +//@} + + +/** + * @name Creation of a PackagedOperation object + */ +//@{ + +namespace +{ + // Poor man's trait class that determines whether type T is a vector: + // FIXME: Implement this as a proper type trait - similar to + // isBlockVector + + template + class has_vector_interface + { + template + static std::false_type test(...); + + template + static std::true_type test(decltype(&C::operator+=), + decltype(&C::operator-=), + decltype(&C::l2_norm)); + + public: + // type is std::true_type if Matrix provides vmult_add and Tvmult_add, + // otherwise it is std::false_type + + typedef decltype(test(0, 0, 0)) type; + }; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object that stores the addition of two vectors. + * + * The PackagedOperation object that is created stores a reference to @p u and @p + * v. Thus, the vectors must remain valid references for the whole lifetime + * of the PackagedOperation object. All changes made on @p u or @p v after the + * creation of the PackagedOperation object are reflected by the operator object. + * + * @ingroup LAOperators + */ + +template ::type::value>::type> +PackagedOperation operator+(const Range &u, const Range &v) +{ + PackagedOperation return_comp; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.reinit_vector = [&u](Range &x, bool fast) + { + x.reinit(u, fast); + }; + + return_comp.apply = [&u, &v](Range &x) + { + x = u; + x += v; + }; + + return_comp.apply_add = [&u, &v](Range &x) + { + x += u; + x += v; + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object that stores the subtraction of two vectors. + * + * The PackagedOperation object that is created stores a reference to @p u + * and @p v. Thus, the vectors must remain valid references for the whole + * lifetime of the PackagedOperation object. All changes made on @p u or @p + * v after the creation of the PackagedOperation object are reflected by + * the operator object. + * + * @ingroup LAOperators + */ + +template ::type::value>::type> +PackagedOperation operator-(const Range &u, const Range &v) +{ + PackagedOperation return_comp; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.reinit_vector = [&u](Range &x, bool fast) + { + x.reinit(u, fast); + }; + + return_comp.apply = [&u, &v](Range &x) + { + x = u; + x -= v; + }; + + return_comp.apply_add = [&u, &v](Range &x) + { + x += u; + x -= v; + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object that stores the scaling of a vector + * with a @p number. + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vectors must remain valid references for the whole lifetime of + * the PackagedOperation object. All changes made on @p u or @p v after the + * creation of the PackagedOperation object are reflected by the operator + * object. + * + * @ingroup LAOperators + */ +template ::type::value>::type> +PackagedOperation operator*(const Range &u, + typename Range::value_type number) +{ + return PackagedOperation(u) * number; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object that stores the scaling of a vector + * with a @p number. + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vectors must remain valid references for the whole lifetime of + * the PackagedOperation object. All changes made on @p u or @p v after the + * creation of the PackagedOperation object are reflected by the operator + * object. + * + * @ingroup LAOperators + */ +template ::type::value>::type> +PackagedOperation operator*(typename Range::value_type number, + const Range &u) +{ + return number * PackagedOperation(u); +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object from a LinearOperator and a reference to a + * vector @p u of the Domain space. The object stores the PackagedOperation + * $\text{op} \,u$ (in matrix notation). return + * (return_add) are implemented with vmult(__1,u) + * (vmult_add(__1,u)). + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vector must remain a valid reference for the whole lifetime of + * the PackagedOperation object. All changes made on @p u after the + * creation of the PackagedOperation object are reflected by the operator + * object. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const LinearOperator &op, + const Domain &u) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_range_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.apply = [op, &u](Range &v) + { + op.vmult(v, u); + }; + + return_comp.apply_add = [op, &u](Range &v) + { + op.vmult_add(v, u); + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object from a LinearOperator and a reference + * to a vector @p u of the Range space. The object stores the + * PackagedOperation $\text{op}^T \,u$ (in matrix notation). + * return (return_add) are implemented with + * Tvmult(__1,u) (Tvmult_add(__1,u)). + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vector must remain a valid reference for the whole lifetime of + * the PackagedOperation object. All changes made on @p u after the + * creation of the PackagedOperation object are reflected by the operator + * object. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const Range &u, + const LinearOperator &op) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_domain_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.apply = [op, &u](Domain &v) + { + op.Tvmult(v, u); + }; + + return_comp.apply_add = [op, &u](Domain &v) + { + op.Tvmult_add(v, u); + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Composition of a PackagedOperation object with a LinearOperator. The + * object stores the computation $\text{op} \,comp$ (in matrix notation). + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const LinearOperator &op, + const PackagedOperation &comp) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_range_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.apply = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_domain_vector(*i, /*bool fast =*/ true); + + comp.apply(*i); + op.vmult(v, *i); + + vector_memory.free(i); + }; + + return_comp.apply_add = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/ true); + + comp.apply(*i); + op.vmult_add(v, *i); + + vector_memory.free(i); + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Composition of a PackagedOperation object with a LinearOperator. The + * object stores the computation $\text{op}^T \,comp$ (in matrix notation). + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const PackagedOperation &comp, + const LinearOperator &op) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_domain_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.apply = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/ true); + + comp.apply(*i); + op.Tvmult(v, *i); + + vector_memory.free(i); + }; + + return_comp.apply_add = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/ true); + + comp.apply(*i); + op.Tvmult_add(v, *i); + + vector_memory.free(i); + }; + + return return_comp; +} + +//@} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_CXX11 +#endif diff --git a/tests/lac/packaged_operation_01.cc b/tests/lac/packaged_operation_01.cc new file mode 100644 index 0000000000..9acbc29b52 --- /dev/null +++ b/tests/lac/packaged_operation_01.cc @@ -0,0 +1,84 @@ +// --------------------------------------------------------------------- +// +// 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 PackagedOperation for +// dealii::Vector +// dealii::SparseMatrix + +#include "../tests.h" + +#include +#include + +using namespace dealii; + + +void test_applies(std::string description, + const PackagedOperation > &expr) +{ + // test apply + Vector tmp = expr; + deallog << description << ": " << tmp << std::endl; + + // test apply_add + for (auto &i : tmp) + i = 100.; + expr.apply_add(tmp); + deallog << "100. * 1_n + " << description << ": " << tmp << std::endl; +} + + +int main() +{ + initlog(); + deallog << std::setprecision(10); + + static const int dim = 2; + + // Tests: + + Vector u(25); + for (unsigned int i = 0; i < u.size(); ++i) { + u[i] = (double)(i+1); + } + + deallog << "u: " << u << std::endl; + + // creation via operator+, operator-, operator* + + test_applies("u + u", u + u); + test_applies("u - u", u - u); + test_applies("3. * u", 3. * u); + test_applies("u * 3.", u * 3.); + + // creation via mixed operator+, operator- + + auto expr = 2. * u; + + test_applies("2. * u + u", expr + u); + test_applies("2. * u - u", expr - u); + + test_applies("u + 2. * u", u + expr); + test_applies("u - 2. * u", u - expr); + + // operator+, operator-, operator* + + PackagedOperation > expr2 = u; + + test_applies("2. * u + u", expr + expr2); + test_applies("2. * u - u", expr - expr2); + test_applies("3. * u", 3. * expr2); + test_applies("u * 3.", expr2 * 3.); +} diff --git a/tests/lac/packaged_operation_01.with_cxx11=on.output b/tests/lac/packaged_operation_01.with_cxx11=on.output new file mode 100644 index 0000000000..8efd0139ea --- /dev/null +++ b/tests/lac/packaged_operation_01.with_cxx11=on.output @@ -0,0 +1,51 @@ + +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::u + u: 2.000000000 4.000000000 6.000000000 8.000000000 10.00000000 12.00000000 14.00000000 16.00000000 18.00000000 20.00000000 22.00000000 24.00000000 26.00000000 28.00000000 30.00000000 32.00000000 34.00000000 36.00000000 38.00000000 40.00000000 42.00000000 44.00000000 46.00000000 48.00000000 50.00000000 +DEAL:: +DEAL::100. * 1_n + u + u: 102.0000000 104.0000000 106.0000000 108.0000000 110.0000000 112.0000000 114.0000000 116.0000000 118.0000000 120.0000000 122.0000000 124.0000000 126.0000000 128.0000000 130.0000000 132.0000000 134.0000000 136.0000000 138.0000000 140.0000000 142.0000000 144.0000000 146.0000000 148.0000000 150.0000000 +DEAL:: +DEAL::u - u: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL:: +DEAL::100. * 1_n + u - u: 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 +DEAL:: +DEAL::3. * u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + 3. * u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: +DEAL::u * 3.: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + u * 3.: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: +DEAL::2. * u + u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + 2. * u + u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: +DEAL::2. * u - 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::100. * 1_n + 2. * u - u: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 +DEAL:: +DEAL::u + 2. * u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + u + 2. * u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: +DEAL::u - 2. * 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::100. * 1_n + u - 2. * u: 99.00000000 98.00000000 97.00000000 96.00000000 95.00000000 94.00000000 93.00000000 92.00000000 91.00000000 90.00000000 89.00000000 88.00000000 87.00000000 86.00000000 85.00000000 84.00000000 83.00000000 82.00000000 81.00000000 80.00000000 79.00000000 78.00000000 77.00000000 76.00000000 75.00000000 +DEAL:: +DEAL::2. * u + u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + 2. * u + u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: +DEAL::2. * u - 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::100. * 1_n + 2. * u - u: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 +DEAL:: +DEAL::3. * u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + 3. * u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: +DEAL::u * 3.: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 +DEAL:: +DEAL::100. * 1_n + u * 3.: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 +DEAL:: diff --git a/tests/lac/packaged_operation_02.cc b/tests/lac/packaged_operation_02.cc new file mode 100644 index 0000000000..070fa4e4e9 --- /dev/null +++ b/tests/lac/packaged_operation_02.cc @@ -0,0 +1,115 @@ +// --------------------------------------------------------------------- +// +// 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 PackagedOperation for +// dealii::SparseMatrix with LinearOperator + +#include "../tests.h" + +#include + +// and a _lot_ of stuff to create a linera oprator +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +void test_applies(std::string description, + const PackagedOperation > &expr) +{ + // test apply + Vector tmp = expr; + deallog << description << ": " << tmp << std::endl; + + // test apply_add + for (auto &i : tmp) + i = 100.; + expr.apply_add(tmp); + deallog << "100. * 1_n + " << description << ": " << tmp << std::endl; +} + + +int main() +{ + initlog(); + deallog << std::setprecision(10); + + static const int dim = 2; + + // Create mass marix M, and an iterative inverse MInv: + + 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 m (sparsity_pattern); + + QGauss quadrature(4); + MatrixCreator::create_mass_matrix(mapping_q1, dof_handler, quadrature, m); + + auto M = linear_operator(m); + + SolverControl solver_control (1000, 1e-12); + SolverCG<> solver (solver_control); + + auto MInv = inverse_operator(M, solver, PreconditionIdentity()); + + // Tests: + + Vector u; + M.reinit_domain_vector(u, true); + for (unsigned int i = 0; i < u.size(); ++i) { + u[i] = (double)(i+1); + } + + deallog << "u: " << u << std::endl; + + deallog.depth_file(0); + Vector b = MInv * u; + deallog.depth_file(3); + deallog << "b: " << b << std::endl; + + test_applies("M * b", M * b); + test_applies("b * M", b * M); + + auto expr = b; + test_applies("M * b", M * expr); + test_applies("b * M", expr * M); +} diff --git a/tests/lac/packaged_operation_02.with_cxx11=on.output b/tests/lac/packaged_operation_02.with_cxx11=on.output new file mode 100644 index 0000000000..eb13437136 --- /dev/null +++ b/tests/lac/packaged_operation_02.with_cxx11=on.output @@ -0,0 +1,21 @@ + +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::b: 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::M * b: 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::100. * 1_n + M * b: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 +DEAL:: +DEAL::b * M: 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::100. * 1_n + b * M: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 +DEAL:: +DEAL::M * b: 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::100. * 1_n + M * b: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 +DEAL:: +DEAL::b * M: 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::100. * 1_n + b * M: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 +DEAL::