From: kronbichler Date: Fri, 22 Aug 2008 21:56:05 +0000 (+0000) Subject: Updated Trilinos AMG preconditioner to interact with Trilinos wrappers for matrices... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2c78dfa8bd9a52baa155f1185927697a103ad1e3;p=dealii-svn.git Updated Trilinos AMG preconditioner to interact with Trilinos wrappers for matrices and vectors. git-svn-id: https://svn.dealii.org/trunk@16658 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_precondition_amg.h b/deal.II/lac/include/lac/trilinos_precondition_amg.h index 3ab5cfbc1b..1ecf569047 100755 --- a/deal.II/lac/include/lac/trilinos_precondition_amg.h +++ b/deal.II/lac/include/lac/trilinos_precondition_amg.h @@ -15,8 +15,11 @@ #include +#include #include #include +#include +#include #include @@ -24,7 +27,6 @@ #ifdef DEAL_II_USE_TRILINOS #ifdef DEAL_II_COMPILER_SUPPORTS_MPI # include -# include #else # include #endif @@ -65,6 +67,12 @@ namespace TrilinosWrappers * fact that some of the entries in the preconditioner matrix are zero * and hence can be neglected. * + * The implementation is able to distinguish between matrices from + * elliptic problems and convection dominated problems. We use the standard + * options provided by Trilinos ML for elliptic problems, except that we use a + * Chebyshev smoother instead of a symmetric Gauss-Seidel smoother. + * For most elliptic problems, Chebyshev is better than Gauss-Seidel (SSOR). + * * @author Martin Kronbichler, 2008 */ class PreconditionAMG : public Subscriptor @@ -84,15 +92,59 @@ namespace TrilinosWrappers * Let Trilinos compute a * multilevel hierarchy for the * solution of a linear system - * with the given matrix. + * with the given matrix. The + * function uses the matrix + * format specified in + * TrilinosWrappers::SparseMatrix. + */ + void initialize (const SparseMatrix &matrix, + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool higher_order_elements, + const bool elliptic, + const bool output_details); + + /** + * Let Trilinos compute a + * multilevel hierarchy for the + * solution of a linear system + * with the given matrix. This + * function takes a deal.ii matrix + * and copies the content into a + * Trilinos matrix, so the function + * can be considered rather + * inefficient. */ void initialize (const dealii::SparseMatrix &matrix, const std::vector &null_space, const unsigned int null_space_dimension, const bool higher_order_elements, const bool elliptic, - const bool output_details, - const double drop_tolerance = 1e-13); + const bool output_details); + + /** + * This function can be used + * for a faster recalculation of + * the preconditioner construction + * when the matrix entries + * underlying the preconditioner + * have changed, + * but the matrix sparsity pattern + * has remained the same. What this + * function does is to take the + * already generated coarsening + * structure, compute the AMG + * prolongation and restriction + * according to a smoothed aggregation + * strategy and then builds the whole + * multilevel hiearchy. This function + * can be considerably faster in that + * case, since the coarsening pattern + * is usually the most difficult thing + * to do when setting up the + * AMG ML preconditioner. + */ + void reinit (); /** * Multiply the source vector @@ -101,6 +153,15 @@ namespace TrilinosWrappers * and return it in the * destination vector. */ + void vmult (Vector &dst, + const Vector &src) const; + + /** + * Do the same as before, but + * now use deal.II vectors instead + * of the ones provided in the + * Trilinos wrapper class. + */ void vmult (dealii::Vector &dst, const dealii::Vector &src) const; @@ -122,7 +183,7 @@ namespace TrilinosWrappers * A copy of the deal.II matrix * into Trilinos format. */ - boost::shared_ptr Matrix; + boost::shared_ptr Matrix; }; } diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 54b078c972..f3e432bc2a 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -16,7 +16,8 @@ #include #include -#include +#include +#include #include #include @@ -414,8 +415,8 @@ namespace TrilinosWrappers * not directly available, use one of the * other reinit functions. */ - void reinit (const CompressedSparsityPattern &sparsity_pattern, - const unsigned int n_max_entries_per_row); + void reinit (const SparsityPattern &sparsity_pattern, + const unsigned int n_max_entries_per_row); /** * This function initializes the @@ -426,7 +427,22 @@ namespace TrilinosWrappers * of nonzeros from the sparsity * pattern internally. */ - void reinit (const CompressedSparsityPattern &sparsity_pattern); + void reinit (const SparsityPattern &sparsity_pattern); + + /** + * 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). + */ + void reinit (const Epetra_Map &input_map, + const ::dealii::SparseMatrix &deal_ii_sparse_matrix, + const double drop_tolerance=1e-15); /** * This function is similar to the @@ -437,8 +453,8 @@ namespace TrilinosWrappers * when the matrix structure changes, * e.g. when the grid is refined. */ - void reinit (const Epetra_Map &input_map, - const CompressedSparsityPattern &sparsity_pattern); + void reinit (const Epetra_Map &input_map, + const SparsityPattern &sparsity_pattern); /** * This function is similar to the @@ -449,9 +465,9 @@ namespace TrilinosWrappers * To be used e.g. for rectangular * matrices after remeshing. */ - void reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const CompressedSparsityPattern &sparsity_pattern); + void reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const SparsityPattern &sparsity_pattern); /** * Release all memory and return diff --git a/deal.II/lac/source/trilinos_precondition_amg.cc b/deal.II/lac/source/trilinos_precondition_amg.cc index 332241451f..ca338ed565 100755 --- a/deal.II/lac/source/trilinos_precondition_amg.cc +++ b/deal.II/lac/source/trilinos_precondition_amg.cc @@ -13,14 +13,10 @@ #include -#include -#include - #ifdef DEAL_II_USE_TRILINOS #include -#include #include #include #include @@ -43,80 +39,30 @@ namespace TrilinosWrappers PreconditionAMG::~PreconditionAMG () {} - - + + void PreconditionAMG:: - initialize (const dealii::SparseMatrix &matrix, - const std::vector &null_space, - const unsigned int null_space_dimension, - const bool elliptic, - const bool higher_order_elements, - const bool output_details, - const double drop_tolerance) + initialize (const SparseMatrix &matrix, + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool elliptic, + const bool higher_order_elements, + const bool output_details) { - Assert (drop_tolerance >= 0, - ExcMessage ("Drop tolerance must be a non-negative number.")); - const unsigned int n_rows = matrix.m(); - - // Init Epetra Matrix, avoid - // storing the nonzero elements. - { - Map.reset (new Epetra_Map(n_rows, 0, communicator)); - std::vector row_lengths (n_rows); - for (dealii::SparseMatrix::const_iterator p = matrix.begin(); - p != matrix.end(); ++p) - if (std::abs(p->value()) > drop_tolerance) - ++row_lengths[p->row()]; - - Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true)); - - const unsigned int max_nonzero_entries - = *std::max_element (row_lengths.begin(), row_lengths.end()); - - std::vector values(max_nonzero_entries, 0); - std::vector row_indices(max_nonzero_entries); - - for (unsigned int row=0; row::const_iterator p = matrix.begin(row); - p != matrix.end(row); ++p) - if (std::abs(p->value()) > drop_tolerance) - { - row_indices[index] = p->column(); - values[index] = p->value(); - ++index; - } - - Assert (index == static_cast(row_lengths[row]), - ExcMessage("Filtering out zeros could not " - "be successfully finished!")); - - Matrix->InsertGlobalValues(row, row_lengths[row], - &values[0], &row_indices[0]); - } - - Matrix->FillComplete(); - } - + if (!Matrix) + { + Matrix = boost::shared_ptr + (const_cast(&matrix)); + Map = boost::shared_ptr + (const_cast(&matrix.row_map)); + } + // Build the AMG preconditioner. Teuchos::ParameterList parameter_list; - // The implementation is able - // to distinguish between - // matrices from elliptic problems - // and convection dominated - // problems. We use the standard - // options for elliptic problems, - // except that we use a - // Chebyshev smoother instead - // of a symmetric Gauss-Seidel - // smoother. For most elliptic - // problems, Chebyshev is better - // than Gauss-Seidel (SSOR). if (elliptic) { ML_Epetra::SetDefaults("SA",parameter_list); @@ -148,20 +94,67 @@ namespace TrilinosWrappers parameter_list.set("null space: dimension", int(null_space_dimension)); parameter_list.set("null space: vectors", (double *)&null_space[0]); } - + multigrid_operator = boost::shared_ptr - (new ML_Epetra::MultiLevelPreconditioner(*Matrix, parameter_list, true)); + (new ML_Epetra::MultiLevelPreconditioner( + *matrix.matrix, parameter_list, true)); if (output_details) multigrid_operator->PrintUnused(0); } + + + void + PreconditionAMG:: + initialize (const dealii::SparseMatrix &deal_ii_sparse_matrix, + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool elliptic, + const bool higher_order_elements, + const bool output_details) + { + const unsigned int n_rows = deal_ii_sparse_matrix.m(); + + // Init Epetra Matrix, avoid + // storing the nonzero elements. + + Map.reset (new Epetra_Map(n_rows, 0, communicator)); + + Matrix.reset(); + Matrix = boost::shared_ptr (new SparseMatrix()); + + Matrix->reinit (*Map, deal_ii_sparse_matrix); + Matrix->compress(); + + initialize (*Matrix, null_space, null_space_dimension, elliptic, + higher_order_elements, output_details); + } + + + void PreconditionAMG:: + reinit () + { + multigrid_operator->ReComputePreconditioner(); + } + + + + void PreconditionAMG::vmult (Vector &dst, + const Vector &src) const + { + const int ierr = multigrid_operator->ApplyInverse (*src.vector, *dst.vector); + + Assert (ierr == 0, ExcTrilinosError(ierr)); + } + + // For the implementation of the // vmult function we // note that invoking a call of // the Trilinos preconditioner // requires us to use Epetra vectors - // as well. Luckily, it is sufficient + // as well. It is faster // to provide a view, i.e., feed // Trilinos with a pointer to the // data, so we avoid copying the @@ -176,6 +169,15 @@ namespace TrilinosWrappers void PreconditionAMG::vmult (dealii::Vector &dst, const dealii::Vector &src) const { + Assert (Map->SameAs(Matrix->row_map), + ExcMessage("The sparse matrix given to the preconditioner uses " + "a map that is not compatible. Check ML preconditioner " + "setup.")); + Assert (Map->SameAs(Matrix->col_map), + ExcMessage("The sparse matrix given to the preconditioner uses " + "a map that is not compatible. Check ML preconditioner " + "setup.")); + Epetra_Vector LHS (View, *Map, dst.begin()); Epetra_Vector RHS (View, *Map, const_cast(src.begin())); diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 0e9789f71c..e5b4862043 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -147,8 +147,8 @@ namespace TrilinosWrappers void - SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern, - const unsigned int n_max_entries_per_row) + SparseMatrix::reinit (const SparsityPattern &sparsity_pattern, + const unsigned int n_max_entries_per_row) { unsigned int n_rows = sparsity_pattern.n_rows(); @@ -187,14 +187,14 @@ namespace TrilinosWrappers void - SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern) + SparseMatrix::reinit (const SparsityPattern &sparsity_pattern) { unsigned int n_rows = sparsity_pattern.n_rows(); - + Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), ExcDimensionMismatch (matrix->NumGlobalRows(), sparsity_pattern.n_rows())); - + // Trilinos seems to have a bug for // rectangular matrices at this point, // so do not check for consistent @@ -204,10 +204,10 @@ namespace TrilinosWrappers // sparsity_pattern.n_cols())); std::vector n_entries_per_row(n_rows); - + for (unsigned int row=0; row n_entries_per_row(n_rows); - + for (unsigned int row=0; row (new Epetra_FECrsMatrix ( Copy, row_map, &n_entries_per_row[0], true)); - + const unsigned int n_max_entries_per_row = *std::max_element ( &n_entries_per_row[0], &n_entries_per_row[n_rows-1]); @@ -250,32 +251,33 @@ namespace TrilinosWrappers void - SparseMatrix::reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const CompressedSparsityPattern &sparsity_pattern) + SparseMatrix::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const SparsityPattern &sparsity_pattern) { + matrix.reset(); unsigned int n_rows = sparsity_pattern.n_rows(); - matrix.reset(); - row_map = input_row_map; - col_map = input_col_map; - + 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, col_map, &n_entries_per_row[0], true)); - + const unsigned int n_max_entries_per_row = *std::max_element ( &n_entries_per_row[0], &n_entries_per_row[n_rows-1]); @@ -284,6 +286,63 @@ namespace TrilinosWrappers + void + SparseMatrix::reinit (const Epetra_Map &input_map, + const ::dealii::SparseMatrix &deal_ii_sparse_matrix, + const double drop_tolerance) + { + matrix.reset(); + + unsigned int n_rows = deal_ii_sparse_matrix.m(); + + Assert (input_map.NumGlobalElements() == (int)n_rows, + ExcDimensionMismatch (input_map.NumGlobalElements(), + n_rows)); + Assert (input_map.NumGlobalElements() == (int)deal_ii_sparse_matrix.n(), + ExcDimensionMismatch (input_map.NumGlobalElements(), + deal_ii_sparse_matrix.n())); + + row_map = input_map; + col_map = row_map; + + std::vector n_entries_per_row(n_rows); + + for (unsigned int row=0; row + (new Epetra_FECrsMatrix(Copy, row_map, col_map, + &n_entries_per_row[0], true)); + + std::vector values; + std::vector row_indices; + + for (unsigned int row=0; row::const_iterator + p = deal_ii_sparse_matrix.begin(row); + p != deal_ii_sparse_matrix.end(row); ++p) + if (std::abs(p->value()) > drop_tolerance) + { + row_indices[index] = p->column(); + values[index] = p->value(); + ++index; + } + const int n_row_entries = index; + + matrix->InsertGlobalValues(row, n_row_entries, + &values[0], &row_indices[0]); + } + + } + + + void SparseMatrix::clear () {