Enable symmetry mode and pick data from the PETSc environment.

git-svn-id: https://svn.dealii.org/trunk@27190 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 49469fe47643743ac79e1bce34c106c0c374c874..02a25131b95a2d189a14fdfc2beb9bedeaf8702a 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -127,6 +127,12 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: The PETScWrappers::SparseDirectMUMPS class now allows to
+exploit symmetry of the matrix, using the
+PETScWrappers::SparseDirectMUMPS::set_symmetric_mode() function.
+<br>
+(Alexander Grayver, 2012/10/23)
+
 <li> 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 984ea7b6b5ebd18c7e2818a4bf38671d9144b326..4636db9c08538c6876dd5d483739bf3b65c531a6 100644 (file)
--- 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<SolverDataMUMPS> 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 d01dae47f093932a1cea97103968cb4552016564..b0713c88b676b8ef81a0ad6ef981396a86240da4 100644 (file)
--- 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