ROL
ROL_Bundle_TT.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #ifndef ROL_BUNDLE_TT_H
45 #define ROL_BUNDLE_TT_H
46 
47 #include "ROL_Types.hpp"
48 #include "ROL_Vector.hpp"
49 #include "ROL_StdVector.hpp"
50 #include "ROL_Bundle.hpp"
51 
52 #include "Teuchos_RCP.hpp"
53 
54 #include <vector>
55 #include <limits.h>
56 #include <stdint.h>
57 #include <float.h>
58 #include <math.h>
59 #include <algorithm> // TT: std::find
60 
61 #include "Teuchos_SerialDenseMatrix.hpp"
62 #include "Teuchos_SerialDenseVector.hpp"
63 #include "Teuchos_LAPACK.hpp"
64 
65 #define DEBUG_TT 0
66 #define EXACT 1
67 #define TABOO_LIST 1
68 #define FIRST_VIOLATED 0
69 
74 namespace ROL {
75 
76 template<class Real>
77 class Bundle_TT : public Bundle<Real> {
78 
79 private:
80  unsigned maxSize_;
81  Teuchos::LAPACK<int, Real> lapack_; // TT
82 
83 public:
84  Bundle_TT (const unsigned maxSize = 10, const Real coeff = 0.0, const unsigned remSize = 2)
85  : Bundle<Real>(maxSize,coeff,remSize), maxSize_(maxSize) {
86  INF = std::numeric_limits<double>::max();
87  maxind = std::numeric_limits<int>::max();
88  Id.reshape(maxSize_,maxSize_);
89  for(unsigned i=0;i<maxSize_;i++)
90  Id(i,i) = 1.0;
91  }
92 
93  Real GiTGj(const int i, const int j){
94  return (this->subgradient(i)).dot(this->subgradient(j));
95  }
96 
97 /***********************************************************************************************/
98 /****************** DUAL CUTTING PLANE SUBPROBLEM ROUTINES *************************************/
99 /***********************************************************************************************/
100 private:
102  Real INF;
103  int maxind;
104  int entering_; // index of entering item
105  std::vector<int> taboo_; // list of "taboo" items
106  bool optimal_; // flag for optimality of restricted solution
107  unsigned dependent_; // number of lin. dependent items in base
108  Real rho_;
109  unsigned currSize_; // current size of base
110  std::vector<int> base; // base
111  Teuchos::SerialDenseMatrix<int, Real> L;
112  Teuchos::SerialDenseMatrix<int, Real> Id;
113  Teuchos::SerialDenseVector<int, Real> tempv;
114  Teuchos::SerialDenseVector<int, Real> tempw1;
115  Teuchos::SerialDenseVector<int, Real> tempw2;
116  Teuchos::SerialDenseVector<int, Real> lh;
117  Teuchos::SerialDenseVector<int, Real> lj;
118  Teuchos::SerialDenseVector<int, Real> z1;
119  Teuchos::SerialDenseVector<int, Real> z2;
121  int LiMax; // index of max element of diag(L)
122  int LiMin; // index of min element of diag(L)
123  Real kappa; // condition number of matrix L ( >= max|L_ii|/min|L_ii| )
124  Real Quad, Lin; // quadratic and linear terms of objective
125  Real objval; // value of objective
126  Real minobjval; // min value of objective (ever reached)
127  Real deltaLh, deltaLj; // needed in case dependent row becomes independent
128 
129  void swapRowsL(unsigned ind1, unsigned ind2, bool trans=false){
130  if( ind1 > ind2){
131  unsigned tmp = ind1;
132  ind2 = ind1;
133  ind1 = tmp;
134  }
135  unsigned dd = ind1;
136  for (unsigned n=ind1+1; n<=ind2; n++){
137  Teuchos::SerialDenseMatrix<int, Real> Id_n(Teuchos::Copy,Id,currSize_,currSize_);
138  Id_n(dd,dd) = 0; Id_n(dd,n) = 1.0;
139  Id_n(n,dd) = 1.0; Id_n(n,n) = 0;
140  Teuchos::SerialDenseMatrix<int, Real> prod(currSize_,currSize_);
141  if( !trans )
142  prod.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0,Id_n,L,0.0);
143  else
144  prod.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0,L,Id_n,0.0);
145  L = prod;
146  dd++;
147  }
148  }
149 
150  void updateK(){
151  if (currSize_ <= dependent_) { // L is empty
152  kappa = 1.0;
153  }
154  else{
155  Real tmpdiagMax = -INF;
156  Real tmpdiagMin = INF;
157  for (unsigned j=0;j<currSize_-dependent_;j++){
158  if( L(j,j) > tmpdiagMax ){
159  tmpdiagMax = L(j,j);
160  LiMax = j;
161  }
162  if( L(j,j) < tmpdiagMin ){
163  tmpdiagMin = L(j,j);
164  LiMin = j;
165  }
166  }
167  kappa = tmpdiagMax/tmpdiagMin;
168  }
169  }
170 
171  void addSubgradToBase(unsigned ind, Real delta){
172  // update z1, z2, kappa
173  // swap rows if: dependent == 1 and we insert independent row (dependent row is always last)
174  // dependent == 2 and Lj has become independent (Lh still dependent)
175  if(dependent_ && (ind == currSize_-1)){
176  swapRowsL(currSize_-2,currSize_-1);
177  int tmp;
178  tmp = base[currSize_-2];
179  base[currSize_-2] = base[currSize_-1];
180  base[currSize_-1] = tmp;
181  ind--;
182 #if( DEBUG_TT )
183  std::cout << "Swapped last two rows of L " << std::endl;
184 #endif
185  } // end if dependent
186 
187  L(ind,ind) = delta;
188 
189  // update z1 and z2
190  unsigned zsize = ind+1;
191  z1.resize(zsize); z2.resize(zsize);
192  z1[ind] = ( 1.0 - lhTz1 ) / delta;
193  z2[ind] = ( this->alpha(base[ind]) - lhTz2 ) / delta;
194  //z2[zsize-1] = ( this->alpha(entering_) - lhTz2 ) / delta;
195 
196  // update kappa
197  if(delta > L(LiMax,LiMax)){
198  LiMax = ind;
199  kappa = delta/L(LiMin,LiMin);
200  }
201  if(delta < L(LiMin,LiMin)){
202  LiMin = ind;
203  kappa = L(LiMax,LiMax)/delta;
204  }
205  }
206 
207  void deleteSubgradFromBase(unsigned ind, Real tol){
208  // update L, currSize, base, z1, z2, dependent, dualVariables_, kappa
209  if (ind >= currSize_-dependent_){
210 #if( DEBUG_TT )
211  std::cout << "Eliminating dependent item " << base[ind] << std::endl;
212 #endif
213  // if dependent > 0, the last one or two rows of L are lin. dependent
214  if (ind < currSize_-1){ // eliminate currSize_-2 but keep currSize_-1
215 #if( DEBUG_TT )
216  std::cout << "Eliminating Lh but keeping Lj" << std::endl;
217 #endif
218  // swap the last row with the second to last
219  swapRowsL(ind,currSize_-1);
220  base[ind] = base[currSize_-1];
221 #if( ! EXACT )
222  lhNorm = ljNorm; // new last row is lh
223 #endif
224  }
225 
226  dependent_--;
227  currSize_--;
228  L.reshape(currSize_,currSize_); // the row to be eliminated is the last row
229  base.resize(currSize_);
230 
231 #if( DEBUG_TT )
232  std::cout << "New base = " << std::endl;
233  for (unsigned kk=0;kk<base.size();kk++){
234  std::cout << base[kk] << std::endl;
235  }
236  std::cout << "\n";
237  std::cout << "New L " << std::endl;
238  std::cout << L << std::endl;
239  std::cout << "\n";
240 #endif
241  // note: z1, z2, kappa need not be updated
242  return;
243  } // end if dependent item
244 
245 #if( DEBUG_TT )
246  std::cout << "Eliminating INdependent item " << base[baseitem] << std::endl;
247 #endif
248  /* currently L_B is lower trapezoidal
249 
250  | L_1 0 0 |
251  L_B = | l d 0 |
252  | Z v L_2 |
253 
254  Apply Givens rotations to transform it to
255 
256  | L_1 0 0 |
257  | l d 0 |
258  | Z 0 L_2' |
259 
260  then delete row and column to obtain factorization of L_B' with B' = B/{i}
261 
262  L_B' = | L_1 0 |
263  | Z L_2' |
264 
265  */
266  for (unsigned j=ind+1;j<currSize_-dependent_;j++){
267  Real ai = L(j,ind);
268  if (std::abs(ai) <= tol*currSize_) // nothing to do
269  continue;
270  Real aj = L(j,j);
271  Real d, Gc, Gs;
272  // Find Givens
273  // Anderson's implementation
274  if (std::abs(aj) <= tol*currSize_){ // aj is zero
275  Gc = 0.0; d = std::abs(ai); Gs = (ai < 0) - (ai > 0); // Gs = -sgn(ai)
276  }
277  else if ( std::abs(ai) > std::abs(aj) ){
278  Real t = aj/ai;
279  Real sgnb = (ai > 0) - (ai < 0);
280  Real u = sgnb * std::sqrt(1.0 + t*t );
281  Gs = -1.0/u;
282  Gc = -Gs*t;
283  d = ai*u;
284  }
285  else{
286  Real t = ai/aj;
287  Real sgna = (aj > 0) - (aj < 0);
288  Real u = sgna * std::sqrt(1.0 + t*t );
289  Gc = 1.0/u;
290  Gs = -Gc*t;
291  d = aj*u;
292  }
293 
294  // // "Naive" implementation
295  // d = hypot(ai,aj);
296  // Gc = aj/d;
297  // Gs = -ai/d;
298 
299  L(j,j) = d; L(j,ind) = 0.0;
300  // apply Givens to columns i,j of L
301  for (unsigned h=j+1;h<currSize_;h++){
302  Real tmp1 = L(h,ind);
303  Real tmp2 = L(h,j);
304  L(h,ind) = Gc*tmp1 + Gs*tmp2;
305  L(h,j) = Gc*tmp2 - Gs*tmp1;
306  }
307  // apply Givens to z1, z2
308  Real tmp1 = z1[ind];
309  Real tmp2 = z1[j];
310  Real tmp3 = z2[ind];
311  Real tmp4 = z2[j];
312  z1[ind] = Gc*tmp1 + Gs*tmp2;
313  z1[j] = Gc*tmp2 - Gs*tmp1;
314  z2[ind] = Gc*tmp3 + Gs*tmp4;
315  z2[j] = Gc*tmp4 - Gs*tmp3;
316  }// end loop over j
317 
318 #if( DEBUG_TT )
319  std::cout << "After Givens: L,z1,z2 " << std::endl;
320  std::cout << L << std::endl;
321  std::cout << z1 << std::endl;
322  std::cout << z2 << std::endl;
323 #endif
324 
325  if(dependent_){
326  deltaLh = L(currSize_-dependent_,ind); // h = currSize_ - dependent
327  if( dependent_ > 1 ) // j = currSize_ - 1, h = currSize_ - 2
328  deltaLj = L(currSize_-1,ind);
329  }
330 
331  // shift rows and columns of L by exchanging i-th row with next row and i-th column with next column until the row to be deleted is the last, then deleting last row and column
332  swapRowsL(ind,currSize_-1);
333  swapRowsL(ind,currSize_-1,true);
334  L.reshape(currSize_-1,currSize_-1);
335 
336  // delete i-th item from z1 and z2
337  // note: z1 and z2 are of size currSize_-dependent
338  unsigned zsize = currSize_-dependent_;
339  for( unsigned m=ind; m<zsize; m++ ){
340  z1[m] = z1[m+1];
341  z2[m] = z2[m+1];
342  }
343  z1.resize(zsize-1);
344  z2.resize(zsize-1);
345 
346 #if( DEBUG_TT )
347  std::cout << "After elimination: L,z1,z2 " << std::endl;
348  std::cout << L << std::endl;
349  std::cout << z1 << std::endl;
350  std::cout << z2 << std::endl;
351 #endif
352 
353  // update base
354  base.erase(base.begin()+ind);
355 
356 #if( DEBUG_TT )
357  std::cout << "New base = " << std::endl;
358  for (unsigned kk=0;kk<base.size();kk++){
359  std::cout << base[kk] << std::endl;
360  }
361 #endif
362 
363  currSize_--; // update size
364 
365  // update kappa
366  updateK();
367 
368  if(dependent_){
369  // if some previously dependent item have become independent
370 #if( DEBUG_TT )
371  std::cout << "deltaLh = " << deltaLh << std::endl;
372 #endif
373  // recompute deltaLh
374  Real ghNorm = GiTGj(base[currSize_-dependent_],base[currSize_-dependent_]);
375  Real lhnrm = 0.0; // exact lhNorm
376 #if( EXACT)
377  for (unsigned ii=0;ii<currSize_-dependent_;ii++){
378  lhnrm += L(currSize_-dependent_,ii)*L(currSize_-dependent_,ii);
379  }
380  deltaLh = std::abs(ghNorm - lhnrm);
381 #else
382  Real sgn1 = (deltaLh > 0) ? 1 : ((deltaLh < 0) ? -1 : 0);
383  Real sgn2 = (deltaLj > 0) ? 1 : ((deltaLj < 0) ? -1 : 0);
384  Real signum = sgn1 * sgn2; // sgn( deltaLh ) * sgn ( deltaLj );
385  deltaLh = std::abs( ghNorm - lhNorm + deltaLh * deltaLh);
386 #endif
387 
388 #if( DEBUG_TT )
389  std::cout << "ghNorm = " << ghNorm << std::endl;
390  std::cout << "lhNorm (exact) = " << lhnrm << std::endl;
391  std::cout << "lhNorm = " << lhNorm << std::endl;
392  std::cout << "deltaLh = " << std::sqrt(deltaLh) << std::endl;
393  std::cout << "kappa = " << kappa << std::endl;
394 #endif
395 
396  if( std::sqrt(deltaLh) > tol*kappa*std::max(1.0,ghNorm) ){ // originally had just deltaLh (without sqrt)
397 #if( DEBUG_TT )
398  std::cout << "Lh has become lin. INdependent" << std::endl;
399 #endif
400  unsigned newind = currSize_-dependent_;
401  dependent_--;
402  // get the last row of L
403  lh.size(newind); // initialize to zeros;
404  lhTz1 = 0.0;
405  lhTz2 = 0.0;
406  for (unsigned ii=0;ii<newind;ii++){
407  lh[ii] = L(newind,ii);
408  lhTz1 += lh[ii]*z1[ii];
409  lhTz2 += lh[ii]*z2[ii];
410  }
411  deltaLh = std::sqrt(deltaLh);
412  addSubgradToBase(newind,deltaLh);
413 
414  if(dependent_){ // dependent was 2
415 #if( ! EXACT )
416  Real gjNorm = GiTGj(base[currSize_-1],base[currSize_-1]);
417  ljNorm -= deltaLj * deltaLj;
418  lhNorm = gjNorm;
419 
420  deltaLj = std::abs(gjNorm - ljNorm);
421  if ( signum < 0 )
422  deltaLj = - std::sqrt( deltaLj );
423  else
424  deltaLj = std::sqrt( deltaLj );
425 #else
426  // recompute deltaLj
427  Real gjTgh = GiTGj(base[currSize_-1],base[currSize_-2]);
428  Real ljTlh = 0.0;
429  for (unsigned ii=0;ii<currSize_;ii++){
430  ljTlh += L(currSize_-1,ii)*L(currSize_-2,ii);
431  }
432  deltaLj = (gjTgh - ljTlh) / deltaLh;
433 #endif
434 
435  L(currSize_-1,currSize_-2) = deltaLj;
436  }
437 
438 #if( DEBUG_TT )
439  std::cout << "Updated L, z1, z2: " << std::endl;
440  std::cout << L << std::endl;
441  std::cout << z1 << std::endl;
442  std::cout << z2 << std::endl;
443  std::cout << "kappa = " << kappa << std::endl;
444 #endif
445 
446  } // end if deltaLh > 0
447 
448  if (dependent_ > 1){ // deltaLh is 0 but deltaLj is not
449  // recompute deltaLj
450  Real gjNorm = GiTGj(base[currSize_-1],base[currSize_-1]);
451  Real ljnrm = 0.0; // exact ljNorm
452 #if( EXACT )
453  for (unsigned ii=0;ii<currSize_;ii++){
454  ljnrm += L(currSize_-1,ii)*L(currSize_-1,ii);
455  }
456  deltaLj = std::abs(gjNorm - ljnrm);
457 #else
458  deltaLj = std::abs( gjNorm - ljNorm + deltaLj * deltaLj);
459 #endif
460 
461 #if( DEBUG_TT )
462  std::cout << "gjNorm = " << gjNorm << std::endl;
463  std::cout << "ljNorm (exact) = " << ljnrm << std::endl;
464  std::cout << "ljNorm = " << ljNorm << std::endl;
465  std::cout << "deltaLj = " << std::sqrt(deltaLj) << std::endl;
466  std::cout << "kappa = " << kappa << std::endl;
467 #endif
468 
469  if( std::sqrt(deltaLj) > tol*kappa*std::max(1.0,gjNorm) ){ // originally just had deltaLj (without sqrt)
470 #if( DEBUG_TT )
471  std::cout << "Lj has become lin. INdependent" << std::endl;
472 #endif
473  unsigned newind = currSize_-1;
474  dependent_--;
475  // get the last row of L
476  lj.size(newind-1); // initialize to zeros;
477  Real ljTz1 = 0.0;
478  Real ljTz2 = 0.0;
479  for (unsigned ii=0;ii<newind-1;ii++){
480  lj[ii] = L(newind,ii);
481  ljTz1 += lj[ii]*z1[ii];
482  ljTz2 += lj[ii]*z2[ii];
483  }
484  deltaLj = std::sqrt(deltaLj);
485  addSubgradToBase(newind,deltaLj);
486 #if( EXACT )
487  deltaLh = GiTGj(base[currSize_-2],base[currSize_-1]);
488  for (unsigned ii=0;ii<currSize_-1;ii++){
489  deltaLh -= L(currSize_-2,ii)*L(currSize_-1,ii);
490  }
491  deltaLh /= deltaLj;
492 #else
493  if ( signum < 0)
494  deltaLh = - std::sqrt( deltaLh );
495  else
496  deltaLh = std::sqrt ( deltaLh );
497 #endif
498 
499  L(currSize_-1,currSize_-2) = deltaLh;
500 
501  } // end if deltaLj > 0
502  } // end if ( dependent > 1 )
503 
504  } // end if(dependent)
505 
506  }// end deleteSubgradFromBase()
507 
508  Real evaluateObjective(std::vector<Real> &g, const std::vector<Real> &x, const Real t) const {
509  this->gx_->zero(); this->eG_->zero();
510  for (unsigned i = 0; i < this->size(); i++) {
511  // Compute Gx using Kahan's compensated sum
512  this->tG_->set(*this->gx_);
513  this->yG_->set(*this->eG_); this->yG_->axpy(x[i],this->subgradient(i));
514  this->gx_->set(*this->tG_); this->gx_->plus(*this->yG_);
515  this->eG_->set(*this->tG_); this->eG_->axpy(-1.0,*this->gx_); this->eG_->plus(*this->yG_);
516  }
517  Real Hx = 0.0, val = 0.0, err = 0.0, tmp = 0.0, y = 0.0;
518  for (unsigned i = 0; i < this->size(); i++) {
519  // Compute < g_i, Gx >
520  Hx = this->gx_->dot(this->subgradient(i));
521  // Add to the objective function value using Kahan's compensated sum
522  tmp = val;
523  y = x[i]*(0.5*Hx + this->alpha(i)/t) + err;
524  val = tmp + y;
525  err = (tmp - val) + y;
526  // Add gradient component
527  g[i] = Hx + this->alpha(i)/t;
528  }
529  return val;
530  }
531 
532  unsigned solveDual_dim1(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
533  this->setDualVariables(0,1.0);
534  return 0;
535  }
536 
537  unsigned solveDual_dim2(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
538  this->gx_->set(this->subgradient(0));
539  //gx_ = this->subgradient(0).clone();
540  this->gx_->axpy(-1.0,this->subgradient(1));
541  Real diffg = this->gx_->dot(*this->gx_);
542  if ( std::abs(diffg) > ROL_EPSILON ) {
543  Real diffa = (this->alpha(0)-this->alpha(1))/t;
544  Real gdiffg = this->subgradient(1).dot(*this->gx_);
545  this->setDualVariables(0,std::min(1.0,std::max(0.0,-(gdiffg+diffa)/diffg)));
546  this->setDualVariables(1,1.0 - this->getDualVariables(0));
547  }
548  else {
549  if ( std::abs(this->alpha(0)-this->alpha(1)) > ROL_EPSILON ) {
550  if ( this->alpha(0) < this->alpha(1) ) {
551  this->setDualVariables(0,1.0); this->setDualVariables(1,0.0);
552  }
553  else if ( this->alpha(0) > this->alpha(1) ) {
554  this->setDualVariables(0,0.0); this->setDualVariables(1,1.0);
555  }
556  }
557  else {
558  this->setDualVariables(0,0.5); this->setDualVariables(1,0.5);
559  }
560  }
561  return 0;
562  }
563 
564  // TT: solving triangular system for TT algorithm
565  void solveSystem(int size, char tran, Teuchos::SerialDenseMatrix<int,Real> &L, Teuchos::SerialDenseVector<int,Real> &v){
566  int info;
567  if( L.numRows()!=size )
568  std::cout << "Error: Wrong size matrix!" << std::endl;
569  else if( v.numRows()!=size )
570  std::cout << "Error: Wrong size vector!" << std::endl;
571  else if( size==0 )
572  return;
573  else{
574  //std::cout << L.stride() << ' ' << size << std::endl;
575  lapack_.TRTRS( 'L', tran, 'N', size, 1, L.values(), L.stride(), v.values(), v.stride(), &info );
576  }
577  }
578 
579  // TT: check that inequality constraints are satisfied for dual variables
580  bool isFeasible(Teuchos::SerialDenseVector<int,Real> &v, const Real &tol){
581  bool feas = true;
582  for(int i=0;i<v.numRows();i++){
583  if(v[i]<-tol){
584  feas = false;
585  }
586  }
587  return feas;
588  }
589 
590  unsigned solveDual_TT(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
591 #if( DEBUG_TT )
592  std::cout << "Calling solveDual_TT" << std::endl;
593  std::cout << "t = " << t << std::endl;
594  std::cout << "maxit = " << maxit << std::endl;
595  std::cout << "tol = " << tol << std::endl;
596 #endif
597  QPStatus_ = 1; // normal status
598  entering_ = maxind;
599 
600  // cold start
601  optimal_ = true;
602  dependent_ = 0;
603  rho_ = INF; // value of rho = -v
604  currSize_ = 1; // current base size
605  base.clear();
606  base.push_back(0); // initial base
607  L.reshape(1,1);
608  L(0,0) = std::sqrt(GiTGj(0,0));
609  this->resetDualVariables();
610  this->setDualVariables(0,1.0);
611  tempv.resize(1);
612  tempw1.resize(1);
613  tempw2.resize(1);
614  lh.resize(1);
615  lj.resize(1);
616  z1.resize(1); z2.resize(1);
617  z1[0] = 1.0/L(0,0);
618  z2[0] = this->alpha(0)/L(0,0);
619  LiMax = 0; // index of max element of diag(L)
620  LiMin = 0; // index of min element of diag(L)
621  kappa = 1.0; // condition number of matrix L ( >= max|L_ii|/min|L_ii| )
622  objval = INF; // value of objective
623  minobjval = INF; // min value of objective (ever reached)
624 
625  unsigned iter;
626  //-------------------------- MAIN LOOP --------------------------------//
627  for (iter=0;iter<maxit;iter++){
628  //---------------------- INNER LOOP -----------------------//
629  while( !optimal_ ){
630 
631  switch( dependent_ ){
632 
633  case(0): // KT system admits unique solution
634  {
635  /*
636  L = L_B'
637  */
638  z1Tz2 = z1.dot(z2);
639  z1Tz1 = z1.dot(z1);
640  rho_ = ( 1 + z1Tz2/t )/z1Tz1;
641  tempv = z1; tempv.scale(rho_);
642  tempw1 = z2; tempw1.scale(1/t);
643  tempv -= tempw1;
644  solveSystem(currSize_,'T',L,tempv); // tempv contains solution
645  optimal_ = true;
646 #if( DEBUG_TT )
647  std::cout << "In case 0" << std::endl;
648  std::cout << "rho_ = " << rho_ << std::endl;
649  std::cout << "Solution tempv = \n" << tempv << std::endl;
650 #endif
651  break;
652  }
653  case(1):
654  {
655  /*
656  L = | L_B' 0 | \ currSize
657  | l_h^T 0 | /
658  */
659  Teuchos::SerialDenseMatrix<int,Real> LBprime( Teuchos::Copy,L,currSize_-1,currSize_-1);
660  lh.size(currSize_-1); // initialize to zeros;
661  lhTz1 = 0.0;
662  lhTz2 = 0.0;
663  for(unsigned i=0;i<currSize_-1;i++){
664  Real tmp = L(currSize_-1,i);
665  lhTz1 += tmp*z1(i);
666  lhTz2 += tmp*z2(i);
667  lh[i] = tmp;
668  }
669 #if( DEBUG_TT )
670  bool singular = false;
671 #endif
672  // Test for singularity
673  if (std::abs(lhTz1-1.0) <= tol*kappa){
674  // tempv is an infinite direction
675 #if( DEBUG_TT )
676  singular = true;
677 #endif
678  tempv = lh;
679  solveSystem(currSize_-1,'T',LBprime,tempv);
680  tempv.resize(currSize_); // add last entry
681  tempv[currSize_-1] = -1.0;
682  optimal_ = false;
683  }
684  else{
685  // system has (unique) solution
686  rho_ = ( (this->alpha(base[currSize_-1]) - lhTz2)/t ) / ( 1.0 - lhTz1 );
687  z1Tz2 = z1.dot(z2);
688  z1Tz1 = z1.dot(z1);
689  Real tmp = ( 1.0 + z1Tz2 / t - rho_ * z1Tz1 )/( 1.0 - lhTz1 );
690  tempw1 = z1; tempw1.scale(rho_);
691  tempw2 = z2; tempw2.scale(1.0/t);
692  tempw1 -= tempw2;
693  tempw2 = lh; tempw2.scale(tmp);
694  tempw1 -= tempw2;
695  solveSystem(currSize_-1,'T',LBprime,tempw1);
696  tempv = tempw1;
697  tempv.resize(currSize_);
698  tempv[currSize_-1] = tmp;
699  optimal_ = true;
700  }
701 #if( DEBUG_TT )
702  std::cout << "In case 1" << std::endl;
703  if (!singular){
704  std::cout << "rho_ = " << rho_ << std::endl;
705  std::cout << "Solution tempv = \n" << tempv << std::endl;
706  }
707  else
708  std::cout << "Direction tempv = \n" << tempv << std::endl;
709 #endif
710  break;
711  } // case(1)
712  case(2):
713  {
714  /* | L_B' 0 0 | \
715  L = | l_h^T 0 0 | | currSize
716  | l_j^T 0 0 | /
717  */
718  Teuchos::SerialDenseMatrix<int,Real> LBprime( Teuchos::Copy,L,currSize_-2,currSize_-2 );
719  lj.size(currSize_-2); // initialize to zeros;
720  lh.size(currSize_-2); // initialize to zeros;
721  ljTz1 = 0.0;
722  lhTz1 = 0.0;
723  for(unsigned i=0;i<currSize_-2;i++){
724  Real tmp1 = L(currSize_-1,i);
725  Real tmp2 = L(currSize_-2,i);
726  ljTz1 += tmp1*z1(i);
727  lhTz1 += tmp2*z1(i);
728  lj[i] = tmp1;
729  lh[i] = tmp2;
730  }
731  if(std::abs(ljTz1-1.0) <= tol*kappa){
732  // tempv is an infinite direction
733  tempv = lj;
734  solveSystem(currSize_-2,'T',LBprime,tempv);
735  tempv.resize(currSize_); // add two last entries
736  tempv[currSize_-2] = 0.0;
737  tempv[currSize_-1] = -1.0;
738  }
739  else{
740  // tempv is an infinite direction
741  Real mu = ( 1.0 - lhTz1 )/( 1.0 - ljTz1 );
742  tempw1 = lj; tempw1.scale(-mu);
743  tempw1 += lh;
744  solveSystem(currSize_-2,'T',LBprime,tempw1);
745  tempv = tempw1;
746  tempv.resize(currSize_);
747  tempv[currSize_-2] = -1.0;
748  tempv[currSize_-1] = mu;
749  }
750  optimal_ = false;
751 #if( DEBUG_TT )
752  std::cout << "In case 2" << std::endl;
753  std::cout << "Direction tempv = \n" << tempv << std::endl;
754 #endif
755  }// case(2)
756  } // end switch(dependent_)
757 
758  // optimal is true if tempv is a solution, otherwise, tempv is an infinite direction
759  if (optimal_){
760  // check feasibility (inequality constraints)
761  if (isFeasible(tempv,tol*currSize_)){
762 #if( DEBUG_TT )
763  std::cout << "Solution tempv is feasible" << std::endl;
764 #endif
765  // set dual variables to values in tempv
766  this->resetDualVariables();
767  for (unsigned i=0;i<currSize_;i++){
768  this->setDualVariables(base[i],tempv[i]);
769  }
770  }
771  else{
772 #if( DEBUG_TT )
773  std::cout << "Solution tempv is NOT feasible" << std::endl;
774 #endif
775  // w_B = /bar{x}_B - x_B
776  for (unsigned i=0;i<currSize_;i++){
777  tempv[i] -= this->getDualVariables(base[i]);
778  }
779  optimal_ = false;
780  }
781  } // if(optimal)
782  else{ // choose the direction of descent
783  if ( ( entering_ < maxind ) && ( this->getDualVariables(entering_) == 0.0 ) ){
784  if ( tempv[currSize_-1] < 0 ) // w_h < 0
785  tempv.scale(-1.0);
786  }
787  else{ // no i such that dualVariables_[i] == 0
788  Real sp = 0.0;
789  for( unsigned kk=0;kk<currSize_;kk++)
790  sp += tempv[kk]*this->alpha(base[kk]);
791  if ( sp > 0 )
792  tempv.scale(-1.0);
793  }
794  }// end else ( not optimal )
795 
796  if(!optimal_){
797  // take a step in direction tempv (possibly infinite)
798  Real myeps = tol*currSize_;
799  Real step = INF;
800  for (unsigned h=0;h<currSize_;h++){
801  if ( (tempv[h] < -myeps) && (-this->getDualVariables(base[h])/tempv[h] < step) )
802  if ( (this->getDualVariables(base[h]) > myeps) || (this->getDualVariables(base[h]) < -myeps) ) // otherwise, consider it 0
803  step = -this->getDualVariables(base[h])/tempv[h];
804 #if( DEBUG_TT )
805  std::cout << "h = " << h << " tempv[h] = " << tempv[h] << " dualV[base[h]] = " << this->getDualVariables(base[h]) << std::endl;
806 #endif
807  }
808 
809 #if( DEBUG_TT )
810  std::cout << "Taking step of size " << step << std::endl;
811 #endif
812 
813  if (step <= 0 || step == INF){
814 #if( DEBUG_TT )
815  std::cout << "Invalid step!" << std::endl;
816 #endif
817  QPStatus_ = -1; // invalid step
818  return iter;
819  }
820 
821  for (unsigned i=0;i<currSize_;i++)
822  this->setDualVariables(base[i],this->getDualVariables(base[i]) + step * tempv[i]);
823  }// if(!optimal)
824 
825  //------------------------- ITEMS ELIMINATION ---------------------------//
826 
827  // Eliminate items with 0 multipliers from base
828  bool deleted = optimal_;
829  for (unsigned baseitem=0;baseitem<currSize_;baseitem++){
830  if (this->getDualVariables(base[baseitem]) <= tol){
831  deleted = true;
832 
833 #if( TABOO_LIST )
834  // item that just entered shouldn't exit; if it does, mark it as taboo
835  if( base[baseitem] == entering_ ){
836 #if( DEBUG_TT )
837  std::cout << "Blocking " << entering_ << " because it just entered" << std::endl;
838 #endif
839  taboo_.push_back(entering_);
840  entering_ = maxind;
841  }
842 #endif
843 
844  // eliminate item from base;
845  deleteSubgradFromBase(baseitem,tol);
846 
847  } // end if(dualVariables_[baseitem] < tol)
848  } // end loop over baseitem
849 
850  if(!deleted){ // nothing deleted and not optimal
851 #if( DEBUG_TT )
852  std::cout << "Returning because nothing deleted and not optimal" << std::endl;
853 #endif
854  QPStatus_ = -2; // loop
855  return iter;
856  }
857  } // end inner loop
858 
859  Real newobjval; // new objective value
860  Lin = 0.0;
861  for (unsigned i=0;i<currSize_;i++){
862  Lin += this->alpha(base[i])*this->getDualVariables(base[i]);
863  }
864 
865  if (rho_ == -INF){
866  Quad = -Lin/t;
867  newobjval = - Quad/2;
868  }
869  else{
870  Quad = rho_ - Lin/t;
871  newobjval = (rho_ + Lin/t)/2;
872  }
873 
874 #if( DEBUG_TT )
875  std::cout << "New Obj value = " << newobjval << std::endl;
876 #endif
877 
878 #if( TABOO_LIST )
879  // -- test for strict decrease -- //
880  // if item didn't provide decrease, move it to taboo list ...
881  if( ( entering_ < maxind ) && ( objval < INF ) ){
882  if( newobjval >= objval - std::max( tol*std::abs(objval), 1.e-16 ) ){
883 #if( DEBUG_TT )
884  std::cout << "Blocking " << entering_ << " because it didn't provide decrease" << std::endl;
885 #endif
886  taboo_.push_back(entering_);
887  }
888  }
889 #endif
890 
891  objval = newobjval;
892 
893  // if sufficient decrease obtained
894  if ( objval + std::max( tol*std::abs(objval), 1.e-16 ) <= minobjval ){
895  taboo_.clear(); // reset taboo list
896 #if( DEBUG_TT )
897  std::cout << "Taboo list has been reset." << std::endl;
898 #endif
899  minobjval = objval;
900  }
901 
902  //---------------------- OPTIMALITY TEST -------------------------//
903  if ( (rho_ >= -INF) && (objval <= -INF) ) // if current x (dualVariables_) is feasible
904  break;
905 
906  entering_ = maxind;
907  Real minro = - std::max( tol*currSize_*std::abs(objval),1.e-16 );
908 #if ( ! FIRST_VIOLATED )
909  Real diff = -INF, olddiff = -INF;
910 #endif
911 
912  for (unsigned bundleitem=0;bundleitem<this->size();bundleitem++){ // loop over items in bundle
913  //for (int bundleitem=size_-1;bundleitem>-1;bundleitem--){ // loop over items in bundle (in reverse order)
914  if ( std::find(taboo_.begin(),taboo_.end(),bundleitem) != taboo_.end() ){
915 #if( DEBUG_TT )
916  std::cout << "Item " << bundleitem << " is blocked." << std::endl;
917 #endif
918  continue; // if item is taboo, move on
919  }
920 
921  if ( std::find(base.begin(),base.end(),bundleitem) == base.end() ){
922  // base does not contain bundleitem
923 #if( DEBUG_TT )
924  std::cout << "Base does not contain index " << bundleitem << std::endl;
925 #endif
926  Real ro = 0.0;
927  for (unsigned j=0;j<currSize_;j++){
928  ro += this->getDualVariables(base[j]) * GiTGj(bundleitem,base[j]);
929  }
930  ro += this->alpha(bundleitem)/t;
931 
932 #if( DEBUG_TT )
933  std::cout << "RO = " << ro << std::endl;
934 #endif
935 
936  if (rho_ >= -INF){
937  ro = ro - rho_; // note: rho = -v
938  }
939  else{
940  ro = -INF;
941  minobjval = INF;
942  objval = INF;
943  }
944 
945  if (ro < minro){
946 #if ( FIRST_VIOLATED )
947  entering_ = bundleitem;
948  break; // skip going through rest of constraints; alternatively, could look for "most violated"
949 #else
950  diff = minro - ro;
951  if ( diff > olddiff ){
952  entering_ = bundleitem;
953  olddiff = diff;
954  }
955 #endif
956  }
957 
958  } // end if item not in base
959  }// end of loop over items in bundle
960 
961 #if( DEBUG_TT )
962  std::cout << "entering_ = " << entering_ << std::endl;
963 #endif
964 
965  //----------------- INSERTING ITEM ------------------------//
966  if (entering_ < maxind){ // dual constraint is violated
967 #if( DEBUG_TT )
968  std::cout << "Inserting " << entering_ << std::endl;
969 #endif
970  optimal_ = false;
971  this->setDualVariables(entering_,0.0);
972  base.push_back(entering_);
973 #if( DEBUG_TT )
974  std::cout << "Current base = " << std::endl;
975  for (unsigned k=0;k<base.size();k++){
976  std::cout << base[k] << std::endl;
977  }
978  std::cout << "dependent_ = " << dependent_ << std::endl;
979 #endif
980 
981  // construct new row of L_B
982  unsigned zsize = currSize_ - dependent_; // zsize is the size of L_Bprime (current one)
983  lh.size(zsize); // initialize to zeros;
984  lhTz1 = 0.0;
985  lhTz2 = 0.0;
986  for (unsigned i=0;i<zsize;i++){
987  lh[i] = GiTGj(entering_,base[i]);
988  }
989  Teuchos::SerialDenseMatrix<int,Real> LBprime( Teuchos::Copy,L,zsize,zsize);
990  solveSystem(zsize,'N',LBprime,lh); // lh = (L_B^{-1})*(G_B^T*g_h)
991  for (unsigned i=0;i<zsize;i++){
992  lhTz1 += lh[i]*z1[i];
993  lhTz2 += lh[i]*z2[i];
994  }
995 
996  Real nrm = lh.dot(lh);
997  Real delta = GiTGj(entering_,entering_) - nrm; // delta squared
998 #if( DEBUG_TT )
999  std::cout << "GiTGj = " << GiTGj(entering_,entering_) << std::endl;
1000  std::cout << "lh_dot_lh = " << nrm << std::endl;
1001  std::cout << "delta = " << delta << std::endl;
1002 #endif
1003 
1004 #if( ! EXACT )
1005  if(dependent_)
1006  ljNorm = nrm; // adding second dependent
1007  else
1008  lhNorm = nrm; // adding first dependent
1009 #endif
1010 
1011  currSize_++; // update base size
1012 
1013  L.reshape(currSize_,currSize_);
1014  zsize = currSize_ - dependent_; // zsize is the size of L_Bprime (new one)
1015  for (unsigned i=0;i<zsize-1;i++){
1016  L(currSize_-1,i) = lh[i];
1017  }
1018 
1019  Real deltaeps = tol*kappa*std::max(1.0,lh.dot(lh));
1020 #if( DEBUG_TT )
1021  std::cout << "kappa = " << kappa << std::endl;
1022  std::cout << "deltaeps = " << deltaeps << std::endl;
1023 #endif
1024  if ( delta > deltaeps ){ // new row is independent
1025  // add subgradient to the base
1026  unsigned ind = currSize_-1;
1027  delta = std::sqrt(delta);
1028  addSubgradToBase(ind,delta);
1029  }
1030  else if(delta < -deltaeps){
1031 #if( DEBUG_TT )
1032  std::cout << "NEGATIVE delta!" << std::endl;
1033 #endif
1034  dependent_++;
1035  QPStatus_ = 0; // negative delta
1036  return iter;
1037  }
1038  else{ // delta zero
1039  dependent_++;
1040  }
1041 
1042 #if( DEBUG_TT )
1043  std::cout << "Current L = " << std::endl;
1044  std::cout << L << std::endl;
1045  std::cout << "Current z1 = " << std::endl;
1046  std::cout << z1 << std::endl;
1047  std::cout << "Current z2 = " << std::endl;
1048  std::cout << z2 << std::endl;
1049 #endif
1050  } // end if(entering_ < maxind)
1051 
1052  else{ // no dual constraint violated
1053  if( objval - std::max( tol*std::abs( objval ), 1.e-16 ) > minobjval ){ // check if f is as good as minf
1054 #if( DEBUG_TT )
1055  std::cout << "Returning because no dual constraint violated and f cannot reach min value " << minobjval << std::endl;
1056 #endif
1057  QPStatus_ = -3; // loop
1058  return iter;
1059  }
1060  }
1061 
1062  if(optimal_)
1063  break;
1064  } // end main loop
1065 
1066  taboo_.clear();
1067  return iter;
1068  }// end solveDual_TT()
1069 
1070 public:
1071 
1072  unsigned solveDual(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
1073  unsigned iter = 0;
1074  if (this->size() == 1) {
1075  iter = solveDual_dim1(t,maxit,tol);
1076  }
1077  else if (this->size() == 2) {
1078  iter = solveDual_dim2(t,maxit,tol);
1079  }
1080  else {
1081  Real mytol = tol;
1082  unsigned outermaxit = 20;
1083  bool increase = false, decrease = false;
1084  iter = 0;
1085  for ( unsigned it=0; it < outermaxit; it++ ){
1086  iter += solveDual_TT(t,maxit,mytol);
1087  if ( QPStatus_ == 1 )
1088  break;
1089  else if ( QPStatus_ == -2 || QPStatus_ == -3 ){
1090  mytol /= 10.0;
1091  decrease = true;
1092  }
1093  else {
1094  mytol *= 10.0;
1095  increase = true;
1096  }
1097  if ( (mytol > 1.e-4) || (mytol < 1.e-16) ){
1098  break;
1099  }
1100  if ( increase && decrease ){
1101  break;
1102  }
1103  }// end outer loop
1104 #if ( DEBUG_TT )
1105  std::cout << "SolveDual returned after " << iter << " iterations with status " << QPStatus_ << " and solution" << std::endl;
1106  std::vector<Real> sol;
1107  for(unsigned i=0;i<this->size();i++){
1108  sol.push_back(this->getDualVariables(i));
1109  std::cout << "x[" << i << "] = " << sol[i] << std::endl;
1110  }
1111  std::cout << std::endl;
1112  std::vector<Real> g(this->size(),0.0);
1113  Real val = evaluateObjective(g,sol,t);
1114  std::cout << "and objective value = " << val << std::endl;
1115  std::cout << "Checking constraints" << std::endl;
1116  bool success = checkPrimary(g,sol,t);
1117  std::cout << success << std::endl;
1118 #endif
1119  }
1120  return iter;
1121  }
1122 
1123  bool checkPrimary(std::vector<Real> &g, const std::vector<Real> &x, const Real t) const {
1124  bool success = true;
1125  this->gx_->zero(); this->eG_->zero();
1126  for (unsigned i = 0; i < this->size(); i++) {
1127  // Compute Gx using Kahan's compensated sum
1128  this->tG_->set(*this->gx_);
1129  this->yG_->set(*this->eG_); this->yG_->axpy(x[i],this->subgradient(i));
1130  this->gx_->set(*this->tG_); this->gx_->plus(*this->yG_);
1131  this->eG_->set(*this->tG_); this->eG_->axpy(-1.0,*this->gx_); this->eG_->plus(*this->yG_);
1132  }
1133  Real Hx = 0.0, v = 0.0, err = 0.0, tmp = 0.0, y = 0.0;
1134  for (unsigned i = 0; i < this->size(); i++) {
1135  // Compute < g_i, Gx > = - < g_i, d >
1136  Hx = this->gx_->dot(this->subgradient(i));
1137  // Add to the objective function value using Kahan's compensated sum
1138  tmp = v;
1139  y = x[i]*(t*Hx + this->alpha(i)) + err;
1140  v = tmp + y;
1141  err = (tmp - v) + y;
1142  // Add gradient component
1143  g[i] = - t*Hx - this->alpha(i);
1144  }
1145  v *= -1.0;
1146  Real myeps = 1.e-8;
1147  for (unsigned i = 0; i < this->size(); i++) {
1148  if ( g[i] > v + myeps ){
1149  std::cout << "Constraint " << i << " is violated!: g[" << i << "] = " << g[i] << ", v = " << v << std::endl;
1150  success = false;
1151  }
1152  else if ( g[i] < v - myeps ){
1153  std::cout << "Constraint " << i << " is inactive" << std::endl;
1154  }
1155  else{
1156  std::cout << "Constraint " << i << " is active" << std::endl;
1157  }
1158  }
1159  return success;
1160  }
1161 
1162 }; // class Bundle_TT
1163 
1164 } // namespace ROL
1165 
1166 #endif
1167 
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Real GiTGj(const int i, const int j)
unsigned solveDual_TT(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void resetDualVariables(void)
Definition: ROL_Bundle.hpp:265
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)
Teuchos::SerialDenseVector< int, Real > tempw2
Teuchos::RCP< Vector< Real > > yG_
Definition: ROL_Bundle.hpp:114
Teuchos::SerialDenseVector< int, Real > tempv
Contains definitions of custom data types in ROL.
Teuchos::SerialDenseMatrix< int, Real > Id
bool isFeasible(Teuchos::SerialDenseVector< int, Real > &v, const Real &tol)
Teuchos::LAPACK< int, Real > lapack_
Teuchos::SerialDenseVector< int, Real > tempw1
const Vector< Real > & subgradient(const unsigned i) const
Definition: ROL_Bundle.hpp:161
Teuchos::SerialDenseVector< int, Real > z1
Teuchos::SerialDenseVector< int, Real > lj
const Real alpha(const unsigned i) const
Definition: ROL_Bundle.hpp:173
bool checkPrimary(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
Bundle_TT(const unsigned maxSize=10, const Real coeff=0.0, const unsigned remSize=2)
Real getDualVariables(const unsigned i)
Definition: ROL_Bundle.hpp:257
Teuchos::SerialDenseMatrix< int, Real > L
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
Teuchos::RCP< Vector< Real > > eG_
Definition: ROL_Bundle.hpp:113
unsigned size(void) const
Definition: ROL_Bundle.hpp:177
std::vector< int > base
void solveSystem(int size, char tran, Teuchos::SerialDenseMatrix< int, Real > &L, Teuchos::SerialDenseVector< int, Real > &v)
Teuchos::RCP< Vector< Real > > tG_
Definition: ROL_Bundle.hpp:112
void deleteSubgradFromBase(unsigned ind, Real tol)
Provides the interface for and implements a bundle. The semidefinite quadratic subproblem is solved u...
void setDualVariables(const unsigned i, const Real val)
Definition: ROL_Bundle.hpp:261
void addSubgradToBase(unsigned ind, Real delta)
std::vector< int > taboo_
Teuchos::SerialDenseVector< int, Real > lh
void swapRowsL(unsigned ind1, unsigned ind2, bool trans=false)
Teuchos::SerialDenseVector< int, Real > z2
Teuchos::RCP< Vector< Real > > gx_
Definition: ROL_Bundle.hpp:273
Provides the interface for and implments a bundle.
Definition: ROL_Bundle.hpp:62
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118