From 9e0bf81a64c1a9c5473c16ebf6782d41721666f4 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 18 Apr 2000 21:57:33 +0000 Subject: [PATCH] power method by von Mises git-svn-id: https://svn.dealii.org/trunk@2745 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/eigen.h | 179 ++++++++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 deal.II/lac/include/lac/eigen.h diff --git a/deal.II/lac/include/lac/eigen.h b/deal.II/lac/include/lac/eigen.h new file mode 100644 index 0000000000..7ccb12a908 --- /dev/null +++ b/deal.II/lac/include/lac/eigen.h @@ -0,0 +1,179 @@ +//---------------------------- eigen.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- eigen.h --------------------------- +#ifndef __deal2__eigen_h +#define __deal2__eigen_h + + +#include +#include +#include +#include + +/** + * Power method (von Mises). + * + * This method determines the largest eigenvalue of a matrix by + * applying increasing powers of this matrix to a vector. If there is + * an eigenvalue $l$ with dominant absolute value, the iteration vectors + * will become aligned to its eigenspace and $Ax = lx$. + * + * A shift parameter allows to shift the spectrum, so it is possible + * to compute the smallest eigenvalue, too. + * + * Convergence of this method is known to be slow. + * + * @author Guido Kanschat, 2000 + */ +template , + class VECTOR = Vector > +class EigenPower : public Solver +{ + public: + /** + * Standardized data struct to + * pipe additional data to the + * solver. This solver does not + * need additional data yet. + */ + struct AdditionalData + { + /** + * Shift parameter. This + * parameter allows to shift + * the spectrum to compute a + * different eigenvalue. + */ + double shift; + /** + * Constructor. Set the shift parameter. + */ + AdditionalData (const double shift): + shift(shift) + {} + + }; + + /** + * Constructor. + */ + EigenPower (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data=AdditionalData()); + + /** + * Virtual destructor. + */ + virtual ~EigenPower (); + + /** + * Power method. @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 +EigenPower::EigenPower (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data): + Solver(cn, mem), + additional_data(data) +{} + + +template +EigenPower::~EigenPower () +{} + + +template +typename Solver::ReturnState +EigenPower::solve (double &value, + const MATRIX &A, + VECTOR &x) +{ + SolverControl::State conv=SolverControl::iterate; + + deallog.push("Power method"); + + VECTOR* Vy = memory.alloc (); VECTOR& y = *Vy; y.reinit (x); + + double length = x.l2_norm (); + double old_length = 0.; + x.scale(1./length); + + + // 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 + old_length = length; + length = y.l2_norm (); + + // do a little trick to compute the sign + // with not too much round-off errors. + double entry = 0.; + unsigned int i = 0; + double thresh = length/x.size(); + do + { + Assert (i