From: Martin Kronbichler Date: Tue, 21 May 2013 13:34:04 +0000 (+0000) Subject: Fix bug in parallel::distributed::Vector::compress(VectorOperation::insert) pointed... X-Git-Tag: v8.0.0~502 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e2a9bd58a34fcd6923ef9a3e3f32c78814a57551;p=dealii.git Fix bug in parallel::distributed::Vector::compress(VectorOperation::insert) pointed out by Martin Steigemann. git-svn-id: https://svn.dealii.org/trunk@29534 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 3159feed5a..1241a2ef19 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -128,6 +128,18 @@ this function.
    +
  1. Fixed: The method parallel::distributed::Vector::compress +(VectorOperation::insert) previously set the elements of ghost elements +unconditionally on the owning processor, even if they had not been touched. +This led to a problem in certain library functions where vector entries became +zero in a spurious way. This is now fixed by discarding the elements in ghost +entries for the VectorOperation::insert operation. This is legitimate since we +assume consistency of set elements across processors, so the owning processor +sets the element already. +
    +(Martin Kronbichler, 2013/05/21) +
    2. +
  2. Improved: DoFTools::make_periodicity_constraints now also works for meshes where the refinement level of the two sides of the domain is not the same, i.e., one side is more refined than the other. diff --git a/deal.II/include/deal.II/lac/parallel_vector.h b/deal.II/include/deal.II/lac/parallel_vector.h index 64cade9c6f..9440ccc1d5 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.h +++ b/deal.II/include/deal.II/lac/parallel_vector.h @@ -85,10 +85,9 @@ namespace parallel { public: /** - * Declare standard types used in all - * containers. These types parallel those in - * the C++ standard libraries - * vector<...> class. + * Declare standard types used in all containers. These types parallel + * those in the C++ standard libraries vector<...> + * class. */ typedef Number value_type; typedef value_type *pointer; @@ -110,35 +109,27 @@ namespace parallel Vector (); /** - * Copy constructor. Uses the parallel - * partitioning of @p in_vector. + * Copy constructor. Uses the parallel partitioning of @p in_vector. */ Vector (const Vector &in_vector); /** - * Constructs a parallel vector of the given - * global size without any actual parallel - * distribution. + * Constructs a parallel vector of the given global size without any + * actual parallel distribution. */ Vector (const unsigned int 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. + * Constructs a parallel vector. The local range is specified by @p + * locally_owned_set (note that this must be a contiguous interval, + * multiple intervals are not possible). The IndexSet @p ghost_indices + * specifies ghost indices, i.e., indices which one might need to read + * data from or accumulate data from. It is allowed that the set of + * ghost indices also contains the local range, but it does not need to. * - * This function involves global - * communication, so it should only be called - * once for a given layout. Use the - * constructor with Vector argument to - * create additional vectors with the same + * This function involves global communication, so it should only be + * called once for a given layout. Use the constructor with + * Vector argument to create additional vectors with the same * parallel layout. */ Vector (const IndexSet &local_range, @@ -146,12 +137,10 @@ namespace parallel 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. + * 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_cxx1x::shared_ptr &partitioner); @@ -161,47 +150,36 @@ namespace parallel ~Vector (); /** - * Sets the global size of the vector to @p - * size without any actual parallel - * distribution. + * Sets the global size of the vector to @p size without any actual + * parallel distribution. */ void reinit (const unsigned int size, const bool fast = 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. + * 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 fast 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). + * If the flag @p fast is set to false, the memory will be initialized + * with zero, otherwise the memory will be untouched (and the user must + * make sure to fill it with reasonable data before using it). */ template void reinit(const Vector &in_vector, const bool fast = 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. + * Initialize the vector. The local range is specified by @p + * locally_owned_set (note that this must be a contiguous interval, + * multiple intervals are not possible). The IndexSet @p ghost_indices + * specifies ghost indices, i.e., indices which one might need to read + * data from or accumulate data from. It is allowed that the set of + * ghost indices also contains the local range, but it does not need to. * - * This function involves global - * communication, so it should only be called - * once for a given layout. Use the @p reinit - * function with Vector argument to - * create additional vectors with the same + * This function involves global communication, so it should only be + * called once for a given layout. Use the @p reinit function with + * Vector argument to create additional vectors with the same * parallel layout. */ void reinit (const IndexSet &local_range, @@ -209,89 +187,81 @@ namespace parallel 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. + * 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_cxx1x::shared_ptr &partitioner); /** - * Swap the contents of this - * vector and the other vector - * @p v. One could do this - * operation with a temporary - * variable and copying over the - * data elements, but this - * function is significantly more - * efficient since it only swaps - * the pointers to the data of - * the two vectors and therefore - * does not need to allocate - * temporary storage and move - * data around. + * Swap the contents of this vector and the other vector @p v. One could + * do this operation with a temporary variable and copying over the data + * elements, but this function is significantly more efficient since it + * only swaps the pointers to the data of the two vectors and therefore + * does not need to allocate temporary storage and move data around. * - * This function is analog to the - * the @p swap function of all C++ - * standard containers. Also, - * there is a global function - * swap(u,v) that simply calls - * u.swap(v), again in analogy - * to standard functions. + * This function is analog to the the @p swap function of all C++ + * standard containers. Also, there is a global function + * swap(u,v) that simply calls u.swap(v), again in + * analogy to standard functions. * - * This function is virtual in - * order to allow for derived - * classes to handle memory - * separately. + * This function is virtual in order to allow for derived classes to + * handle memory separately. */ void swap (Vector &v); /** - * Assigns the vector to the parallel - * partitioning of the input vector @p - * in_vector, and copies all the data. + * Assigns the vector to the parallel partitioning of the input vector + * @p in_vector, and copies all the data. */ Vector & operator = (const Vector &in_vector); /** - * Assigns the vector to the parallel - * partitioning of the input vector @p - * in_vector, and copies all the data. + * Assigns the vector to the parallel partitioning of the input vector + * @p in_vector, and copies all the data. */ template Vector & operator = (const Vector &in_vector); /** - * This method copies the local range from - * another vector with the same local range, - * but possibly different layout of ghost - * indices. + * This method copies the local range from another vector with the same + * local range, but possibly different layout of ghost indices. */ void copy_from (const Vector &in_vector, const bool call_update_ghost_values = false); /** - * 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. + * Sets all elements of the vector to the scalar @p s. If the scalar is + * zero, also ghost elements are set to zero, otherwise they remain + * unchanged. */ Vector &operator = (const Number s); /** - * This function copies the data that has - * accumulated in the data buffer for ghost - * indices to the owning processor. + * This function copies the data that has accumulated in the data buffer + * for ghost indices to the owning processor. For the meaning of the + * argument @p operation, see the entry on @ref GlossCompress + * "Compressing distributed vectors and matrices" in the glossary. + * + * There are two variants for this function. If called with argument @p + * VectorOperation::add adds all the data accumulated in ghost elements + * to the respective elements on the owning processor and clears the + * ghost array afterwards. If called with argument @p + * VectorOperation::insert, a set operation is performed. Since setting + * elements in a vector with ghost elements is ambiguous (as one can set + * both the element on the ghost site as well as the owning site), this + * operation makes the assumption that all data is set correctly on the + * owning processor. Upon call of compress(VectorOperation::insert), all + * ghost entries are therefore simply zeroed out (using + * zero_ghost_values()). In debug mode, a check is performed that makes + * sure that the data set is actually consistent between processors, + * i.e., whenever a non-zero ghost element is found, it is compared to + * the value on the owning processor and an exception is thrown if these + * elements do not agree. * - * For the meaning of this argument, - * see the entry on @ref - * GlossCompress "Compressing - * distributed vectors and matrices" - * in the glossary. */ void compress (::dealii::VectorOperation::values operation); @@ -301,121 +271,97 @@ namespace parallel void compress () DEAL_II_DEPRECATED; /** - * 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. + * Fills the data field for ghost indices with the values stored in the + * respective positions of the owning processor. This function is needed + * before reading from ghosts. The function is @p const even though + * ghost data is changed. This is needed to allow functions with a @p + * const vector to perform the data exchange without creating + * temporaries. */ void update_ghost_values () const; /** - * 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. + * 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(). + * 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); - - /** - * 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 - * whether @p add_ghost_data is @p true or @p - * false) 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(). + * 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(). * - * Must follow a call to the @p compress_start - * function. + * Must follow a call to the @p compress_start function. */ - void compress_finish (const bool add_ghost_data = true); + void compress_finish (::dealii::VectorOperation::values operation); + /** + * @deprecated: use compress(VectorOperation::values) instead. + */ + void compress_finish (const bool add_ghost_data = true) DEAL_II_DEPRECATED; /** - * 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. + * 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(). + * 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. + * 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. + * 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. + * 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. + * This method zeros the entries on ghost dofs, but does not touch + * locally owned DoFs. */ void zero_out_ghosts (); /** - * Return whether the vector contains only - * elements with value zero. This function - * is mainly for internal consistency - * checks and should seldom be used when - * not in debug mode since it uses quite - * some time. + * Return whether the vector contains only elements with value + * zero. This function is mainly for internal consistency checks and + * should seldom be used when not in debug mode since it uses quite some + * time. */ bool all_zero () const; /** - * Return @p true if the vector has no - * negative entries, i.e. all entries are - * zero or positive. This function is - * used, for example, to check whether - * refinement indicators are really all - * positive (or zero). + * 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. + * 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; @@ -438,90 +384,75 @@ namespace parallel Number operator * (const Vector &V) const; /** - * Computes the square of the l2 - * norm of the vector (i.e., the sum of the - * squares of all entries among all - * processors). + * Computes the square of the l2 norm of the vector (i.e., + * the sum of the squares of all entries among all processors). */ real_type norm_sqr () const; /** - * Computes the mean value of all the entries - * in the vector. + * Computes the mean value of all the entries in the vector. */ Number mean_value () const; /** - * Returns the l1 norm of the - * vector (i.e., the sum of the absolute - * values of all entries among all - * processors). + * Returns the l1 norm of the vector (i.e., the sum of the + * absolute values of all entries among all processors). */ real_type l1_norm () const; /** - * Returns the l2 norm of the - * vector (i.e., square root of the sum of the - * square of all entries among all - * processors). + * Returns the l2 norm of the vector (i.e., square root of + * the sum of the square of all entries among all processors). */ real_type l2_norm () const; /** - * Returns the lp norm with real @p - * p of the vector (i.e., the pth root of sum - * of the pth power of all entries among all + * Returns the lp norm with real @p p of the vector (i.e., + * the pth root of sum of the pth power of all entries among all * processors). */ real_type lp_norm (const real_type p) const; /** - * Returns the maximum norm of the vector - * (i.e., maximum absolute value among all - * entries among all processors). + * Returns the maximum norm of the vector (i.e., maximum absolute value + * among all entries among all processors). */ real_type linfty_norm () const; /** - * Returns the global size of the vector, - * equal to the sum of the number of locally - * owned indices among all the processors. + * Returns the global size of the vector, equal to the sum of the number + * of locally owned indices among all the processors. */ types::global_dof_index size () const; /** - * Returns the local size of the vector, i.e., - * the number of indices owned locally. + * Returns the local size of the vector, i.e., the number of indices + * owned locally. */ unsigned int local_size() const; /** - * Returns the half-open interval that - * specifies the locally owned range of the - * vector. Note that local_size() == - * local_range().second - + * Returns the half-open interval that specifies the locally owned range + * of the vector. Note that local_size() == local_range().second - * local_range().first. */ std::pair local_range () const; /** - * Returns true if the given global index is - * in the local range of this processor. + * Returns true if the given global index is in the local range of this + * processor. */ bool in_local_range (const types::global_dof_index global_index) const; /** - * Returns the number of ghost elements - * present on the vector. + * Returns the number of ghost elements present on the vector. */ unsigned int n_ghost_entries () const; /** - * 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. + * 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. */ bool is_ghost_entry (const types::global_dof_index global_index) const; @@ -560,69 +491,62 @@ namespace parallel //@{ /** - * 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. + * Read access to the data in the position corresponding to @p + * global_index. The index must be either in the local range of the + * vector or be specified as a ghost index at construction. + * + * Performance: O(1) for locally owned elements that represent + * a contiguous range and O(log(nranges)) for ghost + * elements (quite fast, but slower than local_element()). */ Number operator () (const types::global_dof_index 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. + * Read and write access to the data in the position corresponding to @p + * global_index. The index must be either in the local range of the + * vector or be specified as a ghost index at construction. + * + * Performance: O(1) for locally owned elements that represent + * a contiguous range and O(log(nranges)) for ghost + * elements (quite fast, but slower than local_element()). */ Number &operator () (const types::global_dof_index 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. + * 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(). + * This function does the same thing as operator(). */ Number operator [] (const types::global_dof_index 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. + * 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(). + * This function does the same thing as operator(). */ Number &operator [] (const types::global_dof_index global_index); /** - * Read access to the data field specified by - * @p local_index. Locally owned indices can - * be accessed with indices - * [0,local_size), and ghost - * indices with indices - * [local_size,local_size+ - * n_ghost_entries]. + * Read access to the data field specified by @p local_index. Locally + * owned indices can be accessed with indices + * [0,local_size), and ghost indices with indices + * [local_size,local_size+ n_ghost_entries]. + * + * Performance: Direct array access (fast). */ Number local_element (const unsigned int local_index) const; /** - * Read and write access to the data field - * specified by @p local_index. Locally owned - * indices can be accessed with indices - * [0,local_size), and ghost - * indices with indices + * Read and write access to the data field specified by @p + * local_index. Locally owned indices can be accessed with indices + * [0,local_size), and ghost indices with indices * [local_size,local_size+n_ghosts]. + * + * Performance: Direct array access (fast). */ Number &local_element (const unsigned int local_index); //@} @@ -634,48 +558,35 @@ namespace parallel //@{ /** - * Add the given vector to the present - * one. + * Add the given vector to the present one. */ Vector &operator += (const Vector &V); /** - * Subtract the given vector from the - * present one. + * Subtract the given vector from the present one. */ Vector &operator -= (const Vector &V); /** - * A collective add operation: - * This funnction adds a whole - * set of values stored in @p - * values to the vector - * components specified by @p - * indices. + * A collective add operation: This funnction adds a whole set of values + * stored in @p values to the vector components specified by @p indices. */ template void add (const std::vector &indices, const std::vector &values); /** - * This is a second collective - * add operation. As a - * difference, this function - * takes a deal.II vector of - * values. + * This is a second collective add operation. As a difference, this + * function takes a deal.II vector of values. */ template void add (const std::vector &indices, const ::dealii::Vector &values); /** - * Take an address where - * n_elements are stored - * contiguously and add them into - * the vector. Handles all cases - * which are not covered by the - * other two add() - * functions above. + * Take an address where n_elements are stored contiguously and + * add them into the vector. Handles all cases which are not covered by + * the other two add() functions above. */ template void add (const unsigned int n_elements, @@ -683,42 +594,37 @@ namespace parallel const OtherNumber *values); /** - * Addition of @p s to all - * components. Note that @p s is a - * scalar and not a vector. + * Addition of @p s to all components. Note that @p s is a scalar and + * not a vector. */ void add (const Number s); /** - * Simple vector addition, equal to the - * operator +=. + * Simple vector addition, equal to the operator +=. */ void add (const Vector &V); /** - * Simple addition of a multiple of a - * vector, i.e. *this += a*V. + * Simple addition of a multiple of a vector, i.e. *this += + * a*V. */ void add (const Number a, const Vector &V); /** - * Multiple addition of scaled vectors, - * i.e. *this += a*V+b*W. + * Multiple addition of scaled vectors, i.e. *this += a*V+b*W. */ void add (const Number a, const Vector &V, const Number b, const Vector &W); /** - * Scaling and simple vector addition, - * i.e. - * *this = s*(*this)+V. + * Scaling and simple vector addition, i.e. *this = + * s*(*this)+V. */ void sadd (const Number s, const Vector &V); /** - * Scaling and simple addition, i.e. - * *this = s*(*this)+a*V. + * Scaling and simple addition, i.e. *this = s*(*this)+a*V. */ void sadd (const Number s, const Number a, @@ -734,8 +640,8 @@ namespace parallel const Vector &W); /** - * Scaling and multiple addition. - * *this = s*(*this)+a*V + b*W + c*X. + * Scaling and multiple addition. *this = s*(*this)+a*V + b*W + + * c*X. */ void sadd (const Number s, const Number a, @@ -746,51 +652,35 @@ namespace parallel const Vector &X); /** - * Scale each element of the - * vector by the given factor. + * Scale each element of the vector by the given factor. * - * @deprecated This function is deprecated - * and will be removed in a - * future version. Use - * operator *= and - * operator /= instead. + * @deprecated This function is deprecated and will be removed in a + * future version. Use operator *= and operator /= + * instead. */ void scale (const Number factor) DEAL_II_DEPRECATED; - /** - * Scale each element of the - * vector by a constant - * value. + * Scale each element of the vector by a constant value. */ Vector &operator *= (const Number factor); /** - * Scale each element of the - * vector by the inverse of the - * given value. + * Scale each element of the vector by the inverse of the given value. */ Vector &operator /= (const Number factor); /** - * Scale each element of this - * vector by the corresponding - * element in the argument. This - * function is mostly meant to - * simulate multiplication (and - * immediate re-assignment) by a - * diagonal scaling matrix. + * Scale each element of this vector by the corresponding element in the + * argument. This function is mostly meant to simulate multiplication + * (and immediate re-assignment) by a diagonal scaling matrix. */ void scale (const Vector &scaling_factors); /** - * Scale each element of this - * vector by the corresponding - * element in the argument. This - * function is mostly meant to - * simulate multiplication (and - * immediate re-assignment) by a - * diagonal scaling matrix. + * Scale each element of this vector by the corresponding element in the + * argument. This function is mostly meant to simulate multiplication + * (and immediate re-assignment) by a diagonal scaling matrix. */ template void scale (const Vector &scaling_factors); @@ -820,20 +710,14 @@ namespace parallel const Number c, const Vector &w); /** - * Compute the elementwise ratio of the - * two given vectors, that is let - * this[i] = a[i]/b[i]. This is - * useful for example if you want to - * compute the cellwise ratio of true to - * estimated error. + * Compute the elementwise ratio of the two given vectors, that is let + * this[i] = a[i]/b[i]. This is useful for example if you want + * to compute the cellwise ratio of true to estimated error. * - * This vector is appropriately - * scaled to hold the result. + * This vector is appropriately scaled to hold the result. * - * If any of the b[i] is - * zero, the result is - * undefined. No attempt is made - * to catch such situations. + * If any of the b[i] is zero, the result is undefined. No + * attempt is made to catch such situations. */ void ratio (const Vector &a, const Vector &b); @@ -845,29 +729,20 @@ namespace parallel */ //@{ /** - * Checks whether the given - * partitioner is compatible with the - * partitioner used for this - * vector. Two partitioners are - * compatible if the have the same - * local size and the same ghost - * indices. They do not necessarily - * need to be the same data - * field. 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. + * Checks whether the given partitioner is compatible with the + * partitioner used for this vector. Two partitioners are compatible if + * the have the same local size and the same ghost indices. They do not + * necessarily need to be the same data field. 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; /** - * Prints the vector to the output stream @p - * out. + * Prints the vector to the output stream @p out. */ void print (std::ostream &out, const unsigned int precision = 3, @@ -875,97 +750,80 @@ namespace parallel const bool across = true) const; /** - * Returns the memory consumption of this - * class in bytes. + * Returns the memory consumption of this class in bytes. */ std::size_t memory_consumption () const; //@} private: /** - * Shared pointer to store the parallel - * partitioning information. This information - * can be shared between several vectors that - * have the same partitioning. + * Shared pointer to store the parallel partitioning information. This + * information can be shared between several vectors that have the same + * partitioning. */ std_cxx1x::shared_ptr partitioner; /** - * The size that is currently allocated in the - * val array. + * The size that is currently allocated in the val array. */ unsigned int allocated_size; /** - * Pointer to the array of - * local elements of this vector. + * 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 + * 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; /** - * Provide this class with all functionality - * of ::dealii::Vector by creating a - * VectorView object. + * Provide this class with all functionality of ::dealii::Vector by + * creating a VectorView object. */ VectorView vector_view; #ifdef DEAL_II_WITH_MPI /** - * A vector that collects all requests from @p - * compress() operations. This class uses - * persistent MPI communicators, i.e., the - * communication channels are stored during - * successive calls to a given function. This - * reduces the overhead involved with setting - * up the MPI machinery, but it does not - * remove the need for a receive operation to - * be posted before the data can actually be - * sent. + * A vector that collects all requests from @p compress() + * operations. This class uses persistent MPI communicators, i.e., the + * communication channels are stored during successive calls to a given + * function. This reduces the overhead involved with setting up the MPI + * machinery, but it does not remove the need for a receive operation to + * be posted before the data can actually be sent. */ std::vector compress_requests; /** - * A vector that collects all requests from @p - * update_ghost_values() operations. This class uses - * persistent MPI communicators. + * A vector that collects all requests from @p update_ghost_values() + * operations. This class uses persistent MPI communicators. */ mutable std::vector update_ghost_values_requests; #endif /** - * A lock that makes sure that - * the @p compress and @p - * update_ghost_values functions - * give reasonable results also - * when used with several - * threads. + * 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. + * 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. + * A helper function that is used to resize the val array. */ void resize_val (const unsigned int new_allocated_size); /* - * Make all other vector types - * friends. + * Make all other vector types friends. */ template friend class Vector; }; @@ -1138,13 +996,12 @@ namespace parallel void Vector::compress (::dealii::VectorOperation::values operation) { - compress_start (); - if (operation == ::dealii::VectorOperation::insert) - compress_finish (false); - else - compress_finish (true); + compress_start (0, operation); + compress_finish(operation); } + + template inline void @@ -1154,6 +1011,20 @@ namespace parallel } + + template + inline + void + Vector::compress_finish (const bool add_value) + { + if (add_value) + compress_finish(VectorOperation::add); + else + compress_finish(VectorOperation::insert); + } + + + template inline void diff --git a/deal.II/include/deal.II/lac/parallel_vector.templates.h b/deal.II/include/deal.II/lac/parallel_vector.templates.h index 5622d55ba2..b29381dbcb 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.templates.h +++ b/deal.II/include/deal.II/lac/parallel_vector.templates.h @@ -1,7 +1,7 @@ //----------%----------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011, 2012 by the deal.II authors +// Copyright (C) 2011, 2012, 2013 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -123,12 +123,10 @@ namespace parallel { delete [] import_data; - // do not reallocate import_data directly, but - // only upon request. It is only used as - // temporary storage for compress() and - // update_ghost_values, and we might have vectors where - // we never call these methods and hence do - // not need to have the storage. + // do not reallocate import_data directly, but only upon request. It + // is only used as temporary storage for compress() and + // update_ghost_values, and we might have vectors where we never + // call these methods and hence do not need to have the storage. import_data = 0; } } @@ -171,12 +169,10 @@ namespace parallel { delete [] import_data; - // do not reallocate import_data directly, but - // only upon request. It is only used as - // temporary storage for compress() and - // update_ghost_values, and we might have vectors where - // we never call these methods and hence do - // not need to have the storage. + // do not reallocate import_data directly, but only upon request. It + // is only used as temporary storage for compress() and + // update_ghost_values, and we might have vectors where we never + // call these methods and hence do not need to have the storage. import_data = 0; } } @@ -200,9 +196,20 @@ namespace parallel template void - Vector::compress_start (const unsigned int counter) + Vector::compress_start (const unsigned int counter, + ::dealii::VectorOperation::values operation) { #ifdef DEAL_II_WITH_MPI + + // nothing to do for insert (only need to zero ghost entries in + // compress_finish(). in debug mode we still want to check consistency + // of the inserted data, therefore the communication is still + // initialized +#ifndef DEBUG + if (operation == VectorOperation::insert) + return; +#endif + const Utilities::MPI::Partitioner &part = *partitioner; // nothing to do when we neither have import @@ -272,18 +279,22 @@ namespace parallel ierr = MPI_Startall(compress_requests.size(),&compress_requests[0]); Assert (ierr == MPI_SUCCESS, ExcInternalError()); } -#else + +#else // ifdef DEAL_II_WITH_MPI (void)counter; + (void)operation; #endif + } template void - Vector::compress_finish (const bool add_ghost_data) + Vector::compress_finish (::dealii::VectorOperation::values operation) { #ifdef DEAL_II_WITH_MPI + const Utilities::MPI::Partitioner &part = *partitioner; // nothing to do when we neither have import @@ -297,11 +308,12 @@ namespace parallel const unsigned int n_import_targets = part.import_targets().size(); const unsigned int n_ghost_targets = part.ghost_targets().size(); - AssertDimension (n_ghost_targets+n_import_targets, - compress_requests.size()); + if (operation != dealii::VectorOperation::insert) + AssertDimension (n_ghost_targets+n_import_targets, + compress_requests.size()); // first wait for the receive to complete - if (n_import_targets > 0) + if (compress_requests.size() > 0 && n_import_targets > 0) { int ierr; ierr = MPI_Waitall (n_import_targets, &compress_requests[0], @@ -312,31 +324,38 @@ namespace parallel std::vector >::const_iterator my_imports = part.import_indices().begin(); - // If add_ghost_data is set, add the imported - // data to the local values. If not, set the - // vector entries. - if (add_ghost_data == true) + // If the operation is no insertion, add the imported data to the + // local values. For insert, nothing is done here (but in debug mode + // we assert that the specified value is either zero or matches with + // the ones already present + if (operation != dealii::VectorOperation::insert) for ( ; my_imports!=part.import_indices().end(); ++my_imports) for (unsigned int j=my_imports->first; jsecond; j++) local_element(j) += *read_position++; else for ( ; my_imports!=part.import_indices().end(); ++my_imports) - for (unsigned int j=my_imports->first; jsecond; j++) - local_element(j) = *read_position++; + for (unsigned int j=my_imports->first; jsecond; + j++, read_position++) + Assert(*read_position == 0. || + std::abs(local_element(j) - *read_position) < + local_element(j) * 100. * + std::numeric_limits::epsilon(), + ExcMessage("Inserted elements do not match.")); AssertDimension(read_position-import_data,part.n_import_indices()); } - if (n_ghost_targets > 0) + if (compress_requests.size() > 0 && n_ghost_targets > 0) { int ierr; ierr = MPI_Waitall (n_ghost_targets, &compress_requests[n_import_targets], MPI_STATUSES_IGNORE); Assert (ierr == MPI_SUCCESS, ExcInternalError()); - zero_out_ghosts (); } else AssertDimension (part.n_ghost_indices(), 0); + + zero_out_ghosts (); #else (void)add_ghost_data; #endif diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h index fa52dc196f..1cee5ba92b 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.h @@ -1851,7 +1851,7 @@ namespace internal inline void compress_finish (parallel::distributed::Vector &vec) { - vec.compress_finish(); + vec.compress_finish(::dealii::VectorOperation::add); } diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 3dec6628b3..e598adf066 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -67,30 +67,6 @@ DEAL_II_NAMESPACE_OPEN namespace VectorTools { - namespace - { - template - void perform_compress_insert (VEC &vec) - { - vec.compress(::dealii::VectorOperation::insert); - } - - template - void perform_compress_insert (::dealii::parallel::distributed::Vector &vec) - { - // should not do compress when inserting - // elements as the ghosts are fixed and some - // will not be set at all. Instead, zero the - // ghost data afterwards to get consistent - // data. Otherwise, at least with the layout - // in r27671 the tests/mpi/step-48/ncpu_10 - // will signal incorrect results because of - // incorrect interpolation. - vec.zero_out_ghosts(); - } - - } - template void interpolate (const Mapping &mapping, @@ -284,7 +260,7 @@ namespace VectorTools = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]]; } } - perform_compress_insert(vec); + vec.compress(VectorOperation::insert); } diff --git a/tests/mpi/p4est_save_01.cc b/tests/mpi/p4est_save_01.cc index b79f172c2c..6ad64226ee 100644 --- a/tests/mpi/p4est_save_01.cc +++ b/tests/mpi/p4est_save_01.cc @@ -48,7 +48,7 @@ void test() deallog << "hyper_cube" << std::endl; std::string filename = - (std::string("p4est_save_01/ncpu_") + Utilities::int_to_string(Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD)) + "/dat"); + (std::string("p4est_save_01/ncpu_") + Utilities::int_to_string(Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD)) + "/dat"); { parallel::distributed::Triangulation tr(MPI_COMM_WORLD); diff --git a/tests/mpi/p4est_save_02.cc b/tests/mpi/p4est_save_02.cc index 9bd7e49e04..e1ecc552dc 100644 --- a/tests/mpi/p4est_save_02.cc +++ b/tests/mpi/p4est_save_02.cc @@ -48,7 +48,7 @@ void test() deallog << "hyper_cube" << std::endl; std::string filename = - (std::string("p4est_save_02/ncpu_") + Utilities::int_to_string(Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD)) + "/dat"); + (std::string("p4est_save_02/ncpu_") + Utilities::int_to_string(Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD)) + "/dat"); { parallel::distributed::Triangulation tr(MPI_COMM_WORLD); diff --git a/tests/mpi/p4est_save_03.cc b/tests/mpi/p4est_save_03.cc index 015d3ef1db..1c09c01562 100644 --- a/tests/mpi/p4est_save_03.cc +++ b/tests/mpi/p4est_save_03.cc @@ -49,7 +49,7 @@ void test() deallog << "hyper_cube" << std::endl; std::string filename = - (std::string ("p4est_save_03/ncpu_") + Utilities::int_to_string (Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD)) + "/dat"); + (std::string ("p4est_save_03/ncpu_") + Utilities::int_to_string (Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD)) + "/dat"); { parallel::distributed::Triangulation tr (MPI_COMM_WORLD); diff --git a/tests/mpi/parallel_partitioner_01.cc b/tests/mpi/parallel_partitioner_01.cc index 898b54129d..8b0f06eca0 100644 --- a/tests/mpi/parallel_partitioner_01.cc +++ b/tests/mpi/parallel_partitioner_01.cc @@ -26,7 +26,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; diff --git a/tests/mpi/parallel_partitioner_02.cc b/tests/mpi/parallel_partitioner_02.cc index 1d016a4bbe..2dcdc26b82 100644 --- a/tests/mpi/parallel_partitioner_02.cc +++ b/tests/mpi/parallel_partitioner_02.cc @@ -25,7 +25,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; diff --git a/tests/mpi/parallel_partitioner_03.cc b/tests/mpi/parallel_partitioner_03.cc index d3b7c2df31..3522564859 100644 --- a/tests/mpi/parallel_partitioner_03.cc +++ b/tests/mpi/parallel_partitioner_03.cc @@ -26,7 +26,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; diff --git a/tests/mpi/parallel_partitioner_04.cc b/tests/mpi/parallel_partitioner_04.cc index defb95d8af..c88d460431 100644 --- a/tests/mpi/parallel_partitioner_04.cc +++ b/tests/mpi/parallel_partitioner_04.cc @@ -26,7 +26,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; diff --git a/tests/mpi/parallel_vector_01.cc b/tests/mpi/parallel_vector_01.cc index 2836b61528..5113a87dd4 100644 --- a/tests/mpi/parallel_vector_01.cc +++ b/tests/mpi/parallel_vector_01.cc @@ -25,7 +25,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -41,7 +41,7 @@ void test () v(myid*2)=myid*2.0; v(myid*2+1)=myid*2.0+1.0; - v.compress(); + v.compress(VectorOperation::insert); v*=2.0; if (myid == 0) diff --git a/tests/mpi/parallel_vector_02.cc b/tests/mpi/parallel_vector_02.cc index 3ad8747346..9252b882d8 100644 --- a/tests/mpi/parallel_vector_02.cc +++ b/tests/mpi/parallel_vector_02.cc @@ -25,7 +25,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -44,7 +44,7 @@ void test () v(myid*2)=myid*2.0; v(myid*2+1)=myid*2.0+1.0; - v.compress(); + v.compress(VectorOperation::add); v*=2.0; Assert(v(myid*2) == myid*4.0, ExcInternalError()); @@ -52,7 +52,7 @@ void test () // set ghost dof, compress v(1) = 7; - v.compress(); + v.compress(VectorOperation::add); if (myid == 0) { diff --git a/tests/mpi/parallel_vector_03.cc b/tests/mpi/parallel_vector_03.cc index c0afd2ecdf..9e42109dc2 100644 --- a/tests/mpi/parallel_vector_03.cc +++ b/tests/mpi/parallel_vector_03.cc @@ -25,7 +25,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -44,16 +44,15 @@ void test () v(myid*2)=myid*2.0; v(myid*2+1)=myid*2.0+1.0; - v.compress(); + v.compress(VectorOperation::insert); v*=2.0; Assert(v(myid*2) == myid*4.0, ExcInternalError()); Assert(v(myid*2+1) == myid*4.0+2.0, ExcInternalError()); - // set ghost dof on remote processors, - // compress (no addition) - if (myid > 0) - v(1) = 7; + // set ghost dof on all processors, compress + // (insert mode) + v(1) = 7; v.compress(VectorOperation::insert); if (myid == 0) diff --git a/tests/mpi/parallel_vector_04.cc b/tests/mpi/parallel_vector_04.cc index 266e9ee3a3..bb605d686c 100644 --- a/tests/mpi/parallel_vector_04.cc +++ b/tests/mpi/parallel_vector_04.cc @@ -26,7 +26,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -45,13 +45,13 @@ void test () v(myid*2)=myid*2.0; v(myid*2+1)=myid*2.0+1.0; - v.compress(); + v.compress(VectorOperation::insert); v*=2.0; Assert(v(myid*2) == myid*4.0, ExcInternalError()); Assert(v(myid*2+1) == myid*4.0+2.0, ExcInternalError()); - // set ghost dof on remote processors, no + // set ghost dof on remote process, no // compress called if (myid > 0) v(1) = 7; @@ -70,7 +70,7 @@ void test () // check that everything remains zero also // after compress - v.compress(); + v.compress(VectorOperation::add); Assert(v(myid*2) == 0., ExcInternalError()); Assert(v(myid*2+1) == 0., ExcInternalError()); diff --git a/tests/mpi/parallel_vector_05.cc b/tests/mpi/parallel_vector_05.cc index 6720c52a26..06f21f7cd7 100644 --- a/tests/mpi/parallel_vector_05.cc +++ b/tests/mpi/parallel_vector_05.cc @@ -25,7 +25,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -44,18 +44,18 @@ void test () v(myid*2)=myid*2.0; v(myid*2+1)=myid*2.0+1.0; - v.compress(); + v.compress(VectorOperation::insert); v*=2.0; Assert(v(myid*2) == myid*4.0, ExcInternalError()); Assert(v(myid*2+1) == myid*4.0+2.0, ExcInternalError()); - // set ghost dof on remote processors, + // set ghost dof to zero on remote processors, // compress if (myid > 0) v(1) = 0; - v.compress(); + v.compress(VectorOperation::add); // check that nothing has changed Assert(v(myid*2) == myid*4.0, ExcInternalError()); diff --git a/tests/mpi/parallel_vector_06.cc b/tests/mpi/parallel_vector_06.cc index d465624c1c..48dcd6aeed 100644 --- a/tests/mpi/parallel_vector_06.cc +++ b/tests/mpi/parallel_vector_06.cc @@ -26,7 +26,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -50,7 +50,7 @@ void test () v(myid*2)=myid*2.0; v(myid*2+1)=myid*2.0+1.0; } - v.compress(); + v.compress(VectorOperation::insert); v*=2.0; if (myid < 8) { diff --git a/tests/mpi/parallel_vector_07.cc b/tests/mpi/parallel_vector_07.cc index 5470bbc4fe..81867be21b 100644 --- a/tests/mpi/parallel_vector_07.cc +++ b/tests/mpi/parallel_vector_07.cc @@ -25,7 +25,7 @@ void test () { unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int numproc = Utilities::System::get_n_mpi_processes (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); if (myid==0) deallog << "numproc=" << numproc << std::endl; @@ -59,7 +59,6 @@ void test () for (unsigned i=0; i +#include +#include +#include +#include +#include + + +void test () +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + + if (myid==0) deallog << "numproc=" << numproc << std::endl; + + + // each processor owns 2 indices and all + // are ghosting element 1 (the second) + IndexSet local_owned(numproc*2); + local_owned.add_range(myid*2,myid*2+2); + IndexSet local_relevant(numproc*2); + local_relevant = local_owned; + local_relevant.add_range(1,2); + local_relevant.add_range(4,5); + + parallel::distributed::Vector v(local_owned, local_relevant, MPI_COMM_WORLD); + + // set local values and check them + v(myid*2)=myid*2.0; + v(myid*2+1)=myid*2.0+1.0; + + v.compress(VectorOperation::insert); + v*=2.0; + v.add(1.0); + + Assert(v(myid*2) == myid*4.0+1, ExcInternalError()); + Assert(v(myid*2+1) == myid*4.0+3.0, ExcInternalError()); + + // set ghost dof on all processors, compress + // (insert mode) + v(1) = 7; + v.compress(VectorOperation::insert); + + if (myid == 0) + { + deallog << myid*2 << ":" << v(myid*2) << std::endl; + deallog << myid*2+1 << ":" << v(myid*2+1) << std::endl; + } + + // import ghosts onto all procs + v.update_ghost_values(); + Assert (v(1) == 7.0, ExcInternalError()); + + // check l2 norm + const double l2_norm = v.l2_norm(); + if (myid == 0) + deallog << "L2 norm: " << l2_norm << std::endl; + + if (myid == 0) + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); + + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile(output_file_for_mpi("parallel_vector_13").c_str()); + deallog.attach(logfile); + deallog << std::setprecision(4); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test(); + } + else + test(); + +} diff --git a/tests/mpi/parallel_vector_13/ncpu_10/cmp/generic b/tests/mpi/parallel_vector_13/ncpu_10/cmp/generic new file mode 100644 index 0000000000..10f1691ab2 --- /dev/null +++ b/tests/mpi/parallel_vector_13/ncpu_10/cmp/generic @@ -0,0 +1,6 @@ + +DEAL:0::numproc=10 +DEAL:0::0:1.000 +DEAL:0::1:7.000 +DEAL:0::L2 norm: 103.4 +DEAL:0::OK diff --git a/tests/mpi/parallel_vector_13/ncpu_4/cmp/generic b/tests/mpi/parallel_vector_13/ncpu_4/cmp/generic new file mode 100644 index 0000000000..3054471d30 --- /dev/null +++ b/tests/mpi/parallel_vector_13/ncpu_4/cmp/generic @@ -0,0 +1,6 @@ + +DEAL:0::numproc=4 +DEAL:0::0:1.000 +DEAL:0::1:7.000 +DEAL:0::L2 norm: 26.83 +DEAL:0::OK diff --git a/tests/mpi/parallel_vector_interpolate.cc b/tests/mpi/parallel_vector_interpolate.cc new file mode 100644 index 0000000000..53b333f8e2 --- /dev/null +++ b/tests/mpi/parallel_vector_interpolate.cc @@ -0,0 +1,95 @@ +//-------------------------- parallel_vector_interpolate.cc ----------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2011, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//-------------------------- parallel_vector_interpolate.cc ----------------------- + +// check FETools::interpolate on parallel vector + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +void test () +{ + const unsigned int dim = 2; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(3); + + FE_Q fe1(1), fe2(2); + DoFHandler dof1(tria), dof2(tria); + dof1.distribute_dofs(fe1); + dof2.distribute_dofs(fe2); + + IndexSet locally_relevant_dofs1, locally_relevant_dofs2; + DoFTools::extract_locally_relevant_dofs (dof1, + locally_relevant_dofs1); + DoFTools::extract_locally_relevant_dofs (dof2, + locally_relevant_dofs2); + + parallel::distributed::Vector v1(dof1.locally_owned_dofs(), + locally_relevant_dofs1, + MPI_COMM_WORLD), + v2(dof2.locally_owned_dofs(), locally_relevant_dofs2, MPI_COMM_WORLD); + + // set first vector to 1 + VectorTools::interpolate(dof1, ConstantFunction(1.), v1); + for (unsigned int i=0; i1e-15) inv_mass_matrix.local_element(k) = 1./inv_mass_matrix.local_element(k);