]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implemented a few more wrappers to Trilinos Ifpack preconditioners (Jacobi, SOR,...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Sep 2008 11:45:20 +0000 (11:45 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Sep 2008 11:45:20 +0000 (11:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@16855 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-31/step-31.cc
deal.II/examples/step-32/step-32.cc
deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/source/trilinos_precondition.cc

index 818bf64b99e1b1743597785b4e6bef40e5847bb2..db3be88046d668337c90c5de5d2de93b14d061c4 100644 (file)
@@ -1144,8 +1144,8 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
                                   // keep the (1,1) block, though
       
   Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionSSOR>
-       (new TrilinosWrappers::PreconditionSSOR(
-          stokes_preconditioner_matrix.block(1,1),1.2));
+                                   (new TrilinosWrappers::PreconditionSSOR());
+  Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1,1),1.2);
 
   std::cout << std::endl;
 
@@ -1841,11 +1841,11 @@ void BoussinesqFlowProblem<dim>::solve ()
                                  1e-8*temperature_rhs.l2_norm());
     SolverCG<TrilinosWrappers::Vector>   cg (solver_control);
 
-    TrilinosWrappers::PreconditionSSOR preconditioner (temperature_matrix,
-                                                      1.2);
+    TrilinosWrappers::PreconditionSSOR preconditioner;
+    preconditioner.initialize (temperature_matrix, 1.2);
+
     cg.solve (temperature_matrix, temperature_solution,
-             temperature_rhs,
-             preconditioner);
+             temperature_rhs, preconditioner);
 
                                     // produce a consistent temperature field
     temperature_constraints.distribute (temperature_solution);
index b72f043873e6ebcc00e44b8446979fe0b26a6bf3..ec451be4a2b7ab615cb0b74fd5e1a2ac2322784e 100644 (file)
@@ -800,8 +800,9 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
                                 true, true, 5e-2, null_space, false);
 
   Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionSSOR>
-       (new TrilinosWrappers::PreconditionSSOR(
-          stokes_preconditioner_matrix.block(1,1),1.2));
+                                   (new TrilinosWrappers::PreconditionSSOR());
+
+  Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1),1.2);
 
   pcout << std::endl;
 
@@ -1286,11 +1287,11 @@ void BoussinesqFlowProblem<dim>::solve ()
                                  1e-8*temperature_rhs.l2_norm());
     SolverCG<TrilinosWrappers::MPI::Vector>   cg (solver_control);
 
-    TrilinosWrappers::PreconditionSSOR preconditioner (temperature_matrix,
-                                                      1.2);
+    TrilinosWrappers::PreconditionSSOR preconditioner;
+    preconditioner.initialize (temperature_matrix, 1.2);
+
     cg.solve (temperature_matrix, temperature_solution,
-             temperature_rhs,
-             preconditioner);
+             temperature_rhs, preconditioner);
 
     TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
     temperature_constraints.distribute (localized_temperature_solution);
index aebf1abc02b83dd34f055ec9fb2e1e967636531d..280d2f436229f90765f5be5c6896157244dc9026 100755 (executable)
 
 
 #include <base/config.h>
-#include <boost/shared_ptr.hpp>
+#include <base/subscriptor.h>
 
 
 #ifdef DEAL_II_USE_TRILINOS
 
+#include <Teuchos_RefCountPtr.hpp>
+
 // forward declarations
 class Ifpack_Preconditioner;
 
@@ -36,43 +38,487 @@ namespace TrilinosWrappers
   class VectorBase;
   class SparseMatrix;
   
+/**  
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionJacobi : public Subscriptor
+  {
+    public:
+
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional flags to the
+                                        * preconditioner. The
+                                        * parameter <tt>omega</tt>
+                                        * specifies the relaxation
+                                        * parameter in the Jacobi
+                                        * preconditioner. The
+                                        * parameter
+                                        * <tt>min_diagonal</tt> can be
+                                        * used to make the application
+                                        * of the preconditioner also
+                                        * possible when some diagonal
+                                        * elements are zero. In a
+                                        * default application this
+                                        * would mean that we divide by
+                                        * zero, so by setting the
+                                        * parameter
+                                        * <tt>min_diagonal</tt> to a
+                                        * small nonzero value the SOR
+                                        * will work on a matrix that
+                                        * is not too far away from the
+                                        * one we want to
+                                        * treat.
+                                        */      
+      struct AdditionalData
+      {
+                                       /**
+                                       * Constructor. By default, set
+                                       * the damping parameter to
+                                       * one, and do not modify the
+                                       * diagonal.
+                                       */
+       AdditionalData (const double       omega = 1,
+                       const double       min_diagonal = 0);
+       
+                                       /**
+                                       * Relaxation parameter and
+                                       * minimal diagonal value.
+                                       */
+       double omega, min_diagonal;
+      };
+
+                                       /**
+                                        * Constructor. Does not do
+                                        * anything. The
+                                        * <tt>initialize</tt> function
+                                        * will have to set create the
+                                        * partitioner.
+                                        */
+      PreconditionJacobi ();
+
+                                       /**
+                                       * Take 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());
+      
+                                      /**
+                                       * Apply the preconditioner.
+                                       */
+      void vmult (VectorBase       &dst,
+                 const VectorBase &src) const;
+
+    private:
+      Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+      //std::auto_ptr<Ifpack_Preconditioner> preconditioner;
+      //boost::shared_ptr<Ifpack_Preconditioner> preconditioner;
+  };
+  
+
+  
+  
 /**  
  * @ingroup TrilinosWrappers
  * @ingroup Preconditioners
  * @author Wolfgang Bangerth, 2008
  */
-  class PreconditionSSOR
+  class PreconditionSSOR : public Subscriptor
   {
     public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional flags to the
+                                        * preconditioner. The
+                                        * parameter <tt>omega</tt>
+                                        * specifies the relaxation
+                                        * parameter in the SSOR
+                                        * preconditioner. The
+                                        * parameter
+                                        * <tt>min_diagonal</tt> can be
+                                        * used to make the application
+                                        * of the preconditioner also
+                                        * possible when some diagonal
+                                        * elements are zero. In a
+                                        * default application this
+                                        * would mean that we divide by
+                                        * zero, so by setting the
+                                        * parameter
+                                        * <tt>min_diagonal</tt> to a
+                                        * small nonzero value the SOR
+                                        * will work on a matrix that
+                                        * is not too far away from the
+                                        * one we want to
+                                        * treat. Finally,
+                                        * <tt>overlap</tt> governs the
+                                        * overlap of the partitions
+                                        * when the preconditioner runs
+                                        * in parallel, forming a
+                                        * so-called additive Schwarz
+                                        * preconditioner.
+                                        */      
+      struct AdditionalData
+      {
+                                       /**
+                                       * Constructor. By default, set
+                                       * the damping parameter to
+                                       * one, we do not modify the
+                                       * diagonal, and there is no
+                                       * overlap (i.e. in parallel,
+                                       * we run a BlockJacobi
+                                       * preconditioner, where each
+                                       * block is inverted
+                                       * approximately by an SSOR.
+                                       */
+       AdditionalData (const double       omega = 1,
+                       const double       min_diagonal = 0,
+                       const unsigned int overlap = 0);
+       
+                                       /**
+                                       * Relaxation parameter,
+                                       * minimal diagonal element,
+                                       * and overlap.
+                                       */
+       double omega, min_diagonal;
+       unsigned int overlap; 
+      };
+
+                                       /**
+                                        * Constructor. Does not do
+                                        * anything. The
+                                        * <tt>initialize</tt> function
+                                        * will have to set create the
+                                        * partitioner.
+                                        */
+      PreconditionSSOR ();
+
+                                       /**
+                                       * Take 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());
+      
+                                      /**
+                                       * Apply the preconditioner.
+                                       */
+      void vmult (VectorBase       &dst,
+                 const VectorBase &src) const;
+
+    private:
+      Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+  };
+  
+
+  
+  
+/**  
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionSOR : public Subscriptor
+  {
     public:
+
                                        /**
                                         * Standardized data struct to
                                         * pipe additional flags to the
+                                        * preconditioner. The
+                                        * parameter <tt>omega</tt>
+                                        * specifies the relaxation
+                                        * parameter in the SOR
+                                        * preconditioner. The
+                                        * parameter
+                                        * <tt>min_diagonal</tt> can be
+                                        * used to make the application
+                                        * of the preconditioner also
+                                        * possible when some diagonal
+                                        * elements are zero. In a
+                                        * default application this
+                                        * would mean that we divide by
+                                        * zero, so by setting the
+                                        * parameter
+                                        * <tt>min_diagonal</tt> to a
+                                        * small nonzero value the SOR
+                                        * will work on a matrix that
+                                        * is not too far away from the
+                                        * one we want to
+                                        * treat. Finally,
+                                        * <tt>overlap</tt> governs the
+                                        * overlap of the partitions
+                                        * when the preconditioner runs
+                                        * in parallel, forming a
+                                        * so-called additive Schwarz
                                         * preconditioner.
                                         */      
       struct AdditionalData
       {
-                                           /**
-                                            * Constructor. By default,
-                                            * set the damping parameter
-                                            * to one.
-                                            */
-          AdditionalData (const double omega = 1);
+                                       /**
+                                       * Constructor. By default, set
+                                       * the damping parameter to
+                                       * one, we do not modify the
+                                       * diagonal, and there is no
+                                       * overlap (i.e. in parallel,
+                                       * we run a BlockJacobi
+                                       * preconditioner, where each
+                                       * block is inverted
+                                       * approximately by an SOR.
+                                       */
+       AdditionalData (const double       omega = 1,
+                       const double       min_diagonal = 0,
+                       const unsigned int overlap = 0);
        
-                                           /**
-                                            * Relaxation parameter.
-                                            */
-          double omega;
+                                       /**
+                                       * Relaxation parameter,
+                                       * minimal diagonal element,
+                                       * and overlap.
+                                       */
+       double omega, min_diagonal;
+       unsigned int overlap; 
+      };
+
+                                       /**
+                                        * Constructor. Does not do
+                                        * anything. The
+                                        * <tt>initialize</tt> function
+                                        * will have to set create the
+                                        * partitioner.
+                                        */
+      PreconditionSOR ();
+
+                                       /**
+                                       * Take 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());
+      
+                                      /**
+                                       * Apply the preconditioner.
+                                       */
+      void vmult (VectorBase       &dst,
+                 const VectorBase &src) const;
+
+    private:
+      Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+  };
+
+
+
+/**  
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionIC : public Subscriptor
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional parameters
+                                        * to the preconditioner. The
+                                        * Trilinos IC decomposition
+                                        * allows for some fill-in, so
+                                        * it actually is a threshold
+                                        * incomplete Cholesky
+                                        * factorization. The amount of
+                                        * fill-in, and hence, the
+                                        * amount of memory used by
+                                        * this preconditioner, is
+                                        * controlled by the parameter
+                                        * <tt>ic_fill</tt>, which
+                                        * specifies this as a
+                                        * double. When forming the
+                                        * preconditioner, for certain
+                                        * problems bad conditioning
+                                        * (or just bad luck) can cause
+                                        * the preconditioner to be
+                                        * very poorly
+                                        * conditioned. Hence it can
+                                        * help to add diagonal
+                                        * perturbations to the
+                                        * original matrix and form the
+                                        * preconditioner for this
+                                        * slightly better
+                                        * matrix. <tt>ic_atol</tt> is
+                                        * an absolute perturbation
+                                        * that is added to the
+                                        * diagonal before forming the
+                                        * prec, and <tt>ic_rtol</tt>
+                                        * is a scaling factor $rtol
+                                        * \geq 1$. The last parameter
+                                        * specifies the overlap of the
+                                        * partitions when the
+                                        * preconditioner runs in
+                                        * parallel.
+                                        */      
+      struct AdditionalData
+      {
+                                       /**
+                                       * Constructor. By default, set
+                                       * the drop tolerance to 0, the
+                                       * level of extra fill-ins is
+                                       * set to be zero (just use the
+                                       * matrix structure, do not
+                                       * generate any additional
+                                       * fill-in), the tolerance
+                                       * level are 0 and 1,
+                                       * respectively, and the
+                                       * overlap in case of a
+                                       * parallel execution is
+                                       * zero. This overlap in a
+                                       * block-application of the IC
+                                       * in the parallel case makes
+                                       * the preconditioner a
+                                       * so-called additive Schwarz
+                                       * preconditioner.
+                                       */
+       AdditionalData (const double       ic_fill = 0.,
+                       const double       ic_atol = 0.,
+                       const double       ic_rtol = 1.,
+                       const unsigned int overlap = 0);
+       
+                                       /**
+                                       * IC parameters and overlap.
+                                       */
+       double ic_fill, ic_atol, ic_rtol;
+       unsigned int overlap;
       };
 
                                        /**
-                                        * Constructor. Take the matrix which
-                                        * is used to form the preconditioner,
-                                        * and additional flags if there are
-                                        * any.
+                                       * (Empty) constructor.
+                                       */
+      PreconditionIC ();
+
+                                       /**
+                                        * 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());
+      
+                                      /**
+                                       * Apply the preconditioner.
+                                       */
+      void vmult (VectorBase       &dst,
+                 const VectorBase &src) const;
+
+    private:
+      Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+  };
+
+
+
+/**  
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionILU : public Subscriptor
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional parameters
+                                        * to the preconditioner. The
+                                        * Trilinos ILU decomposition
+                                        * allows for some fill-in, so
+                                        * it actually is a threshold
+                                        * incomplete LU
+                                        * factorization. The amount of
+                                        * fill-in, and hence, the
+                                        * amount of memory used by
+                                        * this preconditioner, is
+                                        * controlled by the parameter
+                                        * <tt>ilu_fill</tt>, which
+                                        * specifies this as a
+                                        * double. When forming the
+                                        * preconditioner, for certain
+                                        * problems bad conditioning
+                                        * (or just bad luck) can cause
+                                        * the preconditioner to be
+                                        * very poorly
+                                        * conditioned. Hence it can
+                                        * help to add diagonal
+                                        * perturbations to the
+                                        * original matrix and form the
+                                        * preconditioner for this
+                                        * slightly better
+                                        * matrix. <tt>ilu_atol</tt> is
+                                        * an absolute perturbation
+                                        * that is added to the
+                                        * diagonal before forming the
+                                        * prec, and <tt>ilu_rtol</tt>
+                                        * is a scaling factor $rtol
+                                        * \geq 1$. The last parameter
+                                        * specifies the overlap of the
+                                        * partitions when the
+                                        * preconditioner runs in
+                                        * parallel.
+                                        */      
+      struct AdditionalData
+      {
+                                       /**
+                                       * Constructor. By default, the
+                                       * level of extra fill-ins is
+                                       * set to be zero (just use the
+                                       * matrix structure, do not
+                                       * generate any additional
+                                       * fill-in), the tolerance
+                                       * level are 0 and 1,
+                                       * respectively, and the
+                                       * overlap in case of a
+                                       * parallel execution is
+                                       * zero. This overlap in a
+                                       * block-application of the ILU
+                                       * in the parallel case makes
+                                       * the preconditioner a
+                                       * so-called additive Schwarz
+                                       * preconditioner.
+                                       */
+       AdditionalData (const double       ilu_fill = 0.,
+                       const double       ilu_atol = 0.,
+                       const double       ilu_rtol = 1.,
+                       const unsigned int overlap  = 0);
+       
+                                       /**
+                                       * ILU parameters and overlap.
+                                       */
+       double ilu_drop, ilu_fill, ilu_atol, ilu_rtol;
+       unsigned int overlap;
+      };
+
+                                       /**
+                                       * (Empty) constructor.
+                                       */
+      PreconditionILU ();
+
+                                       /**
+                                        * Initialize function. Takes
+                                        * the matrix which is used to
+                                        * form the preconditioner, and
+                                        * additional flags if there
+                                        * are any.
                                         */
-      PreconditionSSOR (const SparseMatrix   &matrix,
-                        const AdditionalData &additional_data = AdditionalData());
+      void initialize (const SparseMatrix   &matrix,
+                      const AdditionalData &additional_data = AdditionalData());
       
                                       /**
                                        * Apply the preconditioner.
@@ -81,7 +527,7 @@ namespace TrilinosWrappers
                  const VectorBase &src) const;
 
     private:
-      boost::shared_ptr<Ifpack_Preconditioner> preconditioner;
+      Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
   };
   
 
@@ -94,7 +540,7 @@ DEAL_II_NAMESPACE_CLOSE
 
 #endif // DEAL_II_USE_TRILINOS
 
-/*----------------------------   trilinos_precondition_base.h     ---------------------------*/
+/*----------------------------   trilinos_precondition.h     ---------------------------*/
 
 #endif
-/*----------------------------   trilinos_precondition_base.h     ---------------------------*/
+/*----------------------------   trilinos_precondition.h     ---------------------------*/
index 0d9eb8768f7fc1567dc307e07e413fc3e49073c3..7a96d6bc3a6d17274d6269430104638fea01ecc8 100755 (executable)
@@ -26,26 +26,96 @@ DEAL_II_NAMESPACE_OPEN
 namespace TrilinosWrappers
 {
 
+  PreconditionJacobi::AdditionalData::
+  AdditionalData (const double omega,
+                 const double min_diagonal)
+                  :
+                  omega (omega),
+                 min_diagonal (min_diagonal)
+  {}
+
+  
+
+  PreconditionJacobi::PreconditionJacobi ()
+  {}
+
+
+
+  void
+  PreconditionJacobi::initialize (const SparseMatrix   &matrix,
+                                 const AdditionalData &additional_data)
+  {
+
+    preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation", &*matrix.matrix, 0));
+    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+                                              "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: type", "Jacobi");
+    parameter_list.set ("relaxation: damping factor", additional_data.omega);
+    parameter_list.set ("relaxation: min diagonal value", 
+                       additional_data.min_diagonal);
+    
+    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));
+  }
+
+
+
+  void
+  PreconditionJacobi::vmult (VectorBase       &dst,
+                            const VectorBase &src) const
+  {
+    const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
   PreconditionSSOR::AdditionalData::
-  AdditionalData (const double omega)
+  AdditionalData (const double       omega,
+                 const double       min_diagonal,
+                 const unsigned int overlap)
                   :
-                  omega (omega)
+                  omega        (omega),
+                 min_diagonal (min_diagonal),
+                 overlap      (overlap)
   {}
 
   
-  PreconditionSSOR::PreconditionSSOR (const SparseMatrix   &matrix,
-                                     const AdditionalData &additional_data)
-                 :
-                 preconditioner (Ifpack().Create ("point relaxation",
-                                                  &*matrix.matrix, 0))
+
+  PreconditionSSOR::PreconditionSSOR ()
+  {}
+
+
+
+  void
+  PreconditionSSOR::initialize (const SparseMatrix   &matrix,
+                               const AdditionalData &additional_data)
   {
-    Assert (preconditioner != 0, ExcMessage ("Trilinos could not create this preconditioner"));
+
+    preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation",
+                                                   &*matrix.matrix, 
+                                                   additional_data.overlap));
+    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+                                              "preconditioner"));
 
     int ierr;
 
     Teuchos::ParameterList parameter_list;
     parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
+    parameter_list.set ("relaxation: min diagonal value", 
+                       additional_data.min_diagonal);
+    parameter_list.set ("schwarz: combine mode", "Add");
     
     ierr = preconditioner->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
@@ -58,11 +128,188 @@ namespace TrilinosWrappers
   }
 
 
+
   void
   PreconditionSSOR::vmult (VectorBase       &dst,
                           const VectorBase &src) const
   {
-    preconditioner->ApplyInverse (*src.vector, *dst.vector);
+    const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  PreconditionSOR::AdditionalData::
+  AdditionalData (const double       omega,
+                 const double       min_diagonal,
+                 const unsigned int overlap)
+                  :
+                  omega        (omega),
+                 min_diagonal (min_diagonal),
+                 overlap      (overlap)
+  {}
+
+  
+
+  PreconditionSOR::PreconditionSOR ()
+  {}
+
+
+
+  void
+  PreconditionSOR::initialize (const SparseMatrix   &matrix,
+                              const AdditionalData &additional_data)
+  {
+
+    preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation",
+                                                   &*matrix.matrix, 
+                                                   additional_data.overlap));
+    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+                                              "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: type", "Gauss-Seidel");
+    parameter_list.set ("relaxation: damping factor", additional_data.omega);
+    parameter_list.set ("relaxation: min diagonal value", 
+                       additional_data.min_diagonal);
+    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));
+  }
+
+
+
+  void
+  PreconditionSOR::vmult (VectorBase       &dst,
+                         const VectorBase &src) const
+  {
+    const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  PreconditionIC::AdditionalData::
+  AdditionalData (const double       ic_fill,
+                 const double       ic_atol,
+                 const double       ic_rtol,
+                 const unsigned int overlap)
+                  :
+                  ic_fill (ic_fill),
+                 ic_atol (ic_atol),
+                 ic_rtol (ic_rtol),
+                 overlap (overlap)
+  {}
+
+
+
+  PreconditionIC::PreconditionIC ()
+  {}
+
+
+
+  void
+  PreconditionIC::initialize (const SparseMatrix   &matrix,
+                             const AdditionalData &additional_data)
+  {
+    preconditioner = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix, 
+                                                   additional_data.overlap));
+    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+                                              "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("fact: level-of-fill",additional_data.ic_fill); 
+    parameter_list.set ("fact: absolute threshold",additional_data.ic_atol); 
+    parameter_list.set ("fact: relative threshold",additional_data.ic_rtol); 
+    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));
+  }
+
+
+
+  void
+  PreconditionIC::vmult (VectorBase       &dst,
+                        const VectorBase &src) const
+  {
+    const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  PreconditionILU::AdditionalData::
+  AdditionalData (const double       ilu_fill,
+                 const double       ilu_atol,
+                 const double       ilu_rtol,
+                 const unsigned int overlap)
+                  :
+                  ilu_fill (ilu_fill),
+                 ilu_atol (ilu_atol),
+                 ilu_rtol (ilu_rtol),
+                 overlap  (overlap)
+  {}
+
+
+
+  PreconditionILU::PreconditionILU ()
+  {}
+
+
+
+  void
+  PreconditionILU::initialize (const SparseMatrix   &matrix,
+                              const AdditionalData &additional_data)
+  {
+    preconditioner = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix, 
+                                                   additional_data.overlap));
+    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+                                              "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("fact: level-of-fill",additional_data.ilu_fill); 
+    parameter_list.set ("fact: absolute threshold",additional_data.ilu_atol); 
+    parameter_list.set ("fact: relative threshold",additional_data.ilu_rtol); 
+    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));
+  }
+
+
+
+  void
+  PreconditionILU::vmult (VectorBase       &dst,
+                         const VectorBase &src) const
+  {
+    const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 }
 

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.