From: kronbichler Date: Wed, 22 Oct 2008 10:15:10 +0000 (+0000) Subject: Trilinos sparse matrices can now also be reinitialized from the compressed sparsity... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cfb5f2984b277b26fa57f76512835e529e8cba21;p=dealii-svn.git Trilinos sparse matrices can now also be reinitialized from the compressed sparsity pattern classes. Change step-31 and step-32 according to this change. git-svn-id: https://svn.dealii.org/trunk@17303 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h index 661f6b6d98..90133e1bdc 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -32,6 +31,11 @@ DEAL_II_NAMESPACE_OPEN + // forward declarations +class BlockSparsityPattern; +class BlockCompressedSparsityPattern; +class BlockCompressedSetSparsityPattern; +class BlockCompressedSimpleSparsityPattern; namespace TrilinosWrappers @@ -181,8 +185,9 @@ namespace TrilinosWrappers * assumes that a quadratic block * matrix is generated. */ + template void reinit (const std::vector &input_maps, - const BlockSparsityPattern &block_sparsity_pattern); + const BlockSparsityType &block_sparsity_pattern); /** * Resize the matrix and initialize it @@ -191,8 +196,9 @@ namespace TrilinosWrappers * result is a block matrix for which * all elements are stored locally. */ - void reinit (const BlockSparsityPattern &block_sparsity_pattern); - + template + void reinit (const BlockSparsityType &block_sparsity_pattern); + /** * This function initializes the * Trilinos matrix using the deal.II @@ -208,7 +214,25 @@ namespace TrilinosWrappers const ::dealii::BlockSparseMatrix &deal_ii_sparse_matrix, 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 whose + * modulus is larger than the + * threshold (so zeros in the + * deal.II matrix can be + * filtered away). Since no + * Epetra_Map is given, all the + * elements will be locally + * stored. + */ + void reinit (const ::dealii::BlockSparseMatrix &deal_ii_sparse_matrix, + const double drop_tolerance=1e-13); + + /** * This function calls the compress() * command of all matrices after * the assembly is diff --git a/deal.II/lac/include/lac/trilinos_precondition_block.h b/deal.II/lac/include/lac/trilinos_precondition_block.h index ae4909388a..7615f9c594 100644 --- a/deal.II/lac/include/lac/trilinos_precondition_block.h +++ b/deal.II/lac/include/lac/trilinos_precondition_block.h @@ -286,8 +286,8 @@ namespace TrilinosWrappers * block, so that the number of * outer iterations is * reduced. See the discussion - * in the @step_22 step-22 - * tutorial program. + * in the @ref step_22 + * "step-22" tutorial program. */ unsigned int Av_do_solve; diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 9ecf706bdf..2bf080ecf0 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -16,7 +16,6 @@ #include #include -#include #include #include #include @@ -41,6 +40,11 @@ DEAL_II_NAMESPACE_OPEN + // forward declarations +class SparsityPattern; +class CompressedSparsityPattern; +class CompressedSetSparsityPattern; +class CompressedSimpleSparsityPattern; namespace TrilinosWrappers { @@ -466,28 +470,33 @@ namespace TrilinosWrappers * Epetra matrix know the * position of nonzero entries * according to the sparsity - * pattern. Note that, when - * using this function, the - * matrix must already be - * initialized with a suitable - * Epetra_Map that describes - * the distribution of the - * matrix among the MPI - * processes. Otherwise, an - * error will be thrown. In a - * parallel run, 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 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. * - * 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 SparsityPattern &sparsity_pattern); + template + void reinit (const SparsityType &sparsity_pattern); /** * This function is initializes @@ -526,8 +535,9 @@ namespace TrilinosWrappers * processors in order to avoid a * dead lock. */ - void reinit (const Epetra_Map &input_map, - const SparsityPattern &sparsity_pattern); + template + void reinit (const Epetra_Map &input_map, + const SparsityType &sparsity_pattern); /** * This function is similar to @@ -545,9 +555,10 @@ namespace TrilinosWrappers * processors in order to avoid a * dead lock. */ - void reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const SparsityPattern &sparsity_pattern); + template + void reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const SparsityType &sparsity_pattern); /** * This function copies the @@ -581,6 +592,34 @@ namespace TrilinosWrappers * processors in order to avoid a * dead lock. */ + void reinit (const ::dealii::SparseMatrix &dealii_sparse_matrix, + 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 + * internally generating one. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. + */ void reinit (const Epetra_Map &input_map, const ::dealii::SparseMatrix &dealii_sparse_matrix, const double drop_tolerance=1e-13); diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc index 0995b26e7b..6c82722086 100644 --- a/deal.II/lac/source/trilinos_block_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_block_sparse_matrix.cc @@ -13,6 +13,7 @@ #include +#include #ifdef DEAL_II_USE_TRILINOS @@ -88,10 +89,11 @@ namespace TrilinosWrappers + template void BlockSparseMatrix:: reinit (const std::vector &input_maps, - const BlockSparsityPattern &block_sparsity_pattern) + const BlockSparsityType &block_sparsity_pattern) { Assert (input_maps.size() == block_sparsity_pattern.n_block_rows(), ExcDimensionMismatch (input_maps.size(), @@ -129,9 +131,10 @@ namespace TrilinosWrappers + template void BlockSparseMatrix:: - reinit (const BlockSparsityPattern &block_sparsity_pattern) + reinit (const BlockSparsityType &block_sparsity_pattern) { Assert (block_sparsity_pattern.n_block_rows() == block_sparsity_pattern.n_block_cols(), @@ -158,7 +161,7 @@ namespace TrilinosWrappers reinit (input_maps, block_sparsity_pattern); } - + void @@ -194,6 +197,39 @@ namespace TrilinosWrappers + void + BlockSparseMatrix:: + reinit (const ::dealii::BlockSparseMatrix &dealii_block_sparse_matrix, + const double drop_tolerance) + { + Assert (dealii_block_sparse_matrix.n_block_rows() == + dealii_block_sparse_matrix.n_block_cols(), + ExcDimensionMismatch (dealii_block_sparse_matrix.n_block_rows(), + dealii_block_sparse_matrix.n_block_cols())); + Assert (dealii_block_sparse_matrix.m() == + dealii_block_sparse_matrix.n(), + ExcDimensionMismatch (dealii_block_sparse_matrix.m(), + dealii_block_sparse_matrix.n())); + + // produce a dummy local map and pass it + // off to the other function +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD); +#else + Epetra_SerialComm trilinos_communicator; +#endif + + std::vector input_maps; + for (unsigned int i=0; i &, + const BlockSparsityPattern &); + template void + BlockSparseMatrix::reinit (const std::vector &, + const BlockCompressedSparsityPattern &); + template void + BlockSparseMatrix::reinit (const std::vector &, + const BlockCompressedSetSparsityPattern &); + template void + BlockSparseMatrix::reinit (const std::vector &, + const BlockCompressedSimpleSparsityPattern &); + } diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 4953d8e3db..36a3844dfd 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -13,6 +13,11 @@ #include +#include +#include +#include +#include + #ifdef DEAL_II_USE_TRILINOS DEAL_II_NAMESPACE_OPEN @@ -172,8 +177,9 @@ namespace TrilinosWrappers + template void - SparseMatrix::reinit (const SparsityPattern &sparsity_pattern) + SparseMatrix::reinit (const SparsityType &sparsity_pattern) { #ifdef DEAL_II_COMPILER_SUPPORTS_MPI Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD); @@ -193,19 +199,21 @@ namespace TrilinosWrappers + template void - SparseMatrix::reinit (const Epetra_Map &input_map, - const SparsityPattern &sparsity_pattern) + SparseMatrix::reinit (const Epetra_Map &input_map, + const SparsityType &sparsity_pattern) { reinit (input_map, input_map, sparsity_pattern); } + template void - SparseMatrix::reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const SparsityPattern &sparsity_pattern) + SparseMatrix::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const SparsityType &sparsity_pattern) { matrix.reset(); @@ -295,6 +303,110 @@ namespace TrilinosWrappers + // The CompressedSetSparsityPattern + // class stores the columns + // differently, so we need to + // explicitly rewrite this function + // in that case. + template<> + void + SparseMatrix::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const CompressedSetSparsityPattern &sparsity_pattern) + { + matrix.reset(); + + unsigned int n_rows = sparsity_pattern.n_rows(); + + if (input_row_map.Comm().MyPID() == 0) + { + Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (input_row_map.NumGlobalElements(), + sparsity_pattern.n_rows())); + Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), + ExcDimensionMismatch (input_col_map.NumGlobalElements(), + sparsity_pattern.n_cols())); + } + + row_map = input_row_map; + col_map = input_col_map; + + std::vector n_entries_per_row(n_rows); + + for (unsigned int row=0; row + (new Epetra_FECrsMatrix(Copy, row_map, + &n_entries_per_row[0], false)); + + + // TODO: As of now, assume that the + // sparsity pattern sits at the all + // processors (completely), and let + // each processor set its rows. Since + // this is wasteful, a better solution + // should be found in the future. + Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (matrix->NumGlobalRows(), + sparsity_pattern.n_rows())); + + // Trilinos has a bug for rectangular + // matrices at this point, so do not + // check for consistent column numbers + // here. + // + // this bug is filed in the Sandia + // bugzilla under #4123 and should be + // fixed for version 9.0 +// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(), +// ExcDimensionMismatch (matrix->NumGlobalCols(), +// sparsity_pattern.n_cols())); + + std::vector values; + std::vector row_indices; + + for (unsigned int row=0; rowInsertGlobalValues (row, row_length, + &values[0], &row_indices[0]); + } + + last_action = Zero; + + // In the end, the matrix needs to + // be compressed in order to be + // really ready. + compress(); + } + + + void SparseMatrix::reinit (const SparseMatrix &sparse_matrix) { @@ -309,6 +421,27 @@ namespace TrilinosWrappers + void + SparseMatrix::reinit (const ::dealii::SparseMatrix &dealii_sparse_matrix, + const double drop_tolerance) + { +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD); +#else + Epetra_SerialComm trilinos_communicator; +#endif + + const Epetra_Map rows (dealii_sparse_matrix.m(), + 0, + trilinos_communicator); + const Epetra_Map columns (dealii_sparse_matrix.n(), + 0, + trilinos_communicator); + reinit (rows, columns, dealii_sparse_matrix, drop_tolerance); + } + + + void SparseMatrix::reinit (const Epetra_Map &input_map, const ::dealii::SparseMatrix &dealii_sparse_matrix, @@ -1104,7 +1237,8 @@ namespace TrilinosWrappers // As of now, no particularly neat // ouput is generated in case of // multiple processors. - void SparseMatrix::print (std::ostream &out) const + void + SparseMatrix::print (std::ostream &out) const { double * values; int * indices; @@ -1121,6 +1255,48 @@ namespace TrilinosWrappers AssertThrow (out, ExcIO()); } + + + + // explicit instantiations + // + template void + SparseMatrix::reinit (const SparsityPattern &); + template void + SparseMatrix::reinit (const CompressedSparsityPattern &); + template void + SparseMatrix::reinit (const CompressedSetSparsityPattern &); + template void + SparseMatrix::reinit (const CompressedSimpleSparsityPattern &); + + + template void + SparseMatrix::reinit (const Epetra_Map &, + const SparsityPattern &); + template void + SparseMatrix::reinit (const Epetra_Map &, + const CompressedSparsityPattern &); + template void + SparseMatrix::reinit (const Epetra_Map &, + const CompressedSetSparsityPattern &); + template void + SparseMatrix::reinit (const Epetra_Map &, + const CompressedSimpleSparsityPattern &); + + + template void + SparseMatrix::reinit (const Epetra_Map &, + const Epetra_Map &, + const SparsityPattern &); + template void + SparseMatrix::reinit (const Epetra_Map &, + const Epetra_Map &, + const CompressedSparsityPattern &); + template void + SparseMatrix::reinit (const Epetra_Map &, + const Epetra_Map &, + const CompressedSimpleSparsityPattern &); + } DEAL_II_NAMESPACE_CLOSE