-// ---------------------------------------------------------------------
+ // ---------------------------------------------------------------------
// $Id$
//
// Copyright (C) 2013 by the deal.II authors
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.
<br>
(Martin Kronbichler, 2013/12/23)
</li>
BlockSparsityPattern (const std::vector<IndexSet> ¶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<IndexSet> &row_parallel_partitioning,
+ const std::vector<IndexSet> &column_parallel_partitioning,
+ const std::vector<IndexSet> &writeable_rows,
+ const MPI_Comm &communicator = MPI_COMM_WORLD);
+
/**
* Resize the matrix to a tensor product of matrices with dimensions
* defined by the arguments.
{
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<Vector> 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;
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<Epetra_Map> &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<IndexSet> &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
- * <tt>num_blocks</tt>
- * 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 <tt>num_blocks</tt> 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 <tt>n.size()</tt> and
- * initialize each block with
- * <tt>n[i]</tt> zero elements.
+ * Constructor. Set the number of blocks to <tt>n.size()</tt> and
+ * initialize each block with <tt>n[i]</tt> zero elements.
*
* References BlockVector.reinit().
*/
explicit BlockVector (const std::vector<size_type> &N);
/**
- * Constructor. Set the number of
- * blocks to
- * <tt>n.size()</tt>. 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
- * <tt>std::vector</tt> 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 <tt>n.size()</tt>. 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
+ * <tt>std::vector</tt> class, but the first argument is needed in order
+ * to know how to subdivide the block vector into different blocks.
*/
template <typename InputIterator>
BlockVector (const std::vector<size_type> &n,
*
* @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;
using BlockVectorBase<Vector>::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 <typename Number>
BlockVector &
operator = (const ::dealii::BlockVector<Number> &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 <tt>fast==false</tt>, the vector
- * is filled with zeros.
+ * If <tt>fast==false</tt>, the vector is filled with zeros.
*/
void reinit (const std::vector<Epetra_Map> &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 <tt>fast==false</tt>, the vector
* is filled with zeros.
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 <tt>fast==false</tt>, the vector
* is filled with zeros.
void reinit (const MPI::BlockVector &V);
/**
- * Change the dimension to that
- * of the vector <tt>V</tt>. The same
- * applies as for the other
- * reinit() function.
+ * Change the dimension to that of the vector <tt>V</tt>. The same applies
+ * as for the other reinit() function.
*
- * The elements of <tt>V</tt> are not
- * copied, i.e. this function is
- * the same as calling <tt>reinit
- * (V.size(), fast)</tt>.
+ * The elements of <tt>V</tt> are not copied, i.e. this function is the
+ * same as calling <tt>reinit (V.size(), fast)</tt>.
*
- * 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
- * <tt>num_blocks</tt>. 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 <tt>collect_sizes</tt>
- * afterwards.
+ * Change the number of blocks to <tt>num_blocks</tt>. 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 <tt>collect_sizes</tt> afterwards.
*/
void reinit (const size_type num_blocks);
/**
- * Swap the contents of this
- * vector and the other vector
- * <tt>v</tt>. 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 <tt>v</tt>. 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
- * <tt>u.swap(v)</tt>, 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 <tt>u.swap(v)</tt>, again in analogy to standard functions.
*/
void swap (BlockVector &v);
{
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<Vector> 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;
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<Epetra_Map> ¶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<IndexSet> ¶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<IndexSet> ¶llel_partitioning,
const std::vector<IndexSet> &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
- * <tt>num_blocks</tt>
- * 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 <tt>num_blocks</tt> 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);
~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 <typename Number>
BlockVector &
operator = (const ::dealii::BlockVector<Number> &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 <tt>fast==false</tt>, the vector
- * is filled with zeros.
+ * If <tt>fast==false</tt>, the vector is filled with zeros.
*/
void reinit (const std::vector<Epetra_Map> ¶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 <tt>fast==false</tt>, the vector
- * is filled with zeros.
+ * If <tt>fast==false</tt>, the vector is filled with zeros.
*/
void reinit (const std::vector<IndexSet> ¶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<IndexSet> &partitioning,
const std::vector<IndexSet> &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 <tt>V</tt>. The same
- * applies as for the other
- * reinit() function.
+ * Change the dimension to that of the vector <tt>V</tt>. The same
+ * applies as for the other reinit() function.
*
- * The elements of <tt>V</tt> are not
- * copied, i.e. this function is
- * the same as calling <tt>reinit
- * (V.size(), fast)</tt>.
+ * The elements of <tt>V</tt> are not copied, i.e. this function is the
+ * same as calling <tt>reinit (V.size(), fast)</tt>.
*
- * 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
- * <tt>num_blocks</tt>. 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 <tt>collect_sizes</tt>
- * afterwards.
+ * Change the number of blocks to <tt>num_blocks</tt>. 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 <tt>collect_sizes</tt> 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<dim>::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<dim>::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);
*
* @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;
/**
- * 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
- * <tt>debug</tt> 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 <tt>debug</tt> mode, since it is quite
+ * expensive to keep track of all operations that lead to the need for
* compress().
*/
bool is_compressed () const;
bool has_ghost_elements() const;
/**
- * Swap the contents of this
- * vector and the other vector
- * <tt>v</tt>. 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 <tt>v</tt>. 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
- * <tt>u.swap(v)</tt>, 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 <tt>u.swap(v)</tt>, again in analogy to standard
+ * functions.
*/
void swap (BlockVector &v);
}
+
inline
BlockVector::BlockVector (const std::vector<IndexSet> ¶llel_partitioning,
const std::vector<IndexSet> &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)
{
}
+
inline
bool
BlockVector::has_ghost_elements() const
return ghosted;
}
+
+
inline
void
BlockVector::swap (BlockVector &v)
*
* 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
- * <b>different</b> 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 <b>different</b>
+ * rows of the matrix from several threads simultaneously under the
+ * following two conditions:
* <ul>
* <li> The matrix uses only one MPI process.
* <li> 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.
* </ul>
*
* @ingroup 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<TrilinosWrappers::types::int_type>(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]));
- // 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
* operations at the same time, for example by placing zero additions if
* necessary.
*
+ * <h3>Thread safety of Trilinos vectors</h3>
+ *
+ * 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
+ * <b>different</b> 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
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;
*/
//@{
/**
- * Default constructor that
- * generates an empty (zero size)
- * vector. The function
- * <tt>reinit()</tt> 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 <tt>reinit()</tt> 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);
~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
- * <tt>v</tt>. If <tt>fast</tt> is
- * not <tt>true</tt>, 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 <tt>v</tt>. If <tt>fast</tt> is not
+ * <tt>true</tt>, 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,
- * <tt>allow_different_maps</tt>,
- * 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
- * <tt>v</tt> is a distributed
- * vector. In this case, the variable
- * <tt>fast</tt> needs to be set to
- * <tt>false</tt>, since it does not
- * make sense to exchange data
- * between differently parallelized
- * vectors without touching the
- * elements.
+ * This function has a third argument, <tt>allow_different_maps</tt>,
+ * 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 <tt>v</tt>
+ * is a distributed vector. In this case, the variable <tt>fast</tt>
+ * needs to be set to <tt>false</tt>, 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 <typename Number>
Vector &
operator = (const ::dealii::Vector<Number> &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<dim>::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<dim>::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,
*/
//@{
/**
- * 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 <typename number>
const dealii::Vector<number> &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 <typename Number>
Vector (const Epetra_Map ¶llel_partitioning,
*/
//@{
/**
- * 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);
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 <typename Number>
Vector (const IndexSet ¶llel_partitioning,
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
- * <tt>fast</tt> 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
+ * <tt>fast</tt> 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);
//@}
};
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
- * <tt>reinit()</tt> will have to
- * give the vector the correct
+ * Default constructor that generates an empty (zero size) vector. The
+ * function <tt>reinit()</tt> 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 <typename Number>
explicit Vector (const dealii::Vector<Number> &v);
/**
- * Reinit function that resizes
- * the vector to the size
- * specified by <tt>n</tt>.
+ * Reinit function that resizes the vector to the size specified by
+ * <tt>n</tt>.
*/
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
- * <tt>fast</tt> 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
+ * <tt>fast</tt> 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
- * <tt>fast</tt> 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
+ * <tt>fast</tt> determines whether the vector should be filled with zero
* (false) or left untouched (true).
*/
void reinit (const IndexSet &input_map,
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 <typename Number>
Vector &
operator = (const ::dealii::Vector<Number> &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;
};
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);
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
- * <tt>v(i)=w(i)</tt>. 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 <tt>v(i)=w(i)</tt>. 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 <tt>s</tt>.
+ * Set the referenced element of the vector to <tt>s</tt>.
*/
const VectorReference &
operator = (const TrilinosScalar &s) const;
/**
- * Add <tt>s</tt> to the
- * referenced element of the
- * vector->
+ * Add <tt>s</tt> to the referenced element of the vector->
*/
const VectorReference &
operator += (const TrilinosScalar &s) const;
/**
- * Subtract <tt>s</tt> from the
- * referenced element of the
- * vector->
+ * Subtract <tt>s</tt> from the referenced element of the vector->
*/
const VectorReference &
operator -= (const TrilinosScalar &s) const;
/**
- * Multiply the referenced
- * element of the vector by
- * <tt>s</tt>.
+ * Multiply the referenced element of the vector by <tt>s</tt>.
*/
const VectorReference &
operator *= (const TrilinosScalar &s) const;
/**
- * Divide the referenced
- * element of the vector by
- * <tt>s</tt>.
+ * Divide the referenced element of the vector by <tt>s</tt>.
*/
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;
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;
{
public:
/**
- * Declare some of the standard
- * types used in all
- * containers. These types
- * parallel those in the
- * <tt>C</tt> standard libraries
+ * Declare some of the standard types used in all containers. These types
+ * parallel those in the <tt>C</tt> standard libraries
* <tt>vector<...></tt> class.
*/
typedef TrilinosScalar value_type;
//@{
/**
- * Default constructor that
- * generates an empty (zero size)
- * vector. The function
- * <tt>reinit()</tt> 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 <tt>reinit()</tt> 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);
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
- * (<code>Add</code> or
- * <code>Insert</code>) 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
+ * (<code>Add</code> or <code>Insert</code>) 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.
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
- * <tt>v=0</tt>. 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
+ * <tt>v=0</tt>. 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 <typename Number>
VectorBase &
operator = (const ::dealii::Vector<Number> &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 <code>(0,N)</code>, otherwise it will
- * be a pair <code>(i,i+n)</code>, where
- * <code>n=local_size()</code> and <code>i</code> 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 <code>(0,N)</code>, otherwise it will be a pair
+ * <code>(i,i+n)</code>, where <code>n=local_size()</code> and
+ * <code>i</code> 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
std::pair<size_type, size_type> 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().
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
- * <i>p</i>th root of the sum of
- * the <i>p</i>th powers of the
- * absolute values of the
- * elements.
+ * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the
+ * <i>p</i>th 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;
//@}
//@{
/**
- * 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 <code>el()</code> command.
+ * Provide read-only access to an element. This is equivalent to the
+ * <code>el()</code> 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().
*/
operator [] (const size_type index);
/**
- * Provide read-only access to an
- * element. This is equivalent to
- * the <code>el()</code> command.
+ * Provide read-only access to an element. This is equivalent to the
+ * <code>el()</code> command.
*
* Exactly the same as operator().
*/
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<size_type> &indices,
std::vector<TrilinosScalar> &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 <typename ForwardIterator, typename OutputIterator>
void extract_subvector_to (ForwardIterator indices_begin,
OutputIterator values_begin) const;
/**
- * Return the value of the vector
- * entry <i>i</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>i</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 <tt>vector<></tt> 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 <tt>vector<></tt> 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().
*/
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<size_type> &indices,
const std::vector<TrilinosScalar> &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<size_type> &indices,
const ::dealii::Vector<TrilinosScalar> &values);
//@{
/**
- * 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<size_type> &indices,
const std::vector<TrilinosScalar> &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<size_type> &indices,
const ::dealii::Vector<TrilinosScalar> &values);
/**
- * Take an address where
- * <tt>n_elements</tt> are stored
- * contiguously and add them into
- * the vector. Handles all cases
- * which are not covered by the
- * other two <tt>add()</tt>
- * functions above.
+ * Take an address where <tt>n_elements</tt> are stored contiguously and
+ * add them into the vector. Handles all cases which are not covered by
+ * the other two <tt>add()</tt> functions above.
*/
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 <tt>operator
- * +=</tt>.
+ * Simple vector addition, equal to the <tt>operator +=</tt>.
*
- * Though, if the second argument
- * <tt>allow_different_maps</tt>
- * is set, then it is possible to
- * add data from a different map.
+ * Though, if the second argument <tt>allow_different_maps</tt> 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. <tt>*this =
- * a*V</tt>.
+ * Simple addition of a multiple of a vector, i.e. <tt>*this = a*V</tt>.
*/
void add (const TrilinosScalar a,
const VectorBase &V);
/**
- * Multiple addition of scaled
- * vectors, i.e. <tt>*this = a*V +
- * b*W</tt>.
+ * Multiple addition of scaled vectors, i.e. <tt>*this = a*V + b*W</tt>.
*/
void add (const TrilinosScalar a,
const VectorBase &V,
const VectorBase &W);
/**
- * Scaling and simple vector
- * addition, i.e. <tt>*this =
- * s*(*this) + V</tt>.
+ * Scaling and simple vector addition, i.e. <tt>*this = s*(*this) +
+ * V</tt>.
*/
void sadd (const TrilinosScalar s,
const VectorBase &V);
/**
- * Scaling and simple addition,
- * i.e. <tt>*this = s*(*this) +
- * a*V</tt>.
+ * Scaling and simple addition, i.e. <tt>*this = s*(*this) + a*V</tt>.
*/
void sadd (const TrilinosScalar s,
const TrilinosScalar a,
const VectorBase &W);
/**
- * Scaling and multiple addition.
- * <tt>*this = s*(*this) + a*V +
- * b*W + c*X</tt>.
+ * Scaling and multiple addition. <tt>*this = s*(*this) + a*V + b*W +
+ * c*X</tt>.
*/
void sadd (const TrilinosScalar s,
const TrilinosScalar a,
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 <tt>*this =
- * a*V</tt>.
+ * Assignment <tt>*this = a*V</tt>.
*/
void equ (const TrilinosScalar a,
const VectorBase &V);
/**
- * Assignment <tt>*this = a*V +
- * b*W</tt>.
+ * Assignment <tt>*this = a*V + b*W</tt>.
*/
void equ (const TrilinosScalar a,
const VectorBase &V,
const VectorBase &W);
/**
- * Compute the elementwise ratio
- * of the two given vectors, that
- * is let <tt>this[i] =
- * a[i]/b[i]</tt>. This is useful
- * for example if you want to
- * compute the cellwise ratio of
- * true to estimated error.
+ * Compute the elementwise ratio of the two given vectors, that is let
+ * <tt>this[i] = a[i]/b[i]</tt>. This is useful for example if you want to
+ * compute the cellwise ratio of true to estimated error.
*
- * This vector is appropriately
- * scaled to hold the result.
+ * This vector is appropriately scaled to hold the result.
*
- * If any of the <tt>b[i]</tt> is
- * zero, the result is
- * undefined. No attempt is made
- * to catch such situations.
+ * If any of the <tt>b[i]</tt> is zero, the result is undefined. No
+ * attempt is made to catch such situations.
*/
void ratio (const VectorBase &a,
const VectorBase &b);
//@{
/**
- * 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<number>
- * class.
+ * Output of vector in user-defined format in analogy to the
+ * dealii::Vector<number> 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,
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
- * <tt>swap(u,v)</tt> that simply
- * calls <tt>u.swap(v)</tt>,
- * 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 <tt>swap(u,v)</tt> that
+ * simply calls <tt>u.swap(v)</tt>, 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;
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<Epetra_FEVector> 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<Epetra_MultiVector> nonlocal_vector;
/**
- * Make the reference class a
- * friend.
+ * Make the reference class a friend.
*/
friend class internal::VectorReference;
friend class Vector;
- 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<const Epetra_MpiComm *>
- (&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 ()
{
const size_type row = indices[i];
const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(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),
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)
}
}
{
const size_type row = indices[i];
const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(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),
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<TrilinosWrappers::types::int_type>(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;
+ }
}
}
+ BlockSparsityPattern::
+ BlockSparsityPattern (const std::vector<IndexSet> &row_parallel_partitioning,
+ const std::vector<IndexSet> &col_parallel_partitioning,
+ const std::vector<IndexSet> &writable_rows,
+ const MPI_Comm &communicator)
+ :
+ BlockSparsityPatternBase<SparsityPattern>
+ (row_parallel_partitioning.size(),
+ col_parallel_partitioning.size())
+ {
+ for (size_type i=0; i<row_parallel_partitioning.size(); ++i)
+ for (size_type j=0; j<col_parallel_partitioning.size(); ++j)
+ this->block(i,j).reinit(row_parallel_partitioning[i],
+ col_parallel_partitioning[j],
+ writable_rows[i],
+ communicator);
+ this->collect_sizes();
+ }
+
+
+
void
BlockSparsityPattern::reinit (const std::vector<size_type> &row_block_sizes,
const std::vector<size_type> &col_block_sizes)
void
BlockVector::reinit (const std::vector<IndexSet> ¶llel_partitioning,
const std::vector<IndexSet> &ghost_values,
- const MPI_Comm &communicator)
+ const MPI_Comm &communicator,
+ const bool vector_writable)
{
const size_type no_blocks = parallel_partitioning.size();
std::vector<size_type> block_sizes (no_blocks);
components.resize(n_blocks());
for (size_type i=0; i<n_blocks(); ++i)
- components[i].reinit(parallel_partitioning[i], ghost_values[i], communicator);
+ components[i].reinit(parallel_partitioning[i], ghost_values[i],
+ communicator, vector_writable);
collect_sizes();
}
void
SparseMatrix::copy_from (const SparseMatrix &m)
{
+ nonlocal_matrix.reset();
// check whether we need to update the
// partitioner or can just copy the data:
matrix.reset (new Epetra_FECrsMatrix(*m.matrix));
}
+ if (m.nonlocal_matrix.get() != 0)
+ nonlocal_matrix.reset(new Epetra_CrsMatrix(Copy, m.nonlocal_matrix->Graph()));
+
compress();
}
{
// 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
(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;
}
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();
}
const Epetra_CrsGraph *graph = &input_matrix.Graph();
+ nonlocal_matrix.reset();
matrix.reset ();
matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
// 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);
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
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)
const MPI_Comm &communicator,
const bool fast)
{
+ nonlocal_vector.reset();
+
Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator,
true);
reinit (map, fast);
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)
}
}
- // 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,
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<v.n_blocks(); ++block)
n_elements += v.block(block).local_size();
}
- void Vector::reinit(const IndexSet &local, const IndexSet &ghost, const MPI_Comm &communicator)
+ void Vector::reinit(const IndexSet &locally_owned_entries,
+ const IndexSet &ghost_entries,
+ const MPI_Comm &communicator,
+ const bool vector_writable)
{
- IndexSet parallel_partitioning = local;
- parallel_partitioning.add_indices(ghost);
- reinit(parallel_partitioning, communicator);
+ nonlocal_vector.reset();
+ if (vector_writable == false)
+ {
+ IndexSet parallel_partitioning = locally_owned_entries;
+ parallel_partitioning.add_indices(ghost_entries);
+ reinit(parallel_partitioning, communicator);
+ }
+ else
+ {
+ Epetra_Map map = locally_owned_entries.make_trilinos_map (communicator,
+ true);
+ Assert (map.IsOneToOne(),
+ ExcMessage("A writable vector must not have ghost entries in "
+ "its parallel partitioning"));
+ reinit (map);
+
+ IndexSet nonlocal_entries(ghost_entries);
+ nonlocal_entries.subtract_set(locally_owned_entries);
+ if (Utilities::MPI::n_mpi_processes(communicator) > 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
has_ghosts = v.has_ghosts;
}
+ if (v.nonlocal_vector.get() != 0)
+ nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
+
return *this;
}
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());
# include <cmath>
# include <Epetra_Import.h>
+# include <Epetra_Export.h>
DEAL_II_NAMESPACE_OPEN
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<TrilinosWrappers::types::int_type>(index));
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
+ 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<const Epetra_MpiComm *>
+ (&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<TrilinosWrappers::types::int_type>(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?
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<TrilinosWrappers::types::int_type>(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,
{
struct Data
{
+ Data(const bool assemble_reference)
+ :
+ assemble_reference(assemble_reference)
+ {}
std::vector<types::global_dof_index> local_dof_indices;
FullMatrix<double> local_matrix;
Vector<double> local_rhs;
+ const bool assemble_reference;
};
}
}
(FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &)> >
(std_cxx1x::bind(&LaplaceProblem<dim>::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);
+ }
}
void
LaplaceProblem<dim>::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);
}
template <int dim>
void LaplaceProblem<dim>::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<dim> assembly_data(fe, quadrature);
for (unsigned int color=0; color<graph.size(); ++color)
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);
}
this,
std_cxx1x::_1),
Assembly::Scratch::Data<dim>(fe, quadrature),
- Assembly::Copy::Data (),
+ Assembly::Copy::Data (false),
2*multithread_info.n_threads(),
1);
test_matrix.compress(VectorOperation::add);
{
if (cycle == 0)
{
- GridGenerator::hyper_cube(triangulation, 0, 1/* ,
- Point<dim>(),
- 0.5, 1., (dim==3) ? 96 : 12, false*/);
+ GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(6);
}
{
struct Data
{
+ Data(const bool assemble_reference)
+ :
+ assemble_reference(assemble_reference)
+ {}
std::vector<types::global_dof_index> local_dof_indices;
FullMatrix<double> local_matrix;
Vector<double> local_rhs;
+ const bool assemble_reference;
};
}
}
(FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &)> >
(std_cxx1x::bind(&LaplaceProblem<dim>::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);
+ }
}
void
LaplaceProblem<dim>::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);
}
template <int dim>
void LaplaceProblem<dim>::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<dim> assembly_data(fe, quadrature);
for (unsigned int color=0; color<graph.size(); ++color)
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);
}
this,
std_cxx1x::_1),
Assembly::Scratch::Data<dim>(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);
{
struct Data
{
+ Data(const bool assemble_reference)
+ :
+ assemble_reference(assemble_reference)
+ {}
std::vector<types::global_dof_index> local_dof_indices;
FullMatrix<double> local_matrix;
Vector<double> local_rhs;
+ const bool assemble_reference;
};
}
}
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);
+ }
}
void
LaplaceProblem<dim>::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);
}
template <int dim>
void LaplaceProblem<dim>::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<dim> assembly_data(fe, quadrature);
for (unsigned int color=0; color<graph.size(); ++color)
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);
}
this,
std_cxx1x::_1),
Assembly::Scratch::Data<dim>(fe, quadrature),
- Assembly::Copy::Data (),
+ Assembly::Copy::Data (false),
2*multithread_info.n_threads(),
1);
test_matrix.compress(VectorOperation::add);
projects / dealii-svn.git / commitdiff