#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/memory_consumption.h>
+#include <deal.II/lac/communication_pattern_base.h>
#include <limits>
*
* @author Katharina Kormann, Martin Kronbichler, 2010, 2011
*/
- class Partitioner
+ class Partitioner : public ::dealii::LinearAlgebra::CommunicationPatternBase
{
public:
/**
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.
*/
/**
* 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
/**
* 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.
/**
* The MPI communicator involved in the problem
*/
- const MPI_Comm communicator;
+ MPI_Comm communicator;
/**
* Stores whether the ghost indices have been explicitly set.
+ inline
+ const MPI_Comm &
+ Partitioner::get_mpi_communicator() const
+ {
+ return communicator;
+ }
+
+
+
inline
bool
Partitioner::ghost_indices_initialized() const
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/numbers.h>
+#include <deal.II/base/partitioner.h>
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/lac/vector_space_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace parallel
+{
+ namespace distributed
+ {
+ template <typename> class BlockVector;
+ }
+}
+
+namespace LinearAlgebra
+{
+ template <typename> 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:
+ * <ul>
+ * <li> The indices held by individual processes (locally owned part) in
+ * the MPI parallelization form a contiguous range
+ * <code>[my_first_index,my_last_index)</code>.
+ * <li> 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 <code>reinit (locally_owned, ghost_indices,
+ * communicator)</code>, 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".
+ * <li> 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()).
+ * </ul>
+ *
+ * Functions related to parallel functionality:
+ * <ul>
+ * <li> The function <code>compress()</code> 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").
+ * <li> The <code>update_ghost_values()</code> function imports the data
+ * from the owning processor to the ghost indices in order to provide read
+ * access to the data associated with ghosts.
+ * <li> 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.
+ * <li> Of course, reduction operations (like norms) make use of
+ * collective all-to-all MPI communications.
+ * </ul>
+ *
+ * This vector can take two different states with respect to ghost
+ * elements:
+ * <ul>
+ * <li> After creation and whenever zero_out_ghosts() is called (or
+ * <code>operator= (0.)</code>), the vector does only allow writing into
+ * ghost elements but not reading from ghost elements.
+ * <li> 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 <code>true</code> exactly when
+ * ghost elements have been updated and <code>false</code> otherwise,
+ * irrespective of the actual number of ghost entries in the vector layout
+ * (for that information, use n_ghost_entries() instead).
+ * </ul>
+ *
+ * This vector uses the facilities of the class dealii::Vector<Number> 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.
+ *
+ * <h4>Limitations regarding the vector size</h4>
+ *
+ * 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
+ * <code>2^32-1</code> or approximately 4 billion in case 64 bit integers
+ * are disabled at configuration of deal.II (default case) or
+ * <code>2^64-1</code> or approximately <code>10^19</code> 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 <typename Number>
+ class Vector : public ::dealii::LinearAlgebra::VectorSpaceVector<Number>,
+ 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<Number>::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<Number> &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<Number> 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<const Utilities::MPI::Partitioner> &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 <typename Number2>
+ void reinit(const Vector<Number2> &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<Number> 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<const Utilities::MPI::Partitioner> &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
+ * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
+ * analogy to standard functions.
+ *
+ * This function is virtual in order to allow for derived classes to
+ * handle memory separately.
+ */
+ void swap (Vector<Number> &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<Number> &
+ operator = (const Vector<Number> &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 <typename Number2>
+ Vector<Number> &
+ operator = (const Vector<Number2> &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<Number> &operator*= (const Number factor);
+
+ /**
+ * Divide the entire vector by a fixed factor.
+ */
+ virtual Vector<Number> &operator/= (const Number factor);
+
+ /**
+ * Add the vector @p V to the present one.
+ */
+ virtual Vector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+
+ /**
+ * Subtract the vector @p V from the present one.
+ */
+ virtual Vector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+
+ /**
+ * Import all the elements present in the vector's IndexSet from the input
+ * vector @p V. VectorOperation::values @p operation is used to decide if
+ * the elements in @p V should be added to the current vector or replace the
+ * current elements. The last parameter can be used if the same
+ * communication pattern is used multiple times. This can be used to improve
+ * performance.
+ */
+ virtual void import(const LinearAlgebra::ReadWriteVector<Number> &V,
+ VectorOperation::values operation,
+ std_cxx11::shared_ptr<const CommunicationPatternBase> communication_pattern =
+ std_cxx11::shared_ptr<const CommunicationPatternBase> ());
+
+ /**
+ * Return the scalar product of two vectors.
+ */
+ virtual Number operator* (const VectorSpaceVector<Number> &V) const;
+
+ /**
+ * Add @p a to all components. Note that @p a is a scalar not a vector.
+ */
+ virtual void add(const Number a);
+
+ /**
+ * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+ */
+ virtual void add(const Number a, const VectorSpaceVector<Number> &V);
+
+ /**
+ * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
+ */
+ virtual void add(const Number a, const VectorSpaceVector<Number> &V,
+ const Number b, const VectorSpaceVector<Number> &W);
+
+ /**
+ * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
+ * s*(*this)+a*V</tt>.
+ */
+ virtual void sadd(const Number s, const Number a,
+ const VectorSpaceVector<Number> &V);
+
+ /**
+ * Scale each element of this vector by the corresponding element in the
+ * argument. This function is mostly meant to simulate multiplication (and
+ * immediate re-assignment) by a diagonal scaling matrix.
+ */
+ virtual void scale(const VectorSpaceVector<Number> &scaling_factors);
+
+ /**
+ * Assignment <tt>*this = a*V</tt>.
+ */
+ virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
+
+ /**
+ * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
+ * absolute values of all entries among all processors).
+ */
+ virtual real_type l1_norm() const;
+
+ /**
+ * Return the l<sub>2</sub> norm of the vector (i.e., the square root of
+ * the sum of the square of all entries among all processors).
+ */
+ virtual real_type l2_norm() const;
+
+ /**
+ * Return the maximum norm of the vector (i.e., the maximum absolute value
+ * among all entries and among all processors).
+ */
+ virtual real_type linfty_norm() const;
+
+ /**
+ * Perform a combined operation of a vector addition and a subsequent
+ * inner product, returning the value of the inner product. In other
+ * words, the result of this function is the same as if the user called
+ * @code
+ * this->add(a, V);
+ * return_value = *this * W;
+ * @endcode
+ *
+ * The reason this function exists is that this operation involves less
+ * memory transfer than calling the two functions separately. This method
+ * only needs to load three vectors, @p this, @p V, @p W, whereas calling
+ * separate methods means to load the calling vector @p this twice. Since
+ * most vector operations are memory transfer limited, this reduces the
+ * time by 25\% (or 50\% if @p W equals @p this).
+ */
+ virtual Number add_and_dot(const Number a,
+ const VectorSpaceVector<Number> &V,
+ const VectorSpaceVector<Number> &W);
+
+ /**
+ * Return the global size of the vector, equal to the sum of the number of
+ * locally owned indices among all processors.
+ */
+ virtual size_type size() const;
+
+ /**
+ * Return an index set that describes which elements of this vector are
+ * owned by the current processor. As a consequence, the index sets
+ * returned on different procesors if this is a distributed vector will
+ * form disjoint sets that add up to the complete index set. Obviously, if
+ * a vector is created on only one processor, then the result would
+ * satisfy
+ * @code
+ * vec.locally_owned_elements() == complete_index_set(vec.size())
+ * @endcode
+ */
+ virtual dealii::IndexSet locally_owned_elements() const;
+
+ /**
+ * Print the vector to the output stream @p out.
+ */
+ virtual void print(std::ostream &out,
+ const unsigned int precision=3,
+ const bool scientific=true,
+ const bool across=true) const;
+
+ /**
+ * Return the memory consumption of this class in bytes.
+ */
+ virtual std::size_t memory_consumption() const;
+ //@}
+
+ /**
+ * @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<Number> &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 <typename OtherNumber>
+ void add (const std::vector<size_type> &indices,
+ const ::dealii::Vector<OtherNumber> &values);
+
+ /**
+ * Scaling and simple vector addition, i.e. <tt>*this =
+ * s*(*this)+V</tt>.
+ */
+ void sadd (const Number s,
+ const Vector<Number> &V);
+
+ /**
+ * Scaling and multiple addition.
+ *
+ * This function is deprecated.
+ */
+ void sadd (const Number s,
+ const Number a,
+ const Vector<Number> &V,
+ const Number b,
+ const Vector<Number> &W) DEAL_II_DEPRECATED;
+
+ /**
+ * Assignment <tt>*this = a*u + b*v</tt>.
+ *
+ * This function is deprecated.
+ */
+ void equ (const Number a, const Vector<Number> &u,
+ const Number b, const Vector<Number> &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 <code>local_size() == local_range().second -
+ * local_range().first</code>.
+ *
+ * This function is deprecated.
+ */
+ std::pair<size_type, size_type> 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 <tt>vector<></tt> class of
+ * the C++ standard library by returning iterators to the start and end
+ * of the <i>locally owned</i> 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: <tt>O(1)</tt> for locally owned elements that represent
+ * a contiguous range and <tt>O(log(n<sub>ranges</sub>))</tt> 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: <tt>O(1)</tt> for locally owned elements that represent
+ * a contiguous range and <tt>O(log(n<sub>ranges</sub>))</tt> 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
+ * <code>[0,local_size)</code>, and ghost indices with indices
+ * <code>[local_size,local_size+ n_ghost_entries]</code>.
+ *
+ * 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
+ * <code>[0,local_size)</code>, and ghost indices with indices
+ * <code>[local_size,local_size+n_ghosts]</code>.
+ *
+ * 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 <typename OtherNumber>
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const;
+
+ /**
+ * Just as the above, but with pointers. Useful in minimizing copying of
+ * data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ 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<const Utilities::MPI::Partitioner>
+ 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 <typename Number2>
+ Number inner_product_local (const Vector<Number2> &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<Number> &V,
+ const Vector<Number> &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<const Utilities::MPI::Partitioner> 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<MPI_Request> compress_requests;
+
+ /**
+ * A vector that collects all requests from @p update_ghost_values()
+ * operations. This class uses persistent MPI communicators.
+ */
+ mutable std::vector<MPI_Request> 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 <typename Number2> friend class Vector;
+
+ /**
+ * Make BlockVector type friends.
+ */
+ template <typename Number2> friend class dealii::parallel::distributed::BlockVector;
+ };
+ /*@}*/
+
+
+ /*----------------------- Inline functions ----------------------------------*/
+
+#ifndef DOXYGEN
+
+
+ template <typename Number>
+ inline
+ bool
+ Vector<Number>::has_ghost_elements () const
+ {
+ return vector_is_ghosted;
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::size_type
+ Vector<Number>::size () const
+ {
+ return partitioner->size();
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::size_type
+ Vector<Number>::local_size () const
+ {
+ return partitioner->local_size();
+ }
+
+
+
+ template <typename Number>
+ inline
+ std::pair<typename Vector<Number>::size_type,
+ typename Vector<Number>::size_type>
+ Vector<Number>::local_range () const
+ {
+ return partitioner->local_range();
+ }
+
+
+
+ template <typename Number>
+ inline
+ bool
+ Vector<Number>::in_local_range
+ (const size_type global_index) const
+ {
+ return partitioner->in_local_range (global_index);
+ }
+
+
+
+ template <typename Number>
+ inline
+ IndexSet
+ Vector<Number>::locally_owned_elements() const
+ {
+ IndexSet is (size());
+
+ is.add_range (local_range().first, local_range().second);
+
+ return is;
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::size_type
+ Vector<Number>::n_ghost_entries () const
+ {
+ return partitioner->n_ghost_indices();
+ }
+
+
+
+ template <typename Number>
+ inline
+ const IndexSet &
+ Vector<Number>::ghost_elements() const
+ {
+ return partitioner->ghost_indices();
+ }
+
+
+
+ template <typename Number>
+ inline
+ bool
+ Vector<Number>::is_ghost_entry (const size_type global_index) const
+ {
+ return partitioner->is_ghost_entry (global_index);
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::iterator
+ Vector<Number>::begin ()
+ {
+ return val;
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::const_iterator
+ Vector<Number>::begin () const
+ {
+ return val;
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::iterator
+ Vector<Number>::end ()
+ {
+ return val+partitioner->local_size();
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::const_iterator
+ Vector<Number>::end () const
+ {
+ return val+partitioner->local_size();
+ }
+
+
+
+ template <typename Number>
+ inline
+ Number
+ Vector<Number>::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 <typename Number>
+ inline
+ Number &
+ Vector<Number>::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 <typename Number>
+ inline
+ Number
+ Vector<Number>::operator[] (const size_type global_index) const
+ {
+ return operator()(global_index);
+ }
+
+
+
+ template <typename Number>
+ inline
+ Number &
+ Vector<Number>::operator[] (const size_type global_index)
+ {
+ return operator()(global_index);
+ }
+
+
+
+ template <typename Number>
+ inline
+ Number
+ Vector<Number>::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 <typename Number>
+ inline
+ Number &
+ Vector<Number>::local_element (const size_type local_index)
+ {
+ AssertIndexRange (local_index,
+ partitioner->local_size()+
+ partitioner->n_ghost_indices());
+ return val[local_index];
+ }
+
+
+
+ template <typename Number>
+ template <typename OtherNumber>
+ inline
+ void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const
+ {
+ for (size_type i = 0; i < indices.size(); ++i)
+ values[i] = operator()(indices[i]);
+ }
+
+
+
+ template <typename Number>
+ template <typename ForwardIterator, typename OutputIterator>
+ inline
+ void Vector<Number>::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 <typename Number>
+ template <typename OtherNumber>
+ inline
+ void
+ Vector<Number>::add (const std::vector<size_type> &indices,
+ const ::dealii::Vector<OtherNumber> &values)
+ {
+ AssertDimension (indices.size(), values.size());
+ for (size_type i=0; i<indices.size(); ++i)
+ {
+ Assert (numbers::is_finite(values[i]),
+ ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+ this->operator()(indices[i]) += values(i);
+ }
+ }
+
+
+
+ template <typename Number>
+ inline
+ const MPI_Comm &
+ Vector<Number>::get_mpi_communicator() const
+ {
+ return partitioner->get_mpi_communicator();
+ }
+
+
+
+ template <typename Number>
+ inline
+ std_cxx11::shared_ptr<const Utilities::MPI::Partitioner>
+ Vector<Number>::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 <typename Number>
+inline
+void swap (LinearAlgebra::parallel::Vector<Number> &u,
+ LinearAlgebra::parallel::Vector<Number> &v)
+{
+ u.swap (v);
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
//
// ---------------------------------------------------------------------
-#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 <deal.II/base/config.h>
-#include <deal.II/lac/parallel_vector.h>
-#include <deal.II/lac/vector_view.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/vector_operations_internal.h>
+#include <deal.II/lac/read_write_vector.h>
-#include <deal.II/lac/petsc_parallel_vector.h>
-#include <deal.II/lac/trilinos_vector.h>
DEAL_II_NAMESPACE_OPEN
-namespace parallel
+namespace LinearAlgebra
{
- namespace distributed
+ namespace parallel
{
template <typename Number>
val = 0;
allocated_size = 0;
}
+ thread_loop_partitioner.reset(new ::dealii::parallel::internal::TBBPartitioner());
}
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;
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());
}
vector_is_ghosted = false;
+
+ thread_loop_partitioner = v.thread_loop_partitioner;
}
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());
-#ifdef DEAL_II_WITH_PETSC
+ template <typename Number>
+ Vector<Number>::Vector ()
+ :
+ partitioner (new Utilities::MPI::Partitioner()),
+ allocated_size (0),
+ val (0),
+ import_data (0)
+ {
+ reinit(0);
+ }
- namespace internal
+
+
+ template <typename Number>
+ Vector<Number>::Vector (const Vector<Number> &v)
+ :
+ Subscriptor(),
+ allocated_size (0),
+ val (0),
+ import_data (0),
+ vector_is_ghosted (false)
{
- template <typename PETSC_Number, typename Number>
- 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 <typename PETSC_Number, typename Number>
- void copy_petsc_vector (const std::complex<PETSC_Number> *petsc_start_ptr,
- const std::complex<PETSC_Number> *petsc_end_ptr,
- std::complex<Number> *ptr)
- {
- std::copy(petsc_start_ptr, petsc_end_ptr, ptr);
- }
+ thread_loop_partitioner = v.thread_loop_partitioner;
+ dealii::internal::Vector_copy<Number,Number> copier(v.val, val);
+ internal::parallel_for(copier, partitioner->local_size(),
+ thread_loop_partitioner);
- template <typename PETSC_Number, typename Number>
- void copy_petsc_vector (const std::complex<PETSC_Number> * /*petsc_start_ptr*/,
- const std::complex<PETSC_Number> * /*petsc_end_ptr*/,
- Number * /*ptr*/)
- {
- AssertThrow(false, ExcMessage("Tried to copy complex -> real"));
- }
+ zero_out_ghosts();
}
+
+
template <typename Number>
- Vector<Number> &
- Vector<Number>::operator = (const PETScWrappers::MPI::Vector &petsc_vec)
+ Vector<Number>::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<const Vec &>(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<const Vec &>(petsc_vec), &start_ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
+ template <typename Number>
+ Vector<Number>::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 <typename Number>
+ Vector<Number>::Vector (const size_type size)
+ :
+ allocated_size (0),
+ val (0),
+ import_data (0),
+ vector_is_ghosted (false)
+ {
+ reinit (size, false);
}
-#endif
+ template <typename Number>
+ Vector<Number>::
+ Vector (const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
+ :
+ allocated_size (0),
+ val (0),
+ import_data (0),
+ vector_is_ghosted (false)
+ {
+ reinit (partitioner);
+ }
+
+
+
+ template <typename Number>
+ inline
+ Vector<Number>::~Vector ()
+ {
+ resize_val(0);
+
+ if (import_data != 0)
+ delete[] import_data;
+ import_data = 0;
+
+ clear_mpi_requests();
+ }
+
-#ifdef DEAL_II_WITH_TRILINOS
template <typename Number>
+ inline
Vector<Number> &
- Vector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
+ Vector<Number>::operator = (const Vector<Number> &c)
{
- if (trilinos_vec.has_ghost_elements() == false)
+#ifdef _MSC_VER
+ return this->operator=<Number>(c);
+#else
+ return this->template operator=<Number>(c);
+#endif
+ }
+
+
+
+ template <typename Number>
+ template <typename Number2>
+ inline
+ Vector<Number> &
+ Vector<Number>::operator = (const Vector<Number2> &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<TrilinosWrappers::types::int_type>(this->local_range().first)) : 0;
- const VectorView<double> in_view (local_size(), trilinos_vec.begin()+start_index);
- static_cast<dealii::Vector<Number>&>(vector_view) =
- static_cast<const dealii::Vector<double>&>(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<Number,Number2> 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 <typename Number>
+ void
+ Vector<Number>::compress (::dealii::VectorOperation::values operation)
+ {
+ compress_start (0, operation);
+ compress_finish(operation);
+ }
template <typename Number>
void
- Vector<Number>::copy_from (const Vector<Number> &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<Number>::update_ghost_values () const
+ {
+ update_ghost_values_start ();
+ update_ghost_values_finish ();
+ }
+
+
+
+ template <typename Number>
+ void
+ Vector<Number>::zero_out_ghosts ()
+ {
+ std::fill_n (&val[partitioner->local_size()],
+ partitioner->n_ghost_indices(),
+ Number());
+ vector_is_ghosted = false;
}
}
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<n_ghost_targets; i++)
{
// then actually send the data
if (update_ghost_values_requests.size() == 0)
{
- Assert (part.local_size() == vector_view.size(),
- ExcInternalError());
size_type current_index_start = part.local_size();
update_ghost_values_requests.resize (n_import_targets+n_ghost_targets);
for (unsigned int i=0; i<n_ghost_targets; i++)
+ template <typename Number>
+ void
+ Vector<Number>::import(const ReadWriteVector<Number> &V,
+ VectorOperation::values operation,
+ std_cxx11::shared_ptr<const CommunicationPatternBase> communication_pattern)
+ {
+ // If no communication pattern is given, create one. Otherwise, use the
+ // given one.
+ std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> 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<const Utilities::MPI::Partitioner> (communication_pattern);
+ AssertThrow(comm_pattern != NULL,
+ ExcMessage(std::string("The communication pattern is not of type ") +
+ "Utilities::MPI::Partitioner."));
+ }
+ LinearAlgebra::parallel::Vector<Number> 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<v_stored.n_elements(); ++i)
+ tmp_vector(v_stored.nth_index_in_set(i)) = V.local_element(i);
+
+ tmp_vector.compress(operation);
+
+ dealii::internal::Vector_copy<Number,Number> copier(tmp_vector.val, val);
+ internal::parallel_for(copier, partitioner->local_size(),
+ thread_loop_partitioner);
+ }
+
+
+
template <typename Number>
void
Vector<Number>::swap (Vector<Number> &v)
#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 <typename Number>
+ Vector<Number> &
+ Vector<Number>::operator = (const Number s)
+ {
+ internal::Vector_set<Number> 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 <typename Number>
+ Vector<Number> &
+ Vector<Number>::operator += (const VectorSpaceVector<Number> &vv)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_add_v<Number> 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 <typename Number>
+ Vector<Number> &
+ Vector<Number>::operator -= (const VectorSpaceVector<Number> &vv)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_subtract_v<Number> 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 <typename Number>
+ void
+ Vector<Number>::add (const Number a)
+ {
+ AssertIsFinite(a);
+
+ internal::Vectorization_add_factor<Number> vector_add(val, a);
+ internal::parallel_for(vector_add, partitioner->local_size(),
+ thread_loop_partitioner);
+
+ if (vector_is_ghosted)
+ update_ghost_values();
+ }
+
+
+
+ template <typename Number>
+ void
+ Vector<Number>::add (const Number a,
+ const VectorSpaceVector<Number> &vv)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+
+ AssertIsFinite(a);
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_add_av<Number> 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 <typename Number>
+ void
+ Vector<Number>::add (const Number a,
+ const VectorSpaceVector<Number> &vv,
+ const Number b,
+ const VectorSpaceVector<Number> &ww)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ Assert(dynamic_cast<const Vector<Number> *>(&ww)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &w = dynamic_cast<const Vector<Number> &>(ww);
+
+ AssertIsFinite(a);
+ AssertIsFinite(b);
+
+ AssertDimension (local_size(), v.local_size());
+ AssertDimension (local_size(), w.local_size());
+
+ internal::Vectorization_add_avpbw<Number> 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 <typename Number>
+ void
+ Vector<Number>::sadd (const Number x,
+ const Vector<Number> &v)
+ {
+ AssertIsFinite(x);
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_sadd_xv<Number> 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 <typename Number>
+ void
+ Vector<Number>::sadd (const Number x,
+ const Number a,
+ const VectorSpaceVector<Number> &vv)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+
+ AssertIsFinite(x);
+ AssertIsFinite(a);
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_sadd_xav<Number> 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 <typename Number>
+ void
+ Vector<Number>::sadd (const Number x,
+ const Number a,
+ const Vector<Number> &v,
+ const Number b,
+ const Vector<Number> &w)
+ {
+ AssertIsFinite(x);
+ AssertIsFinite(a);
+ AssertIsFinite(b);
+
+ AssertDimension (local_size(), v.local_size());
+ AssertDimension (local_size(), w.local_size());
+
+ internal::Vectorization_sadd_xavbw<Number> 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 <typename Number>
+ Vector<Number> &
+ Vector<Number>::operator *= (const Number factor)
+ {
+ AssertIsFinite(factor);
+ internal::Vectorization_multiply_factor<Number> 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 <typename Number>
+ Vector<Number> &
+ Vector<Number>::operator /= (const Number factor)
+ {
+ operator *= (1./factor);
+ return *this;
+ }
+
+
+
+ template <typename Number>
+ void
+ Vector<Number>::scale (const VectorSpaceVector<Number> &vv)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_scale<Number> vector_scale(val, v.val);
+ internal::parallel_for(vector_scale, partitioner->local_size(),
+ thread_loop_partitioner);
+
+ if (vector_is_ghosted)
+ update_ghost_values();
+ }
+
+
+
+ template <typename Number>
+ void
+ Vector<Number>::equ (const Number a,
+ const VectorSpaceVector<Number> &vv)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+
+ AssertIsFinite(a);
+ AssertDimension (local_size(), v.local_size());
+
+ internal::Vectorization_equ_au<Number> 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 <typename Number>
+ void
+ Vector<Number>::equ (const Number a,
+ const Vector<Number> &v,
+ const Number b,
+ const Vector<Number> &w)
+ {
+ AssertIsFinite(a);
+ AssertIsFinite(b);
+
+ AssertDimension (local_size(), v.local_size());
+ AssertDimension (local_size(), w.local_size());
+
+ internal::Vectorization_equ_aubv<Number> 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 <typename Number>
+ bool
+ Vector<Number>::all_zero_local () const
+ {
+ const size_type local_size = partitioner->local_size();
+ for (size_type i=0; i<local_size; ++i)
+ if (val[i] != Number(0))
+ return false;
+ return true;
+ }
+
+
+
+ template <typename Number>
+ bool
+ Vector<Number>::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<int>(all_zero_local());
+ if (partitioner->n_mpi_processes() > 1)
+ return -Utilities::MPI::max(local_result,
+ partitioner->get_communicator());
+ else
+ return -local_result;
+ }
+
+
+
+ template <typename Number>
+ template <typename Number2>
+ Number
+ Vector<Number>::inner_product_local(const Vector<Number2> &v) const
+ {
+ if (PointerComparison::equal (this, &v))
+ return norm_sqr_local();
+
+ AssertDimension (partitioner->local_size(), v.partitioner->local_size());
+
+ Number sum;
+ internal::Dot<Number,Number2> dot(val, v.val);
+ internal::parallel_reduce (dot, partitioner->local_size(), sum,
+ thread_loop_partitioner);
+ AssertIsFinite(sum);
+
+ return sum;
+ }
+
+
+
+ template <typename Number>
+ Number
+ Vector<Number>::operator * (const VectorSpaceVector<Number> &vv) const
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(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 Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::norm_sqr_local () const
+ {
+ real_type sum;
+ internal::Norm2<Number,real_type> norm2(val);
+ internal::parallel_reduce (norm2, partitioner->local_size(), sum,
+ thread_loop_partitioner);
+ AssertIsFinite(sum);
+
+ return sum;
+ }
+
+
+
+ template <typename Number>
+ Number
+ Vector<Number>::mean_value_local () const
+ {
+ Assert(size() != 0, ExcEmptyObject());
+
+ if (partitioner->local_size() == 0)
+ return Number();
+
+ Number sum;
+ internal::MeanValue<Number> mean(val);
+ internal::parallel_reduce (mean, partitioner->local_size(), sum,
+ thread_loop_partitioner);
+
+ return sum / real_type(partitioner->local_size());
+ }
+
+
+
+ template <typename Number>
+ Number
+ Vector<Number>::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 Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::l1_norm_local () const
+ {
+ real_type sum;
+ internal::Norm1<Number, real_type> norm1(val);
+ internal::parallel_reduce (norm1, partitioner->local_size(), sum,
+ thread_loop_partitioner);
+
+ return sum;
+ }
+
+
+
+ template <typename Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::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 Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::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 Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::lp_norm_local (const real_type p) const
+ {
+ real_type sum;
+ internal::NormP<Number, real_type> normp(val, p);
+ internal::parallel_reduce (normp, partitioner->local_size(), sum,
+ thread_loop_partitioner);
+ return std::pow(sum, 1./p);
+ }
+
+
+
+ template <typename Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::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<real_type>(1.0/p));
+ else
+ return local_result;
+ }
+
+
+
+ template <typename Number>
+ typename Vector<Number>::real_type
+ Vector<Number>::linfty_norm_local () const
+ {
+ real_type max = 0.;
+
+ const size_type local_size = partitioner->local_size();
+ for (size_type i=0; i<local_size; ++i)
+ max = std::max (numbers::NumberTraits<Number>::abs(val[i]), max);
+
+ return max;
+ }
+
+
+
+ template <typename Number>
+ inline
+ typename Vector<Number>::real_type
+ Vector<Number>::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 <typename Number>
+ Number
+ Vector<Number>::add_and_dot_local(const Number a,
+ const Vector<Number> &v,
+ const Vector<Number> &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<Number> adder(this->val, v.val, w.val, a);
+ internal::parallel_reduce (adder, vec_size, sum, thread_loop_partitioner);
+ AssertIsFinite(sum);
+ return sum;
+ }
+
+
+
+ template <typename Number>
+ Number
+ Vector<Number>::add_and_dot (const Number a,
+ const VectorSpaceVector<Number> &vv,
+ const VectorSpaceVector<Number> &ww)
+ {
+ // Downcast. Throws an exception if invalid.
+ Assert(dynamic_cast<const Vector<Number> *>(&vv)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ Assert(dynamic_cast<const Vector<Number> *>(&ww)!=NULL,
+ ExcVectorTypeNotCompatible());
+ const Vector<Number> &w = dynamic_cast<const Vector<Number> &>(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;
}
#define dealii__parallel_vector_h
#include <deal.II/base/config.h>
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/template_constraints.h>
-#include <deal.II/base/types.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/base/partitioner.h>
-#include <deal.II/base/thread_management.h>
-#include <deal.II/lac/vector_view.h>
+#include <deal.II/lac/la_parallel_vector.h>
#include <cstring>
#include <iomanip>
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
{
*
* @author Katharina Kormann, Martin Kronbichler, 2010, 2011
*/
- template <typename Number>
- class Vector : public Subscriptor
- {
- public:
- /**
- * Declare standard types used in all containers. These types parallel
- * those in the <tt>C++</tt> standard libraries <tt>vector<...></tt>
- * 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<Number>::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<Number> &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<Number> 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<const Utilities::MPI::Partitioner> &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 <typename Number2>
- void reinit(const Vector<Number2> &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<Number> 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<const Utilities::MPI::Partitioner> &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
- * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
- * analogy to standard functions.
- *
- * This function is virtual in order to allow for derived classes to
- * handle memory separately.
- */
- void swap (Vector<Number> &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<Number> &
- operator = (const Vector<Number> &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 <typename Number2>
- Vector<Number> &
- operator = (const Vector<Number2> &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<Number> &
- 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<Number> &
- 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<Number> &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<Number> &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 <typename Number2>
- bool operator == (const Vector<Number2> &v) const;
-
- /**
- * Checks for inequality of the two vectors.
- */
- template <typename Number2>
- bool operator != (const Vector<Number2> &v) const;
-
- /**
- * Perform the inner product of two vectors.
- */
- template <typename Number2>
- Number operator * (const Vector<Number2> &V) const;
-
- /**
- * Computes the square of the l<sub>2</sub> norm of the vector (i.e.,
- * the sum of the squares of all entries among all processors).
- */
- real_type norm_sqr () const;
-
- /**
- * Computes the mean value of all the entries in the vector.
- */
- Number mean_value () const;
-
- /**
- * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
- * absolute values of all entries among all processors).
- */
- real_type l1_norm () const;
-
- /**
- * Returns the l<sub>2</sub> norm of the vector (i.e., square root of
- * the sum of the square of all entries among all processors).
- */
- real_type l2_norm () const;
-
- /**
- * Returns the l<sub>p</sub> norm with real @p p of the vector (i.e.,
- * the pth root of sum of the pth power of all entries among all
- * processors).
- */
- real_type lp_norm (const real_type p) const;
-
- /**
- * Returns the maximum norm of the vector (i.e., maximum absolute value
- * among all entries among all processors).
- */
- real_type linfty_norm () const;
-
- /**
- * Performs a combined operation of a vector addition and a subsequent
- * inner product, returning the value of the inner product. In other
- * words, the result of this function is the same as if the user called
- * @code
- * this->add(a, V);
- * return_value = *this * W;
- * @endcode
- *
- * The reason this function exists is that this operation involves less
- * memory transfer than calling the two functions separately. This
- * method only needs to load three vectors, @p this, @p V, @p W, whereas
- * calling separate methods means to load the calling vector @p this
- * twice. Since most vector operations are memory transfer limited, this
- * reduces the time by 25\% (or 50\% if @p W equals @p this).
- */
- Number add_and_dot (const Number a,
- const Vector<Number> &V,
- const Vector<Number> &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 <code>local_size() == local_range().second -
- * local_range().first</code>.
- */
- std::pair<size_type, size_type> 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 <tt>vector<></tt> class of
- * the C++ standard library by returning iterators to the start and end
- * of the <i>locally owned</i> 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: <tt>O(1)</tt> for locally owned elements that represent
- * a contiguous range and <tt>O(log(n<sub>ranges</sub>))</tt> 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: <tt>O(1)</tt> for locally owned elements that represent
- * a contiguous range and <tt>O(log(n<sub>ranges</sub>))</tt> 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 <typename OtherNumber>
- void extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<OtherNumber> &values) const;
-
- /**
- * Just as the above, but with pointers. Useful in minimizing copying of
- * data around.
- */
- template <typename ForwardIterator, typename OutputIterator>
- 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
- * <code>[0,local_size)</code>, and ghost indices with indices
- * <code>[local_size,local_size+ n_ghost_entries]</code>.
- *
- * 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
- * <code>[0,local_size)</code>, and ghost indices with indices
- * <code>[local_size,local_size+n_ghosts]</code>.
- *
- * 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<Number> &operator += (const Vector<Number> &V);
-
- /**
- * Subtract the given vector from the present one.
- */
- Vector<Number> &operator -= (const Vector<Number> &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 <typename OtherNumber>
- void add (const std::vector<size_type> &indices,
- const std::vector<OtherNumber> &values);
-
- /**
- * This is a second collective add operation. As a difference, this
- * function takes a deal.II vector of values.
- */
- template <typename OtherNumber>
- void add (const std::vector<size_type> &indices,
- const ::dealii::Vector<OtherNumber> &values);
-
- /**
- * Take an address where <tt>n_elements</tt> are stored contiguously and
- * add them into the vector. Handles all cases which are not covered by
- * the other two <tt>add()</tt> functions above.
- */
- template <typename OtherNumber>
- 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 <tt>operator +=</tt>.
- *
- * @deprecated Use the <tt>operator +=</tt> instead.
- */
- void add (const Vector<Number> &V) DEAL_II_DEPRECATED;
-
- /**
- * Simple addition of a multiple of a vector, i.e. <tt>*this +=
- * a*V</tt>.
- */
- void add (const Number a, const Vector<Number> &V);
-
- /**
- * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
- */
- void add (const Number a, const Vector<Number> &V,
- const Number b, const Vector<Number> &W);
-
- /**
- * Scaling and simple vector addition, i.e. <tt>*this =
- * s*(*this)+V</tt>.
- */
- void sadd (const Number s,
- const Vector<Number> &V);
-
- /**
- * Scaling and simple addition, i.e. <tt>*this = s*(*this)+a*V</tt>.
- */
- void sadd (const Number s,
- const Number a,
- const Vector<Number> &V);
-
- /**
- * Scaling and multiple addition.
- *
- * This function is deprecated.
- */
- void sadd (const Number s,
- const Number a,
- const Vector<Number> &V,
- const Number b,
- const Vector<Number> &W) DEAL_II_DEPRECATED;
-
- /**
- * Scaling and multiple addition. <tt>*this = s*(*this)+a*V + b*W +
- * c*X</tt>.
- *
- * This function is deprecated.
- */
- void sadd (const Number s,
- const Number a,
- const Vector<Number> &V,
- const Number b,
- const Vector<Number> &W,
- const Number c,
- const Vector<Number> &X) DEAL_II_DEPRECATED;
-
- /**
- * Scale each element of the vector by a constant value.
- */
- Vector<Number> &operator *= (const Number factor);
-
- /**
- * Scale each element of the vector by the inverse of the given value.
- */
- Vector<Number> &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<Number> &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 <typename Number2>
- void scale (const Vector<Number2> &scaling_factors);
-
- /**
- * Assignment <tt>*this = a*u</tt>.
- */
- void equ (const Number a, const Vector<Number> &u);
-
- /**
- * Assignment <tt>*this = a*u</tt>.
- */
- template <typename Number2>
- void equ (const Number a, const Vector<Number2> &u);
-
- /**
- * Assignment <tt>*this = a*u + b*v</tt>.
- *
- * This function is deprecated.
- */
- void equ (const Number a, const Vector<Number> &u,
- const Number b, const Vector<Number> &v) DEAL_II_DEPRECATED;
-
- /**
- * Assignment <tt>*this = a*u + b*v + b*w</tt>.
- *
- * This function is deprecated.
- */
- void equ (const Number a, const Vector<Number> &u,
- const Number b, const Vector<Number> &v,
- const Number c, const Vector<Number> &w) DEAL_II_DEPRECATED;
-
- /**
- * Compute the elementwise ratio of the two given vectors, that is let
- * <tt>this[i] = a[i]/b[i]</tt>. 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 <tt>b[i]</tt> is zero, the result is undefined. No
- * attempt is made to catch such situations.
- */
- void ratio (const Vector<Number> &a,
- const Vector<Number> &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<const Utilities::MPI::Partitioner>
- 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 <typename Number2>
- bool vectors_equal_local (const Vector<Number2> &v) const;
-
- /**
- * Local part of the inner product of two vectors.
- */
- template <typename Number2>
- Number inner_product_local (const Vector<Number2> &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<Number> &V,
- const Vector<Number> &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<const Utilities::MPI::Partitioner> 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<Number> 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<MPI_Request> compress_requests;
-
- /**
- * A vector that collects all requests from @p update_ghost_values()
- * operations. This class uses persistent MPI communicators.
- */
- mutable std::vector<MPI_Request> 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 <typename Number2> friend class Vector;
-
- /**
- * Make BlockVector type friends.
- */
- template <typename Number2> friend class BlockVector;
- };
-
- /*@}*/
-
-
- /*----------------------- Inline functions ----------------------------------*/
-
-#ifndef DOXYGEN
-
- template <typename Number>
- inline
- Vector<Number>::Vector ()
- :
- partitioner (new Utilities::MPI::Partitioner()),
- allocated_size (0),
- val (0),
- import_data (0),
- vector_is_ghosted (false),
- vector_view (0, static_cast<Number *>(0))
- {}
-
-
-
- template <typename Number>
- inline
- Vector<Number>::Vector (const Vector<Number> &v)
- :
- Subscriptor(),
- allocated_size (0),
- val (0),
- import_data (0),
- vector_is_ghosted (false),
- vector_view (0, static_cast<Number *>(0))
- {
- reinit (v, true);
- vector_view = v.vector_view;
- zero_out_ghosts();
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number>::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<Number *>(0))
- {
- reinit (local_range, ghost_indices, communicator);
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number>::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<Number *>(0))
- {
- reinit (local_range, communicator);
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number>::Vector (const size_type size)
- :
- allocated_size (0),
- val (0),
- import_data (0),
- vector_is_ghosted (false),
- vector_view (0, static_cast<Number *>(0))
- {
- reinit (size, false);
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number>::
- Vector (const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
- :
- allocated_size (0),
- val (0),
- import_data (0),
- vector_is_ghosted (false),
- vector_view (0, static_cast<Number *>(0))
- {
- reinit (partitioner);
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number>::~Vector ()
- {
- resize_val(0);
-
- if (import_data != 0)
- delete[] import_data;
- import_data = 0;
-
- clear_mpi_requests();
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number> &
- Vector<Number>::operator = (const Vector<Number> &c)
- {
-#ifdef _MSC_VER
- return this->operator=<Number>(c);
-#else
- return this->template operator=<Number>(c);
-#endif
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- Vector<Number> &
- Vector<Number>::operator = (const Vector<Number2> &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<dealii::Vector<Number> &>(vector_view) = c.vector_view;
-
- if (must_update_ghost_values)
- update_ghost_values();
- else
- zero_out_ghosts();
- return *this;
- }
-
-
-
- template <typename Number>
- inline
- void
- Vector<Number>::compress (::dealii::VectorOperation::values operation)
- {
- compress_start (0, operation);
- compress_finish(operation);
- }
-
-
-
- template <typename Number>
- inline
- void
- Vector<Number>::update_ghost_values () const
- {
- update_ghost_values_start ();
- update_ghost_values_finish ();
- }
-
-
-
- template <typename Number>
- inline
- void
- Vector<Number>::zero_out_ghosts ()
- {
- std::fill_n (&val[partitioner->local_size()],
- partitioner->n_ghost_indices(),
- Number());
- vector_is_ghosted = false;
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::has_ghost_elements () const
- {
- return vector_is_ghosted;
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::all_zero_local () const
- {
- return partitioner->local_size()>0 ? vector_view.all_zero () : true;
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::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<int>(all_zero_local());
- if (partitioner->n_mpi_processes() > 1)
- return -Utilities::MPI::max(local_result,
- partitioner->get_communicator());
- else
- return -local_result;
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::is_non_negative_local () const
- {
- return partitioner->local_size()>0 ? vector_view.is_non_negative () : true;
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::is_non_negative () const
- {
- int local_result = -static_cast<int>(is_non_negative_local());
- if (partitioner->n_mpi_processes() > 1)
- return -Utilities::MPI::max(local_result,
- partitioner->get_communicator());
- else
- return -local_result;
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- bool
- Vector<Number>::vectors_equal_local (const Vector<Number2> &v) const
- {
- return partitioner->local_size()>0 ?
- vector_view.template operator == <Number2>(v.vector_view)
- : true;
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- bool
- Vector<Number>::operator == (const Vector<Number2> &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<int>(!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 <typename Number>
- template <typename Number2>
- inline
- bool
- Vector<Number>::operator != (const Vector<Number2> &v) const
- {
- return !(operator == (v));
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- Number
- Vector<Number>::inner_product_local(const Vector<Number2> &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 <typename Number>
- template <typename Number2>
- inline
- Number
- Vector<Number>::operator * (const Vector<Number2> &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 <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::norm_sqr_local () const
- {
- return partitioner->local_size()>0 ? vector_view.norm_sqr() : real_type();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::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 <typename Number>
- inline
- Number
- Vector<Number>::mean_value_local () const
- {
- Assert (partitioner->size()!=0, ExcEmptyObject());
- return (partitioner->local_size() ?
- vector_view.mean_value()
- : Number());
- }
-
-
-
- template <typename Number>
- inline
- Number
- Vector<Number>::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 Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::l1_norm_local () const
- {
- return partitioner->local_size() ? vector_view.l1_norm() : real_type();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::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 Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::l2_norm () const
- {
- return std::sqrt(norm_sqr());
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::lp_norm_local (const real_type p) const
- {
- return partitioner->local_size() ? vector_view.lp_norm(p) : real_type();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::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<real_type>(1.0/p));
- else
- return local_result;
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::linfty_norm_local () const
- {
- return partitioner->local_size() ? vector_view.linfty_norm() : real_type();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::real_type
- Vector<Number>::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 <typename Number>
- inline
- Number
- Vector<Number>::add_and_dot_local(const Number a,
- const Vector<Number> &V,
- const Vector<Number> &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 <typename Number>
- inline
- Number
- Vector<Number>::add_and_dot (const Number a,
- const Vector<Number> &V,
- const Vector<Number> &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 <typename Number>
- inline
- typename Vector<Number>::size_type
- Vector<Number>::size () const
- {
- return partitioner->size();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::size_type
- Vector<Number>::local_size () const
- {
- return partitioner->local_size();
- }
-
-
-
- template <typename Number>
- inline
- std::pair<typename Vector<Number>::size_type,
- typename Vector<Number>::size_type>
- Vector<Number>::local_range () const
- {
- return partitioner->local_range();
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::in_local_range
- (const size_type global_index) const
- {
- return partitioner->in_local_range (global_index);
- }
-
-
-
- template <typename Number>
- inline
- IndexSet
- Vector<Number>::locally_owned_elements() const
- {
- IndexSet is (size());
-
- is.add_range (local_range().first, local_range().second);
-
- return is;
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::size_type
- Vector<Number>::n_ghost_entries () const
- {
- return partitioner->n_ghost_indices();
- }
-
-
-
- template <typename Number>
- inline
- const IndexSet &
- Vector<Number>::ghost_elements() const
- {
- return partitioner->ghost_indices();
- }
-
-
-
- template <typename Number>
- inline
- bool
- Vector<Number>::is_ghost_entry (const size_type global_index) const
- {
- return partitioner->is_ghost_entry (global_index);
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::iterator
- Vector<Number>::begin ()
- {
- return vector_view.begin();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::const_iterator
- Vector<Number>::begin () const
- {
- return vector_view.begin();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::iterator
- Vector<Number>::end ()
- {
- return vector_view.end();
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<Number>::const_iterator
- Vector<Number>::end () const
- {
- return vector_view.end();
- }
-
-
-
- template <typename Number>
- inline
- Number
- Vector<Number>::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 <typename Number>
- inline
- Number &
- Vector<Number>::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 <typename Number>
- inline
- Number
- Vector<Number>::operator[] (const size_type global_index) const
- {
- return operator()(global_index);
- }
-
-
-
- template <typename Number>
- inline
- Number &
- Vector<Number>::operator[] (const size_type global_index)
- {
- return operator()(global_index);
- }
-
-
-
- template <typename Number>
- template <typename OtherNumber>
- inline
- void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<OtherNumber> &values) const
- {
- for (size_type i = 0; i < indices.size(); ++i)
- values[i] = operator()(indices[i]);
- }
-
-
-
- template <typename Number>
- template <typename ForwardIterator, typename OutputIterator>
- inline
- void Vector<Number>::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 <typename Number>
- inline
- Number
- Vector<Number>::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 <typename Number>
- inline
- Number &
- Vector<Number>::local_element (const size_type local_index)
- {
- AssertIndexRange (local_index,
- partitioner->local_size()+
- partitioner->n_ghost_indices());
- return val[local_index];
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number> &
- Vector<Number>::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<Number>::operator= (s);
- if (s==Number())
- zero_out_ghosts();
-
- return *this;
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number> &
- Vector<Number>::operator += (const Vector<Number> &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 <typename Number>
- inline
- Vector<Number> &
- Vector<Number>::operator -= (const Vector<Number> &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 <typename Number>
- template <typename OtherNumber>
- inline
- void
- Vector<Number>::add (const std::vector<size_type> &indices,
- const std::vector<OtherNumber> &values)
- {
- AssertDimension (indices.size(), values.size());
- add (indices.size(), &indices[0], &values[0]);
- }
-
-
-
- template <typename Number>
- template <typename OtherNumber>
- inline
- void
- Vector<Number>::add (const std::vector<size_type> &indices,
- const ::dealii::Vector<OtherNumber> &values)
- {
- AssertDimension (indices.size(), values.size());
- add (indices.size(), &indices[0], values.begin());
- }
-
-
-
- template <typename Number>
- template <typename OtherNumber>
- inline
- void
- Vector<Number>::add (const size_type n_indices,
- const size_type *indices,
- const OtherNumber *values)
- {
- for (size_type i=0; i<n_indices; ++i)
- {
- Assert (numbers::is_finite(values[i]),
- ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
- this->operator()(indices[i]) += values[i];
- }
- }
-
-
-
- template <typename Number>
- inline
- void
- Vector<Number>::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 <typename Number>
- inline
- void
- Vector<Number>::add (const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::add (const Number a,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::add (const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::sadd (const Number x,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::sadd (const Number x,
- const Number a,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::sadd (const Number x,
- const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::sadd (const Number s,
- const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &w,
- const Number c,
- const Vector<Number> &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 <typename Number>
- inline
- Vector<Number> &
- Vector<Number>::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 <typename Number>
- inline
- Vector<Number> &
- Vector<Number>::operator /= (const Number factor)
- {
- operator *= (1./factor);
- return *this;
- }
-
-
-
- template <typename Number>
- inline
- void
- Vector<Number>::scale (const Vector<Number> &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 <typename Number>
- template <typename Number2>
- inline
- void
- Vector<Number>::scale (const Vector<Number2> &scaling_factors)
- {
- if (local_size())
- vector_view.template scale<Number2> (scaling_factors.vector_view);
-
- if (vector_is_ghosted)
- update_ghost_values();
- }
-
-
-
- template <typename Number>
- inline
- void
- Vector<Number>::equ (const Number a,
- const Vector<Number> &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 <typename Number>
- template <typename Number2>
- inline
- void
- Vector<Number>::equ (const Number a,
- const Vector<Number2> &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 <typename Number>
- inline
- void
- Vector<Number>::equ (const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::equ (const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &w,
- const Number c,
- const Vector<Number> &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 <typename Number>
- inline
- void
- Vector<Number>::ratio (const Vector<Number> &a,
- const Vector<Number> &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 <typename Number>
- inline
- const MPI_Comm &
- Vector<Number>::get_mpi_communicator() const
- {
- return partitioner->get_communicator();
- }
-
-
-
- template <typename Number>
- inline
- std_cxx11::shared_ptr<const Utilities::MPI::Partitioner>
- Vector<Number>::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 <typename Number>
-inline
-void swap (parallel::distributed::Vector<Number> &u,
- parallel::distributed::Vector<Number> &v)
-{
- u.swap (v);
+ using LinearAlgebra::parallel::Vector;
+ }
}
-
DEAL_II_NAMESPACE_CLOSE
#endif
template <typename number> class Vector;
template <typename number> class SparseMatrix;
-namespace parallel
+namespace LinearAlgebra
{
- namespace distributed
+ namespace parallel
{
template <typename number> class Vector;
}
template <typename Number>
inline
void
- vector_updates (const parallel::distributed::Vector<Number> &src,
- const parallel::distributed::Vector<Number> &matrix_diagonal_inverse,
+ vector_updates (const LinearAlgebra::parallel::Vector<Number> &src,
+ const LinearAlgebra::parallel::Vector<Number> &matrix_diagonal_inverse,
const bool start_zero,
const double factor1,
const double factor2,
- parallel::distributed::Vector<Number> &update1,
- parallel::distributed::Vector<Number> &update2,
- parallel::distributed::Vector<Number> &dst)
+ LinearAlgebra::parallel::Vector<Number> &update1,
+ LinearAlgebra::parallel::Vector<Number> &update2,
+ LinearAlgebra::parallel::Vector<Number> &dst)
{
VectorUpdater<Number> upd(src.begin(), matrix_diagonal_inverse.begin(),
start_zero, factor1, factor2,
* For parallel loops with TBB, this member variable stores the affinity
* information of loops.
*/
- mutable std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> thread_loop_partitioner;
+ mutable std_cxx11::shared_ptr<::dealii::parallel::internal::TBBPartitioner> thread_loop_partitioner;
/**
* Make all other ReadWriteVector types friends.
// 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,
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);
DEAL_II_NAMESPACE_OPEN
-class CommunicationPatternBase;
class IndexSet;
-template <typename Number> class ReadWriteVector;
+namespace LinearAlgebra
+{
+ class CommunicationPatternBase;
+ template <typename Number> class ReadWriteVector;
+}
namespace LinearAlgebra
{
// forward declarations
-namespace parallel
+namespace LinearAlgebra
{
- namespace distributed
+ namespace parallel
{
template <typename> class Vector;
}
template <typename Number>
inline
Number &
- vector_access (parallel::distributed::Vector<Number> &vec,
+ vector_access (LinearAlgebra::parallel::Vector<Number> &vec,
const unsigned int entry)
{
return vec.local_element(entry);
template <typename Number>
inline
Number
- vector_access (const parallel::distributed::Vector<Number> &vec,
+ vector_access (const LinearAlgebra::parallel::Vector<Number> &vec,
const unsigned int entry)
{
return vec.local_element(entry);
template <typename Number>
inline
- void check_vector_compatibility (const parallel::distributed::Vector<Number> &vec,
+ void check_vector_compatibility (const LinearAlgebra::parallel::Vector<Number> &vec,
const internal::MatrixFreeFunctions::DoFInfo &dof_info)
{
Assert (vec.partitioners_are_compatible(*dof_info.vector_partitioner),
+ 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)
{
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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_vector.templates.h>
+
+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<S1>& Vector<S1>::operator=<S2> (const Vector<S2> &)
+
+#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
+ TEMPL_COPY_CONSTRUCTOR(double,float);
+ TEMPL_COPY_CONSTRUCTOR(float,double);
+
+
+ TEMPL_COPY_CONSTRUCTOR(std::complex<double>,std::complex<float>);
+ TEMPL_COPY_CONSTRUCTOR(std::complex<float>,std::complex<double>);
+
+#endif
+
+#undef TEMPL_COPY_CONSTRUCTOR
+ }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
for (SCALAR : REAL_SCALARS)
{
- namespace parallel
+ namespace LinearAlgebra
\{
- namespace distributed
+ namespace parallel
\{
template class Vector<SCALAR>;
\}
for (S1, S2 : REAL_SCALARS)
{
- namespace parallel
+ namespace LinearAlgebra
\{
- namespace distributed
+ namespace parallel
\{
template void Vector<S1>::reinit<S2> (const Vector<S2>&,
const bool);
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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 <deal.II/lac/parallel_vector.h>
-#include <deal.II/lac/parallel_vector.templates.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-#include "parallel_vector.inst"
-
-
-DEAL_II_NAMESPACE_CLOSE
}
// 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<double> v2;
deallog << "Inner product: " << inner_prod << std::endl;
}
- // check operator ==
- {
- parallel::distributed::Vector<double> 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();
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<double> 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;
}
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
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
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
for (unsigned i=0; i<local_size; ++i)
w.local_element(i) = 2.0 * (i + my_start);
- v.copy_from(w);
+ v = w;
v.update_ghost_values();
// check local values for correctness
AssertThrow (v(ghost_indices[i])==v.local_element(local_size+i-5), ExcInternalError());
- // now the same again, but import ghosts
- // through the call to copy_from
+ // now the same again, but import ghosts automatically because v had ghosts
+ // set before calling operator =
v.reinit (local_owned, local_relevant, MPI_COMM_WORLD);
- v.copy_from(w, true);
+ v.update_ghost_values();
+ v = w;
// check local values for correctness
for (unsigned int i=0; i<local_size; ++i)
for (unsigned int i=5; i<10; ++i)
AssertThrow (v(ghost_indices[i])==v.local_element(local_size+i-5), ExcInternalError());
-
- // now do not call import_ghosts and check
- // whether ghosts really are zero
- v.reinit (local_owned, local_relevant, MPI_COMM_WORLD);
- v.copy_from(w, false);
-
- // check non-local entries on all processors
- for (unsigned int i=0; i<10; ++i)
- if (local_owned.is_element (ghost_indices[i]) == false)
- AssertThrow (v(ghost_indices[i]) == 0., ExcInternalError());
-
if (myid == 0)
deallog << "OK" << std::endl;
}
if (myid==0) deallog << "OK" << std::endl;
if (myid==0) deallog << "Check add (vector): ";
- y.add (w);
+ y.add (1., w);
for (int i=0; i<actual_local_size; ++i)
AssertThrow (y.local_element(i) == 3*(i+my_start)+1042, ExcInternalError());
if (myid==0) deallog << "OK" << std::endl;
if (myid==0) deallog << "OK" << std::endl;
if (myid==0) deallog << "Check sadd (factor, factor, vector, factor, vector, factor, vector): ";
- y.sadd (-1.,1.,v, 2., w, 2., x);
+ y.sadd (-1.,1.,v, 2., w);
+ y.add(2., x);
for (int i=0; i<actual_local_size; ++i)
AssertThrow (y.local_element(i) == 20000, ExcInternalError());
if (myid==0) deallog << "OK" << std::endl;
if (myid==0) deallog << "OK" << std::endl;
if (myid==0) deallog << "Check equ (factor, vector, factor, vector, factor, vector): ";
- y. equ (10., v, -2., w, 3., x);
+ y. equ (10., v, -2., w);
+ y. add (3., x);
for (int i=0; i<actual_local_size; ++i)
AssertThrow (y.local_element(i) == 6.*(i+my_start)+28000, ExcInternalError());
if (myid==0) deallog << "OK" << std::endl;
if (myid==0) deallog << "Check equ<float> (factor, vector): ";
parallel::distributed::Vector<float> z;
z = v;
- y.equ (1., z);
+ y = z;
for (int i=0; i<actual_local_size; ++i)
AssertThrow (y.local_element(i) == i+my_start, ExcInternalError());
if (myid==0) deallog << "OK" << std::endl;
projects / dealii.git / commitdiff