Allow LinOp<LA::dist::Vector>(TrilinosMatrix).

diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h
index 015f184fc4eb5752ce212b308ac7151b2d2a865f..7f4374099a4e8e8946ae47b69b68f3c570317982 100644 (file)
--- a/include/deal.II/lac/la_parallel_vector.h
+++ b/include/deal.II/lac/la_parallel_vector.h
@@ -1759,26 +1759,137 @@ namespace internal
     class ReinitHelper<LinearAlgebra::distributed::Vector<Number>>
     {
     public:
-      template <typename Matrix>
+      // A helper type-trait that leverage SFINAE to figure out if type T has
+      // void T::get_mpi_communicator()
+      template <typename T>
+      struct has_get_mpi_communicator
+      {
+      private:
+        static bool
+        detect(...);
+
+        template <typename U>
+        static decltype(std::declval<U>().get_mpi_communicator())
+        detect(const U &);
+
+      public:
+        static const bool value =
+          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+      };
+
+      // A helper type-trait that leverage SFINAE to figure out if type T has
+      // void T::locally_owned_domain_indices()
+      template <typename T>
+      struct has_locally_owned_domain_indices
+      {
+      private:
+        static bool
+        detect(...);
+
+        template <typename U>
+        static decltype(std::declval<U>().locally_owned_domain_indices())
+        detect(const U &);
+
+      public:
+        static const bool value =
+          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+      };
+
+      // A helper type-trait that leverage SFINAE to figure out if type T has
+      // void T::locally_owned_range_indices()
+      template <typename T>
+      struct has_locally_owned_range_indices
+      {
+      private:
+        static bool
+        detect(...);
+
+        template <typename U>
+        static decltype(std::declval<U>().locally_owned_range_indices())
+        detect(const U &);
+
+      public:
+        static const bool value =
+          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+      };
+
+      // A helper type-trait that leverage SFINAE to figure out if type T has
+      // void T::initialize_dof_vector(VectorType v)
+      template <typename T>
+      struct has_initialize_dof_vector
+      {
+      private:
+        static bool
+        detect(...);
+
+        template <typename U>
+        static decltype(std::declval<U>().initialize_dof_vector(
+          LinearAlgebra::distributed::Vector<double>()))
+        detect(const U &);
+
+      public:
+        static const bool value =
+          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+      };
+
+      // Used for (Trilinos/PETSc)Wrappers::SparseMatrix
+      template <typename MatrixType,
+                typename std::enable_if<
+                  has_get_mpi_communicator<MatrixType>::value &&
+                    has_locally_owned_domain_indices<MatrixType>::value,
+                  MatrixType>::type * = nullptr>
       static void
-      reinit_range_vector(const Matrix &                              matrix,
-                          LinearAlgebra::distributed::Vector<Number> &v,
-                          bool omit_zeroing_entries)
+      reinit_domain_vector(MatrixType &                                mat,
+                           LinearAlgebra::distributed::Vector<double> &vec,
+                           bool /*omit_zeroing_entries*/)
       {
-        matrix.initialize_dof_vector(v);
-        if (!omit_zeroing_entries)
-          v = Number();
+        vec.reinit(mat.locally_owned_domain_indices(),
+                   mat.get_mpi_communicator());
       }
 
-      template <typename Matrix>
+      // Used for MatrixFree and DiagonalMatrix
+      template <
+        typename MatrixType,
+        typename std::enable_if<has_initialize_dof_vector<MatrixType>::value,
+                                MatrixType>::type * = nullptr>
       static void
-      reinit_domain_vector(const Matrix &                              matrix,
-                           LinearAlgebra::distributed::Vector<Number> &v,
+      reinit_domain_vector(MatrixType &                                mat,
+                           LinearAlgebra::distributed::Vector<double> &vec,
                            bool omit_zeroing_entries)
       {
-        matrix.initialize_dof_vector(v);
+        mat.initialize_dof_vector(vec);
+        if (!omit_zeroing_entries)
+          vec = Number();
+      }
+
+      // Used for (Trilinos/PETSc)Wrappers::SparseMatrix
+      template <typename MatrixType,
+                typename std::enable_if<
+                  has_get_mpi_communicator<MatrixType>::value &&
+                    has_locally_owned_range_indices<MatrixType>::value,
+                  MatrixType>::type * = nullptr>
+      static void
+      reinit_range_vector(MatrixType &                                mat,
+                          LinearAlgebra::distributed::Vector<double> &vec,
+                          bool /*omit_zeroing_entries*/)
+      {
+        vec.reinit(mat.locally_owned_range_indices(),
+                   mat.get_mpi_communicator());
+      }
+
+      // Used for MatrixFree and DiagonalMatrix
+      template <
+        typename MatrixType,
+        typename std::enable_if<has_initialize_dof_vector<MatrixType>::value,
+                                MatrixType>::type * = nullptr>
+      static void
+      reinit_range_vector(MatrixType &                                mat,
+                          LinearAlgebra::distributed::Vector<double> &vec,
+                          bool omit_zeroing_entries)
+      {
+        mat.initialize_dof_vector(vec);
         if (!omit_zeroing_entries)
-          v = Number();
+          vec = Number();
       }
     };
 

diff --git a/tests/lac/linear_operator_18.cc b/tests/lac/linear_operator_18.cc
new file mode 100644 (file)
index 0000000..00a962f
--- /dev/null
+++ b/tests/lac/linear_operator_18.cc
@@ -0,0 +1,53 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test that it is possible to use a PackagedOperation created by
+// operator*() of a LinearOperator object for Trilinos matrices and vectors
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+
+#include "../tests.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  initlog();
+  deallog << std::setprecision(10);
+
+  typedef LinearAlgebra::distributed::Vector<double> vector_t;
+  typedef TrilinosWrappers::SparseMatrix             matrix_t;
+
+  matrix_t a(5U, 4U, 3U);
+  a.compress(VectorOperation::add);
+
+  auto     op_a = linear_operator<vector_t>(a);
+  vector_t u, res;
+  op_a.reinit_domain_vector(u, false);
+  res = op_a * u;
+
+  auto op_a_transpose = transpose_operator(op_a);
+  op_a.reinit_range_vector(u, false);
+  res = op_a_transpose * u;
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}

diff --git a/tests/lac/linear_operator_18.with_trilinos=on.output b/tests/lac/linear_operator_18.with_trilinos=on.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
+++ b/tests/lac/linear_operator_18.with_trilinos=on.output
@@ -0,0 +1,2 @@
+
+DEAL::OK