33 :
d_(diag), ev_((calc == WithEigenVector) ?
d_.size() :
34 (calc == WithoutEigenVector) ? 0 :
48 for (
Size k=
n-1; k >=1; --k) {
60 const Real t1 = std::sqrt(
62 - 0.5*
d_[k-1]*
d_[k] + e[k]*e[k]);
66 (std::fabs(t2+t1 -
d_[k]) < std::fabs(t2-t1 -
d_[k]))?
72 q-=((k==
n-1)? 1.25 : 1.0)*lambda;
81 bool recoverUnderflow =
false;
82 for (
Size i=l+1; i <= k && !recoverUnderflow; ++i) {
83 const Real h = cosine*e[i];
84 const Real p = sine*e[i];
86 e[i-1] = std::sqrt(p*p+
q*
q);
92 const Real t = (
d_[i]-g)*sine+2*cosine*h;
100 ev_[j][i-1] = sine*
ev_[j][i] + cosine*tmp;
101 ev_[j][i] = cosine*
ev_[j][i] - sine*tmp;
107 recoverUnderflow =
true;
111 if (!recoverUnderflow) {
121 std::vector<std::pair<Real, std::vector<Real> > > temp(
n);
122 std::vector<Real> eigenVector(
ev_.
rows());
123 for (
Size i=0; i<
n; i++) {
127 temp[i] = std::make_pair(
d_[i], eigenVector);
129 std::sort(temp.begin(), temp.end(), std::greater<>());
131 for (
Size i=0; i<
n; i++) {
132 d_[i] = temp[i].first;
134 if (
ev_.
rows() > 0 && temp[i].second[0]<0.0)
137 ev_[j][i] = sign * temp[i].second[j];
145 return std::fabs(
d_[k-1])+std::fabs(
d_[k])
146 == std::fabs(
d_[k-1])+std::fabs(
d_[k])+std::fabs(e[k]);
1-D array used in linear algebra.
const_iterator end() const
Size size() const
dimension of the array
const_iterator begin() const
const_column_iterator column_begin(Size i) const
const_column_iterator column_end(Size i) const
TqrEigenDecomposition(const Array &diag, const Array &sub, EigenVectorCalculation calc=WithEigenVector, ShiftStrategy strategy=CloseEigenValue)
bool offDiagIsZero(Size k, Array &e)
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
std::size_t Size
size of a container
ext::shared_ptr< YieldTermStructure > q
tridiag. QR eigen decompositions with implicit shift