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