QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
impliciteulerscheme.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) 2009 Andreas Gaida
5 Copyright (C) 2009 Ralph Schreyer
6 Copyright (C) 2009, 2017 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
22#include <ql/functional.hpp>
26#include <utility>
27
28namespace QuantLib {
29
30 ImplicitEulerScheme::ImplicitEulerScheme(ext::shared_ptr<FdmLinearOpComposite> map,
31 const bc_set& bcSet,
32 Real relTol,
33 SolverType solverType)
34 : dt_(Null<Real>()), iterations_(ext::make_shared<Size>(0U)), relTol_(relTol),
35 map_(std::move(map)), bcSet_(bcSet), solverType_(solverType) {}
36
38 return r - (theta*dt_)*map_->apply(r);
39 }
40
42 step(a, t, 1.0);
43 }
44
46 QL_REQUIRE(t-dt_ > -1e-8, "a step towards negative time given");
47 map_->setTime(std::max(0.0, t-dt_), t);
48 bcSet_.setTime(std::max(0.0, t-dt_));
49
51
52 if (map_->size() == 1) {
53 a = map_->solve_splitting(0, a, -theta*dt_);
54 }
55 else {
56 auto preconditioner = [&](const Array& _a){ return map_->preconditioner(_a, -theta*dt_); };
57 auto applyF = [&](const Array& _a){ return apply(_a, theta); };
58
59 if (solverType_ == BiCGstab) {
60 const BiCGStabResult result =
61 QuantLib::BiCGstab(applyF, std::max(Size(10), a.size()),
62 relTol_, preconditioner).solve(a, a);
63
64 (*iterations_) += result.iterations;
65 a = result.x;
66 }
67 else if (solverType_ == GMRES) {
68 const GMRESResult result =
69 QuantLib::GMRES(applyF, std::max(Size(10), a.size() / 10U), relTol_,
70 preconditioner)
71 .solve(a, a);
72
73 (*iterations_) += result.errors.size();
74 a = result.x;
75 }
76 else
77 QL_FAIL("unknown/illegal solver type");
78 }
80 }
81
83 dt_=dt;
84 }
85
87 return *iterations_;
88 }
89}
Biconjugate gradient stabilized method.
1-D array used in linear algebra.
Definition: array.hpp:52
BiCGStabResult solve(const Array &b, const Array &x0=Array()) const
Definition: bicgstab.cpp:37
void applyBeforeSolving(operator_type &op, array_type &a) const
GMRESResult solve(const Array &b, const Array &x0=Array()) const
Definition: gmres.cpp:39
ext::shared_ptr< Size > iterations_
ImplicitEulerScheme(ext::shared_ptr< FdmLinearOpComposite > map, const bc_set &bcSet=bc_set(), Real relTol=1e-8, SolverType solverType=BiCGstab)
const BoundaryConditionSchemeHelper bcSet_
const ext::shared_ptr< FdmLinearOpComposite > map_
Array apply(const Array &r, Real theta) const
void step(array_type &a, Time t)
template class providing a null value for a given type.
Definition: null.hpp:76
const DefaultType & t
#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
Maps function, bind and cref to either the boost or std implementation.
generalized minimal residual method
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Real theta
Implicit-Euler scheme.
Definition: any.hpp:35
STL namespace.
ext::shared_ptr< YieldTermStructure > r
std::list< Real > errors
Definition: gmres.hpp:50