From 7d00ab3caff913b72ad647ef72176167b0b042ea Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 14 Jun 2016 18:54:32 +0200 Subject: [PATCH] Add LinearAlgebra::distributed::BlockVector --- .../deal.II/lac/la_parallel_block_vector.h | 599 +++++++++++ .../lac/la_parallel_block_vector.templates.h | 775 ++++++++++++++ include/deal.II/lac/la_parallel_vector.h | 4 +- include/deal.II/lac/parallel_block_vector.h | 963 +----------------- include/deal.II/lac/parallel_vector.h | 4 +- source/lac/CMakeLists.txt | 2 + source/lac/la_parallel_block_vector.cc | 46 + source/lac/la_parallel_block_vector.inst.in | 41 + source/lac/la_parallel_vector.inst.in | 1 + .../trilinos_epetra_communication_pattern.cc | 2 +- tests/mpi/parallel_block_vector_01.cc | 80 +- .../parallel_block_vector_01.mpirun=1.output | 12 - .../parallel_block_vector_01.mpirun=10.output | 12 - .../parallel_block_vector_01.mpirun=4.output | 12 - 14 files changed, 1478 insertions(+), 1075 deletions(-) create mode 100644 include/deal.II/lac/la_parallel_block_vector.h create mode 100644 include/deal.II/lac/la_parallel_block_vector.templates.h create mode 100644 source/lac/la_parallel_block_vector.cc create mode 100644 source/lac/la_parallel_block_vector.inst.in 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 index 0000000000..11691b6908 --- /dev/null +++ b/include/deal.II/lac/la_parallel_block_vector.h @@ -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 +#include +#include +#include +#include + + +#include +#include + +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 @ + * and @; 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 + class BlockVector : public BlockVectorBase >, + public VectorSpaceVector + { + public: + /** + * Typedef the base class for simpler access to its own typedefs. + */ + typedef BlockVectorBase > 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 num_blocks blocks, but these + * blocks have size zero. The third variant finally initializes all + * blocks to the same size block_size. + * + * 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 &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 OtherNumber to + * Number. 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 + explicit + BlockVector (const BlockVector &v); +#endif + + /** + * Constructor. Set the number of blocks to block_sizes.size() + * and initialize each block with block_sizes[i] zero elements. + */ + BlockVector (const std::vector &block_sizes); + + /** + * Construct a block vector with an IndexSet for the local range and + * ghost entries for each block. + */ + BlockVector (const std::vector &local_ranges, + const std::vector &ghost_indices, + const MPI_Comm communicator); + + /** + * Same as above but the ghost indices are assumed to be empty. + */ + BlockVector (const std::vector &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 + BlockVector & + operator= (const BlockVector &V); + + /** + * Copy a regular vector into a block vector. + */ + BlockVector & + operator= (const Vector &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 & + 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 & + operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec); +#endif + + /** + * Reinitialize the BlockVector to contain num_blocks blocks of + * size block_size 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 omit_zeroing_entries==false, 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 + * block_sizes.size() blocks. Each block is reinitialized to + * dimension block_sizes[i]. + * + * 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 omit_zeroing_entries==false, 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 &N, + const bool omit_zeroing_entries=false); + + /** + * Change the dimension to that of the vector V. The same + * applies as for the other reinit() function. + * + * The elements of V are not copied, i.e. this function is the + * same as calling reinit (V.size(), omit_zeroing_entries). + * + * 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 + void reinit (const BlockVector &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 + void add (const std::vector &indices, + const ::dealii::Vector &values); + + /** + * Scaling and simple vector addition, i.e. *this = + * s*(*this)+V. + */ + void sadd (const Number s, + const BlockVector &V); + + /** + * Assignment *this = a*u + b*v. + * + * This function is deprecated. + */ + void equ (const Number a, const BlockVector &u, + const Number b, const BlockVector &v) DEAL_II_DEPRECATED; + + /** + * Scaling and multiple addition. + * + * This function is deprecated. + */ + void sadd (const Number s, + const Number a, + const BlockVector &V, + const Number b, + const BlockVector &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 v. 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 u.swap(v), again in analogy to standard + * functions. + */ + void swap (BlockVector &v); + //@} + + /** + * @name 2: Implementation of VectorSpaceVector + */ + //@{ + + /** + * Multiply the entire vector by a fixed factor. + */ + virtual BlockVector &operator*= (const Number factor); + + /** + * Divide the entire vector by a fixed factor. + */ + virtual BlockVector &operator/= (const Number factor); + + /** + * Add the vector @p V to the present one. + */ + virtual BlockVector &operator+= (const VectorSpaceVector &V); + + /** + * Subtract the vector @p V from the present one. + */ + virtual BlockVector &operator-= (const VectorSpaceVector &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 &V, + VectorOperation::values operation, + std_cxx11::shared_ptr communication_pattern = + std_cxx11::shared_ptr ()); + + /** + * Return the scalar product of two vectors. + */ + virtual Number operator* (const VectorSpaceVector &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. *this += a*V. + */ + virtual void add(const Number a, const VectorSpaceVector &V); + + /** + * Multiple addition of scaled vectors, i.e. *this += a*V+b*W. + */ + virtual void add(const Number a, const VectorSpaceVector &V, + const Number b, const VectorSpaceVector &W); + + /** + * Scaling and simple addition of a multiple of a vector, i.e. *this = + * s*(*this)+a*V. + */ + virtual void sadd(const Number s, const Number a, + const VectorSpaceVector &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 &scaling_factors); + + /** + * Assignment *this = a*V. + */ + virtual void equ(const Number a, const VectorSpaceVector &V); + + /** + * Return the l1 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 l2 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 &V, + const VectorSpaceVector &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 +inline +void swap (LinearAlgebra::distributed::BlockVector &u, + LinearAlgebra::distributed::BlockVector &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 index 0000000000..a36c5d4cf4 --- /dev/null +++ b/include/deal.II/lac/la_parallel_block_vector.templates.h @@ -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 +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace LinearAlgebra +{ + namespace distributed + { + template + BlockVector::BlockVector (const size_type n_blocks, + const size_type block_size) + { + reinit (n_blocks, block_size); + } + + + + template + BlockVector::BlockVector (const std::vector &n) + { + reinit (n, false); + } + + + template + BlockVector::BlockVector (const std::vector &local_ranges, + const std::vector &ghost_indices, + const MPI_Comm communicator) + { + std::vector sizes(local_ranges.size()); + for (unsigned int i=0; iblock_indices.reinit(sizes); + this->components.resize(this->n_blocks()); + + for (unsigned int i=0; in_blocks(); ++i) + this->block(i).reinit(local_ranges[i], ghost_indices[i], communicator); + } + + + template + BlockVector::BlockVector (const std::vector &local_ranges, + const MPI_Comm communicator) + { + std::vector sizes(local_ranges.size()); + for (unsigned int i=0; iblock_indices.reinit(sizes); + this->components.resize(this->n_blocks()); + + for (unsigned int i=0; in_blocks(); ++i) + this->block(i).reinit(local_ranges[i], communicator); + } + + + + template + BlockVector::BlockVector (const BlockVector &v) + : + BlockVectorBase > () + { + this->components.resize (v.n_blocks()); + this->block_indices = v.block_indices; + + for (size_type i=0; in_blocks(); ++i) + this->components[i] = v.components[i]; + } + + +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + + template + template + BlockVector::BlockVector (const BlockVector &v) + { + reinit (v, true); + *this = v; + } + +#endif + + + + template + void BlockVector::reinit (const size_type n_bl, + const size_type bl_sz, + const bool omit_zeroing_entries) + { + std::vector n(n_bl, bl_sz); + reinit(n, omit_zeroing_entries); + } + + + template + void BlockVector::reinit (const std::vector &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; in_blocks(); ++i) + this->components[i].reinit(n[i], omit_zeroing_entries); + } + + + + template + template + void BlockVector::reinit (const BlockVector &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; in_blocks(); ++i) + this->block(i).reinit(v.block(i), omit_zeroing_entries); + } + + + + template + BlockVector::~BlockVector () + {} + + + + template + BlockVector & + BlockVector::operator = (const value_type s) + { + AssertIsFinite(s); + + BaseClass::operator = (s); + return *this; + } + + + + template + BlockVector & + BlockVector::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; in_blocks(); ++i) + this->components[i] = v.block(i); + + this->collect_sizes(); + return *this; + } + + + + template + BlockVector & + BlockVector::operator = (const Vector &v) + { + BaseClass::operator = (v); + return *this; + } + + + + template + template + BlockVector & + BlockVector::operator = (const BlockVector &v) + { + reinit (v, true); + BaseClass::operator = (v); + return *this; + } + + + +#ifdef DEAL_II_WITH_PETSC + + template + BlockVector & + BlockVector::operator = (const PETScWrappers::MPI::BlockVector &petsc_vec) + { + AssertDimension(this->n_blocks(), petsc_vec.n_blocks()); + for (unsigned int i=0; in_blocks(); ++i) + this->block(i) = petsc_vec.block(i); + + return *this; + } + +#endif + + + +#ifdef DEAL_II_WITH_TRILINOS + + template + BlockVector & + BlockVector::operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec) + { + AssertDimension(this->n_blocks(), trilinos_vec.n_blocks()); + for (unsigned int i=0; in_blocks(); ++i) + this->block(i) = trilinos_vec.block(i); + + return *this; + } + +#endif + + + + template + void + BlockVector::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; blockn_blocks(); ++block) + this->block(block).compress_start(block + 8273, operation); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).compress_finish(operation); + } + + + + template + void + BlockVector::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; blockn_blocks(); ++block) + this->block(block).update_ghost_values_start(block + 9923); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).update_ghost_values_finish(); + } + + + + template + void + BlockVector::zero_out_ghosts () + { + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).zero_out_ghosts(); + } + + + + template + bool + BlockVector::has_ghost_elements () const + { + bool has_ghost_elements = false; + for (unsigned int block=0; blockn_blocks(); ++block) + if (this->block(block).has_ghost_elements() == true) + has_ghost_elements = true; + return has_ghost_elements; + } + + + + template + BlockVector & + BlockVector::operator *= (const Number factor) + { + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block) *= factor; + return *this; + } + + + + template + BlockVector & + BlockVector::operator /= (const Number factor) + { + operator *= (static_cast(1.)/factor); + return *this; + } + + + + template + void + BlockVector::scale (const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).scale(v.block(block)); + } + + + + template + void + BlockVector::equ (const Number a, + const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).equ(a, v.block(block)); + } + + + + template + void + BlockVector::equ (const Number a, + const BlockVector &v, + const Number b, + const BlockVector &w) + { + AssertDimension(this->n_blocks(), v.n_blocks()); + AssertDimension(this->n_blocks(), w.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).equ(a, v.block(block), b, w.block(block)); + } + + + + template + BlockVector & + BlockVector::operator += (const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block) += v.block(block); + + return *this; + } + + + + template + BlockVector & + BlockVector::operator -= (const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block) -= v.block(block); + + return *this; + } + + + + template + void + BlockVector::add (const Number a) + { + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).add(a); + } + + + + template + void + BlockVector::add (const Number a, + const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).add(a, v.block(block)); + } + + + + template + void + BlockVector::add (const Number a, + const VectorSpaceVector &vv, + const Number b, + const VectorSpaceVector &ww) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + Assert(dynamic_cast *>(&ww)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &w = dynamic_cast &>(ww); + AssertDimension(this->n_blocks(), v.n_blocks()); + + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).add(a, v.block(block), b, w.block(block)); + } + + + + template + void + BlockVector::sadd (const Number x, + const Number a, + const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).sadd(x, a, v.block(block)); + } + + + + template + void + BlockVector::sadd (const Number x, + const BlockVector &v) + { + AssertDimension(this->n_blocks(), v.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).sadd(x, v.block(block)); + } + + + + template + void + BlockVector::sadd (const Number x, + const Number a, + const BlockVector &v, + const Number b, + const BlockVector &w) + { + AssertDimension(this->n_blocks(), v.n_blocks()); + AssertDimension(this->n_blocks(), w.n_blocks()); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).sadd(x, a, v.block(block), b, w.block(block)); + } + + + + template + template + void + BlockVector::add (const std::vector &indices, + const ::dealii::Vector &values) + { + for (size_type i=0; i + bool + BlockVector::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; in_blocks(); ++i) + local_result = std::max(local_result, + -static_cast(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 + Number + BlockVector::operator * (const VectorSpaceVector &vv) const + { + Assert (this->n_blocks() > 0, ExcEmptyObject()); + + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + + Number local_result = Number(); + for (unsigned int i=0; in_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 + inline + Number + BlockVector::mean_value () const + { + Assert (this->n_blocks() > 0, ExcEmptyObject()); + + Number local_result = Number(); + for (unsigned int i=0; in_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 + inline + typename BlockVector::real_type + BlockVector::l1_norm () const + { + Assert (this->n_blocks() > 0, ExcEmptyObject()); + + real_type local_result = real_type(); + for (unsigned int i=0; in_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 + inline + typename BlockVector::real_type + BlockVector::l2_norm () const + { + Assert (this->n_blocks() > 0, ExcEmptyObject()); + + real_type local_result = real_type(); + for (unsigned int i=0; in_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 + inline + typename BlockVector::real_type + BlockVector::lp_norm (const real_type p) const + { + Assert (this->n_blocks() > 0, ExcEmptyObject()); + + real_type local_result = real_type(); + for (unsigned int i=0; in_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(1.0/p)); + else + return std::pow (local_result, static_cast(1.0/p)); + } + + + + template + inline + typename BlockVector::real_type + BlockVector::linfty_norm () const + { + Assert (this->n_blocks() > 0, ExcEmptyObject()); + + real_type local_result = real_type(); + for (unsigned int i=0; in_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 + inline + Number + BlockVector::add_and_dot (const Number a, + const VectorSpaceVector &vv, + const VectorSpaceVector &ww) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &v = dynamic_cast &>(vv); + AssertDimension(this->n_blocks(), v.n_blocks()); + + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&ww)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &w = dynamic_cast &>(ww); + AssertDimension(this->n_blocks(), w.n_blocks()); + + Number local_result = Number(); + for (unsigned int i=0; in_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 + inline + void + BlockVector::swap (BlockVector &v) + { + Assert (this->n_blocks() == v.n_blocks(), + ExcDimensionMismatch(this->n_blocks(), v.n_blocks())); + + for (size_type i=0; in_blocks(); ++i) + dealii::swap (this->components[i], v.components[i]); + dealii::swap (this->block_indices, v.block_indices); + } + + + + template + typename BlockVector::size_type + BlockVector::size () const + { + return this->block_indices.total_size(); + } + + + + template + inline + void + BlockVector::import(const LinearAlgebra::ReadWriteVector &, + VectorOperation::values, + std_cxx11::shared_ptr) + { + AssertThrow(false, ExcNotImplemented()); + } + + + + template + IndexSet + BlockVector::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; bn_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 + void + BlockVector::print(std::ostream &out, + const unsigned int precision, + const bool scientific, + const bool across) const + { + for (unsigned int b=0; bn_blocks(); ++b) + this->block(b).print(out, precision, scientific, across); + } + + + + template + std::size_t + BlockVector::memory_consumption () const + { + std::size_t mem = sizeof(this->n_blocks()); + for (size_type i=0; icomponents.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 diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h index 2cc4e0999d..52ef2df1e2 100644 --- a/include/deal.II/lac/la_parallel_vector.h +++ b/include/deal.II/lac/la_parallel_vector.h @@ -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 friend class dealii::parallel::distributed::BlockVector; + template friend class BlockVector; }; /*@}*/ diff --git a/include/deal.II/lac/parallel_block_vector.h b/include/deal.II/lac/parallel_block_vector.h index 2d5928786c..c392776217 100644 --- a/include/deal.II/lac/parallel_block_vector.h +++ b/include/deal.II/lac/parallel_block_vector.h @@ -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. // @@ -16,19 +16,12 @@ #ifndef dealii__parallel_block_vector_h #define dealii__parallel_block_vector_h - #include -#include -#include -#include -#include - -#include -#include +#include +#include +#include -#include -#include 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 - class BlockVector : public BlockVectorBase > - { - public: - /** - * Typedef the base class for simpler access to its own typedefs. - */ - typedef BlockVectorBase > 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 num_blocks blocks, but these - * blocks have size zero. The third variant finally initializes all - * blocks to the same size block_size. - * - * 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 &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 OtherNumber to - * Number. 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 - explicit - BlockVector (const BlockVector &v); -#endif - - /** - * Constructor. Set the number of blocks to block_sizes.size() - * and initialize each block with block_sizes[i] zero elements. - */ - BlockVector (const std::vector &block_sizes); - - /** - * Construct a block vector with an IndexSet for the local range and - * ghost entries for each block. - */ - BlockVector (const std::vector &local_ranges, - const std::vector &ghost_indices, - const MPI_Comm communicator); - - /** - * Same as above but the ghost indices are assumed to be empty. - */ - BlockVector (const std::vector &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 - BlockVector & - operator= (const BlockVector &V); - - /** - * Copy a regular vector into a block vector. - */ - BlockVector & - operator= (const Vector &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 & - 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 & - operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec); -#endif - - /** - * Reinitialize the BlockVector to contain num_blocks blocks of - * size block_size 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 omit_zeroing_entries==false, 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 - * block_sizes.size() blocks. Each block is reinitialized to - * dimension block_sizes[i]. - * - * 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 omit_zeroing_entries==false, 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 &N, - const bool omit_zeroing_entries=false); - - /** - * Change the dimension to that of the vector V. The same - * applies as for the other reinit() function. - * - * The elements of V are not copied, i.e. this function is the - * same as calling reinit (V.size(), omit_zeroing_entries). - * - * 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 - void reinit (const BlockVector &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 - bool operator == (const BlockVector &v) const; - - /** - * Checks for inequality of the two vectors. - */ - template - bool operator != (const BlockVector &v) const; - - /** - * Perform the inner product of two vectors. - */ - template - Number operator * (const BlockVector &V) const; - - /** - * Computes the square of the l2 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 l1 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 l2 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 lp 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 &V, - const BlockVector &W); - - /** - * Multiply each element of this vector by the corresponding element of - * v. - */ - template - void scale (const BlockVector2 &v); - - /** - * Swap the contents of this vector and the other vector v. 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 u.swap(v), again in analogy to standard - * functions. - */ - void swap (BlockVector &v); - - /** - * @addtogroup Exceptions - * @{ - */ - - /** - * Exception - */ - DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize); - //@} - }; + using LinearAlgebra::distributed::BlockVector; /*@}*/ - -#ifndef DOXYGEN - /*----------------------- Inline functions ----------------------------------*/ - - - template - inline - BlockVector::BlockVector (const size_type n_blocks, - const size_type block_size) - { - reinit (n_blocks, block_size); - } - - - - template - inline - BlockVector::BlockVector (const std::vector &n) - { - reinit (n, false); - } - - - template - inline - BlockVector::BlockVector (const std::vector &local_ranges, - const std::vector &ghost_indices, - const MPI_Comm communicator) - { - std::vector sizes(local_ranges.size()); - for (unsigned int i=0; iblock_indices.reinit(sizes); - this->components.resize(this->n_blocks()); - - for (unsigned int i=0; in_blocks(); ++i) - this->block(i).reinit(local_ranges[i], ghost_indices[i], communicator); - } - - - template - inline - BlockVector::BlockVector (const std::vector &local_ranges, - const MPI_Comm communicator) - { - std::vector sizes(local_ranges.size()); - for (unsigned int i=0; iblock_indices.reinit(sizes); - this->components.resize(this->n_blocks()); - - for (unsigned int i=0; in_blocks(); ++i) - this->block(i).reinit(local_ranges[i], communicator); - } - - - - template - inline - BlockVector::BlockVector (const BlockVector &v) - : - BlockVectorBase > () - { - this->components.resize (v.n_blocks()); - this->block_indices = v.block_indices; - - for (size_type i=0; in_blocks(); ++i) - this->components[i] = v.components[i]; - } - - -#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG - - template - template - inline - BlockVector::BlockVector (const BlockVector &v) - { - reinit (v, true); - *this = v; - } - -#endif - - - - template - inline - void BlockVector::reinit (const size_type n_bl, - const size_type bl_sz, - const bool omit_zeroing_entries) - { - std::vector n(n_bl, bl_sz); - reinit(n, omit_zeroing_entries); - } - - - template - inline - void BlockVector::reinit (const std::vector &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; in_blocks(); ++i) - this->components[i].reinit(n[i], omit_zeroing_entries); - } - - - - template - template - inline - void BlockVector::reinit (const BlockVector &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; in_blocks(); ++i) - this->block(i).reinit(v.block(i), omit_zeroing_entries); - } - - - - template - inline - BlockVector::~BlockVector () - {} - - - - template - inline - BlockVector & - BlockVector::operator = (const value_type s) - { - - AssertIsFinite(s); - - BaseClass::operator = (s); - return *this; - } - - - - template - inline - BlockVector & - BlockVector::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; in_blocks(); ++i) - this->components[i] = v.block(i); - - this->collect_sizes(); - return *this; - } - - - - template - inline - BlockVector & - BlockVector::operator = (const Vector &v) - { - BaseClass::operator = (v); - return *this; - } - - - - template - template - inline - BlockVector & - BlockVector::operator = (const BlockVector &v) - { - reinit (v, true); - BaseClass::operator = (v); - return *this; - } - - - -#ifdef DEAL_II_WITH_PETSC - - template - inline - BlockVector & - BlockVector::operator = (const PETScWrappers::MPI::BlockVector &petsc_vec) - { - AssertDimension(this->n_blocks(), petsc_vec.n_blocks()); - for (unsigned int i=0; in_blocks(); ++i) - this->block(i) = petsc_vec.block(i); - - return *this; - } - -#endif - - - -#ifdef DEAL_II_WITH_TRILINOS - - template - inline - BlockVector & - BlockVector::operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec) - { - AssertDimension(this->n_blocks(), trilinos_vec.n_blocks()); - for (unsigned int i=0; in_blocks(); ++i) - this->block(i) = trilinos_vec.block(i); - - return *this; - } - -#endif - - - - template - inline - void - BlockVector::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; blockn_blocks(); ++block) - this->block(block).compress_start(block*10 + 8273, operation); - for (unsigned int block=0; blockn_blocks(); ++block) - this->block(block).compress_finish(operation); - } - - - - template - inline - void - BlockVector::update_ghost_values () const - { - for (unsigned int block=0; blockn_blocks(); ++block) - this->block(block).update_ghost_values_start(block*10 + 9923); - for (unsigned int block=0; blockn_blocks(); ++block) - this->block(block).update_ghost_values_finish(); - } - - - - template - inline - void - BlockVector::zero_out_ghosts () - { - for (unsigned int block=0; blockn_blocks(); ++block) - this->block(block).zero_out_ghosts(); - } - - - - template - inline - bool - BlockVector::has_ghost_elements () const - { - bool has_ghost_elements = false; - for (unsigned int block=0; blockn_blocks(); ++block) - if (this->block(block).has_ghost_elements() == true) - has_ghost_elements = true; - return has_ghost_elements; - } - - - - template - inline - bool - BlockVector::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; in_blocks(); ++i) - local_result = std::max(local_result, - -static_cast(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 - inline - bool - BlockVector::is_non_negative () const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - int local_result = -1; - for (unsigned int i=0; in_blocks(); ++i) - local_result = std::max(local_result, - -static_cast(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 - template - inline - bool - BlockVector::operator == (const BlockVector &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; in_blocks(); ++i) - local_result = std::max(local_result, - static_cast(!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 - template - inline - bool - BlockVector::operator != (const BlockVector &v) const - { - return !(operator == (v)); - } - - - - template - template - inline - Number - BlockVector::operator * (const BlockVector &v) const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - AssertDimension (this->n_blocks(), v.n_blocks()); - - Number local_result = Number(); - for (unsigned int i=0; in_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 - inline - typename BlockVector::real_type - BlockVector::norm_sqr () const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - - real_type local_result = real_type(); - for (unsigned int i=0; in_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 - inline - Number - BlockVector::mean_value () const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - - Number local_result = Number(); - for (unsigned int i=0; in_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 - inline - typename BlockVector::real_type - BlockVector::l1_norm () const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - - real_type local_result = real_type(); - for (unsigned int i=0; in_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 - inline - typename BlockVector::real_type - BlockVector::l2_norm () const - { - return std::sqrt(norm_sqr()); - } - - - - template - inline - typename BlockVector::real_type - BlockVector::lp_norm (const real_type p) const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - - real_type local_result = real_type(); - for (unsigned int i=0; in_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(1.0/p)); - else - return std::pow (local_result, static_cast(1.0/p)); - } - - - - template - inline - typename BlockVector::real_type - BlockVector::linfty_norm () const - { - Assert (this->n_blocks() > 0, ExcEmptyObject()); - - real_type local_result = real_type(); - for (unsigned int i=0; in_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 - inline - Number - BlockVector::add_and_dot (const Number a, - const BlockVector &V, - const BlockVector &W) - { - Number local_result = Number(); - for (unsigned int i=0; in_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 - inline - void - BlockVector::swap (BlockVector &v) - { - Assert (this->n_blocks() == v.n_blocks(), - ExcDimensionMismatch(this->n_blocks(), v.n_blocks())); - - for (size_type i=0; in_blocks(); ++i) - dealii::swap (this->components[i], v.components[i]); - dealii::swap (this->block_indices, v.block_indices); - } - - - - template - template - void BlockVector::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 -inline -void swap (parallel::distributed::BlockVector &u, - parallel::distributed::BlockVector &v) -{ - u.swap (v); + } } DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/lac/parallel_vector.h b/include/deal.II/lac/parallel_vector.h index 146e35085d..b03b0f6121 100644 --- a/include/deal.II/lac/parallel_vector.h +++ b/include/deal.II/lac/parallel_vector.h @@ -30,8 +30,6 @@ namespace parallel { namespace distributed { - template class BlockVector; - /*! @addtogroup Vectors *@{ */ @@ -140,6 +138,8 @@ namespace parallel * @author Katharina Kormann, Martin Kronbichler, 2010, 2011 */ using LinearAlgebra::distributed::Vector; + + /*@}*/ } } diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 9ead26e3fe..924b993baf 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -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 index 0000000000..bbabfc3e7a --- /dev/null +++ b/source/lac/la_parallel_block_vector.cc @@ -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 +#include + +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& BlockVector::operator= (const BlockVector &) + +#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 index 0000000000..adcb86666e --- /dev/null +++ b/source/lac/la_parallel_block_vector.inst.in @@ -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; + \} + \} +} + +for (S1, S2 : REAL_SCALARS) +{ + namespace LinearAlgebra + \{ + namespace distributed + \{ + template void BlockVector::reinit (const BlockVector&, + const bool); + template void BlockVector::add (const std::vector &, + const ::dealii::Vector&); + \} + \} +} diff --git a/source/lac/la_parallel_vector.inst.in b/source/lac/la_parallel_vector.inst.in index e06292f4e9..9b118a8acf 100644 --- a/source/lac/la_parallel_vector.inst.in +++ b/source/lac/la_parallel_vector.inst.in @@ -34,6 +34,7 @@ for (S1, S2 : REAL_SCALARS) \{ template void Vector::reinit (const Vector&, const bool); + template S1 Vector::inner_product_local (const Vector&) const; \} \} } diff --git a/source/lac/trilinos_epetra_communication_pattern.cc b/source/lac/trilinos_epetra_communication_pattern.cc index ee7e9045ae..d84ff99456 100644 --- a/source/lac/trilinos_epetra_communication_pattern.cc +++ b/source/lac/trilinos_epetra_communication_pattern.cc @@ -45,7 +45,7 @@ namespace LinearAlgebra comm = std_cxx11::make_shared(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); diff --git a/tests/mpi/parallel_block_vector_01.cc b/tests/mpi/parallel_block_vector_01.cc index a9cefeee1f..e76fda085d 100644 --- a/tests/mpi/parallel_block_vector_01.cc +++ b/tests/mpi/parallel_block_vector_01.cc @@ -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 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 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 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; } diff --git a/tests/mpi/parallel_block_vector_01.mpirun=1.output b/tests/mpi/parallel_block_vector_01.mpirun=1.output index 2b1e9bf7df..3a3e5186c7 100644 --- a/tests/mpi/parallel_block_vector_01.mpirun=1.output +++ b/tests/mpi/parallel_block_vector_01.mpirun=1.output @@ -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 diff --git a/tests/mpi/parallel_block_vector_01.mpirun=10.output b/tests/mpi/parallel_block_vector_01.mpirun=10.output index 53261a12e3..ba129a4d03 100644 --- a/tests/mpi/parallel_block_vector_01.mpirun=10.output +++ b/tests/mpi/parallel_block_vector_01.mpirun=10.output @@ -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 diff --git a/tests/mpi/parallel_block_vector_01.mpirun=4.output b/tests/mpi/parallel_block_vector_01.mpirun=4.output index 8fd4464a52..9ffa1364ef 100644 --- a/tests/mpi/parallel_block_vector_01.mpirun=4.output +++ b/tests/mpi/parallel_block_vector_01.mpirun=4.output @@ -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 -- 2.39.5