<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
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> ¶llel_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> ¶llel_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> ¶llel_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> ¶llel_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;
};
<< " 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
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;
/**
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;
};
{
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>
private:
/**
- * Make iterator class a
- * friend.
+ * Make iterator class a friend.
*/
template <bool> friend class Iterator;
};
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,
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;
};
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);
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;
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;
private:
/**
- * Store an object of the
- * accessor class.
+ * Store an object of the accessor class.
*/
Accessor<Constness> accessor;
* 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
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;
*/
//@{
/**
- * 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,
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
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);
*/
//@{
/**
- * 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 ¶llel_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 ¶llel_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
- * – 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 – 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 – 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 –
+ * 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 ¶llel_partitioning,
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,
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 ¶llel_partitioning,
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,
*/
//@{
/**
- * 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 ¶llel_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 ¶llel_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,
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 – 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 –
+ * 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,
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 ¶llel_partitioning,
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,
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 ¶llel_partitioning,
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,
//@{
/**
- * 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;
//@{
/**
- * 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.
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,
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,
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,
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,
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,
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,
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,
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);
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 —
- * 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 — 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 ();
//@{
/**
- * 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;
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,
/**
- * 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,
//@{
/**
- * 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;
//@{
/**
- * 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;
//@}
//@{
/**
- * STL-like iterator with the
- * first entry.
+ * STL-like iterator with the first entry.
*/
const_iterator begin () const;
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 ();
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);
//@{
/**
- * 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;
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();
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>;
};
- 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 ()
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,
// 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;
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,
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
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;
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;
};
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,
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;
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;
private:
/**
- * Store an object of the
- * accessor class.
+ * Store an object of the accessor class.
*/
Accessor accessor;
*/
//@{
/**
- * 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,
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,
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 ();
//@}
//@{
/**
- * 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 ¶llel_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 ¶llel_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
- * – 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
+ * – 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,
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 – 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 – 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 ¶llel_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
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 – 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
+ * – 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
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 – 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 – 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,
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
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
//@{
/**
- * 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 ¶llel_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 ¶llel_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 – 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 – 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,
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 – 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 – 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,
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 ¶llel_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
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 – 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
+ * – 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
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,
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
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
//@{
/**
- * 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;
*/
//@{
/**
- * 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,
//@{
/**
- * 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;
//@}
//@{
/**
- * STL-like iterator with the
- * first entry.
+ * STL-like iterator with the first entry.
*/
const_iterator begin () const;
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;
//@{
/**
- * 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;
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;
};
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));
}
return graph->RangeMap().Comm();
}
-
-
- inline
- SparsityPattern::SparsityPattern (const IndexSet ¶llel_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 ¶llel_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 ¶llel_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 ¶llel_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 ¶llel_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
}
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
# 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>
{
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;
}
+ 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 ()
{
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();
# 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
+ SparsityPattern::SparsityPattern (const IndexSet ¶llel_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 ¶llel_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 ()
{}
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;
+ void
+ SparsityPattern::reinit (const IndexSet ¶llel_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 ¶llel_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 ¶llel_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,
*column_space_map, 0));
graph->FillComplete();
+ nonlocal_graph.reset();
+
compressed = true;
}
{
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));
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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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();
+ }
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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();
+ }
+}
+
--- /dev/null
+
+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
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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();
+ }
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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();
+ }
+}
+
--- /dev/null
+
+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
projects / dealii.git / commitdiff