16 #ifndef dealii_lac_utilities_h 17 #define dealii_lac_utilities_h 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> 29 DEAL_II_NAMESPACE_OPEN
65 template<
typename NumberType>
97 template<
typename NumberType>
135 template <
typename OperatorType,
typename VectorType>
137 const VectorType &v0,
138 const unsigned int k,
191 template <
typename OperatorType,
typename VectorType>
193 const OperatorType &H,
194 const unsigned int n,
195 const std::pair<double,double> unwanted_spectrum,
213 template<
typename NumberType>
214 std::array<std::complex<NumberType>,3>
hyperbolic_rotation(
const std::complex<NumberType> &,
215 const std::complex<NumberType> &)
218 std::array<NumberType,3> res;
224 template<
typename NumberType>
229 const NumberType tau = g/f;
231 ExcMessage(
"real-valued Hyperbolic rotation does not exist for ("+
236 const NumberType u = std::copysign(sqrt((1.-tau)*(1.+tau)), f);
237 std::array<NumberType,3> csr;
239 csr[1] = csr[0] * tau;
246 template<
typename NumberType>
247 std::array<std::complex<NumberType>,3>
givens_rotation(
const std::complex<NumberType> &,
248 const std::complex<NumberType> &)
251 std::array<NumberType,3> res;
257 template<
typename NumberType>
261 std::array<NumberType,3> res;
275 if (g == NumberType())
277 res[0] = std::copysign(1., f);
278 res[1] = NumberType();
279 res[2] = std::abs(f);
281 else if (f == NumberType())
283 res[0] = NumberType();
284 res[1] = std::copysign(1., g);
285 res[2] = std::abs(g);
287 else if (std::abs(f) > std::abs(g))
289 const NumberType tau = g/f;
290 const NumberType u = std::copysign(std::sqrt(1.+tau*tau), f);
292 res[1] = res[0] * tau;
297 const NumberType tau = f/g;
298 const NumberType u = std::copysign(std::sqrt(1.+tau*tau), g);
300 res[0] = res[1] * tau;
309 template <
typename OperatorType,
typename VectorType>
311 const VectorType &v0_,
312 const unsigned int k,
329 std::vector<double> subdiagonal;
333 double a = v->l2_norm();
344 for (
unsigned int i = 1; i < k; ++i)
347 const double b = f->l2_norm();
361 subdiagonal.push_back(b);
366 Assert (subdiagonal.size() == k-1,
372 std::vector<double> Z;
374 std::vector<double> work;
378 diagonal.data(), subdiagonal.data(),
379 Z.data(), &ldz, work.data(),
395 return diagonal[k-1] + f->l2_norm();
399 template <
typename OperatorType,
typename VectorType>
401 const OperatorType &op,
402 const unsigned int degree,
403 const std::pair<double,double> unwanted_spectrum,
407 const double a = unwanted_spectrum.first;
408 const double b = unwanted_spectrum.second;
410 ExcMessage (
"Only positive degrees make sense."));
412 const bool scale = (a_L < std::numeric_limits<double>::infinity());
414 ExcMessage(
"Lower bound of the unwanted spectrum should be smaller than the upper bound."));
416 Assert (a_L <= a || a_L >= b || !scale,
417 ExcMessage(
"Scaling point should be outside of the unwanted spectrum."));
427 VectorType &y = *p_y;
428 VectorType &yn = *p_yn;
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;
447 const double sigma1 =
e/(a_L - c);
448 double sigma =
scale ? sigma1 : 1.;
449 const double tau = 2./sigma;
451 y.sadd(alpha*sigma, beta*sigma, x);
453 for (
unsigned int i = 2; i <= degree; ++i)
455 const double sigma_new =
scale ? 1./(tau-sigma) : 1.;
457 yn.sadd(2.*alpha*sigma_new, 2.*beta*sigma_new, y);
458 yn.add(-sigma*sigma_new, x);
474 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)
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)
#define AssertThrow(cond, exc)
std::array< NumberType, 3 > hyperbolic_rotation(const NumberType &x, const NumberType &y)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
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)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInternalError()