]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Provide some more Trilinos preconditioners and solvers.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 5 Jan 2009 07:38:24 +0000 (07:38 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 5 Jan 2009 07:38:24 +0000 (07:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@18069 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_block_vector.h
deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/include/lac/trilinos_solver.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_vector_base.h
deal.II/lac/source/trilinos_block_vector.cc
deal.II/lac/source/trilinos_precondition.cc
deal.II/lac/source/trilinos_solver.cc

index 13e27772f52f3e106ccec83bfa4ad590b0bd3ab0..cf313ad4fef86888c9ac59b4d941cf69aa4b1023 100644 (file)
@@ -237,46 +237,41 @@ namespace TrilinosWrappers
                                           */
         void reinit (const unsigned int num_blocks);
 
-                                      /**
-                                       * This reinit function is
-                                       * meant to be used for
-                                       * parallel calculations where
-                                       * some non-local data has to
-                                       * be used. The typical
-                                       * situation where one needs
-                                       * this function is the call of
-                                       * the
-                                       * FEValues<dim>::get_function_values
-                                       * function (or of some
-                                       * derivatives) in
-                                       * parallel. Since it is
-                                       * usually faster to retrieve
-                                       * the data in advance, this
-                                       * function can be called
-                                       * before the assembly forks
-                                       * out to the different
-                                       * processors. What this
-                                       * function does is the
-                                       * following: It takes the
-                                       * information in the columns
-                                       * of the given matrix and
-                                       * looks which data couples
-                                       * between the different
-                                       * processors. That data is
-                                       * then queried from the input
-                                       * vector. Note that you should
-                                       * not write to the resulting
-                                       * vector any more, since the
-                                       * some data can be stored
-                                       * several times on different
-                                       * processors, leading to
-                                       * unpredictable results. In
-                                       * particular, such a vector
-                                       * cannot be used for
-                                       * matrix-vector products as
-                                       * for example done during the
-                                       * solution of linear systems.
-                                       */
+                                        /**
+                                         * This reinit function is meant to
+                                         * be used for parallel
+                                         * calculations where some
+                                         * non-local data has to be
+                                         * used. The typical situation
+                                         * where one needs this function is
+                                         * the call of the
+                                         * FEValues<dim>::get_function_values
+                                         * function (or of some
+                                         * derivatives) in parallel. Since
+                                         * it is usually faster to retrieve
+                                         * the data in advance, this
+                                         * function can be called before
+                                         * the assembly forks out to the
+                                         * different processors. What this
+                                         * function does is the following:
+                                         * It takes the information in the
+                                         * columns of the given matrix and
+                                         * looks which data couples between
+                                         * the different processors. That
+                                         * data is then queried from the
+                                         * input vector. Note that you
+                                         * should not write to the
+                                         * resulting vector any more, since
+                                         * the some data can be stored
+                                         * several times on different
+                                         * processors, leading to
+                                         * unpredictable results. In
+                                         * particular, such a vector cannot
+                                         * be used for matrix-vector
+                                         * products as for example done
+                                         * during the solution of linear
+                                         * systems.
+                                         */
        void import_nonlocal_data_for_fe (const TrilinosWrappers::BlockSparseMatrix &m,
                                          const BlockVector                         &v);
 
index 04101161ad0e187843273480f440c99851667940..dde9daeff02931ae5f974f7df34133389226ac07 100755 (executable)
@@ -882,6 +882,188 @@ namespace TrilinosWrappers
 
 
 
+
+
+
+/**
+ * A wrapper class for a thresholded incomplete LU factorization (ILU-T)
+ * preconditioner for Trilinos matrices. This preconditioner works both in
+ * serial and in parallel, depending on the matrix it is based on. In
+ * general, an incomplete factorization does not take all fill-in elements
+ * that would appear in a full factorization (that is the basis for a direct
+ * solve). For the ILU-T precondtioner, the parameter <tt>ilut_drop</tt>
+ * lets the user specify which elements should be dropped (i.e., should not
+ * be part of the incomplete decomposition). Trilinos calculates first the
+ * complete factorization for one row, and then skips those elements that
+ * are lower than the threshold. This is the main difference to the
+ * non-thresholded ILU preconditioner, where the parameter
+ * <tt>ilut_fill</tt> governs the incomplete factorization structure. This
+ * parameter is available here as well, but provides only some extra
+ * information here.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. Besides the fill-in arguments, these options are some options
+ * for perturbations (see the documentation of the AdditionalData structure
+ * for details), and a parameter <tt>overlap</tt> that determines if and how
+ * much overlap there should be between the matrix partitions on the various
+ * MPI processes. The default settings are 0 for the additional fill-in, 0
+ * for the absolute augmentation tolerance, 1 for the relative augmentation
+ * tolerance, 0 for the overlap.
+ *
+ * Note that a parallel applicatoin of the ILU-T preconditioner is
+ * actually a block-Jacobi preconditioner with block size equal to the
+ * local matrix size. Spoken more technically, this parallel operation
+ * is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an ILU <em>approximate solve</em> as inner
+ * solver, based on the (outer) parallel partitioning.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2009
+ */
+  class PreconditionILUT : public PreconditionBase
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to pipe
+                                        * additional parameters to the
+                                        * preconditioner. The Trilinos ILU-T
+                                        * 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 parameters <tt>ilut_drop</tt>
+                                        * and <tt>ilut_fill</tt>, which
+                                        * specifies a threshold about which
+                                        * values should form the incomplete
+                                        * factorization and the level of
+                                        * additional fill-in. 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>ilut_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, no
+                                       * element will be dropped, the level
+                                       * of extra fill-ins is set to be
+                                       * zero (just use the matrix
+                                       * structure, do not generate any
+                                       * additional fill-in except the one
+                                       * that results from non-dropping
+                                       * large elements), 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       ilut_drop = 0.,
+                       const unsigned int ilut_fill = 0,
+                       const double       ilut_atol = 0.,
+                       const double       ilut_rtol = 1.,
+                       const unsigned int overlap  = 0);
+
+                                       /**
+                                       * This specifies the relative size
+                                       * of elements which should be
+                                       * dropped when forming an incomplete
+                                       * LU decomposition with threshold.
+                                       */
+       double ilut_drop;
+
+                                       /**
+                                       * 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 ilut_fill;
+
+                                       /**
+                                       * 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 ilut_atol;
+
+                                       /**
+                                       * This specifies the factor by
+                                       * which the diagonal of the
+                                       * matrix will be scaled, which
+                                       * sometimes can help to get
+                                       * better preconditioners.
+                                       */
+       double ilut_rtol;
+
+                                       /**
+                                       * This determines how large
+                                       * the overlap of the local
+                                       * matrix portions on each
+                                       * processor in a parallel
+                                       * application should be.
+                                       */
+       unsigned int overlap; 
+      };
+
+                                       /**
+                                        * 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());
+
+    private:
+                                      /**
+                                       * This is a pointer to the
+                                       * Ifpack data contained in
+                                       * this preconditioner.
+                                       */
+      Teuchos::RCP<Ifpack_Preconditioner> ifpack;
+  };
+
+
+
 /**
  * A wrapper class for a sparse direct LU decomposition on parallel
  * blocks for Trilinos matrices. When run in serial, this corresponds
index f7a01c17f16b398d0f0b24ab5cecd2636aa01821..0ce1af8b7d6adaa78868712848c15d0b75b4455e 100644 (file)
 #include <lac/solver_control.h>
 #include <boost/shared_ptr.hpp>
 
+#include <lac/vector.h>
 
 #ifdef DEAL_II_USE_TRILINOS
 
 #include <Epetra_LinearProblem.h>
 #include <AztecOO.h>
 #include <Epetra_Operator.h>
+#include <Amesos.h>
 
 DEAL_II_NAMESPACE_OPEN
 
-
 namespace TrilinosWrappers
 {
                                    // forward declarations
@@ -54,7 +55,7 @@ namespace TrilinosWrappers
  * deal.II class SolverSelector.
  *
  * @ingroup TrilinosWrappers
- * @author Martin Kronbichler, 2008
+ * @author Martin Kronbichler, 2008, 2009
  */
   class SolverBase
   {
@@ -81,18 +82,36 @@ namespace TrilinosWrappers
                                         * pipe additional data to the
                                         * solver.
                                         */
+
       struct AdditionalData
       {
                                       /**
-                                       * Restart parameter in case
-                                       * the derived class is
-                                       * GMRES.
+                                       * Sets the additional data field to
+                                       * the desired output format and puts
+                                       * the restart parameter in case the
+                                       * derived class is GMRES.
                                        * 
                                        * TODO: Find a better way for
-                                       * doing this.
+                                       * setting the GMRES restart
+                                       * parameter since it is quite
+                                       * inelegant to set a specific option
+                                       * of one solver in the base class
+                                       * for all solvers.
+                                       */
+       AdditionalData (const bool         output_solver_details   = false,
+                       const unsigned int gmres_restart_parameter = 30);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
                                        */
-       AdditionalData (const unsigned int gmres_restart_parameter = 30);
+       const bool output_solver_details;
 
+                                      /**
+                                       * Restart parameter for GMRES
+                                       * solver.
+                                       */
        const unsigned int gmres_restart_parameter;
       };
 
@@ -136,6 +155,30 @@ namespace TrilinosWrappers
              const VectorBase       &b,
              const PreconditionBase &preconditioner);
 
+                                       /**
+                                        * Solve the linear system
+                                        * <tt>Ax=b</tt>. Depending on the
+                                        * information provided by derived
+                                        * classes and the object passed as a
+                                        * preconditioner, one of the linear
+                                        * solvers and preconditioners of
+                                        * Trilinos is chosen. This class
+                                        * works with matrices according to
+                                        * the TrilinosWrappers format, but
+                                        * can take deal.II vectors as
+                                        * argument. Since deal.II are serial
+                                        * vectors (not distributed), this
+                                        * function does only what you expect
+                                        * in case the matrix is locally
+                                        * owned. Otherwise, an exception
+                                        * will be thrown.
+                                        */
+      void
+      solve (const SparseMatrix           &A,
+             dealii::Vector<double>       &x,
+             const dealii::Vector<double> &b,
+             const PreconditionBase       &preconditioner);
+
                                        /**
                                         * Access to object that controls
                                         * convergence.
@@ -210,9 +253,23 @@ namespace TrilinosWrappers
                                         * pipe additional data to the
                                         * solver.
                                         */
+
       struct AdditionalData
-      {};
-       
+      {
+                                      /**
+                                       * Sets the additional data field to
+                                       * the desired output format.
+                                       */
+       AdditionalData (const bool output_solver_details = false);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
+                                       */
+       bool output_solver_details;
+      };
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -252,7 +309,20 @@ namespace TrilinosWrappers
                                         * solver.
                                         */
       struct AdditionalData
-      {};
+      {
+                                      /**
+                                       * Sets the additional data field to
+                                       * the desired output format.
+                                       */
+       AdditionalData (const bool output_solver_details = false);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
+                                       */
+       bool output_solver_details;
+      };
        
                                        /**
                                         * Constructor. In contrast to
@@ -299,12 +369,20 @@ namespace TrilinosWrappers
                                             * 30, i.e. do a restart every 30
                                             * iterations.
                                             */
-          AdditionalData (const unsigned int restart_parameter = 30);
-       
-                                           /**
-                                            * Maximum number of
-                                            * tmp vectors.
-                                            */
+         AdditionalData (const bool         output_solver_details = false,
+                         const unsigned int restart_parameter = 30);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
+                                       */
+         bool output_solver_details;
+
+                                       /**
+                                       * Maximum number of
+                                       * tmp vectors.
+                                       */
           unsigned int restart_parameter;
       };
        
@@ -347,7 +425,20 @@ namespace TrilinosWrappers
                                         * solver.
                                         */
       struct AdditionalData
-      {};
+      {
+                                      /**
+                                       * Sets the additional data field to
+                                       * the desired output format.
+                                       */
+       AdditionalData (const bool output_solver_details = false);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
+                                       */
+       bool output_solver_details;
+      };
        
                                        /**
                                         * Constructor. In contrast to
@@ -388,7 +479,20 @@ namespace TrilinosWrappers
                                         * solver.
                                         */
       struct AdditionalData
-      {};
+      {
+                                      /**
+                                       * Sets the additional data field to
+                                       * the desired output format.
+                                       */
+       AdditionalData (const bool output_solver_details = false);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
+                                       */
+       bool output_solver_details;
+      };
        
                                        /**
                                         * Constructor. In contrast to
@@ -411,6 +515,145 @@ namespace TrilinosWrappers
       const AdditionalData additional_data;
   };
 
+
+
+/**
+ * An implementation of the Trilinos KLU direct solver (using the Amesos
+ * package).
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2009
+ */
+  class SolverDirect
+  {
+    public:
+
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional data to the
+                                        * solver.
+                                        */
+
+      struct AdditionalData
+      {
+                                      /**
+                                       * Sets the additional data field to
+                                       * the desired output format.
+                                       */
+       AdditionalData (const bool output_solver_details = false);
+
+                                      /**
+                                       * Enables/disables the output of
+                                       * solver details (residual in each
+                                       * iterations etc.).
+                                       */
+       bool output_solver_details;
+      };
+
+                                       /**
+                                        * Constructor. Takes the
+                                        * solver control object and
+                                        * creates the solver.
+                                        */
+      SolverDirect (SolverControl  &cn,
+                   const AdditionalData &data = AdditionalData());
+      
+                                       /**
+                                        * Destructor.
+                                        */
+      virtual ~SolverDirect ();
+
+                                       /**
+                                        * Solve the linear system
+                                        * <tt>Ax=b</tt>. Creates a KLU
+                                        * factorization of the matrix and
+                                        * performs the solve. Note that
+                                        * there is no need for a
+                                        * preconditioner here.
+                                        */
+      void
+      solve (const SparseMatrix     &A,
+             VectorBase             &x,
+             const VectorBase       &b);
+
+                                       /**
+                                        * Solve the linear system
+                                        * <tt>Ax=b</tt>. Depending on the
+                                        * information provided by derived
+                                        * classes and the object passed as a
+                                        * preconditioner, one of the linear
+                                        * solvers and preconditioners of
+                                        * Trilinos is chosen. This class
+                                        * works with matrices according to
+                                        * the TrilinosWrappers format, but
+                                        * can take deal.II vectors as
+                                        * argument. Since deal.II are serial
+                                        * vectors (not distributed), this
+                                        * function does only what you expect
+                                        * in case the matrix is locally
+                                        * owned. Otherwise, an exception
+                                        * will be thrown.
+                                        */
+      void
+      solve (const SparseMatrix           &A,
+             dealii::Vector<double>       &x,
+             const dealii::Vector<double> &b);
+
+                                       /**
+                                        * Access to object that controls
+                                        * convergence.
+                                        */
+      SolverControl & control() const;
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException1 (ExcTrilinosError,
+                      int,
+                      << "An error with error number " << arg1
+                      << " occured while calling a Trilinos function");
+
+    private:
+
+                                       /**
+                                        * Reference to the object that
+                                        * controls convergence of the
+                                        * iterative solver. In fact,
+                                        * for these Trilinos wrappers,
+                                        * Trilinos does so itself, but
+                                        * we copy the data from this
+                                        * object before starting the
+                                        * solution process, and copy
+                                        * the data back into it
+                                        * afterwards.
+                                        */
+      SolverControl &solver_control;
+
+                                       /**
+                                       * A structure that collects
+                                       * the Trilinos sparse matrix,
+                                       * the right hand side vector
+                                       * and the solution vector,
+                                       * which is passed down to the
+                                       * Trilinos solver.
+                                       */
+      std::auto_ptr<Epetra_LinearProblem> linear_problem;
+
+                                       /**
+                                        * A structure that contains
+                                        * the Trilinos solver and
+                                        * preconditioner objects.
+                                        */
+      std::auto_ptr<Amesos_BaseSolver> solver;
+
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular solver.
+                                        */
+      const AdditionalData additional_data;
+
+  };
+
 }
 
 DEAL_II_NAMESPACE_CLOSE
index fe8f684854762accbc68d224db8d00eaf8cddfcb..de6f2bbd2fd99b0ac4e11b893ed8bcb3aa37fa42 100755 (executable)
@@ -2316,7 +2316,7 @@ namespace TrilinosWrappers
     Assert (values.m() == values.n(), ExcNotQuadratic());
 
     for (unsigned int i=0; i<indices.size(); ++i)
-      add (indices[i], indices.size(), &indices[0], &values(i,0),
+      add (indices[i], indices.size(), &indices[0], &values(i,0), 
           elide_zero_values);
   }
 
@@ -2335,8 +2335,8 @@ namespace TrilinosWrappers
            ExcDimensionMismatch(col_indices.size(), values.n()));
 
     for (unsigned int i=0; i<row_indices.size(); ++i)
-      add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
-          elide_zero_values);
+      add (row_indices[i], col_indices.size(), &col_indices[0], 
+          &values(i,0), elide_zero_values);
   }
 
 
index 1d12d5a8f4785506ee1970725b8a66af5fc30cb4..5fb50bf98104a3a0d19a5cbfb3d4e6cf2d6bc9a2 100644 (file)
@@ -49,6 +49,7 @@ namespace TrilinosWrappers
 {
                                // forward declaration
   class VectorBase;
+  class SparseMatrix;
 
 
                                       /**
@@ -300,12 +301,11 @@ namespace TrilinosWrappers
       void clear ();
 
                                       /**
-                                       * Reinit functionality, sets
-                                       * the dimension and the
-                                       * possibly parallel
-                                       * partitioning (Epetra_Map) of
-                                       * the calling vector the value
-                                       * of the input vector.
+                                       * Reinit functionality, sets the
+                                       * dimension and possibly the
+                                       * parallel partitioning (Epetra_Map)
+                                       * of the calling vector to the
+                                       * settings of the input vector.
                                        */
       void reinit (const VectorBase &v,
                   const bool        fast = false);
index a2d7bae017c4ba0a57589f8bab6552472d5c1477..3ab97b6054d65044aa0a689ad669a0ae7d342762 100644 (file)
@@ -131,7 +131,7 @@ namespace TrilinosWrappers
         this->components.resize(this->n_blocks());
   
       for (unsigned int i=0;i<this->n_blocks();++i)
-        block(i).clear();
+        components[i].clear();
 
       collect_sizes();
     }
index cdf1e1e9146c33900057ca7fe79865a276f35f2a..5df4665db086752a622b8aaff2fda6fad0177878 100755 (executable)
@@ -371,6 +371,60 @@ namespace TrilinosWrappers
 
 
 
+/* -------------------------- PreconditionILUT -------------------------- */
+
+  PreconditionILUT::AdditionalData::
+  AdditionalData (const double       ilut_drop,
+                 const unsigned int ilut_fill,
+                 const double       ilut_atol,
+                 const double       ilut_rtol,
+                 const unsigned int overlap)
+                  :
+                  ilut_drop (ilut_drop),
+                  ilut_fill (ilut_fill),
+                 ilut_atol (ilut_atol),
+                 ilut_rtol (ilut_rtol),
+                 overlap  (overlap)
+  {}
+
+
+
+  void
+  PreconditionILUT::initialize (const SparseMatrix   &matrix,
+                               const AdditionalData &additional_data)
+  {
+    preconditioner.release();
+    ifpack.release();
+
+    ifpack = Teuchos::rcp (Ifpack().Create ("ILUT", &*matrix.matrix, 
+                                           additional_data.overlap));
+
+    Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                      "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("fact: drop value",additional_data.ilut_drop);  
+    parameter_list.set ("fact: level-of-fill",(int)additional_data.ilut_fill);  
+    parameter_list.set ("fact: absolute threshold",additional_data.ilut_atol); 
+    parameter_list.set ("fact: relative threshold",additional_data.ilut_rtol); 
+    parameter_list.set ("schwarz: combine mode", "Add");
+    
+    ierr = ifpack->SetParameters(parameter_list);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = ifpack->Initialize();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = ifpack->Compute();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    preconditioner = Teuchos::rcp (&*ifpack, false);
+  }
+
+
+
 /* ---------------------- PreconditionBlockDirect --------------------- */
 
   PreconditionBlockwiseDirect::AdditionalData::
index 4ca0f6ed8c88cfb81c16f5cf1e89dd2f6ce0872c..d22462b4b7a380ed37f19bf664da5d99480bdda7 100644 (file)
@@ -12,6 +12,7 @@
 //---------------------------------------------------------------------------
 
 
+#include <base/conditional_ostream.h>
 #include <lac/trilinos_solver.h>
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_vector_base.h>
 
 #ifdef DEAL_II_USE_TRILINOS
 
-#include <Epetra_Operator.h>
-
 DEAL_II_NAMESPACE_OPEN
 
 namespace TrilinosWrappers
 {
 
-  SolverBase::AdditionalData::AdditionalData (const unsigned int gmres_restart_parameter)
+  SolverBase::AdditionalData::AdditionalData (const bool         output_solver_details,
+                                             const unsigned int gmres_restart_parameter)
                   :
+                 output_solver_details (output_solver_details),
                   gmres_restart_parameter (gmres_restart_parameter)
   {}
 
-  
-  
+
+
   SolverBase::SolverBase (SolverControl  &cn)
                   :
                  solver_name    (gmres),
                   solver_control (cn)
   {}
 
-  
+
 
   SolverBase::SolverBase (const enum SolverBase::SolverName  solver_name,
                          SolverControl                     &cn)
@@ -50,12 +51,12 @@ namespace TrilinosWrappers
                   solver_control (cn)
   {}
 
-  
+
 
   SolverBase::~SolverBase ()
   {}
 
-  
+
 
   SolverControl &
   SolverBase::control() const
@@ -119,7 +120,101 @@ namespace TrilinosWrappers
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
                                        // ... set some options, ...
-    solver.SetAztecOption (AZ_output, AZ_none);
+    solver.SetAztecOption (AZ_output, additional_data.output_solver_details ? 
+                                     AZ_all : AZ_none);
+    solver.SetAztecOption (AZ_conv, AZ_noscaled);
+
+                                       // ... and then solve!
+    ierr = solver.Iterate (solver_control.max_steps(), 
+                          solver_control.tolerance());
+    AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+
+                                       // Finally, let the deal.II
+                                       // SolverControl object know
+                                       // what has happened. If the
+                                       // solve succeeded, the status
+                                       // of the solver control will
+                                       // turn into
+                                       // SolverControl::success.
+    solver_control.check (solver.NumIters(), solver.TrueResidual());
+
+    if (solver_control.last_check() != SolverControl::success)
+      throw SolverControl::NoConvergence (solver_control.last_step(),
+                                          solver_control.last_value());
+  }
+
+
+
+  void
+  SolverBase::solve (const SparseMatrix           &A,
+                     dealii::Vector<double>       &x,
+                     const dealii::Vector<double> &b,
+                     const PreconditionBase       &preconditioner)
+  {
+    int ierr;
+    
+    linear_problem.reset();
+
+                                       // In case we call the solver with
+                                       // deal.II vectors, we create views
+                                       // of the vectors in Epetra format.
+    Assert (x.size() == A.n(),
+           ExcDimensionMismatch(x.size(), A.n()));
+    Assert (b.size() == A.m(),
+           ExcDimensionMismatch(b.size(), A.m()));
+    Assert (A.local_range ().second == A.m(),
+           ExcMessage ("Can only work in serial when using deal.II vectors."));
+    Assert (A.matrix->Filled(),
+           ExcMessage ("Matrix is not compressed. Call compress() method."));
+
+    Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin());
+    Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast<double*>(b.begin()));
+
+                                       // We need an
+                                       // Epetra_LinearProblem object
+                                       // to let the AztecOO solver
+                                       // know about the matrix and
+                                       // vectors.
+    linear_problem = std::auto_ptr<Epetra_LinearProblem> 
+      (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b));
+
+                                       // Next we can allocate the
+                                       // AztecOO solver...
+    solver.SetProblem(*linear_problem);
+
+                                       // ... and we can specify the
+                                       // solver to be used.
+    switch (solver_name)
+      {
+      case cg:
+       solver.SetAztecOption(AZ_solver, AZ_cg);
+       break;
+      case cgs:
+       solver.SetAztecOption(AZ_solver, AZ_cgs);
+       break;
+      case gmres:
+       solver.SetAztecOption(AZ_solver, AZ_gmres);
+       solver.SetAztecOption(AZ_kspace, additional_data.gmres_restart_parameter);
+       break;
+      case bicgstab:
+       solver.SetAztecOption(AZ_solver, AZ_bicgstab);
+       break;
+      case tfqmr:
+       solver.SetAztecOption(AZ_solver, AZ_tfqmr);
+       break;
+      default:
+       Assert (false, ExcNotImplemented());
+      }
+
+                                       // Introduce the
+                                       // preconditioner, ...
+    ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>
+                                    (&*preconditioner.preconditioner));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+                                       // ... set some options, ...
+    solver.SetAztecOption (AZ_output, additional_data.output_solver_details ? 
+                                     AZ_all : AZ_none);
     solver.SetAztecOption (AZ_conv, AZ_noscaled);
 
                                        // ... and then solve!
@@ -147,11 +242,19 @@ namespace TrilinosWrappers
   
 /* ---------------------- SolverCG ------------------------ */
 
+  SolverCG::AdditionalData::
+  AdditionalData (const bool output_solver_details)
+                  :
+                  output_solver_details (output_solver_details)
+  {}
+
+
+
   SolverCG::SolverCG (SolverControl        &cn,
                       const AdditionalData &data)
                   :
                   SolverBase (cn),
-                  additional_data (data)
+                  additional_data (data.output_solver_details)
   {
     solver_name = cg;
   }
@@ -160,8 +263,10 @@ namespace TrilinosWrappers
 /* ---------------------- SolverGMRES ------------------------ */
 
   SolverGMRES::AdditionalData::
-  AdditionalData (const unsigned int restart_parameter)
+    AdditionalData (const bool output_solver_details,
+                   const unsigned int restart_parameter)
                   :
+                  output_solver_details (output_solver_details),
                   restart_parameter (restart_parameter)
   {}
 
@@ -171,7 +276,8 @@ namespace TrilinosWrappers
                             const AdditionalData &data)
                   :
                   SolverBase (cn),
-                  additional_data (data.restart_parameter)
+                  additional_data (data.output_solver_details,
+                                  data.restart_parameter)
   {
     solver_name = gmres;
   }
@@ -179,11 +285,20 @@ namespace TrilinosWrappers
 
 /* ---------------------- SolverBicgstab ------------------------ */
 
+  SolverBicgstab::AdditionalData::
+  AdditionalData (const bool output_solver_details)
+                  :
+                  output_solver_details (output_solver_details)
+  {}
+
+
+
+
   SolverBicgstab::SolverBicgstab (SolverControl        &cn,
                                   const AdditionalData &data)
                   :
                   SolverBase (cn),
-                  additional_data (data)
+                  additional_data (data.output_solver_details)
   {
     solver_name = bicgstab;
   }
@@ -191,11 +306,20 @@ namespace TrilinosWrappers
   
 /* ---------------------- SolverCGS ------------------------ */
 
+  SolverCGS::AdditionalData::
+  AdditionalData (const bool output_solver_details)
+                  :
+                  output_solver_details (output_solver_details)
+  {}
+
+
+
+
   SolverCGS::SolverCGS (SolverControl        &cn,
                         const AdditionalData &data)
                   :
                   SolverBase (cn),
-                  additional_data (data)
+                  additional_data (data.output_solver_details)
   {
     solver_name = cgs;
   }
@@ -203,15 +327,177 @@ namespace TrilinosWrappers
 
 /* ---------------------- SolverTFQMR ------------------------ */
 
+  SolverTFQMR::AdditionalData::
+  AdditionalData (const bool output_solver_details)
+                  :
+                  output_solver_details (output_solver_details)
+  {}
+
+
+
   SolverTFQMR::SolverTFQMR (SolverControl        &cn,
                             const AdditionalData &data)
                   :
                   SolverBase (cn),
-                  additional_data (data)
+                  additional_data (data.output_solver_details)
   {
     solver_name = tfqmr;
   }
 
+
+
+/* ---------------------- SolverDirect ------------------------ */
+
+  SolverDirect::AdditionalData::
+  AdditionalData (const bool output_solver_details)
+                  :
+                  output_solver_details (output_solver_details)
+  {}
+
+
+
+
+  SolverDirect::SolverDirect (SolverControl  &cn,
+                             const AdditionalData &data)
+                  :
+                  solver_control (cn),
+                  additional_data (data.output_solver_details)
+  {}
+
+  
+
+  SolverDirect::~SolverDirect ()
+  {}
+
+  
+
+  SolverControl &
+  SolverDirect::control() const
+  {
+    return solver_control;
+  }
+
+
+
+  void
+  SolverDirect::solve (const SparseMatrix     &A,
+                      VectorBase             &x,
+                      const VectorBase       &b)
+  {
+                                      // First set whether we want to print
+                                      // the solver information to screen
+                                      // or not.
+    ConditionalOStream  verbose_cout (std::cout, 
+                                     additional_data.output_solver_details);
+
+    linear_problem.reset();
+    solver.reset();
+
+                                       // We need an
+                                       // Epetra_LinearProblem object
+                                       // to let the AztecOO solver
+                                       // know about the matrix and
+                                       // vectors.
+    linear_problem = std::auto_ptr<Epetra_LinearProblem> (
+                             new Epetra_LinearProblem(&*(A.matrix), &*x.vector,
+                                                      &*b.vector));
+
+                                       // Next we can allocate the
+                                       // AztecOO solver...
+    solver = std::auto_ptr<Amesos_BaseSolver> (Amesos().Create("Amesos_Klu",
+                                                              *linear_problem));
+
+    verbose_cout << "Starting symbolic factorization" << std::endl;
+    solver->SymbolicFactorization();
+
+    verbose_cout << "Starting numeric factorization" << std::endl;
+    solver->NumericFactorization();
+
+    verbose_cout << "Starting solve" << std::endl;
+    solver->Solve();
+
+                                       // Finally, let the deal.II
+                                       // SolverControl object know
+                                       // what has happened. If the
+                                       // solve succeeded, the status
+                                       // of the solver control will
+                                       // turn into
+                                       // SolverControl::success.
+    solver_control.check (0, 0);
+
+    if (solver_control.last_check() != SolverControl::success)
+      throw SolverControl::NoConvergence (solver_control.last_step(),
+                                          solver_control.last_value());
+  }
+
+
+
+  void
+  SolverDirect::solve (const SparseMatrix           &A,
+                      dealii::Vector<double>       &x,
+                      const dealii::Vector<double> &b)
+  {
+                                      // First set whether we want to print
+                                      // the solver information to screen
+                                      // or not.
+    ConditionalOStream  verbose_cout (std::cout, 
+                                     additional_data.output_solver_details);
+
+    linear_problem.reset();
+    solver.reset();
+
+
+                                       // In case we call the solver with
+                                       // deal.II vectors, we create views
+                                       // of the vectors in Epetra format.
+    Assert (x.size() == A.n(),
+           ExcDimensionMismatch(x.size(), A.n()));
+    Assert (b.size() == A.m(),
+           ExcDimensionMismatch(b.size(), A.m()));
+    Assert (A.local_range ().second == A.m(),
+           ExcMessage ("Can only work in serial when using deal.II vectors."));
+    Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin());
+    Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast<double*>(b.begin()));
+
+                                       // We need an
+                                       // Epetra_LinearProblem object
+                                       // to let the AztecOO solver
+                                       // know about the matrix and
+                                       // vectors.
+    linear_problem = std::auto_ptr<Epetra_LinearProblem> 
+      (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b));
+
+                                       // Next we can allocate the
+                                       // AztecOO solver...
+    solver = std::auto_ptr<Amesos_BaseSolver> (Amesos().Create("Amesos_Klu",
+                                                              *linear_problem));
+
+    verbose_cout << "Starting symbolic factorization" << std::endl;
+    solver->SymbolicFactorization();
+
+    verbose_cout << "Starting numeric factorization" << std::endl;
+    solver->NumericFactorization();
+
+    verbose_cout << "Starting solve" << std::endl;
+    solver->Solve();
+
+                                       // Finally, let the deal.II
+                                       // SolverControl object know
+                                       // what has happened. If the
+                                       // solve succeeded, the status
+                                       // of the solver control will
+                                       // turn into
+                                       // SolverControl::success.
+    solver_control.check (0, 0);
+
+    if (solver_control.last_check() != SolverControl::success)
+      throw SolverControl::NoConvergence (solver_control.last_step(),
+                                          solver_control.last_value());
+  }
+  
+    
+
+
 }
 
 DEAL_II_NAMESPACE_CLOSE

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.