From: Sebastian Kinnewig Date: Wed, 31 Jan 2024 08:43:35 +0000 (+0100) Subject: Introduce TpetraWrappers::BlockSparseMatrix X-Git-Tag: relicensing~25^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c567b8451789153cf5e60f3446e6ccfa22bfaf1;p=dealii.git Introduce TpetraWrappers::BlockSparseMatrix --- diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index 0bcc94c9b2..c6d24ec496 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -48,6 +48,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/include/deal.II/lac/trilinos_tpetra_block_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_block_sparse_matrix.h new file mode 100644 index 0000000000..fea9cadf7e --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_block_sparse_matrix.h @@ -0,0 +1,463 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_tpetra_trilinos_block_sparse_matrix_h +#define dealii_tpetra_trilinos_block_sparse_matrix_h + + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# include + +# include +# include +# include +# include +# include +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +// forward declarations +# ifndef DOXYGEN +class BlockSparsityPattern; +template +class BlockSparseMatrix; +# endif + +namespace LinearAlgebra +{ + /** + * @addtogroup TpetraWrappers + * @{ + */ + namespace TpetraWrappers + { + /** + * Blocked sparse matrix based on the + * LinearAlgebra::TpetraWrappers::SparseMatrix class. This class implements + * the functions that are specific to the Trilinos SparseMatrix base objects + * for a blocked sparse matrix, and leaves the actual work relaying most of + * the calls to the individual blocks to the functions implemented in the + * base class. See there also for a description of when this class is + * useful. + * + * In contrast to the deal.II-type SparseMatrix class, the Trilinos matrices + * do not have external objects for the sparsity patterns. Thus, one does + * not determine the size of the individual blocks of a block matrix of this + * type by attaching a block sparsity pattern, but by calling reinit() to + * set the number of blocks and then by setting the size of each block + * separately. In order to fix the data structures of the block matrix, it + * is then necessary to let it know that we have changed the sizes of the + * underlying matrices. For this, one has to call the collect_sizes() + * function, for much the same reason as is documented with the + * BlockSparsityPattern class. + * + * @ingroup Matrix1 @see + * @ref GlossBlockLA "Block (linear algebra)" + */ + template + class BlockSparseMatrix + : public BlockMatrixBase> + { + public: + /** + * Typedef the base class for simpler access to its own alias. + */ + using BaseClass = BlockMatrixBase>; + + /** + * Typedef the type of the underlying matrix. + */ + using BlockType = typename BaseClass::BlockType; + + /** + * Import the alias from the base class. + */ + using value_type = typename BaseClass::value_type; + using pointer = typename BaseClass::pointer; + using const_pointer = typename BaseClass::const_pointer; + using reference = typename BaseClass::reference; + using const_reference = typename BaseClass::const_reference; + using size_type = typename BaseClass::size_type; + using iterator = typename BaseClass::iterator; + using const_iterator = typename BaseClass::const_iterator; + + /** + * Constructor; initializes the matrix to be empty, without any structure, + * i.e. the matrix is not usable at all. This constructor is therefore + * only useful for matrices which are members of a class. All other + * matrices should be created at a point in the data flow where all + * necessary information is available. + * + * You have to initialize the matrix before usage with + * reinit(BlockSparsityPattern). The number of blocks per row and column + * are then determined by that function. + */ + BlockSparseMatrix() = default; + + /** + * Destructor. + */ + ~BlockSparseMatrix() override; + + /** + * Pseudo copy operator only copying empty objects. The sizes of the block + * matrices need to be the same. + */ + BlockSparseMatrix & + operator=(const BlockSparseMatrix &) = default; + + /** + * This operator assigns a scalar to a matrix. Since this does usually not + * make much sense (should we set all matrix entries to this value? Only + * the nonzero entries of the sparsity pattern?), this operation is only + * allowed if the actual value to be assigned is zero. This operator only + * exists to allow for the obvious notation matrix=0, which sets + * all elements of the matrix to zero, but keep the sparsity pattern + * previously used. + */ + BlockSparseMatrix & + operator=(const Number d); + + /** + * Resize the matrix, by setting the number of block rows and columns. + * This deletes all blocks and replaces them with uninitialized ones, i.e. + * ones for which also the sizes are not yet set. You have to do that by + * calling the @p reinit functions of the blocks themselves. Do not forget + * to call collect_sizes() after that on this object. + * + * The reason that you have to set sizes of the blocks yourself is that + * the sizes may be varying, the maximum number of elements per row may be + * varying, etc. It is simpler not to reproduce the interface of the @p + * SparsityPattern class here but rather let the user call whatever + * function they desire. + */ + void + reinit(const size_type n_block_rows, const size_type n_block_columns); + + /** + * Resize the matrix, by using an array of index sets to determine the + * %parallel distribution of the individual matrices. This function + * assumes that a quadratic block matrix is generated. + */ + template + void + reinit(const std::vector &input_maps, + const BlockSparsityPatternType &block_sparsity_pattern, + const MPI_Comm communicator = MPI_COMM_WORLD, + const bool exchange_data = false); + + /** + * Resize the matrix and initialize it by the given sparsity pattern. + * Since no distribution map is given, the result is a block matrix for + * which all elements are stored locally. + */ + template + void + reinit(const BlockSparsityPatternType &block_sparsity_pattern); + + /** + * This function initializes the Trilinos matrix using the deal.II sparse + * matrix and the entries stored therein. It uses a threshold to copy only + * elements whose modulus is larger than the threshold (so zeros in the + * deal.II matrix can be filtered away). + */ + void + reinit( + const std::vector ¶llel_partitioning, + const ::dealii::BlockSparseMatrix &dealii_block_sparse_matrix, + const MPI_Comm communicator = MPI_COMM_WORLD, + const double drop_tolerance = 1e-13); + + /** + * This function initializes the Trilinos matrix using the deal.II sparse + * matrix and the entries stored therein. It uses a threshold to copy only + * elements whose modulus is larger than the threshold (so zeros in the + * deal.II matrix can be filtered away). Since no Epetra_Map is given, all + * the elements will be locally stored. + */ + void + reinit(const ::dealii::BlockSparseMatrix &deal_ii_sparse_matrix, + const double drop_tolerance = 1e-13); + + /** + * Return the state of the matrix, i.e., whether compress() needs to be + * called after an operation requiring data exchange. Does only return + * non-true values when used in debug mode, since it is quite + * expensive to keep track of all operations that lead to the need for + * compress(). + */ + bool + is_compressed() const; + + /** + * This function collects the sizes of the sub-objects and stores them in + * internal arrays, in order to be able to relay global indices into the + * matrix to indices into the subobjects. You *must* call this function + * each time after you have changed the size of the sub-objects. Note that + * this is a @ref GlossCollectiveOperation "collective operation", i.e., + * it needs to be called on all MPI + * processes. This command internally calls the method + * compress(), so you don't need to call that function in case + * you use collect_sizes(). + */ + void + collect_sizes(); + + /** + * Return the total number of nonzero elements of this matrix (summed + * over all MPI processes). + */ + std::uint64_t + n_nonzero_elements() const; + + /** + * Return the underlying MPI communicator. + */ + MPI_Comm + get_mpi_communicator() const; + + /** + * Return the partitioning of the domain space for the individual blocks + * of this matrix, i.e., the partitioning of the block vectors this matrix + * has to be multiplied with. + */ + std::vector + locally_owned_domain_indices() const; + + /** + * Return the partitioning of the range space for the individual blocks of + * this matrix, i.e., the partitioning of the block vectors that result + * from matrix-vector products. + */ + std::vector + locally_owned_range_indices() const; + + /** + * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this + * matrix. The vector types can be block vectors or non-block vectors + * (only if the matrix has only one row or column, respectively), and need + * to define TrilinosWrappers::SparseMatrix::vmult. + */ + template + void + vmult(VectorType1 &dst, const VectorType2 &src) const; + + /** + * Matrix-vector multiplication: let $dst = M^T*src$ with $M$ being this + * matrix. This function does the same as vmult() but takes the transposed + * matrix. + */ + template + void + Tvmult(VectorType1 &dst, const VectorType2 &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. + * + * Note that both vectors have to be distributed vectors generated using + * the same Map as was used for the matrix. + * + * This function only applicable if the matrix only has one block row. + */ + template + Number + residual(VectorType1 &dst, + const VectorType2 &x, + const VectorType3 &b) const; + + /** + * Make the clear() function in the base class visible, though it is + * protected. + */ + using BlockMatrixBase>::clear; + + private: + /** + * Internal version of (T)vmult with two block vectors + */ + template + void + vmult(VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const std::bool_constant, + const std::bool_constant) const; + + /** + * Internal version of (T)vmult where the source vector is a block vector + * but the destination vector is a non-block vector + */ + template + void + vmult(VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const std::bool_constant, + const std::bool_constant) const; + + /** + * Internal version of (T)vmult where the source vector is a non-block + * vector but the destination vector is a block vector + */ + template + void + vmult(VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const std::bool_constant, + const std::bool_constant) const; + + /** + * Internal version of (T)vmult where both source vector and the + * destination vector are non-block vectors (only defined if the matrix + * consists of only one block) + */ + template + void + vmult(VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const std::bool_constant, + const std::bool_constant) const; + }; + + } // namespace TpetraWrappers + + /** @} */ + + // ------------- inline and template functions ----------------- + namespace TpetraWrappers + { + template + template + inline void + BlockSparseMatrix::vmult(VectorType1 &dst, + const VectorType2 &src) const + { + vmult(dst, + src, + false, + std::bool_constant::value>(), + std::bool_constant::value>()); + } + + + + template + template + inline void + BlockSparseMatrix::Tvmult(VectorType1 &dst, + const VectorType2 &src) const + { + vmult(dst, + src, + true, + std::bool_constant::value>(), + std::bool_constant::value>()); + } + + + + template + template + inline void + BlockSparseMatrix::vmult( + VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + std::bool_constant, + std::bool_constant) const + { + if (transpose == true) + BaseClass::Tvmult_block_block(dst, src); + else + BaseClass::vmult_block_block(dst, src); + } + + + + template + template + inline void + BlockSparseMatrix::vmult( + VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + std::bool_constant, + std::bool_constant) const + { + if (transpose == true) + BaseClass::Tvmult_nonblock_block(dst, src); + else + BaseClass::vmult_nonblock_block(dst, src); + } + + + + template + template + inline void + BlockSparseMatrix::vmult( + VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + std::bool_constant, + std::bool_constant) const + { + if (transpose == true) + BaseClass::Tvmult_block_nonblock(dst, src); + else + BaseClass::vmult_block_nonblock(dst, src); + } + + + + template + template + inline void + BlockSparseMatrix::vmult( + VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + std::bool_constant, + std::bool_constant) const + { + if (transpose == true) + BaseClass::Tvmult_nonblock_nonblock(dst, src); + else + BaseClass::vmult_nonblock_nonblock(dst, src); + } + } // namespace TpetraWrappers + +} // namespace LinearAlgebra + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif // dealii_tpetra_trilinos_block_sparse_matrix_h diff --git a/include/deal.II/lac/trilinos_tpetra_block_sparse_matrix.templates.h b/include/deal.II/lac/trilinos_tpetra_block_sparse_matrix.templates.h new file mode 100644 index 0000000000..95305a7e20 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_block_sparse_matrix.templates.h @@ -0,0 +1,340 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_tpetra_trilinos_block_sparse_matrix_templates_h +#define dealii_tpetra_trilinos_block_sparse_matrix_templates_h + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + + namespace TpetraWrappers + { + + template + BlockSparseMatrix::~BlockSparseMatrix() + { + // delete previous content of + // the subobjects array + try + { + clear(); + } + catch (...) + {} + } + + + + template + BlockSparseMatrix & + BlockSparseMatrix::operator=(const Number d) + { + Assert(d == 0, ExcScalarAssignmentOnlyForZeroValue()); + + for (size_type r = 0; r < this->n_block_rows(); ++r) + for (size_type c = 0; c < this->n_block_cols(); ++c) + this->block(r, c) = d; + + return *this; + } + + + + template + void + BlockSparseMatrix::reinit( + const size_type n_block_rows, + const size_type n_block_columns) + { + // first delete previous content of + // the subobjects array + clear(); + + // then resize. set sizes of blocks to + // zero. user will later have to call + // collect_sizes for this + this->sub_objects.reinit(n_block_rows, n_block_columns); + this->row_block_indices.reinit(n_block_rows, 0); + this->column_block_indices.reinit(n_block_columns, 0); + + // and reinitialize the blocks + for (size_type r = 0; r < this->n_block_rows(); ++r) + for (size_type c = 0; c < this->n_block_cols(); ++c) + { + BlockType *p = new BlockType(); + + Assert(this->sub_objects[r][c] == nullptr, ExcInternalError()); + this->sub_objects[r][c] = p; + } + } + + + + template + template + void + BlockSparseMatrix::reinit( + const std::vector ¶llel_partitioning, + const BlockSparsityPatternType &block_sparsity_pattern, + const MPI_Comm communicator, + const bool exchange_data) + { +# ifdef DEBUG + std::vector tpetra_maps; + for (size_type i = 0; i < block_sparsity_pattern.n_block_rows(); ++i) + tpetra_maps.push_back( + parallel_partitioning[i].make_tpetra_map(communicator, false)); + + Assert(tpetra_maps.size() == block_sparsity_pattern.n_block_rows(), + ExcDimensionMismatch(tpetra_maps.size(), + block_sparsity_pattern.n_block_rows())); + Assert(tpetra_maps.size() == block_sparsity_pattern.n_block_cols(), + ExcDimensionMismatch(tpetra_maps.size(), + block_sparsity_pattern.n_block_cols())); + + const size_type n_block_rows = tpetra_maps.size(); + (void)n_block_rows; + + Assert(n_block_rows == block_sparsity_pattern.n_block_rows(), + ExcDimensionMismatch(n_block_rows, + block_sparsity_pattern.n_block_rows())); + Assert(n_block_rows == block_sparsity_pattern.n_block_cols(), + ExcDimensionMismatch(n_block_rows, + block_sparsity_pattern.n_block_cols())); +# endif + + + // Call the other basic reinit function, ... + reinit(block_sparsity_pattern.n_block_rows(), + block_sparsity_pattern.n_block_cols()); + + // ... set the correct sizes, ... + this->row_block_indices = block_sparsity_pattern.get_row_indices(); + this->column_block_indices = block_sparsity_pattern.get_column_indices(); + + // ... and then assign the correct + // data to the blocks. + for (size_type r = 0; r < this->n_block_rows(); ++r) + for (size_type c = 0; c < this->n_block_cols(); ++c) + { + this->sub_objects[r][c]->reinit(parallel_partitioning[r], + parallel_partitioning[c], + block_sparsity_pattern.block(r, c), + communicator, + exchange_data); + } + } + + + + template + template + void + BlockSparseMatrix::reinit( + const BlockSparsityPatternType &block_sparsity_pattern) + { + std::vector parallel_partitioning; + for (size_type i = 0; i < block_sparsity_pattern.n_block_rows(); ++i) + parallel_partitioning.emplace_back( + complete_index_set(block_sparsity_pattern.block(i, 0).n_rows())); + + reinit(parallel_partitioning, block_sparsity_pattern); + } + + + + template + void + BlockSparseMatrix::reinit( + const std::vector ¶llel_partitioning, + const ::dealii::BlockSparseMatrix &dealii_block_sparse_matrix, + const MPI_Comm communicator, + const double drop_tolerance) + { + const size_type n_block_rows = parallel_partitioning.size(); + + Assert(n_block_rows == dealii_block_sparse_matrix.n_block_rows(), + ExcDimensionMismatch(n_block_rows, + dealii_block_sparse_matrix.n_block_rows())); + Assert(n_block_rows == dealii_block_sparse_matrix.n_block_cols(), + ExcDimensionMismatch(n_block_rows, + dealii_block_sparse_matrix.n_block_cols())); + + // Call the other basic reinit function ... + reinit(n_block_rows, n_block_rows); + + // ... and then assign the correct + // data to the blocks. + for (size_type r = 0; r < this->n_block_rows(); ++r) + for (size_type c = 0; c < this->n_block_cols(); ++c) + { + this->sub_objects[r][c]->reinit(parallel_partitioning[r], + parallel_partitioning[c], + dealii_block_sparse_matrix.block(r, + c), + communicator, + drop_tolerance); + } + + collect_sizes(); + } + + + + template + void + BlockSparseMatrix::reinit( + const ::dealii::BlockSparseMatrix &dealii_block_sparse_matrix, + const double drop_tolerance) + { + Assert(dealii_block_sparse_matrix.n_block_rows() == + dealii_block_sparse_matrix.n_block_cols(), + ExcDimensionMismatch(dealii_block_sparse_matrix.n_block_rows(), + dealii_block_sparse_matrix.n_block_cols())); + Assert(dealii_block_sparse_matrix.m() == dealii_block_sparse_matrix.n(), + ExcDimensionMismatch(dealii_block_sparse_matrix.m(), + dealii_block_sparse_matrix.n())); + + std::vector parallel_partitioning; + for (size_type i = 0; i < dealii_block_sparse_matrix.n_block_rows(); ++i) + parallel_partitioning.emplace_back( + complete_index_set(dealii_block_sparse_matrix.block(i, 0).m())); + + reinit(parallel_partitioning, + dealii_block_sparse_matrix, + MPI_COMM_SELF, + drop_tolerance); + } + + + + template + inline bool + BlockSparseMatrix::is_compressed() const + { + bool compressed = true; + for (size_type row = 0; row < this->n_block_rows(); ++row) + for (size_type col = 0; col < this->n_block_cols(); ++col) + if (this->block(row, col).is_compressed() == false) + { + compressed = false; + break; + } + + return compressed; + } + + + + template + void + BlockSparseMatrix::collect_sizes() + { + // simply forward to the (non-public) function of the base class + BaseClass::collect_sizes(); + } + + + + template + std::uint64_t + BlockSparseMatrix::n_nonzero_elements() const + { + std::uint64_t n_nonzero = 0; + for (size_type rows = 0; rows < this->n_block_rows(); ++rows) + for (size_type cols = 0; cols < this->n_block_cols(); ++cols) + n_nonzero += this->block(rows, cols).n_nonzero_elements(); + + return n_nonzero; + } + + + + template + MPI_Comm + BlockSparseMatrix::get_mpi_communicator() const + { + Assert(this->n_block_cols() != 0, ExcNotInitialized()); + Assert(this->n_block_rows() != 0, ExcNotInitialized()); + return this->sub_objects[0][0]->get_mpi_communicator(); + } + + + + template + inline std::vector + BlockSparseMatrix::locally_owned_domain_indices() const + { + Assert(this->n_block_cols() != 0, ExcNotInitialized()); + Assert(this->n_block_rows() != 0, ExcNotInitialized()); + + std::vector domain_indices; + for (size_type c = 0; c < this->n_block_cols(); ++c) + domain_indices.push_back( + this->sub_objects[0][c]->locally_owned_domain_indices()); + + return domain_indices; + } + + + + template + inline std::vector + BlockSparseMatrix::locally_owned_range_indices() const + { + Assert(this->n_block_cols() != 0, ExcNotInitialized()); + Assert(this->n_block_rows() != 0, ExcNotInitialized()); + + std::vector range_indices; + for (size_type r = 0; r < this->n_block_rows(); ++r) + range_indices.push_back( + this->sub_objects[r][0]->locally_owned_range_indices()); + + return range_indices; + } + + + + template + template + Number + BlockSparseMatrix::residual(VectorType1 &dst, + const VectorType2 &x, + const VectorType3 &b) const + { + vmult(dst, x); + dst -= b; + dst *= -1.; + + return dst.l2_norm(); + } + + } // namespace TpetraWrappers + +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif // dealii_tpetra_trilinos_block_sparse_matrix_templates_h diff --git a/include/deal.II/lac/trilinos_tpetra_block_vector.h b/include/deal.II/lac/trilinos_tpetra_block_vector.h index 73da7c105c..b47167fbb2 100644 --- a/include/deal.II/lac/trilinos_tpetra_block_vector.h +++ b/include/deal.II/lac/trilinos_tpetra_block_vector.h @@ -22,6 +22,7 @@ # include # include +# include # include DEAL_II_NAMESPACE_OPEN @@ -39,9 +40,8 @@ namespace LinearAlgebra template class BlockVector; - // TODO - // template - // class BlockSparseMatrix; + template + class BlockSparseMatrix; } // namespace TpetraWrappers } // namespace LinearAlgebra # endif // DOXYGEN diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 624d668cc4..cc9d8da6de 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -144,6 +144,7 @@ if(DEAL_II_WITH_TRILINOS) trilinos_sparse_matrix.cc trilinos_sparsity_pattern.cc trilinos_tpetra_block_vector.cc + trilinos_tpetra_block_sparse_matrix.cc trilinos_tpetra_communication_pattern.cc trilinos_tpetra_solver_direct.cc trilinos_tpetra_sparse_matrix.cc diff --git a/source/lac/affine_constraints.inst.in b/source/lac/affine_constraints.inst.in index e3665527a8..1c232930de 100644 --- a/source/lac/affine_constraints.inst.in +++ b/source/lac/affine_constraints.inst.in @@ -256,8 +256,45 @@ for (S : TRILINOS_SCALARS) LinearAlgebra::TpetraWrappers::Vector &, bool, std::integral_constant) const; - } + // BlockSparseMatrix + template void AffineConstraints::distribute_local_to_global< + LinearAlgebra::TpetraWrappers::BlockSparseMatrix, + LinearAlgebra::TpetraWrappers::Vector>( + const FullMatrix &, + const Vector &, + const std::vector::size_type> &, + LinearAlgebra::TpetraWrappers::BlockSparseMatrix &, + LinearAlgebra::TpetraWrappers::Vector &, + bool, + std::bool_constant) const; + + template void AffineConstraints::distribute_local_to_global< + LinearAlgebra::TpetraWrappers::BlockSparseMatrix, + LinearAlgebra::TpetraWrappers::BlockVector>( + const FullMatrix &, + const Vector &, + const std::vector::size_type> &, + LinearAlgebra::TpetraWrappers::BlockSparseMatrix &, + LinearAlgebra::TpetraWrappers::BlockVector &, + bool, + std::bool_constant) const; + + template void AffineConstraints::distribute_local_to_global< + LinearAlgebra::TpetraWrappers::BlockSparseMatrix>( + const FullMatrix &, + const std::vector::size_type> &, + const std::vector::size_type> &, + LinearAlgebra::TpetraWrappers::BlockSparseMatrix &) const; + + template void AffineConstraints::distribute_local_to_global< + LinearAlgebra::TpetraWrappers::BlockSparseMatrix>( + const FullMatrix &, + const std::vector::size_type> &, + const AffineConstraints &, + const std::vector::size_type> &, + LinearAlgebra::TpetraWrappers::BlockSparseMatrix &) const; + } // --------------------------------------------------------------------- diff --git a/source/lac/trilinos_tpetra_block_sparse_matrix.cc b/source/lac/trilinos_tpetra_block_sparse_matrix.cc new file mode 100644 index 0000000000..fafcebe325 --- /dev/null +++ b/source/lac/trilinos_tpetra_block_sparse_matrix.cc @@ -0,0 +1,73 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# include + +DEAL_II_NAMESPACE_OPEN + +# ifndef DOXYGEN +// explicit instantiations +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + template class BlockSparseMatrix; + + template void + BlockSparseMatrix::reinit( + const ::dealii::BlockDynamicSparsityPattern &); + + template void + BlockSparseMatrix::vmult( + TpetraWrappers::Vector &, + const TpetraWrappers::Vector &) const; + + template void + BlockSparseMatrix::Tvmult( + TpetraWrappers::Vector &, + const TpetraWrappers::Vector &) const; + + template void + BlockSparseMatrix::vmult( + TpetraWrappers::BlockVector &, + const TpetraWrappers::BlockVector &) const; + + template void + BlockSparseMatrix::Tvmult( + TpetraWrappers::BlockVector &, + const TpetraWrappers::BlockVector &) const; + + template void + BlockSparseMatrix::vmult( + ::dealii::BlockVector &, + const ::dealii::BlockVector &) const; + + template void + BlockSparseMatrix::Tvmult( + ::dealii::BlockVector &, + const ::dealii::BlockVector &) const; + + + } // namespace TpetraWrappers +} // namespace LinearAlgebra +# endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA