Reference documentation for deal.II version GIT a189bc2bdf 2022-12-07 02:45:02+00:00
\(\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 - 2021 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 #include <limits>
33 
35 
36 
55 template <typename VectorType = Vector<double>>
56 class EigenPower : private SolverBase<VectorType>
57 {
58 public:
63 
68  {
73  double shift;
77  AdditionalData(const double shift = 0.)
78  : shift(shift)
79  {}
80  };
81 
87  const AdditionalData & data = AdditionalData());
88 
89 
96  template <typename MatrixType>
97  void
98  solve(double &value, const MatrixType &A, VectorType &x);
99 
100 protected:
105 };
106 
128 template <typename VectorType = Vector<double>>
129 class EigenInverse : private SolverBase<VectorType>
130 {
131 public:
136 
141  {
145  double relaxation;
146 
150  unsigned int start_adaption;
159  unsigned int start_adaption = 6,
160  bool use_residual = true)
164  {}
165  };
166 
172  const AdditionalData & data = AdditionalData());
173 
181  template <typename MatrixType>
182  void
183  solve(double &value, const MatrixType &A, VectorType &x);
184 
185 protected:
190 };
191 
193 //---------------------------------------------------------------------------
194 
195 
196 template <class VectorType>
199  const AdditionalData & data)
200  : SolverBase<VectorType>(cn, mem)
201  , additional_data(data)
202 {}
203 
204 
205 
206 template <class VectorType>
207 template <typename MatrixType>
208 void
209 EigenPower<VectorType>::solve(double &value, const MatrixType &A, VectorType &x)
210 {
212 
213  LogStream::Prefix prefix("Power method");
214 
215  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
216  VectorType & y = *Vy;
217  y.reinit(x);
218  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
219  VectorType & r = *Vr;
220  r.reinit(x);
221 
222  double length = x.l2_norm();
223  double old_length = 0.;
224  x *= 1. / length;
225 
226  A.vmult(y, x);
227 
228  // Main loop
229  int iter = 0;
230  for (; conv == SolverControl::iterate; ++iter)
231  {
232  y.add(additional_data.shift, x);
233 
234  // Compute absolute value of eigenvalue
235  old_length = length;
236  length = y.l2_norm();
237 
238  // do a little trick to compute the sign
239  // with not too much effect of round-off errors.
240  double entry = 0.;
241  size_type i = 0;
242  double thresh = length / x.size();
243  do
244  {
245  Assert(i < x.size(), ExcInternalError());
246  entry = y(i++);
247  }
248  while (std::fabs(entry) < thresh);
249 
250  --i;
251 
252  // Compute unshifted eigenvalue
253  value = (entry * x(i) < 0.) ? -length : length;
254  value -= additional_data.shift;
255 
256  // Update normalized eigenvector
257  x.equ(1 / length, y);
258 
259  // Compute residual
260  A.vmult(y, x);
261 
262  // Check the change of the eigenvalue
263  // Brrr, this is not really a good criterion
264  conv = this->iteration_status(iter,
265  std::fabs(1. / length - 1. / old_length),
266  x);
267  }
268 
269  // in case of failure: throw exception
272  iter, std::fabs(1. / length - 1. / old_length)));
273 
274  // otherwise exit as normal
275 }
276 
277 //---------------------------------------------------------------------------
278 
279 template <class VectorType>
282  const AdditionalData & data)
283  : SolverBase<VectorType>(cn, mem)
284  , additional_data(data)
285 {}
286 
287 
288 
289 template <class VectorType>
290 template <typename MatrixType>
291 void
293  const MatrixType &A,
294  VectorType & x)
295 {
296  LogStream::Prefix prefix("Wielandt");
297 
299 
300  // Prepare matrix for solver
301  auto A_op = linear_operator(A);
302  double current_shift = -value;
303  auto A_s = A_op + current_shift * identity_operator(A_op);
304 
305  // Define solver
306  ReductionControl inner_control(5000, 1.e-16, 1.e-5, false, false);
308  SolverGMRES<VectorType> solver(inner_control, this->memory);
309 
310  // Next step for recomputing the shift
311  unsigned int goal = additional_data.start_adaption;
312 
313  // Auxiliary vector
314  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
315  VectorType & y = *Vy;
316  y.reinit(x);
317  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
318  VectorType & r = *Vr;
319  r.reinit(x);
320 
321  double length = x.l2_norm();
322  double old_value = value;
323 
324  x *= 1. / length;
325 
326  // Main loop
327  double res = std::numeric_limits<double>::lowest();
328  size_type iter = 0;
329  for (; conv == SolverControl::iterate; ++iter)
330  {
331  solver.solve(A_s, y, x, prec);
332 
333  // Compute absolute value of eigenvalue
334  length = y.l2_norm();
335 
336  // do a little trick to compute the sign
337  // with not too much effect of round-off errors.
338  double entry = 0.;
339  size_type i = 0;
340  double thresh = length / x.size();
341  do
342  {
343  Assert(i < x.size(), ExcInternalError());
344  entry = y(i++);
345  }
346  while (std::fabs(entry) < thresh);
347 
348  --i;
349 
350  // Compute unshifted eigenvalue
351  value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;
352 
353  if (iter == goal)
354  {
355  const auto & relaxation = additional_data.relaxation;
356  const double new_shift =
357  relaxation * (-value) + (1. - relaxation) * current_shift;
358 
359  A_s = A_op + new_shift * identity_operator(A_op);
360  current_shift = new_shift;
361 
362  ++goal;
363  }
364 
365  // Update normalized eigenvector
366  x.equ(1. / length, y);
367  // Compute residual
368  if (additional_data.use_residual)
369  {
370  y.equ(value, x);
371  A.vmult(r, x);
372  r.sadd(-1., value, x);
373  res = r.l2_norm();
374  // Check the residual
375  conv = this->iteration_status(iter, res, x);
376  }
377  else
378  {
379  res = std::fabs(1. / value - 1. / old_value);
380  conv = this->iteration_status(iter, res, x);
381  }
382  old_value = value;
383  }
384 
385  // in case of failure: throw
386  // exception
388  SolverControl::NoConvergence(iter, res));
389  // otherwise exit as normal
390 }
391 
393 
394 #endif
types::global_dof_index size_type
Definition: eigen.h:135
EigenInverse(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:280
AdditionalData additional_data
Definition: eigen.h:189
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:292
EigenPower(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition: eigen.h:197
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:209
AdditionalData additional_data
Definition: eigen.h:104
types::global_dof_index size_type
Definition: eigen.h:62
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
Expression fabs(const Expression &x)
static const char A
unsigned int global_dof_index
Definition: types.h:81
AdditionalData(double relaxation=1., unsigned int start_adaption=6, bool use_residual=true)
Definition: eigen.h:158
unsigned int start_adaption
Definition: eigen.h:150
AdditionalData(const double shift=0.)
Definition: eigen.h:77