Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
eigen.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_eigen_h
17 #define dealii_eigen_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/lac/linear_operator.h>
23 #include <deal.II/lac/precondition.h>
24 #include <deal.II/lac/solver.h>
25 #include <deal.II/lac/solver_cg.h>
26 #include <deal.II/lac/solver_control.h>
27 #include <deal.II/lac/solver_gmres.h>
28 #include <deal.II/lac/solver_minres.h>
29 #include <deal.II/lac/vector_memory.h>
30 
31 #include <cmath>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 
38 
54 template <typename VectorType = Vector<double>>
55 class EigenPower : private SolverBase<VectorType>
56 {
57 public:
62 
67  {
72  double shift;
76  AdditionalData(const double shift = 0.)
77  : shift(shift)
78  {}
79  };
80 
86  const AdditionalData & data = AdditionalData());
87 
91  virtual ~EigenPower();
92 
99  template <typename MatrixType>
100  void
101  solve(double &value, const MatrixType &A, VectorType &x);
102 
103 protected:
108 };
109 
133 template <typename VectorType = Vector<double>>
134 class EigenInverse : private SolverBase<VectorType>
135 {
136 public:
141 
146  {
150  double relaxation;
151 
155  unsigned int start_adaption;
164  unsigned int start_adaption = 6,
165  bool use_residual = true)
169  {}
170  };
171 
177  const AdditionalData & data = AdditionalData());
178 
179 
183  virtual ~EigenInverse();
184 
192  template <typename MatrixType>
193  void
194  solve(double &value, const MatrixType &A, VectorType &x);
195 
196 protected:
201 };
202 
204 //---------------------------------------------------------------------------
205 
206 
207 template <class VectorType>
210  const AdditionalData & data)
211  : SolverBase<VectorType>(cn, mem)
212  , additional_data(data)
213 {}
214 
215 
216 
217 template <class VectorType>
219 {}
220 
221 
222 
223 template <class VectorType>
224 template <typename MatrixType>
225 void
226 EigenPower<VectorType>::solve(double &value, const MatrixType &A, VectorType &x)
227 {
229 
230  LogStream::Prefix prefix("Power method");
231 
232  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
233  VectorType & y = *Vy;
234  y.reinit(x);
235  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
236  VectorType & r = *Vr;
237  r.reinit(x);
238 
239  double length = x.l2_norm();
240  double old_length = 0.;
241  x *= 1. / length;
242 
243  A.vmult(y, x);
244 
245  // Main loop
246  int iter = 0;
247  for (; conv == SolverControl::iterate; iter++)
248  {
249  y.add(additional_data.shift, x);
250 
251  // Compute absolute value of eigenvalue
252  old_length = length;
253  length = y.l2_norm();
254 
255  // do a little trick to compute the sign
256  // with not too much effect of round-off errors.
257  double entry = 0.;
258  size_type i = 0;
259  double thresh = length / x.size();
260  do
261  {
262  Assert(i < x.size(), ExcInternalError());
263  entry = y(i++);
264  }
265  while (std::fabs(entry) < thresh);
266 
267  --i;
268 
269  // Compute unshifted eigenvalue
270  value = (entry * x(i) < 0.) ? -length : length;
271  value -= additional_data.shift;
272 
273  // Update normalized eigenvector
274  x.equ(1 / length, y);
275 
276  // Compute residual
277  A.vmult(y, x);
278 
279  // Check the change of the eigenvalue
280  // Brrr, this is not really a good criterion
281  conv = this->iteration_status(iter,
282  std::fabs(1. / length - 1. / old_length),
283  x);
284  }
285 
286  // in case of failure: throw exception
289  iter, std::fabs(1. / length - 1. / old_length)));
290 
291  // otherwise exit as normal
292 }
293 
294 //---------------------------------------------------------------------------
295 
296 template <class VectorType>
299  const AdditionalData & data)
300  : SolverBase<VectorType>(cn, mem)
301  , additional_data(data)
302 {}
303 
304 
305 
306 template <class VectorType>
308 {}
309 
310 
311 
312 template <class VectorType>
313 template <typename MatrixType>
314 void
316  const MatrixType &A,
317  VectorType & x)
318 {
319  LogStream::Prefix prefix("Wielandt");
320 
322 
323  // Prepare matrix for solver
324  auto A_op = linear_operator(A);
325  double current_shift = -value;
326  auto A_s = A_op + current_shift * identity_operator(A_op);
327 
328  // Define solver
329  ReductionControl inner_control(5000, 1.e-16, 1.e-5, false, false);
331  SolverGMRES<VectorType> solver(inner_control, this->memory);
332 
333  // Next step for recomputing the shift
334  unsigned int goal = additional_data.start_adaption;
335 
336  // Auxiliary vector
337  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
338  VectorType & y = *Vy;
339  y.reinit(x);
340  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
341  VectorType & r = *Vr;
342  r.reinit(x);
343 
344  double length = x.l2_norm();
345  double old_value = value;
346 
347  x *= 1. / length;
348 
349  // Main loop
350  double res = -std::numeric_limits<double>::max();
351  size_type iter = 0;
352  for (; conv == SolverControl::iterate; iter++)
353  {
354  solver.solve(A_s, y, x, prec);
355 
356  // Compute absolute value of eigenvalue
357  length = y.l2_norm();
358 
359  // do a little trick to compute the sign
360  // with not too much effect of round-off errors.
361  double entry = 0.;
362  size_type i = 0;
363  double thresh = length / x.size();
364  do
365  {
366  Assert(i < x.size(), ExcInternalError());
367  entry = y(i++);
368  }
369  while (std::fabs(entry) < thresh);
370 
371  --i;
372 
373  // Compute unshifted eigenvalue
374  value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;
375 
376  if (iter == goal)
377  {
378  const auto & relaxation = additional_data.relaxation;
379  const double new_shift =
380  relaxation * (-value) + (1. - relaxation) * current_shift;
381 
382  A_s = A_op + new_shift * identity_operator(A_op);
383  current_shift = new_shift;
384 
385  ++goal;
386  }
387 
388  // Update normalized eigenvector
389  x.equ(1. / length, y);
390  // Compute residual
391  if (additional_data.use_residual)
392  {
393  y.equ(value, x);
394  A.vmult(r, x);
395  r.sadd(-1., value, x);
396  res = r.l2_norm();
397  // Check the residual
398  conv = this->iteration_status(iter, res, x);
399  }
400  else
401  {
402  res = std::fabs(1. / value - 1. / old_value);
403  conv = this->iteration_status(iter, res, x);
404  }
405  old_value = value;
406  }
407 
408  // in case of failure: throw
409  // exception
411  SolverControl::NoConvergence(iter, res));
412  // otherwise exit as normal
413 }
414 
415 DEAL_II_NAMESPACE_CLOSE
416 
417 #endif
types::global_dof_index size_type
Definition: eigen.h:61
Continue iteration.
AdditionalData additional_data
Definition: eigen.h:107
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:226
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:315
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1407
AdditionalData additional_data
Definition: eigen.h:200
EigenPower(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:208
AdditionalData(const double shift=0.)
Definition: eigen.h:76
EigenInverse(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:297
virtual ~EigenInverse()
Definition: eigen.h:307
virtual ~EigenPower()
Definition: eigen.h:218
unsigned int start_adaption
Definition: eigen.h:155
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
unsigned int global_dof_index
Definition: types.h:89
types::global_dof_index size_type
Definition: eigen.h:140
AdditionalData(double relaxation=1., unsigned int start_adaption=6, bool use_residual=true)
Definition: eigen.h:163
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
static ::ExceptionBase & ExcInternalError()