Reference documentation for deal.II version 8.1.0
eigen.h
1 // ---------------------------------------------------------------------
2 // @f$Id: eigen.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 2000 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__eigen_h
18 #define __deal2__eigen_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/lac/shifted_matrix.h>
23 #include <deal.II/lac/solver.h>
24 #include <deal.II/lac/solver_control.h>
25 #include <deal.II/lac/solver_cg.h>
26 #include <deal.II/lac/solver_gmres.h>
27 #include <deal.II/lac/solver_minres.h>
28 #include <deal.II/lac/vector_memory.h>
29 #include <deal.II/lac/precondition.h>
30 
31 #include <cmath>
32 
34 
35 
38 
54 template <class VECTOR = Vector<double> >
55 class EigenPower : private Solver<VECTOR>
56 {
57 public:
62 
69  {
76  double shift;
80  AdditionalData (const double shift = 0.)
81  :
82  shift(shift)
83  {}
84 
85  };
86 
92  const AdditionalData &data=AdditionalData());
93 
97  virtual ~EigenPower ();
98 
110  template <class MATRIX>
111  void
112  solve (double &value,
113  const MATRIX &A,
114  VECTOR &x);
115 
116 protected:
121 };
122 
148 template <class VECTOR = Vector<double> >
149 class EigenInverse : private Solver<VECTOR>
150 {
151 public:
156 
163  {
167  double relaxation;
168 
173  unsigned int start_adaption;
181  AdditionalData (double relaxation = 1.,
182  unsigned int start_adaption = 6,
183  bool use_residual = true)
184  :
185  relaxation(relaxation),
186  start_adaption(start_adaption),
187  use_residual(use_residual)
188  {}
189 
190  };
191 
197  const AdditionalData &data=AdditionalData());
198 
199 
203  virtual ~EigenInverse ();
204 
218  template <class MATRIX>
219  void
220  solve (double &value,
221  const MATRIX &A,
222  VECTOR &x);
223 
224 protected:
229 };
230 
232 //---------------------------------------------------------------------------
233 
234 
235 template <class VECTOR>
238  const AdditionalData &data)
239  :
240  Solver<VECTOR>(cn, mem),
241  additional_data(data)
242 {}
243 
244 
245 
246 template <class VECTOR>
248 {}
249 
250 
251 
252 template <class VECTOR>
253 template <class MATRIX>
254 void
256  const MATRIX &A,
257  VECTOR &x)
258 {
260 
261  deallog.push("Power method");
262 
263  VECTOR *Vy = this->memory.alloc ();
264  VECTOR &y = *Vy;
265  y.reinit (x);
266  VECTOR *Vr = this->memory.alloc ();
267  VECTOR &r = *Vr;
268  r.reinit (x);
269 
270  double length = x.l2_norm ();
271  double old_length = 0.;
272  x.scale(1./length);
273 
274  A.vmult (y,x);
275 
276  // Main loop
277  for (int iter=0; conv==SolverControl::iterate; iter++)
278  {
279  y.add(additional_data.shift, x);
280 
281  // Compute absolute value of eigenvalue
282  old_length = length;
283  length = y.l2_norm ();
284 
285  // do a little trick to compute the sign
286  // with not too much effect of round-off errors.
287  double entry = 0.;
288  size_type i = 0;
289  double thresh = length/x.size();
290  do
291  {
292  Assert (i<x.size(), ExcInternalError());
293  entry = y (i++);
294  }
295  while (std::fabs(entry) < thresh);
296 
297  --i;
298 
299  // Compute unshifted eigenvalue
300  value = (entry * x (i) < 0.) ? -length : length;
301  value -= additional_data.shift;
302 
303  // Update normalized eigenvector
304  x.equ (1/length, y);
305 
306  // Compute residual
307  A.vmult (y,x);
308 
309  // Check the change of the eigenvalue
310  // Brrr, this is not really a good criterion
311  conv = this->control().check (iter, std::fabs(1./length-1./old_length));
312  }
313 
314  this->memory.free(Vy);
315  this->memory.free(Vr);
316 
317  deallog.pop();
318 
319  // in case of failure: throw exception
320  if (this->control().last_check() != SolverControl::success)
321  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
322  this->control().last_value()));
323  // otherwise exit as normal
324 }
325 
326 //---------------------------------------------------------------------------
327 
328 template <class VECTOR>
331  const AdditionalData &data)
332  :
333  Solver<VECTOR>(cn, mem),
334  additional_data(data)
335 {}
336 
337 
338 
339 template <class VECTOR>
341 {}
342 
343 
344 
345 template <class VECTOR>
346 template <class MATRIX>
347 void
349  const MATRIX &A,
350  VECTOR &x)
351 {
352  deallog.push("Wielandt");
353 
355 
356  // Prepare matrix for solver
357  ShiftedMatrix <MATRIX> A_s(A, -value);
358 
359  // Define solver
360  ReductionControl inner_control (5000, 1.e-16, 1.e-5, false, false);
363  solver(inner_control, this->memory);
364 
365  // Next step for recomputing the shift
366  unsigned int goal = additional_data.start_adaption;
367 
368  // Auxiliary vector
369  VECTOR *Vy = this->memory.alloc ();
370  VECTOR &y = *Vy;
371  y.reinit (x);
372  VECTOR *Vr = this->memory.alloc ();
373  VECTOR &r = *Vr;
374  r.reinit (x);
375 
376  double length = x.l2_norm ();
377  double old_value = value;
378 
379  x.scale(1./length);
380 
381  // Main loop
382  for (size_type iter=0; conv==SolverControl::iterate; iter++)
383  {
384  solver.solve (A_s, y, x, prec);
385 
386  // Compute absolute value of eigenvalue
387  length = y.l2_norm ();
388 
389  // do a little trick to compute the sign
390  // with not too much effect of round-off errors.
391  double entry = 0.;
392  size_type i = 0;
393  double thresh = length/x.size();
394  do
395  {
396  Assert (i<x.size(), ExcInternalError());
397  entry = y (i++);
398  }
399  while (std::fabs(entry) < thresh);
400 
401  --i;
402 
403  // Compute unshifted eigenvalue
404  value = (entry * x (i) < 0.) ? -length : length;
405  value = 1./value;
406  value -= A_s.shift ();
407 
408  if (iter==goal)
409  {
410  const double new_shift = - additional_data.relaxation * value
411  + (1.-additional_data.relaxation) * A_s.shift();
412  A_s.shift(new_shift);
413  ++goal;
414  }
415 
416  // Update normalized eigenvector
417  x.equ (1./length, y);
418  // Compute residual
420  {
421  y.equ (value, x);
422  A.vmult(r,x);
423  r.sadd(-1., value, x);
424  double res = r.l2_norm();
425  // Check the residual
426  conv = this->control().check (iter, res);
427  }
428  else
429  {
430  conv = this->control().check (iter, std::fabs(1./value-1./old_value));
431  }
432  old_value = value;
433  }
434 
435  this->memory.free(Vy);
436  this->memory.free(Vr);
437 
438  deallog.pop();
439 
440  // in case of failure: throw
441  // exception
442  if (this->control().last_check() != SolverControl::success)
443  throw SolverControl::NoConvergence (this->control().last_step(),
444  this->control().last_value());
445  // otherwise exit as normal
446 }
447 
448 DEAL_II_NAMESPACE_CLOSE
449 
450 #endif
EigenInverse(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:329
VectorMemory< VECTOR > & memory
Definition: solver.h:210
void solve(double &value, const MATRIX &A, VECTOR &x)
Definition: eigen.h:255
void pop()
virtual State check(const unsigned int step, const double check_value)
void vmult(VECTOR &u, const VECTOR &v) const
Continue iteration.
virtual ~EigenPower()
Definition: eigen.h:247
virtual ~EigenInverse()
Definition: eigen.h:340
types::global_dof_index size_type
Definition: eigen.h:61
void shift(const double sigma)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
AdditionalData additional_data
Definition: eigen.h:120
SolverControl & control() const
Definition: solver.h:239
unsigned int global_dof_index
Definition: types.h:100
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:299
AdditionalData additional_data
Definition: eigen.h:228
AdditionalData(const double shift=0.)
Definition: eigen.h:80
void solve(double &value, const MATRIX &A, VECTOR &x)
Definition: eigen.h:348
unsigned int start_adaption
Definition: eigen.h:173
Definition: solver.h:147
void push(const std::string &text)
types::global_dof_index size_type
Definition: eigen.h:155
AdditionalData(double relaxation=1., unsigned int start_adaption=6, bool use_residual=true)
Definition: eigen.h:181
::ExceptionBase & ExcInternalError()
EigenPower(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:236