]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Restructured the internals of PETScWrappers::Precondition* to allow a PETSc PC object...
authorTimo Heister <timo.heister@gmail.com>
Mon, 17 Jan 2011 18:07:20 +0000 (18:07 +0000)
committerTimo Heister <timo.heister@gmail.com>
Mon, 17 Jan 2011 18:07:20 +0000 (18:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@23205 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/petsc_precondition.h
deal.II/source/lac/petsc_precondition.cc
deal.II/source/lac/petsc_solver.cc

index 1f0fd34a7dd53b848bcbb6359c56642b6a68c258..982a50704422ede1b3c2ddeae0d3146670b7b411 100644 (file)
@@ -40,13 +40,22 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Restructured the internals of PETScWrappers::Precondition* to allow a
+PETSc PC object to exist without a solver. New: use Precondition*::vmult() to
+apply the preconditioner once. Preconditioners now have a default constructor
+and an initialize() function and are no longer initialized in the solver call,
+but in the constructor or initialize().
+<br>
+(Timo Heister, 2011/01/17)
+</li>
+
 <li> Fixed: Boundary conditions in the step-23 tutorial program are now
 applied correctly. Matrix columns get eliminated with the used method
 and introduce some contribution to the right hand side coming from
 inhomogeneous boundary values. The old implementation did not reset the
 matrix columns before applying new boundary values.<br>
 (Martin Stoll, Martin Kronbichler, 2011/01/14)
-</ol>
+</li>
 
 <li> Extended: <code>base/tensor.h</code> has an additional collection of
 contractions between three tensors (<i>ie</i>. <code>contract3</code>).
index 6b0b26c2e24a7120a1fcda61fdea975b043f6e9a..47903a6b8e667c71c9e13a3ca9b434997af81d7c 100644 (file)
@@ -47,31 +47,54 @@ namespace PETScWrappers
  * for parallel jobs, such as for example the ILU preconditioner.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionerBase
   {
     public:
                                        /**
-                                        * Constructor. Take a pointer to the
-                                        * matrix from which the preconditioner
-                                        * shall be constructed.
+                                        * Constructor.
                                         */
-      PreconditionerBase (const MatrixBase &matrix);
+      PreconditionerBase ();
 
                                        /**
                                         * Destructor.
                                         */
       virtual ~PreconditionerBase ();
 
+                                      /**
+                                       * Apply the preconditioner once to the
+                                       * given src vector.
+                                       */
+      void vmult (VectorBase       &dst,
+                 const VectorBase &src) const;
 
+
+                                      /**
+                                       * Gives access to the underlying PETSc
+                                       * object.
+                                       */
+      const PC & get_pc () const;
+      
     protected:
+                                      /**
+                                       * the PETSc preconditioner object
+                                       */
+      PC pc;
+      
                                        /**
                                         * A pointer to the matrix that acts as
                                         * a preconditioner.
                                         */
-      const Mat matrix;
+      Mat matrix;
 
+                                      /**
+                                       * Internal function to create the
+                                       * PETSc preconditioner object. Fails
+                                       * if called twice.
+                                       */
+      void create_pc ();
+      
                                        /**
                                         * Conversion operator to get a
                                         * representation of the matrix that
@@ -81,16 +104,7 @@ namespace PETScWrappers
                                         * the PETSc solvers.
                                         */
       operator Mat () const;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is requested by
-                                        * the derived class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const = 0;
-
+      
                                        /**
                                         * Make the solver class a friend,
                                         * since it needs to call the
@@ -109,7 +123,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionJacobi : public PreconditionerBase
   {
@@ -122,6 +136,14 @@ namespace PETScWrappers
       struct AdditionalData
       {};
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionJacobi ();
+      
+      
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -131,21 +153,26 @@ namespace PETScWrappers
       PreconditionJacobi (const MatrixBase     &matrix,
                           const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                      const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -165,7 +192,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionBlockJacobi : public PreconditionerBase
   {
@@ -177,6 +204,13 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
+      
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionBlockJacobi ();
 
                                        /**
                                         * Constructor. Take the matrix which
@@ -187,21 +221,26 @@ namespace PETScWrappers
       PreconditionBlockJacobi (const MatrixBase     &matrix,
                                const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                      const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -214,7 +253,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionSOR : public PreconditionerBase
   {
@@ -239,6 +278,13 @@ namespace PETScWrappers
           double omega;
       };
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionSOR ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -248,21 +294,26 @@ namespace PETScWrappers
       PreconditionSOR (const MatrixBase     &matrix,
                        const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -275,7 +326,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionSSOR : public PreconditionerBase
   {
@@ -300,6 +351,13 @@ namespace PETScWrappers
           double omega;
       };
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionSSOR ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -309,21 +367,26 @@ namespace PETScWrappers
       PreconditionSSOR (const MatrixBase     &matrix,
                         const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -336,7 +399,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionEisenstat : public PreconditionerBase
   {
@@ -361,6 +424,13 @@ namespace PETScWrappers
           double omega;
       };
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionEisenstat ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -370,21 +440,26 @@ namespace PETScWrappers
       PreconditionEisenstat (const MatrixBase     &matrix,
                              const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -397,7 +472,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionICC : public PreconditionerBase
   {
@@ -422,6 +497,13 @@ namespace PETScWrappers
           unsigned int levels;
       };
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionICC ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -431,21 +513,26 @@ namespace PETScWrappers
       PreconditionICC (const MatrixBase     &matrix,
                        const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -458,7 +545,7 @@ namespace PETScWrappers
  * preconditioner may or may not work.
  *
  * @ingroup PETScWrappers
- * @author Wolfgang Bangerth, 2004
+ * @author Wolfgang Bangerth, Timo Heister, 2004, 2011
  */
   class PreconditionILU : public PreconditionerBase
   {
@@ -483,6 +570,13 @@ namespace PETScWrappers
           unsigned int levels;
       };
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionILU ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -492,21 +586,26 @@ namespace PETScWrappers
       PreconditionILU (const MatrixBase     &matrix,
                        const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -568,6 +667,13 @@ namespace PETScWrappers
          double damping;
       };
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionLU ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -577,21 +683,26 @@ namespace PETScWrappers
       PreconditionLU (const MatrixBase     &matrix,
                      const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 
 
@@ -696,6 +807,13 @@ namespace PETScWrappers
 
 
 
+                                      /**
+                                       * Empty Constructor. You need to call
+                                       * initialize() before using this
+                                       * object.
+                                       */
+      PreconditionBoomerAMG ();
+
                                        /**
                                         * Constructor. Take the matrix which
                                         * is used to form the preconditioner,
@@ -705,21 +823,26 @@ namespace PETScWrappers
       PreconditionBoomerAMG (const MatrixBase     &matrix,
                             const AdditionalData &additional_data = AdditionalData());
 
+                                      /**
+                                       * Initializes the preconditioner
+                                       * object and calculate all data that
+                                       * is necessary for applying it in a
+                                       * solver. This function is
+                                       * automatically called when calling
+                                       * the constructor with the same
+                                       * arguments and is only used if you
+                                       * create the preconditioner without
+                                       * arguments.
+                                       */
+      void initialize (const MatrixBase     &matrix,
+                       const AdditionalData &additional_data = AdditionalData());
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
                                         * particular preconditioner.
                                         */
-      const AdditionalData additional_data;
-
-                                       /**
-                                        * Function that takes a Krylov
-                                        * Subspace Preconditioner context
-                                        * object, and sets the type of
-                                        * preconditioner that is appropriate
-                                        * for the present class.
-                                        */
-      virtual void set_preconditioner_type (PC &pc) const;
+      AdditionalData additional_data;
   };
 }
 
index d1b0e496eaef3d92cad1d95e3325ea70a0cceba8..33d54111e1d0a0d029672eee15634f350fb97429 100644 (file)
@@ -19,6 +19,7 @@
 #  include <base/utilities.h>
 #  include <lac/petsc_matrix_base.h>
 #  include <lac/petsc_vector_base.h>
+#  include <lac/petsc_solver.h>
 #  include <petscconf.h>
 #  include <cmath>
 
@@ -26,18 +27,64 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace PETScWrappers
 {
-  PreconditionerBase::PreconditionerBase (const MatrixBase &matrix)
-                  :
-                  matrix (matrix)
+  PreconditionerBase::PreconditionerBase ()
+                 :
+                 pc(NULL), matrix(NULL)
   {}
 
 
-
   PreconditionerBase::~PreconditionerBase ()
-  {}
+  {
+    if (pc!=NULL)
+      {
+       int ierr = PCDestroy(pc);
+       AssertThrow (ierr == 0, ExcPETScError(ierr));
+      }
+  }
+  
+  
+  void
+  PreconditionerBase::vmult (VectorBase       &dst,
+                            const VectorBase &src) const
+  {
+    AssertThrow (pc != NULL, StandardExceptions::ExcInvalidState ());
+
+    int ierr;
+    ierr = PCApply(pc, src, dst);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+  }  
 
+  
+  void
+  PreconditionerBase::create_pc ()
+  {
+                                    // only allow the creation of the
+                                    // preconditioner once
+    AssertThrow (pc == NULL, StandardExceptions::ExcInvalidState ());
+    
+    MPI_Comm comm;
+    int ierr;
+                                    // this ugly cast is necessay because the
+                                    // type Mat and PETScObject are
+                                    // unrelated.
+    ierr = PetscObjectGetComm(reinterpret_cast<PetscObject>(matrix), &comm);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+    
+    ierr = PCCreate(comm, &pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
 
+    ierr = PCSetOperators(pc , matrix, matrix, SAME_PRECONDITIONER);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+  }
 
+
+  const PC &
+  PreconditionerBase::get_pc () const
+  {
+    return pc;
+  }
+
+  
   PreconditionerBase::operator Mat () const
   {
     return matrix;
@@ -46,52 +93,69 @@ namespace PETScWrappers
 
 /* ----------------- PreconditionJacobi -------------------- */
 
-
+  PreconditionJacobi::PreconditionJacobi ()
+  {}
+  
+    
   PreconditionJacobi::PreconditionJacobi (const MatrixBase     &matrix,
                                           const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionJacobi::set_preconditioner_type (PC &pc) const
+  PreconditionJacobi::initialize (const MatrixBase     &matrix_,
+                                 const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+    
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCJACOBI));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
-/* ----------------- PreconditionJacobi -------------------- */
-
+/* ----------------- PreconditionBlockJacobi -------------------- */
 
+  PreconditionBlockJacobi::PreconditionBlockJacobi ()
+  {}
+  
+    
   PreconditionBlockJacobi::
   PreconditionBlockJacobi (const MatrixBase     &matrix,
                            const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
-
+  {
+    initialize(matrix, additional_data);    
+  }
 
   void
-  PreconditionBlockJacobi::set_preconditioner_type (PC &pc) const
+  PreconditionBlockJacobi::initialize (const MatrixBase     &matrix_,
+                                      const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCBJACOBI));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -104,20 +168,26 @@ namespace PETScWrappers
   {}
 
 
+  PreconditionSOR::PreconditionSOR ()
+  {}
+  
 
   PreconditionSOR::PreconditionSOR (const MatrixBase     &matrix,
                                     const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionSOR::set_preconditioner_type (PC &pc) const
+  PreconditionSOR::initialize (const MatrixBase     &matrix_,
+                              const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCSOR));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -128,6 +198,9 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -139,21 +212,27 @@ namespace PETScWrappers
                   omega (omega)
   {}
 
+  
+  PreconditionSSOR::PreconditionSSOR ()
+  {}
 
 
   PreconditionSSOR::PreconditionSSOR (const MatrixBase     &matrix,
                                       const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionSSOR::set_preconditioner_type (PC &pc) const
+  PreconditionSSOR::initialize (const MatrixBase     &matrix_,
+                               const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCSOR));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -168,6 +247,9 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -180,20 +262,26 @@ namespace PETScWrappers
   {}
 
 
+  PreconditionEisenstat::PreconditionEisenstat ()
+  {}
+  
 
   PreconditionEisenstat::PreconditionEisenstat (const MatrixBase     &matrix,
                                                 const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionEisenstat::set_preconditioner_type (PC &pc) const
+  PreconditionEisenstat::initialize (const MatrixBase     &matrix_,
+                                    const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCEISENSTAT));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -204,6 +292,9 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -217,20 +308,26 @@ namespace PETScWrappers
   {}
 
 
+  PreconditionICC::PreconditionICC ()
+  {}
 
+  
   PreconditionICC::PreconditionICC (const MatrixBase     &matrix,
                                     const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionICC::set_preconditioner_type (PC &pc) const
+  PreconditionICC::initialize (const MatrixBase     &matrix_,
+                              const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCICC));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -245,6 +342,9 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -257,20 +357,26 @@ namespace PETScWrappers
   {}
 
 
+  PreconditionILU::PreconditionILU ()
+  {}
+  
 
   PreconditionILU::PreconditionILU (const MatrixBase     &matrix,
                                     const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionILU::set_preconditioner_type (PC &pc) const
+  PreconditionILU::initialize (const MatrixBase     &matrix_,
+                              const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCILU));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -285,6 +391,9 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -306,20 +415,27 @@ namespace PETScWrappers
   {}
 
 
+  PreconditionBoomerAMG::PreconditionBoomerAMG ()
+  {}
+
+  
   PreconditionBoomerAMG::PreconditionBoomerAMG (const MatrixBase     &matrix,
                                                const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionBoomerAMG::set_preconditioner_type (PC &pc) const
+  PreconditionBoomerAMG::initialize (const MatrixBase     &matrix_,
+                                    const AdditionalData &additional_data_)
   {
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+
 #ifdef PETSC_HAVE_HYPRE
-                                     // set the right type for the
-                                     // preconditioner
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -358,6 +474,10 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
 #else // PETSC_HAVE_HYPRE
     (void)pc;
     Assert (false,
@@ -380,20 +500,26 @@ namespace PETScWrappers
   {}
 
 
+  PreconditionLU::PreconditionLU ()
+  {}
 
+  
   PreconditionLU::PreconditionLU (const MatrixBase     &matrix,
                                  const AdditionalData &additional_data)
-                  :
-                  PreconditionerBase (matrix),
-                  additional_data (additional_data)
-  {}
+  {
+    initialize(matrix, additional_data);    
+  }
 
 
   void
-  PreconditionLU::set_preconditioner_type (PC &pc) const
+  PreconditionLU::initialize (const MatrixBase     &matrix_,
+                             const AdditionalData &additional_data_)
   {
-                                     // set the right type for the
-                                     // preconditioner
+    matrix = static_cast<Mat>(matrix_);
+    additional_data = additional_data_;
+    
+    create_pc();
+
     int ierr;
     ierr = PCSetType (pc, const_cast<char *>(PCLU));
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -427,6 +553,9 @@ namespace PETScWrappers
 
     ierr = PCSetFromOptions (pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = PCSetUp (pc);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
index 8f5af44aa60ddd089fb5c4e83b8701d7a00e4fc3..2c8fbdaef60fca5e0dcbfde3fee21f0b481688cb 100644 (file)
@@ -68,7 +68,6 @@ namespace PETScWrappers
                      const PreconditionerBase &preconditioner)
   {
     int ierr;
-
                                      // first create a solver object if this
                                      // is necessary
     if (solver_data.get() == 0)
@@ -92,10 +91,8 @@ namespace PETScWrappers
                                          // preconditioner
         set_solver_type (solver_data->ksp);
 
-        ierr = KSPGetPC (solver_data->ksp, &solver_data->pc);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        preconditioner.set_preconditioner_type (solver_data->pc);
+       ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc());
+       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
                                          // then a convergence monitor
                                          // function. that function simply

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.