16 #ifndef dealii_lac_utilities_h 17 #define dealii_lac_utilities_h 19 #include <deal.II/base/config.h> 21 #include <deal.II/base/exceptions.h> 22 #include <deal.II/base/signaling_nan.h> 24 #include <deal.II/lac/lapack_support.h> 25 #include <deal.II/lac/lapack_templates.h> 26 #include <deal.II/lac/vector_memory.h> 31 DEAL_II_NAMESPACE_OPEN
66 template <
typename NumberType>
67 std::array<NumberType, 3>
98 template <
typename NumberType>
99 std::array<NumberType, 3>
139 template <
typename OperatorType,
typename VectorType>
142 const VectorType & v0,
143 const unsigned int k,
196 template <
typename OperatorType,
typename VectorType>
199 const OperatorType & H,
200 const unsigned int n,
201 const std::pair<double, double> unwanted_spectrum,
218 template <
typename NumberType>
219 std::array<std::complex<NumberType>, 3>
221 const std::complex<NumberType> & )
224 std::array<NumberType, 3> res;
230 template <
typename NumberType>
231 std::array<NumberType, 3>
235 const NumberType tau = g / f;
238 "real-valued Hyperbolic rotation does not exist for (" +
239 std::to_string(f) +
"," + std::to_string(g) +
")"));
241 std::copysign(std::sqrt((1. - tau) * (1. + tau)),
243 std::array<NumberType, 3> csr;
245 csr[1] = csr[0] * tau;
252 template <
typename NumberType>
253 std::array<std::complex<NumberType>, 3>
255 const std::complex<NumberType> & )
258 std::array<NumberType, 3> res;
264 template <
typename NumberType>
265 std::array<NumberType, 3>
268 std::array<NumberType, 3> res;
282 if (g == NumberType())
284 res[0] = std::copysign(1., f);
285 res[1] = NumberType();
286 res[2] = std::abs(f);
288 else if (f == NumberType())
290 res[0] = NumberType();
291 res[1] = std::copysign(1., g);
292 res[2] = std::abs(g);
294 else if (std::abs(f) > std::abs(g))
296 const NumberType tau = g / f;
297 const NumberType u = std::copysign(std::sqrt(1. + tau * tau), f);
299 res[1] = res[0] * tau;
304 const NumberType tau = f / g;
305 const NumberType u = std::copysign(std::sqrt(1. + tau * tau), g);
307 res[0] = res[1] * tau;
316 template <
typename OperatorType,
typename VectorType>
319 const VectorType & v0_,
320 const unsigned int k,
337 std::vector<double> subdiagonal;
341 double a = v->l2_norm();
352 for (
unsigned int i = 1; i < k; ++i)
355 const double b = f->l2_norm();
369 subdiagonal.push_back(b);
378 std::vector<double> Z;
380 std::vector<double> work;
404 return diagonal[k - 1] + f->l2_norm();
408 template <
typename OperatorType,
typename VectorType>
411 const OperatorType & op,
412 const unsigned int degree,
413 const std::pair<double, double> unwanted_spectrum,
417 const double a = unwanted_spectrum.first;
418 const double b = unwanted_spectrum.second;
421 const bool scale = (a_L < std::numeric_limits<double>::infinity());
425 "Lower bound of the unwanted spectrum should be smaller than the upper bound."));
429 "Scaling point should be outside of the unwanted spectrum."));
439 VectorType &y = *p_y;
440 VectorType &yn = *p_yn;
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;
460 const double sigma1 =
462 double sigma =
scale ? sigma1 : 1.;
463 const double tau = 2. / sigma;
465 y.sadd(alpha * sigma, beta * sigma, x);
467 for (
unsigned int i = 2; i <= degree; ++i)
469 const double sigma_new =
scale ? 1. / (tau - sigma) : 1.;
471 yn.sadd(2. * alpha * sigma_new, 2. * beta * sigma_new, y);
472 yn.add(-sigma * sigma_new, x);
488 DEAL_II_NAMESPACE_CLOSE
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)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
static ::ExceptionBase & ExcDivideByZero()
__global__ void scale(Number *val, const Number *V_val, const size_type N)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
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)
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInternalError()