10 #ifndef ROL_BUNDLE_TT_H
11 #define ROL_BUNDLE_TT_H
18 #include "ROL_Ptr.hpp"
19 #include "ROL_LAPACK.hpp"
28 #include "ROL_LinearAlgebra.hpp"
32 #define FIRST_VIOLATED 0
85 Real
GiGj(
const int i,
const int j)
const {
89 Real
sgn(
const Real x)
const {
90 const Real
zero(0), one(1);
91 return ((x <
zero) ? -one :
97 const Real coeff = 0.0,
98 const Real omega = 2.0,
99 const unsigned remSize = 2)
100 :
Bundle<Real>(maxSize,coeff,omega,remSize),
102 maxind_ = std::numeric_limits<int>::max();
104 for(
unsigned i=0; i<
maxSize_; ++i) {
105 Id_(i,i) =
static_cast<Real
>(1);
109 unsigned solveDual(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
127 void swapRowsL(
unsigned ind1,
unsigned ind2,
bool trans=
false) {
128 const Real
zero(0), one(1);
135 for (
unsigned n=ind1+1; n<=ind2; ++n){
137 Id_n(dd,dd) =
zero; Id_n(dd,n) = one;
138 Id_n(n,dd) = one; Id_n(n,n) =
zero;
141 prod.multiply(LA::ETransp::NO_TRANS,LA::ETransp::NO_TRANS,one,Id_n,
L_,
zero);
144 prod.multiply(LA::ETransp::NO_TRANS,LA::ETransp::NO_TRANS,one,
L_,Id_n,
zero);
153 kappa_ =
static_cast<Real
>(1);
156 Real tmpdiagMax = ROL_NINF<Real>();
157 Real tmpdiagMin = ROL_INF<Real>();
159 if(
L_(j,j) > tmpdiagMax ){
160 tmpdiagMax =
L_(j,j);
163 if(
L_(j,j) < tmpdiagMin ){
164 tmpdiagMin =
L_(j,j);
168 kappa_ = tmpdiagMax/tmpdiagMin;
188 unsigned zsize = ind+1;
189 z1_.resize(zsize);
z2_.resize(zsize);
190 z1_[ind] = (
static_cast<Real
>(1) -
lhz1_ ) / delta;
206 const Real
zero(0), one(1);
255 if (std::abs(aj) <= tol*currSize_){
260 else if ( std::abs(ai) > std::abs(aj) ){
263 Real u = sgnb * std::sqrt(one + t*t );
271 Real u = sgna * std::sqrt(one + t*t );
285 Real tmp1 =
L_(h,ind);
287 L_(h,ind) = Gc*tmp1 + Gs*tmp2;
288 L_(h,j) = Gc*tmp2 - Gs*tmp1;
291 Real tmp1 =
z1_[ind];
293 Real tmp3 =
z2_[ind];
295 z1_[ind] = Gc*tmp1 + Gs*tmp2;
296 z1_[j] = Gc*tmp2 - Gs*tmp1;
297 z2_[ind] = Gc*tmp3 + Gs*tmp4;
298 z2_[j] = Gc*tmp4 - Gs*tmp3;
315 for(
unsigned m=ind; m<zsize; m++ ){
338 deltaLh_ = std::abs(ghNorm - lhnrm);
341 Real sgn2 =
sgn(dletaLj);
342 Real signum = sgn1 * sgn2;
346 if( std::sqrt(deltaLh_) > tol*
kappa_*std::max(static_cast<Real>(1),ghNorm) ){
353 for (
unsigned ii=0; ii<newind; ++ii){
354 lh_[ii] =
L_(newind,ii);
358 deltaLh_ = std::sqrt(deltaLh_);
366 deltaLj_ = std::abs(gjNorm -
ljNorm);
368 deltaLj_ = - std::sqrt( deltaLj_ );
370 deltaLj_ = std::sqrt( deltaLj_ );
376 ljTlh +=
L_(currSize_-1,ii)*
L_(currSize_-2,ii);
378 deltaLj_ = (gjTgh - ljTlh) / deltaLh_;
389 for (
unsigned ii=0; ii<
currSize_; ++ii) {
390 ljnrm +=
L_(currSize_-1,ii)*
L_(currSize_-1,ii);
392 deltaLj_ = std::abs(gjNorm - ljnrm);
397 if( std::sqrt(deltaLj_) > tol*
kappa_*std::max(static_cast<Real>(1),gjNorm) ){
398 unsigned newind = currSize_-1;
404 for (
unsigned ii=0;ii<newind-1;ii++){
405 lj_[ii] =
L_(newind,ii);
406 ljz1__ +=
lj_[ii]*
z1_[ii];
409 deltaLj_ = std::sqrt(deltaLj_);
413 for (
unsigned ii=0;ii<currSize_-1;ii++){
414 deltaLh_ -=
L_(currSize_-2,ii)*
L_(currSize_-1,ii);
419 deltaLh_ = - std::sqrt( deltaLh_ );
422 deltaLh_ = std::sqrt ( deltaLh_ );
434 if( L.numRows()!=
size )
435 std::cout <<
"Error: Wrong size matrix!" << std::endl;
436 else if( v.numRows()!=
size )
437 std::cout <<
"Error: Wrong size vector!" << std::endl;
442 lapack_.TRTRS(
'L', tran,
'N', size, 1, L.values(), L.stride(), v.values(), v.stride(), &info );
449 for(
int i=0;i<v.numRows();i++){
457 unsigned solveDual_TT(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
458 const Real
zero(0), half(0.5), one(1);
459 Real z1z2(0), z1z1(0);
466 rho_ = ROL_INF<Real>();
471 L_(0,0) = std::sqrt(
GiGj(0,0));
479 z1_.resize(1);
z2_.resize(1);
480 z1_[0] = one/
L_(0,0);
490 for (iter=0;iter<maxit;iter++){
501 rho_ = ( one + z1z2/t )/z1z1;
515 LA::Matrix<Real> LBprime( LA::Copy,
L_,currSize_-1,currSize_-1);
516 lh_.size(currSize_-1);
519 for(
unsigned i=0; i<currSize_-1; ++i){
520 Real tmp =
L_(currSize_-1,i);
531 tempv_[currSize_-1] = -one;
539 Real tmp = ( one + z1z2 / t -
rho_ * z1z1 )/( one -
lhz1_ );
548 tempv_[currSize_-1] = tmp;
559 LA::Matrix<Real> LBprime( LA::Copy,
L_,currSize_-2,currSize_-2 );
560 lj_.size(currSize_-2);
561 lh_.size(currSize_-2);
564 for(
unsigned i=0; i<currSize_-2; ++i){
565 Real tmp1 =
L_(currSize_-1,i);
566 Real tmp2 =
L_(currSize_-2,i);
578 tempv_[currSize_-1] = -one;
588 tempv_[currSize_-2] = -one;
630 Real step = ROL_INF<Real>();
638 if (step <=
zero || step == ROL_INF<Real>()){
651 for (
unsigned baseitem=0; baseitem<
currSize_; ++baseitem){
675 Real newobjval(0), Lin(0), Quad(0);
679 if (
rho_ == ROL_NINF<Real>()){
681 newobjval = -half*Quad;
685 newobjval = half*(
rho_ + Lin/t);
692 if( newobjval >=
objval_ - std::max( tol*std::abs(
objval_), ROL_EPSILON<Real>() ) ){
701 if (
objval_ + std::max( tol*std::abs(
objval_), ROL_EPSILON<Real>() ) <= minobjval_ ){
707 if ( (
rho_ >= ROL_NINF<Real>()) && (
objval_ <= ROL_NINF<Real>()) )
711 Real minro = - std::max( tol*currSize_*std::abs(
objval_), ROL_EPSILON<Real>() );
712 #if ( ! FIRST_VIOLATED )
713 Real diff = ROL_NINF<Real>(), olddiff = ROL_NINF<Real>();
716 for (
unsigned bundleitem=0; bundleitem<Bundle<Real>::size(); ++bundleitem){
731 if (
rho_ >= ROL_NINF<Real>()){
735 ro = ROL_NINF<Real>();
736 minobjval_ = ROL_INF<Real>();
741 #if ( FIRST_VIOLATED )
746 if ( diff > olddiff ){
766 for (
unsigned i=0; i<zsize; ++i){
769 LA::Matrix<Real> LBprime( LA::Copy,
L_,zsize,zsize);
771 for (
unsigned i=0; i<zsize; ++i){
787 L_.reshape(currSize_,currSize_);
789 for (
unsigned i=0; i<zsize-1; ++i){
790 L_(currSize_-1,i) =
lh_[i];
794 if ( delta > deltaeps ){
796 unsigned ind = currSize_-1;
797 delta = std::sqrt(delta);
800 else if(delta < -deltaeps){
810 if(
objval_ - std::max( tol*std::abs(
objval_ ), ROL_EPSILON<Real>() ) > minobjval_ ){
826 unsigned outermaxit = 20;
827 bool increase =
false, decrease =
false;
829 for (
unsigned it=0; it < outermaxit; ++it ){
835 mytol /=
static_cast<Real
>(10);
839 mytol *=
static_cast<Real
>(10);
842 if ( (mytol > static_cast<Real>(1e-4))
843 || (mytol < static_cast<Real>(1e-16)) ){
846 if ( increase && decrease ) {
Bundle_TT(const unsigned maxSize=10, const Real coeff=0.0, const Real omega=2.0, const unsigned remSize=2)
unsigned solveDual_TT(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void resetDualVariables(void)
bool isFeasible(LA::Vector< Real > &v, const Real &tol)
unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Real sgn(const Real x) const
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
LA::Vector< Real > tempw2_
void setDualVariable(const unsigned i, const Real val)
Contains definitions of custom data types in ROL.
ROL::LAPACK< int, Real > lapack_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
LA::Vector< Real > tempw1_
const Real alpha(const unsigned i) const
LA::Vector< Real > tempv_
void solveSystem(int size, char tran, LA::Matrix< Real > &L, LA::Vector< Real > &v)
const Real getDualVariable(const unsigned i) const
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
unsigned size(void) const
void deleteSubgradFromBase(unsigned ind, Real tol)
Provides the interface for and implements a bundle. The semidefinite quadratic subproblem is solved u...
void addSubgradToBase(unsigned ind, Real delta)
Real GiGj(const int i, const int j) const
std::vector< int > taboo_
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void swapRowsL(unsigned ind1, unsigned ind2, bool trans=false)
Provides the interface for and implements a bundle.