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
utilities.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2017 - 2023 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_lac_utilities_h
17#define dealii_lac_utilities_h
18
19#include <deal.II/base/config.h>
20
23
26
27#include <array>
28#include <complex>
29#include <limits>
30
32
33namespace Utilities
34{
38 namespace LinearAlgebra
39 {
64 template <typename NumberType>
65 std::array<NumberType, 3>
66 givens_rotation(const NumberType &x, const NumberType &y);
67
94 template <typename NumberType>
95 std::array<NumberType, 3>
96 hyperbolic_rotation(const NumberType &x, const NumberType &y);
97
133 template <typename OperatorType, typename VectorType>
134 double
135 lanczos_largest_eigenvalue(const OperatorType & H,
136 const VectorType & v0,
137 const unsigned int k,
138 VectorMemory<VectorType> &vector_memory,
139 std::vector<double> * eigenvalues = nullptr);
140
188 template <typename OperatorType, typename VectorType>
189 void
190 chebyshev_filter(VectorType & x,
191 const OperatorType & H,
192 const unsigned int n,
193 const std::pair<double, double> unwanted_spectrum,
194 const double tau,
195 VectorMemory<VectorType> & vector_memory);
196
197 } // namespace LinearAlgebra
198
199} // namespace Utilities
200
201
202/*------------------------- Implementation ----------------------------*/
203
204#ifndef DOXYGEN
205
206namespace internal
207{
208 namespace UtilitiesImplementation
209 {
210 // We want to avoid including our own LAPACK wrapper header in any external
211 // headers to avoid possible conflicts with other packages that may define
212 // their own such header. At the same time we want to be able to call some
213 // LAPACK functions from the template functions below. To resolve both
214 // problems define some extra wrappers here that can be in the header:
215 template <typename Number>
216 void
217 call_stev(const char jobz,
218 const types::blas_int n,
219 Number * d,
220 Number * e,
221 Number * z,
222 const types::blas_int ldz,
223 Number * work,
224 types::blas_int * info);
225 } // namespace UtilitiesImplementation
226} // namespace internal
227
228namespace Utilities
229{
230 namespace LinearAlgebra
231 {
232 template <typename NumberType>
233 std::array<std::complex<NumberType>, 3>
234 hyperbolic_rotation(const std::complex<NumberType> & /*f*/,
235 const std::complex<NumberType> & /*g*/)
236 {
238 std::array<NumberType, 3> res;
239 return res;
240 }
241
242
243
244 template <typename NumberType>
245 std::array<NumberType, 3>
246 hyperbolic_rotation(const NumberType &f, const NumberType &g)
247 {
248 Assert(f != 0, ExcDivideByZero());
249 const NumberType tau = g / f;
250 AssertThrow(std::abs(tau) < 1.,
252 "real-valued Hyperbolic rotation does not exist for (" +
253 std::to_string(f) + "," + std::to_string(g) + ")"));
254 const NumberType u =
255 std::copysign(std::sqrt((1. - tau) * (1. + tau)),
256 f); // <-- more stable than std::sqrt(1.-tau*tau)
257 std::array<NumberType, 3> csr;
258 csr[0] = 1. / u; // c
259 csr[1] = csr[0] * tau; // s
260 csr[2] = f * u; // r
261 return csr;
262 }
263
264
265
266 template <typename NumberType>
267 std::array<std::complex<NumberType>, 3>
268 givens_rotation(const std::complex<NumberType> & /*f*/,
269 const std::complex<NumberType> & /*g*/)
270 {
272 std::array<NumberType, 3> res;
273 return res;
274 }
275
276
277
278 template <typename NumberType>
279 std::array<NumberType, 3>
280 givens_rotation(const NumberType &f, const NumberType &g)
281 {
282 std::array<NumberType, 3> res;
283 // naive calculation for "r" may overflow or underflow:
284 // c = x / \sqrt{x^2+y^2}
285 // s = -y / \sqrt{x^2+y^2}
286
287 // See Golub 2013, Matrix computations, Chapter 5.1.8
288 // Algorithm 5.1.3
289 // and
290 // Anderson (2000),
291 // Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem.
292 // LAPACK Working Note 150, University of Tennessee, UT-CS-00-454,
293 // December 4, 2000.
294 // Algorithm 4
295 // We implement the latter below:
296 if (g == NumberType())
297 {
298 res[0] = std::copysign(1., f);
299 res[1] = NumberType();
300 res[2] = std::abs(f);
301 }
302 else if (f == NumberType())
303 {
304 res[0] = NumberType();
305 res[1] = std::copysign(1., g);
306 res[2] = std::abs(g);
307 }
308 else if (std::abs(f) > std::abs(g))
309 {
310 const NumberType tau = g / f;
311 const NumberType u = std::copysign(std::sqrt(1. + tau * tau), f);
312 res[0] = 1. / u; // c
313 res[1] = res[0] * tau; // s
314 res[2] = f * u; // r
315 }
316 else
317 {
318 const NumberType tau = f / g;
319 const NumberType u = std::copysign(std::sqrt(1. + tau * tau), g);
320 res[1] = 1. / u; // s
321 res[0] = res[1] * tau; // c
322 res[2] = g * u; // r
323 }
324
325 return res;
326 }
327
328
329
330 template <typename OperatorType, typename VectorType>
331 double
332 lanczos_largest_eigenvalue(const OperatorType & H,
333 const VectorType & v0_,
334 const unsigned int k,
335 VectorMemory<VectorType> &vector_memory,
336 std::vector<double> * eigenvalues)
337 {
338 // Do k-step Lanczos:
339
340 typename VectorMemory<VectorType>::Pointer v(vector_memory);
341 typename VectorMemory<VectorType>::Pointer v0(vector_memory);
342 typename VectorMemory<VectorType>::Pointer f(vector_memory);
343
344 v->reinit(v0_);
345 v0->reinit(v0_);
346 f->reinit(v0_);
347
348 // two vectors to store diagonal and subdiagonal of the Lanczos
349 // matrix
350 std::vector<double> diagonal;
351 std::vector<double> subdiagonal;
352
353 // 1. Normalize input vector
354 (*v) = v0_;
355 double a = v->l2_norm();
356 Assert(a != 0, ExcDivideByZero());
357 (*v) *= 1. / a;
358
359 // 2. Compute f = Hv; a = f*v; f <- f - av; T(0,0)=a;
360 H.vmult(*f, *v);
361 a = (*f) * (*v);
362 f->add(-a, *v);
363 diagonal.push_back(a);
364
365 // 3. Loop over steps
366 for (unsigned int i = 1; i < k; ++i)
367 {
368 // 4. L2 norm of f
369 const double b = f->l2_norm();
370 Assert(b != 0, ExcDivideByZero());
371 // 5. v0 <- v; v <- f/b
372 *v0 = *v;
373 *v = *f;
374 (*v) *= 1. / b;
375 // 6. f = Hv; f <- f - b v0;
376 H.vmult(*f, *v);
377 f->add(-b, *v0);
378 // 7. a = f*v; f <- f - a v;
379 a = (*f) * (*v);
380 f->add(-a, *v);
381 // 8. T(i,i-1) = T(i-1,i) = b; T(i,i) = a;
382 diagonal.push_back(a);
383 subdiagonal.push_back(b);
384 }
385
386 Assert(diagonal.size() == k, ExcInternalError());
387 Assert(subdiagonal.size() == k - 1, ExcInternalError());
388
389 // Use Lapack dstev to get ||T||_2 norm, i.e. the largest eigenvalue
390 // of T
391 const types::blas_int n = k;
392 std::vector<double> Z; // unused for eigenvalues-only ("N") job
393 const types::blas_int ldz = 1; // ^^ (>=1)
394 std::vector<double> work; // ^^
395 types::blas_int info;
396 // call lapack_templates.h wrapper:
398 n,
399 diagonal.data(),
400 subdiagonal.data(),
401 Z.data(),
402 ldz,
403 work.data(),
404 &info);
405
406 Assert(info == 0, LAPACKSupport::ExcErrorCode("dstev", info));
407
408 if (eigenvalues != nullptr)
409 {
410 eigenvalues->resize(diagonal.size());
411 std::copy(diagonal.begin(), diagonal.end(), eigenvalues->begin());
412 }
413
414 // note that the largest eigenvalue of T is below the largest
415 // eigenvalue of the operator.
416 // return ||T||_2 + ||f||_2, although it is not guaranteed to be an upper
417 // bound.
418 return diagonal[k - 1] + f->l2_norm();
419 }
420
421
422 template <typename OperatorType, typename VectorType>
423 void
424 chebyshev_filter(VectorType & x,
425 const OperatorType & op,
426 const unsigned int degree,
427 const std::pair<double, double> unwanted_spectrum,
428 const double a_L,
429 VectorMemory<VectorType> & vector_memory)
430 {
431 const double a = unwanted_spectrum.first;
432 const double b = unwanted_spectrum.second;
433 Assert(degree > 0, ExcMessage("Only positive degrees make sense."));
434
435 const bool scale = numbers::is_finite(a_L);
436 Assert(
437 a < b,
439 "Lower bound of the unwanted spectrum should be smaller than the upper bound."));
440
441 Assert(a_L <= a || a_L >= b || !scale,
443 "Scaling point should be outside of the unwanted spectrum."));
444
445 // Setup auxiliary vectors:
446 typename VectorMemory<VectorType>::Pointer p_y(vector_memory);
447 typename VectorMemory<VectorType>::Pointer p_yn(vector_memory);
448
449 p_y->reinit(x);
450 p_yn->reinit(x);
451
452 // convenience to avoid pointers
453 VectorType &y = *p_y;
454 VectorType &yn = *p_yn;
455
456 // Below is an implementation of
457 // Algorithm 3.2 in Zhou et al, Journal of Computational Physics 274
458 // (2014) 770-782 with **a bugfix for sigma1**. Here is the original
459 // algorithm verbatim:
460 //
461 // [Y]=chebyshev_filter_scaled(X, m, a, b, aL).
462 // e=(b-a)/2; c=(a+b)/2; σ=e/(c-aL); τ=2/σ;
463 // Y=(H∗X-c∗X)∗(σ/e);
464 // for i=2 to m do
465 // σnew =1/(τ - σ);
466 // Yt =(H∗Y - c∗Y)∗(2∗σnew/e)-(σ∗σnew)∗X;
467 // X =Y; Y =Yt; σ =σnew;
468
469 const double e = (b - a) / 2.;
470 const double c = (a + b) / 2.;
471 const double alpha = 1. / e;
472 const double beta = -c / e;
473
474 const double sigma1 =
475 e / (a_L - c); // BUGFIX which is relevant for odd degrees
476 double sigma = scale ? sigma1 : 1.;
477 const double tau = 2. / sigma;
478 op.vmult(y, x);
479 y.sadd(alpha * sigma, beta * sigma, x);
480
481 for (unsigned int i = 2; i <= degree; ++i)
482 {
483 const double sigma_new = scale ? 1. / (tau - sigma) : 1.;
484 op.vmult(yn, y);
485 yn.sadd(2. * alpha * sigma_new, 2. * beta * sigma_new, y);
486 yn.add(-sigma * sigma_new, x);
487 x.swap(y);
488 y.swap(yn);
489 sigma = sigma_new;
490 }
491
492 x.swap(y);
493 }
494
495 } // namespace LinearAlgebra
496} // namespace Utilities
497
498#endif
499
500
501
503
504
505#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
const unsigned int v0
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
@ diagonal
Matrix is diagonal.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void chebyshev_filter(VectorType &x, const OperatorType &H, const unsigned int n, const std::pair< double, double > unwanted_spectrum, const double tau, VectorMemory< VectorType > &vector_memory)
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
std::array< NumberType, 3 > hyperbolic_rotation(const NumberType &x, const NumberType &y)
double lanczos_largest_eigenvalue(const OperatorType &H, const VectorType &v0, const unsigned int k, VectorMemory< VectorType > &vector_memory, std::vector< double > *eigenvalues=nullptr)
void call_stev(const char jobz, const types::blas_int n, Number *d, Number *e, Number *z, const types::blas_int ldz, Number *work, types::blas_int *info)
Definition utilities.cc:32
bool is_finite(const double x)
Definition numbers.h:538
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)