From 7f9921d3787bef9ad85cc28004712afee895f989 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 4 Jul 2014 13:46:44 +0000 Subject: [PATCH] Implement some new Trilinos preconditioner. Tests. git-svn-id: https://svn.dealii.org/trunk@33109 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 11 + .../deal.II/lac/trilinos_precondition.h | 1159 +++++++++-------- deal.II/source/lac/trilinos_precondition.cc | 205 ++- tests/trilinos/precondition.cc | 419 ++++++ tests/trilinos/precondition.output | 45 + 5 files changed, 1260 insertions(+), 579 deletions(-) create mode 100644 tests/trilinos/precondition.cc create mode 100644 tests/trilinos/precondition.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 99ee610c27..974eef62d3 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,17 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Martin Kronbichler, 2014/07/04) +
    2. +
  2. 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 diff --git a/deal.II/include/deal.II/lac/trilinos_precondition.h b/deal.II/include/deal.II/lac/trilinos_precondition.h index 9ab30fb98f..2d1b86302e 100644 --- a/deal.II/include/deal.II/lac/trilinos_precondition.h +++ b/deal.II/include/deal.II/lac/trilinos_precondition.h @@ -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 - * initialize function - * of the derived classes will - * have to create the - * preconditioner from a given + * Constructor. Does not do anything. The initialize 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 &dst, const dealii::Vector &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 &dst, const dealii::Vector &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 &dst, const dealii::parallel::distributed::Vector &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 &dst, const dealii::parallel::distributed::Vector &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 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 vector_distributor; }; @@ -233,69 +212,50 @@ namespace TrilinosWrappers public: /** - * Standardized data struct to - * pipe additional flags to the - * preconditioner. The - * parameter omega - * specifies the relaxation - * parameter in the Jacobi - * preconditioner. The - * parameter - * min_diagonal 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 - * min_diagonal 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 omega specifies the relaxation + * parameter in the Jacobi preconditioner. The parameter + * min_diagonal 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 min_diagonal 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 omega - * specifies the relaxation - * parameter in the SSOR - * preconditioner. The - * parameter - * min_diagonal 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 - * min_diagonal 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, - * overlap 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 omega specifies the relaxation + * parameter in the SSOR preconditioner. The parameter + * min_diagonal 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 min_diagonal 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, overlap 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 omega - * specifies the relaxation - * parameter in the SOR - * preconditioner. The - * parameter - * min_diagonal 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 - * min_diagonal 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, - * overlap 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 omega specifies the relaxation + * parameter in the SOR preconditioner. The parameter + * min_diagonal 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 min_diagonal 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, overlap 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 block + * Jacobi 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 block_size 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 block_creation_type allows to + * pass the strategy for finding the blocks to Ifpack. The parameter + * omega specifies the relaxation parameter in the SOR + * preconditioner. The parameter min_diagonal 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 min_diagonal 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, overlap 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 additive + * Schwarz method with a block SSOR approximate solve 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 block_size 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 block_creation_type allows to + * pass the strategy for finding the blocks to Ifpack. The parameter + * omega specifies the relaxation parameter in the SSOR + * preconditioner. The parameter min_diagonal 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 min_diagonal 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, overlap 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 additive + * Schwarz method with a block SOR approximate solve 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 block_size 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 block_creation_type allows to + * pass the strategy for finding the blocks to Ifpack. The parameter + * omega specifies the relaxation parameter in the SOR + * preconditioner. The parameter min_diagonal 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 min_diagonal 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, overlap 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 - * ic_fill, 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. ic_atol is - * an absolute perturbation - * that is added to the - * diagonal before forming the - * prec, and ic_rtol - * 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 ic_fill, 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. ic_atol is an + * absolute perturbation that is added to the diagonal before forming the + * prec, and ic_rtol 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 - * ic_fill 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 ic_fill 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 - * ilu_fill, 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. ilu_atol is - * an absolute perturbation - * that is added to the - * diagonal before forming the - * prec, and ilu_rtol - * 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 ilu_fill, 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. ilu_atol is an + * absolute perturbation that is added to the diagonal before forming the + * prec, and ilu_rtol 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 - * ilu_fill 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 ilu_fill 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 ilut_drop - * and ilut_fill, 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. ilut_atol - * is an absolute perturbation that - * is added to the diagonal before - * forming the prec, and - * ilu_rtol 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 ilut_drop and + * ilut_fill, 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. ilut_atol is an absolute perturbation that is + * added to the diagonal before forming the prec, and ilu_rtol 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 - * ilu_fill 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 ilu_fill 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 - * true, it enables the - * method vmult(dst, src) to - * use non-zero data in the vector - * dst, 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 - * dst for the - * preconditioner application). The - * user should really know what she - * is doing when touching this flag. + * When this flag is set to true, it enables the method + * vmult(dst, src) to use non-zero data in the vector + * dst, 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 dst + * 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 vmult - // 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 - // const in all - // deal.II calls) to - // non-constant value, as this - // is the way Trilinos wants - // to have them. + // For the implementation of the vmult 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 + // const 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 &dst, const dealii::Vector &src) const diff --git a/deal.II/source/lac/trilinos_precondition.cc b/deal.II/source/lac/trilinos_precondition.cc index e4db7d8632..13d093e154 100644 --- a/deal.II/source/lac/trilinos_precondition.cc +++ b/deal.II/source/lac/trilinos_precondition.cc @@ -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(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(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(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(&matrix.trilinos_matrix()), + 0)); + + Ifpack_Preconditioner *ifpack = static_cast + (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(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(&matrix.trilinos_matrix()), + additional_data.overlap)); + + Ifpack_Preconditioner *ifpack = static_cast + (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(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(&matrix.trilinos_matrix()), + additional_data.overlap)); + + Ifpack_Preconditioner *ifpack = static_cast + (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(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 index 0000000000..2ec84a4251 --- /dev/null +++ b/tests/trilinos/precondition.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + + +template +double RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i=0; i +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4 () + : + fe (1), + dof_handler (triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (6); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + 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 +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + const RightHandSide right_hand_side; + + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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_pointget_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 +void Step4::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 +void Step4::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 index 0000000000..395cc23bff --- /dev/null +++ b/tests/trilinos/precondition.output @@ -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 -- 2.39.5