QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
triplebandlinearop.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) 2008 Andreas Gaida
5 Copyright (C) 2008 Ralph Schreyer
6 Copyright (C) 2008 Klaus Spanderen
7 Copyright (C) 2014 Johannes Göttker-Schnetmann
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
27
28namespace QuantLib {
29
31 Size direction,
32 const ext::shared_ptr<FdmMesher>& mesher)
33 : direction_(direction),
34 i0_ (new Size[mesher->layout()->size()]),
35 i2_ (new Size[mesher->layout()->size()]),
36 reverseIndex_ (new Size[mesher->layout()->size()]),
37 lower_ (new Real[mesher->layout()->size()]),
38 diag_ (new Real[mesher->layout()->size()]),
39 upper_ (new Real[mesher->layout()->size()]),
40 mesher_(mesher) {
41
42 std::vector<Size> newDim(mesher->layout()->dim());
43 std::iter_swap(newDim.begin(), newDim.begin()+direction_);
44 std::vector<Size> newSpacing = FdmLinearOpLayout(newDim).spacing();
45 std::iter_swap(newSpacing.begin(), newSpacing.begin()+direction_);
46
47 for (const auto& iter : *mesher->layout()) {
48 const Size i = iter.index();
49
50 i0_[i] = mesher->layout()->neighbourhood(iter, direction, -1);
51 i2_[i] = mesher->layout()->neighbourhood(iter, direction, 1);
52
53 const std::vector<Size>& coordinates = iter.coordinates();
54 const Size newIndex =
55 std::inner_product(coordinates.begin(), coordinates.end(),
56 newSpacing.begin(), Size(0));
57 reverseIndex_[newIndex] = i;
58 }
59 }
60
63 i0_ (new Size[m.mesher_->layout()->size()]),
64 i2_ (new Size[m.mesher_->layout()->size()]),
65 reverseIndex_(new Size[m.mesher_->layout()->size()]),
66 lower_(new Real[m.mesher_->layout()->size()]),
67 diag_ (new Real[m.mesher_->layout()->size()]),
68 upper_(new Real[m.mesher_->layout()->size()]),
69 mesher_(m.mesher_) {
70 const Size len = m.mesher_->layout()->size();
71 std::copy(m.i0_.get(), m.i0_.get() + len, i0_.get());
72 std::copy(m.i2_.get(), m.i2_.get() + len, i2_.get());
73 std::copy(m.reverseIndex_.get(), m.reverseIndex_.get()+len,
74 reverseIndex_.get());
75 std::copy(m.lower_.get(), m.lower_.get() + len, lower_.get());
76 std::copy(m.diag_.get(), m.diag_.get() + len, diag_.get());
77 std::copy(m.upper_.get(), m.upper_.get() + len, upper_.get());
78 }
79
81 mesher_.swap(m.mesher_);
82 std::swap(direction_, m.direction_);
83
84 i0_.swap(m.i0_); i2_.swap(m.i2_);
85 reverseIndex_.swap(m.reverseIndex_);
86 lower_.swap(m.lower_); diag_.swap(m.diag_); upper_.swap(m.upper_);
87 }
88
90 const TripleBandLinearOp& x,
91 const TripleBandLinearOp& y,
92 const Array& b) {
93 const Size size = mesher_->layout()->size();
94
95 Real *diag(diag_.get());
96 Real *lower(lower_.get());
97 Real *upper(upper_.get());
98
99 const Real *y_diag (y.diag_.get());
100 const Real *y_lower(y.lower_.get());
101 const Real *y_upper(y.upper_.get());
102
103 if (a.empty()) {
104 if (b.empty()) {
105 //#pragma omp parallel for
106 for (Size i=0; i < size; ++i) {
107 diag[i] = y_diag[i];
108 lower[i] = y_lower[i];
109 upper[i] = y_upper[i];
110 }
111 }
112 else {
113 Array::const_iterator bptr(b.begin());
114 const Size binc = (b.size() > 1) ? 1 : 0;
115 //#pragma omp parallel for
116 for (Size i=0; i < size; ++i) {
117 diag[i] = y_diag[i] + bptr[i*binc];
118 lower[i] = y_lower[i];
119 upper[i] = y_upper[i];
120 }
121 }
122 }
123 else if (b.empty()) {
125 const Size ainc = (a.size() > 1) ? 1 : 0;
126
127 const Real *x_diag (x.diag_.get());
128 const Real *x_lower(x.lower_.get());
129 const Real *x_upper(x.upper_.get());
130
131 //#pragma omp parallel for
132 for (Size i=0; i < size; ++i) {
133 const Real s = aptr[i*ainc];
134 diag[i] = y_diag[i] + s*x_diag[i];
135 lower[i] = y_lower[i] + s*x_lower[i];
136 upper[i] = y_upper[i] + s*x_upper[i];
137 }
138 }
139 else {
140 Array::const_iterator bptr(b.begin());
141 const Size binc = (b.size() > 1) ? 1 : 0;
142
144 const Size ainc = (a.size() > 1) ? 1 : 0;
145
146 const Real *x_diag (x.diag_.get());
147 const Real *x_lower(x.lower_.get());
148 const Real *x_upper(x.upper_.get());
149
150 //#pragma omp parallel for
151 for (Size i=0; i < size; ++i) {
152 const Real s = aptr[i*ainc];
153 diag[i] = y_diag[i] + s*x_diag[i] + bptr[i*binc];
154 lower[i] = y_lower[i] + s*x_lower[i];
155 upper[i] = y_upper[i] + s*x_upper[i];
156 }
157 }
158 }
159
161
163 const Size size = mesher_->layout()->size();
164 //#pragma omp parallel for
165 for (Size i=0; i < size; ++i) {
166 retVal.lower_[i]= lower_[i] + m.lower_[i];
167 retVal.diag_[i] = diag_[i] + m.diag_[i];
168 retVal.upper_[i]= upper_[i] + m.upper_[i];
169 }
170
171 return retVal;
172 }
173
174
176
178
179 const Size size = mesher_->layout()->size();
180 //#pragma omp parallel for
181 for (Size i=0; i < size; ++i) {
182 const Real s = u[i];
183 retVal.lower_[i]= lower_[i]*s;
184 retVal.diag_[i] = diag_[i]*s;
185 retVal.upper_[i]= upper_[i]*s;
186 }
187
188 return retVal;
189 }
190
192 const Size size = mesher_->layout()->size();
193 QL_REQUIRE(u.size() == size, "inconsistent size of rhs");
195
196 #pragma omp parallel for
197 for (long i=0; i < (long)size; ++i) {
198 const Real sm1 = i > 0? u[i-1] : 1.0;
199 const Real s0 = u[i];
200 const Real sp1 = i < (long)size-1? u[i+1] : 1.0;
201 retVal.lower_[i]= lower_[i]*sm1;
202 retVal.diag_[i] = diag_[i]*s0;
203 retVal.upper_[i]= upper_[i]*sp1;
204 }
205
206 return retVal;
207 }
208
210
212
213 const Size size = mesher_->layout()->size();
214 //#pragma omp parallel for
215 for (Size i=0; i < size; ++i) {
216 retVal.lower_[i]= lower_[i];
217 retVal.upper_[i]= upper_[i];
218 retVal.diag_[i] = diag_[i]+u[i];
219 }
220
221 return retVal;
222 }
223
225 QL_REQUIRE(r.size() == mesher_->layout()->size(), "inconsistent length of r");
226
227 const Real* lptr = lower_.get();
228 const Real* dptr = diag_.get();
229 const Real* uptr = upper_.get();
230 const Size* i0ptr = i0_.get();
231 const Size* i2ptr = i2_.get();
232
233 array_type retVal(r.size());
234 //#pragma omp parallel for
235 for (Size i=0; i < mesher_->layout()->size(); ++i) {
236 retVal[i] = r[i0ptr[i]]*lptr[i]+r[i]*dptr[i]+r[i2ptr[i]]*uptr[i];
237 }
238
239 return retVal;
240 }
241
243 const Size n = mesher_->layout()->size();
244
245 SparseMatrix retVal(n, n, 3*n);
246 for (Size i=0; i < n; ++i) {
247 retVal(i, i0_[i]) += lower_[i];
248 retVal(i, i ) += diag_[i];
249 retVal(i, i2_[i]) += upper_[i];
250 }
251
252 return retVal;
253 }
254
255
257 QL_REQUIRE(r.size() == mesher_->layout()->size(), "inconsistent size of rhs");
258
259#ifdef QL_EXTRA_SAFETY_CHECKS
260 for (const auto& iter : *mesher_->layout()) {
261 const std::vector<Size>& coordinates = iter.coordinates();
262 QL_REQUIRE( coordinates[direction_] != 0
263 || lower_[iter.index()] == 0,"removing non zero entry!");
264 QL_REQUIRE( coordinates[direction_] != mesher_->layout()->dim()[direction_]-1
265 || upper_[iter.index()] == 0,"removing non zero entry!");
266 }
267#endif
268
269 Array retVal(r.size()), tmp(r.size());
270
271 const Real* lptr = lower_.get();
272 const Real* dptr = diag_.get();
273 const Real* uptr = upper_.get();
274
275 // Thomson algorithm to solve a tridiagonal system.
276 // Example code taken from Tridiagonalopertor and
277 // changed to fit for the triple band operator.
278 Size rim1 = reverseIndex_[0];
279 Real bet=1.0/(a*dptr[rim1]+b);
280 QL_REQUIRE(bet != 0.0, "division by zero");
281 retVal[reverseIndex_[0]] = r[rim1]*bet;
282
283 for (Size j=1; j<=mesher_->layout()->size()-1; j++){
284 const Size ri = reverseIndex_[j];
285 tmp[j] = a*uptr[rim1]*bet;
286
287 bet=b+a*(dptr[ri]-tmp[j]*lptr[ri]);
288 QL_ENSURE(bet != 0.0, "division by zero");
289 bet=1.0/bet;
290
291 retVal[ri] = (r[ri]-a*lptr[ri]*retVal[rim1])*bet;
292 rim1 = ri;
293 }
294 // cannot be j>=0 with Size j
295 for (Size j=mesher_->layout()->size()-2; j>0; --j)
296 retVal[reverseIndex_[j]] -= tmp[j+1]*retVal[reverseIndex_[j+1]];
297 retVal[reverseIndex_[0]] -= tmp[1]*retVal[reverseIndex_[1]];
298
299 return retVal;
300 }
301}
1-D array used in linear algebra.
Definition: array.hpp:52
const Real * const_iterator
Definition: array.hpp:124
bool empty() const
whether the array is empty
Definition: array.hpp:499
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
const std::vector< Size > & spacing() const
std::unique_ptr< Real[]> diag_
void swap(TripleBandLinearOp &m) noexcept
std::unique_ptr< Real[]> lower_
std::unique_ptr< Size[]> i0_
SparseMatrix toMatrix() const override
TripleBandLinearOp multR(const Array &u) const
Array solve_splitting(const Array &r, Real a, Real b=1.0) const
std::unique_ptr< Real[]> upper_
void axpyb(const Array &a, const TripleBandLinearOp &x, const TripleBandLinearOp &y, const Array &b)
TripleBandLinearOp mult(const Array &u) const
Array apply(const Array &r) const override
ext::shared_ptr< FdmMesher > mesher_
std::unique_ptr< Size[]> reverseIndex_
TripleBandLinearOp add(const TripleBandLinearOp &m) const
std::unique_ptr< Size[]> i2_
#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
memory layout of a fdm linear operator
mesher for a fdm grid
const ext::shared_ptr< FdmMesher > mesher_
const Size direction_
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
boost::numeric::ublas::compressed_matrix< Real > SparseMatrix
ext::shared_ptr< YieldTermStructure > r
tridiagonal operator
general triple band linear operator