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 
33 namespace 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 Disposable<TridiagonalOperator>
45         operator+(const TridiagonalOperator&);
46         friend Disposable<TridiagonalOperator>
47         operator-(const TridiagonalOperator&);
48         // binary operators
49         friend Disposable<TridiagonalOperator>
50         operator+(const TridiagonalOperator&,
51                   const TridiagonalOperator&);
52         friend Disposable<TridiagonalOperator>
53         operator-(const TridiagonalOperator&,
54                   const TridiagonalOperator&);
55         friend Disposable<TridiagonalOperator>
56         operator*(Real,
57                   const TridiagonalOperator&);
58         friend Disposable<TridiagonalOperator>
59         operator*(const TridiagonalOperator&,
60                   Real);
61         friend Disposable<TridiagonalOperator>
62         operator/(const TridiagonalOperator&,
63                   Real);
64       public:
65         typedef Array array_type;
66         // constructors
67         explicit TridiagonalOperator(Size size = 0);
68         TridiagonalOperator(const Array& low,
69                             const Array& mid,
70                             const Array& high);
71         TridiagonalOperator(const Disposable<TridiagonalOperator>&);
72         TridiagonalOperator& operator=(const Disposable<TridiagonalOperator>&);
73         //! \name Operator interface
74         //@{
75         //! apply operator to a given array
76         Disposable<Array> applyTo(const Array& v) const;
77         //! solve linear system for a given right-hand side
78         Disposable<Array> solveFor(const Array& rhs) const;
79         /*! solve linear system for a given right-hand side
80             without result Array allocation. The rhs and result parameters
81             can be the same Array, in which case rhs will be changed
82         */
83         void solveFor(const Array& rhs,
84                       Array& result) const;
85         //! solve linear system with SOR approach
86         Disposable<Array> SOR(const Array& rhs,
87                               Real tol) const;
88         //! identity instance
89         static Disposable<TridiagonalOperator> identity(Size size);
90         //@}
91         //! \name Inspectors
92         //@{
size() const93         Size size() const { return n_; }
isTimeDependent() const94         bool isTimeDependent() const { return !!timeSetter_; }
lowerDiagonal() const95         const Array& lowerDiagonal() const { return lowerDiagonal_; }
diagonal() const96         const Array& diagonal() const { return diagonal_; }
upperDiagonal() const97         const Array& upperDiagonal() const { return upperDiagonal_; }
98         //@}
99         //! \name Modifiers
100         //@{
101         void setFirstRow(Real, Real);
102         void setMidRow(Size, Real, Real, Real);
103         void setMidRows(Real, Real, Real);
104         void setLastRow(Real, Real);
105         void setTime(Time t);
106         //@}
107         //! \name Utilities
108         //@{
109         void swap(TridiagonalOperator&);
110         //@}
111         //! encapsulation of time-setting logic
112         class TimeSetter {
113           public:
~TimeSetter()114             virtual ~TimeSetter() {}
115             virtual void setTime(Time t,
116                                  TridiagonalOperator& L) const = 0;
117         };
118       protected:
119         Size n_;
120         Array diagonal_, lowerDiagonal_, upperDiagonal_;
121         mutable Array temp_;
122         ext::shared_ptr<TimeSetter> timeSetter_;
123     };
124 
125     /* \relates TridiagonalOperator */
126     void swap(TridiagonalOperator&, TridiagonalOperator&);
127 
128 
129     // inline definitions
130 
operator =(const Disposable<TridiagonalOperator> & from)131     inline TridiagonalOperator& TridiagonalOperator::operator=(
132                                 const Disposable<TridiagonalOperator>& from) {
133         swap(const_cast<Disposable<TridiagonalOperator>&>(from));
134         return *this;
135     }
136 
setFirstRow(Real valB,Real valC)137     inline void TridiagonalOperator::setFirstRow(Real valB,
138                                                  Real valC) {
139         diagonal_[0]      = valB;
140         upperDiagonal_[0] = valC;
141     }
142 
setMidRow(Size i,Real valA,Real valB,Real valC)143     inline void TridiagonalOperator::setMidRow(Size i,
144                                                Real valA,
145                                                Real valB,
146                                                Real valC) {
147         QL_REQUIRE(i>=1 && i<=n_-2,
148                    "out of range in TridiagonalSystem::setMidRow");
149         lowerDiagonal_[i-1] = valA;
150         diagonal_[i]        = valB;
151         upperDiagonal_[i]   = valC;
152     }
153 
setMidRows(Real valA,Real valB,Real valC)154     inline void TridiagonalOperator::setMidRows(Real valA,
155                                                 Real valB,
156                                                 Real valC) {
157         for (Size i=1; i<=n_-2; i++) {
158             lowerDiagonal_[i-1] = valA;
159             diagonal_[i]        = valB;
160             upperDiagonal_[i]   = valC;
161         }
162     }
163 
setLastRow(Real valA,Real valB)164     inline void TridiagonalOperator::setLastRow(Real valA,
165                                                 Real valB) {
166         lowerDiagonal_[n_-2] = valA;
167         diagonal_[n_-1]      = valB;
168     }
169 
setTime(Time t)170     inline void TridiagonalOperator::setTime(Time t) {
171         if (timeSetter_ != 0)
172             timeSetter_->setTime(t, *this);
173     }
174 
swap(TridiagonalOperator & from)175     inline void TridiagonalOperator::swap(TridiagonalOperator& from) {
176         using std::swap;
177         swap(n_, from.n_);
178         diagonal_.swap(from.diagonal_);
179         lowerDiagonal_.swap(from.lowerDiagonal_);
180         upperDiagonal_.swap(from.upperDiagonal_);
181         temp_.swap(from.temp_);
182         swap(timeSetter_, from.timeSetter_);
183     }
184 
185 
186     // Time constant algebra
187 
188     inline Disposable<TridiagonalOperator>
operator +(const TridiagonalOperator & D)189     operator+(const TridiagonalOperator& D) {
190         TridiagonalOperator D1 = D;
191         return D1;
192     }
193 
194     inline Disposable<TridiagonalOperator>
operator -(const TridiagonalOperator & D)195     operator-(const TridiagonalOperator& D) {
196         Array low = -D.lowerDiagonal_,
197             mid = -D.diagonal_,
198             high = -D.upperDiagonal_;
199         TridiagonalOperator result(low, mid, high);
200         return result;
201     }
202 
203     inline Disposable<TridiagonalOperator>
operator +(const TridiagonalOperator & D1,const TridiagonalOperator & D2)204     operator+(const TridiagonalOperator& D1,
205               const TridiagonalOperator& D2) {
206         Array low = D1.lowerDiagonal_ + D2.lowerDiagonal_,
207             mid = D1.diagonal_ + D2.diagonal_,
208             high = D1.upperDiagonal_ + D2.upperDiagonal_;
209         TridiagonalOperator result(low, mid, high);
210         return result;
211     }
212 
213     inline Disposable<TridiagonalOperator>
operator -(const TridiagonalOperator & D1,const TridiagonalOperator & D2)214     operator-(const TridiagonalOperator& D1,
215               const TridiagonalOperator& D2) {
216         Array low = D1.lowerDiagonal_ - D2.lowerDiagonal_,
217             mid = D1.diagonal_ - D2.diagonal_,
218             high = D1.upperDiagonal_ - D2.upperDiagonal_;
219         TridiagonalOperator result(low, mid, high);
220         return result;
221     }
222 
223     inline Disposable<TridiagonalOperator>
operator *(Real a,const TridiagonalOperator & D)224     operator*(Real a,
225               const TridiagonalOperator& D) {
226         Array low = D.lowerDiagonal_ * a,
227             mid = D.diagonal_ * a,
228             high = D.upperDiagonal_ * a;
229         TridiagonalOperator result(low, mid, high);
230         return result;
231     }
232 
233     inline Disposable<TridiagonalOperator>
operator *(const TridiagonalOperator & D,Real a)234     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 Disposable<TridiagonalOperator>
operator /(const TridiagonalOperator & D,Real a)244     operator/(const TridiagonalOperator& D,
245               Real a) {
246         Array low = D.lowerDiagonal_ / a,
247             mid = D.diagonal_ / a,
248             high = D.upperDiagonal_ / a;
249         TridiagonalOperator result(low, mid, high);
250         return result;
251     }
252 
swap(TridiagonalOperator & L1,TridiagonalOperator & L2)253     inline void swap(TridiagonalOperator& L1,
254                      TridiagonalOperator& L2) {
255         L1.swap(L2);
256     }
257 
258 }
259 
260 #endif
261