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>
30#include <ql/math/comparison.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 */
42 class TridiagonalOperator {
43 // unary operators
44 friend TridiagonalOperator operator+(const TridiagonalOperator&);
45 friend TridiagonalOperator operator-(const TridiagonalOperator&);
46 // binary operators
47 friend TridiagonalOperator operator+(const TridiagonalOperator&,
48 const TridiagonalOperator&);
49 friend TridiagonalOperator operator-(const TridiagonalOperator&,
50 const TridiagonalOperator&);
51 friend TridiagonalOperator operator*(Real,
52 const TridiagonalOperator&);
53 friend TridiagonalOperator operator*(const TridiagonalOperator&,
54 Real);
55 friend TridiagonalOperator operator/(const TridiagonalOperator&,
56 Real);
57 public:
58 typedef Array array_type;
59 // constructors
60 explicit TridiagonalOperator(Size size = 0);
61 TridiagonalOperator(const Array& low,
62 const Array& mid,
63 const Array& high);
64 TridiagonalOperator(const TridiagonalOperator&) = default;
65 TridiagonalOperator(TridiagonalOperator&&) noexcept;
66 TridiagonalOperator& operator=(const TridiagonalOperator&);
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
85 static TridiagonalOperator identity(Size size);
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
108 class TimeSetter {
109 public:
110 virtual ~TimeSetter() = default;
111 virtual void setTime(Time t,
112 TridiagonalOperator& L) const = 0;
113 };
114 protected:
115 Size n_;
116 Array diagonal_, lowerDiagonal_, upperDiagonal_;
117 mutable Array temp_;
118 ext::shared_ptr<TimeSetter> timeSetter_;
119 };
120
121 /* \relates TridiagonalOperator */
122 void swap(TridiagonalOperator&, TridiagonalOperator&) noexcept;
123
124
125 // inline definitions
126
127 inline TridiagonalOperator::TridiagonalOperator(TridiagonalOperator&& from) noexcept : n_(0) {
128 swap(from);
129 }
130
131 inline TridiagonalOperator& TridiagonalOperator::operator=(
132 const TridiagonalOperator& from) {
133 TridiagonalOperator temp(from);
134 swap(temp);
135 return *this;
136 }
137
138 inline TridiagonalOperator&
139 TridiagonalOperator::operator=(TridiagonalOperator&& from) noexcept {
140 swap(from);
141 return *this;
142 }
143
144 inline void TridiagonalOperator::setFirstRow(Real valB,
145 Real valC) {
146 diagonal_[0] = valB;
147 upperDiagonal_[0] = valC;
148 }
149
150 inline void TridiagonalOperator::setMidRow(Size i,
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
161 inline void TridiagonalOperator::setMidRows(Real valA,
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
171 inline void TridiagonalOperator::setLastRow(Real valA,
172 Real valB) {
173 lowerDiagonal_[n_-2] = valA;
174 diagonal_[n_-1] = valB;
175 }
176
177 inline void TridiagonalOperator::setTime(Time t) {
178 if (timeSetter_ != nullptr)
179 timeSetter_->setTime(t, L&: *this);
180 }
181
182 inline void TridiagonalOperator::swap(TridiagonalOperator& from) noexcept {
183 std::swap(a&: n_, b&: from.n_);
184 diagonal_.swap(from&: from.diagonal_);
185 lowerDiagonal_.swap(from&: from.lowerDiagonal_);
186 upperDiagonal_.swap(from&: from.upperDiagonal_);
187 temp_.swap(from&: from.temp_);
188 timeSetter_.swap(other&: from.timeSetter_);
189 }
190
191
192 // Time constant algebra
193
194 inline TridiagonalOperator operator+(const TridiagonalOperator& D) {
195 TridiagonalOperator D1 = D;
196 return D1;
197 }
198
199 inline TridiagonalOperator operator-(const TridiagonalOperator& D) {
200 Array low = -D.lowerDiagonal_,
201 mid = -D.diagonal_,
202 high = -D.upperDiagonal_;
203 TridiagonalOperator result(low, mid, high);
204 return result;
205 }
206
207 inline TridiagonalOperator operator+(const TridiagonalOperator& D1,
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
216 inline TridiagonalOperator operator-(const TridiagonalOperator& D1,
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
225 inline TridiagonalOperator operator*(Real a,
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
234 inline TridiagonalOperator operator*(const TridiagonalOperator& D,
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
243 inline TridiagonalOperator operator/(const TridiagonalOperator& D,
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(from&: L2);
255 }
256
257}
258
259#endif
260

source code of quantlib/ql/methods/finitedifferences/tridiagonaloperator.hpp