]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make TpetraWrapps::SparseMatrix compatible with BlockMatrixBase.
authorSebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Fri, 2 Feb 2024 10:35:09 +0000 (11:35 +0100)
committerSebastian Kinnewig <sebastian@kinnewig.org>
Sun, 18 Feb 2024 15:20:42 +0000 (16:20 +0100)
include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

index 47d61170cbdaf70c4cbb8edc1e003a6bb0706362..543ceb8ec68f04d3ee64b826fcc86eae1abedfb0 100644 (file)
@@ -40,6 +40,9 @@ template <typename Number>
 class SparseMatrix;
 
 #  ifndef DOXYGEN
+template <typename MatrixType>
+class BlockMatrixBase;
+
 namespace LinearAlgebra
 {
   namespace TpetraWrappers
@@ -102,6 +105,22 @@ namespace LinearAlgebra
        */
       using size_type = dealii::types::global_dof_index;
 
+      /**
+       * A structure that describes some of the traits of this class in terms of
+       * its run-time behavior. Some other classes (such as the block matrix
+       * classes) that take one or other of the matrix classes as its template
+       * parameters can tune their behavior based on the variables in this
+       * class.
+       */
+      struct Traits
+      {
+        /**
+         * It is safe to elide additions of zeros to individual elements of this
+         * matrix.
+         */
+        static const bool zero_addition_can_be_elided = true;
+      };
+
       /**
        * Declare an alias for the type used to store matrix elements, in analogy
        * to all the other container classes.
@@ -442,6 +461,12 @@ namespace LinearAlgebra
       SparseMatrix &
       operator/=(const Number factor);
 
+      /**
+       * Copy the given (Trilinos) matrix (sparsity pattern and entries).
+       */
+      void
+      copy_from(const SparseMatrix<Number, MemorySpace> &source);
+
       /**
        * Add @p value to the element (<i>i,j</i>).
        * Just as the respective call in deal.II SparseMatrix<Number,
@@ -475,6 +500,16 @@ namespace LinearAlgebra
           const bool       elide_zero_values      = true,
           const bool       col_indices_are_sorted = false);
 
+      /**
+       * Add <tt>matrix</tt> scaled by <tt>factor</tt> to this matrix, i.e. the
+       * matrix <tt>factor*matrix</tt> is added to <tt>this</tt>. If the
+       * sparsity pattern of the calling matrix does not contain all the
+       * elements in the sparsity pattern of the input matrix, this function
+       * will throw an exception.
+       */
+      void
+      add(const Number factor, const SparseMatrix<Number, MemorySpace> &matrix);
+
       /**
        * Set the element (<i>i,j</i>) to @p value.
        *
@@ -614,6 +649,16 @@ namespace LinearAlgebra
           const Number    *values,
           const bool       elide_zero_values = false);
 
+      /**
+       * Release all memory and return to a state just like after having called
+       * the default constructor.
+       *
+       * This is a @ref GlossCollectiveOperation "collective operation" that needs to be called on all
+       * processors in order to avoid a dead lock.
+       */
+      void
+      clear();
+
       /** @} */
       /**
        * @name Entry Access
@@ -710,6 +755,28 @@ namespace LinearAlgebra
       void
       Tvmult_add(Vector<Number, MemorySpace>       &dst,
                  const Vector<Number, MemorySpace> &src) const;
+
+      /**
+       * Compute the residual of an equation <i>Mx=b</i>, where the residual is
+       * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The
+       * <i>l<sub>2</sub></i> norm of the residual vector is returned.
+       *
+       * Source <i>x</i> and destination <i>dst</i> must not be the same vector.
+       *
+       * The vectors @p dst and @p b have to be initialized with the same
+       * IndexSet that was used for the row indices of the matrix and the vector
+       * @p x has to be initialized with the same IndexSet that was used for the
+       * column indices of the matrix.
+       *
+       * In case of a localized Vector, this function will only work when
+       * running on one processor, since the matrix object is inherently
+       * distributed. Otherwise, an exception will be thrown.
+       */
+      Number
+      residual(Vector<Number, MemorySpace>       &dst,
+               const Vector<Number, MemorySpace> &x,
+               const Vector<Number, MemorySpace> &b) const;
+
       /** @} */
 
       /**
@@ -944,11 +1011,49 @@ namespace LinearAlgebra
        */
       bool compressed;
 
+      /**
+       * For some matrix storage formats, in particular for the PETSc
+       * distributed blockmatrices, set and add operations on individual
+       * elements can not be freely mixed. Rather, one has to synchronize
+       * operations when one wants to switch from setting elements to adding to
+       * elements.  BlockMatrixBase automatically synchronizes the access by
+       * calling this helper function for each block.  This function ensures
+       * that the matrix is in a state that allows adding elements; if it
+       * previously already was in this state, the function does nothing.
+       *
+       * This function is called from BlockMatrixBase.
+       */
+      void
+      prepare_add();
+
+      /**
+       * Same as prepare_add() but prepare the matrix for setting elements if
+       * the representation of elements in this class requires such an
+       * operation.
+       *
+       * This function is called from BlockMatrixBase.
+       */
+      void
+      prepare_set();
+
+      // To allow calling protected prepare_add() and prepare_set().
+      friend class BlockMatrixBase<SparseMatrix<Number, MemorySpace>>;
     }; // class SparseMatrix
 
 
     /* ------------------------- Inline functions ---------------------- */
 
+    template <typename Number, typename MemorySpace>
+    inline void
+    SparseMatrix<Number, MemorySpace>::set(const size_type i,
+                                           const size_type j,
+                                           const Number    value)
+    {
+      set(i, 1, &j, &value, false);
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     inline void
     SparseMatrix<Number, MemorySpace>::add(const size_type i,
@@ -960,6 +1065,22 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    inline Number
+    SparseMatrix<Number, MemorySpace>::residual(
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &x,
+      const Vector<Number, MemorySpace> &b) const
+    {
+      vmult(dst, x);
+      dst -= b;
+      dst *= -1.;
+
+      return dst.l2_norm();
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     inline dealii::types::signed_global_dof_index
     SparseMatrix<Number, MemorySpace>::m() const
@@ -991,6 +1112,24 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    inline void
+    SparseMatrix<Number, MemorySpace>::prepare_add()
+    {
+      // nothing to do here
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    SparseMatrix<Number, MemorySpace>::prepare_set()
+    {
+      // nothing to do here
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     inline const Tpetra::CrsMatrix<
       Number,
index 1562ab627e4ab7b54702cfd0160c49f2efd5af76..50deeda6de6bd45f621fac0dae39101e2a17478a 100644 (file)
@@ -657,6 +657,9 @@ namespace LinearAlgebra
       // 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 the matrix is marked as compressed, we need to
+      // call resumeFill() first.
       if (compressed || matrix->isFillComplete())
         {
           matrix->resumeFill();
@@ -677,19 +680,6 @@ namespace LinearAlgebra
 
 
 
-    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(
@@ -776,6 +766,84 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    void
+    SparseMatrix<Number, MemorySpace>::add(
+      const Number                             factor,
+      const SparseMatrix<Number, MemorySpace> &source)
+    {
+      AssertDimension(source.m(), m());
+      AssertDimension(source.n(), n());
+      AssertDimension(source.local_range().first, local_range().first);
+      AssertDimension(source.local_range().second, local_range().second);
+      Assert(matrix->getRowMap()->isSameAs(*source.matrix->getRowMap()),
+             ExcMessage(
+               "Can only add matrices with same distribution of rows"));
+      Assert(matrix->isFillComplete() && source.matrix->isFillComplete(),
+             ExcMessage("Addition of matrices only allowed if matrices are "
+                        "filled, i.e., compress() has been called"));
+
+      matrix->add(factor,
+                  *source.matrix,
+                  1.0,
+                  matrix->getDomainMap(),
+                  matrix->getRangeMap(),
+                  Teuchos::null);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    SparseMatrix<Number, MemorySpace>::copy_from(
+      const SparseMatrix<Number, MemorySpace> &source)
+    {
+      if (this == &source)
+        return;
+
+      // release memory before reallocation
+      matrix.reset();
+      column_space_map.reset();
+
+      // TODO:
+      // If the source and the target matrix have the same structure, we do
+      // not need to perform a deep copy.
+
+      // Perform a deep copy
+      matrix =
+        Utilities::Trilinos::internal::make_rcp<MatrixType>(*source.matrix,
+                                                            Teuchos::Copy);
+      column_space_map = Teuchos::rcp_const_cast<MapType>(matrix->getColMap());
+      compressed       = source.compressed;
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    SparseMatrix<Number, MemorySpace>::clear()
+    {
+      // When we clear the matrix, reset
+      // the pointer and generate an
+      // empty matrix.
+      column_space_map = Utilities::Trilinos::internal::make_rcp<MapType>(
+        0, 0, Utilities::Trilinos::tpetra_comm_self());
+
+      // Prepare the graph
+      Teuchos::RCP<GraphType> graph =
+        Utilities::Trilinos::internal::make_rcp<GraphType>(column_space_map,
+                                                           column_space_map,
+                                                           0);
+      graph->fillComplete();
+
+      // Create the matrix from the graph
+      matrix = Utilities::Trilinos::internal::make_rcp<MatrixType>(graph);
+
+      compressed = true;
+    }
+
+
+
     // Multiplications
 
     template <typename Number, typename MemorySpace>
@@ -855,6 +923,7 @@ namespace LinearAlgebra
     }
 
 
+
     template <typename Number, typename MemorySpace>
     void
     SparseMatrix<Number, MemorySpace>::print(

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.