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