Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20:02+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\}}\)
Loading...
Searching...
No Matches
utilities.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_lac_utilities_h
16#define dealii_lac_utilities_h
17
18#include <deal.II/base/config.h>
19
22
25
26#include <array>
27#include <complex>
28#include <limits>
29
31
32namespace 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
205namespace 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
227namespace 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 {
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.,
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 {
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 = numbers::is_finite(a_L);
435 Assert(
436 a < b,
438 "Lower bound of the unwanted spectrum should be smaller than the upper bound."));
439
440 Assert(a_L <= a || a_L >= b || !scale,
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:31
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)