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