15#ifndef dealii_lac_utilities_h
16#define dealii_lac_utilities_h
63 template <
typename NumberType>
64 std::array<NumberType, 3>
93 template <
typename NumberType>
94 std::array<NumberType, 3>
120 template <
typename OperatorType,
typename VectorType>
123 const VectorType &
v0,
124 const unsigned int k,
164 template <
typename OperatorType,
typename VectorType>
168 const unsigned int n,
184 namespace UtilitiesImplementation
191 template <
typename Number>
208 template <
typename NumberType>
209 std::array<std::complex<NumberType>, 3>
211 const std::complex<NumberType> & )
214 std::array<NumberType, 3>
res;
220 template <
typename NumberType>
221 std::array<NumberType, 3>
225 const NumberType
tau =
g / f;
228 "real-valued Hyperbolic rotation does not exist for (" +
229 std::to_string(f) +
"," + std::to_string(
g) +
")"));
233 std::array<NumberType, 3>
csr;
242 template <
typename NumberType>
243 std::array<std::complex<NumberType>, 3>
245 const std::complex<NumberType> & )
248 std::array<NumberType, 3>
res;
254 template <
typename NumberType>
255 std::array<NumberType, 3>
258 std::array<NumberType, 3>
res;
272 if (
g == NumberType())
274 res[0] = std::copysign(1., f);
275 res[1] = NumberType();
278 else if (f == NumberType())
280 res[0] = NumberType();
281 res[1] = std::copysign(1.,
g);
286 const NumberType
tau =
g / f;
294 const NumberType
tau = f /
g;
306 template <
typename OperatorType,
typename VectorType>
309 const VectorType &
v0_,
310 const unsigned int k,
331 double a = v->l2_norm();
342 for (
unsigned int i = 1; i <
k; ++i)
345 const double b = f->l2_norm();
368 std::vector<double>
Z;
370 std::vector<double> work;
398 template <
typename OperatorType,
typename VectorType>
402 const unsigned int degree,
415 "Lower bound of the unwanted spectrum should be smaller than the upper bound."));
419 "Scaling point should be outside of the unwanted spectrum."));
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;
457 for (
unsigned int i = 2; i <= degree; ++i)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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)
@ diagonal
Matrix is diagonal.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
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)
bool is_finite(const double x)
::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)