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