From 427753a8cceb14fcfd8b147366f995cdc5351d82 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sat, 27 Dec 2008 20:42:06 +0000 Subject: [PATCH] Added a function TrilinosWrappers::SparsityPattern that can generate distributed sparsity patterns for use with Trilinos sparse matrices. git-svn-id: https://svn.dealii.org/trunk@18021 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/source/dofs/dof_constraints.cc | 24 +- deal.II/deal.II/source/dofs/dof_tools.cc | 124 +- .../lac/include/lac/block_sparsity_pattern.h | 136 +- .../lac/include/lac/trilinos_sparse_matrix.h | 441 +++--- .../include/lac/trilinos_sparsity_pattern.h | 1272 +++++++++++++++++ deal.II/lac/source/block_sparsity_pattern.cc | 119 +- .../source/trilinos_block_sparse_matrix.cc | 29 +- deal.II/lac/source/trilinos_solver_block.cc | 2 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 80 +- .../lac/source/trilinos_sparsity_pattern.cc | 735 ++++++++++ 10 files changed, 2687 insertions(+), 275 deletions(-) create mode 100755 deal.II/lac/include/lac/trilinos_sparsity_pattern.h create mode 100755 deal.II/lac/source/trilinos_sparsity_pattern.cc diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 2421304424..80c251d5ee 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -2090,16 +2090,23 @@ add_entries_local_to_global (const std::vector &) const; - template void ConstraintMatrix:: add_entries_local_to_global (const std::vector &, CompressedSimpleSparsityPattern &, const bool, const Table<2,bool> &) const; +#ifdef DEAL_II_USE_TRILINOS +template void ConstraintMatrix:: +add_entries_local_to_global + (const std::vector &, + TrilinosWrappers::SparsityPattern &, + const bool, + const Table<2,bool> &) const; +#endif template void ConstraintMatrix:: add_entries_local_to_global (const std::vector &, - BlockSparsityPattern &, + BlockSparsityPattern &, const bool, const Table<2,bool> &) const; template void ConstraintMatrix:: @@ -2109,14 +2116,21 @@ add_entries_local_to_global (const std::vector &) const; template void ConstraintMatrix:: add_entries_local_to_global (const std::vector &, - BlockCompressedSetSparsityPattern &, + BlockCompressedSetSparsityPattern &, const bool, const Table<2,bool> &) const; - template void ConstraintMatrix:: add_entries_local_to_global (const std::vector &, - BlockCompressedSimpleSparsityPattern &, + BlockCompressedSimpleSparsityPattern &, const bool, const Table<2,bool> &) const; +#ifdef DEAL_II_USE_TRILINOS +template void ConstraintMatrix:: +add_entries_local_to_global + (const std::vector &, + TrilinosWrappers::BlockSparsityPattern &, + const bool, + const Table<2,bool> &) const; +#endif DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 9bba2f39f0..7fae1afa9f 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -33,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -66,7 +68,24 @@ DoFTools::make_sparsity_pattern (const DH &dof, dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); + + // In case we work with a distributed + // sparsity pattern of Trilinos type, we + // only have to do the work if the + // current cell is owned by the calling + // processor. Otherwise, just continue. for (; cell!=endc; ++cell) +#ifdef DEAL_II_USE_TRILINOS + if ((types_are_equal::value + || + types_are_equal::value) + && + cell->subdomain_id() != + Utilities::Trilinos::get_this_mpi_process(Utilities::Trilinos::comm_world())) + continue; + else + +#endif { const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; dofs_on_this_cell.resize (dofs_per_cell); @@ -156,7 +175,24 @@ DoFTools::make_sparsity_pattern ( std::vector dofs_on_this_cell(fe_collection.max_dofs_per_cell()); typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); - for (; cell!=endc; ++cell) + + // In case we work with a distributed + // sparsity pattern of Trilinos type, we + // only have to do the work if the + // current cell is owned by the calling + // processor. Otherwise, just continue. + for (; cell!=endc; ++cell) +#ifdef DEAL_II_USE_TRILINOS + if ((types_are_equal::value + || + types_are_equal::value) + && + cell->subdomain_id() != + Utilities::Trilinos::get_this_mpi_process(Utilities::Trilinos::comm_world())) + continue; + else + +#endif { const unsigned int fe_index = cell->active_fe_index(); @@ -5394,6 +5430,16 @@ DoFTools::make_sparsity_pattern, CompressedSimpleSparsityPattern &sparsity, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::SparsityPattern> +(const DoFHandler &dof, + TrilinosWrappers::SparsityPattern &sparsity, + const ConstraintMatrix &, + const bool); +#endif + template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> @@ -5422,6 +5468,15 @@ DoFTools::make_sparsity_pattern, BlockCompressedSimpleSparsityPattern &sparsity, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::BlockSparsityPattern> +(const DoFHandler &dof, + TrilinosWrappers::BlockSparsityPattern &sparsity, + const ConstraintMatrix &, + const bool); +#endif template void DoFTools::make_sparsity_pattern, @@ -5430,8 +5485,6 @@ DoFTools::make_sparsity_pattern, SparsityPattern &sparsity, const ConstraintMatrix &, const bool); - - template void DoFTools::make_sparsity_pattern, CompressedSparsityPattern> @@ -5439,7 +5492,6 @@ DoFTools::make_sparsity_pattern, CompressedSparsityPattern &sparsity, const ConstraintMatrix &, const bool); - template void DoFTools::make_sparsity_pattern, CompressedSetSparsityPattern> @@ -5447,7 +5499,6 @@ DoFTools::make_sparsity_pattern, CompressedSetSparsityPattern &sparsity, const ConstraintMatrix &, const bool); - template void DoFTools::make_sparsity_pattern, CompressedSimpleSparsityPattern> @@ -5455,7 +5506,15 @@ DoFTools::make_sparsity_pattern, CompressedSimpleSparsityPattern &sparsity, const ConstraintMatrix &, const bool); - +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::SparsityPattern> +(const hp::DoFHandler &dof, + TrilinosWrappers::SparsityPattern &sparsity, + const ConstraintMatrix &, + const bool); +#endif template void DoFTools::make_sparsity_pattern, @@ -5485,6 +5544,16 @@ DoFTools::make_sparsity_pattern, BlockCompressedSimpleSparsityPattern &sparsity, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::BlockSparsityPattern> +(const hp::DoFHandler &dof, + TrilinosWrappers::BlockSparsityPattern &sparsity, + const ConstraintMatrix &, + const bool); +#endif + template void @@ -5519,6 +5588,17 @@ DoFTools::make_sparsity_pattern, CompressedSimpleSparsityPattern&, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::SparsityPattern> +(const DoFHandler&, + const Table<2,Coupling>&, + TrilinosWrappers::SparsityPattern&, + const ConstraintMatrix &, + const bool); +#endif + template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> @@ -5551,6 +5631,16 @@ DoFTools::make_sparsity_pattern, BlockCompressedSimpleSparsityPattern&, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::BlockSparsityPattern> +(const DoFHandler&, + const Table<2,Coupling>&, + TrilinosWrappers::BlockSparsityPattern&, + const ConstraintMatrix &, + const bool); +#endif template void DoFTools::make_sparsity_pattern, @@ -5584,6 +5674,17 @@ DoFTools::make_sparsity_pattern, CompressedSimpleSparsityPattern&, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::SparsityPattern> +(const hp::DoFHandler&, + const Table<2,Coupling>&, + TrilinosWrappers::SparsityPattern&, + const ConstraintMatrix &, + const bool); +#endif + template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> @@ -5616,6 +5717,17 @@ DoFTools::make_sparsity_pattern, BlockCompressedSimpleSparsityPattern&, const ConstraintMatrix &, const bool); +#ifdef DEAL_II_USE_TRILINOS +template void +DoFTools::make_sparsity_pattern, + TrilinosWrappers::BlockSparsityPattern> +(const hp::DoFHandler&, + const Table<2,Coupling>&, + TrilinosWrappers::BlockSparsityPattern&, + const ConstraintMatrix &, + const bool); +#endif + template void diff --git a/deal.II/lac/include/lac/block_sparsity_pattern.h b/deal.II/lac/include/lac/block_sparsity_pattern.h index 2952413051..45b914cba8 100644 --- a/deal.II/lac/include/lac/block_sparsity_pattern.h +++ b/deal.II/lac/include/lac/block_sparsity_pattern.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -33,6 +34,12 @@ class BlockSparsityPattern; class BlockCompressedSparsityPattern; class BlockCompressedSimpleSparsityPattern; class BlockCompressedSetSparsityPattern; +#ifdef DEAL_II_USE_TRILINOS +namespace TrilinosWrappers +{ + class BlockSparsityPattern; +} +#endif /*! @addtogroup Sparsity *@{ @@ -687,8 +694,6 @@ typedef BlockCompressedSparsityPattern CompressedBlockSparsityPattern; * There are several, exchangeable variations of this class, see @ref Sparsity, * section 'Dynamic block sparsity patterns' for more information. * - * This class is used in @ref step_22 "step-22". - * * @author Wolfgang Bangerth, 2007 */ class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase @@ -776,6 +781,8 @@ class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase @@ -803,7 +810,7 @@ class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase + { + 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. + */ + 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. + */ + BlockSparsityPattern (const unsigned int n_rows, + const unsigned int 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. + */ + 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. + */ + BlockSparsityPattern (const std::vector& input_maps); + + /** + * 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]. + */ + void reinit (const std::vector< unsigned int > &row_block_sizes, + const std::vector< unsigned int > &col_block_sizes); + + /** + * 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& input_maps); + + + /** + * Allow the use of the reinit + * functions of the base class as + * well. + */ + using BlockSparsityPatternBase::reinit; + }; +} + +#endif + + /*@}*/ /*---------------------- Template functions -----------------------------------*/ diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 21e9cf65e2..0a646f8325 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -41,17 +41,13 @@ DEAL_II_NAMESPACE_OPEN - // forward declarations -class SparsityPattern; -class CompressedSparsityPattern; -class CompressedSetSparsityPattern; -class CompressedSimpleSparsityPattern; namespace TrilinosWrappers { // forward declarations class VectorBase; class SparseMatrix; + class SparsityPattern; namespace MatrixIterators { @@ -60,10 +56,11 @@ namespace TrilinosWrappers * over the elements of Trilinos matrices. The implementation of this * class is similar to the one for PETSc matrices. * - * Note that Trilinos does not give any guarantees as to the order of - * elements within each row. Note also that accessing the elements of - * a full matrix surprisingly only shows the nonzero elements of the - * matrix, not all elements. + * Note that Trilinos stores the elements within each row in ascending + * order. This is opposed to the deal.II sparse matrix style where the + * diagonal element (if it exists) is stored before all other values, and + * the PETSc sparse matrices, where one can't guarantee a certain order of + * the elements. * * @ingroup TrilinosWrappers * @author Martin Kronbichler, Wolfgang Bangerth, 2008 @@ -83,9 +80,9 @@ namespace TrilinosWrappers * access, a const matrix * pointer is sufficient. */ - Accessor (const SparseMatrix *matrix, - const unsigned int row, - const unsigned int index); + Accessor (const SparseMatrix *matrix, + const unsigned int row, + const unsigned int index); /** * Row number of the element @@ -199,7 +196,7 @@ namespace TrilinosWrappers * matrix for the given row and * the index within it. */ - const_iterator (const SparseMatrix *matrix, + const_iterator (const SparseMatrix *matrix, const unsigned int row, const unsigned int index); @@ -344,22 +341,23 @@ namespace TrilinosWrappers SparseMatrix (); /** - * Constructor using an - * Epetra_Map and a maximum - * number of nonzero matrix - * entries. 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. + * Constructor using an Epetra_Map + * and a maximum number of nonzero + * matrix entries. 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 should be no + * effect in the performance of + * matrix-vector products, since + * Trilinos wants to reorganize the + * matrix memory prior to use. */ SparseMatrix (const Epetra_Map &InputMap, const unsigned int n_max_entries_per_row); @@ -383,29 +381,28 @@ namespace TrilinosWrappers 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 takes - * care of the columns and the - * other one of the rows. 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 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 arragements when + * doing matrix-vector products with + * vectors based on that column map. * * The number of columns entries * per row is specified as the @@ -424,16 +421,22 @@ namespace TrilinosWrappers * to be used for generating * rectangular matrices, where one * map specifies the parallel - * distribution of rows and the - * second one specifies the number of - * columns in the total matrix. It - * also provides information for the + * 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, but is not used - * for the distribution of the - * columns – rather, all column - * elements of a row are stored on - * the same processor. The vector + * 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 @@ -472,6 +475,12 @@ namespace TrilinosWrappers const unsigned int n, const std::vector &n_entries_per_row); + /** + * Generate a matrix from a Trilinos + * sparsity pattern object. + */ + SparseMatrix (const SparsityPattern &InputSparsityPattern); + /** * Copy constructor. Sets the * calling matrix to be the same @@ -489,75 +498,64 @@ namespace TrilinosWrappers 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 + * 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 the - * 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. + * 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 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 + * 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. In a - * parallel run, it is currently - * necessary that each processor - * holds the sparsity_pattern - * structure because each - * processor sets its - * rows. + * internally, so it won't be needed + * any more after this call, in + * contrast to the deal.II own + * object. In a parallel run, 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 + * This is a collective operation + * that needs to be called on all * processors in order to avoid a * dead lock. */ @@ -566,18 +564,16 @@ namespace TrilinosWrappers const SparsityType &sparsity_pattern); /** - * 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. + * 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. * - * This is a - * collective operation that - * needs to be called on all + * This is a collective operation + * that needs to be called on all * processors in order to avoid a * dead lock. */ @@ -587,34 +583,42 @@ namespace TrilinosWrappers const SparsityType &sparsity_pattern); /** - * This function copies the - * content in - * sparse_matrix to - * the calling matrix. + * 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. + */ + void reinit (const SparsityPattern &sparsity_pattern); + + /** + * This function copies the content + * in sparse_matrix to the + * calling matrix. * - * This is a - * collective operation that - * needs to be called on all + * 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). + * 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). * - * This is a - * collective operation that - * needs to be called on all + * This is a collective operation + * that needs to be called on all * processors in order to avoid a * dead lock. */ @@ -622,27 +626,23 @@ namespace TrilinosWrappers const double drop_tolerance=1e-13); /** - * 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 + * 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 one. * - * This is a - * collective operation that - * needs to be called on all + * This is a collective operation + * that needs to be called on all * processors in order to avoid a * dead lock. */ @@ -651,18 +651,16 @@ namespace TrilinosWrappers const double drop_tolerance=1e-13); /** - * 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. + * 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. * - * This is a - * collective operation that - * needs to be called on all + * This is a collective operation + * that needs to be called on all * processors in order to avoid a * dead lock. */ @@ -679,58 +677,51 @@ namespace TrilinosWrappers void reinit (const Epetra_CrsMatrix &input_matrix); /** - * 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. + * 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 + * This is a collective operation + * that needs to be called on all * processors in order to avoid a * dead lock. */ void clear (); /** - * Trilinos matrices store their - * own sparsity patterns. So, in - * analogy to our own - * SparsityPattern class, this - * function compresses the - * sparsity pattern and allows - * the resulting matrix to be - * used in all other operations - * where before only assembly - * functions were allowed. This - * function must therefore be - * called once you have assembled - * the matrix. This is a + * Trilinos matrices store their own + * sparsity patterns. So, in analogy + * to our own SparsityPattern class, + * this function compresses the + * sparsity pattern and allows the + * resulting matrix to be used in all + * other operations where before only + * assembly functions were + * allowed. This function must + * therefore be called once you have + * assembled the matrix. This is a * collective operation, i.e., it - * needs to be run on all - * processors when used in - * parallel. + * needs to be run on all processors + * when used in parallel. */ void compress (); @@ -739,9 +730,9 @@ namespace TrilinosWrappers * i.e., whether compress() needs to * be called after an operation * requiring data exchange. A call to - * compress() is also after the - * method set() is called (even when - * working in serial). + * compress() is also needed when the + * method set() has been called (even + * when working in serial). */ bool is_compressed () const; //@} @@ -1741,7 +1732,7 @@ namespace TrilinosWrappers * to the individual processes. */ Epetra_Map col_map; - + /** * Trilinos doesn't allow to * mix additions to matrix @@ -2476,7 +2467,7 @@ DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_USE_TRILINOS -/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/ +/*----------------------- trilinos_sparse_matrix.h --------------------*/ #endif -/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/ +/*----------------------- trilinos_sparse_matrix.h --------------------*/ diff --git a/deal.II/lac/include/lac/trilinos_sparsity_pattern.h b/deal.II/lac/include/lac/trilinos_sparsity_pattern.h new file mode 100755 index 0000000000..11707a96db --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_sparsity_pattern.h @@ -0,0 +1,1272 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__trilinos_sparsity_pattern_h +#define __deal2__trilinos_sparsity_pattern_h + + +#include +#include +#include + +#include +#include +#include +#include + +#ifdef DEAL_II_USE_TRILINOS + +# include +# include +# ifdef DEAL_II_COMPILER_SUPPORTS_MPI +# include +# include "mpi.h" +# else +# include "Epetra_SerialComm.h" +# endif + +DEAL_II_NAMESPACE_OPEN + + // forward declarations +class SparsityPattern; +class CompressedSparsityPattern; +class CompressedSetSparsityPattern; +class CompressedSimpleSparsityPattern; + +namespace TrilinosWrappers +{ + // forward declarations + class SparsityPattern; + + namespace SparsityPatternIterators + { +/** + * STL conforming iterator. This class acts as an iterator walking + * over the elements of Trilinos sparsity pattern. + * + * @ingroup TrilinosWrappers + * @author Martin Kronbichler, Wolfgang Bangerth, 2008 + */ + class const_iterator + { + private: + /** + * Accessor class for iterators + */ + class Accessor + { + public: + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const SparsityPattern *sparsity_pattern, + const unsigned int row, + const unsigned int index); + + /** + * Row number of the element + * represented by this object. + */ + unsigned int row() const; + + /** + * Index in row of the element + * represented by this object. + */ + unsigned int index() const; + + /** + * Column number of the element + * represented by this object. + */ + unsigned int column() const; + + /** + * Exception + */ + DeclException0 (ExcBeyondEndOfSparsityPattern); + + /** + * Exception + */ + DeclException3 (ExcAccessToNonlocalRow, + int, int, int, + << "You tried to access row " << arg1 + << " of a distributed sparsity pattern, " + << " but only rows " << arg2 << " through " << arg3 + << " are stored locally and can be accessed."); + + private: + /** + * The matrix accessed. + */ + mutable SparsityPattern *sparsity_pattern; + + /** + * Current row number. + */ + unsigned int a_row; + + /** + * Current index in row. + */ + unsigned int 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. + * + * 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. + */ + boost::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 + * this accessor. + */ + void visit_present_row (); + + /** + * Make enclosing class a + * friend. + */ + friend class const_iterator; + }; + + public: + + /** + * Constructor. Create an + * iterator into the matrix @p + * matrix for the given row and + * the index within it. + */ + const_iterator (const SparsityPattern *sparsity_pattern, + const unsigned int row, + const unsigned int index); + + /** + * Prefix increment. + */ + const_iterator& operator++ (); + + /** + * Postfix increment. + */ + const_iterator operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor& operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor* operator-> () const; + + /** + * Comparison. True, if both + * iterators point to the same + * matrix position. + */ + bool operator == (const const_iterator&) const; + + /** + * Inverse of ==. + */ + bool operator != (const 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 + * smaller. + */ + bool operator < (const const_iterator&) const; + + /** + * Exception + */ + DeclException2 (ExcInvalidIndexWithinRow, + int, int, + << "Attempt to access element " << arg2 + << " of row " << arg1 + << " which doesn't have that many elements."); + + private: + /** + * Store an object of the + * accessor class. + */ + Accessor accessor; + }; + + } + + +/** + * This class implements a wrapper class to use the Trilinos distributed + * sparsity pattern class Epetra_FECrsGraph. This class is designed to be + * used for construction of parallel Trilinos matrices. The functionality of + * this class is modeled after the existing sparsity pattern classes, with + * the difference that this class can work fully in parallel according to a + * partitioning of the sparsity pattern rows. + * + * This class has many similarities to the compressed sparsity pattern + * classes of deal.II (i.e., the classes CompressedSparsityPattern, + * CompressedSetSparsityPattern, and CompressedSimpleSparsityPattern), since + * it can dynamically add elements to the pattern without any memory being + * previously reserved for it. However, it also has a method + * SparsityPattern::compress(), that finalizes the pattern and enables its + * use with Trilinos sparse matrices. + * + * @ingroup TrilinosWrappers + * @ingroup Sparsity + * @author Martin Kronbichler, 2008 + */ + class SparsityPattern : public Subscriptor + { + public: + + /** + * Declare a typedef for the + * iterator class. + */ + typedef SparsityPatternIterators::const_iterator const_iterator; + +/** + * @name Constructors and initalization. + */ +//@{ + /** + * Default constructor. Generates an + * empty (zero-size) sparsity + * pattern. + */ + SparsityPattern (); + + /** + * Constructor for a square sparsity + * pattern using an Epetra_Map 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. + */ + SparsityPattern (const Epetra_Map &InputMap, + const unsigned int n_entries_per_row = 1); + + /** + * 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 &InputMap, + 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 + * 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 number of entries + * argument. + */ + SparsityPattern (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const unsigned int n_entries_per_row = 1); + + /** + * 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 &InputRowMap, + const Epetra_Map &InputColMap, + const std::vector &n_entries_per_row); + + /** + * Generate 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. As + * above, 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. + */ + SparsityPattern (const unsigned int m, + const unsigned int n, + const unsigned int n_entries_per_row = 1); + + /** + * Generate 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. + */ + SparsityPattern (const unsigned int m, + const unsigned int n, + const std::vector &n_entries_per_row); + + /** + * Copy constructor. Sets the calling + * sparsity pattern to be the same as + * the input sparsity pattern. + */ + SparsityPattern (const SparsityPattern &SP); + + /** + * Destructor. Made virtual so that + * one can use pointers to this + * class. + */ + virtual ~SparsityPattern (); + + /** + * Reinitialization function for + * generating a square sparsity + * pattern using an Epetra_Map 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. + */ + void + reinit (const Epetra_Map &InputMap, + const unsigned int n_entries_per_row = 1); + + /** + * 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 + reinit (const Epetra_Map &InputMap, + 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. + * + * The number of columns entries per + * row is specified by the argument + * n_entries_per_row. + */ + void + reinit (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const unsigned int n_entries_per_row = 1); + + /** + * 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 &InputRowMap, + const Epetra_Map &InputColMap, + const std::vector &n_entries_per_row); + + /** + * 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. As + * above, 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 unsigned int m, + const unsigned int n, + const unsigned int n_entries_per_row = 1); + + /** + * 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. + */ + void + reinit (const unsigned int m, + const unsigned int n, + 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. + */ + template + void + reinit (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const SparsityType &SP); + + /** + * 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. + */ + template + void + reinit (const Epetra_Map &InputMap, + const SparsityType &SP); + + /** + * Copy function. Sets the calling + * sparsity pattern to be the same as + * the input sparsity pattern. + */ + void + copy_from (const SparsityPattern &SP); + + /** + * 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 &SP); + + /** + * 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. + */ + 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 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 (); +//@} +/** + * @name Information on the sparsity pattern + */ +//@{ + + /** + * 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. + */ + unsigned int max_entries_per_row () const; + + /** + * Return the number of rows in this + * sparsity pattern. + */ + unsigned int n_rows () const; + + /** + * Return the number of columns in + * this sparsity pattern. + */ + unsigned int 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. + * + * 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 + * n=local_size(). + */ + std::pair + local_range () const; + + /** + * Return whether @p index is + * in the local range or not, + * see also local_range(). + */ + bool in_local_range (const unsigned int index) const; + + /** + * Return the number of nonzero + * elements of this sparsity pattern. + */ + unsigned int n_nonzero_elements () const; + + /** + * Number of entries in a + * specific row. + */ + unsigned int row_length (const unsigned int row) 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. + */ + bool empty () const; + + /** + * Return whether the index + * (i,j) exists in the + * sparsity pattern (i.e., it may be + * non-zero) or not. + */ + bool exists (const unsigned int i, + const unsigned int j) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. Currently not + * implemented for this class. + */ + unsigned int memory_consumption () const; + +//@} +/** + * @name Adding entries + */ +//@{ + /** + * Add the element (i,j) to + * the sparsity pattern. + */ + void add (const unsigned int i, + const unsigned int j); + + + /** + * Add several elements in one row to + * the sparsity pattern. + */ + void add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices); +//@} +/** + * @name Iterators + */ +//@{ + + /** + * STL-like iterator with the + * first entry. + */ + const_iterator begin () const; + + /** + * Final iterator. + */ + const_iterator end () const; + + /** + * 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. + */ + const_iterator begin (const unsigned int 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. + * + * 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 unsigned int r) const; + +//@} +/** + * @name Input/Output + */ +//@{ + + /** + * 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 sparsity pattern to the + * given stream, using the format + * (line,col). + */ + void print (std::ostream &out) const; + + // TODO: Write an overloading + // of the operator << for output. + // Since the underlying Trilinos + // object supports it, this should + // be very easy. + +//@} + /** @addtogroup Exceptions + * @{ */ + /** + * Exception + */ + DeclException1 (ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occured while calling a Trilinos function"); + + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The entry with index <" << arg1 << ',' << arg2 + << "> does not exist."); + + /** + * Exception + */ + DeclException0 (ExcSourceEqualsDestination); + + /** + * Exception + */ + DeclException4 (ExcAccessToNonLocalElement, + int, int, int, int, + << "You tried to access element (" << arg1 + << "/" << arg2 << ")" + << " of a distributed matrix, but only rows " + << arg3 << " through " << arg4 + << " are stored locally and can be accessed."); + + /** + * Exception + */ + DeclException2 (ExcAccessToNonPresentElement, + int, int, + << "You tried to access element (" << arg1 + << "/" << arg2 << ")" + << " of a sparse matrix, but it appears to not" + << " exist in the Trilinos sparsity pattern."); + //@} + private: + /** + * Epetra Trilinos + * mapping of the matrix rows + * that assigns parts of the + * matrix to the individual + * processes. This map is + * provided either via the + * constructor or in a reinit + * function. + */ + Epetra_Map row_map; + + /** + * Pointer to the user-supplied + * Epetra Trilinos mapping of + * the matrix columns that + * assigns parts of the matrix + * to the individual processes. + */ + Epetra_Map col_map; + + /** + * 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. + */ + std::auto_ptr graph; + + friend class SparseMatrix; + friend class SparsityPatternIterators::const_iterator; + }; + + + +// -------------------------- inline and template functions ---------------------- + + +#ifndef DOXYGEN + + namespace SparsityPatternIterators + { + + inline + const_iterator::Accessor:: + Accessor (const SparsityPattern *sp, + const unsigned int row, + const unsigned int index) + : + sparsity_pattern(const_cast(sp)), + a_row(row), + a_index(index) + { + visit_present_row (); + } + + + inline + unsigned int + const_iterator::Accessor::row() const + { + Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern()); + return a_row; + } + + + + inline + unsigned int + const_iterator::Accessor::column() const + { + Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern()); + return (*colnum_cache)[a_index]; + } + + + + inline + unsigned int + const_iterator::Accessor::index() const + { + Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern()); + return a_index; + } + + + + inline + const_iterator:: + const_iterator(const SparsityPattern *sp, + const unsigned int row, + const unsigned int index) + : + accessor(sp, row, index) + {} + + + + inline + const_iterator & + const_iterator::operator++ () + { + Assert (accessor.a_row < accessor.sparsity_pattern->n_rows(), + ExcIteratorPastEnd()); + + ++accessor.a_index; + + // If at end of line: do one + // step, then cycle until we + // find a row with a nonzero + // number of entries. + if (accessor.a_index >= accessor.colnum_cache->size()) + { + accessor.a_index = 0; + ++accessor.a_row; + + while ((accessor.a_row < accessor.sparsity_pattern->n_rows()) + && + (accessor.sparsity_pattern->row_length(accessor.a_row) == 0)) + ++accessor.a_row; + + accessor.visit_present_row(); + } + return *this; + } + + + + inline + const_iterator + const_iterator::operator++ (int) + { + const const_iterator old_state = *this; + ++(*this); + return old_state; + } + + + + inline + const const_iterator::Accessor & + const_iterator::operator* () const + { + return accessor; + } + + + + inline + const const_iterator::Accessor * + const_iterator::operator-> () const + { + return &accessor; + } + + + + inline + bool + const_iterator:: + operator == (const const_iterator& other) const + { + return (accessor.a_row == other.accessor.a_row && + accessor.a_index == other.accessor.a_index); + } + + + + inline + bool + const_iterator:: + operator != (const const_iterator& other) const + { + return ! (*this == other); + } + + + + inline + bool + const_iterator:: + operator < (const const_iterator& other) const + { + return (accessor.row() < other.accessor.row() || + (accessor.row() == other.accessor.row() && + accessor.index() < other.accessor.index())); + } + + } + + + + inline + SparsityPattern::const_iterator + SparsityPattern::begin() const + { + return const_iterator(this, 0, 0); + } + + + + inline + SparsityPattern::const_iterator + SparsityPattern::end() const + { + return const_iterator(this, n_rows(), 0); + } + + + + inline + SparsityPattern::const_iterator + SparsityPattern::begin(const unsigned int r) const + { + Assert (r < n_rows(), ExcIndexRange(r, 0, n_rows())); + if (row_length(r) > 0) + return const_iterator(this, r, 0); + else + return end (r); + } + + + + inline + SparsityPattern::const_iterator + SparsityPattern::end(const unsigned int r) const + { + Assert (r < n_rows(), ExcIndexRange(r, 0, n_rows())); + + // place the iterator on the first entry + // past this line, or at the end of the + // matrix + for (unsigned int i=r+1; i 0) + return const_iterator(this, i, 0); + + // if there is no such line, then take the + // end iterator of the matrix + return end(); + } + + + + inline + bool + SparsityPattern::in_local_range (const unsigned int index) const + { + int begin, end; + begin = graph->RowMap().MinMyGID(); + end = graph->RowMap().MaxMyGID()+1; + + return ((index >= static_cast(begin)) && + (index < static_cast(end))); + } + + + + inline + bool + SparsityPattern::is_compressed () const + { + return compressed; + } + + + + inline + bool + SparsityPattern::empty () const + { + return ((n_rows() == 0) && (n_cols() == 0)); + } + + + + inline + void + SparsityPattern::add (const unsigned int i, + const unsigned int j) + { + add (i, 1, &j); + } + + + + inline + void + SparsityPattern::add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices) + { + int * col_index_ptr = (int*)col_indices; + compressed = false; + + int ierr; + + // If the calling sparsity pattern owns + // the row to which we want to add + // values, we can directly call the + // Epetra_CrsGraph input function, which + // is much faster than the + // Epetra_FECrsGraph function. + if (row_map.MyGID(row) == true) + ierr = graph->Epetra_CrsGraph::InsertGlobalIndices(row, + n_cols, + col_index_ptr); + else + { + // When we're at off-processor data, we + // have to stick with the standard + // SumIntoGlobalValues + // 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 a time). + + ierr = graph->InsertGlobalIndices (1, (int*)&row, n_cols, + col_index_ptr); + } + + //Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0])); + AssertThrow (ierr >= 0, ExcTrilinosError(ierr)); + } + + +#endif // DOXYGEN +} + + +DEAL_II_NAMESPACE_CLOSE + + +#endif // DEAL_II_USE_TRILINOS + + +/*-------------------- trilinos_sparsity_pattern.h --------------------*/ + +#endif +/*-------------------- trilinos_sparsity_pattern.h --------------------*/ diff --git a/deal.II/lac/source/block_sparsity_pattern.cc b/deal.II/lac/source/block_sparsity_pattern.cc index 6f6ea0be85..37481370fc 100644 --- a/deal.II/lac/source/block_sparsity_pattern.cc +++ b/deal.II/lac/source/block_sparsity_pattern.cc @@ -309,7 +309,9 @@ template class BlockSparsityPatternBase; template class BlockSparsityPatternBase; template class BlockSparsityPatternBase; template class BlockSparsityPatternBase; - +#ifdef DEAL_II_USE_TRILINOS +template class BlockSparsityPatternBase; +#endif @@ -473,25 +475,27 @@ BlockCompressedSparsityPattern::reinit ( this->collect_sizes(); } -BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern () + + +BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern () {} -BlockCompressedSimpleSparsityPattern:: -BlockCompressedSimpleSparsityPattern (const unsigned int n_rows, +BlockCompressedSetSparsityPattern:: +BlockCompressedSetSparsityPattern (const unsigned int n_rows, const unsigned int n_columns) : - BlockSparsityPatternBase(n_rows, + BlockSparsityPatternBase(n_rows, n_columns) {} -BlockCompressedSimpleSparsityPattern:: -BlockCompressedSimpleSparsityPattern (const std::vector& row_indices, +BlockCompressedSetSparsityPattern:: +BlockCompressedSetSparsityPattern (const std::vector& row_indices, const std::vector& col_indices) : - BlockSparsityPatternBase(row_indices.size(), + BlockSparsityPatternBase(row_indices.size(), col_indices.size()) { for (unsigned int i=0;i& row_indic void -BlockCompressedSimpleSparsityPattern::reinit ( +BlockCompressedSetSparsityPattern::reinit ( const std::vector< unsigned int > &row_block_sizes, const std::vector< unsigned int > &col_block_sizes) { - BlockSparsityPatternBase::reinit(row_block_sizes.size(), col_block_sizes.size()); + BlockSparsityPatternBase::reinit(row_block_sizes.size(), col_block_sizes.size()); for (unsigned int i=0;iblock(i,j).reinit(row_block_sizes[i],col_block_sizes[j]); @@ -515,25 +519,26 @@ BlockCompressedSimpleSparsityPattern::reinit ( -BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern () +BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern () {} -BlockCompressedSetSparsityPattern:: -BlockCompressedSetSparsityPattern (const unsigned int n_rows, +BlockCompressedSimpleSparsityPattern:: +BlockCompressedSimpleSparsityPattern (const unsigned int n_rows, const unsigned int n_columns) : - BlockSparsityPatternBase(n_rows, + BlockSparsityPatternBase(n_rows, n_columns) {} -BlockCompressedSetSparsityPattern:: -BlockCompressedSetSparsityPattern (const std::vector& row_indices, + +BlockCompressedSimpleSparsityPattern:: +BlockCompressedSimpleSparsityPattern (const std::vector& row_indices, const std::vector& col_indices) : - BlockSparsityPatternBase(row_indices.size(), + BlockSparsityPatternBase(row_indices.size(), col_indices.size()) { for (unsigned int i=0;i& row_indices, } + void -BlockCompressedSetSparsityPattern::reinit ( +BlockCompressedSimpleSparsityPattern::reinit ( const std::vector< unsigned int > &row_block_sizes, const std::vector< unsigned int > &col_block_sizes) { - BlockSparsityPatternBase::reinit(row_block_sizes.size(), col_block_sizes.size()); + BlockSparsityPatternBase::reinit(row_block_sizes.size(), col_block_sizes.size()); for (unsigned int i=0;iblock(i,j).reinit(row_block_sizes[i],col_block_sizes[j]); @@ -557,6 +563,81 @@ BlockCompressedSetSparsityPattern::reinit ( +#ifdef DEAL_II_USE_TRILINOS +namespace TrilinosWrappers +{ + + BlockSparsityPattern::BlockSparsityPattern () + {} + + + + BlockSparsityPattern:: + BlockSparsityPattern (const unsigned int n_rows, + const unsigned int n_columns) + : + dealii::BlockSparsityPatternBase(n_rows, + n_columns) + {} + + + + BlockSparsityPattern:: + BlockSparsityPattern (const std::vector& row_indices, + const std::vector& col_indices) + : + BlockSparsityPatternBase(row_indices.size(), + col_indices.size()) + { + for (unsigned int i=0;iblock(i,j).reinit(row_indices[i],col_indices[j]); + this->collect_sizes(); + } + + + + BlockSparsityPattern:: + BlockSparsityPattern (const std::vector& input_maps) + : + BlockSparsityPatternBase(input_maps.size(), + input_maps.size()) + { + for (unsigned int i=0;iblock(i,j).reinit(input_maps[i],input_maps[j]); + this->collect_sizes(); + } + + + + void + BlockSparsityPattern::reinit (const std::vector &row_block_sizes, + const std::vector &col_block_sizes) + { + dealii::BlockSparsityPatternBase::reinit(row_block_sizes.size(), col_block_sizes.size()); + for (unsigned int i=0;iblock(i,j).reinit(row_block_sizes[i],col_block_sizes[j]); + this->collect_sizes(); + } + + + + void + BlockSparsityPattern::reinit (const std::vector &input_maps) + { + dealii::BlockSparsityPatternBase::reinit(input_maps.size(), + input_maps.size()); + for (unsigned int i=0;iblock(i,j).reinit(input_maps[i],input_maps[j]); + this->collect_sizes(); + } + +} + +#endif // Remark: The explicit instantiations for "BlockSparsityPatternBase" were moved // to the top of this source file. The reason is a slightly buggy version diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc index 8fc6101880..eea3d942dc 100644 --- a/deal.II/lac/source/trilinos_block_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_block_sparse_matrix.cc @@ -164,6 +164,31 @@ namespace TrilinosWrappers + template <> + void + BlockSparseMatrix:: + reinit (const BlockSparsityPattern &block_sparsity_pattern) + { + + // Call the other basic reinit function, ... + reinit (block_sparsity_pattern.n_block_rows(), + block_sparsity_pattern.n_block_cols()); + + // ... set the correct sizes, ... + this->row_block_indices = block_sparsity_pattern.get_row_indices(); + this->column_block_indices = block_sparsity_pattern.get_column_indices(); + + // ... and then assign the correct + // data to the blocks. + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + this->sub_objects[r][c]->reinit (block_sparsity_pattern.block(r,c)); + } + } + + + void BlockSparseMatrix:: reinit (const std::vector &input_maps, @@ -369,7 +394,7 @@ namespace TrilinosWrappers // -------------------- explicit instantiations ----------------------- // template void - BlockSparseMatrix::reinit (const BlockSparsityPattern &); + BlockSparseMatrix::reinit (const dealii::BlockSparsityPattern &); template void BlockSparseMatrix::reinit (const BlockCompressedSparsityPattern &); template void @@ -380,7 +405,7 @@ namespace TrilinosWrappers template void BlockSparseMatrix::reinit (const std::vector &, - const BlockSparsityPattern &); + const dealii::BlockSparsityPattern &); template void BlockSparseMatrix::reinit (const std::vector &, const BlockCompressedSparsityPattern &); diff --git a/deal.II/lac/source/trilinos_solver_block.cc b/deal.II/lac/source/trilinos_solver_block.cc index 5021328f0d..4b5fa34598 100644 --- a/deal.II/lac/source/trilinos_solver_block.cc +++ b/deal.II/lac/source/trilinos_solver_block.cc @@ -276,7 +276,7 @@ namespace TrilinosWrappers break; case gmres: aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Aztec Solver", "GMRES"); + .sublist("AztecOO Settings").set("Aztec Solver", "MinRes"); aztecBlockParams->sublist("Forward Solve") .sublist("AztecOO Settings").set("Size of Krylov Subspace", (int)additional_data.gmres_restart_parameter); diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 3d4f2dedaf..31ced248c3 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -13,6 +13,7 @@ #include +#include #include #include #include @@ -97,6 +98,8 @@ namespace TrilinosWrappers matrix->FillComplete(); } + + SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, const unsigned int n_max_entries_per_row) : @@ -109,6 +112,8 @@ namespace TrilinosWrappers int(n_max_entries_per_row), false))) {} + + SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, const std::vector &n_entries_per_row) : @@ -122,6 +127,8 @@ namespace TrilinosWrappers false))) {} + + SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap, const Epetra_Map &InputColMap, const unsigned int n_max_entries_per_row) @@ -135,6 +142,8 @@ namespace TrilinosWrappers int(n_max_entries_per_row), false))) {} + + SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap, const Epetra_Map &InputColMap, const std::vector &n_entries_per_row) @@ -149,6 +158,8 @@ namespace TrilinosWrappers false))) {} + + SparseMatrix::SparseMatrix (const unsigned int m, const unsigned int n, const unsigned int n_max_entries_per_row) @@ -167,6 +178,8 @@ namespace TrilinosWrappers int(n_max_entries_per_row), false))) {} + + SparseMatrix::SparseMatrix (const unsigned int m, const unsigned int n, const std::vector &n_entries_per_row) @@ -186,6 +199,25 @@ namespace TrilinosWrappers false))) {} + + + SparseMatrix::SparseMatrix (const SparsityPattern &InputSP) + : + Subscriptor(), + row_map (InputSP.row_map), + col_map (InputSP.col_map), + last_action (Zero), + compressed (true), + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, *InputSP.graph, false))) + { + Assert(InputSP.graph->Filled() == true, + ExcMessage("The Trilinos sparsity pattern has not been compressed.")); + compress(); + } + + + SparseMatrix::SparseMatrix (const SparseMatrix &InputMatrix) : Subscriptor(), @@ -474,6 +506,25 @@ namespace TrilinosWrappers + void + SparseMatrix::reinit (const SparsityPattern &sparsity_pattern) + { + matrix.reset(); + + row_map = sparsity_pattern.row_map; + col_map = sparsity_pattern.col_map; + + Assert (sparsity_pattern.graph->Filled() == true, + ExcMessage("The Trilinos sparsity pattern has not been compressed")); + + matrix = std::auto_ptr + (new Epetra_FECrsMatrix(Copy, *sparsity_pattern.graph, false)); + + compress(); + } + + + void SparseMatrix::reinit (const SparseMatrix &sparse_matrix) { @@ -748,20 +799,20 @@ namespace TrilinosWrappers // If the data is not on the // present processor, we throw - // an exception. This is on of + // an exception. This is one of // the two tiny differences to // the el(i,j) call, which does // not throw any assertions. - if ((trilinos_i == -1 ) || (trilinos_j == -1)) + if (trilinos_i == -1) { Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first, local_range().second)); } else { - // Check whether the matrix - // already is transformed to - // local indices. + // Check whether the matrix has + // already been transformed to local + // indices. if (matrix->Filled() == false) matrix->GlobalAssemble(col_map, row_map, true); @@ -914,7 +965,7 @@ namespace TrilinosWrappers unsigned int SparseMatrix::n () const { - int n_cols = matrix -> NumGlobalCols(); + unsigned int n_cols = matrix -> NumGlobalCols(); return n_cols; } @@ -923,7 +974,7 @@ namespace TrilinosWrappers unsigned int SparseMatrix::local_size () const { - int n_rows = matrix -> NumMyRows(); + unsigned int n_rows = matrix -> NumMyRows(); return n_rows; } @@ -933,7 +984,7 @@ namespace TrilinosWrappers std::pair SparseMatrix::local_range () const { - int begin, end; + unsigned int begin, end; begin = matrix -> RowMap().MinMyGID(); end = matrix -> RowMap().MaxMyGID()+1; @@ -945,9 +996,9 @@ namespace TrilinosWrappers unsigned int SparseMatrix::n_nonzero_elements () const { - int nnz = matrix->NumGlobalNonzeros(); + unsigned int nnz = matrix->NumGlobalNonzeros(); - return static_cast(nnz); + return nnz; } @@ -960,7 +1011,7 @@ namespace TrilinosWrappers // get a representation of the // present row int ncols = -1; - int local_row = matrix->RowMap().LID(row); + int local_row = matrix->LRID(row); // on the processor who owns this // row, we'll have a non-negative @@ -975,6 +1026,7 @@ namespace TrilinosWrappers } + TrilinosScalar SparseMatrix::l1_norm () const { @@ -1374,7 +1426,7 @@ namespace TrilinosWrappers // explicit instantiations // template void - SparseMatrix::reinit (const SparsityPattern &); + SparseMatrix::reinit (const dealii::SparsityPattern &); template void SparseMatrix::reinit (const CompressedSparsityPattern &); template void @@ -1385,7 +1437,7 @@ namespace TrilinosWrappers template void SparseMatrix::reinit (const Epetra_Map &, - const SparsityPattern &); + const dealii::SparsityPattern &); template void SparseMatrix::reinit (const Epetra_Map &, const CompressedSparsityPattern &); @@ -1400,7 +1452,7 @@ namespace TrilinosWrappers template void SparseMatrix::reinit (const Epetra_Map &, const Epetra_Map &, - const SparsityPattern &); + const dealii::SparsityPattern &); template void SparseMatrix::reinit (const Epetra_Map &, const Epetra_Map &, diff --git a/deal.II/lac/source/trilinos_sparsity_pattern.cc b/deal.II/lac/source/trilinos_sparsity_pattern.cc new file mode 100755 index 0000000000..5de4418625 --- /dev/null +++ b/deal.II/lac/source/trilinos_sparsity_pattern.cc @@ -0,0 +1,735 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include + +#include +#include +#include +#include + +#ifdef DEAL_II_USE_TRILINOS + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + namespace SparsityPatternIterators + { + void + SparsityPattern::const_iterator::Accessor:: + visit_present_row () + { + // if we are asked to visit the + // past-the-end line, then simply + // release all our caches and go on + // with life + if (this->a_row == sparsity_pattern->n_rows()) + { + colnum_cache.reset (); + + return; + } + + // otherwise first flush Trilinos caches + sparsity_pattern->compress (); + + // get a representation of the present + // row + int ncols; + int colnums = sparsity_pattern->n_cols(); + + int ierr; + ierr = sparsity_pattern->graph->ExtractGlobalRowCopy((int)this->a_row, + colnums, + ncols, + (int*)&(*colnum_cache)[0]); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + // copy it into our caches if the + // line isn't empty. if it is, then + // we've done something wrong, since + // we shouldn't have initialized an + // iterator for an empty line (what + // would it point to?) + Assert (ncols != 0, ExcInternalError()); + colnum_cache.reset (new std::vector (colnums, + colnums+ncols)); + } + } + + + // The constructor is actually the + // only point where we have to check + // whether we build a serial or a + // parallel Trilinos matrix. + // Actually, it does not even matter + // how many threads there are, but + // only if we use an MPI compiler or + // a standard compiler. So, even one + // thread on a configuration with + // MPI will still get a parallel + // interface. + SparsityPattern::SparsityPattern () + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)), +#else + row_map (0, 0, Epetra_SerialComm()), +#endif + col_map (row_map), + compressed (true), + graph (std::auto_ptr + (new Epetra_FECrsGraph(View, row_map, 0))) + { + graph->FillComplete(); + } + + SparsityPattern::SparsityPattern (const Epetra_Map &InputMap, + const unsigned int n_entries_per_row) + : + row_map (InputMap), + col_map (row_map), + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + int(n_entries_per_row), false))) + {} + + SparsityPattern::SparsityPattern (const Epetra_Map &InputMap, + const std::vector &n_entries_per_row) + : + row_map (InputMap), + col_map (row_map), + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + (int*)const_cast(&(n_entries_per_row[0])), + false))) + {} + + SparsityPattern::SparsityPattern (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const unsigned int n_entries_per_row) + : + row_map (InputRowMap), + col_map (InputColMap), + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + int(n_entries_per_row), false))) + {} + + SparsityPattern::SparsityPattern (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const std::vector &n_entries_per_row) + : + row_map (InputRowMap), + col_map (InputColMap), + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + (int*)const_cast(&(n_entries_per_row[0])), + false))) + {} + + SparsityPattern::SparsityPattern (const unsigned int m, + const unsigned int n, + const unsigned int n_entries_per_row) + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)), + col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)), +#else + row_map (m, 0, Epetra_SerialComm()), + col_map (n, 0, Epetra_SerialComm()), +#endif + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + int(n_entries_per_row), false))) + {} + + SparsityPattern::SparsityPattern (const unsigned int m, + const unsigned int n, + const std::vector &n_entries_per_row) + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)), + col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)), +#else + row_map (m, 0, Epetra_SerialComm()), + col_map (n, 0, Epetra_SerialComm()), +#endif + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + (int*)const_cast(&(n_entries_per_row[0])), + false))) + {} + + // Copy function is currently not working + // because the Trilinos Epetra_FECrsGraph + // does not implement a reinit function + // from another graph. + /* + SparsityPattern::SparsityPattern (const SparsityPattern &InputSP) + : + Subscriptor(), + row_map (InputSP.row_map), + col_map (InputSP.col_map), + compressed (false), + graph (std::auto_ptr + (new Epetra_FECrsGraph(*InputSP.graph))) + {} + */ + + + + SparsityPattern::~SparsityPattern () + {} + + + + void + SparsityPattern::reinit (const Epetra_Map &input_map, + const unsigned int n_entries_per_row) + { + reinit (input_map, input_map, n_entries_per_row); + } + + + void + SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const unsigned int n_entries_per_row) + { +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD); +#else + Epetra_SerialComm trilinos_communicator; +#endif + + const Epetra_Map rows (m, 0, trilinos_communicator); + const Epetra_Map columns (n, 0, trilinos_communicator); + + reinit (rows, columns, n_entries_per_row); + } + + + void + SparsityPattern::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const unsigned int n_entries_per_row) + { + graph.reset(); + + row_map = input_row_map; + col_map = input_col_map; + + graph = std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, n_entries_per_row, false)); + } + + + + void + SparsityPattern::reinit (const Epetra_Map &input_map, + const std::vector &n_entries_per_row) + { + reinit (input_map, input_map, n_entries_per_row); + } + + + + void + SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const std::vector &n_entries_per_row) + { +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD); +#else + Epetra_SerialComm trilinos_communicator; +#endif + + const Epetra_Map rows (m, 0, trilinos_communicator); + const Epetra_Map columns (n, 0, trilinos_communicator); + + reinit (rows, columns, n_entries_per_row); + } + + + + void + SparsityPattern::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const std::vector &n_entries_per_row) + { + graph.reset(); + + Assert (n_entries_per_row.size() == + static_cast(input_row_map.NumGlobalElements()), + ExcDimensionMismatch (n_entries_per_row.size(), + input_row_map.NumGlobalElements())); + row_map = input_row_map; + col_map = input_col_map; + + graph = std::auto_ptr + (new Epetra_FECrsGraph(Copy, row_map, + n_entries_per_row[input_row_map.MinMyGID()], + false)); + } + + + + template + void + SparsityPattern::reinit (const Epetra_Map &input_map, + const SparsityType &sp) + { + reinit (input_map, input_map, sp); + } + + + + template + void + SparsityPattern::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const SparsityType &sp) + { + graph.reset(); + + Assert (sp.n_rows() == + static_cast(input_row_map.NumGlobalElements()), + ExcDimensionMismatch (sp.n_rows(), + input_row_map.NumGlobalElements())); + Assert (sp.n_cols() == + static_cast(input_col_map.NumGlobalElements()), + ExcDimensionMismatch (sp.n_cols(), + input_col_map.NumGlobalElements())); + + row_map = input_row_map; + col_map = input_col_map; + + const unsigned int n_rows = sp.n_rows(); + + std::vector n_entries_per_row(n_rows); + + for (unsigned int row=0; row + (new Epetra_FECrsGraph(Copy, row_map, + n_entries_per_row[input_row_map.MinMyGID()], + false)); + + Assert (graph->NumGlobalRows() == (int)sp.n_rows(), + ExcDimensionMismatch (graph->NumGlobalRows(), + sp.n_rows())); + + + std::vector row_indices; + + for (unsigned int row=0; rowEpetra_CrsGraph::InsertGlobalIndices (row, row_length, + &row_indices[0]); + } + + compress(); + } + + + + template<> + void + SparsityPattern::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const CompressedSetSparsityPattern &sp) + { + graph.reset(); + + Assert (sp.n_rows() == + static_cast(input_row_map.NumGlobalElements()), + ExcDimensionMismatch (sp.n_rows(), + input_row_map.NumGlobalElements())); + Assert (sp.n_cols() == + static_cast(input_col_map.NumGlobalElements()), + ExcDimensionMismatch (sp.n_cols(), + input_col_map.NumGlobalElements())); + + row_map = input_row_map; + col_map = input_col_map; + + const unsigned int n_rows = sp.n_rows(); + + std::vector n_entries_per_row(n_rows); + + for (unsigned int row=0; row + (new Epetra_FECrsGraph(Copy, row_map, + n_entries_per_row[input_row_map.MinMyGID()], + false)); + + Assert (graph->NumGlobalRows() == (int)sp.n_rows(), + ExcDimensionMismatch (graph->NumGlobalRows(), + sp.n_rows())); + + + std::vector row_indices; + + for (unsigned int row=0; rowEpetra_CrsGraph::InsertGlobalIndices (row, row_length, + &row_indices[0]); + } + + compress(); + } + + + + /* void + SparsityPattern::copy_from (const SparsityPattern &sp) + { + graph.reset(); + row_map = sp.row_map; + col_map = sp.col_map; + + graph = std::auto_ptr (new Epetra_FECrsGraph(*sp.graph)); + } + */ + + + template + void + SparsityPattern::copy_from (const SparsityType &sp) + { +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD); +#else + Epetra_SerialComm trilinos_communicator; +#endif + + const Epetra_Map rows (sp.n_rows(), 0, trilinos_communicator); + const Epetra_Map columns (sp.n_cols(), 0, trilinos_communicator); + + reinit (rows, columns, sp); + } + + + + void + SparsityPattern::clear () + { + // When we clear the matrix, reset + // the pointer and generate an + // empty matrix. + graph.reset(); + +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + row_map = Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)); +#else + row_map = Epetra_Map (0, 0, Epetra_SerialComm()); +#endif + + col_map = row_map; + + graph = std::auto_ptr + (new Epetra_FECrsGraph(View, row_map, 0)); + + graph->FillComplete(); + + compressed = true; + } + + + + void + SparsityPattern::compress () + { + // flush buffers + int ierr; + ierr = graph->GlobalAssemble (col_map, row_map, true); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = graph->OptimizeStorage (); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + compressed = true; + } + + + + bool + SparsityPattern::exists (const unsigned int i, + const unsigned int j) const + { + // Extract local indices in + // the matrix. + int trilinos_i = graph->LRID(i), trilinos_j = graph->LRID(j); + + // If the data is not on the + // present processor, we throw + // an exception. This is on of + // the two tiny differences to + // the el(i,j) call, which does + // not throw any assertions. + if (trilinos_i == -1) + { + return false; + } + else + { + // Check whether the matrix + // already is transformed to + // local indices. + if (graph->Filled() == false) + { + int nnz_present = graph->NumGlobalIndices(i); + int nnz_extracted; + int *col_indices; + + // Generate the view and make + // sure that we have not generated + // an error. + int ierr = graph->ExtractGlobalRowView(trilinos_i, nnz_extracted, + col_indices); + Assert (ierr==0, ExcTrilinosError(ierr)); + Assert (nnz_present == nnz_extracted, + ExcDimensionMismatch(nnz_present, nnz_extracted)); + + // Search the index + int* el_find = std::find(col_indices, col_indices + nnz_present, + trilinos_j); + + int local_col_index = (int)(el_find - col_indices); + + if (local_col_index == nnz_present) + return false; + } + else + { + // Prepare pointers for extraction + // of a view of the row. + int nnz_present = graph->NumGlobalIndices(i); + int nnz_extracted; + int *col_indices; + + // Generate the view and make + // sure that we have not generated + // an error. + int ierr = graph->ExtractMyRowView(trilinos_i, nnz_extracted, + col_indices); + Assert (ierr==0, ExcTrilinosError(ierr)); + + Assert (nnz_present == nnz_extracted, + ExcDimensionMismatch(nnz_present, nnz_extracted)); + + // Search the index + int* el_find = std::find(col_indices, col_indices + nnz_present, + trilinos_j); + + int local_col_index = (int)(el_find - col_indices); + + if (local_col_index == nnz_present) + return false; + } + } + + return true; + } + + + + unsigned int + SparsityPattern::n_rows () const + { + int n_rows = graph -> NumGlobalRows(); + + return n_rows; + } + + + + unsigned int + SparsityPattern::n_cols () const + { + int n_cols = graph -> NumGlobalCols(); + return n_cols; + } + + + + unsigned int + SparsityPattern::local_size () const + { + int n_rows = graph -> NumMyRows(); + + return n_rows; + } + + + + std::pair + SparsityPattern::local_range () const + { + unsigned int begin, end; + begin = graph -> RowMap().MinMyGID(); + end = graph -> RowMap().MaxMyGID()+1; + + return std::make_pair (begin, end); + } + + + + unsigned int + SparsityPattern::n_nonzero_elements () const + { + int nnz = graph->NumGlobalEntries(); + + return static_cast(nnz); + } + + + + unsigned int + SparsityPattern::max_entries_per_row () const + { + int nnz = graph->MaxRowDim(); + + return static_cast(nnz); + } + + + + unsigned int + SparsityPattern::row_length (const unsigned int row) const + { + Assert (row < n_rows(), ExcInternalError()); + + // get a representation of the + // present row + int ncols = -1; + int local_row = graph->LRID(row); + + // on the processor who owns this + // row, we'll have a non-negative + // value. + if (local_row >= 0) + ncols = graph->NumMyIndices (local_row); + + return static_cast(ncols); + } + + + + void + SparsityPattern::write_ascii () + { + Assert (false, ExcNotImplemented()); + } + + + + // As of now, no particularly neat + // ouput is generated in case of + // multiple processors. + void + SparsityPattern::print (std::ostream &out) const + { + int * indices; + int num_entries; + + for (int i=0; iNumMyRows(); ++i) + { + graph->ExtractMyRowView (i, num_entries, indices); + for (int j=0; jGRID(j)] << ") " + << std::endl; + } + + AssertThrow (out, ExcIO()); + } + + + + + // explicit instantiations + // + template void + SparsityPattern::copy_from (const dealii::SparsityPattern &); + template void + SparsityPattern::copy_from (const CompressedSparsityPattern &); + template void + SparsityPattern::copy_from (const CompressedSetSparsityPattern &); + template void + SparsityPattern::copy_from (const CompressedSimpleSparsityPattern &); + + + template void + SparsityPattern::reinit (const Epetra_Map &, + const dealii::SparsityPattern &); + template void + SparsityPattern::reinit (const Epetra_Map &, + const CompressedSparsityPattern &); + template void + SparsityPattern::reinit (const Epetra_Map &, + const CompressedSetSparsityPattern &); + template void + SparsityPattern::reinit (const Epetra_Map &, + const CompressedSimpleSparsityPattern &); + + + template void + SparsityPattern::reinit (const Epetra_Map &, + const Epetra_Map &, + const dealii::SparsityPattern &); + template void + SparsityPattern::reinit (const Epetra_Map &, + const Epetra_Map &, + const CompressedSparsityPattern &); + template void + SparsityPattern::reinit (const Epetra_Map &, + const Epetra_Map &, + const CompressedSimpleSparsityPattern &); + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS -- 2.39.5