12#ifndef TLAPACK_LAHQR_HH
13#define TLAPACK_LAHQR_HH
98 const idx_t n = ncols(
A);
109 if (
nh <= 0)
return 0;
115 const idx_t
itmax = 30 * std::max<idx_t>(10,
nh);
134 "The QR algorithm failed to compute all the eigenvalues"
135 " in a total of 30 iterations per eigenvalue. Elements"
136 " i:ihi of w contain those eigenvalues which have been"
137 " successfully computed.");
171 tst =
tst + abs(
A(i - 1, i - 2));
306 conj(
v[0]) *
A(i, i - 1) +
conj(
v[1]) *
A(i + 1, i - 1);
310 abs1(
A(i + 1, i + 2)))) {
318 const idx_t
nr = std::min<idx_t>(3,
istop - i);
333 v[1] =
A(i + 1, i - 1);
334 if (
nr == 3)
v[2] =
A(i + 2, i - 1);
340 A(i + 1, i - 1) =
zero;
341 if (
nr == 3)
A(i + 2, i - 1) =
zero;
357 for (idx_t j = i; j <
istop_m; ++j) {
361 A(i + 1, j) =
A(i + 1, j) -
sum *
t2;
362 A(i + 2, j) =
A(i + 2, j) -
sum *
t3;
366 sum =
A(j, i) +
v2 *
A(j, i + 1) +
v3 *
A(j, i + 2);
373 for (idx_t j = 0; j < n; ++j) {
374 sum =
Z(j, i) +
v2 *
Z(j, i + 1) +
v3 *
Z(j, i + 2);
383 for (idx_t j = i; j <
istop_m; ++j) {
386 A(i + 1, j) =
A(i + 1, j) -
sum *
t2;
390 sum =
A(j, i) +
v2 *
A(j, i + 1);
396 for (idx_t j = 0; j < n; ++j) {
397 sum =
Z(j, i) +
v2 *
Z(j, i + 1);
439 const idx_t n = ncols(
A);
450 if (
nh <= 0)
return 0;
456 const idx_t
itmax = 30 * std::max<idx_t>(10,
nh);
506 tst =
tst + abs(
A(i - 1, i - 2));
600 if (i >
istart)
A(i, i - 1) =
A(i, i - 1) *
cs;
607 A(i + 1, i - 1) =
zero;
611 for (idx_t j = i; j <
istop_m; ++j) {
613 A(i + 1, j) = -
conj(
sn) *
A(i, j) +
cs *
A(i + 1, j);
619 A(j, i + 1) = -
sn *
A(j, i) +
cs *
A(j, i + 1);
624 for (idx_t j = 0; j < n; ++j) {
626 Z(j, i + 1) = -
sn *
Z(j, i) +
cs *
Z(j, i + 1);
#define TLAPACK_CSMATRIX
Macro for tlapack::concepts::ConstructableAndSliceableMatrix compatible with C++17.
Definition concepts.hpp:961
#define TLAPACK_VECTOR
Macro for tlapack::concepts::Vector compatible with C++17.
Definition concepts.hpp:906
void lahqr_eig22(T a00, T a01, T a10, T a11, complex_type< T > &s1, complex_type< T > &s2)
Computes the eigenvalues of a 2x2 matrix A.
Definition lahqr_eig22.hpp:34
int lahqr(bool want_t, bool want_z, size_type< matrix_t > ilo, size_type< matrix_t > ihi, matrix_t &A, vector_t &w, matrix_t &Z)
lahqr computes the eigenvalues and optionally the Schur factorization of an upper Hessenberg matrix,...
Definition lahqr.hpp:73
void larfg(storage_t storeMode, type_t< vector_t > &alpha, vector_t &x, type_t< vector_t > &tau)
Generates a elementary Householder reflection.
Definition larfg.hpp:73
void lahqr_schur22(T &a, T &b, T &c, T &d, complex_type< T > &s1, complex_type< T > &s2, T &cs, T &sn)
Computes the Schur factorization of a 2x2 matrix A.
Definition lahqr_schur22.hpp:44
int lahqr_shiftcolumn(const matrix_t &H, vector_t &v, complex_type< type_t< matrix_t > > s1, complex_type< type_t< matrix_t > > s2)
Given a 2-by-2 or 3-by-3 matrix H, lahqr_shiftcolumn calculates a multiple of the product: (H - s1*I)...
Definition lahqr_shiftcolumn.hpp:41
void rotg(T &a, T &b, T &c, T &s)
Construct plane rotation that eliminates b, such that:
Definition rotg.hpp:39
void rot(vectorX_t &x, vectorY_t &y, const c_type &c, const s_type &s)
Apply plane rotation:
Definition rot.hpp:44
#define tlapack_error(info, detailedInfo)
Error handler.
Definition exceptionHandling.hpp:142
#define tlapack_check_false(cond)
Throw an error if cond is true.
Definition exceptionHandling.hpp:113
Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ).
Definition arrayTraits.hpp:15
typename traits::real_type_traits< Types..., int >::type real_type
The common real type of the list of types.
Definition scalar_type_traits.hpp:113
constexpr T conj(const T &x) noexcept
Extends std::conj() to real datatypes.
Definition utils.hpp:100
constexpr real_type< T > abs1(const T &x)
1-norm absolute value, |Re(x)| + |Im(x)|
Definition utils.hpp:133
constexpr internal::Forward FORWARD
Forward direction.
Definition types.hpp:381
constexpr real_type< T > imag(const T &x) noexcept
Extends std::imag() to real datatypes.
Definition utils.hpp:86
constexpr internal::ColumnwiseStorage COLUMNWISE_STORAGE
Columnwise storage.
Definition types.hpp:414