27 : diagonal_(
s.rows()), eigenVectors_(
s.rows(),
s.columns(), 0.0) {
29 QL_REQUIRE(
s.rows() > 0 &&
s.columns() > 0,
"null matrix given");
30 QL_REQUIRE(
s.rows()==
s.columns(),
"input matrix must be square");
40 std::vector<Real> tmpAccumulate(size, 0.0);
41 Real threshold, epsPrec = 1e-15;
42 bool keeplooping =
true;
43 Size maxIterations = 100, ite = 1;
47 for (
Size a=0; a<size-1; a++) {
49 sum += std::fabs(ss[a][
b]);
59 if (ite<5) threshold = 0.2*sum/(size*size);
63 for (j=0; j<size-1; j++) {
64 for (k=j+1; k<size; k++) {
66 Real smll = std::fabs(ss[j][k]);
71 }
else if (std::fabs(ss[j][k])>threshold) {
73 if (smll<epsPrec*std::fabs(heig)) {
76 beta = 0.5*heig/ss[j][k];
77 tang = 1.0/(std::fabs(
beta)+
82 cosin = 1/std::sqrt(1+tang*tang);
86 tmpAccumulate[j] -= heig;
87 tmpAccumulate[k] += heig;
91 for (l=0; l+1<=j; l++)
93 for (l=j+1; l<=k-1; l++)
95 for (l=k+1; l<size; l++)
97 for (l=0; l<size; l++)
99 rho, sine, l, j, l, k);
103 for (k=0; k<size; k++) {
104 tmpDiag[k] += tmpAccumulate[k];
106 tmpAccumulate[k] = 0.0;
109 }
while (++ite<=maxIterations && keeplooping);
112 "Too many iterations (" << maxIterations <<
") reached");
116 std::vector<std::pair<Real, std::vector<Real> > > temp(size);
117 std::vector<Real> eigenVector(size);
119 for (col=0; col<size; col++) {
122 temp[col] = std::make_pair(
diagonal_[col], eigenVector);
124 std::sort(temp.begin(), temp.end(), std::greater<>());
125 Real maxEv = temp[0].first;
126 for (col=0; col<size; col++) {
129 (std::fabs(temp[col].first/maxEv)<1e-16 ? 0.0 :
132 if (temp[col].second[0]<0.0)
134 for (row=0; row<size; row++) {
const_iterator end() const
const_iterator begin() const
Matrix used in linear algebra.
const_column_iterator column_begin(Size i) const
const_column_iterator column_end(Size i) const
SymmetricSchurDecomposition(const Matrix &s)
void jacobiRotate_(Matrix &m, Real rot, Real dil, Size j1, Size k1, Size j2, Size k2) const
This routines implements the Jacobi, a.k.a. Givens, rotation.
#define QL_ENSURE(condition, message)
throw an error if the given post-condition is not verified
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
ext::function< Real(Real)> b
std::size_t Size
size of a container
ext::shared_ptr< YieldTermStructure > q
Eigenvalues/eigenvectors of a real symmetric matrix.