Reference documentation for deal.II version 9.5.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// Copyright (C) 2000 - 2022 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>
30
31#include <cmath>
32#include <limits>
33
35
36
55template <typename VectorType = Vector<double>>
56class EigenPower : private SolverBase<VectorType>
57{
58public:
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
100protected:
105};
106
128template <typename VectorType = Vector<double>>
129class EigenInverse : private SolverBase<VectorType>
130{
131public:
136
141 {
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
185protected:
190};
191
193//---------------------------------------------------------------------------
194
195
196template <class VectorType>
199 const AdditionalData & data)
200 : SolverBase<VectorType>(cn, mem)
201 , additional_data(data)
202{}
203
204
205
206template <class VectorType>
207template <typename MatrixType>
208void
209EigenPower<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
279template <class VectorType>
282 const AdditionalData & data)
283 : SolverBase<VectorType>(cn, mem)
284 , additional_data(data)
285{}
286
287
288
289template <class VectorType>
290template <typename MatrixType>
291void
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
389 // otherwise exit as normal
390}
391
393
394#endif
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
@ 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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:82
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