block_matrix_array.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2016 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/pointer_matrix.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/block_vector.h>

#include <vector>
#include <map>
#include <string>
#include <memory>
#include <sstream>

DEAL_II_NAMESPACE_OPEN

template <typename number = double, typename BlockVectorType=BlockVector<number> >
class BlockMatrixArray : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

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

  std::vector<Entry> entries;

private:
  unsigned int block_rows;
  unsigned int block_cols;

  Entry &operator= (const Entry &);
};

template <typename number = double, typename BlockVectorType = BlockVector<number> >
class BlockTrianglePrecondition
  : private BlockMatrixArray<number,BlockVectorType>
{
public:
  typedef types::global_dof_index size_type;

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

  typedef std::map<const PointerMatrixBase<typename BlockVectorType::BlockType > *, std::string> NameMap;
  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) != "" && 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