From: Wolfgang Bangerth 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, 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(result.get()), ExcInternalError()); + matrix.reset(dynamic_cast(result.release().get())); }