From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Thu, 7 Mar 2024 21:51:35 +0000 (-0700)
Subject: Fix Tpetra SparseMatrix::add().
X-Git-Tag: v9.6.0-rc1~498^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b735a06ff287265c937bdf8ff8da0f6dd5e45af5;p=dealii.git

Fix Tpetra SparseMatrix::add().
---

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 f4ead0de68..1c7a911970 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
@@ -1137,17 +1137,28 @@ namespace LinearAlgebra
       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"));
+               "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);
+             ExcMessage("Addition of matrices is only allowed if the matrices "
+                        "are filled, i.e., if compress() has been called."));
+
+      // Tpetra does not have a function that adds a matrix in-place to the
+      // we're currently storing. But it (inefficiently) has one that returns a
+      // new matrix with the result. As a consequence, replace the current one
+      // by the one that is returned.
+      //
+      // Inconveniently, however, the Tpetra::CrsMatrix::add() function returns
+      // a RCP<Tpetra::RowMatrix>, though it guarantees that the stored object
+      // is actually a Tpetra::CrsMatrix for our combination of arguments. So,
+      // we have to do some casting around to assign things back to 'matrix':
+      auto result = matrix->add(factor,
+                                *source.matrix,
+                                1.0,
+                                matrix->getDomainMap(),
+                                matrix->getRangeMap(),
+                                Teuchos::null);
+      Assert(dynamic_cast<MatrixType *>(result.get()), ExcInternalError());
+      matrix.reset(dynamic_cast<MatrixType *>(result.release().get()));
     }