Add reinit() method to create a LA::TpetraWrappers::SparseMatrix from a dealii::Spars...

diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
index 543ceb8ec68f04d3ee64b826fcc86eae1abedfb0..ad91ed457a65947698d4e6c4b84ec94309900497 100644 (file)
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
@@ -24,6 +24,8 @@
 #  include <deal.II/base/subscriptor.h>
 #  include <deal.II/base/trilinos_utilities.h>
 
+#  include <deal.II/lac/sparse_matrix.h>
+#  include <deal.II/lac/sparsity_pattern.h>
 #  include <deal.II/lac/trilinos_tpetra_sparsity_pattern.h>
 #  include <deal.II/lac/trilinos_tpetra_vector.h>
 
@@ -369,6 +371,32 @@ namespace LinearAlgebra
              const SparsityPatternType &sparsity_pattern,
              const MPI_Comm             communicator  = MPI_COMM_WORLD,
              const bool                 exchange_data = false);
+
+      /**
+       * This function initializes the Trilinos matrix using the deal.II sparse
+       * matrix and the entries stored therein. It uses a threshold to copy only
+       * elements with modulus larger than the threshold (so zeros in the
+       * deal.II matrix can be filtered away). In contrast to the other reinit
+       * function with deal.II sparse matrix argument, this function takes a
+       * %parallel partitioning specified by the user instead of internally
+       * generating it.
+       *
+       * The optional parameter <tt>copy_values</tt> decides whether only the
+       * sparsity structure of the input matrix should be used or the matrix
+       * entries should be copied, too.
+       *
+       * This is a @ref GlossCollectiveOperation "collective operation" that needs to be called on all
+       * processors in order to avoid a dead lock.
+       */
+      void
+      reinit(const IndexSet                     &row_parallel_partitioning,
+             const IndexSet                     &col_parallel_partitioning,
+             const dealii::SparseMatrix<Number> &dealii_sparse_matrix,
+             const MPI_Comm                      communicator = MPI_COMM_WORLD,
+             const double                        drop_tolerance    = 1e-13,
+             const bool                          copy_values       = true,
+             const dealii::SparsityPattern      *use_this_sparsity = nullptr);
+
       /** @} */
 
       /**

diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
index 650c1feb20ff9cdff27fd5156627d0027ca5db9f..48b836e320dd2b132562bc1006a8f5ca7f9388d5 100644 (file)
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
@@ -583,6 +583,108 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    void
+    SparseMatrix<Number, MemorySpace>::reinit(
+      const IndexSet                     &row_parallel_partitioning,
+      const IndexSet                     &col_parallel_partitioning,
+      const dealii::SparseMatrix<Number> &dealii_sparse_matrix,
+      const MPI_Comm                      communicator,
+      const double                        drop_tolerance,
+      const bool                          copy_values,
+      const dealii::SparsityPattern      *use_this_sparsity)
+    {
+      dealii::types::signed_global_dof_index n_rows = dealii_sparse_matrix.m();
+      AssertDimension(row_parallel_partitioning.size(), n_rows);
+      AssertDimension(col_parallel_partitioning.size(),
+                      dealii_sparse_matrix.n());
+
+      const dealii::SparsityPattern &sparsity_pattern =
+        (use_this_sparsity != nullptr) ?
+          *use_this_sparsity :
+          dealii_sparse_matrix.get_sparsity_pattern();
+
+      if (matrix.is_null() || m() != n_rows ||
+          n_nonzero_elements() != sparsity_pattern.n_nonzero_elements() ||
+          copy_values)
+        if (use_this_sparsity == nullptr)
+          reinit(row_parallel_partitioning,
+                 col_parallel_partitioning,
+                 sparsity_pattern,
+                 communicator,
+                 false);
+
+      // in case we do not copy values, we are done
+      if (copy_values == false)
+        return;
+
+        // fill the values: go through all rows of the
+        // matrix, and then all columns. since the sparsity patterns of the
+        // input matrix and the specified sparsity pattern might be different,
+        // need to go through the row for both these sparsity structures
+        // simultaneously in order to really set the correct values.
+#  if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
+      size_type maximum_row_length = matrix->getLocalMaxNumRowEntries();
+#  else
+      size_type maximum_row_length = matrix->getNodeMaxNumRowEntries();
+#  endif
+      std::vector<size_type> row_indices(maximum_row_length);
+      std::vector<Number>    values(maximum_row_length);
+
+      for (dealii::types::signed_global_dof_index row = 0; row < n_rows; ++row)
+        // see if the row is locally stored on this processor
+        if (row_parallel_partitioning.is_element(row) == true)
+          {
+            dealii::SparsityPattern::iterator select_index =
+              sparsity_pattern.begin(row);
+            typename dealii::SparseMatrix<Number>::const_iterator it =
+              dealii_sparse_matrix.begin(row);
+            size_type col = 0;
+            if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols())
+              {
+                // optimized diagonal
+                AssertDimension(it->column(), row);
+                if (std::fabs(it->value()) > drop_tolerance)
+                  {
+                    values[col]        = it->value();
+                    row_indices[col++] = it->column();
+                  }
+                ++select_index;
+                ++it;
+              }
+
+            while (it != dealii_sparse_matrix.end(row) &&
+                   select_index != sparsity_pattern.end(row))
+              {
+                while (select_index->column() < it->column() &&
+                       select_index != sparsity_pattern.end(row))
+                  ++select_index;
+                while (it->column() < select_index->column() &&
+                       it != dealii_sparse_matrix.end(row))
+                  ++it;
+
+                if (it == dealii_sparse_matrix.end(row))
+                  break;
+                if (std::fabs(it->value()) > drop_tolerance)
+                  {
+                    values[col]        = it->value();
+                    row_indices[col++] = it->column();
+                  }
+                ++select_index;
+                ++it;
+              }
+            set(row,
+                col,
+                reinterpret_cast<size_type *>(row_indices.data()),
+                values.data(),
+                false);
+          }
+
+      compress(VectorOperation::insert);
+    }
+
+
+
     // Information on the matrix
 
     template <typename Number, typename MemorySpace>