Reference documentation for deal.II version 9.6.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\}}\)
Loading...
Searching...
No Matches
eigen.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_eigen_h
16#define dealii_eigen_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/lac/solver.h>
29
30#include <cmath>
31#include <limits>
32
34
35
54template <typename VectorType = Vector<double>>
55class EigenPower : private SolverBase<VectorType>
56{
57public:
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
99protected:
104};
105
127template <typename VectorType = Vector<double>>
128class EigenInverse : private SolverBase<VectorType>
129{
130public:
135
140 {
145
149 unsigned int start_adaption;
158 unsigned int start_adaption = 6,
159 bool use_residual = true)
163 {}
164 };
165
171 const AdditionalData &data = AdditionalData());
172
180 template <typename MatrixType>
181 void
182 solve(double &value, const MatrixType &A, VectorType &x);
183
184protected:
189};
190
192//---------------------------------------------------------------------------
193
194
195template <typename VectorType>
198 const AdditionalData &data)
199 : SolverBase<VectorType>(cn, mem)
200 , additional_data(data)
201{}
202
203
204
205template <typename VectorType>
206template <typename MatrixType>
207void
208EigenPower<VectorType>::solve(double &value, const MatrixType &A, VectorType &x)
209{
211
212 LogStream::Prefix prefix("Power method");
213
214 typename VectorMemory<VectorType>::Pointer Vy(this->memory);
215 VectorType &y = *Vy;
216 y.reinit(x);
217 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
218 VectorType &r = *Vr;
219 r.reinit(x);
220
221 double length = x.l2_norm();
222 double old_length = 0.;
223 x *= 1. / length;
224
225 A.vmult(y, x);
226
227 // Main loop
228 int iter = 0;
229 for (; conv == SolverControl::iterate; ++iter)
230 {
231 y.add(additional_data.shift, x);
232
233 // Compute absolute value of eigenvalue
234 old_length = length;
235 length = y.l2_norm();
236
237 // do a little trick to compute the sign
238 // with not too much effect of round-off errors.
239 double entry = 0.;
240 size_type i = 0;
241 double thresh = length / x.size();
242 do
243 {
244 Assert(i < x.size(), ExcInternalError());
245 entry = y(i++);
246 }
247 while (std::fabs(entry) < thresh);
248
249 --i;
250
251 // Compute unshifted eigenvalue
252 value = (entry * x(i) < 0.) ? -length : length;
253 value -= additional_data.shift;
254
255 // Update normalized eigenvector
256 x.equ(1 / length, y);
257
258 // Compute residual
259 A.vmult(y, x);
260
261 // Check the change of the eigenvalue
262 // Brrr, this is not really a good criterion
263 conv = this->iteration_status(iter,
264 std::fabs(1. / length - 1. / old_length),
265 x);
266 }
267
268 // in case of failure: throw exception
271 iter, std::fabs(1. / length - 1. / old_length)));
272
273 // otherwise exit as normal
274}
275
276//---------------------------------------------------------------------------
277
278template <typename VectorType>
281 const AdditionalData &data)
282 : SolverBase<VectorType>(cn, mem)
283 , additional_data(data)
284{}
285
286
287
288template <typename VectorType>
289template <typename MatrixType>
290void
292 const MatrixType &A,
293 VectorType &x)
294{
295 LogStream::Prefix prefix("Wielandt");
296
298
299 // Prepare matrix for solver
300 auto A_op = linear_operator(A);
301 double current_shift = -value;
302 auto A_s = A_op + current_shift * identity_operator(A_op);
303
304 // Define solver
305 ReductionControl inner_control(5000, 1.e-16, 1.e-5, false, false);
307 SolverGMRES<VectorType> solver(inner_control, this->memory);
308
309 // Next step for recomputing the shift
310 unsigned int goal = additional_data.start_adaption;
311
312 // Auxiliary vector
313 typename VectorMemory<VectorType>::Pointer Vy(this->memory);
314 VectorType &y = *Vy;
315 y.reinit(x);
316 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
317 VectorType &r = *Vr;
318 r.reinit(x);
319
320 double length = x.l2_norm();
321 double old_value = value;
322
323 x *= 1. / length;
324
325 // Main loop
326 double res = std::numeric_limits<double>::lowest();
327 size_type iter = 0;
328 for (; conv == SolverControl::iterate; ++iter)
329 {
330 solver.solve(A_s, y, x, prec);
331
332 // Compute absolute value of eigenvalue
333 length = y.l2_norm();
334
335 // do a little trick to compute the sign
336 // with not too much effect of round-off errors.
337 double entry = 0.;
338 size_type i = 0;
339 double thresh = length / x.size();
340 do
341 {
342 Assert(i < x.size(), ExcInternalError());
343 entry = y(i++);
344 }
345 while (std::fabs(entry) < thresh);
346
347 --i;
348
349 // Compute unshifted eigenvalue
350 value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;
351
352 if (iter == goal)
353 {
354 const auto &relaxation = additional_data.relaxation;
355 const double new_shift =
356 relaxation * (-value) + (1. - relaxation) * current_shift;
357
358 A_s = A_op + new_shift * identity_operator(A_op);
359 current_shift = new_shift;
360
361 ++goal;
362 }
363
364 // Update normalized eigenvector
365 x.equ(1. / length, y);
366 // Compute residual
367 if (additional_data.use_residual)
368 {
369 y.equ(value, x);
370 A.vmult(r, x);
371 r.sadd(-1., value, x);
372 res = r.l2_norm();
373 // Check the residual
374 conv = this->iteration_status(iter, res, x);
375 }
376 else
377 {
378 res = std::fabs(1. / value - 1. / old_value);
379 conv = this->iteration_status(iter, res, x);
380 }
381 old_value = value;
382 }
383
384 // in case of failure: throw
385 // exception
388 // otherwise exit as normal
389}
390
392
393#endif
EigenInverse(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition eigen.h:279
AdditionalData additional_data
Definition eigen.h:188
void solve(double &value, const MatrixType &A, VectorType &x)
Definition eigen.h:291
EigenPower(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Definition eigen.h:196
void solve(double &value, const MatrixType &A, VectorType &x)
Definition eigen.h:208
AdditionalData additional_data
Definition eigen.h:103
@ 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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &)
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:157
unsigned int start_adaption
Definition eigen.h:149
AdditionalData(const double shift=0.)
Definition eigen.h:76