Reference documentation for deal.II version 9.4.1
Searching...
No Matches
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
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
34
35
38
52template <typename VectorType = Vector<double>>
53class EigenPower : private SolverBase<VectorType>
54{
55public:
60
65 {
70 double shift;
74 AdditionalData(const double shift = 0.)
75 : shift(shift)
76 {}
77 };
78
85
86
93 template <typename MatrixType>
94 void
95 solve(double &value, const MatrixType &A, VectorType &x);
96
97protected:
102};
103
125template <typename VectorType = Vector<double>>
126class EigenInverse : private SolverBase<VectorType>
127{
128public:
133
138 {
143
156 unsigned int start_adaption = 6,
157 bool use_residual = true)
161 {}
162 };
163
170
178 template <typename MatrixType>
179 void
180 solve(double &value, const MatrixType &A, VectorType &x);
181
182protected:
187};
188
190//---------------------------------------------------------------------------
191
192
193template <class VectorType>
197 : SolverBase<VectorType>(cn, mem)
199{}
200
201
202
203template <class VectorType>
204template <typename MatrixType>
205void
206EigenPower<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 {
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;
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
276template <class VectorType>
280 : SolverBase<VectorType>(cn, mem)
282{}
283
284
285
286template <class VectorType>
287template <typename MatrixType>
288void
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
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>::lowest();
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);
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
386 // otherwise exit as normal
387}
388
390
391#endif
Definition: eigen.h:277
Definition: eigen.h:186
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:289
Definition: eigen.h:194
void solve(double &value, const MatrixType &A, VectorType &x)
Definition: eigen.h:206
Definition: eigen.h:101
@ 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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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:76