Logo Search packages:      
Sourcecode: quantlib version File versions  Download package

triplebandlinearop.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2008 Andreas Gaida
 Copyright (C) 2008 Ralph Schreyer
 Copyright (C) 2008 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/methods/finitedifferences/tridiagonaloperator.hpp>
#include <ql/experimental/finitedifferences/triplebandlinearop.hpp>

namespace QuantLib {

    TripleBandLinearOp::TripleBandLinearOp(
        Size direction,
        const boost::shared_ptr<FdmMesher>& mesher)
    : direction_(direction),
      i0_       (new Size[mesher->layout()->size()]),
      i2_       (new Size[mesher->layout()->size()]),
      reverseIndex_ (new Size[mesher->layout()->size()]),
      lower_    (new Real[mesher->layout()->size()]),
      diag_     (new Real[mesher->layout()->size()]),
      upper_    (new Real[mesher->layout()->size()]),
      mesher_(mesher) {

        const boost::shared_ptr<FdmLinearOpLayout> layout = mesher->layout();
        const FdmLinearOpIterator endIter = layout->end();

        std::vector<Size> newDim(layout->dim());
        std::iter_swap(newDim.begin(), newDim.begin()+direction_);
        std::vector<Size> newSpacing = FdmLinearOpLayout(newDim).spacing();
        std::iter_swap(newSpacing.begin(), newSpacing.begin()+direction_);

        for (FdmLinearOpIterator iter = layout->begin(); iter!=endIter; ++iter) {
            const Size i = iter.index();

            i0_[i] = layout->neighbourhood(iter, direction, -1);
            i2_[i] = layout->neighbourhood(iter, direction,  1);

            const std::vector<Size>& coordinates = iter.coordinates();
            const Size newIndex =
                  std::inner_product(coordinates.begin(), coordinates.end(),
                                     newSpacing.begin(), Size(0));
            reverseIndex_[newIndex] = i;

            //lowerEdge_[newIndex] = (coordinates[direction_] != 0);
            //upperEdge_[newIndex] =
            //    (coordinates[direction_] != layout->dim()[direction_]-1);
        }
    }

    TripleBandLinearOp::TripleBandLinearOp(const TripleBandLinearOp& m)
    : direction_(m.direction_),
      i0_   (new Size[m.mesher_->layout()->size()]),
      i2_   (new Size[m.mesher_->layout()->size()]),
      reverseIndex_(new Size[m.mesher_->layout()->size()]),
      lower_(new Real[m.mesher_->layout()->size()]),
      diag_ (new Real[m.mesher_->layout()->size()]),
      upper_(new Real[m.mesher_->layout()->size()]),
      mesher_(m.mesher_) {
        const Size len = m.mesher_->layout()->size();
        std::copy(m.i0_.get(), m.i0_.get() + len, i0_.get());
        std::copy(m.i2_.get(), m.i2_.get() + len, i2_.get());
        std::copy(m.reverseIndex_.get(), m.reverseIndex_.get()+len,
                  reverseIndex_.get());
        std::copy(m.lower_.get(), m.lower_.get() + len, lower_.get());
        std::copy(m.diag_.get(),  m.diag_.get() + len,  diag_.get());
        std::copy(m.upper_.get(), m.upper_.get() + len, upper_.get());
    }


    TripleBandLinearOp::TripleBandLinearOp(
        const Disposable<TripleBandLinearOp>& from) {
        swap(const_cast<Disposable<TripleBandLinearOp>&>(from));
    }

    TripleBandLinearOp& TripleBandLinearOp::operator=(
        const TripleBandLinearOp& m) {
        TripleBandLinearOp tmp(m);
        swap(tmp);
        return *this;
    }

    TripleBandLinearOp& TripleBandLinearOp::operator=(
        const Disposable<TripleBandLinearOp>& m) {
        swap(const_cast<Disposable<TripleBandLinearOp>&>(m));
        return *this;
    }

    void TripleBandLinearOp::swap(TripleBandLinearOp& m) {
        std::swap(mesher_, m.mesher_);
        std::swap(direction_, m.direction_);

        i0_.swap(m.i0_); i2_.swap(m.i2_);
        reverseIndex_.swap(m.reverseIndex_);
        lower_.swap(m.lower_); diag_.swap(m.diag_); upper_.swap(m.upper_);
    }

    void TripleBandLinearOp::axpyb(const Array& a,
                                   const TripleBandLinearOp& x,
                                   const TripleBandLinearOp& y,
                                   const Array& b) {
        const Size size = mesher_->layout()->size();

        Real *diag(diag_.get());
        Real *lower(lower_.get());
        Real *upper(upper_.get());

        const Real *y_diag (y.diag_.get());
        const Real *y_lower(y.lower_.get());
        const Real *y_upper(y.upper_.get());

        if (a.empty()) {
            if (b.empty()) {
                for (Size i=0; i < size; ++i) {
                    diag[i]  = y_diag[i];
                    lower[i] = y_lower[i];
                    upper[i] = y_upper[i];
                }
            }
            else {
                Array::const_iterator bptr(b.begin());
                const Size binc = (b.size() > 1) ? 1 : 0;

                for (Size i=0; i < size; ++i) {
                    diag[i]  = y_diag[i] + bptr[i*binc];
                    lower[i] = y_lower[i];
                    upper[i] = y_upper[i];
                }
            }
        }
        else if (b.empty()) {
            Array::const_iterator aptr(a.begin());
            const Size ainc = (a.size() > 1) ? 1 : 0;

            const Real *x_diag (x.diag_.get());
            const Real *x_lower(x.lower_.get());
            const Real *x_upper(x.upper_.get());

            for (Size i=0; i < size; ++i) {
                const Real s = aptr[i*ainc];
                diag[i]  = y_diag[i]  + s*x_diag[i];
                lower[i] = y_lower[i] + s*x_lower[i];
                upper[i] = y_upper[i] + s*x_upper[i];
            }
        }
        else {
            Array::const_iterator bptr(b.begin());
            const Size binc = (b.size() > 1) ? 1 : 0;

            Array::const_iterator aptr(a.begin());
            const Size ainc = (a.size() > 1) ? 1 : 0;

            const Real *x_diag (x.diag_.get());
            const Real *x_lower(x.lower_.get());
            const Real *x_upper(x.upper_.get());

            for (Size i=0; i < size; ++i) {
                const Real s = aptr[i*ainc];
                diag[i]  = y_diag[i]  + s*x_diag[i] + bptr[i*binc];
                lower[i] = y_lower[i] + s*x_lower[i];
                upper[i] = y_upper[i] + s*x_upper[i];
            }
        }
    }

    Disposable<TripleBandLinearOp>
    TripleBandLinearOp::add(const TripleBandLinearOp& m) const {

        TripleBandLinearOp retVal(direction_, mesher_);
        const Size size = mesher_->layout()->size();
        for (Size i=0; i < size; ++i) {
            retVal.lower_[i]= lower_[i] + m.lower_[i];
            retVal.diag_[i] = diag_[i]  + m.diag_[i];
            retVal.upper_[i]= upper_[i] + m.upper_[i];
        }

        return retVal;
    }


    Disposable<TripleBandLinearOp> TripleBandLinearOp::mult(const Array& u) const {

        TripleBandLinearOp retVal(direction_, mesher_);

        const Size size = mesher_->layout()->size();
        for (Size i=0; i < size; ++i) {
            const Real s = u[i];
            retVal.lower_[i]= lower_[i]*s;
            retVal.diag_[i] = diag_[i]*s;
            retVal.upper_[i]= upper_[i]*s;
        }

        return retVal;
    }

    Disposable<TripleBandLinearOp> TripleBandLinearOp::add(const Array& u) const {

        TripleBandLinearOp retVal(direction_, mesher_);

        const Size size = mesher_->layout()->size();
        for (Size i=0; i < size; ++i) {
            retVal.lower_[i]= lower_[i];
            retVal.upper_[i]= upper_[i];
            retVal.diag_[i] = diag_[i]+u[i];
        }

        return retVal;
    }

    Disposable<Array> TripleBandLinearOp::apply(const Array& r) const {
        const boost::shared_ptr<FdmLinearOpLayout> index = mesher_->layout();

        QL_REQUIRE(r.size() == index->size(), "inconsistent length of r");

        const Real* lptr = lower_.get();
        const Real* dptr = diag_.get();
        const Real* uptr = upper_.get();
        const Size* i0ptr = i0_.get();
        const Size* i2ptr = i2_.get();

        array_type retVal(r.size());
        for (Size i=0; i < index->size(); ++i) {
            retVal[i] = r[i0ptr[i]]*lptr[i]+r[i]*dptr[i]+r[i2ptr[i]]*uptr[i];
        }

        return retVal;
    }

    Disposable<Array>
    TripleBandLinearOp::solve_splitting(const Array& r, Real a, Real b) const {
        const boost::shared_ptr<FdmLinearOpLayout> layout = mesher_->layout();
        QL_REQUIRE(r.size() == layout->size(), "inconsistent size of rhs");

#ifdef QL_EXTRA_SAFETY_CHECKS
        for (FdmLinearOpIterator iter = layout->begin();
             iter!=layout->end(); ++iter) {
            const std::vector<Size>& coordinates = iter.coordinates();
            QL_REQUIRE(   coordinates[direction_] != 0
                       || lower_[iter.index()] == 0,"removing non zero entry!");
            QL_REQUIRE(   coordinates[direction_] != layout->dim()[direction_]-1
                       || upper_[iter.index()] == 0,"removing non zero entry!");
        }
#endif

        Array retVal(r.size()), tmp(r.size());

        const Real* lptr = lower_.get();
        const Real* dptr = diag_.get();
        const Real* uptr = upper_.get();

        // Thomson algorithm to solve a tridiagonal system.
        // Example code taken from Tridiagonalopertor and
        // changed to fit for the triple band operator.
        Size rim1 = reverseIndex_[0];
        Real bet=1.0/(a*dptr[rim1]+b);
        QL_REQUIRE(bet != 0.0, "division by zero");
        retVal[reverseIndex_[0]] = r[rim1]*bet;

        for (Size j=1; j<=layout->size()-1; j++){
            const Size ri = reverseIndex_[j];
            tmp[j] = a*uptr[rim1]*bet;

            bet=b+a*(dptr[ri]-tmp[j]*lptr[ri]);
            QL_ENSURE(bet != 0.0, "division by zero");
            bet=1.0/bet;

            retVal[ri] = (r[ri]-a*lptr[ri]*retVal[rim1])*bet;
            rim1 = ri;
        }
        // cannot be j>=0 with Size j
        for (Size j=layout->size()-2; j>0; --j)
            retVal[reverseIndex_[j]] -= tmp[j+1]*retVal[reverseIndex_[j+1]];
        retVal[reverseIndex_[0]] -= tmp[1]*retVal[reverseIndex_[1]];

        return retVal;
    }
}

Generated by  Doxygen 1.6.0   Back to index