From 019e0dafac15584270030d35f9fe0213dcf79471 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 6 Oct 2008 17:02:37 +0000 Subject: [PATCH] Added a new function for data exchange to Trilinos vector. Added a block direct solve as Trilinos preconditioner. git-svn-id: https://svn.dealii.org/trunk@17113 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/trilinos_block_sparse_matrix.h | 2 +- .../lac/include/lac/trilinos_block_vector.h | 48 +++++++++++++- .../lac/include/lac/trilinos_precondition.h | 66 +++++++++++++++++-- deal.II/lac/source/trilinos_block_vector.cc | 24 +++++++ deal.II/lac/source/trilinos_precondition.cc | 41 ++++++++++++ deal.II/lac/source/trilinos_vector.cc | 27 ++++++++ 6 files changed, 200 insertions(+), 8 deletions(-) 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 4b9d3a761a..afb5e31cc3 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -116,7 +116,7 @@ namespace TrilinosWrappers */ ~BlockSparseMatrix (); - /** + /** * Pseudo copy operator only copying * empty objects. The sizes of the block * matrices need to be the same. diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index 30ca2aec75..daee99f3a3 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -39,6 +39,7 @@ namespace TrilinosWrappers class BlockVector; } class BlockVector; + class BlockSparseMatrix; namespace MPI @@ -236,7 +237,50 @@ namespace TrilinosWrappers */ void reinit (const unsigned int num_blocks); - /** + /** + * This reinit function is + * meant to be used for + * parallel calculations where + * some non-local data has to + * be used. The typical + * situation where one needs + * this function is the call of + * the + * FEValues::get_function_values + * function (or of some + * derivatives) in + * parallel. Since it is + * usually faster to retrieve + * the data in advance, this + * function can be called + * before the assembly forks + * out to the different + * processors. What this + * function does is the + * following: It takes the + * information in the columns + * of the given matrix and + * looks which data couples + * between the different + * processors. That data is + * then queried from the input + * vector. Note that you should + * not write to the resulting + * vector any more, since the + * some data can be stored + * several times on different + * processors, leading to + * unpredictable results. In + * particular, such a vector + * cannot be used for + * matrix-vector products as + * for example done during the + * solution of linear systems. + */ + void do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m, + const BlockVector &v); + + /** * Compresses all the components * after assembling together all * elements. @@ -245,7 +289,7 @@ namespace TrilinosWrappers /** * Returns the state of the - * matrix, i.e., whether + * vector, i.e., whether * compress() needs to be * called after an operation * requiring data diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index ac8def9bca..2093233864 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -88,7 +88,7 @@ namespace TrilinosWrappers friend class SolverBase; friend class SolverSaddlePoint; - protected: + //protected: /** * This is a pointer to the * preconditioner object. @@ -407,8 +407,9 @@ namespace TrilinosWrappers * AdditionalData structure for details), and a parameter * overlap that determines if and how much overlap there * should be between the matrix partitions on the various MPI - * processes. The default settings are 1 for the relaxation parameter, - * 0 for the diagonal augmentation and 0 for the overlap. + * processes. The default settings are 0 for the additional fill-in, 0 + * for the absolute augmentation tolerance, 1 for the relative + * augmentation tolerance, 0 for the overlap. * * Note that a parallel applicatoin of the IC preconditioner is * actually a block-Jacobi preconditioner with block size equal to the @@ -533,8 +534,9 @@ namespace TrilinosWrappers * AdditionalData structure for details), and a parameter * overlap that determines if and how much overlap there * should be between the matrix partitions on the various MPI - * processes. The default settings are 1 for the relaxation parameter, - * 0 for the diagonal augmentation and 0 for the overlap. + * processes. The default settings are 0 for the additional fill-in, 0 + * for the absolute augmentation tolerance, 1 for the relative + * augmentation tolerance, 0 for the overlap. * * Note that a parallel applicatoin of the ILU preconditioner is * actually a block-Jacobi preconditioner with block size equal to the @@ -633,6 +635,60 @@ namespace TrilinosWrappers void initialize (const SparseMatrix &matrix, const AdditionalData &additional_data = AdditionalData()); }; + + + +/** + * A wrapper class for a sparse direct LU decomposition on parallel + * blocks for Trilinos matrices. When run in serial, this corresponds + * to a direct solve on the matrix. + * + * The AdditionalData data structure allows to set preconditioner + * options. + * + * Note that a parallel applicatoin of the block direct solve + * preconditioner is actually a block-Jacobi preconditioner with block + * size equal to the local matrix size. Spoken more technically, this + * parallel operation is an additive + * Schwarz method with an exact solve as inner solver, + * based on the (outer) parallel partitioning. + * + * @ingroup TrilinosWrappers + * @ingroup Preconditioners + * @author Martin Kronbichler, 2008 + */ + class PreconditionBlockDirect : public PreconditionBase + { + public: + /** + * Standardized data struct to + * pipe additional parameters + * to the preconditioner. + */ + struct AdditionalData + { + /** + * Constructor. + */ + AdditionalData (const unsigned int overlap = 0); + + /** + * Block direct parameters and overlap. + */ + unsigned int overlap; + }; + + /** + * Initialize function. Takes + * the matrix which is used to + * form the preconditioner, and + * additional flags if there + * are any. + */ + void initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data = AdditionalData()); + }; } diff --git a/deal.II/lac/source/trilinos_block_vector.cc b/deal.II/lac/source/trilinos_block_vector.cc index 27f52188f4..acaa66ffbf 100644 --- a/deal.II/lac/source/trilinos_block_vector.cc +++ b/deal.II/lac/source/trilinos_block_vector.cc @@ -12,6 +12,7 @@ //--------------------------------------------------------------------------- #include +#include #ifdef DEAL_II_USE_TRILINOS @@ -135,6 +136,29 @@ namespace TrilinosWrappers + void + BlockVector::do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m, + const BlockVector &v) + { + Assert (m.n_block_rows() == v.n_blocks(), + ExcDimensionMismatch(m.n_block_rows(),v.n_blocks())); + Assert (m.n_block_cols() == v.n_blocks(), + ExcDimensionMismatch(m.n_block_cols(),v.n_blocks())); + + if (v.n_blocks() != n_blocks()) + { + block_indices = v.get_block_indices(); + components.resize(v.n_blocks()); + } + + for (unsigned int i=0; in_blocks(); ++i) + components[i].do_data_exchange(m.block(i,i), v.block(i)); + + collect_sizes(); + } + + + void BlockVector::compress () { diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 55bc5c1bcd..9b4c4ae03c 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -72,6 +72,7 @@ namespace TrilinosWrappers int ierr; Teuchos::ParameterList parameter_list; + parameter_list.set ("relaxation: sweeps", 1); parameter_list.set ("relaxation: type", "Jacobi"); parameter_list.set ("relaxation: damping factor", additional_data.omega); parameter_list.set ("relaxation: min diagonal value", @@ -118,6 +119,7 @@ namespace TrilinosWrappers int ierr; Teuchos::ParameterList parameter_list; + parameter_list.set ("relaxation: sweeps", 1); parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel"); parameter_list.set ("relaxation: damping factor", additional_data.omega); parameter_list.set ("relaxation: min diagonal value", @@ -165,6 +167,7 @@ namespace TrilinosWrappers int ierr; Teuchos::ParameterList parameter_list; + parameter_list.set ("relaxation: sweeps", 1); parameter_list.set ("relaxation: type", "Gauss-Seidel"); parameter_list.set ("relaxation: damping factor", additional_data.omega); parameter_list.set ("relaxation: min diagonal value", @@ -275,6 +278,44 @@ namespace TrilinosWrappers AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } + + +/* ---------------------- PreconditionBlockDirect --------------------- */ + + PreconditionBlockDirect::AdditionalData:: + AdditionalData (const unsigned int overlap) + : + overlap (overlap) + {} + + + + void + PreconditionBlockDirect::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) + { + preconditioner.release(); + + preconditioner = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix, + additional_data.overlap)); + Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); + + int ierr; + + Teuchos::ParameterList parameter_list; + parameter_list.set ("schwarz: combine mode", "Add"); + + ierr = preconditioner->SetParameters(parameter_list); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Initialize(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Compute(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 3501c88181..5e4ad408fe 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -215,6 +215,33 @@ namespace TrilinosWrappers return *this; } + + + void + Vector::do_data_exchange (const TrilinosWrappers::SparseMatrix &m, + const Vector &v) + { + Assert (m.matrix->Filled() == true, + ExcMessage ("Matrix is not compressed. " + "Cannot find exchange information!")); + Assert (v.vector->Map().UniqueGIDs() == true, + ExcMessage ("The input vector has overlapping data, " + "which is not allowed.")); + + if (vector->Map().SameAs(m.matrix->ColMap()) == false) + { + map = m.matrix->ColMap(); + vector = std::auto_ptr (new Epetra_FEVector(map)); + } + + Epetra_Import data_exchange (vector->Map(), v.vector->Map()); + const int ierr = vector->Import(*v.vector, data_exchange, Insert); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + last_action = Insert; + } + } /* end of namespace MPI */ -- 2.39.5