]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix Trilinos vectors
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 29 Dec 2013 19:52:11 +0000 (19:52 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 29 Dec 2013 19:52:11 +0000 (19:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@32126 0785d39b-7218-0410-832d-ea1e28bc413d

16 files changed:
deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/block_sparsity_pattern.h
deal.II/include/deal.II/lac/trilinos_block_vector.h
deal.II/include/deal.II/lac/trilinos_parallel_block_vector.h
deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
deal.II/include/deal.II/lac/trilinos_vector.h
deal.II/include/deal.II/lac/trilinos_vector_base.h
deal.II/source/lac/block_sparsity_pattern.cc
deal.II/source/lac/trilinos_block_vector.cc
deal.II/source/lac/trilinos_sparse_matrix.cc
deal.II/source/lac/trilinos_sparsity_pattern.cc
deal.II/source/lac/trilinos_vector.cc
deal.II/source/lac/trilinos_vector_base.cc
tests/trilinos/assemble_matrix_parallel_02.cc
tests/trilinos/assemble_matrix_parallel_03.cc
tests/trilinos/assemble_matrix_parallel_04.cc

index e398c01f87cf74c1fc2d2d8e9491a2a389a863e7..8b705acd64638047a5739a5365ef347914720c57 100644 (file)
@@ -1,4 +1,4 @@
-// ---------------------------------------------------------------------
+ // ---------------------------------------------------------------------
 // $Id$
 //
 // Copyright (C) 2013 by the deal.II authors
@@ -106,7 +106,9 @@ inconvenience this causes.
   written into as an optional argument. This allows for pre-allocating all
   possible entries right away, which makes writing into the matrix from
   several threads possible (otherwise, only one processor at a time can write
-  off-processor data).
+  off-processor data). Similarly, TrilinosWrappers::MPI::Vector objects can
+  be initialized with hints to ghost elements for a writable vector that can
+  be added into from multiple threads.
   <br>
   (Martin Kronbichler, 2013/12/23)
   </li>
index e22741f97c971659a2df15d6a874524b901a8ea7..8082826f243a580dfbf5b9a2ae6255c13d19225b 100644 (file)
@@ -1060,6 +1060,22 @@ namespace TrilinosWrappers
     BlockSparsityPattern (const std::vector<IndexSet> &parallel_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.
index 9006b134c22d983fb62b9fd8dc8c44de4dc0c840..276f778902637c1292d34b091a64ad272fffcbca 100644 (file)
@@ -71,20 +71,17 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Typedef the base class for simpler
-     * access to its own typedefs.
+     * Typedef the base class for simpler access to its own typedefs.
      */
     typedef BlockVectorBase<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;
@@ -96,90 +93,59 @@ namespace TrilinosWrappers
     typedef BaseClass::const_iterator  const_iterator;
 
     /**
-     * Default constructor. Generate an
-     * empty vector without any blocks.
+     * Default constructor. Generate an empty vector without any blocks.
      */
     BlockVector ();
 
     /**
-     * Constructor. Generate a block
-     * vector with as many blocks as
-     * there are entries in Input_Maps.
-     * For this non-distributed vector,
-     * the %parallel partitioning is not
-     * used, just the global size of the
-     * partitioner.
+     * Constructor. Generate a block vector with as many blocks as there are
+     * entries in Input_Maps.  For this non-distributed vector, the %parallel
+     * partitioning is not used, just the global size of the partitioner.
      */
     explicit BlockVector (const std::vector<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,
@@ -196,8 +162,7 @@ namespace TrilinosWrappers
     *
     * @deprecated
     *
-    * See @ref GlossCompress "Compressing
-    * distributed objects" for more
+    * See @ref GlossCompress "Compressing distributed objects" for more
     * information.
     */
     void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
@@ -208,91 +173,60 @@ namespace TrilinosWrappers
     using BlockVectorBase<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.
@@ -302,13 +236,9 @@ namespace TrilinosWrappers
                  const bool                   fast = false);
 
     /**
-     * Reinitialize the BlockVector to
-     * contain as many blocks as there
-     * are elements in the first
-     * argument, and with the respective
-     * sizes. Since no distribution map
-     * is given, all vectors are local
-     * vectors.
+     * Reinitialize the BlockVector to contain as many blocks as there are
+     * elements in the first argument, and with the respective sizes. Since no
+     * distribution map is given, all vectors are local vectors.
      *
      * If <tt>fast==false</tt>, the vector
      * is filled with zeros.
@@ -324,76 +254,44 @@ namespace TrilinosWrappers
     void reinit (const MPI::BlockVector &V);
 
     /**
-     * Change the dimension to that
-     * of the vector <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);
 
index 7ec7a0e217891eac190dc548609ee3366ccd4f63..de56aa813d3f01ba91dd30c1a2811f54ffeb628b 100644 (file)
@@ -73,20 +73,17 @@ namespace TrilinosWrappers
     {
     public:
       /**
-       * Typedef the base class for simpler
-       * access to its own typedefs.
+       * Typedef the base class for simpler access to its own typedefs.
        */
       typedef BlockVectorBase<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;
@@ -98,62 +95,47 @@ namespace TrilinosWrappers
       typedef BaseClass::const_iterator  const_iterator;
 
       /**
-       * Default constructor. Generate an
-       * empty vector without any blocks.
+       * Default constructor. Generate an empty vector without any blocks.
        */
       BlockVector ();
 
       /**
-       * Constructor. Generate a block
-       * vector with as many blocks as
-       * there are entries in @p
-       * partitioning. Each Epetra_Map
-       * contains the layout of the
-       * distribution of data among the MPI
-       * processes.
+       * Constructor. Generate a block vector with as many blocks as there are
+       * entries in @p partitioning. Each Epetra_Map contains the layout of
+       * the distribution of data among the MPI processes.
        */
       explicit BlockVector (const std::vector<Epetra_Map> &parallel_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> &parallel_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> &parallel_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);
 
@@ -163,168 +145,124 @@ namespace TrilinosWrappers
       ~BlockVector ();
 
       /**
-       * Copy operator: fill all
-       * components of the vector that
-       * are locally stored with the
-       * given scalar value.
+       * Copy operator: fill all components of the vector that are locally
+       * stored with the given scalar value.
        */
       BlockVector &
       operator = (const value_type s);
 
       /**
-       * Copy operator for arguments of
-       * the same type.
+       * Copy operator for arguments of the same type.
        */
       BlockVector &
       operator = (const BlockVector &V);
 
       /**
-       * Copy operator for arguments of
-       * the localized Trilinos vector
-       * type.
+       * Copy operator for arguments of the localized Trilinos vector type.
        */
       BlockVector &
       operator = (const ::dealii::TrilinosWrappers::BlockVector &V);
 
       /**
-       * Another copy function. This
-       * one takes a deal.II block
-       * vector and copies it into a
-       * TrilinosWrappers block
-       * vector. Note that the number
-       * of blocks has to be the same
-       * in the vector as in the input
-       * vector. Use the reinit()
-       * command for resizing the
-       * BlockVector or for changing
-       * the internal structure of the
-       * block components.
+       * Another copy function. This one takes a deal.II block vector and
+       * copies it into a TrilinosWrappers block vector. Note that the number
+       * of blocks has to be the same in the vector as in the input
+       * vector. Use the reinit() command for resizing the BlockVector or for
+       * changing the internal structure of the block components.
        *
-       * Since Trilinos only works on
-       * doubles, this function is
-       * limited to accept only one
-       * possible number type in the
-       * deal.II vector.
+       * Since Trilinos only works on doubles, this function is limited to
+       * accept only one possible number type in the deal.II vector.
        */
       template <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> &parallel_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> &parallel_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);
@@ -335,8 +273,7 @@ namespace TrilinosWrappers
        *
        * @deprecated
        *
-       * See @ref GlossCompress "Compressing
-       * distributed objects" for more
+       * See @ref GlossCompress "Compressing distributed objects" for more
        * information.
        */
       void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
@@ -348,17 +285,10 @@ namespace TrilinosWrappers
 
 
       /**
-       * Returns the state of the
-       * vector, i.e., whether
-       * compress() needs to be
-       * called after an operation
-       * requiring data
-       * exchange. Does only return
-       * non-true values when used in
-       * <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;
@@ -369,34 +299,21 @@ namespace TrilinosWrappers
       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);
 
@@ -446,15 +363,19 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     BlockVector::BlockVector (const std::vector<IndexSet> &parallel_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)
     {
@@ -533,6 +454,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     bool
     BlockVector::has_ghost_elements() const
@@ -545,6 +467,8 @@ namespace TrilinosWrappers
       return ghosted;
     }
 
+
+
     inline
     void
     BlockVector::swap (BlockVector &v)
index 8babcea0032da0b64826169943ee87ce9833e834..2a21d5daa9d8c10508ca50b7dfe66cb4c83dab27 100644 (file)
@@ -452,23 +452,23 @@ namespace TrilinosWrappers
    *
    * When writing into Trilinos matrices from several threads in shared
    * memory, several things must be kept in mind as there is no built-in locks
-   * in this class to prevent data races. Therefore, simultaneous access to
-   * the same matrix row at the same time leads to data races in general and
-   * must be explicitly avoided by the user. However, it is possible to access
-   * <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
@@ -2499,25 +2499,17 @@ namespace TrilinosWrappers
             if (ierr > 0)
               ierr = 0;
           }
-        else if (nonlocal_matrix.get() != 0)
-          {
-            // this is the case when we have explicitly set the off-processor
-            // rows and want to create a separate matrix object for them (to
-            // retain thread-safety)
-            Assert (nonlocal_matrix->RowMap().LID(static_cast<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]));
@@ -2762,8 +2754,7 @@ namespace TrilinosWrappers
 
 
 
-  // inline "simple" functions that are
-  // called frequently and do only involve
+  // inline "simple" functions that are called frequently and do only involve
   // a call to some Trilinos function.
   inline
   SparseMatrix::size_type
index 1c4762e18cd6a0e4b9eb063c96d073ef65bb7558..9e7896532b2066bfd985661cb04ae211693c1ee1 100644 (file)
@@ -186,6 +186,18 @@ namespace TrilinosWrappers
      * 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
@@ -199,15 +211,14 @@ namespace TrilinosWrappers
       typedef dealii::types::global_dof_index size_type;
 
       /**
-       * A variable that indicates whether this vector
-       * supports distributed data storage. If true, then
-       * this vector also needs an appropriate compress()
-       * function that allows communicating recent set or
-       * add operations to individual elements to be communicated
-       * to other processors.
+       * A variable that indicates whether this vector supports distributed
+       * data storage. If true, then this vector also needs an appropriate
+       * compress() function that allows communicating recent set or add
+       * operations to individual elements to be communicated to other
+       * processors.
        *
-       * For the current class, the variable equals
-       * true, since it does support parallel data storage.
+       * For the current class, the variable equals true, since it does
+       * support parallel data storage.
        */
       static const bool supports_distributed_data = true;
 
@@ -216,20 +227,14 @@ namespace TrilinosWrappers
        */
       //@{
       /**
-       * Default constructor that
-       * generates an empty (zero size)
-       * vector. The function
-       * <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);
 
@@ -239,156 +244,88 @@ namespace TrilinosWrappers
       ~Vector ();
 
       /**
-       * Reinit functionality. This
-       * function sets the calling vector
-       * to the dimension and the parallel
-       * distribution of the input vector,
-       * but does not copy the elements in
-       * <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,
@@ -399,35 +336,24 @@ namespace TrilinosWrappers
        */
 //@{
       /**
-       * This constructor takes an
-       * Epetra_Map that already knows
-       * how to distribute the
-       * individual components among
-       * the MPI processors. Since it
-       * also includes information
-       * about the size of the vector,
-       * this is all we need to
-       * generate a parallel vector.
+       * This constructor takes an Epetra_Map that already knows how to
+       * distribute the individual components among the MPI processors. Since
+       * it also includes information about the size of the vector, this is
+       * all we need to generate a parallel vector.
        */
       explicit Vector (const Epetra_Map &parallel_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 &parallel_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>
@@ -435,20 +361,15 @@ namespace TrilinosWrappers
                    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 &parallel_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             &parallel_partitioning,
@@ -459,14 +380,10 @@ namespace TrilinosWrappers
        */
 //@{
       /**
-       * This constructor takes an IndexSet
-       * that defines how to distribute the
-       * individual components among the
-       * MPI processors. Since it also
-       * includes information about the
-       * size of the vector, this is all we
-       * need to generate a %parallel
-       * vector.
+       * This constructor takes an IndexSet that defines how to distribute the
+       * individual components among the MPI processors. Since it also
+       * includes information about the size of the vector, this is all we
+       * need to generate a %parallel vector.
        */
       explicit Vector (const IndexSet &parallel_partitioning,
                        const MPI_Comm &communicator = MPI_COMM_WORLD);
@@ -479,25 +396,18 @@ namespace TrilinosWrappers
               const MPI_Comm &communicator = MPI_COMM_WORLD);
 
       /**
-       * Copy constructor from the
-       * TrilinosWrappers vector
-       * class. Since a vector of this
-       * class does not necessarily need to
-       * be distributed among processes,
-       * the user needs to supply us with
-       * an IndexSet and an MPI
-       * communicator that set the
-       * partitioning details.
+       * Copy constructor from the TrilinosWrappers vector class. Since a
+       * vector of this class does not necessarily need to be distributed
+       * among processes, the user needs to supply us with an IndexSet and an
+       * MPI communicator that set the partitioning details.
        */
       Vector (const IndexSet   &parallel_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               &parallel_partitioning,
@@ -505,17 +415,37 @@ namespace TrilinosWrappers
               const MPI_Comm               &communicator = MPI_COMM_WORLD);
 
       /**
-       * Reinit functionality. This function
-       * destroys the old vector content and
-       * generates a new one based on the
-       * input partitioning.  The flag
-       * <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 &parallel_partitioning,
                    const MPI_Comm &communicator = MPI_COMM_WORLD,
                    const bool      fast = false);
+
+      /**
+       * Reinit functionality. This function destroys the old vector content
+       * and generates a new one based on the input partitioning. In addition
+       * to just specifying one index set as in all the other methods above,
+       * this method allows to supply an additional set of ghost
+       * entries. There are two different versions of a vector that can be
+       * created. If the flag @p vector_writable is set to @p false, the
+       * vector only allows read access to the joint set of @p
+       * parallel_partitioning and @p ghost_entries. The effect of the reinit
+       * method is then equivalent to calling the other reinit method with an
+       * index set containing both the locally owned entries and the ghost
+       * entries.
+       *
+       * If the flag @p vector_writable is set to true, this creates an
+       * alternative storage scheme for ghost elements that allows multiple
+       * threads to write into the vector (for the other reinit methods, only
+       * one thread is allowed to write into the ghost entries at a time).
+       */
+      void reinit (const IndexSet &locally_owned_entries,
+                   const IndexSet &ghost_entries,
+                   const MPI_Comm &communicator = MPI_COMM_WORLD,
+                   const bool      vector_writable = false);
 //@}
     };
 
@@ -639,124 +569,86 @@ namespace TrilinosWrappers
     typedef dealii::types::global_dof_index size_type;
 
     /**
-     * A variable that indicates whether this vector
-     * supports distributed data storage. If true, then
-     * this vector also needs an appropriate compress()
-     * function that allows communicating recent set or
-     * add operations to individual elements to be communicated
-     * to other processors.
+     * A variable that indicates whether this vector supports distributed data
+     * storage. If true, then this vector also needs an appropriate compress()
+     * function that allows communicating recent set or add operations to
+     * individual elements to be communicated to other processors.
      *
-     * For the current class, the variable equals
-     * false, since it does not support parallel data storage.
-     * If you do need parallel data storage, use
-     * TrilinosWrappers::MPI::Vector.
+     * For the current class, the variable equals false, since it does not
+     * support parallel data storage.  If you do need parallel data storage,
+     * use TrilinosWrappers::MPI::Vector.
      */
     static const bool supports_distributed_data = false;
 
     /**
-     * Default constructor that
-     * generates an empty (zero size)
-     * vector. The function
-     * <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,
@@ -764,66 +656,52 @@ namespace TrilinosWrappers
                  const bool        fast = false);
 
     /**
-     * Reinit function. Takes the
-     * information of a Vector and copies
-     * everything to the calling vector,
-     * now also allowing different maps.
+     * Reinit function. Takes the information of a Vector and copies
+     * everything to the calling vector, now also allowing different maps.
      */
     void reinit (const VectorBase &V,
                  const bool        fast = false,
                  const bool        allow_different_maps = false);
 
     /**
-     * Set all components of the
-     * vector to the given number @p
-     * s. Simply pass this down to
-     * the base class, but we still
-     * need to declare this function
-     * to make the example given in
-     * the discussion about making
-     * the constructor explicit work.
+     * Set all components of the vector to the given number @p s. Simply pass
+     * this down to the base class, but we still need to declare this function
+     * to make the example given in the discussion about making the
+     * constructor explicit work.
      */
     Vector &operator = (const TrilinosScalar s);
 
     /**
-     * Sets the left hand argument to
-     * the (parallel) Trilinos
-     * Vector. Equivalent to the @p
-     * reinit function.
+     * Sets the left hand argument to the (parallel) Trilinos
+     * Vector. Equivalent to the @p reinit function.
      */
     Vector &
     operator = (const MPI::Vector &V);
 
     /**
-     * Sets the left hand argument to
-     * the deal.II vector.
+     * Sets the left hand argument to the deal.II vector.
      */
     template <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;
   };
index 66995f19cdc04203431d936c21129d4340f30e0e..649513392adb6d5d6d346509e7133f80225b403e 100644 (file)
@@ -75,29 +75,20 @@ namespace TrilinosWrappers
     typedef dealii::types::global_dof_index size_type;
 
     /**
-     * This class implements a
-     * wrapper for accessing the
-     * Trilinos vector in the same
-     * way as we access deal.II
-     * objects: it is initialized
-     * with a vector and an element
-     * within it, and has a
-     * conversion operator to
-     * extract the scalar value of
-     * this element. It also has a
-     * variety of assignment
-     * operator for writing to this
-     * one element.  @ingroup
-     * TrilinosWrappers
+     * This class implements a wrapper for accessing the Trilinos vector in
+     * the same way as we access deal.II objects: it is initialized with a
+     * vector and an element within it, and has a conversion operator to
+     * extract the scalar value of this element. It also has a variety of
+     * assignment operator for writing to this one element.
+     *
+     * @ingroup TrilinosWrappers
      */
     class VectorReference
     {
     private:
       /**
-       * Constructor. It is made
-       * private so as to only allow
-       * the actual vector class to
-       * create it.
+       * Constructor. It is made private so as to only allow the actual vector
+       * class to create it.
        */
       VectorReference (VectorBase     &vector,
                        const size_type  index);
@@ -105,83 +96,58 @@ namespace TrilinosWrappers
     public:
 
       /**
-       * This looks like a copy
-       * operator, but does something
-       * different than usual. In
-       * particular, it does not copy
-       * the member variables of this
-       * reference. Rather, it
-       * handles the situation where
-       * we have two vectors @p v and
-       * @p w, and assign elements
-       * like in
-       * <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;
 
@@ -209,21 +175,17 @@ namespace TrilinosWrappers
 
     private:
       /**
-       * Point to the vector we are
-       * referencing.
+       * Point to the vector we are referencing.
        */
       VectorBase   &vector;
 
       /**
-       * Index of the referenced element
-       * of the vector.
+       * Index of the referenced element of the vector.
        */
       const size_type  index;
 
       /**
-       * Make the vector class a
-       * friend, so that it can
-       * create objects of the
+       * Make the vector class a friend, so that it can create objects of the
        * present type.
        */
       friend class ::dealii::TrilinosWrappers::VectorBase;
@@ -270,11 +232,8 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Declare some of the standard
-     * types used in all
-     * containers. These types
-     * parallel those in the
-     * <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;
@@ -291,22 +250,15 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * 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);
 
@@ -316,47 +268,30 @@ namespace TrilinosWrappers
     virtual ~VectorBase ();
 
     /**
-     * Release all memory and return
-     * to a state just like after
-     * having called the default
-     * constructor.
+     * Release all memory and return to a state just like after having called
+     * the default constructor.
      */
     void clear ();
 
     /**
-     * Reinit functionality, sets the
-     * dimension and possibly the
-     * parallel partitioning (Epetra_Map)
-     * of the calling vector to the
-     * settings of the input vector.
+     * Reinit functionality, sets the dimension and possibly the parallel
+     * partitioning (Epetra_Map) of the calling vector to the settings of the
+     * input vector.
      */
     void reinit (const VectorBase &v,
                  const bool        fast = false);
 
     /**
-     * Compress the underlying
-     * representation of the Trilinos
-     * object, i.e. flush the buffers
-     * of the vector object if it has
-     * any. This function is
-     * necessary after writing into a
-     * vector element-by-element and
-     * before anything else can be
-     * done on it.
+     * Compress the underlying representation of the Trilinos object,
+     * i.e. flush the buffers of the vector object if it has any. This
+     * function is necessary after writing into a vector element-by-element
+     * and before anything else can be done on it.
      *
-     * The (defaulted) argument can
-     * be used to specify the
-     * compress mode
-     * (<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.
@@ -375,136 +310,87 @@ namespace TrilinosWrappers
     void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
 
     /**
-     * Returns the state of the
-     * vector, i.e., whether
-     * compress() has already been
-     * called after an operation
-     * requiring data exchange.
+     * Returns the state of the vector, i.e., whether compress() has already
+     * been called after an operation requiring data exchange.
      */
     bool is_compressed () const;
 
     /**
-     * Set all components of the
-     * vector to the given number @p
-     * s. Simply pass this down to
-     * the Trilinos Epetra object,
-     * but we still need to declare
-     * this function to make the
-     * example given in the
-     * discussion about making the
-     * constructor explicit work.
+     * Set all components of the vector to the given number @p s. Simply pass
+     * this down to the Trilinos Epetra object, but we still need to declare
+     * this function to make the example given in the discussion about making
+     * the constructor explicit work.
      *
-     * Since the semantics of
-     * assigning a scalar to a vector
-     * are not immediately clear,
-     * this operator should really
-     * only be used if you want to
-     * set the entire vector to
-     * zero. This allows the
-     * intuitive notation
-     * <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
@@ -522,9 +408,8 @@ namespace TrilinosWrappers
     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().
@@ -548,88 +433,65 @@ namespace TrilinosWrappers
     IndexSet locally_owned_elements () const;
 
     /**
-     * Return if the vector contains ghost
-     * elements. This answer is true if there
-     * are ghost elements on at least one
-     * process.
+     * Return if the vector contains ghost elements. This answer is true if
+     * there are ghost elements on at least one process.
      */
     bool has_ghost_elements() const;
 
     /**
-     * Return the scalar (inner)
-     * product of two vectors. The
-     * vectors must have the same
-     * size.
+     * Return the scalar (inner) product of two vectors. The vectors must have
+     * the same size.
      */
     TrilinosScalar operator * (const VectorBase &vec) const;
 
     /**
-     * Return square of the
-     * $l_2$-norm.
+     * Return square of the $l_2$-norm.
      */
     real_type norm_sqr () const;
 
     /**
-     * Mean value of the elements of
-     * this vector.
+     * Mean value of the elements of this vector.
      */
     TrilinosScalar mean_value () const;
 
     /**
-     * Compute the minimal value of
-     * the elements of this vector.
+     * Compute the minimal value of the elements of this vector.
      */
     TrilinosScalar minimal_value () const;
 
     /**
-     * $l_1$-norm of the vector.  The
-     * sum of the absolute values.
+     * $l_1$-norm of the vector.  The sum of the absolute values.
      */
     real_type l1_norm () const;
 
     /**
-     * $l_2$-norm of the vector.  The
-     * square root of the sum of the
-     * squares of the elements.
+     * $l_2$-norm of the vector.  The square root of the sum of the squares of
+     * the elements.
      */
     real_type l2_norm () const;
 
     /**
-     * $l_p$-norm of the vector. The
-     * <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;
     //@}
@@ -641,23 +503,20 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * Provide access to a given
-     * element, both read and write.
+     * Provide access to a given element, both read and write.
      */
     reference
     operator () (const size_type index);
 
     /**
-     * Provide read-only access to an
-     * element. This is equivalent to
-     * the <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().
      */
@@ -665,9 +524,8 @@ namespace TrilinosWrappers
     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().
      */
@@ -675,21 +533,17 @@ namespace TrilinosWrappers
     operator [] (const size_type index) const;
 
     /**
-     * A collective get operation: instead
-     * of getting individual elements of a
-     * vector, this function allows to get
-     * a whole set of elements at once. The
-     * indices of the elements to be read
-     * are stated in the first argument,
-     * the corresponding values are returned in the
-     * second.
+     * A collective get operation: instead of getting individual elements of a
+     * vector, this function allows to get a whole set of elements at
+     * once. The indices of the elements to be read are stated in the first
+     * argument, the corresponding values are returned in the second.
      */
     void extract_subvector_to (const std::vector<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,
@@ -697,22 +551,20 @@ namespace TrilinosWrappers
                                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().
      */
@@ -731,31 +583,23 @@ namespace TrilinosWrappers
     iterator end ();
 
     /**
-     * Return a constant iterator pointing to the element past the end of
-     * the array of the locally owned entries.
+     * Return a constant iterator pointing to the element past the end of the
+     * array of the locally owned entries.
      */
     const_iterator end () const;
 
     /**
-     * A collective set operation:
-     * instead of setting individual
-     * elements of a vector, this
-     * function allows to set a whole
-     * set of elements at once. The
-     * indices of the elements to be
-     * set are stated in the first
-     * argument, the corresponding
-     * values in the second.
+     * A collective set operation: instead of setting individual elements of a
+     * vector, this function allows to set a whole set of elements at
+     * once. The indices of the elements to be set are stated in the first
+     * argument, the corresponding values in the second.
      */
     void set (const std::vector<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);
@@ -768,109 +612,80 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * This collective set operation
-     * is of lower level and can
-     * handle anything else &mdash;
-     * 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 &mdash; 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,
@@ -878,17 +693,14 @@ namespace TrilinosWrappers
               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,
@@ -904,9 +716,8 @@ namespace TrilinosWrappers
                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,
@@ -917,26 +728,20 @@ namespace TrilinosWrappers
                const VectorBase     &X);
 
     /**
-     * Scale each element of this
-     * vector by the corresponding
-     * element in the argument. This
-     * function is mostly meant to
-     * simulate multiplication (and
-     * immediate re-assignment) by a
-     * diagonal scaling matrix.
+     * Scale each element of this vector by the corresponding element in the
+     * argument. This function is mostly meant to simulate multiplication (and
+     * immediate re-assignment) by a diagonal scaling matrix.
      */
     void scale (const VectorBase &scaling_factors);
 
     /**
-     * Assignment <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,
@@ -944,21 +749,14 @@ namespace TrilinosWrappers
               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);
@@ -971,46 +769,34 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos
-     * Epetra_MultiVector class.
+     * Return a const reference to the underlying Trilinos Epetra_MultiVector
+     * class.
      */
     const Epetra_MultiVector &trilinos_vector () const;
 
     /**
-     * Return a (modifyable) reference to
-     * the underlying Trilinos
+     * Return a (modifyable) reference to the underlying Trilinos
      * Epetra_FEVector class.
      */
     Epetra_FEVector &trilinos_vector ();
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the parallel
-     * partitioning of the vector.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the parallel partitioning of the vector.
      */
     const Epetra_Map &vector_partitioner () const;
 
     /**
-     *  Output of vector in
-     *  user-defined format in analogy
-     *  to the dealii::Vector<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,
@@ -1019,42 +805,27 @@ namespace TrilinosWrappers
                 const bool          across     = true) const;
 
     /**
-     * Swap the contents of this
-     * vector and the other vector @p
-     * v. One could do this operation
-     * with a temporary variable and
-     * copying over the data
-     * elements, but this function is
-     * significantly more efficient
-     * since it only swaps the
-     * pointers to the data of the
-     * two vectors and therefore does
-     * not need to allocate temporary
-     * storage and move data
-     * around. Note that the vectors
-     * need to be of the same size
-     * and base on the same map.
+     * Swap the contents of this vector and the other vector @p v. One could
+     * do this operation with a temporary variable and copying over the data
+     * elements, but this function is significantly more efficient since it
+     * only swaps the pointers to the data of the two vectors and therefore
+     * does not need to allocate temporary storage and move data around. Note
+     * that the vectors need to be of the same size and base on the same map.
      *
-     * This function is analog to the
-     * the @p swap function of all C
-     * standard containers. Also,
-     * there is a global function
-     * <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;
@@ -1091,57 +862,45 @@ namespace TrilinosWrappers
 
   private:
     /**
-     * Trilinos doesn't allow to
-     * mix additions to matrix
-     * entries and overwriting them
-     * (to make synchronisation of
-     * parallel computations
-     * simpler). The way we do it
-     * is to, for each access
-     * operation, store whether it
-     * is an insertion or an
-     * addition. If the previous
-     * one was of different type,
-     * then we first have to flush
-     * the Trilinos buffers;
-     * otherwise, we can simply go
-     * on.  Luckily, Trilinos has
-     * an object for this which
-     * does already all the
-     * parallel communications in
-     * such a case, so we simply
-     * use their model, which
-     * stores whether the last
-     * operation was an addition or
-     * an insertion.
+     * Trilinos doesn't allow to mix additions to matrix entries and
+     * overwriting them (to make synchronisation of parallel computations
+     * simpler). The way we do it is to, for each access operation, store
+     * whether it is an insertion or an addition. If the previous one was of
+     * different type, then we first have to flush the Trilinos buffers;
+     * otherwise, we can simply go on.  Luckily, Trilinos has an object for
+     * this which does already all the parallel communications in such a case,
+     * so we simply use their model, which stores whether the last operation
+     * was an addition or an insertion.
      */
     Epetra_CombineMode last_action;
 
     /**
-     * A boolean variable to hold
-     * information on whether the
-     * vector is compressed or not.
+     * A boolean variable to hold information on whether the vector is
+     * compressed or not.
      */
     bool compressed;
 
     /**
-     * Whether this vector has ghost elements. This is true
-     * on all processors even if only one of them has any
-     * ghost elements.
+     * Whether this vector has ghost elements. This is true on all processors
+     * even if only one of them has any ghost elements.
      */
     bool has_ghosts;
 
     /**
-     * An Epetra distibuted vector
-     * type. Requires an existing
-     * Epetra_Map for storing data.
+     * An Epetra distibuted vector type. Requires an existing Epetra_Map for
+     * storing data.
      */
     std_cxx1x::shared_ptr<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;
@@ -1440,58 +1199,6 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  void
-  VectorBase::compress (::dealii::VectorOperation::values given_last_action)
-  {
-    //Select which mode to send to
-    //Trilinos. Note that we use last_action
-    //if available and ignore what the user
-    //tells us to detect wrongly mixed
-    //operations. Typically given_last_action
-    //is only used on machines that do not
-    //execute an operation (because they have
-    //no own cells for example).
-    Epetra_CombineMode mode = last_action;
-    if (last_action == Zero)
-      {
-        if (given_last_action==::dealii::VectorOperation::add)
-          mode = Add;
-        else if (given_last_action==::dealii::VectorOperation::insert)
-          mode = Insert;
-      }
-
-#ifdef DEBUG
-#  ifdef DEAL_II_WITH_MPI
-    // check that every process has decided
-    // to use the same mode. This will
-    // otherwise result in undefined
-    // behaviour in the call to
-    // GlobalAssemble().
-    double double_mode = mode;
-    Utilities::MPI::MinMaxAvg result
-      = Utilities::MPI::min_max_avg (double_mode,
-                                     dynamic_cast<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 ()
@@ -1574,7 +1281,9 @@ namespace TrilinosWrappers
       {
         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),
@@ -1582,8 +1291,11 @@ namespace TrilinosWrappers
             AssertThrow (ierr == 0, ExcTrilinosError(ierr));
             compressed = false;
           }
-        else
-          (*vector)[0][local_row] = values[i];
+        // in set operation, do not use the pre-allocated vector for nonlocal
+        // entries even if it exists. This is to ensure that we really only
+        // set the elements touched by the set() method and not all contained
+        // in the nonlocal entries vector (there is no way to distinguish them
+        // on the receiving processor)
       }
   }
 
@@ -1642,7 +1354,9 @@ namespace TrilinosWrappers
       {
         const size_type row = indices[i];
         const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<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),
@@ -1651,7 +1365,17 @@ namespace TrilinosWrappers
             compressed = false;
           }
         else
-          (*vector)[0][local_row] += values[i];
+          {
+            // use pre-allocated vector for non-local entries if it exists for
+            // addition operation
+            const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<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;
+          }
       }
   }
 
index 0536414ef5f7da6492603d0f392e533685b261c9..d07c32c912e3ca7fc17c47fdd9b6489e301d5467 100644 (file)
@@ -741,6 +741,27 @@ namespace TrilinosWrappers
 
 
 
+  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)
index 4f21ce6136401bba2f25f5d8bba86e8c2fcde249..d1b5a6f5369e2b6428c6cec1f5413721f085cb49 100644 (file)
@@ -146,7 +146,8 @@ namespace TrilinosWrappers
     void
     BlockVector::reinit (const std::vector<IndexSet> &parallel_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);
@@ -161,7 +162,8 @@ namespace TrilinosWrappers
         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();
     }
index 665c7ba1513884dfa8a454e23b2a50f38a0e406a..c486406f7fd7a92c1a65a1d3b08c8869521e3aaf 100644 (file)
@@ -389,6 +389,7 @@ namespace TrilinosWrappers
   void
   SparseMatrix::copy_from (const SparseMatrix &m)
   {
+    nonlocal_matrix.reset();
 
     // check whether we need to update the
     // partitioner or can just copy the data:
@@ -405,6 +406,9 @@ namespace TrilinosWrappers
         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();
   }
 
@@ -475,6 +479,7 @@ namespace TrilinosWrappers
   {
     // release memory before reallocation
     matrix.reset();
+    nonlocal_matrix.reset();
 
     // if we want to exchange data, build a usual Trilinos sparsity pattern
     // and let that handle the exchange. otherwise, manually create a
@@ -582,9 +587,10 @@ namespace TrilinosWrappers
                   (Copy, sparsity_pattern.trilinos_sparsity_pattern(), false));
 
     if (sparsity_pattern.nonlocal_graph.get() != 0)
-      {
-        nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph));
-      }
+      nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph));
+    else
+      nonlocal_matrix.reset ();
+
     compress();
     last_action = Zero;
   }
@@ -598,6 +604,11 @@ namespace TrilinosWrappers
     matrix.reset ();
     matrix.reset (new Epetra_FECrsMatrix
                   (Copy, sparse_matrix.trilinos_sparsity_pattern(), false));
+    if (sparse_matrix.nonlocal_matrix != 0)
+      nonlocal_matrix.reset (new Epetra_CrsMatrix
+                             (Copy, sparse_matrix.nonlocal_matrix->Graph()));
+    else
+      nonlocal_matrix.reset();
 
     compress();
   }
@@ -755,6 +766,7 @@ namespace TrilinosWrappers
 
     const Epetra_CrsGraph *graph = &input_matrix.Graph();
 
+    nonlocal_matrix.reset();
     matrix.reset ();
     matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
 
@@ -801,8 +813,10 @@ namespace TrilinosWrappers
 
     // flush buffers
     int ierr;
-    if (nonlocal_matrix.get() != 0)
+    if (nonlocal_matrix.get() != 0 && mode == Add)
       {
+        // do only export in case of an add() operation, otherwise the owning
+        // processor must have set the correct entry
         nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap());
         Epetra_Export exporter(nonlocal_matrix->RowMap(), matrix->RowMap());
         ierr = matrix->Export(*nonlocal_matrix, exporter, mode);
index 85ca1111fcb3981a5710ef7cf35322dcec16e16e..b406c7300f20e1cc5b971e2e86a230337a80d345 100644 (file)
@@ -544,7 +544,7 @@ namespace TrilinosWrappers
     if (Utilities::MPI::n_mpi_processes(communicator) > 1)
       {
         Epetra_Map nonlocal_map =
-          nonlocal_partitioner.make_trilinos_map(communicator, false);
+          nonlocal_partitioner.make_trilinos_map(communicator, true);
         nonlocal_graph.reset(new Epetra_CrsGraph(Copy, nonlocal_map, 0));
       }
     else
index ea330e88e3a893acd1fb2254a17b809a41cc7ae5..4ff1d598e869a21466cbf2a8dbe67afba2af7ded 100644 (file)
@@ -186,6 +186,8 @@ namespace TrilinosWrappers
     Vector::reinit (const Epetra_Map &input_map,
                     const bool        fast)
     {
+      nonlocal_vector.reset();
+
       if (vector->Map().SameAs(input_map)==false)
         vector.reset (new Epetra_FEVector(input_map));
       else if (fast == false)
@@ -205,6 +207,8 @@ namespace TrilinosWrappers
                     const MPI_Comm &communicator,
                     const bool      fast)
     {
+      nonlocal_vector.reset();
+
       Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator,
                                                                true);
       reinit (map, fast);
@@ -217,13 +221,11 @@ namespace TrilinosWrappers
                     const bool        fast,
                     const bool        allow_different_maps)
     {
-      // In case we do not allow to
-      // have different maps, this
-      // call means that we have to
-      // reset the vector. So clear
-      // the vector, initialize our
-      // map with the map in v, and
-      // generate the vector.
+      nonlocal_vector.reset();
+
+      // In case we do not allow to have different maps, this call means that
+      // we have to reset the vector. So clear the vector, initialize our map
+      // with the map in v, and generate the vector.
       if (allow_different_maps == false)
         {
           if (vector->Map().SameAs(v.vector->Map()) == false)
@@ -250,14 +252,10 @@ namespace TrilinosWrappers
             }
         }
 
-      // Otherwise, we have to check
-      // that the two vectors are
-      // already of the same size,
-      // create an object for the data
-      // exchange and then insert all
-      // the data. The first assertion
-      // is only a check whether the
-      // user knows what she is doing.
+      // Otherwise, we have to check that the two vectors are already of the
+      // same size, create an object for the data exchange and then insert all
+      // the data. The first assertion is only a check whether the user knows
+      // what she is doing.
       else
         {
           Assert (fast == false,
@@ -284,19 +282,16 @@ namespace TrilinosWrappers
     Vector::reinit (const BlockVector &v,
                     const bool         import_data)
     {
-      // In case we do not allow to
-      // have different maps, this
-      // call means that we have to
-      // reset the vector. So clear
-      // the vector, initialize our
-      // map with the map in v, and
-      // generate the vector.
+      nonlocal_vector.reset();
+
+      // In case we do not allow to have different maps, this call means that
+      // we have to reset the vector. So clear the vector, initialize our map
+      // with the map in v, and generate the vector.
       if (v.n_blocks() == 0)
         return;
 
-      // create a vector that holds all the elements
-      // contained in the block vector. need to
-      // manually create an Epetra_Map.
+      // create a vector that holds all the elements contained in the block
+      // vector. need to manually create an Epetra_Map.
       size_type n_elements = 0, added_elements = 0, block_offset = 0;
       for (size_type block=0; block<v.n_blocks(); ++block)
         n_elements += v.block(block).local_size();
@@ -349,26 +344,51 @@ namespace TrilinosWrappers
     }
 
 
-    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
@@ -389,6 +409,9 @@ namespace TrilinosWrappers
           has_ghosts = v.has_ghosts;
         }
 
+      if (v.nonlocal_vector.get() != 0)
+        nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
+      
       return *this;
     }
 
@@ -397,6 +420,8 @@ namespace TrilinosWrappers
     Vector &
     Vector::operator = (const TrilinosWrappers::Vector &v)
     {
+      nonlocal_vector.reset();
+
       Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size()));
 
       Epetra_Import data_exchange (vector->Map(), v.vector->Map());
index 48b0afcd0cef06b832a5e7062513aa82c5c6d615..95a9b1387baa949084d9cc0666fe08f2e8cf7c69 100644 (file)
@@ -21,6 +21,7 @@
 
 #  include <cmath>
 #  include <Epetra_Import.h>
+#  include <Epetra_Export.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -55,13 +56,9 @@ namespace TrilinosWrappers
       Assert (index < vector.size(),
               ExcIndexRange (index, 0, vector.size()));
 
-      // Trilinos allows for vectors
-      // to be referenced by the [] or
-      // () operators but only ()
-      // checks index bounds. We check
-      // these bounds by ourselves, so
-      // we can use []. Note that we
-      // can only get local values.
+      // Trilinos allows for vectors to be referenced by the [] or ()
+      // operators but only () checks index bounds. We check these bounds by
+      // ourselves, so we can use []. Note that we can only get local values.
 
       const TrilinosWrappers::types::int_type local_index =
         vector.vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
@@ -112,9 +109,8 @@ namespace TrilinosWrappers
   void
   VectorBase::clear ()
   {
-    // When we clear the vector,
-    // reset the pointer and generate
-    // an empty vector.
+    // When we clear the vector, reset the pointer and generate an empty
+    // vector.
 #ifdef DEAL_II_WITH_MPI
     Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF));
 #else
@@ -186,18 +182,67 @@ namespace TrilinosWrappers
 
 
 
+  void
+  VectorBase::compress (::dealii::VectorOperation::values given_last_action)
+  {
+    //Select which mode to send to Trilinos. Note that we use last_action if
+    //available and ignore what the user tells us to detect wrongly mixed
+    //operations. Typically given_last_action is only used on machines that do
+    //not execute an operation (because they have no own cells for example).
+    Epetra_CombineMode mode = last_action;
+    if (last_action == Zero)
+      {
+        if (given_last_action==::dealii::VectorOperation::add)
+          mode = Add;
+        else if (given_last_action==::dealii::VectorOperation::insert)
+          mode = Insert;
+      }
+
+#ifdef DEBUG
+#  ifdef DEAL_II_WITH_MPI
+    // check that every process has decided to use the same mode. This will
+    // otherwise result in undefined behaviour in the call to
+    // GlobalAssemble().
+    double double_mode = mode;
+    Utilities::MPI::MinMaxAvg result
+      = Utilities::MPI::min_max_avg (double_mode,
+                                     dynamic_cast<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?
@@ -218,17 +263,13 @@ namespace TrilinosWrappers
   TrilinosScalar
   VectorBase::operator () (const size_type index) const
   {
-    // Extract local indices in
-    // the vector.
+    // Extract local indices in the vector.
     TrilinosWrappers::types::int_type trilinos_i =
       vector->Map().LID(static_cast<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,
index 954afca20aace38595b242e6d5b3600a4f629195..96cabf815a29b23aa75a2939a267d1077025d17a 100644 (file)
@@ -88,9 +88,14 @@ namespace Assembly
   {
     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;
     };
   }
 }
@@ -246,18 +251,27 @@ void LaplaceProblem<dim>::setup_system ()
         (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);
+  }
 }
 
 
@@ -310,9 +324,14 @@ template <int dim>
 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);
 }
 
 
@@ -320,10 +339,10 @@ LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
 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)
@@ -333,11 +352,8 @@ void LaplaceProblem<dim>::assemble_reference ()
         local_assemble(*p, assembly_data, copy_data);
         copy_local_to_global(copy_data);
       }
-  test_matrix.compress(VectorOperation::add);
-  test_rhs.compress(VectorOperation::add);
-
-  reference_matrix.add(1., test_matrix);
-  reference_rhs = test_rhs;
+  reference_matrix.compress(VectorOperation::add);
+  reference_rhs.compress(VectorOperation::add);
 }
 
 
@@ -361,7 +377,7 @@ void LaplaceProblem<dim>::assemble_test ()
                           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);
@@ -401,9 +417,7 @@ void LaplaceProblem<dim>::run ()
     {
       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);
         }
 
index e61846e1b5e0f5567195cb14ecc37b2f5ea5c3a1..ce114164bd66f8918d86a196e6cfcf2c23f8a9c1 100644 (file)
@@ -89,9 +89,14 @@ namespace Assembly
   {
     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;
     };
   }
 }
@@ -247,18 +252,27 @@ void LaplaceProblem<dim>::setup_system ()
         (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);
+  }
 }
 
 
@@ -311,9 +325,14 @@ template <int dim>
 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);
 }
 
 
@@ -321,10 +340,10 @@ LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
 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)
@@ -334,11 +353,8 @@ void LaplaceProblem<dim>::assemble_reference ()
         local_assemble(*p, assembly_data, copy_data);
         copy_local_to_global(copy_data);
       }
-  test_matrix.compress(VectorOperation::add);
-  test_rhs.compress(VectorOperation::add);
-
-  reference_matrix.add(1., test_matrix);
-  reference_rhs = test_rhs;
+  reference_matrix.compress(VectorOperation::add);
+  reference_rhs.compress(VectorOperation::add);
 }
 
 
@@ -362,12 +378,13 @@ void LaplaceProblem<dim>::assemble_test ()
                           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);
index b7745d6e2a2e66a6e5889436b0a5613bec0c60d3..970137a13311fac71687cf5fa6096be7a4b79a79 100644 (file)
@@ -90,9 +90,14 @@ namespace Assembly
   {
     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;
     };
   }
 }
@@ -265,15 +270,22 @@ void LaplaceProblem<dim>::setup_system ()
   relevant_set[1] = relevant_total.get_view(dofs_per_block[0],
                                             dof_handler.n_dofs());
 
-  csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD);
-  DoFTools::make_sparsity_pattern (dof_handler, csp,
-                                   constraints, false);
-  csp.compress();
-  reference_matrix.reinit (csp);
-  test_matrix.reinit (csp);
-
-  reference_rhs.reinit (locally_owned, MPI_COMM_WORLD);
-  test_rhs.reinit (reference_rhs);
+  {
+    csp.reinit(locally_owned, locally_owned, MPI_COMM_WORLD);
+    DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                     constraints, false);
+    csp.compress();
+    reference_matrix.reinit (csp);
+    reference_rhs.reinit (locally_owned, MPI_COMM_WORLD);
+  }
+  {
+    csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD);
+    DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                     constraints, false);
+    csp.compress();
+    test_matrix.reinit (csp);
+    test_rhs.reinit (locally_owned, relevant_set, MPI_COMM_WORLD, true);
+  }
 }
 
 
@@ -328,9 +340,14 @@ template <int dim>
 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);
 }
 
 
@@ -338,10 +355,10 @@ LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
 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)
@@ -351,11 +368,8 @@ void LaplaceProblem<dim>::assemble_reference ()
         local_assemble(*p, assembly_data, copy_data);
         copy_local_to_global(copy_data);
       }
-  test_matrix.compress(VectorOperation::add);
-  test_rhs.compress(VectorOperation::add);
-
-  reference_matrix.add(1., test_matrix);
-  reference_rhs = test_rhs;
+  reference_matrix.compress(VectorOperation::add);
+  reference_rhs.compress(VectorOperation::add);
 }
 
 
@@ -379,7 +393,7 @@ void LaplaceProblem<dim>::assemble_test ()
                           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);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.