From 02baeecf521dabaf21f53e30e679affc04548fa1 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 23 Apr 2003 17:25:45 +0000 Subject: [PATCH] Wielandt iteration improved git-svn-id: https://svn.dealii.org/trunk@7460 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/eigen.h | 61 +++++++++++++++++++++++++-------- 1 file changed, 46 insertions(+), 15 deletions(-) diff --git a/deal.II/lac/include/lac/eigen.h b/deal.II/lac/include/lac/eigen.h index 85db6bbd77..ab474523bd 100644 --- a/deal.II/lac/include/lac/eigen.h +++ b/deal.II/lac/include/lac/eigen.h @@ -19,6 +19,8 @@ #include #include #include +#include +#include #include #include @@ -111,9 +113,20 @@ class EigenPower : private Solver * converge to zero for non-symmetric matrices with non-trivial Jordan * blocks, it can be replaced by checking the difference of successive * eigenvalues. Use @p{AdditionalData::use_residual} for switching - * these options. + * this option. * - * @author Guido Kanschat, 2000 + * Usually, the initial guess entering this method is updated after + * each step, replacing it with the new approximation of the + * eigenvalue. Using a parameter @p{AdditionalData::relaxation} + * between 0 and 1, this update can be damped. With relaxation + * parameter 0, no update is performed. This damping allows for slower + * adaption of the shift value to make sure that the method converges + * to the eigenvalue closest to the initial guess. This can be aided + * by the parameter @p{AdditionalData::start_adaption}, which + * indicates the first iteration step in which the shift value should + * be adapted. + * + * @author Guido Kanschat, 2000, 2003 */ template > class EigenInverse : private Solver @@ -126,16 +139,30 @@ class EigenInverse : private Solver */ struct AdditionalData { + /** + * Damping of the updated shift value. + */ + double relaxation; + + /** + * Start step of adaptive + * shift parameter. + */ + unsigned int start_adaption; /** * Flag for the stopping criterion. */ - bool use_residual; + bool use_residual; /** * Constructor. */ - AdditionalData (bool use_residual = true): - use_residual(use_residual) - {} + AdditionalData (double relaxation = 1., + unsigned int start_adaption = 6, + bool use_residual = true): + relaxation(relaxation), + start_adaption(start_adaption), + use_residual(use_residual) + {} }; @@ -173,7 +200,7 @@ class EigenInverse : private Solver protected: /** - * Shift parameter. + * Flags for execution. */ AdditionalData additional_data; }; @@ -270,8 +297,8 @@ EigenPower::solve (double &value, template EigenInverse::EigenInverse (SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data): + VectorMemory &mem, + const AdditionalData &data): Solver(cn, mem), additional_data(data) {} @@ -289,7 +316,7 @@ EigenInverse::solve (double &value, const MATRIX &A, VECTOR &x) { - deallog.push("Wieland"); + deallog.push("Wielandt"); SolverControl::State conv=SolverControl::iterate; @@ -297,13 +324,13 @@ EigenInverse::solve (double &value, ShiftedMatrix A_s(A, -value); // Define solver - ReductionControl inner_control (A.m(), 1.e-16, 1.e-8, false, false); + ReductionControl inner_control (5000, 1.e-16, 1.e-5, false, false); PreconditionIdentity prec; - SolverCG + SolverGMRES solver(inner_control, this->memory); // Next step for recomputing the shift - unsigned int goal = 10; + unsigned int goal = additional_data.start_adaption; // Auxiliary vector VECTOR* Vy = this->memory.alloc (); VECTOR& y = *Vy; y.reinit (x); @@ -345,7 +372,9 @@ EigenInverse::solve (double &value, if (iter==goal) { - A_s.shift(-value); + const double new_shift = - additional_data.relaxation * value + + (1.-additional_data.relaxation) * A_s.shift(); + A_s.shift(new_shift); ++goal; } @@ -355,7 +384,9 @@ EigenInverse::solve (double &value, if (additional_data.use_residual) { y.equ (value, x); - double res = A.residual (r,x,y); + A.vmult(r,x); + r.sadd(-1., value, x); + double res = r.l2_norm(); // Check the residual conv = this->control().check (iter, res); } else { -- 2.39.5