From 498f0caef5aa3c1054ddfb1e48b6535b72add9d4 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 18 May 2016 17:46:32 +0200 Subject: [PATCH] Make parallel vector derive from VectorSpaceVector --- include/deal.II/base/partitioner.h | 34 +- include/deal.II/lac/la_parallel_vector.h | 1383 ++++++++++ ...lates.h => la_parallel_vector.templates.h} | 894 ++++++- include/deal.II/lac/parallel_vector.h | 2295 +---------------- include/deal.II/lac/precondition.h | 14 +- include/deal.II/lac/read_write_vector.h | 2 +- .../deal.II/lac/vector_operations_internal.h | 7 +- include/deal.II/lac/vector_space_vector.h | 7 +- include/deal.II/matrix_free/fe_evaluation.h | 10 +- source/base/partitioner.cc | 13 + source/lac/CMakeLists.txt | 4 +- source/lac/la_parallel_vector.cc | 50 + ...tor.inst.in => la_parallel_vector.inst.in} | 8 +- source/lac/parallel_vector.cc | 24 - tests/mpi/parallel_vector_06.cc | 76 +- tests/mpi/parallel_vector_06.mpirun=1.output | 12 - tests/mpi/parallel_vector_06.mpirun=10.output | 12 - tests/mpi/parallel_vector_06.mpirun=4.output | 12 - tests/mpi/parallel_vector_08.cc | 20 +- tests/mpi/parallel_vector_11.cc | 10 +- 20 files changed, 2311 insertions(+), 2576 deletions(-) create mode 100644 include/deal.II/lac/la_parallel_vector.h rename include/deal.II/lac/{parallel_vector.templates.h => la_parallel_vector.templates.h} (50%) create mode 100644 source/lac/la_parallel_vector.cc rename source/lac/{parallel_vector.inst.in => la_parallel_vector.inst.in} (90%) delete mode 100644 source/lac/parallel_vector.cc diff --git a/include/deal.II/base/partitioner.h b/include/deal.II/base/partitioner.h index 415782de2d..0ae382e178 100644 --- a/include/deal.II/base/partitioner.h +++ b/include/deal.II/base/partitioner.h @@ -22,6 +22,7 @@ #include #include #include +#include #include @@ -58,7 +59,7 @@ namespace Utilities * * @author Katharina Kormann, Martin Kronbichler, 2010, 2011 */ - class Partitioner + class Partitioner : public ::dealii::LinearAlgebra::CommunicationPatternBase { public: /** @@ -93,6 +94,17 @@ namespace Utilities Partitioner (const IndexSet &locally_owned_indices, const MPI_Comm communicator_in); + /** + * Reinitialize the communication pattern. The first argument @p + * vector_space_vector_index_set is the index set associated to a + * VectorSpaceVector object. The second argument @p + * read_write_vector_index_set is the index set associated to a + * ReadWriteVector object. + */ + virtual void reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator); + /** * Sets the locally owned indices. Used in the constructor. */ @@ -250,7 +262,12 @@ namespace Utilities /** * Returns the MPI communicator underlying the partitioner object. */ - const MPI_Comm &get_communicator() const; + const MPI_Comm &get_communicator() const DEAL_II_DEPRECATED; + + /** + * Returns the MPI communicator underlying the partitioner object. + */ + virtual const MPI_Comm &get_mpi_communicator() const; /** * Returns whether ghost indices have been explicitly added as a @p @@ -277,7 +294,7 @@ namespace Utilities /** * The global size of the vector over all processors */ - const types::global_dof_index global_size; + types::global_dof_index global_size; /** * The range of the vector that is stored locally. @@ -341,7 +358,7 @@ namespace Utilities /** * The MPI communicator involved in the problem */ - const MPI_Comm communicator; + MPI_Comm communicator; /** * Stores whether the ghost indices have been explicitly set. @@ -533,6 +550,15 @@ namespace Utilities + inline + const MPI_Comm & + Partitioner::get_mpi_communicator() const + { + return communicator; + } + + + inline bool Partitioner::ghost_indices_initialized() const diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h new file mode 100644 index 0000000000..e9bf996890 --- /dev/null +++ b/include/deal.II/lac/la_parallel_vector.h @@ -0,0 +1,1383 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2011 - 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_vector_h +#define dealii__la_parallel_vector_h + +#include +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace parallel +{ + namespace distributed + { + template class BlockVector; + } +} + +namespace LinearAlgebra +{ + template class ReadWriteVector; +} + +namespace LinearAlgebra +{ + namespace parallel + { + /*! @addtogroup Vectors + *@{ + */ + + /** + * Implementation of a parallel vector class. The design of this class is + * similar to the standard ::dealii::Vector class in deal.II, with the + * exception that storage is distributed with MPI. + * + * The vector is designed for the following scheme of parallel + * partitioning: + *
    + *
  • The indices held by individual processes (locally owned part) in + * the MPI parallelization form a contiguous range + * [my_first_index,my_last_index). + *
  • Ghost indices residing on arbitrary positions of other processors + * are allowed. It is in general more efficient if ghost indices are + * clustered, since they are stored as a set of intervals. The + * communication pattern of the ghost indices is determined when calling + * the function reinit (locally_owned, ghost_indices, + * communicator), and retained until the partitioning is changed. + * This allows for efficient parallel communication of indices. In + * particular, it stores the communication pattern, rather than having to + * compute it again for every communication. For more information on ghost + * vectors, see also the + * @ref GlossGhostedVector "glossary entry on vectors with ghost elements". + *
  • Besides the usual global access operator() it is also possible to + * access vector entries in the local index space with the function @p + * local_element(). Locally owned indices are placed first, [0, + * local_size()), and then all ghost indices follow after them + * contiguously, [local_size(), local_size()+n_ghost_entries()). + *
+ * + * Functions related to parallel functionality: + *
    + *
  • The function compress() goes through the data + * associated with ghost indices and communicates it to the owner process, + * which can then add it to the correct position. This can be used e.g. + * after having run an assembly routine involving ghosts that fill this + * vector. Note that the @p insert mode of @p compress() does not set the + * elements included in ghost entries but simply discards them, assuming + * that the owning processor has set them to the desired value already + * (See also the + * @ref GlossCompress "glossary entry on compress"). + *
  • The update_ghost_values() function imports the data + * from the owning processor to the ghost indices in order to provide read + * access to the data associated with ghosts. + *
  • It is possible to split the above functions into two phases, where + * the first initiates the communication and the second one finishes it. + * These functions can be used to overlap communication with computations + * in other parts of the code. + *
  • Of course, reduction operations (like norms) make use of + * collective all-to-all MPI communications. + *
+ * + * This vector can take two different states with respect to ghost + * elements: + *
    + *
  • After creation and whenever zero_out_ghosts() is called (or + * operator= (0.)), the vector does only allow writing into + * ghost elements but not reading from ghost elements. + *
  • After a call to update_ghost_values(), the vector does not allow + * writing into ghost elements but only reading from them. This is to + * avoid undesired ghost data artifacts when calling compress() after + * modifying some vector entries. The current status of the ghost entries + * (read mode or write mode) can be queried by the method + * has_ghost_elements(), which returns true exactly when + * ghost elements have been updated and false otherwise, + * irrespective of the actual number of ghost entries in the vector layout + * (for that information, use n_ghost_entries() instead). + *
+ * + * This vector uses the facilities of the class dealii::Vector for + * implementing the operations on the local range of the vector. In + * particular, it also inherits thread parallelism that splits most + * vector-vector operations into smaller chunks if the program uses + * multiple threads. This may or may not be desired when working also with + * MPI. + * + *

Limitations regarding the vector size

+ * + * This vector class is based on two different number types for indexing. + * The so-called global index type encodes the overall size of the vector. + * Its type is types::global_dof_index. The largest possible value is + * 2^32-1 or approximately 4 billion in case 64 bit integers + * are disabled at configuration of deal.II (default case) or + * 2^64-1 or approximately 10^19 if 64 bit + * integers are enabled (see the glossary entry on + * @ref GlobalDoFIndex + * for further information). + * + * The second relevant index type is the local index used within one MPI + * rank. As opposed to the global index, the implementation assumes 32-bit + * unsigned integers unconditionally. In other words, to actually use a + * vector with more than four billion entries, you need to use MPI with + * more than one rank (which in general is a safe assumption since four + * billion entries consume at least 16 GB of memory for floats or 32 GB of + * memory for doubles) and enable 64-bit indices. If more than 4 billion + * local elements are present, the implementation tries to detect that, + * which triggers an exception and aborts the code. Note, however, that + * the detection of overflow is tricky and the detection mechanism might + * fail in some circumstances. Therefore, it is strongly recommended to + * not rely on this class to automatically detect the unsupported case. + * + * @author Katharina Kormann, Martin Kronbichler, 2010, 2011, 2016 + */ + template + class Vector : public ::dealii::LinearAlgebra::VectorSpaceVector, + public Subscriptor + { + public: + typedef Number value_type; + typedef value_type *pointer; + typedef const value_type *const_pointer; + typedef value_type *iterator; + typedef const value_type *const_iterator; + typedef value_type &reference; + typedef const value_type &const_reference; + typedef types::global_dof_index size_type; + typedef typename numbers::NumberTraits::real_type real_type; + + /** + * A variable that indicates whether this vector supports distributed + * data storage. If true, then this vector also needs an appropriate + * compress() function that allows communicating recent set or add + * operations to individual elements to be communicated to other + * processors. + * + * For the current class, the variable equals true, since it does + * support parallel data storage. + */ + static const bool supports_distributed_data = true; + + /** + * @name 1: Basic Object-handling + */ + //@{ + /** + * Empty constructor. + */ + Vector (); + + /** + * Copy constructor. Uses the parallel partitioning of @p in_vector. + */ + Vector (const Vector &in_vector); + + /** + * Constructs a parallel vector of the given global size without any + * actual parallel distribution. + */ + Vector (const size_type size); + + /** + * Constructs a parallel vector. The local range is specified by @p + * locally_owned_set (note that this must be a contiguous interval, + * multiple intervals are not possible). The IndexSet @p ghost_indices + * specifies ghost indices, i.e., indices which one might need to read + * data from or accumulate data from. It is allowed that the set of + * ghost indices also contains the local range, but it does not need to. + * + * This function involves global communication, so it should only be + * called once for a given layout. Use the constructor with + * Vector argument to create additional vectors with the same + * parallel layout. + * + * @see + * @ref GlossGhostedVector "vectors with ghost elements" + */ + Vector (const IndexSet &local_range, + const IndexSet &ghost_indices, + const MPI_Comm communicator); + + /** + * Same constructor as above but without any ghost indices. + */ + Vector (const IndexSet &local_range, + const MPI_Comm communicator); + + /** + * Create the vector based on the parallel partitioning described in @p + * partitioner. The input argument is a shared pointer, which store the + * partitioner data only once and share it between several vectors with + * the same layout. + */ + Vector (const std_cxx11::shared_ptr &partitioner); + + /** + * Destructor. + */ + virtual ~Vector (); + + /** + * Sets the global size of the vector to @p size without any actual + * parallel distribution. + */ + void reinit (const size_type size, + const bool omit_zeroing_entries = false); + + /** + * Uses the parallel layout of the input vector @p in_vector and + * allocates memory for this vector. Recommended initialization function + * when several vectors with the same layout should be created. + * + * If the flag @p omit_zeroing_entries is set to false, the memory will + * be initialized with zero, otherwise the memory will be untouched (and + * the user must make sure to fill it with reasonable data before using + * it). + */ + template + void reinit(const Vector &in_vector, + const bool omit_zeroing_entries = false); + + /** + * Initialize the vector. The local range is specified by @p + * locally_owned_set (note that this must be a contiguous interval, + * multiple intervals are not possible). The IndexSet @p ghost_indices + * specifies ghost indices, i.e., indices which one might need to read + * data from or accumulate data from. It is allowed that the set of + * ghost indices also contains the local range, but it does not need to. + * + * This function involves global communication, so it should only be + * called once for a given layout. Use the @p reinit function with + * Vector argument to create additional vectors with the same + * parallel layout. + * + * @see + * @ref GlossGhostedVector "vectors with ghost elements" + */ + void reinit (const IndexSet &local_range, + const IndexSet &ghost_indices, + const MPI_Comm communicator); + + /** + * Same as above, but without ghost entries. + */ + void reinit (const IndexSet &local_range, + const MPI_Comm communicator); + + /** + * Initialize the vector given to the parallel partitioning described in + * @p partitioner. The input argument is a shared pointer, which store + * the partitioner data only once and share it between several vectors + * with the same layout. + */ + void reinit (const std_cxx11::shared_ptr &partitioner); + + /** + * Swap the contents of this vector and the other vector @p 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. + * + * This function is analog to the the @p 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. + * + * This function is virtual in order to allow for derived classes to + * handle memory separately. + */ + void swap (Vector &v); + + /** + * Assigns the vector to the parallel partitioning of the input vector + * @p in_vector, and copies all the data. + * + * If one of the input vector or the calling vector (to the left of the + * assignment operator) had ghost elements set before this operation, + * the calling vector will have ghost values set. Otherwise, it will be + * in write mode. If the input vector does not have any ghost elements + * at all, the vector will also update its ghost values in analogy to + * the respective setting the Trilinos and PETSc vectors. + */ + Vector & + operator = (const Vector &in_vector); + + /** + * Assigns the vector to the parallel partitioning of the input vector + * @p in_vector, and copies all the data. + * + * If one of the input vector or the calling vector (to the left of the + * assignment operator) had ghost elements set before this operation, + * the calling vector will have ghost values set. Otherwise, it will be + * in write mode. If the input vector does not have any ghost elements + * at all, the vector will also update its ghost values in analogy to + * the respective setting the Trilinos and PETSc vectors. + */ + template + Vector & + operator = (const Vector &in_vector); + //@} + + /** + * @name 2: Parallel data exchange + */ + //@{ + /** + * 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 thus simply zeroed out (using zero_ghost_values()). + * In debug mode, a check is performed for whether 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. + * + * After calling this method, write access to ghost elements of the + * vector is forbidden and an exception is thrown. Only read access to + * ghost elements is allowed in this state. Note that all subsequent + * operations on this vector, like global vector addition, etc., will + * also update the ghost values by a call to this method after the + * operation. However, global reduction operations like norms or the + * inner product will always ignore ghost elements in order to avoid + * counting the ghost data more than once. To allow writing to ghost + * elements again, call zero_out_ghosts(). + * + * @see + * @ref GlossGhostedVector "vectors with ghost elements" + */ + void update_ghost_values () const; + + /** + * Initiates communication for the @p compress() function with non- + * blocking communication. This function does not wait for the transfer + * to finish, in order to allow for other computations during the time + * it takes until all data arrives. + * + * Before the data is actually exchanged, the function must be followed + * by a call to @p compress_finish(). + * + * In case this function is called for more than one vector before @p + * compress_finish() is invoked, it is mandatory to specify a unique + * communication channel to each such call, in order to avoid several + * messages with the same ID that will corrupt this operation. + */ + void compress_start (const unsigned int communication_channel = 0, + ::dealii::VectorOperation::values operation = VectorOperation::add); + + /** + * For all requests that have been initiated in compress_start, wait for + * the communication to finish. Once it is finished, add or set the data + * (depending on the flag operation) to the respective positions in the + * owning processor, and clear the contents in the ghost data fields. + * The meaning of this argument is the same as in compress(). + * + * This function should be called exactly once per vector after calling + * compress_start, otherwise the result is undefined. In particular, it + * is not well-defined to call compress_start on the same vector again + * before compress_finished has been called. However, there is no + * warning to prevent this situation. + * + * Must follow a call to the @p compress_start function. + */ + void compress_finish (::dealii::VectorOperation::values operation); + + /** + * Initiates communication for the @p update_ghost_values() function + * with non-blocking communication. This function does not wait for the + * transfer to finish, in order to allow for other computations during + * the time it takes until all data arrives. + * + * Before the data is actually exchanged, the function must be followed + * by a call to @p update_ghost_values_finish(). + * + * In case this function is called for more than one vector before @p + * update_ghost_values_finish() is invoked, it is mandatory to specify a + * unique communication channel to each such call, in order to avoid + * several messages with the same ID that will corrupt this operation. + */ + void update_ghost_values_start (const unsigned int communication_channel = 0) const; + + + /** + * For all requests that have been started in update_ghost_values_start, + * wait for the communication to finish. + * + * Must follow a call to the @p update_ghost_values_start function + * before reading data from ghost indices. + */ + void update_ghost_values_finish () 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 whether the vector currently is in a state where ghost values + * can be read or not. This is the same functionality as other parallel + * vectors have. If this method returns false, this only means that + * read-access to ghost elements is prohibited whereas write access is + * still possible (to those entries specified as ghosts during + * initialization), not that there are no ghost elements at all. + * + * @see + * @ref GlossGhostedVector "vectors with ghost elements" + */ + bool has_ghost_elements() const; + //@} + + /** + * @name 3: Implementation of VectorSpaceVector + */ + //@{ + + /** + * Multiply the entire vector by a fixed factor. + */ + virtual Vector &operator*= (const Number factor); + + /** + * Divide the entire vector by a fixed factor. + */ + virtual Vector &operator/= (const Number factor); + + /** + * Add the vector @p V to the present one. + */ + virtual Vector &operator+= (const VectorSpaceVector &V); + + /** + * Subtract the vector @p V from the present one. + */ + virtual Vector &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; + //@} + + /** + * @name 4: Other vector operations not included in VectorSpaceVector + */ + //@{ + + /** + * Sets all elements of the vector to the scalar @p s. If the scalar is + * zero, also ghost elements are set to zero, otherwise they remain + * unchanged. + */ + Vector &operator = (const Number s); + + /** + * 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 Vector &V); + + /** + * Scaling and multiple addition. + * + * This function is deprecated. + */ + void sadd (const Number s, + const Number a, + const Vector &V, + const Number b, + const Vector &W) DEAL_II_DEPRECATED; + + /** + * Assignment *this = a*u + b*v. + * + * This function is deprecated. + */ + void equ (const Number a, const Vector &u, + const Number b, const Vector &v) DEAL_II_DEPRECATED; + + //@} + + + /** + * @name 5: Entry access and local data representation + */ + //@{ + + /** + * Returns the local size of the vector, i.e., the number of indices + * owned locally. + */ + size_type local_size() const; + + /** + * Returns the half-open interval that specifies the locally owned range + * of the vector. Note that local_size() == local_range().second - + * local_range().first. + * + * This function is deprecated. + */ + std::pair local_range () const DEAL_II_DEPRECATED; + + /** + * Returns true if the given global index is in the local range of this + * processor. + * + * This function is deprecated. + */ + bool in_local_range (const size_type global_index) const DEAL_II_DEPRECATED; + + /** + * Returns the number of ghost elements present on the vector. + * + * This function is deprecated. + */ + size_type n_ghost_entries () const DEAL_II_DEPRECATED; + + /** + * Return an index set that describes which elements of this vector are + * not owned by the current processor but can be written into or read + * from locally (ghost elements). + * + * This function is deprecated. + */ + const IndexSet &ghost_elements() const DEAL_II_DEPRECATED; + + /** + * Returns whether the given global index is a ghost index on the + * present processor. Returns false for indices that are owned locally + * and for indices not present at all. + * + * This function is deprecated. + */ + bool is_ghost_entry (const types::global_dof_index global_index) const DEAL_II_DEPRECATED; + + /** + * Make the @p Vector class a bit like the vector<> class of + * the C++ standard library by returning iterators to the start and end + * of the locally owned elements of this vector. + * + * It holds that end() - begin() == local_size(). + */ + iterator begin (); + + /** + * Return constant iterator to the start of the locally owned elements + * of the vector. + */ + const_iterator begin () const; + + /** + * Return an iterator pointing to the element past the end of the array + * of locally owned entries. + */ + iterator end (); + + /** + * Return a constant iterator pointing to the element past the end of + * the array of the locally owned entries. + */ + const_iterator end () const; + + /** + * Read access to the data in the position corresponding to @p + * global_index. The index must be either in the local range of the + * vector or be specified as a ghost index at construction. + * + * Performance: O(1) for locally owned elements that represent + * a contiguous range and O(log(nranges)) for ghost + * elements (quite fast, but slower than local_element()). + */ + Number operator () (const size_type global_index) const; + + /** + * Read and write access to the data in the position corresponding to @p + * global_index. The index must be either in the local range of the + * vector or be specified as a ghost index at construction. + * + * Performance: O(1) for locally owned elements that represent + * a contiguous range and O(log(nranges)) for ghost + * elements (quite fast, but slower than local_element()). + */ + Number &operator () (const size_type global_index); + + /** + * Read access to the data in the position corresponding to @p + * global_index. The index must be either in the local range of the + * vector or be specified as a ghost index at construction. + * + * This function does the same thing as operator(). + */ + Number operator [] (const size_type global_index) const; + /** + * Read and write access to the data in the position corresponding to @p + * global_index. The index must be either in the local range of the + * vector or be specified as a ghost index at construction. + * + * This function does the same thing as operator(). + */ + Number &operator [] (const size_type global_index); + + /** + * Read access to the data field specified by @p local_index. Locally + * owned indices can be accessed with indices + * [0,local_size), and ghost indices with indices + * [local_size,local_size+ n_ghost_entries]. + * + * Performance: Direct array access (fast). + */ + Number local_element (const size_type local_index) const; + + /** + * Read and write access to the data field specified by @p local_index. + * Locally owned indices can be accessed with indices + * [0,local_size), and ghost indices with indices + * [local_size,local_size+n_ghosts]. + * + * Performance: Direct array access (fast). + */ + Number &local_element (const size_type local_index); + + /** + * A collective get operation: instead of getting individual elements of + * a vector, this function allows to get a whole set of elements at + * once. The indices of the elements to be read are stated in the first + * argument, the corresponding values are returned in the second. + */ + template + void extract_subvector_to (const std::vector &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. Useful in minimizing copying of + * data around. + */ + template + void extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; + /** + * Return whether the vector contains only elements with value zero. + * This is a collective operation. This function is expensive, because + * potentially all elements have to be checked. + */ + 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; + //@} + + /** + * @name 6: Mixed stuff + */ + //@{ + + /** + * Return a reference to the MPI communicator object in use with this + * vector. + */ + const MPI_Comm &get_mpi_communicator () const; + + /** + * Return the MPI partitioner that describes the parallel layout of the + * vector. This object can be used to initialize another vector with the + * respective reinit() call, for additional queries regarding the + * parallel communication, or the compatibility of partitioners. + */ + std_cxx11::shared_ptr + get_partitioner () const; + + /** + * Checks whether the given partitioner is compatible with the + * partitioner used for this vector. Two partitioners are compatible if + * they have the same local size and the same ghost indices. They do not + * necessarily need to be the same data field of the shared pointer. + * This is a local operation only, i.e., if only some processors decide + * that the partitioning is not compatible, only these processors will + * return @p false, whereas the other processors will return @p true. + */ + bool + partitioners_are_compatible (const Utilities::MPI::Partitioner &part) const; + + /** + * Checks whether the given partitioner is compatible with the + * partitioner used for this vector. Two partitioners are compatible if + * they have the same local size and the same ghost indices. They do not + * necessarily need to be the same data field. As opposed to + * partitioners_are_compatible(), this method checks for compatibility + * among all processors and the method only returns @p true if the + * partitioner is the same on all processors. + * + * This method performs global communication, so make sure to use it + * only in a context where all processors call it the same number of + * times. + */ + bool + partitioners_are_globally_compatible (const Utilities::MPI::Partitioner &part) const; + //@} + + /** + * Attempt to perform an operation between two incompatible vector types. + * + * @ingroup Exceptions + */ + DeclException0(ExcVectorTypeNotCompatible); + + /** + * Exception + */ + DeclException3 (ExcNonMatchingElements, + double, double, unsigned int, + << "Called compress(VectorOperation::insert), but" + << " the element received from a remote processor, value " + << std::setprecision(16) << arg1 + << ", does not match with the value " + << std::setprecision(16) << arg2 + << " on the owner processor " << arg3); + + /** + * Exception + */ + DeclException4 (ExcAccessToNonLocalElement, + size_type, size_type, size_type, size_type, + << "You tried to access element " << arg1 + << " of a distributed vector, but this element is not " + << "stored on the current processor. Note: The range of " + << "locally owned elements is " << arg2 << " to " + << arg3 << ", and there are " << arg4 << " ghost elements " + << "that this vector can access."); + + private: + /** + * Local part of all_zero(). + */ + bool all_zero_local () const; + + /** + * Local part of the inner product of two vectors. + */ + template + Number inner_product_local (const Vector &V) const; + + /** + * Local part of norm_sqr(). + */ + real_type norm_sqr_local () const; + + /** + * Local part of mean_value(). + */ + Number mean_value_local () const; + + /** + * Local part of l1_norm(). + */ + real_type l1_norm_local () const; + + /** + * Local part of lp_norm(). + */ + real_type lp_norm_local (const real_type p) const; + + /** + * Local part of linfty_norm(). + */ + real_type linfty_norm_local () const; + + /** + * Local part of the addition followed by an inner product of two + * vectors. + */ + Number add_and_dot_local (const Number a, + const Vector &V, + const Vector &W); + + /** + * Shared pointer to store the parallel partitioning information. This + * information can be shared between several vectors that have the same + * partitioning. + */ + std_cxx11::shared_ptr partitioner; + + /** + * The size that is currently allocated in the val array. + */ + size_type allocated_size; + + /** + * Pointer to the array of local elements of this vector. + */ + Number *val; + + /** + * For parallel loops with TBB, this member variable stores the affinity + * information of loops. + */ + mutable std_cxx11::shared_ptr<::dealii::parallel::internal::TBBPartitioner> thread_loop_partitioner; + + /** + * Temporary storage that holds the data that is sent to this processor + * in @p compress() or sent from this processor in @p + * update_ghost_values. + */ + mutable Number *import_data; + + /** + * Stores whether the vector currently allows for reading ghost elements + * or not. Note that this is to ensure consistent ghost data and does + * not indicate whether the vector actually can store ghost elements. In + * particular, when assembling a vector we do not allow reading + * elements, only writing them. + */ + mutable bool vector_is_ghosted; + +#ifdef DEAL_II_WITH_MPI + /** + * A vector that collects all requests from @p compress() operations. + * This class uses persistent MPI communicators, i.e., the communication + * channels are stored during successive calls to a given function. This + * reduces the overhead involved with setting up the MPI machinery, but + * it does not remove the need for a receive operation to be posted + * before the data can actually be sent. + */ + std::vector compress_requests; + + /** + * A vector that collects all requests from @p update_ghost_values() + * operations. This class uses persistent MPI communicators. + */ + mutable std::vector update_ghost_values_requests; +#endif + + /** + * A lock that makes sure that the @p compress and @p + * update_ghost_values functions give reasonable results also when used + * with several threads. + */ + mutable Threads::Mutex mutex; + + /** + * A helper function that clears the compress_requests and + * update_ghost_values_requests field. Used in reinit functions. + */ + void clear_mpi_requests (); + + /** + * A helper function that is used to resize the val array. + */ + void resize_val (const size_type new_allocated_size); + + /* + * Make all other vector types friends. + */ + template friend class Vector; + + /** + * Make BlockVector type friends. + */ + template friend class dealii::parallel::distributed::BlockVector; + }; + /*@}*/ + + + /*----------------------- Inline functions ----------------------------------*/ + +#ifndef DOXYGEN + + + template + inline + bool + Vector::has_ghost_elements () const + { + return vector_is_ghosted; + } + + + + template + inline + typename Vector::size_type + Vector::size () const + { + return partitioner->size(); + } + + + + template + inline + typename Vector::size_type + Vector::local_size () const + { + return partitioner->local_size(); + } + + + + template + inline + std::pair::size_type, + typename Vector::size_type> + Vector::local_range () const + { + return partitioner->local_range(); + } + + + + template + inline + bool + Vector::in_local_range + (const size_type global_index) const + { + return partitioner->in_local_range (global_index); + } + + + + template + inline + IndexSet + Vector::locally_owned_elements() const + { + IndexSet is (size()); + + is.add_range (local_range().first, local_range().second); + + return is; + } + + + + template + inline + typename Vector::size_type + Vector::n_ghost_entries () const + { + return partitioner->n_ghost_indices(); + } + + + + template + inline + const IndexSet & + Vector::ghost_elements() const + { + return partitioner->ghost_indices(); + } + + + + template + inline + bool + Vector::is_ghost_entry (const size_type global_index) const + { + return partitioner->is_ghost_entry (global_index); + } + + + + template + inline + typename Vector::iterator + Vector::begin () + { + return val; + } + + + + template + inline + typename Vector::const_iterator + Vector::begin () const + { + return val; + } + + + + template + inline + typename Vector::iterator + Vector::end () + { + return val+partitioner->local_size(); + } + + + + template + inline + typename Vector::const_iterator + Vector::end () const + { + return val+partitioner->local_size(); + } + + + + template + inline + Number + Vector::operator() (const size_type global_index) const + { + Assert (in_local_range (global_index) || + partitioner->ghost_indices().is_element(global_index), + ExcAccessToNonLocalElement(global_index, local_range().first, + local_range().second, + partitioner->ghost_indices().n_elements())); + // do not allow reading a vector which is not in ghost mode + Assert (in_local_range (global_index) || vector_is_ghosted == true, + ExcMessage("You tried to read a ghost element of this vector, " + "but it has not imported its ghost values.")); + return val[partitioner->global_to_local(global_index)]; + } + + + + template + inline + Number & + Vector::operator() (const size_type global_index) + { + Assert (in_local_range (global_index) || + partitioner->ghost_indices().is_element(global_index), + ExcAccessToNonLocalElement(global_index, local_range().first, + local_range().second, + partitioner->ghost_indices().n_elements())); + // we would like to prevent reading ghosts from a vector that does not + // have them imported, but this is not possible because we might be in a + // part of the code where the vector has enabled ghosts but is non-const + // (then, the compiler picks this method according to the C++ rule book + // even if a human would pick the const method when this subsequent use + // is just a read) + return val[partitioner->global_to_local (global_index)]; + } + + + + template + inline + Number + Vector::operator[] (const size_type global_index) const + { + return operator()(global_index); + } + + + + template + inline + Number & + Vector::operator[] (const size_type global_index) + { + return operator()(global_index); + } + + + + template + inline + Number + Vector::local_element (const size_type local_index) const + { + AssertIndexRange (local_index, + partitioner->local_size()+ + partitioner->n_ghost_indices()); + // do not allow reading a vector which is not in ghost mode + Assert (local_index < local_size() || vector_is_ghosted == true, + ExcMessage("You tried to read a ghost element of this vector, " + "but it has not imported its ghost values.")); + return val[local_index]; + } + + + + template + inline + Number & + Vector::local_element (const size_type local_index) + { + AssertIndexRange (local_index, + partitioner->local_size()+ + partitioner->n_ghost_indices()); + return val[local_index]; + } + + + + template + template + inline + void Vector::extract_subvector_to (const std::vector &indices, + std::vector &values) const + { + for (size_type i = 0; i < indices.size(); ++i) + values[i] = operator()(indices[i]); + } + + + + template + template + inline + void Vector::extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const + { + while (indices_begin != indices_end) + { + *values_begin = operator()(*indices_begin); + indices_begin++; + values_begin++; + } + } + + + + template + template + inline + void + Vector::add (const std::vector &indices, + const ::dealii::Vector &values) + { + AssertDimension (indices.size(), values.size()); + for (size_type i=0; ioperator()(indices[i]) += values(i); + } + } + + + + template + inline + const MPI_Comm & + Vector::get_mpi_communicator() const + { + return partitioner->get_mpi_communicator(); + } + + + + template + inline + std_cxx11::shared_ptr + Vector::get_partitioner () const + { + return partitioner; + } + +#endif + + } +} + + + +/** + * Global function @p swap 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 Vector + * @author Katharina Kormann, Martin Kronbichler, 2011 + */ +template +inline +void swap (LinearAlgebra::parallel::Vector &u, + LinearAlgebra::parallel::Vector &v) +{ + u.swap (v); +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/lac/parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h similarity index 50% rename from include/deal.II/lac/parallel_vector.templates.h rename to include/deal.II/lac/la_parallel_vector.templates.h index eca52c03d3..2e6b7cf7ea 100644 --- a/include/deal.II/lac/parallel_vector.templates.h +++ b/include/deal.II/lac/la_parallel_vector.templates.h @@ -13,23 +13,22 @@ // // --------------------------------------------------------------------- -#ifndef dealii__parallel_vector_templates_h -#define dealii__parallel_vector_templates_h +#ifndef dealii__la_parallel_vector_templates_h +#define dealii__la_parallel_vector_templates_h #include -#include -#include +#include +#include +#include -#include -#include DEAL_II_NAMESPACE_OPEN -namespace parallel +namespace LinearAlgebra { - namespace distributed + namespace parallel { template @@ -70,6 +69,7 @@ namespace parallel val = 0; allocated_size = 0; } + thread_loop_partitioner.reset(new ::dealii::parallel::internal::TBBPartitioner()); } @@ -80,12 +80,10 @@ namespace parallel const bool omit_zeroing_entries) { clear_mpi_requests(); + // check whether we need to reallocate resize_val (size); - // reset vector view - vector_view.reinit (size, val); - // delete previous content in import data if (import_data != 0) delete[] import_data; @@ -122,11 +120,7 @@ namespace parallel const size_type new_allocated_size = partitioner->local_size() + partitioner->n_ghost_indices(); resize_val (new_allocated_size); - vector_view.reinit (partitioner->local_size(), val); } - else - Assert (vector_view.size() == partitioner->local_size(), - ExcInternalError()); if (omit_zeroing_entries == false) this->operator= (Number()); @@ -143,6 +137,8 @@ namespace parallel } vector_is_ghosted = false; + + thread_loop_partitioner = v.thread_loop_partitioner; } @@ -187,7 +183,6 @@ namespace parallel const size_type new_allocated_size = partitioner->local_size() + partitioner->n_ghost_indices(); resize_val (new_allocated_size); - vector_view.reinit (partitioner->local_size(), val); // initialize to zero this->operator= (Number()); @@ -208,122 +203,215 @@ namespace parallel -#ifdef DEAL_II_WITH_PETSC + template + Vector::Vector () + : + partitioner (new Utilities::MPI::Partitioner()), + allocated_size (0), + val (0), + import_data (0) + { + reinit(0); + } - namespace internal + + + template + Vector::Vector (const Vector &v) + : + Subscriptor(), + allocated_size (0), + val (0), + import_data (0), + vector_is_ghosted (false) { - template - void copy_petsc_vector (const PETSC_Number *petsc_start_ptr, - const PETSC_Number *petsc_end_ptr, - Number *ptr) - { - std::copy(petsc_start_ptr, petsc_end_ptr, ptr); - } + reinit (v, true); - template - void copy_petsc_vector (const std::complex *petsc_start_ptr, - const std::complex *petsc_end_ptr, - std::complex *ptr) - { - std::copy(petsc_start_ptr, petsc_end_ptr, ptr); - } + thread_loop_partitioner = v.thread_loop_partitioner; + dealii::internal::Vector_copy copier(v.val, val); + internal::parallel_for(copier, partitioner->local_size(), + thread_loop_partitioner); - template - void copy_petsc_vector (const std::complex * /*petsc_start_ptr*/, - const std::complex * /*petsc_end_ptr*/, - Number * /*ptr*/) - { - AssertThrow(false, ExcMessage("Tried to copy complex -> real")); - } + zero_out_ghosts(); } + + template - Vector & - Vector::operator = (const PETScWrappers::MPI::Vector &petsc_vec) + Vector::Vector (const IndexSet &local_range, + const IndexSet &ghost_indices, + const MPI_Comm communicator) + : + allocated_size (0), + val (0), + import_data (0), + vector_is_ghosted (false) { - Assert(petsc_vec.locally_owned_elements() == locally_owned_elements(), - StandardExceptions::ExcInvalidState()); + reinit (local_range, ghost_indices, communicator); + } - // get a representation of the vector and copy it - PetscScalar *start_ptr; - int ierr = VecGetArray (static_cast(petsc_vec), &start_ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - const size_type vec_size = local_size(); - internal::copy_petsc_vector (start_ptr, start_ptr + vec_size, begin()); - // restore the representation of the vector - ierr = VecRestoreArray (static_cast(petsc_vec), &start_ptr); - AssertThrow (ierr == 0, ExcPETScError(ierr)); + template + Vector::Vector (const IndexSet &local_range, + const MPI_Comm communicator) + : + allocated_size (0), + val (0), + import_data (0), + vector_is_ghosted (false) + { + reinit (local_range, communicator); + } + - // spread ghost values between processes? - if (vector_is_ghosted || petsc_vec.has_ghost_elements()) - update_ghost_values(); - // return a pointer to this object per normal c++ operator overloading - // semantics - return *this; + template + Vector::Vector (const size_type size) + : + allocated_size (0), + val (0), + import_data (0), + vector_is_ghosted (false) + { + reinit (size, false); } -#endif + template + Vector:: + Vector (const std_cxx11::shared_ptr &partitioner) + : + allocated_size (0), + val (0), + import_data (0), + vector_is_ghosted (false) + { + reinit (partitioner); + } + + + + template + inline + Vector::~Vector () + { + resize_val(0); + + if (import_data != 0) + delete[] import_data; + import_data = 0; + + clear_mpi_requests(); + } + -#ifdef DEAL_II_WITH_TRILINOS template + inline Vector & - Vector::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec) + Vector::operator = (const Vector &c) { - if (trilinos_vec.has_ghost_elements() == false) +#ifdef _MSC_VER + return this->operator=(c); +#else + return this->template operator=(c); +#endif + } + + + + template + template + inline + Vector & + Vector::operator = (const Vector &c) + { + Assert (c.partitioner.get() != 0, ExcNotInitialized()); + + // we update ghost values whenever one of the input or output vector + // already held ghost values or when we import data from a vector with + // the same local range but different ghost layout + bool must_update_ghost_values = c.vector_is_ghosted; + + // check whether the two vectors use the same parallel partitioner. if + // not, check if all local ranges are the same (that way, we can + // exchange data between different parallel layouts). One variant which + // is included here and necessary for compatibility with the other + // distributed vector classes (Trilinos, PETSc) is the case when vector + // c does not have any ghosts (constructed without ghost elements given) + // but the current vector does: In that case, we need to exchange data + // also when none of the two vector had updated its ghost values before. + if (partitioner.get() == 0) + reinit (c, true); + else if (partitioner.get() != c.partitioner.get()) { - Assert(trilinos_vec.locally_owned_elements() == locally_owned_elements(), - StandardExceptions::ExcInvalidState()); + // local ranges are also the same if both partitioners are empty + // (even if they happen to define the empty range as [0,0) or [c,c) + // for some c!=0 in a different way). + int local_ranges_are_identical = + (local_range() == c.local_range() || + (local_range().second == local_range().first && + c.local_range().second == c.local_range().first)); + if ((c.partitioner->n_mpi_processes() > 1 && + Utilities::MPI::min(local_ranges_are_identical, + c.partitioner->get_communicator()) == 0) + || + !local_ranges_are_identical) + reinit (c, true); + else + must_update_ghost_values |= vector_is_ghosted; + + must_update_ghost_values |= + (c.partitioner->ghost_indices_initialized() == false && + partitioner->ghost_indices_initialized() == true); } else - // ghosted trilinos vector must contain the local range of this vector - // which is contiguous - { - Assert((trilinos_vec.locally_owned_elements() & locally_owned_elements()) - == locally_owned_elements(), - StandardExceptions::ExcInvalidState()); - } + must_update_ghost_values |= vector_is_ghosted; - // create on trilinos data - const std::size_t start_index = - trilinos_vec.vector_partitioner().NumMyElements() > 0 ? - trilinos_vec.vector_partitioner(). - LID(static_cast(this->local_range().first)) : 0; - const VectorView in_view (local_size(), trilinos_vec.begin()+start_index); - static_cast&>(vector_view) = - static_cast&>(in_view); - - // spread ghost values between processes? - if (vector_is_ghosted || trilinos_vec.has_ghost_elements()) - update_ghost_values(); + thread_loop_partitioner = c.thread_loop_partitioner; + dealii::internal::Vector_copy copier(c.val, val); + internal::parallel_for(copier, partitioner->local_size(), + thread_loop_partitioner); - // return a pointer to this object per normal c++ operator overloading - // semantics + if (must_update_ghost_values) + update_ghost_values(); + else + zero_out_ghosts(); return *this; } -#endif + + + template + void + Vector::compress (::dealii::VectorOperation::values operation) + { + compress_start (0, operation); + compress_finish(operation); + } template void - Vector::copy_from (const Vector &c, - const bool call_update_ghost_values) - { - AssertDimension (local_range().first, c.local_range().first); - AssertDimension (local_range().second, c.local_range().second); - AssertDimension (vector_view.size(), c.vector_view.size()); - vector_view = c.vector_view; - if (call_update_ghost_values == true) - update_ghost_values(); - else - vector_is_ghosted = false; + Vector::update_ghost_values () const + { + update_ghost_values_start (); + update_ghost_values_finish (); + } + + + + template + void + Vector::zero_out_ghosts () + { + std::fill_n (&val[partitioner->local_size()], + partitioner->n_ghost_indices(), + Number()); + vector_is_ghosted = false; } @@ -396,7 +484,6 @@ namespace parallel } AssertDimension(current_index_start, part.n_import_indices()); - Assert (part.local_size() == vector_view.size(), ExcInternalError()); current_index_start = part.local_size(); for (unsigned int i=0; i + void + Vector::import(const ReadWriteVector &V, + VectorOperation::values operation, + std_cxx11::shared_ptr communication_pattern) + { + // If no communication pattern is given, create one. Otherwise, use the + // given one. + std_cxx11::shared_ptr comm_pattern; + if (communication_pattern.get() == NULL) + { + comm_pattern.reset(new Utilities::MPI::Partitioner(locally_owned_elements(), + V.get_stored_elements(), + get_mpi_communicator())); + } + else + { + comm_pattern = + std_cxx11::dynamic_pointer_cast (communication_pattern); + AssertThrow(comm_pattern != NULL, + ExcMessage(std::string("The communication pattern is not of type ") + + "Utilities::MPI::Partitioner.")); + } + LinearAlgebra::parallel::Vector tmp_vector(comm_pattern); + + // fill entries from ReadWriteVector into the distributed vector, + // including ghost entries. this is not really efficient right now + // because indices are translated twice, once for + const IndexSet &v_stored = V.get_stored_elements(); + for (size_type i=0; i copier(tmp_vector.val, val); + internal::parallel_for(copier, partitioner->local_size(), + thread_loop_partitioner); + } + + + template void Vector::swap (Vector &v) @@ -672,15 +798,575 @@ namespace parallel #endif std::swap (partitioner, v.partitioner); + std::swap (thread_loop_partitioner, v.thread_loop_partitioner); std::swap (allocated_size, v.allocated_size); std::swap (val, v.val); std::swap (import_data, v.import_data); std::swap (vector_is_ghosted, v.vector_is_ghosted); + } + + + + template + Vector & + Vector::operator = (const Number s) + { + internal::Vector_set setter(s, val); + + internal::parallel_for(setter, partitioner->local_size(), + thread_loop_partitioner); + + // if we call Vector::operator=0, we want to zero out all the entries + // plus ghosts. + if (s==Number()) + zero_out_ghosts(); + + return *this; + } + + + + template + Vector & + Vector::operator += (const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_add_v vector_add(val, v.val); + internal::parallel_for(vector_add, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + + return *this; + } + + + + template + Vector & + Vector::operator -= (const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_subtract_v vector_subtract(val, v.val); + internal::parallel_for(vector_subtract, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + + return *this; + } + + + + template + void + Vector::add (const Number a) + { + AssertIsFinite(a); + + internal::Vectorization_add_factor vector_add(val, a); + internal::parallel_for(vector_add, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + void + Vector::add (const Number a, + const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + AssertIsFinite(a); + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_add_av vector_add(val, v.val, a); + internal::parallel_for(vector_add, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + void + Vector::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 Vector &v = dynamic_cast &>(vv); + Assert(dynamic_cast *>(&ww)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &w = dynamic_cast &>(ww); + + AssertIsFinite(a); + AssertIsFinite(b); + + AssertDimension (local_size(), v.local_size()); + AssertDimension (local_size(), w.local_size()); + + internal::Vectorization_add_avpbw vector_add(val, v.val, w.val, a, b); + internal::parallel_for(vector_add, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + void + Vector::sadd (const Number x, + const Vector &v) + { + AssertIsFinite(x); + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_sadd_xv vector_sadd(val, v.val, x); + internal::parallel_for(vector_sadd, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + void + Vector::sadd (const Number x, + const Number a, + const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + AssertIsFinite(x); + AssertIsFinite(a); + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_sadd_xav vector_sadd(val, v.val, a, x); + internal::parallel_for(vector_sadd, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } - // vector view cannot be swapped so reset it manually (without touching - // the vector elements) - vector_view.reinit (partitioner->local_size(), val); - v.vector_view.reinit (v.partitioner->local_size(), v.val); + + + template + void + Vector::sadd (const Number x, + const Number a, + const Vector &v, + const Number b, + const Vector &w) + { + AssertIsFinite(x); + AssertIsFinite(a); + AssertIsFinite(b); + + AssertDimension (local_size(), v.local_size()); + AssertDimension (local_size(), w.local_size()); + + internal::Vectorization_sadd_xavbw vector_sadd(val, v.val, w.val, + x, a, b); + internal::parallel_for(vector_sadd, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + Vector & + Vector::operator *= (const Number factor) + { + AssertIsFinite(factor); + internal::Vectorization_multiply_factor vector_multiply(val, factor); + + internal::parallel_for(vector_multiply, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + + return *this; + } + + + + template + Vector & + Vector::operator /= (const Number factor) + { + operator *= (1./factor); + return *this; + } + + + + template + void + Vector::scale (const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_scale vector_scale(val, v.val); + internal::parallel_for(vector_scale, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + void + Vector::equ (const Number a, + const VectorSpaceVector &vv) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + AssertIsFinite(a); + AssertDimension (local_size(), v.local_size()); + + internal::Vectorization_equ_au vector_equ(val, v.val, a); + internal::parallel_for(vector_equ, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + void + Vector::equ (const Number a, + const Vector &v, + const Number b, + const Vector &w) + { + AssertIsFinite(a); + AssertIsFinite(b); + + AssertDimension (local_size(), v.local_size()); + AssertDimension (local_size(), w.local_size()); + + internal::Vectorization_equ_aubv vector_equ(val, v.val, w.val, a, b); + internal::parallel_for(vector_equ, partitioner->local_size(), + thread_loop_partitioner); + + if (vector_is_ghosted) + update_ghost_values(); + } + + + + template + bool + Vector::all_zero_local () const + { + const size_type local_size = partitioner->local_size(); + for (size_type i=0; i + bool + Vector::all_zero () const + { + // 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 reduce functions below + int local_result = -static_cast(all_zero_local()); + if (partitioner->n_mpi_processes() > 1) + return -Utilities::MPI::max(local_result, + partitioner->get_communicator()); + else + return -local_result; + } + + + + template + template + Number + Vector::inner_product_local(const Vector &v) const + { + if (PointerComparison::equal (this, &v)) + return norm_sqr_local(); + + AssertDimension (partitioner->local_size(), v.partitioner->local_size()); + + Number sum; + internal::Dot dot(val, v.val); + internal::parallel_reduce (dot, partitioner->local_size(), sum, + thread_loop_partitioner); + AssertIsFinite(sum); + + return sum; + } + + + + template + Number + Vector::operator * (const VectorSpaceVector &vv) const + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&vv)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &v = dynamic_cast &>(vv); + + Number local_result = inner_product_local(v); + if (partitioner->n_mpi_processes() > 1) + return Utilities::MPI::sum (local_result, + partitioner->get_communicator()); + else + return local_result; + } + + + + template + typename Vector::real_type + Vector::norm_sqr_local () const + { + real_type sum; + internal::Norm2 norm2(val); + internal::parallel_reduce (norm2, partitioner->local_size(), sum, + thread_loop_partitioner); + AssertIsFinite(sum); + + return sum; + } + + + + template + Number + Vector::mean_value_local () const + { + Assert(size() != 0, ExcEmptyObject()); + + if (partitioner->local_size() == 0) + return Number(); + + Number sum; + internal::MeanValue mean(val); + internal::parallel_reduce (mean, partitioner->local_size(), sum, + thread_loop_partitioner); + + return sum / real_type(partitioner->local_size()); + } + + + + template + Number + Vector::mean_value () const + { + Number local_result = mean_value_local(); + if (partitioner->n_mpi_processes() > 1) + return Utilities::MPI::sum (local_result * + (real_type)partitioner->local_size(), + partitioner->get_communicator()) + /(real_type)partitioner->size(); + else + return local_result; + } + + + + template + typename Vector::real_type + Vector::l1_norm_local () const + { + real_type sum; + internal::Norm1 norm1(val); + internal::parallel_reduce (norm1, partitioner->local_size(), sum, + thread_loop_partitioner); + + return sum; + } + + + + template + typename Vector::real_type + Vector::l1_norm () const + { + real_type local_result = l1_norm_local(); + if (partitioner->n_mpi_processes() > 1) + return Utilities::MPI::sum(local_result, + partitioner->get_communicator()); + else + return local_result; + } + + + + template + typename Vector::real_type + Vector::l2_norm () const + { + real_type local_result = norm_sqr_local(); + if (partitioner->n_mpi_processes() > 1) + return std::sqrt(Utilities::MPI::sum(local_result, + partitioner->get_communicator())); + else + return std::sqrt(local_result); + } + + + + template + typename Vector::real_type + Vector::lp_norm_local (const real_type p) const + { + real_type sum; + internal::NormP normp(val, p); + internal::parallel_reduce (normp, partitioner->local_size(), sum, + thread_loop_partitioner); + return std::pow(sum, 1./p); + } + + + + template + typename Vector::real_type + Vector::lp_norm (const real_type p) const + { + const real_type local_result = lp_norm_local(p); + if (partitioner->n_mpi_processes() > 1) + return std::pow (Utilities::MPI::sum(std::pow(local_result,p), + partitioner->get_communicator()), + static_cast(1.0/p)); + else + return local_result; + } + + + + template + typename Vector::real_type + Vector::linfty_norm_local () const + { + real_type max = 0.; + + const size_type local_size = partitioner->local_size(); + for (size_type i=0; i::abs(val[i]), max); + + return max; + } + + + + template + inline + typename Vector::real_type + Vector::linfty_norm () const + { + const real_type local_result = linfty_norm_local(); + if (partitioner->n_mpi_processes() > 1) + return Utilities::MPI::max (local_result, + partitioner->get_communicator()); + else + return local_result; + } + + + + template + Number + Vector::add_and_dot_local(const Number a, + const Vector &v, + const Vector &w) + { + const size_type vec_size = partitioner->local_size(); + AssertDimension (vec_size, v.local_size()); + AssertDimension (vec_size, w.local_size()); + + Number sum; + internal::AddAndDot adder(this->val, v.val, w.val, a); + internal::parallel_reduce (adder, vec_size, sum, thread_loop_partitioner); + AssertIsFinite(sum); + return sum; + } + + + + template + Number + Vector::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 Vector &v = dynamic_cast &>(vv); + Assert(dynamic_cast *>(&ww)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &w = dynamic_cast &>(ww); + + Number local_result = add_and_dot_local(a, v, w); + if (partitioner->n_mpi_processes() > 1) + return Utilities::MPI::sum (local_result, + partitioner->get_communicator()); + else + return local_result; } diff --git a/include/deal.II/lac/parallel_vector.h b/include/deal.II/lac/parallel_vector.h index 7d472ae2df..8cb7113f42 100644 --- a/include/deal.II/lac/parallel_vector.h +++ b/include/deal.II/lac/parallel_vector.h @@ -17,14 +17,7 @@ #define dealii__parallel_vector_h #include -#include -#include -#include -#include -#include -#include -#include -#include +#include #include #include @@ -32,26 +25,6 @@ DEAL_II_NAMESPACE_OPEN -#ifdef DEAL_II_WITH_PETSC -namespace PETScWrappers -{ - namespace MPI - { - class Vector; - } -} -#endif - -#ifdef DEAL_II_WITH_TRILINOS -namespace TrilinosWrappers -{ - namespace MPI - { - class Vector; - } -} -#endif - namespace parallel { @@ -166,2272 +139,10 @@ namespace parallel * * @author Katharina Kormann, Martin Kronbichler, 2010, 2011 */ - template - class Vector : public Subscriptor - { - public: - /** - * Declare standard types used in all containers. These types parallel - * those in the C++ standard libraries vector<...> - * class. - */ - typedef Number value_type; - typedef value_type *pointer; - typedef const value_type *const_pointer; - typedef value_type *iterator; - typedef const value_type *const_iterator; - typedef value_type &reference; - typedef const value_type &const_reference; - typedef types::global_dof_index size_type; - typedef typename numbers::NumberTraits::real_type real_type; - - /** - * A variable that indicates whether this vector supports distributed - * data storage. If true, then this vector also needs an appropriate - * compress() function that allows communicating recent set or add - * operations to individual elements to be communicated to other - * processors. - * - * For the current class, the variable equals true, since it does - * support parallel data storage. - */ - static const bool supports_distributed_data = true; - - /** - * @name 1: Basic Object-handling - */ - //@{ - /** - * Empty constructor. - */ - Vector (); - - /** - * Copy constructor. Uses the parallel partitioning of @p in_vector. - */ - Vector (const Vector &in_vector); - - /** - * Constructs a parallel vector of the given global size without any - * actual parallel distribution. - */ - Vector (const size_type size); - - /** - * Constructs a parallel vector. The local range is specified by @p - * locally_owned_set (note that this must be a contiguous interval, - * multiple intervals are not possible). The IndexSet @p ghost_indices - * specifies ghost indices, i.e., indices which one might need to read - * data from or accumulate data from. It is allowed that the set of - * ghost indices also contains the local range, but it does not need to. - * - * This function involves global communication, so it should only be - * called once for a given layout. Use the constructor with - * Vector argument to create additional vectors with the same - * parallel layout. - * - * @see - * @ref GlossGhostedVector "vectors with ghost elements" - */ - Vector (const IndexSet &local_range, - const IndexSet &ghost_indices, - const MPI_Comm communicator); - - /** - * Same constructor as above but without any ghost indices. - */ - Vector (const IndexSet &local_range, - const MPI_Comm communicator); - - /** - * Create the vector based on the parallel partitioning described in @p - * partitioner. The input argument is a shared pointer, which store the - * partitioner data only once and share it between several vectors with - * the same layout. - */ - Vector (const std_cxx11::shared_ptr &partitioner); - - /** - * Destructor. - */ - ~Vector (); - - /** - * Sets the global size of the vector to @p size without any actual - * parallel distribution. - */ - void reinit (const size_type size, - const bool omit_zeroing_entries = false); - - /** - * Uses the parallel layout of the input vector @p in_vector and - * allocates memory for this vector. Recommended initialization function - * when several vectors with the same layout should be created. - * - * If the flag @p omit_zeroing_entries is set to false, the memory will - * be initialized with zero, otherwise the memory will be untouched (and - * the user must make sure to fill it with reasonable data before using - * it). - */ - template - void reinit(const Vector &in_vector, - const bool omit_zeroing_entries = false); - - /** - * Initialize the vector. The local range is specified by @p - * locally_owned_set (note that this must be a contiguous interval, - * multiple intervals are not possible). The IndexSet @p ghost_indices - * specifies ghost indices, i.e., indices which one might need to read - * data from or accumulate data from. It is allowed that the set of - * ghost indices also contains the local range, but it does not need to. - * - * This function involves global communication, so it should only be - * called once for a given layout. Use the @p reinit function with - * Vector argument to create additional vectors with the same - * parallel layout. - * - * @see - * @ref GlossGhostedVector "vectors with ghost elements" - */ - void reinit (const IndexSet &local_range, - const IndexSet &ghost_indices, - const MPI_Comm communicator); - - /** - * Same as above, but without ghost entries. - */ - void reinit (const IndexSet &local_range, - const MPI_Comm communicator); - - /** - * Initialize the vector given to the parallel partitioning described in - * @p partitioner. The input argument is a shared pointer, which store - * the partitioner data only once and share it between several vectors - * with the same layout. - */ - void reinit (const std_cxx11::shared_ptr &partitioner); - - /** - * Swap the contents of this vector and the other vector @p 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. - * - * This function is analog to the the @p 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. - * - * This function is virtual in order to allow for derived classes to - * handle memory separately. - */ - void swap (Vector &v); - - /** - * Assigns the vector to the parallel partitioning of the input vector - * @p in_vector, and copies all the data. - * - * If one of the input vector or the calling vector (to the left of the - * assignment operator) had ghost elements set before this operation, - * the calling vector will have ghost values set. Otherwise, it will be - * in write mode. If the input vector does not have any ghost elements - * at all, the vector will also update its ghost values in analogy to - * the respective setting the Trilinos and PETSc vectors. - */ - Vector & - operator = (const Vector &in_vector); - - /** - * Assigns the vector to the parallel partitioning of the input vector - * @p in_vector, and copies all the data. - * - * If one of the input vector or the calling vector (to the left of the - * assignment operator) had ghost elements set before this operation, - * the calling vector will have ghost values set. Otherwise, it will be - * in write mode. If the input vector does not have any ghost elements - * at all, the vector will also update its ghost values in analogy to - * the respective setting the Trilinos and PETSc vectors. - */ - template - Vector & - operator = (const Vector &in_vector); - -#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. - */ - Vector & - operator = (const PETScWrappers::MPI::Vector &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. - */ - Vector & - operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec); -#endif - - /** - * This method copies the local range from another vector with the same - * local range, but possibly different layout of ghost indices. - * - * This function is deprecated. - */ - void copy_from (const Vector &in_vector, - const bool call_update_ghost_values = false) DEAL_II_DEPRECATED; - - /** - * Sets all elements of the vector to the scalar @p s. If the scalar is - * zero, also ghost elements are set to zero, otherwise they remain - * unchanged. - */ - Vector &operator = (const Number s); - - /** - * 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 thus simply zeroed out (using zero_ghost_values()). - * In debug mode, a check is performed for whether 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. - * - * After calling this method, write access to ghost elements of the - * vector is forbidden and an exception is thrown. Only read access to - * ghost elements is allowed in this state. Note that all subsequent - * operations on this vector, like global vector addition, etc., will - * also update the ghost values by a call to this method after the - * operation. However, global reduction operations like norms or the - * inner product will always ignore ghost elements in order to avoid - * counting the ghost data more than once. To allow writing to ghost - * elements again, call zero_out_ghosts(). - * - * @see - * @ref GlossGhostedVector "vectors with ghost elements" - */ - void update_ghost_values () const; - - /** - * Initiates communication for the @p compress() function with non- - * blocking communication. This function does not wait for the transfer - * to finish, in order to allow for other computations during the time - * it takes until all data arrives. - * - * Before the data is actually exchanged, the function must be followed - * by a call to @p compress_finish(). - * - * In case this function is called for more than one vector before @p - * compress_finish() is invoked, it is mandatory to specify a unique - * communication channel to each such call, in order to avoid several - * messages with the same ID that will corrupt this operation. - */ - void compress_start (const unsigned int communication_channel = 0, - ::dealii::VectorOperation::values operation = VectorOperation::add); - - /** - * For all requests that have been initiated in compress_start, wait for - * the communication to finish. Once it is finished, add or set the data - * (depending on the flag operation) to the respective positions in the - * owning processor, and clear the contents in the ghost data fields. - * The meaning of this argument is the same as in compress(). - * - * This function should be called exactly once per vector after calling - * compress_start, otherwise the result is undefined. In particular, it - * is not well-defined to call compress_start on the same vector again - * before compress_finished has been called. However, there is no - * warning to prevent this situation. - * - * Must follow a call to the @p compress_start function. - */ - void compress_finish (::dealii::VectorOperation::values operation); - - /** - * Initiates communication for the @p update_ghost_values() function - * with non-blocking communication. This function does not wait for the - * transfer to finish, in order to allow for other computations during - * the time it takes until all data arrives. - * - * Before the data is actually exchanged, the function must be followed - * by a call to @p update_ghost_values_finish(). - * - * In case this function is called for more than one vector before @p - * update_ghost_values_finish() is invoked, it is mandatory to specify a - * unique communication channel to each such call, in order to avoid - * several messages with the same ID that will corrupt this operation. - */ - void update_ghost_values_start (const unsigned int communication_channel = 0) const; - - - /** - * For all requests that have been started in update_ghost_values_start, - * wait for the communication to finish. - * - * Must follow a call to the @p update_ghost_values_start function - * before reading data from ghost indices. - */ - void update_ghost_values_finish () 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 whether the vector currently is in a state where ghost values - * can be read or not. This is the same functionality as other parallel - * vectors have. If this method returns false, this only means that - * read-access to ghost elements is prohibited whereas write access is - * still possible (to those entries specified as ghosts during - * initialization), not that there are no ghost elements at all. - * - * @see - * @ref GlossGhostedVector "vectors with ghost elements" - */ - bool has_ghost_elements() const; - - /** - * Return whether the vector contains only elements with value zero. - * This is a collective operation. This function is expensive, because - * potentially all elements have to be checked. - */ - 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 Vector &v) const; - - /** - * Checks for inequality of the two vectors. - */ - template - bool operator != (const Vector &v) const; - - /** - * Perform the inner product of two vectors. - */ - template - Number operator * (const Vector &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 Vector &V, - const Vector &W); - - /** - * Returns the global size of the vector, equal to the sum of the number - * of locally owned indices among all the processors. - */ - size_type size () const; - - /** - * Returns the local size of the vector, i.e., the number of indices - * owned locally. - */ - size_type local_size() const; - - /** - * Returns the half-open interval that specifies the locally owned range - * of the vector. Note that local_size() == local_range().second - - * local_range().first. - */ - std::pair local_range () const; - - /** - * Returns true if the given global index is in the local range of this - * processor. - */ - bool in_local_range (const size_type global_index) const; - - /** - * Return an index set that describes which elements of this vector are - * owned by the current processor. Note that this index set does not - * include elements this vector may store locally as ghost elements but - * that are in fact owned by another processor. As a consequence, the - * index sets returned on different processors 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 - */ - IndexSet locally_owned_elements () const; - - /** - * Returns the number of ghost elements present on the vector. - * - * This function is deprecated. - */ - size_type n_ghost_entries () const DEAL_II_DEPRECATED; - - /** - * Return an index set that describes which elements of this vector are - * not owned by the current processor but can be written into or read - * from locally (ghost elements). - * - * This function is deprecated. - */ - const IndexSet &ghost_elements() const DEAL_II_DEPRECATED; - - /** - * Returns whether the given global index is a ghost index on the - * present processor. Returns false for indices that are owned locally - * and for indices not present at all. - * - * This function is deprecated. - */ - bool is_ghost_entry (const types::global_dof_index global_index) const DEAL_II_DEPRECATED; - - /** - * Make the @p Vector class a bit like the vector<> class of - * the C++ standard library by returning iterators to the start and end - * of the locally owned elements of this vector. - * - * It holds that end() - begin() == local_size(). - */ - iterator begin (); - - /** - * Return constant iterator to the start of the locally owned elements - * of the vector. - */ - const_iterator begin () const; - - /** - * Return an iterator pointing to the element past the end of the array - * of locally owned entries. - */ - iterator end (); - - /** - * Return a constant iterator pointing to the element past the end of - * the array of the locally owned entries. - */ - const_iterator end () const; - //@} - - - /** - * @name 2: Data-Access - */ - //@{ - - /** - * Read access to the data in the position corresponding to @p - * global_index. The index must be either in the local range of the - * vector or be specified as a ghost index at construction. - * - * Performance: O(1) for locally owned elements that represent - * a contiguous range and O(log(nranges)) for ghost - * elements (quite fast, but slower than local_element()). - */ - Number operator () (const size_type global_index) const; - - /** - * Read and write access to the data in the position corresponding to @p - * global_index. The index must be either in the local range of the - * vector or be specified as a ghost index at construction. - * - * Performance: O(1) for locally owned elements that represent - * a contiguous range and O(log(nranges)) for ghost - * elements (quite fast, but slower than local_element()). - */ - Number &operator () (const size_type global_index); - - /** - * Read access to the data in the position corresponding to @p - * global_index. The index must be either in the local range of the - * vector or be specified as a ghost index at construction. - * - * This function does the same thing as operator(). - */ - Number operator [] (const size_type global_index) const; - - /** - * Read and write access to the data in the position corresponding to @p - * global_index. The index must be either in the local range of the - * vector or be specified as a ghost index at construction. - * - * This function does the same thing as operator(). - */ - Number &operator [] (const size_type global_index); - - /** - * A collective get operation: instead of getting individual elements of - * a vector, this function allows to get a whole set of elements at - * once. The indices of the elements to be read are stated in the first - * argument, the corresponding values are returned in the second. - */ - template - void extract_subvector_to (const std::vector &indices, - std::vector &values) const; - - /** - * Just as the above, but with pointers. Useful in minimizing copying of - * data around. - */ - template - void extract_subvector_to (ForwardIterator indices_begin, - const ForwardIterator indices_end, - OutputIterator values_begin) const; - - /** - * Read access to the data field specified by @p local_index. Locally - * owned indices can be accessed with indices - * [0,local_size), and ghost indices with indices - * [local_size,local_size+ n_ghost_entries]. - * - * Performance: Direct array access (fast). - */ - Number local_element (const size_type local_index) const; - - /** - * Read and write access to the data field specified by @p local_index. - * Locally owned indices can be accessed with indices - * [0,local_size), and ghost indices with indices - * [local_size,local_size+n_ghosts]. - * - * Performance: Direct array access (fast). - */ - Number &local_element (const size_type local_index); - //@} - - - /** - * @name 3: Modification of vectors - */ - //@{ - - /** - * Add the given vector to the present one. - */ - Vector &operator += (const Vector &V); - - /** - * Subtract the given vector from the present one. - */ - Vector &operator -= (const Vector &V); - - /** - * A collective add operation: This function 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 std::vector &values); - - /** - * This is a second collective add operation. As a difference, this - * function takes a deal.II vector of values. - */ - template - void add (const std::vector &indices, - const ::dealii::Vector &values); - - /** - * Take an address where n_elements are stored contiguously and - * add them into the vector. Handles all cases which are not covered by - * the other two add() functions above. - */ - template - void add (const size_type n_elements, - const size_type *indices, - const OtherNumber *values); - - /** - * Addition of @p s to all components. Note that @p s is a scalar and - * not a vector. - */ - void add (const Number s); - - /** - * Simple vector addition, equal to the operator +=. - * - * @deprecated Use the operator += instead. - */ - void add (const Vector &V) DEAL_II_DEPRECATED; - - /** - * Simple addition of a multiple of a vector, i.e. *this += - * a*V. - */ - void add (const Number a, const Vector &V); - - /** - * Multiple addition of scaled vectors, i.e. *this += a*V+b*W. - */ - void add (const Number a, const Vector &V, - const Number b, const Vector &W); - - /** - * Scaling and simple vector addition, i.e. *this = - * s*(*this)+V. - */ - void sadd (const Number s, - const Vector &V); - - /** - * Scaling and simple addition, i.e. *this = s*(*this)+a*V. - */ - void sadd (const Number s, - const Number a, - const Vector &V); - - /** - * Scaling and multiple addition. - * - * This function is deprecated. - */ - void sadd (const Number s, - const Number a, - const Vector &V, - const Number b, - const Vector &W) DEAL_II_DEPRECATED; - - /** - * Scaling and multiple addition. *this = s*(*this)+a*V + b*W + - * c*X. - * - * This function is deprecated. - */ - void sadd (const Number s, - const Number a, - const Vector &V, - const Number b, - const Vector &W, - const Number c, - const Vector &X) DEAL_II_DEPRECATED; - - /** - * Scale each element of the vector by a constant value. - */ - Vector &operator *= (const Number factor); - - /** - * Scale each element of the vector by the inverse of the given value. - */ - Vector &operator /= (const Number factor); - - /** - * 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. - */ - void scale (const Vector &scaling_factors); - - /** - * 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. - */ - template - void scale (const Vector &scaling_factors); - - /** - * Assignment *this = a*u. - */ - void equ (const Number a, const Vector &u); - - /** - * Assignment *this = a*u. - */ - template - void equ (const Number a, const Vector &u); - - /** - * Assignment *this = a*u + b*v. - * - * This function is deprecated. - */ - void equ (const Number a, const Vector &u, - const Number b, const Vector &v) DEAL_II_DEPRECATED; - - /** - * Assignment *this = a*u + b*v + b*w. - * - * This function is deprecated. - */ - void equ (const Number a, const Vector &u, - const Number b, const Vector &v, - const Number c, const Vector &w) DEAL_II_DEPRECATED; - - /** - * Compute the elementwise ratio of the two given vectors, that is let - * this[i] = a[i]/b[i]. This is useful for example if you want - * to compute the cellwise ratio of true to estimated error. - * - * This vector is appropriately scaled to hold the result. - * - * If any of the b[i] is zero, the result is undefined. No - * attempt is made to catch such situations. - */ - void ratio (const Vector &a, - const Vector &b) DEAL_II_DEPRECATED; - //@} - - - /** - * @name 4: Mixed stuff - */ - //@{ - /** - * Return a reference to the MPI communicator object in use with this - * vector. - */ - const MPI_Comm &get_mpi_communicator () const; - - /** - * Return the MPI partitioner that describes the parallel layout of the - * vector. This object can be used to initialize another vector with the - * respective reinit() call, for additional queries regarding the - * parallel communication, or the compatibility of partitioners. - */ - std_cxx11::shared_ptr - get_partitioner () const; - - /** - * Checks whether the given partitioner is compatible with the - * partitioner used for this vector. Two partitioners are compatible if - * they have the same local size and the same ghost indices. They do not - * necessarily need to be the same data field of the shared pointer. - * This is a local operation only, i.e., if only some processors decide - * that the partitioning is not compatible, only these processors will - * return @p false, whereas the other processors will return @p true. - */ - bool - partitioners_are_compatible (const Utilities::MPI::Partitioner &part) const; - - /** - * Checks whether the given partitioner is compatible with the - * partitioner used for this vector. Two partitioners are compatible if - * they have the same local size and the same ghost indices. They do not - * necessarily need to be the same data field. As opposed to - * partitioners_are_compatible(), this method checks for compatibility - * among all processors and the method only returns @p true if the - * partitioner is the same on all processors. - * - * This method performs global communication, so make sure to use it - * only in a context where all processors call it the same number of - * times. - */ - bool - partitioners_are_globally_compatible (const Utilities::MPI::Partitioner &part) const; - - /** - * Prints the vector to the output stream @p out. - */ - void print (std::ostream &out, - const unsigned int precision = 3, - const bool scientific = true, - const bool across = true) const; - - /** - * Returns the memory consumption of this class in bytes. - */ - std::size_t memory_consumption () const; - //@} - - /** - * Exception - */ - DeclException3 (ExcNonMatchingElements, - double, double, unsigned int, - << "Called compress(VectorOperation::insert), but" - << " the element received from a remote processor, value " - << std::setprecision(16) << arg1 - << ", does not match with the value " - << std::setprecision(16) << arg2 - << " on the owner processor " << arg3); - - /** - * Exception - */ - DeclException4 (ExcAccessToNonLocalElement, - size_type, size_type, size_type, size_type, - << "You tried to access element " << arg1 - << " of a distributed vector, but this element is not " - << "stored on the current processor. Note: The range of " - << "locally owned elements is " << arg2 << " to " - << arg3 << ", and there are " << arg4 << " ghost elements " - << "that this vector can access."); - - private: - /** - * Local part of all_zero(). - */ - bool all_zero_local () const; - - /** - * Local part of is_non_negative(). - */ - bool is_non_negative_local () const; - - /** - * Local part of operator==. - */ - template - bool vectors_equal_local (const Vector &v) const; - - /** - * Local part of the inner product of two vectors. - */ - template - Number inner_product_local (const Vector &V) const; - - /** - * Local part of norm_sqr(). - */ - real_type norm_sqr_local () const; - - /** - * Local part of mean_value(). - */ - Number mean_value_local () const; - - /** - * Local part of l1_norm(). - */ - real_type l1_norm_local () const; - - /** - * Local part of lp_norm(). - */ - real_type lp_norm_local (const real_type p) const; - - /** - * Local part of linfty_norm(). - */ - real_type linfty_norm_local () const; - - /** - * Local part of the addition followed by an inner product of two - * vectors. - */ - Number add_and_dot_local (const Number a, - const Vector &V, - const Vector &W); - - /** - * Shared pointer to store the parallel partitioning information. This - * information can be shared between several vectors that have the same - * partitioning. - */ - std_cxx11::shared_ptr partitioner; - - /** - * The size that is currently allocated in the val array. - */ - size_type allocated_size; - - /** - * Pointer to the array of local elements of this vector. - */ - Number *val; - - /** - * Temporary storage that holds the data that is sent to this processor - * in @p compress() or sent from this processor in @p - * update_ghost_values. - */ - mutable Number *import_data; - - /** - * Stores whether the vector currently allows for reading ghost elements - * or not. Note that this is to ensure consistent ghost data and does - * not indicate whether the vector actually can store ghost elements. In - * particular, when assembling a vector we do not allow reading - * elements, only writing them. - */ - mutable bool vector_is_ghosted; - - /** - * Provide this class with all functionality of ::dealii::Vector by - * creating a VectorView object. - */ - VectorView vector_view; - -#ifdef DEAL_II_WITH_MPI - /** - * A vector that collects all requests from @p compress() operations. - * This class uses persistent MPI communicators, i.e., the communication - * channels are stored during successive calls to a given function. This - * reduces the overhead involved with setting up the MPI machinery, but - * it does not remove the need for a receive operation to be posted - * before the data can actually be sent. - */ - std::vector compress_requests; - - /** - * A vector that collects all requests from @p update_ghost_values() - * operations. This class uses persistent MPI communicators. - */ - mutable std::vector update_ghost_values_requests; -#endif - - /** - * A lock that makes sure that the @p compress and @p - * update_ghost_values functions give reasonable results also when used - * with several threads. - */ - mutable Threads::Mutex mutex; - - /** - * A helper function that clears the compress_requests and - * update_ghost_values_requests field. Used in reinit functions. - */ - void clear_mpi_requests (); - - /** - * A helper function that is used to resize the val array. - */ - void resize_val (const size_type new_allocated_size); - - /* - * Make all other vector types friends. - */ - template friend class Vector; - - /** - * Make BlockVector type friends. - */ - template friend class BlockVector; - }; - - /*@}*/ - - - /*----------------------- Inline functions ----------------------------------*/ - -#ifndef DOXYGEN - - template - inline - Vector::Vector () - : - partitioner (new Utilities::MPI::Partitioner()), - allocated_size (0), - val (0), - import_data (0), - vector_is_ghosted (false), - vector_view (0, static_cast(0)) - {} - - - - template - inline - Vector::Vector (const Vector &v) - : - Subscriptor(), - allocated_size (0), - val (0), - import_data (0), - vector_is_ghosted (false), - vector_view (0, static_cast(0)) - { - reinit (v, true); - vector_view = v.vector_view; - zero_out_ghosts(); - } - - - - template - inline - Vector::Vector (const IndexSet &local_range, - const IndexSet &ghost_indices, - const MPI_Comm communicator) - : - allocated_size (0), - val (0), - import_data (0), - vector_is_ghosted (false), - vector_view (0, static_cast(0)) - { - reinit (local_range, ghost_indices, communicator); - } - - - - template - inline - Vector::Vector (const IndexSet &local_range, - const MPI_Comm communicator) - : - allocated_size (0), - val (0), - import_data (0), - vector_is_ghosted (false), - vector_view (0, static_cast(0)) - { - reinit (local_range, communicator); - } - - - - template - inline - Vector::Vector (const size_type size) - : - allocated_size (0), - val (0), - import_data (0), - vector_is_ghosted (false), - vector_view (0, static_cast(0)) - { - reinit (size, false); - } - - - - template - inline - Vector:: - Vector (const std_cxx11::shared_ptr &partitioner) - : - allocated_size (0), - val (0), - import_data (0), - vector_is_ghosted (false), - vector_view (0, static_cast(0)) - { - reinit (partitioner); - } - - - - template - inline - Vector::~Vector () - { - resize_val(0); - - if (import_data != 0) - delete[] import_data; - import_data = 0; - - clear_mpi_requests(); - } - - - - template - inline - Vector & - Vector::operator = (const Vector &c) - { -#ifdef _MSC_VER - return this->operator=(c); -#else - return this->template operator=(c); -#endif - } - - - - template - template - inline - Vector & - Vector::operator = (const Vector &c) - { - Assert (c.partitioner.get() != 0, ExcNotInitialized()); - - // we update ghost values whenever one of the input or output vector - // already held ghost values or when we import data from a vector with - // the same local range but different ghost layout - bool must_update_ghost_values = c.vector_is_ghosted; - - // check whether the two vectors use the same parallel partitioner. if - // not, check if all local ranges are the same (that way, we can - // exchange data between different parallel layouts). One variant which - // is included here and necessary for compatibility with the other - // distributed vector classes (Trilinos, PETSc) is the case when vector - // c does not have any ghosts (constructed without ghost elements given) - // but the current vector does: In that case, we need to exchange data - // also when none of the two vector had updated its ghost values before. - if (partitioner.get() == 0) - reinit (c, true); - else if (partitioner.get() != c.partitioner.get()) - { - // local ranges are also the same if both partitioners are empty - // (even if they happen to define the empty range as [0,0) or [c,c) - // for some c!=0 in a different way). - int local_ranges_are_identical = - (local_range() == c.local_range() || - (local_range().second == local_range().first && - c.local_range().second == c.local_range().first)); - if ((c.partitioner->n_mpi_processes() > 1 && - Utilities::MPI::min(local_ranges_are_identical, - c.partitioner->get_communicator()) == 0) - || - !local_ranges_are_identical) - reinit (c, true); - else - must_update_ghost_values |= vector_is_ghosted; - - must_update_ghost_values |= - (c.partitioner->ghost_indices_initialized() == false && - partitioner->ghost_indices_initialized() == true); - } - else - must_update_ghost_values |= vector_is_ghosted; - - // Need to explicitly downcast to dealii::Vector to make templated - // operator= available. - AssertDimension(vector_view.size(), c.vector_view.size()); - static_cast &>(vector_view) = c.vector_view; - - if (must_update_ghost_values) - update_ghost_values(); - else - zero_out_ghosts(); - return *this; - } - - - - template - inline - void - Vector::compress (::dealii::VectorOperation::values operation) - { - compress_start (0, operation); - compress_finish(operation); - } - - - - template - inline - void - Vector::update_ghost_values () const - { - update_ghost_values_start (); - update_ghost_values_finish (); - } - - - - template - inline - void - Vector::zero_out_ghosts () - { - std::fill_n (&val[partitioner->local_size()], - partitioner->n_ghost_indices(), - Number()); - vector_is_ghosted = false; - } - - - - template - inline - bool - Vector::has_ghost_elements () const - { - return vector_is_ghosted; - } - - - - template - inline - bool - Vector::all_zero_local () const - { - return partitioner->local_size()>0 ? vector_view.all_zero () : true; - } - - - - template - inline - bool - Vector::all_zero () const - { - // 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 = -static_cast(all_zero_local()); - if (partitioner->n_mpi_processes() > 1) - return -Utilities::MPI::max(local_result, - partitioner->get_communicator()); - else - return -local_result; - } - - - - template - inline - bool - Vector::is_non_negative_local () const - { - return partitioner->local_size()>0 ? vector_view.is_non_negative () : true; - } - - - - template - inline - bool - Vector::is_non_negative () const - { - int local_result = -static_cast(is_non_negative_local()); - if (partitioner->n_mpi_processes() > 1) - return -Utilities::MPI::max(local_result, - partitioner->get_communicator()); - else - return -local_result; - } - - - - template - template - inline - bool - Vector::vectors_equal_local (const Vector &v) const - { - return partitioner->local_size()>0 ? - vector_view.template operator == (v.vector_view) - : true; - } - - - - template - template - inline - bool - Vector::operator == (const Vector &v) const - { - // 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 = static_cast(!vectors_equal_local(v)); - unsigned int result = - partitioner->n_mpi_processes() > 1 - ? - Utilities::MPI::max(local_result, partitioner->get_communicator()) - : - local_result; - return result==0; - } - - - - template - template - inline - bool - Vector::operator != (const Vector &v) const - { - return !(operator == (v)); - } - - - - template - template - inline - Number - Vector::inner_product_local(const Vector &V) const - { - // on some processors, the size might be zero, which is not allowed by - // the dealii::Vector class. Therefore, insert a check here - return (partitioner->local_size()>0 ? - vector_view.operator* (V.vector_view) - : Number()); - } - - - - template - template - inline - Number - Vector::operator * (const Vector &V) const - { - Number local_result = inner_product_local(V); - if (partitioner->n_mpi_processes() > 1) - return Utilities::MPI::sum (local_result, - partitioner->get_communicator()); - else - return local_result; - } - - - - template - inline - typename Vector::real_type - Vector::norm_sqr_local () const - { - return partitioner->local_size()>0 ? vector_view.norm_sqr() : real_type(); - } - - - - template - inline - typename Vector::real_type - Vector::norm_sqr () const - { - real_type local_result = norm_sqr_local(); - if (partitioner->n_mpi_processes() > 1) - return Utilities::MPI::sum(local_result, - partitioner->get_communicator()); - else - return local_result; - } - - - - template - inline - Number - Vector::mean_value_local () const - { - Assert (partitioner->size()!=0, ExcEmptyObject()); - return (partitioner->local_size() ? - vector_view.mean_value() - : Number()); - } - - - - template - inline - Number - Vector::mean_value () const - { - Number local_result = mean_value_local(); - if (partitioner->n_mpi_processes() > 1) - return Utilities::MPI::sum (local_result * - (real_type)partitioner->local_size(), - partitioner->get_communicator()) - /(real_type)partitioner->size(); - else - return local_result; - } - - - - template - inline - typename Vector::real_type - Vector::l1_norm_local () const - { - return partitioner->local_size() ? vector_view.l1_norm() : real_type(); - } - - - - template - inline - typename Vector::real_type - Vector::l1_norm () const - { - real_type local_result = l1_norm_local(); - if (partitioner->n_mpi_processes() > 1) - return Utilities::MPI::sum(local_result, - partitioner->get_communicator()); - else - return local_result; - } - - - - template - inline - typename Vector::real_type - Vector::l2_norm () const - { - return std::sqrt(norm_sqr()); - } - - - - template - inline - typename Vector::real_type - Vector::lp_norm_local (const real_type p) const - { - return partitioner->local_size() ? vector_view.lp_norm(p) : real_type(); - } - - - - template - inline - typename Vector::real_type - Vector::lp_norm (const real_type p) const - { - const real_type local_result = lp_norm_local(p); - if (partitioner->n_mpi_processes() > 1) - return std::pow (Utilities::MPI::sum(std::pow(local_result,p), - partitioner->get_communicator()), - static_cast(1.0/p)); - else - return local_result; - } - - - - template - inline - typename Vector::real_type - Vector::linfty_norm_local () const - { - return partitioner->local_size() ? vector_view.linfty_norm() : real_type(); - } - - - - template - inline - typename Vector::real_type - Vector::linfty_norm () const - { - const real_type local_result = linfty_norm_local(); - if (partitioner->n_mpi_processes() > 1) - return Utilities::MPI::max (local_result, - partitioner->get_communicator()); - else - return local_result; - } - - - - template - inline - Number - Vector::add_and_dot_local(const Number a, - const Vector &V, - const Vector &W) - { - // on some processors, the size might be zero, which is not allowed by - // the dealii::Vector class. Therefore, insert a check here - return (partitioner->local_size()>0 ? - vector_view.add_and_dot(a, V.vector_view, W.vector_view) - : Number()); - } - - - - template - inline - Number - Vector::add_and_dot (const Number a, - const Vector &V, - const Vector &W) - { - Number local_result = add_and_dot_local(a, V, W); - if (partitioner->n_mpi_processes() > 1) - return Utilities::MPI::sum (local_result, - partitioner->get_communicator()); - else - return local_result; - } - - - - template - inline - typename Vector::size_type - Vector::size () const - { - return partitioner->size(); - } - - - - template - inline - typename Vector::size_type - Vector::local_size () const - { - return partitioner->local_size(); - } - - - - template - inline - std::pair::size_type, - typename Vector::size_type> - Vector::local_range () const - { - return partitioner->local_range(); - } - - - - template - inline - bool - Vector::in_local_range - (const size_type global_index) const - { - return partitioner->in_local_range (global_index); - } - - - - template - inline - IndexSet - Vector::locally_owned_elements() const - { - IndexSet is (size()); - - is.add_range (local_range().first, local_range().second); - - return is; - } - - - - template - inline - typename Vector::size_type - Vector::n_ghost_entries () const - { - return partitioner->n_ghost_indices(); - } - - - - template - inline - const IndexSet & - Vector::ghost_elements() const - { - return partitioner->ghost_indices(); - } - - - - template - inline - bool - Vector::is_ghost_entry (const size_type global_index) const - { - return partitioner->is_ghost_entry (global_index); - } - - - - template - inline - typename Vector::iterator - Vector::begin () - { - return vector_view.begin(); - } - - - - template - inline - typename Vector::const_iterator - Vector::begin () const - { - return vector_view.begin(); - } - - - - template - inline - typename Vector::iterator - Vector::end () - { - return vector_view.end(); - } - - - - template - inline - typename Vector::const_iterator - Vector::end () const - { - return vector_view.end(); - } - - - - template - inline - Number - Vector::operator() (const size_type global_index) const - { - Assert (in_local_range (global_index) || - partitioner->ghost_indices().is_element(global_index), - ExcAccessToNonLocalElement(global_index, local_range().first, - local_range().second, - partitioner->ghost_indices().n_elements())); - // do not allow reading a vector which is not in ghost mode - Assert (in_local_range (global_index) || vector_is_ghosted == true, - ExcMessage("You tried to read a ghost element of this vector, " - "but it has not imported its ghost values.")); - return val[partitioner->global_to_local(global_index)]; - } - - - - template - inline - Number & - Vector::operator() (const size_type global_index) - { - Assert (in_local_range (global_index) || - partitioner->ghost_indices().is_element(global_index), - ExcAccessToNonLocalElement(global_index, local_range().first, - local_range().second, - partitioner->ghost_indices().n_elements())); - // we would like to prevent reading ghosts from a vector that does not - // have them imported, but this is not possible because we might be in a - // part of the code where the vector has enabled ghosts but is non-const - // (then, the compiler picks this method according to the C++ rule book - // even if a human would pick the const method when this subsequent use - // is just a read) - return val[partitioner->global_to_local (global_index)]; - } - - - - template - inline - Number - Vector::operator[] (const size_type global_index) const - { - return operator()(global_index); - } - - - - template - inline - Number & - Vector::operator[] (const size_type global_index) - { - return operator()(global_index); - } - - - - template - template - inline - void Vector::extract_subvector_to (const std::vector &indices, - std::vector &values) const - { - for (size_type i = 0; i < indices.size(); ++i) - values[i] = operator()(indices[i]); - } - - - - template - template - inline - void Vector::extract_subvector_to (ForwardIterator indices_begin, - const ForwardIterator indices_end, - OutputIterator values_begin) const - { - while (indices_begin != indices_end) - { - *values_begin = operator()(*indices_begin); - indices_begin++; - values_begin++; - } - } - - - - template - inline - Number - Vector::local_element (const size_type local_index) const - { - AssertIndexRange (local_index, - partitioner->local_size()+ - partitioner->n_ghost_indices()); - // do not allow reading a vector which is not in ghost mode - Assert (local_index < local_size() || vector_is_ghosted == true, - ExcMessage("You tried to read a ghost element of this vector, " - "but it has not imported its ghost values.")); - return val[local_index]; - } - - - - template - inline - Number & - Vector::local_element (const size_type local_index) - { - AssertIndexRange (local_index, - partitioner->local_size()+ - partitioner->n_ghost_indices()); - return val[local_index]; - } - - - - template - inline - Vector & - Vector::operator = (const Number s) - { - // if we call Vector::operator=0, we want to zero out all the entries - // plus ghosts. - if (partitioner->local_size() > 0) - vector_view.dealii::template Vector::operator= (s); - if (s==Number()) - zero_out_ghosts(); - - return *this; - } - - - - template - inline - Vector & - Vector::operator += (const Vector &v) - { - AssertDimension (local_size(), v.local_size()); - - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()>0) - vector_view += v.vector_view; - - if (vector_is_ghosted) - update_ghost_values(); - - return *this; - } - - - - template - inline - Vector & - Vector::operator -= (const Vector &v) - { - AssertDimension (local_size(), v.local_size()); - - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()>0) - vector_view -= v.vector_view; - - if (vector_is_ghosted) - update_ghost_values(); - - return *this; - } - - - - template - template - inline - void - Vector::add (const std::vector &indices, - const std::vector &values) - { - AssertDimension (indices.size(), values.size()); - add (indices.size(), &indices[0], &values[0]); - } - - - - template - template - inline - void - Vector::add (const std::vector &indices, - const ::dealii::Vector &values) - { - AssertDimension (indices.size(), values.size()); - add (indices.size(), &indices[0], values.begin()); - } - - - - template - template - inline - void - Vector::add (const size_type n_indices, - const size_type *indices, - const OtherNumber *values) - { - for (size_type i=0; ioperator()(indices[i]) += values[i]; - } - } - - - - template - inline - void - Vector::add (const Number a) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.add (a); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::add (const Vector &v) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.add (v.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::add (const Number a, - const Vector &v) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.add (a, v.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::add (const Number a, - const Vector &v, - const Number b, - const Vector &w) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.add (a, v.vector_view, b, w.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::sadd (const Number x, - const Vector &v) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.sadd (x, v.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::sadd (const Number x, - const Number a, - const Vector &v) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.sadd (x, a, v.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::sadd (const Number x, - const Number a, - const Vector &v, - const Number b, - const Vector &w) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.sadd (x, a, v.vector_view, b, w.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::sadd (const Number s, - const Number a, - const Vector &v, - const Number b, - const Vector &w, - const Number c, - const Vector &x) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.sadd (s, a, v.vector_view, b, w.vector_view, - c, x.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - Vector & - Vector::operator *= (const Number factor) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view *= factor; - - if (vector_is_ghosted) - update_ghost_values(); - - return *this; - } - - - - template - inline - Vector & - Vector::operator /= (const Number factor) - { - operator *= (1./factor); - return *this; - } - - - - template - inline - void - Vector::scale (const Vector &scaling_factors) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.scale (scaling_factors.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - template - inline - void - Vector::scale (const Vector &scaling_factors) - { - if (local_size()) - vector_view.template scale (scaling_factors.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::equ (const Number a, - const Vector &v) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.equ (a, v.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - template - inline - void - Vector::equ (const Number a, - const Vector &v) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.equ (a, v.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::equ (const Number a, - const Vector &v, - const Number b, - const Vector &w) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.equ (a, v.vector_view, b, w.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::equ (const Number a, - const Vector &v, - const Number b, - const Vector &w, - const Number c, - const Vector &x) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.equ (a, v.vector_view, b, w.vector_view, - c, x.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - void - Vector::ratio (const Vector &a, - const Vector &b) - { - // dealii::Vector does not allow empty fields but this might happen on - // some processors for parallel implementation - if (local_size()) - vector_view.ratio (a.vector_view, b.vector_view); - - if (vector_is_ghosted) - update_ghost_values(); - } - - - - template - inline - const MPI_Comm & - Vector::get_mpi_communicator() const - { - return partitioner->get_communicator(); - } - - - - template - inline - std_cxx11::shared_ptr - Vector::get_partitioner () const - { - return partitioner; - } - -#endif // ifndef DOXYGEN - - } // end of namespace distributed - -} // end of namespace parallel - - - -/** - * Global function @p swap 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 Vector - * @author Katharina Kormann, Martin Kronbichler, 2011 - */ -template -inline -void swap (parallel::distributed::Vector &u, - parallel::distributed::Vector &v) -{ - u.swap (v); + using LinearAlgebra::parallel::Vector; + } } - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 5cf413b32d..b633050b42 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -33,9 +33,9 @@ DEAL_II_NAMESPACE_OPEN template class Vector; template class SparseMatrix; -namespace parallel +namespace LinearAlgebra { - namespace distributed + namespace parallel { template class Vector; } @@ -1715,14 +1715,14 @@ namespace internal template inline void - vector_updates (const parallel::distributed::Vector &src, - const parallel::distributed::Vector &matrix_diagonal_inverse, + vector_updates (const LinearAlgebra::parallel::Vector &src, + const LinearAlgebra::parallel::Vector &matrix_diagonal_inverse, const bool start_zero, const double factor1, const double factor2, - parallel::distributed::Vector &update1, - parallel::distributed::Vector &update2, - parallel::distributed::Vector &dst) + LinearAlgebra::parallel::Vector &update1, + LinearAlgebra::parallel::Vector &update2, + LinearAlgebra::parallel::Vector &dst) { VectorUpdater upd(src.begin(), matrix_diagonal_inverse.begin(), start_zero, factor1, factor2, diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index 00c1d7d31d..d42be66c2f 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -529,7 +529,7 @@ namespace LinearAlgebra * For parallel loops with TBB, this member variable stores the affinity * information of loops. */ - mutable std_cxx11::shared_ptr thread_loop_partitioner; + mutable std_cxx11::shared_ptr<::dealii::parallel::internal::TBBPartitioner> thread_loop_partitioner; /** * Make all other ReadWriteVector types friends. diff --git a/include/deal.II/lac/vector_operations_internal.h b/include/deal.II/lac/vector_operations_internal.h index 7fadf02cda..593256c0eb 100644 --- a/include/deal.II/lac/vector_operations_internal.h +++ b/include/deal.II/lac/vector_operations_internal.h @@ -902,6 +902,11 @@ namespace internal // order to obtain known loop bounds for most of the work. size_type index = first; ResultType outer_results [vector_accumulation_recursion_threshold]; + + // set the zeroth element to zero to correctly handle the case where + // vec_size == 0 + outer_results[0] = ResultType(); + size_type n_chunks = vec_size / 32; const size_type remainder = vec_size % 32; Assert (remainder == 0 || n_chunks < vector_accumulation_recursion_threshold, @@ -1214,7 +1219,7 @@ namespace internal partitioner->release_one_partitioner(tbb_partitioner); result = generic_functor.do_sum(); } - else if (vec_size > 0) + else accumulate_recursive(op,0,vec_size,result); #else accumulate_recursive(op,0,vec_size,result); diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index 5decb1574e..de05d3e1e4 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -24,9 +24,12 @@ DEAL_II_NAMESPACE_OPEN -class CommunicationPatternBase; class IndexSet; -template class ReadWriteVector; +namespace LinearAlgebra +{ + class CommunicationPatternBase; + template class ReadWriteVector; +} namespace LinearAlgebra { diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index b8dcb9fc3a..0903bc2f37 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -34,9 +34,9 @@ DEAL_II_NAMESPACE_OPEN // forward declarations -namespace parallel +namespace LinearAlgebra { - namespace distributed + namespace parallel { template class Vector; } @@ -2117,7 +2117,7 @@ namespace internal template inline Number & - vector_access (parallel::distributed::Vector &vec, + vector_access (LinearAlgebra::parallel::Vector &vec, const unsigned int entry) { return vec.local_element(entry); @@ -2131,7 +2131,7 @@ namespace internal template inline Number - vector_access (const parallel::distributed::Vector &vec, + vector_access (const LinearAlgebra::parallel::Vector &vec, const unsigned int entry) { return vec.local_element(entry); @@ -2155,7 +2155,7 @@ namespace internal template inline - void check_vector_compatibility (const parallel::distributed::Vector &vec, + void check_vector_compatibility (const LinearAlgebra::parallel::Vector &vec, const internal::MatrixFreeFunctions::DoFInfo &dof_info) { Assert (vec.partitioners_are_compatible(*dof_info.vector_partitioner), diff --git a/source/base/partitioner.cc b/source/base/partitioner.cc index 3e44d8b2a3..abeae873e4 100644 --- a/source/base/partitioner.cc +++ b/source/base/partitioner.cc @@ -89,6 +89,19 @@ namespace Utilities + void + Partitioner::reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator_in) + { + have_ghost_indices = false; + communicator = communicator_in; + set_owned_indices (vector_space_vector_index_set); + set_ghost_indices (read_write_vector_index_set); + } + + + void Partitioner::set_owned_indices (const IndexSet &locally_owned_indices) { diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 9dc145c2c5..9ead26e3fe 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -28,9 +28,9 @@ SET(_src full_matrix.cc lapack_full_matrix.cc la_vector.cc + la_parallel_vector.cc matrix_lib.cc matrix_out.cc - parallel_vector.cc precondition_block.cc precondition_block_ez.cc relaxation_block.cc @@ -62,7 +62,7 @@ SET(_inst full_matrix.inst.in lapack_full_matrix.inst.in la_vector.inst.in - parallel_vector.inst.in + la_parallel_vector.inst.in precondition_block.inst.in relaxation_block.inst.in read_write_vector.inst.in diff --git a/source/lac/la_parallel_vector.cc b/source/lac/la_parallel_vector.cc new file mode 100644 index 0000000000..f62ac75049 --- /dev/null +++ b/source/lac/la_parallel_vector.cc @@ -0,0 +1,50 @@ +// --------------------------------------------------------------------- +// +// 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_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 parallel + { +#define TEMPL_COPY_CONSTRUCTOR(S1,S2) \ + template Vector& Vector::operator= (const Vector &) + +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + TEMPL_COPY_CONSTRUCTOR(double,float); + TEMPL_COPY_CONSTRUCTOR(float,double); + + + TEMPL_COPY_CONSTRUCTOR(std::complex,std::complex); + TEMPL_COPY_CONSTRUCTOR(std::complex,std::complex); + +#endif + +#undef TEMPL_COPY_CONSTRUCTOR + } +} + + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/parallel_vector.inst.in b/source/lac/la_parallel_vector.inst.in similarity index 90% rename from source/lac/parallel_vector.inst.in rename to source/lac/la_parallel_vector.inst.in index 1291f9eeeb..99b7760bd4 100644 --- a/source/lac/parallel_vector.inst.in +++ b/source/lac/la_parallel_vector.inst.in @@ -17,9 +17,9 @@ for (SCALAR : REAL_SCALARS) { - namespace parallel + namespace LinearAlgebra \{ - namespace distributed + namespace parallel \{ template class Vector; \} @@ -28,9 +28,9 @@ for (SCALAR : REAL_SCALARS) for (S1, S2 : REAL_SCALARS) { - namespace parallel + namespace LinearAlgebra \{ - namespace distributed + namespace parallel \{ template void Vector::reinit (const Vector&, const bool); diff --git a/source/lac/parallel_vector.cc b/source/lac/parallel_vector.cc deleted file mode 100644 index fd01c6db2d..0000000000 --- a/source/lac/parallel_vector.cc +++ /dev/null @@ -1,24 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1999 - 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. -// -// --------------------------------------------------------------------- - -#include -#include - -DEAL_II_NAMESPACE_OPEN - -#include "parallel_vector.inst" - - -DEAL_II_NAMESPACE_CLOSE diff --git a/tests/mpi/parallel_vector_06.cc b/tests/mpi/parallel_vector_06.cc index 2460d8226b..0ad8c1e430 100644 --- a/tests/mpi/parallel_vector_06.cc +++ b/tests/mpi/parallel_vector_06.cc @@ -105,7 +105,7 @@ void test () } // check inner product { - const double norm_sqr = v.norm_sqr(); + const double norm_sqr = v.l2_norm() * v.l2_norm(); AssertThrow (std::fabs(v * v - norm_sqr) < 1e-15, ExcInternalError()); parallel::distributed::Vector v2; @@ -120,48 +120,6 @@ void test () deallog << "Inner product: " << inner_prod << std::endl; } - // check operator == - { - parallel::distributed::Vector v2 (v); - bool equal = (v2 == v); - if (myid == 0) - deallog << " v==v2 ? " << equal << std::endl; - - bool not_equal = (v2 != v); - if (myid == 0) - deallog << " v!=v2 ? " << not_equal << std::endl; - - // change v2 on one proc only - if (myid == 0) - v2.local_element(1) = 2.2212; - - equal = (v2 == v); - if (myid == 0) - deallog << " v==v2 ? " << equal << std::endl; - not_equal = (v2 != v); - if (myid == 0) - deallog << " v!=v2 ? " << not_equal << std::endl; - - // reset - v2 = v; - equal = (v2 == v); - if (myid == 0) - deallog << " v==v2 ? " << equal << std::endl; - not_equal = (v2 != v); - if (myid == 0) - deallog << " v!=v2 ? " << not_equal << std::endl; - - // change some value on all procs - if (myid < 8) - v2.local_element(0) = -1; - equal = (v2 == v); - if (myid == 0) - deallog << " v==v2 ? " << equal << std::endl; - not_equal = (v2 != v); - if (myid == 0) - deallog << " v!=v2 ? " << not_equal << std::endl; - } - // check all_zero { bool allzero = v.all_zero(); @@ -181,38 +139,6 @@ void test () deallog << " v2==0 ? " << allzero << std::endl; } - - // check all_non_negative - { - bool allnonneg = v.is_non_negative(); - if (myid == 0) - deallog << " v>=0 ? " << allnonneg << std::endl; - parallel::distributed::Vector v2, v3; - - // vector where all processors have - // non-negative entries - v2 = v; - if (myid < 8) - v2.local_element(0) = -1; - allnonneg = v2.is_non_negative(); - if (myid == 0) - deallog << " v2>=0 ? " << allnonneg << std::endl; - - // zero vector - v3.reinit (v2); - allnonneg = v3.is_non_negative(); - if (myid == 0) - deallog << " v3>=0 ? " << allnonneg << std::endl; - - // only one processor has non-negative entry - v3 = v; - if (myid == 1 || numproc==1) - v3.local_element(0) = -1; - allnonneg = v3.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_vector_06.mpirun=1.output b/tests/mpi/parallel_vector_06.mpirun=1.output index 8ea2068ba2..b86e4d10e6 100644 --- a/tests/mpi/parallel_vector_06.mpirun=1.output +++ b/tests/mpi/parallel_vector_06.mpirun=1.output @@ -6,19 +6,7 @@ DEAL:0::linfty norm: 2.000 DEAL:0::l2.2 norm: 2.000 DEAL:0::Mean value: 1.000 DEAL:0::Inner product: 4.000 -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_vector_06.mpirun=10.output b/tests/mpi/parallel_vector_06.mpirun=10.output index d147bb71a6..9684a3ed6e 100644 --- a/tests/mpi/parallel_vector_06.mpirun=10.output +++ b/tests/mpi/parallel_vector_06.mpirun=10.output @@ -6,19 +6,7 @@ DEAL:0::linfty norm: 30.00 DEAL:0::l2.2 norm: 63.47 DEAL:0::Mean value: 15.00 DEAL:0::Inner product: 2608. -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_vector_06.mpirun=4.output b/tests/mpi/parallel_vector_06.mpirun=4.output index d16831617e..2c6ee43693 100644 --- a/tests/mpi/parallel_vector_06.mpirun=4.output +++ b/tests/mpi/parallel_vector_06.mpirun=4.output @@ -6,19 +6,7 @@ DEAL:0::linfty norm: 14.00 DEAL:0::l2.2 norm: 22.04 DEAL:0::Mean value: 7.000 DEAL:0::Inner product: 312.0 -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_vector_08.cc b/tests/mpi/parallel_vector_08.cc index 97865e5540..931cc4926c 100644 --- a/tests/mpi/parallel_vector_08.cc +++ b/tests/mpi/parallel_vector_08.cc @@ -68,7 +68,7 @@ void test () for (unsigned i=0; i (factor, vector): "; parallel::distributed::Vector z; z = v; - y.equ (1., z); + y = z; for (int i=0; i