ROL
ROL_Bundle_U_AS_Def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_BUNDLE_U_AS_DEF_H
11 #define ROL_BUNDLE_U_AS_DEF_H
12 
13 namespace ROL {
14 
15 template<typename Real>
16 Bundle_U_AS<Real>::Bundle_U_AS(const unsigned maxSize,
17  const Real coeff,
18  const Real omega,
19  const unsigned remSize)
20  : Bundle_U<Real>(maxSize,coeff,omega,remSize), isInitialized_(false) {}
21 
22 template<typename Real>
25  if ( !isInitialized_ ) {
26  tG_ = g.clone();
27  yG_ = g.clone();
28  eG_ = g.clone();
29  gx_ = g.clone();
30  ge_ = g.clone();
31  isInitialized_ = true;
32  }
33 }
34 
35 template<typename Real>
36 unsigned Bundle_U_AS<Real>::solveDual(const Real t, const unsigned maxit, const Real tol) {
37  unsigned iter = 0;
38  if (Bundle_U<Real>::size() == 1) {
39  iter = Bundle_U<Real>::solveDual_dim1(t,maxit,tol);
40  }
41  else if (Bundle_U<Real>::size() == 2) {
42  iter = Bundle_U<Real>::solveDual_dim2(t,maxit,tol);
43  }
44  else {
45  iter = solveDual_arbitrary(t,maxit,tol);
46  }
47  return iter;
48 }
49 
50 template<typename Real>
52  const Real zero(0);
53  Real sum(0), err(0), tmp(0), y(0);
54  for (unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
55  // Compute sum of dualVariables_ using Kahan's compensated sum
56  //sum += Bundle_U<Real>::getDualVariable(i);
58  tmp = sum + y;
59  err = (tmp - sum) - y;
60  sum = tmp;
61  }
62  for (unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
65  }
66  nworkingSet_.clear();
67  workingSet_.clear();
68  for (unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
70  workingSet_.insert(i);
71  }
72  else {
73  nworkingSet_.insert(i);
74  }
75  }
76 }
77 
78 template<typename Real>
79 void Bundle_U_AS<Real>::computeLagMult(std::vector<Real> &lam, const Real mu, const std::vector<Real> &g) const {
80  const Real zero(0);
81  unsigned n = workingSet_.size();
82  if ( n > 0 ) {
83  lam.resize(n,zero);
84  typename std::set<unsigned>::iterator it = workingSet_.begin();
85  for (unsigned i = 0; i < n; ++i) {
86  lam[i] = g[*it] - mu; it++;
87  }
88  }
89  else {
90  lam.clear();
91  }
92 }
93 
94 template<typename Real>
95 bool Bundle_U_AS<Real>::isNonnegative(unsigned &ind, const std::vector<Real> &x) const {
96  bool flag = true;
97  unsigned n = workingSet_.size();
98  ind = Bundle_U<Real>::size();
99  if ( n > 0 ) {
100  Real min = ROL_OVERFLOW<Real>();
101  typename std::set<unsigned>::iterator it = workingSet_.begin();
102  for (unsigned i = 0; i < n; ++i) {
103  if ( x[i] < min ) {
104  ind = *it;
105  min = x[i];
106  }
107  it++;
108  }
109  flag = ((min < -ROL_EPSILON<Real>()) ? false : true);
110  }
111  return flag;
112 }
113 
114 template<typename Real>
115 Real Bundle_U_AS<Real>::computeStepSize(unsigned &ind, const std::vector<Real> &x, const std::vector<Real> &p) const {
116  const Real zero(0);
117  Real alpha(1), tmp(0);
118  ind = Bundle_U<Real>::size();
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 ) {
124  alpha = tmp;
125  ind = *it;
126  }
127  }
128  }
129  return std::max(zero,alpha);
130 }
131 
132 template<typename Real>
133 unsigned Bundle_U_AS<Real>::solveEQPsubproblem(std::vector<Real> &s, Real &mu,
134  const std::vector<Real> &g, const Real tol) const {
135  // Build reduced QP information
136  const Real zero(0);
137  unsigned n = nworkingSet_.size(), cnt = 0;
138  mu = zero;
139  s.assign(Bundle_U<Real>::size(),zero);
140  if ( n > 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++;
145  }
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++;
151  }
152  }
153  return cnt;
154 }
155 
156 template<typename Real>
157 void Bundle_U_AS<Real>::applyPreconditioner(std::vector<Real> &Px, const std::vector<Real> &x) const {
158  const Real zero(0);
159  int type = 0;
160  std::vector<Real> tmp(Px.size(),zero);
161  switch (type) {
162  case 0: applyPreconditioner_Identity(tmp,x); break;
163  case 1: applyPreconditioner_Jacobi(tmp,x); break;
164  case 2: applyPreconditioner_SymGS(tmp,x); break;
165  }
166  applyPreconditioner_Identity(Px,tmp);
167 }
168 
169 template<typename Real>
170 void Bundle_U_AS<Real>::applyG(std::vector<Real> &Gx, const std::vector<Real> &x) const {
171  int type = 0;
172  switch (type) {
173  case 0: applyG_Identity(Gx,x); break;
174  case 1: applyG_Jacobi(Gx,x); break;
175  case 2: applyG_SymGS(Gx,x); break;
176  }
177 }
178 
179 template<typename Real>
180 void Bundle_U_AS<Real>::applyPreconditioner_Identity(std::vector<Real> &Px, const std::vector<Real> &x) const {
181  unsigned dim = nworkingSet_.size();
182  Real sum(0), err(0), tmp(0), y(0);
183  for (unsigned i = 0; i < dim; ++i) {
184  // Compute sum of x using Kahan's compensated sum
185  //sum += x[i];
186  y = x[i] - err;
187  tmp = sum + y;
188  err = (tmp - sum) - y;
189  sum = tmp;
190  }
191  sum /= (Real)dim;
192  for (unsigned i = 0; i < dim; ++i) {
193  Px[i] = x[i] - sum;
194  }
195 }
196 
197 template<typename Real>
198 void Bundle_U_AS<Real>::applyG_Identity(std::vector<Real> &Gx, const std::vector<Real> &x) const {
199  Gx.assign(x.begin(),x.end());
200 }
201 
202 template<typename Real>
203 void Bundle_U_AS<Real>::applyPreconditioner_Jacobi(std::vector<Real> &Px, const std::vector<Real> &x) const {
204  const Real zero(0), one(1);
205  unsigned dim = nworkingSet_.size();
206  Real eHe(0), sum(0);
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) {
211  gg[i] = one/std::abs(Bundle_U<Real>::GiGj(*it,*it)); it++;
212  // Compute sum of inv(D)x using Kahan's aggregated sum
213  //sum += x[i]*gg[i];
214  yX = x[i]*gg[i] - errX;
215  tmpX = sum + yX;
216  errX = (tmpX - sum) - yX;
217  sum = tmpX;
218  // Compute sum of inv(D)e using Kahan's aggregated sum
219  //eHe += gg[i];
220  yE = gg[i] - errE;
221  tmpE = eHe + yE;
222  errE = (tmpE - eHe) - yE;
223  eHe = tmpE;
224  }
225  sum /= eHe;
226  for (unsigned i = 0; i < dim; ++i) {
227  Px[i] = (x[i]-sum)*gg[i];
228  }
229 }
230 
231 template<typename Real>
232 void Bundle_U_AS<Real>::applyG_Jacobi(std::vector<Real> &Gx, const std::vector<Real> &x) const {
233  unsigned dim = nworkingSet_.size();
234  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
235  for (unsigned i = 0; i < dim; ++i) {
236  Gx[i] = std::abs(Bundle_U<Real>::GiGj(*it,*it))*x[i]; it++;
237  }
238 }
239 
240 template<typename Real>
241 void Bundle_U_AS<Real>::applyPreconditioner_SymGS(std::vector<Real> &Px, const std::vector<Real> &x) const {
242  int dim = nworkingSet_.size();
243  //unsigned cnt = 0;
244  gx_->zero(); ge_->zero();
245  Real eHx(0), eHe(0), one(1);
246  // Forward substitution
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) {
253  Bundle_U<Real>::addGi(*jt,x1[j],*gx_);
254  Bundle_U<Real>::addGi(*jt,e1[j],*ge_);
255  jt++;
256  }
257  gg[i] = Bundle_U<Real>::GiGj(*it,*it);
258  x1[i] = (x[i] - Bundle_U<Real>::dotGi(*it,*gx_))/gg[i];
259  e1[i] = (one - Bundle_U<Real>::dotGi(*it,*ge_))/gg[i];
260  it++;
261  }
262  // Apply diagonal
263  for (int i = 0; i < dim; ++i) {
264  x1[i] *= gg[i];
265  e1[i] *= gg[i];
266  }
267  // Back substitution
268  std::vector<Real> Hx(dim,0), He(dim,0); it = nworkingSet_.end();
269  for (int i = dim-1; i >= 0; --i) {
270  it--;
271  gx_->zero(); ge_->zero(); jt = nworkingSet_.end();
272  for (int j = dim-1; j >= i+1; --j) {
273  jt--;
274  Bundle_U<Real>::addGi(*jt,Hx[j],*gx_);
275  Bundle_U<Real>::addGi(*jt,He[j],*ge_);
276  }
277  Hx[i] = (x1[i] - Bundle_U<Real>::dotGi(*it,*gx_))/gg[i];
278  He[i] = (e1[i] - Bundle_U<Real>::dotGi(*it,*ge_))/gg[i];
279  // Compute sums
280  eHx += Hx[i];
281  eHe += He[i];
282  }
283  // Accumulate preconditioned vector
284  for (int i = 0; i < dim; ++i) {
285  Px[i] = Hx[i] - (eHx/eHe)*He[i];
286  }
287 }
288 
289 template<typename Real>
290 void Bundle_U_AS<Real>::applyG_SymGS(std::vector<Real> &Gx, const std::vector<Real> &x) const {
291  unsigned dim = nworkingSet_.size();
292  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
293  for (unsigned i = 0; i < dim; ++i) {
294  Gx[i] = std::abs(Bundle_U<Real>::GiGj(*it,*it))*x[i]; it++;
295  }
296 }
297 
298 template<typename Real>
299 void Bundle_U_AS<Real>::computeResidualUpdate(std::vector<Real> &r, std::vector<Real> &g) const {
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);
304  applyG(Gg,g);
305  // Compute multiplier using Kahan's compensated sum
306  for (unsigned i = 0; i < n; ++i) {
307  //y += (r[i] - Gg[i]);
308  yprt = (r[i] - Gg[i]) - yerr;
309  ytmp = y + yprt;
310  yerr = (ytmp - y) - yprt;
311  y = ytmp;
312  }
313  y /= (Real)n;
314  for (unsigned i = 0; i < n; ++i) {
315  r[i] -= y;
316  }
317  applyPreconditioner(g,r);
318 }
319 
320 template<typename Real>
321 void Bundle_U_AS<Real>::applyFullMatrix(std::vector<Real> &Hx, const std::vector<Real> &x) const {
322  const Real one(1);
323  gx_->zero(); eG_->zero();
324  for (unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
325  // Compute Gx using Kahan's compensated sum
326  //gx_->axpy(x[i],Bundle_U<Real>::subgradient(i));
327  yG_->set(Bundle_U<Real>::subgradient(i)); yG_->scale(x[i]); yG_->axpy(-one,*eG_);
328  tG_->set(*gx_); tG_->plus(*yG_);
329  eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
330  gx_->set(*tG_);
331  }
332  for (unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
333  // Compute < g_i, Gx >
334  Hx[i] = Bundle_U<Real>::dotGi(i,*gx_);
335  }
336 }
337 
338 template<typename Real>
339 void Bundle_U_AS<Real>::applyMatrix(std::vector<Real> &Hx, const std::vector<Real> &x) const {
340  const Real one(1);
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) {
345  // Compute Gx using Kahan's compensated sum
346  //gx_->axpy(x[i],Bundle_U<Real>::subgradient(*it));
347  yG_->set(Bundle_U<Real>::subgradient(*it)); yG_->scale(x[i]); yG_->axpy(-one,*eG_);
348  tG_->set(*gx_); tG_->plus(*yG_);
349  eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
350  gx_->set(*tG_);
351  it++;
352  }
353  it = nworkingSet_.begin();
354  for (unsigned i = 0; i < n; ++i) {
355  // Compute < g_i, Gx >
356  Hx[i] = Bundle_U<Real>::dotGi(*it,*gx_); it++;
357  }
358 }
359 
360 template<typename Real>
361 unsigned Bundle_U_AS<Real>::projectedCG(std::vector<Real> &x, Real &mu, const std::vector<Real> &b, const Real tol) const {
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);
365  // Compute residual Hx + g = g with x = 0
366  x.assign(n,0);
367  scale(r,one,b);
368  r0.assign(r.begin(),r.end());
369  // Precondition residual
370  computeResidualUpdate(r,g);
371  Real rg = dot(r,g), rg0(0);
372  // Get search direction
373  scale(d,-one,g);
374  Real alpha(0), kappa(0), beta(0), TOL(1.e-2);
375  Real CGtol = std::min(tol,TOL*rg);
376  unsigned cnt = 0;
377  while (rg > CGtol && cnt < 2*n+1) {
378  applyMatrix(Ad,d);
379  kappa = dot(d,Ad);
380  alpha = rg/kappa;
381  axpy(alpha,d,x);
382  axpy(alpha,Ad,r);
383  axpy(alpha,Ad,r0);
384  computeResidualUpdate(r,g);
385  rg0 = rg;
386  rg = dot(r,g);
387  beta = rg/rg0;
388  scale(d,beta);
389  axpy(-one,g,d);
390  cnt++;
391  }
392  // Compute multiplier for equality constraint using Kahan's compensated sum
393  mu = zero;
394  Real err(0), tmp(0), y(0);
395  for (unsigned i = 0; i < n; ++i) {
396  //mu += r0[i];
397  y = r0[i] - err;
398  tmp = mu + y;
399  err = (tmp - mu) - y;
400  mu = tmp;
401  }
402  mu /= static_cast<Real>(n);
403  // Return iteration count
404  return cnt;
405 }
406 
407 template<typename Real>
408 Real Bundle_U_AS<Real>::dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
409  // Compute dot product of two vectors using Kahan's compensated sum
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) {
413  //val += x[i]*y[i];
414  y0 = x[i]*y[i] - err;
415  tmp = val + y0;
416  err = (tmp - val) - y0;
417  val = tmp;
418  }
419  return val;
420 }
421 
422 template<typename Real>
423 Real Bundle_U_AS<Real>::norm(const std::vector<Real> &x) const {
424  return std::sqrt(dot(x,x));
425 }
426 
427 template<typename Real>
428 void Bundle_U_AS<Real>::axpy(const Real a, const std::vector<Real> &x, std::vector<Real> &y) const {
429  unsigned n = std::min(y.size(),x.size());
430  for (unsigned i = 0; i < n; ++i) {
431  y[i] += a*x[i];
432  }
433 }
434 
435 template<typename Real>
436 void Bundle_U_AS<Real>::scale(std::vector<Real> &x, const Real a) const {
437  for (unsigned i = 0; i < x.size(); ++i) {
438  x[i] *= a;
439  }
440 }
441 
442 template<typename Real>
443 void Bundle_U_AS<Real>::scale(std::vector<Real> &x, const Real a, const std::vector<Real> &y) const {
444  unsigned n = std::min(x.size(),y.size());
445  for (unsigned i = 0; i < n; ++i) {
446  x[i] = a*y[i];
447  }
448 }
449 
450 template<typename Real>
451 unsigned Bundle_U_AS<Real>::solveDual_arbitrary(const Real t, const unsigned maxit, const Real tol) {
452  const Real zero(0), one(1);
453  initializeDualSolver();
454  bool nonneg = false;
455  unsigned ind = 0, i = 0, CGiter = 0;
456  Real snorm(0), alpha(0), mu(0);
457  std::vector<Real> s(Bundle_U<Real>::size(),0), Hs(Bundle_U<Real>::size(),0);
458  std::vector<Real> g(Bundle_U<Real>::size(),0), lam(Bundle_U<Real>::size()+1,0);
459  std::vector<Real> dualVariables(Bundle_U<Real>::size(),0);
460  for (unsigned j = 0; j < Bundle_U<Real>::size(); ++j) {
461  dualVariables[j] = Bundle_U<Real>::getDualVariable(j);
462  }
463  //Real val = Bundle_U<Real>::evaluateObjective(g,dualVariables,t);
464  Bundle_U<Real>::evaluateObjective(g,dualVariables,t);
465  for (i = 0; i < maxit; ++i) {
466  CGiter += solveEQPsubproblem(s,mu,g,tol);
467  snorm = norm(s);
468  if ( snorm < ROL_EPSILON<Real>() ) {
469  computeLagMult(lam,mu,g);
470  nonneg = isNonnegative(ind,lam);
471  if ( nonneg ) {
472  break;
473  }
474  else {
475  alpha = one;
476  if ( ind < Bundle_U<Real>::size() ) {
477  workingSet_.erase(ind);
478  nworkingSet_.insert(ind);
479  }
480  }
481  }
482  else {
483  alpha = computeStepSize(ind,dualVariables,s);
484  if ( alpha > zero ) {
485  axpy(alpha,s,dualVariables);
486  applyFullMatrix(Hs,s);
487  axpy(alpha,Hs,g);
488  }
489  if (ind < Bundle_U<Real>::size()) {
490  workingSet_.insert(ind);
491  nworkingSet_.erase(ind);
492  }
493  }
494  //std::cout << "iter = " << i << " snorm = " << snorm << " alpha = " << alpha << "\n";
495  }
496  //Real crit = computeCriticality(g,dualVariables);
497  //std::cout << "Criticality Measure: " << crit << "\n";
498  //std::cout << "dim = " << Bundle_U<Real>::size() << " iter = " << i << " CGiter = " << CGiter << " CONVERGED!\n";
499  for (unsigned j = 0; j < Bundle_U<Real>::size(); ++j) {
500  Bundle_U<Real>::setDualVariable(j,dualVariables[j]);
501  }
502  return i;
503 }
504 
505 template<typename Real>
506 void Bundle_U_AS<Real>::project(std::vector<Real> &x, const std::vector<Real> &v) const {
507  const Real zero(0), one(1);
508  std::vector<Real> vsort(Bundle_U<Real>::size(),0);
509  vsort.assign(v.begin(),v.end());
510  std::sort(vsort.begin(),vsort.end());
511  Real sum(-1), lam(0);
512  for (unsigned i = Bundle_U<Real>::size()-1; i > 0; i--) {
513  sum += vsort[i];
514  if ( sum >= static_cast<Real>(Bundle_U<Real>::size()-i)*vsort[i-1] ) {
515  lam = sum/static_cast<Real>(Bundle_U<Real>::size()-i);
516  break;
517  }
518  }
519  if (lam == zero) {
520  lam = (sum+vsort[0])/static_cast<Real>(Bundle_U<Real>::size());
521  }
522  for (unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
523  x[i] = std::max(zero,v[i] - lam);
524  }
525 }
526 
527 template<typename Real>
528 Real Bundle_U_AS<Real>::computeCriticality(const std::vector<Real> &g, const std::vector<Real> &sol) const {
529  const Real zero(0), one(1);
530  std::vector<Real> x(Bundle_U<Real>::size(),0), Px(Bundle_U<Real>::size(),0);
531  axpy(one,sol,x);
532  axpy(-one,g,x);
533  project(Px,x);
534  scale(x,zero);
535  axpy(one,sol,x);
536  axpy(-one,Px,x);
537  return norm(x);
538 }
539 
540 } // namespace ROL
541 
542 #endif
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.
Definition: ROL_Vector.hpp:46
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
constexpr auto dim
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