QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
gmres.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) 2017 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 license0/0 iee 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 gmres.cpp
21 \brief generalized minimal residual method
22*/
23
24
28#include <numeric>
29#include <utility>
30
31namespace QuantLib {
32
33 GMRES::GMRES(GMRES::MatrixMult A, Size maxIter, Real relTol, GMRES::MatrixMult preConditioner)
34 : A_(std::move(A)), M_(std::move(preConditioner)), maxIter_(maxIter), relTol_(relTol) {
35
36 QL_REQUIRE(maxIter_ > 0, "maxIter must be greater than zero");
37 }
38
39 GMRESResult GMRES::solve(const Array& b, const Array& x0) const {
40 GMRESResult result = solveImpl(b, x0);
41
42 QL_REQUIRE(result.errors.back() < relTol_, "could not converge");
43
44 return result;
45 }
46
48 Size restart, const Array& b, const Array& x0) const {
49
50 GMRESResult result = solveImpl(b, x0);
51
52 std::list<Real> errors = result.errors;
53
54 for (Size i=0; i < restart-1 && result.errors.back() >= relTol_;++i) {
55 result = solveImpl(b, result.x);
56 errors.insert(
57 errors.end(), result.errors.begin(), result.errors.end());
58 }
59
60 QL_REQUIRE(errors.back() < relTol_, "could not converge");
61
62 result.errors = errors;
63 return result;
64 }
65
66 GMRESResult GMRES::solveImpl(const Array& b, const Array& x0) const {
67 const Real bn = Norm2(b);
68 if (bn == 0.0) {
69 GMRESResult result = { std::list<Real>(1, 0.0), b };
70 return result;
71 }
72
73 Array x = ((!x0.empty()) ? x0 : Array(b.size(), 0.0));
74 Array r = b - A_(x);
75
76 const Real g = Norm2(r);
77 if (g/bn < relTol_) {
78 GMRESResult result = { std::list<Real>(1, g/bn), x };
79 return result;
80 }
81
82 std::vector<Array> v(1, r/g);
83 std::vector<Array> h(1, Array(maxIter_, 0.0));
84 std::vector<Real> c(maxIter_+1), s(maxIter_+1), z(maxIter_+1);
85
86 z[0] = g;
87
88 std::list<Real> errors(1, g/bn);
89
90 for (Size j=0; j < maxIter_ && errors.back() >= relTol_; ++j) {
91 h.emplace_back(maxIter_, 0.0);
92 Array w = A_(!M_ ? v[j] : M_(v[j]));
93
94 for (Size i=0; i <= j; ++i) {
95 h[i][j] = DotProduct(w, v[i]);
96 w -= h[i][j] * v[i];
97 }
98
99 h[j+1][j] = Norm2(w);
100
101 if (h[j+1][j] < QL_EPSILON*QL_EPSILON)
102 break;
103
104 v.push_back(w / h[j+1][j]);
105
106 for (Size i=0; i < j; ++i) {
107 const Real h0 = c[i]*h[i][j] + s[i]*h[i+1][j];
108 const Real h1 =-s[i]*h[i][j] + c[i]*h[i+1][j];
109
110 h[i][j] = h0;
111 h[i+1][j] = h1;
112 }
113
114 const Real nu = std::sqrt(squared(h[j][j]) + squared(h[j+1][j]));
115
116 c[j] = h[j][j]/nu;
117 s[j] = h[j+1][j]/nu;
118
119 h[j][j] = nu;
120 h[j+1][j] = 0.0;
121
122 z[j+1] = -s[j]*z[j];
123 z[j] = c[j] * z[j];
124
125 errors.push_back(std::fabs(z[j+1]/bn));
126 }
127
128 const Size k = v.size()-1;
129
130 Array y(k, 0.0);
131 y[k-1]=z[k-1]/h[k-1][k-1];
132
133 for (Integer i=k-2; i >= 0; --i) {
134 y[i] = (z[i] - std::inner_product(
135 h[i].begin()+i+1, h[i].begin()+k, y.begin()+i+1, Real(0.0)))/h[i][i];
136 }
137
138 Array xm = std::inner_product(
139 v.begin(), v.begin()+k, y.begin(), Array(x.size(), Real(0.0)));
140
141 xm = x + (!M_ ? xm : M_(xm));
142
143 GMRESResult result = { errors, xm };
144 return result;
145 }
146
147}
1-D array used in linear algebra.
Definition: array.hpp:52
bool empty() const
whether the array is empty
Definition: array.hpp:499
Size size() const
dimension of the array
Definition: array.hpp:495
GMRESResult solveWithRestart(Size restart, const Array &b, const Array &x0=Array()) const
Definition: gmres.cpp:47
const MatrixMult A_
Definition: gmres.hpp:67
GMRESResult solve(const Array &b, const Array &x0=Array()) const
Definition: gmres.cpp:39
const Real relTol_
Definition: gmres.hpp:69
const MatrixMult M_
Definition: gmres.hpp:67
GMRES(MatrixMult A, Size maxIter, Real relTol, MatrixMult preConditioner=MatrixMult())
Definition: gmres.cpp:33
const Size maxIter_
Definition: gmres.hpp:68
GMRESResult solveImpl(const Array &b, const Array &x0) const
Definition: gmres.cpp:66
ext::function< Array(const Array &)> MatrixMult
Definition: gmres.hpp:56
#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
generalized minimal residual method
#define QL_EPSILON
Definition: qldefines.hpp:178
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
functionals and combinators not included in the STL
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
Real Norm2(const Array &v)
Definition: array.hpp:560
Real DotProduct(const Array &v1, const Array &v2)
Definition: array.hpp:553
STL namespace.
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
QR decomposition.
Real nu
Definition: sabr.cpp:200
std::list< Real > errors
Definition: gmres.hpp:50