Skip to content

Commit

Permalink
address comments, add documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
tjhei committed Mar 23, 2020
1 parent c559680 commit 2664fb2
Showing 1 changed file with 32 additions and 27 deletions.
59 changes: 32 additions & 27 deletions include/deal.II/lac/precondition.h
Original file line number Diff line number Diff line change
Expand Up @@ -914,7 +914,7 @@ class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
* 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
* requirements. The eigenvalue algorithm can be controlled by
* required. 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
Expand All @@ -929,16 +929,20 @@ class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
* 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. This is because
* temporary vectors of the same layout as the source and destination vectors
* are necessary for these computations and this information gets only
* available through vmult().
* 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.
*
* Due to the cost of the eigenvalue estimate in the first vmult(), this class
* is most appropriate if it is applied repeatedly, e.g. in a smoother for a
* geometric multigrid solver, that can in turn be used to solve several
* linear systems.
* 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
* factor of 1.2.
*
* Due to the cost of the eigenvalue estimate, this class is most appropriate
* if it is applied repeatedly, e.g. in a smoother for a geometric multigrid
* solver, that can in turn be used to solve several linear systems.
*
* <h4>Bypassing the eigenvalue computation</h4>
*
Expand Down Expand Up @@ -1161,25 +1165,26 @@ class PreconditionChebyshev : public Subscriptor
n() const;

/**
* Struct that contains information about the eigenvalue estimation performed
* by this class.
* A struct that contains information about the eigenvalue estimation
* performed by the PreconditionChebychev class.
*/
struct EigenvalueInformation
{
/**
* Estimate for the smallest eigenvalue.
*/
double min_eigenvalue;
double min_eigenvalue_estimate;
/**
* Estimate for the largest eigenvalue.
*/
double max_eigenvalue;
double max_eigenvalue_estimate;
/**
* Number of CG iterations performed or 0.
*/
unsigned int cg_iterations;
/**
* The degree (either as set or estimated).
* The degree of the Chebyshev polynomial (either as set using
* AdditionalData::degree or estimated as described there).
*/
unsigned int degree;
};
Expand Down Expand Up @@ -2325,28 +2330,28 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::

// read the eigenvalues from the attached eigenvalue tracker
if (eigenvalue_tracker.values.empty())
info.min_eigenvalue = info.max_eigenvalue = 1.;
info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
else
{
info.min_eigenvalue = eigenvalue_tracker.values.front();
info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();

// include a safety factor since the CG method will in general not
// be converged
info.max_eigenvalue = 1.2 * eigenvalue_tracker.values.back();
info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
}

info.cg_iterations = control.last_step();
}
else
{
info.max_eigenvalue = data.max_eigenvalue;
info.min_eigenvalue = data.max_eigenvalue / data.smoothing_range;
info.max_eigenvalue_estimate = data.max_eigenvalue;
info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
}

const double alpha =
(data.smoothing_range > 1. ?
info.max_eigenvalue / data.smoothing_range :
std::min(0.9 * info.max_eigenvalue, info.min_eigenvalue));
const double alpha = (data.smoothing_range > 1. ?
info.max_eigenvalue_estimate / data.smoothing_range :
std::min(0.9 * info.max_eigenvalue_estimate,
info.min_eigenvalue_estimate));

// in case the user set the degree to invalid unsigned int, we have to
// determine the number of necessary iterations from the Chebyshev error
Expand All @@ -2355,7 +2360,7 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
// R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
if (data.degree == numbers::invalid_unsigned_int)
{
const double actual_range = info.max_eigenvalue / alpha;
const double actual_range = info.max_eigenvalue_estimate / alpha;
const double sigma = (1. - std::sqrt(1. / actual_range)) /
(1. + std::sqrt(1. / actual_range));
const double eps = data.smoothing_range;
Expand All @@ -2372,10 +2377,10 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::

const_cast<
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
->delta = (info.max_eigenvalue - alpha) * 0.5;
->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
const_cast<
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
->theta = (info.max_eigenvalue + alpha) * 0.5;
->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;

// We do not need the second temporary vector in case we have a
// DiagonalMatrix as preconditioner and use deal.II's own vectors
Expand Down

0 comments on commit 2664fb2

Please sign in to comment.