45 #ifndef ROL_DOUGLASRACHFORDPROJECTION_DEF_H
46 #define ROL_DOUGLASRACHFORDPROJECTION_DEF_H
50 template<
typename Real>
59 DEFAULT_rtol_ (1e-2*std::sqrt(
ROL_EPSILON<Real>())),
60 DEFAULT_maxit_ (10000),
61 DEFAULT_verbosity_ (0),
62 DEFAULT_alpha1_ (0.5),
63 DEFAULT_gamma_ (10.0),
65 atol_ (DEFAULT_atol_),
66 rtol_ (DEFAULT_rtol_),
67 maxit_ (DEFAULT_maxit_),
68 verbosity_ (DEFAULT_verbosity_),
69 alpha1_ (DEFAULT_alpha1_),
70 alpha2_ (1.0-alpha1_),
71 gamma_ (DEFAULT_gamma_),
80 Real tol(std::sqrt(ROL_EPSILON<Real>()));
85 mul_->setScalar(static_cast<Real>(1));
86 con_->applyAdjointJacobian(*
z_,*
mul_,xprim,tol);
93 template<
typename Real>
102 atol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Absolute Tolerance",
DEFAULT_atol_);
103 rtol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Relative Tolerance",
DEFAULT_rtol_);
104 maxit_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Iteration Limit",
DEFAULT_maxit_);
106 alpha1_ = list.sublist(
"General").sublist(
"Polyhedral Projection").sublist(
"Douglas-Rachford").get(
"Constraint Weight",
DEFAULT_alpha1_);
108 gamma_ = list.sublist(
"General").sublist(
"Polyhedral Projection").sublist(
"Douglas-Rachford").get(
"Penalty Parameter",
DEFAULT_gamma_);
109 t0_ = list.sublist(
"General").sublist(
"Polyhedral Projection").sublist(
"Douglas-Rachford").get(
"Relaxation Parameter",
DEFAULT_t0_);
112 template<
typename Real>
114 if (con_ == nullPtr) {
118 project_DouglasRachford(x, stream);
122 template<
typename Real>
124 return xprim_->dot(x) + b_;
127 template<
typename Real>
129 Real tol(std::sqrt(ROL_EPSILON<Real>()));
131 con_->value(r,y,tol);
134 template<
typename Real>
140 template<
typename Real>
143 Real rhs = residual_1d(y);
144 Real lam = -rhs/cdot_;
149 Real tol = std::sqrt(ROL_EPSILON<Real>());
150 residual_nd(*res_,y);
151 con_->solveAugmentedSystem(x,*mul_,*z_,*res_,y,tol);
152 x.
scale(static_cast<Real>(-1));
157 template<
typename Real>
159 const Real one(1), two(2), xnorm(x.
norm()), ctol(std::min(atol_,rtol_*xnorm));
161 p_->zero(); q_->zero(); y_->set(x);
162 std::ios_base::fmtflags streamFlags(stream.flags());
163 if (verbosity_ > 2) {
164 stream << std::scientific << std::setprecision(6);
166 stream <<
" Polyhedral Projection using Douglas Rachford Splitting" << std::endl;
168 stream << std::setw(6) << std::left <<
"iter";
169 stream << std::setw(15) << std::left <<
"error";
170 stream << std::setw(15) << std::left <<
"tol";
173 for (
int cnt=0; cnt < maxit_; ++cnt) {
176 tmp_->axpy(alpha1_*gamma_,x);
177 tmp_->scale(one/(alpha1_*gamma_+one));
178 project_con(*p_,*tmp_);
182 tmp_->axpy(-one,*y_);
183 tmp_->axpy(alpha2_*gamma_,x);
184 tmp_->scale(one/(alpha2_*gamma_+one));
185 project_bnd(*q_,*tmp_);
188 tmp_->axpy(-one,*p_);
189 rnorm = tmp_->norm();
190 if (verbosity_ > 2) {
192 stream << std::setw(6) << std::left << cnt;
193 stream << std::setw(15) << std::left << rnorm;
194 stream << std::setw(15) << std::left << ctol;
197 if (rnorm <= ctol)
break;
200 if (verbosity_ > 2) stream << std::endl;
204 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
205 stream << rnorm <<
" rtol = " << ctol << std::endl;
207 stream.flags(streamFlags);
void project_DouglasRachford(Vector< Real > &x, std::ostream &stream=std::cout) const
Ptr< Vector< Real > > xprim_
void project_bnd(Vector< Real > &x, const Vector< Real > &y) const
virtual void scale(const Real alpha)=0
Compute where .
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_
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
DouglasRachfordProjection(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)
Real residual_1d(const Vector< Real > &x) const
Defines the linear algebra or vector space interface.
Ptr< Vector< Real > > tmp_
const Ptr< Constraint< Real > > con_
Ptr< Vector< Real > > mul_
void residual_nd(Vector< Real > &r, const Vector< Real > &y) const
void project_con(Vector< Real > &x, const Vector< Real > &y) const
Provides the interface to apply upper and lower bound constraints.
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Defines the general constraint operator interface.