QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
symmetricschurdecomposition.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) 2003 Ferdinando Ametrano
5 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
22#include <vector>
23
24namespace QuantLib {
25
27 : diagonal_(s.rows()), eigenVectors_(s.rows(), s.columns(), 0.0) {
28
29 QL_REQUIRE(s.rows() > 0 && s.columns() > 0, "null matrix given");
30 QL_REQUIRE(s.rows()==s.columns(), "input matrix must be square");
31
32 Size size = s.rows();
33 for (Size q=0; q<size; q++) {
34 diagonal_[q] = s[q][q];
35 eigenVectors_[q][q] = 1.0;
36 }
37 Matrix ss = s;
38
39 std::vector<Real> tmpDiag(diagonal_.begin(), diagonal_.end());
40 std::vector<Real> tmpAccumulate(size, 0.0);
41 Real threshold, epsPrec = 1e-15;
42 bool keeplooping = true;
43 Size maxIterations = 100, ite = 1;
44 do {
45 //main loop
46 Real sum = 0;
47 for (Size a=0; a<size-1; a++) {
48 for (Size b=a+1; b<size; b++) {
49 sum += std::fabs(ss[a][b]);
50 }
51 }
52
53 if (sum==0) {
54 keeplooping = false;
55 } else {
56 /* To speed up computation a threshold is introduced to
57 make sure it is worthy to perform the Jacobi rotation
58 */
59 if (ite<5) threshold = 0.2*sum/(size*size);
60 else threshold = 0.0;
61
62 Size j, k, l;
63 for (j=0; j<size-1; j++) {
64 for (k=j+1; k<size; k++) {
65 Real sine, rho, cosin, heig, tang, beta;
66 Real smll = std::fabs(ss[j][k]);
67 if(ite> 5 &&
68 smll<epsPrec*std::fabs(diagonal_[j]) &&
69 smll<epsPrec*std::fabs(diagonal_[k])) {
70 ss[j][k] = 0;
71 } else if (std::fabs(ss[j][k])>threshold) {
72 heig = diagonal_[k]-diagonal_[j];
73 if (smll<epsPrec*std::fabs(heig)) {
74 tang = ss[j][k]/heig;
75 } else {
76 beta = 0.5*heig/ss[j][k];
77 tang = 1.0/(std::fabs(beta)+
78 std::sqrt(1+beta*beta));
79 if (beta<0)
80 tang = -tang;
81 }
82 cosin = 1/std::sqrt(1+tang*tang);
83 sine = tang*cosin;
84 rho = sine/(1+cosin);
85 heig = tang*ss[j][k];
86 tmpAccumulate[j] -= heig;
87 tmpAccumulate[k] += heig;
88 diagonal_[j] -= heig;
89 diagonal_[k] += heig;
90 ss[j][k] = 0.0;
91 for (l=0; l+1<=j; l++)
92 jacobiRotate_(ss, rho, sine, l, j, l, k);
93 for (l=j+1; l<=k-1; l++)
94 jacobiRotate_(ss, rho, sine, j, l, l, k);
95 for (l=k+1; l<size; l++)
96 jacobiRotate_(ss, rho, sine, j, l, k, l);
97 for (l=0; l<size; l++)
99 rho, sine, l, j, l, k);
100 }
101 }
102 }
103 for (k=0; k<size; k++) {
104 tmpDiag[k] += tmpAccumulate[k];
105 diagonal_[k] = tmpDiag[k];
106 tmpAccumulate[k] = 0.0;
107 }
108 }
109 } while (++ite<=maxIterations && keeplooping);
110
111 QL_ENSURE(ite<=maxIterations,
112 "Too many iterations (" << maxIterations << ") reached");
113
114
115 // sort (eigenvalues, eigenvectors)
116 std::vector<std::pair<Real, std::vector<Real> > > temp(size);
117 std::vector<Real> eigenVector(size);
118 Size row, col;
119 for (col=0; col<size; col++) {
120 std::copy(eigenVectors_.column_begin(col),
121 eigenVectors_.column_end(col), eigenVector.begin());
122 temp[col] = std::make_pair(diagonal_[col], eigenVector);
123 }
124 std::sort(temp.begin(), temp.end(), std::greater<>());
125 Real maxEv = temp[0].first;
126 for (col=0; col<size; col++) {
127 // check for round-off errors
128 diagonal_[col] =
129 (std::fabs(temp[col].first/maxEv)<1e-16 ? 0.0 :
130 temp[col].first);
131 Real sign = 1.0;
132 if (temp[col].second[0]<0.0)
133 sign = -1.0;
134 for (row=0; row<size; row++) {
135 eigenVectors_[row][col] = sign * temp[col].second[row];
136 }
137 }
138 }
139
140}
const_iterator end() const
Definition: array.hpp:511
const_iterator begin() const
Definition: array.hpp:503
Matrix used in linear algebra.
Definition: matrix.hpp:41
const_column_iterator column_begin(Size i) const
Definition: matrix.hpp:415
const_column_iterator column_end(Size i) const
Definition: matrix.hpp:434
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
Definition: errors.hpp:130
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
ext::function< Real(Real)> b
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Real rho
Definition: any.hpp:35
ext::shared_ptr< YieldTermStructure > q
Real beta
Definition: sabr.cpp:200
Eigenvalues/eigenvectors of a real symmetric matrix.