ROL
ROL_NullSpaceOperator.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_NULL_SPACE_OPERATOR_H
11 #define ROL_NULL_SPACE_OPERATOR_H
12 
13 #include "ROL_LinearOperator.hpp"
14 #include "ROL_Constraint.hpp"
15 
23 namespace ROL {
24 
25 template <class Real>
26 class NullSpaceOperator : public LinearOperator<Real> {
27 private:
28  const Ptr<Constraint<Real>> con_;
29  const Ptr<Vector<Real>> x_;
30  const bool useAugSys_;
31 
32  mutable Ptr<Vector<Real>> b1_;
33  mutable Ptr<Vector<Real>> b1dual_;
34  mutable Ptr<Vector<Real>> b2_;
35  mutable Ptr<Vector<Real>> mul_;
36 
37  int dim_;
38  Real b1sqr_;
39 
40 public:
41  virtual ~NullSpaceOperator() {}
43  const Vector<Real> &dom,
44  const Vector<Real> &ran,
45  const bool useAugSys = false)
46  : con_(con), x_(dom.clone()), useAugSys_(useAugSys) {
47  x_->set(dom);
48  dim_ = ran.dimension();
49  if (dim_==1 && !useAugSys_) {
50  Real tol = std::sqrt(ROL_EPSILON<Real>());
51  b1_ = dom.dual().clone();
52  b1dual_ = dom.clone();
53  b2_ = ran.dual().clone(); b2_->setScalar(1.0);
54  con_->applyAdjointJacobian(*b1_,*b2_,dom,tol);
55  b1dual_->set(b1_->dual());
56  b1sqr_ = b1_->dot(*b1_);
57  }
58  else {
59  b1_ = dom.dual().clone();
60  b2_ = ran.clone();
61  mul_ = ran.dual().clone();
62  }
63  }
64 
66  const Ptr<const Vector<Real>> &dom,
67  const Ptr<const Vector<Real>> &ran)
68  : NullSpaceOperator(con,*dom,*ran) {}
69 
70  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
71  x_->set(x);
72  if (dim_==1 && !useAugSys_) {
73  Real tol = std::sqrt(ROL_EPSILON<Real>());
74  con_->applyAdjointJacobian(*b1_,*b2_,x,tol);
75  b1dual_->set(b1_->dual());
76  b1sqr_ = b1_->dot(*b1_);
77  }
78  }
79 
80  virtual void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
81  if (dim_==1 && !useAugSys_) {
82  Real dot = v.dot(*b1dual_);
83  Hv.set(v);
84  Hv.axpy(-dot/b1sqr_,*b1dual_);
85  }
86  else {
87  b1_->set(v.dual()); b2_->zero();
88  con_->solveAugmentedSystem(Hv,*mul_,*b1_,*b2_,*x_,tol);
89  }
90  }
91 
92  void applyAdjoint( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
93  apply(Hv,v,tol);
94  }
95 
96  void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
97  throw Exception::NotImplemented(">>> NullSpaceOperator::applyInverse : Not Implemented!");
98  }
99 
100  void applyAdjointInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
101  throw Exception::NotImplemented(">>> NullSpaceOperator::applyAdjointInverse : Not Implemented!");
102  }
103 
104 }; // class NullSpaceOperator
105 
106 } // namespace ROL
107 
108 #endif
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:162
Ptr< Vector< Real > > b1_
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
Ptr< Vector< Real > > mul_
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual Real dot(const Vector &x) const =0
Compute where .
void applyAdjointInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply adjoint of the inverse linear operator.
Ptr< Vector< Real > > b1dual_
Projects on to the null space of a linear constraint.
const Ptr< Vector< Real > > x_
const Ptr< Constraint< Real > > con_
NullSpaceOperator(const Ptr< Constraint< Real >> &con, const Vector< Real > &dom, const Vector< Real > &ran, const bool useAugSys=false)
Provides the interface to apply a linear operator.
NullSpaceOperator(const Ptr< Constraint< Real >> &con, const Ptr< const Vector< Real >> &dom, const Ptr< const Vector< Real >> &ran)
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
void applyAdjoint(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply adjoint of linear operator.
Ptr< Vector< Real > > b2_
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update linear operator.
Defines the general constraint operator interface.