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