]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Made the interface to Trilinos preconditioners more generic. Now AMG can be controlle...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 Oct 2008 13:09:28 +0000 (13:09 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 Oct 2008 13:09:28 +0000 (13:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@17145 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
deal.II/lac/source/trilinos_solver.cc

index b17f68bf530f08dabd127cd47de762060b29f989..6ef4629c47bbac6f0fa4fa97dcae6b6172639d14 100644 (file)
                                 // classes. In particular, we will need
                                 // interfaces to the matrix and vector
                                 // classes based on Trilinos as well as
-                                // generic preconditioners and the Trilinos
-                                // Algebraic Multigrid (AMG) preconditioner
-                                // that we will use for the $A$ block of the
-                                // Stokes matrix:
+                                // Trilinos preconditioners:
 #include <lac/trilinos_block_vector.h>
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_block_sparse_matrix.h>
 #include <lac/trilinos_precondition.h>
-#include <lac/trilinos_precondition_amg.h>
 
                                 // Finally, here are two C++ headers that
                                 // haven't been included yet by one of the
@@ -1231,7 +1227,8 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
   DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components, 
                                    null_space);
   Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
-                                true, true, 5e-2, null_space, false);
+                                TrilinosWrappers::PreconditionAMG::AdditionalData
+                                  (true, true, 5e-2, null_space, 0, false));
 
                                   // TODO: we could throw away the (0,0)
                                   // block here since things have been
index 853aecf245d402184786b818f971c1d9bd704207..d1b125f44ecbd6e3e14c2b6d17afa5fde511c119 100644 (file)
@@ -28,7 +28,6 @@
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_block_sparse_matrix.h>
 #include <lac/trilinos_precondition.h>
-#include <lac/trilinos_precondition_amg.h>
 
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
@@ -805,7 +804,8 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
                                    null_space);
   
   Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
-                                true, true, 5e-2, null_space, false);
+                                TrilinosWrappers::PreconditionAMG::AdditionalData
+                                  (true, true, 5e-2, null_space, 0, false));
 
   Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionIC>
                                    (new TrilinosWrappers::PreconditionIC());
index 20932338642f47b2b06b60d14780cfd5ecfba59b..54db649b4d172ae4c82d867912d4bfdada364b8c 100755 (executable)
 
 #include <base/config.h>
 #include <base/subscriptor.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+
+#include <boost/shared_ptr.hpp>
 
 
 #ifdef DEAL_II_USE_TRILINOS
 
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+#  include <Epetra_MpiComm.h>
+#else
+#  include <Epetra_SerialComm.h>
+#endif
+#include <Epetra_Map.h>
+
 #include <Teuchos_RefCountPtr.hpp>
+#include <Epetra_Operator.h>
 
 // forward declarations
 class Ifpack_Preconditioner;
+namespace ML_Epetra
+{
+  class MultiLevelPreconditioner;
+}
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -68,6 +85,17 @@ namespace TrilinosWrappers
                                         * sparse matrix.
                                         */
       PreconditionBase ();
+
+                                       /**
+                                        * Constructor. Does not do
+                                        * anything. The
+                                        * <tt>initialize</tt> function
+                                        * of the derived classes will
+                                        * have to create the
+                                        * preconditioner from a given
+                                        * sparse matrix.
+                                        */
+      ~PreconditionBase ();
       
                                       /**
                                        * Apply the preconditioner.
@@ -75,6 +103,16 @@ namespace TrilinosWrappers
       void vmult (VectorBase       &dst,
                  const VectorBase &src) const;
 
+                                      /**
+                                       * Apply the preconditioner on
+                                       * deal.II data structures
+                                       * instead of the ones provided
+                                       * in the Trilinos wrapper
+                                       * class.
+                                       */
+      void vmult (dealii::Vector<double>       &dst,
+                 const dealii::Vector<double> &src) const;
+
                                        /**
                                        * Exception.
                                        */
@@ -88,13 +126,33 @@ namespace TrilinosWrappers
       friend class SolverBase;
       friend class SolverSaddlePoint;
 
-      //protected:
+    protected:
                                       /**
                                        * This is a pointer to the
-                                       * preconditioner object.
+                                       * preconditioner object that
+                                       * is used when applying the
+                                       * preconditioner.
+                                       */
+      Teuchos::RCP<const Epetra_Operator> preconditioner;
+
+                                       /**
+                                       * Internal communication
+                                       * pattern in case the matrix
+                                       * needs to be copied from
+                                       * deal.II format.
                                        */
-        Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+      Epetra_MpiComm     communicator;
+#else
+      Epetra_SerialComm  communicator;
+#endif
 
+                                       /**
+                                       * Internal Trilinos map in
+                                       * case the matrix needs to be
+                                       * copied from deal.II format.
+                                       */
+      std::auto_ptr<Epetra_Map>   map;
   };
 
   
@@ -156,10 +214,25 @@ namespace TrilinosWrappers
                        const double       min_diagonal = 0);
        
                                        /**
-                                       * Relaxation parameter and
-                                       * minimal diagonal value.
+                                       * This specifies the
+                                       * relaxation parameter in the
+                                       * Jacobi preconditioner.
                                        */
-       double omega, min_diagonal;
+       double omega;
+
+                                       /**
+                                       * This specifies the minimum
+                                       * value the diagonal elements
+                                       * should have. This might be
+                                       * necessary when the Jacobi
+                                       * preconditioner is used on
+                                       * matrices with zero diagonal
+                                       * elements. In that case, a
+                                       * straight-forward application
+                                       * would not be possible since
+                                       * we would divide by zero.
+                                       */
+       double min_diagonal;
       };
 
                                        /**
@@ -171,6 +244,14 @@ namespace TrilinosWrappers
                                        */
       void initialize (const SparseMatrix   &matrix,
                       const AdditionalData &additional_data = AdditionalData());
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
   };
   
 
@@ -257,11 +338,34 @@ namespace TrilinosWrappers
                        const unsigned int overlap = 0);
        
                                        /**
-                                       * Relaxation parameter,
-                                       * minimal diagonal element,
-                                       * and overlap.
+                                       * This specifies the (over-)
+                                       * relaxation parameter in the
+                                       * SSOR preconditioner.
+                                       */
+       double omega;
+
+                                       /**
+                                       * This specifies the minimum
+                                       * value the diagonal elements
+                                       * should have. This might be
+                                       * necessary when the SSOR
+                                       * preconditioner is used on
+                                       * matrices with zero diagonal
+                                       * elements. In that case, a
+                                       * straight-forward application
+                                       * would not be possible since
+                                       * we divide by the diagonal
+                                       * element.
+                                       */
+       double min_diagonal;
+
+                                       /**
+                                       * This determines how large
+                                       * the overlap of the local
+                                       * matrix portions on each
+                                       * processor in a parallel
+                                       * application should be.
                                        */
-       double omega, min_diagonal;
        unsigned int overlap; 
       };
 
@@ -276,6 +380,14 @@ namespace TrilinosWrappers
                                        */
       void initialize (const SparseMatrix   &matrix,
                       const AdditionalData &additional_data = AdditionalData());
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
   };
   
 
@@ -360,13 +472,36 @@ namespace TrilinosWrappers
        AdditionalData (const double       omega = 1,
                        const double       min_diagonal = 0,
                        const unsigned int overlap = 0);
-       
+
                                        /**
-                                       * Relaxation parameter,
-                                       * minimal diagonal element,
-                                       * and overlap.
+                                       * This specifies the (over-)
+                                       * relaxation parameter in the
+                                       * SOR preconditioner.
+                                       */
+       double omega;
+
+                                       /**
+                                       * This specifies the minimum
+                                       * value the diagonal elements
+                                       * should have. This might be
+                                       * necessary when the SOR
+                                       * preconditioner is used on
+                                       * matrices with zero diagonal
+                                       * elements. In that case, a
+                                       * straight-forward application
+                                       * would not be possible since
+                                       * we divide by the diagonal
+                                       * element.
+                                       */
+       double min_diagonal;
+
+                                       /**
+                                       * This determines how large
+                                       * the overlap of the local
+                                       * matrix portions on each
+                                       * processor in a parallel
+                                       * application should be.
                                        */
-       double omega, min_diagonal;
        unsigned int overlap; 
       };
 
@@ -380,7 +515,15 @@ namespace TrilinosWrappers
                                         * are any.
                                        */
       void initialize (const SparseMatrix   &matrix,
-                      const AdditionalData &additional_data = AdditionalData());
+                      const AdditionalData &additional_data = AdditionalData()); 
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
   };
 
 
@@ -490,13 +633,55 @@ namespace TrilinosWrappers
                        const double       ic_atol = 0.,
                        const double       ic_rtol = 1.,
                        const unsigned int overlap = 0);
-       
+
                                        /**
-                                       * IC parameters and overlap.
+                                       * This specifies the amount of
+                                       * additional fill-in elements
+                                       * besides the sparse matrix
+                                       * structure. When
+                                       * <tt>ic_fill</tt> is large,
+                                       * this means that many
+                                       * fill-ins will be added, so
+                                       * that the IC preconditioner
+                                       * comes closer to a direct
+                                       * sparse Cholesky
+                                       * decomposition. Note,
+                                       * however, that this will
+                                       * drastically increase the
+                                       * memory requirement,
+                                       * especially when the
+                                       * preconditioner is used in
+                                       * 3D.
                                        */
        unsigned int ic_fill;
-       double ic_atol, ic_rtol;
-       unsigned int overlap;
+
+                                       /**
+                                       * This specifies the amount of
+                                       * an absolute perturbation
+                                       * that will be added to the
+                                       * diagonal of the matrix,
+                                       * which sometimes can help to
+                                       * get better preconditioners.
+                                       */
+       double ic_atol;
+
+                                       /**
+                                       * This specifies the factor by
+                                       * which the diagonal of the
+                                       * matrix will be scaled, which
+                                       * sometimes can help to get
+                                       * better preconditioners.
+                                       */
+       double ic_rtol;
+
+                                       /**
+                                       * This determines how large
+                                       * the overlap of the local
+                                       * matrix portions on each
+                                       * processor in a parallel
+                                       * application should be.
+                                       */
+       unsigned int overlap; 
       };
 
                                        /**
@@ -508,6 +693,14 @@ namespace TrilinosWrappers
                                         */
       void initialize (const SparseMatrix   &matrix,
                       const AdditionalData &additional_data = AdditionalData());
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
   };
 
 
@@ -616,13 +809,55 @@ namespace TrilinosWrappers
                        const double       ilu_atol = 0.,
                        const double       ilu_rtol = 1.,
                        const unsigned int overlap  = 0);
-       
+
                                        /**
-                                       * ILU parameters and overlap.
+                                       * This specifies the amount of
+                                       * additional fill-in elements
+                                       * besides the sparse matrix
+                                       * structure. When
+                                       * <tt>ilu_fill</tt> is large,
+                                       * this means that many
+                                       * fill-ins will be added, so
+                                       * that the ILU preconditioner
+                                       * comes closer to a (direct)
+                                       * sparse LU
+                                       * decomposition. Note,
+                                       * however, that this will
+                                       * drastically increase the
+                                       * memory requirement,
+                                       * especially when the
+                                       * preconditioner is used in
+                                       * 3D.
                                        */
        unsigned int ilu_fill;
-       double ilu_atol, ilu_rtol;
-       unsigned int overlap;
+
+                                       /**
+                                       * This specifies the amount of
+                                       * an absolute perturbation
+                                       * that will be added to the
+                                       * diagonal of the matrix,
+                                       * which sometimes can help to
+                                       * get better preconditioners.
+                                       */
+       double ilu_atol;
+
+                                       /**
+                                       * This specifies the factor by
+                                       * which the diagonal of the
+                                       * matrix will be scaled, which
+                                       * sometimes can help to get
+                                       * better preconditioners.
+                                       */
+       double ilu_rtol;
+
+                                       /**
+                                       * This determines how large
+                                       * the overlap of the local
+                                       * matrix portions on each
+                                       * processor in a parallel
+                                       * application should be.
+                                       */
+       unsigned int overlap; 
       };
 
                                        /**
@@ -634,6 +869,14 @@ namespace TrilinosWrappers
                                         */
       void initialize (const SparseMatrix   &matrix,
                       const AdditionalData &additional_data = AdditionalData());
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
   };
 
 
@@ -673,8 +916,13 @@ namespace TrilinosWrappers
                                        */
        AdditionalData (const unsigned int overlap  = 0);
        
+
                                        /**
-                                       * Block direct parameters and overlap.
+                                       * This determines how large
+                                       * the overlap of the local
+                                       * matrix portions on each
+                                       * processor in a parallel
+                                       * application should be.
                                        */
        unsigned int overlap;
       };
@@ -688,9 +936,219 @@ namespace TrilinosWrappers
                                         */
       void initialize (const SparseMatrix   &matrix,
                       const AdditionalData &additional_data = AdditionalData());
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
   };
   
 
+
+/**
+ * This class implements an algebraic multigrid (AMG) preconditioner
+ * based on the Trilinos ML implementation, which is a black-box
+ * preconditioner that works well for many PDE-based linear problems.
+ * What this class does is twofold.  When the initialize() function is
+ * invoked, a ML preconditioner object is created based on the matrix
+ * that we want the preconditioner to be based on. A call of the
+ * respective <code>vmult</code> function does call the respective
+ * operation in the Trilinos package, where it is called
+ * <code>ApplyInverse</code>. Use of this class is explained in the
+ * @ref step_31 "step-31" tutorial program.
+ *
+ * There are a few pecularities in initialize(). Since the Trilinos
+ * objects we want to use are heavily dependent on Epetra objects, the
+ * fundamental construction routines for vectors and matrices in
+ * Trilinos, we do a copy of our deal.II preconditioner matrix to a
+ * Epetra matrix. This is of course not optimal, but for the time
+ * being there is no direct support for our data interface.  When
+ * doing this time-consuming operation, we can still profit from the
+ * fact that some of the entries in the preconditioner matrix are zero
+ * and hence can be neglected.
+ *
+ * The implementation is able to distinguish between matrices from
+ * elliptic problems and convection dominated problems. We use the
+ * standard options provided by Trilinos ML for elliptic problems,
+ * except that we use a Chebyshev smoother instead of a symmetric
+ * Gauss-Seidel smoother.  For most elliptic problems, Chebyshev
+ * provides a better damping of high frequencies (in the algebraic
+ * sense) than Gauss-Seidel (SSOR).
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionAMG : public PreconditionBase
+  {
+    public:
+
+    struct AdditionalData
+      {
+                                       /**
+                                       * Constructor. By default, we
+                                       * pretend to work on elliptic
+                                       * problems with linear finite
+                                       * elements on a scalar
+                                       * equation.
+                                       */
+       AdditionalData (const bool                             elliptic = true,
+                       const bool                             higher_order_elements = false,
+                       const double                           aggregation_threshold = 1e-4,
+                       const std::vector<std::vector<bool> > &null_space = std::vector<std::vector<bool> > (1),
+                       const unsigned int                     smoother_overlap = 0,
+                       const bool                             output_details = false);
+
+                                      /**
+                                       * Determines whether the AMG
+                                       * preconditioner should be
+                                       * optimized for elliptic
+                                       * problems (ML option smoothed
+                                       * aggregation SA, using a
+                                       * Chebyshev smoother) or for
+                                       * non-elliptic problems (ML
+                                       * option non-symmetric
+                                       * smoothed aggregation NSSA,
+                                       * smoother is SSOR with
+                                       * underrelaxation).
+                                       */
+       bool elliptic;
+
+                                      /**
+                                       * Determines whether the
+                                       * matrix that the
+                                       * preconditioner is built upon
+                                       * is generated from linear or
+                                       * higher-order elements.
+                                       */
+       bool higher_order_elements;
+
+                                      /**
+                                       * This threshold tells the AMG
+                                       * setup how the coarsening
+                                       * should be performed. In the
+                                       * AMG used by ML, all points
+                                       * that strongly couple with
+                                       * the tentative coarse-level
+                                       * point form one
+                                       * aggregate. The term
+                                       * <em>strong coupling</em> is
+                                       * controlled by the variable
+                                       * <tt>aggregation_threshold</tt>,
+                                       * meaning that all elements
+                                       * that are not smaller than
+                                       * <tt>aggregation_threshold</tt>
+                                       * times the diagonal element
+                                       * do couple strongly.
+                                       */
+       double aggregation_threshold;
+
+                                      /**
+                                       * Specifies the constant modes
+                                       * (near null space) of the
+                                       * matrix. This parameter tells
+                                       * AMG whether we work on a
+                                       * scalar equation (where the
+                                       * near null space only
+                                       * consists of ones) or on a
+                                       * vector-valued equation.
+                                       */
+       std::vector<std::vector<bool> > null_space;
+
+                                      /**
+                                       * Determines the overlap in
+                                       * the SSOR/Chebyshev error
+                                       * smoother when run in
+                                       * parallel.
+                                       */
+       unsigned int smoother_overlap;
+
+                                      /**
+                                       * If this flag is set to
+                                       * <tt>true</tt>, then internal
+                                       * information from the ML
+                                       * preconditioner is printed to
+                                       * screen. This can be useful
+                                       * when debugging the
+                                       * preconditioner.
+                                       */
+       bool output_details;
+      };
+
+                                      /**
+                                       * Let Trilinos compute a
+                                       * multilevel hierarchy for the
+                                       * solution of a linear system
+                                       * with the given matrix. The
+                                       * function uses the matrix
+                                       * format specified in
+                                       * TrilinosWrappers::SparseMatrix.
+                                       */
+      void initialize (const SparseMatrix                    &matrix,
+                      const AdditionalData &additional_data = AdditionalData());
+
+                                      /**
+                                       * Let Trilinos compute a
+                                       * multilevel hierarchy for the
+                                       * solution of a linear system
+                                       * with the given matrix. This
+                                       * function takes a deal.ii
+                                       * matrix and copies the
+                                       * content into a Trilinos
+                                       * matrix, so the function can
+                                       * be considered rather
+                                       * inefficient.
+                                       */
+      void initialize (const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+                      const AdditionalData                 &additional_data = AdditionalData(),
+                      const double                          drop_tolerance = 1e-13);
+
+                                      /**
+                                       * This function can be used
+                                       * for a faster recalculation
+                                       * of the preconditioner
+                                       * construction when the matrix
+                                       * entries underlying the
+                                       * preconditioner have changed,
+                                       * but the matrix sparsity
+                                       * pattern has remained the
+                                       * same. What this function
+                                       * does is taking the already
+                                       * generated coarsening
+                                       * structure, computing the AMG
+                                       * prolongation and restriction
+                                       * according to a smoothed
+                                       * aggregation strategy and
+                                       * then building the whole
+                                       * multilevel hiearchy. This
+                                       * function can be considerably
+                                       * faster than the initialize
+                                       * function, since the
+                                       * coarsening pattern is
+                                       * usually the most difficult
+                                       * thing to do when setting up
+                                       * the AMG ML preconditioner.
+                                       */
+      void reinit ();
+
+    private:
+
+                                      /**
+                                       * A pointer to the
+                                       * preconditioner object.
+                                       */
+      Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> multilevel_operator;
+
+                                      /**
+                                       * A copy of the deal.II matrix
+                                       * into Trilinos format.
+                                       */
+      boost::shared_ptr<SparseMatrix> Matrix;
+  };
+
 }
 
 /*@}*/
index 9b4c4ae03c664fb11cd20179c96d545e864a2945..b09cae729504a988fbf2761cfb49de6582fd478f 100755 (executable)
 
 #include <Ifpack.h>
 #include <Teuchos_ParameterList.hpp>
+#include <Epetra_Vector.h>
+#include <Epetra_MultiVector.h>
+#include <Epetra_Vector.h>
+#include <ml_include.h>
+#include <ml_MultiLevelPreconditioner.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -28,10 +33,21 @@ namespace TrilinosWrappers
 {
 
   PreconditionBase::PreconditionBase()
+                    :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                    communicator (MPI_COMM_WORLD)
+#endif
   {}
 
 
 
+  PreconditionBase::~PreconditionBase()
+  {
+    preconditioner.release();
+  }
+
+
+
   void
   PreconditionBase::vmult (VectorBase       &dst,
                           const VectorBase &src) const
@@ -47,6 +63,45 @@ namespace TrilinosWrappers
 
 
 
+                                       // For the implementation of
+                                       // the <code>vmult</code>
+                                       // function we note that
+                                       // invoking a call of the
+                                       // Trilinos preconditioner
+                                       // requires us to use Epetra
+                                       // vectors as well. It is
+                                       // faster to provide a view,
+                                       // i.e., feed Trilinos with a
+                                       // pointer to the data, so we
+                                       // avoid copying the content
+                                       // of the vectors during the
+                                       // iteration. In the
+                                       // declaration of the right
+                                       // hand side, we need to cast
+                                       // the source vector (that is
+                                       // <code>const</code> in all
+                                       // deal.II calls) to
+                                       // non-constant value, as this
+                                       // is the way Trilinos wants
+                                       // to have them.
+  void PreconditionBase::vmult (dealii::Vector<double>       &dst,
+                               const dealii::Vector<double> &src) const
+  {
+    
+    Epetra_Vector LHS (View, preconditioner->OperatorDomainMap(),
+                      dst.begin());
+    Epetra_Vector RHS (View, preconditioner->OperatorRangeMap(),
+                      const_cast<double*>(src.begin()));
+  
+    const int res = preconditioner->ApplyInverse (RHS, LHS);
+  
+    Assert (res == 0,
+           ExcMessage ("Trilinos AMG MultiLevel preconditioner returned "
+                       "with an error!"));
+  }
+
+
+
 /* -------------------------- PreconditionJacobi -------------------------- */
 
   PreconditionJacobi::AdditionalData::
@@ -64,10 +119,13 @@ namespace TrilinosWrappers
                                  const AdditionalData &additional_data)
   {
     preconditioner.release();
+    ifpack.release();
 
-    preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation", &*matrix.matrix, 0));
-    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
-                                              "preconditioner"));
+    ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation", 
+                                           &*matrix.matrix, 0));
+
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                      "preconditioner"));
 
     int ierr;
 
@@ -78,14 +136,16 @@ namespace TrilinosWrappers
     parameter_list.set ("relaxation: min diagonal value", 
                        additional_data.min_diagonal);
     
-    ierr = preconditioner->SetParameters(parameter_list);
+    ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Initialize();
+    ierr = ifpack->Initialize();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Compute();
+    ierr = ifpack->Compute();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
   }
 
 
@@ -109,12 +169,14 @@ namespace TrilinosWrappers
                                const AdditionalData &additional_data)
   {
     preconditioner.release();
+    ifpack.release();
 
-    preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation",
-                                                   &*matrix.matrix, 
-                                                   additional_data.overlap));
-    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
-                                              "preconditioner"));
+    ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
+                                           &*matrix.matrix, 
+                                           additional_data.overlap));
+
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                      "preconditioner"));
 
     int ierr;
 
@@ -126,14 +188,16 @@ namespace TrilinosWrappers
                        additional_data.min_diagonal);
     parameter_list.set ("schwarz: combine mode", "Add");
     
-    ierr = preconditioner->SetParameters(parameter_list);
+    ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Initialize();
+    ierr = ifpack->Initialize();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Compute();
+    ierr = ifpack->Compute();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
   }
 
 
@@ -157,12 +221,14 @@ namespace TrilinosWrappers
                               const AdditionalData &additional_data)
   {
     preconditioner.release();
+    ifpack.release();
 
-    preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation",
+    ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
                                                    &*matrix.matrix, 
                                                    additional_data.overlap));
-    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
-                                              "preconditioner"));
+
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                      "preconditioner"));
 
     int ierr;
 
@@ -174,14 +240,16 @@ namespace TrilinosWrappers
                        additional_data.min_diagonal);
     parameter_list.set ("schwarz: combine mode", "Add");
     
-    ierr = preconditioner->SetParameters(parameter_list);
+    ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Initialize();
+    ierr = ifpack->Initialize();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Compute();
+    ierr = ifpack->Compute();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
   }
 
 
@@ -207,10 +275,12 @@ namespace TrilinosWrappers
                              const AdditionalData &additional_data)
   {
     preconditioner.release();
+    ifpack.release();
 
-    preconditioner = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix, 
-                                                   additional_data.overlap));
-    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+    ifpack = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix, 
+                                           additional_data.overlap));
+
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
                                               "preconditioner"));
 
     int ierr;
@@ -221,14 +291,16 @@ namespace TrilinosWrappers
     parameter_list.set ("fact: relative threshold",additional_data.ic_rtol); 
     parameter_list.set ("schwarz: combine mode", "Add");
     
-    ierr = preconditioner->SetParameters(parameter_list);
+    ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Initialize();
+    ierr = ifpack->Initialize();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Compute();
+    ierr = ifpack->Compute();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
   }
 
 
@@ -254,10 +326,12 @@ namespace TrilinosWrappers
                               const AdditionalData &additional_data)
   {
     preconditioner.release();
+    ifpack.release();
 
-    preconditioner = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix, 
-                                                   additional_data.overlap));
-    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+    ifpack = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix, 
+                                           additional_data.overlap));
+
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
                                               "preconditioner"));
 
     int ierr;
@@ -268,14 +342,16 @@ namespace TrilinosWrappers
     parameter_list.set ("fact: relative threshold",additional_data.ilu_rtol); 
     parameter_list.set ("schwarz: combine mode", "Add");
     
-    ierr = preconditioner->SetParameters(parameter_list);
+    ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Initialize();
+    ierr = ifpack->Initialize();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Compute();
+    ierr = ifpack->Compute();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
   }
 
 
@@ -295,25 +371,164 @@ namespace TrilinosWrappers
                                       const AdditionalData &additional_data)
   {
     preconditioner.release();
+    ifpack.release();
 
-    preconditioner = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix, 
-                                                   additional_data.overlap));
-    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
-                                              "preconditioner"));
+    ifpack = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix, 
+                                           additional_data.overlap));
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                      "preconditioner"));
 
     int ierr;
 
     Teuchos::ParameterList parameter_list;
     parameter_list.set ("schwarz: combine mode", "Add");
     
-    ierr = preconditioner->SetParameters(parameter_list);
+    ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Initialize();
+    ierr = ifpack->Initialize();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    ierr = preconditioner->Compute();
+    ierr = ifpack->Compute();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
+  }
+
+
+
+/* -------------------------- PreconditionAMG -------------------------- */
+
+  PreconditionAMG::AdditionalData::
+  AdditionalData (const bool                             elliptic,
+                 const bool                             higher_order_elements,
+                 const double                           aggregation_threshold,
+                 const std::vector<std::vector<bool> > &null_space,
+                 const unsigned int                     smoother_overlap,
+                 const bool                             output_details)
+                  :
+                  elliptic (elliptic),
+                 higher_order_elements (higher_order_elements),
+                 aggregation_threshold (aggregation_threshold),
+                 null_space (null_space),
+                 smoother_overlap (smoother_overlap),
+                 output_details (output_details)
+  {}
+
+
+
+  void
+  PreconditionAMG:: initialize (const SparseMatrix   &matrix,
+                               const AdditionalData &additional_data)
+  {
+    preconditioner.release();
+    multilevel_operator.release();
+
+    const unsigned int n_rows = matrix.m();
+    const unsigned int null_space_dimension = additional_data.null_space.size();
+
+                                       // Build the AMG preconditioner.
+    Teuchos::ParameterList parameter_list;
+  
+    if (additional_data.elliptic == true)
+      {
+       ML_Epetra::SetDefaults("SA",parameter_list);
+       parameter_list.set("smoother: type", "Chebyshev");
+       parameter_list.set("smoother: sweeps", 4);
+      }
+    else
+      {
+       ML_Epetra::SetDefaults("NSSA",parameter_list);
+       parameter_list.set("aggregation: type", "Uncoupled");
+       parameter_list.set("aggregation: block scaling", true);
+      }
+  
+    parameter_list.set("smoother: ifpack overlap", 
+                      (int)additional_data.smoother_overlap);
+    parameter_list.set("aggregation: threshold", 
+                      additional_data.aggregation_threshold);
+    
+    if (additional_data.output_details)
+      parameter_list.set("ML output", 10);
+    else
+      parameter_list.set("ML output", 0);
+  
+    if (additional_data.higher_order_elements)
+      parameter_list.set("aggregation: type", "MIS");
+
+    const Epetra_Map * domain_map = &(matrix.matrix->DomainMap());
+    
+    Epetra_MultiVector null_space_modes (*domain_map, null_space_dimension);
+  
+    if (null_space_dimension > 1)
+      {
+       Assert (n_rows == additional_data.null_space[0].size(),
+               ExcDimensionMismatch(n_rows,
+                                    additional_data.null_space[0].size()));
+       Assert (n_rows == (unsigned int)null_space_modes.GlobalLength(),
+               ExcDimensionMismatch(n_rows,
+                                    null_space_modes.GlobalLength()));
+
+       const unsigned int my_size = domain_map->NumMyElements();
+       Assert (my_size == (unsigned int)domain_map->MaxLID()+1,
+               ExcDimensionMismatch (my_size, domain_map->MaxLID()+1));
+       
+                                       // Reshape null space as a
+                                       // contiguous vector of
+                                       // doubles so that Trilinos
+                                       // can read from it.
+       for (unsigned int d=0; d<null_space_dimension; ++d)
+         for (unsigned int row=0; row<my_size; ++row)
+           {
+             int global_row_id = domain_map->GID(row);
+             null_space_modes.ReplaceMyValue(row, d, 
+                                (double)additional_data.null_space[d][global_row_id]);
+           }
+  
+       parameter_list.set("null space: type", "pre-computed");
+       parameter_list.set("null space: dimension", null_space_modes.NumVectors());
+       parameter_list.set("null space: vectors", null_space_modes.Values());
+      }
+
+    multilevel_operator = Teuchos::rcp (new ML_Epetra::MultiLevelPreconditioner(
+                                       *matrix.matrix, parameter_list, true));
+
+    if (additional_data.output_details)
+      multilevel_operator->PrintUnused(0);
+
+    preconditioner = Teuchos::rcp (&*multilevel_operator);
+  }
+
+
+
+  void
+  PreconditionAMG::
+  initialize (const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+             const AdditionalData                 &additional_data,
+             const double                          drop_tolerance)
+  {
+    const unsigned int n_rows = deal_ii_sparse_matrix.m();
+  
+                                       // Init Epetra Matrix, avoid
+                                       // storing the nonzero
+                                       // elements.
+
+    map.reset (new Epetra_Map(n_rows, 0, communicator));
+
+    Matrix.reset();
+    Matrix = boost::shared_ptr<SparseMatrix> (new SparseMatrix());
+
+    Matrix->reinit (*map, deal_ii_sparse_matrix, drop_tolerance);
+    Matrix->compress();
+
+    initialize (*Matrix, additional_data);
+  }
+  
+  
+  void PreconditionAMG::
+  reinit ()
+  {
+    multilevel_operator->ReComputePreconditioner();
   }
 
 }
index 7f865051f33e4c1d13f57690160ef1d5f1791efb..57ba88ac89cfbd6fc17a3b61de9c23e85d765a87 100644 (file)
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_vector_base.h>
 #include <lac/trilinos_precondition.h>
-#include <lac/trilinos_precondition_amg.h>
 
 #include <cmath>
 
 #ifdef DEAL_II_USE_TRILINOS
 
-#include <Ifpack.h>
-#include <ml_MultiLevelPreconditioner.h>
+#include <Epetra_Operator.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -50,28 +48,6 @@ namespace TrilinosWrappers
 
   
 
-  void
-  SolverBase::solve (const SparseMatrix     &A,
-                     VectorBase             &x,
-                     const VectorBase       &b,
-                     const PreconditionBase &preconditioner)
-  {
-    solve (A, x, b, &*preconditioner.preconditioner);
-  }
-
-
-
-  void
-  SolverBase::solve (const SparseMatrix     &A,
-                     VectorBase             &x,
-                     const VectorBase       &b,
-                     const PreconditionAMG  &preconditioner)
-  {
-    solve (A, x, b, &*preconditioner.multilevel_operator);
-  }
-
-
-
   SolverControl &
   SolverBase::control() const
   {
@@ -84,7 +60,7 @@ namespace TrilinosWrappers
   SolverBase::solve (const SparseMatrix     &A,
                      VectorBase             &x,
                      const VectorBase       &b,
-                     const Epetra_Operator*  preconditioner)
+                     const PreconditionBase &preconditioner)
   {
     int ierr;
     
@@ -129,7 +105,8 @@ namespace TrilinosWrappers
 
                                        // Introduce the
                                        // preconditioner, ...
-    ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>(preconditioner));
+    ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>
+                                    (&*preconditioner.preconditioner));
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
                                        // ... set some options, ...

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.