block_sparsity_pattern.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2018 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_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/trilinos_sparsity_pattern.h>

DEAL_II_NAMESPACE_OPEN


template <typename number>
class BlockSparseMatrix;
class BlockDynamicSparsityPattern;

template <typename SparsityPatternType>
class BlockSparsityPatternBase : public Subscriptor
{
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);

  ~BlockSparsityPatternBase() override;

  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);

  size_type
  n_rows() const;

  size_type
  n_cols() const;

  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;

  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 rows;

  size_type columns;

  Table<2,
        SmartPointer<SparsityPatternType,
                     BlockSparsityPatternBase<SparsityPatternType>>>
    sub_objects;

  BlockIndices row_indices;

  BlockIndices column_indices;

private:
  std::vector<size_type> counter_within_block;

  std::vector<std::vector<size_type>> block_column_indices;

  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;

  std::size_t
  memory_consumption() 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);

    DEAL_II_DEPRECATED
    BlockSparsityPattern(const std::vector<Epetra_Map> &parallel_partitioning);

    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);

    DEAL_II_DEPRECATED
    void
    reinit(const std::vector<Epetra_Map> &parallel_partitioning);

    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)
{
  Assert(row < rows, ExcIndexRange(row, 0, rows));
  Assert(column < columns, ExcIndexRange(column, 0, columns));
  return *sub_objects[row][column];
}



template <typename SparsityPatternType>
inline const SparsityPatternType &
BlockSparsityPatternBase<SparsityPatternType>::block(
  const size_type row,
  const size_type column) const
{
  Assert(row < rows, ExcIndexRange(row, 0, rows));
  Assert(column < columns, ExcIndexRange(column, 0, columns));
  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)
{
  // Resize scratch arrays
  if (block_column_indices.size() < this->n_block_cols())
    {
      block_column_indices.resize(this->n_block_cols());
      counter_within_block.resize(this->n_block_cols());
    }

  const size_type n_cols = static_cast<size_type>(end - begin);

  // 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;
    }

#ifdef DEBUG
  // If in debug mode, do a check whether
  // the right length has been obtained.
  size_type length = 0;
  for (size_type i = 0; i < this->n_block_cols(); ++i)
    length += counter_within_block[i];
  Assert(length == n_cols, ExcInternalError());
#endif

  // 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>
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 < 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 columns;
}



template <typename SparsityPatternType>
inline typename BlockSparsityPatternBase<SparsityPatternType>::size_type
BlockSparsityPatternBase<SparsityPatternType>::n_block_rows() const
{
  return 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);

  Assert(index < row_length(row), ExcIndexRange(index, 0, 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 < columns; ++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();
    }

  Assert(false, ExcInternalError());
  return 0;
}


inline void
BlockSparsityPattern::reinit(const size_type n_block_rows,
                             const size_type n_block_columns)
{
  BlockSparsityPatternBase<SparsityPattern>::reinit(n_block_rows,
                                                    n_block_columns);
}


DEAL_II_NAMESPACE_CLOSE

#endif