Reference documentation for deal.II version Git e8f99f7e31 2020-07-08 12:50:43 +0200
\(\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 
52 template <typename VectorType = Vector<double>>
53 class EigenPower : private SolverBase<VectorType>
54 {
55 public:
60 
65  {
70  double shift;
74  AdditionalData(const double shift = 0.)
75  : shift(shift)
76  {}
77  };
78 
84  const AdditionalData & data = AdditionalData());
85 
86 
93  template <typename MatrixType>
94  void
95  solve(double &value, const MatrixType &A, VectorType &x);
96 
97 protected:
102 };
103 
125 template <typename VectorType = Vector<double>>
126 class EigenInverse : private SolverBase<VectorType>
127 {
128 public:
133 
138  {
142  double relaxation;
143 
147  unsigned int start_adaption;
155  AdditionalData(double relaxation = 1.,
156  unsigned int start_adaption = 6,
157  bool use_residual = true)
158  : relaxation(relaxation)
159  , start_adaption(start_adaption)
160  , use_residual(use_residual)
161  {}
162  };
163 
169  const AdditionalData & data = AdditionalData());
170 
178  template <typename MatrixType>
179  void
180  solve(double &value, const MatrixType &A, VectorType &x);
181 
182 protected:
187 };
188 
190 //---------------------------------------------------------------------------
191 
192 
193 template <class VectorType>
196  const AdditionalData & data)
197  : SolverBase<VectorType>(cn, mem)
198  , additional_data(data)
199 {}
200 
201 
202 
203 template <class VectorType>
204 template <typename MatrixType>
205 void
206 EigenPower<VectorType>::solve(double &value, const MatrixType &A, VectorType &x)
207 {
209 
210  LogStream::Prefix prefix("Power method");
211 
212  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
213  VectorType & y = *Vy;
214  y.reinit(x);
215  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
216  VectorType & r = *Vr;
217  r.reinit(x);
218 
219  double length = x.l2_norm();
220  double old_length = 0.;
221  x *= 1. / length;
222 
223  A.vmult(y, x);
224 
225  // Main loop
226  int iter = 0;
227  for (; conv == SolverControl::iterate; iter++)
228  {
229  y.add(additional_data.shift, x);
230 
231  // Compute absolute value of eigenvalue
232  old_length = length;
233  length = y.l2_norm();
234 
235  // do a little trick to compute the sign
236  // with not too much effect of round-off errors.
237  double entry = 0.;
238  size_type i = 0;
239  double thresh = length / x.size();
240  do
241  {
242  Assert(i < x.size(), ExcInternalError());
243  entry = y(i++);
244  }
245  while (std::fabs(entry) < thresh);
246 
247  --i;
248 
249  // Compute unshifted eigenvalue
250  value = (entry * x(i) < 0.) ? -length : length;
251  value -= additional_data.shift;
252 
253  // Update normalized eigenvector
254  x.equ(1 / length, y);
255 
256  // Compute residual
257  A.vmult(y, x);
258 
259  // Check the change of the eigenvalue
260  // Brrr, this is not really a good criterion
261  conv = this->iteration_status(iter,
262  std::fabs(1. / length - 1. / old_length),
263  x);
264  }
265 
266  // in case of failure: throw exception
269  iter, std::fabs(1. / length - 1. / old_length)));
270 
271  // otherwise exit as normal
272 }
273 
274 //---------------------------------------------------------------------------
275 
276 template <class VectorType>
279  const AdditionalData & data)
280  : SolverBase<VectorType>(cn, mem)
281  , additional_data(data)
282 {}
283 
284 
285 
286 template <class VectorType>
287 template <typename MatrixType>
288 void
290  const MatrixType &A,
291  VectorType & x)
292 {
293  LogStream::Prefix prefix("Wielandt");
294 
296 
297  // Prepare matrix for solver
298  auto A_op = linear_operator(A);
299  double current_shift = -value;
300  auto A_s = A_op + current_shift * identity_operator(A_op);
301 
302  // Define solver
303  ReductionControl inner_control(5000, 1.e-16, 1.e-5, false, false);
305  SolverGMRES<VectorType> solver(inner_control, this->memory);
306 
307  // Next step for recomputing the shift
308  unsigned int goal = additional_data.start_adaption;
309 
310  // Auxiliary vector
311  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
312  VectorType & y = *Vy;
313  y.reinit(x);
314  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
315  VectorType & r = *Vr;
316  r.reinit(x);
317 
318  double length = x.l2_norm();
319  double old_value = value;
320 
321  x *= 1. / length;
322 
323  // Main loop
324  double res = -std::numeric_limits<double>::max();
325  size_type iter = 0;
326  for (; conv == SolverControl::iterate; iter++)
327  {
328  solver.solve(A_s, y, x, prec);
329 
330  // Compute absolute value of eigenvalue
331  length = y.l2_norm();
332 
333  // do a little trick to compute the sign
334  // with not too much effect of round-off errors.
335  double entry = 0.;
336  size_type i = 0;
337  double thresh = length / x.size();
338  do
339  {
340  Assert(i < x.size(), ExcInternalError());
341  entry = y(i++);
342  }
343  while (std::fabs(entry) < thresh);
344 
345  --i;
346 
347  // Compute unshifted eigenvalue
348  value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;
349 
350  if (iter == goal)
351  {
352  const auto & relaxation = additional_data.relaxation;
353  const double new_shift =
354  relaxation * (-value) + (1. - relaxation) * current_shift;
355 
356  A_s = A_op + new_shift * identity_operator(A_op);
357  current_shift = new_shift;
358 
359  ++goal;
360  }
361 
362  // Update normalized eigenvector
363  x.equ(1. / length, y);
364  // Compute residual
366  {
367  y.equ(value, x);
368  A.vmult(r, x);
369  r.sadd(-1., value, x);
370  res = r.l2_norm();
371  // Check the residual
372  conv = this->iteration_status(iter, res, x);
373  }
374  else
375  {
376  res = std::fabs(1. / value - 1. / old_value);
377  conv = this->iteration_status(iter, res, x);
378  }
379  old_value = value;
380  }
381 
382  // in case of failure: throw
383  // exception
385  SolverControl::NoConvergence(iter, res));
386  // otherwise exit as normal
387 }
388 
390 
391 #endif
Continue iteration.
AdditionalData additional_data
Definition: eigen.h:101
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:471
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:206
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:289
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1403
AdditionalData additional_data
Definition: eigen.h:186
EigenPower(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:194
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
AdditionalData(const double shift=0.)
Definition: eigen.h:74
EigenInverse(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:277
Expression fabs(const Expression &x)
static const char A
unsigned int start_adaption
Definition: eigen.h:147
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:155
T max(const T &t, const MPI_Comm &mpi_communicator)
VectorMemory< VectorType > & memory
Definition: solver.h:420
static ::ExceptionBase & ExcInternalError()