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