44 #ifndef ROL_BUNDLE_AS_H
45 #define ROL_BUNDLE_AS_H
62 ROL::Ptr<Vector<Real> >
tG_;
63 ROL::Ptr<Vector<Real> >
eG_;
64 ROL::Ptr<Vector<Real> >
yG_;
65 ROL::Ptr<Vector<Real> >
gx_;
66 ROL::Ptr<Vector<Real> >
ge_;
78 const Real coeff = 0.0,
79 const Real omega = 2.0,
80 const unsigned remSize = 2)
95 unsigned solveDual(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
114 Real sum(0), err(0), tmp(0), y(0),
zero(0);
115 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
120 err = (tmp - sum) - y;
123 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
129 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
139 void computeLagMult(std::vector<Real> &lam,
const Real mu,
const std::vector<Real> &g)
const {
144 typename std::set<unsigned>::iterator it =
workingSet_.begin();
145 for (
unsigned i = 0; i < n; ++i) {
146 lam[i] = g[*it] - mu; it++;
158 Real min = ROL_OVERFLOW<Real>();
159 typename std::set<unsigned>::iterator it =
workingSet_.begin();
160 for (
unsigned i = 0; i < n; ++i) {
167 flag = ((min < -ROL_EPSILON<Real>()) ?
false :
true);
172 Real
computeStepSize(
unsigned &ind,
const std::vector<Real> &x,
const std::vector<Real> &p)
const {
174 typename std::set<unsigned>::iterator it;
176 if ( p[*it] < -ROL_EPSILON<Real>() ) {
177 tmp = -x[*it]/p[*it];
178 if (
alpha >= tmp ) {
184 return std::max(zero,
alpha);
188 const std::vector<Real> &g,
const Real tol)
const {
195 std::vector<Real> gk(n,zero);
196 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
197 for (
unsigned i = 0; i < n; ++i) {
198 gk[i] = g[*it]; it++;
200 std::vector<Real> sk(n,zero);
203 for (
unsigned i = 0; i < n; ++i) {
204 s[*it] = sk[i]; it++;
213 std::vector<Real> tmp(Px.size(),
zero);
222 void applyG(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
233 Real sum(0), err(0), tmp(0), y(0);
234 for (
unsigned i = 0; i < dim; ++i) {
239 err = (tmp - sum) - y;
243 for (
unsigned i = 0; i < dim; ++i) {
249 Gx.assign(x.begin(),x.end());
254 Real eHe(0), sum(0), one(1),
zero(0);
255 Real errX(0), tmpX(0), yX(0), errE(0), tmpE(0), yE(0);
256 std::vector<Real> gg(dim,zero);
257 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
258 for (
unsigned i = 0; i < dim; ++i) {
262 yX = x[i]*gg[i] - errX;
264 errX = (tmpX - sum) - yX;
270 errE = (tmpE - eHe) - yE;
274 for (
unsigned i = 0; i < dim; ++i) {
275 Px[i] = (x[i]-sum)*gg[i];
279 void applyG_Jacobi(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
281 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
282 for (
unsigned i = 0; i < dim; ++i) {
291 Real eHx(0), eHe(0), one(1);
293 std::vector<Real> x1(dim,0), e1(dim,0),gg(dim,0);
294 typename std::set<unsigned>::iterator it, jt;
296 for (
int i = 0; i < dim; ++i) {
298 for (
int j = 0; j < i; ++j) {
309 for (
int i = 0; i < dim; ++i) {
314 std::vector<Real> Hx(dim,0), He(dim,0); it =
nworkingSet_.end();
315 for (
int i = dim-1; i >= 0; --i) {
318 for (
int j = dim-1; j >= i+1; --j) {
330 for (
int i = 0; i < dim; ++i) {
331 Px[i] = Hx[i] - (eHx/eHe)*He[i];
335 void applyG_SymGS(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
337 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
338 for (
unsigned i = 0; i < dim; ++i) {
344 unsigned n = g.size();
345 std::vector<Real> Gg(n,0);
346 Real y(0), ytmp(0), yprt(0), yerr(0);
350 for (
unsigned i = 0; i < n; ++i) {
352 yprt = (r[i] - Gg[i]) - yerr;
354 yerr = (ytmp - y) - yprt;
358 for (
unsigned i = 0; i < n; ++i) {
367 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
375 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
381 void applyMatrix(std::vector<Real> &Hx,
const std::vector<Real> &x)
const {
385 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
386 for (
unsigned i = 0; i < n; ++i) {
396 for (
unsigned i = 0; i < n; ++i) {
402 unsigned projectedCG(std::vector<Real> &x, Real &mu,
const std::vector<Real> &b,
const Real tol)
const {
403 Real one(1),
zero(0);
405 std::vector<Real> r(n,0), r0(n,0), g(n,0), d(n,0), Ad(n,0);
409 r0.assign(r.begin(),r.end());
412 Real rg =
dot(r,g), rg0(0);
415 Real
alpha(0), kappa(0), beta(0), TOL(1.e-2);
416 Real CGtol = std::min(tol,TOL*rg);
418 while (rg > CGtol && cnt < 2*n+1) {
435 Real err(0), tmp(0), y(0);
436 for (
unsigned i = 0; i < n; ++i) {
440 err = (tmp - mu) - y;
443 mu /=
static_cast<Real
>(n);
448 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
450 Real val(0), err(0), tmp(0), y0(0);
451 unsigned n = std::min(x.size(),y.size());
452 for (
unsigned i = 0; i < n; ++i) {
454 y0 = x[i]*y[i] - err;
456 err = (tmp - val) - y0;
462 Real
norm(
const std::vector<Real> &x)
const {
463 return std::sqrt(
dot(x,x));
466 void axpy(
const Real a,
const std::vector<Real> &x, std::vector<Real> &y)
const {
467 unsigned n = std::min(y.size(),x.size());
468 for (
unsigned i = 0; i < n; ++i) {
473 void scale(std::vector<Real> &x,
const Real a)
const {
474 for (
unsigned i = 0; i < x.size(); ++i) {
479 void scale(std::vector<Real> &x,
const Real a,
const std::vector<Real> &y)
const {
480 unsigned n = std::min(x.size(),y.size());
481 for (
unsigned i = 0; i < n; ++i) {
489 unsigned ind = 0, i = 0, CGiter = 0;
490 Real snorm(0),
alpha(0), mu(0), one(1),
zero(0);
494 for (
unsigned j = 0; j < Bundle<Real>::size(); ++j) {
499 for (i = 0; i < maxit; ++i) {
502 if ( snorm < ROL_EPSILON<Real>() ) {
518 if (
alpha > zero ) {
533 for (
unsigned j = 0; j < Bundle<Real>::size(); ++j) {
542 void project(std::vector<Real> &x,
const std::vector<Real> &v)
const {
544 vsort.assign(v.begin(),v.end());
545 std::sort(vsort.begin(),vsort.end());
546 Real sum(-1), lam(0),
zero(0), one(1);
557 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
558 x[i] = std::max(
zero,v[i] - lam);
563 Real
zero(0), one(1);
std::set< unsigned > nworkingSet_
const Real dotGi(const unsigned i, const Vector< Real > &x) const
ROL::Ptr< Vector< Real > > gx_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void applyPreconditioner(std::vector< Real > &Px, const std::vector< Real > &x) const
void applyG_Jacobi(std::vector< Real > &Gx, const std::vector< Real > &x) const
Real dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void applyPreconditioner_SymGS(std::vector< Real > &Px, const std::vector< Real > &x) const
void applyPreconditioner_Jacobi(std::vector< Real > &Px, const std::vector< Real > &x) const
ROL::Ptr< Vector< Real > > eG_
void axpy(const Real a, const std::vector< Real > &x, std::vector< Real > &y) const
void applyG_Identity(std::vector< Real > &Gx, const std::vector< Real > &x) const
void setDualVariable(const unsigned i, const Real val)
std::set< unsigned > workingSet_
Real norm(const std::vector< Real > &x) const
Defines the linear algebra or vector space interface.
unsigned projectedCG(std::vector< Real > &x, Real &mu, const std::vector< Real > &b, const Real tol) const
void scale(std::vector< Real > &x, const Real a) const
Real computeCriticality(const std::vector< Real > &g, const std::vector< Real > &sol)
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Real computeStepSize(unsigned &ind, const std::vector< Real > &x, const std::vector< Real > &p) const
unsigned solveEQPsubproblem(std::vector< Real > &s, Real &mu, const std::vector< Real > &g, const Real tol) const
void computeLagMult(std::vector< Real > &lam, const Real mu, const std::vector< Real > &g) const
const Real alpha(const unsigned i) const
void initialize(const Vector< Real > &g)
void computeResidualUpdate(std::vector< Real > &r, std::vector< Real > &g) const
ROL::Ptr< Vector< Real > > tG_
void scale(std::vector< Real > &x, const Real a, const std::vector< Real > &y) const
void applyFullMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
bool isNonnegative(unsigned &ind, const std::vector< Real > &x) const
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
const Real getDualVariable(const unsigned i) const
virtual void initialize(const Vector< Real > &g)
void initializeDualSolver(void)
Provides the interface for and implements an active set bundle.
void addGi(const unsigned i, const Real a, Vector< Real > &x) const
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
ROL::Ptr< Vector< Real > > ge_
unsigned size(void) const
void applyPreconditioner_Identity(std::vector< Real > &Px, const std::vector< Real > &x) const
void applyG(std::vector< Real > &Gx, const std::vector< Real > &x) const
const Real GiGj(const unsigned i, const unsigned j) const
void applyMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void applyG_SymGS(std::vector< Real > &Gx, const std::vector< Real > &x) const
unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
ROL::Ptr< Vector< Real > > yG_
void project(std::vector< Real > &x, const std::vector< Real > &v) const
Provides the interface for and implements a bundle.
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Bundle_AS(const unsigned maxSize=10, const Real coeff=0.0, const Real omega=2.0, const unsigned remSize=2)