From 1e7bc6cad7b819a20d12f663f9aab7e36a881ea1 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 24 Oct 2012 02:44:37 +0000 Subject: [PATCH] Enable symmetry mode and pick data from the PETSc environment. git-svn-id: https://svn.dealii.org/trunk@27190 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 +++ deal.II/include/deal.II/lac/petsc_solver.h | 61 +++++++++++++++++----- deal.II/source/lac/petsc_solver.cc | 37 ++++++++++--- 3 files changed, 82 insertions(+), 22 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 49469fe476..02a25131b9 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -127,6 +127,12 @@ DoFHandler, in particular removal of specializations.

Specific improvements

    +
  1. New: The PETScWrappers::SparseDirectMUMPS class now allows to +exploit symmetry of the matrix, using the +PETScWrappers::SparseDirectMUMPS::set_symmetric_mode() function. +
    +(Alexander Grayver, 2012/10/23) +
  2. Fixed: Several static const member variables of the Accessor classes were not properly instantiated. This only rarely created trouble because they are typically only used as template arguments diff --git a/deal.II/include/deal.II/lac/petsc_solver.h b/deal.II/include/deal.II/lac/petsc_solver.h index 984ea7b6b5..4636db9c08 100644 --- a/deal.II/include/deal.II/lac/petsc_solver.h +++ b/deal.II/include/deal.II/lac/petsc_solver.h @@ -197,7 +197,6 @@ namespace PETScWrappers */ virtual void set_solver_type (KSP &ksp) const = 0; - private: /** * Solver prefix name to qualify options * specific to the PETSc KSP object in the @@ -209,6 +208,7 @@ namespace PETScWrappers */ std::string prefix_name; + private: /** * A function that is used in PETSc as * a callback to check on @@ -1143,40 +1143,64 @@ namespace PETScWrappers }; /** - * An implementation of the solver interface using the MUMPS lu solver - * through PETSc. + * An implementation of the solver interface using the sparse direct MUMPS + * solver through PETSc. This class has the usual interface of all other + * solver classes but it is of course different in that it doesn't implement + * an iterative solver. As a consequence, things like the SolverControl object + * have no particular meaning here. + * + * MUMPS allows to make use of symmetry in this matrix. In this class this is + * made possible by the set_symmetric_mode() function. If your matrix is + * symmetric, you can use this class as follows: + * @code + * SolverControl cn; + * PETScWrappers::SparseDirectMUMPS solver(cn, mpi_communicator); + * solver.set_symmetric_mode(true); + * solver.solve(system_matrix, solutions, system_rhs); + * @endcode + * + * @note The class internally calls KSPSetFromOptions thus you are + * able to use all the PETSc parameters for MATSOLVERMUMPS package. + * See http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATSOLVERMUMPS.html * * @ingroup PETScWrappers - * @author Daniel Brauss, 2012 + * @author Daniel Brauss, Alexander Grayver, 2012 */ class SparseDirectMUMPS : public SolverBase { public: /** * Standardized data structure - * to pipe additional data to + * to pipe additional data to * the solver. */ struct AdditionalData {}; /** - * constructor + * Constructor */ SparseDirectMUMPS (SolverControl &cn, const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); /** - * the method to solve the - * linear system. SolverBase has a method - * with this name already. However, the - * different function calls and objects used here - * were not easily incorporated + * The method to solve the + * linear system. */ void solve (const MatrixBase &A, VectorBase &x, const VectorBase &b); + /** + * The method allows to take advantage + * if the system matrix is symmetric by + * using LDL^T decomposition unstead of + * more expensive LU. The argument + * indicates whether the matrix is + * symmetric or not. + */ + void set_symmetric_mode (const bool flag); + protected: /** * Store a copy of flags for this @@ -1188,7 +1212,7 @@ namespace PETScWrappers private: /** - * A function that is used in PETSc + * A function that is used in PETSc * as a callback to check convergence. * It takes the information provided * from PETSc and checks it against @@ -1211,12 +1235,12 @@ namespace PETScWrappers KSPConvergedReason *reason, void *solver_control); /** - * A structure that contains the + * A structure that contains the * PETSc solver and preconditioner * objects. Since the solve member * function in the base is not used * here, the private SolverData struct - * located in the base could not be used + * located in the base could not be used * either */ struct SolverDataMUMPS @@ -1226,6 +1250,15 @@ namespace PETScWrappers }; std_cxx1x::shared_ptr solver_data; + + /** + * Flag specifies whether matrix + * being factorized is symmetric + * or not. It influences the type + * of the used preconditioner + * (PCLU or PCCHOLESKY) + */ + bool symmetric_mode; }; } diff --git a/deal.II/source/lac/petsc_solver.cc b/deal.II/source/lac/petsc_solver.cc index d01dae47f0..b0713c88b6 100644 --- a/deal.II/source/lac/petsc_solver.cc +++ b/deal.II/source/lac/petsc_solver.cc @@ -599,7 +599,8 @@ namespace PETScWrappers const AdditionalData &data) : SolverBase (cn, mpi_communicator), - additional_data (data) + additional_data (data), + symmetric_mode(false) {} @@ -723,13 +724,14 @@ namespace PETScWrappers /** * build PETSc PC for particular - * PCLU preconditioner - * (note: if the matrix is - * symmetric, then a cholesky - * decomposition PCCHOLESKY - * could be used) - */ - ierr = PCSetType (solver_data->pc, PCLU); + * PCLU or PCCHOLESKY preconditioner + * depending on whether the + * symmetric mode has been set + */ + if (symmetric_mode) + ierr = PCSetType (solver_data->pc, PCCHOLESKY); + else + ierr = PCSetType (solver_data->pc, PCLU); AssertThrow (ierr == 0, ExcPETScError(ierr)); /** @@ -771,6 +773,19 @@ namespace PETScWrappers * to MUMPS */ ierr = MatMumpsSetIcntl (F, icntl, ival); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + /** + * set the command line option prefix name + */ + ierr = KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str()); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + /** + * set the command line options provided + * by the user to override the defaults + */ + ierr = KSPSetFromOptions (solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -862,6 +877,12 @@ namespace PETScWrappers return 0; } + void + SparseDirectMUMPS::set_symmetric_mode(const bool flag) + { + symmetric_mode = flag; + } + } DEAL_II_NAMESPACE_CLOSE -- 2.39.5