projectedgradientstep.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#include <dune/solvers/computeenergy.hh>

template<class MatrixType, class VectorType>
inline void ProjectedGradientStep<MatrixType, VectorType>::computeGeneralizedCP(const VectorType& dir)
{

    // Compute the first local minimum along the piecewise path of the projected direction
    VectorType projGrad = dir;

    ////////////////////////////////////////////////
    // First compute the projection of the direction
    // ////////////////////////////////////////////

    // Compute all the break points and sort them
    std::set<field_type> bPValues;
    field_type infinity = std::numeric_limits<field_type>::max();

    VectorType bPCoords(rhs_->size());
    bPCoords = 0;

    const VectorType& oldIt = *this->x_;

    for (size_t i=0; i<obstacles_->size(); i++)
        for (int j=0; j<blocksize; j++) {

            if (this->ignore()[i][j]) {
                projGrad[i][j]=0;
                continue;
            }

            bPCoords[i][j] = infinity;

            if (projGrad[i][j]>0 && (*obstacles_)[i].upper(j)<infinity)
                bPCoords[i][j] = ((*obstacles_)[i].upper(j)-oldIt[i][j])/projGrad[i][j];

            if (projGrad[i][j]<0 && (*obstacles_)[i].lower(j)>-infinity)
                bPCoords[i][j] = ((*obstacles_)[i].lower(j)-oldIt[i][j])/projGrad[i][j];

            if (bPCoords[i][j] <1e-14)
                projGrad[i][j] = 0;
            else if (bPCoords[i][j] < infinity)
                bPValues.insert(bPCoords[i][j]);
        }

    // Walk along the piecewise linear path until we find a minimizer
    VectorType cauchyPoint(rhs_->size());
    cauchyPoint = *this->x_;

    field_type oldBreakPoint(0);
    for ( auto bP : bPValues) {

        // if there are multiple breakpoints at the same position then skip them
        if (bP-oldBreakPoint<1e-14)
            continue;

        // Compute minimal value on the line segment
        VectorType dummy(projGrad.size());
        mat_->mv(projGrad,dummy);

        field_type linTerm = -((*rhs_)*projGrad) + (dummy*cauchyPoint);
        field_type localMin = -linTerm/(dummy*projGrad);

        if (dummy*projGrad <0)
            DUNE_THROW(Dune::Exception,"Not positive definite!\n");

        // Check if we really have a local minimum in this segment
        if (localMin < 0) {
            //field_type min = computeEnergy(*mat_,cauchyPoint,*rhs_);
            //std::cout<<"Generalized Cauchy Point "<<min<<" at previous breakpoint "<<oldBreakPoint<<std::endl;

            *this->x_ = cauchyPoint;

            return;

        } else if (localMin<=(bP-oldBreakPoint)) {

            cauchyPoint.axpy(localMin,projGrad);
            //field_type min = computeEnergy(*mat_,cauchyPoint,*rhs_);
            //std::cout<<"Generalized Cauchy Point found  "<<min<<" at local min "<<localMin
            //         <<" BP "<<bP+localMin<<std::endl;

            *this->x_ = cauchyPoint;

            return;
        }

        // Move to next interval
        cauchyPoint.axpy(bP-oldBreakPoint,projGrad);
        oldBreakPoint = bP;

        // Truncate projected gradient
        for (size_t i=0; i<bPCoords.size(); i++)
            for (int j=0; j<blocksize; j++)
                if (std::fabs(bPCoords[i][j]-bP)<1e-14)
                    projGrad[i][j] = 0;
    }
    *this->x_ += cauchyPoint;
}