From 340fdb1c58ef55b329e822c568883b7e0d04b3f4 Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Fri, 2 Feb 2024 11:35:09 +0100 Subject: [PATCH] Make TpetraWrapps::SparseMatrix compatible with BlockMatrixBase. --- .../lac/trilinos_tpetra_sparse_matrix.h | 139 ++++++++++++++++++ .../trilinos_tpetra_sparse_matrix.templates.h | 95 ++++++++++-- 2 files changed, 221 insertions(+), 13 deletions(-) diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h index 47d61170cb..543ceb8ec6 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h @@ -40,6 +40,9 @@ template class SparseMatrix; # ifndef DOXYGEN +template +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 &source); + /** * Add @p value to the element (i,j). * Just as the respective call in deal.II SparseMatrixmatrix scaled by factor to this matrix, i.e. the + * matrix factor*matrix is added to this. 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 &matrix); + /** * Set the element (i,j) 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 &dst, const Vector &src) const; + + /** + * Compute the residual of an equation Mx=b, where the residual is + * defined to be r=b-Mx. Write the residual into @p dst. The + * l2 norm of the residual vector is returned. + * + * Source x and destination dst 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 &dst, + const Vector &x, + const Vector &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>; }; // class SparseMatrix /* ------------------------- Inline functions ---------------------- */ + template + inline void + SparseMatrix::set(const size_type i, + const size_type j, + const Number value) + { + set(i, 1, &j, &value, false); + } + + + template inline void SparseMatrix::add(const size_type i, @@ -960,6 +1065,22 @@ namespace LinearAlgebra + template + inline Number + SparseMatrix::residual( + Vector &dst, + const Vector &x, + const Vector &b) const + { + vmult(dst, x); + dst -= b; + dst *= -1.; + + return dst.l2_norm(); + } + + + template inline dealii::types::signed_global_dof_index SparseMatrix::m() const @@ -991,6 +1112,24 @@ namespace LinearAlgebra + template + inline void + SparseMatrix::prepare_add() + { + // nothing to do here + } + + + + template + inline void + SparseMatrix::prepare_set() + { + // nothing to do here + } + + + template inline const Tpetra::CrsMatrix< Number, 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 1562ab627e..50deeda6de 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h @@ -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 - inline void - SparseMatrix::set(const size_type i, - const size_type j, - const Number value) - { - AssertIsFinite(value); - - set(i, 1, &j, &value, false); - } - - - template inline void SparseMatrix::set( @@ -776,6 +766,84 @@ namespace LinearAlgebra + template + void + SparseMatrix::add( + const Number factor, + const SparseMatrix &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 + void + SparseMatrix::copy_from( + const SparseMatrix &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(*source.matrix, + Teuchos::Copy); + column_space_map = Teuchos::rcp_const_cast(matrix->getColMap()); + compressed = source.compressed; + } + + + + template + void + SparseMatrix::clear() + { + // When we clear the matrix, reset + // the pointer and generate an + // empty matrix. + column_space_map = Utilities::Trilinos::internal::make_rcp( + 0, 0, Utilities::Trilinos::tpetra_comm_self()); + + // Prepare the graph + Teuchos::RCP graph = + Utilities::Trilinos::internal::make_rcp(column_space_map, + column_space_map, + 0); + graph->fillComplete(); + + // Create the matrix from the graph + matrix = Utilities::Trilinos::internal::make_rcp(graph); + + compressed = true; + } + + + // Multiplications template @@ -855,6 +923,7 @@ namespace LinearAlgebra } + template void SparseMatrix::print( -- 2.39.5