50 #include "Teuchos_RCP.hpp" 81 void remove(
const std::vector<unsigned> &ind) {
83 for (
unsigned j = ind.back()+1; j <
size_; j++) {
93 for (
unsigned i = ind.size()-1; i > 0; --i) {
94 for (
unsigned j = ind[i-1]+1; j <
size_; j++) {
114 Teuchos::RCP<Vector<Real> >
tG_;
115 Teuchos::RCP<Vector<Real> >
eG_;
116 Teuchos::RCP<Vector<Real> >
yG_;
123 Bundle(
const unsigned maxSize = 10,
const Real coeff = 0.0,
const unsigned remSize = 2)
139 Real zero(0), one(1);
140 for (
unsigned i = 0; i <
maxSize_; i++) {
170 Real
alpha = le, two(2);
171 if (
coeff_ > ROL_EPSILON<Real>() ) {
177 const Real
alpha(
const unsigned i)
const {
186 Real zero(0), one(1);
187 aggSubGrad.
zero(); aggLinErr = zero; aggDistMeas = zero;
eG_->zero();
188 Real eLE(0), eDM(0), yLE(0), yDM(0), tLE(0), tDM(0);
189 for (
unsigned i = 0; i <
size_; i++) {
191 tG_->set(aggSubGrad);
198 aggLinErr = tLE + yLE;
199 eLE = (tLE - aggLinErr) + yLE;
203 aggDistMeas = tDM + yDM;
204 eDM = (tDM - aggDistMeas) + yDM;
211 unsigned loc =
size_, cnt = 0;
212 std::vector<unsigned> ind(
remSize_,0);
213 for (
unsigned i =
size_; i > 0; --i) {
219 for (
unsigned i = 0; i <
size_; i++) {
220 if ( i < loc || i > loc ) {
235 void update(
const bool flag,
const Real linErr,
const Real distMeas,
240 for (
unsigned i = 0; i <
size_; i++) {
280 Teuchos::RCP<Vector<Real> >
gx_;
281 Teuchos::RCP<Vector<Real> >
ge_;
300 Real sum(0), err(0), tmp(0), y(0), zero(0);
301 for (
unsigned i = 0; i <
size_; i++) {
306 err = (tmp - sum) + y;
308 for (
unsigned i = 0; i <
size_; i++) {
313 for (
unsigned i = 0; i <
size_; i++) {
324 Real one(1), half(0.5);
326 for (
unsigned i = 0; i <
size_; i++) {
333 Real Hx(0), val(0), err(0), tmp(0), y(0);
334 for (
unsigned i = 0; i <
size_; i++) {
339 y = x[i]*(half*Hx +
alpha(i)/t) + err;
341 err = (tmp - val) + y;
343 g[i] = Hx +
alpha(i)/t;
351 for (
unsigned i = 0; i <
size_; i++) {
358 for (
unsigned i = 0; i <
size_; i++) {
364 void applyMatrix(std::vector<Real> &Hx,
const std::vector<Real> &x)
const {
368 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
369 for (
unsigned i = 0; i < n; i++) {
378 for (
unsigned i = 0; i < n; i++) {
384 void computeLagMult(std::vector<Real> &lam,
const Real mu,
const std::vector<Real> &g)
const {
389 typename std::set<unsigned>::iterator it =
workingSet_.begin();
390 for (
unsigned i = 0; i < n; i++) {
391 lam[i] = g[*it] - mu; it++;
403 Real min = ROL_OVERFLOW<Real>();
404 typename std::set<unsigned>::iterator it =
workingSet_.begin();
405 for (
unsigned i = 0; i < n; i++) {
412 flag = ((min < -ROL_EPSILON<Real>()) ? false :
true);
417 Real
computeAlpha(
unsigned &ind,
const std::vector<Real> &x,
const std::vector<Real> &p)
const {
419 typename std::set<unsigned>::iterator it;
421 if ( p[*it] < -ROL_EPSILON<Real>() ) {
422 tmp = -x[*it]/p[*it];
423 if (
alpha >= tmp ) {
429 return std::max(zero,
alpha);
433 const std::vector<Real> &g,
const Real tol)
const {
438 s.assign(
size_,zero);
440 std::vector<Real> gk(n,zero);
441 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
442 for (
unsigned i = 0; i < n; i++) {
443 gk[i] = g[*it]; it++;
445 std::vector<Real> sk(n,zero);
448 for (
unsigned i = 0; i < n; i++) {
449 s[*it] = sk[i]; it++;
458 std::vector<Real> tmp(Px.size(),zero);
467 void applyG(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
478 Real sum(0), err(0), tmp(0), y(0);
479 for (
unsigned i = 0; i < dim; i++) {
484 err = (tmp - sum) + y;
487 for (
unsigned i = 0; i < dim; i++) {
493 Gx.assign(x.begin(),x.end());
498 Real eHe(0), sum(0), one(1), zero(0);
499 Real errX(0), tmpX(0), yX(0), errE(0), tmpE(0), yE(0);
500 std::vector<Real> gg(dim,zero);
501 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
502 for (
unsigned i = 0; i < dim; i++) {
506 yX = x[i]*gg[i] + errX;
508 errX = (tmpX - sum) + yX;
513 errE = (tmpE - eHe) + yE;
516 for (
unsigned i = 0; i < dim; i++) {
517 Px[i] = (x[i]-sum)*gg[i];
521 void applyG_Jacobi(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
523 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
524 for (
unsigned i = 0; i < dim; i++) {
533 Real eHx(0), eHe(0), one(1);
535 std::vector<Real> x1(dim,0), e1(dim,0),gg(dim,0);
536 typename std::set<unsigned>::iterator it, jt;
538 for (
int i = 0; i < dim; i++) {
540 for (
int j = 0; j < i; j++) {
551 for (
int i = 0; i < dim; i++) {
556 std::vector<Real> Hx(dim,0), He(dim,0); it =
nworkingSet_.end();
557 for (
int i = dim-1; i >= 0; --i) {
560 for (
int j = dim-1; j >= i+1; --j) {
572 for (
int i = 0; i < dim; i++) {
573 Px[i] = Hx[i] - (eHx/eHe)*He[i];
577 void applyG_SymGS(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
579 typename std::set<unsigned>::iterator it =
nworkingSet_.begin();
580 for (
unsigned i = 0; i < dim; i++) {
586 unsigned n = g.size();
587 std::vector<Real> Gg(n,0);
588 Real y(0), ytmp(0), yprt(0), yerr(0);
592 for (
unsigned i = 0; i < n; i++) {
594 yprt = (r[i] - Gg[i]) + yerr;
596 yerr = (ytmp - y) + yprt;
599 for (
unsigned i = 0; i < n; i++) {
605 unsigned projectedCG(std::vector<Real> &x, Real &mu,
const std::vector<Real> &b,
const Real tol)
const {
606 Real one(1), zero(0);
608 std::vector<Real> r(n,0), r0(n,0), g(n,0), d(n,0), Ad(n,0);
612 r0.assign(r.begin(),r.end());
615 Real rg =
dot(r,g), rg0(0);
618 Real
alpha(0), kappa(0), beta(0), TOL(1.e-2);
619 Real CGtol = std::min(tol,TOL*rg);
621 while (rg > CGtol && cnt < 2*n+1) {
638 Real err(0), tmp(0), y(0);
639 for (
unsigned i = 0; i < n; i++) {
643 err = (tmp - mu) + y;
650 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
652 Real val(0), err(0), tmp(0), y0(0);
653 unsigned n = std::min(x.size(),y.size());
654 for (
unsigned i = 0; i < n; i++) {
656 y0 = x[i]*y[i] + err;
658 err = (tmp - val) + y0;
663 Real
norm(
const std::vector<Real> &x)
const {
664 return std::sqrt(
dot(x,x));
667 void axpy(
const Real a,
const std::vector<Real> &x, std::vector<Real> &y)
const {
668 unsigned n = std::min(y.size(),x.size());
669 for (
unsigned i = 0; i < n; i++) {
674 void scale(std::vector<Real> &x,
const Real a)
const {
675 for (
unsigned i = 0; i < x.size(); i++) {
680 void scale(std::vector<Real> &x,
const Real a,
const std::vector<Real> &y)
const {
681 unsigned n = std::min(x.size(),y.size());
682 for (
unsigned i = 0; i < n; i++) {
687 unsigned solveDual_dim1(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
693 unsigned solveDual_dim2(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
694 Real diffg =
gx_->dot(*
gx_), zero(0), one(1), half(0.5);
696 if ( std::abs(diffg) > ROL_EPSILON<Real>() ) {
699 dualVariables_[0] = std::min(one,std::max(zero,-(gdiffg+diffa)/diffg));
703 if ( std::abs(
alpha(0)-
alpha(1)) > ROL_EPSILON<Real>() ) {
722 unsigned ind = 0, i = 0, CGiter = 0;
723 Real snorm(0),
alpha(0), mu(0), one(1), zero(0);
727 for (i = 0; i < maxit; i++) {
730 if ( snorm < ROL_EPSILON<Real>() ) {
746 if (
alpha > zero ) {
765 virtual unsigned solveDual(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
770 else if (
size_ == 2) {
783 void project(std::vector<Real> &x,
const std::vector<Real> &v)
const {
784 std::vector<Real> vsort(
size_,0);
785 vsort.assign(v.begin(),v.end());
786 std::sort(vsort.begin(),vsort.end());
787 Real sum(-1), lam(0), zero(0), one(1);
788 for (
int i =
size_-1; i > 0; i--) {
790 if ( sum >= ((Real)(
size_-i))*vsort[i-1] ) {
791 lam = sum/(Real)(
size_-i);
796 lam = (sum+vsort[0])/(Real)
size_;
798 for (
int i = 0; i <
size_; i++) {
799 x[i] = std::max(zero,v[i] - lam);
804 Real zero(0), one(1);
void initializeDualSolver(void)
void applyG_Jacobi(std::vector< Real > &Gx, const std::vector< Real > &x) const
void resetDualVariables(void)
void applyG(std::vector< Real > &Gx, const std::vector< Real > &x) const
unsigned projectedCG(std::vector< Real > &x, Real &mu, const std::vector< Real > &b, const Real tol) const
void computeResidualUpdate(std::vector< Real > &r, std::vector< Real > &g) const
void applyFullMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
std::set< unsigned > workingSet_
virtual void plus(const Vector &x)=0
Compute , where .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void update(const bool flag, const Real linErr, const Real distMeas, const Vector< Real > &g, const Vector< Real > &s)
Teuchos::RCP< Vector< Real > > yG_
const Real alpha(const unsigned i) const
Contains definitions of custom data types in ROL.
Real computeCriticality(const std::vector< Real > &g)
const Vector< Real > & subgradient(const unsigned i) const
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void reset(const Vector< Real > &g, const Real le, const Real dm)
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
void applyPreconditioner_SymGS(std::vector< Real > &Px, const std::vector< Real > &x) const
std::set< unsigned > nworkingSet_
std::vector< Real > dualVariables_
const Real distanceMeasure(const unsigned i) const
Teuchos::RCP< Vector< Real > > ge_
void applyG_SymGS(std::vector< Real > &Gx, const std::vector< Real > &x) const
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void scale(std::vector< Real > &x, const Real a, const std::vector< Real > &y) const
std::vector< Teuchos::RCP< Vector< Real > > > subgradients_
void computeLagMult(std::vector< Real > &lam, const Real mu, const std::vector< Real > &g) const
Real norm(const std::vector< Real > &x) const
Real computeAlpha(unsigned &ind, const std::vector< Real > &x, const std::vector< Real > &p) const
Real getDualVariables(const unsigned i)
std::vector< Real > distanceMeasures_
void applyPreconditioner_Jacobi(std::vector< Real > &Px, const std::vector< Real > &x) const
unsigned size(void) const
const Real computeAlpha(const Real dm, const Real le) const
void applyPreconditioner(std::vector< Real > &Px, const std::vector< Real > &x) const
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void axpy(const Real a, const std::vector< Real > &x, std::vector< Real > &y) const
Teuchos::RCP< Vector< Real > > eG_
Teuchos::RCP< Vector< Real > > tG_
virtual void set(const Vector &x)
Set where .
void applyPreconditioner_Identity(std::vector< Real > &Px, const std::vector< Real > &x) const
void applyG_Identity(std::vector< Real > &Gx, const std::vector< Real > &x) const
void applyMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
unsigned solveEQPsubproblem(std::vector< Real > &s, Real &mu, const std::vector< Real > &g, const Real tol) const
void setDualVariables(const unsigned i, const Real val)
virtual unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
std::vector< Real > linearizationErrors_
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void initialize(const Vector< Real > &g)
void scale(std::vector< Real > &x, const Real a) const
void add(const Vector< Real > &g, const Real le, const Real dm)
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
Teuchos::RCP< Vector< Real > > gx_
Bundle(const unsigned maxSize=10, const Real coeff=0.0, const unsigned remSize=2)
bool isNonnegative(unsigned &ind, const std::vector< Real > &x) const
const Real linearizationError(const unsigned i) const
Provides the interface for and implments a bundle.
void project(std::vector< Real > &x, const std::vector< Real > &v) const
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void aggregate(Vector< Real > &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const