45 #ifndef ROL_DYKSTRAPROJECTION_DEF_H
46 #define ROL_DYKSTRAPROJECTION_DEF_H
50 template<
typename Real>
60 DEFAULT_maxit_ (10000),
61 DEFAULT_verbosity_ (0),
62 atol_ (DEFAULT_atol_),
63 rtol_ (DEFAULT_rtol_),
64 maxit_ (DEFAULT_maxit_),
65 verbosity_ (DEFAULT_verbosity_) {
73 Real tol(std::sqrt(ROL_EPSILON<Real>()));
78 mul_->setScalar(static_cast<Real>(1));
79 con_->applyAdjointJacobian(*
z_,*
mul_,xprim,tol);
86 template<
typename Real>
95 atol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Absolute Tolerance",
DEFAULT_atol_);
96 rtol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Relative Tolerance",
DEFAULT_rtol_);
97 maxit_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Iteration Limit",
DEFAULT_maxit_);
101 template<
typename Real>
103 if (con_ == nullPtr) {
107 project_Dykstra(x, stream);
111 template<
typename Real>
113 return xprim_->dot(x) + b_;
116 template<
typename Real>
118 Real tol(std::sqrt(ROL_EPSILON<Real>()));
120 con_->value(r,y,tol);
123 template<
typename Real>
129 template<
typename Real>
132 Real rhs = residual_1d(y);
133 Real lam = -rhs/cdot_;
138 Real tol = std::sqrt(ROL_EPSILON<Real>());
139 residual_nd(*res_,y);
140 con_->solveAugmentedSystem(x,*mul_,*z_,*res_,y,tol);
141 x.
scale(static_cast<Real>(-1));
146 template<
typename Real>
148 const Real one(1), xnorm(x.
norm()), ctol(std::min(atol_,rtol_*xnorm));
149 Real norm1(0), norm2(0), rnorm(0);
150 p_->zero(); q_->zero();
151 std::ios_base::fmtflags streamFlags(stream.flags());
152 if (verbosity_ > 2) {
153 stream << std::scientific << std::setprecision(6);
155 stream <<
" Polyhedral Projection using Dykstra's Algorithm" << std::endl;
157 stream << std::setw(6) << std::left <<
"iter";
158 stream << std::setw(15) << std::left <<
"con norm";
159 stream << std::setw(15) << std::left <<
"bnd norm";
160 stream << std::setw(15) << std::left <<
"error";
161 stream << std::setw(15) << std::left <<
"tol";
164 for (
int cnt=0; cnt < maxit_; ++cnt) {
166 tmp_->set(x); tmp_->plus(*p_);
167 project_con(*y_,*tmp_);
168 p_->set(*tmp_); p_->axpy(-one,*y_);
170 tmp_->set(x); tmp_->axpy(-one,*y_);
171 norm1 = tmp_->norm();
173 tmp_->set(*y_); tmp_->plus(*q_);
174 project_bnd(x,*tmp_);
175 q_->set(*tmp_); q_->axpy(-one,x);
177 tmp_->set(x); tmp_->axpy(-one,*y_);
178 norm2 = tmp_->norm();
182 rnorm = std::sqrt(norm1*norm1 + norm2*norm2);
183 if (verbosity_ > 2) {
185 stream << std::setw(6) << std::left << cnt;
186 stream << std::setw(15) << std::left << norm1;
187 stream << std::setw(15) << std::left << norm2;
188 stream << std::setw(15) << std::left << rnorm;
189 stream << std::setw(15) << std::left << ctol;
192 if (rnorm <= ctol)
break;
194 if (verbosity_ > 2) {
199 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
200 stream << rnorm <<
" rtol = " << ctol << std::endl;
202 stream.flags(streamFlags);
Ptr< Vector< Real > > xprim_
virtual void scale(const Real alpha)=0
Compute where .
Real residual_1d(const Vector< Real > &x) const
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
Ptr< Vector< Real > > res_
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
void residual_nd(Vector< Real > &r, const Vector< Real > &y) const
Defines the linear algebra or vector space interface.
Ptr< Vector< Real > > tmp_
void project_con(Vector< Real > &x, const Vector< Real > &y) const
const Ptr< Constraint< Real > > con_
Ptr< Vector< Real > > mul_
void project_bnd(Vector< Real > &x, const Vector< Real > &y) const
Provides the interface to apply upper and lower bound constraints.
DykstraProjection(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)
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
void project_Dykstra(Vector< Real > &x, std::ostream &stream=std::cout) const
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Defines the general constraint operator interface.