Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

PreconditionChebyshev: Implement power iteration for eigenvalue estimate #13792

Merged
merged 3 commits into from
May 23, 2022
Merged

PreconditionChebyshev: Implement power iteration for eigenvalue estimate #13792

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3227411
PreconditionChebyshev: Implement power iteration for eigenvalue estimate
kronbichler May 23, 2022
359e3cd
Add test case
kronbichler May 23, 2022
ff356b5
Changelog
kronbichler May 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
7 changes: 7 additions & 0 deletions doc/news/changes/minor/20220523Kronbichler
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
New: PreconditionChebyshev can now also estimate eigenvalues with the power
method, rather than the Lanczos algorithm underlying the SolverCG
estimation. This is useful for matrices that are slightly non-symmetric such
that the eigenvalue estimate of SolverCG breaks down, but not too much so that
the Chebyshev iteration is still a useful iteration.
<br>
(Martin Kronbichler, 2022/05/23)
212 changes: 149 additions & 63 deletions include/deal.II/lac/precondition.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1435,33 +1435,37 @@ class PreconditionPSOR
* <h4>Estimation of the eigenvalues</h4>
*
* 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
* estimate_eigenvalues() is called directly. In the latter case, it is
* 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
Expand All @@ -1471,12 +1475,11 @@ class PreconditionPSOR
* <h4>Bypassing the eigenvalue computation</h4>
*
* 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`.

* <h4>Using the PreconditionChebyshev as a solver</h4>
Expand All @@ -1489,8 +1492,10 @@ class PreconditionPSOR
*
* 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).
*
Expand Down Expand Up @@ -1597,14 +1602,39 @@ class PreconditionChebyshev : public Subscriptor
*/
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.
Expand Down Expand Up @@ -1670,6 +1700,11 @@ class PreconditionChebyshev : public Subscriptor
* Stores the preconditioner object that the Chebyshev is wrapped around.
*/
std::shared_ptr<PreconditionerType> preconditioner;

/**
* Specifies the underlying eigenvalue estimation algorithm.
*/
EigenvalueAlgorithm eigenvalue_algorithm;
};


Expand Down Expand Up @@ -2766,23 +2801,60 @@ namespace internal

std::vector<double> values;
};



template <typename MatrixType,
typename VectorType,
typename PreconditionerType>
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<PreconditionerType, PreconditionIdentity>::value)
vector2.reinit(eigenvector, true);
for (unsigned int i = 0; i < n_iterations; ++i)
{
if (!std::is_same<PreconditionerType, PreconditionIdentity>::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



template <typename MatrixType, class VectorType, typename PreconditionerType>
inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
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)
{}


Expand All @@ -2794,12 +2866,13 @@ inline typename PreconditionChebyshev<MatrixType,
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
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;
Expand Down Expand Up @@ -2878,22 +2951,8 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
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<typename VectorType::value_type>::epsilon()),
1e-10,
false,
false);

internal::PreconditionChebyshevImplementation::EigenvalueTracker
eigenvalue_tracker;
SolverCG<VectorType> solver(control);
solver.connect_eigenvalues_slot(
[&eigenvalue_tracker](const std::vector<double> &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)
Expand All @@ -2902,15 +2961,44 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
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<VectorType> solver(control);
solver.connect_eigenvalues_slot(
[&eigenvalue_tracker](const std::vector<double> &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())
Expand All @@ -2923,8 +3011,6 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
// be converged
info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
}

info.cg_iterations = control.last_step();
}
else
{
Expand Down
6 changes: 3 additions & 3 deletions tests/lac/precondition_chebyshev_03.with_lapack=true.output
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand Down