projectedgradientstep.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_PROJECTED_GRADIENT_STEP_HH
#define DUNE_SOLVERS_PROJECTED_GRADIENT_STEP_HH

/** \file
    \brief Computes the exact generalized Cauchy point for a possibly nonconvex,
           box-constraints quadratic problem.
*/

#include <memory>

#include <dune/common/shared_ptr.hh>

#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
#include <dune/solvers/common/boxconstraint.hh>

/*
   \brief Computes the exact generalized Cauchy point for a possibly nonconvex,
   box-constraints quadratic problem, which is the first local minimizer along the
   projected gradient path.
 */
template<class MatrixType, class VectorType>
class ProjectedGradientStep : public IterationStep<VectorType,Dune::BitSetVector<VectorType::block_type::dimension> >
{

    typedef typename VectorType::block_type VectorBlock;

    enum {blocksize = VectorBlock::dimension};

    typedef typename VectorType::field_type field_type;
    typedef IterationStep<VectorType,Dune::BitSetVector<blocksize> > Base;

public:

    //! Default constructor.  Doesn't init anything
    ProjectedGradientStep() {}

    //! Constructor that gets the quadratic energy functional x^T mat x - rhs^T*x, and the initial iterate
    ProjectedGradientStep(const MatrixType& mat, VectorType& x, const VectorType& rhs)
        : Base(x), rhs_(&rhs)
    {
        mat_     = Dune::stackobject_to_shared_ptr(mat);
    }

    /** \brief Compute the generalized Cauchy point along a given direction. This is sometimes needed in
     *         nonlinear programming where this direction then corresponds to the gradient of the nonlinear function
     *         and not the gradient of the quadratic approximation. This can also be used as a criticality measure.
     **/
    void computeGeneralizedCP(const VectorType& dir);

    //! Perform one iteration
    virtual void iterate() {
        VectorType negativeGradient=*rhs_;
        mat_->mmv(*this->x_,negativeGradient);
        computeGeneralizedCP(negativeGradient);
    }

    //! The obstacles
    const std::vector<BoxConstraint<field_type,blocksize> >* obstacles_;
private:
    //! The container for the right hand side
    const VectorType* rhs_;

    //! The linear operator
    std::shared_ptr<const MatrixType> mat_;
};

#include "projectedgradientstep.cc"

#endif