block_matrix_array.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2019 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_matrix_array_h
#define dealii_block_matrix_array_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table.h>
 
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/pointer_matrix.h>
#include <deal.II/lac/vector_memory.h>
 
#include <map>
#include <memory>
#include <sstream>
#include <string>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename number          = double,
          typename BlockVectorType = BlockVector<number>>
class DEAL_II_DEPRECATED BlockMatrixArray : public Subscriptor
{
public:
  using size_type = types::global_dof_index;
 
  BlockMatrixArray();
 
  BlockMatrixArray(const unsigned int n_block_rows,
                   const unsigned int n_block_cols);
 
  void
  initialize(const unsigned int n_block_rows, const unsigned int n_block_cols);
 
  void
  reinit(const unsigned int n_block_rows, const unsigned int n_block_cols);
 
  template <typename MatrixType>
  void
  enter(const MatrixType & matrix,
        const unsigned int row,
        const unsigned int col,
        const number       prefix    = 1.,
        const bool         transpose = false);
 
  void
  clear();
 
  unsigned int
  n_block_rows() const;
 
  unsigned int
  n_block_cols() const;
 
  void
  vmult(BlockVectorType &dst, const BlockVectorType &src) const;
 
  void
  vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
 
  void
  Tvmult(BlockVectorType &dst, const BlockVectorType &src) const;
 
  void
  Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
 
  number
  matrix_scalar_product(const BlockVectorType &u,
                        const BlockVectorType &v) const;
 
  number
  matrix_norm_square(const BlockVectorType &u) const;
 
  template <class StreamType>
  void
  print_latex(StreamType &out) const;
 
protected:
  class Entry
  {
  public:
    template <typename MatrixType>
    Entry(const MatrixType &matrix,
          size_type         row,
          size_type         col,
          number            prefix,
          bool              transpose);
 
    Entry(const Entry &);
 
    ~Entry();
 
    size_type row;
 
    size_type col;
 
    number prefix;
 
    bool transpose;
 
    PointerMatrixBase<typename BlockVectorType::BlockType> *matrix;
 
    Entry &
    operator=(const Entry &) = delete;
  };
 
  std::vector<Entry> entries;
 
private:
  unsigned int block_rows;
  unsigned int block_cols;
};
 
template <typename number          = double,
          typename BlockVectorType = BlockVector<number>>
class DEAL_II_DEPRECATED BlockTrianglePrecondition
  : private BlockMatrixArray<number, BlockVectorType>
{
public:
  using size_type = types::global_dof_index;
 
  BlockTrianglePrecondition();
 
  BlockTrianglePrecondition(const unsigned int n_blocks);
 
  void
  reinit(const unsigned int n_block_rows);
 
 
  template <typename MatrixType>
  void
  enter(const MatrixType &matrix,
        const size_type   row,
        const size_type   col,
        const number      prefix    = 1.,
        const bool        transpose = false);
 
  void
  vmult(BlockVectorType &dst, const BlockVectorType &src) const;
 
  void
  vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
 
  void
  Tvmult(BlockVectorType &dst, const BlockVectorType &src) const;
 
  void
  Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
 
  using BlockMatrixArray<number, BlockVectorType>::print_latex;
 
  using BlockMatrixArray<number, BlockVectorType>::n_block_rows;
 
  using BlockMatrixArray<number, BlockVectorType>::n_block_cols;
  using BlockMatrixArray<number, BlockVectorType>::clear;
  using BlockMatrixArray<number, BlockVectorType>::Subscriptor::subscribe;
  using BlockMatrixArray<number, BlockVectorType>::Subscriptor::unsubscribe;
 
  DeclException1(ExcNoDiagonal,
                 size_type,
                 << "No diagonal entry was added for block " << arg1);
 
  DeclException1(ExcMultipleDiagonal,
                 size_type,
                 << "Inverse diagonal entries may not be added in block "
                 << arg1);
private:
  void
  do_row(BlockVectorType &dst, size_type row_num) const;
 
  bool backward;
};
 
 
#ifndef DOXYGEN
//---------------------------------------------------------------------------
 
template <typename number, typename BlockVectorType>
template <typename MatrixType>
inline BlockMatrixArray<number, BlockVectorType>::Entry::Entry(
  const MatrixType &m,
  size_type         row,
  size_type         col,
  number            prefix,
  bool              transpose)
  : row(row)
  , col(col)
  , prefix(prefix)
  , transpose(transpose)
  , matrix(new_pointer_matrix_base(m,
                                   typename BlockVectorType::BlockType(),
                                   typeid(*this).name()))
{}
 
 
 
template <typename number, typename BlockVectorType>
template <typename MatrixType>
inline void
BlockMatrixArray<number, BlockVectorType>::enter(const MatrixType &matrix,
                                                 unsigned int      row,
                                                 unsigned int      col,
                                                 number            prefix,
                                                 bool              transpose)
{
  Assert(row < n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
  Assert(col < n_block_cols(), ExcIndexRange(col, 0, n_block_cols()));
  entries.push_back(Entry(matrix, row, col, prefix, transpose));
}
 
 
template <typename number, typename BlockVectorType>
template <class StreamType>
inline void
BlockMatrixArray<number, BlockVectorType>::print_latex(StreamType &out) const
{
  out << "\\begin{array}{" << std::string(n_block_cols(), 'c') << "}"
      << std::endl;
 
  Table<2, std::string> array(n_block_rows(), n_block_cols());
 
  using NameMap =
    std::map<const PointerMatrixBase<typename BlockVectorType::BlockType> *,
             std::string>;
  NameMap matrix_names;
 
  typename std::vector<Entry>::const_iterator m   = entries.begin();
  typename std::vector<Entry>::const_iterator end = entries.end();
 
  size_type matrix_number = 0;
  for (; m != end; ++m)
    {
      if (matrix_names.find(m->matrix) == matrix_names.end())
        {
          std::pair<typename NameMap::iterator, bool> x = matrix_names.insert(
            std::pair<
              const PointerMatrixBase<typename BlockVectorType::BlockType> *,
              std::string>(m->matrix, std::string("M")));
          std::ostringstream stream;
          stream << matrix_number++;
 
          x.first->second += stream.str();
        }
 
      std::ostringstream stream;
 
      if (!array(m->row, m->col).empty() && m->prefix >= 0)
        stream << "+";
      if (m->prefix != 1.)
        stream << m->prefix << 'x';
      stream << matrix_names.find(m->matrix)->second;
      //      stream << '(' << m->matrix << ')';
      if (m->transpose)
        stream << "^T";
 
      array(m->row, m->col) += stream.str();
    }
  for (unsigned int i = 0; i < n_block_rows(); ++i)
    for (unsigned int j = 0; j < n_block_cols(); ++j)
      {
        out << '\t' << array(i, j);
        if (j == n_block_cols() - 1)
          {
            if (i != n_block_rows() - 1)
              out << "\\\\" << std::endl;
            else
              out << std::endl;
          }
        else
          out << " &";
      }
  out << "\\end{array}" << std::endl;
}
 
template <typename number, typename BlockVectorType>
template <typename MatrixType>
inline void
BlockTrianglePrecondition<number, BlockVectorType>::enter(
  const MatrixType &matrix,
  size_type         row,
  size_type         col,
  number            prefix,
  bool              transpose)
{
  BlockMatrixArray<number, BlockVectorType>::enter(
    matrix, row, col, prefix, transpose);
}
 
 
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif