From: Sebastian Kinnewig Date: Mon, 29 Jan 2024 20:55:26 +0000 (+0100) Subject: Introduce TpetraWrappers::BlockVector X-Git-Tag: relicensing~25^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2630dc65cc6b98257b30658c001d9fe16d2540e5;p=dealii.git Introduce TpetraWrappers::BlockVector --- diff --git a/include/deal.II/lac/trilinos_tpetra_block_vector.h b/include/deal.II/lac/trilinos_tpetra_block_vector.h new file mode 100644 index 0000000000..73da7c105c --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_block_vector.h @@ -0,0 +1,249 @@ +// --------------------------------------------------------------------- +// +// 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_trilinos_tpetra_block_vector_h +#define dealii_trilinos_tpetra_block_vector_h + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# include +# include +# include + +DEAL_II_NAMESPACE_OPEN + +// forward declaration +# ifndef DOXYGEN +template +class BlockVectorBase; + +namespace LinearAlgebra +{ + // forward declaration + namespace TpetraWrappers + { + template + class BlockVector; + + // TODO + // template + // class BlockSparseMatrix; + } // namespace TpetraWrappers +} // namespace LinearAlgebra +# endif // DOXYGEN + +namespace LinearAlgebra +{ + /** + * @addtogroup TpetraWrappers + * @{ + */ + namespace TpetraWrappers + { + /** + * An implementation of block vectors based on the vector class + * implemented in LinearAlgebra::TpetraWrappers. + * While the base class provides for most of the interface, this class + * handles the actual allocation of vectors and provides functions that + * are specific to the underlying vector type. + * + * The model of distribution of data is such that each of the blocks is + * distributed across all MPI processes named in the MPI communicator. + * I.e. we don't just distribute the whole vector, but each component. In + * the constructors and reinit() functions, one therefore not only has to + * specify the sizes of the individual blocks, but also the number of + * elements of each of these blocks to be stored on the local process. + * + * @ingroup Vectors + * @ingroup TrilinosWrappers @see + * @ref GlossBlockLA "Block (linear algebra)" + */ + template + class BlockVector : public dealii::BlockVectorBase< + TpetraWrappers::Vector> + { + public: + /** + * Typedef the base class for simpler access to its own alias. + */ + using BaseClass = + dealii::BlockVectorBase>; + + /** + * Typedef the type of the underlying vector. + */ + 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; + + /** + * Default constructor. Generate an empty vector without any blocks. + */ + BlockVector() = default; + + /** + * Constructor. Generate a block vector with as many blocks as there are + * entries in @p partitioning. Each IndexSet together with the MPI + * communicator contains the layout of the distribution of data among + * the MPI processes. + */ + BlockVector(const std::vector ¶llel_partitioning, + const MPI_Comm communicator = MPI_COMM_WORLD); + + /** + * Creates a BlockVector with ghost elements. See the respective + * reinit() method for more details. @p ghost_values may contain any + * elements in @p parallel_partitioning, they will be ignored. + */ + BlockVector(const std::vector ¶llel_partitioning, + const std::vector &ghost_values, + const MPI_Comm communicator, + const bool vector_writable = false); + + /** + * Copy-Constructor. Set all the properties of the parallel vector to + * those of the given argument and copy the elements. + */ + BlockVector(const BlockVector &v); + + /** + * Creates a block vector consisting of num_blocks components, + * but there is no content in the individual components and the user has + * to fill appropriate data using a reinit of the blocks. + */ + explicit BlockVector(const size_type num_blocks); + + /** + * Destructor. Clears memory + */ + ~BlockVector() override = default; + + /** + * Copy operator: fill all components of the vector that are locally + * stored with the given scalar value. + */ + BlockVector & + operator=(const Number s); + + /** + * Copy operator for arguments of the same type. + */ + BlockVector & + operator=(const BlockVector &v); + + /** + * Reinitialize the BlockVector to contain as many blocks as there are + * index sets given in the input argument, according to the parallel + * distribution of the individual components described in the maps. + * + * If omit_zeroing_entries==false, the vector is filled with + * zeros. + */ + void + reinit(const std::vector ¶llel_partitioning, + const MPI_Comm communicator = MPI_COMM_WORLD, + const bool omit_zeroing_entries = false); + + /** + * Reinit functionality. This function destroys the old vector content + * and generates a new one based on the input partitioning. In addition + * to just specifying one index set as in all the other methods above, + * this method allows to supply an additional set of ghost entries. + * There are two different versions of a vector that can be created. If + * the flag @p vector_writable is set to @p false, the vector only + * allows read access to the joint set of @p parallel_partitioning and + * @p ghost_entries. The effect of the reinit method is then equivalent + * to calling the other reinit method with an index set containing both + * the locally owned entries and the ghost entries. + * + * If the flag @p vector_writable is set to true, this creates an + * alternative storage scheme for ghost elements that allows multiple + * threads to write into the vector (for the other reinit methods, only + * one thread is allowed to write into the ghost entries at a time). + */ + void + reinit(const std::vector &partitioning, + const std::vector &ghost_values, + const MPI_Comm communicator = MPI_COMM_WORLD, + const bool vector_writable = false); + + /** + * Change the dimension to that of the vector V. The same + * applies as for the other reinit() function. + * + * The elements of V are not copied, i.e. this function is the + * same as calling reinit (V.size(), omit_zeroing_entries). + * + * Note that you must call this (or the other reinit() functions) + * function, rather than calling the reinit() functions of an individual + * block, to allow the block vector to update its caches of vector + * sizes. If you call reinit() on one of the blocks, then subsequent + * actions on this object may yield unpredictable results since they may + * be routed to the wrong block. + */ + void + reinit(const BlockVector &V, + const bool omit_zeroing_entries = false); + + /** + * Change the number of blocks to num_blocks. The individual + * blocks will get initialized with zero size, so it is assumed that the + * user resizes the individual blocks by herself in an appropriate way, + * and calls collect_sizes afterwards. + */ + void + reinit(const size_type num_blocks); + + /** + * Return if this Vector contains ghost elements. + * + * @see + * @ref GlossGhostedVector "vectors with ghost elements" + */ + bool + has_ghost_elements() const; + + /** + * Print to a stream. + */ + void + print(std::ostream &out, + const unsigned int precision = 3, + const bool scientific = true, + const bool across = true) const; + }; + } // namespace TpetraWrappers + + /** @} */ + +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif // dealii_trilinos_tpetra_block_vector_h diff --git a/include/deal.II/lac/trilinos_tpetra_block_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_block_vector.templates.h new file mode 100644 index 0000000000..4ed8f65992 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_block_vector.templates.h @@ -0,0 +1,237 @@ +// --------------------------------------------------------------------- +// +// 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_trilinos_tpetra_block_vector_templates_h +#define dealii_trilinos_tpetra_block_vector_templates_h + +#include + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + template + BlockVector::BlockVector( + const std::vector ¶llel_partitioning, + const MPI_Comm communicator) + { + reinit(parallel_partitioning, communicator, false); + } + + + + template + BlockVector::BlockVector( + const std::vector ¶llel_partitioning, + const std::vector &ghost_values, + const MPI_Comm communicator, + const bool vector_writable) + { + reinit(parallel_partitioning, + ghost_values, + communicator, + vector_writable); + } + + + + template + BlockVector::BlockVector( + const BlockVector &v) + : dealii::BlockVectorBase>() + { + this->block_indices = v.block_indices; + + this->components.resize(this->n_blocks()); + for (unsigned int i = 0; i < this->n_blocks(); ++i) + this->components[i] = v.components[i]; + + this->collect_sizes(); + } + + + + template + BlockVector::BlockVector(const size_type num_blocks) + : dealii::BlockVectorBase>() + { + reinit(num_blocks); + } + + + + template + BlockVector & + BlockVector::operator=(const Number s) + { + BaseClass::operator=(s); + return *this; + } + + + + template + BlockVector & + BlockVector::operator=(const BlockVector &v) + { + // we only allow assignment to vectors with the same number of blocks + // or to an empty BlockVector + Assert(this->n_blocks() == 0 || this->n_blocks() == v.n_blocks(), + ExcDimensionMismatch(this->n_blocks(), v.n_blocks())); + + if (this->n_blocks() != v.n_blocks()) + this->block_indices = v.block_indices; + + this->components.resize(this->n_blocks()); + for (unsigned int i = 0; i < this->n_blocks(); ++i) + this->components[i] = v.components[i]; + + this->collect_sizes(); + + return *this; + } + + + + template + void + BlockVector::reinit( + const std::vector ¶llel_partitioning, + const MPI_Comm communicator, + const bool omit_zeroing_entries) + { + // update the number of blocks + this->block_indices.reinit(parallel_partitioning.size(), 0); + + // initialize each block + this->components.resize(this->n_blocks()); + for (unsigned int i = 0; i < this->n_blocks(); ++i) + this->components[i].reinit(parallel_partitioning[i], + communicator, + omit_zeroing_entries); + + // update block_indices content + this->collect_sizes(); + } + + + + template + void + BlockVector::reinit( + const std::vector ¶llel_partitioning, + const std::vector &ghost_values, + const MPI_Comm communicator, + const bool vector_writable) + { + AssertDimension(parallel_partitioning.size(), ghost_values.size()); + + // update the number of blocks + this->block_indices.reinit(parallel_partitioning.size(), 0); + + // initialize each block + this->components.resize(this->n_blocks()); + for (unsigned int i = 0; i < this->n_blocks(); ++i) + this->components[i].reinit(parallel_partitioning[i], + ghost_values[i], + communicator, + vector_writable); + + // update block_indices content + this->collect_sizes(); + } + + + + template + void + BlockVector::reinit( + const BlockVector &v, + const bool omit_zeroing_entries) + { + if (this->n_blocks() != v.n_blocks()) + this->block_indices = v.get_block_indices(); + + this->components.resize(this->n_blocks()); + for (unsigned int i = 0; i < this->n_blocks(); ++i) + this->components[i].reinit(v.components[i], omit_zeroing_entries); + + this->collect_sizes(); + } + + + + template + void + BlockVector::reinit(const size_type num_blocks) + { + this->block_indices.reinit(num_blocks, 0); + + this->components.resize(this->n_blocks()); + for (auto &c : this->components) + c.clear(); + } + + + + template + bool + BlockVector::has_ghost_elements() const + { + bool ghosted = this->block(0).has_ghost_elements(); +# ifdef DEBUG + for (unsigned int i = 0; i < this->n_blocks(); ++i) + Assert(this->block(i).has_ghost_elements() == ghosted, + ExcInternalError()); +# endif + return ghosted; + } + + + + template + void + BlockVector::print(std::ostream &out, + const unsigned int precision, + const bool scientific, + const bool across) const + { + for (unsigned int i = 0; i < this->n_blocks(); ++i) + { + if (across) + out << 'C' << i << ':'; + else + out << "Component " << i << std::endl; + this->components[i].print(out, precision, scientific, across); + } + } + + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif // dealii_trilinos_tpetra_block_vector_templates_h diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index f2548f3924..624d668cc4 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -143,6 +143,7 @@ if(DEAL_II_WITH_TRILINOS) trilinos_solver.cc trilinos_sparse_matrix.cc trilinos_sparsity_pattern.cc + trilinos_tpetra_block_vector.cc trilinos_tpetra_communication_pattern.cc trilinos_tpetra_solver_direct.cc trilinos_tpetra_sparse_matrix.cc diff --git a/source/lac/trilinos_tpetra_block_vector.cc b/source/lac/trilinos_tpetra_block_vector.cc new file mode 100644 index 0000000000..8eccf1971f --- /dev/null +++ b/source/lac/trilinos_tpetra_block_vector.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// 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 + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + // Instantiate these vectors types for specific scalar types. + // + // While there: + // Check that the class we declare here satisfies the + // vector-space-vector concept. If we catch it here, + // any mistake in the vector class declaration would + // show up in uses of this class later on as well. + +# ifdef HAVE_TPETRA_INST_FLOAT +# ifdef DEAL_II_HAVE_CXX20 + static_assert(concepts::is_vector_space_vector>); +# endif + template class BlockVector; +# endif +# ifdef HAVE_TPETRA_INST_DOUBLE +# ifdef DEAL_II_HAVE_CXX20 + static_assert(concepts::is_vector_space_vector>); +# endif + template class BlockVector; +# endif +# ifdef DEAL_II_WITH_COMPLEX_VALUES +# ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT +# ifdef DEAL_II_HAVE_CXX20 + static_assert(concepts::is_vector_space_vector < + BlockVector>); +# endif + template class BlockVector>; +# endif +# ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE +# ifdef DEAL_II_HAVE_CXX20 + static_assert(concepts::is_vector_space_vector < + BlockVector>); +# endif + template class BlockVector>; +# endif +# endif + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA