block_sparsity_pattern.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2000 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_block_sparsity_pattern_h
#define dealii_block_sparsity_pattern_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table.h>
 
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/sparsity_pattern_base.h>
#include <deal.II/lac/trilinos_sparsity_pattern.h>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
template <typename number>
class BlockSparseMatrix;
class BlockDynamicSparsityPattern;
#endif
 
template <typename SparsityPatternType>
class BlockSparsityPatternBase : public SparsityPatternBase
{
public:
  using size_type = types::global_dof_index;
 
  static const size_type invalid_entry = SparsityPattern::invalid_entry;
 
  BlockSparsityPatternBase();
 
  BlockSparsityPatternBase(const size_type n_block_rows,
                           const size_type n_block_columns);
 
  BlockSparsityPatternBase(const BlockSparsityPatternBase &bsp);
 
  void
  reinit(const size_type n_block_rows, const size_type n_block_columns);
 
  BlockSparsityPatternBase &
  operator=(const BlockSparsityPatternBase &);
 
  void
  collect_sizes();
 
  SparsityPatternType &
  block(const size_type row, const size_type column);
 
  const SparsityPatternType &
  block(const size_type row, const size_type column) const;
 
  const BlockIndices &
  get_row_indices() const;
 
  const BlockIndices &
  get_column_indices() const;
 
  void
  compress();
 
  size_type
  n_block_rows() const;
 
  size_type
  n_block_cols() const;
 
  bool
  empty() const;
 
  size_type
  max_entries_per_row() const;
 
  void
  add(const size_type i, const size_type j);
 
  template <typename ForwardIterator>
  void
  add_entries(const size_type row,
              ForwardIterator begin,
              ForwardIterator end,
              const bool      indices_are_sorted = false);
 
  virtual void
  add_row_entries(const size_type                  &row,
                  const ArrayView<const size_type> &columns,
                  const bool indices_are_sorted = false) override;
 
  using SparsityPatternBase::add_entries;
 
  using SparsityPatternBase::n_rows;
 
  using SparsityPatternBase::n_cols;
 
  bool
  exists(const size_type i, const size_type j) const;
 
  unsigned int
  row_length(const size_type row) const;
 
  size_type
  n_nonzero_elements() const;
 
  void
  print(std::ostream &out) const;
 
  void
  print_gnuplot(std::ostream &out) const;
 
  void
  print_svg(std::ostream &out) const;
 
  std::size_t
  memory_consumption() const;
 
  DeclExceptionMsg(
    ExcNeedsCollectSizes,
    "The number of rows and columns (returned by n_rows() and n_cols()) does "
    "not match their directly computed values. This typically means that a "
    "call to collect_sizes() is missing.");
 
  DeclException4(ExcIncompatibleRowNumbers,
                 int,
                 int,
                 int,
                 int,
                 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
                 << ',' << arg4 << "] have differing row numbers.");
  DeclException4(ExcIncompatibleColNumbers,
                 int,
                 int,
                 int,
                 int,
                 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
                 << ',' << arg4 << "] have differing column numbers.");
protected:
  size_type block_rows;
 
  size_type block_columns;
 
  Table<2, std::unique_ptr<SparsityPatternType>> sub_objects;
 
  BlockIndices row_indices;
 
  BlockIndices column_indices;
 
private:
  size_type
  compute_n_rows() const;
 
  size_type
  compute_n_cols() const;
 
  std::vector<size_type> counter_within_block;
 
  std::vector<std::vector<size_type>> block_column_indices;
 
  // Make the block sparse matrix a friend, so that it can use our
  // #row_indices and #column_indices objects.
  template <typename number>
  friend class BlockSparseMatrix;
};
 
 
 
class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
{
public:
  BlockSparsityPattern() = default;
 
  BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
 
  void
  reinit(const size_type n_block_rows, const size_type n_block_columns);
 
  void
  reinit(const BlockIndices                           &row_indices,
         const BlockIndices                           &col_indices,
         const std::vector<std::vector<unsigned int>> &row_lengths);
 
 
  bool
  is_compressed() const;
 
  void
  copy_from(const BlockDynamicSparsityPattern &dsp);
};
 
 
 
class BlockDynamicSparsityPattern
  : public BlockSparsityPatternBase<DynamicSparsityPattern>
{
public:
  BlockDynamicSparsityPattern() = default;
 
  BlockDynamicSparsityPattern(const size_type n_rows,
                              const size_type n_columns);
 
  BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
                              const std::vector<size_type> &col_block_sizes);
 
  BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
 
  BlockDynamicSparsityPattern(const BlockIndices &row_indices,
                              const BlockIndices &col_indices);
 
 
  void
  reinit(const std::vector<size_type> &row_block_sizes,
         const std::vector<size_type> &col_block_sizes);
 
  void
  reinit(const std::vector<IndexSet> &partitioning);
 
  void
  reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
 
  size_type
  column_number(const size_type row, const unsigned int index) const;
 
  using BlockSparsityPatternBase<DynamicSparsityPattern>::reinit;
};
 
#ifdef DEAL_II_WITH_TRILINOS
 
 
namespace TrilinosWrappers
{
  class BlockSparsityPattern
    : public ::BlockSparsityPatternBase<SparsityPattern>
  {
  public:
    BlockSparsityPattern() = default;
 
    BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
 
    BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
                         const std::vector<size_type> &col_block_sizes);
 
    BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
                         const MPI_Comm communicator = MPI_COMM_WORLD);
 
    BlockSparsityPattern(
      const std::vector<IndexSet> &row_parallel_partitioning,
      const std::vector<IndexSet> &column_parallel_partitioning,
      const std::vector<IndexSet> &writeable_rows,
      const MPI_Comm               communicator = MPI_COMM_WORLD);
 
    void
    reinit(const std::vector<size_type> &row_block_sizes,
           const std::vector<size_type> &col_block_sizes);
 
    void
    reinit(const std::vector<IndexSet> &parallel_partitioning,
           const MPI_Comm               communicator = MPI_COMM_WORLD);
 
    void
    reinit(const std::vector<IndexSet> &row_parallel_partitioning,
           const std::vector<IndexSet> &column_parallel_partitioning,
           const MPI_Comm               communicator = MPI_COMM_WORLD);
 
    void
    reinit(const std::vector<IndexSet> &row_parallel_partitioning,
           const std::vector<IndexSet> &column_parallel_partitioning,
           const std::vector<IndexSet> &writeable_rows,
           const MPI_Comm               communicator = MPI_COMM_WORLD);
 
    using BlockSparsityPatternBase<SparsityPattern>::reinit;
  };
 
} /* namespace TrilinosWrappers */
 
#endif
 
/*--------------------- Template functions ----------------------------------*/
 
 
 
template <typename SparsityPatternType>
inline SparsityPatternType &
BlockSparsityPatternBase<SparsityPatternType>::block(const size_type row,
                                                     const size_type column)
{
  AssertIndexRange(row, n_block_rows());
  AssertIndexRange(column, n_block_cols());
  return *sub_objects(row, column);
}
 
 
 
template <typename SparsityPatternType>
inline const SparsityPatternType &
BlockSparsityPatternBase<SparsityPatternType>::block(
  const size_type row,
  const size_type column) const
{
  AssertIndexRange(row, n_block_rows());
  AssertIndexRange(column, n_block_cols());
  return *sub_objects(row, column);
}
 
 
 
template <typename SparsityPatternType>
inline const BlockIndices &
BlockSparsityPatternBase<SparsityPatternType>::get_row_indices() const
{
  return row_indices;
}
 
 
 
template <typename SparsityPatternType>
inline const BlockIndices &
BlockSparsityPatternBase<SparsityPatternType>::get_column_indices() const
{
  return column_indices;
}
 
 
 
template <typename SparsityPatternType>
inline void
BlockSparsityPatternBase<SparsityPatternType>::add(const size_type i,
                                                   const size_type j)
{
  // if you get an error here, are
  // you sure you called
  // <tt>collect_sizes()</tt> before?
  const std::pair<size_type, size_type> row_index =
                                          row_indices.global_to_local(i),
                                        col_index =
                                          column_indices.global_to_local(j);
  sub_objects[row_index.first][col_index.first]->add(row_index.second,
                                                     col_index.second);
}
 
 
 
template <typename SparsityPatternType>
template <typename ForwardIterator>
void
BlockSparsityPatternBase<SparsityPatternType>::add_entries(
  const size_type row,
  ForwardIterator begin,
  ForwardIterator end,
  const bool      indices_are_sorted)
{
  // In debug mode, verify that collect_sizes() was called by redoing the
  // calculation
  Assert(n_rows() == compute_n_rows(), ExcNeedsCollectSizes());
  Assert(n_cols() == compute_n_cols(), ExcNeedsCollectSizes());
 
  const size_type n_cols = static_cast<size_type>(end - begin);
 
  if (indices_are_sorted && n_cols > 0)
    {
      block_column_indices[0].resize(0);
 
      const std::pair<size_type, size_type> row_index =
        this->row_indices.global_to_local(row);
      const auto n_blocks = column_indices.size();
 
      // Assume we start with the first block: since we assemble sparsity
      // patterns one cell at a time this should always be true
      size_type current_block       = 0;
      size_type current_start_index = column_indices.block_start(current_block);
      size_type next_start_index =
        current_block == n_blocks - 1 ?
          numbers::invalid_dof_index :
          column_indices.block_start(current_block + 1);
 
      for (auto it = begin; it < end; ++it)
        {
          // BlockIndices::global_to_local() is a bit slow so instead we just
          // keep track of which block we are in - as the indices are sorted we
          // know that the block number can only increase.
          if (*it >= next_start_index)
            {
              // we found a column outside the present block: write the present
              // block entries and continue to the next block
              sub_objects[row_index.first][current_block]->add_entries(
                row_index.second,
                block_column_indices[0].begin(),
                block_column_indices[0].end(),
                true);
              block_column_indices[0].clear();
 
              auto block_and_col  = column_indices.global_to_local(*it);
              current_block       = block_and_col.first;
              current_start_index = column_indices.block_start(current_block);
              next_start_index =
                current_block == n_blocks - 1 ?
                  numbers::invalid_dof_index :
                  column_indices.block_start(current_block + 1);
            }
          const size_type local_index = *it - current_start_index;
          block_column_indices[0].push_back(local_index);
 
          // Check that calculation:
#ifdef DEBUG
          {
            auto check_block_and_col = column_indices.global_to_local(*it);
            Assert(current_block == check_block_and_col.first,
                   ExcInternalError());
            Assert(local_index == check_block_and_col.second,
                   ExcInternalError());
          }
#endif
        }
      // add whatever is left over:
      sub_objects[row_index.first][current_block]->add_entries(
        row_index.second,
        block_column_indices[0].begin(),
        block_column_indices[0].end(),
        true);
 
      return;
    }
  else
    {
      // Resize sub-arrays to n_cols. This
      // is a bit wasteful, but we resize
      // only a few times (then the maximum
      // row length won't increase that
      // much any more). At least we know
      // that all arrays are going to be of
      // the same size, so we can check
      // whether the size of one is large
      // enough before actually going
      // through all of them.
      if (block_column_indices[0].size() < n_cols)
        for (size_type i = 0; i < this->n_block_cols(); ++i)
          block_column_indices[i].resize(n_cols);
 
      // Reset the number of added elements
      // in each block to zero.
      for (size_type i = 0; i < this->n_block_cols(); ++i)
        counter_within_block[i] = 0;
 
      // Go through the column indices to
      // find out which portions of the
      // values should be set in which
      // block of the matrix. We need to
      // touch all the data, since we can't
      // be sure that the data of one block
      // is stored contiguously (in fact,
      // indices will be intermixed when it
      // comes from an element matrix).
      for (ForwardIterator it = begin; it != end; ++it)
        {
          const size_type col = *it;
 
          const std::pair<size_type, size_type> col_index =
            this->column_indices.global_to_local(col);
 
          const size_type local_index = counter_within_block[col_index.first]++;
 
          block_column_indices[col_index.first][local_index] = col_index.second;
        }
 
      // Now we found out about where the
      // individual columns should start and
      // where we should start reading out
      // data. Now let's write the data into
      // the individual blocks!
      const std::pair<size_type, size_type> row_index =
        this->row_indices.global_to_local(row);
      for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
        {
          if (counter_within_block[block_col] == 0)
            continue;
          sub_objects[row_index.first][block_col]->add_entries(
            row_index.second,
            block_column_indices[block_col].begin(),
            block_column_indices[block_col].begin() +
              counter_within_block[block_col],
            indices_are_sorted);
        }
    }
}
 
 
 
template <typename SparsityPatternType>
void
BlockSparsityPatternBase<SparsityPatternType>::add_row_entries(
  const size_type                  &row,
  const ArrayView<const size_type> &columns,
  const bool                        indices_are_sorted)
{
  add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
}
 
 
 
template <typename SparsityPatternType>
inline bool
BlockSparsityPatternBase<SparsityPatternType>::exists(const size_type i,
                                                      const size_type j) const
{
  // if you get an error here, are
  // you sure you called
  // <tt>collect_sizes()</tt> before?
  const std::pair<size_type, size_type> row_index =
                                          row_indices.global_to_local(i),
                                        col_index =
                                          column_indices.global_to_local(j);
  return sub_objects[row_index.first][col_index.first]->exists(
    row_index.second, col_index.second);
}
 
 
 
template <typename SparsityPatternType>
inline unsigned int
BlockSparsityPatternBase<SparsityPatternType>::row_length(
  const size_type row) const
{
  const std::pair<size_type, size_type> row_index =
    row_indices.global_to_local(row);
 
  unsigned int c = 0;
 
  for (size_type b = 0; b < n_block_rows(); ++b)
    c += sub_objects[row_index.first][b]->row_length(row_index.second);
 
  return c;
}
 
 
 
template <typename SparsityPatternType>
inline typename BlockSparsityPatternBase<SparsityPatternType>::size_type
BlockSparsityPatternBase<SparsityPatternType>::n_block_cols() const
{
  return block_columns;
}
 
 
 
template <typename SparsityPatternType>
inline typename BlockSparsityPatternBase<SparsityPatternType>::size_type
BlockSparsityPatternBase<SparsityPatternType>::n_block_rows() const
{
  return block_rows;
}
 
 
inline BlockDynamicSparsityPattern::size_type
BlockDynamicSparsityPattern::column_number(const size_type    row,
                                           const unsigned int index) const
{
  // .first= ith block, .second = jth row in that block
  const std::pair<size_type, size_type> row_index =
    row_indices.global_to_local(row);
 
  AssertIndexRange(index, row_length(row));
 
  size_type c             = 0;
  size_type block_columns = 0; // sum of n_cols for all blocks to the left
  for (unsigned int b = 0; b < this->n_block_cols(); ++b)
    {
      unsigned int rowlen =
        sub_objects[row_index.first][b]->row_length(row_index.second);
      if (index < c + rowlen)
        return block_columns +
               sub_objects[row_index.first][b]->column_number(row_index.second,
                                                              index - c);
      c += rowlen;
      block_columns += sub_objects[row_index.first][b]->n_cols();
    }
 
  DEAL_II_ASSERT_UNREACHABLE();
  return 0;
}
 
 
inline void
BlockSparsityPattern::reinit(const size_type new_block_rows,
                             const size_type new_block_columns)
{
  BlockSparsityPatternBase<SparsityPattern>::reinit(new_block_rows,
                                                    new_block_columns);
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif