ROL_SemismoothNewtonProjection.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef ROL_SEMISMOOTHNEWTONPROJECTION_H
#define ROL_SEMISMOOTHNEWTONPROJECTION_H

#include "ROL_PolyhedralProjection.hpp"
#include "ROL_ParameterList.hpp"
#include "ROL_KrylovFactory.hpp"

namespace ROL {

template<typename Real>
class SemismoothNewtonProjection : public PolyhedralProjection<Real> {
private:
  int dim_;
  Ptr<Krylov<Real>> krylov_;
  Ptr<Vector<Real>> xnew_, lnew_, dlam_;

  Real DEFAULT_atol_, DEFAULT_rtol_, DEFAULT_stol_;
  Real DEFAULT_decr_, DEFAULT_factor_, DEFAULT_regscale_, DEFAULT_errscale_;
  int DEFAULT_maxit_, DEFAULT_lstype_, DEFAULT_verbosity_;
  bool DEFAULT_useproj_;

  Real atol_, rtol_, stol_, decr_, factor_, regscale_, errscale_, ctol_;
  int maxit_, lstype_, verbosity_;
  bool useproj_;

  using PolyhedralProjection<Real>::bnd_;
  using PolyhedralProjection<Real>::con_;
  using PolyhedralProjection<Real>::xprim_;
  using PolyhedralProjection<Real>::xdual_;
  using PolyhedralProjection<Real>::mul_;
  using PolyhedralProjection<Real>::res_;

public:

  SemismoothNewtonProjection(const Vector<Real>               &xprim,
                             const Vector<Real>               &xdual,
                             const Ptr<BoundConstraint<Real>> &bnd,
                             const Ptr<Constraint<Real>>      &con,
                             const Vector<Real>               &mul,
                             const Vector<Real>               &res);

  SemismoothNewtonProjection(const Vector<Real>               &xprim,
                             const Vector<Real>               &xdual,
                             const Ptr<BoundConstraint<Real>> &bnd,
                             const Ptr<Constraint<Real>>      &con,
                             const Vector<Real>               &mul,
                             const Vector<Real>               &res,
                             ParameterList                    &list);

  void project(Vector<Real> &x, std::ostream &stream = std::cout) override;

private:

  // Jv = inv(A DP(y) A^T) v
  // y  = P(x + A^T lam)
  class Jacobian : public LinearOperator<Real> {
    private:
      const Ptr<Constraint<Real>>      con_;
      const Ptr<BoundConstraint<Real>> bnd_;
      const Ptr<const Vector<Real>>    y_;
      const Ptr<Vector<Real>> xdual_, xprim_;
      const Real alpha_;
    public:
      Jacobian(const Ptr<Constraint<Real>>      &con,
               const Ptr<BoundConstraint<Real>> &bnd,
               const Ptr<const Vector<Real>>    &y,
               const Ptr<Vector<Real>>          &xdual,
               const Ptr<Vector<Real>>          &xprim,
               const Real                        alpha = 1e-4)
        : con_(con), bnd_(bnd), y_(y), xdual_(xdual), xprim_(xprim), alpha_(alpha) {}
      void apply(Vector<Real> &Jx, const Vector<Real> &x, Real &tol) const {
        con_->applyAdjointJacobian(*xdual_,x.dual(),*y_,tol);
        xprim_->set(xdual_->dual());
        bnd_->pruneActive(*xprim_,*y_);
        con_->applyJacobian(Jx,*xprim_,*y_,tol);
        // This is a hack to make the Jacobian invertible
        Jx.axpy(alpha_,x);
      }
  };

  class Precond : public LinearOperator<Real> {
  private:
    const Real alpha_;
  public:
    Precond(Real alpha = 1e-4) : alpha_(alpha) {}
    void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
      const Real one(1);
      Hv.set(v.dual());
      Hv.scale(one+alpha_);
    }
    void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
      const Real one(1);
      Hv.set(v.dual());
      Hv.scale(one/(one+alpha_));
    }
  };

  Real residual(Vector<Real> &r, const Vector<Real> &y) const;

  void solve_newton_system(Vector<Real>       &s,
                           const Vector<Real> &r,
                           const Vector<Real> &y,
                           const Real          mu,
                           const Real          rho,
                           int                &iter,
                           int                &flag) const;

  void update_primal(Vector<Real>       &y,
                     const Vector<Real> &x,
                     const Vector<Real> &lam) const;

  void project_ssn(Vector<Real> &x,
                   Vector<Real> &lam,
                   Vector<Real> &dlam,
                   std::ostream &stream = std::cout) const;

   Real compute_tolerance() const;

}; // class SemismoothNewtonProjection

} // namespace ROL

#include "ROL_SemismoothNewtonProjection_Def.hpp"

#endif