#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_tpetra_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_tpetra_block_vector.h>
#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
#include <deal.II/lac/trilinos_tpetra_vector.h>
#include <deal.II/lac/trilinos_vector.h>
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+
+# include <deal.II/base/template_constraints.h>
+
+# include <deal.II/lac/block_matrix_base.h>
+# include <deal.II/lac/block_sparse_matrix.h>
+# include <deal.II/lac/exceptions.h>
+# include <deal.II/lac/full_matrix.h>
+# include <deal.II/lac/trilinos_tpetra_block_vector.h>
+# include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+# include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+// forward declarations
+# ifndef DOXYGEN
+class BlockSparsityPattern;
+template <typename number>
+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 <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+ class BlockSparseMatrix
+ : public BlockMatrixBase<SparseMatrix<Number, MemorySpace>>
+ {
+ public:
+ /**
+ * Typedef the base class for simpler access to its own alias.
+ */
+ using BaseClass = BlockMatrixBase<SparseMatrix<Number, MemorySpace>>;
+
+ /**
+ * 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<Number, MemorySpace> &
+ operator=(const BlockSparseMatrix<Number, MemorySpace> &) = 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 <tt>matrix=0</tt>, which sets
+ * all elements of the matrix to zero, but keep the sparsity pattern
+ * previously used.
+ */
+ BlockSparseMatrix<Number, MemorySpace> &
+ 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 <typename BlockSparsityPatternType>
+ void
+ reinit(const std::vector<IndexSet> &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 <typename BlockSparsityPatternType>
+ 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<IndexSet> ¶llel_partitioning,
+ const ::dealii::BlockSparseMatrix<double> &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<double> &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 <tt>debug</tt> 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
+ * <tt>compress()</tt>, so you don't need to call that function in case
+ * you use <tt>collect_sizes()</tt>.
+ */
+ 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<IndexSet>
+ 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<IndexSet>
+ 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 <typename VectorType1, typename VectorType2>
+ 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 <typename VectorType1, typename VectorType2>
+ void
+ Tvmult(VectorType1 &dst, const VectorType2 &src) const;
+
+ /**
+ * Compute the residual of an equation <i>Mx=b</i>, where the residual is
+ * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The
+ * <i>l<sub>2</sub></i> norm of the residual vector is returned.
+ *
+ * Source <i>x</i> and destination <i>dst</i> 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 <typename VectorType1,
+ typename VectorType2,
+ typename VectorType3>
+ 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<SparseMatrix<Number, MemorySpace>>::clear;
+
+ private:
+ /**
+ * Internal version of (T)vmult with two block vectors
+ */
+ template <typename VectorType1, typename VectorType2>
+ void
+ vmult(VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const std::bool_constant<true>,
+ const std::bool_constant<true>) const;
+
+ /**
+ * Internal version of (T)vmult where the source vector is a block vector
+ * but the destination vector is a non-block vector
+ */
+ template <typename VectorType1, typename VectorType2>
+ void
+ vmult(VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const std::bool_constant<false>,
+ const std::bool_constant<true>) const;
+
+ /**
+ * Internal version of (T)vmult where the source vector is a non-block
+ * vector but the destination vector is a block vector
+ */
+ template <typename VectorType1, typename VectorType2>
+ void
+ vmult(VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const std::bool_constant<true>,
+ const std::bool_constant<false>) 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 <typename VectorType1, typename VectorType2>
+ void
+ vmult(VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const std::bool_constant<false>,
+ const std::bool_constant<false>) const;
+ };
+
+ } // namespace TpetraWrappers
+
+ /** @} */
+
+ // ------------- inline and template functions -----------------
+ namespace TpetraWrappers
+ {
+ template <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2>
+ inline void
+ BlockSparseMatrix<Number, MemorySpace>::vmult(VectorType1 &dst,
+ const VectorType2 &src) const
+ {
+ vmult(dst,
+ src,
+ false,
+ std::bool_constant<IsBlockVector<VectorType1>::value>(),
+ std::bool_constant<IsBlockVector<VectorType2>::value>());
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2>
+ inline void
+ BlockSparseMatrix<Number, MemorySpace>::Tvmult(VectorType1 &dst,
+ const VectorType2 &src) const
+ {
+ vmult(dst,
+ src,
+ true,
+ std::bool_constant<IsBlockVector<VectorType1>::value>(),
+ std::bool_constant<IsBlockVector<VectorType2>::value>());
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2>
+ inline void
+ BlockSparseMatrix<Number, MemorySpace>::vmult(
+ VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ std::bool_constant<true>,
+ std::bool_constant<true>) const
+ {
+ if (transpose == true)
+ BaseClass::Tvmult_block_block(dst, src);
+ else
+ BaseClass::vmult_block_block(dst, src);
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2>
+ inline void
+ BlockSparseMatrix<Number, MemorySpace>::vmult(
+ VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ std::bool_constant<false>,
+ std::bool_constant<true>) const
+ {
+ if (transpose == true)
+ BaseClass::Tvmult_nonblock_block(dst, src);
+ else
+ BaseClass::vmult_nonblock_block(dst, src);
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2>
+ inline void
+ BlockSparseMatrix<Number, MemorySpace>::vmult(
+ VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ std::bool_constant<true>,
+ std::bool_constant<false>) const
+ {
+ if (transpose == true)
+ BaseClass::Tvmult_block_nonblock(dst, src);
+ else
+ BaseClass::vmult_block_nonblock(dst, src);
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2>
+ inline void
+ BlockSparseMatrix<Number, MemorySpace>::vmult(
+ VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ std::bool_constant<false>,
+ std::bool_constant<false>) 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+
+# include <deal.II/lac/trilinos_tpetra_block_sparse_matrix.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+
+ namespace TpetraWrappers
+ {
+
+ template <typename Number, typename MemorySpace>
+ BlockSparseMatrix<Number, MemorySpace>::~BlockSparseMatrix()
+ {
+ // delete previous content of
+ // the subobjects array
+ try
+ {
+ clear();
+ }
+ catch (...)
+ {}
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ BlockSparseMatrix<Number, MemorySpace> &
+ BlockSparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+ void
+ BlockSparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+ template <typename BlockSparsityPatternType>
+ void
+ BlockSparseMatrix<Number, MemorySpace>::reinit(
+ const std::vector<IndexSet> ¶llel_partitioning,
+ const BlockSparsityPatternType &block_sparsity_pattern,
+ const MPI_Comm communicator,
+ const bool exchange_data)
+ {
+# ifdef DEBUG
+ std::vector<typename BlockType::MapType> 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 <typename Number, typename MemorySpace>
+ template <typename BlockSparsityPatternType>
+ void
+ BlockSparseMatrix<Number, MemorySpace>::reinit(
+ const BlockSparsityPatternType &block_sparsity_pattern)
+ {
+ std::vector<IndexSet> 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 <typename Number, typename MemorySpace>
+ void
+ BlockSparseMatrix<Number, MemorySpace>::reinit(
+ const std::vector<IndexSet> ¶llel_partitioning,
+ const ::dealii::BlockSparseMatrix<double> &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 <typename Number, typename MemorySpace>
+ void
+ BlockSparseMatrix<Number, MemorySpace>::reinit(
+ const ::dealii::BlockSparseMatrix<double> &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<IndexSet> 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 <typename Number, typename MemorySpace>
+ inline bool
+ BlockSparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+ void
+ BlockSparseMatrix<Number, MemorySpace>::collect_sizes()
+ {
+ // simply forward to the (non-public) function of the base class
+ BaseClass::collect_sizes();
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ std::uint64_t
+ BlockSparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+ MPI_Comm
+ BlockSparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+ inline std::vector<IndexSet>
+ BlockSparseMatrix<Number, MemorySpace>::locally_owned_domain_indices() const
+ {
+ Assert(this->n_block_cols() != 0, ExcNotInitialized());
+ Assert(this->n_block_rows() != 0, ExcNotInitialized());
+
+ std::vector<IndexSet> 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 <typename Number, typename MemorySpace>
+ inline std::vector<IndexSet>
+ BlockSparseMatrix<Number, MemorySpace>::locally_owned_range_indices() const
+ {
+ Assert(this->n_block_cols() != 0, ExcNotInitialized());
+ Assert(this->n_block_rows() != 0, ExcNotInitialized());
+
+ std::vector<IndexSet> 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 <typename Number, typename MemorySpace>
+ template <typename VectorType1, typename VectorType2, typename VectorType3>
+ Number
+ BlockSparseMatrix<Number, MemorySpace>::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
# include <deal.II/lac/block_indices.h>
# include <deal.II/lac/block_vector_base.h>
+# include <deal.II/lac/trilinos_tpetra_block_sparse_matrix.h>
# include <deal.II/lac/trilinos_tpetra_vector.h>
DEAL_II_NAMESPACE_OPEN
template <typename Number, typename MemorySpace>
class BlockVector;
- // TODO
- // template <typename Number, typename MemorySpace>
- // class BlockSparseMatrix;
+ template <typename Number, typename MemorySpace>
+ class BlockSparseMatrix;
} // namespace TpetraWrappers
} // namespace LinearAlgebra
# endif // DOXYGEN
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
LinearAlgebra::TpetraWrappers::Vector<S> &,
bool,
std::integral_constant<bool, false>) const;
- }
+ // BlockSparseMatrix
+ template void AffineConstraints<S>::distribute_local_to_global<
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S>,
+ LinearAlgebra::TpetraWrappers::Vector<S>>(
+ const FullMatrix<S> &,
+ const Vector<S> &,
+ const std::vector<AffineConstraints<S>::size_type> &,
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S> &,
+ LinearAlgebra::TpetraWrappers::Vector<S> &,
+ bool,
+ std::bool_constant<true>) const;
+
+ template void AffineConstraints<S>::distribute_local_to_global<
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S>,
+ LinearAlgebra::TpetraWrappers::BlockVector<S>>(
+ const FullMatrix<S> &,
+ const Vector<S> &,
+ const std::vector<AffineConstraints<S>::size_type> &,
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S> &,
+ LinearAlgebra::TpetraWrappers::BlockVector<S> &,
+ bool,
+ std::bool_constant<true>) const;
+
+ template void AffineConstraints<S>::distribute_local_to_global<
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S>>(
+ const FullMatrix<S> &,
+ const std::vector<AffineConstraints<S>::size_type> &,
+ const std::vector<AffineConstraints<S>::size_type> &,
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S> &) const;
+
+ template void AffineConstraints<S>::distribute_local_to_global<
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S>>(
+ const FullMatrix<S> &,
+ const std::vector<AffineConstraints<S>::size_type> &,
+ const AffineConstraints<S> &,
+ const std::vector<AffineConstraints<S>::size_type> &,
+ LinearAlgebra::TpetraWrappers::BlockSparseMatrix<S> &) const;
+ }
// ---------------------------------------------------------------------
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+
+# include <deal.II/lac/trilinos_tpetra_block_sparse_matrix.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+# ifndef DOXYGEN
+// explicit instantiations
+namespace LinearAlgebra
+{
+ namespace TpetraWrappers
+ {
+ template class BlockSparseMatrix<double>;
+
+ template void
+ BlockSparseMatrix<double>::reinit(
+ const ::dealii::BlockDynamicSparsityPattern &);
+
+ template void
+ BlockSparseMatrix<double>::vmult(
+ TpetraWrappers::Vector<double> &,
+ const TpetraWrappers::Vector<double> &) const;
+
+ template void
+ BlockSparseMatrix<double>::Tvmult(
+ TpetraWrappers::Vector<double> &,
+ const TpetraWrappers::Vector<double> &) const;
+
+ template void
+ BlockSparseMatrix<double>::vmult(
+ TpetraWrappers::BlockVector<double> &,
+ const TpetraWrappers::BlockVector<double> &) const;
+
+ template void
+ BlockSparseMatrix<double>::Tvmult(
+ TpetraWrappers::BlockVector<double> &,
+ const TpetraWrappers::BlockVector<double> &) const;
+
+ template void
+ BlockSparseMatrix<double>::vmult(
+ ::dealii::BlockVector<double> &,
+ const ::dealii::BlockVector<double> &) const;
+
+ template void
+ BlockSparseMatrix<double>::Tvmult(
+ ::dealii::BlockVector<double> &,
+ const ::dealii::BlockVector<double> &) const;
+
+
+ } // namespace TpetraWrappers
+} // namespace LinearAlgebra
+# endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_TRILINOS_WITH_TPETRA
projects / dealii.git / commitdiff