From: Martin Kronbichler Date: Sun, 29 Dec 2013 19:52:11 +0000 (+0000) Subject: Fix Trilinos vectors X-Git-Tag: v8.2.0-rc1~1103 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3185724e0d895ff5ad925def23a5ac70ba8d2edb;p=dealii.git Fix Trilinos vectors git-svn-id: https://svn.dealii.org/trunk@32126 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index e398c01f87..8b705acd64 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -1,4 +1,4 @@ -// --------------------------------------------------------------------- + // --------------------------------------------------------------------- // $Id$ // // Copyright (C) 2013 by the deal.II authors @@ -106,7 +106,9 @@ inconvenience this causes. written into as an optional argument. This allows for pre-allocating all possible entries right away, which makes writing into the matrix from several threads possible (otherwise, only one processor at a time can write - off-processor data). + off-processor data). Similarly, TrilinosWrappers::MPI::Vector objects can + be initialized with hints to ghost elements for a writable vector that can + be added into from multiple threads.
(Martin Kronbichler, 2013/12/23) diff --git a/deal.II/include/deal.II/lac/block_sparsity_pattern.h b/deal.II/include/deal.II/lac/block_sparsity_pattern.h index e22741f97c..8082826f24 100644 --- a/deal.II/include/deal.II/lac/block_sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/block_sparsity_pattern.h @@ -1060,6 +1060,22 @@ namespace TrilinosWrappers BlockSparsityPattern (const std::vector ¶llel_partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD); + /** + * Initialize the pattern with two arrays of index sets that specify rows + * and columns of the matrix, where the size() of the IndexSets specifies + * the size of the blocks and the values in each IndexSet denotes the rows + * that are going to be saved in each block. The additional index set + * writable_rows is used to set all rows that we allow to write + * locally. This constructor is used to create matrices that allow several + * threads to write simultaneously into the matrix (to different rows, of + * course), see the method TrilinosWrappers::SparsityPattern::reinit + * method with three index set arguments for more details. + */ + BlockSparsityPattern (const std::vector &row_parallel_partitioning, + const std::vector &column_parallel_partitioning, + const std::vector &writeable_rows, + const MPI_Comm &communicator = MPI_COMM_WORLD); + /** * Resize the matrix to a tensor product of matrices with dimensions * defined by the arguments. diff --git a/deal.II/include/deal.II/lac/trilinos_block_vector.h b/deal.II/include/deal.II/lac/trilinos_block_vector.h index 9006b134c2..276f778902 100644 --- a/deal.II/include/deal.II/lac/trilinos_block_vector.h +++ b/deal.II/include/deal.II/lac/trilinos_block_vector.h @@ -71,20 +71,17 @@ namespace TrilinosWrappers { public: /** - * Typedef the base class for simpler - * access to its own typedefs. + * Typedef the base class for simpler access to its own typedefs. */ typedef BlockVectorBase BaseClass; /** - * Typedef the type of the underlying - * vector. + * Typedef the type of the underlying vector. */ typedef BaseClass::BlockType BlockType; /** - * Import the typedefs from the base - * class. + * Import the typedefs from the base class. */ typedef BaseClass::value_type value_type; typedef BaseClass::pointer pointer; @@ -96,90 +93,59 @@ namespace TrilinosWrappers typedef BaseClass::const_iterator const_iterator; /** - * Default constructor. Generate an - * empty vector without any blocks. + * Default constructor. Generate an empty vector without any blocks. */ BlockVector (); /** - * Constructor. Generate a block - * vector with as many blocks as - * there are entries in Input_Maps. - * For this non-distributed vector, - * the %parallel partitioning is not - * used, just the global size of the - * partitioner. + * Constructor. Generate a block vector with as many blocks as there are + * entries in Input_Maps. For this non-distributed vector, the %parallel + * partitioning is not used, just the global size of the partitioner. */ explicit BlockVector (const std::vector &partitioner); /** - * Constructor. Generate a block - * vector with as many blocks as - * there are entries in Input_Maps. - * For this non-distributed vector, - * the %parallel partitioning is not - * used, just the global size of the - * partitioner. + * Constructor. Generate a block vector with as many blocks as there are + * entries in Input_Maps. For this non-distributed vector, the %parallel + * partitioning is not used, just the global size of the partitioner. */ explicit BlockVector (const std::vector &partitioner, const MPI_Comm &communicator = MPI_COMM_WORLD); /** - * Copy-Constructor. Set all the - * properties of the non-%parallel - * vector to those of the given - * %parallel vector and import the - * elements. + * Copy-Constructor. Set all the properties of the non-%parallel vector to + * those of the given %parallel vector and import the elements. */ BlockVector (const MPI::BlockVector &V); /** - * Copy-Constructor. Set all the - * properties of the vector to those - * of the given input vector and copy - * the elements. + * Copy-Constructor. Set all the properties of the vector to those of the + * given input vector and copy the elements. */ BlockVector (const BlockVector &V); /** - * Creates a block vector - * consisting of - * num_blocks - * components, but there is no - * content in the individual - * components and the user has to - * fill appropriate data using a - * reinit of the blocks. + * Creates a block vector consisting of num_blocks components, + * but there is no content in the individual components and the user has + * to fill appropriate data using a reinit of the blocks. */ explicit BlockVector (const size_type num_blocks); /** - * Constructor. Set the number of - * blocks to n.size() and - * initialize each block with - * n[i] zero elements. + * Constructor. Set the number of blocks to n.size() and + * initialize each block with n[i] zero elements. * * References BlockVector.reinit(). */ explicit BlockVector (const std::vector &N); /** - * Constructor. Set the number of - * blocks to - * n.size(). Initialize the - * vector with the elements - * pointed to by the range of - * iterators given as second and - * third argument. Apart from the - * first argument, this - * constructor is in complete - * analogy to the respective - * constructor of the - * std::vector class, but the - * first argument is needed in - * order to know how to subdivide - * the block vector into - * different blocks. + * Constructor. Set the number of blocks to n.size(). Initialize + * the vector with the elements pointed to by the range of iterators given + * as second and third argument. Apart from the first argument, this + * constructor is in complete analogy to the respective constructor of the + * std::vector class, but the first argument is needed in order + * to know how to subdivide the block vector into different blocks. */ template BlockVector (const std::vector &n, @@ -196,8 +162,7 @@ namespace TrilinosWrappers * * @deprecated * - * See @ref GlossCompress "Compressing - * distributed objects" for more + * See @ref GlossCompress "Compressing distributed objects" for more * information. */ void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED; @@ -208,91 +173,60 @@ namespace TrilinosWrappers using BlockVectorBase::compress; /** - * Copy operator: fill all - * components of the vector that - * are locally stored with the - * given scalar value. + * Copy operator: fill all components of the vector that are locally + * stored with the given scalar value. */ BlockVector & operator = (const value_type s); /** - * Copy operator for a - * distributed Trilinos vector to - * a localized one. + * Copy operator for a distributed Trilinos vector to a localized one. */ BlockVector & operator = (const MPI::BlockVector &V); /** - * Copy operator for arguments of - * the same type. + * Copy operator for arguments of the same type. */ BlockVector & operator = (const BlockVector &V); /** - * Another copy function. This - * one takes a deal.II block - * vector and copies it into a - * TrilinosWrappers block - * vector. Note that the number - * of blocks has to be the same - * in the vector as in the input - * vector. Use the reinit() - * command for resizing the - * BlockVector or for changing - * the internal structure of the - * block components. + * Another copy function. This one takes a deal.II block vector and copies + * it into a TrilinosWrappers block vector. Note that the number of blocks + * has to be the same in the vector as in the input vector. Use the + * reinit() command for resizing the BlockVector or for changing the + * internal structure of the block components. * - * Since Trilinos only works on - * doubles, this function is - * limited to accept only one - * possible number type in the - * deal.II vector. + * Since Trilinos only works on doubles, this function is limited to + * accept only one possible number type in the deal.II vector. */ template BlockVector & operator = (const ::dealii::BlockVector &V); /** - * Reinitialize the BlockVector to - * contain as many blocks as there - * are Epetra_Maps given in the - * input argument, according to the - * global size of the individual - * components described in the - * maps. Note that the resulting - * vector will be stored completely - * on each process. The Epetra_Map - * is useful when data exchange - * with a distributed vector based - * on the same Epetra_map is - * intended. In that case, the same - * communicator is used for data - * exchange. + * Reinitialize the BlockVector to contain as many blocks as there are + * Epetra_Maps given in the input argument, according to the global size + * of the individual components described in the maps. Note that the + * resulting vector will be stored completely on each process. The + * Epetra_Map is useful when data exchange with a distributed vector based + * on the same Epetra_map is intended. In that case, the same communicator + * is used for data exchange. * - * If fast==false, the vector - * is filled with zeros. + * If fast==false, the vector is filled with zeros. */ void reinit (const std::vector &partitioning, const bool fast = false); /** - * Reinitialize the BlockVector to - * contain as many blocks as there - * are index sets given in the - * input argument, according to the - * global size of the individual - * components described in the - * index set, and using a given MPI - * communicator. The MPI - * communicator is useful when data - * exchange with a distributed - * vector based on the same - * initialization is intended. In - * that case, the same communicator - * is used for data exchange. + * Reinitialize the BlockVector to contain as many blocks as there are + * index sets given in the input argument, according to the global size of + * the individual components described in the index set, and using a given + * MPI communicator. The MPI communicator is useful when data exchange + * with a distributed vector based on the same initialization is + * intended. In that case, the same communicator is used for data + * exchange. * * If fast==false, the vector * is filled with zeros. @@ -302,13 +236,9 @@ namespace TrilinosWrappers const bool fast = false); /** - * Reinitialize the BlockVector to - * contain as many blocks as there - * are elements in the first - * argument, and with the respective - * sizes. Since no distribution map - * is given, all vectors are local - * vectors. + * Reinitialize the BlockVector to contain as many blocks as there are + * elements in the first argument, and with the respective sizes. Since no + * distribution map is given, all vectors are local vectors. * * If fast==false, the vector * is filled with zeros. @@ -324,76 +254,44 @@ namespace TrilinosWrappers void reinit (const MPI::BlockVector &V); /** - * Change the dimension to that - * of the vector V. The same - * applies as for the other - * reinit() function. + * Change the dimension to that of the vector V. The same applies + * as for the other reinit() function. * - * The elements of V are not - * copied, i.e. this function is - * the same as calling reinit - * (V.size(), fast). + * The elements of V are not copied, i.e. this function is the + * same as calling reinit (V.size(), fast). * - * Note that you must call this - * (or the other reinit() - * functions) function, rather - * than calling the reinit() - * functions of an individual - * block, to allow the block - * vector to update its caches of - * vector sizes. If you call - * reinit() on one of the - * blocks, then subsequent - * actions on this object may - * yield unpredictable results - * since they may be routed to - * the wrong block. + * Note that you must call this (or the other reinit() functions) + * function, rather than calling the reinit() functions of an individual + * block, to allow the block vector to update its caches of vector + * sizes. If you call reinit() on one of the blocks, then subsequent + * actions on this object may yield unpredictable results since they may + * be routed to the wrong block. */ void reinit (const BlockVector &V, const bool fast = false); /** - * Change the number of blocks to - * num_blocks. The individual - * blocks will get initialized with - * zero size, so it is assumed that - * the user resizes the - * individual blocks by herself - * in an appropriate way, and - * calls collect_sizes - * afterwards. + * Change the number of blocks to num_blocks. The individual + * blocks will get initialized with zero size, so it is assumed that the + * user resizes the individual blocks by herself in an appropriate way, + * and calls collect_sizes afterwards. */ void reinit (const size_type num_blocks); /** - * Swap the contents of this - * vector and the other vector - * v. One could do this - * operation with a temporary - * variable and copying over the - * data elements, but this - * function is significantly more - * efficient since it only swaps - * the pointers to the data of - * the two vectors and therefore - * does not need to allocate - * temporary storage and move - * data around. + * Swap the contents of this vector and the other vector v. One + * could do this operation with a temporary variable and copying over the + * data elements, but this function is significantly more efficient since + * it only swaps the pointers to the data of the two vectors and therefore + * does not need to allocate temporary storage and move data around. * - * Limitation: right now this - * function only works if both - * vectors have the same number - * of blocks. If needed, the - * numbers of blocks should be + * Limitation: right now this function only works if both vectors have the + * same number of blocks. If needed, the numbers of blocks should be * exchanged, too. * - * This function is analog to the - * the swap() function of all C++ - * standard containers. Also, - * there is a global function - * swap(u,v) that simply calls - * u.swap(v), again in analogy - * to standard functions. + * This function is analog to the the swap() function of all C++ standard + * containers. Also, there is a global function swap(u,v) that simply + * calls u.swap(v), again in analogy to standard functions. */ void swap (BlockVector &v); diff --git a/deal.II/include/deal.II/lac/trilinos_parallel_block_vector.h b/deal.II/include/deal.II/lac/trilinos_parallel_block_vector.h index 7ec7a0e217..de56aa813d 100644 --- a/deal.II/include/deal.II/lac/trilinos_parallel_block_vector.h +++ b/deal.II/include/deal.II/lac/trilinos_parallel_block_vector.h @@ -73,20 +73,17 @@ namespace TrilinosWrappers { public: /** - * Typedef the base class for simpler - * access to its own typedefs. + * Typedef the base class for simpler access to its own typedefs. */ typedef BlockVectorBase BaseClass; /** - * Typedef the type of the underlying - * vector. + * Typedef the type of the underlying vector. */ typedef BaseClass::BlockType BlockType; /** - * Import the typedefs from the base - * class. + * Import the typedefs from the base class. */ typedef BaseClass::value_type value_type; typedef BaseClass::pointer pointer; @@ -98,62 +95,47 @@ namespace TrilinosWrappers typedef BaseClass::const_iterator const_iterator; /** - * Default constructor. Generate an - * empty vector without any blocks. + * Default constructor. Generate an empty vector without any blocks. */ BlockVector (); /** - * Constructor. Generate a block - * vector with as many blocks as - * there are entries in @p - * partitioning. Each Epetra_Map - * contains the layout of the - * distribution of data among the MPI - * processes. + * Constructor. Generate a block vector with as many blocks as there are + * entries in @p partitioning. Each Epetra_Map contains the layout of + * the distribution of data among the MPI processes. */ explicit BlockVector (const std::vector ¶llel_partitioning); /** - * Constructor. Generate a block - * vector with as many blocks as - * there are entries in - * @p partitioning. Each IndexSet - * together with the MPI communicator - * contains the layout of the - * distribution of data among the MPI - * processes. + * Constructor. Generate a block vector with as many blocks as there are + * entries in @p partitioning. Each IndexSet together with the MPI + * communicator contains the layout of the distribution of data among + * the MPI processes. */ explicit BlockVector (const std::vector ¶llel_partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD); /** - * Creates a BlockVector with ghost elements. @p ghost_values - * may contain any elements in @p parallel_partitioning, they will - * be ignored. + * Creates a BlockVector with ghost elements. See the respective + * reinit() method for more details. @p ghost_values may contain any + * elements in @p parallel_partitioning, they will be ignored. */ BlockVector (const std::vector ¶llel_partitioning, const std::vector &ghost_values, - const MPI_Comm &communicator); + const MPI_Comm &communicator, + const bool vector_writable); /** - * Copy-Constructor. Set all the - * properties of the parallel vector - * to those of the given argument and - * copy the elements. + * Copy-Constructor. Set all the properties of the parallel vector to + * those of the given argument and copy the elements. */ BlockVector (const BlockVector &V); /** - * Creates a block vector - * consisting of - * num_blocks - * components, but there is no - * content in the individual - * components and the user has to - * fill appropriate data using a - * reinit of the blocks. + * Creates a block vector consisting of num_blocks components, + * but there is no content in the individual components and the user has + * to fill appropriate data using a reinit of the blocks. */ explicit BlockVector (const size_type num_blocks); @@ -163,168 +145,124 @@ namespace TrilinosWrappers ~BlockVector (); /** - * Copy operator: fill all - * components of the vector that - * are locally stored with the - * given scalar value. + * Copy operator: fill all components of the vector that are locally + * stored with the given scalar value. */ BlockVector & operator = (const value_type s); /** - * Copy operator for arguments of - * the same type. + * Copy operator for arguments of the same type. */ BlockVector & operator = (const BlockVector &V); /** - * Copy operator for arguments of - * the localized Trilinos vector - * type. + * Copy operator for arguments of the localized Trilinos vector type. */ BlockVector & operator = (const ::dealii::TrilinosWrappers::BlockVector &V); /** - * Another copy function. This - * one takes a deal.II block - * vector and copies it into a - * TrilinosWrappers block - * vector. Note that the number - * of blocks has to be the same - * in the vector as in the input - * vector. Use the reinit() - * command for resizing the - * BlockVector or for changing - * the internal structure of the - * block components. + * Another copy function. This one takes a deal.II block vector and + * copies it into a TrilinosWrappers block vector. Note that the number + * of blocks has to be the same in the vector as in the input + * vector. Use the reinit() command for resizing the BlockVector or for + * changing the internal structure of the block components. * - * Since Trilinos only works on - * doubles, this function is - * limited to accept only one - * possible number type in the - * deal.II vector. + * Since Trilinos only works on doubles, this function is limited to + * accept only one possible number type in the deal.II vector. */ template BlockVector & operator = (const ::dealii::BlockVector &V); /** - * Reinitialize the BlockVector to - * contain as many blocks as there - * are Epetra_Maps given in the input - * argument, according to the - * parallel distribution of the - * individual components described - * in the maps. + * Reinitialize the BlockVector to contain as many blocks as there are + * Epetra_Maps given in the input argument, according to the parallel + * distribution of the individual components described in the maps. * - * If fast==false, the vector - * is filled with zeros. + * If fast==false, the vector is filled with zeros. */ void reinit (const std::vector ¶llel_partitioning, const bool fast = false); /** - * Reinitialize the BlockVector to - * contain as many blocks as there - * are index sets given in the input - * argument, according to the - * parallel distribution of the - * individual components described - * in the maps. + * Reinitialize the BlockVector to contain as many blocks as there are + * index sets given in the input argument, according to the parallel + * distribution of the individual components described in the maps. * - * If fast==false, the vector - * is filled with zeros. + * If fast==false, the vector is filled with zeros. */ void reinit (const std::vector ¶llel_partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD, const bool fast = false); + /** - * like above, but with a second set of indices for - * ghost entries. + * Reinit functionality. This function destroys the old vector content + * and generates a new one based on the input partitioning. In addition + * to just specifying one index set as in all the other methods above, + * this method allows to supply an additional set of ghost + * entries. There are two different versions of a vector that can be + * created. If the flag @p vector_writable is set to @p false, the + * vector only allows read access to the joint set of @p + * parallel_partitioning and @p ghost_entries. The effect of the reinit + * method is then equivalent to calling the other reinit method with an + * index set containing both the locally owned entries and the ghost + * entries. + * + * If the flag @p vector_writable is set to true, this creates an + * alternative storage scheme for ghost elements that allows multiple + * threads to write into the vector (for the other reinit methods, only + * one thread is allowed to write into the ghost entries at a time). */ void reinit (const std::vector &partitioning, const std::vector &ghost_values, - const MPI_Comm &communicator = MPI_COMM_WORLD); + const MPI_Comm &communicator = MPI_COMM_WORLD, + const bool vector_writable = false); /** - * Change the dimension to that - * of the vector V. The same - * applies as for the other - * reinit() function. + * Change the dimension to that of the vector V. The same + * applies as for the other reinit() function. * - * The elements of V are not - * copied, i.e. this function is - * the same as calling reinit - * (V.size(), fast). + * The elements of V are not copied, i.e. this function is the + * same as calling reinit (V.size(), fast). * - * Note that you must call this - * (or the other reinit() - * functions) function, rather - * than calling the reinit() - * functions of an individual - * block, to allow the block - * vector to update its caches of - * vector sizes. If you call - * reinit() on one of the - * blocks, then subsequent - * actions on this object may - * yield unpredictable results - * since they may be routed to - * the wrong block. + * Note that you must call this (or the other reinit() functions) + * function, rather than calling the reinit() functions of an individual + * block, to allow the block vector to update its caches of vector + * sizes. If you call reinit() on one of the blocks, then subsequent + * actions on this object may yield unpredictable results since they may + * be routed to the wrong block. */ void reinit (const BlockVector &V, const bool fast = false); /** - * Change the number of blocks to - * num_blocks. The individual - * blocks will get initialized with - * zero size, so it is assumed that - * the user resizes the - * individual blocks by herself - * in an appropriate way, and - * calls collect_sizes - * afterwards. + * Change the number of blocks to num_blocks. The individual + * blocks will get initialized with zero size, so it is assumed that the + * user resizes the individual blocks by herself in an appropriate way, + * and calls collect_sizes afterwards. */ void reinit (const size_type num_blocks); /** - * This reinit function is meant to - * be used for parallel - * calculations where some - * non-local data has to be - * used. The typical situation - * where one needs this function is - * the call of the - * FEValues::get_function_values - * function (or of some - * derivatives) in parallel. Since - * it is usually faster to retrieve - * the data in advance, this - * function can be called before - * the assembly forks out to the - * different processors. What this - * function does is the following: - * It takes the information in the - * columns of the given matrix and - * looks which data couples between - * the different processors. That - * data is then queried from the - * input vector. Note that you - * should not write to the - * resulting vector any more, since - * the some data can be stored - * several times on different - * processors, leading to - * unpredictable results. In - * particular, such a vector cannot - * be used for matrix-vector - * products as for example done - * during the solution of linear - * systems. + * This reinit function is meant to be used for parallel calculations + * where some non-local data has to be used. The typical situation where + * one needs this function is the call of the + * FEValues::get_function_values function (or of some derivatives) + * in parallel. Since it is usually faster to retrieve the data in + * advance, this function can be called before the assembly forks out to + * the different processors. What this function does is the following: + * It takes the information in the columns of the given matrix and looks + * which data couples between the different processors. That data is + * then queried from the input vector. Note that you should not write to + * the resulting vector any more, since the some data can be stored + * several times on different processors, leading to unpredictable + * results. In particular, such a vector cannot be used for + * matrix-vector products as for example done during the solution of + * linear systems. */ void import_nonlocal_data_for_fe (const TrilinosWrappers::BlockSparseMatrix &m, const BlockVector &v); @@ -335,8 +273,7 @@ namespace TrilinosWrappers * * @deprecated * - * See @ref GlossCompress "Compressing - * distributed objects" for more + * See @ref GlossCompress "Compressing distributed objects" for more * information. */ void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED; @@ -348,17 +285,10 @@ namespace TrilinosWrappers /** - * Returns the state of the - * vector, i.e., whether - * compress() needs to be - * called after an operation - * requiring data - * exchange. Does only return - * non-true values when used in - * debug mode, since - * it is quite expensive to - * keep track of all operations - * that lead to the need for + * Returns the state of the vector, i.e., whether compress() needs to be + * called after an operation requiring data exchange. Does only return + * non-true values when used in debug mode, since it is quite + * expensive to keep track of all operations that lead to the need for * compress(). */ bool is_compressed () const; @@ -369,34 +299,21 @@ namespace TrilinosWrappers bool has_ghost_elements() const; /** - * Swap the contents of this - * vector and the other vector - * v. One could do this - * operation with a temporary - * variable and copying over the - * data elements, but this - * function is significantly more - * efficient since it only swaps - * the pointers to the data of - * the two vectors and therefore - * does not need to allocate - * temporary storage and move - * data around. + * Swap the contents of this vector and the other vector v. One + * could do this operation with a temporary variable and copying over + * the data elements, but this function is significantly more efficient + * since it only swaps the pointers to the data of the two vectors and + * therefore does not need to allocate temporary storage and move data + * around. * - * Limitation: right now this - * function only works if both - * vectors have the same number - * of blocks. If needed, the - * numbers of blocks should be + * Limitation: right now this function only works if both vectors have + * the same number of blocks. If needed, the numbers of blocks should be * exchanged, too. * - * This function is analog to the - * the swap() function of all C++ - * standard containers. Also, - * there is a global function - * swap(u,v) that simply calls - * u.swap(v), again in analogy - * to standard functions. + * This function is analog to the the swap() function of all C++ + * standard containers. Also, there is a global function swap(u,v) that + * simply calls u.swap(v), again in analogy to standard + * functions. */ void swap (BlockVector &v); @@ -446,15 +363,19 @@ namespace TrilinosWrappers } + inline BlockVector::BlockVector (const std::vector ¶llel_partitioning, const std::vector &ghost_values, - const MPI_Comm &communicator) + const MPI_Comm &communicator, + const bool vector_writable) { - reinit(parallel_partitioning, ghost_values, communicator); + reinit(parallel_partitioning, ghost_values, communicator, + vector_writable); } + inline BlockVector::BlockVector (const size_type num_blocks) { @@ -533,6 +454,7 @@ namespace TrilinosWrappers } + inline bool BlockVector::has_ghost_elements() const @@ -545,6 +467,8 @@ namespace TrilinosWrappers return ghosted; } + + inline void BlockVector::swap (BlockVector &v) diff --git a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h index 8babcea003..2a21d5daa9 100644 --- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h @@ -452,23 +452,23 @@ namespace TrilinosWrappers * * When writing into Trilinos matrices from several threads in shared * memory, several things must be kept in mind as there is no built-in locks - * in this class to prevent data races. Therefore, simultaneous access to - * the same matrix row at the same time leads to data races in general and - * must be explicitly avoided by the user. However, it is possible to access - * different rows of the matrix from several threads simultaneously - * under the following two conditions: + * in this class to prevent data races. Simultaneous access to the same + * matrix row at the same time can lead to data races and must be explicitly + * avoided by the user. However, it is possible to access different + * rows of the matrix from several threads simultaneously under the + * following two conditions: *
    *
  • The matrix uses only one MPI process. *
  • The matrix has been initialized from a * TrilinosWrappers::SparsityPattern object that in turn has been * initialized with the reinit function specifying three index sets, one * for the rows, one for the columns and for the larger set of @p - * writeable_rows. Note that all other reinit methods and constructors of - * TrilinosWrappers::SparsityPattern will result in a matrix that needs to - * allocate off-processor entries on demand, which breaks - * thread-safety. Of course, using the respective reinit method for the - * block Trilinos sparsity pattern and block matrix also results in - * thread-safety. + * writeable_rows, and the operation is an addition. Note that all other + * reinit methods and constructors of TrilinosWrappers::SparsityPattern + * will result in a matrix that needs to allocate off-processor entries on + * demand, which breaks thread-safety. Of course, using the respective + * reinit method for the block Trilinos sparsity pattern and block matrix + * also results in thread-safety. *
* * @ingroup TrilinosWrappers @@ -2499,25 +2499,17 @@ namespace TrilinosWrappers if (ierr > 0) ierr = 0; } - else if (nonlocal_matrix.get() != 0) - { - // this is the case when we have explicitly set the off-processor - // rows and want to create a separate matrix object for them (to - // retain thread-safety) - Assert (nonlocal_matrix->RowMap().LID(static_cast(row)) != -1, - ExcMessage("Attempted to write into off-processor matrix row " - "that has not be specified as being writable upon " - "initialization")); - ierr = nonlocal_matrix->ReplaceGlobalValues(row, n_columns, - col_value_ptr, - col_index_ptr); - } else ierr = matrix->ReplaceGlobalValues (1, (TrilinosWrappers::types::int_type *)&row, n_columns, col_index_ptr, &col_value_ptr, Epetra_FECrsMatrix::ROW_MAJOR); + // use the FECrsMatrix facilities for set even in the case when we + // have explicitly set the off-processor rows because that only works + // properly when adding elements, not when setting them (since we want + // to only touch elements that have been set explicitly, and there is + // no way on the receiving processor to identify them otherwise) } Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0])); @@ -2762,8 +2754,7 @@ namespace TrilinosWrappers - // inline "simple" functions that are - // called frequently and do only involve + // inline "simple" functions that are called frequently and do only involve // a call to some Trilinos function. inline SparseMatrix::size_type diff --git a/deal.II/include/deal.II/lac/trilinos_vector.h b/deal.II/include/deal.II/lac/trilinos_vector.h index 1c4762e18c..9e7896532b 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector.h +++ b/deal.II/include/deal.II/lac/trilinos_vector.h @@ -186,6 +186,18 @@ namespace TrilinosWrappers * operations at the same time, for example by placing zero additions if * necessary. * + *

Thread safety of Trilinos vectors

+ * + * When writing into Trilinos vectors from several threads in shared + * memory, several things must be kept in mind as there is no built-in + * locks in this class to prevent data races. Simultaneous access to the + * same vector entry at the same time results in data races and must be + * explicitly avoided by the user. However, it is possible to access + * different entries of the vector from several threads + * simultaneously when only one MPI process is present or the vector has + * been constructed with an additional index set for ghost entries in + * write mode. + * * @ingroup TrilinosWrappers * @ingroup Vectors * @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009 @@ -199,15 +211,14 @@ namespace TrilinosWrappers typedef dealii::types::global_dof_index size_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. + * 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. + * For the current class, the variable equals true, since it does + * support parallel data storage. */ static const bool supports_distributed_data = true; @@ -216,20 +227,14 @@ namespace TrilinosWrappers */ //@{ /** - * Default constructor that - * generates an empty (zero size) - * vector. The function - * reinit() will have to - * give the vector the correct - * size and distribution among - * processes in case of an MPI - * run. + * Default constructor that generates an empty (zero size) vector. The + * function reinit() will have to give the vector the correct + * size and distribution among processes in case of an MPI run. */ Vector (); /** - * Copy constructor using the - * given vector. + * Copy constructor using the given vector. */ Vector (const Vector &V); @@ -239,156 +244,88 @@ namespace TrilinosWrappers ~Vector (); /** - * Reinit functionality. This - * function sets the calling vector - * to the dimension and the parallel - * distribution of the input vector, - * but does not copy the elements in - * v. If fast is - * not true, the elements in - * the vector are initialized with - * zero, otherwise the content will - * be left unchanged and the user has - * to set all elements. + * Reinit functionality. This function sets the calling vector to the + * dimension and the parallel distribution of the input vector, but does + * not copy the elements in v. If fast is not + * true, the elements in the vector are initialized with zero, + * otherwise the content will be left unchanged and the user has to set + * all elements. * - * This function has a third argument, - * allow_different_maps, - * that allows for an exchange of - * data between two equal-sized - * vectors (but being distributed - * differently among the - * processors). A trivial application - * of this function is to generate a - * replication of a whole vector on - * each machine, when the calling - * vector is built according to the - * localized vector class - * TrilinosWrappers::Vector, and - * v is a distributed - * vector. In this case, the variable - * fast needs to be set to - * false, since it does not - * make sense to exchange data - * between differently parallelized - * vectors without touching the - * elements. + * This function has a third argument, allow_different_maps, + * that allows for an exchange of data between two equal-sized vectors + * (but being distributed differently among the processors). A trivial + * application of this function is to generate a replication of a whole + * vector on each machine, when the calling vector is built according to + * the localized vector class TrilinosWrappers::Vector, and v + * is a distributed vector. In this case, the variable fast + * needs to be set to false, since it does not make sense to + * exchange data between differently parallelized vectors without + * touching the elements. */ void reinit (const VectorBase &v, const bool fast = false, const bool allow_different_maps = false); - void reinit (const BlockVector &v, - const bool import_data = false); - /** - * Reinit function. Creates a vector - * using the IndexSet local as our - * own unknowns, add optional ghost - * values ghost. + * Create vector by merging components from a block vector. */ - void reinit (const IndexSet &local, - const IndexSet &ghost, - const MPI_Comm &communicator = MPI_COMM_WORLD); + void reinit (const BlockVector &v, + const bool import_data = false); /** - * Set all components of the - * vector to the given number @p - * s. Simply pass this down to - * the base class, but we still - * need to declare this function - * to make the example given in - * the discussion about making - * the constructor explicit work. + * Set all components of the vector to the given number @p s. Simply + * pass this down to the base class, but we still need to declare this + * function to make the example given in the discussion about making the + * constructor explicit work. */ Vector &operator = (const TrilinosScalar s); /** - * Copy the given vector. Resize - * the present vector if - * necessary. In this case, also - * the Epetra_Map that designs - * the parallel partitioning is - * taken from the input vector. + * Copy the given vector. Resize the present vector if necessary. In + * this case, also the Epetra_Map that designs the parallel partitioning + * is taken from the input vector. */ Vector & operator = (const Vector &V); /** - * Copy operator from a given - * localized vector (present on - * all processes) in - * TrilinosWrappers format to the - * current distributed - * vector. This function assumes - * that the calling vector (left - * hand object) already is of the - * same size as the right hand - * side vector. Otherwise, an - * exception will be thrown. + * Copy operator from a given localized vector (present on all + * processes) in TrilinosWrappers format to the current distributed + * vector. This function assumes that the calling vector (left hand + * object) already is of the same size as the right hand side + * vector. Otherwise, an exception will be thrown. */ Vector & operator = (const ::dealii::TrilinosWrappers::Vector &V); /** - * Another copy function. This - * one takes a deal.II vector and - * copies it into a - * TrilinosWrapper vector. Note - * that since we do not provide - * any Epetra_map that tells - * about the partitioning of the - * vector among the MPI - * processes, the size of the - * TrilinosWrapper vector has to - * be the same as the size of the - * input vector. In order to - * change the map, use the - * reinit(const Epetra_Map - * &input_map) function. + * Another copy function. This one takes a deal.II vector and copies it + * into a TrilinosWrapper vector. Note that since we do not provide any + * Epetra_map that tells about the partitioning of the vector among the + * MPI processes, the size of the TrilinosWrapper vector has to be the + * same as the size of the input vector. In order to change the map, use + * the reinit(const Epetra_Map &input_map) function. */ template Vector & operator = (const ::dealii::Vector &v); /** - * This reinit function is - * meant to be used for - * parallel calculations where - * some non-local data has to - * be used. The typical - * situation where one needs - * this function is the call of - * the - * FEValues::get_function_values - * function (or of some - * derivatives) in - * parallel. Since it is - * usually faster to retrieve - * the data in advance, this - * function can be called - * before the assembly forks - * out to the different - * processors. What this - * function does is the - * following: It takes the - * information in the columns - * of the given matrix and - * looks which data couples - * between the different - * processors. That data is - * then queried from the input - * vector. Note that you should - * not write to the resulting - * vector any more, since the - * some data can be stored - * several times on different - * processors, leading to - * unpredictable results. In - * particular, such a vector - * cannot be used for - * matrix-vector products as - * for example done during the - * solution of linear systems. + * This reinit function is meant to be used for parallel calculations + * where some non-local data has to be used. The typical situation where + * one needs this function is the call of the + * FEValues::get_function_values function (or of some derivatives) + * in parallel. Since it is usually faster to retrieve the data in + * advance, this function can be called before the assembly forks out to + * the different processors. What this function does is the following: + * It takes the information in the columns of the given matrix and looks + * which data couples between the different processors. That data is + * then queried from the input vector. Note that you should not write to + * the resulting vector any more, since the some data can be stored + * several times on different processors, leading to unpredictable + * results. In particular, such a vector cannot be used for + * matrix-vector products as for example done during the solution of + * linear systems. */ void import_nonlocal_data_for_fe (const dealii::TrilinosWrappers::SparseMatrix &matrix, @@ -399,35 +336,24 @@ namespace TrilinosWrappers */ //@{ /** - * This constructor takes an - * Epetra_Map that already knows - * how to distribute the - * individual components among - * the MPI processors. Since it - * also includes information - * about the size of the vector, - * this is all we need to - * generate a parallel vector. + * This constructor takes an Epetra_Map that already knows how to + * distribute the individual components among the MPI processors. Since + * it also includes information about the size of the vector, this is + * all we need to generate a parallel vector. */ explicit Vector (const Epetra_Map ¶llel_partitioning); /** - * Copy constructor from the - * TrilinosWrappers vector - * class. Since a vector of this - * class does not necessarily - * need to be distributed among - * processes, the user needs to - * supply us with an Epetra_Map - * that sets the partitioning - * details. + * Copy constructor from the TrilinosWrappers vector class. Since a + * vector of this class does not necessarily need to be distributed + * among processes, the user needs to supply us with an Epetra_Map that + * sets the partitioning details. */ Vector (const Epetra_Map ¶llel_partitioning, const VectorBase &v); /** - * Reinitialize from a deal.II - * vector. The Epetra_Map specifies the + * Reinitialize from a deal.II vector. The Epetra_Map specifies the * %parallel partitioning. */ template @@ -435,20 +361,15 @@ namespace TrilinosWrappers const dealii::Vector &v); /** - * Reinit functionality. This - * function destroys the old - * vector content and generates a - * new one based on the input - * map. + * Reinit functionality. This function destroys the old vector content + * and generates a new one based on the input map. */ void reinit (const Epetra_Map ¶llel_partitioning, const bool fast = false); /** - * Copy-constructor from deal.II - * vectors. Sets the dimension to that - * of the given vector, and copies all - * elements. + * Copy-constructor from deal.II vectors. Sets the dimension to that of + * the given vector, and copies all elements. */ template Vector (const Epetra_Map ¶llel_partitioning, @@ -459,14 +380,10 @@ namespace TrilinosWrappers */ //@{ /** - * This constructor takes an IndexSet - * that defines how to distribute the - * individual components among the - * MPI processors. Since it also - * includes information about the - * size of the vector, this is all we - * need to generate a %parallel - * vector. + * This constructor takes an IndexSet that defines how to distribute the + * individual components among the MPI processors. Since it also + * includes information about the size of the vector, this is all we + * need to generate a %parallel vector. */ explicit Vector (const IndexSet ¶llel_partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD); @@ -479,25 +396,18 @@ namespace TrilinosWrappers const MPI_Comm &communicator = MPI_COMM_WORLD); /** - * Copy constructor from the - * TrilinosWrappers vector - * class. Since a vector of this - * class does not necessarily need to - * be distributed among processes, - * the user needs to supply us with - * an IndexSet and an MPI - * communicator that set the - * partitioning details. + * Copy constructor from the TrilinosWrappers vector class. Since a + * vector of this class does not necessarily need to be distributed + * among processes, the user needs to supply us with an IndexSet and an + * MPI communicator that set the partitioning details. */ Vector (const IndexSet ¶llel_partitioning, const VectorBase &v, const MPI_Comm &communicator = MPI_COMM_WORLD); /** - * Copy-constructor from deal.II - * vectors. Sets the dimension to - * that of the given vector, and - * copies all the elements. + * Copy-constructor from deal.II vectors. Sets the dimension to that of + * the given vector, and copies all the elements. */ template Vector (const IndexSet ¶llel_partitioning, @@ -505,17 +415,37 @@ namespace TrilinosWrappers const MPI_Comm &communicator = MPI_COMM_WORLD); /** - * Reinit functionality. This function - * destroys the old vector content and - * generates a new one based on the - * input partitioning. The flag - * fast determines whether the - * vector should be filled with zero - * (false) or left untouched (true). + * Reinit functionality. This function destroys the old vector content + * and generates a new one based on the input partitioning. The flag + * fast determines whether the vector should be filled with + * zero (false) or left untouched (true). */ void reinit (const IndexSet ¶llel_partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD, const bool fast = false); + + /** + * Reinit functionality. This function destroys the old vector content + * and generates a new one based on the input partitioning. In addition + * to just specifying one index set as in all the other methods above, + * this method allows to supply an additional set of ghost + * entries. There are two different versions of a vector that can be + * created. If the flag @p vector_writable is set to @p false, the + * vector only allows read access to the joint set of @p + * parallel_partitioning and @p ghost_entries. The effect of the reinit + * method is then equivalent to calling the other reinit method with an + * index set containing both the locally owned entries and the ghost + * entries. + * + * If the flag @p vector_writable is set to true, this creates an + * alternative storage scheme for ghost elements that allows multiple + * threads to write into the vector (for the other reinit methods, only + * one thread is allowed to write into the ghost entries at a time). + */ + void reinit (const IndexSet &locally_owned_entries, + const IndexSet &ghost_entries, + const MPI_Comm &communicator = MPI_COMM_WORLD, + const bool vector_writable = false); //@} }; @@ -639,124 +569,86 @@ namespace TrilinosWrappers typedef dealii::types::global_dof_index size_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. + * 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 - * false, since it does not support parallel data storage. - * If you do need parallel data storage, use - * TrilinosWrappers::MPI::Vector. + * For the current class, the variable equals false, since it does not + * support parallel data storage. If you do need parallel data storage, + * use TrilinosWrappers::MPI::Vector. */ static const bool supports_distributed_data = false; /** - * Default constructor that - * generates an empty (zero size) - * vector. The function - * reinit() will have to - * give the vector the correct + * Default constructor that generates an empty (zero size) vector. The + * function reinit() will have to give the vector the correct * size. */ Vector (); /** - * This constructor takes as - * input the number of elements - * in the vector. + * This constructor takes as input the number of elements in the vector. */ explicit Vector (const size_type n); /** - * This constructor takes as - * input the number of elements - * in the vector. If the map is - * not localized, i.e., if there - * are some elements that are not - * present on all processes, only - * the global size of the map - * will be taken and a localized - * map will be generated - * internally. + * This constructor takes as input the number of elements in the + * vector. If the map is not localized, i.e., if there are some elements + * that are not present on all processes, only the global size of the map + * will be taken and a localized map will be generated internally. */ explicit Vector (const Epetra_Map &partitioning); /** - * This constructor takes as input - * the number of elements in the - * vector. If the index set is not - * localized, i.e., if there are some - * elements that are not present on - * all processes, only the global - * size of the index set will be - * taken and a localized version will - * be generated internally. + * This constructor takes as input the number of elements in the + * vector. If the index set is not localized, i.e., if there are some + * elements that are not present on all processes, only the global size of + * the index set will be taken and a localized version will be generated + * internally. */ explicit Vector (const IndexSet &partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD); /** - * This constructor takes a - * (possibly parallel) Trilinos - * Vector and generates a - * localized version of the whole - * content on each processor. + * This constructor takes a (possibly parallel) Trilinos Vector and + * generates a localized version of the whole content on each processor. */ explicit Vector (const VectorBase &V); /** - * Copy-constructor from deal.II - * vectors. Sets the dimension to that - * of the given vector, and copies all - * elements. + * Copy-constructor from deal.II vectors. Sets the dimension to that of + * the given vector, and copies all elements. */ template explicit Vector (const dealii::Vector &v); /** - * Reinit function that resizes - * the vector to the size - * specified by n. + * Reinit function that resizes the vector to the size specified by + * n. */ void reinit (const size_type n, const bool fast = false); /** - * Initialization with an - * Epetra_Map. Similar to the call in - * the other class MPI::Vector, with - * the difference that now a copy on - * all processes is generated. This - * initialization function is - * appropriate when the data in the - * localized vector should be - * imported from a distributed vector - * that has been initialized with the - * same communicator. The variable - * fast determines whether - * the vector should be filled with - * zero or left untouched. + * Initialization with an Epetra_Map. Similar to the call in the other + * class MPI::Vector, with the difference that now a copy on all processes + * is generated. This initialization function is appropriate when the data + * in the localized vector should be imported from a distributed vector + * that has been initialized with the same communicator. The variable + * fast determines whether the vector should be filled with zero + * or left untouched. */ void reinit (const Epetra_Map &input_map, const bool fast = false); /** - * Initialization with an - * IndexSet. Similar to the call in the - * other class MPI::Vector, with the - * difference that now a copy on all - * processes is generated. This - * initialization function is - * appropriate in case the data in the - * localized vector should be imported - * from a distributed vector that has - * been initialized with the same - * communicator. The variable - * fast determines whether the - * vector should be filled with zero + * Initialization with an IndexSet. Similar to the call in the other class + * MPI::Vector, with the difference that now a copy on all processes is + * generated. This initialization function is appropriate in case the data + * in the localized vector should be imported from a distributed vector + * that has been initialized with the same communicator. The variable + * fast determines whether the vector should be filled with zero * (false) or left untouched (true). */ void reinit (const IndexSet &input_map, @@ -764,66 +656,52 @@ namespace TrilinosWrappers const bool fast = false); /** - * Reinit function. Takes the - * information of a Vector and copies - * everything to the calling vector, - * now also allowing different maps. + * Reinit function. Takes the information of a Vector and copies + * everything to the calling vector, now also allowing different maps. */ void reinit (const VectorBase &V, const bool fast = false, const bool allow_different_maps = false); /** - * Set all components of the - * vector to the given number @p - * s. Simply pass this down to - * the base class, but we still - * need to declare this function - * to make the example given in - * the discussion about making - * the constructor explicit work. + * Set all components of the vector to the given number @p s. Simply pass + * this down to the base class, but we still need to declare this function + * to make the example given in the discussion about making the + * constructor explicit work. */ Vector &operator = (const TrilinosScalar s); /** - * Sets the left hand argument to - * the (parallel) Trilinos - * Vector. Equivalent to the @p - * reinit function. + * Sets the left hand argument to the (parallel) Trilinos + * Vector. Equivalent to the @p reinit function. */ Vector & operator = (const MPI::Vector &V); /** - * Sets the left hand argument to - * the deal.II vector. + * Sets the left hand argument to the deal.II vector. */ template Vector & operator = (const ::dealii::Vector &V); /** - * Copy operator. Copies both the - * dimension and the content in - * the right hand argument. + * Copy operator. Copies both the dimension and the content in the right + * hand argument. */ Vector & operator = (const Vector &V); /** - * This function does nothing but is - * there for compatibility with the - * @p PETScWrappers::Vector class. + * This function does nothing but is there for compatibility with the @p + * PETScWrappers::Vector class. * - * For the PETSc vector wrapper class, - * this function updates the ghost - * values of the PETSc vector. This - * is necessary after any modification + * For the PETSc vector wrapper class, this function updates the ghost + * values of the PETSc vector. This is necessary after any modification * before reading ghost values. * - * However, for the implementation of - * this class, it is immaterial and thus - * an empty function. + * However, for the implementation of this class, it is immaterial and + * thus an empty function. */ void update_ghost_values () const; }; diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h index 66995f19cd..649513392a 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector_base.h +++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h @@ -75,29 +75,20 @@ namespace TrilinosWrappers typedef dealii::types::global_dof_index size_type; /** - * This class implements a - * wrapper for accessing the - * Trilinos vector in the same - * way as we access deal.II - * objects: it is initialized - * with a vector and an element - * within it, and has a - * conversion operator to - * extract the scalar value of - * this element. It also has a - * variety of assignment - * operator for writing to this - * one element. @ingroup - * TrilinosWrappers + * This class implements a wrapper for accessing the Trilinos vector in + * the same way as we access deal.II objects: it is initialized with a + * vector and an element within it, and has a conversion operator to + * extract the scalar value of this element. It also has a variety of + * assignment operator for writing to this one element. + * + * @ingroup TrilinosWrappers */ class VectorReference { private: /** - * Constructor. It is made - * private so as to only allow - * the actual vector class to - * create it. + * Constructor. It is made private so as to only allow the actual vector + * class to create it. */ VectorReference (VectorBase &vector, const size_type index); @@ -105,83 +96,58 @@ namespace TrilinosWrappers public: /** - * This looks like a copy - * operator, but does something - * different than usual. In - * particular, it does not copy - * the member variables of this - * reference. Rather, it - * handles the situation where - * we have two vectors @p v and - * @p w, and assign elements - * like in - * v(i)=w(i). Here, - * both left and right hand - * side of the assignment have - * data type VectorReference, - * but what we really mean is - * to assign the vector - * elements represented by the - * two references. This - * operator implements this - * operation. Note also that - * this allows us to make the - * assignment operator const. + * This looks like a copy operator, but does something different than + * usual. In particular, it does not copy the member variables of this + * reference. Rather, it handles the situation where we have two vectors + * @p v and @p w, and assign elements like in v(i)=w(i). Here, + * both left and right hand side of the assignment have data type + * VectorReference, but what we really mean is to assign the vector + * elements represented by the two references. This operator implements + * this operation. Note also that this allows us to make the assignment + * operator const. */ const VectorReference & operator = (const VectorReference &r) const; /** - * Same as above but for non-const - * reference objects. + * Same as above but for non-const reference objects. */ const VectorReference & operator = (const VectorReference &r); /** - * Set the referenced element of the - * vector to s. + * Set the referenced element of the vector to s. */ const VectorReference & operator = (const TrilinosScalar &s) const; /** - * Add s to the - * referenced element of the - * vector-> + * Add s to the referenced element of the vector-> */ const VectorReference & operator += (const TrilinosScalar &s) const; /** - * Subtract s from the - * referenced element of the - * vector-> + * Subtract s from the referenced element of the vector-> */ const VectorReference & operator -= (const TrilinosScalar &s) const; /** - * Multiply the referenced - * element of the vector by - * s. + * Multiply the referenced element of the vector by s. */ const VectorReference & operator *= (const TrilinosScalar &s) const; /** - * Divide the referenced - * element of the vector by - * s. + * Divide the referenced element of the vector by s. */ const VectorReference & operator /= (const TrilinosScalar &s) const; /** - * Convert the reference to an - * actual value, i.e. return - * the value of the referenced - * element of the vector. + * Convert the reference to an actual value, i.e. return the value of + * the referenced element of the vector. */ operator TrilinosScalar () const; @@ -209,21 +175,17 @@ namespace TrilinosWrappers private: /** - * Point to the vector we are - * referencing. + * Point to the vector we are referencing. */ VectorBase &vector; /** - * Index of the referenced element - * of the vector. + * Index of the referenced element of the vector. */ const size_type index; /** - * Make the vector class a - * friend, so that it can - * create objects of the + * Make the vector class a friend, so that it can create objects of the * present type. */ friend class ::dealii::TrilinosWrappers::VectorBase; @@ -270,11 +232,8 @@ namespace TrilinosWrappers { public: /** - * Declare some of the standard - * types used in all - * containers. These types - * parallel those in the - * C standard libraries + * Declare some of the standard types used in all containers. These types + * parallel those in the C standard libraries * vector<...> class. */ typedef TrilinosScalar value_type; @@ -291,22 +250,15 @@ namespace TrilinosWrappers //@{ /** - * Default constructor that - * generates an empty (zero size) - * vector. The function - * reinit() will have to - * give the vector the correct - * size and distribution among - * processes in case of an MPI - * run. + * Default constructor that generates an empty (zero size) vector. The + * function reinit() will have to give the vector the correct + * size and distribution among processes in case of an MPI run. */ VectorBase (); /** - * Copy constructor. Sets the - * dimension to that of the given - * vector, and copies all the - * elements. + * Copy constructor. Sets the dimension to that of the given vector, and + * copies all the elements. */ VectorBase (const VectorBase &v); @@ -316,47 +268,30 @@ namespace TrilinosWrappers virtual ~VectorBase (); /** - * Release all memory and return - * to a state just like after - * having called the default - * constructor. + * Release all memory and return to a state just like after having called + * the default constructor. */ void clear (); /** - * Reinit functionality, sets the - * dimension and possibly the - * parallel partitioning (Epetra_Map) - * of the calling vector to the - * settings of the input vector. + * Reinit functionality, sets the dimension and possibly the parallel + * partitioning (Epetra_Map) of the calling vector to the settings of the + * input vector. */ void reinit (const VectorBase &v, const bool fast = false); /** - * Compress the underlying - * representation of the Trilinos - * object, i.e. flush the buffers - * of the vector object if it has - * any. This function is - * necessary after writing into a - * vector element-by-element and - * before anything else can be - * done on it. + * Compress the underlying representation of the Trilinos object, + * i.e. flush the buffers of the vector object if it has any. This + * function is necessary after writing into a vector element-by-element + * and before anything else can be done on it. * - * The (defaulted) argument can - * be used to specify the - * compress mode - * (Add or - * Insert) in case - * the vector has not been - * written to since the last - * time this function was - * called. The argument is - * ignored if the vector has - * been added or written to - * since the last time - * compress() was called. + * The (defaulted) argument can be used to specify the compress mode + * (Add or Insert) in case the vector has not + * been written to since the last time this function was called. The + * argument is ignored if the vector has been added or written to since + * the last time compress() was called. * * See @ref GlossCompress "Compressing distributed objects" * for more information. @@ -375,136 +310,87 @@ namespace TrilinosWrappers void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED; /** - * Returns the state of the - * vector, i.e., whether - * compress() has already been - * called after an operation - * requiring data exchange. + * Returns the state of the vector, i.e., whether compress() has already + * been called after an operation requiring data exchange. */ bool is_compressed () const; /** - * Set all components of the - * vector to the given number @p - * s. Simply pass this down to - * the Trilinos Epetra object, - * but we still need to declare - * this function to make the - * example given in the - * discussion about making the - * constructor explicit work. + * Set all components of the vector to the given number @p s. Simply pass + * this down to the Trilinos Epetra object, but we still need to declare + * this function to make the example given in the discussion about making + * the constructor explicit work. * - * Since the semantics of - * assigning a scalar to a vector - * are not immediately clear, - * this operator should really - * only be used if you want to - * set the entire vector to - * zero. This allows the - * intuitive notation - * v=0. Assigning other - * values is deprecated and may - * be disallowed in the future. + * Since the semantics of assigning a scalar to a vector are not + * immediately clear, this operator should really only be used if you want + * to set the entire vector to zero. This allows the intuitive notation + * v=0. Assigning other values is deprecated and may be + * disallowed in the future. */ VectorBase & operator = (const TrilinosScalar s); /** - * Copy function. This function takes - * a VectorBase vector and copies all - * the elements. The target vector - * will have the same parallel - * distribution as the calling - * vector. + * Copy function. This function takes a VectorBase vector and copies all + * the elements. The target vector will have the same parallel + * distribution as the calling vector. */ VectorBase & operator = (const VectorBase &v); /** - * Another copy function. This - * one takes a deal.II vector and - * copies it into a - * TrilinosWrapper vector. Note - * that since we do not provide - * any Epetra_map that tells - * about the partitioning of the - * vector among the MPI - * processes, the size of the - * TrilinosWrapper vector has to - * be the same as the size of the - * input vector. In order to - * change the map, use the - * reinit(const Epetra_Map - * &input_map) function. + * Another copy function. This one takes a deal.II vector and copies it + * into a TrilinosWrapper vector. Note that since we do not provide any + * Epetra_map that tells about the partitioning of the vector among the + * MPI processes, the size of the TrilinosWrapper vector has to be the + * same as the size of the input vector. In order to change the map, use + * the reinit(const Epetra_Map &input_map) function. */ template VectorBase & operator = (const ::dealii::Vector &v); /** - * Test for equality. This - * function assumes that the - * present vector and the one to - * compare with have the same - * size already, since comparing - * vectors of different sizes - * makes not much sense anyway. + * Test for equality. This function assumes that the present vector and + * the one to compare with have the same size already, since comparing + * vectors of different sizes makes not much sense anyway. */ bool operator == (const VectorBase &v) const; /** - * Test for inequality. This - * function assumes that the - * present vector and the one to - * compare with have the same - * size already, since comparing - * vectors of different sizes - * makes not much sense anyway. + * Test for inequality. This function assumes that the present vector and + * the one to compare with have the same size already, since comparing + * vectors of different sizes makes not much sense anyway. */ bool operator != (const VectorBase &v) const; /** - * Return the global dimension of - * the vector. + * Return the global dimension of the vector. */ size_type size () const; /** - * Return the local dimension of - * the vector, i.e. the number of - * elements stored on the present - * MPI process. For sequential - * vectors, this number is the - * same as size(), but for - * parallel vectors it may be - * smaller. + * Return the local dimension of the vector, i.e. the number of elements + * stored on the present MPI process. For sequential vectors, this number + * is the same as size(), but for parallel vectors it may be smaller. * - * To figure out which elements - * exactly are stored locally, - * use local_range(). + * To figure out which elements exactly are stored locally, use + * local_range(). * - * If the vector contains ghost - * elements, they are included in - * this number. + * If the vector contains ghost elements, they are included in this + * number. */ size_type local_size () const; /** - * Return a pair of indices - * indicating which elements of - * this vector are stored - * locally. The first number is - * the index of the first element - * stored, the second the index - * of the one past the last one - * that is stored locally. If - * this is a sequential vector, - * then the result will be the - * pair (0,N), otherwise it will - * be a pair (i,i+n), where - * n=local_size() and i is the first - * element of the vector stored on this processor, corresponding - * to the half open interval $[i,i+n)$ + * Return a pair of indices indicating which elements of this vector are + * stored locally. The first number is the index of the first element + * stored, the second the index of the one past the last one that is + * stored locally. If this is a sequential vector, then the result will be + * the pair (0,N), otherwise it will be a pair + * (i,i+n), where n=local_size() and + * i is the first element of the vector stored on this + * processor, corresponding to the half open interval $[i,i+n)$ * * @note The description above is true most of the time, but * not always. In particular, Trilinos vectors need not store @@ -522,9 +408,8 @@ namespace TrilinosWrappers std::pair local_range () const; /** - * Return whether @p index is in - * the local range or not, see - * also local_range(). + * Return whether @p index is in the local range or not, see also + * local_range(). * * @note The same limitation for the applicability of this * function applies as listed in the documentation of local_range(). @@ -548,88 +433,65 @@ namespace TrilinosWrappers IndexSet locally_owned_elements () const; /** - * Return if the vector contains ghost - * elements. This answer is true if there - * are ghost elements on at least one - * process. + * Return if the vector contains ghost elements. This answer is true if + * there are ghost elements on at least one process. */ bool has_ghost_elements() const; /** - * Return the scalar (inner) - * product of two vectors. The - * vectors must have the same - * size. + * Return the scalar (inner) product of two vectors. The vectors must have + * the same size. */ TrilinosScalar operator * (const VectorBase &vec) const; /** - * Return square of the - * $l_2$-norm. + * Return square of the $l_2$-norm. */ real_type norm_sqr () const; /** - * Mean value of the elements of - * this vector. + * Mean value of the elements of this vector. */ TrilinosScalar mean_value () const; /** - * Compute the minimal value of - * the elements of this vector. + * Compute the minimal value of the elements of this vector. */ TrilinosScalar minimal_value () const; /** - * $l_1$-norm of the vector. The - * sum of the absolute values. + * $l_1$-norm of the vector. The sum of the absolute values. */ real_type l1_norm () const; /** - * $l_2$-norm of the vector. The - * square root of the sum of the - * squares of the elements. + * $l_2$-norm of the vector. The square root of the sum of the squares of + * the elements. */ real_type l2_norm () const; /** - * $l_p$-norm of the vector. The - * pth root of the sum of - * the pth powers of the - * absolute values of the - * elements. + * $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 TrilinosScalar p) const; /** - * Maximum absolute value of the - * elements. + * Maximum absolute value of the elements. */ real_type linfty_norm () const; /** - * Return whether the vector - * contains only elements with - * value zero. This function is - * mainly for internal - * consistency checks and should - * seldom be used when not in - * debug mode since it uses quite - * some time. + * 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). */ bool is_non_negative () const; //@} @@ -641,23 +503,20 @@ namespace TrilinosWrappers //@{ /** - * Provide access to a given - * element, both read and write. + * Provide access to a given element, both read and write. */ reference operator () (const size_type index); /** - * Provide read-only access to an - * element. This is equivalent to - * the el() command. + * Provide read-only access to an element. This is equivalent to the + * el() command. */ TrilinosScalar operator () (const size_type index) const; /** - * Provide access to a given - * element, both read and write. + * Provide access to a given element, both read and write. * * Exactly the same as operator(). */ @@ -665,9 +524,8 @@ namespace TrilinosWrappers operator [] (const size_type index); /** - * Provide read-only access to an - * element. This is equivalent to - * the el() command. + * Provide read-only access to an element. This is equivalent to the + * el() command. * * Exactly the same as operator(). */ @@ -675,21 +533,17 @@ namespace TrilinosWrappers operator [] (const size_type index) const; /** - * 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. + * 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. */ void extract_subvector_to (const std::vector &indices, std::vector &values) const; /** - * Just as the above, but with pointers. - * Useful in minimizing copying of data around. + * Just as the above, but with pointers. Useful in minimizing copying of + * data around. */ template void extract_subvector_to (ForwardIterator indices_begin, @@ -697,22 +551,20 @@ namespace TrilinosWrappers OutputIterator values_begin) const; /** - * Return the value of the vector - * entry i. Note that this - * function does only work - * properly when we request a - * data stored on the local - * processor. The function will - * throw an exception in case the - * elements sits on another - * process. + * Return the value of the vector entry i. Note that this function + * does only work properly when we request a data stored on the local + * processor. The function will throw an exception in case the elements + * sits on another process. */ TrilinosScalar el (const size_type index) const; /** - * Make the Vector class a bit like the vector<> class of - * the C++ standard library by returning iterators to the start and end - * of the locally owned elements of this vector. The ordering of local elements corresponds to the one given + * Make the Vector class a bit like the vector<> class of the C++ + * standard library by returning iterators to the start and end of the + * locally owned elements of this vector. The ordering of local elements + * corresponds to the one given by the global indices in case the vector + * is constructed from an IndexSet or other methods in deal.II (note that + * an Epetra_Map can contain elements in arbitrary orders, though). * * It holds that end() - begin() == local_size(). */ @@ -731,31 +583,23 @@ namespace TrilinosWrappers iterator end (); /** - * Return a constant iterator pointing to the element past the end of - * the array of the locally owned entries. + * Return a constant iterator pointing to the element past the end of the + * array of the locally owned entries. */ const_iterator end () const; /** - * A collective set operation: - * instead of setting individual - * elements of a vector, this - * function allows to set a whole - * set of elements at once. The - * indices of the elements to be - * set are stated in the first - * argument, the corresponding - * values in the second. + * A collective set operation: instead of setting individual elements of a + * vector, this function allows to set a whole set of elements at + * once. The indices of the elements to be set are stated in the first + * argument, the corresponding values in the second. */ void set (const std::vector &indices, const std::vector &values); /** - * This is a second collective - * set operation. As a - * difference, this function - * takes a deal.II vector of - * values. + * This is a second collective set operation. As a difference, this + * function takes a deal.II vector of values. */ void set (const std::vector &indices, const ::dealii::Vector &values); @@ -768,109 +612,80 @@ namespace TrilinosWrappers //@{ /** - * This collective set operation - * is of lower level and can - * handle anything else — - * the only thing you have to - * provide is an address where - * all the indices are stored and - * the number of elements to be - * set. + * This collective set operation is of lower level and can handle anything + * else — the only thing you have to provide is an address where all + * the indices are stored and the number of elements to be set. */ void set (const size_type n_elements, const size_type *indices, const TrilinosScalar *values); /** - * 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. */ 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. */ 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. */ void add (const size_type n_elements, const size_type *indices, const TrilinosScalar *values); /** - * Multiply the entire vector by - * a fixed factor. + * Multiply the entire vector by a fixed factor. */ VectorBase &operator *= (const TrilinosScalar factor); /** - * Divide the entire vector by a - * fixed factor. + * Divide the entire vector by a fixed factor. */ VectorBase &operator /= (const TrilinosScalar factor); /** - * Add the given vector to the - * present one. + * Add the given vector to the present one. */ VectorBase &operator += (const VectorBase &V); /** - * Subtract the given vector from - * the present one. + * Subtract the given vector from the present one. */ VectorBase &operator -= (const VectorBase &V); /** - * 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 TrilinosScalar s); /** - * Simple vector addition, equal - * to the operator - * +=. + * Simple vector addition, equal to the operator +=. * - * Though, if the second argument - * allow_different_maps - * is set, then it is possible to - * add data from a different map. + * Though, if the second argument allow_different_maps is set, + * then it is possible to add data from a different map. */ void add (const VectorBase &V, const bool allow_different_maps = false); /** - * 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 TrilinosScalar a, const VectorBase &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 TrilinosScalar a, const VectorBase &V, @@ -878,17 +693,14 @@ namespace TrilinosWrappers const VectorBase &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 TrilinosScalar s, const VectorBase &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 TrilinosScalar s, const TrilinosScalar a, @@ -904,9 +716,8 @@ namespace TrilinosWrappers const VectorBase &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 TrilinosScalar s, const TrilinosScalar a, @@ -917,26 +728,20 @@ namespace TrilinosWrappers const VectorBase &X); /** - * 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 VectorBase &scaling_factors); /** - * Assignment *this = - * a*V. + * Assignment *this = a*V. */ void equ (const TrilinosScalar a, const VectorBase &V); /** - * Assignment *this = a*V + - * b*W. + * Assignment *this = a*V + b*W. */ void equ (const TrilinosScalar a, const VectorBase &V, @@ -944,21 +749,14 @@ namespace TrilinosWrappers const VectorBase &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 VectorBase &a, const VectorBase &b); @@ -971,46 +769,34 @@ namespace TrilinosWrappers //@{ /** - * Return a const reference to the - * underlying Trilinos - * Epetra_MultiVector class. + * Return a const reference to the underlying Trilinos Epetra_MultiVector + * class. */ const Epetra_MultiVector &trilinos_vector () const; /** - * Return a (modifyable) reference to - * the underlying Trilinos + * Return a (modifyable) reference to the underlying Trilinos * Epetra_FEVector class. */ Epetra_FEVector &trilinos_vector (); /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the parallel - * partitioning of the vector. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the parallel partitioning of the vector. */ const Epetra_Map &vector_partitioner () const; /** - * Output of vector in - * user-defined format in analogy - * to the dealii::Vector - * class. + * Output of vector in user-defined format in analogy to the + * dealii::Vector class. */ void print (const char *format = 0) const; /** - * Print to a stream. @p - * precision denotes the desired - * precision with which values - * shall be printed, @p - * scientific whether scientific - * notation shall be used. If @p - * across is @p true then the - * vector is printed in a line, - * while if @p false then the - * elements are printed on a + * Print to a stream. @p precision denotes the desired precision with + * which values shall be printed, @p scientific whether scientific + * notation shall be used. If @p across is @p true then the vector is + * printed in a line, while if @p false then the elements are printed on a * separate line each. */ void print (std::ostream &out, @@ -1019,42 +805,27 @@ namespace TrilinosWrappers const bool across = true) const; /** - * 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. Note that the vectors - * need to be of the same size - * and base on the same map. + * 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. Note + * that the vectors need to be of the same size and base on the same map. * - * 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 + * 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. */ void swap (VectorBase &v); /** - * Estimate for the memory - * consumption in bytes. + * Estimate for the memory consumption in bytes. */ std::size_t memory_consumption () const; /** - * Return a reference to the MPI - * communicator object in use with this + * Return a reference to the MPI communicator object in use with this * object. */ const MPI_Comm &get_mpi_communicator () const; @@ -1091,57 +862,45 @@ namespace TrilinosWrappers private: /** - * Trilinos doesn't allow to - * mix additions to matrix - * entries and overwriting them - * (to make synchronisation of - * parallel computations - * simpler). The way we do it - * is to, for each access - * operation, store whether it - * is an insertion or an - * addition. If the previous - * one was of different type, - * then we first have to flush - * the Trilinos buffers; - * otherwise, we can simply go - * on. Luckily, Trilinos has - * an object for this which - * does already all the - * parallel communications in - * such a case, so we simply - * use their model, which - * stores whether the last - * operation was an addition or - * an insertion. + * Trilinos doesn't allow to mix additions to matrix entries and + * overwriting them (to make synchronisation of parallel computations + * simpler). The way we do it is to, for each access operation, store + * whether it is an insertion or an addition. If the previous one was of + * different type, then we first have to flush the Trilinos buffers; + * otherwise, we can simply go on. Luckily, Trilinos has an object for + * this which does already all the parallel communications in such a case, + * so we simply use their model, which stores whether the last operation + * was an addition or an insertion. */ Epetra_CombineMode last_action; /** - * A boolean variable to hold - * information on whether the - * vector is compressed or not. + * A boolean variable to hold information on whether the vector is + * compressed or not. */ bool compressed; /** - * Whether this vector has ghost elements. This is true - * on all processors even if only one of them has any - * ghost elements. + * Whether this vector has ghost elements. This is true on all processors + * even if only one of them has any ghost elements. */ bool has_ghosts; /** - * An Epetra distibuted vector - * type. Requires an existing - * Epetra_Map for storing data. + * An Epetra distibuted vector type. Requires an existing Epetra_Map for + * storing data. */ std_cxx1x::shared_ptr vector; + /** + * A vector object in Trilinos to be used for collecting the non-local + * elements if the vector was constructed with an additional IndexSet + * describing ghost elements. + */ + std_cxx1x::shared_ptr nonlocal_vector; /** - * Make the reference class a - * friend. + * Make the reference class a friend. */ friend class internal::VectorReference; friend class Vector; @@ -1440,58 +1199,6 @@ namespace TrilinosWrappers - inline - void - VectorBase::compress (::dealii::VectorOperation::values given_last_action) - { - //Select which mode to send to - //Trilinos. Note that we use last_action - //if available and ignore what the user - //tells us to detect wrongly mixed - //operations. Typically given_last_action - //is only used on machines that do not - //execute an operation (because they have - //no own cells for example). - Epetra_CombineMode mode = last_action; - if (last_action == Zero) - { - if (given_last_action==::dealii::VectorOperation::add) - mode = Add; - else if (given_last_action==::dealii::VectorOperation::insert) - mode = Insert; - } - -#ifdef DEBUG -# ifdef DEAL_II_WITH_MPI - // check that every process has decided - // to use the same mode. This will - // otherwise result in undefined - // behaviour in the call to - // GlobalAssemble(). - double double_mode = mode; - Utilities::MPI::MinMaxAvg result - = Utilities::MPI::min_max_avg (double_mode, - dynamic_cast - (&vector_partitioner().Comm())->GetMpiComm()); - Assert(result.max-result.min<1e-5, - ExcMessage ("Not all processors agree whether the last operation on " - "this vector was an addition or a set operation. This will " - "prevent the compress() operation from succeeding.")); - -# endif -#endif - - // Now pass over the information about - // what we did last to the vector. - const int ierr = vector->GlobalAssemble(mode); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - last_action = Zero; - - compressed = true; - } - - - inline void VectorBase::compress () @@ -1574,7 +1281,9 @@ namespace TrilinosWrappers { const size_type row = indices[i]; const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast(row)); - if (local_row == -1) + if (local_row != -1) + (*vector)[0][local_row] = values[i]; + else { const int ierr = vector->ReplaceGlobalValues (1, (const TrilinosWrappers::types::int_type *)(&row), @@ -1582,8 +1291,11 @@ namespace TrilinosWrappers AssertThrow (ierr == 0, ExcTrilinosError(ierr)); compressed = false; } - else - (*vector)[0][local_row] = values[i]; + // in set operation, do not use the pre-allocated vector for nonlocal + // entries even if it exists. This is to ensure that we really only + // set the elements touched by the set() method and not all contained + // in the nonlocal entries vector (there is no way to distinguish them + // on the receiving processor) } } @@ -1642,7 +1354,9 @@ namespace TrilinosWrappers { const size_type row = indices[i]; const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast(row)); - if (local_row == -1) + if (local_row != -1) + (*vector)[0][local_row] += values[i]; + else if (nonlocal_vector.get() == 0) { const int ierr = vector->SumIntoGlobalValues (1, (const TrilinosWrappers::types::int_type *)(&row), @@ -1651,7 +1365,17 @@ namespace TrilinosWrappers compressed = false; } else - (*vector)[0][local_row] += values[i]; + { + // use pre-allocated vector for non-local entries if it exists for + // addition operation + const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast(row)); + Assert(my_row != -1, + ExcMessage("Attempted to write into off-processor vector entry " + "that has not be specified as being writable upon " + "initialization")); + (*nonlocal_vector)[0][my_row] += values[i]; + compressed = false; + } } } diff --git a/deal.II/source/lac/block_sparsity_pattern.cc b/deal.II/source/lac/block_sparsity_pattern.cc index 0536414ef5..d07c32c912 100644 --- a/deal.II/source/lac/block_sparsity_pattern.cc +++ b/deal.II/source/lac/block_sparsity_pattern.cc @@ -741,6 +741,27 @@ namespace TrilinosWrappers + BlockSparsityPattern:: + BlockSparsityPattern (const std::vector &row_parallel_partitioning, + const std::vector &col_parallel_partitioning, + const std::vector &writable_rows, + const MPI_Comm &communicator) + : + BlockSparsityPatternBase + (row_parallel_partitioning.size(), + col_parallel_partitioning.size()) + { + for (size_type i=0; iblock(i,j).reinit(row_parallel_partitioning[i], + col_parallel_partitioning[j], + writable_rows[i], + communicator); + this->collect_sizes(); + } + + + void BlockSparsityPattern::reinit (const std::vector &row_block_sizes, const std::vector &col_block_sizes) diff --git a/deal.II/source/lac/trilinos_block_vector.cc b/deal.II/source/lac/trilinos_block_vector.cc index 4f21ce6136..d1b5a6f536 100644 --- a/deal.II/source/lac/trilinos_block_vector.cc +++ b/deal.II/source/lac/trilinos_block_vector.cc @@ -146,7 +146,8 @@ namespace TrilinosWrappers void BlockVector::reinit (const std::vector ¶llel_partitioning, const std::vector &ghost_values, - const MPI_Comm &communicator) + const MPI_Comm &communicator, + const bool vector_writable) { const size_type no_blocks = parallel_partitioning.size(); std::vector block_sizes (no_blocks); @@ -161,7 +162,8 @@ namespace TrilinosWrappers components.resize(n_blocks()); for (size_type i=0; iGraph())); + compress(); } @@ -475,6 +479,7 @@ namespace TrilinosWrappers { // release memory before reallocation matrix.reset(); + nonlocal_matrix.reset(); // if we want to exchange data, build a usual Trilinos sparsity pattern // and let that handle the exchange. otherwise, manually create a @@ -582,9 +587,10 @@ namespace TrilinosWrappers (Copy, sparsity_pattern.trilinos_sparsity_pattern(), false)); if (sparsity_pattern.nonlocal_graph.get() != 0) - { - nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph)); - } + nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph)); + else + nonlocal_matrix.reset (); + compress(); last_action = Zero; } @@ -598,6 +604,11 @@ namespace TrilinosWrappers matrix.reset (); matrix.reset (new Epetra_FECrsMatrix (Copy, sparse_matrix.trilinos_sparsity_pattern(), false)); + if (sparse_matrix.nonlocal_matrix != 0) + nonlocal_matrix.reset (new Epetra_CrsMatrix + (Copy, sparse_matrix.nonlocal_matrix->Graph())); + else + nonlocal_matrix.reset(); compress(); } @@ -755,6 +766,7 @@ namespace TrilinosWrappers const Epetra_CrsGraph *graph = &input_matrix.Graph(); + nonlocal_matrix.reset(); matrix.reset (); matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false)); @@ -801,8 +813,10 @@ namespace TrilinosWrappers // flush buffers int ierr; - if (nonlocal_matrix.get() != 0) + if (nonlocal_matrix.get() != 0 && mode == Add) { + // do only export in case of an add() operation, otherwise the owning + // processor must have set the correct entry nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap()); Epetra_Export exporter(nonlocal_matrix->RowMap(), matrix->RowMap()); ierr = matrix->Export(*nonlocal_matrix, exporter, mode); diff --git a/deal.II/source/lac/trilinos_sparsity_pattern.cc b/deal.II/source/lac/trilinos_sparsity_pattern.cc index 85ca1111fc..b406c7300f 100644 --- a/deal.II/source/lac/trilinos_sparsity_pattern.cc +++ b/deal.II/source/lac/trilinos_sparsity_pattern.cc @@ -544,7 +544,7 @@ namespace TrilinosWrappers if (Utilities::MPI::n_mpi_processes(communicator) > 1) { Epetra_Map nonlocal_map = - nonlocal_partitioner.make_trilinos_map(communicator, false); + nonlocal_partitioner.make_trilinos_map(communicator, true); nonlocal_graph.reset(new Epetra_CrsGraph(Copy, nonlocal_map, 0)); } else diff --git a/deal.II/source/lac/trilinos_vector.cc b/deal.II/source/lac/trilinos_vector.cc index ea330e88e3..4ff1d598e8 100644 --- a/deal.II/source/lac/trilinos_vector.cc +++ b/deal.II/source/lac/trilinos_vector.cc @@ -186,6 +186,8 @@ namespace TrilinosWrappers Vector::reinit (const Epetra_Map &input_map, const bool fast) { + nonlocal_vector.reset(); + if (vector->Map().SameAs(input_map)==false) vector.reset (new Epetra_FEVector(input_map)); else if (fast == false) @@ -205,6 +207,8 @@ namespace TrilinosWrappers const MPI_Comm &communicator, const bool fast) { + nonlocal_vector.reset(); + Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator, true); reinit (map, fast); @@ -217,13 +221,11 @@ namespace TrilinosWrappers const bool fast, const bool allow_different_maps) { - // In case we do not allow to - // have different maps, this - // call means that we have to - // reset the vector. So clear - // the vector, initialize our - // map with the map in v, and - // generate the vector. + nonlocal_vector.reset(); + + // In case we do not allow to have different maps, this call means that + // we have to reset the vector. So clear the vector, initialize our map + // with the map in v, and generate the vector. if (allow_different_maps == false) { if (vector->Map().SameAs(v.vector->Map()) == false) @@ -250,14 +252,10 @@ namespace TrilinosWrappers } } - // Otherwise, we have to check - // that the two vectors are - // already of the same size, - // create an object for the data - // exchange and then insert all - // the data. The first assertion - // is only a check whether the - // user knows what she is doing. + // Otherwise, we have to check that the two vectors are already of the + // same size, create an object for the data exchange and then insert all + // the data. The first assertion is only a check whether the user knows + // what she is doing. else { Assert (fast == false, @@ -284,19 +282,16 @@ namespace TrilinosWrappers Vector::reinit (const BlockVector &v, const bool import_data) { - // In case we do not allow to - // have different maps, this - // call means that we have to - // reset the vector. So clear - // the vector, initialize our - // map with the map in v, and - // generate the vector. + nonlocal_vector.reset(); + + // In case we do not allow to have different maps, this call means that + // we have to reset the vector. So clear the vector, initialize our map + // with the map in v, and generate the vector. if (v.n_blocks() == 0) return; - // create a vector that holds all the elements - // contained in the block vector. need to - // manually create an Epetra_Map. + // create a vector that holds all the elements contained in the block + // vector. need to manually create an Epetra_Map. size_type n_elements = 0, added_elements = 0, block_offset = 0; for (size_type block=0; block 1) + { + Epetra_Map nonlocal_map = + nonlocal_entries.make_trilinos_map(communicator, true); + nonlocal_vector.reset(new Epetra_MultiVector(nonlocal_map, 1)); + } + } } Vector & Vector::operator = (const Vector &v) { - // distinguish three cases. First case: both - // vectors have the same layout (just need to - // copy the local data, not reset the memory - // and the underlying Epetra_Map). The third - // case means that we have to rebuild the - // calling vector. + // distinguish three cases. First case: both vectors have the same + // layout (just need to copy the local data, not reset the memory and + // the underlying Epetra_Map). The third case means that we have to + // rebuild the calling vector. if (vector->Map().SameAs(v.vector->Map())) { *vector = *v.vector; + if (v.nonlocal_vector.get() != 0) + nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1)); last_action = Zero; } // Second case: vectors have the same global @@ -389,6 +409,9 @@ namespace TrilinosWrappers has_ghosts = v.has_ghosts; } + if (v.nonlocal_vector.get() != 0) + nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1)); + return *this; } @@ -397,6 +420,8 @@ namespace TrilinosWrappers Vector & Vector::operator = (const TrilinosWrappers::Vector &v) { + nonlocal_vector.reset(); + Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size())); Epetra_Import data_exchange (vector->Map(), v.vector->Map()); diff --git a/deal.II/source/lac/trilinos_vector_base.cc b/deal.II/source/lac/trilinos_vector_base.cc index 48b0afcd0c..95a9b1387b 100644 --- a/deal.II/source/lac/trilinos_vector_base.cc +++ b/deal.II/source/lac/trilinos_vector_base.cc @@ -21,6 +21,7 @@ # include # include +# include DEAL_II_NAMESPACE_OPEN @@ -55,13 +56,9 @@ namespace TrilinosWrappers Assert (index < vector.size(), ExcIndexRange (index, 0, vector.size())); - // Trilinos allows for vectors - // to be referenced by the [] or - // () operators but only () - // checks index bounds. We check - // these bounds by ourselves, so - // we can use []. Note that we - // can only get local values. + // Trilinos allows for vectors to be referenced by the [] or () + // operators but only () checks index bounds. We check these bounds by + // ourselves, so we can use []. Note that we can only get local values. const TrilinosWrappers::types::int_type local_index = vector.vector->Map().LID(static_cast(index)); @@ -112,9 +109,8 @@ namespace TrilinosWrappers void VectorBase::clear () { - // When we clear the vector, - // reset the pointer and generate - // an empty vector. + // When we clear the vector, reset the pointer and generate an empty + // vector. #ifdef DEAL_II_WITH_MPI Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF)); #else @@ -186,18 +182,67 @@ namespace TrilinosWrappers + void + VectorBase::compress (::dealii::VectorOperation::values given_last_action) + { + //Select which mode to send to Trilinos. Note that we use last_action if + //available and ignore what the user tells us to detect wrongly mixed + //operations. Typically given_last_action is only used on machines that do + //not execute an operation (because they have no own cells for example). + Epetra_CombineMode mode = last_action; + if (last_action == Zero) + { + if (given_last_action==::dealii::VectorOperation::add) + mode = Add; + else if (given_last_action==::dealii::VectorOperation::insert) + mode = Insert; + } + +#ifdef DEBUG +# ifdef DEAL_II_WITH_MPI + // check that every process has decided to use the same mode. This will + // otherwise result in undefined behaviour in the call to + // GlobalAssemble(). + double double_mode = mode; + Utilities::MPI::MinMaxAvg result + = Utilities::MPI::min_max_avg (double_mode, + dynamic_cast + (&vector_partitioner().Comm())->GetMpiComm()); + Assert(result.max-result.min<1e-5, + ExcMessage ("Not all processors agree whether the last operation on " + "this vector was an addition or a set operation. This will " + "prevent the compress() operation from succeeding.")); + +# endif +#endif + + // Now pass over the information about what we did last to the vector. + int ierr = 0; + if (nonlocal_vector.get() == 0 || mode != Add) + ierr = vector->GlobalAssemble(mode); + else + { + Epetra_Export exporter(nonlocal_vector->Map(), vector->Map()); + ierr = vector->Export(*nonlocal_vector, exporter, mode); + nonlocal_vector->PutScalar(0.); + } + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + last_action = Zero; + + compressed = true; + } + + + TrilinosScalar VectorBase::el (const size_type index) const { - // Extract local indices in - // the vector. + // Extract local indices in the vector. TrilinosWrappers::types::int_type trilinos_i = vector->Map().LID(static_cast(index)); TrilinosScalar value = 0.; - // If the element is not - // present on the current - // processor, we can't + // If the element is not present on the current processor, we can't // continue. Just print out 0. // TODO: Is this reasonable? @@ -218,17 +263,13 @@ namespace TrilinosWrappers TrilinosScalar VectorBase::operator () (const size_type index) const { - // Extract local indices in - // the vector. + // Extract local indices in the vector. TrilinosWrappers::types::int_type trilinos_i = vector->Map().LID(static_cast(index)); TrilinosScalar value = 0.; - // If the element is not present - // on the current processor, we - // can't continue. This is the - // main difference to the el() - // function. + // If the element is not present on the current processor, we can't + // continue. This is the main difference to the el() function. if (trilinos_i == -1 ) { Assert (false, ExcAccessToNonlocalElement(index, local_range().first, diff --git a/tests/trilinos/assemble_matrix_parallel_02.cc b/tests/trilinos/assemble_matrix_parallel_02.cc index 954afca20a..96cabf815a 100644 --- a/tests/trilinos/assemble_matrix_parallel_02.cc +++ b/tests/trilinos/assemble_matrix_parallel_02.cc @@ -88,9 +88,14 @@ namespace Assembly { struct Data { + Data(const bool assemble_reference) + : + assemble_reference(assemble_reference) + {} std::vector local_dof_indices; FullMatrix local_matrix; Vector local_rhs; + const bool assemble_reference; }; } } @@ -246,18 +251,27 @@ void LaplaceProblem::setup_system () (FilteredIterator::active_cell_iterator> const &)> > (std_cxx1x::bind(&LaplaceProblem::get_conflict_indices, this,std_cxx1x::_1))); - TrilinosWrappers::SparsityPattern csp; - IndexSet locally_owned = dof_handler.locally_owned_dofs(), relevant_set; - DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set); - csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); - DoFTools::make_sparsity_pattern (dof_handler, csp, - constraints, false); - csp.compress(); - reference_matrix.reinit (csp); - test_matrix.reinit (csp); - - reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); - test_rhs.reinit (reference_rhs); + IndexSet locally_owned = dof_handler.locally_owned_dofs(); + { + TrilinosWrappers::SparsityPattern csp; + csp.reinit(locally_owned, locally_owned, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + reference_matrix.reinit (csp); + reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); + } + { + TrilinosWrappers::SparsityPattern csp; + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set); + csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + test_matrix.reinit (csp); + test_rhs.reinit (locally_owned, relevant_set, MPI_COMM_WORLD, true); + } } @@ -310,9 +324,14 @@ template void LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) { - constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, - data.local_dof_indices, - test_matrix, test_rhs); + if (data.assemble_reference) + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + reference_matrix, reference_rhs); + else + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); } @@ -320,10 +339,10 @@ LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) template void LaplaceProblem::assemble_reference () { - test_matrix = 0; - test_rhs = 0; + reference_matrix = 0; + reference_rhs = 0; - Assembly::Copy::Data copy_data; + Assembly::Copy::Data copy_data(true); Assembly::Scratch::Data assembly_data(fe, quadrature); for (unsigned int color=0; color::assemble_reference () local_assemble(*p, assembly_data, copy_data); copy_local_to_global(copy_data); } - test_matrix.compress(VectorOperation::add); - test_rhs.compress(VectorOperation::add); - - reference_matrix.add(1., test_matrix); - reference_rhs = test_rhs; + reference_matrix.compress(VectorOperation::add); + reference_rhs.compress(VectorOperation::add); } @@ -361,7 +377,7 @@ void LaplaceProblem::assemble_test () this, std_cxx1x::_1), Assembly::Scratch::Data(fe, quadrature), - Assembly::Copy::Data (), + Assembly::Copy::Data (false), 2*multithread_info.n_threads(), 1); test_matrix.compress(VectorOperation::add); @@ -401,9 +417,7 @@ void LaplaceProblem::run () { if (cycle == 0) { - GridGenerator::hyper_cube(triangulation, 0, 1/* , - Point(), - 0.5, 1., (dim==3) ? 96 : 12, false*/); + GridGenerator::hyper_cube(triangulation, 0, 1); triangulation.refine_global(6); } diff --git a/tests/trilinos/assemble_matrix_parallel_03.cc b/tests/trilinos/assemble_matrix_parallel_03.cc index e61846e1b5..ce114164bd 100644 --- a/tests/trilinos/assemble_matrix_parallel_03.cc +++ b/tests/trilinos/assemble_matrix_parallel_03.cc @@ -89,9 +89,14 @@ namespace Assembly { struct Data { + Data(const bool assemble_reference) + : + assemble_reference(assemble_reference) + {} std::vector local_dof_indices; FullMatrix local_matrix; Vector local_rhs; + const bool assemble_reference; }; } } @@ -247,18 +252,27 @@ void LaplaceProblem::setup_system () (FilteredIterator::active_cell_iterator> const &)> > (std_cxx1x::bind(&LaplaceProblem::get_conflict_indices, this,std_cxx1x::_1))); - TrilinosWrappers::SparsityPattern csp; - IndexSet locally_owned = dof_handler.locally_owned_dofs(), relevant_set; - DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set); - csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); - DoFTools::make_sparsity_pattern (dof_handler, csp, - constraints, false); - csp.compress(); - reference_matrix.reinit (csp); - test_matrix.reinit (csp); - - reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); - test_rhs.reinit (reference_rhs); + IndexSet locally_owned = dof_handler.locally_owned_dofs(); + { + TrilinosWrappers::SparsityPattern csp; + csp.reinit(locally_owned, locally_owned, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + reference_matrix.reinit (csp); + reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); + } + { + TrilinosWrappers::SparsityPattern csp; + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set); + csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + test_matrix.reinit (csp); + test_rhs.reinit (locally_owned, relevant_set, MPI_COMM_WORLD, true); + } } @@ -311,9 +325,14 @@ template void LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) { - constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, - data.local_dof_indices, - test_matrix, test_rhs); + if (data.assemble_reference) + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + reference_matrix, reference_rhs); + else + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); } @@ -321,10 +340,10 @@ LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) template void LaplaceProblem::assemble_reference () { - test_matrix = 0; - test_rhs = 0; + reference_matrix = 0; + reference_rhs = 0; - Assembly::Copy::Data copy_data; + Assembly::Copy::Data copy_data(true); Assembly::Scratch::Data assembly_data(fe, quadrature); for (unsigned int color=0; color::assemble_reference () local_assemble(*p, assembly_data, copy_data); copy_local_to_global(copy_data); } - test_matrix.compress(VectorOperation::add); - test_rhs.compress(VectorOperation::add); - - reference_matrix.add(1., test_matrix); - reference_rhs = test_rhs; + reference_matrix.compress(VectorOperation::add); + reference_rhs.compress(VectorOperation::add); } @@ -362,12 +378,13 @@ void LaplaceProblem::assemble_test () this, std_cxx1x::_1), Assembly::Scratch::Data(fe, quadrature), - Assembly::Copy::Data (), + Assembly::Copy::Data (false), 2*multithread_info.n_threads(), 1); test_matrix.compress(VectorOperation::add); test_rhs.compress(VectorOperation::add); test_matrix.compress(VectorOperation::add); + test_rhs.compress(VectorOperation::add); test_matrix.add(-1, reference_matrix); test_matrix.compress(VectorOperation::add); diff --git a/tests/trilinos/assemble_matrix_parallel_04.cc b/tests/trilinos/assemble_matrix_parallel_04.cc index b7745d6e2a..970137a133 100644 --- a/tests/trilinos/assemble_matrix_parallel_04.cc +++ b/tests/trilinos/assemble_matrix_parallel_04.cc @@ -90,9 +90,14 @@ namespace Assembly { struct Data { + Data(const bool assemble_reference) + : + assemble_reference(assemble_reference) + {} std::vector local_dof_indices; FullMatrix local_matrix; Vector local_rhs; + const bool assemble_reference; }; } } @@ -265,15 +270,22 @@ void LaplaceProblem::setup_system () relevant_set[1] = relevant_total.get_view(dofs_per_block[0], dof_handler.n_dofs()); - csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); - DoFTools::make_sparsity_pattern (dof_handler, csp, - constraints, false); - csp.compress(); - reference_matrix.reinit (csp); - test_matrix.reinit (csp); - - reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); - test_rhs.reinit (reference_rhs); + { + csp.reinit(locally_owned, locally_owned, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + reference_matrix.reinit (csp); + reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); + } + { + csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + test_matrix.reinit (csp); + test_rhs.reinit (locally_owned, relevant_set, MPI_COMM_WORLD, true); + } } @@ -328,9 +340,14 @@ template void LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) { - constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, - data.local_dof_indices, - test_matrix, test_rhs); + if (data.assemble_reference) + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + reference_matrix, reference_rhs); + else + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); } @@ -338,10 +355,10 @@ LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) template void LaplaceProblem::assemble_reference () { - test_matrix = 0; - test_rhs = 0; + reference_matrix = 0; + reference_rhs = 0; - Assembly::Copy::Data copy_data; + Assembly::Copy::Data copy_data(true); Assembly::Scratch::Data assembly_data(fe, quadrature); for (unsigned int color=0; color::assemble_reference () local_assemble(*p, assembly_data, copy_data); copy_local_to_global(copy_data); } - test_matrix.compress(VectorOperation::add); - test_rhs.compress(VectorOperation::add); - - reference_matrix.add(1., test_matrix); - reference_rhs = test_rhs; + reference_matrix.compress(VectorOperation::add); + reference_rhs.compress(VectorOperation::add); } @@ -379,7 +393,7 @@ void LaplaceProblem::assemble_test () this, std_cxx1x::_1), Assembly::Scratch::Data(fe, quadrature), - Assembly::Copy::Data (), + Assembly::Copy::Data (false), 2*multithread_info.n_threads(), 1); test_matrix.compress(VectorOperation::add);