179 const Real
zero(0), one(1);
181 if (ind >= currSize_-dependent_){
183 if (ind < currSize_-1){
185 swapRowsL(ind,currSize_-1);
186 base_[ind] = base_[currSize_-1];
194 L_.reshape(currSize_,currSize_);
195 base_.resize(currSize_);
219 for (
unsigned j=ind+1; j<currSize_-dependent_; ++j){
221 if (std::abs(ai) <= tol*currSize_) {
228 if (std::abs(aj) <= tol*currSize_){
233 else if ( std::abs(ai) > std::abs(aj) ){
236 Real u = sgnb * std::sqrt(one + t*t );
244 Real u = sgna * std::sqrt(one + t*t );
255 L_(j,j) = d; L_(j,ind) =
zero;
257 for (
unsigned h=j+1; h<currSize_; ++h){
258 Real tmp1 = L_(h,ind);
260 L_(h,ind) = Gc*tmp1 + Gs*tmp2;
261 L_(h,j) = Gc*tmp2 - Gs*tmp1;
264 Real tmp1 = z1_[ind];
266 Real tmp3 = z2_[ind];
268 z1_[ind] = Gc*tmp1 + Gs*tmp2;
269 z1_[j] = Gc*tmp2 - Gs*tmp1;
270 z2_[ind] = Gc*tmp3 + Gs*tmp4;
271 z2_[j] = Gc*tmp4 - Gs*tmp3;
275 deltaLh_ = L_(currSize_-dependent_,ind);
277 deltaLj_ = L_(currSize_-1,ind);
281 swapRowsL(ind,currSize_-1);
282 swapRowsL(ind,currSize_-1,
true);
283 L_.reshape(currSize_-1,currSize_-1);
287 unsigned zsize = currSize_-dependent_;
288 for(
unsigned m=ind; m<zsize; m++ ){
296 base_.erase(base_.begin()+ind);
305 Real ghNorm = GiGj(base_[currSize_-dependent_],base_[currSize_-dependent_]);
308 for (
unsigned ii=0; ii<currSize_-dependent_; ++ii){
309 lhnrm += L_(currSize_-dependent_,ii)*L_(currSize_-dependent_,ii);
311 deltaLh_ = std::abs(ghNorm - lhnrm);
313 Real sgn1 = sgn(deltaLh_);
314 Real sgn2 = sgn(dletaLj);
315 Real signum = sgn1 * sgn2;
316 deltaLh_ = std::abs( ghNorm - lhNorm + deltaLh_ * deltaLh_);
319 if( std::sqrt(deltaLh_) > tol*kappa_*std::max(
static_cast<Real
>(1),ghNorm) ){
320 unsigned newind = currSize_-dependent_;
326 for (
unsigned ii=0; ii<newind; ++ii){
327 lh_[ii] = L_(newind,ii);
328 lhz1_ += lh_[ii]*z1_[ii];
329 lhz2_ += lh_[ii]*z2_[ii];
331 deltaLh_ = std::sqrt(deltaLh_);
332 addSubgradToBase(newind,deltaLh_);
336 Real gjNorm = GiGj(base_[currSize_-1],base_[currSize_-1]);
337 ljNorm -= deltaLj_ * deltaLj_;
339 deltaLj_ = std::abs(gjNorm - ljNorm);
341 deltaLj_ = - std::sqrt( deltaLj_ );
343 deltaLj_ = std::sqrt( deltaLj_ );
346 Real gjTgh = GiGj(base_[currSize_-1],base_[currSize_-2]);
348 for (
unsigned ii=0;ii<currSize_;ii++){
349 ljTlh += L_(currSize_-1,ii)*L_(currSize_-2,ii);
351 deltaLj_ = (gjTgh - ljTlh) / deltaLh_;
353 L_(currSize_-1,currSize_-2) = deltaLj_;
359 Real gjNorm = GiGj(base_[currSize_-1],base_[currSize_-1]);
362 for (
unsigned ii=0; ii<currSize_; ++ii) {
363 ljnrm += L_(currSize_-1,ii)*L_(currSize_-1,ii);
365 deltaLj_ = std::abs(gjNorm - ljnrm);
367 deltaLj_ = std::abs( gjNorm - ljNorm + deltaLj_ * deltaLj_);
370 if( std::sqrt(deltaLj_) > tol*kappa_*std::max(
static_cast<Real
>(1),gjNorm) ){
371 unsigned newind = currSize_-1;
377 for (
unsigned ii=0;ii<newind-1;ii++){
378 lj_[ii] = L_(newind,ii);
379 ljz1_ += lj_[ii]*z1_[ii];
380 ljTz2 += lj_[ii]*z2_[ii];
382 deltaLj_ = std::sqrt(deltaLj_);
383 addSubgradToBase(newind,deltaLj_);
385 deltaLh_ = GiGj(base_[currSize_-2],base_[currSize_-1]);
386 for (
unsigned ii=0;ii<currSize_-1;ii++){
387 deltaLh_ -= L_(currSize_-2,ii)*L_(currSize_-1,ii);
389 deltaLh_ /= deltaLj_;
392 deltaLh_ = - std::sqrt( deltaLh_ );
395 deltaLh_ = std::sqrt ( deltaLh_ );
398 L_(currSize_-1,currSize_-2) = deltaLh_;
432 const Real
zero(0), half(0.5), one(1);
433 Real z1z2(0), z1z1(0);
440 rho_ = ROL_INF<Real>();
445 L_(0,0) = std::sqrt(GiGj(0,0));
453 z1_.resize(1); z2_.resize(1);
454 z1_[0] = one/L_(0,0);
459 objval_ = ROL_INF<Real>();
460 minobjval_ = ROL_INF<Real>();
464 for (iter=0;iter<maxit;iter++){
467 switch( dependent_ ){
475 rho_ = ( one + z1z2/t )/z1z1;
476 tempv_ = z1_; tempv_.scale(rho_);
477 tempw1_ = z2_; tempw1_.scale(one/t);
479 solveSystem(currSize_,
'T',L_,tempv_);
489 LA::Matrix<Real> LBprime( LA::Copy,L_,currSize_-1,currSize_-1);
490 lh_.size(currSize_-1);
493 for(
unsigned i=0; i<currSize_-1; ++i){
494 Real tmp = L_(currSize_-1,i);
500 if (std::abs(lhz1_-one) <= tol*kappa_){
503 solveSystem(currSize_-1,
'T',LBprime,tempv_);
504 tempv_.resize(currSize_);
505 tempv_[currSize_-1] = -one;
513 Real tmp = ( one + z1z2 / t - rho_ * z1z1 )/( one - lhz1_ );
514 tempw1_ = z1_; tempw1_.scale(rho_);
515 tempw2_ = z2_; tempw2_.scale(one/t);
517 tempw2_ = lh_; tempw2_.scale(tmp);
519 solveSystem(currSize_-1,
'T',LBprime,tempw1_);
521 tempv_.resize(currSize_);
522 tempv_[currSize_-1] = tmp;
533 LA::Matrix<Real> LBprime( LA::Copy,L_,currSize_-2,currSize_-2 );
534 lj_.size(currSize_-2);
535 lh_.size(currSize_-2);
538 for(
unsigned i=0; i<currSize_-2; ++i){
539 Real tmp1 = L_(currSize_-1,i);
540 Real tmp2 = L_(currSize_-2,i);
541 ljz1_ += tmp1*z1_(i);
542 lhz1_ += tmp2*z1_(i);
546 if(std::abs(ljz1_-one) <= tol*kappa_){
549 solveSystem(currSize_-2,
'T',LBprime,tempv_);
550 tempv_.resize(currSize_);
551 tempv_[currSize_-2] =
zero;
552 tempv_[currSize_-1] = -one;
556 Real mu = ( one - lhz1_ )/( one - ljz1_ );
557 tempw1_ = lj_; tempw1_.scale(-mu);
559 solveSystem(currSize_-2,
'T',LBprime,tempw1_);
561 tempv_.resize(currSize_);
562 tempv_[currSize_-2] = -one;
563 tempv_[currSize_-1] = mu;
572 if (isFeasible(tempv_,tol*currSize_)){
575 for (
unsigned i=0; i<currSize_; ++i){
581 for (
unsigned i=0; i<currSize_; ++i){
589 if ( tempv_[currSize_-1] <
zero )
594 for(
unsigned kk=0; kk<currSize_; ++kk)
603 Real myeps = tol*currSize_;
604 Real step = ROL_INF<Real>();
605 for (
unsigned h=0; h<currSize_; ++h){
612 if (step <=
zero || step == ROL_INF<Real>()){
617 for (
unsigned i=0; i<currSize_; ++i)
624 bool deleted = optimal_;
625 for (
unsigned baseitem=0; baseitem<currSize_; ++baseitem){
631 if( base_[baseitem] == entering_ ){
632 taboo_.push_back(entering_);
638 deleteSubgradFromBase(baseitem,tol);
649 Real newobjval(0), Lin(0), Quad(0);
650 for (
unsigned i=0; i<currSize_; ++i){
653 if (rho_ == ROL_NINF<Real>()){
655 newobjval = -half*Quad;
659 newobjval = half*(rho_ + Lin/t);
665 if( ( entering_ < maxind_ ) && ( objval_ < ROL_INF<Real>() ) ){
666 if( newobjval >= objval_ - std::max( tol*std::abs(objval_), ROL_EPSILON<Real>() ) ){
667 taboo_.push_back(entering_);
675 if ( objval_ + std::max( tol*std::abs(objval_), ROL_EPSILON<Real>() ) <= minobjval_ ){
677 minobjval_ = objval_;
681 if ( (rho_ >= ROL_NINF<Real>()) && (objval_ <= ROL_NINF<Real>()) )
685 Real minro = - std::max( tol*currSize_*std::abs(objval_), ROL_EPSILON<Real>() );
686#if ( ! FIRST_VIOLATED )
687 Real diff = ROL_NINF<Real>(), olddiff = ROL_NINF<Real>();
690 for (
unsigned bundleitem=0; bundleitem<Bundle_U<Real>::size(); ++bundleitem){
692 if ( std::find(taboo_.begin(),taboo_.end(),bundleitem) != taboo_.end() ){
696 if ( std::find(base_.begin(),base_.end(),bundleitem) == base_.end() ){
699 for (
unsigned j=0;j<currSize_;j++){
705 if (rho_ >= ROL_NINF<Real>()){
709 ro = ROL_NINF<Real>();
710 minobjval_ = ROL_INF<Real>();
711 objval_ = ROL_INF<Real>();
715#if ( FIRST_VIOLATED )
716 entering_ = bundleitem;
720 if ( diff > olddiff ){
721 entering_ = bundleitem;
731 if (entering_ < maxind_){
734 base_.push_back(entering_);
736 unsigned zsize = currSize_ - dependent_;
740 for (
unsigned i=0; i<zsize; ++i){
741 lh_[i] = GiGj(entering_,base_[i]);
743 LA::Matrix<Real> LBprime( LA::Copy,L_,zsize,zsize);
744 solveSystem(zsize,
'N',LBprime,lh_);
745 for (
unsigned i=0; i<zsize; ++i){
746 lhz1_ += lh_[i]*z1_[i];
747 lhz2_ += lh_[i]*z2_[i];
750 Real nrm = lh_.dot(lh_);
751 Real delta = GiGj(entering_,entering_) - nrm;
761 L_.reshape(currSize_,currSize_);
762 zsize = currSize_ - dependent_;
763 for (
unsigned i=0; i<zsize-1; ++i){
764 L_(currSize_-1,i) = lh_[i];
767 Real deltaeps = tol*kappa_*std::max(one,lh_.dot(lh_));
768 if ( delta > deltaeps ){
770 unsigned ind = currSize_-1;
771 delta = std::sqrt(delta);
772 addSubgradToBase(ind,delta);
774 else if(delta < -deltaeps){
784 if( objval_ - std::max( tol*std::abs( objval_ ), ROL_EPSILON<Real>() ) > minobjval_ ){