From 3227411fba5ba5c285975fcb052adf0be37fefc3 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 23 May 2022 14:13:07 +0200 Subject: [PATCH] PreconditionChebyshev: Implement power iteration for eigenvalue estimate --- include/deal.II/lac/precondition.h | 212 ++++++++++++------ ...ition_chebyshev_03.with_lapack=true.output | 6 +- 2 files changed, 152 insertions(+), 66 deletions(-) diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 4303f33d45..fb4d54b67b 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -1435,23 +1435,26 @@ public: *

Estimation of the eigenvalues

* * The Chebyshev method relies on an estimate of the eigenvalues of the matrix - * which are computed during the first invocation of vmult(). The algorithm - * invokes a conjugate gradient solver (i.e., Lanczos iteration) so symmetry - * and positive definiteness of the (preconditioned) matrix system are - * required. The eigenvalue algorithm can be controlled by + * which are computed during the first invocation of vmult(). This class + * offers several algorithms to this end, see + * PreconditionChebyshev::AdditionalData::EigenvalueAlgorithm. The default + * algorithm invokes the Lanczos method via the SolverCG class, which requires + * symmetry and positive definiteness of the (preconditioned) matrix system + * are required. Also note that deal.II needs to be configured with LAPACK + * support to use this option. The eigenvalue algorithm can be controlled by * PreconditionChebyshev::AdditionalData::eig_cg_n_iterations specifying how - * many iterations should be performed. The iterations are started from an - * initial vector that depends on the vector type. For the classes - * dealii::Vector or dealii::LinearAlgebra::distributed::Vector, which have - * fast element access, it is a vector with entries `(-5.5, -4.5, -3.5, - * -2.5, ..., 3.5, 4.5, 5.5)` with appropriate epilogue and adjusted such that - * its mean is always zero, which works well for the Laplacian. This setup is - * stable in parallel in the sense that for a different number of processors - * but the same ordering of unknowns, the same initial vector and thus - * eigenvalue distribution will be computed, apart from roundoff errors. For - * other vector types, the initial vector contains all ones, scaled by the - * length of the vector, except for the very first entry that is zero, - * triggering high-frequency content again. + * many iterations should be performed. For all algorithms, the iterative + * process is started from an initial vector that depends on the vector + * type. For the classes dealii::Vector or + * dealii::LinearAlgebra::distributed::Vector, which have fast element access, + * it is a vector with entries `(-5.5, -4.5, -3.5, -2.5, ..., 3.5, 4.5, 5.5)` + * with appropriate epilogue and adjusted such that its mean is always zero, + * which works well for the Laplacian. This setup is stable in parallel in the + * sense that for a different number of processors but the same ordering of + * unknowns, the same initial vector and thus eigenvalue distribution will be + * computed, apart from roundoff errors. For other vector types, the initial + * vector contains all ones, scaled by the length of the vector, except for + * the very first entry that is zero, triggering high-frequency content again. * * The computation of eigenvalues happens the first time one of the vmult(), * Tvmult(), step() or Tstep() functions is called or when @@ -1459,9 +1462,10 @@ public: * necessary to provide a temporary vector of the same layout as the source * and destination vectors used during application of the preconditioner. * - * The estimates for minimum and maximum eigenvalue are taken from SolverCG - * (even if the solver did not converge in the requested number of - * iterations). Finally, the maximum eigenvalue is multiplied by a safety + * The estimates for minimum and maximum eigenvalue are taken from the + * underlying solver or eigenvalue algorithm in the given number of + * iterations, even if the solver did not converge in the requested number of + * iterations. Finally, the maximum eigenvalue is multiplied by a safety * factor of 1.2. * * Due to the cost of the eigenvalue estimate, this class is most appropriate @@ -1471,12 +1475,11 @@ public: *

Bypassing the eigenvalue computation

* * In some contexts, the automatic eigenvalue computation of this class may - * result in bad quality, or it may be unstable when used in parallel with - * different enumerations of the degrees of freedom, making computations - * strongly dependent on the parallel configuration. It is possible to bypass - * the automatic eigenvalue computation by setting - * AdditionalData::eig_cg_n_iterations to zero, and provide the variable - * AdditionalData::max_eigenvalue instead. The minimal eigenvalue is + * result in a bad quality, e.g. when the polynomial basis or numbering of + * unknowns is such that the initial vector described above is a bad + * choice. It is possible to bypass the automatic eigenvalue computation by + * setting AdditionalData::eig_cg_n_iterations to zero, and provide the + * variable AdditionalData::max_eigenvalue instead. The minimal eigenvalue is * implicitly specified via `max_eigenvalue/smoothing_range`. *

Using the PreconditionChebyshev as a solver

@@ -1489,8 +1492,10 @@ public: * * In order to use Chebyshev as a solver, set the degree to * numbers::invalid_unsigned_int to force the automatic computation of the - * number of iterations needed to reach a given target tolerance. In this - * case, the target tolerance is read from the variable + * number of iterations needed to reach a given target tolerance. Note that + * this currently only works for symmetric positive definite matrices with the + * eigenvalue algorithm set to the conjugate gradient algorithm. In this case, + * the target tolerance is read from the variable * PreconditionChebyshev::AdditionalData::smoothing_range (it needs to be a * number less than one to force any iterations obviously). * @@ -1597,14 +1602,39 @@ public: */ struct AdditionalData { + /** + * An enum to define the available types of eigenvalue estimation + * algorithms. + */ + enum class EigenvalueAlgorithm + { + /** + * This option runs the conjugate gradient solver and computes an + * eigenvalue estimation from the underlying Lanczos space. This only + * works for symmetric positive definite matrices. + */ + lanczos, + /** + * This option runs a power iteration to estimate the largest + * eigenvalue. This algorithm also works for non-symmetric matrices, + * but typically gives less accurate estimates than the option 'lanczos' + * because it does not take the relation between vectors in the iterations + * into account (roughly speaking the off-diagonal entries in the + * tri-diagonal matrix of the Lanczos iteration). + */ + power_iteration + }; + /** * Constructor. */ - AdditionalData(const unsigned int degree = 1, - const double smoothing_range = 0., - const unsigned int eig_cg_n_iterations = 8, - const double eig_cg_residual = 1e-2, - const double max_eigenvalue = 1); + AdditionalData(const unsigned int degree = 1, + const double smoothing_range = 0., + const unsigned int eig_cg_n_iterations = 8, + const double eig_cg_residual = 1e-2, + const double max_eigenvalue = 1, + const EigenvalueAlgorithm eigenvalue_algorithm = + EigenvalueAlgorithm::lanczos); /** * Copy assignment operator. @@ -1670,6 +1700,11 @@ public: * Stores the preconditioner object that the Chebyshev is wrapped around. */ std::shared_ptr preconditioner; + + /** + * Specifies the underlying eigenvalue estimation algorithm. + */ + EigenvalueAlgorithm eigenvalue_algorithm; }; @@ -2766,6 +2801,41 @@ namespace internal std::vector values; }; + + + + template + double + power_iteration(const MatrixType & matrix, + VectorType & eigenvector, + const PreconditionerType &preconditioner, + const unsigned int n_iterations) + { + double eigenvalue_estimate = 0.; + eigenvector /= eigenvector.l2_norm(); + VectorType vector1, vector2; + vector1.reinit(eigenvector, true); + if (!std::is_same::value) + vector2.reinit(eigenvector, true); + for (unsigned int i = 0; i < n_iterations; ++i) + { + if (!std::is_same::value) + { + matrix.vmult(vector2, eigenvector); + preconditioner.vmult(vector1, vector2); + } + else + matrix.vmult(vector1, eigenvector); + + eigenvalue_estimate = eigenvector * vector1; + + vector1 /= vector1.l2_norm(); + eigenvector.swap(vector1); + } + return eigenvalue_estimate; + } } // namespace PreconditionChebyshevImplementation } // namespace internal @@ -2773,16 +2843,18 @@ namespace internal template inline PreconditionChebyshev:: - AdditionalData::AdditionalData(const unsigned int degree, - const double smoothing_range, - const unsigned int eig_cg_n_iterations, - const double eig_cg_residual, - const double max_eigenvalue) + AdditionalData::AdditionalData(const unsigned int degree, + const double smoothing_range, + const unsigned int eig_cg_n_iterations, + const double eig_cg_residual, + const double max_eigenvalue, + const EigenvalueAlgorithm eigenvalue_algorithm) : degree(degree) , smoothing_range(smoothing_range) , eig_cg_n_iterations(eig_cg_n_iterations) , eig_cg_residual(eig_cg_residual) , max_eigenvalue(max_eigenvalue) + , eigenvalue_algorithm(eigenvalue_algorithm) {} @@ -2794,12 +2866,13 @@ inline typename PreconditionChebyshev:: AdditionalData::operator=(const AdditionalData &other_data) { - degree = other_data.degree; - smoothing_range = other_data.smoothing_range; - eig_cg_n_iterations = other_data.eig_cg_n_iterations; - eig_cg_residual = other_data.eig_cg_residual; - max_eigenvalue = other_data.max_eigenvalue; - preconditioner = other_data.preconditioner; + degree = other_data.degree; + smoothing_range = other_data.smoothing_range; + eig_cg_n_iterations = other_data.eig_cg_n_iterations; + eig_cg_residual = other_data.eig_cg_residual; + max_eigenvalue = other_data.max_eigenvalue; + preconditioner = other_data.preconditioner; + eigenvalue_algorithm = other_data.eigenvalue_algorithm; constraints.copy_from(other_data.constraints); return *this; @@ -2878,22 +2951,8 @@ PreconditionChebyshev:: ExcMessage( "Need to set at least two iterations to find eigenvalues.")); - // set a very strict tolerance to force at least two iterations - ReductionControl control( - data.eig_cg_n_iterations, - std::sqrt( - std::numeric_limits::epsilon()), - 1e-10, - false, - false); - internal::PreconditionChebyshevImplementation::EigenvalueTracker - eigenvalue_tracker; - SolverCG solver(control); - solver.connect_eigenvalues_slot( - [&eigenvalue_tracker](const std::vector &eigenvalues) { - eigenvalue_tracker.slot(eigenvalues); - }); + eigenvalue_tracker; // set an initial guess that contains some high-frequency parts (to the // extent possible without knowing the discretization and the numbering) @@ -2902,15 +2961,44 @@ PreconditionChebyshev:: temp_vector1); data.constraints.set_zero(temp_vector1); - try + if (data.eigenvalue_algorithm == + AdditionalData::EigenvalueAlgorithm::lanczos) { + // set a very strict tolerance to force at least two iterations + IterationNumberControl control(data.eig_cg_n_iterations, + 1e-10, + false, + false); + + SolverCG solver(control); + solver.connect_eigenvalues_slot( + [&eigenvalue_tracker](const std::vector &eigenvalues) { + eigenvalue_tracker.slot(eigenvalues); + }); + solver.solve(*matrix_ptr, solution_old, temp_vector1, *data.preconditioner); + + info.cg_iterations = control.last_step(); } - catch (SolverControl::NoConvergence &) - {} + else if (data.eigenvalue_algorithm == + AdditionalData::EigenvalueAlgorithm::power_iteration) + { + Assert(data.degree != numbers::invalid_unsigned_int, + ExcMessage("Cannot estimate the minimal eigenvalue with the " + "power iteration")); + + eigenvalue_tracker.values.push_back( + internal::PreconditionChebyshevImplementation::power_iteration( + *matrix_ptr, + temp_vector1, + *data.preconditioner, + data.eig_cg_n_iterations)); + } + else + Assert(false, ExcNotImplemented()); // read the eigenvalues from the attached eigenvalue tracker if (eigenvalue_tracker.values.empty()) @@ -2923,8 +3011,6 @@ PreconditionChebyshev:: // be converged info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back(); } - - info.cg_iterations = control.last_step(); } else { diff --git a/tests/lac/precondition_chebyshev_03.with_lapack=true.output b/tests/lac/precondition_chebyshev_03.with_lapack=true.output index d8ef70a6e4..48a71cb206 100644 --- a/tests/lac/precondition_chebyshev_03.with_lapack=true.output +++ b/tests/lac/precondition_chebyshev_03.with_lapack=true.output @@ -1,8 +1,8 @@ DEAL::Size 4 Unknowns 9 -DEAL::Residual step i=0: ilu=0.3321, cheby=0.0789 -DEAL::Residual step i=1: ilu=0.3592, cheby=0.0789 -DEAL::Residual step i=2: ilu=0.1717, cheby=0.0507 +DEAL::Residual step i=0: ilu=0.3321, cheby=0.0794 +DEAL::Residual step i=1: ilu=0.3592, cheby=0.0794 +DEAL::Residual step i=2: ilu=0.1717, cheby=0.0509 DEAL::Size 8 Unknowns 49 DEAL::Residual step i=0: ilu=1.9436, cheby=0.1549 DEAL::Residual step i=1: ilu=2.1487, cheby=0.1583 -- 2.39.5