From ca103e81cd6956debce03641324042329dbf3811 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 1 Feb 2017 10:52:38 +0100 Subject: [PATCH] Refactor identity_operator and null_operator for non-trivial payloads When performing compound operations with LinearOperators, it is necessary for the Payload type of all operations to be the same. identity_operator did not take in a operator, so we used to define a variant of it in the TrilinosWrappers namespace in order to set up the underlying Payload. This lead to some compiler confusion as it could not deduce which version of the function to call unless the namespace was used to explicitly remove this ambiguity. We've now removed the second identity_operator in TrilinosWrappers, and have instead defined an alternative identity_operator that takes in an exemplar LinearOperator. This version of the function will automatically extract the range of operations from the exemplar, deduce the Payload type AND initialise it correctly (this would not have been done correctly before). This required the introduction of the inverse_payload function to the Payload classes. In anticipation of similar problems with null_operator, this LinearOperator is also slightly refactored as to ensure that the payload is correctly configured. This required the introduction of the null_payload function to the Payload classes. Fixes #3846 --- include/deal.II/lac/linear_operator.h | 57 +++++++++++- .../deal.II/lac/trilinos_linear_operator.h | 22 ----- include/deal.II/lac/trilinos_sparse_matrix.h | 11 ++- source/lac/trilinos_sparse_matrix.cc | 70 ++++++++++++++ tests/lac/linear_operator_15.cc | 92 +++++++++++++++++++ ...r_15.with_cxx11=on.with_trilinos=on.output | 2 + tests/lac/linear_operator_16.cc | 89 ++++++++++++++++++ ...r_16.with_cxx11=on.with_trilinos=on.output | 2 + 8 files changed, 321 insertions(+), 24 deletions(-) create mode 100644 tests/lac/linear_operator_15.cc create mode 100644 tests/lac/linear_operator_15.with_cxx11=on.with_trilinos=on.output create mode 100644 tests/lac/linear_operator_16.cc create mode 100644 tests/lac/linear_operator_16.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 index 4c67f9225d..d0fa73dbfc 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -788,6 +788,39 @@ identity_operator(const std::function &reinit_vector) } +/** + * @relates LinearOperator + * + * Return a LinearOperator that is the identity of the vector space @p Range. + * + * The function takes a LinearOperator @p op and uses its range initializer + * to create an identiy operator. In contrast to the function above, this + * function also ensures that the underlying Payload matches that of the + * input @p op. + * + * @ingroup LAOperators + */ +template +LinearOperator +identity_operator(const LinearOperator &op) +{ + const LinearOperator id_op + = identity_operator(op.reinit_range_vector); + LinearOperator return_op ( + op.identity_payload() + ); + + return_op.reinit_range_vector = id_op.reinit_range_vector; + return_op.reinit_domain_vector = id_op.reinit_domain_vector; + return_op.vmult = id_op.vmult; + return_op.vmult_add = id_op.vmult_add; + return_op.Tvmult = id_op.Tvmult; + return_op.Tvmult_add = id_op.Tvmult_add; + + return return_op; +} + + /** * @relates LinearOperator * @@ -801,7 +834,9 @@ template LinearOperator null_operator(const LinearOperator &op) { - auto return_op = op; + LinearOperator return_op ( + op.null_payload() + ); return_op.is_null_operator = true; @@ -910,6 +945,26 @@ namespace internal { } + /** + * Returns a payload configured for identity operations + */ + EmptyPayload + identity_payload () const + { + return *this; + } + + + /** + * Returns a payload configured for null operations + */ + EmptyPayload + null_payload () const + { + return *this; + } + + /** * Returns a payload configured for transpose operations */ diff --git a/include/deal.II/lac/trilinos_linear_operator.h b/include/deal.II/lac/trilinos_linear_operator.h index 9b788a2e90..c5cb9eecc8 100644 --- a/include/deal.II/lac/trilinos_linear_operator.h +++ b/include/deal.II/lac/trilinos_linear_operator.h @@ -54,28 +54,6 @@ namespace TrilinosWrappers //@{ - /** - * @relates LinearOperator - * - * Return a LinearOperator that is the identity of the vector space @p Range. - * - * This function is the equivalent of the dealii::identity_operator, but - * ensures full compatibility with Trilinos operations by preselecting the - * appropriate template parameters. - * - * @author Jean-Paul Pelteret, 2016 - * - * @ingroup TrilinosWrappers - */ - template - inline LinearOperator - identity_operator(const std::function &reinit_vector) - { - typedef TrilinosWrappers::internal::LinearOperator::TrilinosPayload Payload; - return dealii::identity_operator(reinit_vector); - } - - /** * @relates LinearOperator * diff --git a/include/deal.II/lac/trilinos_sparse_matrix.h b/include/deal.II/lac/trilinos_sparse_matrix.h index ab74e2f6ac..ba1b2923d0 100644 --- a/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_sparse_matrix.h @@ -2200,12 +2200,21 @@ namespace TrilinosWrappers */ virtual ~TrilinosPayload() {} + /** + * Returns a payload configured for identity operations + */ + TrilinosPayload identity_payload () const; + + /** + * Returns a payload configured for null operations + */ + TrilinosPayload null_payload () const; + /** * Returns a payload configured for transpose operations */ TrilinosPayload transpose_payload () const; - /** * Returns a payload configured for inverse operations * diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc index e67368cade..7b598a4dc5 100644 --- a/source/lac/trilinos_sparse_matrix.cc +++ b/source/lac/trilinos_sparse_matrix.cc @@ -2763,6 +2763,76 @@ namespace TrilinosWrappers + TrilinosPayload + TrilinosPayload::identity_payload () const + { + TrilinosPayload return_op (*this); + + return_op.vmult = [](Range &tril_dst, + const Range &tril_src) + { + tril_dst = tril_src; + }; + + return_op.Tvmult = [](Range &tril_dst, + const Range &tril_src) + { + tril_dst = tril_src; + }; + + return_op.inv_vmult = [](Range &tril_dst, + const Range &tril_src) + { + tril_dst = tril_src; + }; + + return_op.inv_Tvmult = [](Range &tril_dst, + const Range &tril_src) + { + tril_dst = tril_src; + }; + + return return_op; + } + + + + TrilinosPayload + TrilinosPayload::null_payload () const + { + TrilinosPayload return_op (*this); + + return_op.vmult = [](Range &tril_dst, + const Domain &) + { + tril_dst.Scale(0); + }; + + return_op.Tvmult = [](Domain &tril_dst, + const Range &) + { + tril_dst.Scale(0); + }; + + return_op.inv_vmult = [](Domain &tril_dst, + const Range &) + { + AssertThrow(false, ExcMessage("Cannot compute inverse of null operator")); + tril_dst.Scale(0); + }; + + return_op.inv_Tvmult = [](Range &tril_dst, + const Domain &) + { + AssertThrow(false, ExcMessage("Cannot compute inverse of null operator")); + tril_dst.Scale(0); + }; + + return return_op; + } + + + TrilinosPayload TrilinosPayload::transpose_payload () const { diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_15.cc new file mode 100644 index 0000000000..34a90a2916 --- /dev/null +++ b/tests/lac/linear_operator_15.cc @@ -0,0 +1,92 @@ + +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2017 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 identity_operator LinearOperator and its payload initialization + + +#include "../tests.h" +#include "../testmatrix.h" + +#include +#include +#include +#include +#include + +using namespace dealii; + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + initlog(); + + // Need to get the linear operator to do some nonsense work to verify + // that it can indeed be used as expected + const unsigned int size = 32; + unsigned int dim = (size-1)*(size-1); + + // Make matrix + FDMatrix testproblem(size, size); + + { + SparsityPattern sp (dim, dim); + testproblem.five_point_structure(sp); + sp.compress(); + + SparseMatrix A; + A.reinit(sp); + testproblem.five_point(A); + + Vector f(dim); + Vector u(dim); + + const auto lo_A = linear_operator >(A); + const auto lo_id = identity_operator >(lo_A.reinit_range_vector); + const auto lo_A_plus_id = lo_A + lo_id; + + u = lo_A_plus_id * f; + } + + { + DynamicSparsityPattern csp (dim, dim); + testproblem.five_point_structure(csp); + + TrilinosWrappers::SparseMatrix A; + A.reinit(csp); + testproblem.five_point(A); + + TrilinosWrappers::Vector f(dim); + TrilinosWrappers::Vector u(dim); + + A.compress (VectorOperation::insert); + f.compress (VectorOperation::insert); + u.compress (VectorOperation::insert); + + const auto lo_A = linear_operator(A); + const auto lo_id_1 = identity_operator(lo_A.reinit_range_vector); + const auto lo_id_2 = identity_operator(lo_A); + const auto lo_A_plus_id_1 = lo_A + lo_id_1; // Not a good idea. See below. + const auto lo_A_plus_id_2 = lo_A + lo_id_2; + + // u = lo_A_plus_id_1 * f; // This is not allowed -- different payloads + u = lo_A_plus_id_2 * f; + } + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/lac/linear_operator_15.with_cxx11=on.with_trilinos=on.output b/tests/lac/linear_operator_15.with_cxx11=on.with_trilinos=on.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/linear_operator_15.with_cxx11=on.with_trilinos=on.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/lac/linear_operator_16.cc b/tests/lac/linear_operator_16.cc new file mode 100644 index 0000000000..b7ac9c2d11 --- /dev/null +++ b/tests/lac/linear_operator_16.cc @@ -0,0 +1,89 @@ + +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2017 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 null_operator LinearOperator and its payload initialization + + +#include "../tests.h" +#include "../testmatrix.h" + +#include +#include +#include +#include +#include + +using namespace dealii; + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + initlog(); + + // Need to get the linear operator to do some nonsense work to verify + // that it can indeed be used as expected + const unsigned int size = 32; + unsigned int dim = (size-1)*(size-1); + + // Make matrix + FDMatrix testproblem(size, size); + + { + SparsityPattern sp (dim, dim); + testproblem.five_point_structure(sp); + sp.compress(); + + SparseMatrix A; + A.reinit(sp); + testproblem.five_point(A); + + Vector f(dim); + Vector u(dim); + + const auto lo_A = linear_operator >(A); + const auto lo_null = null_operator(lo_A); + const auto lo_A_plus_null = lo_A + lo_null; + + u = lo_A_plus_null * f; + } + + { + DynamicSparsityPattern csp (dim, dim); + testproblem.five_point_structure(csp); + + TrilinosWrappers::SparseMatrix A; + A.reinit(csp); + testproblem.five_point(A); + + TrilinosWrappers::Vector f(dim); + TrilinosWrappers::Vector u(dim); + + A.compress (VectorOperation::insert); + f.compress (VectorOperation::insert); + u.compress (VectorOperation::insert); + + const auto lo_A = linear_operator(A); + const auto lo_null = null_operator(lo_A); + const auto lo_A_plus_null = lo_A + lo_null; + + u = lo_A_plus_null * f; + } + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/lac/linear_operator_16.with_cxx11=on.with_trilinos=on.output b/tests/lac/linear_operator_16.with_cxx11=on.with_trilinos=on.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/linear_operator_16.with_cxx11=on.with_trilinos=on.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5