]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement new reinit method to Trilinos sparsity pattern that makes writing into...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 23 Dec 2013 22:14:25 +0000 (22:14 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 23 Dec 2013 22:14:25 +0000 (22:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@32106 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_sparse_matrix.h
deal.II/include/deal.II/lac/trilinos_sparsity_pattern.h
deal.II/source/lac/block_sparsity_pattern.cc
deal.II/source/lac/trilinos_sparse_matrix.cc
deal.II/source/lac/trilinos_sparsity_pattern.cc
tests/trilinos/assemble_matrix_parallel_01.cc [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_01.output [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_02.cc [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_03.cc [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_04.cc [new file with mode: 0644]
tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output [new file with mode: 0644]

index be811da6063c3b227d723e40d4bc38f1b5f06fa2..e2c2969b4b329b582d1e670752ee1f0f4bbd05fb 100644 (file)
@@ -84,6 +84,17 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+
+  <li> New: A new reinit() method has been introduced to
+  TrilinosWrappers::SparsityPattern that takes all rows that are possibly
+  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).
+  <br>
+  (Martin Kronbichler, 2013/12/23)
+  </li>
+
   <li> New: The TableBase::fill function has become more powerful in that
   it now doesn't just take pointers to initializing elements but can deal
   with arbitrary input iterators. It now also takes a flag that denotes the
index 392eaa4fca35564736737c7a2f14496a060cedf7..e22741f97c971659a2df15d6a874524b901a8ea7 100644 (file)
@@ -1016,113 +1016,100 @@ namespace TrilinosWrappers
   public:
 
     /**
-     * Initialize the matrix empty,
-     * that is with no memory
-     * allocated. This is useful if
-     * you want such objects as
-     * member variables in other
-     * classes. You can make the
-     * structure usable by calling
-     * the reinit() function.
+     * Initialize the matrix empty, that is with no memory allocated. This is
+     * useful if you want such objects as member variables in other
+     * classes. You can make the structure usable by calling the reinit()
+     * function.
      */
     BlockSparsityPattern ();
 
     /**
-     * Initialize the matrix with the
-     * given number of block rows and
-     * columns. The blocks themselves
-     * are still empty, and you have
-     * to call collect_sizes() after
-     * you assign them sizes.
+     * Initialize the matrix with the given number of block rows and
+     * columns. The blocks themselves are still empty, and you have to call
+     * collect_sizes() after you assign them sizes.
      */
     BlockSparsityPattern (const size_type n_rows,
                           const size_type n_columns);
 
     /**
-     * Initialize the pattern with
-     * two BlockIndices for the block
-     * structures of matrix rows and
-     * columns. This function is
-     * equivalent to calling the
-     * previous constructor with the
-     * length of the two index vector
-     * and then entering the index
-     * values.
+     * Initialize the pattern with two BlockIndices for the block structures
+     * of matrix rows and columns. This function is equivalent to calling the
+     * previous constructor with the length of the two index vector and then
+     * entering the index values.
      */
     BlockSparsityPattern (const std::vector<size_type> &row_block_sizes,
                           const std::vector<size_type> &col_block_sizes);
 
     /**
-     * Initialize the pattern with an array
-     * Epetra_Map that specifies both rows
-     * and columns of the matrix (so the
-     * final matrix will be a square
-     * matrix), where the Epetra_Map
-     * specifies the parallel distribution
-     * of the degrees of freedom on the
-     * individual block.  This function is
-     * equivalent to calling the second
-     * constructor with the length of the
-     * mapping vector and then entering the
-     * index values.
+     * Initialize the pattern with an array Epetra_Map that specifies both
+     * rows and columns of the matrix (so the final matrix will be a square
+     * matrix), where the Epetra_Map specifies the parallel distribution of
+     * the degrees of freedom on the individual block.  This function is
+     * equivalent to calling the second constructor with the length of the
+     * mapping vector and then entering the index values.
      */
     BlockSparsityPattern (const std::vector<Epetra_Map> &parallel_partitioning);
 
     /**
-     * Initialize the pattern with an array
-     * of index sets that specifies both rows
-     * and columns of the matrix (so the
-     * final matrix will be a square 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.
+     * Initialize the pattern with an array of index sets that specifies both
+     * rows and columns of the matrix (so the final matrix will be a square
+     * 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.
      */
     BlockSparsityPattern (const std::vector<IndexSet> &parallel_partitioning,
                           const MPI_Comm &communicator = MPI_COMM_WORLD);
 
     /**
-     * Resize the matrix to a tensor
-     * product of matrices with
-     * dimensions defined by the
-     * arguments.
+     * Resize the matrix to a tensor product of matrices with dimensions
+     * defined by the arguments.
      *
-     * The matrix will have as many
-     * block rows and columns as
-     * there are entries in the two
-     * arguments. The block at
-     * position (<i>i,j</i>) will
-     * have the dimensions
-     * <tt>row_block_sizes[i]</tt>
-     * times <tt>col_block_sizes[j]</tt>.
+     * The matrix will have as many block rows and columns as there are
+     * entries in the two arguments. The block at position (<i>i,j</i>) will
+     * have the dimensions <tt>row_block_sizes[i]</tt> times
+     * <tt>col_block_sizes[j]</tt>.
      */
     void reinit (const std::vector<size_type> &row_block_sizes,
                  const std::vector<size_type> &col_block_sizes);
 
     /**
-     * Resize the matrix to a square tensor
-     * product of matrices with parallel
-     * distribution according to the
-     * specifications in the array of
+     * Resize the matrix to a square tensor product of matrices with parallel
+     * distribution according to the specifications in the array of
      * Epetra_Maps.
      */
     void reinit (const std::vector<Epetra_Map> &parallel_partitioning);
 
     /**
-     * Resize the matrix to a square tensor
-     * product of matrices. See the
-     * constructor that takes a vector of
-     * IndexSets for details.
+     * Resize the matrix to a square tensor product of matrices. See the
+     * constructor that takes a vector of IndexSets for details.
      */
     void reinit (const std::vector<IndexSet> &parallel_partitioning,
                  const MPI_Comm              &communicator = MPI_COMM_WORLD);
 
+    /**
+     * Resize the matrix to a rectangular block matrices. This method allows
+     * rows and columns to be different, both in the outer block structure and
+     * within the blocks.
+     */
+    void reinit (const std::vector<IndexSet> &row_parallel_partitioning,
+                 const std::vector<IndexSet> &column_parallel_partitioning,
+                 const MPI_Comm              &communicator = MPI_COMM_WORLD);
+
+    /**
+     * Resize the matrix to a rectangular block matrices that furthermore
+     * explicitly specify the writable rows in each of the blocks. This method
+     * 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.
+     */
+    void reinit (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);
 
     /**
-     * Allow the use of the reinit
-     * functions of the base class as
-     * well.
+     * Allow the use of the reinit functions of the base class as well.
      */
     using BlockSparsityPatternBase<SparsityPattern>::reinit;
   };
index ed72eb0060b676f8321a94f489ab69697ec81c9b..8babcea0032da0b64826169943ee87ce9833e834 100644 (file)
@@ -84,17 +84,11 @@ namespace TrilinosWrappers
                     << " are stored locally and can be accessed.");
 
     /**
-     * Handling of indices for both
-     * constant and non constant
-     * Accessor objects
-     *
-     * For a regular
-     * dealii::SparseMatrix, we would
-     * use an accessor for the sparsity
-     * pattern. For Trilinos matrices,
-     * this does not seem so simple,
-     * therefore, we write a little
-     * base class here.
+     * Handling of indices for both constant and non constant Accessor objects
+     *
+     * For a regular dealii::SparseMatrix, we would use an accessor for the
+     * sparsity pattern. For Trilinos matrices, this does not seem so simple,
+     * therefore, we write a little base class here.
      *
      * @author Guido Kanschat
      * @date 2012
@@ -115,34 +109,26 @@ namespace TrilinosWrappers
                     const size_type  index);
 
       /**
-       * Row number of the element
-       * represented by this object.
+       * Row number of the element represented by this object.
        */
       size_type row() const;
 
       /**
-       * Index in row of the element
-       * represented by this object.
+       * Index in row of the element represented by this object.
        */
       size_type index() const;
 
       /**
-       * Column number of the element
-       * represented by this object.
+       * Column number of the element represented by this object.
        */
       size_type column() const;
 
     protected:
       /**
-       * Pointer to the matrix
-       * object. This object should
-       * be handled as a const
-       * pointer or non-const by the
-       * appropriate derived
-       * classes. In order to be able
-       * to implement both, it is not
-       * const here, so handle with
-       * care!
+       * Pointer to the matrix object. This object should be handled as a
+       * const pointer or non-const by the appropriate derived classes. In
+       * order to be able to implement both, it is not const here, so handle
+       * with care!
        */
       mutable SparseMatrix *matrix;
       /**
@@ -156,45 +142,28 @@ namespace TrilinosWrappers
       size_type a_index;
 
       /**
-       * Discard the old row caches
-       * (they may still be used by
-       * other accessors) and
-       * generate new ones for the
-       * row pointed to presently by
+       * Discard the old row caches (they may still be used by other
+       * accessors) and generate new ones for the row pointed to presently by
        * this accessor.
        */
       void visit_present_row ();
 
       /**
-       * Cache where we store the
-       * column indices of the
-       * present row. This is
-       * necessary, since Trilinos
-       * makes access to the elements
-       * of its matrices rather hard,
-       * and it is much more
-       * efficient to copy all column
-       * entries of a row once when
-       * we enter it than repeatedly
-       * asking Trilinos for
-       * individual ones. This also
-       * makes some sense since it is
-       * likely that we will access
-       * them sequentially anyway.
+       * Cache where we store the column indices of the present row. This is
+       * necessary, since Trilinos makes access to the elements of its
+       * matrices rather hard, and it is much more efficient to copy all
+       * column entries of a row once when we enter it than repeatedly asking
+       * Trilinos for individual ones. This also makes some sense since it is
+       * likely that we will access them sequentially anyway.
        *
-       * In order to make copying of
-       * iterators/accessor of
-       * acceptable performance, we
-       * keep a shared pointer to
-       * these entries so that more
-       * than one accessor can access
-       * this data if necessary.
+       * In order to make copying of iterators/accessor of acceptable
+       * performance, we keep a shared pointer to these entries so that more
+       * than one accessor can access this data if necessary.
        */
       std_cxx1x::shared_ptr<std::vector<size_type> > colnum_cache;
 
       /**
-       * Cache for the values
-       * of this row.
+       * Cache for the values of this row.
        */
       std_cxx1x::shared_ptr<std::vector<TrilinosScalar> > value_cache;
     };
@@ -231,25 +200,21 @@ namespace TrilinosWrappers
     {
     public:
       /**
-       * Typedef for the type (including
-       * constness) of the matrix to be
-       * used here.
+       * Typedef for the type (including constness) of the matrix to be used
+       * here.
        */
       typedef const SparseMatrix MatrixType;
 
       /**
-       * Constructor. Since we use
-       * accessors only for read
-       * access, a const matrix
-       * pointer is sufficient.
+       * Constructor. Since we use accessors only for read access, a const
+       * matrix pointer is sufficient.
        */
       Accessor (MatrixType *matrix,
                 const size_type  row,
                 const size_type  index);
 
       /**
-       * Copy constructor to get from a
-       * const or non-const accessor to a const
+       * Copy constructor to get from a const or non-const accessor to a const
        * accessor.
        */
       template <bool Other>
@@ -262,8 +227,7 @@ namespace TrilinosWrappers
 
     private:
       /**
-       * Make iterator class a
-       * friend.
+       * Make iterator class a friend.
        */
       template <bool> friend class Iterator;
     };
@@ -283,66 +247,53 @@ namespace TrilinosWrappers
         Reference (const Accessor<false> &accessor);
 
         /**
-         * Conversion operator to the
-         * data type of the matrix.
+         * Conversion operator to the data type of the matrix.
          */
         operator TrilinosScalar () const;
 
         /**
-         * Set the element of the matrix
-         * we presently point to to @p n.
+         * Set the element of the matrix we presently point to to @p n.
          */
         const Reference &operator = (const TrilinosScalar n) const;
 
         /**
-         * Add @p n to the element of the
-         * matrix we presently point to.
+         * Add @p n to the element of the matrix we presently point to.
          */
         const Reference &operator += (const TrilinosScalar n) const;
 
         /**
-         * Subtract @p n from the element
-         * of the matrix we presently
-         * point to.
+         * Subtract @p n from the element of the matrix we presently point to.
          */
         const Reference &operator -= (const TrilinosScalar n) const;
 
         /**
-         * Multiply the element of the
-         * matrix we presently point to
-         * by @p n.
+         * Multiply the element of the matrix we presently point to by @p n.
          */
         const Reference &operator *= (const TrilinosScalar n) const;
 
         /**
-         * Divide the element of the
-         * matrix we presently point to
-         * by @p n.
+         * Divide the element of the matrix we presently point to by @p n.
          */
         const Reference &operator /= (const TrilinosScalar n) const;
 
       private:
         /**
-         * Pointer to the accessor that
-         * denotes which element we
-         * presently point to.
+         * Pointer to the accessor that denotes which element we presently
+         * point to.
          */
         Accessor &accessor;
       };
 
     public:
       /**
-       * Typedef for the type (including
-       * constness) of the matrix to be
-       * used here.
+       * Typedef for the type (including constness) of the matrix to be used
+       * here.
        */
       typedef SparseMatrix MatrixType;
 
       /**
-       * Constructor. Since we use
-       * accessors only for read
-       * access, a const matrix
-       * pointer is sufficient.
+       * Constructor. Since we use accessors only for read access, a const
+       * matrix pointer is sufficient.
        */
       Accessor (MatrixType *matrix,
                 const size_type  row,
@@ -355,13 +306,11 @@ namespace TrilinosWrappers
 
     private:
       /**
-       * Make iterator class a
-       * friend.
+       * Make iterator class a friend.
        */
       template <bool> friend class Iterator;
       /**
-       * Make Reference object a
-       * friend.
+       * Make Reference object a friend.
        */
       friend class Reference;
     };
@@ -390,25 +339,21 @@ namespace TrilinosWrappers
       typedef dealii::types::global_dof_index size_type;
 
       /**
-       * Typedef for the matrix type
-       * (including constness) we are to
-       * operate on.
+       * Typedef for the matrix type (including constness) we are to operate
+       * on.
        */
       typedef typename Accessor<Constness>::MatrixType MatrixType;
 
       /**
-       * Constructor. Create an
-       * iterator into the matrix @p
-       * matrix for the given row and
-       * the index within it.
+       * Constructor. Create an iterator into the matrix @p matrix for the
+       * given row and the index within it.
        */
       Iterator (MatrixType *matrix,
                 const size_type  row,
                 const size_type  index);
 
       /**
-       * Copy constructor with
-       * optional change of constness.
+       * Copy constructor with optional change of constness.
        */
       template <bool Other>
       Iterator(const Iterator<Other> &other);
@@ -434,9 +379,8 @@ namespace TrilinosWrappers
       const Accessor<Constness> *operator-> () const;
 
       /**
-       * Comparison. True, if both
-       * iterators point to the same
-       * matrix position.
+       * Comparison. True, if both iterators point to the same matrix
+       * position.
        */
       bool operator == (const Iterator<Constness> &) const;
 
@@ -446,11 +390,8 @@ namespace TrilinosWrappers
       bool operator != (const Iterator<Constness> &) const;
 
       /**
-       * Comparison operator. Result
-       * is true if either the first
-       * row number is smaller or if
-       * the row numbers are equal
-       * and the first index is
+       * Comparison operator. Result is true if either the first row number is
+       * smaller or if the row numbers are equal and the first index is
        * smaller.
        */
       bool operator < (const Iterator<Constness> &) const;
@@ -471,8 +412,7 @@ namespace TrilinosWrappers
 
     private:
       /**
-       * Store an object of the
-       * accessor class.
+       * Store an object of the accessor class.
        */
       Accessor<Constness> accessor;
 
@@ -508,6 +448,29 @@ namespace TrilinosWrappers
    * unused elements. Trilinos allows to continue with assembling the
    * matrix after calls to these functions, though.
    *
+   * <h3>Thread safety of Trilinos matrices</h3>
+   *
+   * 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:
+   * <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.
+   * </ul>
+   *
    * @ingroup TrilinosWrappers
    * @ingroup Matrix1
    * @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009
@@ -521,43 +484,33 @@ namespace TrilinosWrappers
     typedef dealii::types::global_dof_index size_type;
 
     /**
-     * A structure that describes
-     * some of the traits of this
-     * class in terms of its run-time
-     * behavior. Some other classes
-     * (such as the block matrix
-     * classes) that take one or
-     * other of the matrix classes as
-     * its template parameters can
-     * tune their behavior based on
-     * the variables in this class.
+     * A structure that describes some of the traits of this class in terms of
+     * its run-time behavior. Some other classes (such as the block matrix
+     * classes) that take one or other of the matrix classes as its template
+     * parameters can tune their behavior based on the variables in this
+     * class.
      */
     struct Traits
     {
       /**
-       * It is safe to elide additions
-       * of zeros to individual
-       * elements of this matrix.
+       * It is safe to elide additions of zeros to individual elements of this
+       * matrix.
        */
       static const bool zero_addition_can_be_elided = true;
     };
 
     /**
-     * Declare a typedef for the
-     * iterator class.
+     * Declare a typedef for the iterator class.
      */
     typedef SparseMatrixIterators::Iterator<false> iterator;
 
     /**
-     * Declare a typedef for the
-     * const iterator class.
+     * Declare a typedef for the const iterator class.
      */
     typedef SparseMatrixIterators::Iterator<true> const_iterator;
 
     /**
-     * Declare a typedef in analogy
-     * to all the other container
-     * classes.
+     * Declare a typedef in analogy to all the other container classes.
      */
     typedef TrilinosScalar value_type;
 
@@ -566,18 +519,15 @@ namespace TrilinosWrappers
      */
 //@{
     /**
-     * Default constructor. Generates
-     * an empty (zero-size) matrix.
+     * Default constructor. Generates an empty (zero-size) matrix.
      */
     SparseMatrix ();
 
     /**
-     * Generate a matrix that is completely
-     * stored locally, having #m rows and
+     * Generate a matrix that is completely stored locally, having #m rows and
      * #n columns.
      *
-     * The number of columns entries per
-     * row is specified as the maximum
+     * The number of columns entries per row is specified as the maximum
      * number of entries argument.
      */
     SparseMatrix (const size_type  m,
@@ -585,118 +535,85 @@ namespace TrilinosWrappers
                   const unsigned int  n_max_entries_per_row);
 
     /**
-     * Generate a matrix that is completely
-     * stored locally, having #m rows and
+     * Generate a matrix that is completely stored locally, having #m rows and
      * #n columns.
      *
-     * The vector
-     * <tt>n_entries_per_row</tt>
-     * specifies the number of entries in
-     * each row.
+     * The vector <tt>n_entries_per_row</tt> specifies the number of entries
+     * in each row.
      */
     SparseMatrix (const size_type                  m,
                   const size_type                  n,
                   const std::vector<unsigned int> &n_entries_per_row);
 
     /**
-     * Generate a matrix from a Trilinos
-     * sparsity pattern object.
+     * Generate a matrix from a Trilinos sparsity pattern object.
      */
     SparseMatrix (const SparsityPattern &InputSparsityPattern);
 
     /**
-     * Copy constructor. Sets the
-     * calling matrix to be the same
-     * as the input matrix, i.e.,
-     * using the same sparsity
-     * pattern and entries.
+     * Copy constructor. Sets the calling matrix to be the same as the input
+     * matrix, i.e., using the same sparsity pattern and entries.
      */
     SparseMatrix (const SparseMatrix &InputMatrix);
 
     /**
-     * Destructor. Made virtual so
-     * that one can use pointers to
-     * this class.
+     * Destructor. Made virtual so that one can use pointers to this class.
      */
     virtual ~SparseMatrix ();
 
     /**
-     * This function initializes the
-     * Trilinos matrix with a deal.II
-     * sparsity pattern, i.e. it makes
-     * the Trilinos Epetra matrix know
-     * the position of nonzero entries
-     * according to the sparsity
-     * pattern. This function is meant
-     * for use in serial programs, where
-     * there is no need to specify how
-     * the matrix is going to be
-     * distributed among different
-     * processors. This function works in
-     * %parallel, too, but it is
-     * recommended to manually specify
-     * the %parallel partioning of the
-     * matrix using an Epetra_Map. When
-     * run in %parallel, it is currently
-     * necessary that each processor
-     * holds the sparsity_pattern
-     * structure because each processor
+     * This function initializes the Trilinos matrix with a deal.II sparsity
+     * pattern, i.e. it makes the Trilinos Epetra matrix know the position of
+     * nonzero entries according to the sparsity pattern. This function is
+     * meant for use in serial programs, where there is no need to specify how
+     * the matrix is going to be distributed among different processors. This
+     * function works in %parallel, too, but it is recommended to manually
+     * specify the %parallel partioning of the matrix using an
+     * Epetra_Map. When run in %parallel, it is currently necessary that each
+     * processor holds the sparsity_pattern structure because each processor
      * sets its rows.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template<typename SparsityType>
     void reinit (const SparsityType &sparsity_pattern);
 
     /**
-     * This function reinitializes the
-     * Trilinos sparse matrix from a
-     * (possibly distributed) Trilinos
-     * sparsity pattern.
+     * This function reinitializes the Trilinos sparse matrix from a (possibly
+     * distributed) Trilinos sparsity pattern.
+     *
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * If you want to write to the matrix from several threads and use MPI,
+     * you need to use this reinit method with a sparsity pattern that has
+     * been created with explicitly stating writeable rows. In all other
+     * cases, you cannot mix MPI with multithreaded writing into the matrix.
      */
     void reinit (const SparsityPattern &sparsity_pattern);
 
     /**
-     * This function copies the content
-     * in <tt>sparse_matrix</tt> to the
+     * This function copies the content in <tt>sparse_matrix</tt> to the
      * calling matrix.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     void reinit (const SparseMatrix &sparse_matrix);
 
     /**
-     * This function initializes the
-     * Trilinos matrix using the deal.II
-     * sparse matrix and the entries
-     * stored therein. It uses a
-     * threshold to copy only elements
-     * with modulus larger than the
-     * threshold (so zeros in the deal.II
-     * matrix can be filtered away).
-     *
-     * The optional parameter
-     * <tt>copy_values</tt> decides
-     * whether only the sparsity
-     * structure of the input matrix
-     * should be used or the matrix
+     * This function initializes the Trilinos matrix using the deal.II sparse
+     * matrix and the entries stored therein. It uses a threshold to copy only
+     * elements with modulus larger than the threshold (so zeros in the
+     * deal.II matrix can be filtered away).
+     *
+     * The optional parameter <tt>copy_values</tt> decides whether only the
+     * sparsity structure of the input matrix should be used or the matrix
      * entries should be copied, too.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * deadlock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a deadlock.
      *
      * @note If a different sparsity pattern is given in the last argument
      * (i.e., one that differs from the one used in the sparse matrix given
@@ -712,11 +629,9 @@ namespace TrilinosWrappers
                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 
     /**
-     * This reinit function takes as
-     * input a Trilinos Epetra_CrsMatrix
-     * and copies its sparsity
-     * pattern. If so requested, even the
-     * content (values) will be copied.
+     * This reinit function takes as input a Trilinos Epetra_CrsMatrix and
+     * copies its sparsity pattern. If so requested, even the content (values)
+     * will be copied.
      */
     void reinit (const Epetra_CrsMatrix &input_matrix,
                  const bool              copy_values = true);
@@ -726,150 +641,87 @@ namespace TrilinosWrappers
      */
 //@{
     /**
-     * Constructor using an Epetra_Map to
-     * describe the %parallel
-     * partitioning. The parameter @p
-     * n_max_entries_per_row sets the
-     * number of nonzero entries in each
-     * row that will be allocated. Note
-     * that this number does not need to
-     * be exact, and it is even allowed
-     * that the actual matrix structure
-     * has more nonzero entries than
-     * specified in the
-     * constructor. However it is still
-     * advantageous to provide good
-     * estimates here since this will
-     * considerably increase the
-     * performance of the matrix
-     * setup. However, there is no effect
-     * in the performance of
-     * matrix-vector products, since
-     * Trilinos reorganizes the matrix
-     * memory prior to use (in the
-     * compress() step).
+     * Constructor using an Epetra_Map to describe the %parallel
+     * partitioning. The parameter @p n_max_entries_per_row sets the number of
+     * nonzero entries in each row that will be allocated. Note that this
+     * number does not need to be exact, and it is even allowed that the
+     * actual matrix structure has more nonzero entries than specified in the
+     * constructor. However it is still advantageous to provide good estimates
+     * here since this will considerably increase the performance of the
+     * matrix setup. However, there is no effect in the performance of
+     * matrix-vector products, since Trilinos reorganizes the matrix memory
+     * prior to use (in the compress() step).
      */
     SparseMatrix (const Epetra_Map  &parallel_partitioning,
                   const size_type    n_max_entries_per_row = 0);
 
     /**
-     * Same as before, but now set a
-     * value of nonzeros for each matrix
-     * row. Since we know the number of
-     * elements in the matrix exactly in
-     * this case, we can already allocate
-     * the right amount of memory, which
-     * makes the creation process
-     * including the insertion of nonzero
-     * elements by the respective
-     * SparseMatrix::reinit call
-     * considerably faster.
+     * Same as before, but now set a value of nonzeros for each matrix
+     * row. Since we know the number of elements in the matrix exactly in this
+     * case, we can already allocate the right amount of memory, which makes
+     * the creation process including the insertion of nonzero elements by the
+     * respective SparseMatrix::reinit call considerably faster.
      */
     SparseMatrix (const Epetra_Map                &parallel_partitioning,
                   const std::vector<unsigned int> &n_entries_per_row);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different Epetra maps for rows and
-     * columns. This interface is meant
-     * to be used for generating
-     * rectangular matrices, where one
-     * map describes the %parallel
-     * partitioning of the dofs
-     * associated with the matrix rows
-     * and the other one the partitioning
-     * of dofs in the matrix
-     * columns. Note that there is no
-     * real parallelism along the columns
-     * &ndash; the processor that owns a
-     * certain row always owns all the
-     * column elements, no matter how far
-     * they might be spread out. The
-     * second Epetra_Map is only used to
-     * specify the number of columns and
-     * for internal arrangements when
-     * doing matrix-vector products with
-     * vectors based on that column map.
-     *
-     * The integer input @p
-     * n_max_entries_per_row defines the
-     * number of columns entries per row
-     * that will be allocated.
+     * This constructor is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular matrices, where one map describes
+     * the %parallel partitioning of the dofs associated with the matrix rows
+     * and the other one the partitioning of dofs in the matrix columns. Note
+     * that there is no real parallelism along the columns &ndash; the
+     * processor that owns a certain row always owns all the column elements,
+     * no matter how far they might be spread out. The second Epetra_Map is
+     * only used to specify the number of columns and for internal
+     * arrangements when doing matrix-vector products with vectors based on
+     * that column map.
+     *
+     * The integer input @p n_max_entries_per_row defines the number of
+     * columns entries per row that will be allocated.
      */
     SparseMatrix (const Epetra_Map &row_parallel_partitioning,
                   const Epetra_Map &col_parallel_partitioning,
                   const size_type   n_max_entries_per_row = 0);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different Epetra maps for rows and
-     * columns. This interface is meant
-     * to be used for generating
-     * rectangular matrices, where one
-     * map specifies the %parallel
-     * distribution of degrees of freedom
-     * associated with matrix rows and
-     * the second one specifies the
-     * %parallel distribution the dofs
-     * associated with columns in the
-     * matrix. The second map also
-     * provides information for the
-     * internal arrangement in matrix
-     * vector products (i.e., the
-     * distribution of vector this matrix
-     * is to be multiplied with), but is
-     * not used for the distribution of
-     * the columns &ndash; rather, all
-     * column elements of a row are
-     * stored on the same processor in
-     * any case. The vector
-     * <tt>n_entries_per_row</tt>
-     * specifies the number of entries in
-     * each row of the newly generated
-     * matrix.
+     * This constructor is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular matrices, where one map specifies
+     * the %parallel distribution of degrees of freedom associated with matrix
+     * rows and the second one specifies the %parallel distribution the dofs
+     * associated with columns in the matrix. The second map also provides
+     * information for the internal arrangement in matrix vector products
+     * (i.e., the distribution of vector this matrix is to be multiplied
+     * with), but is not used for the distribution of the columns &ndash;
+     * rather, all column elements of a row are stored on the same processor
+     * in any case. The vector <tt>n_entries_per_row</tt> specifies the number
+     * of entries in each row of the newly generated matrix.
      */
     SparseMatrix (const Epetra_Map                &row_parallel_partitioning,
                   const Epetra_Map                &col_parallel_partitioning,
                   const std::vector<unsigned int> &n_entries_per_row);
 
     /**
-     * This function is initializes the
-     * Trilinos Epetra matrix according to
-     * the specified sparsity_pattern, and
-     * also reassigns the matrix rows to
-     * different processes according to a
-     * user-supplied Epetra map. In
-     * programs following the style of the
-     * tutorial programs, this function
-     * (and the respective call for a
-     * rectangular matrix) are the natural
-     * way to initialize the matrix size,
-     * its distribution among the MPI
-     * processes (if run in %parallel) as
-     * well as the locatoin of non-zero
-     * elements. Trilinos stores the
-     * sparsity pattern internally, so it
-     * won't be needed any more after this
-     * call, in contrast to the deal.II own
-     * object. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
-     * CompressedSimpleSparsityPattern. If
-     * the flag is not set, each processor
-     * just sets the elements in the
-     * sparsity pattern that belong to its
-     * rows.
-     *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This function is initializes the Trilinos Epetra matrix according to
+     * the specified sparsity_pattern, and also reassigns the matrix rows to
+     * different processes according to a user-supplied Epetra map. In
+     * programs following the style of the tutorial programs, this function
+     * (and the respective call for a rectangular matrix) are the natural way
+     * to initialize the matrix size, its distribution among the MPI processes
+     * (if run in %parallel) as well as the locatoin of non-zero
+     * elements. Trilinos stores the sparsity pattern internally, so it won't
+     * be needed any more after this call, in contrast to the deal.II own
+     * object. The optional argument @p exchange_data can be used for
+     * reinitialization with a sparsity pattern that is not fully
+     * constructed. This feature is only implemented for input sparsity
+     * patterns of type CompressedSimpleSparsityPattern. If the flag is not
+     * set, each processor just sets the elements in the sparsity pattern that
+     * belong to its rows.
+     *
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template<typename SparsityType>
     void reinit (const Epetra_Map    &parallel_partitioning,
@@ -877,25 +729,16 @@ namespace TrilinosWrappers
                  const bool          exchange_data = false);
 
     /**
-     * This function is similar to the
-     * other initialization function
-     * above, but now also reassigns the
-     * matrix rows and columns according
-     * to two user-supplied Epetra maps.
-     * To be used for rectangular
-     * matrices. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
+     * This function is similar to the other initialization function above,
+     * but now also reassigns the matrix rows and columns according to two
+     * user-supplied Epetra maps.  To be used for rectangular matrices. The
+     * optional argument @p exchange_data can be used for reinitialization
+     * with a sparsity pattern that is not fully constructed. This feature is
+     * only implemented for input sparsity patterns of type
      * CompressedSimpleSparsityPattern.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template<typename SparsityType>
     void reinit (const Epetra_Map    &row_parallel_partitioning,
@@ -904,32 +747,20 @@ namespace TrilinosWrappers
                  const bool          exchange_data = false);
 
     /**
-     * This function initializes the
-     * Trilinos matrix using the deal.II
-     * sparse matrix and the entries
-     * stored therein. It uses a
-     * threshold to copy only elements
-     * with modulus larger than the
-     * threshold (so zeros in the deal.II
-     * matrix can be filtered away). In
-     * contrast to the other reinit
-     * function with deal.II sparse
-     * matrix argument, this function
-     * takes a %parallel partitioning
-     * specified by the user instead of
-     * internally generating it.
-     *
-     * The optional parameter
-     * <tt>copy_values</tt> decides
-     * whether only the sparsity
-     * structure of the input matrix
-     * should be used or the matrix
+     * This function initializes the Trilinos matrix using the deal.II sparse
+     * matrix and the entries stored therein. It uses a threshold to copy only
+     * elements with modulus larger than the threshold (so zeros in the
+     * deal.II matrix can be filtered away). In contrast to the other reinit
+     * function with deal.II sparse matrix argument, this function takes a
+     * %parallel partitioning specified by the user instead of internally
+     * generating it.
+     *
+     * The optional parameter <tt>copy_values</tt> decides whether only the
+     * sparsity structure of the input matrix should be used or the matrix
      * entries should be copied, too.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template <typename number>
     void reinit (const Epetra_Map                     &parallel_partitioning,
@@ -939,25 +770,17 @@ namespace TrilinosWrappers
                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 
     /**
-     * This function is similar to the
-     * other initialization function with
-     * deal.II sparse matrix input above,
-     * but now takes Epetra maps for both
-     * the rows and the columns of the
-     * matrix. Chosen for rectangular
+     * This function is similar to the other initialization function with
+     * deal.II sparse matrix input above, but now takes Epetra maps for both
+     * the rows and the columns of the matrix. Chosen for rectangular
      * matrices.
      *
-     * The optional parameter
-     * <tt>copy_values</tt> decides
-     * whether only the sparsity
-     * structure of the input matrix
-     * should be used or the matrix
+     * The optional parameter <tt>copy_values</tt> decides whether only the
+     * sparsity structure of the input matrix should be used or the matrix
      * entries should be copied, too.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template <typename number>
     void reinit (const Epetra_Map                      &row_parallel_partitioning,
@@ -972,77 +795,45 @@ namespace TrilinosWrappers
      */
 //@{
     /**
-     * Constructor using an IndexSet and
-     * an MPI communicator to describe
-     * the %parallel partitioning. The
-     * parameter @p n_max_entries_per_row
-     * sets the number of nonzero entries
-     * in each row that will be
-     * allocated. Note that this number
-     * does not need to be exact, and it
-     * is even allowed that the actual
-     * matrix structure has more nonzero
-     * entries than specified in the
-     * constructor. However it is still
-     * advantageous to provide good
-     * estimates here since this will
-     * considerably increase the
-     * performance of the matrix
-     * setup. However, there is no effect
-     * in the performance of
-     * matrix-vector products, since
-     * Trilinos reorganizes the matrix
-     * memory prior to use (in the
-     * compress() step).
+     * Constructor using an IndexSet and an MPI communicator to describe the
+     * %parallel partitioning. The parameter @p n_max_entries_per_row sets the
+     * number of nonzero entries in each row that will be allocated. Note that
+     * this number does not need to be exact, and it is even allowed that the
+     * actual matrix structure has more nonzero entries than specified in the
+     * constructor. However it is still advantageous to provide good estimates
+     * here since this will considerably increase the performance of the
+     * matrix setup. However, there is no effect in the performance of
+     * matrix-vector products, since Trilinos reorganizes the matrix memory
+     * prior to use (in the compress() step).
      */
     SparseMatrix (const IndexSet    &parallel_partitioning,
                   const MPI_Comm    &communicator = MPI_COMM_WORLD,
                   const unsigned int n_max_entries_per_row = 0);
 
     /**
-     * Same as before, but now set the
-     * number of nonzeros in each matrix
-     * row separately. Since we know the
-     * number of elements in the matrix
-     * exactly in this case, we can
-     * already allocate the right amount
-     * of memory, which makes the
-     * creation process including the
-     * insertion of nonzero elements by
-     * the respective
-     * SparseMatrix::reinit call
-     * considerably faster.
+     * Same as before, but now set the number of nonzeros in each matrix row
+     * separately. Since we know the number of elements in the matrix exactly
+     * in this case, we can already allocate the right amount of memory, which
+     * makes the creation process including the insertion of nonzero elements
+     * by the respective SparseMatrix::reinit call considerably faster.
      */
     SparseMatrix (const IndexSet                  &parallel_partitioning,
                   const MPI_Comm                  &communicator,
                   const std::vector<unsigned int> &n_entries_per_row);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different IndexSet partitions for
-     * row and columns. This interface is
-     * meant to be used for generating
-     * rectangular matrices, where the
-     * first index set describes the
-     * %parallel partitioning of the
-     * degrees of freedom associated with
-     * the matrix rows and the second one
-     * the partitioning of the matrix
-     * columns. The second index set
-     * specifies the partitioning of the
-     * vectors this matrix is to be
-     * multiplied with, not the
-     * distribution of the elements that
-     * actually appear in the matrix.
-     *
-     * The parameter @p
-     * n_max_entries_per_row defines how
-     * much memory will be allocated for
-     * each row. This number does not
-     * need to be accurate, as the
-     * structure is reorganized in the
-     * compress() call.
+     * This constructor is similar to the one above, but it now takes two
+     * different IndexSet partitions for row and columns. This interface is
+     * meant to be used for generating rectangular matrices, where the first
+     * index set describes the %parallel partitioning of the degrees of
+     * freedom associated with the matrix rows and the second one the
+     * partitioning of the matrix columns. The second index set specifies the
+     * partitioning of the vectors this matrix is to be multiplied with, not
+     * the distribution of the elements that actually appear in the matrix.
+     *
+     * The parameter @p n_max_entries_per_row defines how much memory will be
+     * allocated for each row. This number does not need to be accurate, as
+     * the structure is reorganized in the compress() call.
      */
     SparseMatrix (const IndexSet  &row_parallel_partitioning,
                   const IndexSet  &col_parallel_partitioning,
@@ -1050,33 +841,18 @@ namespace TrilinosWrappers
                   const size_type  n_max_entries_per_row = 0);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different Epetra maps for rows and
-     * columns. This interface is meant
-     * to be used for generating
-     * rectangular matrices, where one
-     * map specifies the %parallel
-     * distribution of degrees of freedom
-     * associated with matrix rows and
-     * the second one specifies the
-     * %parallel distribution the dofs
-     * associated with columns in the
-     * matrix. The second map also
-     * provides information for the
-     * internal arrangement in matrix
-     * vector products (i.e., the
-     * distribution of vector this matrix
-     * is to be multiplied with), but is
-     * not used for the distribution of
-     * the columns &ndash; rather, all
-     * column elements of a row are
-     * stored on the same processor in
-     * any case. The vector
-     * <tt>n_entries_per_row</tt>
-     * specifies the number of entries in
-     * each row of the newly generated
-     * matrix.
+     * This constructor is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular matrices, where one map specifies
+     * the %parallel distribution of degrees of freedom associated with matrix
+     * rows and the second one specifies the %parallel distribution the dofs
+     * associated with columns in the matrix. The second map also provides
+     * information for the internal arrangement in matrix vector products
+     * (i.e., the distribution of vector this matrix is to be multiplied
+     * with), but is not used for the distribution of the columns &ndash;
+     * rather, all column elements of a row are stored on the same processor
+     * in any case. The vector <tt>n_entries_per_row</tt> specifies the number
+     * of entries in each row of the newly generated matrix.
      */
     SparseMatrix (const IndexSet                  &row_parallel_partitioning,
                   const IndexSet                  &col_parallel_partitioning,
@@ -1084,43 +860,24 @@ namespace TrilinosWrappers
                   const std::vector<unsigned int> &n_entries_per_row);
 
     /**
-     * This function is initializes the
-     * Trilinos Epetra matrix according
-     * to the specified sparsity_pattern,
-     * and also reassigns the matrix rows
-     * to different processes according
-     * to a user-supplied index set and
-     * %parallel communicator. In
-     * programs following the style of
-     * the tutorial programs, this
-     * function (and the respective call
-     * for a rectangular matrix) are the
-     * natural way to initialize the
-     * matrix size, its distribution
-     * among the MPI processes (if run in
-     * %parallel) as well as the locatoin
-     * of non-zero elements. Trilinos
-     * stores the sparsity pattern
-     * internally, so it won't be needed
-     * any more after this call, in
-     * contrast to the deal.II own
-     * object. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
-     * CompressedSimpleSparsityPattern. If
-     * the flag is not set, each
-     * processor just sets the elements
-     * in the sparsity pattern that
+     * This function is initializes the Trilinos Epetra matrix according to
+     * the specified sparsity_pattern, and also reassigns the matrix rows to
+     * different processes according to a user-supplied index set and
+     * %parallel communicator. In programs following the style of the tutorial
+     * programs, this function (and the respective call for a rectangular
+     * matrix) are the natural way to initialize the matrix size, its
+     * distribution among the MPI processes (if run in %parallel) as well as
+     * the locatoin of non-zero elements. Trilinos stores the sparsity pattern
+     * internally, so it won't be needed any more after this call, in contrast
+     * to the deal.II own object. The optional argument @p exchange_data can
+     * be used for reinitialization with a sparsity pattern that is not fully
+     * constructed. This feature is only implemented for input sparsity
+     * patterns of type CompressedSimpleSparsityPattern. If the flag is not
+     * set, each processor just sets the elements in the sparsity pattern that
      * belong to its rows.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template<typename SparsityType>
     void reinit (const IndexSet      &parallel_partitioning,
@@ -1129,25 +886,16 @@ namespace TrilinosWrappers
                  const bool           exchange_data = false);
 
     /**
-     * This function is similar to the
-     * other initialization function
-     * above, but now also reassigns the
-     * matrix rows and columns according
-     * to two user-supplied index sets.
-     * To be used for rectangular
-     * matrices. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
+     * This function is similar to the other initialization function above,
+     * but now also reassigns the matrix rows and columns according to two
+     * user-supplied index sets.  To be used for rectangular matrices. The
+     * optional argument @p exchange_data can be used for reinitialization
+     * with a sparsity pattern that is not fully constructed. This feature is
+     * only implemented for input sparsity patterns of type
      * CompressedSimpleSparsityPattern.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template<typename SparsityType>
     void reinit (const IndexSet      &row_parallel_partitioning,
@@ -1157,32 +905,20 @@ namespace TrilinosWrappers
                  const bool           exchange_data = false);
 
     /**
-     * This function initializes the
-     * Trilinos matrix using the deal.II
-     * sparse matrix and the entries
-     * stored therein. It uses a
-     * threshold to copy only elements
-     * with modulus larger than the
-     * threshold (so zeros in the deal.II
-     * matrix can be filtered away). In
-     * contrast to the other reinit
-     * function with deal.II sparse
-     * matrix argument, this function
-     * takes a %parallel partitioning
-     * specified by the user instead of
-     * internally generating it.
-     *
-     * The optional parameter
-     * <tt>copy_values</tt> decides
-     * whether only the sparsity
-     * structure of the input matrix
-     * should be used or the matrix
+     * This function initializes the Trilinos matrix using the deal.II sparse
+     * matrix and the entries stored therein. It uses a threshold to copy only
+     * elements with modulus larger than the threshold (so zeros in the
+     * deal.II matrix can be filtered away). In contrast to the other reinit
+     * function with deal.II sparse matrix argument, this function takes a
+     * %parallel partitioning specified by the user instead of internally
+     * generating it.
+     *
+     * The optional parameter <tt>copy_values</tt> decides whether only the
+     * sparsity structure of the input matrix should be used or the matrix
      * entries should be copied, too.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template <typename number>
     void reinit (const IndexSet                       &parallel_partitioning,
@@ -1193,25 +929,17 @@ namespace TrilinosWrappers
                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 
     /**
-     * This function is similar to the
-     * other initialization function with
-     * deal.II sparse matrix input above,
-     * but now takes index sets for both
-     * the rows and the columns of the
-     * matrix. Chosen for rectangular
+     * This function is similar to the other initialization function with
+     * deal.II sparse matrix input above, but now takes index sets for both
+     * the rows and the columns of the matrix. Chosen for rectangular
      * matrices.
      *
-     * The optional parameter
-     * <tt>copy_values</tt> decides
-     * whether only the sparsity
-     * structure of the input matrix
-     * should be used or the matrix
+     * The optional parameter <tt>copy_values</tt> decides whether only the
+     * sparsity structure of the input matrix should be used or the matrix
      * entries should be copied, too.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     template <typename number>
     void reinit (const IndexSet                        &row_parallel_partitioning,
@@ -1228,89 +956,64 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Return the number of rows in
-     * this matrix.
+     * Return the number of rows in this matrix.
      */
     size_type m () const;
 
     /**
-     * Return the number of columns
-     * in this matrix.
+     * Return the number of columns in this matrix.
      */
     size_type n () const;
 
     /**
-     * Return the local dimension
-     * of the matrix, i.e. the
-     * number of rows stored on the
-     * present MPI process. For
-     * sequential matrices, this
-     * number is the same as m(),
-     * but for %parallel matrices it
-     * may be smaller.
+     * Return the local dimension of the matrix, i.e. the number of rows
+     * stored on the present MPI process. For sequential matrices, this number
+     * is the same as m(), but for %parallel matrices 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().
      */
     unsigned int local_size () const;
 
     /**
-     * Return a pair of indices
-     * indicating which rows of
-     * this matrix are stored
-     * locally. The first number is
-     * the index of the first row
-     * stored, the second the index
-     * of the one past the last one
-     * that is stored locally. If
-     * this is a sequential matrix,
-     * then the result will be the
-     * pair (0,m()), otherwise it
-     * will be a pair (i,i+n),
-     * where
+     * Return a pair of indices indicating which rows of this matrix are
+     * stored locally. The first number is the index of the first row stored,
+     * the second the index of the one past the last one that is stored
+     * locally. If this is a sequential matrix, then the result will be the
+     * pair (0,m()), otherwise it will be a pair (i,i+n), where
      * <tt>n=local_size()</tt>.
      */
     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().
      */
     bool in_local_range (const size_type index) const;
 
     /**
-     * Return the number of nonzero
-     * elements of this matrix.
+     * Return the number of nonzero elements of this matrix.
      */
     size_type n_nonzero_elements () const;
 
     /**
-     * Number of entries in a
-     * specific row.
+     * Number of entries in a specific row.
      */
     unsigned int row_length (const size_type row) const;
 
     /**
-     * Returns the state of the matrix,
-     * i.e., whether compress() needs to
-     * be called after an operation
-     * requiring data exchange. A call to
-     * compress() is also needed when the
-     * method set() has been called (even
-     * when working in serial).
+     * Returns the state of the matrix, i.e., whether compress() needs to be
+     * called after an operation requiring data exchange. A call to compress()
+     * is also needed when the method set() has been called (even when working
+     * in serial).
      */
     bool is_compressed () const;
 
     /**
-     * Determine an estimate for the memory
-     * consumption (in bytes) of this
-     * object. Note that only the memory
-     * reserved on the current processor is
-     * returned in case this is called in
-     * an MPI-based program.
+     * Determine an estimate for the memory consumption (in bytes) of this
+     * object. Note that only the memory reserved on the current processor is
+     * returned in case this is called in an MPI-based program.
      */
     size_type memory_consumption () const;
 
@@ -1321,75 +1024,48 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * This operator assigns a scalar to
-     * a matrix. Since this does usually
-     * not make much sense (should we set
-     * all matrix entries to this value?
-     * Only the nonzero entries of the
-     * sparsity pattern?), this operation
-     * is only allowed if the actual
-     * value to be assigned is zero. This
-     * operator only exists to allow for
-     * the obvious notation
-     * <tt>matrix=0</tt>, which sets all
-     * elements of the matrix to zero,
-     * but keeps the sparsity pattern
+     * This operator assigns a scalar to a matrix. Since this does usually not
+     * make much sense (should we set all matrix entries to this value?  Only
+     * the nonzero entries of the sparsity pattern?), this operation is only
+     * allowed if the actual value to be assigned is zero. This operator only
+     * exists to allow for the obvious notation <tt>matrix=0</tt>, which sets
+     * all elements of the matrix to zero, but keeps the sparsity pattern
      * previously used.
      */
     SparseMatrix &
     operator = (const double d);
 
     /**
-     * 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.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     void clear ();
 
     /**
      * This command does two things:
      * <ul>
-     * <li> If the matrix was initialized
-     * without a sparsity pattern,
-     * elements have been added manually
-     * using the set() command. When this
-     * process is completed, a call to
-     * compress() reorganizes the
-     * internal data structures (aparsity
-     * pattern) so that a fast access to
-     * data is possible in matrix-vector
-     * products.
-     * <li> If the matrix structure has
-     * already been fixed (either by
-     * initialization with a sparsity
-     * pattern or by calling compress()
-     * during the setup phase), this
-     * command does the %parallel
-     * exchange of data. This is
-     * necessary when we perform assembly
-     * on more than one (MPI) process,
-     * because then some non-local row
-     * data will accumulate on nodes that
-     * belong to the current's processor
-     * element, but are actually held by
-     * another. This command is usually
-     * called after all elements have
-     * been traversed.
+     * <li> If the matrix was initialized without a sparsity pattern, elements
+     * have been added manually using the set() command. When this process is
+     * completed, a call to compress() reorganizes the internal data
+     * structures (aparsity pattern) so that a fast access to data is possible
+     * in matrix-vector products.
+     * <li> If the matrix structure has already been fixed (either by
+     * initialization with a sparsity pattern or by calling compress() during
+     * the setup phase), this command does the %parallel exchange of
+     * data. This is necessary when we perform assembly on more than one (MPI)
+     * process, because then some non-local row data will accumulate on nodes
+     * that belong to the current's processor element, but are actually held
+     * by another. This command is usually called after all elements have been
+     * traversed.
      * </ul>
      *
-     * In both cases, this function
-     * compresses the data structures and
-     * allows the resulting matrix to be
-     * used in all other operations like
-     * matrix-vector products. This is a
-     * collective operation, i.e., it
-     * needs to be run on all processors
-     * when used in %parallel.
+     * In both cases, this function compresses the data structures and allows
+     * the resulting matrix to be used in all other operations like
+     * matrix-vector products. This is a collective operation, i.e., it needs
+     * to be run on all processors when used in %parallel.
      *
      * See @ref GlossCompress "Compressing distributed objects"
      * for more information.
@@ -1402,69 +1078,42 @@ namespace TrilinosWrappers
     void compress () DEAL_II_DEPRECATED;
 
     /**
-     * Set the element (<i>i,j</i>)
-     * to @p value.
-     *
-     * This function is able to insert new
-     * elements into the matrix as long as
-     * compress() has not been called, so
-     * the sparsity pattern will be
-     * extended. When compress() is called
-     * for the first time, then this is no
-     * longer possible and an insertion of
-     * elements at positions which have not
-     * been initialized will throw an
-     * exception. Note that in case
-     * elements need to be inserted, it is
-     * mandatory that elements are inserted
-     * only once. Otherwise, the elements
-     * will actually be added in the end
-     * (since it is not possible to
-     * efficiently find values to the same
-     * entry before compress() has been
-     * called). In the case that an element
-     * is set more than once, initialize
-     * the matrix with a sparsity pattern
-     * first.
+     * Set the element (<i>i,j</i>) to @p value.
+     *
+     * This function is able to insert new elements into the matrix as long as
+     * compress() has not been called, so the sparsity pattern will be
+     * extended. When compress() is called for the first time, then this is no
+     * longer possible and an insertion of elements at positions which have
+     * not been initialized will throw an exception. Note that in case
+     * elements need to be inserted, it is mandatory that elements are
+     * inserted only once. Otherwise, the elements will actually be added in
+     * the end (since it is not possible to efficiently find values to the
+     * same entry before compress() has been called). In the case that an
+     * element is set more than once, initialize the matrix with a sparsity
+     * pattern first.
      */
     void set (const size_type i,
               const size_type j,
               const TrilinosScalar value);
 
     /**
-     * Set all elements given in a
-     * FullMatrix<double> into the sparse
-     * matrix locations given by
-     * <tt>indices</tt>. In other words,
-     * this function writes the elements
-     * in <tt>full_matrix</tt> into the
-     * calling matrix, using the
-     * local-to-global indexing specified
-     * by <tt>indices</tt> for both the
-     * rows and the columns of the
-     * matrix. This function assumes a
-     * quadratic sparse matrix and a
-     * quadratic full_matrix, the usual
+     * Set all elements given in a FullMatrix<double> into the sparse matrix
+     * locations given by <tt>indices</tt>. In other words, this function
+     * writes the elements in <tt>full_matrix</tt> into the calling matrix,
+     * using the local-to-global indexing specified by <tt>indices</tt> for
+     * both the rows and the columns of the matrix. This function assumes a
+     * quadratic sparse matrix and a quadratic full_matrix, the usual
      * situation in FE calculations.
      *
-     * This function is able to insert
-     * new elements into the matrix as
-     * long as compress() has not been
-     * called, so the sparsity pattern
-     * will be extended. When compress()
-     * is called for the first time, then
-     * this is no longer possible and an
-     * insertion of elements at positions
-     * which have not been initialized
-     * will throw an exception.
+     * This function is able to insert new elements into the matrix as long as
+     * compress() has not been called, so the sparsity pattern will be
+     * extended. When compress() is called for the first time, then this is no
+     * longer possible and an insertion of elements at positions which have
+     * not been initialized will throw an exception.
      *
-     * The optional parameter
-     * <tt>elide_zero_values</tt> can be
-     * used to specify whether zero
-     * values should be inserted anyway
-     * or they should be filtered
-     * away. The default value is
-     * <tt>false</tt>, i.e., even zero
+     * The optional parameter <tt>elide_zero_values</tt> can be used to
+     * specify whether zero values should be inserted anyway or they should be
+     * filtered away. The default value is <tt>false</tt>, i.e., even zero
      * values are inserted/replaced.
      */
     void set (const std::vector<size_type>     &indices,
@@ -1472,11 +1121,9 @@ namespace TrilinosWrappers
               const bool                        elide_zero_values = false);
 
     /**
-     * Same function as before, but now
-     * including the possibility to use
-     * rectangular full_matrices and
-     * different local-to-global indexing
-     * on rows and columns, respectively.
+     * Same function as before, but now including the possibility to use
+     * rectangular full_matrices and different local-to-global indexing on
+     * rows and columns, respectively.
      */
     void set (const std::vector<size_type>     &row_indices,
               const std::vector<size_type>     &col_indices,
@@ -1484,30 +1131,18 @@ namespace TrilinosWrappers
               const bool                        elide_zero_values = false);
 
     /**
-     * Set several elements in the
-     * specified row of the matrix with
-     * column indices as given by
-     * <tt>col_indices</tt> to the
-     * respective value.
-     *
-     * This function is able to insert
-     * new elements into the matrix as
-     * long as compress() has not been
-     * called, so the sparsity pattern
-     * will be extended. When compress()
-     * is called for the first time, then
-     * this is no longer possible and an
-     * insertion of elements at positions
-     * which have not been initialized
-     * will throw an exception.
+     * Set several elements in the specified row of the matrix with column
+     * indices as given by <tt>col_indices</tt> to the respective value.
      *
-     * The optional parameter
-     * <tt>elide_zero_values</tt> can be
-     * used to specify whether zero
-     * values should be inserted anyway
-     * or they should be filtered
-     * away. The default value is
-     * <tt>false</tt>, i.e., even zero
+     * This function is able to insert new elements into the matrix as long as
+     * compress() has not been called, so the sparsity pattern will be
+     * extended. When compress() is called for the first time, then this is no
+     * longer possible and an insertion of elements at positions which have
+     * not been initialized will throw an exception.
+     *
+     * The optional parameter <tt>elide_zero_values</tt> can be used to
+     * specify whether zero values should be inserted anyway or they should be
+     * filtered away. The default value is <tt>false</tt>, i.e., even zero
      * values are inserted/replaced.
      */
     void set (const size_type                    row,
@@ -1516,30 +1151,18 @@ namespace TrilinosWrappers
               const bool                         elide_zero_values = false);
 
     /**
-     * Set several elements to values
-     * given by <tt>values</tt> in a
-     * given row in columns given by
-     * col_indices into the sparse
-     * matrix.
+     * Set several elements to values given by <tt>values</tt> in a given row
+     * in columns given by col_indices into the sparse matrix.
      *
-     * This function is able to insert
-     * new elements into the matrix as
-     * long as compress() has not been
-     * called, so the sparsity pattern
-     * will be extended. When compress()
-     * is called for the first time, then
-     * this is no longer possible and an
-     * insertion of elements at positions
-     * which have not been initialized
-     * will throw an exception.
+     * This function is able to insert new elements into the matrix as long as
+     * compress() has not been called, so the sparsity pattern will be
+     * extended. When compress() is called for the first time, then this is no
+     * longer possible and an insertion of elements at positions which have
+     * not been initialized will throw an exception.
      *
-     * The optional parameter
-     * <tt>elide_zero_values</tt> can be
-     * used to specify whether zero
-     * values should be inserted anyway
-     * or they should be filtered
-     * away. The default value is
-     * <tt>false</tt>, i.e., even zero
+     * The optional parameter <tt>elide_zero_values</tt> can be used to
+     * specify whether zero values should be inserted anyway or they should be
+     * filtered away. The default value is <tt>false</tt>, i.e., even zero
      * values are inserted/replaced.
      */
     void set (const size_type       row,
@@ -1549,70 +1172,44 @@ namespace TrilinosWrappers
               const bool            elide_zero_values = false);
 
     /**
-     * Add @p value to the element
-     * (<i>i,j</i>).
+     * Add @p value to the element (<i>i,j</i>).
      *
-     * Just as the respective call in
-     * deal.II SparseMatrix<Number>
-     * class (but in contrast to the
-     * situation for PETSc based
-     * matrices), this function
-     * throws an exception if an
-     * entry does not exist in the
-     * sparsity pattern. Moreover, if
-     * <tt>value</tt> is not a finite
-     * number an exception is thrown.
+     * Just as the respective call in deal.II SparseMatrix<Number> class (but
+     * in contrast to the situation for PETSc based matrices), this function
+     * throws an exception if an entry does not exist in the sparsity
+     * pattern. Moreover, if <tt>value</tt> is not a finite number an
+     * exception is thrown.
      */
     void add (const size_type      i,
               const size_type      j,
               const TrilinosScalar value);
 
     /**
-     * Add all elements given in a
-     * FullMatrix<double> into sparse
-     * matrix locations given by
-     * <tt>indices</tt>. In other words,
-     * this function adds the elements in
-     * <tt>full_matrix</tt> to the
-     * respective entries in calling
-     * matrix, using the local-to-global
-     * indexing specified by
-     * <tt>indices</tt> for both the rows
-     * and the columns of the
-     * matrix. This function assumes a
-     * quadratic sparse matrix and a
-     * quadratic full_matrix, the usual
-     * situation in FE calculations.
+     * Add all elements given in a FullMatrix<double> into sparse matrix
+     * locations given by <tt>indices</tt>. In other words, this function adds
+     * the elements in <tt>full_matrix</tt> to the respective entries in
+     * calling matrix, using the local-to-global indexing specified by
+     * <tt>indices</tt> for both the rows and the columns of the matrix. This
+     * function assumes a quadratic sparse matrix and a quadratic full_matrix,
+     * the usual situation in FE calculations.
+     *
+     * Just as the respective call in deal.II SparseMatrix<Number> class (but
+     * in contrast to the situation for PETSc based matrices), this function
+     * throws an exception if an entry does not exist in the sparsity pattern.
      *
-     * Just as the respective call in
-     * deal.II SparseMatrix<Number>
-     * class (but in contrast to the
-     * situation for PETSc based
-     * matrices), this function
-     * throws an exception if an
-     * entry does not exist in the
-     * sparsity pattern.
-     *
-     * The optional parameter
-     * <tt>elide_zero_values</tt> can be
-     * used to specify whether zero
-     * values should be added anyway or
-     * these should be filtered away and
-     * only non-zero data is added. The
-     * default value is <tt>true</tt>,
-     * i.e., zero values won't be added
-     * into the matrix.
+     * The optional parameter <tt>elide_zero_values</tt> can be used to
+     * specify whether zero values should be added anyway or these should be
+     * filtered away and only non-zero data is added. The default value is
+     * <tt>true</tt>, i.e., zero values won't be added into the matrix.
      */
     void add (const std::vector<size_type>  &indices,
               const FullMatrix<TrilinosScalar> &full_matrix,
               const bool                        elide_zero_values = true);
 
     /**
-     * Same function as before, but now
-     * including the possibility to use
-     * rectangular full_matrices and
-     * different local-to-global indexing
-     * on rows and columns, respectively.
+     * Same function as before, but now including the possibility to use
+     * rectangular full_matrices and different local-to-global indexing on
+     * rows and columns, respectively.
      */
     void add (const std::vector<size_type>     &row_indices,
               const std::vector<size_type>     &col_indices,
@@ -1620,30 +1217,17 @@ namespace TrilinosWrappers
               const bool                        elide_zero_values = true);
 
     /**
-     * Set several elements in the
-     * specified row of the matrix with
-     * column indices as given by
-     * <tt>col_indices</tt> to the
-     * respective value.
-     *
-     * Just as the respective call in
-     * deal.II SparseMatrix<Number>
-     * class (but in contrast to the
-     * situation for PETSc based
-     * matrices), this function
-     * throws an exception if an
-     * entry does not exist in the
-     * sparsity pattern.
-     *
-     * The optional parameter
-     * <tt>elide_zero_values</tt> can be
-     * used to specify whether zero
-     * values should be added anyway or
-     * these should be filtered away and
-     * only non-zero data is added. The
-     * default value is <tt>true</tt>,
-     * i.e., zero values won't be added
-     * into the matrix.
+     * Set several elements in the specified row of the matrix with column
+     * indices as given by <tt>col_indices</tt> to the respective value.
+     *
+     * Just as the respective call in deal.II SparseMatrix<Number> class (but
+     * in contrast to the situation for PETSc based matrices), this function
+     * throws an exception if an entry does not exist in the sparsity pattern.
+     *
+     * The optional parameter <tt>elide_zero_values</tt> can be used to
+     * specify whether zero values should be added anyway or these should be
+     * filtered away and only non-zero data is added. The default value is
+     * <tt>true</tt>, i.e., zero values won't be added into the matrix.
      */
     void add (const size_type                    row,
               const std::vector<size_type>      &col_indices,
@@ -1651,29 +1235,17 @@ namespace TrilinosWrappers
               const bool                         elide_zero_values = true);
 
     /**
-     * Add an array of values given by
-     * <tt>values</tt> in the given
-     * global matrix row at columns
-     * specified by col_indices in the
-     * sparse matrix.
-     *
-     * Just as the respective call in
-     * deal.II SparseMatrix<Number> class
-     * (but in contrast to the situation
-     * for PETSc based matrices), this
-     * function throws an exception if an
-     * entry does not exist in the
-     * sparsity pattern.
-     *
-     * The optional parameter
-     * <tt>elide_zero_values</tt> can be
-     * used to specify whether zero
-     * values should be added anyway or
-     * these should be filtered away and
-     * only non-zero data is added. The
-     * default value is <tt>true</tt>,
-     * i.e., zero values won't be added
-     * into the matrix.
+     * Add an array of values given by <tt>values</tt> in the given global
+     * matrix row at columns specified by col_indices in the sparse matrix.
+     *
+     * Just as the respective call in deal.II SparseMatrix<Number> class (but
+     * in contrast to the situation for PETSc based matrices), this function
+     * throws an exception if an entry does not exist in the sparsity pattern.
+     *
+     * The optional parameter <tt>elide_zero_values</tt> can be used to
+     * specify whether zero values should be added anyway or these should be
+     * filtered away and only non-zero data is added. The default value is
+     * <tt>true</tt>, i.e., zero values won't be added into the matrix.
      */
     void add (const size_type       row,
               const size_type       n_cols,
@@ -1683,14 +1255,12 @@ namespace TrilinosWrappers
               const bool            col_indices_are_sorted = false);
 
     /**
-     * Multiply the entire matrix
-     * by a fixed factor.
+     * Multiply the entire matrix by a fixed factor.
      */
     SparseMatrix &operator *= (const TrilinosScalar factor);
 
     /**
-     * Divide the entire matrix by
-     * a fixed factor.
+     * Divide the entire matrix by a fixed factor.
      */
     SparseMatrix &operator /= (const TrilinosScalar factor);
 
@@ -1701,84 +1271,48 @@ namespace TrilinosWrappers
     void copy_from (const SparseMatrix &source);
 
     /**
-     * Add <tt>matrix</tt> scaled by
-     * <tt>factor</tt> to this matrix,
-     * i.e. the matrix
-     * <tt>factor*matrix</tt> is added to
-     * <tt>this</tt>. If the sparsity
-     * pattern of the calling matrix does
-     * not contain all the elements in
-     * the sparsity pattern of the input
-     * matrix, this function will throw
-     * an exception.
+     * Add <tt>matrix</tt> scaled by <tt>factor</tt> to this matrix, i.e. the
+     * matrix <tt>factor*matrix</tt> is added to <tt>this</tt>. If the
+     * sparsity pattern of the calling matrix does not contain all the
+     * elements in the sparsity pattern of the input matrix, this function
+     * will throw an exception.
      */
     void add (const TrilinosScalar  factor,
               const SparseMatrix   &matrix);
 
     /**
-     * Remove all elements from
-     * this <tt>row</tt> by setting
-     * them to zero. The function
-     * does not modify the number
-     * of allocated nonzero
-     * entries, it only sets some
-     * entries to zero. It may drop
-     * them from the sparsity
-     * pattern, though (but retains
-     * the allocated memory in case
-     * new entries are again added
-     * later). Note that this is a
-     * global operation, so this
-     * needs to be done on all MPI
-     * processes.
-     *
-     * This operation is used in
-     * eliminating constraints
-     * (e.g. due to hanging nodes)
-     * and makes sure that we can
-     * write this modification to
-     * the matrix without having to
-     * read entries (such as the
-     * locations of non-zero
-     * elements) from it &mdash;
-     * without this operation,
-     * removing constraints on
-     * %parallel matrices is a
-     * rather complicated
-     * procedure.
-     *
-     * The second parameter can be
-     * used to set the diagonal
-     * entry of this row to a value
-     * different from zero. The
-     * default is to set it to
-     * zero.
+     * Remove all elements from this <tt>row</tt> by setting them to zero. The
+     * function does not modify the number of allocated nonzero entries, it
+     * only sets some entries to zero. It may drop them from the sparsity
+     * pattern, though (but retains the allocated memory in case new entries
+     * are again added later). Note that this is a global operation, so this
+     * needs to be done on all MPI processes.
+     *
+     * This operation is used in eliminating constraints (e.g. due to hanging
+     * nodes) and makes sure that we can write this modification to the matrix
+     * without having to read entries (such as the locations of non-zero
+     * elements) from it &mdash; without this operation, removing constraints
+     * on %parallel matrices is a rather complicated procedure.
+     *
+     * The second parameter can be used to set the diagonal entry of this row
+     * to a value different from zero. The default is to set it to zero.
      */
     void clear_row (const size_type      row,
                     const TrilinosScalar new_diag_value = 0);
 
     /**
-     * Same as clear_row(), except
-     * that it works on a number of
-     * rows at once.
+     * Same as clear_row(), except that it works on a number of rows at once.
      *
-     * The second parameter can be
-     * used to set the diagonal
-     * entries of all cleared rows
-     * to something different from
-     * zero. Note that all of these
-     * diagonal entries get the
-     * same value -- if you want
-     * different values for the
-     * diagonal entries, you have
-     * to set them by hand.
+     * The second parameter can be used to set the diagonal entries of all
+     * cleared rows to something different from zero. Note that all of these
+     * diagonal entries get the same value -- if you want different values for
+     * the diagonal entries, you have to set them by hand.
      */
     void clear_rows (const std::vector<size_type> &rows,
                      const TrilinosScalar          new_diag_value = 0);
 
     /**
-     * Make an in-place transpose
-     * of a matrix.
+     * Make an in-place transpose of a matrix.
      */
     void transpose ();
 
@@ -1789,77 +1323,40 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Return the value of the
-     * entry (<i>i,j</i>).  This
-     * may be an expensive
-     * operation and you should
-     * always take care where to
-     * call this function. As in
-     * the deal.II sparse matrix
-     * class, we throw an exception
-     * if the respective entry
-     * doesn't exist in the
-     * sparsity pattern of this
-     * class, which is requested
-     * from Trilinos. Moreover, an
-     * exception will be thrown
-     * when the requested element
-     * is not saved on the calling
-     * process.
+     * Return the value of the entry (<i>i,j</i>).  This may be an expensive
+     * operation and you should always take care where to call this
+     * function. As in the deal.II sparse matrix class, we throw an exception
+     * if the respective entry doesn't exist in the sparsity pattern of this
+     * class, which is requested from Trilinos. Moreover, an exception will be
+     * thrown when the requested element is not saved on the calling process.
      */
     TrilinosScalar operator () (const size_type i,
                                 const size_type j) const;
 
     /**
-     * Return the value of the
-     * matrix entry
-     * (<i>i,j</i>). If this entry
-     * does not exist in the
-     * sparsity pattern, then zero
-     * is returned. While this may
-     * be convenient in some cases,
-     * note that it is simple to
-     * write algorithms that are
-     * slow compared to an optimal
-     * solution, since the sparsity
-     * of the matrix is not used.
-     * On the other hand, if you
-     * want to be sure the entry
-     * exists, you should use
-     * operator() instead.
-     *
-     * The lack of error checking
-     * in this function can also
-     * yield surprising results if
-     * you have a parallel
-     * matrix. In that case, just
-     * because you get a zero
-     * result from this function
-     * does not mean that either
-     * the entry does not exist in
-     * the sparsity pattern or that
-     * it does but has a value of
-     * zero. Rather, it could also
-     * be that it simply isn't
-     * stored on the current
-     * processor; in that case, it
-     * may be stored on a different
-     * processor, and possibly so
-     * with a nonzero value.
+     * Return the value of the matrix entry (<i>i,j</i>). If this entry does
+     * not exist in the sparsity pattern, then zero is returned. While this
+     * may be convenient in some cases, note that it is simple to write
+     * algorithms that are slow compared to an optimal solution, since the
+     * sparsity of the matrix is not used.  On the other hand, if you want to
+     * be sure the entry exists, you should use operator() instead.
+     *
+     * The lack of error checking in this function can also yield surprising
+     * results if you have a parallel matrix. In that case, just because you
+     * get a zero result from this function does not mean that either the
+     * entry does not exist in the sparsity pattern or that it does but has a
+     * value of zero. Rather, it could also be that it simply isn't stored on
+     * the current processor; in that case, it may be stored on a different
+     * processor, and possibly so with a nonzero value.
      */
     TrilinosScalar el (const size_type i,
                        const size_type j) const;
 
     /**
-     * Return the main diagonal
-     * element in the <i>i</i>th
-     * row. This function throws an
-     * error if the matrix is not
-     * quadratic and it also throws
-     * an error if <i>(i,i)</i> is not
-     * element of the local matrix.
-     * See also the comment in
-     * trilinos_sparse_matrix.cc.
+     * Return the main diagonal element in the <i>i</i>th row. This function
+     * throws an error if the matrix is not quadratic and it also throws an
+     * error if <i>(i,i)</i> is not element of the local matrix.  See also the
+     * comment in trilinos_sparse_matrix.cc.
      */
     TrilinosScalar diag_element (const size_type i) const;
 
@@ -1975,155 +1472,86 @@ namespace TrilinosWrappers
                      const VectorType &src) const;
 
     /**
-     * Return the square of the norm
-     * of the vector $v$ with respect
-     * to the norm induced by this
-     * matrix, i.e.,
-     * $\left(v,Mv\right)$. This is
-     * useful, e.g. in the finite
-     * element context, where the
-     * $L_2$ norm of a function
-     * equals the matrix norm with
-     * respect to the mass matrix of
-     * the vector representing the
-     * nodal values of the finite
-     * element function.
-     *
-     * Obviously, the matrix needs to
-     * be quadratic for this
-     * operation.
-     *
-     * The implementation of this
-     * function is not as efficient
-     * as the one in the @p
-     * SparseMatrix class used in
-     * deal.II (i.e. the original
-     * one, not the Trilinos wrapper
-     * class) since Trilinos doesn't
-     * support this operation and
-     * needs a temporary vector.
-     *
-     * Note that both vectors have to
-     * be distributed vectors
-     * generated using the same Map
-     * as was used for the matrix in
-     * case you work on a distributed
-     * memory architecture, using the
-     * interface in the
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::Vector).
-     *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
+     * Return the square of the norm of the vector $v$ with respect to the
+     * norm induced by this matrix, i.e., $\left(v,Mv\right)$. This is useful,
+     * e.g. in the finite element context, where the $L_2$ norm of a function
+     * equals the matrix norm with respect to the mass matrix of the vector
+     * representing the nodal values of the finite element function.
+     *
+     * Obviously, the matrix needs to be quadratic for this operation.
+     *
+     * The implementation of this function is not as efficient as the one in
+     * the @p SparseMatrix class used in deal.II (i.e. the original one, not
+     * the Trilinos wrapper class) since Trilinos doesn't support this
+     * operation and needs a temporary vector.
+     *
+     * Note that both vectors have to be distributed vectors generated using
+     * the same Map as was used for the matrix in case you work on a
+     * distributed memory architecture, using the interface in the
+     * TrilinosWrappers::VectorBase class (or one of the two derived classes
+     * Vector and MPI::Vector).
+     *
+     * In case of a localized Vector, this function will only work when
+     * running on one processor, since the matrix object is inherently
+     * distributed. Otherwise, and exception will be thrown.
      */
     TrilinosScalar matrix_norm_square (const VectorBase &v) const;
 
     /**
-     * Compute the matrix scalar
-     * product $\left(u,Mv\right)$.
-     *
-     * The implementation of this
-     * function is not as efficient
-     * as the one in the @p
-     * SparseMatrix class used in
-     * deal.II (i.e. the original
-     * one, not the Trilinos
-     * wrapper class) since
-     * Trilinos doesn't support
-     * this operation and needs a
-     * temporary vector.
-     *
-     * Note that both vectors have to
-     * be distributed vectors
-     * generated using the same Map
-     * as was used for the matrix in
-     * case you work on a distributed
-     * memory architecture, using the
-     * interface in the
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::Vector).
-     *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
+     * Compute the matrix scalar product $\left(u,Mv\right)$.
+     *
+     * The implementation of this function is not as efficient as the one in
+     * the @p SparseMatrix class used in deal.II (i.e. the original one, not
+     * the Trilinos wrapper class) since Trilinos doesn't support this
+     * operation and needs a temporary vector.
+     *
+     * Note that both vectors have to be distributed vectors generated using
+     * the same Map as was used for the matrix in case you work on a
+     * distributed memory architecture, using the interface in the
+     * TrilinosWrappers::VectorBase class (or one of the two derived classes
+     * Vector and MPI::Vector).
+     *
+     * In case of a localized Vector, this function will only work when
+     * running on one processor, since the matrix object is inherently
+     * distributed. Otherwise, and exception will be thrown.
      */
     TrilinosScalar matrix_scalar_product (const VectorBase &u,
                                           const VectorBase &v) const;
 
     /**
-     * Compute the residual of an
-     * equation <i>Mx=b</i>, where
-     * the residual is defined to
-     * be <i>r=b-Mx</i>. Write the
-     * residual into @p dst. The
-     * <i>l<sub>2</sub></i> norm of
-     * the residual vector is
-     * returned.
-     *
-     * Source <i>x</i> and
-     * destination <i>dst</i> must
-     * not be the same vector.
-     *
-     * Note that both vectors have to
-     * be distributed vectors
-     * generated using the same Map
-     * as was used for the matrix in
-     * case you work on a distributed
-     * memory architecture, using the
-     * interface in the
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::Vector).
-     *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
+     * Compute the residual of an equation <i>Mx=b</i>, where the residual is
+     * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The
+     * <i>l<sub>2</sub></i> norm of the residual vector is returned.
+     *
+     * Source <i>x</i> and destination <i>dst</i> must not be the same vector.
+     *
+     * Note that both vectors have to be distributed vectors generated using
+     * the same Map as was used for the matrix in case you work on a
+     * distributed memory architecture, using the interface in the
+     * TrilinosWrappers::VectorBase class (or one of the two derived classes
+     * Vector and MPI::Vector).
+     *
+     * In case of a localized Vector, this function will only work when
+     * running on one processor, since the matrix object is inherently
+     * distributed. Otherwise, and exception will be thrown.
      */
     TrilinosScalar residual (VectorBase       &dst,
                              const VectorBase &x,
                              const VectorBase &b) const;
 
     /**
-    * Perform the matrix-matrix
-    * multiplication <tt>C = A * B</tt>,
-    * or, if an optional vector argument
-    * is given, <tt>C = A * diag(V) *
-    * B</tt>, where <tt>diag(V)</tt>
-    * defines a diagonal matrix with the
-    * vector entries.
+    * Perform the matrix-matrix multiplication <tt>C = A * B</tt>, or, if an
+    * optional vector argument is given, <tt>C = A * diag(V) * B</tt>, where
+    * <tt>diag(V)</tt> defines a diagonal matrix with the vector entries.
     *
-    * This function assumes that the
-    * calling matrix <tt>A</tt> and
-    * <tt>B</tt> have compatible
-    * sizes. The size of <tt>C</tt> will
-    * be set within this function.
+    * This function assumes that the calling matrix <tt>A</tt> and <tt>B</tt>
+    * have compatible sizes. The size of <tt>C</tt> will be set within this
+    * function.
     *
-    * The content as well as the sparsity
-    * pattern of the matrix C will be
-    * changed by this function, so make
-    * sure that the sparsity pattern is
-    * not used somewhere else in your
-    * program. This is an expensive
-    * operation, so think twice before you
-    * use this function.
+    * The content as well as the sparsity pattern of the matrix C will be
+    * changed by this function, so make sure that the sparsity pattern is not
+    * used somewhere else in your program. This is an expensive operation, so
+    * think twice before you use this function.
     */
     void mmult (SparseMatrix       &C,
                 const SparseMatrix &B,
@@ -2131,30 +1559,20 @@ namespace TrilinosWrappers
 
 
     /**
-    * Perform the matrix-matrix
-    * multiplication with the transpose of
-    * <tt>this</tt>, i.e., <tt>C =
-    * A<sup>T</sup> * B</tt>, or, if an
-    * optional vector argument is given,
-    * <tt>C = A<sup>T</sup> * diag(V) *
-    * B</tt>, where <tt>diag(V)</tt>
-    * defines a diagonal matrix with the
-    * vector entries.
+    * Perform the matrix-matrix multiplication with the transpose of
+    * <tt>this</tt>, i.e., <tt>C = A<sup>T</sup> * B</tt>, or, if an optional
+    * vector argument is given, <tt>C = A<sup>T</sup> * diag(V) * B</tt>,
+    * where <tt>diag(V)</tt> defines a diagonal matrix with the vector
+    * entries.
     *
-    * This function assumes that the
-    * calling matrix <tt>A</tt> and
-    * <tt>B</tt> have compatible
-    * sizes. The size of <tt>C</tt> will
-    * be set within this function.
+    * This function assumes that the calling matrix <tt>A</tt> and <tt>B</tt>
+    * have compatible sizes. The size of <tt>C</tt> will be set within this
+    * function.
     *
-    * The content as well as the sparsity
-    * pattern of the matrix C will be
-    * changed by this function, so make
-    * sure that the sparsity pattern is
-    * not used somewhere else in your
-    * program. This is an expensive
-    * operation, so think twice before you
-    * use this function.
+    * The content as well as the sparsity pattern of the matrix C will be
+    * changed by this function, so make sure that the sparsity pattern is not
+    * used somewhere else in your program. This is an expensive operation, so
+    * think twice before you use this function.
     */
     void Tmmult (SparseMatrix       &C,
                  const SparseMatrix &B,
@@ -2167,46 +1585,27 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Return the
-     * <i>l</i><sub>1</sub>-norm of
-     * the matrix, that is
-     * $|M|_1=
-     * \max_{\mathrm{all\ columns\ } j}
-     * \sum_{\mathrm{all\ rows\ } i}
-     * |M_{ij}|$, (max. sum
-     * of columns).  This is the
-     * natural matrix norm that is
-     * compatible to the l1-norm for
-     * vectors, i.e.  $|Mv|_1 \leq
-     * |M|_1 |v|_1$.
-     * (cf. Haemmerlin-Hoffmann:
-     * Numerische Mathematik)
+     * Return the <i>l</i><sub>1</sub>-norm of the matrix, that is $|M|_1=
+     * \max_{\mathrm{all\ columns\ } j} \sum_{\mathrm{all\ rows\ } i}
+     * |M_{ij}|$, (max. sum of columns).  This is the natural matrix norm that
+     * is compatible to the l1-norm for vectors, i.e.  $|Mv|_1 \leq |M|_1
+     * |v|_1$.  (cf. Haemmerlin-Hoffmann: Numerische Mathematik)
      */
     TrilinosScalar l1_norm () const;
 
     /**
-     * Return the linfty-norm of the
-     * matrix, that is
-     * $|M|_\infty=\max_{\mathrm{all\
-     * rows\ } i}\sum_{\mathrm{all\
-     * columns\ } j} |M_{ij}|$,
-     * (max. sum of rows).  This is
-     * the natural matrix norm that
-     * is compatible to the
-     * linfty-norm of vectors, i.e.
-     * $|Mv|_\infty \leq |M|_\infty
-     * |v|_\infty$.
-     * (cf. Haemmerlin-Hoffmann:
-     * Numerische Mathematik)
+     * Return the linfty-norm of the matrix, that is
+     * $|M|_\infty=\max_{\mathrm{all\ rows\ } i}\sum_{\mathrm{all\ columns\ }
+     * j} |M_{ij}|$, (max. sum of rows).  This is the natural matrix norm that
+     * is compatible to the linfty-norm of vectors, i.e.  $|Mv|_\infty \leq
+     * |M|_\infty |v|_\infty$.  (cf. Haemmerlin-Hoffmann: Numerische
+     * Mathematik)
      */
     TrilinosScalar linfty_norm () const;
 
     /**
-     * Return the frobenius norm of
-     * the matrix, i.e. the square
-     * root of the sum of squares
-     * of all entries in the
-     * matrix.
+     * Return the frobenius norm of the matrix, i.e. the square root of the
+     * sum of squares of all entries in the matrix.
      */
     TrilinosScalar frobenius_norm () const;
 
@@ -2217,60 +1616,44 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos
-     * Epetra_CrsMatrix data.
+     * Return a const reference to the underlying Trilinos Epetra_CrsMatrix
+     * data.
      */
     const Epetra_CrsMatrix &trilinos_matrix () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos
-     * Epetra_CrsGraph data that stores
-     * the sparsity pattern of the
-     * matrix.
+     * Return a const reference to the underlying Trilinos Epetra_CrsGraph
+     * data that stores the sparsity pattern of the matrix.
      */
     const Epetra_CrsGraph &trilinos_sparsity_pattern () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * domain space of this matrix, i.e.,
-     * the partitioning of the vectors
-     * this matrix has to be multiplied
-     * with.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the domain space of this matrix, i.e., the
+     * partitioning of the vectors this matrix has to be multiplied with.
      */
     const Epetra_Map &domain_partitioner () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * range space of this matrix, i.e.,
-     * the partitioning of the vectors
-     * that are result from matrix-vector
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the range space of this matrix, i.e., the
+     * partitioning of the vectors that are result from matrix-vector
      * products.
      */
     const Epetra_Map &range_partitioner () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * matrix rows. Equal to the
-     * partitioning of the range.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the matrix rows. Equal to the partitioning of
+     * the range.
      */
     const Epetra_Map &row_partitioner () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * matrix columns. This is in general
-     * not equal to the partitioner
-     * Epetra_Map for the domain because
-     * of overlap in the matrix.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the matrix columns. This is in general not
+     * equal to the partitioner Epetra_Map for the domain because of overlap
+     * in the matrix.
      */
     const Epetra_Map &col_partitioner () const;
 //@}
@@ -2280,8 +1663,7 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * STL-like iterator with the
-     * first entry.
+     * STL-like iterator with the first entry.
      */
     const_iterator begin () const;
 
@@ -2291,40 +1673,27 @@ namespace TrilinosWrappers
     const_iterator end () const;
 
     /**
-     * STL-like iterator with the
-     * first entry of row @p r.
+     * STL-like iterator with the first entry of row @p r.
      *
-     * Note that if the given row
-     * is empty, i.e. does not
-     * contain any nonzero entries,
-     * then the iterator returned
-     * by this function equals
-     * <tt>end(r)</tt>. Note also
-     * that the iterator may not be
-     * dereferencable in that case.
+     * Note that if the given row is empty, i.e. does not contain any nonzero
+     * entries, then the iterator returned by this function equals
+     * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable
+     * in that case.
      */
     const_iterator begin (const size_type r) const;
 
     /**
-     * Final iterator of row
-     * <tt>r</tt>. It points to the
-     * first element past the end
-     * of line @p r, or past the
-     * end of the entire sparsity
-     * pattern.
+     * Final iterator of row <tt>r</tt>. It points to the first element past
+     * the end of line @p r, or past the end of the entire sparsity pattern.
      *
-     * Note that the end iterator
-     * is not necessarily
-     * dereferencable. This is in
-     * particular the case if it is
-     * the end iterator for the
-     * last row of a matrix.
+     * Note that the end iterator is not necessarily dereferencable. This is
+     * in particular the case if it is the end iterator for the last row of a
+     * matrix.
      */
     const_iterator end (const size_type r) const;
 
     /**
-     * STL-like iterator with the
-     * first entry.
+     * STL-like iterator with the first entry.
      */
     iterator begin ();
 
@@ -2334,34 +1703,22 @@ namespace TrilinosWrappers
     iterator end ();
 
     /**
-     * STL-like iterator with the
-     * first entry of row @p r.
+     * STL-like iterator with the first entry of row @p r.
      *
-     * Note that if the given row
-     * is empty, i.e. does not
-     * contain any nonzero entries,
-     * then the iterator returned
-     * by this function equals
-     * <tt>end(r)</tt>. Note also
-     * that the iterator may not be
-     * dereferencable in that case.
+     * Note that if the given row is empty, i.e. does not contain any nonzero
+     * entries, then the iterator returned by this function equals
+     * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable
+     * in that case.
      */
     iterator begin (const size_type r);
 
     /**
-     * Final iterator of row
-     * <tt>r</tt>. It points to the
-     * first element past the end
-     * of line @p r, or past the
-     * end of the entire sparsity
-     * pattern.
+     * Final iterator of row <tt>r</tt>. It points to the first element past
+     * the end of line @p r, or past the end of the entire sparsity pattern.
      *
-     * Note that the end iterator
-     * is not necessarily
-     * dereferencable. This is in
-     * particular the case if it is
-     * the end iterator for the
-     * last row of a matrix.
+     * Note that the end iterator is not necessarily dereferencable. This is
+     * in particular the case if it is the end iterator for the last row of a
+     * matrix.
      */
     iterator end (const size_type r);
 
@@ -2372,27 +1729,18 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Abstract Trilinos object
-     * that helps view in ASCII
-     * other Trilinos
-     * objects. Currently this
-     * function is not
-     * implemented.  TODO: Not
+     * Abstract Trilinos object that helps view in ASCII other Trilinos
+     * objects. Currently this function is not implemented.  TODO: Not
      * implemented.
      */
     void write_ascii ();
 
     /**
-     * Print the matrix to the given
-     * stream, using the format
-     * <tt>(line,col) value</tt>, i.e. one
-     * nonzero entry of the matrix per
-     * line. The optional flag outputs the
-     * sparsity pattern in Trilinos style,
-     * where the data is sorted according
-     * to the processor number when printed
-     * to the stream, as well as a summary
-     * of the matrix like the global size.
+     * Print the matrix to the given stream, using the format <tt>(line,col)
+     * value</tt>, i.e. one nonzero entry of the matrix per line. The optional
+     * flag outputs the sparsity pattern in Trilinos style, where the data is
+     * sorted according to the processor number when printed to the stream, as
+     * well as a summary of the matrix like the global size.
      */
     void print (std::ostream &out,
                 const bool    write_extended_trilinos_info = false) const;
@@ -2455,35 +1803,20 @@ namespace TrilinosWrappers
   protected:
 
     /**
-    * For some matrix storage
-    * formats, in particular for the
-    * PETSc distributed blockmatrices,
-    * set and add operations on
-    * individual elements can not be
-    * freely mixed. Rather, one has
-    * to synchronize operations when
-    * one wants to switch from
-    * setting elements to adding to
-    * elements.
-    * BlockMatrixBase automatically
-    * synchronizes the access by
-    * calling this helper function
-    * for each block.
-    * This function ensures that the
-    * matrix is in a state that
-    * allows adding elements; if it
-    * previously already was in this
-    * state, the function does
-    * nothing.
+    * For some matrix storage formats, in particular for the PETSc distributed
+    * blockmatrices, set and add operations on individual elements can not be
+    * freely mixed. Rather, one has to synchronize operations when one wants
+    * to switch from setting elements to adding to elements.  BlockMatrixBase
+    * automatically synchronizes the access by calling this helper function
+    * for each block.  This function ensures that the matrix is in a state
+    * that allows adding elements; if it previously already was in this state,
+    * the function does nothing.
     */
     void prepare_add();
 
     /**
-    * Same as prepare_add() but
-    * prepare the matrix for setting
-    * elements if the representation
-    * of elements in this class
-    * requires such an operation.
+    * Same as prepare_add() but prepare the matrix for setting elements if the
+    * representation of elements in this class requires such an operation.
     */
     void prepare_set();
 
@@ -2492,60 +1825,46 @@ namespace TrilinosWrappers
   private:
 
     /**
-     * Pointer to the user-supplied
-     * Epetra Trilinos mapping of
-     * the matrix columns that
-     * assigns parts of the matrix
-     * to the individual processes.
+     * Pointer to the user-supplied Epetra Trilinos mapping of the matrix
+     * columns that assigns parts of the matrix to the individual processes.
      */
     std_cxx1x::shared_ptr<Epetra_Map> column_space_map;
 
     /**
-     * A sparse matrix object in
-     * Trilinos to be used for
-     * finite element based
-     * problems which allows for
-     * assembling into non-local
-     * elements.  The actual type,
-     * a sparse matrix, is set in
-     * the constructor.
+     * A sparse matrix object in Trilinos to be used for finite element based
+     * problems which allows for assembling into non-local elements.  The
+     * actual type, a sparse matrix, is set in the constructor.
      */
     std_cxx1x::shared_ptr<Epetra_FECrsMatrix> matrix;
 
     /**
-     * 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.
+     * A sparse matrix object in Trilinos to be used for collecting the
+     * non-local elements if the matrix was constructed from a Trilinos
+     * sparsity pattern with the respective option.
+     */
+    std_cxx1x::shared_ptr<Epetra_CrsMatrix> nonlocal_matrix;
+
+    /**
+     * 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;
 
     /**
-     *  To allow calling protected
-     *  prepare_add() and
-     *  prepare_set().
+     *  To allow calling protected prepare_add() and prepare_set().
      */
     friend class BlockMatrixBase<SparseMatrix>;
   };
@@ -2965,46 +2284,6 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  void
-  SparseMatrix::compress (::dealii::VectorOperation::values operation)
-  {
-
-    Epetra_CombineMode mode = last_action;
-    if (last_action == Zero)
-      {
-        if ((operation==::dealii::VectorOperation::add) ||
-            (operation==::dealii::VectorOperation::unknown))
-          mode = Add;
-        else if (operation==::dealii::VectorOperation::insert)
-          mode = Insert;
-      }
-    else
-      {
-        Assert(
-          ((last_action == Add) && (operation!=::dealii::VectorOperation::insert))
-          ||
-          ((last_action == Insert) && (operation!=::dealii::VectorOperation::add)),
-          ExcMessage("operation and argument to compress() do not match"));
-      }
-
-    // flush buffers
-    int ierr;
-    ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(),
-                                   true, mode);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    ierr = matrix->OptimizeStorage ();
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    last_action = Zero;
-
-    compressed = true;
-  }
-
-
-
   inline
   void
   SparseMatrix::compress ()
@@ -3023,19 +2302,17 @@ namespace TrilinosWrappers
 
     const int ierr = matrix->PutScalar(d);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+    if (nonlocal_matrix.get() != 0)
+      nonlocal_matrix->PutScalar(d);
 
     return *this;
   }
 
 
 
-  // Inline the set() and add()
-  // functions, since they will be
-  // called frequently, and the
-  // compiler can optimize away
-  // some unnecessary loops when
-  // the sizes are given at
-  // compile time.
+  // Inline the set() and add() functions, since they will be called
+  // frequently, and the compiler can optimize away some unnecessary loops
+  // when the sizes are given at compile time.
   inline
   void
   SparseMatrix::set (const size_type      i,
@@ -3208,7 +2485,7 @@ namespace TrilinosWrappers
         // standard Insert/ReplaceGlobalValues function. Nevertheless, the way
         // we call it is the fastest one (any other will lead to repeated
         // allocation and deallocation of memory in order to call the function
-        // we already use, which is very unefficient if writing one element at
+        // we already use, which is very inefficient if writing one element at
         // a time).
         compressed = false;
 
@@ -3222,6 +2499,19 @@ 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,
@@ -3413,6 +2703,20 @@ namespace TrilinosWrappers
                                                              col_value_ptr,
                                                              col_index_ptr);
       }
+    else if (nonlocal_matrix.get() != 0)
+      {
+        compressed = false;
+        // 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->SumIntoGlobalValues(row, n_columns,
+                                                    col_value_ptr,
+                                                    col_index_ptr);
+      }
     else
       {
         // When we're at off-processor data, we have to stick with the
index 3519bf8af1810c6fd0cec3658342afec3dbd748f..01d6d16676600c06a878ff43313b3da793c1a575 100644 (file)
@@ -93,20 +93,17 @@ namespace TrilinosWrappers
       Accessor (const Accessor &a);
 
       /**
-       * Row number of the element
-       * represented by this object.
+       * Row number of the element represented by this object.
        */
       size_type row() const;
 
       /**
-       * Index in row of the element
-       * represented by this object.
+       * Index in row of the element represented by this object.
        */
       size_type index() const;
 
       /**
-       * Column number of the element
-       * represented by this object.
+       * Column number of the element represented by this object.
        */
       size_type column() const;
 
@@ -142,45 +139,28 @@ namespace TrilinosWrappers
       size_type a_index;
 
       /**
-       * Cache where we store the
-       * column indices of the
-       * present row. This is
-       * necessary, since Trilinos
-       * makes access to the elements
-       * of its matrices rather hard,
-       * and it is much more
-       * efficient to copy all column
-       * entries of a row once when
-       * we enter it than repeatedly
-       * asking Trilinos for
-       * individual ones. This also
-       * makes some sense since it is
-       * likely that we will access
-       * them sequentially anyway.
+       * Cache where we store the column indices of the present row. This is
+       * necessary, since Trilinos makes access to the elements of its
+       * matrices rather hard, and it is much more efficient to copy all
+       * column entries of a row once when we enter it than repeatedly asking
+       * Trilinos for individual ones. This also makes some sense since it is
+       * likely that we will access them sequentially anyway.
        *
-       * In order to make copying of
-       * iterators/accessor of
-       * acceptable performance, we
-       * keep a shared pointer to
-       * these entries so that more
-       * than one accessor can access
-       * this data if necessary.
+       * In order to make copying of iterators/accessor of acceptable
+       * performance, we keep a shared pointer to these entries so that more
+       * than one accessor can access this data if necessary.
        */
       std_cxx1x::shared_ptr<const std::vector<size_type> > colnum_cache;
 
       /**
-       * Discard the old row caches
-       * (they may still be used by
-       * other accessors) and
-       * generate new ones for the
-       * row pointed to presently by
+       * Discard the old row caches (they may still be used by other
+       * accessors) and generate new ones for the row pointed to presently by
        * this accessor.
        */
       void visit_present_row ();
 
       /**
-       * Make enclosing class a
-       * friend.
+       * Make enclosing class a friend.
        */
       friend class Iterator;
     };
@@ -199,10 +179,8 @@ namespace TrilinosWrappers
       typedef dealii::types::global_dof_index size_type;
 
       /**
-       * Constructor. Create an
-       * iterator into the matrix @p
-       * matrix for the given row and
-       * the index within it.
+       * Constructor. Create an iterator into the matrix @p matrix for the
+       * given row and the index within it.
        */
       Iterator (const SparsityPattern *sparsity_pattern,
                 const size_type        row,
@@ -234,9 +212,8 @@ namespace TrilinosWrappers
       const Accessor *operator-> () const;
 
       /**
-       * Comparison. True, if both
-       * iterators point to the same
-       * matrix position.
+       * Comparison. True, if both iterators point to the same matrix
+       * position.
        */
       bool operator == (const Iterator &) const;
 
@@ -246,11 +223,8 @@ namespace TrilinosWrappers
       bool operator != (const Iterator &) const;
 
       /**
-       * Comparison operator. Result
-       * is true if either the first
-       * row number is smaller or if
-       * the row numbers are equal
-       * and the first index is
+       * Comparison operator. Result is true if either the first row number is
+       * smaller or if the row numbers are equal and the first index is
        * smaller.
        */
       bool operator < (const Iterator &) const;
@@ -266,8 +240,7 @@ namespace TrilinosWrappers
 
     private:
       /**
-       * Store an object of the
-       * accessor class.
+       * Store an object of the accessor class.
        */
       Accessor accessor;
 
@@ -317,92 +290,61 @@ namespace TrilinosWrappers
      */
 //@{
     /**
-     * Default constructor. Generates an
-     * empty (zero-size) sparsity
-     * pattern.
+     * Default constructor. Generates an empty (zero-size) sparsity pattern.
      */
     SparsityPattern ();
 
     /**
-     * Generate a sparsity pattern that is
-     * completely stored locally, having
-     * $m$ rows and $n$ columns. The
-     * resulting matrix will be completely
+     * Generate a sparsity pattern that is completely stored locally, having
+     * $m$ rows and $n$ columns. The resulting matrix will be completely
      * stored locally, too.
      *
-     * It is possible to specify the
-     * number of columns entries per row
-     * using the optional @p
-     * n_entries_per_row
-     * argument. However, this value does
-     * not need to be accurate or even
-     * given at all, since one does
-     * usually not have this kind of
-     * information before building the
-     * sparsity pattern (the usual case
-     * when the function
-     * DoFTools::make_sparsity_pattern()
-     * is called). The entries are
-     * allocated dynamically in a similar
-     * manner as for the deal.II
-     * CompressedSparsityPattern
-     * classes. However, a good estimate
-     * will reduce the setup time of the
-     * sparsity pattern.
+     * It is possible to specify the number of columns entries per row using
+     * the optional @p n_entries_per_row argument. However, this value does
+     * not need to be accurate or even given at all, since one does usually
+     * not have this kind of information before building the sparsity pattern
+     * (the usual case when the function DoFTools::make_sparsity_pattern() is
+     * called). The entries are allocated dynamically in a similar manner as
+     * for the deal.II CompressedSparsityPattern classes. However, a good
+     * estimate will reduce the setup time of the sparsity pattern.
      */
     SparsityPattern (const size_type  m,
                      const size_type  n,
                      const size_type  n_entries_per_row = 0);
 
     /**
-     * Generate a sparsity pattern that is
-     * completely stored locally, having
-     * $m$ rows and $n$ columns. The
-     * resulting matrix will be completely
+     * Generate a sparsity pattern that is completely stored locally, having
+     * $m$ rows and $n$ columns. The resulting matrix will be completely
      * stored locally, too.
      *
-     * The vector
-     * <tt>n_entries_per_row</tt>
-     * specifies the number of entries in
-     * each row (an information usually
-     * not available, though).
+     * The vector <tt>n_entries_per_row</tt> specifies the number of entries
+     * in each row (an information usually not available, though).
      */
     SparsityPattern (const size_type               m,
                      const size_type               n,
                      const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * Copy constructor. Sets the calling
-     * sparsity pattern to be the same as
+     * Copy constructor. Sets the calling sparsity pattern to be the same as
      * the input sparsity pattern.
      */
     SparsityPattern (const SparsityPattern &input_sparsity_pattern);
 
     /**
-     * Destructor. Made virtual so that
-     * one can use pointers to this
-     * class.
+     * Destructor. Made virtual so that one can use pointers to this class.
      */
     virtual ~SparsityPattern ();
 
     /**
-     * Initialize a sparsity pattern that
-     * is completely stored locally,
-     * having $m$ rows and $n$
-     * columns. The resulting matrix will
-     * be completely stored locally.
+     * Initialize a sparsity pattern that is completely stored locally, having
+     * $m$ rows and $n$ columns. The resulting matrix will be completely
+     * stored locally.
      *
-     * The number of columns entries per
-     * row is specified as the maximum
-     * number of entries argument.  This
-     * does not need to be an accurate
-     * number since the entries are
-     * allocated dynamically in a similar
-     * manner as for the deal.II
-     * CompressedSparsityPattern classes,
-     * but a good estimate will reduce
-     * the setup time of the sparsity
-     * pattern.
+     * The number of columns entries per row is specified as the maximum
+     * number of entries argument.  This does not need to be an accurate
+     * number since the entries are allocated dynamically in a similar manner
+     * as for the deal.II CompressedSparsityPattern classes, but a good
+     * estimate will reduce the setup time of the sparsity pattern.
      */
     void
     reinit (const size_type  m,
@@ -410,16 +352,12 @@ namespace TrilinosWrappers
             const size_type  n_entries_per_row = 0);
 
     /**
-     * Initialize a sparsity pattern that
-     * is completely stored locally,
-     * having $m$ rows and $n$ columns. The
-     * resulting matrix will be
-     * completely stored locally.
+     * Initialize a sparsity pattern that is completely stored locally, having
+     * $m$ rows and $n$ columns. The resulting matrix will be completely
+     * stored locally.
      *
-     * The vector
-     * <tt>n_entries_per_row</tt>
-     * specifies the number of entries in
-     * each row.
+     * The vector <tt>n_entries_per_row</tt> specifies the number of entries
+     * in each row.
      */
     void
     reinit (const size_type               m,
@@ -427,65 +365,46 @@ namespace TrilinosWrappers
             const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * Copy function. Sets the calling
-     * sparsity pattern to be the same as
-     * the input sparsity pattern.
+     * Copy function. Sets the calling sparsity pattern to be the same as the
+     * input sparsity pattern.
      */
     void
     copy_from (const SparsityPattern &input_sparsity_pattern);
 
     /**
-     * Copy function from one of the
-     * deal.II sparsity patterns. If used
-     * in parallel, this function uses an
-     * ad-hoc partitioning of the rows
-     * and columns.
+     * Copy function from one of the deal.II sparsity patterns. If used in
+     * parallel, this function uses an ad-hoc partitioning of the rows and
+     * columns.
      */
     template<typename SparsityType>
     void
     copy_from (const SparsityType &nontrilinos_sparsity_pattern);
 
     /**
-     * Copy operator. This operation is
-     * only allowed for empty objects, to
-     * avoid potentially very costly
-     * operations automatically
-     * synthesized by the compiler. Use
-     * copy_from() instead if you know
-     * that you really want to copy a
-     * sparsity pattern with non-trivial
-     * content.
+     * Copy operator. This operation is only allowed for empty objects, to
+     * avoid potentially very costly operations automatically synthesized by
+     * the compiler. Use copy_from() instead if you know that you really want
+     * to copy a sparsity pattern with non-trivial content.
      */
     SparsityPattern &operator = (const SparsityPattern &input_sparsity_pattern);
 
     /**
-     * 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.
      *
-     * This is a collective operation
-     * that needs to be called on all
-     * processors in order to avoid a
-     * dead lock.
+     * This is a collective operation that needs to be called on all
+     * processors in order to avoid a dead lock.
      */
     void clear ();
 
     /**
-     * In analogy to our own
-     * SparsityPattern class, this
-     * function compresses the sparsity
-     * pattern and allows the resulting
-     * pattern to be used for actually
-     * generating a (Trilinos-based)
-     * matrix. This function also
-     * exchanges non-local data that
-     * might have accumulated during the
-     * addition of new elements. This
-     * function must therefore be called
-     * once the structure is fixed. This
-     * is a collective operation, i.e.,
-     * it needs to be run on all
-     * processors when used in parallel.
+     * In analogy to our own SparsityPattern class, this function compresses
+     * the sparsity pattern and allows the resulting pattern to be used for
+     * actually generating a (Trilinos-based) matrix. This function also
+     * exchanges non-local data that might have accumulated during the
+     * addition of new elements. This function must therefore be called once
+     * the structure is fixed. This is a collective operation, i.e., it needs
+     * to be run on all processors when used in parallel.
      */
     void compress ();
 //@}
@@ -495,79 +414,47 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Constructor for a square sparsity
-     * pattern using an Epetra_map for
-     * the description of the %parallel
-     * partitioning. Moreover, the number
-     * of nonzero entries in the rows of
-     * the sparsity pattern can be
-     * specified. Note that this number
-     * does not need to be exact, and it
-     * is allowed that the actual
-     * sparsity structure has more
-     * nonzero entries than specified in
-     * the constructor (the usual case
-     * when the function
-     * DoFTools::make_sparsity_pattern()
-     * is called). However it is still
-     * advantageous to provide good
-     * estimates here since a good value
-     * will avoid repeated allocation of
-     * memory, which considerably
-     * increases the performance when
-     * creating the sparsity pattern.
+     * Constructor for a square sparsity pattern using an Epetra_map for the
+     * description of the %parallel partitioning. Moreover, the number of
+     * nonzero entries in the rows of the sparsity pattern can be
+     * specified. Note that this number does not need to be exact, and it is
+     * allowed that the actual sparsity structure has more nonzero entries
+     * than specified in the constructor (the usual case when the function
+     * DoFTools::make_sparsity_pattern() is called). However it is still
+     * advantageous to provide good estimates here since a good value will
+     * avoid repeated allocation of memory, which considerably increases the
+     * performance when creating the sparsity pattern.
      */
     SparsityPattern (const Epetra_Map &parallel_partitioning,
                      const size_type   n_entries_per_row = 0);
 
     /**
-     * Same as before, but now use the
-     * exact number of nonzeros in each m
-     * row. Since we know the number of
-     * elements in the sparsity pattern
-     * exactly in this case, we can
-     * already allocate the right amount
-     * of memory, which makes the
-     * creation process by the respective
-     * SparsityPattern::reinit call
-     * considerably faster. However, this
-     * is a rather unusual situation,
-     * since knowing the number of
-     * entries in each row is usually
-     * connected to knowing the indices
-     * of nonzero entries, which the
-     * sparsity pattern is designed to
-     * describe.
+     * Same as before, but now use the exact number of nonzeros in each m
+     * row. Since we know the number of elements in the sparsity pattern
+     * exactly in this case, we can already allocate the right amount of
+     * memory, which makes the creation process by the respective
+     * SparsityPattern::reinit call considerably faster. However, this is a
+     * rather unusual situation, since knowing the number of entries in each
+     * row is usually connected to knowing the indices of nonzero entries,
+     * which the sparsity pattern is designed to describe.
      */
     SparsityPattern (const Epetra_Map             &parallel_partitioning,
                      const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different Epetra maps for rows and
-     * columns. This interface is meant to
-     * be used for generating rectangular
-     * sparsity pattern, where one map
-     * describes the %parallel partitioning
-     * of the dofs associated with the
-     * sparsity pattern rows and the other
-     * one of the sparsity pattern
-     * columns. Note that there is no real
-     * parallelism along the columns
-     * &ndash; the processor that owns a
-     * certain row always owns all the
-     * column elements, no matter how far
-     * they might be spread out. The second
-     * Epetra_Map is only used to specify
-     * the number of columns and for
-     * specifying the correct domain space
-     * when performing matrix-vector
-     * products with vectors based on the
-     * same column map.
+     * This constructor is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular sparsity pattern, where one map
+     * describes the %parallel partitioning of the dofs associated with the
+     * sparsity pattern rows and the other one of the sparsity pattern
+     * columns. Note that there is no real parallelism along the columns
+     * &ndash; the processor that owns a certain row always owns all the
+     * column elements, no matter how far they might be spread out. The second
+     * Epetra_Map is only used to specify the number of columns and for
+     * specifying the correct domain space when performing matrix-vector
+     * products with vectors based on the same column map.
      *
-     * The number of columns entries per
-     * row is specified as the maximum
+     * The number of columns entries per row is specified as the maximum
      * number of entries argument.
      */
     SparsityPattern (const Epetra_Map   &row_parallel_partitioning,
@@ -575,74 +462,47 @@ namespace TrilinosWrappers
                      const size_type     n_entries_per_row = 0);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different Epetra maps for rows and
-     * columns. This interface is meant to
-     * be used for generating rectangular
-     * matrices, where one map specifies
-     * the %parallel distribution of rows
-     * and the second one specifies the
-     * distribution of degrees of freedom
-     * associated with matrix columns. This
-     * second map is however not used for
-     * the distribution of the columns
-     * themselves &ndash; rather, all
-     * column elements of a row are stored
-     * on the same processor. The vector
-     * <tt>n_entries_per_row</tt> specifies
-     * the number of entries in each row of
-     * the newly generated matrix.
+     * This constructor is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular matrices, where one map specifies
+     * the %parallel distribution of rows and the second one specifies the
+     * distribution of degrees of freedom associated with matrix columns. This
+     * second map is however not used for the distribution of the columns
+     * themselves &ndash; rather, all column elements of a row are stored on
+     * the same processor. The vector <tt>n_entries_per_row</tt> specifies the
+     * number of entries in each row of the newly generated matrix.
      */
     SparsityPattern (const Epetra_Map             &row_parallel_partitioning,
                      const Epetra_Map             &col_parallel_partitioning,
                      const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * Reinitialization function for
-     * generating a square sparsity pattern
-     * using an Epetra_Map for the
-     * description of the %parallel
-     * partitioning and the number of
-     * nonzero entries in the rows of the
-     * sparsity pattern. Note that this
-     * number does not need to be exact,
-     * and it is even allowed that the
-     * actual sparsity structure has more
-     * nonzero entries than specified in
-     * the constructor. However it is still
-     * advantageous to provide good
-     * estimates here since this will
-     * considerably increase the
-     * performance when creating the
-     * sparsity pattern.
+     * Reinitialization function for generating a square sparsity pattern
+     * using an Epetra_Map for the description of the %parallel partitioning
+     * and the number of nonzero entries in the rows of the sparsity
+     * pattern. Note that this number does not need to be exact, and it is
+     * even allowed that the actual sparsity structure has more nonzero
+     * entries than specified in the constructor. However it is still
+     * advantageous to provide good estimates here since this will
+     * considerably increase the performance when creating the sparsity
+     * pattern.
      *
-     * This function does not create any
-     * entries by itself, but provides
-     * the correct data structures that
-     * can be used by the respective
-     * add() function.
+     * This function does not create any entries by itself, but provides the
+     * correct data structures that can be used by the respective add()
+     * function.
      */
     void
     reinit (const Epetra_Map &parallel_partitioning,
             const size_type   n_entries_per_row = 0);
 
     /**
-     * Same as before, but now use the
-     * exact number of nonzeros in each m
-     * row. Since we know the number of
-     * elements in the sparsity pattern
-     * exactly in this case, we can
-     * already allocate the right amount
-     * of memory, which makes process of
-     * adding entries to the sparsity
-     * pattern considerably
-     * faster. However, this is a rather
-     * unusual situation, since knowing
-     * the number of entries in each row
-     * is usually connected to knowing
-     * the indices of nonzero entries,
-     * which the sparsity pattern is
+     * Same as before, but now use the exact number of nonzeros in each m
+     * row. Since we know the number of elements in the sparsity pattern
+     * exactly in this case, we can already allocate the right amount of
+     * memory, which makes process of adding entries to the sparsity pattern
+     * considerably faster. However, this is a rather unusual situation, since
+     * knowing the number of entries in each row is usually connected to
+     * knowing the indices of nonzero entries, which the sparsity pattern is
      * designed to describe.
      */
     void
@@ -650,31 +510,19 @@ namespace TrilinosWrappers
             const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * This reinit function is similar to
-     * the one above, but it now takes
-     * two different Epetra maps for rows
-     * and columns. This interface is
-     * meant to be used for generating
-     * rectangular sparsity pattern,
-     * where one map describes the
-     * %parallel partitioning of the dofs
-     * associated with the sparsity
-     * pattern rows and the other one of
-     * the sparsity pattern columns. Note
-     * that there is no real parallelism
-     * along the columns &ndash; the
-     * processor that owns a certain row
-     * always owns all the column
-     * elements, no matter how far they
-     * might be spread out. The second
-     * Epetra_Map is only used to specify
-     * the number of columns and for
-     * internal arragements when doing
-     * matrix-vector products with
-     * vectors based on that column map.
+     * This reinit function is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular sparsity pattern, where one map
+     * describes the %parallel partitioning of the dofs associated with the
+     * sparsity pattern rows and the other one of the sparsity pattern
+     * columns. Note that there is no real parallelism along the columns
+     * &ndash; the processor that owns a certain row always owns all the
+     * column elements, no matter how far they might be spread out. The second
+     * Epetra_Map is only used to specify the number of columns and for
+     * internal arragements when doing matrix-vector products with vectors
+     * based on that column map.
      *
-     * The number of columns entries per
-     * row is specified by the argument
+     * The number of columns entries per row is specified by the argument
      * <tt>n_entries_per_row</tt>.
      */
     void
@@ -683,27 +531,15 @@ namespace TrilinosWrappers
             const size_type     n_entries_per_row = 0);
 
     /**
-     * This reinit function is similar to
-     * the one above, but it now takes
-     * two different Epetra maps for rows
-     * and columns. This interface is
-     * meant to be used for generating
-     * rectangular matrices, where one
-     * map specifies the %parallel
-     * distribution of rows and the
-     * second one specifies the
-     * distribution of degrees of freedom
-     * associated with matrix
-     * columns. This second map is
-     * however not used for the
-     * distribution of the columns
-     * themselves &ndash; rather, all
-     * column elements of a row are
-     * stored on the same processor. The
-     * vector <tt>n_entries_per_row</tt>
-     * specifies the number of entries in
-     * each row of the newly generated
-     * matrix.
+     * This reinit function is similar to the one above, but it now takes two
+     * different Epetra maps for rows and columns. This interface is meant to
+     * be used for generating rectangular matrices, where one map specifies
+     * the %parallel distribution of rows and the second one specifies the
+     * distribution of degrees of freedom associated with matrix columns. This
+     * second map is however not used for the distribution of the columns
+     * themselves &ndash; rather, all column elements of a row are stored on
+     * the same processor. The vector <tt>n_entries_per_row</tt> specifies the
+     * number of entries in each row of the newly generated matrix.
      */
     void
     reinit (const Epetra_Map             &row_parallel_partitioning,
@@ -711,19 +547,12 @@ namespace TrilinosWrappers
             const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * Reinit function. Takes one of the
-     * deal.II sparsity patterns and a
-     * %parallel partitioning of the rows
-     * and columns for initializing the
-     * current Trilinos sparsity
-     * pattern. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
-     * CompressedSimpleSparsityPattern.
+     * Reinit function. Takes one of the deal.II sparsity patterns and a
+     * %parallel partitioning of the rows and columns for initializing the
+     * current Trilinos sparsity pattern. The optional argument @p
+     * exchange_data can be used for reinitialization with a sparsity pattern
+     * that is not fully constructed. This feature is only implemented for
+     * input sparsity patterns of type CompressedSimpleSparsityPattern.
      */
     template<typename SparsityType>
     void
@@ -733,19 +562,12 @@ namespace TrilinosWrappers
             const bool          exchange_data = false);
 
     /**
-     * Reinit function. Takes one of the
-     * deal.II sparsity patterns and a
-     * %parallel partitioning of the rows
-     * and columns for initializing the
-     * current Trilinos sparsity
-     * pattern. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
-     * CompressedSimpleSparsityPattern.
+     * Reinit function. Takes one of the deal.II sparsity patterns and a
+     * %parallel partitioning of the rows and columns for initializing the
+     * current Trilinos sparsity pattern. The optional argument @p
+     * exchange_data can be used for reinitialization with a sparsity pattern
+     * that is not fully constructed. This feature is only implemented for
+     * input sparsity patterns of type CompressedSimpleSparsityPattern.
      */
     template<typename SparsityType>
     void
@@ -759,74 +581,47 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Constructor for a square sparsity
-     * pattern using an IndexSet and an
-     * MPI communicator for the
-     * description of the %parallel
-     * partitioning. Moreover, the number
-     * of nonzero entries in the rows of
-     * the sparsity pattern can be
-     * specified. Note that this number
-     * does not need to be exact, and it
-     * is even allowed that the actual
-     * sparsity structure has more
-     * nonzero entries than specified in
-     * the constructor. However it is
-     * still advantageous to provide good
-     * estimates here since a good value
-     * will avoid repeated allocation of
-     * memory, which considerably
-     * increases the performance when
-     * creating the sparsity pattern.
+     * Constructor for a square sparsity pattern using an IndexSet and an MPI
+     * communicator for the description of the %parallel
+     * partitioning. Moreover, the number of nonzero entries in the rows of
+     * the sparsity pattern can be specified. Note that this number does not
+     * need to be exact, and it is even allowed that the actual sparsity
+     * structure has more nonzero entries than specified in the
+     * constructor. However it is still advantageous to provide good estimates
+     * here since a good value will avoid repeated allocation of memory, which
+     * considerably increases the performance when creating the sparsity
+     * pattern.
      */
     SparsityPattern (const IndexSet  &parallel_partitioning,
                      const MPI_Comm  &communicator = MPI_COMM_WORLD,
                      const size_type  n_entries_per_row = 0);
 
     /**
-     * Same as before, but now use the
-     * exact number of nonzeros in each m
-     * row. Since we know the number of
-     * elements in the sparsity pattern
-     * exactly in this case, we can
-     * already allocate the right amount
-     * of memory, which makes the
-     * creation process by the respective
-     * SparsityPattern::reinit call
-     * considerably faster. However, this
-     * is a rather unusual situation,
-     * since knowing the number of
-     * entries in each row is usually
-     * connected to knowing the indices
-     * of nonzero entries, which the
-     * sparsity pattern is designed to
-     * describe.
+     * Same as before, but now use the exact number of nonzeros in each m
+     * row. Since we know the number of elements in the sparsity pattern
+     * exactly in this case, we can already allocate the right amount of
+     * memory, which makes the creation process by the respective
+     * SparsityPattern::reinit call considerably faster. However, this is a
+     * rather unusual situation, since knowing the number of entries in each
+     * row is usually connected to knowing the indices of nonzero entries,
+     * which the sparsity pattern is designed to describe.
      */
     SparsityPattern (const IndexSet                  &parallel_partitioning,
                      const MPI_Comm                  &communicator,
                      const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different index sets to describe the
-     * %parallel partitioning of rows and
-     * columns. This interface is meant to
-     * be used for generating rectangular
-     * sparsity pattern. Note that there is
-     * no real parallelism along the
-     * columns &ndash; the processor that
-     * owns a certain row always owns all
-     * the column elements, no matter how
-     * far they might be spread out. The
-     * second Epetra_Map is only used to
-     * specify the number of columns and
-     * for internal arragements when doing
-     * matrix-vector products with vectors
+     * This constructor is similar to the one above, but it now takes two
+     * different index sets to describe the %parallel partitioning of rows and
+     * columns. This interface is meant to be used for generating rectangular
+     * sparsity pattern. Note that there is no real parallelism along the
+     * columns &ndash; the processor that owns a certain row always owns all
+     * the column elements, no matter how far they might be spread out. The
+     * second Epetra_Map is only used to specify the number of columns and for
+     * internal arragements when doing matrix-vector products with vectors
      * based on that column map.
      *
-     * The number of columns entries per
-     * row is specified as the maximum
+     * The number of columns entries per row is specified as the maximum
      * number of entries argument.
      */
     SparsityPattern (const IndexSet  &row_parallel_partitioning,
@@ -835,24 +630,15 @@ namespace TrilinosWrappers
                      const size_type  n_entries_per_row = 0);
 
     /**
-     * This constructor is similar to the
-     * one above, but it now takes two
-     * different index sets for rows and
-     * columns. This interface is meant to
-     * be used for generating rectangular
-     * matrices, where one map specifies
-     * the %parallel distribution of rows
-     * and the second one specifies the
-     * distribution of degrees of freedom
-     * associated with matrix columns. This
-     * second map is however not used for
-     * the distribution of the columns
-     * themselves &ndash; rather, all
-     * column elements of a row are stored
-     * on the same processor. The vector
-     * <tt>n_entries_per_row</tt> specifies
-     * the number of entries in each row of
-     * the newly generated matrix.
+     * This constructor is similar to the one above, but it now takes two
+     * different index sets for rows and columns. This interface is meant to
+     * be used for generating rectangular matrices, where one map specifies
+     * the %parallel distribution of rows and the second one specifies the
+     * distribution of degrees of freedom associated with matrix columns. This
+     * second map is however not used for the distribution of the columns
+     * themselves &ndash; rather, all column elements of a row are stored on
+     * the same processor. The vector <tt>n_entries_per_row</tt> specifies the
+     * number of entries in each row of the newly generated matrix.
      */
     SparsityPattern (const IndexSet               &row_parallel_partitioning,
                      const IndexSet               &col_parallel_partitioning,
@@ -860,30 +646,51 @@ namespace TrilinosWrappers
                      const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * Reinitialization function for
-     * generating a square sparsity
-     * pattern using an IndexSet and an
-     * MPI communicator for the
-     * description of the %parallel
-     * partitioning and the number of
-     * nonzero entries in the rows of the
-     * sparsity pattern. Note that this
-     * number does not need to be exact,
-     * and it is even allowed that the
-     * actual sparsity structure has more
-     * nonzero entries than specified in
-     * the constructor. However it is
-     * still advantageous to provide good
-     * estimates here since this will
-     * considerably increase the
-     * performance when creating the
-     * sparsity pattern.
+     * This constructor constructs general sparsity patterns, possible
+     * non-square ones. Constructing a sparsity pattern this way allows the
+     * user to explicitly specify the rows into which we are going to add
+     * elements. This set is required to be a superset of the first index set
+     * @p row_parallel_partitioning that includes also rows that are owned by
+     * another processor (ghost rows). Note that elements can only be added to
+     * rows specified by @p writable_rows.
+     *
+     * This method is beneficial when the rows to which a processor is going
+     * to write can be determined before actually inserting elements into the
+     * matrix. For the typical parallel::distributed::Triangulation class used
+     * in deal.II, we know that a processor only will add row elements for
+     * what we call the locally relevant dofs (see
+     * DoFTools::extract_locally_relevant_dofs). The other constructors
+     * methods use general Trilinos facilities that allow to add elements to
+     * arbitrary rows (as done by all the other reinit functions). However,
+     * this flexbility come at a cost, the most prominent being that adding
+     * elements into the same matrix from multiple threads in shared memory is
+     * not safe whenever MPI is used. For these settings, the current method
+     * is the one to choose: It will store the off-processor data as an
+     * additional sparsity pattern (that is then passed to the Trilinos matrix
+     * via the reinit mehtod) which can be organized in such a way that
+     * thread-safety can be ensured (as long as the user makes sure to never
+     * write into the same matrix row simultaneously, of course).
+     */
+    SparsityPattern (const IndexSet  &row_parallel_partitioning,
+                     const IndexSet  &col_parallel_partitioning,
+                     const IndexSet  &writable_rows,
+                     const MPI_Comm  &communicator = MPI_COMM_WORLD,
+                     const size_type  n_entries_per_row = 0);
+
+    /**
+     * Reinitialization function for generating a square sparsity pattern
+     * using an IndexSet and an MPI communicator for the description of the
+     * %parallel partitioning and the number of nonzero entries in the rows of
+     * the sparsity pattern. Note that this number does not need to be exact,
+     * and it is even allowed that the actual sparsity structure has more
+     * nonzero entries than specified in the constructor. However it is still
+     * advantageous to provide good estimates here since this will
+     * considerably increase the performance when creating the sparsity
+     * pattern.
      *
-     * This function does not create any
-     * entries by itself, but provides
-     * the correct data structures that
-     * can be used by the respective
-     * add() function.
+     * This function does not create any entries by itself, but provides the
+     * correct data structures that can be used by the respective add()
+     * function.
      */
     void
     reinit (const IndexSet  &parallel_partitioning,
@@ -891,21 +698,13 @@ namespace TrilinosWrappers
             const size_type  n_entries_per_row = 0);
 
     /**
-     * Same as before, but now use the
-     * exact number of nonzeros in each m
-     * row. Since we know the number of
-     * elements in the sparsity pattern
-     * exactly in this case, we can
-     * already allocate the right amount
-     * of memory, which makes process of
-     * adding entries to the sparsity
-     * pattern considerably
-     * faster. However, this is a rather
-     * unusual situation, since knowing
-     * the number of entries in each row
-     * is usually connected to knowing
-     * the indices of nonzero entries,
-     * which the sparsity pattern is
+     * Same as before, but now use the exact number of nonzeros in each m
+     * row. Since we know the number of elements in the sparsity pattern
+     * exactly in this case, we can already allocate the right amount of
+     * memory, which makes process of adding entries to the sparsity pattern
+     * considerably faster. However, this is a rather unusual situation, since
+     * knowing the number of entries in each row is usually connected to
+     * knowing the indices of nonzero entries, which the sparsity pattern is
      * designed to describe.
      */
     void
@@ -914,32 +713,19 @@ namespace TrilinosWrappers
             const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * This reinit function is similar to
-     * the one above, but it now takes
-     * two different index sets for rows
-     * and columns. This interface is
-     * meant to be used for generating
-     * rectangular sparsity pattern,
-     * where one index set describes the
-     * %parallel partitioning of the dofs
-     * associated with the sparsity
-     * pattern rows and the other one of
-     * the sparsity pattern columns. Note
-     * that there is no real parallelism
-     * along the columns &ndash; the
-     * processor that owns a certain row
-     * always owns all the column
-     * elements, no matter how far they
-     * might be spread out. The second
-     * IndexSet is only used to specify
-     * the number of columns and for
-     * internal arragements when doing
-     * matrix-vector products with
-     * vectors based on an EpetraMap
-     * based on that IndexSet.
+     * This reinit function is similar to the one above, but it now takes two
+     * different index sets for rows and columns. This interface is meant to
+     * be used for generating rectangular sparsity pattern, where one index
+     * set describes the %parallel partitioning of the dofs associated with
+     * the sparsity pattern rows and the other one of the sparsity pattern
+     * columns. Note that there is no real parallelism along the columns
+     * &ndash; the processor that owns a certain row always owns all the
+     * column elements, no matter how far they might be spread out. The second
+     * IndexSet is only used to specify the number of columns and for internal
+     * arragements when doing matrix-vector products with vectors based on an
+     * EpetraMap based on that IndexSet.
      *
-     * The number of columns entries per
-     * row is specified by the argument
+     * The number of columns entries per row is specified by the argument
      * <tt>n_entries_per_row</tt>.
      */
     void
@@ -949,11 +735,40 @@ namespace TrilinosWrappers
             const size_type  n_entries_per_row = 0);
 
     /**
-     * Same as before, but now using a
-     * vector <tt>n_entries_per_row</tt>
-     * for specifying the number of
-     * entries in each row of the
-     * sparsity pattern.
+     * This reinit function is used to specify general matrices, possibly
+     * non-square ones. In addition to the arguments of the other reinit
+     * method above, it allows the user to explicitly specify the rows into
+     * which we are going to add elements. This set is a superset of the first
+     * index set @p row_parallel_partitioning that includes also rows that are
+     * owned by another processor (ghost rows).
+     *
+     * This method is beneficial when the rows to which a processor is going
+     * to write can be determined before actually inserting elements into the
+     * matrix. For the typical parallel::distributed::Triangulation class used
+     * in deal.II, we know that a processor only will add row elements for
+     * what we call the locally relevant dofs (see
+     * DoFTools::extract_locally_relevant_dofs). Trilinos matrices allow to
+     * add elements to arbitrary rows (as done by all the other reinit
+     * functions) and this is what all the other reinit methods do,
+     * too. However, this flexbility come at a cost, the most prominent being
+     * that adding elements into the same matrix from multiple threads in
+     * shared memory is not safe whenever MPI is used. For these settings, the
+     * current method is the one to choose: It will store the off-processor
+     * data as an additional sparsity pattern (that is then passed to the
+     * Trilinos matrix via the reinit mehtod) which can be organized in such a
+     * way that thread-safety can be ensured (as long as the user makes sure
+     * to never write into the same matrix row simultaneously, of course).
+     */
+    void
+    reinit (const IndexSet  &row_parallel_partitioning,
+            const IndexSet  &col_parallel_partitioning,
+            const IndexSet  &writeable_rows,
+            const MPI_Comm  &communicator = MPI_COMM_WORLD,
+            const size_type  n_entries_per_row = 0);
+
+    /**
+     * Same as before, but now using a vector <tt>n_entries_per_row</tt> for
+     * specifying the number of entries in each row of the sparsity pattern.
      */
     void
     reinit (const IndexSet               &row_parallel_partitioning,
@@ -962,20 +777,13 @@ namespace TrilinosWrappers
             const std::vector<size_type> &n_entries_per_row);
 
     /**
-     * Reinit function. Takes one of the
-     * deal.II sparsity patterns and the
-     * %parallel partitioning of the rows
-     * and columns specified by two index
-     * sets and a %parallel communicator
-     * for initializing the current
-     * Trilinos sparsity pattern. The
-     * optional argument @p exchange_data
-     * can be used for reinitialization
-     * with a sparsity pattern that is
-     * not fully constructed. This
-     * feature is only implemented for
-     * input sparsity patterns of type
-     * CompressedSimpleSparsityPattern.
+     * Reinit function. Takes one of the deal.II sparsity patterns and the
+     * %parallel partitioning of the rows and columns specified by two index
+     * sets and a %parallel communicator for initializing the current Trilinos
+     * sparsity pattern. The optional argument @p exchange_data can be used
+     * for reinitialization with a sparsity pattern that is not fully
+     * constructed. This feature is only implemented for input sparsity
+     * patterns of type CompressedSimpleSparsityPattern.
      */
     template<typename SparsityType>
     void
@@ -986,19 +794,12 @@ namespace TrilinosWrappers
             const bool          exchange_data = false);
 
     /**
-     * Reinit function. Takes one of the
-     * deal.II sparsity patterns and a
-     * %parallel partitioning of the rows
-     * and columns for initializing the
-     * current Trilinos sparsity
-     * pattern. The optional argument @p
-     * exchange_data can be used for
-     * reinitialization with a sparsity
-     * pattern that is not fully
-     * constructed. This feature is only
-     * implemented for input sparsity
-     * patterns of type
-     * CompressedSimpleSparsityPattern.
+     * Reinit function. Takes one of the deal.II sparsity patterns and a
+     * %parallel partitioning of the rows and columns for initializing the
+     * current Trilinos sparsity pattern. The optional argument @p
+     * exchange_data can be used for reinitialization with a sparsity pattern
+     * that is not fully constructed. This feature is only implemented for
+     * input sparsity patterns of type CompressedSimpleSparsityPattern.
      */
     template<typename SparsityType>
     void
@@ -1013,119 +814,88 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Returns the state of the sparsity
-     * pattern, i.e., whether compress()
-     * needs to be called after an
-     * operation requiring data
-     * exchange.
+     * Returns the state of the sparsity pattern, i.e., whether compress()
+     * needs to be called after an operation requiring data exchange.
      */
     bool is_compressed () const;
 
     /**
-     * Gives the maximum number of
-     * entries per row on the current
-     * processor.
+     * Gives the maximum number of entries per row on the current processor.
      */
     unsigned int max_entries_per_row () const;
 
     /**
-     * Return the number of rows in this
-     * sparsity pattern.
+     * Return the number of rows in this sparsity pattern.
      */
     size_type n_rows () const;
 
     /**
-     * Return the number of columns in
-     * this sparsity pattern.
+     * Return the number of columns in this sparsity pattern.
      */
     size_type n_cols () const;
 
     /**
-     * Return the local dimension of the
-     * sparsity pattern, i.e. the number
-     * of rows stored on the present MPI
-     * process. In the sequential case,
-     * this number is the same as
-     * n_rows(), but for parallel
-     * matrices it may be smaller.
+     * Return the local dimension of the sparsity pattern, i.e. the number of
+     * rows stored on the present MPI process. In the sequential case, this
+     * number is the same as n_rows(), but for parallel matrices 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().
      */
     unsigned int local_size () const;
 
     /**
-     * Return a pair of indices
-     * indicating which rows of this
-     * sparsity pattern are stored
-     * locally. The first number is the
-     * index of the first row stored, the
-     * second the index of the one past
-     * the last one that is stored
-     * locally. If this is a sequential
-     * matrix, then the result will be
-     * the pair (0,n_rows()), otherwise
-     * it will be a pair (i,i+n), where
+     * Return a pair of indices indicating which rows of this sparsity pattern
+     * are stored locally. The first number is the index of the first row
+     * stored, the second the index of the one past the last one that is
+     * stored locally. If this is a sequential matrix, then the result will be
+     * the pair (0,n_rows()), otherwise it will be a pair (i,i+n), where
      * <tt>n=local_size()</tt>.
      */
     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().
      */
     bool in_local_range (const size_type index) const;
 
     /**
-     * Return the number of nonzero
-     * elements of this sparsity pattern.
+     * Return the number of nonzero elements of this sparsity pattern.
      */
     size_type n_nonzero_elements () const;
 
     /**
-     * Number of entries in a
-     * specific row.
+     * Number of entries in a specific row.
      */
     size_type row_length (const size_type row) const;
 
     /**
-     * Compute the bandwidth of the
-     * matrix represented by this
-     * structure. The bandwidth is the
-     * maximum of $|i-j|$ for which the
-     * index pair $(i,j)$ represents a
-     * nonzero entry of the
-     * matrix. Consequently, the maximum
-     * bandwidth a $n\times m$ matrix can
-     * have is $\max\{n-1,m-1\}$.
+     * Compute the bandwidth of the matrix represented by this structure. The
+     * bandwidth is the maximum of $|i-j|$ for which the index pair $(i,j)$
+     * represents a nonzero entry of the matrix. Consequently, the maximum
+     * bandwidth a $n\times m$ matrix can have is $\max\{n-1,m-1\}$.
      */
     size_type bandwidth () const;
 
     /**
-     * Return whether the object is
-     * empty. It is empty if no memory is
-     * allocated, which is the same as
-     * when both dimensions are zero.
+     * Return whether the object is empty. It is empty if no memory is
+     * allocated, which is the same as when both dimensions are zero.
      */
     bool empty () const;
 
     /**
-     * Return whether the index
-     * (<i>i,j</i>) exists in the
-     * sparsity pattern (i.e., it may be
-     * non-zero) or not.
+     * Return whether the index (<i>i,j</i>) exists in the sparsity pattern
+     * (i.e., it may be non-zero) or not.
      */
     bool exists (const size_type i,
                  const size_type j) const;
 
     /**
-     * Determine an estimate for the
-     * memory consumption (in bytes)
-     * of this object. Currently not
-     * implemented for this class.
+     * Determine an estimate for the memory consumption (in bytes) of this
+     * object. Currently not implemented for this class.
      */
     std::size_t memory_consumption () const;
 
@@ -1135,16 +905,14 @@ namespace TrilinosWrappers
      */
 //@{
     /**
-     * Add the element (<i>i,j</i>) to
-     * the sparsity pattern.
+     * Add the element (<i>i,j</i>) to the sparsity pattern.
      */
     void add (const size_type i,
               const size_type j);
 
 
     /**
-     * Add several elements in one row to
-     * the sparsity pattern.
+     * Add several elements in one row to the sparsity pattern.
      */
     template <typename ForwardIterator>
     void add_entries (const size_type  row,
@@ -1158,61 +926,44 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos
-     * Epetra_CrsGraph data that stores
-     * the sparsity pattern.
+     * Return a const reference to the underlying Trilinos Epetra_CrsGraph
+     * data that stores the sparsity pattern.
      */
     const Epetra_FECrsGraph &trilinos_sparsity_pattern () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the parallel
-     * partitioning of the domain space
-     * of this sparsity pattern, i.e.,
-     * the partitioning of the vectors
-     * matrices based on this sparsity
-     * pattern are multiplied with.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the parallel partitioning of the domain space of this sparsity
+     * pattern, i.e., the partitioning of the vectors matrices based on this
+     * sparsity pattern are multiplied with.
      */
     const Epetra_Map &domain_partitioner () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * range space of this sparsity
-     * pattern, i.e., the partitioning of
-     * the vectors that are result from
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the range space of this sparsity pattern,
+     * i.e., the partitioning of the vectors that are result from
      * matrix-vector products.
      */
     const Epetra_Map &range_partitioner () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * sparsity pattern rows. Equal to
-     * the partitioning of the range.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the sparsity pattern rows. Equal to the
+     * partitioning of the range.
      */
     const Epetra_Map &row_partitioner () const;
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the partitioning of the
-     * sparsity pattern columns. This is
-     * in general not equal to the
-     * partitioner Epetra_Map for the
-     * domain because of overlap in the
-     * matrix.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the partitioning of the sparsity pattern columns. This is in
+     * general not equal to the partitioner Epetra_Map for the domain because
+     * of overlap in the matrix.
      */
     const Epetra_Map &col_partitioner () const;
 
     /**
-     * Return a const reference to
-     * the communicator used for
-     * this object.
+     * Return a const reference to the communicator used for this object.
      */
     const Epetra_Comm &trilinos_communicator () const;
 //@}
@@ -1222,8 +973,7 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * STL-like iterator with the
-     * first entry.
+     * STL-like iterator with the first entry.
      */
     const_iterator begin () const;
 
@@ -1233,34 +983,22 @@ namespace TrilinosWrappers
     const_iterator end () const;
 
     /**
-     * STL-like iterator with the
-     * first entry of row @p r.
+     * STL-like iterator with the first entry of row @p r.
      *
-     * Note that if the given row
-     * is empty, i.e. does not
-     * contain any nonzero entries,
-     * then the iterator returned
-     * by this function equals
-     * <tt>end(r)</tt>. Note also
-     * that the iterator may not be
-     * dereferencable in that case.
+     * Note that if the given row is empty, i.e. does not contain any nonzero
+     * entries, then the iterator returned by this function equals
+     * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable
+     * in that case.
      */
     const_iterator begin (const size_type r) const;
 
     /**
-     * Final iterator of row
-     * <tt>r</tt>. It points to the
-     * first element past the end
-     * of line @p r, or past the
-     * end of the entire sparsity
-     * pattern.
+     * Final iterator of row <tt>r</tt>. It points to the first element past
+     * the end of line @p r, or past the end of the entire sparsity pattern.
      *
-     * Note that the end iterator
-     * is not necessarily
-     * dereferencable. This is in
-     * particular the case if it is
-     * the end iterator for the
-     * last row of a matrix.
+     * Note that the end iterator is not necessarily dereferencable. This is
+     * in particular the case if it is the end iterator for the last row of a
+     * matrix.
      */
     const_iterator end (const size_type r) const;
 
@@ -1271,55 +1009,35 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Abstract Trilinos object
-     * that helps view in ASCII
-     * other Trilinos
-     * objects. Currently this
-     * function is not
-     * implemented.  TODO: Not
+     * Abstract Trilinos object that helps view in ASCII other Trilinos
+     * objects. Currently this function is not implemented.  TODO: Not
      * implemented.
      */
     void write_ascii ();
 
     /**
-     * Print (the locally owned part of)
-     * the sparsity pattern to the given
-     * stream, using the format
-     * <tt>(line,col)</tt>. The optional
-     * flag outputs the sparsity pattern
-     * in Trilinos style, where even the
-     * according processor number is
-     * printed to the stream, as well as
-     * a summary before actually writing
-     * the entries.
+     * Print (the locally owned part of) the sparsity pattern to the given
+     * stream, using the format <tt>(line,col)</tt>. The optional flag outputs
+     * the sparsity pattern in Trilinos style, where even the according
+     * processor number is printed to the stream, as well as a summary before
+     * actually writing the entries.
      */
     void print (std::ostream &out,
                 const bool    write_extended_trilinos_info = false) const;
 
     /**
-     * Print the sparsity of the matrix
-     * in a format that <tt>gnuplot</tt>
-     * understands and which can be used
-     * to plot the sparsity pattern in a
-     * graphical way. The format consists
-     * of pairs <tt>i j</tt> of nonzero
-     * elements, each representing one
-     * entry of this matrix, one per line
-     * of the output file. Indices are
-     * counted from zero on, as
-     * usual. Since sparsity patterns are
-     * printed in the same way as
-     * matrices are displayed, we print
-     * the negative of the column index,
-     * which means that the
-     * <tt>(0,0)</tt> element is in the
-     * top left rather than in the bottom
+     * Print the sparsity of the matrix in a format that <tt>gnuplot</tt>
+     * understands and which can be used to plot the sparsity pattern in a
+     * graphical way. The format consists of pairs <tt>i j</tt> of nonzero
+     * elements, each representing one entry of this matrix, one per line of
+     * the output file. Indices are counted from zero on, as usual. Since
+     * sparsity patterns are printed in the same way as matrices are
+     * displayed, we print the negative of the column index, which means that
+     * the <tt>(0,0)</tt> element is in the top left rather than in the bottom
      * left corner.
      *
-     * Print the sparsity pattern in
-     * gnuplot by setting the data style
-     * to dots or points and use the
-     * <tt>plot</tt> command.
+     * Print the sparsity pattern in gnuplot by setting the data style to dots
+     * or points and use the <tt>plot</tt> command.
      */
     void print_gnuplot (std::ostream &out) const;
 
@@ -1371,30 +1089,31 @@ namespace TrilinosWrappers
   private:
 
     /**
-     * Pointer to the user-supplied
-     * Epetra Trilinos mapping of
-     * the matrix columns that
-     * assigns parts of the matrix
-     * to the individual processes.
+     * Pointer to the user-supplied Epetra Trilinos mapping of the matrix
+     * columns that assigns parts of the matrix to the individual processes.
      */
     std_cxx1x::shared_ptr<Epetra_Map> column_space_map;
 
     /**
-     * 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;
 
     /**
-     * A sparsity pattern object in
-     * Trilinos to be used for finite
-     * element based problems which
-     * allows for adding non-local
-     * elements to the pattern.
+     * A sparsity pattern object in Trilinos to be used for finite element
+     * based problems which allows for adding non-local elements to the
+     * pattern.
      */
     std_cxx1x::shared_ptr<Epetra_FECrsGraph> graph;
 
+    /**
+     * A sparsity pattern object for the non-local part of the sparsity
+     * pattern that is going to be sent to the owning processor. Only used when the particular constructor or reinit method with writable_rows argument is set
+     */
+    std_cxx1x::shared_ptr<Epetra_CrsGraph> nonlocal_graph;
+
+    friend class SparseMatrix;
     friend class SparsityPatternIterators::Accessor;
     friend class SparsityPatternIterators::Iterator;
   };
@@ -1696,9 +1415,23 @@ namespace TrilinosWrappers
     const int n_cols = static_cast<int>(end - begin);
     compressed = false;
 
-    const int ierr = graph->InsertGlobalIndices (1,
-                                                 (TrilinosWrappers::types::int_type *)&row,
-                                                 n_cols, col_index_ptr);
+    int ierr;
+    if ( graph->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row)) != -1)
+      ierr = graph->InsertGlobalIndices (row, n_cols, col_index_ptr);
+    else if (nonlocal_graph.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_graph->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_graph->InsertGlobalIndices (row, n_cols, col_index_ptr);
+      }
+    else
+      ierr = graph->InsertGlobalIndices
+        (1, (TrilinosWrappers::types::int_type *)&row, n_cols, col_index_ptr);
 
     AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
   }
@@ -1757,157 +1490,6 @@ namespace TrilinosWrappers
     return graph->RangeMap().Comm();
   }
 
-
-
-  inline
-  SparsityPattern::SparsityPattern  (const IndexSet  &parallel_partitioning,
-                                     const MPI_Comm  &communicator,
-                                     const size_type  n_entries_per_row)
-    :
-    compressed (false)
-  {
-    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
-                                                              false);
-    reinit (map, map, n_entries_per_row);
-  }
-
-
-
-  inline
-  SparsityPattern::SparsityPattern  (const IndexSet     &parallel_partitioning,
-                                     const MPI_Comm     &communicator,
-                                     const std::vector<size_type> &n_entries_per_row)
-    :
-    compressed (false)
-  {
-    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
-                                                              false);
-    reinit (map, map, n_entries_per_row);
-  }
-
-
-
-  inline
-  SparsityPattern::SparsityPattern  (const IndexSet  &row_parallel_partitioning,
-                                     const IndexSet  &col_parallel_partitioning,
-                                     const MPI_Comm  &communicator,
-                                     const size_type  n_entries_per_row)
-    :
-    compressed (false)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, n_entries_per_row);
-  }
-
-
-
-  inline
-  SparsityPattern::
-  SparsityPattern  (const IndexSet     &row_parallel_partitioning,
-                    const IndexSet     &col_parallel_partitioning,
-                    const MPI_Comm     &communicator,
-                    const std::vector<size_type> &n_entries_per_row)
-    :
-    compressed (false)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, n_entries_per_row);
-  }
-
-
-
-  inline
-  void
-  SparsityPattern::reinit (const IndexSet  &parallel_partitioning,
-                           const MPI_Comm  &communicator,
-                           const size_type  n_entries_per_row)
-  {
-    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
-                                                              false);
-    reinit (map, map, n_entries_per_row);
-  }
-
-
-
-  inline
-  void SparsityPattern::reinit (const IndexSet     &parallel_partitioning,
-                                const MPI_Comm     &communicator,
-                                const std::vector<size_type> &n_entries_per_row)
-  {
-    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
-                                                              false);
-    reinit (map, map, n_entries_per_row);
-  }
-
-
-
-  inline
-  void SparsityPattern::reinit (const IndexSet &row_parallel_partitioning,
-                                const IndexSet &col_parallel_partitioning,
-                                const MPI_Comm &communicator,
-                                const size_type  n_entries_per_row)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, n_entries_per_row);
-  }
-
-
-  inline
-  void
-  SparsityPattern::reinit (const IndexSet     &row_parallel_partitioning,
-                           const IndexSet     &col_parallel_partitioning,
-                           const MPI_Comm     &communicator,
-                           const std::vector<size_type> &n_entries_per_row)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, n_entries_per_row);
-  }
-
-
-
-  template<typename SparsityType>
-  inline
-  void
-  SparsityPattern::reinit (const IndexSet     &row_parallel_partitioning,
-                           const IndexSet     &col_parallel_partitioning,
-                           const SparsityType &nontrilinos_sparsity_pattern,
-                           const MPI_Comm     &communicator,
-                           const bool          exchange_data)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, nontrilinos_sparsity_pattern, exchange_data);
-  }
-
-
-
-  template<typename SparsityType>
-  inline
-  void
-  SparsityPattern::reinit (const IndexSet     &parallel_partitioning,
-                           const SparsityType &nontrilinos_sparsity_pattern,
-                           const MPI_Comm     &communicator,
-                           const bool          exchange_data)
-  {
-    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
-                                                              false);
-    reinit (map, map, nontrilinos_sparsity_pattern, exchange_data);
-  }
-
 #endif // DOXYGEN
 }
 
index 695ae57da6c659fd2a678ba52f93987b5cbcfb7f..0536414ef5f7da6492603d0f392e533685b261c9 100644 (file)
@@ -785,6 +785,45 @@ namespace TrilinosWrappers
     this->collect_sizes();
   }
 
+
+
+  void
+  BlockSparsityPattern::reinit (const std::vector<IndexSet> &row_parallel_partitioning,
+                                const std::vector<IndexSet> &col_parallel_partitioning,
+                                const MPI_Comm &communicator)
+  {
+    dealii::BlockSparsityPatternBase<SparsityPattern>::
+    reinit(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],
+                                communicator);
+    this->collect_sizes();
+  }
+
+
+
+  void
+  BlockSparsityPattern::reinit (const std::vector<IndexSet> &row_parallel_partitioning,
+                                const std::vector<IndexSet> &col_parallel_partitioning,
+                                const std::vector<IndexSet> &writable_rows,
+                                const MPI_Comm &communicator)
+  {
+    AssertDimension(writable_rows.size(), row_parallel_partitioning.size());
+    dealii::BlockSparsityPatternBase<SparsityPattern>::
+    reinit(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();
+  }
+
 }
 
 #endif
index b8f2956a7449951f5d7c4577eab2d6774f88bf1e..665c7ba1513884dfa8a454e23b2a50f38a0e406a 100644 (file)
@@ -26,6 +26,7 @@
 #  include <deal.II/lac/compressed_set_sparsity_pattern.h>
 #  include <deal.II/lac/compressed_simple_sparsity_pattern.h>
 
+#  include <Epetra_Export.h>
 #  include <ml_epetra_utils.h>
 #  include <ml_struct.h>
 #  include <Teuchos_RCP.hpp>
@@ -574,13 +575,18 @@ namespace TrilinosWrappers
   {
     matrix.reset ();
 
-    // reinit with a (parallel) Trilinos
-    // sparsity pattern.
+    // reinit with a (parallel) Trilinos sparsity pattern.
     column_space_map.reset (new Epetra_Map
                             (sparsity_pattern.domain_partitioner()));
     matrix.reset (new Epetra_FECrsMatrix
                   (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));
+      }
     compress();
+    last_action = Zero;
   }
 
 
@@ -770,6 +776,56 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  void
+  SparseMatrix::compress (::dealii::VectorOperation::values operation)
+  {
+
+    Epetra_CombineMode mode = last_action;
+    if (last_action == Zero)
+      {
+        if ((operation==::dealii::VectorOperation::add) ||
+            (operation==::dealii::VectorOperation::unknown))
+          mode = Add;
+        else if (operation==::dealii::VectorOperation::insert)
+          mode = Insert;
+      }
+    else
+      {
+        Assert(
+          ((last_action == Add) && (operation!=::dealii::VectorOperation::insert))
+          ||
+          ((last_action == Insert) && (operation!=::dealii::VectorOperation::add)),
+          ExcMessage("operation and argument to compress() do not match"));
+      }
+
+    // flush buffers
+    int ierr;
+    if (nonlocal_matrix.get() != 0)
+      {
+        nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap());
+        Epetra_Export exporter(nonlocal_matrix->RowMap(), matrix->RowMap());
+        ierr = matrix->Export(*nonlocal_matrix, exporter, mode);
+        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
+        ierr = matrix->FillComplete(*column_space_map, matrix->RowMap());
+        nonlocal_matrix->PutScalar(0);
+      }
+    else
+      ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(),
+                                     true, mode);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = matrix->OptimizeStorage ();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    last_action = Zero;
+
+    compressed = true;
+  }
+
+
+
   void
   SparseMatrix::clear ()
   {
@@ -779,6 +835,7 @@ namespace TrilinosWrappers
     column_space_map.reset (new Epetra_Map (0, 0,
                                             Utilities::Trilinos::comm_self()));
     matrix.reset (new Epetra_FECrsMatrix(View, *column_space_map, 0));
+    nonlocal_matrix.reset();
 
     matrix->FillComplete();
 
index 24f46180920d18e70bafe2660aa59e08a80f01ab..85ca1111fcb3981a5710ef7cf35322dcec16e16e 100644 (file)
@@ -24,6 +24,8 @@
 #  include <deal.II/lac/compressed_set_sparsity_pattern.h>
 #  include <deal.II/lac/compressed_simple_sparsity_pattern.h>
 
+#  include <Epetra_Export.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace TrilinosWrappers
@@ -250,6 +252,80 @@ namespace TrilinosWrappers
 
 
 
+  SparsityPattern::SparsityPattern  (const IndexSet  &parallel_partitioning,
+                                     const MPI_Comm  &communicator,
+                                     const size_type  n_entries_per_row)
+    :
+    compressed (false)
+  {
+    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
+                                                              false);
+    reinit (map, map, n_entries_per_row);
+  }
+
+
+
+  SparsityPattern::SparsityPattern  (const IndexSet     &parallel_partitioning,
+                                     const MPI_Comm     &communicator,
+                                     const std::vector<size_type> &n_entries_per_row)
+    :
+    compressed (false)
+  {
+    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
+                                                              false);
+    reinit (map, map, n_entries_per_row);
+  }
+
+
+
+  SparsityPattern::SparsityPattern  (const IndexSet  &row_parallel_partitioning,
+                                     const IndexSet  &col_parallel_partitioning,
+                                     const MPI_Comm  &communicator,
+                                     const size_type  n_entries_per_row)
+    :
+    compressed (false)
+  {
+    Epetra_Map row_map =
+      row_parallel_partitioning.make_trilinos_map (communicator, false);
+    Epetra_Map col_map =
+      col_parallel_partitioning.make_trilinos_map (communicator, false);
+    reinit (row_map, col_map, n_entries_per_row);
+  }
+
+
+
+  SparsityPattern::
+  SparsityPattern  (const IndexSet     &row_parallel_partitioning,
+                    const IndexSet     &col_parallel_partitioning,
+                    const MPI_Comm     &communicator,
+                    const std::vector<size_type> &n_entries_per_row)
+    :
+    compressed (false)
+  {
+    Epetra_Map row_map =
+      row_parallel_partitioning.make_trilinos_map (communicator, false);
+    Epetra_Map col_map =
+      col_parallel_partitioning.make_trilinos_map (communicator, false);
+    reinit (row_map, col_map, n_entries_per_row);
+  }
+
+
+
+  SparsityPattern::
+  SparsityPattern  (const IndexSet     &row_parallel_partitioning,
+                    const IndexSet     &col_parallel_partitioning,
+                    const IndexSet     &writable_rows,
+                    const MPI_Comm     &communicator,
+                    const size_type     n_max_entries_per_row)
+    :
+    compressed (false)
+  {
+    reinit (row_parallel_partitioning, col_parallel_partitioning,
+            writable_rows, communicator, n_max_entries_per_row);
+  }
+
+
+
   SparsityPattern::~SparsityPattern ()
   {}
 
@@ -282,6 +358,12 @@ namespace TrilinosWrappers
                            const Epetra_Map  &input_col_map,
                            const size_type    n_entries_per_row)
   {
+    Assert(input_row_map.IsOneToOne(),
+           ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
+                      "the maps of different processors."));
+    Assert(input_col_map.IsOneToOne(),
+           ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
+                      "the maps of different processors."));
     graph.reset ();
     column_space_map.reset (new Epetra_Map (input_col_map));
     compressed = false;
@@ -355,6 +437,122 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparsityPattern::reinit (const IndexSet  &parallel_partitioning,
+                           const MPI_Comm  &communicator,
+                           const size_type  n_entries_per_row)
+  {
+    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
+                                                              false);
+    reinit (map, map, n_entries_per_row);
+  }
+
+
+
+  void SparsityPattern::reinit (const IndexSet     &parallel_partitioning,
+                                const MPI_Comm     &communicator,
+                                const std::vector<size_type> &n_entries_per_row)
+  {
+    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
+                                                              false);
+    reinit (map, map, n_entries_per_row);
+  }
+
+
+
+  void SparsityPattern::reinit (const IndexSet &row_parallel_partitioning,
+                                const IndexSet &col_parallel_partitioning,
+                                const MPI_Comm &communicator,
+                                const size_type  n_entries_per_row)
+  {
+    Epetra_Map row_map =
+      row_parallel_partitioning.make_trilinos_map (communicator, false);
+    Epetra_Map col_map =
+      col_parallel_partitioning.make_trilinos_map (communicator, false);
+    reinit (row_map, col_map, n_entries_per_row);
+  }
+
+
+
+  void
+  SparsityPattern::reinit (const IndexSet     &row_parallel_partitioning,
+                           const IndexSet     &col_parallel_partitioning,
+                           const MPI_Comm     &communicator,
+                           const std::vector<size_type> &n_entries_per_row)
+  {
+    Epetra_Map row_map =
+      row_parallel_partitioning.make_trilinos_map (communicator, false);
+    Epetra_Map col_map =
+      col_parallel_partitioning.make_trilinos_map (communicator, false);
+    reinit (row_map, col_map, n_entries_per_row);
+  }
+
+
+
+  template<typename SparsityType>
+  void
+  SparsityPattern::reinit (const IndexSet     &row_parallel_partitioning,
+                           const IndexSet     &col_parallel_partitioning,
+                           const SparsityType &nontrilinos_sparsity_pattern,
+                           const MPI_Comm     &communicator,
+                           const bool          exchange_data)
+  {
+    Epetra_Map row_map =
+      row_parallel_partitioning.make_trilinos_map (communicator, false);
+    Epetra_Map col_map =
+      col_parallel_partitioning.make_trilinos_map (communicator, false);
+    reinit (row_map, col_map, nontrilinos_sparsity_pattern, exchange_data);
+  }
+
+
+
+  template<typename SparsityType>
+  void
+  SparsityPattern::reinit (const IndexSet     &parallel_partitioning,
+                           const SparsityType &nontrilinos_sparsity_pattern,
+                           const MPI_Comm     &communicator,
+                           const bool          exchange_data)
+  {
+    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator,
+                                                              false);
+    reinit (map, map, nontrilinos_sparsity_pattern, exchange_data);
+  }
+
+
+
+  void
+  SparsityPattern::reinit (const IndexSet  &row_parallel_partitioning,
+                           const IndexSet  &col_parallel_partitioning,
+                           const IndexSet  &writable_rows,
+                           const MPI_Comm  &communicator,
+                           const size_type  n_entries_per_row)
+  {
+    reinit(row_parallel_partitioning, col_parallel_partitioning,
+           communicator,n_entries_per_row);
+
+    IndexSet nonlocal_partitioner = writable_rows;
+    AssertDimension(nonlocal_partitioner.size(), row_parallel_partitioning.size());
+#ifdef DEBUG
+    {
+      IndexSet tmp = writable_rows & row_parallel_partitioning;
+      Assert (tmp == row_parallel_partitioning,
+              ExcMessage("The set of writable rows passed to this method does not "
+                         "contain the locally owned rows, which is not allowed."));
+    }
+#endif
+    nonlocal_partitioner.subtract_set(row_parallel_partitioning);
+    if (Utilities::MPI::n_mpi_processes(communicator) > 1)
+      {
+        Epetra_Map nonlocal_map =
+          nonlocal_partitioner.make_trilinos_map(communicator, false);
+        nonlocal_graph.reset(new Epetra_CrsGraph(Copy, nonlocal_map, 0));
+      }
+    else
+      Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
+  }
+
+
+
   template <typename SparsityType>
   void
   SparsityPattern::reinit (const Epetra_Map   &input_map,
@@ -509,6 +707,8 @@ namespace TrilinosWrappers
                                        *column_space_map, 0));
     graph->FillComplete();
 
+    nonlocal_graph.reset();
+
     compressed = true;
   }
 
@@ -519,9 +719,31 @@ namespace TrilinosWrappers
   {
     int ierr;
     Assert (column_space_map.get() != 0, ExcInternalError());
-    ierr = graph->GlobalAssemble (*column_space_map,
-                                  static_cast<const Epetra_Map &>(graph->RangeMap()),
-                                  true);
+    if (nonlocal_graph.get() != 0)
+      {
+        if (nonlocal_graph->IndicesAreGlobal() == false &&
+            nonlocal_graph->RowMap().NumMyElements() > 0)
+          {
+            // insert dummy element
+            TrilinosWrappers::types::int_type row = nonlocal_graph->RowMap().MyGID(0);
+            nonlocal_graph->InsertGlobalIndices(row, 1, &row);
+          }
+        Assert(nonlocal_graph->IndicesAreGlobal() == true,
+               ExcInternalError());
+        nonlocal_graph->FillComplete(*column_space_map,
+                                     static_cast<const Epetra_Map &>(graph->RangeMap()));
+        nonlocal_graph->OptimizeStorage();
+        Epetra_Export exporter(nonlocal_graph->RowMap(), graph->RowMap());
+        ierr = graph->Export(*nonlocal_graph, exporter, Add);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+        ierr =
+          graph->FillComplete(*column_space_map,
+                              static_cast<const Epetra_Map &>(graph->RangeMap()));
+      }
+    else
+      ierr = graph->GlobalAssemble (*column_space_map,
+                                    static_cast<const Epetra_Map &>(graph->RangeMap()),
+                                    true);
 
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
@@ -849,6 +1071,55 @@ namespace TrilinosWrappers
                            const dealii::CompressedSimpleSparsityPattern &,
                            bool);
 
+
+
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const dealii::SparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const dealii::CompressedSparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const dealii::CompressedSetSparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const dealii::CompressedSimpleSparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+
+
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const IndexSet &,
+                           const dealii::SparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const IndexSet &,
+                           const dealii::CompressedSparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const IndexSet &,
+                           const dealii::CompressedSetSparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+  template void
+  SparsityPattern::reinit (const IndexSet &,
+                           const IndexSet &,
+                           const dealii::CompressedSimpleSparsityPattern &,
+                           const MPI_Comm &,
+                           bool);
+
 }
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/trilinos/assemble_matrix_parallel_01.cc b/tests/trilinos/assemble_matrix_parallel_01.cc
new file mode 100644 (file)
index 0000000..209fe9c
--- /dev/null
@@ -0,0 +1,451 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same as deal.II/assemble_matrix_parallel_01, but for trilinos matrices
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/graph_coloring.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim>
+    struct Data
+    {
+      Data (const hp::FECollection<dim> &fe,
+            const hp::QCollection<dim>  &quadrature)
+        :
+        hp_fe_values(fe,
+                     quadrature,
+                     update_values    |  update_gradients |
+                     update_quadrature_points  |  update_JxW_values)
+      {}
+
+      Data (const Data &data)
+        :
+        hp_fe_values(data.hp_fe_values.get_mapping_collection(),
+                     data.hp_fe_values.get_fe_collection(),
+                     data.hp_fe_values.get_quadrature_collection(),
+                     data.hp_fe_values.get_update_flags())
+      {}
+
+      hp::FEValues<dim>               hp_fe_values;
+    };
+  }
+
+  namespace Copy
+  {
+    struct Data
+    {
+      std::vector<types::global_dof_index> local_dof_indices;
+      FullMatrix<double> local_matrix;
+      Vector<double> local_rhs;
+    };
+  }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem ();
+  ~LaplaceProblem ();
+
+  void run ();
+
+private:
+  void setup_system ();
+  void test_equality ();
+  void assemble_reference ();
+  void assemble_test ();
+  void solve ();
+  void create_coarse_grid ();
+  void postprocess ();
+
+  void local_assemble (const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                       Assembly::Scratch::Data<dim>  &scratch,
+                       Assembly::Copy::Data          &data);
+  void copy_local_to_global (const Assembly::Copy::Data &data);
+
+  std::vector<types::global_dof_index>
+  get_conflict_indices (typename hp::DoFHandler<dim>::active_cell_iterator const &cell) const;
+
+  Triangulation<dim>   triangulation;
+
+  hp::DoFHandler<dim>      dof_handler;
+  hp::FECollection<dim>    fe_collection;
+  hp::QCollection<dim>     quadrature_collection;
+  hp::QCollection<dim-1>   face_quadrature_collection;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::SparseMatrix reference_matrix;
+  TrilinosWrappers::SparseMatrix test_matrix;
+
+  Vector<double>       reference_rhs;
+  Vector<double>       test_rhs;
+
+  std::vector<std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> > graph;
+
+  const unsigned int max_degree;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                            const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                           const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+  :
+  dof_handler (triangulation),
+  max_degree (5)
+{
+  if (dim == 2)
+    for (unsigned int degree=2; degree<=max_degree; ++degree)
+      {
+        fe_collection.push_back (FE_Q<dim>(degree));
+        quadrature_collection.push_back (QGauss<dim>(degree+1));
+        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+  else
+    for (unsigned int degree=1; degree<max_degree-1; ++degree)
+      {
+        fe_collection.push_back (FE_Q<dim>(degree));
+        quadrature_collection.push_back (QGauss<dim>(degree+1));
+        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (typename hp::DoFHandler<dim>::active_cell_iterator const &cell) const
+{
+  std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+  cell->get_dof_indices(local_dof_indices);
+
+  constraints.resolve_indices(local_dof_indices);
+  return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  reference_matrix.clear();
+  test_matrix.clear();
+  dof_handler.distribute_dofs (fe_collection);
+
+  reference_rhs.reinit (dof_handler.n_dofs());
+  test_rhs.reinit (dof_handler.n_dofs());
+
+  constraints.clear ();
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+  // add boundary conditions as inhomogeneous constraints here, do it after
+  // having added the hanging node constraints in order to be consistent and
+  // skip dofs that are already constrained (i.e., are hanging nodes on the
+  // boundary in 3D). In contrast to step-27, we choose a sine function.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  graph = GraphColoring::make_graph_coloring(dof_handler.begin_active(),dof_handler.end(),
+                                             static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+                                             (typename hp::DoFHandler<dim>::active_cell_iterator const &)> >
+                                             (std_cxx1x::bind(&LaplaceProblem<dim>::get_conflict_indices, this,std_cxx1x::_1)));
+
+
+  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+                                       dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                   constraints, false);
+  reference_matrix.reinit (dof_handler.locally_owned_dofs(), csp, MPI_COMM_WORLD);
+  test_matrix.reinit (reference_matrix);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                                     Assembly::Scratch::Data<dim>  &scratch,
+                                     Assembly::Copy::Data          &data)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+  data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+  data.local_matrix = 0;
+
+  data.local_rhs.reinit (dofs_per_cell);
+  data.local_rhs = 0;
+
+  scratch.hp_fe_values.reinit (cell);
+
+  const FEValues<dim> &fe_values = scratch.hp_fe_values.get_present_fe_values ();
+
+  const RightHandSide<dim> rhs_function;
+
+  for (unsigned int q_point=0;
+       q_point<fe_values.n_quadrature_points;
+       ++q_point)
+    {
+      const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            data.local_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                       fe_values.shape_grad(j,q_point) *
+                                       fe_values.JxW(q_point));
+
+            data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                  rhs_value *
+                                  fe_values.JxW(q_point));
+        }
+    }
+
+  data.local_dof_indices.resize (dofs_per_cell);
+  cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+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);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  Assembly::Copy::Data copy_data;
+  Assembly::Scratch::Data<dim> assembly_data(fe_collection, quadrature_collection);
+
+  for (unsigned int color=0; color<graph.size(); ++color)
+    for (typename std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::const_iterator p = graph[color].begin();
+         p != graph[color].end(); ++p)
+      {
+        local_assemble(*p, assembly_data, copy_data);
+        copy_local_to_global(copy_data);
+      }
+
+  reference_matrix.add(1., test_matrix);
+  reference_rhs = test_rhs;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe_collection, quadrature_collection),
+         Assembly::Copy::Data ());
+
+  test_matrix.add(-1, reference_matrix);
+
+  // there should not even be roundoff difference between matrices
+  deallog.threshold_double(1.e-30);
+  deallog << "error in matrix: " << test_matrix.frobenius_norm() << std::endl;
+  test_rhs.add(-1., reference_rhs);
+  deallog << "error in vector: " << test_rhs.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+    estimated_error_per_cell(i) = i;
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+
+  for (typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+    cell->set_active_fe_index (rand() % fe_collection.size());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube(triangulation, 0, 1/* ,
+                                     Point<dim>(),
+                                     0.5, 1., (dim==3) ? 96 : 12, false*/);
+          triangulation.refine_global(2);
+        }
+
+      setup_system ();
+
+      assemble_reference ();
+      assemble_test ();
+
+      if (cycle < 2)
+        postprocess ();
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int);
+
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+
+  {
+    deallog.push("3d");
+    LaplaceProblem<3> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/trilinos/assemble_matrix_parallel_01.output b/tests/trilinos/assemble_matrix_parallel_01.output
new file mode 100644 (file)
index 0000000..3b8bf47
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:3d::error in matrix: 0
+DEAL:3d::error in vector: 0
+DEAL:3d::error in matrix: 0
+DEAL:3d::error in vector: 0
+DEAL:3d::error in matrix: 0
+DEAL:3d::error in vector: 0
diff --git a/tests/trilinos/assemble_matrix_parallel_02.cc b/tests/trilinos/assemble_matrix_parallel_02.cc
new file mode 100644 (file)
index 0000000..5dfdc8a
--- /dev/null
@@ -0,0 +1,437 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests parallel assembly of Trilinos matrices when also using MPI
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/graph_coloring.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim>
+    struct Data
+    {
+      Data (const FiniteElement<dim> &fe,
+            const Quadrature<dim>    &quadrature)
+        :
+        fe_values(fe,
+                  quadrature,
+                  update_values    |  update_gradients |
+                  update_quadrature_points  |  update_JxW_values)
+      {}
+
+      Data (const Data &data)
+        :
+        fe_values(data.fe_values.get_mapping(),
+                  data.fe_values.get_fe(),
+                  data.fe_values.get_quadrature(),
+                  data.fe_values.get_update_flags())
+      {}
+
+      FEValues<dim> fe_values;
+    };
+  }
+
+  namespace Copy
+  {
+    struct Data
+    {
+      std::vector<types::global_dof_index> local_dof_indices;
+      FullMatrix<double> local_matrix;
+      Vector<double> local_rhs;
+    };
+  }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem ();
+  ~LaplaceProblem ();
+
+  void run ();
+
+private:
+  void setup_system ();
+  void test_equality ();
+  void assemble_reference ();
+  void assemble_test ();
+  void solve ();
+  void create_coarse_grid ();
+  void postprocess ();
+
+  void local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+                       Assembly::Scratch::Data<dim>  &scratch,
+                       Assembly::Copy::Data          &data);
+  void copy_local_to_global (const Assembly::Copy::Data &data);
+
+  std::vector<types::global_dof_index>
+  get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const;
+
+  parallel::distributed::Triangulation<dim>   triangulation;
+
+  DoFHandler<dim>      dof_handler;
+  FE_Q<dim>            fe;
+  QGauss<dim>          quadrature;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::SparseMatrix reference_matrix;
+  TrilinosWrappers::SparseMatrix test_matrix;
+
+  TrilinosWrappers::MPI::Vector       reference_rhs;
+  TrilinosWrappers::MPI::Vector       test_rhs;
+
+  std::vector<std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > graph;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                            const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                           const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+  :
+  triangulation (MPI_COMM_WORLD),
+  dof_handler (triangulation),
+  fe (1),
+  quadrature(fe.degree+1)
+{
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const
+{
+  std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+  cell->get_dof_indices(local_dof_indices);
+
+  constraints.resolve_indices(local_dof_indices);
+  return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  reference_matrix.clear();
+  test_matrix.clear();
+  dof_handler.distribute_dofs (fe);
+
+  constraints.clear ();
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+  // add boundary conditions as inhomogeneous constraints here, do it after
+  // having added the hanging node constraints in order to be consistent and
+  // skip dofs that are already constrained (i.e., are hanging nodes on the
+  // boundary in 3D). In contrast to step-27, we choose a sine function.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  typedef FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> CellFilter;
+  CellFilter begin(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active());
+  CellFilter end(IteratorFilters::LocallyOwnedCell(),dof_handler.end());
+  graph = GraphColoring::make_graph_coloring(begin,end,
+        static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+        (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);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+                                     Assembly::Scratch::Data<dim>  &scratch,
+                                     Assembly::Copy::Data          &data)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+  data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+  data.local_matrix = 0;
+
+  data.local_rhs.reinit (dofs_per_cell);
+  data.local_rhs = 0;
+
+  scratch.fe_values.reinit (cell);
+
+  const FEValues<dim> &fe_values = scratch.fe_values;
+
+  const RightHandSide<dim> rhs_function;
+
+  for (unsigned int q_point=0;
+       q_point<fe_values.n_quadrature_points;
+       ++q_point)
+    {
+      const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            data.local_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                       fe_values.shape_grad(j,q_point) *
+                                       fe_values.JxW(q_point));
+
+            data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                  rhs_value *
+                                  fe_values.JxW(q_point));
+        }
+    }
+
+  data.local_dof_indices.resize (dofs_per_cell);
+  cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+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);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  Assembly::Copy::Data copy_data;
+  Assembly::Scratch::Data<dim> assembly_data(fe, quadrature);
+
+  for (unsigned int color=0; color<graph.size(); ++color)
+    for (typename std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >::const_iterator p = graph[color].begin();
+         p != graph[color].end(); ++p)
+      {
+        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;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe, quadrature),
+         Assembly::Copy::Data (),
+         1);
+  test_matrix.compress(VectorOperation::add);
+  test_rhs.compress(VectorOperation::add);
+
+  test_matrix.add(-1, reference_matrix);
+
+  // there should not even be roundoff difference between matrices
+  deallog.threshold_double(1.e-30);
+  deallog << "error in matrix: " << test_matrix.frobenius_norm() << std::endl;
+  test_rhs.add(-1., reference_rhs);
+  deallog << "error in vector: " << test_rhs.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+    estimated_error_per_cell(i) = i;
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube(triangulation, 0, 1/* ,
+                                     Point<dim>(),
+                                     0.5, 1., (dim==3) ? 96 : 12, false*/);
+          triangulation.refine_global(6);
+        }
+
+      setup_system ();
+
+      assemble_reference ();
+      assemble_test ();
+
+      if (cycle < 2)
+        postprocess ();
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int);
+
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output b/tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output
new file mode 100644 (file)
index 0000000..6e5eabc
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
diff --git a/tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output b/tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output
new file mode 100644 (file)
index 0000000..6e5eabc
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
diff --git a/tests/trilinos/assemble_matrix_parallel_03.cc b/tests/trilinos/assemble_matrix_parallel_03.cc
new file mode 100644 (file)
index 0000000..65a82e1
--- /dev/null
@@ -0,0 +1,440 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Same test as assemble_matrix_parallel_02, but calling compress() several
+// times on test matrix
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/graph_coloring.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim>
+    struct Data
+    {
+      Data (const FiniteElement<dim> &fe,
+            const Quadrature<dim>    &quadrature)
+        :
+        fe_values(fe,
+                  quadrature,
+                  update_values    |  update_gradients |
+                  update_quadrature_points  |  update_JxW_values)
+      {}
+
+      Data (const Data &data)
+        :
+        fe_values(data.fe_values.get_mapping(),
+                  data.fe_values.get_fe(),
+                  data.fe_values.get_quadrature(),
+                  data.fe_values.get_update_flags())
+      {}
+
+      FEValues<dim> fe_values;
+    };
+  }
+
+  namespace Copy
+  {
+    struct Data
+    {
+      std::vector<types::global_dof_index> local_dof_indices;
+      FullMatrix<double> local_matrix;
+      Vector<double> local_rhs;
+    };
+  }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem ();
+  ~LaplaceProblem ();
+
+  void run ();
+
+private:
+  void setup_system ();
+  void test_equality ();
+  void assemble_reference ();
+  void assemble_test ();
+  void solve ();
+  void create_coarse_grid ();
+  void postprocess ();
+
+  void local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+                       Assembly::Scratch::Data<dim>  &scratch,
+                       Assembly::Copy::Data          &data);
+  void copy_local_to_global (const Assembly::Copy::Data &data);
+
+  std::vector<types::global_dof_index>
+  get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const;
+
+  parallel::distributed::Triangulation<dim>   triangulation;
+
+  DoFHandler<dim>      dof_handler;
+  FE_Q<dim>            fe;
+  QGauss<dim>          quadrature;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::SparseMatrix reference_matrix;
+  TrilinosWrappers::SparseMatrix test_matrix;
+
+  TrilinosWrappers::MPI::Vector       reference_rhs;
+  TrilinosWrappers::MPI::Vector       test_rhs;
+
+  std::vector<std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > graph;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                            const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                           const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+  :
+  triangulation (MPI_COMM_WORLD),
+  dof_handler (triangulation),
+  fe (1),
+  quadrature(fe.degree+1)
+{
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const
+{
+  std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+  cell->get_dof_indices(local_dof_indices);
+
+  constraints.resolve_indices(local_dof_indices);
+  return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  reference_matrix.clear();
+  test_matrix.clear();
+  dof_handler.distribute_dofs (fe);
+
+  constraints.clear ();
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+  // add boundary conditions as inhomogeneous constraints here, do it after
+  // having added the hanging node constraints in order to be consistent and
+  // skip dofs that are already constrained (i.e., are hanging nodes on the
+  // boundary in 3D). In contrast to step-27, we choose a sine function.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  typedef FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> CellFilter;
+  CellFilter begin(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active());
+  CellFilter end(IteratorFilters::LocallyOwnedCell(),dof_handler.end());
+  graph = GraphColoring::make_graph_coloring(begin,end,
+        static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+        (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);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+                                     Assembly::Scratch::Data<dim>  &scratch,
+                                     Assembly::Copy::Data          &data)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+  data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+  data.local_matrix = 0;
+
+  data.local_rhs.reinit (dofs_per_cell);
+  data.local_rhs = 0;
+
+  scratch.fe_values.reinit (cell);
+
+  const FEValues<dim> &fe_values = scratch.fe_values;
+
+  const RightHandSide<dim> rhs_function;
+
+  for (unsigned int q_point=0;
+       q_point<fe_values.n_quadrature_points;
+       ++q_point)
+    {
+      const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            data.local_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                       fe_values.shape_grad(j,q_point) *
+                                       fe_values.JxW(q_point));
+
+            data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                  rhs_value *
+                                  fe_values.JxW(q_point));
+        }
+    }
+
+  data.local_dof_indices.resize (dofs_per_cell);
+  cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+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);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  Assembly::Copy::Data copy_data;
+  Assembly::Scratch::Data<dim> assembly_data(fe, quadrature);
+
+  for (unsigned int color=0; color<graph.size(); ++color)
+    for (typename std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >::const_iterator p = graph[color].begin();
+         p != graph[color].end(); ++p)
+      {
+        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;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe, quadrature),
+         Assembly::Copy::Data (),
+         1);
+  test_matrix.compress(VectorOperation::add);
+  test_rhs.compress(VectorOperation::add);
+  test_matrix.compress(VectorOperation::add);
+
+  test_matrix.add(-1, reference_matrix);
+  test_matrix.compress(VectorOperation::add);
+
+  // there should not even be roundoff difference between matrices
+  deallog.threshold_double(1.e-30);
+  deallog << "error in matrix: " << test_matrix.frobenius_norm() << std::endl;
+  test_rhs.add(-1., reference_rhs);
+  deallog << "error in vector: " << test_rhs.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+    estimated_error_per_cell(i) = i;
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube(triangulation,
+                                     Point<dim>(),
+                                     0.5, 1., (dim==3) ? 96 : 12, false);
+          triangulation.refine_global(3);
+        }
+
+      setup_system ();
+
+      assemble_reference ();
+      assemble_test ();
+
+      if (cycle < 2)
+        postprocess ();
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int);
+
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output b/tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output
new file mode 100644 (file)
index 0000000..6e5eabc
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
diff --git a/tests/trilinos/assemble_matrix_parallel_04.cc b/tests/trilinos/assemble_matrix_parallel_04.cc
new file mode 100644 (file)
index 0000000..90ffbde
--- /dev/null
@@ -0,0 +1,463 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Same as assemble_matrix_parallel_02, but using block matrices
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/graph_coloring.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim>
+    struct Data
+    {
+      Data (const FiniteElement<dim> &fe,
+            const Quadrature<dim>    &quadrature)
+        :
+        fe_values(fe,
+                  quadrature,
+                  update_values    |  update_gradients |
+                  update_quadrature_points  |  update_JxW_values)
+      {}
+
+      Data (const Data &data)
+        :
+        fe_values(data.fe_values.get_mapping(),
+                  data.fe_values.get_fe(),
+                  data.fe_values.get_quadrature(),
+                  data.fe_values.get_update_flags())
+      {}
+
+      FEValues<dim> fe_values;
+    };
+  }
+
+  namespace Copy
+  {
+    struct Data
+    {
+      std::vector<types::global_dof_index> local_dof_indices;
+      FullMatrix<double> local_matrix;
+      Vector<double> local_rhs;
+    };
+  }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem ();
+  ~LaplaceProblem ();
+
+  void run ();
+
+private:
+  void setup_system ();
+  void test_equality ();
+  void assemble_reference ();
+  void assemble_test ();
+  void solve ();
+  void create_coarse_grid ();
+  void postprocess ();
+
+  void local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+                       Assembly::Scratch::Data<dim>  &scratch,
+                       Assembly::Copy::Data          &data);
+  void copy_local_to_global (const Assembly::Copy::Data &data);
+
+  std::vector<types::global_dof_index>
+  get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const;
+
+  parallel::distributed::Triangulation<dim>   triangulation;
+
+  DoFHandler<dim>      dof_handler;
+  FESystem<dim>        fe;
+  QGauss<dim>          quadrature;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::BlockSparseMatrix reference_matrix;
+  TrilinosWrappers::BlockSparseMatrix test_matrix;
+
+  TrilinosWrappers::MPI::BlockVector       reference_rhs;
+  TrilinosWrappers::MPI::BlockVector       test_rhs;
+
+  std::vector<std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > graph;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim> (2) {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                            const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                           const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+  :
+  triangulation (MPI_COMM_WORLD),
+  dof_handler (triangulation),
+  fe (FE_Q<dim>(1),1, FE_Q<dim>(2),1),
+  quadrature(3)
+{
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const
+{
+  std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+  cell->get_dof_indices(local_dof_indices);
+
+  constraints.resolve_indices(local_dof_indices);
+  return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  reference_matrix.clear();
+  test_matrix.clear();
+  dof_handler.distribute_dofs (fe);
+  std::vector<unsigned int> blocks(2,0);
+  blocks[1] = 1;
+  DoFRenumbering::component_wise(dof_handler, blocks);
+
+  constraints.clear ();
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+  // add boundary conditions as inhomogeneous constraints here, do it after
+  // having added the hanging node constraints in order to be consistent and
+  // skip dofs that are already constrained (i.e., are hanging nodes on the
+  // boundary in 3D). In contrast to step-27, we choose a sine function.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  typedef FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> CellFilter;
+  CellFilter begin(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active());
+  CellFilter end(IteratorFilters::LocallyOwnedCell(),dof_handler.end());
+  graph = GraphColoring::make_graph_coloring(begin,end,
+        static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+        (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &)> >
+        (std_cxx1x::bind(&LaplaceProblem<dim>::get_conflict_indices, this,std_cxx1x::_1)));
+
+  TrilinosWrappers::BlockSparsityPattern csp(2,2);
+  std::vector<IndexSet> locally_owned(2), relevant_set(2);
+  IndexSet locally_owned_total = dof_handler.locally_owned_dofs(), relevant_total;
+  DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_total);
+
+  std::vector<types::global_dof_index> dofs_per_block (2);
+  DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, blocks);
+  locally_owned[0] = locally_owned_total.get_view(0, dofs_per_block[0]);
+  locally_owned[1] = locally_owned_total.get_view(dofs_per_block[0],
+                                                  dof_handler.n_dofs());
+  relevant_set[0] = relevant_total.get_view(0, dofs_per_block[0]);
+  relevant_set[1] = relevant_total.get_view(dofs_per_block[0],
+                                            dof_handler.n_dofs());
+
+  csp.reinit(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);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+                                     Assembly::Scratch::Data<dim>  &scratch,
+                                     Assembly::Copy::Data          &data)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+  data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+  data.local_matrix = 0;
+
+  data.local_rhs.reinit (dofs_per_cell);
+  data.local_rhs = 0;
+
+  scratch.fe_values.reinit (cell);
+
+  const FEValues<dim> &fe_values = scratch.fe_values;
+
+  const RightHandSide<dim> rhs_function;
+
+  // this does not make a lot of sense physically but it serves the purpose of
+  // the test well
+  for (unsigned int q_point=0;
+       q_point<fe_values.n_quadrature_points;
+       ++q_point)
+    {
+      const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            data.local_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                       fe_values.shape_grad(j,q_point) *
+                                       fe_values.JxW(q_point));
+
+            data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                  rhs_value *
+                                  fe_values.JxW(q_point));
+        }
+    }
+
+  data.local_dof_indices.resize (dofs_per_cell);
+  cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+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);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  Assembly::Copy::Data copy_data;
+  Assembly::Scratch::Data<dim> assembly_data(fe, quadrature);
+
+  for (unsigned int color=0; color<graph.size(); ++color)
+    for (typename std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >::const_iterator p = graph[color].begin();
+         p != graph[color].end(); ++p)
+      {
+        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;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe, quadrature),
+         Assembly::Copy::Data (),
+         1);
+  test_matrix.compress(VectorOperation::add);
+  test_rhs.compress(VectorOperation::add);
+
+  test_matrix.add(-1, reference_matrix);
+
+  // there should not even be roundoff difference between matrices
+  deallog.threshold_double(1.e-30);
+  double frobenius_norm = 0;
+  for (unsigned int i=0; i<2; ++i)
+    for (unsigned int j=0; j<2; ++j)
+      frobenius_norm += numbers::NumberTraits<double>::abs_square(test_matrix.block(i,j).frobenius_norm());
+  deallog << "error in matrix: " << std::sqrt(frobenius_norm) << std::endl;
+  test_rhs.add(-1., reference_rhs);
+  deallog << "error in vector: " << test_rhs.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+    estimated_error_per_cell(i) = i;
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_shell(triangulation,
+                                     Point<dim>(),
+                                     0.5, 1., (dim==3) ? 96 : 12, false);
+#ifdef DEBUG
+          triangulation.refine_global(3);
+#else
+          triangulation.refine_global(5);
+#endif
+        }
+
+      setup_system ();
+
+      assemble_reference ();
+      assemble_test ();
+
+      if (cycle < 2)
+        postprocess ();
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int);
+
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output b/tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output
new file mode 100644 (file)
index 0000000..6e5eabc
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0

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.