Tpetra: Enable lac/linear_operator_15

include/deal.II/lac/trilinos_tpetra_sparse_matrix.h

@@ -446,6 +446,146 @@ namespace LinearAlgebra
           const TrilinosScalar *values,
           const bool            elide_zero_values      = true,
           const bool            col_indices_are_sorted = false);
+
+      /**
+       * Set the element (<i>i,j</i>) to @p value.
+       *
+       * This function is able to insert new elements into the matrix as long as
+       * compress() has not been called, so the sparsity pattern will be
+       * extended. When compress() is called for the first time (or in case the
+       * matrix is initialized from a sparsity pattern), no new elements can be
+       * added and an insertion of elements at positions which have not been
+       * initialized will throw an exception.
+       *
+       * For the case that the matrix is constructed without a sparsity pattern
+       * and new matrix entries are added on demand, please note the following
+       * behavior imposed by the underlying Tpetra::CrsMatrix data structure:
+       * If the same matrix entry is inserted more than once, the matrix entries
+       * will be added upon calling compress() (since Tpetra does not track
+       * values to the same entry before the final compress() is called), even
+       * if VectorOperation::insert is specified as argument to compress(). In
+       * the case you cannot make sure that matrix entries are only set once,
+       * initialize the matrix with a sparsity pattern to fix the matrix
+       * structure before inserting elements.
+       */
+      void
+      set(const size_type i, const size_type j, const Number value);
+
+      /**
+       * Set all elements given in a FullMatrix<double> into the sparse matrix
+       * locations given by <tt>indices</tt>. In other words, this function
+       * writes the elements in <tt>full_matrix</tt> into the calling matrix,
+       * using the local-to-global indexing specified by <tt>indices</tt> for
+       * both the rows and the columns of the matrix. This function assumes a
+       * quadratic sparse matrix and a quadratic full_matrix, the usual
+       * situation in FE calculations.
+       *
+       * This function is able to insert new elements into the matrix as long as
+       * compress() has not been called, so the sparsity pattern will be
+       * extended. After compress() has been called for the first time or the
+       * matrix has been initialized from a sparsity pattern, extending the
+       * sparsity pattern is no longer possible and an insertion of elements at
+       * positions which have not been initialized will throw an exception.
+       *
+       * The optional parameter <tt>elide_zero_values</tt> can be used to
+       * specify whether zero values should be inserted anyway or they should be
+       * filtered away. The default value is <tt>false</tt>, i.e., even zero
+       * values are inserted/replaced.
+       *
+       * For the case that the matrix is constructed without a sparsity pattern
+       * and new matrix entries are added on demand, please note the following
+       * behavior imposed by the underlying Tpetra::CrsMatrix data structure:
+       * If the same matrix entry is inserted more than once, the matrix entries
+       * will be added upon calling compress() (since Epetra does not track
+       * values to the same entry before the final compress() is called), even
+       * if VectorOperation::insert is specified as argument to compress(). In
+       * the case you cannot make sure that matrix entries are only set once,
+       * initialize the matrix with a sparsity pattern to fix the matrix
+       * structure before inserting elements.
+       */
+      void
+      set(const std::vector<size_type> &indices,
+          const FullMatrix<Number>     &full_matrix,
+          const bool                    elide_zero_values = false);
+
+      /**
+       * Same function as before, but now including the possibility to use
+       * rectangular full_matrices and different local-to-global indexing on
+       * rows and columns, respectively.
+       */
+      void
+      set(const std::vector<size_type> &row_indices,
+          const std::vector<size_type> &col_indices,
+          const FullMatrix<Number>     &full_matrix,
+          const bool                    elide_zero_values = false);
+
+      /**
+       * Set several elements in the specified row of the matrix with column
+       * indices as given by <tt>col_indices</tt> to the respective value.
+       *
+       * This function is able to insert new elements into the matrix as long as
+       * compress() has not been called, so the sparsity pattern will be
+       * extended. After compress() has been called for the first time or the
+       * matrix has been initialized from a sparsity pattern, extending the
+       * sparsity pattern is no longer possible and an insertion of elements at
+       * positions which have not been initialized will throw an exception.
+       *
+       * The optional parameter <tt>elide_zero_values</tt> can be used to
+       * specify whether zero values should be inserted anyway or they should be
+       * filtered away. The default value is <tt>false</tt>, i.e., even zero
+       * values are inserted/replaced.
+       *
+       * For the case that the matrix is constructed without a sparsity pattern
+       * and new matrix entries are added on demand, please note the following
+       * behavior imposed by the underlying Tpetra::CrsMatrix data structure:
+       * If the same matrix entry is inserted more than once, the matrix entries
+       * will be added upon calling compress() (since Epetra does not track
+       * values to the same entry before the final compress() is called), even
+       * if VectorOperation::insert is specified as argument to compress(). In
+       * the case you cannot make sure that matrix entries are only set once,
+       * initialize the matrix with a sparsity pattern to fix the matrix
+       * structure before inserting elements.
+       */
+      void
+      set(const size_type               row,
+          const std::vector<size_type> &col_indices,
+          const std::vector<Number>    &values,
+          const bool                    elide_zero_values = false);
+
+      /**
+       * Set several elements to values given by <tt>values</tt> in a given row
+       * in columns given by col_indices into the sparse matrix.
+       *
+       * This function is able to insert new elements into the matrix as long as
+       * compress() has not been called, so the sparsity pattern will be
+       * extended. After compress() has been called for the first time or the
+       * matrix has been initialized from a sparsity pattern, extending the
+       * sparsity pattern is no longer possible and an insertion of elements at
+       * positions which have not been initialized will throw an exception.
+       *
+       * The optional parameter <tt>elide_zero_values</tt> can be used to
+       * specify whether zero values should be inserted anyway or they should be
+       * filtered away. The default value is <tt>false</tt>, i.e., even zero
+       * values are inserted/replaced.
+       *
+       * For the case that the matrix is constructed without a sparsity pattern
+       * and new matrix entries are added on demand, please note the following
+       * behavior imposed by the underlying Tpetra::CrsMatrix data structure:
+       * If the same matrix entry is inserted more than once, the matrix entries
+       * will be added upon calling compress() (since Epetra does not track
+       * values to the same entry before the final compress() is called), even
+       * if VectorOperation::insert is specified as argument to compress(). In
+       * the case you cannot make sure that matrix entries are only set once,
+       * initialize the matrix with a sparsity pattern to fix the matrix
+       * structure before inserting elements.
+       */
+      void
+      set(const size_type  row,
+          const size_type  n_cols,
+          const size_type *col_indices,
+          const Number    *values,
+          const bool       elide_zero_values = false);
+
       /** @} */
 
       /**

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

@@ -21,6 +21,7 @@
 #ifdef DEAL_II_TRILINOS_WITH_TPETRA
 
 #  include <deal.II/lac/dynamic_sparsity_pattern.h>
+#  include <deal.II/lac/full_matrix.h>
 #  include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
 #  include <deal.II/lac/trilinos_tpetra_sparsity_pattern.h>
 
@@ -516,6 +517,118 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    void
+    SparseMatrix<Number, MemorySpace>::set(const size_type  row,
+                                           const size_type  n_cols,
+                                           const size_type *col_indices,
+                                           const Number    *values,
+                                           const bool       elide_zero_values)
+    {
+      AssertIndexRange(row, this->m());
+      const types::signed_global_dof_index *col_index_ptr;
+      const Number                         *col_value_ptr;
+      const types::signed_global_dof_index  trilinos_row = row;
+      types::signed_global_dof_index        n_columns;
+
+      boost::container::small_vector<Number, 200> local_value_array(
+        elide_zero_values ? n_cols : 0);
+      boost::container::small_vector<types::signed_global_dof_index, 200>
+        local_index_array(elide_zero_values ? n_cols : 0);
+
+      // If we don't elide zeros, the pointers are already available... need to
+      // cast to non-const pointers as that is the format taken by Trilinos (but
+      // we will not modify const data)
+      if (elide_zero_values == false)
+        {
+          col_index_ptr =
+            reinterpret_cast<const dealii::types::signed_global_dof_index *>(
+              col_indices);
+          col_value_ptr = values;
+          n_columns     = n_cols;
+        }
+      else
+        {
+          // Otherwise, extract nonzero values in each row and get the
+          // respective indices.
+          col_index_ptr = local_index_array.data();
+          col_value_ptr = local_value_array.data();
+
+          n_columns = 0;
+          for (size_type j = 0; j < n_cols; ++j)
+            {
+              const double value = values[j];
+              AssertIsFinite(value);
+              if (value != 0)
+                {
+                  local_index_array[n_columns] = col_indices[j];
+                  local_value_array[n_columns] = value;
+                  ++n_columns;
+                }
+            }
+
+          AssertIndexRange(n_columns, n_cols + 1);
+        }
+
+      // We distinguish between two cases: the first one is when the matrix is
+      // not filled (i.e., it is possible to add new elements to the sparsity
+      // pattern), and the second one is when the pattern is already fixed. In
+      // the former case, we add the possibility to insert new values, and in
+      // the second we just replace data.
+      if (compressed || matrix->isFillComplete())
+        {
+          matrix->resumeFill();
+          compressed = false;
+        }
+
+      if (!matrix->isStaticGraph())
+        matrix->insertGlobalValues(trilinos_row,
+                                   n_columns,
+                                   col_value_ptr,
+                                   col_index_ptr);
+      else
+        matrix->replaceGlobalValues(trilinos_row,
+                                    n_columns,
+                                    col_value_ptr,
+                                    col_index_ptr);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    SparseMatrix<Number, MemorySpace>::set(const size_type i,
+                                           const size_type j,
+                                           const Number    value)
+    {
+      AssertIsFinite(value);
+
+      set(i, 1, &j, &value, false);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    SparseMatrix<Number, MemorySpace>::set(
+      const std::vector<size_type> &indices,
+      const FullMatrix<Number>     &values,
+      const bool                    elide_zero_values)
+    {
+      Assert(indices.size() == values.m(),
+             ExcDimensionMismatch(indices.size(), values.m()));
+      Assert(values.m() == values.n(), ExcNotQuadratic());
+
+      for (size_type i = 0; i < indices.size(); ++i)
+        set(indices[i],
+            indices.size(),
+            indices.data(),
+            &values(i, 0),
+            elide_zero_values);
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     void
     SparseMatrix<Number, MemorySpace>::add(

include/deal.II/lac/trilinos_tpetra_vector.h

@@ -335,7 +335,7 @@ namespace LinearAlgebra
        */
       void
       reinit(const IndexSet &parallel_partitioner,
-             const MPI_Comm  communicator,
+             const MPI_Comm  communicator         = MPI_COMM_WORLD,
              const bool      omit_zeroing_entries = false);
 
       /**

source/lac/trilinos_tpetra_sparse_matrix.cc

@@ -37,6 +37,9 @@ namespace LinearAlgebra
       const MPI_Comm                        communicator,
       const bool                            exchange_data);
 
+    template void
+    SparseMatrix<double>::reinit(const dealii::DynamicSparsityPattern &);
+
   } // namespace TpetraWrappers
 } // namespace LinearAlgebra
 #  endif // DOXYGEN

tests/lac/linear_operator_15.cc

@@ -22,6 +22,10 @@
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/trilinos_sparse_matrix.h>
 #include <deal.II/lac/trilinos_vector.h>
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+#  include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#  include <deal.II/lac/trilinos_tpetra_vector.h>
+#endif
 #include <deal.II/lac/vector.h>
 
 #include "../tests.h"
@@ -91,6 +95,39 @@ main(int argc, char *argv[])
     u = lo_A_plus_id_2 * f;
   }
 
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+  {
+    DynamicSparsityPattern csp(dim, dim);
+    testproblem.five_point_structure(csp);
+
+    LinearAlgebra::TpetraWrappers::SparseMatrix<double> A;
+    A.reinit(csp);
+    testproblem.five_point(A);
+
+    LinearAlgebra::TpetraWrappers::Vector<double> f;
+    f.reinit(complete_index_set(dim));
+    LinearAlgebra::TpetraWrappers::Vector<double> u;
+    u.reinit(complete_index_set(dim));
+
+    A.compress(VectorOperation::insert);
+    f.compress(VectorOperation::insert);
+    u.compress(VectorOperation::insert);
+
+    const auto lo_A =
+      linear_operator<LinearAlgebra::TpetraWrappers::Vector<double>>(A);
+    const auto lo_id_1 =
+      identity_operator<LinearAlgebra::TpetraWrappers::Vector<double>>(
+        lo_A.reinit_range_vector);
+    const auto lo_id_2 =
+      identity_operator<LinearAlgebra::TpetraWrappers::Vector<double>>(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;
+  }
+#endif
+
   deallog << "OK" << std::endl;
 
   return 0;