trilinos_block_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 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__trilinos_block_vector_h
#define dealii__trilinos_block_vector_h


#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_TRILINOS

#  include <deal.II/lac/trilinos_vector.h>
#  include <deal.II/lac/trilinos_parallel_block_vector.h>
#  include <deal.II/lac/block_indices.h>
#  include <deal.II/lac/block_vector_base.h>
#  include <deal.II/lac/exceptions.h>

DEAL_II_NAMESPACE_OPEN

// forward declaration
template <typename Number> class BlockVector;

namespace TrilinosWrappers
{
  // forward declaration
  namespace MPI
  {
    class BlockVector;
  }
  class BlockVector;
  class BlockSparseMatrix;


  class BlockVector : public BlockVectorBase<Vector>
  {
  public:
    typedef BlockVectorBase<Vector> BaseClass;

    typedef BaseClass::BlockType  BlockType;

    typedef BaseClass::value_type      value_type;
    typedef BaseClass::pointer         pointer;
    typedef BaseClass::const_pointer   const_pointer;
    typedef BaseClass::reference       reference;
    typedef BaseClass::const_reference const_reference;
    typedef BaseClass::size_type       size_type;
    typedef BaseClass::iterator        iterator;
    typedef BaseClass::const_iterator  const_iterator;

    BlockVector ();

    explicit BlockVector (const std::vector<Epetra_Map> &partitioner) DEAL_II_DEPRECATED;

    explicit BlockVector (const std::vector<IndexSet> &partitioner,
                          const MPI_Comm              &communicator = MPI_COMM_WORLD) DEAL_II_DEPRECATED;

    BlockVector (const MPI::BlockVector &V) DEAL_II_DEPRECATED;

    BlockVector (const BlockVector  &V) DEAL_II_DEPRECATED;

    explicit BlockVector (const size_type num_blocks) DEAL_II_DEPRECATED;

    explicit BlockVector (const std::vector<size_type> &N) DEAL_II_DEPRECATED;

    template <typename InputIterator>
    BlockVector (const std::vector<size_type> &n,
                 const InputIterator           first,
                 const InputIterator           end) DEAL_II_DEPRECATED;

    ~BlockVector ();

    BlockVector &
    operator = (const value_type s);

    BlockVector &
    operator = (const MPI::BlockVector &V);

    BlockVector &
    operator = (const BlockVector &V);

    template <typename Number>
    BlockVector &
    operator = (const ::BlockVector<Number> &V);

    void reinit (const std::vector<Epetra_Map> &partitioning,
                 const bool                     omit_zeroing_entries = false);

    void reinit (const std::vector<IndexSet> &partitioning,
                 const MPI_Comm              &communicator = MPI_COMM_WORLD,
                 const bool                   omit_zeroing_entries = false);

    void reinit (const std::vector<size_type> &N,
                 const bool                    omit_zeroing_entries=false);

    void reinit (const MPI::BlockVector &V);

    void reinit (const BlockVector &V,
                 const bool omit_zeroing_entries = false);

    void reinit (const size_type num_blocks);

    void swap (BlockVector &v);

    void print (std::ostream       &out,
                const unsigned int  precision = 3,
                const bool          scientific = true,
                const bool          across = true) const;

    DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);

    DeclException0 (ExcNonMatchingBlockVectors);

    DeclException2 (ExcNonLocalizedMap,
                    int, int,
                    << "For the generation of a localized vector the map has "
                    << "to assign all elements to all vectors! "
                    << "local_size = global_size is a necessary condition, but"
                    << arg1 << " != " << arg2 << " was given!");

  };



  /*----------------------- Inline functions ----------------------------------*/



  inline
  BlockVector::BlockVector ()
  {}



  inline
  BlockVector::BlockVector (const std::vector<Epetra_Map> &partitioning)
  {
    reinit (partitioning);
  }



  inline
  BlockVector::BlockVector (const std::vector<IndexSet> &partitioning,
                            const MPI_Comm              &communicator)
  {
    reinit (partitioning, communicator);
  }



  inline
  BlockVector::BlockVector (const std::vector<size_type> &N)
  {
    reinit (N);
  }



  template <typename InputIterator>
  BlockVector::BlockVector (const std::vector<size_type> &n,
                            const InputIterator           first,
                            const InputIterator           end)
  {
    // first set sizes of blocks, but
    // don't initialize them as we will
    // copy elements soon
    (void)end;
    reinit (n, true);
    InputIterator start = first;
    for (size_type b=0; b<n.size(); ++b)
      {
        InputIterator end = start;
        std::advance (end, static_cast<size_type>(n[b]));

        for (size_type i=0; i<n[b]; ++i, ++start)
          this->block(b)(i) = *start;
      }
    Assert (start == end, ExcIteratorRangeDoesNotMatchVectorSize());
  }



  inline
  BlockVector::BlockVector (const size_type num_blocks)
  {
    reinit (num_blocks);
  }



  inline
  BlockVector::~BlockVector()
  {}



  inline
  BlockVector::BlockVector (const MPI::BlockVector &v)
  {
    reinit (v);
  }



  inline
  BlockVector::BlockVector (const BlockVector &v)
    :
    BlockVectorBase<Vector > ()
  {
    this->components.resize (v.n_blocks());
    this->block_indices = v.block_indices;

    for (size_type i=0; i<this->n_blocks(); ++i)
      this->components[i] = v.components[i];
  }


  inline
  void
  BlockVector::swap (BlockVector &v)
  {
    Assert (n_blocks() == v.n_blocks(),
            ExcDimensionMismatch(n_blocks(),v.n_blocks()));

    for (unsigned int row=0; row<n_blocks(); ++row)
      block(row).swap (v.block(row));
  }


  template <typename Number>
  BlockVector &
  BlockVector::operator = (const ::BlockVector<Number> &v)
  {
    if (n_blocks() != v.n_blocks())
      {
        std::vector<size_type> block_sizes (v.n_blocks(), 0);
        block_indices.reinit (block_sizes);
        if (components.size() != n_blocks())
          components.resize(n_blocks());
      }

    for (size_type i=0; i<this->n_blocks(); ++i)
      this->components[i] = v.block(i);

    collect_sizes();

    return *this;
  }


  inline
  void swap (BlockVector &u,
             BlockVector &v)
  {
    u.swap (v);
  }

} /* namespace TrilinosWrappers */

namespace internal
{
  namespace LinearOperator
  {
    template <typename> class ReinitHelper;

    template<>
    class ReinitHelper<TrilinosWrappers::BlockVector>
    {
    public:
      template <typename Matrix>
      static
      void reinit_range_vector (const Matrix &matrix,
                                TrilinosWrappers::BlockVector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit(matrix.range_partitioner(), omit_zeroing_entries);
      }

      template <typename Matrix>
      static
      void reinit_domain_vector(const Matrix &matrix,
                                TrilinosWrappers::BlockVector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit(matrix.domain_partitioner(), omit_zeroing_entries);
      }
    };

  } /* namespace LinearOperator */
} /* namespace internal */


DEAL_II_NAMESPACE_CLOSE

#endif  // DEAL_II_WITH_TRILINOS

#endif