10 #ifndef ROL_BUNDLE_U_AS_DEF_H
11 #define ROL_BUNDLE_U_AS_DEF_H
15 template<
typename Real>
19 const unsigned remSize)
20 :
Bundle_U<Real>(maxSize,coeff,omega,remSize), isInitialized_(false) {}
22 template<
typename Real>
25 if ( !isInitialized_ ) {
31 isInitialized_ =
true;
35 template<
typename Real>
45 iter = solveDual_arbitrary(t,maxit,tol);
50 template<
typename Real>
53 Real sum(0), err(0), tmp(0), y(0);
54 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
59 err = (tmp - sum) - y;
62 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
68 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
70 workingSet_.insert(i);
73 nworkingSet_.insert(i);
78 template<
typename Real>
81 unsigned n = workingSet_.size();
84 typename std::set<unsigned>::iterator it = workingSet_.begin();
85 for (
unsigned i = 0; i < n; ++i) {
86 lam[i] = g[*it] - mu; it++;
94 template<
typename Real>
97 unsigned n = workingSet_.size();
100 Real min = ROL_OVERFLOW<Real>();
101 typename std::set<unsigned>::iterator it = workingSet_.begin();
102 for (
unsigned i = 0; i < n; ++i) {
109 flag = ((min < -ROL_EPSILON<Real>()) ?
false :
true);
114 template<
typename Real>
117 Real alpha(1), tmp(0);
119 typename std::set<unsigned>::iterator it;
120 for (it = nworkingSet_.begin(); it != nworkingSet_.end(); it++) {
121 if ( p[*it] < -ROL_EPSILON<Real>() ) {
122 tmp = -x[*it]/p[*it];
123 if ( alpha >= tmp ) {
129 return std::max(zero,alpha);
132 template<
typename Real>
134 const std::vector<Real> &g,
const Real tol)
const {
137 unsigned n = nworkingSet_.size(), cnt = 0;
141 std::vector<Real> gk(n,zero);
142 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
143 for (
unsigned i = 0; i < n; ++i) {
144 gk[i] = g[*it]; it++;
146 std::vector<Real> sk(n,zero);
147 cnt = projectedCG(sk,mu,gk,tol);
148 it = nworkingSet_.begin();
149 for (
unsigned i = 0; i < n; ++i) {
150 s[*it] = sk[i]; it++;
156 template<
typename Real>
160 std::vector<Real> tmp(Px.size(),
zero);
162 case 0: applyPreconditioner_Identity(tmp,x);
break;
163 case 1: applyPreconditioner_Jacobi(tmp,x);
break;
164 case 2: applyPreconditioner_SymGS(tmp,x);
break;
166 applyPreconditioner_Identity(Px,tmp);
169 template<
typename Real>
173 case 0: applyG_Identity(Gx,x);
break;
174 case 1: applyG_Jacobi(Gx,x);
break;
175 case 2: applyG_SymGS(Gx,x);
break;
179 template<
typename Real>
181 unsigned dim = nworkingSet_.size();
182 Real sum(0), err(0), tmp(0), y(0);
183 for (
unsigned i = 0; i <
dim; ++i) {
188 err = (tmp - sum) - y;
192 for (
unsigned i = 0; i <
dim; ++i) {
197 template<
typename Real>
199 Gx.assign(x.begin(),x.end());
202 template<
typename Real>
204 const Real
zero(0), one(1);
205 unsigned dim = nworkingSet_.size();
207 Real errX(0), tmpX(0), yX(0), errE(0), tmpE(0), yE(0);
208 std::vector<Real> gg(dim,
zero);
209 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
210 for (
unsigned i = 0; i <
dim; ++i) {
214 yX = x[i]*gg[i] - errX;
216 errX = (tmpX - sum) - yX;
222 errE = (tmpE - eHe) - yE;
226 for (
unsigned i = 0; i <
dim; ++i) {
227 Px[i] = (x[i]-sum)*gg[i];
231 template<
typename Real>
233 unsigned dim = nworkingSet_.size();
234 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
235 for (
unsigned i = 0; i <
dim; ++i) {
240 template<
typename Real>
242 int dim = nworkingSet_.size();
244 gx_->zero(); ge_->zero();
245 Real eHx(0), eHe(0), one(1);
247 std::vector<Real> x1(dim,0), e1(dim,0),gg(dim,0);
248 typename std::set<unsigned>::iterator it, jt;
249 it = nworkingSet_.begin();
250 for (
int i = 0; i <
dim; ++i) {
251 gx_->zero(); ge_->zero(); jt = nworkingSet_.begin();
252 for (
int j = 0; j < i; ++j) {
263 for (
int i = 0; i <
dim; ++i) {
268 std::vector<Real> Hx(dim,0), He(dim,0); it = nworkingSet_.end();
269 for (
int i = dim-1; i >= 0; --i) {
271 gx_->zero(); ge_->zero(); jt = nworkingSet_.end();
272 for (
int j = dim-1; j >= i+1; --j) {
284 for (
int i = 0; i <
dim; ++i) {
285 Px[i] = Hx[i] - (eHx/eHe)*He[i];
289 template<
typename Real>
291 unsigned dim = nworkingSet_.size();
292 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
293 for (
unsigned i = 0; i <
dim; ++i) {
298 template<
typename Real>
300 unsigned n = g.size();
301 std::vector<Real> Gg(n,0);
302 Real y(0), ytmp(0), yprt(0), yerr(0);
303 applyPreconditioner(g,r);
306 for (
unsigned i = 0; i < n; ++i) {
308 yprt = (r[i] - Gg[i]) - yerr;
310 yerr = (ytmp - y) - yprt;
314 for (
unsigned i = 0; i < n; ++i) {
317 applyPreconditioner(g,r);
320 template<
typename Real>
323 gx_->zero(); eG_->zero();
324 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
328 tG_->set(*gx_); tG_->plus(*yG_);
329 eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
332 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
338 template<
typename Real>
341 gx_->zero(); eG_->zero();
342 unsigned n = nworkingSet_.size();
343 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
344 for (
unsigned i = 0; i < n; ++i) {
348 tG_->set(*gx_); tG_->plus(*yG_);
349 eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
353 it = nworkingSet_.begin();
354 for (
unsigned i = 0; i < n; ++i) {
360 template<
typename Real>
362 const Real one(1),
zero(0);
363 unsigned n = nworkingSet_.size();
364 std::vector<Real> r(n,0), r0(n,0), g(n,0), d(n,0), Ad(n,0);
368 r0.assign(r.begin(),r.end());
370 computeResidualUpdate(r,g);
371 Real rg = dot(r,g), rg0(0);
374 Real alpha(0), kappa(0), beta(0), TOL(1.e-2);
375 Real CGtol = std::min(tol,TOL*rg);
377 while (rg > CGtol && cnt < 2*n+1) {
384 computeResidualUpdate(r,g);
394 Real err(0), tmp(0), y(0);
395 for (
unsigned i = 0; i < n; ++i) {
399 err = (tmp - mu) - y;
402 mu /=
static_cast<Real
>(n);
407 template<
typename Real>
410 Real val(0), err(0), tmp(0), y0(0);
411 unsigned n = std::min(x.size(),y.size());
412 for (
unsigned i = 0; i < n; ++i) {
414 y0 = x[i]*y[i] - err;
416 err = (tmp - val) - y0;
422 template<
typename Real>
424 return std::sqrt(dot(x,x));
427 template<
typename Real>
429 unsigned n = std::min(y.size(),x.size());
430 for (
unsigned i = 0; i < n; ++i) {
435 template<
typename Real>
437 for (
unsigned i = 0; i < x.size(); ++i) {
442 template<
typename Real>
444 unsigned n = std::min(x.size(),y.size());
445 for (
unsigned i = 0; i < n; ++i) {
450 template<
typename Real>
452 const Real
zero(0), one(1);
453 initializeDualSolver();
455 unsigned ind = 0, i = 0, CGiter = 0;
456 Real snorm(0), alpha(0), mu(0);
460 for (
unsigned j = 0; j < Bundle_U<Real>::size(); ++j) {
465 for (i = 0; i < maxit; ++i) {
466 CGiter += solveEQPsubproblem(s,mu,g,tol);
468 if ( snorm < ROL_EPSILON<Real>() ) {
469 computeLagMult(lam,mu,g);
470 nonneg = isNonnegative(ind,lam);
477 workingSet_.erase(ind);
478 nworkingSet_.insert(ind);
483 alpha = computeStepSize(ind,dualVariables,s);
484 if ( alpha >
zero ) {
485 axpy(alpha,s,dualVariables);
486 applyFullMatrix(Hs,s);
490 workingSet_.insert(ind);
491 nworkingSet_.erase(ind);
499 for (
unsigned j = 0; j < Bundle_U<Real>::size(); ++j) {
505 template<
typename Real>
507 const Real
zero(0), one(1);
509 vsort.assign(v.begin(),v.end());
510 std::sort(vsort.begin(),vsort.end());
511 Real sum(-1), lam(0);
522 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
523 x[i] = std::max(
zero,v[i] - lam);
527 template<
typename Real>
529 const Real
zero(0), one(1);
bool isNonnegative(unsigned &ind, const std::vector< Real > &x) const
Real norm(const std::vector< Real > &x) const
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void scale(std::vector< Real > &x, const Real a) const
void applyPreconditioner(std::vector< Real > &Px, const std::vector< Real > &x) const
void setDualVariable(const unsigned i, const Real val)
unsigned solveEQPsubproblem(std::vector< Real > &s, Real &mu, const std::vector< Real > &g, const Real tol) const
const Real dotGi(const unsigned i, const Vector< Real > &x) const
void applyPreconditioner_Jacobi(std::vector< Real > &Px, const std::vector< Real > &x) const
const Real GiGj(const unsigned i, const unsigned j) const
Real dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void applyG_Identity(std::vector< Real > &Gx, const std::vector< Real > &x) const
Provides the interface for and implements a bundle.
Defines the linear algebra or vector space interface.
virtual void initialize(const Vector< Real > &g)
void applyPreconditioner_SymGS(std::vector< Real > &Px, const std::vector< Real > &x) const
void applyMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
unsigned size(void) const
Real computeCriticality(const std::vector< Real > &g, const std::vector< Real > &sol) const
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
unsigned projectedCG(std::vector< Real > &x, Real &mu, const std::vector< Real > &b, const Real tol) const
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void computeResidualUpdate(std::vector< Real > &r, std::vector< Real > &g) const
void applyG(std::vector< Real > &Gx, const std::vector< Real > &x) const
void initialize(const Vector< Real > &g)
Real computeStepSize(unsigned &ind, const std::vector< Real > &x, const std::vector< Real > &p) const
void addGi(const unsigned i, const Real a, Vector< Real > &x) const
void applyG_Jacobi(std::vector< Real > &Gx, const std::vector< Real > &x) const
void axpy(const Real a, const std::vector< Real > &x, std::vector< Real > &y) const
void applyPreconditioner_Identity(std::vector< Real > &Px, const std::vector< Real > &x) const
void initializeDualSolver(void)
void applyFullMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
const Real getDualVariable(const unsigned i) const
void computeLagMult(std::vector< Real > &lam, const Real mu, const std::vector< Real > &g) const
void project(std::vector< Real > &x, const std::vector< Real > &v) const
Bundle_U_AS(const unsigned maxSize=10, const Real coeff=0.0, const Real omega=2.0, const unsigned remSize=2)
void applyG_SymGS(std::vector< Real > &Gx, const std::vector< Real > &x) const