44 #ifndef ROL_BUNDLE_TT_H 45 #define ROL_BUNDLE_TT_H 52 #include "Teuchos_RCP.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_LAPACK.hpp" 68 #define FIRST_VIOLATED 0 84 Bundle_TT (
const unsigned maxSize = 10,
const Real coeff = 0.0,
const unsigned remSize = 2)
86 INF = std::numeric_limits<double>::max();
87 maxind = std::numeric_limits<int>::max();
93 Real
GiTGj(
const int i,
const int j){
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;
129 void swapRowsL(
unsigned ind1,
unsigned ind2,
bool trans=
false){
136 for (
unsigned n=ind1+1; n<=ind2; n++){
138 Id_n(dd,dd) = 0; Id_n(dd,n) = 1.0;
139 Id_n(n,dd) = 1.0; Id_n(n,n) = 0;
142 prod.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0,Id_n,
L,0.0);
144 prod.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0,
L,Id_n,0.0);
155 Real tmpdiagMax = -
INF;
156 Real tmpdiagMin =
INF;
158 if(
L(j,j) > tmpdiagMax ){
162 if(
L(j,j) < tmpdiagMin ){
167 kappa = tmpdiagMax/tmpdiagMin;
183 std::cout <<
"Swapped last two rows of L " << std::endl;
190 unsigned zsize = ind+1;
191 z1.resize(zsize);
z2.resize(zsize);
192 z1[ind] = ( 1.0 -
lhTz1 ) / delta;
211 std::cout <<
"Eliminating dependent item " <<
base[ind] << std::endl;
216 std::cout <<
"Eliminating Lh but keeping Lj" << std::endl;
232 std::cout <<
"New base = " << std::endl;
233 for (
unsigned kk=0;kk<
base.size();kk++){
234 std::cout <<
base[kk] << std::endl;
237 std::cout <<
"New L " << std::endl;
238 std::cout <<
L << std::endl;
246 std::cout <<
"Eliminating INdependent item " <<
base[baseitem] << std::endl;
275 Gc = 0.0; d = std::abs(ai); Gs = (ai < 0) - (ai > 0);
277 else if ( std::abs(ai) > std::abs(aj) ){
279 Real sgnb = (ai > 0) - (ai < 0);
280 Real u = sgnb * std::sqrt(1.0 + t*t );
287 Real sgna = (aj > 0) - (aj < 0);
288 Real u = sgna * std::sqrt(1.0 + t*t );
299 L(j,j) = d;
L(j,ind) = 0.0;
302 Real tmp1 =
L(h,ind);
304 L(h,ind) = Gc*tmp1 + Gs*tmp2;
305 L(h,j) = Gc*tmp2 - Gs*tmp1;
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;
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;
339 for(
unsigned m=ind; m<zsize; m++ ){
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;
357 std::cout <<
"New base = " << std::endl;
358 for (
unsigned kk=0;kk<
base.size();kk++){
359 std::cout <<
base[kk] << std::endl;
371 std::cout <<
"deltaLh = " <<
deltaLh << std::endl;
380 deltaLh = std::abs(ghNorm - lhnrm);
384 Real signum = sgn1 * sgn2;
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;
396 if( std::sqrt(
deltaLh) > tol*
kappa*std::max(1.0,ghNorm) ){
398 std::cout <<
"Lh has become lin. INdependent" << std::endl;
406 for (
unsigned ii=0;ii<newind;ii++){
407 lh[ii] =
L(newind,ii);
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;
456 deltaLj = std::abs(gjNorm - ljnrm);
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;
469 if( std::sqrt(
deltaLj) > tol*
kappa*std::max(1.0,gjNorm) ){
471 std::cout <<
"Lj has become lin. INdependent" << std::endl;
479 for (
unsigned ii=0;ii<newind-1;ii++){
480 lj[ii] =
L(newind,ii);
482 ljTz2 +=
lj[ii]*
z2[ii];
509 this->
gx_->zero(); this->
eG_->zero();
510 for (
unsigned i = 0; i < this->
size(); i++) {
512 this->
tG_->set(*this->
gx_);
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_);
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++) {
523 y = x[i]*(0.5*Hx + this->
alpha(i)/t) + err;
525 err = (tmp - val) + y;
527 g[i] = Hx + this->
alpha(i)/t;
532 unsigned solveDual_dim1(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
537 unsigned solveDual_dim2(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
541 Real diffg = this->
gx_->dot(*this->
gx_);
542 if ( std::abs(diffg) > ROL_EPSILON<Real>() ) {
543 Real diffa = (this->
alpha(0)-this->
alpha(1))/t;
545 this->
setDualVariables(0,std::min(1.0,std::max(0.0,-(gdiffg+diffa)/diffg)));
549 if ( std::abs(this->
alpha(0)-this->
alpha(1)) > ROL_EPSILON<Real>() ) {
565 void solveSystem(
int size,
char tran, Teuchos::SerialDenseMatrix<int,Real> &
L, Teuchos::SerialDenseVector<int,Real> &v){
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;
575 lapack_.TRTRS(
'L', tran,
'N',
size, 1,
L.values(),
L.stride(), v.values(), v.stride(), &info );
580 bool isFeasible(Teuchos::SerialDenseVector<int,Real> &v,
const Real &tol){
582 for(
int i=0;i<v.numRows();i++){
590 unsigned solveDual_TT(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
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;
608 L(0,0) = std::sqrt(
GiTGj(0,0));
616 z1.resize(1);
z2.resize(1);
627 for (iter=0;iter<maxit;iter++){
647 std::cout <<
"In case 0" << std::endl;
648 std::cout <<
"rho_ = " <<
rho_ << std::endl;
649 std::cout <<
"Solution tempv = \n" <<
tempv << std::endl;
670 bool singular =
false;
702 std::cout <<
"In case 1" << std::endl;
704 std::cout <<
"rho_ = " <<
rho_ << std::endl;
705 std::cout <<
"Solution tempv = \n" <<
tempv << std::endl;
708 std::cout <<
"Direction tempv = \n" <<
tempv << std::endl;
752 std::cout <<
"In case 2" << std::endl;
753 std::cout <<
"Direction tempv = \n" <<
tempv << std::endl;
763 std::cout <<
"Solution tempv is feasible" << std::endl;
773 std::cout <<
"Solution tempv is NOT feasible" << std::endl;
805 std::cout <<
"h = " << h <<
" tempv[h] = " <<
tempv[h] <<
" dualV[base[h]] = " << this->
getDualVariables(base[h]) << std::endl;
810 std::cout <<
"Taking step of size " << step << std::endl;
813 if (step <= 0 || step ==
INF){
815 std::cout <<
"Invalid step!" << std::endl;
829 for (
unsigned baseitem=0;baseitem<
currSize_;baseitem++){
837 std::cout <<
"Blocking " <<
entering_ <<
" because it just entered" << std::endl;
852 std::cout <<
"Returning because nothing deleted and not optimal" << std::endl;
867 newobjval = -
Quad/2;
875 std::cout <<
"New Obj value = " << newobjval << std::endl;
882 if( newobjval >=
objval - std::max( tol*std::abs(
objval), 1.e-16 ) ){
884 std::cout <<
"Blocking " <<
entering_ <<
" because it didn't provide decrease" << std::endl;
897 std::cout <<
"Taboo list has been reset." << std::endl;
908 #if ( ! FIRST_VIOLATED ) 909 Real diff = -
INF, olddiff = -
INF;
912 for (
unsigned bundleitem=0;bundleitem<this->
size();bundleitem++){
916 std::cout <<
"Item " << bundleitem <<
" is blocked." << std::endl;
921 if ( std::find(
base.begin(),
base.end(),bundleitem) ==
base.end() ){
924 std::cout <<
"Base does not contain index " << bundleitem << std::endl;
930 ro += this->
alpha(bundleitem)/t;
933 std::cout <<
"RO = " << ro << std::endl;
946 #if ( FIRST_VIOLATED ) 951 if ( diff > olddiff ){
962 std::cout <<
"entering_ = " <<
entering_ << std::endl;
968 std::cout <<
"Inserting " <<
entering_ << std::endl;
974 std::cout <<
"Current base = " << std::endl;
975 for (
unsigned k=0;k<
base.size();k++){
976 std::cout <<
base[k] << std::endl;
978 std::cout <<
"dependent_ = " <<
dependent_ << std::endl;
986 for (
unsigned i=0;i<zsize;i++){
989 Teuchos::SerialDenseMatrix<int,Real> LBprime( Teuchos::Copy,
L,zsize,zsize);
991 for (
unsigned i=0;i<zsize;i++){
996 Real nrm =
lh.dot(
lh);
1000 std::cout <<
"lh_dot_lh = " << nrm << std::endl;
1001 std::cout <<
"delta = " << delta << std::endl;
1015 for (
unsigned i=0;i<zsize-1;i++){
1019 Real deltaeps = tol*
kappa*std::max(1.0,
lh.dot(
lh));
1021 std::cout <<
"kappa = " <<
kappa << std::endl;
1022 std::cout <<
"deltaeps = " << deltaeps << std::endl;
1024 if ( delta > deltaeps ){
1027 delta = std::sqrt(delta);
1030 else if(delta < -deltaeps){
1032 std::cout <<
"NEGATIVE delta!" << std::endl;
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;
1055 std::cout <<
"Returning because no dual constraint violated and f cannot reach min value " <<
minobjval << std::endl;
1072 unsigned solveDual(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
1074 if (this->
size() == 1) {
1077 else if (this->
size() == 2) {
1082 unsigned outermaxit = 20;
1083 bool increase =
false, decrease =
false;
1085 for (
unsigned it=0; it < outermaxit; it++ ){
1097 if ( (mytol > 1.e-4) || (mytol < 1.e-16) ){
1100 if ( increase && decrease ){
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++){
1109 std::cout <<
"x[" << i <<
"] = " << sol[i] << std::endl;
1111 std::cout << std::endl;
1112 std::vector<Real> g(this->
size(),0.0);
1114 std::cout <<
"and objective value = " << val << std::endl;
1115 std::cout <<
"Checking constraints" << std::endl;
1117 std::cout << success << std::endl;
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++) {
1128 this->
tG_->set(*this->
gx_);
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_);
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++) {
1139 y = x[i]*(t*Hx + this->
alpha(i)) + err;
1141 err = (tmp - v) + y;
1143 g[i] = - t*Hx - this->
alpha(i);
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;
1152 else if ( g[i] < v - myeps ){
1153 std::cout <<
"Constraint " << i <<
" is inactive" << std::endl;
1156 std::cout <<
"Constraint " << i <<
" is active" << std::endl;
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)
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_
const Real alpha(const unsigned i) const
Teuchos::SerialDenseVector< int, Real > tempv
Contains definitions of custom data types in ROL.
const Vector< Real > & subgradient(const unsigned i) const
Teuchos::SerialDenseMatrix< int, Real > Id
bool checkPrimary(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
bool isFeasible(Teuchos::SerialDenseVector< int, Real > &v, const Real &tol)
Teuchos::LAPACK< int, Real > lapack_
Teuchos::SerialDenseVector< int, Real > tempw1
Teuchos::SerialDenseVector< int, Real > z1
Teuchos::SerialDenseVector< int, Real > lj
Bundle_TT(const unsigned maxSize=10, const Real coeff=0.0, const unsigned remSize=2)
Real getDualVariables(const unsigned i)
Teuchos::SerialDenseMatrix< int, Real > L
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
unsigned size(void) const
Teuchos::RCP< Vector< Real > > eG_
void solveSystem(int size, char tran, Teuchos::SerialDenseMatrix< int, Real > &L, Teuchos::SerialDenseVector< int, Real > &v)
Teuchos::RCP< Vector< Real > > tG_
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)
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_
Provides the interface for and implments a bundle.