QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
tridiagonaloperator.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) 2002, 2003, 2011 Ferdinando Ametrano
5 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
22
23namespace QuantLib {
24
26 if (size>=2) {
27 n_ = size;
31 temp_ = Array(size);
32 } else if (size==0) {
33 n_ = 0;
34 diagonal_ = Array(0);
37 temp_ = Array(0);
38 } else {
39 QL_FAIL("invalid size (" << size << ") for tridiagonal operator "
40 "(must be null or >= 2)");
41 }
42 }
43
45 const Array& mid,
46 const Array& high)
47 : n_(mid.size()),
48 diagonal_(mid), lowerDiagonal_(low), upperDiagonal_(high), temp_(n_) {
49 QL_REQUIRE(low.size() == n_-1,
50 "low diagonal vector of size " << low.size() <<
51 " instead of " << n_-1);
52 QL_REQUIRE(high.size() == n_-1,
53 "high diagonal vector of size " << high.size() <<
54 " instead of " << n_-1);
55 }
56
58 QL_REQUIRE(n_!=0,
59 "uninitialized TridiagonalOperator");
60 QL_REQUIRE(v.size()==n_,
61 "vector of the wrong size " << v.size() <<
62 " instead of " << n_);
63 Array result(n_);
64 std::transform(diagonal_.begin(), diagonal_.end(),
65 v.begin(),
66 result.begin(),
67 std::multiplies<>());
68
69 // matricial product
70 result[0] += upperDiagonal_[0]*v[1];
71 for (Size j=1; j<=n_-2; j++)
72 result[j] += lowerDiagonal_[j-1]*v[j-1]+
73 upperDiagonal_[j]*v[j+1];
74 result[n_-1] += lowerDiagonal_[n_-2]*v[n_-2];
75
76 return result;
77 }
78
80 Array result(rhs.size());
81 solveFor(rhs, result);
82 return result;
83 }
84
86 Array& result) const {
87
88 QL_REQUIRE(n_!=0,
89 "uninitialized TridiagonalOperator");
90 QL_REQUIRE(rhs.size()==n_,
91 "rhs vector of size " << rhs.size() <<
92 " instead of " << n_);
93
94 Real bet = diagonal_[0];
95 QL_REQUIRE(!close(bet, 0.0),
96 "diagonal's first element (" << bet <<
97 ") cannot be close to zero");
98 result[0] = rhs[0]/bet;
99 for (Size j=1; j<=n_-1; ++j) {
100 temp_[j] = upperDiagonal_[j-1]/bet;
101 bet = diagonal_[j]-lowerDiagonal_[j-1]*temp_[j];
102 QL_ENSURE(!close(bet, 0.0), "division by zero");
103 result[j] = (rhs[j] - lowerDiagonal_[j-1]*result[j-1])/bet;
104 }
105 // cannot be j>=0 with Size j
106 for (Size j=n_-2; j>0; --j)
107 result[j] -= temp_[j+1]*result[j+1];
108 result[0] -= temp_[1]*result[1];
109 }
110
112 Real tol) const {
113 QL_REQUIRE(n_!=0,
114 "uninitialized TridiagonalOperator");
115 QL_REQUIRE(rhs.size()==n_,
116 "rhs vector of size " << rhs.size() <<
117 " instead of " << n_);
118
119 // initial guess
120 Array result = rhs;
121
122 // solve tridiagonal system with SOR technique
123 Real omega = 1.5;
124 Real err = 2.0*tol;
125 Real temp;
126 for (Size sorIteration=0; err>tol ; ++sorIteration) {
127 QL_REQUIRE(sorIteration<100000,
128 "tolerance (" << tol << ") not reached in " <<
129 sorIteration << " iterations. " <<
130 "The error still is " << err);
131
132 temp = omega * (rhs[0] -
133 upperDiagonal_[0] * result[1]-
134 diagonal_[0] * result[0])/diagonal_[0];
135 err = temp*temp;
136 result[0] += temp;
137 Size i;
138 for (i=1; i<n_-1 ; ++i) {
139 temp = omega *(rhs[i] -
140 upperDiagonal_[i] * result[i+1]-
141 diagonal_[i] * result[i] -
142 lowerDiagonal_[i-1] * result[i-1])/diagonal_[i];
143 err += temp * temp;
144 result[i] += temp;
145 }
146
147 temp = omega * (rhs[i] -
148 diagonal_[i] * result[i] -
149 lowerDiagonal_[i-1] * result[i-1])/diagonal_[i];
150 err += temp*temp;
151 result[i] += temp;
152 }
153 return result;
154 }
155
157 return TridiagonalOperator(Array(size-1, 0.0), // lower diagonal
158 Array(size, 1.0), // diagonal
159 Array(size-1, 0.0)); // upper diagonal
160 }
161
162}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
Base implementation for tridiagonal operator.
Array solveFor(const Array &rhs) const
solve linear system for a given right-hand side
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
static TridiagonalOperator identity(Size size)
identity instance
#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
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
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
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
ext::shared_ptr< BlackVolTermStructure > v
bool lowerDiagonal_
Definition: pseudosqrt.cpp:77
tridiagonal operator