Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
eigen.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2020 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 
24 #include <deal.II/lac/solver.h>
25 #include <deal.II/lac/solver_cg.h>
30 
31 #include <cmath>
32 
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 
88 
95  template <typename MatrixType>
96  void
97  solve(double &value, const MatrixType &A, VectorType &x);
98 
99 protected:
104 };
105 
129 template <typename VectorType = Vector<double>>
130 class EigenInverse : private SolverBase<VectorType>
131 {
132 public:
137 
142  {
146  double relaxation;
147 
151  unsigned int start_adaption;
159  AdditionalData(double relaxation = 1.,
160  unsigned int start_adaption = 6,
161  bool use_residual = true)
162  : relaxation(relaxation)
163  , start_adaption(start_adaption)
164  , use_residual(use_residual)
165  {}
166  };
167 
173  const AdditionalData & data = AdditionalData());
174 
182  template <typename MatrixType>
183  void
184  solve(double &value, const MatrixType &A, VectorType &x);
185 
186 protected:
191 };
192 
194 //---------------------------------------------------------------------------
195 
196 
197 template <class VectorType>
200  const AdditionalData & data)
201  : SolverBase<VectorType>(cn, mem)
202  , additional_data(data)
203 {}
204 
205 
206 
207 template <class VectorType>
208 template <typename MatrixType>
209 void
210 EigenPower<VectorType>::solve(double &value, const MatrixType &A, VectorType &x)
211 {
213 
214  LogStream::Prefix prefix("Power method");
215 
216  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
217  VectorType & y = *Vy;
218  y.reinit(x);
219  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
220  VectorType & r = *Vr;
221  r.reinit(x);
222 
223  double length = x.l2_norm();
224  double old_length = 0.;
225  x *= 1. / length;
226 
227  A.vmult(y, x);
228 
229  // Main loop
230  int iter = 0;
231  for (; conv == SolverControl::iterate; iter++)
232  {
233  y.add(additional_data.shift, x);
234 
235  // Compute absolute value of eigenvalue
236  old_length = length;
237  length = y.l2_norm();
238 
239  // do a little trick to compute the sign
240  // with not too much effect of round-off errors.
241  double entry = 0.;
242  size_type i = 0;
243  double thresh = length / x.size();
244  do
245  {
246  Assert(i < x.size(), ExcInternalError());
247  entry = y(i++);
248  }
249  while (std::fabs(entry) < thresh);
250 
251  --i;
252 
253  // Compute unshifted eigenvalue
254  value = (entry * x(i) < 0.) ? -length : length;
255  value -= additional_data.shift;
256 
257  // Update normalized eigenvector
258  x.equ(1 / length, y);
259 
260  // Compute residual
261  A.vmult(y, x);
262 
263  // Check the change of the eigenvalue
264  // Brrr, this is not really a good criterion
265  conv = this->iteration_status(iter,
266  std::fabs(1. / length - 1. / old_length),
267  x);
268  }
269 
270  // in case of failure: throw exception
273  iter, std::fabs(1. / length - 1. / old_length)));
274 
275  // otherwise exit as normal
276 }
277 
278 //---------------------------------------------------------------------------
279 
280 template <class VectorType>
283  const AdditionalData & data)
284  : SolverBase<VectorType>(cn, mem)
285  , additional_data(data)
286 {}
287 
288 
289 
290 template <class VectorType>
291 template <typename MatrixType>
292 void
294  const MatrixType &A,
295  VectorType & x)
296 {
297  LogStream::Prefix prefix("Wielandt");
298 
300 
301  // Prepare matrix for solver
302  auto A_op = linear_operator(A);
303  double current_shift = -value;
304  auto A_s = A_op + current_shift * identity_operator(A_op);
305 
306  // Define solver
307  ReductionControl inner_control(5000, 1.e-16, 1.e-5, false, false);
309  SolverGMRES<VectorType> solver(inner_control, this->memory);
310 
311  // Next step for recomputing the shift
312  unsigned int goal = additional_data.start_adaption;
313 
314  // Auxiliary vector
315  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
316  VectorType & y = *Vy;
317  y.reinit(x);
318  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
319  VectorType & r = *Vr;
320  r.reinit(x);
321 
322  double length = x.l2_norm();
323  double old_value = value;
324 
325  x *= 1. / length;
326 
327  // Main loop
328  double res = -std::numeric_limits<double>::max();
329  size_type iter = 0;
330  for (; conv == SolverControl::iterate; iter++)
331  {
332  solver.solve(A_s, y, x, prec);
333 
334  // Compute absolute value of eigenvalue
335  length = y.l2_norm();
336 
337  // do a little trick to compute the sign
338  // with not too much effect of round-off errors.
339  double entry = 0.;
340  size_type i = 0;
341  double thresh = length / x.size();
342  do
343  {
344  Assert(i < x.size(), ExcInternalError());
345  entry = y(i++);
346  }
347  while (std::fabs(entry) < thresh);
348 
349  --i;
350 
351  // Compute unshifted eigenvalue
352  value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;
353 
354  if (iter == goal)
355  {
356  const auto & relaxation = additional_data.relaxation;
357  const double new_shift =
358  relaxation * (-value) + (1. - relaxation) * current_shift;
359 
360  A_s = A_op + new_shift * identity_operator(A_op);
361  current_shift = new_shift;
362 
363  ++goal;
364  }
365 
366  // Update normalized eigenvector
367  x.equ(1. / length, y);
368  // Compute residual
370  {
371  y.equ(value, x);
372  A.vmult(r, x);
373  r.sadd(-1., value, x);
374  res = r.l2_norm();
375  // Check the residual
376  conv = this->iteration_status(iter, res, x);
377  }
378  else
379  {
380  res = std::fabs(1. / value - 1. / old_value);
381  conv = this->iteration_status(iter, res, x);
382  }
383  old_value = value;
384  }
385 
386  // in case of failure: throw
387  // exception
389  SolverControl::NoConvergence(iter, res));
390  // otherwise exit as normal
391 }
392 
394 
395 #endif
Continue iteration.
AdditionalData additional_data
Definition: eigen.h:103
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:463
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:210
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:293
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1419
AdditionalData additional_data
Definition: eigen.h:190
EigenPower(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:198
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
AdditionalData(const double shift=0.)
Definition: eigen.h:76
EigenInverse(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:281
Expression fabs(const Expression &x)
static const char A
unsigned int start_adaption
Definition: eigen.h:151
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &)
AdditionalData(double relaxation=1., unsigned int start_adaption=6, bool use_residual=true)
Definition: eigen.h:159
static const bool value
T max(const T &t, const MPI_Comm &mpi_communicator)
VectorMemory< VectorType > & memory
Definition: solver.h:412
static ::ExceptionBase & ExcInternalError()