From 0c69608b21f8adf03c6bee7e8cc7751c09798f2d Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 23 Dec 2013 22:14:25 +0000 Subject: [PATCH] Implement new reinit method to Trilinos sparsity pattern that makes writing into matrices thread-safe. git-svn-id: https://svn.dealii.org/trunk@32106 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 11 + .../deal.II/lac/block_sparsity_pattern.h | 121 +- .../deal.II/lac/trilinos_sparse_matrix.h | 2096 ++++++----------- .../deal.II/lac/trilinos_sparsity_pattern.h | 1252 ++++------ deal.II/source/lac/block_sparsity_pattern.cc | 39 + deal.II/source/lac/trilinos_sparse_matrix.cc | 61 +- .../source/lac/trilinos_sparsity_pattern.cc | 277 ++- tests/trilinos/assemble_matrix_parallel_01.cc | 451 ++++ .../assemble_matrix_parallel_01.output | 13 + tests/trilinos/assemble_matrix_parallel_02.cc | 437 ++++ ...ssemble_matrix_parallel_02.mpirun=1.output | 7 + ...ssemble_matrix_parallel_02.mpirun=4.output | 7 + tests/trilinos/assemble_matrix_parallel_03.cc | 440 ++++ ...ssemble_matrix_parallel_03.mpirun=4.output | 7 + tests/trilinos/assemble_matrix_parallel_04.cc | 463 ++++ ...ssemble_matrix_parallel_04.mpirun=4.output | 7 + 16 files changed, 3386 insertions(+), 2303 deletions(-) create mode 100644 tests/trilinos/assemble_matrix_parallel_01.cc create mode 100644 tests/trilinos/assemble_matrix_parallel_01.output create mode 100644 tests/trilinos/assemble_matrix_parallel_02.cc create mode 100644 tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output create mode 100644 tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output create mode 100644 tests/trilinos/assemble_matrix_parallel_03.cc create mode 100644 tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output create mode 100644 tests/trilinos/assemble_matrix_parallel_04.cc create mode 100644 tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index be811da606..e2c2969b4b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -84,6 +84,17 @@ inconvenience this causes.

Specific improvements

    + +
  1. 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). +
    + (Martin Kronbichler, 2013/12/23) +
    2. +
  2. 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 diff --git a/deal.II/include/deal.II/lac/block_sparsity_pattern.h b/deal.II/include/deal.II/lac/block_sparsity_pattern.h index 392eaa4fca..e22741f97c 100644 --- a/deal.II/include/deal.II/lac/block_sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/block_sparsity_pattern.h @@ -1016,113 +1016,100 @@ namespace TrilinosWrappers public: /** - * Initialize the matrix empty, - * that is with no memory - * allocated. This is useful if - * you want such objects as - * member variables in other - * classes. You can make the - * structure usable by calling - * the reinit() function. + * Initialize the matrix empty, that is with no memory allocated. This is + * useful if you want such objects as member variables in other + * classes. You can make the structure usable by calling the reinit() + * function. */ BlockSparsityPattern (); /** - * Initialize the matrix with the - * given number of block rows and - * columns. The blocks themselves - * are still empty, and you have - * to call collect_sizes() after - * you assign them sizes. + * Initialize the matrix with the given number of block rows and + * columns. The blocks themselves are still empty, and you have to call + * collect_sizes() after you assign them sizes. */ BlockSparsityPattern (const size_type n_rows, const size_type n_columns); /** - * Initialize the pattern with - * two BlockIndices for the block - * structures of matrix rows and - * columns. This function is - * equivalent to calling the - * previous constructor with the - * length of the two index vector - * and then entering the index - * values. + * Initialize the pattern with two BlockIndices for the block structures + * of matrix rows and columns. This function is equivalent to calling the + * previous constructor with the length of the two index vector and then + * entering the index values. */ BlockSparsityPattern (const std::vector &row_block_sizes, const std::vector &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 ¶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 ¶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,j) will - * have the dimensions - * row_block_sizes[i] - * times col_block_sizes[j]. + * The matrix will have as many block rows and columns as there are + * entries in the two arguments. The block at position (i,j) will + * have the dimensions row_block_sizes[i] times + * col_block_sizes[j]. */ void reinit (const std::vector &row_block_sizes, const std::vector &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 ¶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 ¶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 &row_parallel_partitioning, + const std::vector &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 &row_parallel_partitioning, + const std::vector &column_parallel_partitioning, + const std::vector &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::reinit; }; diff --git a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h index ed72eb0060..8babcea003 100644 --- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h @@ -84,17 +84,11 @@ namespace TrilinosWrappers << " are stored locally and can be accessed."); /** - * Handling of indices for both - * constant and non constant - * Accessor objects - * - * For a regular - * dealii::SparseMatrix, we would - * use an accessor for the sparsity - * pattern. For Trilinos matrices, - * this does not seem so simple, - * therefore, we write a little - * base class here. + * Handling of indices for both constant and non constant Accessor objects + * + * For a regular dealii::SparseMatrix, we would use an accessor for the + * sparsity pattern. For Trilinos matrices, this does not seem so simple, + * therefore, we write a little base class here. * * @author Guido Kanschat * @date 2012 @@ -115,34 +109,26 @@ namespace TrilinosWrappers const size_type index); /** - * Row number of the element - * represented by this object. + * Row number of the element represented by this object. */ size_type row() const; /** - * Index in row of the element - * represented by this object. + * Index in row of the element represented by this object. */ size_type index() const; /** - * Column number of the element - * represented by this object. + * Column number of the element represented by this object. */ size_type column() const; protected: /** - * Pointer to the matrix - * object. This object should - * be handled as a const - * pointer or non-const by the - * appropriate derived - * classes. In order to be able - * to implement both, it is not - * const here, so handle with - * care! + * Pointer to the matrix object. This object should be handled as a + * const pointer or non-const by the appropriate derived classes. In + * order to be able to implement both, it is not const here, so handle + * with care! */ mutable SparseMatrix *matrix; /** @@ -156,45 +142,28 @@ namespace TrilinosWrappers size_type a_index; /** - * Discard the old row caches - * (they may still be used by - * other accessors) and - * generate new ones for the - * row pointed to presently by + * Discard the old row caches (they may still be used by other + * accessors) and generate new ones for the row pointed to presently by * this accessor. */ void visit_present_row (); /** - * Cache where we store the - * column indices of the - * present row. This is - * necessary, since Trilinos - * makes access to the elements - * of its matrices rather hard, - * and it is much more - * efficient to copy all column - * entries of a row once when - * we enter it than repeatedly - * asking Trilinos for - * individual ones. This also - * makes some sense since it is - * likely that we will access - * them sequentially anyway. + * Cache where we store the column indices of the present row. This is + * necessary, since Trilinos makes access to the elements of its + * matrices rather hard, and it is much more efficient to copy all + * column entries of a row once when we enter it than repeatedly asking + * Trilinos for individual ones. This also makes some sense since it is + * likely that we will access them sequentially anyway. * - * In order to make copying of - * iterators/accessor of - * acceptable performance, we - * keep a shared pointer to - * these entries so that more - * than one accessor can access - * this data if necessary. + * In order to make copying of iterators/accessor of acceptable + * performance, we keep a shared pointer to these entries so that more + * than one accessor can access this data if necessary. */ std_cxx1x::shared_ptr > colnum_cache; /** - * Cache for the values - * of this row. + * Cache for the values of this row. */ std_cxx1x::shared_ptr > value_cache; }; @@ -231,25 +200,21 @@ namespace TrilinosWrappers { public: /** - * Typedef for the type (including - * constness) of the matrix to be - * used here. + * Typedef for the type (including constness) of the matrix to be used + * here. */ typedef const SparseMatrix MatrixType; /** - * Constructor. Since we use - * accessors only for read - * access, a const matrix - * pointer is sufficient. + * Constructor. Since we use accessors only for read access, a const + * matrix pointer is sufficient. */ Accessor (MatrixType *matrix, const size_type row, const size_type index); /** - * Copy constructor to get from a - * const or non-const accessor to a const + * Copy constructor to get from a const or non-const accessor to a const * accessor. */ template @@ -262,8 +227,7 @@ namespace TrilinosWrappers private: /** - * Make iterator class a - * friend. + * Make iterator class a friend. */ template friend class Iterator; }; @@ -283,66 +247,53 @@ namespace TrilinosWrappers Reference (const Accessor &accessor); /** - * Conversion operator to the - * data type of the matrix. + * Conversion operator to the data type of the matrix. */ operator TrilinosScalar () const; /** - * Set the element of the matrix - * we presently point to to @p n. + * Set the element of the matrix we presently point to to @p n. */ const Reference &operator = (const TrilinosScalar n) const; /** - * Add @p n to the element of the - * matrix we presently point to. + * Add @p n to the element of the matrix we presently point to. */ const Reference &operator += (const TrilinosScalar n) const; /** - * Subtract @p n from the element - * of the matrix we presently - * point to. + * Subtract @p n from the element of the matrix we presently point to. */ const Reference &operator -= (const TrilinosScalar n) const; /** - * Multiply the element of the - * matrix we presently point to - * by @p n. + * Multiply the element of the matrix we presently point to by @p n. */ const Reference &operator *= (const TrilinosScalar n) const; /** - * Divide the element of the - * matrix we presently point to - * by @p n. + * Divide the element of the matrix we presently point to by @p n. */ const Reference &operator /= (const TrilinosScalar n) const; private: /** - * Pointer to the accessor that - * denotes which element we - * presently point to. + * Pointer to the accessor that denotes which element we presently + * point to. */ Accessor &accessor; }; public: /** - * Typedef for the type (including - * constness) of the matrix to be - * used here. + * Typedef for the type (including constness) of the matrix to be used + * here. */ typedef SparseMatrix MatrixType; /** - * Constructor. Since we use - * accessors only for read - * access, a const matrix - * pointer is sufficient. + * Constructor. Since we use accessors only for read access, a const + * matrix pointer is sufficient. */ Accessor (MatrixType *matrix, const size_type row, @@ -355,13 +306,11 @@ namespace TrilinosWrappers private: /** - * Make iterator class a - * friend. + * Make iterator class a friend. */ template friend class Iterator; /** - * Make Reference object a - * friend. + * Make Reference object a friend. */ friend class Reference; }; @@ -390,25 +339,21 @@ namespace TrilinosWrappers typedef dealii::types::global_dof_index size_type; /** - * Typedef for the matrix type - * (including constness) we are to - * operate on. + * Typedef for the matrix type (including constness) we are to operate + * on. */ typedef typename Accessor::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 Iterator(const Iterator &other); @@ -434,9 +379,8 @@ namespace TrilinosWrappers const Accessor *operator-> () const; /** - * Comparison. True, if both - * iterators point to the same - * matrix position. + * Comparison. True, if both iterators point to the same matrix + * position. */ bool operator == (const Iterator &) const; @@ -446,11 +390,8 @@ namespace TrilinosWrappers bool operator != (const Iterator &) const; /** - * Comparison operator. Result - * is true if either the first - * row number is smaller or if - * the row numbers are equal - * and the first index is + * Comparison operator. Result is true if either the first row number is + * smaller or if the row numbers are equal and the first index is * smaller. */ bool operator < (const Iterator &) const; @@ -471,8 +412,7 @@ namespace TrilinosWrappers private: /** - * Store an object of the - * accessor class. + * Store an object of the accessor class. */ Accessor accessor; @@ -508,6 +448,29 @@ namespace TrilinosWrappers * unused elements. Trilinos allows to continue with assembling the * matrix after calls to these functions, though. * + *

    Thread safety of Trilinos matrices

    + * + * When writing into Trilinos matrices from several threads in shared + * memory, several things must be kept in mind as there is no built-in locks + * in this class to prevent data races. Therefore, simultaneous access to + * the same matrix row at the same time leads to data races in general and + * must be explicitly avoided by the user. However, it is possible to access + * different rows of the matrix from several threads simultaneously + * under the following two conditions: + *
      + *
    • The matrix uses only one MPI process. + *
    • The matrix has been initialized from a + * TrilinosWrappers::SparsityPattern object that in turn has been + * initialized with the reinit function specifying three index sets, one + * for the rows, one for the columns and for the larger set of @p + * writeable_rows. Note that all other reinit methods and constructors of + * TrilinosWrappers::SparsityPattern will result in a matrix that needs to + * allocate off-processor entries on demand, which breaks + * thread-safety. Of course, using the respective reinit method for the + * block Trilinos sparsity pattern and block matrix also results in + * thread-safety. + *
    + * * @ingroup TrilinosWrappers * @ingroup Matrix1 * @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009 @@ -521,43 +484,33 @@ namespace TrilinosWrappers typedef dealii::types::global_dof_index size_type; /** - * A structure that describes - * some of the traits of this - * class in terms of its run-time - * behavior. Some other classes - * (such as the block matrix - * classes) that take one or - * other of the matrix classes as - * its template parameters can - * tune their behavior based on - * the variables in this class. + * A structure that describes some of the traits of this class in terms of + * its run-time behavior. Some other classes (such as the block matrix + * classes) that take one or other of the matrix classes as its template + * parameters can tune their behavior based on the variables in this + * class. */ struct Traits { /** - * It is safe to elide additions - * of zeros to individual - * elements of this matrix. + * It is safe to elide additions of zeros to individual elements of this + * matrix. */ static const bool zero_addition_can_be_elided = true; }; /** - * Declare a typedef for the - * iterator class. + * Declare a typedef for the iterator class. */ typedef SparseMatrixIterators::Iterator iterator; /** - * Declare a typedef for the - * const iterator class. + * Declare a typedef for the const iterator class. */ typedef SparseMatrixIterators::Iterator const_iterator; /** - * Declare a typedef in analogy - * to all the other container - * classes. + * Declare a typedef in analogy to all the other container classes. */ typedef TrilinosScalar value_type; @@ -566,18 +519,15 @@ namespace TrilinosWrappers */ //@{ /** - * Default constructor. Generates - * an empty (zero-size) matrix. + * Default constructor. Generates an empty (zero-size) matrix. */ SparseMatrix (); /** - * Generate a matrix that is completely - * stored locally, having #m rows and + * Generate a matrix that is completely stored locally, having #m rows and * #n columns. * - * The number of columns entries per - * row is specified as the maximum + * The number of columns entries per row is specified as the maximum * number of entries argument. */ SparseMatrix (const size_type m, @@ -585,118 +535,85 @@ namespace TrilinosWrappers const unsigned int n_max_entries_per_row); /** - * Generate a matrix that is completely - * stored locally, having #m rows and + * Generate a matrix that is completely stored locally, having #m rows and * #n columns. * - * The vector - * n_entries_per_row - * specifies the number of entries in - * each row. + * The vector n_entries_per_row specifies the number of entries + * in each row. */ SparseMatrix (const size_type m, const size_type n, const std::vector &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 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 sparse_matrix to the + * This function copies the content in sparse_matrix 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 - * copy_values 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 copy_values decides whether only the + * sparsity structure of the input matrix should be used or the matrix * entries should be copied, too. * - * This is a collective operation - * that needs to be called on all - * processors in order to avoid a - * deadlock. + * This is a collective operation that needs to be called on all + * processors in order to avoid a deadlock. * * @note If a different sparsity pattern is given in the last argument * (i.e., one that differs from the one used in the sparse matrix given @@ -712,11 +629,9 @@ namespace TrilinosWrappers const ::dealii::SparsityPattern *use_this_sparsity=0); /** - * This reinit function takes as - * input a Trilinos Epetra_CrsMatrix - * and copies its sparsity - * pattern. If so requested, even the - * content (values) will be copied. + * This reinit function takes as input a Trilinos Epetra_CrsMatrix and + * copies its sparsity pattern. If so requested, even the content (values) + * will be copied. */ void reinit (const Epetra_CrsMatrix &input_matrix, const bool copy_values = true); @@ -726,150 +641,87 @@ namespace TrilinosWrappers */ //@{ /** - * Constructor using an Epetra_Map to - * describe the %parallel - * partitioning. The parameter @p - * n_max_entries_per_row sets the - * number of nonzero entries in each - * row that will be allocated. Note - * that this number does not need to - * be exact, and it is even allowed - * that the actual matrix structure - * has more nonzero entries than - * specified in the - * constructor. However it is still - * advantageous to provide good - * estimates here since this will - * considerably increase the - * performance of the matrix - * setup. However, there is no effect - * in the performance of - * matrix-vector products, since - * Trilinos reorganizes the matrix - * memory prior to use (in the - * compress() step). + * Constructor using an Epetra_Map to describe the %parallel + * partitioning. The parameter @p n_max_entries_per_row sets the number of + * nonzero entries in each row that will be allocated. Note that this + * number does not need to be exact, and it is even allowed that the + * actual matrix structure has more nonzero entries than specified in the + * constructor. However it is still advantageous to provide good estimates + * here since this will considerably increase the performance of the + * matrix setup. However, there is no effect in the performance of + * matrix-vector products, since Trilinos reorganizes the matrix memory + * prior to use (in the compress() step). */ SparseMatrix (const Epetra_Map ¶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 &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 - * n_entries_per_row - * 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 n_entries_per_row 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 &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 void reinit (const Epetra_Map ¶llel_partitioning, @@ -877,25 +729,16 @@ namespace TrilinosWrappers const bool exchange_data = false); /** - * This function is similar to the - * other initialization function - * above, but now also reassigns the - * matrix rows and columns according - * to two user-supplied Epetra maps. - * To be used for rectangular - * matrices. The optional argument @p - * exchange_data can be used for - * reinitialization with a sparsity - * pattern that is not fully - * constructed. This feature is only - * implemented for input sparsity - * patterns of type + * This function is similar to the other initialization function above, + * but now also reassigns the matrix rows and columns according to two + * user-supplied Epetra maps. To be used for rectangular matrices. The + * optional argument @p exchange_data can be used for reinitialization + * with a sparsity pattern that is not fully constructed. This feature is + * only implemented for input sparsity patterns of type * CompressedSimpleSparsityPattern. * - * This is a collective operation - * that needs to be called on all - * processors in order to avoid a - * dead lock. + * This is a collective operation that needs to be called on all + * processors in order to avoid a dead lock. */ template void reinit (const Epetra_Map &row_parallel_partitioning, @@ -904,32 +747,20 @@ namespace TrilinosWrappers const bool exchange_data = false); /** - * This function initializes the - * Trilinos matrix using the deal.II - * sparse matrix and the entries - * stored therein. It uses a - * threshold to copy only elements - * with modulus larger than the - * threshold (so zeros in the deal.II - * matrix can be filtered away). In - * contrast to the other reinit - * function with deal.II sparse - * matrix argument, this function - * takes a %parallel partitioning - * specified by the user instead of - * internally generating it. - * - * The optional parameter - * copy_values 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 copy_values 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 void reinit (const Epetra_Map ¶llel_partitioning, @@ -939,25 +770,17 @@ namespace TrilinosWrappers const ::dealii::SparsityPattern *use_this_sparsity=0); /** - * This function is similar to the - * other initialization function with - * deal.II sparse matrix input above, - * but now takes Epetra maps for both - * the rows and the columns of the - * matrix. Chosen for rectangular + * This function is similar to the other initialization function with + * deal.II sparse matrix input above, but now takes Epetra maps for both + * the rows and the columns of the matrix. Chosen for rectangular * matrices. * - * The optional parameter - * copy_values decides - * whether only the sparsity - * structure of the input matrix - * should be used or the matrix + * The optional parameter copy_values 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 void reinit (const Epetra_Map &row_parallel_partitioning, @@ -972,77 +795,45 @@ namespace TrilinosWrappers */ //@{ /** - * Constructor using an IndexSet and - * an MPI communicator to describe - * the %parallel partitioning. The - * parameter @p n_max_entries_per_row - * sets the number of nonzero entries - * in each row that will be - * allocated. Note that this number - * does not need to be exact, and it - * is even allowed that the actual - * matrix structure has more nonzero - * entries than specified in the - * constructor. However it is still - * advantageous to provide good - * estimates here since this will - * considerably increase the - * performance of the matrix - * setup. However, there is no effect - * in the performance of - * matrix-vector products, since - * Trilinos reorganizes the matrix - * memory prior to use (in the - * compress() step). + * Constructor using an IndexSet and an MPI communicator to describe the + * %parallel partitioning. The parameter @p n_max_entries_per_row sets the + * number of nonzero entries in each row that will be allocated. Note that + * this number does not need to be exact, and it is even allowed that the + * actual matrix structure has more nonzero entries than specified in the + * constructor. However it is still advantageous to provide good estimates + * here since this will considerably increase the performance of the + * matrix setup. However, there is no effect in the performance of + * matrix-vector products, since Trilinos reorganizes the matrix memory + * prior to use (in the compress() step). */ SparseMatrix (const IndexSet ¶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 &n_entries_per_row); /** - * This constructor is similar to the - * one above, but it now takes two - * different IndexSet partitions for - * row and columns. This interface is - * meant to be used for generating - * rectangular matrices, where the - * first index set describes the - * %parallel partitioning of the - * degrees of freedom associated with - * the matrix rows and the second one - * the partitioning of the matrix - * columns. The second index set - * specifies the partitioning of the - * vectors this matrix is to be - * multiplied with, not the - * distribution of the elements that - * actually appear in the matrix. - * - * The parameter @p - * n_max_entries_per_row defines how - * much memory will be allocated for - * each row. This number does not - * need to be accurate, as the - * structure is reorganized in the - * compress() call. + * This constructor is similar to the one above, but it now takes two + * different IndexSet partitions for row and columns. This interface is + * meant to be used for generating rectangular matrices, where the first + * index set describes the %parallel partitioning of the degrees of + * freedom associated with the matrix rows and the second one the + * partitioning of the matrix columns. The second index set specifies the + * partitioning of the vectors this matrix is to be multiplied with, not + * the distribution of the elements that actually appear in the matrix. + * + * The parameter @p n_max_entries_per_row defines how much memory will be + * allocated for each row. This number does not need to be accurate, as + * the structure is reorganized in the compress() call. */ SparseMatrix (const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, @@ -1050,33 +841,18 @@ namespace TrilinosWrappers const size_type n_max_entries_per_row = 0); /** - * This constructor is similar to the - * one above, but it now takes two - * different Epetra maps for rows and - * columns. This interface is meant - * to be used for generating - * rectangular matrices, where one - * map specifies the %parallel - * distribution of degrees of freedom - * associated with matrix rows and - * the second one specifies the - * %parallel distribution the dofs - * associated with columns in the - * matrix. The second map also - * provides information for the - * internal arrangement in matrix - * vector products (i.e., the - * distribution of vector this matrix - * is to be multiplied with), but is - * not used for the distribution of - * the columns – rather, all - * column elements of a row are - * stored on the same processor in - * any case. The vector - * n_entries_per_row - * 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 n_entries_per_row specifies the number + * of entries in each row of the newly generated matrix. */ SparseMatrix (const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, @@ -1084,43 +860,24 @@ namespace TrilinosWrappers const std::vector &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 void reinit (const IndexSet ¶llel_partitioning, @@ -1129,25 +886,16 @@ namespace TrilinosWrappers const bool exchange_data = false); /** - * This function is similar to the - * other initialization function - * above, but now also reassigns the - * matrix rows and columns according - * to two user-supplied index sets. - * To be used for rectangular - * matrices. The optional argument @p - * exchange_data can be used for - * reinitialization with a sparsity - * pattern that is not fully - * constructed. This feature is only - * implemented for input sparsity - * patterns of type + * This function is similar to the other initialization function above, + * but now also reassigns the matrix rows and columns according to two + * user-supplied index sets. To be used for rectangular matrices. The + * optional argument @p exchange_data can be used for reinitialization + * with a sparsity pattern that is not fully constructed. This feature is + * only implemented for input sparsity patterns of type * CompressedSimpleSparsityPattern. * - * This is a collective operation - * that needs to be called on all - * processors in order to avoid a - * dead lock. + * This is a collective operation that needs to be called on all + * processors in order to avoid a dead lock. */ template void reinit (const IndexSet &row_parallel_partitioning, @@ -1157,32 +905,20 @@ namespace TrilinosWrappers const bool exchange_data = false); /** - * This function initializes the - * Trilinos matrix using the deal.II - * sparse matrix and the entries - * stored therein. It uses a - * threshold to copy only elements - * with modulus larger than the - * threshold (so zeros in the deal.II - * matrix can be filtered away). In - * contrast to the other reinit - * function with deal.II sparse - * matrix argument, this function - * takes a %parallel partitioning - * specified by the user instead of - * internally generating it. - * - * The optional parameter - * copy_values 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 copy_values 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 void reinit (const IndexSet ¶llel_partitioning, @@ -1193,25 +929,17 @@ namespace TrilinosWrappers const ::dealii::SparsityPattern *use_this_sparsity=0); /** - * This function is similar to the - * other initialization function with - * deal.II sparse matrix input above, - * but now takes index sets for both - * the rows and the columns of the - * matrix. Chosen for rectangular + * This function is similar to the other initialization function with + * deal.II sparse matrix input above, but now takes index sets for both + * the rows and the columns of the matrix. Chosen for rectangular * matrices. * - * The optional parameter - * copy_values decides - * whether only the sparsity - * structure of the input matrix - * should be used or the matrix + * The optional parameter copy_values 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 void reinit (const IndexSet &row_parallel_partitioning, @@ -1228,89 +956,64 @@ namespace TrilinosWrappers //@{ /** - * Return the number of rows in - * this matrix. + * Return the number of rows in this matrix. */ size_type m () const; /** - * Return the number of columns - * in this matrix. + * Return the number of columns in this matrix. */ size_type n () const; /** - * Return the local dimension - * of the matrix, i.e. the - * number of rows stored on the - * present MPI process. For - * sequential matrices, this - * number is the same as m(), - * but for %parallel matrices it - * may be smaller. + * Return the local dimension of the matrix, i.e. the number of rows + * stored on the present MPI process. For sequential matrices, this number + * is the same as m(), but for %parallel matrices it may be smaller. * - * To figure out which elements - * exactly are stored locally, - * use local_range(). + * To figure out which elements exactly are stored locally, use + * local_range(). */ unsigned int local_size () const; /** - * Return a pair of indices - * indicating which rows of - * this matrix are stored - * locally. The first number is - * the index of the first row - * stored, the second the index - * of the one past the last one - * that is stored locally. If - * this is a sequential matrix, - * then the result will be the - * pair (0,m()), otherwise it - * will be a pair (i,i+n), - * where + * Return a pair of indices indicating which rows of this matrix are + * stored locally. The first number is the index of the first row stored, + * the second the index of the one past the last one that is stored + * locally. If this is a sequential matrix, then the result will be the + * pair (0,m()), otherwise it will be a pair (i,i+n), where * n=local_size(). */ std::pair local_range () const; /** - * Return whether @p index is - * in the local range or not, - * see also local_range(). + * Return whether @p index is in the local range or not, see also + * local_range(). */ bool in_local_range (const size_type index) const; /** - * Return the number of nonzero - * elements of this matrix. + * Return the number of nonzero elements of this matrix. */ size_type n_nonzero_elements () const; /** - * Number of entries in a - * specific row. + * Number of entries in a specific row. */ unsigned int row_length (const size_type row) const; /** - * Returns the state of the matrix, - * i.e., whether compress() needs to - * be called after an operation - * requiring data exchange. A call to - * compress() is also needed when the - * method set() has been called (even - * when working in serial). + * Returns the state of the matrix, i.e., whether compress() needs to be + * called after an operation requiring data exchange. A call to compress() + * is also needed when the method set() has been called (even when working + * in serial). */ bool is_compressed () const; /** - * Determine an estimate for the memory - * consumption (in bytes) of this - * object. Note that only the memory - * reserved on the current processor is - * returned in case this is called in - * an MPI-based program. + * Determine an estimate for the memory consumption (in bytes) of this + * object. Note that only the memory reserved on the current processor is + * returned in case this is called in an MPI-based program. */ size_type memory_consumption () const; @@ -1321,75 +1024,48 @@ namespace TrilinosWrappers //@{ /** - * This operator assigns a scalar to - * a matrix. Since this does usually - * not make much sense (should we set - * all matrix entries to this value? - * Only the nonzero entries of the - * sparsity pattern?), this operation - * is only allowed if the actual - * value to be assigned is zero. This - * operator only exists to allow for - * the obvious notation - * matrix=0, 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 matrix=0, 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: *
      - *
    • 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. - *
    • 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. + *
    • 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. + *
    • 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. *
    * - * In both cases, this function - * compresses the data structures and - * allows the resulting matrix to be - * used in all other operations like - * matrix-vector products. This is a - * collective operation, i.e., it - * needs to be run on all processors - * when used in %parallel. + * In both cases, this function compresses the data structures and allows + * the resulting matrix to be used in all other operations like + * matrix-vector products. This is a collective operation, i.e., it needs + * to be run on all processors when used in %parallel. * * See @ref GlossCompress "Compressing distributed objects" * for more information. @@ -1402,69 +1078,42 @@ namespace TrilinosWrappers void compress () DEAL_II_DEPRECATED; /** - * Set the element (i,j) - * 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,j) 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 into the sparse - * matrix locations given by - * indices. In other words, - * this function writes the elements - * in full_matrix into the - * calling matrix, using the - * local-to-global indexing specified - * by indices 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 into the sparse matrix + * locations given by indices. In other words, this function + * writes the elements in full_matrix into the calling matrix, + * using the local-to-global indexing specified by indices 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 - * elide_zero_values can be - * used to specify whether zero - * values should be inserted anyway - * or they should be filtered - * away. The default value is - * false, i.e., even zero + * The optional parameter elide_zero_values can be used to + * specify whether zero values should be inserted anyway or they should be + * filtered away. The default value is false, i.e., even zero * values are inserted/replaced. */ void set (const std::vector &indices, @@ -1472,11 +1121,9 @@ namespace TrilinosWrappers const bool elide_zero_values = false); /** - * Same function as before, but now - * including the possibility to use - * rectangular full_matrices and - * different local-to-global indexing - * on rows and columns, respectively. + * Same function as before, but now including the possibility to use + * rectangular full_matrices and different local-to-global indexing on + * rows and columns, respectively. */ void set (const std::vector &row_indices, const std::vector &col_indices, @@ -1484,30 +1131,18 @@ namespace TrilinosWrappers const bool elide_zero_values = false); /** - * Set several elements in the - * specified row of the matrix with - * column indices as given by - * col_indices 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 col_indices to the respective value. * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be inserted anyway - * or they should be filtered - * away. The default value is - * false, 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 elide_zero_values can be used to + * specify whether zero values should be inserted anyway or they should be + * filtered away. The default value is false, i.e., even zero * values are inserted/replaced. */ void set (const size_type row, @@ -1516,30 +1151,18 @@ namespace TrilinosWrappers const bool elide_zero_values = false); /** - * Set several elements to values - * given by values in a - * given row in columns given by - * col_indices into the sparse - * matrix. + * Set several elements to values given by values 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 - * elide_zero_values can be - * used to specify whether zero - * values should be inserted anyway - * or they should be filtered - * away. The default value is - * false, i.e., even zero + * The optional parameter elide_zero_values can be used to + * specify whether zero values should be inserted anyway or they should be + * filtered away. The default value is false, i.e., even zero * values are inserted/replaced. */ void set (const size_type row, @@ -1549,70 +1172,44 @@ namespace TrilinosWrappers const bool elide_zero_values = false); /** - * Add @p value to the element - * (i,j). + * Add @p value to the element (i,j). * - * Just as the respective call in - * deal.II SparseMatrix - * 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 - * value is not a finite - * number an exception is thrown. + * Just as the respective call in deal.II SparseMatrix 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 value 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 into sparse - * matrix locations given by - * indices. In other words, - * this function adds the elements in - * full_matrix to the - * respective entries in calling - * matrix, using the local-to-global - * indexing specified by - * indices 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 into sparse matrix + * locations given by indices. In other words, this function adds + * the elements in full_matrix to the respective entries in + * calling matrix, using the local-to-global indexing specified by + * indices 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 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 - * 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 - * elide_zero_values 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 true, - * i.e., zero values won't be added - * into the matrix. + * The optional parameter elide_zero_values 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 + * true, i.e., zero values won't be added into the matrix. */ void add (const std::vector &indices, const FullMatrix &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 &row_indices, const std::vector &col_indices, @@ -1620,30 +1217,17 @@ namespace TrilinosWrappers const bool elide_zero_values = true); /** - * Set several elements in the - * specified row of the matrix with - * column indices as given by - * col_indices to the - * respective value. - * - * Just as the respective call in - * deal.II SparseMatrix - * 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 - * elide_zero_values 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 true, - * 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 col_indices to the respective value. + * + * Just as the respective call in deal.II SparseMatrix 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 elide_zero_values 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 + * true, i.e., zero values won't be added into the matrix. */ void add (const size_type row, const std::vector &col_indices, @@ -1651,29 +1235,17 @@ namespace TrilinosWrappers const bool elide_zero_values = true); /** - * Add an array of values given by - * values 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 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 - * elide_zero_values 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 true, - * i.e., zero values won't be added - * into the matrix. + * Add an array of values given by values 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 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 elide_zero_values 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 + * true, i.e., zero values won't be added into the matrix. */ void add (const size_type row, const size_type n_cols, @@ -1683,14 +1255,12 @@ namespace TrilinosWrappers const bool col_indices_are_sorted = false); /** - * Multiply the entire matrix - * by a fixed factor. + * Multiply the entire matrix by a fixed factor. */ SparseMatrix &operator *= (const TrilinosScalar factor); /** - * Divide the entire matrix by - * a fixed factor. + * Divide the entire matrix by a fixed factor. */ SparseMatrix &operator /= (const TrilinosScalar factor); @@ -1701,84 +1271,48 @@ namespace TrilinosWrappers void copy_from (const SparseMatrix &source); /** - * Add matrix scaled by - * factor to this matrix, - * i.e. the matrix - * factor*matrix is added to - * this. 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 matrix scaled by factor to this matrix, i.e. the + * matrix factor*matrix is added to this. 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 row 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 row 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 &rows, const TrilinosScalar new_diag_value = 0); /** - * Make an in-place transpose - * of a matrix. + * Make an in-place transpose of a matrix. */ void transpose (); @@ -1789,77 +1323,40 @@ namespace TrilinosWrappers //@{ /** - * Return the value of the - * entry (i,j). 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,j). 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,j). 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,j). 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 ith - * row. This function throws an - * error if the matrix is not - * quadratic and it also throws - * an error if (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 ith row. This function + * throws an error if the matrix is not quadratic and it also throws an + * error if (i,i) is not element of the local matrix. See also the + * comment in trilinos_sparse_matrix.cc. */ TrilinosScalar diag_element (const size_type i) const; @@ -1975,155 +1472,86 @@ namespace TrilinosWrappers const VectorType &src) const; /** - * Return the square of the norm - * of the vector $v$ with respect - * to the norm induced by this - * matrix, i.e., - * $\left(v,Mv\right)$. This is - * useful, e.g. in the finite - * element context, where the - * $L_2$ norm of a function - * equals the matrix norm with - * respect to the mass matrix of - * the vector representing the - * nodal values of the finite - * element function. - * - * Obviously, the matrix needs to - * be quadratic for this - * operation. - * - * The implementation of this - * function is not as efficient - * as the one in the @p - * SparseMatrix class used in - * deal.II (i.e. the original - * one, not the Trilinos wrapper - * class) since Trilinos doesn't - * support this operation and - * needs a temporary vector. - * - * Note that both vectors have to - * be distributed vectors - * generated using the same Map - * as was used for the matrix in - * case you work on a distributed - * memory architecture, using the - * interface in the - * TrilinosWrappers::VectorBase - * class (or one of the two - * derived classes Vector and - * MPI::Vector). - * - * In case of a localized Vector, - * this function will only work - * when running on one processor, - * since the matrix object is - * inherently - * distributed. Otherwise, and - * exception will be thrown. + * Return the square of the norm of the vector $v$ with respect to the + * norm induced by this matrix, i.e., $\left(v,Mv\right)$. This is useful, + * e.g. in the finite element context, where the $L_2$ norm of a function + * equals the matrix norm with respect to the mass matrix of the vector + * representing the nodal values of the finite element function. + * + * Obviously, the matrix needs to be quadratic for this operation. + * + * The implementation of this function is not as efficient as the one in + * the @p SparseMatrix class used in deal.II (i.e. the original one, not + * the Trilinos wrapper class) since Trilinos doesn't support this + * operation and needs a temporary vector. + * + * Note that both vectors have to be distributed vectors generated using + * the same Map as was used for the matrix in case you work on a + * distributed memory architecture, using the interface in the + * TrilinosWrappers::VectorBase class (or one of the two derived classes + * Vector and MPI::Vector). + * + * In case of a localized Vector, this function will only work when + * running on one processor, since the matrix object is inherently + * distributed. Otherwise, and exception will be thrown. */ TrilinosScalar matrix_norm_square (const VectorBase &v) const; /** - * Compute the matrix scalar - * product $\left(u,Mv\right)$. - * - * The implementation of this - * function is not as efficient - * as the one in the @p - * SparseMatrix class used in - * deal.II (i.e. the original - * one, not the Trilinos - * wrapper class) since - * Trilinos doesn't support - * this operation and needs a - * temporary vector. - * - * Note that both vectors have to - * be distributed vectors - * generated using the same Map - * as was used for the matrix in - * case you work on a distributed - * memory architecture, using the - * interface in the - * TrilinosWrappers::VectorBase - * class (or one of the two - * derived classes Vector and - * MPI::Vector). - * - * In case of a localized Vector, - * this function will only work - * when running on one processor, - * since the matrix object is - * inherently - * distributed. Otherwise, and - * exception will be thrown. + * Compute the matrix scalar product $\left(u,Mv\right)$. + * + * The implementation of this function is not as efficient as the one in + * the @p SparseMatrix class used in deal.II (i.e. the original one, not + * the Trilinos wrapper class) since Trilinos doesn't support this + * operation and needs a temporary vector. + * + * Note that both vectors have to be distributed vectors generated using + * the same Map as was used for the matrix in case you work on a + * distributed memory architecture, using the interface in the + * TrilinosWrappers::VectorBase class (or one of the two derived classes + * Vector and MPI::Vector). + * + * In case of a localized Vector, this function will only work when + * running on one processor, since the matrix object is inherently + * distributed. Otherwise, and exception will be thrown. */ TrilinosScalar matrix_scalar_product (const VectorBase &u, const VectorBase &v) const; /** - * Compute the residual of an - * equation Mx=b, where - * the residual is defined to - * be r=b-Mx. Write the - * residual into @p dst. The - * l2 norm of - * the residual vector is - * returned. - * - * Source x and - * destination dst 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 Mx=b, where the residual is + * defined to be r=b-Mx. Write the residual into @p dst. The + * l2 norm of the residual vector is returned. + * + * Source x and destination dst 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 C = A * B, - * or, if an optional vector argument - * is given, C = A * diag(V) * - * B, where diag(V) - * defines a diagonal matrix with the - * vector entries. + * Perform the matrix-matrix multiplication C = A * B, or, if an + * optional vector argument is given, C = A * diag(V) * B, where + * diag(V) defines a diagonal matrix with the vector entries. * - * This function assumes that the - * calling matrix A and - * B have compatible - * sizes. The size of C will - * be set within this function. + * This function assumes that the calling matrix A and B + * have compatible sizes. The size of C will be set within this + * function. * - * The content as well as the sparsity - * pattern of the matrix C will be - * changed by this function, so make - * sure that the sparsity pattern is - * not used somewhere else in your - * program. This is an expensive - * operation, so think twice before you - * use this function. + * The content as well as the sparsity pattern of the matrix C will be + * changed by this function, so make sure that the sparsity pattern is not + * used somewhere else in your program. This is an expensive operation, so + * think twice before you use this function. */ void mmult (SparseMatrix &C, const SparseMatrix &B, @@ -2131,30 +1559,20 @@ namespace TrilinosWrappers /** - * Perform the matrix-matrix - * multiplication with the transpose of - * this, i.e., C = - * AT * B, or, if an - * optional vector argument is given, - * C = AT * diag(V) * - * B, where diag(V) - * defines a diagonal matrix with the - * vector entries. + * Perform the matrix-matrix multiplication with the transpose of + * this, i.e., C = AT * B, or, if an optional + * vector argument is given, C = AT * diag(V) * B, + * where diag(V) defines a diagonal matrix with the vector + * entries. * - * This function assumes that the - * calling matrix A and - * B have compatible - * sizes. The size of C will - * be set within this function. + * This function assumes that the calling matrix A and B + * have compatible sizes. The size of C will be set within this + * function. * - * The content as well as the sparsity - * pattern of the matrix C will be - * changed by this function, so make - * sure that the sparsity pattern is - * not used somewhere else in your - * program. This is an expensive - * operation, so think twice before you - * use this function. + * The content as well as the sparsity pattern of the matrix C will be + * changed by this function, so make sure that the sparsity pattern is not + * used somewhere else in your program. This is an expensive operation, so + * think twice before you use this function. */ void Tmmult (SparseMatrix &C, const SparseMatrix &B, @@ -2167,46 +1585,27 @@ namespace TrilinosWrappers //@{ /** - * Return the - * l1-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 l1-norm of the matrix, that is $|M|_1= + * \max_{\mathrm{all\ columns\ } j} \sum_{\mathrm{all\ rows\ } i} + * |M_{ij}|$, (max. sum of columns). This is the natural matrix norm that + * is compatible to the l1-norm for vectors, i.e. $|Mv|_1 \leq |M|_1 + * |v|_1$. (cf. Haemmerlin-Hoffmann: Numerische Mathematik) */ TrilinosScalar l1_norm () const; /** - * Return the linfty-norm of the - * matrix, that is - * $|M|_\infty=\max_{\mathrm{all\ - * rows\ } i}\sum_{\mathrm{all\ - * columns\ } j} |M_{ij}|$, - * (max. sum of rows). This is - * the natural matrix norm that - * is compatible to the - * linfty-norm of vectors, i.e. - * $|Mv|_\infty \leq |M|_\infty - * |v|_\infty$. - * (cf. Haemmerlin-Hoffmann: - * Numerische Mathematik) + * Return the linfty-norm of the matrix, that is + * $|M|_\infty=\max_{\mathrm{all\ rows\ } i}\sum_{\mathrm{all\ columns\ } + * j} |M_{ij}|$, (max. sum of rows). This is the natural matrix norm that + * is compatible to the linfty-norm of vectors, i.e. $|Mv|_\infty \leq + * |M|_\infty |v|_\infty$. (cf. Haemmerlin-Hoffmann: Numerische + * Mathematik) */ TrilinosScalar linfty_norm () const; /** - * Return the frobenius norm of - * the matrix, i.e. the square - * root of the sum of squares - * of all entries in the - * matrix. + * Return the frobenius norm of the matrix, i.e. the square root of the + * sum of squares of all entries in the matrix. */ TrilinosScalar frobenius_norm () const; @@ -2217,60 +1616,44 @@ namespace TrilinosWrappers //@{ /** - * Return a const reference to the - * underlying Trilinos - * Epetra_CrsMatrix data. + * Return a const reference to the underlying Trilinos Epetra_CrsMatrix + * data. */ const Epetra_CrsMatrix &trilinos_matrix () const; /** - * Return a const reference to the - * underlying Trilinos - * Epetra_CrsGraph data that stores - * the sparsity pattern of the - * matrix. + * Return a const reference to the underlying Trilinos Epetra_CrsGraph + * data that stores the sparsity pattern of the matrix. */ const Epetra_CrsGraph &trilinos_sparsity_pattern () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * domain space of this matrix, i.e., - * the partitioning of the vectors - * this matrix has to be multiplied - * with. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the domain space of this matrix, i.e., the + * partitioning of the vectors this matrix has to be multiplied with. */ const Epetra_Map &domain_partitioner () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * range space of this matrix, i.e., - * the partitioning of the vectors - * that are result from matrix-vector + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the range space of this matrix, i.e., the + * partitioning of the vectors that are result from matrix-vector * products. */ const Epetra_Map &range_partitioner () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * matrix rows. Equal to the - * partitioning of the range. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the matrix rows. Equal to the partitioning of + * the range. */ const Epetra_Map &row_partitioner () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * matrix columns. This is in general - * not equal to the partitioner - * Epetra_Map for the domain because - * of overlap in the matrix. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the matrix columns. This is in general not + * equal to the partitioner Epetra_Map for the domain because of overlap + * in the matrix. */ const Epetra_Map &col_partitioner () const; //@} @@ -2280,8 +1663,7 @@ namespace TrilinosWrappers //@{ /** - * STL-like iterator with the - * first entry. + * STL-like iterator with the first entry. */ const_iterator begin () const; @@ -2291,40 +1673,27 @@ namespace TrilinosWrappers const_iterator end () const; /** - * STL-like iterator with the - * first entry of row @p r. + * STL-like iterator with the first entry of row @p r. * - * Note that if the given row - * is empty, i.e. does not - * contain any nonzero entries, - * then the iterator returned - * by this function equals - * end(r). 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 + * end(r). Note also that the iterator may not be dereferencable + * in that case. */ const_iterator begin (const size_type r) const; /** - * Final iterator of row - * r. 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 r. It points to the first element past + * the end of line @p r, or past the end of the entire sparsity pattern. * - * Note that the end iterator - * is not necessarily - * dereferencable. This is in - * particular the case if it is - * the end iterator for the - * last row of a matrix. + * Note that the end iterator is not necessarily dereferencable. This is + * in particular the case if it is the end iterator for the last row of a + * matrix. */ const_iterator end (const size_type r) const; /** - * STL-like iterator with the - * first entry. + * STL-like iterator with the first entry. */ iterator begin (); @@ -2334,34 +1703,22 @@ namespace TrilinosWrappers iterator end (); /** - * STL-like iterator with the - * first entry of row @p r. + * STL-like iterator with the first entry of row @p r. * - * Note that if the given row - * is empty, i.e. does not - * contain any nonzero entries, - * then the iterator returned - * by this function equals - * end(r). 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 + * end(r). Note also that the iterator may not be dereferencable + * in that case. */ iterator begin (const size_type r); /** - * Final iterator of row - * r. 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 r. It points to the first element past + * the end of line @p r, or past the end of the entire sparsity pattern. * - * Note that the end iterator - * is not necessarily - * dereferencable. This is in - * particular the case if it is - * the end iterator for the - * last row of a matrix. + * Note that the end iterator is not necessarily dereferencable. This is + * in particular the case if it is the end iterator for the last row of a + * matrix. */ iterator end (const size_type r); @@ -2372,27 +1729,18 @@ namespace TrilinosWrappers //@{ /** - * Abstract Trilinos object - * that helps view in ASCII - * other Trilinos - * objects. Currently this - * function is not - * implemented. TODO: Not + * Abstract Trilinos object that helps view in ASCII other Trilinos + * objects. Currently this function is not implemented. TODO: Not * implemented. */ void write_ascii (); /** - * Print the matrix to the given - * stream, using the format - * (line,col) value, 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 (line,col) + * value, i.e. one nonzero entry of the matrix per line. The optional + * flag outputs the sparsity pattern in Trilinos style, where the data is + * sorted according to the processor number when printed to the stream, as + * well as a summary of the matrix like the global size. */ void print (std::ostream &out, const bool write_extended_trilinos_info = false) const; @@ -2455,35 +1803,20 @@ namespace TrilinosWrappers protected: /** - * For some matrix storage - * formats, in particular for the - * PETSc distributed blockmatrices, - * set and add operations on - * individual elements can not be - * freely mixed. Rather, one has - * to synchronize operations when - * one wants to switch from - * setting elements to adding to - * elements. - * BlockMatrixBase automatically - * synchronizes the access by - * calling this helper function - * for each block. - * This function ensures that the - * matrix is in a state that - * allows adding elements; if it - * previously already was in this - * state, the function does - * nothing. + * For some matrix storage formats, in particular for the PETSc distributed + * blockmatrices, set and add operations on individual elements can not be + * freely mixed. Rather, one has to synchronize operations when one wants + * to switch from setting elements to adding to elements. BlockMatrixBase + * automatically synchronizes the access by calling this helper function + * for each block. This function ensures that the matrix is in a state + * that allows adding elements; if it previously already was in this state, + * the function does nothing. */ void prepare_add(); /** - * Same as prepare_add() but - * prepare the matrix for setting - * elements if the representation - * of elements in this class - * requires such an operation. + * Same as prepare_add() but prepare the matrix for setting elements if the + * representation of elements in this class requires such an operation. */ void prepare_set(); @@ -2492,60 +1825,46 @@ namespace TrilinosWrappers private: /** - * Pointer to the user-supplied - * Epetra Trilinos mapping of - * the matrix columns that - * assigns parts of the matrix - * to the individual processes. + * Pointer to the user-supplied Epetra Trilinos mapping of the matrix + * columns that assigns parts of the matrix to the individual processes. */ std_cxx1x::shared_ptr 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 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 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; }; @@ -2965,46 +2284,6 @@ namespace TrilinosWrappers - inline - void - SparseMatrix::compress (::dealii::VectorOperation::values operation) - { - - Epetra_CombineMode mode = last_action; - if (last_action == Zero) - { - if ((operation==::dealii::VectorOperation::add) || - (operation==::dealii::VectorOperation::unknown)) - mode = Add; - else if (operation==::dealii::VectorOperation::insert) - mode = Insert; - } - else - { - Assert( - ((last_action == Add) && (operation!=::dealii::VectorOperation::insert)) - || - ((last_action == Insert) && (operation!=::dealii::VectorOperation::add)), - ExcMessage("operation and argument to compress() do not match")); - } - - // flush buffers - int ierr; - ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(), - true, mode); - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - - ierr = matrix->OptimizeStorage (); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - - last_action = Zero; - - compressed = true; - } - - - inline void SparseMatrix::compress () @@ -3023,19 +2302,17 @@ namespace TrilinosWrappers const int ierr = matrix->PutScalar(d); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + if (nonlocal_matrix.get() != 0) + nonlocal_matrix->PutScalar(d); return *this; } - // Inline the set() and add() - // functions, since they will be - // called frequently, and the - // compiler can optimize away - // some unnecessary loops when - // the sizes are given at - // compile time. + // Inline the set() and add() functions, since they will be called + // frequently, and the compiler can optimize away some unnecessary loops + // when the sizes are given at compile time. inline void SparseMatrix::set (const size_type i, @@ -3208,7 +2485,7 @@ namespace TrilinosWrappers // standard Insert/ReplaceGlobalValues function. Nevertheless, the way // we call it is the fastest one (any other will lead to repeated // allocation and deallocation of memory in order to call the function - // we already use, which is very unefficient if writing one element at + // we already use, which is very inefficient if writing one element at // a time). compressed = false; @@ -3222,6 +2499,19 @@ namespace TrilinosWrappers if (ierr > 0) ierr = 0; } + else if (nonlocal_matrix.get() != 0) + { + // this is the case when we have explicitly set the off-processor + // rows and want to create a separate matrix object for them (to + // retain thread-safety) + Assert (nonlocal_matrix->RowMap().LID(static_cast(row)) != -1, + ExcMessage("Attempted to write into off-processor matrix row " + "that has not be specified as being writable upon " + "initialization")); + ierr = nonlocal_matrix->ReplaceGlobalValues(row, n_columns, + col_value_ptr, + col_index_ptr); + } else ierr = matrix->ReplaceGlobalValues (1, (TrilinosWrappers::types::int_type *)&row, @@ -3413,6 +2703,20 @@ namespace TrilinosWrappers col_value_ptr, col_index_ptr); } + else if (nonlocal_matrix.get() != 0) + { + compressed = false; + // this is the case when we have explicitly set the off-processor rows + // and want to create a separate matrix object for them (to retain + // thread-safety) + Assert (nonlocal_matrix->RowMap().LID(static_cast(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 diff --git a/deal.II/include/deal.II/lac/trilinos_sparsity_pattern.h b/deal.II/include/deal.II/lac/trilinos_sparsity_pattern.h index 3519bf8af1..01d6d16676 100644 --- a/deal.II/include/deal.II/lac/trilinos_sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/trilinos_sparsity_pattern.h @@ -93,20 +93,17 @@ namespace TrilinosWrappers Accessor (const Accessor &a); /** - * Row number of the element - * represented by this object. + * Row number of the element represented by this object. */ size_type row() const; /** - * Index in row of the element - * represented by this object. + * Index in row of the element represented by this object. */ size_type index() const; /** - * Column number of the element - * represented by this object. + * Column number of the element represented by this object. */ size_type column() const; @@ -142,45 +139,28 @@ namespace TrilinosWrappers size_type a_index; /** - * Cache where we store the - * column indices of the - * present row. This is - * necessary, since Trilinos - * makes access to the elements - * of its matrices rather hard, - * and it is much more - * efficient to copy all column - * entries of a row once when - * we enter it than repeatedly - * asking Trilinos for - * individual ones. This also - * makes some sense since it is - * likely that we will access - * them sequentially anyway. + * Cache where we store the column indices of the present row. This is + * necessary, since Trilinos makes access to the elements of its + * matrices rather hard, and it is much more efficient to copy all + * column entries of a row once when we enter it than repeatedly asking + * Trilinos for individual ones. This also makes some sense since it is + * likely that we will access them sequentially anyway. * - * In order to make copying of - * iterators/accessor of - * acceptable performance, we - * keep a shared pointer to - * these entries so that more - * than one accessor can access - * this data if necessary. + * In order to make copying of iterators/accessor of acceptable + * performance, we keep a shared pointer to these entries so that more + * than one accessor can access this data if necessary. */ std_cxx1x::shared_ptr > colnum_cache; /** - * Discard the old row caches - * (they may still be used by - * other accessors) and - * generate new ones for the - * row pointed to presently by + * Discard the old row caches (they may still be used by other + * accessors) and generate new ones for the row pointed to presently by * this accessor. */ void visit_present_row (); /** - * Make enclosing class a - * friend. + * Make enclosing class a friend. */ friend class Iterator; }; @@ -199,10 +179,8 @@ namespace TrilinosWrappers typedef dealii::types::global_dof_index size_type; /** - * Constructor. Create an - * iterator into the matrix @p - * matrix for the given row and - * the index within it. + * Constructor. Create an iterator into the matrix @p matrix for the + * given row and the index within it. */ Iterator (const SparsityPattern *sparsity_pattern, const size_type row, @@ -234,9 +212,8 @@ namespace TrilinosWrappers const Accessor *operator-> () const; /** - * Comparison. True, if both - * iterators point to the same - * matrix position. + * Comparison. True, if both iterators point to the same matrix + * position. */ bool operator == (const Iterator &) const; @@ -246,11 +223,8 @@ namespace TrilinosWrappers bool operator != (const Iterator &) const; /** - * Comparison operator. Result - * is true if either the first - * row number is smaller or if - * the row numbers are equal - * and the first index is + * Comparison operator. Result is true if either the first row number is + * smaller or if the row numbers are equal and the first index is * smaller. */ bool operator < (const Iterator &) const; @@ -266,8 +240,7 @@ namespace TrilinosWrappers private: /** - * Store an object of the - * accessor class. + * Store an object of the accessor class. */ Accessor accessor; @@ -317,92 +290,61 @@ namespace TrilinosWrappers */ //@{ /** - * Default constructor. Generates an - * empty (zero-size) sparsity - * pattern. + * Default constructor. Generates an empty (zero-size) sparsity pattern. */ SparsityPattern (); /** - * Generate a sparsity pattern that is - * completely stored locally, having - * $m$ rows and $n$ columns. The - * resulting matrix will be completely + * Generate a sparsity pattern that is completely stored locally, having + * $m$ rows and $n$ columns. The resulting matrix will be completely * stored locally, too. * - * It is possible to specify the - * number of columns entries per row - * using the optional @p - * n_entries_per_row - * argument. However, this value does - * not need to be accurate or even - * given at all, since one does - * usually not have this kind of - * information before building the - * sparsity pattern (the usual case - * when the function - * DoFTools::make_sparsity_pattern() - * is called). The entries are - * allocated dynamically in a similar - * manner as for the deal.II - * CompressedSparsityPattern - * classes. However, a good estimate - * will reduce the setup time of the - * sparsity pattern. + * It is possible to specify the number of columns entries per row using + * the optional @p n_entries_per_row argument. However, this value does + * not need to be accurate or even given at all, since one does usually + * not have this kind of information before building the sparsity pattern + * (the usual case when the function DoFTools::make_sparsity_pattern() is + * called). The entries are allocated dynamically in a similar manner as + * for the deal.II CompressedSparsityPattern classes. However, a good + * estimate will reduce the setup time of the sparsity pattern. */ SparsityPattern (const size_type m, const size_type n, const size_type n_entries_per_row = 0); /** - * Generate a sparsity pattern that is - * completely stored locally, having - * $m$ rows and $n$ columns. The - * resulting matrix will be completely + * Generate a sparsity pattern that is completely stored locally, having + * $m$ rows and $n$ columns. The resulting matrix will be completely * stored locally, too. * - * The vector - * n_entries_per_row - * specifies the number of entries in - * each row (an information usually - * not available, though). + * The vector n_entries_per_row 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 &n_entries_per_row); /** - * Copy constructor. Sets the calling - * sparsity pattern to be the same as + * Copy constructor. Sets the calling sparsity pattern to be the same as * the input sparsity pattern. */ SparsityPattern (const SparsityPattern &input_sparsity_pattern); /** - * Destructor. Made virtual so that - * one can use pointers to this - * class. + * Destructor. Made virtual so that one can use pointers to this class. */ virtual ~SparsityPattern (); /** - * Initialize a sparsity pattern that - * is completely stored locally, - * having $m$ rows and $n$ - * columns. The resulting matrix will - * be completely stored locally. + * Initialize a sparsity pattern that is completely stored locally, having + * $m$ rows and $n$ columns. The resulting matrix will be completely + * stored locally. * - * The number of columns entries per - * row is specified as the maximum - * number of entries argument. This - * does not need to be an accurate - * number since the entries are - * allocated dynamically in a similar - * manner as for the deal.II - * CompressedSparsityPattern classes, - * but a good estimate will reduce - * the setup time of the sparsity - * pattern. + * The number of columns entries per row is specified as the maximum + * number of entries argument. This does not need to be an accurate + * number since the entries are allocated dynamically in a similar manner + * as for the deal.II CompressedSparsityPattern classes, but a good + * estimate will reduce the setup time of the sparsity pattern. */ void reinit (const size_type m, @@ -410,16 +352,12 @@ namespace TrilinosWrappers const size_type n_entries_per_row = 0); /** - * Initialize a sparsity pattern that - * is completely stored locally, - * having $m$ rows and $n$ columns. The - * resulting matrix will be - * completely stored locally. + * Initialize a sparsity pattern that is completely stored locally, having + * $m$ rows and $n$ columns. The resulting matrix will be completely + * stored locally. * - * The vector - * n_entries_per_row - * specifies the number of entries in - * each row. + * The vector n_entries_per_row specifies the number of entries + * in each row. */ void reinit (const size_type m, @@ -427,65 +365,46 @@ namespace TrilinosWrappers const std::vector &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 void copy_from (const SparsityType &nontrilinos_sparsity_pattern); /** - * Copy operator. This operation is - * only allowed for empty objects, to - * avoid potentially very costly - * operations automatically - * synthesized by the compiler. Use - * copy_from() instead if you know - * that you really want to copy a - * sparsity pattern with non-trivial - * content. + * Copy operator. This operation is only allowed for empty objects, to + * avoid potentially very costly operations automatically synthesized by + * the compiler. Use copy_from() instead if you know that you really want + * to copy a sparsity pattern with non-trivial content. */ SparsityPattern &operator = (const SparsityPattern &input_sparsity_pattern); /** - * Release all memory and return to a - * state just like after having - * called the default constructor. + * Release all memory and return to a state just like after having called + * the default constructor. * - * This is a collective operation - * that needs to be called on all - * processors in order to avoid a - * dead lock. + * This is a collective operation that needs to be called on all + * processors in order to avoid a dead lock. */ void clear (); /** - * In analogy to our own - * SparsityPattern class, this - * function compresses the sparsity - * pattern and allows the resulting - * pattern to be used for actually - * generating a (Trilinos-based) - * matrix. This function also - * exchanges non-local data that - * might have accumulated during the - * addition of new elements. This - * function must therefore be called - * once the structure is fixed. This - * is a collective operation, i.e., - * it needs to be run on all - * processors when used in parallel. + * In analogy to our own SparsityPattern class, this function compresses + * the sparsity pattern and allows the resulting pattern to be used for + * actually generating a (Trilinos-based) matrix. This function also + * exchanges non-local data that might have accumulated during the + * addition of new elements. This function must therefore be called once + * the structure is fixed. This is a collective operation, i.e., it needs + * to be run on all processors when used in parallel. */ void compress (); //@} @@ -495,79 +414,47 @@ namespace TrilinosWrappers //@{ /** - * Constructor for a square sparsity - * pattern using an Epetra_map for - * the description of the %parallel - * partitioning. Moreover, the number - * of nonzero entries in the rows of - * the sparsity pattern can be - * specified. Note that this number - * does not need to be exact, and it - * is allowed that the actual - * sparsity structure has more - * nonzero entries than specified in - * the constructor (the usual case - * when the function - * DoFTools::make_sparsity_pattern() - * is called). However it is still - * advantageous to provide good - * estimates here since a good value - * will avoid repeated allocation of - * memory, which considerably - * increases the performance when - * creating the sparsity pattern. + * Constructor for a square sparsity pattern using an Epetra_map for the + * description of the %parallel partitioning. Moreover, the number of + * nonzero entries in the rows of the sparsity pattern can be + * specified. Note that this number does not need to be exact, and it is + * allowed that the actual sparsity structure has more nonzero entries + * than specified in the constructor (the usual case when the function + * DoFTools::make_sparsity_pattern() is called). However it is still + * advantageous to provide good estimates here since a good value will + * avoid repeated allocation of memory, which considerably increases the + * performance when creating the sparsity pattern. */ SparsityPattern (const Epetra_Map ¶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 &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, @@ -575,74 +462,47 @@ namespace TrilinosWrappers const size_type n_entries_per_row = 0); /** - * This constructor is similar to the - * one above, but it now takes two - * different Epetra maps for rows and - * columns. This interface is meant to - * be used for generating rectangular - * matrices, where one map specifies - * the %parallel distribution of rows - * and the second one specifies the - * distribution of degrees of freedom - * associated with matrix columns. This - * second map is however not used for - * the distribution of the columns - * themselves – rather, all - * column elements of a row are stored - * on the same processor. The vector - * n_entries_per_row 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 n_entries_per_row 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 &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 @@ -650,31 +510,19 @@ namespace TrilinosWrappers const std::vector &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 * n_entries_per_row. */ void @@ -683,27 +531,15 @@ namespace TrilinosWrappers const size_type n_entries_per_row = 0); /** - * This reinit function is similar to - * the one above, but it now takes - * two different Epetra maps for rows - * and columns. This interface is - * meant to be used for generating - * rectangular matrices, where one - * map specifies the %parallel - * distribution of rows and the - * second one specifies the - * distribution of degrees of freedom - * associated with matrix - * columns. This second map is - * however not used for the - * distribution of the columns - * themselves – rather, all - * column elements of a row are - * stored on the same processor. The - * vector n_entries_per_row - * 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 n_entries_per_row specifies the + * number of entries in each row of the newly generated matrix. */ void reinit (const Epetra_Map &row_parallel_partitioning, @@ -711,19 +547,12 @@ namespace TrilinosWrappers const std::vector &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 void @@ -733,19 +562,12 @@ namespace TrilinosWrappers const bool exchange_data = false); /** - * Reinit function. Takes one of the - * deal.II sparsity patterns and a - * %parallel partitioning of the rows - * and columns for initializing the - * current Trilinos sparsity - * pattern. The optional argument @p - * exchange_data can be used for - * reinitialization with a sparsity - * pattern that is not fully - * constructed. This feature is only - * implemented for input sparsity - * patterns of type - * CompressedSimpleSparsityPattern. + * Reinit function. Takes one of the deal.II sparsity patterns and a + * %parallel partitioning of the rows and columns for initializing the + * current Trilinos sparsity pattern. The optional argument @p + * exchange_data can be used for reinitialization with a sparsity pattern + * that is not fully constructed. This feature is only implemented for + * input sparsity patterns of type CompressedSimpleSparsityPattern. */ template void @@ -759,74 +581,47 @@ namespace TrilinosWrappers //@{ /** - * Constructor for a square sparsity - * pattern using an IndexSet and an - * MPI communicator for the - * description of the %parallel - * partitioning. Moreover, the number - * of nonzero entries in the rows of - * the sparsity pattern can be - * specified. Note that this number - * does not need to be exact, and it - * is even allowed that the actual - * sparsity structure has more - * nonzero entries than specified in - * the constructor. However it is - * still advantageous to provide good - * estimates here since a good value - * will avoid repeated allocation of - * memory, which considerably - * increases the performance when - * creating the sparsity pattern. + * Constructor for a square sparsity pattern using an IndexSet and an MPI + * communicator for the description of the %parallel + * partitioning. Moreover, the number of nonzero entries in the rows of + * the sparsity pattern can be specified. Note that this number does not + * need to be exact, and it is even allowed that the actual sparsity + * structure has more nonzero entries than specified in the + * constructor. However it is still advantageous to provide good estimates + * here since a good value will avoid repeated allocation of memory, which + * considerably increases the performance when creating the sparsity + * pattern. */ SparsityPattern (const IndexSet ¶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 &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, @@ -835,24 +630,15 @@ namespace TrilinosWrappers const size_type n_entries_per_row = 0); /** - * This constructor is similar to the - * one above, but it now takes two - * different index sets for rows and - * columns. This interface is meant to - * be used for generating rectangular - * matrices, where one map specifies - * the %parallel distribution of rows - * and the second one specifies the - * distribution of degrees of freedom - * associated with matrix columns. This - * second map is however not used for - * the distribution of the columns - * themselves – rather, all - * column elements of a row are stored - * on the same processor. The vector - * n_entries_per_row 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 n_entries_per_row specifies the + * number of entries in each row of the newly generated matrix. */ SparsityPattern (const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, @@ -860,30 +646,51 @@ namespace TrilinosWrappers const std::vector &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, @@ -891,21 +698,13 @@ namespace TrilinosWrappers const size_type n_entries_per_row = 0); /** - * Same as before, but now use the - * exact number of nonzeros in each m - * row. Since we know the number of - * elements in the sparsity pattern - * exactly in this case, we can - * already allocate the right amount - * of memory, which makes process of - * adding entries to the sparsity - * pattern considerably - * faster. However, this is a rather - * unusual situation, since knowing - * the number of entries in each row - * is usually connected to knowing - * the indices of nonzero entries, - * which the sparsity pattern is + * Same as before, but now use the exact number of nonzeros in each m + * row. Since we know the number of elements in the sparsity pattern + * exactly in this case, we can already allocate the right amount of + * memory, which makes process of adding entries to the sparsity pattern + * considerably faster. However, this is a rather unusual situation, since + * knowing the number of entries in each row is usually connected to + * knowing the indices of nonzero entries, which the sparsity pattern is * designed to describe. */ void @@ -914,32 +713,19 @@ namespace TrilinosWrappers const std::vector &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 * n_entries_per_row. */ void @@ -949,11 +735,40 @@ namespace TrilinosWrappers const size_type n_entries_per_row = 0); /** - * Same as before, but now using a - * vector n_entries_per_row - * 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 n_entries_per_row for + * specifying the number of entries in each row of the sparsity pattern. */ void reinit (const IndexSet &row_parallel_partitioning, @@ -962,20 +777,13 @@ namespace TrilinosWrappers const std::vector &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 void @@ -986,19 +794,12 @@ namespace TrilinosWrappers const bool exchange_data = false); /** - * Reinit function. Takes one of the - * deal.II sparsity patterns and a - * %parallel partitioning of the rows - * and columns for initializing the - * current Trilinos sparsity - * pattern. The optional argument @p - * exchange_data can be used for - * reinitialization with a sparsity - * pattern that is not fully - * constructed. This feature is only - * implemented for input sparsity - * patterns of type - * CompressedSimpleSparsityPattern. + * Reinit function. Takes one of the deal.II sparsity patterns and a + * %parallel partitioning of the rows and columns for initializing the + * current Trilinos sparsity pattern. The optional argument @p + * exchange_data can be used for reinitialization with a sparsity pattern + * that is not fully constructed. This feature is only implemented for + * input sparsity patterns of type CompressedSimpleSparsityPattern. */ template void @@ -1013,119 +814,88 @@ namespace TrilinosWrappers //@{ /** - * Returns the state of the sparsity - * pattern, i.e., whether compress() - * needs to be called after an - * operation requiring data - * exchange. + * Returns the state of the sparsity pattern, i.e., whether compress() + * needs to be called after an operation requiring data exchange. */ bool is_compressed () const; /** - * Gives the maximum number of - * entries per row on the current - * processor. + * Gives the maximum number of entries per row on the current processor. */ unsigned int max_entries_per_row () const; /** - * Return the number of rows in this - * sparsity pattern. + * Return the number of rows in this sparsity pattern. */ size_type n_rows () const; /** - * Return the number of columns in - * this sparsity pattern. + * Return the number of columns in this sparsity pattern. */ size_type n_cols () const; /** - * Return the local dimension of the - * sparsity pattern, i.e. the number - * of rows stored on the present MPI - * process. In the sequential case, - * this number is the same as - * n_rows(), but for parallel - * matrices it may be smaller. + * Return the local dimension of the sparsity pattern, i.e. the number of + * rows stored on the present MPI process. In the sequential case, this + * number is the same as n_rows(), but for parallel matrices it may be + * smaller. * - * To figure out which elements - * exactly are stored locally, - * use local_range(). + * To figure out which elements exactly are stored locally, use + * local_range(). */ unsigned int local_size () const; /** - * Return a pair of indices - * indicating which rows of this - * sparsity pattern are stored - * locally. The first number is the - * index of the first row stored, the - * second the index of the one past - * the last one that is stored - * locally. If this is a sequential - * matrix, then the result will be - * the pair (0,n_rows()), otherwise - * it will be a pair (i,i+n), where + * Return a pair of indices indicating which rows of this sparsity pattern + * are stored locally. The first number is the index of the first row + * stored, the second the index of the one past the last one that is + * stored locally. If this is a sequential matrix, then the result will be + * the pair (0,n_rows()), otherwise it will be a pair (i,i+n), where * n=local_size(). */ std::pair local_range () const; /** - * Return whether @p index is - * in the local range or not, - * see also local_range(). + * Return whether @p index is in the local range or not, see also + * local_range(). */ 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,j) exists in the - * sparsity pattern (i.e., it may be - * non-zero) or not. + * Return whether the index (i,j) exists in the sparsity pattern + * (i.e., it may be non-zero) or not. */ bool exists (const size_type i, const size_type j) const; /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. Currently not - * implemented for this class. + * Determine an estimate for the memory consumption (in bytes) of this + * object. Currently not implemented for this class. */ std::size_t memory_consumption () const; @@ -1135,16 +905,14 @@ namespace TrilinosWrappers */ //@{ /** - * Add the element (i,j) to - * the sparsity pattern. + * Add the element (i,j) 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 void add_entries (const size_type row, @@ -1158,61 +926,44 @@ namespace TrilinosWrappers //@{ /** - * Return a const reference to the - * underlying Trilinos - * Epetra_CrsGraph data that stores - * the sparsity pattern. + * Return a const reference to the underlying Trilinos Epetra_CrsGraph + * data that stores the sparsity pattern. */ const Epetra_FECrsGraph &trilinos_sparsity_pattern () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the parallel - * partitioning of the domain space - * of this sparsity pattern, i.e., - * the partitioning of the vectors - * matrices based on this sparsity - * pattern are multiplied with. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the parallel partitioning of the domain space of this sparsity + * pattern, i.e., the partitioning of the vectors matrices based on this + * sparsity pattern are multiplied with. */ const Epetra_Map &domain_partitioner () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * range space of this sparsity - * pattern, i.e., the partitioning of - * the vectors that are result from + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the range space of this sparsity pattern, + * i.e., the partitioning of the vectors that are result from * matrix-vector products. */ const Epetra_Map &range_partitioner () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * sparsity pattern rows. Equal to - * the partitioning of the range. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the sparsity pattern rows. Equal to the + * partitioning of the range. */ const Epetra_Map &row_partitioner () const; /** - * Return a const reference to the - * underlying Trilinos Epetra_Map - * that sets the partitioning of the - * sparsity pattern columns. This is - * in general not equal to the - * partitioner Epetra_Map for the - * domain because of overlap in the - * matrix. + * Return a const reference to the underlying Trilinos Epetra_Map that + * sets the partitioning of the sparsity pattern columns. This is in + * general not equal to the partitioner Epetra_Map for the domain because + * of overlap in the matrix. */ const Epetra_Map &col_partitioner () const; /** - * Return a const reference to - * the communicator used for - * this object. + * Return a const reference to the communicator used for this object. */ const Epetra_Comm &trilinos_communicator () const; //@} @@ -1222,8 +973,7 @@ namespace TrilinosWrappers //@{ /** - * STL-like iterator with the - * first entry. + * STL-like iterator with the first entry. */ const_iterator begin () const; @@ -1233,34 +983,22 @@ namespace TrilinosWrappers const_iterator end () const; /** - * STL-like iterator with the - * first entry of row @p r. + * STL-like iterator with the first entry of row @p r. * - * Note that if the given row - * is empty, i.e. does not - * contain any nonzero entries, - * then the iterator returned - * by this function equals - * end(r). 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 + * end(r). Note also that the iterator may not be dereferencable + * in that case. */ const_iterator begin (const size_type r) const; /** - * Final iterator of row - * r. 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 r. It points to the first element past + * the end of line @p r, or past the end of the entire sparsity pattern. * - * Note that the end iterator - * is not necessarily - * dereferencable. This is in - * particular the case if it is - * the end iterator for the - * last row of a matrix. + * Note that the end iterator is not necessarily dereferencable. This is + * in particular the case if it is the end iterator for the last row of a + * matrix. */ const_iterator end (const size_type r) const; @@ -1271,55 +1009,35 @@ namespace TrilinosWrappers //@{ /** - * Abstract Trilinos object - * that helps view in ASCII - * other Trilinos - * objects. Currently this - * function is not - * implemented. TODO: Not + * Abstract Trilinos object that helps view in ASCII other Trilinos + * objects. Currently this function is not implemented. TODO: Not * implemented. */ void write_ascii (); /** - * Print (the locally owned part of) - * the sparsity pattern to the given - * stream, using the format - * (line,col). 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 (line,col). 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 gnuplot - * understands and which can be used - * to plot the sparsity pattern in a - * graphical way. The format consists - * of pairs i j 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 - * (0,0) element is in the - * top left rather than in the bottom + * Print the sparsity of the matrix in a format that gnuplot + * understands and which can be used to plot the sparsity pattern in a + * graphical way. The format consists of pairs i j 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 (0,0) 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 - * plot command. + * Print the sparsity pattern in gnuplot by setting the data style to dots + * or points and use the plot command. */ void print_gnuplot (std::ostream &out) const; @@ -1371,30 +1089,31 @@ namespace TrilinosWrappers private: /** - * Pointer to the user-supplied - * Epetra Trilinos mapping of - * the matrix columns that - * assigns parts of the matrix - * to the individual processes. + * Pointer to the user-supplied Epetra Trilinos mapping of the matrix + * columns that assigns parts of the matrix to the individual processes. */ std_cxx1x::shared_ptr 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 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 nonlocal_graph; + + friend class SparseMatrix; friend class SparsityPatternIterators::Accessor; friend class SparsityPatternIterators::Iterator; }; @@ -1696,9 +1415,23 @@ namespace TrilinosWrappers const int n_cols = static_cast(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(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(row)) != -1, + ExcMessage("Attempted to write into off-processor matrix row " + "that has not be specified as being writable upon " + "initialization")); + ierr = nonlocal_graph->InsertGlobalIndices (row, n_cols, col_index_ptr); + } + else + ierr = graph->InsertGlobalIndices + (1, (TrilinosWrappers::types::int_type *)&row, n_cols, col_index_ptr); AssertThrow (ierr >= 0, ExcTrilinosError(ierr)); } @@ -1757,157 +1490,6 @@ namespace TrilinosWrappers return graph->RangeMap().Comm(); } - - - inline - SparsityPattern::SparsityPattern (const IndexSet ¶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 &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 &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 &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 &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 - 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 - 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 } diff --git a/deal.II/source/lac/block_sparsity_pattern.cc b/deal.II/source/lac/block_sparsity_pattern.cc index 695ae57da6..0536414ef5 100644 --- a/deal.II/source/lac/block_sparsity_pattern.cc +++ b/deal.II/source/lac/block_sparsity_pattern.cc @@ -785,6 +785,45 @@ namespace TrilinosWrappers this->collect_sizes(); } + + + void + BlockSparsityPattern::reinit (const std::vector &row_parallel_partitioning, + const std::vector &col_parallel_partitioning, + const MPI_Comm &communicator) + { + dealii::BlockSparsityPatternBase:: + reinit(row_parallel_partitioning.size(), + col_parallel_partitioning.size()); + for (size_type i=0; iblock(i,j).reinit(row_parallel_partitioning[i], + col_parallel_partitioning[j], + communicator); + this->collect_sizes(); + } + + + + void + BlockSparsityPattern::reinit (const std::vector &row_parallel_partitioning, + const std::vector &col_parallel_partitioning, + const std::vector &writable_rows, + const MPI_Comm &communicator) + { + AssertDimension(writable_rows.size(), row_parallel_partitioning.size()); + dealii::BlockSparsityPatternBase:: + reinit(row_parallel_partitioning.size(), + col_parallel_partitioning.size()); + for (size_type i=0; iblock(i,j).reinit(row_parallel_partitioning[i], + col_parallel_partitioning[j], + writable_rows[i], + communicator); + this->collect_sizes(); + } + } #endif diff --git a/deal.II/source/lac/trilinos_sparse_matrix.cc b/deal.II/source/lac/trilinos_sparse_matrix.cc index b8f2956a74..665c7ba151 100644 --- a/deal.II/source/lac/trilinos_sparse_matrix.cc +++ b/deal.II/source/lac/trilinos_sparse_matrix.cc @@ -26,6 +26,7 @@ # include # include +# include # include # include # include @@ -574,13 +575,18 @@ namespace TrilinosWrappers { matrix.reset (); - // reinit with a (parallel) Trilinos - // sparsity pattern. + // reinit with a (parallel) Trilinos sparsity pattern. column_space_map.reset (new Epetra_Map (sparsity_pattern.domain_partitioner())); matrix.reset (new Epetra_FECrsMatrix (Copy, sparsity_pattern.trilinos_sparsity_pattern(), false)); + + if (sparsity_pattern.nonlocal_graph.get() != 0) + { + nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph)); + } compress(); + last_action = Zero; } @@ -770,6 +776,56 @@ namespace TrilinosWrappers + inline + void + SparseMatrix::compress (::dealii::VectorOperation::values operation) + { + + Epetra_CombineMode mode = last_action; + if (last_action == Zero) + { + if ((operation==::dealii::VectorOperation::add) || + (operation==::dealii::VectorOperation::unknown)) + mode = Add; + else if (operation==::dealii::VectorOperation::insert) + mode = Insert; + } + else + { + Assert( + ((last_action == Add) && (operation!=::dealii::VectorOperation::insert)) + || + ((last_action == Insert) && (operation!=::dealii::VectorOperation::add)), + ExcMessage("operation and argument to compress() do not match")); + } + + // flush buffers + int ierr; + if (nonlocal_matrix.get() != 0) + { + nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap()); + Epetra_Export exporter(nonlocal_matrix->RowMap(), matrix->RowMap()); + ierr = matrix->Export(*nonlocal_matrix, exporter, mode); + AssertThrow(ierr == 0, ExcTrilinosError(ierr)); + ierr = matrix->FillComplete(*column_space_map, matrix->RowMap()); + nonlocal_matrix->PutScalar(0); + } + else + ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(), + true, mode); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = matrix->OptimizeStorage (); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + last_action = Zero; + + compressed = true; + } + + + void SparseMatrix::clear () { @@ -779,6 +835,7 @@ namespace TrilinosWrappers column_space_map.reset (new Epetra_Map (0, 0, Utilities::Trilinos::comm_self())); matrix.reset (new Epetra_FECrsMatrix(View, *column_space_map, 0)); + nonlocal_matrix.reset(); matrix->FillComplete(); diff --git a/deal.II/source/lac/trilinos_sparsity_pattern.cc b/deal.II/source/lac/trilinos_sparsity_pattern.cc index 24f4618092..85ca1111fc 100644 --- a/deal.II/source/lac/trilinos_sparsity_pattern.cc +++ b/deal.II/source/lac/trilinos_sparsity_pattern.cc @@ -24,6 +24,8 @@ # include # include +# include + DEAL_II_NAMESPACE_OPEN namespace TrilinosWrappers @@ -250,6 +252,80 @@ 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 &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 &n_entries_per_row) + : + compressed (false) + { + Epetra_Map row_map = + row_parallel_partitioning.make_trilinos_map (communicator, false); + Epetra_Map col_map = + col_parallel_partitioning.make_trilinos_map (communicator, false); + reinit (row_map, col_map, n_entries_per_row); + } + + + + SparsityPattern:: + SparsityPattern (const IndexSet &row_parallel_partitioning, + const IndexSet &col_parallel_partitioning, + const IndexSet &writable_rows, + const MPI_Comm &communicator, + const size_type n_max_entries_per_row) + : + compressed (false) + { + reinit (row_parallel_partitioning, col_parallel_partitioning, + writable_rows, communicator, n_max_entries_per_row); + } + + + SparsityPattern::~SparsityPattern () {} @@ -282,6 +358,12 @@ namespace TrilinosWrappers const Epetra_Map &input_col_map, const size_type n_entries_per_row) { + Assert(input_row_map.IsOneToOne(), + ExcMessage("Row map must be 1-to-1, i.e., no overlap between " + "the maps of different processors.")); + Assert(input_col_map.IsOneToOne(), + ExcMessage("Column map must be 1-to-1, i.e., no overlap between " + "the maps of different processors.")); graph.reset (); column_space_map.reset (new Epetra_Map (input_col_map)); compressed = false; @@ -355,6 +437,122 @@ namespace TrilinosWrappers + void + SparsityPattern::reinit (const IndexSet ¶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 &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 &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 + 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 + 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 void SparsityPattern::reinit (const Epetra_Map &input_map, @@ -509,6 +707,8 @@ namespace TrilinosWrappers *column_space_map, 0)); graph->FillComplete(); + nonlocal_graph.reset(); + compressed = true; } @@ -519,9 +719,31 @@ namespace TrilinosWrappers { int ierr; Assert (column_space_map.get() != 0, ExcInternalError()); - ierr = graph->GlobalAssemble (*column_space_map, - static_cast(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(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(graph->RangeMap())); + } + else + ierr = graph->GlobalAssemble (*column_space_map, + static_cast(graph->RangeMap()), + true); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); @@ -849,6 +1071,55 @@ namespace TrilinosWrappers const dealii::CompressedSimpleSparsityPattern &, bool); + + + template void + SparsityPattern::reinit (const IndexSet &, + const dealii::SparsityPattern &, + const MPI_Comm &, + bool); + template void + SparsityPattern::reinit (const IndexSet &, + const dealii::CompressedSparsityPattern &, + const MPI_Comm &, + bool); + template void + SparsityPattern::reinit (const IndexSet &, + const dealii::CompressedSetSparsityPattern &, + const MPI_Comm &, + bool); + template void + SparsityPattern::reinit (const IndexSet &, + const dealii::CompressedSimpleSparsityPattern &, + const MPI_Comm &, + bool); + + + template void + SparsityPattern::reinit (const IndexSet &, + const IndexSet &, + const dealii::SparsityPattern &, + const MPI_Comm &, + bool); + template void + SparsityPattern::reinit (const IndexSet &, + const IndexSet &, + const dealii::CompressedSparsityPattern &, + const MPI_Comm &, + bool); + template void + SparsityPattern::reinit (const IndexSet &, + const IndexSet &, + const dealii::CompressedSetSparsityPattern &, + const MPI_Comm &, + bool); + template void + SparsityPattern::reinit (const IndexSet &, + const IndexSet &, + const dealii::CompressedSimpleSparsityPattern &, + const MPI_Comm &, + bool); + } DEAL_II_NAMESPACE_CLOSE diff --git a/tests/trilinos/assemble_matrix_parallel_01.cc b/tests/trilinos/assemble_matrix_parallel_01.cc new file mode 100644 index 0000000000..209fe9c025 --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_01.cc @@ -0,0 +1,451 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2009 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// same as deal.II/assemble_matrix_parallel_01, but for trilinos matrices + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("output"); + +using namespace dealii; + + +namespace Assembly +{ + namespace Scratch + { + template + struct Data + { + Data (const hp::FECollection &fe, + const hp::QCollection &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 hp_fe_values; + }; + } + + namespace Copy + { + struct Data + { + std::vector local_dof_indices; + FullMatrix local_matrix; + Vector local_rhs; + }; + } +} + +template +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::active_cell_iterator &cell, + Assembly::Scratch::Data &scratch, + Assembly::Copy::Data &data); + void copy_local_to_global (const Assembly::Copy::Data &data); + + std::vector + get_conflict_indices (typename hp::DoFHandler::active_cell_iterator const &cell) const; + + Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix reference_matrix; + TrilinosWrappers::SparseMatrix test_matrix; + + Vector reference_rhs; + Vector test_rhs; + + std::vector::active_cell_iterator> > graph; + + const unsigned int max_degree; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + double sum = 0; + for (unsigned int d=0; d +class RightHandSide : public Function +{ +public: + RightHandSide () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d +LaplaceProblem::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(degree)); + quadrature_collection.push_back (QGauss(degree+1)); + face_quadrature_collection.push_back (QGauss(degree+1)); + } + else + for (unsigned int degree=1; degree(degree)); + quadrature_collection.push_back (QGauss(degree+1)); + face_quadrature_collection.push_back (QGauss(degree+1)); + } +} + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + + + +template +std::vector +LaplaceProblem:: +get_conflict_indices (typename hp::DoFHandler::active_cell_iterator const &cell) const +{ + std::vector 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 +void LaplaceProblem::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(), + constraints); + constraints.close (); + + graph = GraphColoring::make_graph_coloring(dof_handler.begin_active(),dof_handler.end(), + static_cast + (typename hp::DoFHandler::active_cell_iterator const &)> > + (std_cxx1x::bind(&LaplaceProblem::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 +void +LaplaceProblem::local_assemble (const typename hp::DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::Data &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 &fe_values = scratch.hp_fe_values.get_present_fe_values (); + + const RightHandSide rhs_function; + + for (unsigned int q_point=0; + q_pointget_dof_indices (data.local_dof_indices); +} + + + +template +void +LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) +{ + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); +} + + + +template +void LaplaceProblem::assemble_reference () +{ + test_matrix = 0; + test_rhs = 0; + + Assembly::Copy::Data copy_data; + Assembly::Scratch::Data assembly_data(fe_collection, quadrature_collection); + + for (unsigned int color=0; color::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 +void LaplaceProblem::assemble_test () +{ + test_matrix = 0; + test_rhs = 0; + + WorkStream:: + run (graph, + std_cxx1x::bind (&LaplaceProblem:: + local_assemble, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&LaplaceProblem:: + copy_local_to_global, + this, + std_cxx1x::_1), + Assembly::Scratch::Data(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 +void LaplaceProblem::postprocess () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + for (unsigned int i=0; i::active_cell_iterator cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + cell->set_active_fe_index (rand() % fe_collection.size()); +} + + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation, 0, 1/* , + Point(), + 0.5, 1., (dim==3) ? 96 : 12, false*/); + triangulation.refine_global(2); + } + + setup_system (); + + assemble_reference (); + assemble_test (); + + if (cycle < 2) + postprocess (); + } +} + + + +int main (int argc, char **argv) +{ + deallog << std::setprecision (2); + logfile << std::setprecision (2); + deallog.attach(logfile); + deallog.depth_console(0); + + Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int); + + { + deallog.push("2d"); + LaplaceProblem<2> laplace_problem; + laplace_problem.run (); + deallog.pop(); + } + + { + deallog.push("3d"); + LaplaceProblem<3> laplace_problem; + laplace_problem.run (); + deallog.pop(); + } +} + diff --git a/tests/trilinos/assemble_matrix_parallel_01.output b/tests/trilinos/assemble_matrix_parallel_01.output new file mode 100644 index 0000000000..3b8bf4746d --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_01.output @@ -0,0 +1,13 @@ + +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:3d::error in matrix: 0 +DEAL:3d::error in vector: 0 +DEAL:3d::error in matrix: 0 +DEAL:3d::error in vector: 0 +DEAL:3d::error in matrix: 0 +DEAL:3d::error in vector: 0 diff --git a/tests/trilinos/assemble_matrix_parallel_02.cc b/tests/trilinos/assemble_matrix_parallel_02.cc new file mode 100644 index 0000000000..5dfdc8ad2c --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_02.cc @@ -0,0 +1,437 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2009 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// tests parallel assembly of Trilinos matrices when also using MPI + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("output"); + +using namespace dealii; + + +namespace Assembly +{ + namespace Scratch + { + template + struct Data + { + Data (const FiniteElement &fe, + const Quadrature &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 fe_values; + }; + } + + namespace Copy + { + struct Data + { + std::vector local_dof_indices; + FullMatrix local_matrix; + Vector local_rhs; + }; + } +} + +template +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::active_cell_iterator> &cell, + Assembly::Scratch::Data &scratch, + Assembly::Copy::Data &data); + void copy_local_to_global (const Assembly::Copy::Data &data); + + std::vector + get_conflict_indices (FilteredIterator::active_cell_iterator> const &cell) const; + + parallel::distributed::Triangulation triangulation; + + DoFHandler dof_handler; + FE_Q fe; + QGauss quadrature; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix reference_matrix; + TrilinosWrappers::SparseMatrix test_matrix; + + TrilinosWrappers::MPI::Vector reference_rhs; + TrilinosWrappers::MPI::Vector test_rhs; + + std::vector::active_cell_iterator> > > graph; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + double sum = 0; + for (unsigned int d=0; d +class RightHandSide : public Function +{ +public: + RightHandSide () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d +LaplaceProblem::LaplaceProblem () + : + triangulation (MPI_COMM_WORLD), + dof_handler (triangulation), + fe (1), + quadrature(fe.degree+1) +{ +} + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + + + +template +std::vector +LaplaceProblem:: +get_conflict_indices (FilteredIterator::active_cell_iterator> const &cell) const +{ + std::vector 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 +void LaplaceProblem::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(), + constraints); + constraints.close (); + + typedef FilteredIterator::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 + (FilteredIterator::active_cell_iterator> const &)> > + (std_cxx1x::bind(&LaplaceProblem::get_conflict_indices, this,std_cxx1x::_1))); + + TrilinosWrappers::SparsityPattern csp; + IndexSet locally_owned = dof_handler.locally_owned_dofs(), relevant_set; + DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set); + csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + reference_matrix.reinit (csp); + test_matrix.reinit (csp); + + reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); + test_rhs.reinit (reference_rhs); +} + + + +template +void +LaplaceProblem::local_assemble (const FilteredIterator::active_cell_iterator> &cell, + Assembly::Scratch::Data &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 &fe_values = scratch.fe_values; + + const RightHandSide rhs_function; + + for (unsigned int q_point=0; + q_pointget_dof_indices (data.local_dof_indices); +} + + + +template +void +LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) +{ + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); +} + + + +template +void LaplaceProblem::assemble_reference () +{ + test_matrix = 0; + test_rhs = 0; + + Assembly::Copy::Data copy_data; + Assembly::Scratch::Data assembly_data(fe, quadrature); + + for (unsigned int color=0; color::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 +void LaplaceProblem::assemble_test () +{ + test_matrix = 0; + test_rhs = 0; + + WorkStream:: + run (graph, + std_cxx1x::bind (&LaplaceProblem:: + local_assemble, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&LaplaceProblem:: + copy_local_to_global, + this, + std_cxx1x::_1), + Assembly::Scratch::Data(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 +void LaplaceProblem::postprocess () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + for (unsigned int i=0; i +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation, 0, 1/* , + Point(), + 0.5, 1., (dim==3) ? 96 : 12, false*/); + triangulation.refine_global(6); + } + + setup_system (); + + assemble_reference (); + assemble_test (); + + if (cycle < 2) + postprocess (); + } +} + + + +int main (int argc, char **argv) +{ + deallog << std::setprecision (2); + logfile << std::setprecision (2); + deallog.attach(logfile); + deallog.depth_console(0); + + Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int); + + { + deallog.push("2d"); + LaplaceProblem<2> laplace_problem; + laplace_problem.run (); + deallog.pop(); + } +} + diff --git a/tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output b/tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output new file mode 100644 index 0000000000..6e5eabc7ba --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_02.mpirun=1.output @@ -0,0 +1,7 @@ + +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 diff --git a/tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output b/tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output new file mode 100644 index 0000000000..6e5eabc7ba --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_02.mpirun=4.output @@ -0,0 +1,7 @@ + +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 diff --git a/tests/trilinos/assemble_matrix_parallel_03.cc b/tests/trilinos/assemble_matrix_parallel_03.cc new file mode 100644 index 0000000000..65a82e152e --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_03.cc @@ -0,0 +1,440 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2009 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Same test as assemble_matrix_parallel_02, but calling compress() several +// times on test matrix + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("output"); + +using namespace dealii; + + +namespace Assembly +{ + namespace Scratch + { + template + struct Data + { + Data (const FiniteElement &fe, + const Quadrature &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 fe_values; + }; + } + + namespace Copy + { + struct Data + { + std::vector local_dof_indices; + FullMatrix local_matrix; + Vector local_rhs; + }; + } +} + +template +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::active_cell_iterator> &cell, + Assembly::Scratch::Data &scratch, + Assembly::Copy::Data &data); + void copy_local_to_global (const Assembly::Copy::Data &data); + + std::vector + get_conflict_indices (FilteredIterator::active_cell_iterator> const &cell) const; + + parallel::distributed::Triangulation triangulation; + + DoFHandler dof_handler; + FE_Q fe; + QGauss quadrature; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix reference_matrix; + TrilinosWrappers::SparseMatrix test_matrix; + + TrilinosWrappers::MPI::Vector reference_rhs; + TrilinosWrappers::MPI::Vector test_rhs; + + std::vector::active_cell_iterator> > > graph; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + double sum = 0; + for (unsigned int d=0; d +class RightHandSide : public Function +{ +public: + RightHandSide () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d +LaplaceProblem::LaplaceProblem () + : + triangulation (MPI_COMM_WORLD), + dof_handler (triangulation), + fe (1), + quadrature(fe.degree+1) +{ +} + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + + + +template +std::vector +LaplaceProblem:: +get_conflict_indices (FilteredIterator::active_cell_iterator> const &cell) const +{ + std::vector 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 +void LaplaceProblem::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(), + constraints); + constraints.close (); + + typedef FilteredIterator::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 + (FilteredIterator::active_cell_iterator> const &)> > + (std_cxx1x::bind(&LaplaceProblem::get_conflict_indices, this,std_cxx1x::_1))); + + TrilinosWrappers::SparsityPattern csp; + IndexSet locally_owned = dof_handler.locally_owned_dofs(), relevant_set; + DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set); + csp.reinit(locally_owned, locally_owned, relevant_set, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + csp.compress(); + reference_matrix.reinit (csp); + test_matrix.reinit (csp); + + reference_rhs.reinit (locally_owned, MPI_COMM_WORLD); + test_rhs.reinit (reference_rhs); +} + + + +template +void +LaplaceProblem::local_assemble (const FilteredIterator::active_cell_iterator> &cell, + Assembly::Scratch::Data &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 &fe_values = scratch.fe_values; + + const RightHandSide rhs_function; + + for (unsigned int q_point=0; + q_pointget_dof_indices (data.local_dof_indices); +} + + + +template +void +LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) +{ + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); +} + + + +template +void LaplaceProblem::assemble_reference () +{ + test_matrix = 0; + test_rhs = 0; + + Assembly::Copy::Data copy_data; + Assembly::Scratch::Data assembly_data(fe, quadrature); + + for (unsigned int color=0; color::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 +void LaplaceProblem::assemble_test () +{ + test_matrix = 0; + test_rhs = 0; + + WorkStream:: + run (graph, + std_cxx1x::bind (&LaplaceProblem:: + local_assemble, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&LaplaceProblem:: + copy_local_to_global, + this, + std_cxx1x::_1), + Assembly::Scratch::Data(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 +void LaplaceProblem::postprocess () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + for (unsigned int i=0; i +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation, + Point(), + 0.5, 1., (dim==3) ? 96 : 12, false); + triangulation.refine_global(3); + } + + setup_system (); + + assemble_reference (); + assemble_test (); + + if (cycle < 2) + postprocess (); + } +} + + + +int main (int argc, char **argv) +{ + deallog << std::setprecision (2); + logfile << std::setprecision (2); + deallog.attach(logfile); + deallog.depth_console(0); + + Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int); + + { + deallog.push("2d"); + LaplaceProblem<2> laplace_problem; + laplace_problem.run (); + deallog.pop(); + } +} + diff --git a/tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output b/tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output new file mode 100644 index 0000000000..6e5eabc7ba --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_03.mpirun=4.output @@ -0,0 +1,7 @@ + +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 diff --git a/tests/trilinos/assemble_matrix_parallel_04.cc b/tests/trilinos/assemble_matrix_parallel_04.cc new file mode 100644 index 0000000000..90ffbde2ff --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_04.cc @@ -0,0 +1,463 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2009 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Same as assemble_matrix_parallel_02, but using block matrices + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("output"); + +using namespace dealii; + + +namespace Assembly +{ + namespace Scratch + { + template + struct Data + { + Data (const FiniteElement &fe, + const Quadrature &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 fe_values; + }; + } + + namespace Copy + { + struct Data + { + std::vector local_dof_indices; + FullMatrix local_matrix; + Vector local_rhs; + }; + } +} + +template +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::active_cell_iterator> &cell, + Assembly::Scratch::Data &scratch, + Assembly::Copy::Data &data); + void copy_local_to_global (const Assembly::Copy::Data &data); + + std::vector + get_conflict_indices (FilteredIterator::active_cell_iterator> const &cell) const; + + parallel::distributed::Triangulation triangulation; + + DoFHandler dof_handler; + FESystem fe; + QGauss quadrature; + + ConstraintMatrix constraints; + + TrilinosWrappers::BlockSparseMatrix reference_matrix; + TrilinosWrappers::BlockSparseMatrix test_matrix; + + TrilinosWrappers::MPI::BlockVector reference_rhs; + TrilinosWrappers::MPI::BlockVector test_rhs; + + std::vector::active_cell_iterator> > > graph; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function (2) {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + double sum = 0; + for (unsigned int d=0; d +class RightHandSide : public Function +{ +public: + RightHandSide () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d +LaplaceProblem::LaplaceProblem () + : + triangulation (MPI_COMM_WORLD), + dof_handler (triangulation), + fe (FE_Q(1),1, FE_Q(2),1), + quadrature(3) +{ +} + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + + + +template +std::vector +LaplaceProblem:: +get_conflict_indices (FilteredIterator::active_cell_iterator> const &cell) const +{ + std::vector 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 +void LaplaceProblem::setup_system () +{ + reference_matrix.clear(); + test_matrix.clear(); + dof_handler.distribute_dofs (fe); + std::vector 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(), + constraints); + constraints.close (); + + typedef FilteredIterator::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 + (FilteredIterator::active_cell_iterator> const &)> > + (std_cxx1x::bind(&LaplaceProblem::get_conflict_indices, this,std_cxx1x::_1))); + + TrilinosWrappers::BlockSparsityPattern csp(2,2); + std::vector 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 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 +void +LaplaceProblem::local_assemble (const FilteredIterator::active_cell_iterator> &cell, + Assembly::Scratch::Data &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 &fe_values = scratch.fe_values; + + const RightHandSide 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_pointget_dof_indices (data.local_dof_indices); +} + + + +template +void +LaplaceProblem::copy_local_to_global (const Assembly::Copy::Data &data) +{ + constraints.distribute_local_to_global(data.local_matrix, data.local_rhs, + data.local_dof_indices, + test_matrix, test_rhs); +} + + + +template +void LaplaceProblem::assemble_reference () +{ + test_matrix = 0; + test_rhs = 0; + + Assembly::Copy::Data copy_data; + Assembly::Scratch::Data assembly_data(fe, quadrature); + + for (unsigned int color=0; color::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 +void LaplaceProblem::assemble_test () +{ + test_matrix = 0; + test_rhs = 0; + + WorkStream:: + run (graph, + std_cxx1x::bind (&LaplaceProblem:: + local_assemble, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&LaplaceProblem:: + copy_local_to_global, + this, + std_cxx1x::_1), + Assembly::Scratch::Data(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::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 +void LaplaceProblem::postprocess () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + for (unsigned int i=0; i +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_shell(triangulation, + Point(), + 0.5, 1., (dim==3) ? 96 : 12, false); +#ifdef DEBUG + triangulation.refine_global(3); +#else + triangulation.refine_global(5); +#endif + } + + setup_system (); + + assemble_reference (); + assemble_test (); + + if (cycle < 2) + postprocess (); + } +} + + + +int main (int argc, char **argv) +{ + deallog << std::setprecision (2); + logfile << std::setprecision (2); + deallog.attach(logfile); + deallog.depth_console(0); + + Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int); + + { + deallog.push("2d"); + LaplaceProblem<2> laplace_problem; + laplace_problem.run (); + deallog.pop(); + } +} + diff --git a/tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output b/tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output new file mode 100644 index 0000000000..6e5eabc7ba --- /dev/null +++ b/tests/trilinos/assemble_matrix_parallel_04.mpirun=4.output @@ -0,0 +1,7 @@ + +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 +DEAL:2d::error in matrix: 0 +DEAL:2d::error in vector: 0 -- 2.39.5