QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
tqreigendecomposition.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2005 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file tqreigendecomposition.hpp
21 \brief tridiag. QR eigen decompositions with implicit shift
22*/
23
25#include <vector>
26
27namespace QuantLib {
28
30 const Array& sub,
32 ShiftStrategy strategy)
33 : d_(diag), ev_((calc == WithEigenVector) ? d_.size() :
34 (calc == WithoutEigenVector) ? 0 :
35 1,
36 d_.size(),
37 0) {
38 Size n = diag.size();
39
40 QL_REQUIRE(n == sub.size()+1, "Wrong dimensions");
41
42 Array e(n, 0.0);
43 std::copy(sub.begin(),sub.end(),e.begin()+1);
44 for (Size i=0; i < ev_.rows(); ++i) {
45 ev_[i][i] = 1.0;
46 }
47
48 for (Size k=n-1; k >=1; --k) {
49 while (!offDiagIsZero(k, e)) {
50 Size l = k;
51 while (--l > 0 && !offDiagIsZero(l,e));
52 iter_++;
53
54 Real q = d_[l];
55 if (strategy != NoShift) {
56 // calculated eigenvalue of 2x2 sub matrix of
57 // [ d_[k-1] e_[k] ]
58 // [ e_[k] d_[k] ]
59 // which is closer to d_[k+1].
60 const Real t1 = std::sqrt(
61 0.25*(d_[k]*d_[k] + d_[k-1]*d_[k-1])
62 - 0.5*d_[k-1]*d_[k] + e[k]*e[k]);
63 const Real t2 = 0.5*(d_[k]+d_[k-1]);
64
65 const Real lambda =
66 (std::fabs(t2+t1 - d_[k]) < std::fabs(t2-t1 - d_[k]))?
67 Real(t2+t1) : Real(t2-t1);
68
69 if (strategy == CloseEigenValue) {
70 q-=lambda;
71 } else {
72 q-=((k==n-1)? 1.25 : 1.0)*lambda;
73 }
74 }
75
76 // the QR transformation
77 Real sine = 1.0;
78 Real cosine = 1.0;
79 Real u = 0.0;
80
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];
85
86 e[i-1] = std::sqrt(p*p+q*q);
87 if (e[i-1] != 0.0) {
88 sine = p/e[i-1];
89 cosine = q/e[i-1];
90
91 const Real g = d_[i-1]-u;
92 const Real t = (d_[i]-g)*sine+2*cosine*h;
93
94 u = sine*t;
95 d_[i-1] = g + u;
96 q = cosine*t - h;
97
98 for (Size j=0; j < ev_.rows(); ++j) {
99 const Real tmp = ev_[j][i-1];
100 ev_[j][i-1] = sine*ev_[j][i] + cosine*tmp;
101 ev_[j][i] = cosine*ev_[j][i] - sine*tmp;
102 }
103 } else {
104 // recover from underflow
105 d_[i-1] -= u;
106 e[l] = 0.0;
107 recoverUnderflow = true;
108 }
109 }
110
111 if (!recoverUnderflow) {
112 d_[k] -= u;
113 e[k] = q;
114 e[l] = 0.0;
115 }
116 }
117 }
118
119 // sort (eigenvalues, eigenvectors),
120 // code taken from symmetricSchureDecomposition.cpp
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++) {
124 if (ev_.rows() > 0)
125 std::copy(ev_.column_begin(i),
126 ev_.column_end(i), eigenVector.begin());
127 temp[i] = std::make_pair(d_[i], eigenVector);
128 }
129 std::sort(temp.begin(), temp.end(), std::greater<>());
130 // first element is positive
131 for (Size i=0; i<n; i++) {
132 d_[i] = temp[i].first;
133 Real sign = 1.0;
134 if (ev_.rows() > 0 && temp[i].second[0]<0.0)
135 sign = -1.0;
136 for (Size j=0; j<ev_.rows(); ++j) {
137 ev_[j][i] = sign * temp[i].second[j];
138 }
139 }
140 }
141
142 // see NR for abort assumption as it is
143 // not part of the original Wilkinson algorithm
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]);
147 }
148
149}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
Size rows() const
Definition: matrix.hpp:504
const_column_iterator column_begin(Size i) const
Definition: matrix.hpp:415
const_column_iterator column_end(Size i) const
Definition: matrix.hpp:434
TqrEigenDecomposition(const Array &diag, const Array &sub, EigenVectorCalculation calc=WithEigenVector, ShiftStrategy strategy=CloseEigenValue)
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
ext::shared_ptr< YieldTermStructure > q
tridiag. QR eigen decompositions with implicit shift