]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add LinearAlgebra::distributed::BlockVector
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 14 Jun 2016 16:54:32 +0000 (18:54 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 15 Jun 2016 08:46:01 +0000 (10:46 +0200)
14 files changed:
include/deal.II/lac/la_parallel_block_vector.h [new file with mode: 0644]
include/deal.II/lac/la_parallel_block_vector.templates.h [new file with mode: 0644]
include/deal.II/lac/la_parallel_vector.h
include/deal.II/lac/parallel_block_vector.h
include/deal.II/lac/parallel_vector.h
source/lac/CMakeLists.txt
source/lac/la_parallel_block_vector.cc [new file with mode: 0644]
source/lac/la_parallel_block_vector.inst.in [new file with mode: 0644]
source/lac/la_parallel_vector.inst.in
source/lac/trilinos_epetra_communication_pattern.cc
tests/mpi/parallel_block_vector_01.cc
tests/mpi/parallel_block_vector_01.mpirun=1.output
tests/mpi/parallel_block_vector_01.mpirun=10.output
tests/mpi/parallel_block_vector_01.mpirun=4.output

diff --git a/include/deal.II/lac/la_parallel_block_vector.h b/include/deal.II/lac/la_parallel_block_vector.h
new file mode 100644 (file)
index 0000000..11691b6
--- /dev/null
@@ -0,0 +1,599 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1999 - 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__la_parallel_block_vector_h
+#define dealii__la_parallel_block_vector_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/lac/block_indices.h>
+#include <deal.II/lac/block_vector_base.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+
+#include <cstdio>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+#ifdef DEAL_II_WITH_PETSC
+namespace PETScWrappers
+{
+  namespace MPI
+  {
+    class BlockVector;
+  }
+}
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+namespace TrilinosWrappers
+{
+  namespace MPI
+  {
+    class BlockVector;
+  }
+}
+#endif
+
+namespace LinearAlgebra
+{
+  namespace distributed
+  {
+
+    /*! @addtogroup Vectors
+     *@{
+     */
+
+
+    /**
+     * An implementation of block vectors based on distributed deal.II
+     * vectors. While the base class provides for most of the interface, this
+     * class handles the actual allocation of vectors and provides functions
+     * that are specific to the underlying vector type.
+     *
+     * @note Instantiations for this template are provided for <tt>@<float@>
+     * and @<double@></tt>; others can be generated in application programs
+     * (see the section on
+     * @ref Instantiations
+     * in the manual).
+     *
+     * @see
+     * @ref GlossBlockLA "Block (linear algebra)"
+     * @author Katharina Kormann, Martin Kronbichler, 2011
+     */
+    template <typename Number>
+    class BlockVector : public BlockVectorBase<Vector<Number> >,
+      public VectorSpaceVector<Number>
+    {
+    public:
+      /**
+       * Typedef the base class for simpler access to its own typedefs.
+       */
+      typedef BlockVectorBase<Vector<Number> > BaseClass;
+
+      /**
+       * Typedef the type of the underlying vector.
+       */
+      typedef typename BaseClass::BlockType  BlockType;
+
+      /**
+       * Import the typedefs from the base class.
+       */
+      typedef typename BaseClass::value_type      value_type;
+      typedef typename BaseClass::real_type       real_type;
+      typedef typename BaseClass::pointer         pointer;
+      typedef typename BaseClass::const_pointer   const_pointer;
+      typedef typename BaseClass::reference       reference;
+      typedef typename BaseClass::const_reference const_reference;
+      typedef typename BaseClass::size_type       size_type;
+      typedef typename BaseClass::iterator        iterator;
+      typedef typename BaseClass::const_iterator  const_iterator;
+
+      /**
+       * @name 1: Basic operations
+       */
+      //@{
+
+      /**
+       * Constructor. There are three ways to use this constructor. First,
+       * without any arguments, it generates an object with no blocks. Given
+       * one argument, it initializes <tt>num_blocks</tt> blocks, but these
+       * blocks have size zero. The third variant finally initializes all
+       * blocks to the same size <tt>block_size</tt>.
+       *
+       * Confer the other constructor further down if you intend to use blocks
+       * of different sizes.
+       */
+      explicit BlockVector (const size_type num_blocks = 0,
+                            const size_type block_size = 0);
+
+      /**
+       * Copy-Constructor. Dimension set to that of V, all components are
+       * copied from V
+       */
+      BlockVector (const BlockVector<Number> &V);
+
+
+#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
+      /**
+       * Copy constructor taking a BlockVector of another data type. This will
+       * fail if there is no conversion path from <tt>OtherNumber</tt> to
+       * <tt>Number</tt>. Note that you may lose accuracy when copying to a
+       * BlockVector with data elements with less accuracy.
+       *
+       * Older versions of gcc did not honor the @p explicit keyword on
+       * template constructors. In such cases, it is easy to accidentally
+       * write code that can be very inefficient, since the compiler starts
+       * performing hidden conversions. To avoid this, this function is
+       * disabled if we have detected a broken compiler during configuration.
+       */
+      template <typename OtherNumber>
+      explicit
+      BlockVector (const BlockVector<OtherNumber> &v);
+#endif
+
+      /**
+       * Constructor. Set the number of blocks to <tt>block_sizes.size()</tt>
+       * and initialize each block with <tt>block_sizes[i]</tt> zero elements.
+       */
+      BlockVector (const std::vector<size_type> &block_sizes);
+
+      /**
+       * Construct a block vector with an IndexSet for the local range and
+       * ghost entries for each block.
+       */
+      BlockVector (const std::vector<IndexSet> &local_ranges,
+                   const std::vector<IndexSet> &ghost_indices,
+                   const MPI_Comm  communicator);
+
+      /**
+       * Same as above but the ghost indices are assumed to be empty.
+       */
+      BlockVector (const std::vector<IndexSet> &local_ranges,
+                   const MPI_Comm  communicator);
+
+      /**
+       * Destructor. Clears memory.
+       */
+      ~BlockVector ();
+
+      /**
+       * Copy operator: fill all components of the vector with the given
+       * scalar value.
+       */
+      BlockVector &operator = (const value_type s);
+
+      /**
+       * Copy operator for arguments of the same type. Resize the present
+       * vector if necessary.
+       */
+      BlockVector &
+      operator= (const BlockVector &V);
+
+      /**
+       * Copy operator for template arguments of different types. Resize the
+       * present vector if necessary.
+       */
+      template <class Number2>
+      BlockVector &
+      operator= (const BlockVector<Number2> &V);
+
+      /**
+       * Copy a regular vector into a block vector.
+       */
+      BlockVector &
+      operator= (const Vector<Number> &V);
+
+#ifdef DEAL_II_WITH_PETSC
+      /**
+       * Copy the content of a PETSc vector into the calling vector. This
+       * function assumes that the vectors layouts have already been
+       * initialized to match.
+       *
+       * This operator is only available if deal.II was configured with PETSc.
+       */
+      BlockVector<Number> &
+      operator = (const PETScWrappers::MPI::BlockVector &petsc_vec);
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+      /**
+       * Copy the content of a Trilinos vector into the calling vector. This
+       * function assumes that the vectors layouts have already been
+       * initialized to match.
+       *
+       * This operator is only available if deal.II was configured with
+       * Trilinos.
+       */
+      BlockVector<Number> &
+      operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec);
+#endif
+
+      /**
+       * Reinitialize the BlockVector to contain <tt>num_blocks</tt> blocks of
+       * size <tt>block_size</tt> each.
+       *
+       * If the second argument is left at its default value, then the block
+       * vector allocates the specified number of blocks but leaves them at
+       * zero size. You then need to later reinitialize the individual blocks,
+       * and call collect_sizes() to update the block system's knowledge of
+       * its individual block's sizes.
+       *
+       * If <tt>omit_zeroing_entries==false</tt>, the vector is filled with
+       * zeros.
+       */
+      void reinit (const size_type num_blocks,
+                   const size_type block_size = 0,
+                   const bool omit_zeroing_entries = false);
+
+      /**
+       * Reinitialize the BlockVector such that it contains
+       * <tt>block_sizes.size()</tt> blocks. Each block is reinitialized to
+       * dimension <tt>block_sizes[i]</tt>.
+       *
+       * If the number of blocks is the same as before this function was
+       * called, all vectors remain the same and reinit() is called for each
+       * vector.
+       *
+       * If <tt>omit_zeroing_entries==false</tt>, the vector is filled with
+       * zeros.
+       *
+       * Note that you must call this (or the other reinit() functions)
+       * function, rather than calling the reinit() functions of an individual
+       * block, to allow the block vector to update its caches of vector
+       * sizes. If you call reinit() on one of the blocks, then subsequent
+       * actions on this object may yield unpredictable results since they may
+       * be routed to the wrong block.
+       */
+      void reinit (const std::vector<size_type> &N,
+                   const bool                    omit_zeroing_entries=false);
+
+      /**
+       * Change the dimension to that of the vector <tt>V</tt>. The same
+       * applies as for the other reinit() function.
+       *
+       * The elements of <tt>V</tt> are not copied, i.e.  this function is the
+       * same as calling <tt>reinit (V.size(), omit_zeroing_entries)</tt>.
+       *
+       * Note that you must call this (or the other reinit() functions)
+       * function, rather than calling the reinit() functions of an individual
+       * block, to allow the block vector to update its caches of vector
+       * sizes. If you call reinit() of one of the blocks, then subsequent
+       * actions of this object may yield unpredictable results since they may
+       * be routed to the wrong block.
+       */
+      template <typename Number2>
+      void reinit (const BlockVector<Number2> &V,
+                   const bool                 omit_zeroing_entries=false);
+
+      /**
+       * This function copies the data that has accumulated in the data buffer
+       * for ghost indices to the owning processor. For the meaning of the
+       * argument @p operation, see the entry on
+       * @ref GlossCompress "Compressing distributed vectors and matrices"
+       * in the glossary.
+       *
+       * There are two variants for this function. If called with argument @p
+       * VectorOperation::add adds all the data accumulated in ghost elements
+       * to the respective elements on the owning processor and clears the
+       * ghost array afterwards. If called with argument @p
+       * VectorOperation::insert, a set operation is performed. Since setting
+       * elements in a vector with ghost elements is ambiguous (as one can set
+       * both the element on the ghost site as well as the owning site), this
+       * operation makes the assumption that all data is set correctly on the
+       * owning processor. Upon call of compress(VectorOperation::insert), all
+       * ghost entries are therefore simply zeroed out (using
+       * zero_ghost_values()). In debug mode, a check is performed that makes
+       * sure that the data set is actually consistent between processors,
+       * i.e., whenever a non-zero ghost element is found, it is compared to
+       * the value on the owning processor and an exception is thrown if these
+       * elements do not agree.
+       *
+       */
+      void compress (::dealii::VectorOperation::values operation);
+
+      /**
+       * Fills the data field for ghost indices with the values stored in the
+       * respective positions of the owning processor. This function is needed
+       * before reading from ghosts. The function is @p const even though
+       * ghost data is changed. This is needed to allow functions with a @p
+       * const vector to perform the data exchange without creating
+       * temporaries.
+       */
+      void update_ghost_values () const;
+
+      /**
+       * This method zeros the entries on ghost dofs, but does not touch
+       * locally owned DoFs.
+       *
+       * After calling this method, read access to ghost elements of the
+       * vector is forbidden and an exception is thrown. Only write access to
+       * ghost elements is allowed in this state.
+       */
+      void zero_out_ghosts ();
+
+      /**
+       * Returns if this Vector contains ghost elements.
+       */
+      bool has_ghost_elements() const;
+
+      /**
+       * This is a collective add operation that adds a whole set of values
+       * stored in @p values to the vector components specified by @p indices.
+       */
+      template <typename OtherNumber>
+      void add (const std::vector<size_type>        &indices,
+                const ::dealii::Vector<OtherNumber> &values);
+
+      /**
+       * Scaling and simple vector addition, i.e.  <tt>*this =
+       * s*(*this)+V</tt>.
+       */
+      void sadd (const Number               s,
+                 const BlockVector<Number> &V);
+
+      /**
+       * Assignment <tt>*this = a*u + b*v</tt>.
+       *
+       * This function is deprecated.
+       */
+      void equ (const Number a, const BlockVector<Number> &u,
+                const Number b, const BlockVector<Number> &v) DEAL_II_DEPRECATED;
+
+      /**
+       * Scaling and multiple addition.
+       *
+       * This function is deprecated.
+       */
+      void sadd (const Number               s,
+                 const Number               a,
+                 const BlockVector<Number> &V,
+                 const Number               b,
+                 const BlockVector<Number> &W) DEAL_II_DEPRECATED;
+
+      /**
+       * Return whether the vector contains only elements with value zero.
+       * This function is mainly for internal consistency checks and should
+       * seldom be used when not in debug mode since it uses quite some time.
+       */
+      bool all_zero () const;
+
+      /**
+       * Computes the mean value of all the entries in the vector.
+       */
+      Number mean_value () const;
+
+      /**
+       * $l_p$-norm of the vector. The pth root of the sum of the pth powers
+       * of the absolute values of the elements.
+       */
+      real_type lp_norm (const real_type p) const;
+
+      /**
+       * Swap the contents of this vector and the other vector <tt>v</tt>. One
+       * could do this operation with a temporary variable and copying over
+       * the data elements, but this function is significantly more efficient
+       * since it only swaps the pointers to the data of the two vectors and
+       * therefore does not need to allocate temporary storage and move data
+       * around.
+       *
+       * Limitation: right now this function only works if both vectors have
+       * the same number of blocks. If needed, the numbers of blocks should be
+       * exchanged, too.
+       *
+       * This function is analog to the the swap() function of all C++
+       * standard containers. Also, there is a global function swap(u,v) that
+       * simply calls <tt>u.swap(v)</tt>, again in analogy to standard
+       * functions.
+       */
+      void swap (BlockVector<Number> &v);
+      //@}
+
+      /**
+       * @name 2: Implementation of VectorSpaceVector
+       */
+      //@{
+
+      /**
+       * Multiply the entire vector by a fixed factor.
+       */
+      virtual BlockVector<Number> &operator*= (const Number factor);
+
+      /**
+       * Divide the entire vector by a fixed factor.
+       */
+      virtual BlockVector<Number> &operator/= (const Number factor);
+
+      /**
+       * Add the vector @p V to the present one.
+       */
+      virtual BlockVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+
+      /**
+       * Subtract the vector @p V from the present one.
+       */
+      virtual BlockVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+
+      /**
+       * Import all the elements present in the vector's IndexSet from the input
+       * vector @p V. VectorOperation::values @p operation is used to decide if
+       * the elements in @p V should be added to the current vector or replace the
+       * current elements. The last parameter can be used if the same
+       * communication pattern is used multiple times. This can be used to improve
+       * performance.
+       */
+      virtual void import(const LinearAlgebra::ReadWriteVector<Number> &V,
+                          VectorOperation::values operation,
+                          std_cxx11::shared_ptr<const CommunicationPatternBase> communication_pattern =
+                            std_cxx11::shared_ptr<const CommunicationPatternBase> ());
+
+      /**
+       * Return the scalar product of two vectors.
+       */
+      virtual Number operator* (const VectorSpaceVector<Number> &V) const;
+
+      /**
+       * Add @p a to all components. Note that @p a is a scalar not a vector.
+       */
+      virtual void add(const Number a);
+
+      /**
+       * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+       */
+      virtual void add(const Number a, const VectorSpaceVector<Number> &V);
+
+      /**
+       * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
+       */
+      virtual void add(const Number a, const VectorSpaceVector<Number> &V,
+                       const Number b, const VectorSpaceVector<Number> &W);
+
+      /**
+       * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
+       * s*(*this)+a*V</tt>.
+       */
+      virtual void sadd(const Number s, const Number a,
+                        const VectorSpaceVector<Number> &V);
+
+      /**
+       * Scale each element of this vector by the corresponding element in the
+       * argument. This function is mostly meant to simulate multiplication (and
+       * immediate re-assignment) by a diagonal scaling matrix.
+       */
+      virtual void scale(const VectorSpaceVector<Number> &scaling_factors);
+
+      /**
+       * Assignment <tt>*this = a*V</tt>.
+       */
+      virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
+
+      /**
+       * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
+       * absolute values of all entries among all processors).
+       */
+      virtual real_type l1_norm() const;
+
+      /**
+       * Return the l<sub>2</sub> norm of the vector (i.e., the square root of
+       * the sum of the square of all entries among all processors).
+       */
+      virtual real_type l2_norm() const;
+
+      /**
+       * Return the maximum norm of the vector (i.e., the maximum absolute value
+       * among all entries and among all processors).
+       */
+      virtual real_type linfty_norm() const;
+
+      /**
+       * Perform a combined operation of a vector addition and a subsequent
+       * inner product, returning the value of the inner product. In other
+       * words, the result of this function is the same as if the user called
+       * @code
+       * this->add(a, V);
+       * return_value = *this * W;
+       * @endcode
+       *
+       * The reason this function exists is that this operation involves less
+       * memory transfer than calling the two functions separately. This method
+       * only needs to load three vectors, @p this, @p V, @p W, whereas calling
+       * separate methods means to load the calling vector @p this twice. Since
+       * most vector operations are memory transfer limited, this reduces the
+       * time by 25\% (or 50\% if @p W equals @p this).
+       */
+      virtual Number add_and_dot(const Number a,
+                                 const VectorSpaceVector<Number> &V,
+                                 const VectorSpaceVector<Number> &W);
+
+      /**
+       * Return the global size of the vector, equal to the sum of the number of
+       * locally owned indices among all processors.
+       */
+      virtual size_type size() const;
+
+      /**
+       * Return an index set that describes which elements of this vector are
+       * owned by the current processor. As a consequence, the index sets
+       * returned on different procesors if this is a distributed vector will
+       * form disjoint sets that add up to the complete index set. Obviously, if
+       * a vector is created on only one processor, then the result would
+       * satisfy
+       * @code
+       *  vec.locally_owned_elements() == complete_index_set(vec.size())
+       * @endcode
+       */
+      virtual dealii::IndexSet locally_owned_elements() const;
+
+      /**
+       * Print the vector to the output stream @p out.
+       */
+      virtual void print(std::ostream &out,
+                         const unsigned int precision=3,
+                         const bool scientific=true,
+                         const bool across=true) const;
+
+      /**
+       * Return the memory consumption of this class in bytes.
+       */
+      virtual std::size_t memory_consumption() const;
+      //@}
+
+      /**
+       * @addtogroup Exceptions
+       * @{
+       */
+
+      /**
+       * Attempt to perform an operation between two incompatible vector types.
+       *
+       * @ingroup Exceptions
+       */
+      DeclException0(ExcVectorTypeNotCompatible);
+
+      /**
+       * Exception
+       */
+      DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
+      //@}
+    };
+
+    /*@}*/
+
+  } // end of namespace distributed
+
+} // end of namespace parallel
+
+/**
+ * Global function which overloads the default implementation of the C++
+ * standard library which uses a temporary object. The function simply
+ * exchanges the data of the two vectors.
+ *
+ * @relates BlockVector
+ * @author Katharina Kormann, Martin Kronbichler, 2011
+ */
+template <typename Number>
+inline
+void swap (LinearAlgebra::distributed::BlockVector<Number> &u,
+           LinearAlgebra::distributed::BlockVector<Number> &v)
+{
+  u.swap (v);
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/lac/la_parallel_block_vector.templates.h b/include/deal.II/lac/la_parallel_block_vector.templates.h
new file mode 100644 (file)
index 0000000..a36c5d4
--- /dev/null
@@ -0,0 +1,775 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1999 - 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__parallel_block_vector_templates_h
+#define dealii__parallel_block_vector_templates_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace LinearAlgebra
+{
+  namespace distributed
+  {
+    template <typename Number>
+    BlockVector<Number>::BlockVector (const size_type n_blocks,
+                                      const size_type block_size)
+    {
+      reinit (n_blocks, block_size);
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number>::BlockVector (const std::vector<size_type> &n)
+    {
+      reinit (n, false);
+    }
+
+
+    template <typename Number>
+    BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
+                                      const std::vector<IndexSet> &ghost_indices,
+                                      const MPI_Comm  communicator)
+    {
+      std::vector<size_type> sizes(local_ranges.size());
+      for (unsigned int i=0; i<local_ranges.size(); ++i)
+        sizes[i] = local_ranges[i].size();
+
+      this->block_indices.reinit(sizes);
+      this->components.resize(this->n_blocks());
+
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        this->block(i).reinit(local_ranges[i], ghost_indices[i], communicator);
+    }
+
+
+    template <typename Number>
+    BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
+                                      const MPI_Comm  communicator)
+    {
+      std::vector<size_type> sizes(local_ranges.size());
+      for (unsigned int i=0; i<local_ranges.size(); ++i)
+        sizes[i] = local_ranges[i].size();
+
+      this->block_indices.reinit(sizes);
+      this->components.resize(this->n_blocks());
+
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        this->block(i).reinit(local_ranges[i], communicator);
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number>::BlockVector (const BlockVector<Number> &v)
+      :
+      BlockVectorBase<Vector<Number> > ()
+    {
+      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];
+    }
+
+
+#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
+
+    template <typename Number>
+    template <typename OtherNumber>
+    BlockVector<Number>::BlockVector (const BlockVector<OtherNumber> &v)
+    {
+      reinit (v, true);
+      *this = v;
+    }
+
+#endif
+
+
+
+    template <typename Number>
+    void BlockVector<Number>::reinit (const size_type n_bl,
+                                      const size_type bl_sz,
+                                      const bool         omit_zeroing_entries)
+    {
+      std::vector<size_type> n(n_bl, bl_sz);
+      reinit(n, omit_zeroing_entries);
+    }
+
+
+    template <typename Number>
+    void BlockVector<Number>::reinit (const std::vector<size_type> &n,
+                                      const bool                    omit_zeroing_entries)
+    {
+      this->block_indices.reinit (n);
+      if (this->components.size() != this->n_blocks())
+        this->components.resize(this->n_blocks());
+
+      for (size_type i=0; i<this->n_blocks(); ++i)
+        this->components[i].reinit(n[i], omit_zeroing_entries);
+    }
+
+
+
+    template <typename Number>
+    template <typename Number2>
+    void BlockVector<Number>::reinit (const BlockVector<Number2> &v,
+                                      const bool omit_zeroing_entries)
+    {
+      this->block_indices = v.get_block_indices();
+      if (this->components.size() != this->n_blocks())
+        this->components.resize(this->n_blocks());
+
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        this->block(i).reinit(v.block(i), omit_zeroing_entries);
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number>::~BlockVector ()
+    {}
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator = (const value_type s)
+    {
+      AssertIsFinite(s);
+
+      BaseClass::operator = (s);
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator = (const BlockVector &v)
+    {
+      // we only allow assignment to vectors with the same number of blocks
+      // or to an empty BlockVector
+      Assert (this->n_blocks() == 0 || this->n_blocks() == v.n_blocks(),
+              ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
+
+      if (this->n_blocks() != v.n_blocks())
+        reinit(v.n_blocks(), true);
+
+      for (size_type i=0; i<this->n_blocks(); ++i)
+        this->components[i] = v.block(i);
+
+      this->collect_sizes();
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator = (const Vector<Number> &v)
+    {
+      BaseClass::operator = (v);
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    template <typename Number2>
+    BlockVector<Number> &
+    BlockVector<Number>::operator = (const BlockVector<Number2> &v)
+    {
+      reinit (v, true);
+      BaseClass::operator = (v);
+      return *this;
+    }
+
+
+
+#ifdef DEAL_II_WITH_PETSC
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator = (const PETScWrappers::MPI::BlockVector &petsc_vec)
+    {
+      AssertDimension(this->n_blocks(), petsc_vec.n_blocks());
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        this->block(i) = petsc_vec.block(i);
+
+      return *this;
+    }
+
+#endif
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec)
+    {
+      AssertDimension(this->n_blocks(), trilinos_vec.n_blocks());
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        this->block(i) = trilinos_vec.block(i);
+
+      return *this;
+    }
+
+#endif
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::compress (::dealii::VectorOperation::values operation)
+    {
+      // start all requests for all blocks before finishing the transfers as
+      // this saves repeated synchronizations. In order to avoid conflict with
+      // possible other ongoing communication requests (from
+      // LA::distributed::Vector that supports unfinished requests), add an
+      // arbitrary number 8273 to the communication tag
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).compress_start(block + 8273, operation);
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).compress_finish(operation);
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::update_ghost_values () const
+    {
+      // In order to avoid conflict with possible other ongoing communication
+      // requests (from LA::distributed::Vector that supports unfinished
+      // requests), add an arbitrary number 9923 to the communication tag
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).update_ghost_values_start(block + 9923);
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).update_ghost_values_finish();
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::zero_out_ghosts ()
+    {
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).zero_out_ghosts();
+    }
+
+
+
+    template <typename Number>
+    bool
+    BlockVector<Number>::has_ghost_elements () const
+    {
+      bool has_ghost_elements = false;
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        if (this->block(block).has_ghost_elements() == true)
+          has_ghost_elements = true;
+      return has_ghost_elements;
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator *= (const Number factor)
+    {
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block) *= factor;
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator /= (const Number factor)
+    {
+      operator *= (static_cast<Number>(1.)/factor);
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::scale (const VectorSpaceVector<Number> &vv)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).scale(v.block(block));
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::equ (const Number a,
+                              const VectorSpaceVector<Number> &vv)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).equ(a, v.block(block));
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::equ (const Number a,
+                              const BlockVector<Number> &v,
+                              const Number b,
+                              const BlockVector<Number> &w)
+    {
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      AssertDimension(this->n_blocks(), w.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).equ(a, v.block(block), b, w.block(block));
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator += (const VectorSpaceVector<Number> &vv)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block) += v.block(block);
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    BlockVector<Number> &
+    BlockVector<Number>::operator -= (const VectorSpaceVector<Number> &vv)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block) -= v.block(block);
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::add (const Number a)
+    {
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).add(a);
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::add (const Number a,
+                              const VectorSpaceVector<Number> &vv)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).add(a, v.block(block));
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::add (const Number a,
+                              const VectorSpaceVector<Number> &vv,
+                              const Number b,
+                              const VectorSpaceVector<Number> &ww)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      Assert(dynamic_cast<const BlockVector<Number> *>(&ww)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &w = dynamic_cast<const BlockVector<Number> &>(ww);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).add(a, v.block(block), b, w.block(block));
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::sadd (const Number x,
+                               const Number a,
+                               const VectorSpaceVector<Number> &vv)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).sadd(x, a, v.block(block));
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::sadd (const Number x,
+                               const BlockVector<Number> &v)
+    {
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).sadd(x, v.block(block));
+    }
+
+
+
+    template <typename Number>
+    void
+    BlockVector<Number>::sadd (const Number x,
+                               const Number a,
+                               const BlockVector<Number> &v,
+                               const Number b,
+                               const BlockVector<Number> &w)
+    {
+      AssertDimension(this->n_blocks(), v.n_blocks());
+      AssertDimension(this->n_blocks(), w.n_blocks());
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).sadd(x, a, v.block(block), b, w.block(block));
+    }
+
+
+
+    template <typename Number>
+    template <typename OtherNumber>
+    void
+    BlockVector<Number>::add (const std::vector<size_type>        &indices,
+                              const ::dealii::Vector<OtherNumber> &values)
+    {
+      for (size_type i=0; i<indices.size(); ++i)
+        (*this)(indices[i]) += values[i];
+    }
+
+
+
+    template <typename Number>
+    bool
+    BlockVector<Number>::all_zero () const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      // use int instead of bool. in order to make global reduction operations
+      // work also when MPI_Init was not called, only call MPI_Allreduce
+      // commands when there is more than one processor (note that reinit()
+      // functions handle this case correctly through the job_supports_mpi()
+      // query). this is the same in all the functions below
+      int local_result = -1;
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result = std::max(local_result,
+                                -static_cast<int>(this->block(i).all_zero_local()));
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return -Utilities::MPI::max(local_result,
+                                    this->block(0).partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
+    template <typename Number>
+    Number
+    BlockVector<Number>::operator * (const VectorSpaceVector<Number> &vv) const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+
+      Number local_result = Number();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += this->block(i).inner_product_local(v.block(i));
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum (local_result,
+                                    this->block(0).partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
+    template <typename Number>
+    inline
+    Number
+    BlockVector<Number>::mean_value () const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      Number local_result = Number();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += this->block(i).mean_value_local()*(real_type)this->block(i).partitioner->local_size();
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum (local_result,
+                                    this->block(0).partitioner->get_communicator())/
+               (real_type)this->size();
+      else
+        return local_result/(real_type)this->size();
+    }
+
+
+
+    template <typename Number>
+    inline
+    typename BlockVector<Number>::real_type
+    BlockVector<Number>::l1_norm () const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      real_type local_result = real_type();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += this->block(i).l1_norm_local();
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum (local_result,
+                                    this->block(0).partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
+    template <typename Number>
+    inline
+    typename BlockVector<Number>::real_type
+    BlockVector<Number>::l2_norm () const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      real_type local_result = real_type();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += this->block(i).norm_sqr_local();
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return std::sqrt(Utilities::MPI::sum (local_result,
+                                              this->block(0).partitioner->get_communicator()));
+      else
+        return std::sqrt(local_result);
+    }
+
+
+
+    template <typename Number>
+    inline
+    typename BlockVector<Number>::real_type
+    BlockVector<Number>::lp_norm (const real_type p) const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      real_type local_result = real_type();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += std::pow(this->block(i).lp_norm_local(p), p);
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return std::pow (Utilities::MPI::sum(local_result,
+                                             this->block(0).partitioner->get_communicator()),
+                         static_cast<real_type>(1.0/p));
+      else
+        return std::pow (local_result, static_cast<real_type>(1.0/p));
+    }
+
+
+
+    template <typename Number>
+    inline
+    typename BlockVector<Number>::real_type
+    BlockVector<Number>::linfty_norm () const
+    {
+      Assert (this->n_blocks() > 0, ExcEmptyObject());
+
+      real_type local_result = real_type();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result = std::max(local_result, this->block(i).linfty_norm_local());
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::max (local_result,
+                                    this->block(0).partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
+    template <typename Number>
+    inline
+    Number
+    BlockVector<Number>::add_and_dot (const Number                     a,
+                                      const VectorSpaceVector<Number> &vv,
+                                      const VectorSpaceVector<Number> &ww)
+    {
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&vv)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &v = dynamic_cast<const BlockVector<Number> &>(vv);
+      AssertDimension(this->n_blocks(), v.n_blocks());
+
+      // Downcast. Throws an exception if invalid.
+      Assert(dynamic_cast<const BlockVector<Number> *>(&ww)!=NULL,
+             ExcVectorTypeNotCompatible());
+      const BlockVector<Number> &w = dynamic_cast<const BlockVector<Number> &>(ww);
+      AssertDimension(this->n_blocks(), w.n_blocks());
+
+      Number local_result = Number();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += this->block(i).add_and_dot_local(a, v.block(i), w.block(i));
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum (local_result,
+                                    this->block(0).partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
+    template <typename Number>
+    inline
+    void
+    BlockVector<Number>::swap (BlockVector<Number> &v)
+    {
+      Assert (this->n_blocks() == v.n_blocks(),
+              ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
+
+      for (size_type i=0; i<this->n_blocks(); ++i)
+        dealii::swap (this->components[i], v.components[i]);
+      dealii::swap (this->block_indices, v.block_indices);
+    }
+
+
+
+    template <typename Number>
+    typename BlockVector<Number>::size_type
+    BlockVector<Number>::size () const
+    {
+      return this->block_indices.total_size();
+    }
+
+
+
+    template <typename Number>
+    inline
+    void
+    BlockVector<Number>::import(const LinearAlgebra::ReadWriteVector<Number> &,
+                                VectorOperation::values,
+                                std_cxx11::shared_ptr<const CommunicationPatternBase>)
+    {
+      AssertThrow(false, ExcNotImplemented());
+    }
+
+
+
+    template <typename Number>
+    IndexSet
+    BlockVector<Number>::locally_owned_elements () const
+    {
+      IndexSet is (size());
+
+      // copy index sets from blocks into the global one, shifted by the
+      // appropriate amount for each block
+      for (unsigned int b=0; b<this->n_blocks(); ++b)
+        {
+          IndexSet x = this->block(b).locally_owned_elements();
+          is.add_indices(x, this->block_indices.block_start(b));
+        }
+
+      is.compress();
+
+      return is;
+    }
+
+    template <typename Number>
+    void
+    BlockVector<Number>::print(std::ostream &out,
+                               const unsigned int precision,
+                               const bool scientific,
+                               const bool across) const
+    {
+      for (unsigned int b=0; b<this->n_blocks(); ++b)
+        this->block(b).print(out, precision, scientific, across);
+    }
+
+
+
+    template <typename Number>
+    std::size_t
+    BlockVector<Number>::memory_consumption () const
+    {
+      std::size_t mem = sizeof(this->n_blocks());
+      for (size_type i=0; i<this->components.size(); ++i)
+        mem += MemoryConsumption::memory_consumption (this->components[i]);
+      mem += MemoryConsumption::memory_consumption (this->block_indices);
+      return mem;
+    }
+
+  } // end of namespace distributed
+
+} // end of namespace parallel
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 2cc4e0999dab049497080f805d54dcfd577dde02..52ef2df1e2410040d2f94420c396eb01efabc52b 100644 (file)
@@ -26,7 +26,7 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace parallel
+namespace LinearAlgebra
 {
   namespace distributed
   {
@@ -1105,7 +1105,7 @@ namespace LinearAlgebra
       /**
        * Make BlockVector type friends.
        */
-      template <typename Number2> friend class dealii::parallel::distributed::BlockVector;
+      template <typename Number2> friend class BlockVector;
     };
     /*@}*/
 
index 2d5928786c64e8646f4b1bd5edc37d541ed8314e..c39277621768d5fc27267f06a0402d5a4024b879 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1999 - 2016 by the deal.II authors
+// Copyright (C) 2011 - 2016 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
 #ifndef dealii__parallel_block_vector_h
 #define dealii__parallel_block_vector_h
 
-
 #include <deal.II/base/config.h>
-#include <deal.II/base/exceptions.h>
-#include <deal.II/lac/block_indices.h>
-#include <deal.II/lac/block_vector_base.h>
-#include <deal.II/lac/parallel_vector.h>
-
-#include <deal.II/lac/petsc_parallel_block_vector.h>
-#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
 
+#include <cstring>
+#include <iomanip>
 
-#include <cstdio>
-#include <vector>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -37,12 +30,10 @@ namespace parallel
 {
   namespace distributed
   {
-
     /*! @addtogroup Vectors
      *@{
      */
 
-
     /**
      * An implementation of block vectors based on distributed deal.II
      * vectors. While the base class provides for most of the interface, this
@@ -59,952 +50,10 @@ namespace parallel
      * @ref GlossBlockLA "Block (linear algebra)"
      * @author Katharina Kormann, Martin Kronbichler, 2011
      */
-    template <typename Number>
-    class BlockVector : public BlockVectorBase<Vector<Number> >
-    {
-    public:
-      /**
-       * Typedef the base class for simpler access to its own typedefs.
-       */
-      typedef BlockVectorBase<Vector<Number> > BaseClass;
-
-      /**
-       * Typedef the type of the underlying vector.
-       */
-      typedef typename BaseClass::BlockType  BlockType;
-
-      /**
-       * Import the typedefs from the base class.
-       */
-      typedef typename BaseClass::value_type      value_type;
-      typedef typename BaseClass::real_type       real_type;
-      typedef typename BaseClass::pointer         pointer;
-      typedef typename BaseClass::const_pointer   const_pointer;
-      typedef typename BaseClass::reference       reference;
-      typedef typename BaseClass::const_reference const_reference;
-      typedef typename BaseClass::size_type       size_type;
-      typedef typename BaseClass::iterator        iterator;
-      typedef typename BaseClass::const_iterator  const_iterator;
-
-      /**
-       * Constructor. There are three ways to use this constructor. First,
-       * without any arguments, it generates an object with no blocks. Given
-       * one argument, it initializes <tt>num_blocks</tt> blocks, but these
-       * blocks have size zero. The third variant finally initializes all
-       * blocks to the same size <tt>block_size</tt>.
-       *
-       * Confer the other constructor further down if you intend to use blocks
-       * of different sizes.
-       */
-      explicit BlockVector (const size_type num_blocks = 0,
-                            const size_type block_size = 0);
-
-      /**
-       * Copy-Constructor. Dimension set to that of V, all components are
-       * copied from V
-       */
-      BlockVector (const BlockVector<Number> &V);
-
-
-#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
-      /**
-       * Copy constructor taking a BlockVector of another data type. This will
-       * fail if there is no conversion path from <tt>OtherNumber</tt> to
-       * <tt>Number</tt>. Note that you may lose accuracy when copying to a
-       * BlockVector with data elements with less accuracy.
-       *
-       * Older versions of gcc did not honor the @p explicit keyword on
-       * template constructors. In such cases, it is easy to accidentally
-       * write code that can be very inefficient, since the compiler starts
-       * performing hidden conversions. To avoid this, this function is
-       * disabled if we have detected a broken compiler during configuration.
-       */
-      template <typename OtherNumber>
-      explicit
-      BlockVector (const BlockVector<OtherNumber> &v);
-#endif
-
-      /**
-       * Constructor. Set the number of blocks to <tt>block_sizes.size()</tt>
-       * and initialize each block with <tt>block_sizes[i]</tt> zero elements.
-       */
-      BlockVector (const std::vector<size_type> &block_sizes);
-
-      /**
-       * Construct a block vector with an IndexSet for the local range and
-       * ghost entries for each block.
-       */
-      BlockVector (const std::vector<IndexSet> &local_ranges,
-                   const std::vector<IndexSet> &ghost_indices,
-                   const MPI_Comm  communicator);
-
-      /**
-       * Same as above but the ghost indices are assumed to be empty.
-       */
-      BlockVector (const std::vector<IndexSet> &local_ranges,
-                   const MPI_Comm  communicator);
-
-      /**
-       * Destructor. Clears memory.
-       */
-      ~BlockVector ();
-
-      /**
-       * Copy operator: fill all components of the vector with the given
-       * scalar value.
-       */
-      BlockVector &operator = (const value_type s);
-
-      /**
-       * Copy operator for arguments of the same type. Resize the present
-       * vector if necessary.
-       */
-      BlockVector &
-      operator= (const BlockVector &V);
-
-      /**
-       * Copy operator for template arguments of different types. Resize the
-       * present vector if necessary.
-       */
-      template <class Number2>
-      BlockVector &
-      operator= (const BlockVector<Number2> &V);
-
-      /**
-       * Copy a regular vector into a block vector.
-       */
-      BlockVector &
-      operator= (const Vector<Number> &V);
-
-#ifdef DEAL_II_WITH_PETSC
-      /**
-       * Copy the content of a PETSc vector into the calling vector. This
-       * function assumes that the vectors layouts have already been
-       * initialized to match.
-       *
-       * This operator is only available if deal.II was configured with PETSc.
-       */
-      BlockVector<Number> &
-      operator = (const PETScWrappers::MPI::BlockVector &petsc_vec);
-#endif
-
-#ifdef DEAL_II_WITH_TRILINOS
-      /**
-       * Copy the content of a Trilinos vector into the calling vector. This
-       * function assumes that the vectors layouts have already been
-       * initialized to match.
-       *
-       * This operator is only available if deal.II was configured with
-       * Trilinos.
-       */
-      BlockVector<Number> &
-      operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec);
-#endif
-
-      /**
-       * Reinitialize the BlockVector to contain <tt>num_blocks</tt> blocks of
-       * size <tt>block_size</tt> each.
-       *
-       * If the second argument is left at its default value, then the block
-       * vector allocates the specified number of blocks but leaves them at
-       * zero size. You then need to later reinitialize the individual blocks,
-       * and call collect_sizes() to update the block system's knowledge of
-       * its individual block's sizes.
-       *
-       * If <tt>omit_zeroing_entries==false</tt>, the vector is filled with
-       * zeros.
-       */
-      void reinit (const size_type num_blocks,
-                   const size_type block_size = 0,
-                   const bool omit_zeroing_entries = false);
-
-      /**
-       * Reinitialize the BlockVector such that it contains
-       * <tt>block_sizes.size()</tt> blocks. Each block is reinitialized to
-       * dimension <tt>block_sizes[i]</tt>.
-       *
-       * If the number of blocks is the same as before this function was
-       * called, all vectors remain the same and reinit() is called for each
-       * vector.
-       *
-       * If <tt>omit_zeroing_entries==false</tt>, the vector is filled with
-       * zeros.
-       *
-       * Note that you must call this (or the other reinit() functions)
-       * function, rather than calling the reinit() functions of an individual
-       * block, to allow the block vector to update its caches of vector
-       * sizes. If you call reinit() on one of the blocks, then subsequent
-       * actions on this object may yield unpredictable results since they may
-       * be routed to the wrong block.
-       */
-      void reinit (const std::vector<size_type> &N,
-                   const bool                    omit_zeroing_entries=false);
-
-      /**
-       * Change the dimension to that of the vector <tt>V</tt>. The same
-       * applies as for the other reinit() function.
-       *
-       * The elements of <tt>V</tt> are not copied, i.e.  this function is the
-       * same as calling <tt>reinit (V.size(), omit_zeroing_entries)</tt>.
-       *
-       * Note that you must call this (or the other reinit() functions)
-       * function, rather than calling the reinit() functions of an individual
-       * block, to allow the block vector to update its caches of vector
-       * sizes. If you call reinit() of one of the blocks, then subsequent
-       * actions of this object may yield unpredictable results since they may
-       * be routed to the wrong block.
-       */
-      template <typename Number2>
-      void reinit (const BlockVector<Number2> &V,
-                   const bool                 omit_zeroing_entries=false);
-
-      /**
-       * This function copies the data that has accumulated in the data buffer
-       * for ghost indices to the owning processor. For the meaning of the
-       * argument @p operation, see the entry on
-       * @ref GlossCompress "Compressing distributed vectors and matrices"
-       * in the glossary.
-       *
-       * There are two variants for this function. If called with argument @p
-       * VectorOperation::add adds all the data accumulated in ghost elements
-       * to the respective elements on the owning processor and clears the
-       * ghost array afterwards. If called with argument @p
-       * VectorOperation::insert, a set operation is performed. Since setting
-       * elements in a vector with ghost elements is ambiguous (as one can set
-       * both the element on the ghost site as well as the owning site), this
-       * operation makes the assumption that all data is set correctly on the
-       * owning processor. Upon call of compress(VectorOperation::insert), all
-       * ghost entries are therefore simply zeroed out (using
-       * zero_ghost_values()). In debug mode, a check is performed that makes
-       * sure that the data set is actually consistent between processors,
-       * i.e., whenever a non-zero ghost element is found, it is compared to
-       * the value on the owning processor and an exception is thrown if these
-       * elements do not agree.
-       *
-       */
-      void compress (::dealii::VectorOperation::values operation);
-
-      /**
-       * Fills the data field for ghost indices with the values stored in the
-       * respective positions of the owning processor. This function is needed
-       * before reading from ghosts. The function is @p const even though
-       * ghost data is changed. This is needed to allow functions with a @p
-       * const vector to perform the data exchange without creating
-       * temporaries.
-       */
-      void update_ghost_values () const;
-
-      /**
-       * This method zeros the entries on ghost dofs, but does not touch
-       * locally owned DoFs.
-       *
-       * After calling this method, read access to ghost elements of the
-       * vector is forbidden and an exception is thrown. Only write access to
-       * ghost elements is allowed in this state.
-       */
-      void zero_out_ghosts ();
-
-      /**
-       * Returns if this Vector contains ghost elements.
-       */
-      bool has_ghost_elements() const;
-
-      /**
-       * Return whether the vector contains only elements with value zero.
-       * This function is mainly for internal consistency checks and should
-       * seldom be used when not in debug mode since it uses quite some time.
-       */
-      bool all_zero () const;
-
-      /**
-       * Return @p true if the vector has no negative entries, i.e. all
-       * entries are zero or positive. This function is used, for example, to
-       * check whether refinement indicators are really all positive (or
-       * zero).
-       *
-       * The function obviously only makes sense if the template argument of
-       * this class is a real type. If it is a complex type, then an exception
-       * is thrown.
-       */
-      bool is_non_negative () const;
-
-      /**
-       * Checks for equality of the two vectors.
-       */
-      template <typename Number2>
-      bool operator == (const BlockVector<Number2> &v) const;
-
-      /**
-       * Checks for inequality of the two vectors.
-       */
-      template <typename Number2>
-      bool operator != (const BlockVector<Number2> &v) const;
-
-      /**
-       * Perform the inner product of two vectors.
-       */
-      template <typename Number2>
-      Number operator * (const BlockVector<Number2> &V) const;
-
-      /**
-       * Computes the square of the l<sub>2</sub> norm of the vector (i.e.,
-       * the sum of the squares of all entries among all processors).
-       */
-      real_type norm_sqr () const;
-
-      /**
-       * Computes the mean value of all the entries in the vector.
-       */
-      Number mean_value () const;
-
-      /**
-       * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
-       * absolute values of all entries among all processors).
-       */
-      real_type l1_norm () const;
-
-      /**
-       * Returns the l<sub>2</sub> norm of the vector (i.e., square root of
-       * the sum of the square of all entries among all processors).
-       */
-      real_type l2_norm () const;
-
-      /**
-       * Returns the l<sub>p</sub> norm with real @p p of the vector (i.e.,
-       * the pth root of sum of the pth power of all entries among all
-       * processors).
-       */
-      real_type lp_norm (const real_type p) const;
-
-      /**
-       * Returns the maximum norm of the vector (i.e., maximum absolute value
-       * among all entries among all processors).
-       */
-      real_type linfty_norm () const;
-
-      /**
-       * Performs a combined operation of a vector addition and a subsequent
-       * inner product, returning the value of the inner product. In other
-       * words, the result of this function is the same as if the user called
-       * @code
-       * this->add(a, V);
-       * return_value = *this * W;
-       * @endcode
-       *
-       * The reason this function exists is that this operation involves less
-       * memory transfer than calling the two functions separately. This
-       * method only needs to load three vectors, @p this, @p V, @p W, whereas
-       * calling separate methods means to load the calling vector @p this
-       * twice. Since most vector operations are memory transfer limited, this
-       * reduces the time by 25\% (or 50\% if @p W equals @p this).
-       */
-      Number add_and_dot (const Number               a,
-                          const BlockVector<Number> &V,
-                          const BlockVector<Number> &W);
-
-      /**
-       * Multiply each element of this vector by the corresponding element of
-       * <tt>v</tt>.
-       */
-      template <class BlockVector2>
-      void scale (const BlockVector2 &v);
-
-      /**
-       * Swap the contents of this vector and the other vector <tt>v</tt>. One
-       * could do this operation with a temporary variable and copying over
-       * the data elements, but this function is significantly more efficient
-       * since it only swaps the pointers to the data of the two vectors and
-       * therefore does not need to allocate temporary storage and move data
-       * around.
-       *
-       * Limitation: right now this function only works if both vectors have
-       * the same number of blocks. If needed, the numbers of blocks should be
-       * exchanged, too.
-       *
-       * This function is analog to the the swap() function of all C++
-       * standard containers. Also, there is a global function swap(u,v) that
-       * simply calls <tt>u.swap(v)</tt>, again in analogy to standard
-       * functions.
-       */
-      void swap (BlockVector<Number> &v);
-
-      /**
-       * @addtogroup Exceptions
-       * @{
-       */
-
-      /**
-       * Exception
-       */
-      DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
-      //@}
-    };
+    using LinearAlgebra::distributed::BlockVector;
 
     /*@}*/
-
-#ifndef DOXYGEN
-    /*----------------------- Inline functions ----------------------------------*/
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number>::BlockVector (const size_type n_blocks,
-                                      const size_type block_size)
-    {
-      reinit (n_blocks, block_size);
-    }
-
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number>::BlockVector (const std::vector<size_type> &n)
-    {
-      reinit (n, false);
-    }
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
-                                      const std::vector<IndexSet> &ghost_indices,
-                                      const MPI_Comm  communicator)
-    {
-      std::vector<size_type> sizes(local_ranges.size());
-      for (unsigned int i=0; i<local_ranges.size(); ++i)
-        sizes[i] = local_ranges[i].size();
-
-      this->block_indices.reinit(sizes);
-      this->components.resize(this->n_blocks());
-
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        this->block(i).reinit(local_ranges[i], ghost_indices[i], communicator);
-    }
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
-                                      const MPI_Comm  communicator)
-    {
-      std::vector<size_type> sizes(local_ranges.size());
-      for (unsigned int i=0; i<local_ranges.size(); ++i)
-        sizes[i] = local_ranges[i].size();
-
-      this->block_indices.reinit(sizes);
-      this->components.resize(this->n_blocks());
-
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        this->block(i).reinit(local_ranges[i], communicator);
-    }
-
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number>::BlockVector (const BlockVector<Number> &v)
-      :
-      BlockVectorBase<Vector<Number> > ()
-    {
-      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];
-    }
-
-
-#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
-
-    template <typename Number>
-    template <typename OtherNumber>
-    inline
-    BlockVector<Number>::BlockVector (const BlockVector<OtherNumber> &v)
-    {
-      reinit (v, true);
-      *this = v;
-    }
-
-#endif
-
-
-
-    template <typename Number>
-    inline
-    void BlockVector<Number>::reinit (const size_type n_bl,
-                                      const size_type bl_sz,
-                                      const bool         omit_zeroing_entries)
-    {
-      std::vector<size_type> n(n_bl, bl_sz);
-      reinit(n, omit_zeroing_entries);
-    }
-
-
-    template <typename Number>
-    inline
-    void BlockVector<Number>::reinit (const std::vector<size_type> &n,
-                                      const bool                    omit_zeroing_entries)
-    {
-      this->block_indices.reinit (n);
-      if (this->components.size() != this->n_blocks())
-        this->components.resize(this->n_blocks());
-
-      for (size_type i=0; i<this->n_blocks(); ++i)
-        this->components[i].reinit(n[i], omit_zeroing_entries);
-    }
-
-
-
-    template <typename Number>
-    template <typename Number2>
-    inline
-    void BlockVector<Number>::reinit (const BlockVector<Number2> &v,
-                                      const bool omit_zeroing_entries)
-    {
-      this->block_indices = v.get_block_indices();
-      if (this->components.size() != this->n_blocks())
-        this->components.resize(this->n_blocks());
-
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        this->block(i).reinit(v.block(i), omit_zeroing_entries);
-    }
-
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number>::~BlockVector ()
-    {}
-
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number> &
-    BlockVector<Number>::operator = (const value_type s)
-    {
-
-      AssertIsFinite(s);
-
-      BaseClass::operator = (s);
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number> &
-    BlockVector<Number>::operator = (const BlockVector &v)
-    {
-      // we only allow assignment to vectors with the same number of blocks
-      // or to an empty BlockVector
-      Assert (this->n_blocks() == 0 || this->n_blocks() == v.n_blocks(),
-              ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
-
-      if (this->n_blocks() != v.n_blocks())
-        reinit(v.n_blocks(), true);
-
-      for (size_type i=0; i<this->n_blocks(); ++i)
-        this->components[i] = v.block(i);
-
-      this->collect_sizes();
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    inline
-    BlockVector<Number> &
-    BlockVector<Number>::operator = (const Vector<Number> &v)
-    {
-      BaseClass::operator = (v);
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    template <typename Number2>
-    inline
-    BlockVector<Number> &
-    BlockVector<Number>::operator = (const BlockVector<Number2> &v)
-    {
-      reinit (v, true);
-      BaseClass::operator = (v);
-      return *this;
-    }
-
-
-
-#ifdef DEAL_II_WITH_PETSC
-
-    template <typename Number>
-    inline
-    BlockVector<Number> &
-    BlockVector<Number>::operator = (const PETScWrappers::MPI::BlockVector &petsc_vec)
-    {
-      AssertDimension(this->n_blocks(), petsc_vec.n_blocks());
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        this->block(i) = petsc_vec.block(i);
-
-      return *this;
-    }
-
-#endif
-
-
-
-#ifdef DEAL_II_WITH_TRILINOS
-
-    template <typename Number>
-    inline
-    BlockVector<Number> &
-    BlockVector<Number>::operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec)
-    {
-      AssertDimension(this->n_blocks(), trilinos_vec.n_blocks());
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        this->block(i) = trilinos_vec.block(i);
-
-      return *this;
-    }
-
-#endif
-
-
-
-    template <typename Number>
-    inline
-    void
-    BlockVector<Number>::compress (::dealii::VectorOperation::values operation)
-    {
-      // start all requests for all blocks before finishing the transfers as
-      // this saves repeated synchronizations
-      for (unsigned int block=0; block<this->n_blocks(); ++block)
-        this->block(block).compress_start(block*10 + 8273, operation);
-      for (unsigned int block=0; block<this->n_blocks(); ++block)
-        this->block(block).compress_finish(operation);
-    }
-
-
-
-    template <typename Number>
-    inline
-    void
-    BlockVector<Number>::update_ghost_values () const
-    {
-      for (unsigned int block=0; block<this->n_blocks(); ++block)
-        this->block(block).update_ghost_values_start(block*10 + 9923);
-      for (unsigned int block=0; block<this->n_blocks(); ++block)
-        this->block(block).update_ghost_values_finish();
-    }
-
-
-
-    template <typename Number>
-    inline
-    void
-    BlockVector<Number>::zero_out_ghosts ()
-    {
-      for (unsigned int block=0; block<this->n_blocks(); ++block)
-        this->block(block).zero_out_ghosts();
-    }
-
-
-
-    template <typename Number>
-    inline
-    bool
-    BlockVector<Number>::has_ghost_elements () const
-    {
-      bool has_ghost_elements = false;
-      for (unsigned int block=0; block<this->n_blocks(); ++block)
-        if (this->block(block).has_ghost_elements() == true)
-          has_ghost_elements = true;
-      return has_ghost_elements;
-    }
-
-
-
-    template <typename Number>
-    inline
-    bool
-    BlockVector<Number>::all_zero () const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-
-      // use int instead of bool. in order to make global reduction operations
-      // work also when MPI_Init was not called, only call MPI_Allreduce
-      // commands when there is more than one processor (note that reinit()
-      // functions handle this case correctly through the job_supports_mpi()
-      // query). this is the same in all the functions below
-      int local_result = -1;
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result = std::max(local_result,
-                                -static_cast<int>(this->block(i).all_zero_local()));
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return -Utilities::MPI::max(local_result,
-                                    this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    inline
-    bool
-    BlockVector<Number>::is_non_negative () const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-      int local_result = -1;
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result = std::max(local_result,
-                                -static_cast<int>(this->block(i).is_non_negative_local()));
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::max(local_result,
-                                   this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    template <typename Number2>
-    inline
-    bool
-    BlockVector<Number>::operator == (const BlockVector<Number2> &v) const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-      AssertDimension (this->n_blocks(), v.n_blocks());
-
-      // MPI does not support bools, so use unsigned int instead. Two vectors
-      // are equal if the check for non-equal fails on all processors
-      unsigned int local_result = 0;
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result = std::max(local_result,
-                                static_cast<unsigned int>(!this->block(i).vectors_equal_local(v.block(i))));
-      unsigned int result =
-        this->block(0).partitioner->n_mpi_processes() > 1
-        ?
-        Utilities::MPI::max(local_result, this->block(0).partitioner->get_communicator())
-        :
-        local_result;
-      return result==0;
-    }
-
-
-
-    template <typename Number>
-    template <typename Number2>
-    inline
-    bool
-    BlockVector<Number>::operator != (const BlockVector<Number2> &v) const
-    {
-      return !(operator == (v));
-    }
-
-
-
-    template <typename Number>
-    template <typename Number2>
-    inline
-    Number
-    BlockVector<Number>::operator * (const BlockVector<Number2> &v) const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-      AssertDimension (this->n_blocks(), v.n_blocks());
-
-      Number local_result = Number();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result += this->block(i).inner_product_local(v.block(i));
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::sum (local_result,
-                                    this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    inline
-    typename BlockVector<Number>::real_type
-    BlockVector<Number>::norm_sqr () const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-
-      real_type local_result = real_type();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result += this->block(i).norm_sqr_local();
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::sum (local_result,
-                                    this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    inline
-    Number
-    BlockVector<Number>::mean_value () const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-
-      Number local_result = Number();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result += this->block(i).mean_value_local()*(real_type)this->block(i).partitioner->local_size();
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::sum (local_result,
-                                    this->block(0).partitioner->get_communicator())/
-               (real_type)this->size();
-      else
-        return local_result/(real_type)this->size();
-    }
-
-
-
-    template <typename Number>
-    inline
-    typename BlockVector<Number>::real_type
-    BlockVector<Number>::l1_norm () const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-
-      real_type local_result = real_type();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result += this->block(i).l1_norm_local();
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::sum (local_result,
-                                    this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    inline
-    typename BlockVector<Number>::real_type
-    BlockVector<Number>::l2_norm () const
-    {
-      return std::sqrt(norm_sqr());
-    }
-
-
-
-    template <typename Number>
-    inline
-    typename BlockVector<Number>::real_type
-    BlockVector<Number>::lp_norm (const real_type p) const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-
-      real_type local_result = real_type();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result += std::pow(this->block(i).lp_norm_local(p), p);
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return std::pow (Utilities::MPI::sum(local_result,
-                                             this->block(0).partitioner->get_communicator()),
-                         static_cast<real_type>(1.0/p));
-      else
-        return std::pow (local_result, static_cast<real_type>(1.0/p));
-    }
-
-
-
-    template <typename Number>
-    inline
-    typename BlockVector<Number>::real_type
-    BlockVector<Number>::linfty_norm () const
-    {
-      Assert (this->n_blocks() > 0, ExcEmptyObject());
-
-      real_type local_result = real_type();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result = std::max(local_result, this->block(i).linfty_norm_local());
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::max (local_result,
-                                    this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    inline
-    Number
-    BlockVector<Number>::add_and_dot (const Number               a,
-                                      const BlockVector<Number> &V,
-                                      const BlockVector<Number> &W)
-    {
-      Number local_result = Number();
-      for (unsigned int i=0; i<this->n_blocks(); ++i)
-        local_result += this->block(i).add_and_dot_local(a, V.block(i), W.block(i));
-
-      if (this->block(0).partitioner->n_mpi_processes() > 1)
-        return Utilities::MPI::sum (local_result,
-                                    this->block(0).partitioner->get_communicator());
-      else
-        return local_result;
-    }
-
-
-
-    template <typename Number>
-    inline
-    void
-    BlockVector<Number>::swap (BlockVector<Number> &v)
-    {
-      Assert (this->n_blocks() == v.n_blocks(),
-              ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
-
-      for (size_type i=0; i<this->n_blocks(); ++i)
-        dealii::swap (this->components[i], v.components[i]);
-      dealii::swap (this->block_indices, v.block_indices);
-    }
-
-
-
-    template <typename Number>
-    template <class BlockVector2>
-    void BlockVector<Number>::scale (const BlockVector2 &v)
-    {
-      BaseClass::scale (v);
-    }
-
-#endif // DOXYGEN
-
-  } // end of namespace distributed
-
-} // end of namespace parallel
-
-/**
- * Global function which overloads the default implementation of the C++
- * standard library which uses a temporary object. The function simply
- * exchanges the data of the two vectors.
- *
- * @relates BlockVector
- * @author Katharina Kormann, Martin Kronbichler, 2011
- */
-template <typename Number>
-inline
-void swap (parallel::distributed::BlockVector<Number> &u,
-           parallel::distributed::BlockVector<Number> &v)
-{
-  u.swap (v);
+  }
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 146e35085d84f02a47edc8c640d37ed7b13fcbe9..b03b0f612146e357cd0382e7c6d907c57c07484f 100644 (file)
@@ -30,8 +30,6 @@ namespace parallel
 {
   namespace distributed
   {
-    template <typename Number> class BlockVector;
-
     /*! @addtogroup Vectors
      *@{
      */
@@ -140,6 +138,8 @@ namespace parallel
      * @author Katharina Kormann, Martin Kronbichler, 2010, 2011
      */
     using LinearAlgebra::distributed::Vector;
+
+    /*@}*/
   }
 }
 
index 9ead26e3fef7a520500424a64f8a896f2b4f37b6..924b993baf9b84f1168290ec1510ec54dc739314 100644 (file)
@@ -29,6 +29,7 @@ SET(_src
   lapack_full_matrix.cc
   la_vector.cc
   la_parallel_vector.cc
+  la_parallel_block_vector.cc
   matrix_lib.cc
   matrix_out.cc
   precondition_block.cc
@@ -63,6 +64,7 @@ SET(_inst
   lapack_full_matrix.inst.in
   la_vector.inst.in
   la_parallel_vector.inst.in
+  la_parallel_block_vector.inst.in
   precondition_block.inst.in
   relaxation_block.inst.in
   read_write_vector.inst.in
diff --git a/source/lac/la_parallel_block_vector.cc b/source/lac/la_parallel_block_vector.cc
new file mode 100644 (file)
index 0000000..bbabfc3
--- /dev/null
@@ -0,0 +1,46 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#include "la_parallel_block_vector.inst"
+
+// do a few functions that currently don't fit the scheme because they have
+// two template arguments that need to be different (the case of same
+// arguments is covered by the default copy constructor and copy operator that
+// is declared separately)
+
+namespace LinearAlgebra
+{
+  namespace distributed
+  {
+#define TEMPL_COPY_CONSTRUCTOR(S1,S2)                                   \
+  template BlockVector<S1>& BlockVector<S1>::operator=<S2> (const BlockVector<S2> &)
+
+#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
+    TEMPL_COPY_CONSTRUCTOR(double,float);
+    TEMPL_COPY_CONSTRUCTOR(float,double);
+
+#endif
+
+#undef TEMPL_COPY_CONSTRUCTOR
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/lac/la_parallel_block_vector.inst.in b/source/lac/la_parallel_block_vector.inst.in
new file mode 100644 (file)
index 0000000..adcb866
--- /dev/null
@@ -0,0 +1,41 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2011 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (SCALAR : REAL_SCALARS)
+{
+  namespace LinearAlgebra
+  \{
+    namespace distributed
+    \{
+      template class BlockVector<SCALAR>;
+    \}
+  \}
+}
+
+for (S1, S2 : REAL_SCALARS)
+{
+  namespace LinearAlgebra
+  \{
+    namespace distributed
+    \{
+      template void BlockVector<S1>::reinit<S2> (const BlockVector<S2>&,
+                                                 const bool);
+      template void BlockVector<S1>::add<S2> (const std::vector<size_type> &,
+                                              const ::dealii::Vector<S2>&);
+    \}
+  \}
+}
index e06292f4e9baea3c1f001ffff9d80fd0027fd807..9b118a8acfb30aef2cf313979f6546b81d8addb9 100644 (file)
@@ -34,6 +34,7 @@ for (S1, S2 : REAL_SCALARS)
     \{
       template void Vector<S1>::reinit<S2> (const Vector<S2>&,
                                             const bool);
+      template S1 Vector<S1>::inner_product_local<S2> (const Vector<S2>&) const;
     \}
   \}
 }
index ee7e9045ae84ba9584480464c79decac61b4f30c..d84ff994563069cf7062006b8507fe6389507adf 100644 (file)
@@ -45,7 +45,7 @@ namespace LinearAlgebra
       comm = std_cxx11::make_shared<const MPI_Comm>(communicator);
 
       Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm,
-                                           true);
+                                           false);
       Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm,
                                          true);
 
index a9cefeee1f244677b0e9e895c287ab91efab464a..e76fda085d682df85c4639953562888fd6ac249e 100644 (file)
@@ -118,12 +118,12 @@ void test ()
   }
   // check inner product
   {
-    const double norm_sqr = w.norm_sqr();
-    AssertThrow (std::fabs(w * w - norm_sqr) < 1e-15,
+    const double norm_sqr = w.l2_norm() * w.l2_norm();
+    AssertThrow (std::fabs(w * w - norm_sqr) < 1e-12,
                  ExcInternalError());
     parallel::distributed::BlockVector<double> w2;
     w2 = w;
-    AssertThrow (std::fabs(w2 * w - norm_sqr) < 1e-15,
+    AssertThrow (std::fabs(w2 * w - norm_sqr) < 1e-12,
                  ExcInternalError());
 
     if (myid<8)
@@ -133,48 +133,6 @@ void test ()
       deallog << "Inner product: " << inner_prod << std::endl;
   }
 
-  // check operator ==
-  {
-    parallel::distributed::BlockVector<double> w2 (w);
-    bool equal = (w2 == w);
-    if (myid == 0)
-      deallog << " v==v2 ? " << equal << std::endl;
-
-    bool not_equal = (w2 != w);
-    if (myid == 0)
-      deallog << " v!=v2 ? " << not_equal << std::endl;
-
-    // change v2 on one proc only
-    if (myid == 0)
-      w2.block(0).local_element(1) = 2.2212;
-
-    equal = (w2 == w);
-    if (myid == 0)
-      deallog << " v==v2 ? " << equal << std::endl;
-    not_equal = (w2 != w);
-    if (myid == 0)
-      deallog << " v!=v2 ? " << not_equal << std::endl;
-
-    // reset
-    w2 = w;
-    equal = (w2 == w);
-    if (myid == 0)
-      deallog << " v==v2 ? " << equal << std::endl;
-    not_equal = (w2 != w);
-    if (myid == 0)
-      deallog << " v!=v2 ? " << not_equal << std::endl;
-
-    // change some value on all procs
-    if (myid < 8)
-      w2.block(1).local_element(0) = -1;
-    equal = (w2 == w);
-    if (myid == 0)
-      deallog << " v==v2 ? " << equal << std::endl;
-    not_equal = (w2 != w);
-    if (myid == 0)
-      deallog << " v!=v2 ? " << not_equal << std::endl;
-  }
-
   // check all_zero
   {
     bool allzero = w.all_zero();
@@ -194,38 +152,6 @@ void test ()
       deallog << " v2==0 ? " << allzero << std::endl;
   }
 
-
-  // check all_non_negative
-  {
-    bool allnonneg = w.is_non_negative();
-    if (myid == 0)
-      deallog << " v>=0 ? " << allnonneg << std::endl;
-    parallel::distributed::BlockVector<double> w2, w3;
-
-    // vector where all processors have
-    // non-negative entries
-    w2 = w;
-    if (myid < 8)
-      w2.block(0).local_element(0) = -1;
-    allnonneg = w2.is_non_negative();
-    if (myid == 0)
-      deallog << " v2>=0 ? " << allnonneg << std::endl;
-
-    // zero vector
-    w3.reinit (w2);
-    allnonneg = w3.is_non_negative();
-    if (myid == 0)
-      deallog << " v3>=0 ? " << allnonneg << std::endl;
-
-    // only one processor has non-negative entry
-    w3 = w;
-    if (myid == 1 || numproc==1)
-      w3.block(0).local_element(0) = -1;
-    allnonneg = w3.is_non_negative();
-    if (myid == 0)
-      deallog << " v3>=0 ? " << allnonneg << std::endl;
-  }
-
   if (myid == 0)
     deallog << "OK" << std::endl;
 }
index 2b1e9bf7df7cf2094205fc3b5a0242101b7318d4..3a3e5186c73c3e2ce1a56fdb5d32ce6e023f3aa0 100644 (file)
@@ -6,19 +6,7 @@ DEAL:0::linfty norm: 2.000
 DEAL:0::l2.2 norm: 3.295
 DEAL:0::Mean value: 1.000
 DEAL:0::Inner product: 12.00
-DEAL:0:: v==v2 ? 1
-DEAL:0:: v!=v2 ? 0
-DEAL:0:: v==v2 ? 0
-DEAL:0:: v!=v2 ? 1
-DEAL:0:: v==v2 ? 1
-DEAL:0:: v!=v2 ? 0
-DEAL:0:: v==v2 ? 0
-DEAL:0:: v!=v2 ? 1
 DEAL:0:: v==0 ? 0
 DEAL:0:: v2==0 ? 1
 DEAL:0:: v2==0 ? 0
-DEAL:0:: v>=0 ? 1
-DEAL:0:: v2>=0 ? 0
-DEAL:0:: v3>=0 ? 1
-DEAL:0:: v3>=0 ? 0
 DEAL:0::OK
index 53261a12e39d9599cd796858500df1e2d051dd73..ba129a4d0397c7a84bd101e5ace08cb636d83a25 100644 (file)
@@ -6,19 +6,7 @@ DEAL:0::linfty norm: 30.00
 DEAL:0::l2.2 norm: 104.6
 DEAL:0::Mean value: 15.00
 DEAL:0::Inner product: 1.253e+04
-DEAL:0:: v==v2 ? 1
-DEAL:0:: v!=v2 ? 0
-DEAL:0:: v==v2 ? 0
-DEAL:0:: v!=v2 ? 1
-DEAL:0:: v==v2 ? 1
-DEAL:0:: v!=v2 ? 0
-DEAL:0:: v==v2 ? 0
-DEAL:0:: v!=v2 ? 1
 DEAL:0:: v==0 ? 0
 DEAL:0:: v2==0 ? 1
 DEAL:0:: v2==0 ? 0
-DEAL:0:: v>=0 ? 1
-DEAL:0:: v2>=0 ? 0
-DEAL:0:: v3>=0 ? 1
-DEAL:0:: v3>=0 ? 0
 DEAL:0::OK
index 8fd4464a522ab0fd30d6d62f14278e87ac8745b8..9ffa1364ef6190066c571388a289c1bd6eaf7eea 100644 (file)
@@ -6,19 +6,7 @@ DEAL:0::linfty norm: 14.00
 DEAL:0::l2.2 norm: 36.31
 DEAL:0::Mean value: 7.000
 DEAL:0::Inner product: 1432.
-DEAL:0:: v==v2 ? 1
-DEAL:0:: v!=v2 ? 0
-DEAL:0:: v==v2 ? 0
-DEAL:0:: v!=v2 ? 1
-DEAL:0:: v==v2 ? 1
-DEAL:0:: v!=v2 ? 0
-DEAL:0:: v==v2 ? 0
-DEAL:0:: v!=v2 ? 1
 DEAL:0:: v==0 ? 0
 DEAL:0:: v2==0 ? 1
 DEAL:0:: v2==0 ? 0
-DEAL:0:: v>=0 ? 1
-DEAL:0:: v2>=0 ? 0
-DEAL:0:: v3>=0 ? 1
-DEAL:0:: v3>=0 ? 0
 DEAL:0::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.