QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
ninepointlinearop.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) 2008 Andreas Gaida
5 Copyright (C) 2008 Ralph Schreyer
6 Copyright (C) 2008 Klaus Spanderen
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
25
26namespace QuantLib {
27
29 Size d0, Size d1,
30 const ext::shared_ptr<FdmMesher>& mesher)
31 : d0_(d0), d1_(d1),
32 i00_(new Size[mesher->layout()->size()]),
33 i10_(new Size[mesher->layout()->size()]),
34 i20_(new Size[mesher->layout()->size()]),
35 i01_(new Size[mesher->layout()->size()]),
36 i21_(new Size[mesher->layout()->size()]),
37 i02_(new Size[mesher->layout()->size()]),
38 i12_(new Size[mesher->layout()->size()]),
39 i22_(new Size[mesher->layout()->size()]),
40 a00_(new Real[mesher->layout()->size()]),
41 a10_(new Real[mesher->layout()->size()]),
42 a20_(new Real[mesher->layout()->size()]),
43 a01_(new Real[mesher->layout()->size()]),
44 a11_(new Real[mesher->layout()->size()]),
45 a21_(new Real[mesher->layout()->size()]),
46 a02_(new Real[mesher->layout()->size()]),
47 a12_(new Real[mesher->layout()->size()]),
48 a22_(new Real[mesher->layout()->size()]),
49 mesher_(mesher) {
50
52 && d0_ < mesher->layout()->dim().size()
53 && d1_ < mesher->layout()->dim().size(),
54 "inconsistent derivative directions");
55
56 for (const auto& iter : *mesher->layout()) {
57 const Size i = iter.index();
58
59 i10_[i] = mesher->layout()->neighbourhood(iter, d1_, -1);
60 i01_[i] = mesher->layout()->neighbourhood(iter, d0_, -1);
61 i21_[i] = mesher->layout()->neighbourhood(iter, d0_, 1);
62 i12_[i] = mesher->layout()->neighbourhood(iter, d1_, 1);
63 i00_[i] = mesher->layout()->neighbourhood(iter, d0_, -1, d1_, -1);
64 i20_[i] = mesher->layout()->neighbourhood(iter, d0_, 1, d1_, -1);
65 i02_[i] = mesher->layout()->neighbourhood(iter, d0_, -1, d1_, 1);
66 i22_[i] = mesher->layout()->neighbourhood(iter, d0_, 1, d1_, 1);
67 }
68 }
69
71 : i00_(new Size[m.mesher_->layout()->size()]),
72 i10_(new Size[m.mesher_->layout()->size()]),
73 i20_(new Size[m.mesher_->layout()->size()]),
74 i01_(new Size[m.mesher_->layout()->size()]),
75 i21_(new Size[m.mesher_->layout()->size()]),
76 i02_(new Size[m.mesher_->layout()->size()]),
77 i12_(new Size[m.mesher_->layout()->size()]),
78 i22_(new Size[m.mesher_->layout()->size()]),
79 a00_(new Real[m.mesher_->layout()->size()]),
80 a10_(new Real[m.mesher_->layout()->size()]),
81 a20_(new Real[m.mesher_->layout()->size()]),
82 a01_(new Real[m.mesher_->layout()->size()]),
83 a11_(new Real[m.mesher_->layout()->size()]),
84 a21_(new Real[m.mesher_->layout()->size()]),
85 a02_(new Real[m.mesher_->layout()->size()]),
86 a12_(new Real[m.mesher_->layout()->size()]),
87 a22_(new Real[m.mesher_->layout()->size()]),
88 mesher_(m.mesher_) {
89
90 const Size size = mesher_->layout()->size();
91 std::copy(m.i00_.get(), m.i00_.get()+size, i00_.get());
92 std::copy(m.i10_.get(), m.i10_.get()+size, i10_.get());
93 std::copy(m.i20_.get(), m.i20_.get()+size, i20_.get());
94 std::copy(m.i01_.get(), m.i01_.get()+size, i01_.get());
95 std::copy(m.i21_.get(), m.i21_.get()+size, i21_.get());
96 std::copy(m.i02_.get(), m.i02_.get()+size, i02_.get());
97 std::copy(m.i12_.get(), m.i12_.get()+size, i12_.get());
98 std::copy(m.i22_.get(), m.i22_.get()+size, i22_.get());
99 std::copy(m.a00_.get(), m.a00_.get()+size, a00_.get());
100 std::copy(m.a10_.get(), m.a10_.get()+size, a10_.get());
101 std::copy(m.a20_.get(), m.a20_.get()+size, a20_.get());
102 std::copy(m.a01_.get(), m.a01_.get()+size, a01_.get());
103 std::copy(m.a11_.get(), m.a11_.get()+size, a11_.get());
104 std::copy(m.a21_.get(), m.a21_.get()+size, a21_.get());
105 std::copy(m.a02_.get(), m.a02_.get()+size, a02_.get());
106 std::copy(m.a12_.get(), m.a12_.get()+size, a12_.get());
107 std::copy(m.a22_.get(), m.a22_.get()+size, a22_.get());
108 }
109
111
112 QL_REQUIRE(u.size() == mesher_->layout()->size(),"inconsistent length of r "
113 << u.size() << " vs " << mesher_->layout()->size());
114
115 Array retVal(u.size());
116 // direct access to make the following code faster.
117 const Real *a00(a00_.get()), *a01(a01_.get()), *a02(a02_.get());
118 const Real *a10(a10_.get()), *a11(a11_.get()), *a12(a12_.get());
119 const Real *a20(a20_.get()), *a21(a21_.get()), *a22(a22_.get());
120 const Size *i00(i00_.get()), *i01(i01_.get()), *i02(i02_.get());
121 const Size *i10(i10_.get()), *i12(i12_.get());
122 const Size *i20(i20_.get()), *i21(i21_.get()), *i22(i22_.get());
123
124 //#pragma omp parallel for
125 for (Size i=0; i < retVal.size(); ++i) {
126 retVal[i] = a00[i]*u[i00[i]]
127 + a01[i]*u[i01[i]]
128 + a02[i]*u[i02[i]]
129 + a10[i]*u[i10[i]]
130 + a11[i]*u[i]
131 + a12[i]*u[i12[i]]
132 + a20[i]*u[i20[i]]
133 + a21[i]*u[i21[i]]
134 + a22[i]*u[i22[i]];
135 }
136 return retVal;
137 }
138
140 const Size n = mesher_->layout()->size();
141
142 SparseMatrix retVal(n, n, 9*n);
143 for (Size i=0; i < mesher_->layout()->size(); ++i) {
144 retVal(i, i00_[i]) += a00_[i];
145 retVal(i, i01_[i]) += a01_[i];
146 retVal(i, i02_[i]) += a02_[i];
147 retVal(i, i10_[i]) += a10_[i];
148 retVal(i, i ) += a11_[i];
149 retVal(i, i12_[i]) += a12_[i];
150 retVal(i, i20_[i]) += a20_[i];
151 retVal(i, i21_[i]) += a21_[i];
152 retVal(i, i22_[i]) += a22_[i];
153 }
154
155 return retVal;
156 }
157
158
160
162 const Size size = mesher_->layout()->size();
163
164 //#pragma omp parallel for
165 for (Size i=0; i < size; ++i) {
166 const Real s = u[i];
167 retVal.a11_[i]=a11_[i]*s; retVal.a00_[i]=a00_[i]*s;
168 retVal.a01_[i]=a01_[i]*s; retVal.a02_[i]=a02_[i]*s;
169 retVal.a10_[i]=a10_[i]*s; retVal.a20_[i]=a20_[i]*s;
170 retVal.a21_[i]=a21_[i]*s; retVal.a12_[i]=a12_[i]*s;
171 retVal.a22_[i]=a22_[i]*s;
172 }
173
174 return retVal;
175 }
176
178 std::swap(d0_, m.d0_);
179 std::swap(d1_, m.d1_);
180
181 i00_.swap(m.i00_); i10_.swap(m.i10_); i20_.swap(m.i20_);
182 i01_.swap(m.i01_); i21_.swap(m.i21_); i02_.swap(m.i02_);
183 i12_.swap(m.i12_); i22_.swap(m.i22_);
184 a00_.swap(m.a00_); a10_.swap(m.a10_); a20_.swap(m.a20_);
185 a01_.swap(m.a01_); a21_.swap(m.a21_); a02_.swap(m.a02_);
186 a12_.swap(m.a12_); a22_.swap(m.a22_); a11_.swap(m.a11_);
187
188 mesher_.swap(m.mesher_);
189 }
190}
1-D array used in linear algebra.
Definition: array.hpp:52
Size size() const
dimension of the array
Definition: array.hpp:495
std::unique_ptr< Size[]> i12_
SparseMatrix toMatrix() const override
std::unique_ptr< Size[]> i20_
std::unique_ptr< Real[]> a20_
std::unique_ptr< Size[]> i02_
std::unique_ptr< Size[]> i21_
std::unique_ptr< Real[]> a11_
std::unique_ptr< Size[]> i01_
std::unique_ptr< Real[]> a00_
std::unique_ptr< Real[]> a02_
std::unique_ptr< Size[]> i10_
Array apply(const Array &r) const override
NinePointLinearOp mult(const Array &u) const
ext::shared_ptr< FdmMesher > mesher_
std::unique_ptr< Real[]> a22_
std::unique_ptr< Real[]> a10_
void swap(NinePointLinearOp &m) noexcept
std::unique_ptr< Real[]> a12_
std::unique_ptr< Real[]> a21_
std::unique_ptr< Size[]> i00_
std::unique_ptr< Size[]> i22_
std::unique_ptr< Real[]> a01_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
memory layout of a fdm linear operator
mesher for a fdm grid
const ext::shared_ptr< FdmMesher > mesher_
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
boost::numeric::ublas::compressed_matrix< Real > SparseMatrix
nine point linear operator