]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor identity_operator and null_operator for non-trivial payloads
authorJean-Paul Pelteret <jppelteret@gmail.com>
Wed, 1 Feb 2017 09:52:38 +0000 (10:52 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Wed, 1 Feb 2017 09:52:38 +0000 (10:52 +0100)
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
include/deal.II/lac/trilinos_linear_operator.h
include/deal.II/lac/trilinos_sparse_matrix.h
source/lac/trilinos_sparse_matrix.cc
tests/lac/linear_operator_15.cc [new file with mode: 0644]
tests/lac/linear_operator_15.with_cxx11=on.with_trilinos=on.output [new file with mode: 0644]
tests/lac/linear_operator_16.cc [new file with mode: 0644]
tests/lac/linear_operator_16.with_cxx11=on.with_trilinos=on.output [new file with mode: 0644]

index 4c67f9225d05822347f17e5d25c7ee512a5109de..d0fa73dbfc530a67abf0df4b2f7525da7d9127f3 100644 (file)
@@ -788,6 +788,39 @@ identity_operator(const std::function<void(Range &, bool)> &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 <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+identity_operator(const LinearOperator<Range, Domain, Payload> &op)
+{
+  const LinearOperator<Range, Domain, Payload> id_op
+    = identity_operator<Range,Payload>(op.reinit_range_vector);
+  LinearOperator<Range, Domain, Payload> 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 <typename Range, typename Domain, typename Payload>
 LinearOperator<Range, Domain, Payload>
 null_operator(const LinearOperator<Range, Domain, Payload> &op)
 {
-  auto return_op = op;
+  LinearOperator<Range, Domain, Payload> 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
       */
index 9b788a2e9003d9b74e55b916dd4097aa8712fc88..c5cb9eecc80ba76ab2359c5222243df68c120279 100644 (file)
@@ -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 <typename Range>
-  inline LinearOperator<Range, Range, TrilinosWrappers::internal::LinearOperator::TrilinosPayload>
-  identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
-  {
-    typedef TrilinosWrappers::internal::LinearOperator::TrilinosPayload Payload;
-    return dealii::identity_operator<Range, Payload>(reinit_vector);
-  }
-
-
   /**
    * @relates LinearOperator
    *
index ab74e2f6acc9481ef4cdc21c0ea6e042129362d1..ba1b2923d08908a4838cbe32fa2cea44f53677c2 100644 (file)
@@ -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
         *
index e67368cade0717210c9d005760a45d5061c26d3a..7b598a4dc5486ed43feff460591227ea1b887111 100644 (file)
@@ -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 (file)
index 0000000..34a90a2
--- /dev/null
@@ -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 <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/linear_operator_tools.h>
+
+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<double> A;
+    A.reinit(sp);
+    testproblem.five_point(A);
+
+    Vector<double>       f(dim);
+    Vector<double>       u(dim);
+
+    const auto lo_A = linear_operator<Vector<double> >(A);
+    const auto lo_id = identity_operator<Vector<double> >(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<TrilinosWrappers::Vector>(A);
+    const auto lo_id_1 = identity_operator<TrilinosWrappers::Vector>(lo_A.reinit_range_vector);
+    const auto lo_id_2 = identity_operator<TrilinosWrappers::Vector>(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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -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 (file)
index 0000000..b7ac9c2
--- /dev/null
@@ -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 <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/linear_operator_tools.h>
+
+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<double> A;
+    A.reinit(sp);
+    testproblem.five_point(A);
+
+    Vector<double>       f(dim);
+    Vector<double>       u(dim);
+
+    const auto lo_A = linear_operator<Vector<double> >(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<TrilinosWrappers::Vector>(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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.