From 498b476ae0c3caf945fe85c38fe7d06c0b804596 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Thu, 20 Apr 2000 17:11:39 +0000 Subject: [PATCH] inverse iteration git-svn-id: https://svn.dealii.org/trunk@2763 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/eigen.h | 293 ++++++++++++++++++++++++- deal.II/lac/include/lac/precondition.h | 95 ++++++++ 2 files changed, 381 insertions(+), 7 deletions(-) diff --git a/deal.II/lac/include/lac/eigen.h b/deal.II/lac/include/lac/eigen.h index 7ccb12a908..e12dcdea72 100644 --- a/deal.II/lac/include/lac/eigen.h +++ b/deal.II/lac/include/lac/eigen.h @@ -17,7 +17,62 @@ #include #include #include +#include #include +#include + +/** + * Matrix with shifted diagonal values. + * + * Given a matrix $A$, this class implements a matrix-vector product with + * $A+\sigma I$, where sigma is a provided shift parameter. + * + * @author Guido Kanschat, 2000 + */ +template +class ShiftedMatrix +{ + public: + /** + * Constructor. + * Provide the base matrix and a shift parameter. + */ + ShiftedMatrix (const MATRIX& A, const double sigma); + + /** + * Set the shift parameter. + */ + void shift (const double sigma); + + /** + * Access to the shift parameter. + */ + double shift () const; + + /** + * Matrix-vector-product. + */ + template + void vmult (VECTOR& dst, const VECTOR& src) const; + + /** + * Residual. + */ + template + double residual (VECTOR& dst, const VECTOR& src, const VECTOR& rhs) const; + + private: + /** + * Storage for base matrix. + */ + const MATRIX& A; + + /** + * Shift parameter. + */ + double sigma; +}; + /** * Power method (von Mises). @@ -36,14 +91,13 @@ */ template , class VECTOR = Vector > -class EigenPower : public Solver +class EigenPower : private Solver { public: /** * Standardized data struct to * pipe additional data to the - * solver. This solver does not - * need additional data yet. + * solver. */ struct AdditionalData { @@ -57,7 +111,7 @@ class EigenPower : public Solver /** * Constructor. Set the shift parameter. */ - AdditionalData (const double shift): + AdditionalData (const double shift = 0.): shift(shift) {} @@ -97,8 +151,122 @@ class EigenPower : public Solver AdditionalData additional_data; }; +//TODO: Finish documentation after fixing parameters. + +/** + * Inverse iteration (Wieland). + * + * This class implements an adaptive version of the inverse iteration by Wieland. + * + * @author Guido Kanschat, 2000 + */ +template , + class VECTOR = Vector > +class EigenInverse : private Solver +{ + public: + /** + * Standardized data struct to + * pipe additional data to the + * solver. This solver does not + * need additional data yet. + */ + struct AdditionalData + {}; + + /** + * Constructor. + */ + EigenInverse (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data=AdditionalData()); + + + /** + * Virtual destructor. + */ + virtual ~EigenInverse (); + + /** + * Inverse method. @p value is + * the start guess for the + * eigenvalue and @p x is the + * (not necessarily normalized) + * start vector for the power + * method. After the iteration, + * @p value is the approximated + * eigenvalue and @p x is the + * corresponding eigenvector, + * normalized with respect to the + * l2-norm. + */ + typename Solver::ReturnState + solve (double &value, + const MATRIX &A, + VECTOR &x); + + protected: + /** + * Shift parameter. + */ + AdditionalData additional_data; +}; + //----------------------------------------------------------------------// +template +inline +ShiftedMatrix::ShiftedMatrix (const MATRIX& A, const double sigma) + : + A(A), sigma(sigma) +{} + + + +template +inline void +ShiftedMatrix::shift (const double s) +{ + sigma = s; +} + + +template +inline double +ShiftedMatrix::shift () const +{ + return sigma; +} + + + +template +template +inline void +ShiftedMatrix::vmult (VECTOR& dst, const VECTOR& src) const +{ + A.vmult(dst, src); + dst.add(sigma, src); +} + + +template +template +inline double +ShiftedMatrix::residual (VECTOR& dst, + const VECTOR& src, + const VECTOR& rhs) const +{ + A.vmult(dst, src); + dst.add(sigma, src); + dst.sadd(-1.,1.,rhs); + return dst.l2_norm (); +} + + +//---------------------------------------------------------------------- + + template EigenPower::EigenPower (SolverControl &cn, VectorMemory &mem, @@ -124,16 +292,17 @@ EigenPower::solve (double &value, deallog.push("Power method"); VECTOR* Vy = memory.alloc (); VECTOR& y = *Vy; y.reinit (x); + VECTOR* Vr = memory.alloc (); VECTOR& r = *Vr; r.reinit (x); double length = x.l2_norm (); double old_length = 0.; x.scale(1./length); + A.vmult (y,x); // Main loop for(int iter=0; conv==SolverControl::iterate; iter++) { - A.vmult (y,x); y.add(additional_data.shift, x); // Compute absolute value of eigenvalue @@ -141,7 +310,7 @@ EigenPower::solve (double &value, length = y.l2_norm (); // do a little trick to compute the sign - // with not too much round-off errors. + // with not too much effect of round-off errors. double entry = 0.; unsigned int i = 0; double thresh = length/x.size(); @@ -161,6 +330,11 @@ EigenPower::solve (double &value, // Update normalized eigenvector x.equ (1/length, y); + // Compute residual + A.vmult (y,x); + + // Check the change of the eigenvalue + // Brrr, this is not really a good criterion conv = control().check (iter, fabs(1.-length/old_length)); } @@ -172,8 +346,113 @@ EigenPower::solve (double &value, return exceeded; else return success; - } +//----------------------------------------------------------------------// + +template +EigenInverse::EigenInverse (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data): + Solver(cn, mem), + additional_data(data) +{} + + +template +EigenInverse::~EigenInverse () +{} + + +template +typename Solver::ReturnState +EigenInverse::solve (double &value, + const MATRIX &A, + VECTOR &x) +{ + deallog.push("Wieland"); + + SolverControl::State conv=SolverControl::iterate; + + // Prepare matrix for solver + ShiftedMatrix A_s(A, -value); + + // Define solver + ReductionControl inner_control (100, 1.e-16, 1.e-8, false, false); + PreconditionIdentity prec; + SolverQMRS, VECTOR> + solver(inner_control, memory); + + // Next step for recomputing the shift + unsigned int goal = 10; + + // Auxiliary vector + VECTOR* Vy = memory.alloc (); VECTOR& y = *Vy; y.reinit (x); + VECTOR* Vr = memory.alloc (); VECTOR& r = *Vr; r.reinit (x); + + double length = x.l2_norm (); + double old_length = 0.; + double old_value = value; + + x.scale(1./length); + + // Main loop + for (unsigned int iter=0; conv==SolverControl::iterate; iter++) + { + solver.solve (A_s, y, x, prec); + + // Compute absolute value of eigenvalue + old_length = length; + length = y.l2_norm (); + + // do a little trick to compute the sign + // with not too much effect of round-off errors. + double entry = 0.; + unsigned int i = 0; + double thresh = length/x.size(); + do + { + Assert (i /** * No preconditioning. @@ -254,6 +255,59 @@ class PreconditionLACSolver }; +/** + * Matrix with preconditioner. + * Given a matrix $A$ and a preconditioner $P$, this class implements a new matrix + * with the matrix-vector product $PA$. It needs an auxiliary vector for that. + * + * By this time, this is considered a temporary object to be plugged + * into eigenvalue solvers. Therefore, no @p SmartPointer is used for + * @p A and @p P. + * + * @author Guido Kanschat, 2000 + */ +template +class PreconditionedMatrix +{ + public: + /** + * Constructor. Provide matrix, + * preconditioner and a memory + * pool to obtain the auxiliary + * vector. + */ + PreconditionedMatrix (const MATRIX& A, + const PRECOND& P, + VectorMemory& mem); + + /** + * Preconditioned matrix-vector-product. + */ + void vmult (VECTOR& dst, const VECTOR& src) const; + + /** + * Residual $b-PAx$. + */ + double residual (VECTOR& dst, const VECTOR& src, const VECTOR& rhs) const; + + private: + /** + * Storage for the matrix. + */ + const MATRIX& A; + /** + * Storage for preconditioner. + */ + const PRECOND& P; + /** + * Memory pool for vectors. + */ + VectorMemory& mem; +}; + + + + /* ---------------------------------- Inline functions ------------------- */ template @@ -316,5 +370,46 @@ PreconditionLACSolver::operator() (VECTOR& dst, solver.solve(matrix, dst, src, precondition); } +////////////////////////////////////////////////////////////////////// + + +template +inline +PreconditionedMatrix +::PreconditionedMatrix (const MATRIX& A, + const PRECOND& P, + VectorMemory& mem): + A(A), P(P), mem(mem) +{} + + +template +inline void +PreconditionedMatrix +::vmult (VECTOR& dst, + const VECTOR& src) const +{ + VECTOR* h = mem.alloc(); + h->reinit(src); + A.vmult(*h, src); + P(dst, *h); + mem.free(h); +} + +template +inline double +PreconditionedMatrix +::residual (VECTOR& dst, + const VECTOR& src, + const VECTOR& rhs) const +{ + VECTOR* h = mem.alloc(); + h->reinit(src); + A.vmult(*h, src); + P(dst, *h); + mem.free(h); + dst.sadd(-1.,1.,rhs); + return dst.l2_norm (); +} #endif -- 2.39.5