32 const ext::shared_ptr<FdmMesher>& mesher)
34 i0_ (new
Size[mesher->layout()->size()]),
35 i2_ (new
Size[mesher->layout()->size()]),
36 reverseIndex_ (new
Size[mesher->layout()->size()]),
37 lower_ (new
Real[mesher->layout()->size()]),
38 diag_ (new
Real[mesher->layout()->size()]),
39 upper_ (new
Real[mesher->layout()->size()]),
42 std::vector<Size> newDim(mesher->layout()->dim());
43 std::iter_swap(newDim.begin(), newDim.begin()+
direction_);
45 std::iter_swap(newSpacing.begin(), newSpacing.begin()+
direction_);
47 for (
const auto& iter : *mesher->layout()) {
48 const Size i = iter.index();
50 i0_[i] = mesher->layout()->neighbourhood(iter, direction, -1);
51 i2_[i] = mesher->layout()->neighbourhood(iter, direction, 1);
53 const std::vector<Size>& coordinates = iter.coordinates();
55 std::inner_product(coordinates.begin(), coordinates.end(),
56 newSpacing.begin(),
Size(0));
65 reverseIndex_(new
Size[m.
mesher_->layout()->size()]),
71 std::copy(m.
i0_.get(), m.
i0_.get() + len,
i0_.get());
72 std::copy(m.
i2_.get(), m.
i2_.get() + len,
i2_.get());
84 i0_.swap(m.i0_); i2_.swap(m.i2_);
85 reverseIndex_.swap(m.reverseIndex_);
86 lower_.swap(m.lower_); diag_.swap(m.diag_); upper_.swap(m.upper_);
99 const Real *y_diag (
y.diag_.get());
100 const Real *y_lower(
y.lower_.get());
101 const Real *y_upper(
y.upper_.get());
106 for (
Size i=0; i < size; ++i) {
108 lower[i] = y_lower[i];
109 upper[i] = y_upper[i];
114 const Size binc = (
b.size() > 1) ? 1 : 0;
116 for (
Size i=0; i < size; ++i) {
117 diag[i] = y_diag[i] + bptr[i*binc];
118 lower[i] = y_lower[i];
119 upper[i] = y_upper[i];
123 else if (
b.empty()) {
125 const Size ainc = (a.
size() > 1) ? 1 : 0;
132 for (
Size i=0; i < size; ++i) {
133 const Real s = aptr[i*ainc];
134 diag[i] = y_diag[i] +
s*x_diag[i];
135 lower[i] = y_lower[i] +
s*x_lower[i];
136 upper[i] = y_upper[i] +
s*x_upper[i];
141 const Size binc = (
b.size() > 1) ? 1 : 0;
144 const Size ainc = (a.
size() > 1) ? 1 : 0;
151 for (
Size i=0; i < size; ++i) {
152 const Real s = aptr[i*ainc];
153 diag[i] = y_diag[i] +
s*x_diag[i] + bptr[i*binc];
154 lower[i] = y_lower[i] +
s*x_lower[i];
155 upper[i] = y_upper[i] +
s*x_upper[i];
165 for (
Size i=0; i < size; ++i) {
181 for (
Size i=0; i < size; ++i) {
196 #pragma omp parallel for
197 for (
long i=0; i < (long)size; ++i) {
198 const Real sm1 = i > 0? u[i-1] : 1.0;
199 const Real s0 = u[i];
200 const Real sp1 = i < (long)size-1? u[i+1] : 1.0;
215 for (
Size i=0; i < size; ++i) {
235 for (
Size i=0; i <
mesher_->layout()->size(); ++i) {
236 retVal[i] =
r[i0ptr[i]]*lptr[i]+
r[i]*dptr[i]+
r[i2ptr[i]]*uptr[i];
246 for (
Size i=0; i <
n; ++i) {
248 retVal(i, i ) +=
diag_[i];
259#ifdef QL_EXTRA_SAFETY_CHECKS
260 for (
const auto& iter : *
mesher_->layout()) {
261 const std::vector<Size>& coordinates = iter.coordinates();
263 ||
lower_[iter.index()] == 0,
"removing non zero entry!");
265 ||
upper_[iter.index()] == 0,
"removing non zero entry!");
269 Array retVal(
r.size()), tmp(
r.size());
279 Real bet=1.0/(a*dptr[rim1]+
b);
283 for (
Size j=1; j<=
mesher_->layout()->size()-1; j++){
285 tmp[j] = a*uptr[rim1]*bet;
287 bet=
b+a*(dptr[ri]-tmp[j]*lptr[ri]);
288 QL_ENSURE(bet != 0.0,
"division by zero");
291 retVal[ri] = (
r[ri]-a*lptr[ri]*retVal[rim1])*bet;
1-D array used in linear algebra.
const Real * const_iterator
bool empty() const
whether the array is empty
Size size() const
dimension of the array
const_iterator begin() const
const std::vector< Size > & spacing() const
std::unique_ptr< Real[]> diag_
void swap(TripleBandLinearOp &m) noexcept
std::unique_ptr< Real[]> lower_
std::unique_ptr< Size[]> i0_
SparseMatrix toMatrix() const override
TripleBandLinearOp multR(const Array &u) const
Array solve_splitting(const Array &r, Real a, Real b=1.0) const
std::unique_ptr< Real[]> upper_
void axpyb(const Array &a, const TripleBandLinearOp &x, const TripleBandLinearOp &y, const Array &b)
TripleBandLinearOp mult(const Array &u) const
Array apply(const Array &r) const override
ext::shared_ptr< FdmMesher > mesher_
std::unique_ptr< Size[]> reverseIndex_
TripleBandLinearOp()=default
TripleBandLinearOp add(const TripleBandLinearOp &m) const
std::unique_ptr< Size[]> i2_
#define QL_ENSURE(condition, message)
throw an error if the given post-condition is not verified
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
ext::function< Real(Real)> b
memory layout of a fdm linear operator
const ext::shared_ptr< FdmMesher > mesher_
std::size_t Size
size of a container
boost::numeric::ublas::compressed_matrix< Real > SparseMatrix
ext::shared_ptr< YieldTermStructure > r
general triple band linear operator