]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement some new Trilinos preconditioner. Tests.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 4 Jul 2014 13:46:44 +0000 (13:46 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 4 Jul 2014 13:46:44 +0000 (13:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@33109 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/trilinos_precondition.h
deal.II/source/lac/trilinos_precondition.cc
tests/trilinos/precondition.cc [new file with mode: 0644]
tests/trilinos/precondition.output [new file with mode: 0644]

index 99ee610c27f95f4ca9ced743b40427d95b2586fc..974eef62d323188401c808bec1949847fa722075 100644 (file)
@@ -148,6 +148,17 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: There are now three new preconditioner classes
+  TrilinosWrappers::PreconditionBlockJacobi,
+  TrilinosWrappers::PreconditionBlockSSOR, and
+  TrilinosWrappers::PreconditionBlockSOR that work on small dense blocks of
+  the global matrix instead of the point-wise relaxation methods in
+  TrilinosWrappers::Precondition{Jacobi,SSOR,SOR} that work on each row
+  separately.
+  <br>
+  (Martin Kronbichler, 2014/07/04)
+  </li>
+
   <li> Fixed: Some versions of DoFTools::extract_boundary_dofs() were
   not instantiated for some combinations of arguments. This could lead
   to missing symbol errors during linking of applications on some
index 9ab30fb98f92292488c6fc5661febcb3f4ca03c7..2d1b86302e15b27cbc2b010f80f1d4200389c709 100644 (file)
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 2008 - 2013 by the deal.II authors
+// Copyright (C) 2008 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -83,20 +83,15 @@ namespace TrilinosWrappers
     typedef dealii::types::global_dof_index size_type;
 
     /**
-     * Standardized data struct to
-     * pipe additional flags to the
+     * Standardized data struct to pipe additional flags to the
      * preconditioner.
      */
     struct AdditionalData
     {};
 
     /**
-     * Constructor. Does not do
-     * anything. The
-     * <tt>initialize</tt> function
-     * of the derived classes will
-     * have to create the
-     * preconditioner from a given
+     * 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 ();
@@ -112,8 +107,7 @@ namespace TrilinosWrappers
     ~PreconditionBase ();
 
     /**
-     * Destroys the preconditioner, leaving
-     * an object like just after having
+     * Destroys the preconditioner, leaving an object like just after having
      * called the constructor.
      */
     void clear ();
@@ -131,39 +125,29 @@ namespace TrilinosWrappers
                          const VectorBase &src) const;
 
     /**
-     * Apply the preconditioner on
-     * deal.II data structures
-     * instead of the ones provided
-     * in the Trilinos wrapper
-     * class.
+     * Apply the preconditioner on deal.II data structures instead of the ones
+     * provided in the Trilinos wrapper class.
      */
     virtual void vmult (dealii::Vector<double>       &dst,
                         const dealii::Vector<double> &src) const;
 
     /**
-     * Apply the transpose preconditioner on
-     * deal.II data structures
-     * instead of the ones provided
-     * in the Trilinos wrapper
-     * class.
+     * Apply the transpose preconditioner on deal.II data structures instead
+     * of the ones provided in the Trilinos wrapper class.
      */
     virtual void Tvmult (dealii::Vector<double>       &dst,
                          const dealii::Vector<double> &src) const;
 
     /**
-     * Apply the preconditioner on deal.II
-     * parallel data structures instead of
-     * the ones provided in the Trilinos
-     * wrapper class.
+     * Apply the preconditioner on deal.II parallel data structures instead of
+     * the ones provided in the Trilinos wrapper class.
      */
     virtual void vmult (dealii::parallel::distributed::Vector<double>       &dst,
                         const dealii::parallel::distributed::Vector<double> &src) const;
 
     /**
-     * Apply the transpose preconditioner on deal.II
-     * parallel data structures instead of
-     * the ones provided in the Trilinos
-     * wrapper class.
+     * Apply the transpose preconditioner on deal.II parallel data structures
+     * instead of the ones provided in the Trilinos wrapper class.
      */
     virtual void Tvmult (dealii::parallel::distributed::Vector<double>       &dst,
                          const dealii::parallel::distributed::Vector<double> &src) const;
@@ -183,18 +167,14 @@ namespace TrilinosWrappers
 
   protected:
     /**
-     * This is a pointer to the
-     * preconditioner object that
-     * is used when applying the
-     * preconditioner.
+     * This is a pointer to the preconditioner object that is used when
+     * applying the preconditioner.
      */
     std_cxx1x::shared_ptr<Epetra_Operator> preconditioner;
 
     /**
-     * Internal communication
-     * pattern in case the matrix
-     * needs to be copied from
-     * deal.II format.
+     * Internal communication pattern in case the matrix needs to be copied
+     * from deal.II format.
      */
 #ifdef DEAL_II_WITH_MPI
     Epetra_MpiComm     communicator;
@@ -203,9 +183,8 @@ namespace TrilinosWrappers
 #endif
 
     /**
-     * Internal Trilinos map in
-     * case the matrix needs to be
-     * copied from deal.II format.
+     * Internal Trilinos map in case the matrix needs to be copied from
+     * deal.II format.
      */
     std_cxx1x::shared_ptr<Epetra_Map>   vector_distributor;
   };
@@ -233,69 +212,50 @@ namespace TrilinosWrappers
   public:
 
     /**
-     * Standardized data struct to
-     * pipe additional flags to the
-     * preconditioner. The
-     * parameter <tt>omega</tt>
-     * specifies the relaxation
-     * parameter in the Jacobi
-     * preconditioner. The
-     * parameter
-     * <tt>min_diagonal</tt> can be
-     * used to make the application
-     * of the preconditioner also
-     * possible when some diagonal
-     * elements are zero. In a
-     * default application this
-     * would mean that we divide by
-     * zero, so by setting the
-     * parameter
-     * <tt>min_diagonal</tt> to a
-     * small nonzero value the SOR
-     * will work on a matrix that
-     * is not too far away from the
-     * one we want to
-     * treat.
+     * Standardized data struct to pipe additional flags to the
+     * preconditioner. The parameter <tt>omega</tt> specifies the relaxation
+     * parameter in the Jacobi preconditioner. The parameter
+     * <tt>min_diagonal</tt> can be used to make the application of the
+     * preconditioner also possible when some diagonal elements are zero. In a
+     * default application this would mean that we divide by zero, so by
+     * setting the parameter <tt>min_diagonal</tt> to a small nonzero value
+     * the SOR will work on a matrix that is not too far away from the one we
+     * want to treat.
      */
     struct AdditionalData
     {
       /**
-       * Constructor. By default, set
-       * the damping parameter to
-       * one, and do not modify the
-       * diagonal.
+       * Constructor. By default, set the damping parameter to one, and do not
+       * modify the diagonal.
        */
       AdditionalData (const double       omega = 1,
-                      const double       min_diagonal = 0);
+                      const double       min_diagonal = 0,
+                      const unsigned int n_sweeps = 1);
 
       /**
-       * This specifies the
-       * relaxation parameter in the
-       * Jacobi preconditioner.
+       * This specifies the relaxation parameter in the Jacobi preconditioner.
        */
       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.
+       * 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;
+
+      /**
+       * Sets how many times the given operation should be applied during the
+       * vmult() operation.
+       */
+      unsigned int n_sweeps;
     };
 
     /**
-     * Take the sparse matrix the
-     * preconditioner object should
-     * be built of, and additional
-     * flags (damping parameter,
-     * etc.) if there are any.
+     * Take the sparse matrix the preconditioner object should be built of,
+     * and additional flags (damping parameter, etc.) if there are any.
      */
     void initialize (const SparseMatrix   &matrix,
                      const AdditionalData &additional_data = AdditionalData());
@@ -337,93 +297,63 @@ namespace TrilinosWrappers
   public:
 
     /**
-     * Standardized data struct to
-     * pipe additional flags to the
-     * preconditioner. The
-     * parameter <tt>omega</tt>
-     * specifies the relaxation
-     * parameter in the SSOR
-     * preconditioner. The
-     * parameter
-     * <tt>min_diagonal</tt> can be
-     * used to make the application
-     * of the preconditioner also
-     * possible when some diagonal
-     * elements are zero. In a
-     * default application this
-     * would mean that we divide by
-     * zero, so by setting the
-     * parameter
-     * <tt>min_diagonal</tt> to a
-     * small nonzero value the SOR
-     * will work on a matrix that
-     * is not too far away from the
-     * one we want to
-     * treat. Finally,
-     * <tt>overlap</tt> governs the
-     * overlap of the partitions
-     * when the preconditioner runs
-     * in parallel, forming a
-     * so-called additive Schwarz
-     * preconditioner.
+     * Standardized data struct to pipe additional flags to the
+     * preconditioner. The parameter <tt>omega</tt> specifies the relaxation
+     * parameter in the SSOR preconditioner. The parameter
+     * <tt>min_diagonal</tt> can be used to make the application of the
+     * preconditioner also possible when some diagonal elements are zero. In a
+     * default application this would mean that we divide by zero, so by
+     * setting the parameter <tt>min_diagonal</tt> to a small nonzero value
+     * the SOR will work on a matrix that is not too far away from the one we
+     * want to treat. Finally, <tt>overlap</tt> governs the overlap of the
+     * partitions when the preconditioner runs in parallel, forming a
+     * so-called additive Schwarz preconditioner.
      */
     struct AdditionalData
     {
       /**
-       * Constructor. By default, set
-       * the damping parameter to
-       * one, we do not modify the
-       * diagonal, and there is no
-       * overlap (i.e. in parallel,
-       * we run a BlockJacobi
-       * preconditioner, where each
-       * block is inverted
-       * approximately by an SSOR.
+       * Constructor. By default, set the damping parameter to one, we do not
+       * modify the diagonal, and there is no overlap (i.e. in parallel, we
+       * run a BlockJacobi preconditioner, where each block is inverted
+       * approximately by an SSOR).
        */
       AdditionalData (const double       omega = 1,
                       const double       min_diagonal = 0,
-                      const unsigned int overlap = 0);
+                      const unsigned int overlap = 0,
+                      const unsigned int n_sweeps = 1);
 
       /**
-       * This specifies the (over-)
-       * relaxation parameter in the
-       * SSOR preconditioner.
+       * 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.
+       * 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.
+       * This determines how large the overlap of the local matrix portions on
+       * each processor in a parallel application should be.
        */
       unsigned int overlap;
+
+      /**
+       * Sets how many times the given operation should be applied during the
+       * vmult() operation.
+       */
+      unsigned int n_sweeps;
     };
 
     /**
-     * Take the sparse matrix the
-     * preconditioner object should
-     * be built of, and additional
-     * flags (damping parameter,
-     * overlap in parallel
-     * computations, etc.) if there
-     * are any.
+     * Take the sparse matrix the preconditioner object should be built of,
+     * and additional flags (damping parameter, overlap in parallel
+     * computations, etc.) if there are any.
      */
     void initialize (const SparseMatrix   &matrix,
                      const AdditionalData &additional_data = AdditionalData());
@@ -465,93 +395,390 @@ namespace TrilinosWrappers
   public:
 
     /**
-     * Standardized data struct to
-     * pipe additional flags to the
-     * preconditioner. The
-     * parameter <tt>omega</tt>
-     * specifies the relaxation
-     * parameter in the SOR
-     * preconditioner. The
-     * parameter
-     * <tt>min_diagonal</tt> can be
-     * used to make the application
-     * of the preconditioner also
-     * possible when some diagonal
-     * elements are zero. In a
-     * default application this
-     * would mean that we divide by
-     * zero, so by setting the
-     * parameter
-     * <tt>min_diagonal</tt> to a
-     * small nonzero value the SOR
-     * will work on a matrix that
-     * is not too far away from the
-     * one we want to
-     * treat. Finally,
-     * <tt>overlap</tt> governs the
-     * overlap of the partitions
-     * when the preconditioner runs
-     * in parallel, forming a
-     * so-called additive Schwarz
-     * preconditioner.
+     * Standardized data struct to pipe additional flags to the
+     * preconditioner. The parameter <tt>omega</tt> specifies the relaxation
+     * parameter in the SOR preconditioner. The parameter
+     * <tt>min_diagonal</tt> can be used to make the application of the
+     * preconditioner also possible when some diagonal elements are zero. In a
+     * default application this would mean that we divide by zero, so by
+     * setting the parameter <tt>min_diagonal</tt> to a small nonzero value
+     * the SOR will work on a matrix that is not too far away from the one we
+     * want to treat. Finally, <tt>overlap</tt> governs the overlap of the
+     * partitions when the preconditioner runs in parallel, forming a
+     * so-called additive Schwarz preconditioner.
      */
     struct AdditionalData
     {
       /**
-       * Constructor. By default, set
-       * the damping parameter to
-       * one, we do not modify the
-       * diagonal, and there is no
-       * overlap (i.e. in parallel,
-       * we run a BlockJacobi
-       * preconditioner, where each
-       * block is inverted
+       * Constructor. By default, set the damping parameter to one, we do not
+       * modify the diagonal, and there is no overlap (i.e. in parallel, we
+       * run a BlockJacobi preconditioner, where each block is inverted
        * approximately by an SOR.
        */
       AdditionalData (const double       omega = 1,
                       const double       min_diagonal = 0,
-                      const unsigned int overlap = 0);
+                      const unsigned int overlap = 0,
+                      const unsigned int n_sweeps = 1);
+
+      /**
+       * 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.
+       */
+      unsigned int overlap;
+
+      /**
+       * Sets how many times the given operation should be applied during the
+       * vmult() operation.
+       */
+      unsigned int n_sweeps;
+    };
+
+    /**
+     * Take the sparse matrix the preconditioner object should be built of,
+     * and additional flags (damping parameter, overlap in parallel
+     * computations etc.) if there are any.
+     */
+    void initialize (const SparseMatrix   &matrix,
+                     const AdditionalData &additional_data = AdditionalData());
+  };
+
+
+
+  /**
+   * A wrapper class for a block Jacobi preconditioner for Trilinos
+   * matrices. As opposed to PreconditionSOR where each row is treated
+   * separately, this scheme collects block of a given size and inverts a full
+   * matrix for all these rows simultaneously. Trilinos allows to select
+   * several strategies for selecting which rows form a block, including
+   * "linear" (i.e., divide the local range of the matrix in slices of the
+   * block size), "greedy" or "metis". Note that the term <em>block
+   * Jacobi</em> does not relate to possible blocks in the MPI setting, but
+   * small blocks of dense matrices extracted from the sparse matrix local to
+   * each processor.
+   *
+   * The AdditionalData data structure allows to set preconditioner
+   * options.
+   *
+   * @ingroup TrilinosWrappers
+   * @ingroup Preconditioners
+   * @author Martin Kronbichler, 2014
+   */
+  class PreconditionBlockJacobi : public PreconditionBase
+  {
+  public:
+
+    /**
+     * Standardized data struct to pipe additional flags to the
+     * preconditioner. The parameter <tt>block_size</tt> sets the size of
+     * small blocks. It is recommended to choose this parameter not too large
+     * (a few hundreds at most) since this implementation uses a dense matrix
+     * for the block. The parameter <tt>block_creation_type</tt> allows to
+     * pass the strategy for finding the blocks to Ifpack. The parameter
+     * <tt>omega</tt> specifies the relaxation parameter in the SOR
+     * preconditioner. The parameter <tt>min_diagonal</tt> can be used to make
+     * the application of the preconditioner also possible when some diagonal
+     * elements are zero. In a default application this would mean that we
+     * divide by zero, so by setting the parameter <tt>min_diagonal</tt> to a
+     * small nonzero value the SOR will work on a matrix that is not too far
+     * away from the one we want to treat. Finally, <tt>overlap</tt> governs
+     * the overlap of the partitions when the preconditioner runs in parallel,
+     * forming a so-called additive Schwarz preconditioner.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor. By default, use a block size of 1, use linear
+       * subdivision of the rows, set the damping parameter to one, and do not
+       * modify the diagonal.
+       */
+      AdditionalData (const unsigned int block_size = 1,
+                      const std::string  block_creation_type = "linear",
+                      const double       omega = 1,
+                      const double       min_diagonal = 0,
+                      const unsigned int n_sweeps = 1);
+
+      /**
+       * This specifies the size of blocks.
+       */
+      unsigned int block_size;
+
+      /**
+       * Strategy for creation of blocks passed on to Ifpack block relaxation
+       * (variable 'partitioner: type') with this string as the given
+       * value. Available types in Ifpack include "linear" (i.e., divide the
+       * local range of the matrix in slices of the block size), "greedy"
+       * "metis". For a full list, see the documentation of Ifpack.
+       */
+      std::string block_creation_type;
+
+      /**
+       * This specifies the (over-) relaxation parameter in the Jacobi
+       * preconditioner.
+       */
+      double omega;
+
+      /**
+       * This specifies the minimum value the diagonal elements should
+       * have. This might be necessary when the block Jacobi 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;
+
+      /**
+       * Sets how many times the given operation should be applied during the
+       * vmult() operation.
+       */
+      unsigned int n_sweeps;
+    };
+
+    /**
+     * Take the sparse matrix the preconditioner object should be built of,
+     * and additional flags (damping parameter, etc.) if there are any.
+     */
+    void initialize (const SparseMatrix   &matrix,
+                     const AdditionalData &additional_data = AdditionalData());
+  };
+
+
+
+
+  /**
+   * A wrapper class for a (pointwise) SSOR preconditioner for Trilinos
+   * matrices. As opposed to PreconditionSSOR where each row is treated
+   * separately, this scheme collects block of a given size and inverts a full
+   * matrix for all these rows simultaneously. Trilinos allows to select
+   * several strategies for selecting which rows form a block, including
+   * "linear" (i.e., divide the local range of the matrix in slices of the
+   * block size), "greedy" or "metis".
+   *
+   * The AdditionalData data structure allows to set preconditioner
+   * options.
+   *
+   * Note that a parallel application of this preconditioner is actually a
+   * block-Jacobi preconditioner with (outer) 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 a block SSOR <em>approximate solve</em> as inner
+   * solver, based on the outer parallel partitioning.
+   *
+   * @ingroup TrilinosWrappers
+   * @ingroup Preconditioners
+   * @author Martin Kronbichler, 2014
+   */
+  class PreconditionBlockSSOR : public PreconditionBase
+  {
+  public:
+
+    /**
+     * Standardized data struct to pipe additional flags to the
+     * preconditioner. The parameter <tt>block_size</tt> sets the size of
+     * small blocks. It is recommended to choose this parameter not too large
+     * (a few hundreds at most) since this implementation uses a dense matrix
+     * for the block. The parameter <tt>block_creation_type</tt> allows to
+     * pass the strategy for finding the blocks to Ifpack. The parameter
+     * <tt>omega</tt> specifies the relaxation parameter in the SSOR
+     * preconditioner. The parameter <tt>min_diagonal</tt> can be used to make
+     * the application of the preconditioner also possible when some diagonal
+     * elements are zero. In a default application this would mean that we
+     * divide by zero, so by setting the parameter <tt>min_diagonal</tt> to a
+     * small nonzero value the SOR will work on a matrix that is not too far
+     * away from the one we want to treat. Finally, <tt>overlap</tt> governs
+     * the overlap of the partitions when the preconditioner runs in parallel,
+     * forming a so-called additive Schwarz preconditioner.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor. By default, use a block size of 1, use linear
+       * subdivision of the rows, set the damping parameter to one, we do not
+       * modify the diagonal, and there is no overlap (i.e. in parallel, we
+       * run a BlockJacobi preconditioner, where each block is inverted
+       * approximately by a block SOR).
+       */
+      AdditionalData (const unsigned int block_size = 1,
+                      const std::string  block_creation_type = "linear",
+                      const double       omega = 1,
+                      const double       min_diagonal = 0,
+                      const unsigned int overlap = 0,
+                      const unsigned int n_sweeps = 1);
+
+      /**
+       * This specifies the size of blocks.
+       */
+      unsigned int block_size;
+
+      /**
+       * Strategy for creation of blocks passed on to Ifpack block relaxation
+       * (variable 'partitioner: type') with this string as the given
+       * value. Available types in Ifpack include "linear" (i.e., divide the
+       * local range of the matrix in slices of the block size), "greedy"
+       * "metis". For a full list, see the documentation of Ifpack.
+       */
+      std::string block_creation_type;
 
       /**
-       * This specifies the (over-)
-       * relaxation parameter in the
-       * SOR preconditioner.
+       * 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.
+       * 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.
+       * This determines how large the overlap of the local matrix portions on
+       * each processor in a parallel application should be.
        */
       unsigned int overlap;
+
+      /**
+       * Sets how many times the given operation should be applied during the
+       * vmult() operation.
+       */
+      unsigned int n_sweeps;
     };
 
     /**
-     * Take the sparse matrix the
-     * preconditioner object should
-     * be built of, and additional
-     * flags (damping parameter,
-     * overlap in parallel
-     * computations etc.) if there
-     * are any.
+     * Take the sparse matrix the preconditioner object should be built of,
+     * and additional flags (damping parameter, overlap in parallel
+     * computations, etc.) if there are any.
+     */
+    void initialize (const SparseMatrix   &matrix,
+                     const AdditionalData &additional_data = AdditionalData());
+  };
+
+
+
+
+  /**
+   * A wrapper class for a block SOR preconditioner for Trilinos matrices. As
+   * opposed to PreconditionSOR where each row is treated separately, this
+   * scheme collects block of a given size and inverts a full matrix for all
+   * these rows simultaneously. Trilinos allows to select several strategies
+   * for selecting which rows form a block, including "linear" (i.e., divide
+   * the local range of the matrix in slices of the block size), "greedy" or
+   * "metis".
+   *
+   * The AdditionalData data structure allows to set preconditioner
+   * options.
+   *
+   * Note that a parallel application of this preconditioner is actually a
+   * block-Jacobi preconditioner with (outer) 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 a block SOR <em>approximate solve</em> as inner
+   * solver, based on the outer parallel partitioning.
+   *
+   * @ingroup TrilinosWrappers
+   * @ingroup Preconditioners
+   * @author Martin Kronbichler, 2014
+   */
+  class PreconditionBlockSOR : public PreconditionBase
+  {
+  public:
+
+    /**
+     * Standardized data struct to pipe additional flags to the
+     * preconditioner. The parameter <tt>block_size</tt> sets the size of
+     * small blocks. It is recommended to choose this parameter not too large
+     * (a few hundreds at most) since this implementation uses a dense matrix
+     * for the block. The parameter <tt>block_creation_type</tt> allows to
+     * pass the strategy for finding the blocks to Ifpack. The parameter
+     * <tt>omega</tt> specifies the relaxation parameter in the SOR
+     * preconditioner. The parameter <tt>min_diagonal</tt> can be used to make
+     * the application of the preconditioner also possible when some diagonal
+     * elements are zero. In a default application this would mean that we
+     * divide by zero, so by setting the parameter <tt>min_diagonal</tt> to a
+     * small nonzero value the SOR will work on a matrix that is not too far
+     * away from the one we want to treat. Finally, <tt>overlap</tt> governs
+     * the overlap of the partitions when the preconditioner runs in parallel,
+     * forming a so-called additive Schwarz preconditioner.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor. By default, use a block size of 1, use linear
+       * subdivision of the rows, set the damping parameter to one, we do not
+       * modify the diagonal, and there is no overlap (i.e. in parallel, we
+       * run a BlockJacobi preconditioner, where each block is inverted
+       * approximately by a block SOR).
+       */
+      AdditionalData (const unsigned int block_size = 1,
+                      const std::string  block_creation_type = "linear",
+                      const double       omega = 1,
+                      const double       min_diagonal = 0,
+                      const unsigned int overlap = 0,
+                      const unsigned int n_sweeps = 1);
+
+      /**
+       * This specifies the size of blocks.
+       */
+      unsigned int block_size;
+
+      /**
+       * Strategy for creation of blocks passed on to Ifpack block relaxation
+       * (variable 'partitioner: type') with this string as the given
+       * value. Available types in Ifpack include "linear" (i.e., divide the
+       * local range of the matrix in slices of the block size), "greedy"
+       * "metis". For a full list, see the documentation of Ifpack.
+       */
+      std::string block_creation_type;
+
+      /**
+       * 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.
+       */
+      unsigned int overlap;
+
+      /**
+       * Sets how many times the given operation should be applied during the
+       * vmult() operation.
+       */
+      unsigned int n_sweeps;
+    };
+
+    /**
+     * Take the sparse matrix the preconditioner object should be built of,
+     * and additional flags (damping parameter, overlap in parallel
+     * computations etc.) if there are any.
      */
     void initialize (const SparseMatrix   &matrix,
                      const AdditionalData &additional_data = AdditionalData());
@@ -601,63 +828,30 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional parameters
-     * to the preconditioner. The
-     * Trilinos IC decomposition
-     * allows for some fill-in, so
-     * it actually is a threshold
-     * incomplete Cholesky
-     * factorization. The amount of
-     * fill-in, and hence, the
-     * amount of memory used by
-     * this preconditioner, is
-     * controlled by the parameter
-     * <tt>ic_fill</tt>, which
-     * specifies this as a
-     * double. When forming the
-     * preconditioner, for certain
-     * problems bad conditioning
-     * (or just bad luck) can cause
-     * the preconditioner to be
-     * very poorly
-     * conditioned. Hence it can
-     * help to add diagonal
-     * perturbations to the
-     * original matrix and form the
-     * preconditioner for this
-     * slightly better
-     * matrix. <tt>ic_atol</tt> is
-     * an absolute perturbation
-     * that is added to the
-     * diagonal before forming the
-     * prec, and <tt>ic_rtol</tt>
-     * is a scaling factor $rtol
-     * \geq 1$. The last parameter
-     * specifies the overlap of the
-     * partitions when the
-     * preconditioner runs in
-     * parallel.
+     * Standardized data struct to pipe additional parameters to the
+     * preconditioner. The Trilinos IC decomposition allows for some fill-in,
+     * so it actually is a threshold incomplete Cholesky factorization. The
+     * amount of fill-in, and hence, the amount of memory used by this
+     * preconditioner, is controlled by the parameter <tt>ic_fill</tt>, which
+     * specifies this as a double. When forming the preconditioner, for
+     * certain problems bad conditioning (or just bad luck) can cause the
+     * preconditioner to be very poorly conditioned. Hence it can help to add
+     * diagonal perturbations to the original matrix and form the
+     * preconditioner for this slightly better matrix. <tt>ic_atol</tt> is an
+     * absolute perturbation that is added to the diagonal before forming the
+     * prec, and <tt>ic_rtol</tt> is a scaling factor $rtol \geq 1$. The last
+     * parameter specifies the overlap of the partitions when the
+     * preconditioner runs in parallel.
      */
     struct AdditionalData
     {
       /**
-       * Constructor. By default, set
-       * the drop tolerance to 0, the
-       * level of extra fill-ins is
-       * set to be zero (just use the
-       * matrix structure, do not
-       * generate any additional
-       * fill-in), the tolerance
-       * level are 0 and 1,
-       * respectively, and the
-       * overlap in case of a
-       * parallel execution is
-       * zero. This overlap in a
-       * block-application of the IC
-       * in the parallel case makes
-       * the preconditioner a
-       * so-called additive Schwarz
+       * Constructor. By default, set the drop tolerance to 0, the level of
+       * extra fill-ins is set to be zero (just use the matrix structure, do
+       * not generate any additional fill-in), the tolerance level are 0 and
+       * 1, respectively, and the overlap in case of a parallel execution is
+       * zero. This overlap in a block-application of the IC in the parallel
+       * case makes the preconditioner a so-called additive Schwarz
        * preconditioner.
        */
       AdditionalData (const unsigned int ic_fill = 0,
@@ -666,61 +860,38 @@ namespace TrilinosWrappers
                       const unsigned int overlap = 0);
 
       /**
-       * 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.
+       * 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;
 
       /**
-       * 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.
+       * 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.
+       * 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.
+       * 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 the
-     * preconditioner should be
-     * computed of, and additional
-     * flags if there are any.
+     * Initialize function. Takes the matrix the preconditioner should be
+     * computed of, and additional flags if there are any.
      */
     void initialize (const SparseMatrix   &matrix,
                      const AdditionalData &additional_data = AdditionalData());
@@ -770,63 +941,30 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional parameters
-     * to the preconditioner. The
-     * Trilinos ILU decomposition
-     * allows for some fill-in, so
-     * it actually is a threshold
-     * incomplete LU
-     * factorization. The amount of
-     * fill-in, and hence, the
-     * amount of memory used by
-     * this preconditioner, is
-     * controlled by the parameter
-     * <tt>ilu_fill</tt>, which
-     * specifies this as a
-     * double. When forming the
-     * preconditioner, for certain
-     * problems bad conditioning
-     * (or just bad luck) can cause
-     * the preconditioner to be
-     * very poorly
-     * conditioned. Hence it can
-     * help to add diagonal
-     * perturbations to the
-     * original matrix and form the
-     * preconditioner for this
-     * slightly better
-     * matrix. <tt>ilu_atol</tt> is
-     * an absolute perturbation
-     * that is added to the
-     * diagonal before forming the
-     * prec, and <tt>ilu_rtol</tt>
-     * is a scaling factor $rtol
-     * \geq 1$. The last parameter
-     * specifies the overlap of the
-     * partitions when the
-     * preconditioner runs in
-     * parallel.
+     * Standardized data struct to pipe additional parameters to the
+     * preconditioner. The Trilinos ILU decomposition allows for some fill-in,
+     * so it actually is a threshold incomplete LU factorization. The amount
+     * of fill-in, and hence, the amount of memory used by this
+     * preconditioner, is controlled by the parameter <tt>ilu_fill</tt>, which
+     * specifies this as a double. When forming the preconditioner, for
+     * certain problems bad conditioning (or just bad luck) can cause the
+     * preconditioner to be very poorly conditioned. Hence it can help to add
+     * diagonal perturbations to the original matrix and form the
+     * preconditioner for this slightly better matrix. <tt>ilu_atol</tt> is an
+     * absolute perturbation that is added to the diagonal before forming the
+     * prec, and <tt>ilu_rtol</tt> is a scaling factor $rtol \geq 1$. The last
+     * parameter specifies the overlap of the partitions when the
+     * preconditioner runs in parallel.
      */
     struct AdditionalData
     {
       /**
-       * Constructor. By default, the
-       * level of extra fill-ins is
-       * set to be zero (just use the
-       * matrix structure, do not
-       * generate any additional
-       * fill-in), the tolerance
-       * level are 0 and 1,
-       * respectively, and the
-       * overlap in case of a
-       * parallel execution is
-       * zero. This overlap in a
-       * block-application of the ILU
-       * in the parallel case makes
-       * the preconditioner a
-       * so-called additive Schwarz
-       * preconditioner.
+       * Constructor. By default, the level of extra fill-ins is set to be
+       * zero (just use the matrix structure, do not generate any additional
+       * fill-in), the tolerance level are 0 and 1, respectively, and the
+       * overlap in case of a parallel execution is zero. This overlap in a
+       * block-application of the ILU in the parallel case makes the
+       * preconditioner a so-called additive Schwarz preconditioner.
        */
       AdditionalData (const unsigned int ilu_fill = 0,
                       const double       ilu_atol = 0.,
@@ -834,61 +972,38 @@ namespace TrilinosWrappers
                       const unsigned int overlap  = 0);
 
       /**
-       * 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.
+       * 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;
 
       /**
-       * 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.
+       * 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.
+       * 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.
+       * 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.
+     * 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());
@@ -940,58 +1055,33 @@ namespace TrilinosWrappers
   {
   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.
+     * 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.
+       * 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,
@@ -1000,69 +1090,44 @@ namespace TrilinosWrappers
                       const unsigned int overlap  = 0);
 
       /**
-       * This specifies the relative size
-       * of elements which should be
-       * dropped when forming an incomplete
-       * LU decomposition with threshold.
+       * 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.
+       * 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.
+       * 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.
+       * 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.
+       * 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.
+     * 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());
@@ -1094,9 +1159,8 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional parameters
-     * to the preconditioner.
+     * Standardized data struct to pipe additional parameters to the
+     * preconditioner.
      */
     struct AdditionalData
     {
@@ -1107,21 +1171,15 @@ namespace TrilinosWrappers
 
 
       /**
-       * This determines how large
-       * the overlap of the local
-       * matrix portions on each
-       * processor in a parallel
-       * application should be.
+       * 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.
+     * 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());
@@ -1146,9 +1204,8 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional parameters
-     * to the preconditioner.
+     * Standardized data struct to pipe additional parameters to the
+     * preconditioner.
      */
     struct AdditionalData
     {
@@ -1163,74 +1220,51 @@ namespace TrilinosWrappers
                       const bool         nonzero_starting = false);
 
       /**
-       * This determines the degree of the
-       * Chebyshev polynomial. The degree
-       * of the polynomial gives the number
-       * of matrix-vector products to be
-       * performed for one application of
-       * the vmult() operation.
+       * This determines the degree of the Chebyshev polynomial. The degree of
+       * the polynomial gives the number of matrix-vector products to be
+       * performed for one application of the vmult() operation.
        */
       unsigned int degree;
 
       /**
-       * This sets the maximum eigenvalue
-       * of the matrix, which needs to be
-       * set properly for appropriate
-       * performance of the Chebyshev
-       * preconditioner.
+       * This sets the maximum eigenvalue of the matrix, which needs to be set
+       * properly for appropriate performance of the Chebyshev preconditioner.
        */
       double max_eigenvalue;
 
       /**
-       * This sets the ratio between the
-       * maximum and the minimum
-       * eigenvalue.
+       * This sets the ratio between the maximum and the minimum eigenvalue.
        */
       double eigenvalue_ratio;
 
       /**
-       * This sets the minimum eigenvalue,
-       * which is an optional parameter
-       * only used internally for checking
-       * whether we use an identity matrix.
+       * This sets the minimum eigenvalue, which is an optional parameter only
+       * used internally for checking whether we use an identity matrix.
        */
       double min_eigenvalue;
 
       /**
-       * This sets a threshold below which
-       * the diagonal element will not be
-       * inverted in the Chebyshev
-       * algorithm.
+       * This sets a threshold below which the diagonal element will not be
+       * inverted in the Chebyshev algorithm.
        */
       double min_diagonal;
 
       /**
-       * When this flag is set to
-       * <tt>true</tt>, it enables the
-       * method <tt>vmult(dst, src)</tt> to
-       * use non-zero data in the vector
-       * <tt>dst</tt>, appending to it the
-       * Chebyshev corrections. This can be
-       * useful in some situations
-       * (e.g. when used for high-frequency
-       * error smoothing), but not the way
-       * the solver classes expect a
-       * preconditioner to work (where one
-       * ignores the content in
-       * <tt>dst</tt> for the
-       * preconditioner application). The
-       * user should really know what she
-       * is doing when touching this flag.
+       * When this flag is set to <tt>true</tt>, it enables the method
+       * <tt>vmult(dst, src)</tt> to use non-zero data in the vector
+       * <tt>dst</tt>, appending to it the Chebyshev corrections. This can be
+       * useful in some situations (e.g. when used for high-frequency error
+       * smoothing), but not the way the solver classes expect a
+       * preconditioner to work (where one ignores the content in <tt>dst</tt>
+       * for the preconditioner application). The user should really know what
+       * she is doing when touching this flag.
        */
       bool nonzero_starting;
     };
 
     /**
-     * Initialize function. Takes
-     * the matrix which is used to
-     * form the preconditioner, and
-     * additional flags if there
-     * are any.
+     * 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());
@@ -1591,30 +1625,15 @@ namespace TrilinosWrappers
   }
 
 
-  // For the implementation of
-  // the <code>vmult</code>
-  // function with deal.II data
-  // structures we note that
-  // invoking a call of the
-  // Trilinos preconditioner
-  // requires us to use Epetra
-  // vectors as well. We do this
-  // by providing a view, i.e.,
-  // feed Trilinos with a
-  // pointer to the data, so we
-  // avoid copying the content
-  // of the vectors during the
-  // iteration (this function is
-  // only useful when used in
-  // serial anyway). 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.
+  // For the implementation of the <code>vmult</code> function with deal.II
+  // data structures we note that invoking a call of the Trilinos
+  // preconditioner requires us to use Epetra vectors as well. We do this by
+  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
+  // avoid copying the content of the vectors during the iteration (this
+  // function is only useful when used in serial anyway). 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.
   inline
   void PreconditionBase::vmult (dealii::Vector<double>       &dst,
                                 const dealii::Vector<double> &src) const
index e4db7d8632d07e46aae86a140065c38a36c61a20..13d093e154645883fad5a4ea1045089b2800f769 100644 (file)
@@ -100,10 +100,12 @@ namespace TrilinosWrappers
 
   PreconditionJacobi::AdditionalData::
   AdditionalData (const double omega,
-                  const double min_diagonal)
+                  const double min_diagonal,
+                  const unsigned int n_sweeps)
     :
     omega (omega),
-    min_diagonal (min_diagonal)
+    min_diagonal (min_diagonal),
+    n_sweeps     (n_sweeps)
   {}
 
 
@@ -127,7 +129,7 @@ namespace TrilinosWrappers
     int ierr;
 
     Teuchos::ParameterList parameter_list;
-    parameter_list.set ("relaxation: sweeps", 1);
+    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
     parameter_list.set ("relaxation: type", "Jacobi");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
     parameter_list.set ("relaxation: min diagonal value",
@@ -150,11 +152,13 @@ namespace TrilinosWrappers
   PreconditionSSOR::AdditionalData::
   AdditionalData (const double       omega,
                   const double       min_diagonal,
-                  const unsigned int overlap)
+                  const unsigned int overlap,
+                  const unsigned int n_sweeps)
     :
     omega        (omega),
     min_diagonal (min_diagonal),
-    overlap      (overlap)
+    overlap      (overlap),
+    n_sweeps     (n_sweeps)
   {}
 
 
@@ -177,7 +181,7 @@ namespace TrilinosWrappers
     int ierr;
 
     Teuchos::ParameterList parameter_list;
-    parameter_list.set ("relaxation: sweeps", 1);
+    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
     parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
     parameter_list.set ("relaxation: min diagonal value",
@@ -201,11 +205,13 @@ namespace TrilinosWrappers
   PreconditionSOR::AdditionalData::
   AdditionalData (const double       omega,
                   const double       min_diagonal,
-                  const unsigned int overlap)
+                  const unsigned int overlap,
+                  const unsigned int n_sweeps)
     :
     omega        (omega),
     min_diagonal (min_diagonal),
-    overlap      (overlap)
+    overlap      (overlap),
+    n_sweeps     (n_sweeps)
   {}
 
 
@@ -228,12 +234,193 @@ namespace TrilinosWrappers
     int ierr;
 
     Teuchos::ParameterList parameter_list;
-    parameter_list.set ("relaxation: sweeps", 1);
+    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
+    parameter_list.set ("relaxation: type", "Gauss-Seidel");
+    parameter_list.set ("relaxation: damping factor", additional_data.omega);
+    parameter_list.set ("relaxation: min diagonal value",
+                        additional_data.min_diagonal);
+    parameter_list.set ("schwarz: combine mode", "Add");
+
+    ierr = 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));
+  }
+
+
+
+  /* ----------------------- PreconditionBlockJacobi ---------------------- */
+
+  PreconditionBlockJacobi::AdditionalData::
+  AdditionalData (const unsigned int block_size,
+                  const std::string  block_creation_type,
+                  const double omega,
+                  const double min_diagonal,
+                  const unsigned int n_sweeps)
+    :
+    block_size(block_size),
+    block_creation_type(block_creation_type),
+    omega (omega),
+    min_diagonal (min_diagonal),
+    n_sweeps     (n_sweeps)
+  {}
+
+
+
+  void
+  PreconditionBlockJacobi::initialize (const SparseMatrix   &matrix,
+                                       const AdditionalData &additional_data)
+  {
+    // release memory before reallocation
+    preconditioner.reset ();
+    preconditioner.reset (Ifpack().Create
+                          ("block relaxation",
+                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                           0));
+
+    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
+                                    (preconditioner.get());
+    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                     "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
+    parameter_list.set ("relaxation: type", "Jacobi");
+    parameter_list.set ("relaxation: damping factor", additional_data.omega);
+    parameter_list.set ("relaxation: min diagonal value",
+                        additional_data.min_diagonal);
+    parameter_list.set ("partitioner: type", additional_data.block_creation_type);
+    int n_local_parts = (matrix.trilinos_matrix().NumMyRows()+additional_data.
+                         block_size-1)/additional_data.block_size;
+    parameter_list.set ("partitioner: local parts", n_local_parts);
+
+    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));
+  }
+
+
+
+  /* ----------------------- PreconditionBlockSSOR ------------------------ */
+
+  PreconditionBlockSSOR::AdditionalData::
+  AdditionalData (const unsigned int block_size,
+                  const std::string  block_creation_type,
+                  const double       omega,
+                  const double       min_diagonal,
+                  const unsigned int overlap,
+                  const unsigned int n_sweeps)
+    :
+    block_size(block_size),
+    block_creation_type(block_creation_type),
+    omega        (omega),
+    min_diagonal (min_diagonal),
+    overlap      (overlap),
+    n_sweeps     (n_sweeps)
+  {}
+
+
+
+  void
+  PreconditionBlockSSOR::initialize (const SparseMatrix   &matrix,
+                                     const AdditionalData &additional_data)
+  {
+    preconditioner.reset ();
+    preconditioner.reset (Ifpack().Create
+                          ("block relaxation",
+                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                           additional_data.overlap));
+
+    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
+                                    (preconditioner.get());
+    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                     "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
+    parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
+    parameter_list.set ("relaxation: damping factor", additional_data.omega);
+    parameter_list.set ("relaxation: min diagonal value",
+                        additional_data.min_diagonal);
+    parameter_list.set ("schwarz: combine mode", "Add");
+    parameter_list.set ("partitioner: type", additional_data.block_creation_type);
+    int n_local_parts = (matrix.trilinos_matrix().NumMyRows()+additional_data.
+                         block_size-1)/additional_data.block_size;
+    parameter_list.set ("partitioner: local parts", n_local_parts);
+
+    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));
+  }
+
+
+
+  /* ------------------------ PreconditionBlockSOR ------------------------ */
+
+  PreconditionBlockSOR::AdditionalData::
+  AdditionalData (const unsigned int block_size,
+                  const std::string  block_creation_type,
+                  const double       omega,
+                  const double       min_diagonal,
+                  const unsigned int overlap,
+                  const unsigned int n_sweeps)
+    :
+    block_size(block_size),
+    block_creation_type(block_creation_type),
+    omega        (omega),
+    min_diagonal (min_diagonal),
+    overlap      (overlap),
+    n_sweeps     (n_sweeps)
+  {}
+
+
+
+  void
+  PreconditionBlockSOR::initialize (const SparseMatrix   &matrix,
+                                    const AdditionalData &additional_data)
+  {
+    preconditioner.reset ();
+    preconditioner.reset (Ifpack().Create
+                          ("block relaxation",
+                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                           additional_data.overlap));
+
+    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
+                                    (preconditioner.get());
+    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                     "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
     parameter_list.set ("relaxation: type", "Gauss-Seidel");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
     parameter_list.set ("relaxation: min diagonal value",
                         additional_data.min_diagonal);
     parameter_list.set ("schwarz: combine mode", "Add");
+    parameter_list.set ("partitioner: type", additional_data.block_creation_type);
+    int n_local_parts = (matrix.trilinos_matrix().NumMyRows()+additional_data.
+                         block_size-1)/additional_data.block_size;
+    parameter_list.set ("partitioner: local parts", n_local_parts);
 
     ierr = ifpack->SetParameters(parameter_list);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
diff --git a/tests/trilinos/precondition.cc b/tests/trilinos/precondition.cc
new file mode 100644 (file)
index 0000000..2ec84a4
--- /dev/null
@@ -0,0 +1,419 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// solves a 2D Poisson equation for linear elements with the simple
+// (Ifpack-based) Trilinos preconditioners
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/base/function.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4 ();
+  void run ();
+
+private:
+  void make_grid ();
+  void setup_system();
+  void assemble_system ();
+  void solve ();
+
+  Triangulation<dim>   triangulation;
+  FE_Q<dim>            fe;
+  DoFHandler<dim>      dof_handler;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                  const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    return_value += 4*std::pow(p(i), 4);
+
+  return return_value;
+}
+
+
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                   const unsigned int /*component*/) const
+{
+  return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4 ()
+  :
+  fe (1),
+  dof_handler (triangulation)
+{}
+
+
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (6);
+}
+
+
+
+template <int dim>
+void Step4<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  constraints.clear();
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close();
+
+  CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false);
+  system_matrix.reinit (c_sparsity);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void Step4<dim>::assemble_system ()
+{
+  QGauss<dim>  quadrature_formula(fe.degree+1);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values   | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
+                                   fe_values.shape_grad (j, q_point) *
+                                   fe_values.JxW (q_point));
+
+            cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                            right_hand_side.value (fe_values.quadrature_point (q_point)) *
+                            fe_values.JxW (q_point));
+          }
+
+      cell->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global(cell_matrix, cell_rhs,
+                                             local_dof_indices,
+                                             system_matrix, system_rhs);
+    }
+  system_matrix.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void Step4<dim>::solve ()
+{
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5));
+
+  {
+    deallog.push("Jacobi");
+    TrilinosWrappers::PreconditionJacobi preconditioner;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("SSOR");
+    TrilinosWrappers::PreconditionSSOR preconditioner;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("SOR");
+    TrilinosWrappers::PreconditionSOR preconditioner;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverBicgstab<>        solver (solver_control);
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("BlockJacobi");
+    TrilinosWrappers::PreconditionBlockJacobi preconditioner;
+    TrilinosWrappers::PreconditionBlockJacobi::AdditionalData data;
+    data.block_size = 4;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix, data);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("BlockSSOR");
+    TrilinosWrappers::PreconditionBlockSSOR preconditioner;
+    TrilinosWrappers::PreconditionBlockSSOR::AdditionalData data;
+    data.block_size = 4;
+    data.omega = 1.2;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix, data);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("BlockSOR");
+    TrilinosWrappers::PreconditionBlockSOR preconditioner;
+    TrilinosWrappers::PreconditionBlockSOR::AdditionalData data;
+    data.block_size = 4;
+    data.omega = 0.8;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverBicgstab<>        solver (solver_control);
+    preconditioner.initialize(system_matrix, data);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("IC");
+    TrilinosWrappers::PreconditionIC preconditioner;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("ILU");
+    TrilinosWrappers::PreconditionILU preconditioner;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("ILUT");
+    TrilinosWrappers::PreconditionILUT preconditioner;
+    TrilinosWrappers::PreconditionILUT::AdditionalData data;
+    data.ilut_drop = 1e-6;
+    data.ilut_fill = 3;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverBicgstab<>        solver (solver_control);
+    preconditioner.initialize(system_matrix, data);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("Chebyshev");
+    TrilinosWrappers::PreconditionChebyshev preconditioner;
+    TrilinosWrappers::PreconditionChebyshev::AdditionalData data;
+    data.max_eigenvalue = 2.5;
+    data.degree = 3;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix, data);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("Direct");
+    TrilinosWrappers::PreconditionBlockwiseDirect preconditioner;
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-10);
+    SolverCG<>              solver (solver_control);
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    deallog.pop();
+  }
+
+
+  deallog.pop();
+}
+
+
+
+template <int dim>
+void Step4<dim>::run()
+{
+  for (unsigned int cycle = 0; cycle < 2; ++cycle)
+    {
+      if (cycle == 0)
+        make_grid();
+      else
+        triangulation.refine_global(1);
+
+      setup_system();
+      assemble_system();
+      solve();
+    }
+}
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+  try
+    {
+      Step4<2> test;
+      test.run();
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/precondition.output b/tests/trilinos/precondition.output
new file mode 100644 (file)
index 0000000..395cc23
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL:04225:Jacobi:cg::Starting value 21.8299
+DEAL:04225:Jacobi:cg::Convergence step 101 value 0
+DEAL:04225:SSOR:cg::Starting value 21.8299
+DEAL:04225:SSOR:cg::Convergence step 78 value 0
+DEAL:04225:SOR:Bicgstab::Starting value 21.8299
+DEAL:04225:SOR:Bicgstab::Convergence step 98 value 0
+DEAL:04225:BlockJacobi:cg::Starting value 21.8299
+DEAL:04225:BlockJacobi:cg::Convergence step 163 value 0
+DEAL:04225:BlockSSOR:cg::Starting value 21.8299
+DEAL:04225:BlockSSOR:cg::Convergence step 60 value 0
+DEAL:04225:BlockSOR:Bicgstab::Starting value 21.8299
+DEAL:04225:BlockSOR:Bicgstab::Convergence step 85 value 0
+DEAL:04225:IC:cg::Starting value 21.8299
+DEAL:04225:IC:cg::Convergence step 70 value 0
+DEAL:04225:ILU:cg::Starting value 21.8299
+DEAL:04225:ILU:cg::Convergence step 57 value 0
+DEAL:04225:ILUT:Bicgstab::Starting value 21.8299
+DEAL:04225:ILUT:Bicgstab::Convergence step 32 value 0
+DEAL:04225:Chebyshev:cg::Starting value 21.8299
+DEAL:04225:Chebyshev:cg::Convergence step 47 value 0
+DEAL:04225:Direct:cg::Starting value 21.8299
+DEAL:04225:Direct:cg::Convergence step 1 value 0
+DEAL:16641:Jacobi:cg::Starting value 30.8770
+DEAL:16641:Jacobi:cg::Convergence step 201 value 0
+DEAL:16641:SSOR:cg::Starting value 30.8770
+DEAL:16641:SSOR:cg::Convergence step 141 value 0
+DEAL:16641:SOR:Bicgstab::Starting value 30.8770
+DEAL:16641:SOR:Bicgstab::Convergence step 209 value 0
+DEAL:16641:BlockJacobi:cg::Starting value 30.8770
+DEAL:16641:BlockJacobi:cg::Convergence step 321 value 0
+DEAL:16641:BlockSSOR:cg::Starting value 30.8770
+DEAL:16641:BlockSSOR:cg::Convergence step 117 value 0
+DEAL:16641:BlockSOR:Bicgstab::Starting value 30.8770
+DEAL:16641:BlockSOR:Bicgstab::Convergence step 173 value 0
+DEAL:16641:IC:cg::Starting value 30.8770
+DEAL:16641:IC:cg::Convergence step 134 value 0
+DEAL:16641:ILU:cg::Starting value 30.8770
+DEAL:16641:ILU:cg::Convergence step 110 value 0
+DEAL:16641:ILUT:Bicgstab::Starting value 30.8770
+DEAL:16641:ILUT:Bicgstab::Convergence step 59 value 0
+DEAL:16641:Chebyshev:cg::Starting value 30.8770
+DEAL:16641:Chebyshev:cg::Convergence step 94 value 0
+DEAL:16641:Direct:cg::Starting value 30.8770
+DEAL:16641:Direct:cg::Convergence step 1 value 0

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.