QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
tridiagonaloperator.hpp
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) 2000, 2001, 2002, 2003 RiskMap srl
5 Copyright (C) 2003, 2004, 2005, 2006 StatPro Italia srl
6 Copyright (C) 2011 Ferdinando Ametrano
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22/*! \file tridiagonaloperator.hpp
23 \brief tridiagonal operator
24*/
25
26#ifndef quantlib_tridiagonal_operator_hpp
27#define quantlib_tridiagonal_operator_hpp
28
29#include <ql/math/array.hpp>
31#include <ql/shared_ptr.hpp>
32
33namespace QuantLib {
34
35 //! Base implementation for tridiagonal operator
36 /*! \warning to use real time-dependant algebra, you must overload
37 the corresponding operators in the inheriting
38 time-dependent class.
39
40 \ingroup findiff
41 */
43 // unary operators
46 // binary operators
48 const TridiagonalOperator&);
50 const TridiagonalOperator&);
52 const TridiagonalOperator&);
54 Real);
56 Real);
57 public:
59 // constructors
60 explicit TridiagonalOperator(Size size = 0);
61 TridiagonalOperator(const Array& low,
62 const Array& mid,
63 const Array& high);
67 TridiagonalOperator& operator=(TridiagonalOperator&&) noexcept;
68 ~TridiagonalOperator() = default;
69 //! \name Operator interface
70 //@{
71 //! apply operator to a given array
72 Array applyTo(const Array& v) const;
73 //! solve linear system for a given right-hand side
74 Array solveFor(const Array& rhs) const;
75 /*! solve linear system for a given right-hand side
76 without result Array allocation. The rhs and result parameters
77 can be the same Array, in which case rhs will be changed
78 */
79 void solveFor(const Array& rhs,
80 Array& result) const;
81 //! solve linear system with SOR approach
82 Array SOR(const Array& rhs,
83 Real tol) const;
84 //! identity instance
86 //@}
87 //! \name Inspectors
88 //@{
89 Size size() const { return n_; }
90 bool isTimeDependent() const { return !!timeSetter_; }
91 const Array& lowerDiagonal() const { return lowerDiagonal_; }
92 const Array& diagonal() const { return diagonal_; }
93 const Array& upperDiagonal() const { return upperDiagonal_; }
94 //@}
95 //! \name Modifiers
96 //@{
97 void setFirstRow(Real, Real);
98 void setMidRow(Size, Real, Real, Real);
99 void setMidRows(Real, Real, Real);
100 void setLastRow(Real, Real);
101 void setTime(Time t);
102 //@}
103 //! \name Utilities
104 //@{
105 void swap(TridiagonalOperator&) noexcept;
106 //@}
107 //! encapsulation of time-setting logic
109 public:
110 virtual ~TimeSetter() = default;
111 virtual void setTime(Time t,
112 TridiagonalOperator& L) const = 0;
113 };
114 protected:
117 mutable Array temp_;
118 ext::shared_ptr<TimeSetter> timeSetter_;
119 };
120
121 /* \relates TridiagonalOperator */
123
124
125 // inline definitions
126
128 swap(from);
129 }
130
132 const TridiagonalOperator& from) {
133 TridiagonalOperator temp(from);
134 swap(temp);
135 return *this;
136 }
137
138 inline TridiagonalOperator&
140 swap(from);
141 return *this;
142 }
143
145 Real valC) {
146 diagonal_[0] = valB;
147 upperDiagonal_[0] = valC;
148 }
149
151 Real valA,
152 Real valB,
153 Real valC) {
154 QL_REQUIRE(i>=1 && i<=n_-2,
155 "out of range in TridiagonalSystem::setMidRow");
156 lowerDiagonal_[i-1] = valA;
157 diagonal_[i] = valB;
158 upperDiagonal_[i] = valC;
159 }
160
162 Real valB,
163 Real valC) {
164 for (Size i=1; i<=n_-2; i++) {
165 lowerDiagonal_[i-1] = valA;
166 diagonal_[i] = valB;
167 upperDiagonal_[i] = valC;
168 }
169 }
170
172 Real valB) {
173 lowerDiagonal_[n_-2] = valA;
174 diagonal_[n_-1] = valB;
175 }
176
178 if (timeSetter_ != nullptr)
179 timeSetter_->setTime(t, *this);
180 }
181
182 inline void TridiagonalOperator::swap(TridiagonalOperator& from) noexcept {
183 std::swap(n_, from.n_);
184 diagonal_.swap(from.diagonal_);
185 lowerDiagonal_.swap(from.lowerDiagonal_);
186 upperDiagonal_.swap(from.upperDiagonal_);
187 temp_.swap(from.temp_);
188 timeSetter_.swap(from.timeSetter_);
189 }
190
191
192 // Time constant algebra
193
195 TridiagonalOperator D1 = D;
196 return D1;
197 }
198
200 Array low = -D.lowerDiagonal_,
201 mid = -D.diagonal_,
202 high = -D.upperDiagonal_;
203 TridiagonalOperator result(low, mid, high);
204 return result;
205 }
206
208 const TridiagonalOperator& D2) {
209 Array low = D1.lowerDiagonal_ + D2.lowerDiagonal_,
210 mid = D1.diagonal_ + D2.diagonal_,
211 high = D1.upperDiagonal_ + D2.upperDiagonal_;
212 TridiagonalOperator result(low, mid, high);
213 return result;
214 }
215
217 const TridiagonalOperator& D2) {
218 Array low = D1.lowerDiagonal_ - D2.lowerDiagonal_,
219 mid = D1.diagonal_ - D2.diagonal_,
220 high = D1.upperDiagonal_ - D2.upperDiagonal_;
221 TridiagonalOperator result(low, mid, high);
222 return result;
223 }
224
226 const TridiagonalOperator& D) {
227 Array low = D.lowerDiagonal_ * a,
228 mid = D.diagonal_ * a,
229 high = D.upperDiagonal_ * a;
230 TridiagonalOperator result(low, mid, high);
231 return result;
232 }
233
235 Real a) {
236 Array low = D.lowerDiagonal_ * a,
237 mid = D.diagonal_ * a,
238 high = D.upperDiagonal_ * a;
239 TridiagonalOperator result(low, mid, high);
240 return result;
241 }
242
244 Real a) {
245 Array low = D.lowerDiagonal_ / a,
246 mid = D.diagonal_ / a,
247 high = D.upperDiagonal_ / a;
248 TridiagonalOperator result(low, mid, high);
249 return result;
250 }
251
252 inline void swap(TridiagonalOperator& L1,
253 TridiagonalOperator& L2) noexcept {
254 L1.swap(L2);
255 }
256
257}
258
259#endif
1-D array used in linear algebra.
1-D array used in linear algebra.
Definition: array.hpp:52
encapsulation of time-setting logic
virtual void setTime(Time t, TridiagonalOperator &L) const =0
Base implementation for tridiagonal operator.
ext::shared_ptr< TimeSetter > timeSetter_
friend TridiagonalOperator operator-(const TridiagonalOperator &)
friend TridiagonalOperator operator+(const TridiagonalOperator &)
Array solveFor(const Array &rhs) const
solve linear system for a given right-hand side
void swap(TridiagonalOperator &) noexcept
Array applyTo(const Array &v) const
apply operator to a given array
Array SOR(const Array &rhs, Real tol) const
solve linear system with SOR approach
const Array & lowerDiagonal() const
const Array & upperDiagonal() const
static TridiagonalOperator identity(Size size)
identity instance
friend TridiagonalOperator operator*(Real, const TridiagonalOperator &)
TridiagonalOperator & operator=(const TridiagonalOperator &)
void setMidRow(Size, Real, Real, Real)
friend TridiagonalOperator operator/(const TridiagonalOperator &, Real)
TridiagonalOperator(const TridiagonalOperator &)=default
floating-point comparisons
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
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
Quantity operator-(const Quantity &m1, const Quantity &m2)
Definition: quantity.hpp:171
Quantity operator*(const Quantity &m, Real x)
Definition: quantity.hpp:177
Quantity operator+(const Quantity &m1, const Quantity &m2)
Definition: quantity.hpp:165
void swap(Array &v, Array &w) noexcept
Definition: array.hpp:903
Real operator/(const Quantity &m1, const Quantity &m2)
Definition: quantity.cpp:86
ext::shared_ptr< BlackVolTermStructure > v
bool lowerDiagonal_
Definition: pseudosqrt.cpp:77
Maps shared_ptr to either the boost or std implementation.