10 #ifndef ROL_BUNDLE_AS_H
11 #define ROL_BUNDLE_AS_H
28 ROL::Ptr<Vector<Real> >
tG_;
29 ROL::Ptr<Vector<Real> >
eG_;
30 ROL::Ptr<Vector<Real> >
yG_;
31 ROL::Ptr<Vector<Real> >
gx_;
32 ROL::Ptr<Vector<Real> >
ge_;
44 const Real coeff = 0.0,
45 const Real omega = 2.0,
46 const unsigned remSize = 2)
61 unsigned solveDual(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
80 Real sum(0), err(0), tmp(0), y(0),
zero(0);
81 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
86 err = (tmp - sum) - y;
89 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
95 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
105 void computeLagMult(std::vector<Real> &lam,
const Real mu,
const std::vector<Real> &g)
const {
110 typename std::set<unsigned>::iterator it =
workingSet_.begin();
111 for (
unsigned i = 0; i < n; ++i) {
112 lam[i] = g[*it] - mu; it++;
124 Real min = ROL_OVERFLOW<Real>();
125 typename std::set<unsigned>::iterator it =
workingSet_.begin();
126 for (
unsigned i = 0; i < n; ++i) {
133 flag = ((min < -ROL_EPSILON<Real>()) ?
false :
true);
138 Real
computeStepSize(
unsigned &ind,
const std::vector<Real> &x,
const std::vector<Real> &p)
const {
140 typename std::set<unsigned>::iterator it;
142 if ( p[*it] < -ROL_EPSILON<Real>() ) {
143 tmp = -x[*it]/p[*it];
144 if (
alpha >= tmp ) {
150 return std::max(zero,
alpha);
154 const std::vector<Real> &g,
const Real tol)
const {
161 std::vector<Real> gk(n,zero);
162 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
163 for (
unsigned i = 0; i < n; ++i) {
164 gk[i] = g[*it]; it++;
166 std::vector<Real> sk(n,zero);
169 for (
unsigned i = 0; i < n; ++i) {
170 s[*it] = sk[i]; it++;
179 std::vector<Real> tmp(Px.size(),
zero);
188 void applyG(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
199 Real sum(0), err(0), tmp(0), y(0);
200 for (
unsigned i = 0; i <
dim; ++i) {
205 err = (tmp - sum) - y;
209 for (
unsigned i = 0; i <
dim; ++i) {
215 Gx.assign(x.begin(),x.end());
220 Real eHe(0), sum(0), one(1),
zero(0);
221 Real errX(0), tmpX(0), yX(0), errE(0), tmpE(0), yE(0);
222 std::vector<Real> gg(dim,zero);
223 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
224 for (
unsigned i = 0; i <
dim; ++i) {
228 yX = x[i]*gg[i] - errX;
230 errX = (tmpX - sum) - yX;
236 errE = (tmpE - eHe) - yE;
240 for (
unsigned i = 0; i <
dim; ++i) {
241 Px[i] = (x[i]-sum)*gg[i];
245 void applyG_Jacobi(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
247 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
248 for (
unsigned i = 0; i <
dim; ++i) {
257 Real eHx(0), eHe(0), one(1);
259 std::vector<Real> x1(dim,0), e1(dim,0),gg(dim,0);
260 typename std::set<unsigned>::iterator it, jt;
262 for (
int i = 0; i <
dim; ++i) {
264 for (
int j = 0; j < i; ++j) {
275 for (
int i = 0; i <
dim; ++i) {
280 std::vector<Real> Hx(dim,0), He(dim,0); it =
nworkingSet_.end();
281 for (
int i = dim-1; i >= 0; --i) {
284 for (
int j = dim-1; j >= i+1; --j) {
296 for (
int i = 0; i <
dim; ++i) {
297 Px[i] = Hx[i] - (eHx/eHe)*He[i];
301 void applyG_SymGS(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
303 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
304 for (
unsigned i = 0; i <
dim; ++i) {
310 unsigned n = g.size();
311 std::vector<Real> Gg(n,0);
312 Real y(0), ytmp(0), yprt(0), yerr(0);
316 for (
unsigned i = 0; i < n; ++i) {
318 yprt = (r[i] - Gg[i]) - yerr;
320 yerr = (ytmp - y) - yprt;
324 for (
unsigned i = 0; i < n; ++i) {
333 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
341 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
347 void applyMatrix(std::vector<Real> &Hx,
const std::vector<Real> &x)
const {
351 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
352 for (
unsigned i = 0; i < n; ++i) {
362 for (
unsigned i = 0; i < n; ++i) {
368 unsigned projectedCG(std::vector<Real> &x, Real &mu,
const std::vector<Real> &b,
const Real tol)
const {
369 Real one(1),
zero(0);
371 std::vector<Real> r(n,0), r0(n,0), g(n,0), d(n,0), Ad(n,0);
375 r0.assign(r.begin(),r.end());
378 Real rg =
dot(r,g), rg0(0);
381 Real
alpha(0), kappa(0), beta(0), TOL(1.e-2);
382 Real CGtol = std::min(tol,TOL*rg);
384 while (rg > CGtol && cnt < 2*n+1) {
401 Real err(0), tmp(0), y(0);
402 for (
unsigned i = 0; i < n; ++i) {
406 err = (tmp - mu) - y;
409 mu /=
static_cast<Real
>(n);
414 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
416 Real val(0), err(0), tmp(0), y0(0);
417 unsigned n = std::min(x.size(),y.size());
418 for (
unsigned i = 0; i < n; ++i) {
420 y0 = x[i]*y[i] - err;
422 err = (tmp - val) - y0;
428 Real
norm(
const std::vector<Real> &x)
const {
429 return std::sqrt(
dot(x,x));
432 void axpy(
const Real a,
const std::vector<Real> &x, std::vector<Real> &y)
const {
433 unsigned n = std::min(y.size(),x.size());
434 for (
unsigned i = 0; i < n; ++i) {
439 void scale(std::vector<Real> &x,
const Real a)
const {
440 for (
unsigned i = 0; i < x.size(); ++i) {
445 void scale(std::vector<Real> &x,
const Real a,
const std::vector<Real> &y)
const {
446 unsigned n = std::min(x.size(),y.size());
447 for (
unsigned i = 0; i < n; ++i) {
455 unsigned ind = 0, i = 0, CGiter = 0;
456 Real snorm(0),
alpha(0), mu(0), one(1),
zero(0);
460 for (
unsigned j = 0; j < Bundle<Real>::size(); ++j) {
465 for (i = 0; i < maxit; ++i) {
468 if ( snorm < ROL_EPSILON<Real>() ) {
484 if (
alpha > zero ) {
499 for (
unsigned j = 0; j < Bundle<Real>::size(); ++j) {
508 void project(std::vector<Real> &x,
const std::vector<Real> &v)
const {
510 vsort.assign(v.begin(),v.end());
511 std::sort(vsort.begin(),vsort.end());
512 Real sum(-1), lam(0),
zero(0), one(1);
523 for (
unsigned i = 0; i < Bundle<Real>::size(); ++i) {
524 x[i] = std::max(
zero,v[i] - lam);
529 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)