60 #include "ConstrainedOptPack_QPSchur.hpp"
61 #include "ConstrainedOptPack_ComputeMinMult.hpp"
62 #include "AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
63 #include "AbstractLinAlgPack_SpVectorOp.hpp"
64 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
65 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
66 #include "AbstractLinAlgPack_GenPermMatrixSliceOut.hpp"
67 #include "AbstractLinAlgPack_SpVectorOut.hpp"
68 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
69 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
70 #include "AbstractLinAlgPack_VectorStdOps.hpp"
71 #include "AbstractLinAlgPack_EtaVector.hpp"
72 #include "DenseLinAlgPack_AssertOp.hpp"
73 #include "DenseLinAlgPack_DVectorClass.hpp"
74 #include "DenseLinAlgPack_DVectorClassExt.hpp"
75 #include "DenseLinAlgPack_DVectorOp.hpp"
76 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
77 #include "DenseLinAlgPack_DVectorOut.hpp"
78 #include "DenseLinAlgPack_DMatrixOut.hpp"
79 #include "Teuchos_Workspace.hpp"
80 #include "Teuchos_Assert.hpp"
82 namespace LinAlgOpPack {
95 T my_max(
const T& v1,
const T& v2 ) {
return v1 > v2 ? v1 : v2; }
101 const char* bnd_str( ConstrainedOptPack::EBounds bnd )
104 case ConstrainedOptPack::FREE:
106 case ConstrainedOptPack::UPPER:
108 case ConstrainedOptPack::LOWER:
110 case ConstrainedOptPack::EQUALITY:
121 const char* bool_str(
bool b )
123 return b ?
"true" :
"false";
130 std::string error_msg(
131 const char file_name[],
const int line_num,
const char err_msg[]
134 std::ostringstream omsg;
136 << file_name <<
":" << line_num <<
":" << err_msg;
143 void deincrement_indices(
145 ,std::vector<DenseLinAlgPack::size_type> *indice_vector
150 typedef std::vector<DenseLinAlgPack::size_type> vec_t;
152 for( vec_t::iterator itr = indice_vector->begin(); itr != indice_vector->begin() + len_vector; ++itr ) {
153 if( *itr > k_remove )
161 void insert_pair_sorted(
165 ,std::vector<DenseLinAlgPack::size_type> *r
166 ,std::vector<DenseLinAlgPack::size_type> *c
169 typedef std::vector<DenseLinAlgPack::size_type> rc_t;
173 itr = std::lower_bound( r->begin(), r->begin() + len_vector-1, r_v );
176 {
for( rc_t::iterator itr_last = r->begin() + len_vector-1;
177 itr_last > r->begin() + p; --itr_last )
179 *itr_last = *(itr_last-1);
181 {
for( rc_t::iterator itr_last = c->begin() + len_vector-1;
182 itr_last > c->begin() + p; --itr_last )
184 *itr_last = *(itr_last-1);
196 ,
const DenseLinAlgPack::DVectorSlice &d_hat
198 ,
const DenseLinAlgPack::DVectorSlice *vo
199 ,DenseLinAlgPack::DVectorSlice *z_hat
207 Workspace<DenseLinAlgPack::value_type> t_ws(wss,d_hat.dim());
208 DenseLinAlgPack::DVectorSlice t(&t_ws[0],t_ws.size());
220 ,
const DenseLinAlgPack::DVectorSlice *fo
222 ,
const DenseLinAlgPack::DVectorSlice &z_hat
223 ,DenseLinAlgPack::DVectorSlice *v
226 using DenseLinAlgPack::norm_inf;
232 Workspace<DenseLinAlgPack::value_type> t_ws(wss,v->dim());
233 DenseLinAlgPack::DVectorSlice t(&t_ws[0],t_ws.size());
242 if( norm_inf(t) > 0.0 )
257 ,
const DenseLinAlgPack::DVectorSlice &x
258 ,
const DenseLinAlgPack::DVectorSlice &v
259 ,DenseLinAlgPack::DVectorSlice *mu_D
264 using LinAlgOpPack::V_MtV;
265 using LinAlgOpPack::V_StMtV;
267 using AbstractLinAlgPack::V_MtV;
268 using AbstractLinAlgPack::Vp_MtV;
280 const DenseLinAlgPack::DVectorSlice g = qp.
g();
283 V_StMtV( mu_D, -1.0, Q_XD_hat,
trans, g );
292 const DenseLinAlgPack::DVectorSlice z_hat = act_set.
z_hat();
293 AbstractLinAlgPack::SpVector P_plus_hat_z_hat;
309 ,
const DenseLinAlgPack::DVectorSlice &p_v
310 ,
const DenseLinAlgPack::DVectorSlice &p_z_hat
312 ,DenseLinAlgPack::DVectorSlice *p_mu_D
317 using AbstractLinAlgPack::V_MtV;
318 using AbstractLinAlgPack::Vp_MtV;
336 AbstractLinAlgPack::SpVector Q_R_p_v1;
337 V_MtV( &Q_R_p_v1, qp.
Q_R(),
no_trans, p_v(1,n_R) );
342 AbstractLinAlgPack::SpVector P_XF_hat_p_z_hat;
352 AbstractLinAlgPack::SpVector p_lambda_bar(
357 p_lambda_bar.insert_element(AbstractLinAlgPack::SpVector::element_type(*ja,1.0));
426 template<
class val_type>
429 ,
const DenseLinAlgPack::DVectorSlice &v
430 ,
const DenseLinAlgPack::DVectorSlice &z_hat
431 ,
const DenseLinAlgPack::value_type ao
432 ,
const DenseLinAlgPack::DVectorSlice *bo
434 ,DenseLinAlgPack::value_type *roR_scaling
435 ,DenseLinAlgPack::value_type *rom_scaling
436 ,
const DenseLinAlgPack::value_type aa
437 ,
const DenseLinAlgPack::DVectorSlice *ba
439 ,DenseLinAlgPack::value_type *ra_scaling
444 using DenseLinAlgPack::norm_inf;
445 using DenseLinAlgPack::DVectorSlice;
447 using AbstractLinAlgPack::V_MtV;
448 using AbstractLinAlgPack::SpVector;
452 using LinAlgOpPack::V_MtV;
453 using LinAlgOpPack::Vp_MtV;
456 namespace COP = ConstrainedOptPack;
460 const COP::QPSchurPack::QP
462 const COP::QPSchurPack::Constraints
463 &constraints = qp.constraints();
464 const GenPermMatrixSlice
472 m_breve = constraints.m_breve(),
473 q_hat = act_set.
q_hat(),
480 lambda = ( m ? v(n_R+1,n_R+m): DVectorSlice() ),
481 boR = ( bo ? (*bo)(1,n_R) : DVectorSlice() ),
482 bom = ( bo && m ? (*bo)(n_R+1,n_R+m) : DVectorSlice() );
488 Workspace<DenseLinAlgPack::value_type>
490 DenseLinAlgPack::DVectorSlice
491 x_free(&x_free_ws[0],x_free_ws.size());
500 t1(&t1_ws[0],t1_ws.size()),
501 t2(&t2_ws[0],t2_ws.size()),
502 t3(&t3_ws[0],t3_ws.size()),
503 tR(&tR_ws[0],tR_ws.size()),
504 tm(tm_ws.size()?&tm_ws[0]:NULL,tm_ws.size()),
505 ta(ta_ws.size()?&ta_ws[0]:NULL,ta_ws.size());
514 LinAlgOpPack::V_MtV( &x_free, Q_R,
no_trans, x_R );
516 LinAlgOpPack::Vp_MtV( &x_free, P_XF_hat,
no_trans, z_hat );
520 V_MtV( &lambda_bar, P_plus_hat,
no_trans, z_hat );
530 LinAlgOpPack::V_MtV( &tR, Q_R,
trans, t1 );
531 *roR_scaling += norm_inf(tR);
534 LinAlgOpPack::V_MtV( &tR, Q_R,
trans, t2 );
535 *roR_scaling += norm_inf(tR);
539 LinAlgOpPack::V_MtV( &tR, Q_R,
trans, t3 );
540 *roR_scaling += norm_inf(tR);
545 *roR_scaling += std::fabs(ao) * norm_inf(boR);
550 *rom_scaling += norm_inf(tm);
554 *rom_scaling += std::fabs(ao)*norm_inf(bom);
561 V_StV( ra, aa, *ba );
562 *ra_scaling += std::fabs(aa) * norm_inf(*ba);
568 LinAlgOpPack::V_MtV( &ta, P_XF_hat,
trans, t1 );
569 *ra_scaling += norm_inf(ta);
574 *ra_scaling += norm_inf(ta);
578 LinAlgOpPack::V_MtV( &ta, P_XF_hat,
trans, t2 );
579 *ra_scaling += norm_inf(ta);
582 if( q_F_hat && q_plus_hat ) {
583 LinAlgOpPack::V_MtV( &ta, P_XF_hat,
trans, t3 );
584 *ra_scaling += norm_inf(ta);
588 const GenPermMatrixSlice
591 for( GenPermMatrixSlice::const_iterator itr = P_FC_hat.begin(); itr != P_FC_hat.end(); ++itr ) {
592 ta(itr->row_i()) = z_hat(itr->col_j());
593 ta(itr->col_j()) = z_hat(itr->row_i());
595 *ra_scaling += norm_inf(ta);
616 int correct_dual_infeas(
618 ,
const ConstrainedOptPack::EBounds bnd_j
619 ,
const DenseLinAlgPack::value_type t_P
620 ,
const DenseLinAlgPack::value_type scale
621 ,
const DenseLinAlgPack::value_type dual_infeas_tol
622 ,
const DenseLinAlgPack::value_type degen_mult_val
625 ,
const bool print_dual_infeas
627 ,DenseLinAlgPack::value_type *nu_j
628 ,DenseLinAlgPack::value_type *scaled_viol
629 ,
const char p_nu_j_n[] = NULL
630 ,DenseLinAlgPack::value_type *p_nu_j = NULL
631 ,
const char nu_j_plus_n[] = NULL
632 ,DenseLinAlgPack::value_type *nu_j_plus = NULL
635 typedef DenseLinAlgPack::value_type value_type;
636 namespace COP = ConstrainedOptPack;
638 value_type nu_j_max = (*scaled_viol) = scale * (*nu_j) * (bnd_j == COP::UPPER ? +1.0 : -1.0);
639 if( nu_j_max > 0.0 || bnd_j == COP::EQUALITY )
642 nu_j_max = std::fabs(nu_j_max);
643 if( nu_j_max < dual_infeas_tol ) {
645 value_type degen_val = degen_mult_val * ( bnd_j == COP::UPPER ? +1.0 : -1.0 );
646 if( out && (
int)olevel >= (int)COP::QPSchur::OUTPUT_BASIC_INFO ) {
648 <<
"\nWarning, the constriant a(" << j <<
") currently at its "
649 << (bnd_j == COP::UPPER ?
"UPPER" :
"LOWER") <<
" bound"
650 <<
" has the wrong Lagrange multiplier value but\n"
651 <<
"scale*|"<<nu_j_n<<
"| = " << scale <<
" * |" << (*nu_j)
652 <<
"| = " << nu_j_max <<
" < dual_infeas_tol = " << dual_infeas_tol
653 <<
"\nTherefore, this is considered a degenerate constraint and this "
654 <<
"multiplier is set to " << degen_val << std::endl;
657 nu_j_max += std::fabs( t_P * (*p_nu_j) ) * scale;
658 if( nu_j_max < dual_infeas_tol ) {
660 if( out && (
int)olevel >= (
int)COP::QPSchur::OUTPUT_BASIC_INFO ) {
662 <<
"Also, the maximum full step scale*(|"<<nu_j_n<<
"|+|t_P*"<<p_nu_j_n<<
"|) = "
663 << scale <<
" * (|" << (*nu_j) <<
"| + |" << t_P <<
" * " << (*p_nu_j) <<
"|) = "
664 << nu_j_max <<
" < dual_infeas_tol = " << dual_infeas_tol
665 <<
"\nTherefore, this is considered degenerate and therefore "
666 <<
"seting "<<p_nu_j_n<<
" = 0";
668 *out <<
" and "<< nu_j_plus_n <<
" = " << degen_val;
673 *nu_j_plus = degen_val;
678 *scaled_viol = scale * (*nu_j) * (bnd_j == COP::UPPER ? +1.0 : -1.0);
682 if( print_dual_infeas && out && (
int)olevel >= (
int)COP::QPSchur::OUTPUT_BASIC_INFO ) {
684 <<
"\nError, the constriant a(" << j <<
") currently at its "
685 << (bnd_j == COP::UPPER ?
"UPPER" :
"LOWER") <<
" bound"
686 <<
" has the wrong Lagrange multiplier value and\n"
687 <<
"scale*|"<<nu_j_n<<
"| = " << scale <<
" * |" << (*nu_j)
688 <<
"| = " << nu_j_max <<
" > dual_infeas_tol = " << dual_infeas_tol
689 <<
"\nThis is an indication of instability in the calculations.\n"
690 <<
"The QP algorithm is terminated!\n";
706 void calc_obj_grad_norm_inf(
708 ,
const DenseLinAlgPack::DVectorSlice &x
709 ,DenseLinAlgPack::value_type *qp_grad_norm_inf
717 namespace ConstrainedOptPack {
733 ,
const MatrixOp *A_bar
734 ,
const GenPermMatrixSlice *Q_R
735 ,
const GenPermMatrixSlice *P_XF_hat
736 ,
const GenPermMatrixSlice *P_plus_hat
743 P_XF_hat_ = P_XF_hat;
744 P_plus_hat_ = P_plus_hat;
749 return Q_R_->cols() + ( A_ ? A_->cols() : 0 );
754 return P_plus_hat_->cols();
807 ,
const DVectorSlice& x, value_type b
813 using AbstractLinAlgPack::V_MtV;
815 using LinAlgOpPack::V_MtV;
821 LinAlgOpPack::Vp_MtV_assert_sizes(y->dim(),
rows(),
cols(),M_trans,x.dim());
830 m = A() ? A()->cols() : 0;
846 y2 = m ? (*y)(n_R+1,n_R+m) : DVectorSlice();
851 if( P_XF_hat().nz() )
852 V_MtV( &P_XF_hat_x, P_XF_hat(),
no_trans, x );
854 if(P_plus_hat().nz())
855 V_MtV( &P_plus_hat_x, P_plus_hat(),
no_trans, x );
858 else if(b!=1.0)
Vt_S(&y1,b);
863 if(P_plus_hat().nz())
868 else if(b!=1.0)
Vt_S(&y2,b);
870 if( P_XF_hat().nz() )
887 x2 = m ? x(n_R+1,n_R+m) : DVectorSlice();
891 V_MtV( &Q_R_x1, Q_R(),
no_trans, x1 );
894 else if(b!=1.0)
Vt_S( y, b );
899 if(P_plus_hat().nz())
902 if( m && P_XF_hat().nz() )
912 ,
const SpVectorSlice& x, value_type b
921 using LinAlgOpPack::V_MtV;
922 using AbstractLinAlgPack::V_MtV;
924 using LinAlgOpPack::V_MtV;
930 LinAlgOpPack::Vp_MtV_assert_sizes(y->dim(),
rows(),
cols(),M_trans,x.dim());
939 m = A() ? A()->cols() : 0;
955 y2 = m ? (*y)(n_R+1,n_R+m) : DVectorSlice();
960 if( P_XF_hat().nz() )
961 V_MtV( &P_XF_hat_x, P_XF_hat(),
no_trans, x );
963 if(P_plus_hat().nz())
964 V_MtV( &P_plus_hat_x, P_plus_hat(),
no_trans, x );
967 else if(b!=1.0)
Vt_S(&y1,b);
972 if(P_plus_hat().nz())
977 else if(b!=1.0)
Vt_S(&y2,b);
996 x2 = m ? x(n_R+1,n_R+m) : SpVectorSlice(NULL,0,0,0);
1000 V_MtV( &Q_R_x1, Q_R(),
no_trans, x1 );
1002 if(b ==0.0) *y = 0.0;
1003 else if(b!=1.0)
Vt_S( y, b );
1008 if(P_plus_hat().nz())
1011 if( m && P_XF_hat().nz() )
1021 QPSchur::ActiveSet::ActiveSet(
1022 const schur_comp_ptr_t& schur_comp
1025 :schur_comp_(schur_comp)
1026 ,pivot_tols_(pivot_tols)
1027 ,initialized_(false)
1040 QP& qp, size_type num_act_change,
const int ij_act_change[]
1041 ,
const EBounds bnds[],
bool test,
bool salvage_init_schur_comp
1045 using LinAlgOpPack::V_MtV;
1048 using AbstractLinAlgPack::V_MtV;
1051 using DenseLinAlgPack::sym;
1052 typedef MatrixSymAddDelUpdateable
MSADU;
1053 namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
1077 m_breve = constraints.
m_breve();
1086 if( num_act_change ) {
1087 for( size_type k = 1; k <= num_act_change; ++k ) {
1088 const int ij = ij_act_change[k-1];
1089 const EBounds bnd = bnds[k-1];
1092 if( x_init(-ij) == FREE ) {
1093 std::ostringstream omsg;
1095 <<
"QPSchur::ActiveSet::initialize(...) : Error, "
1096 <<
"The variable x(" << -ij <<
") is not initially fixed and can not "
1097 <<
"be freed by ij_act_change[" << k-1 <<
"]\n";
1098 throw std::invalid_argument( omsg.str() );
1100 if( x_init(-ij) == EQUALITY ) {
1101 std::ostringstream omsg;
1103 <<
"QPSchur::ActiveSet::initialize(...) : Error, "
1104 <<
"The variable x(" << -ij <<
") is equality fixed and therefore can not "
1105 <<
"be freed by ij_act_change[" << k-1 <<
"]\n";
1106 throw std::invalid_argument( omsg.str() );
1114 EBounds x_init_bnd = x_init(ij);
1115 if( x_init_bnd == FREE ) {
1119 else if ( x_init_bnd == EQUALITY ) {
1123 else if( x_init_bnd == bnd ) {
1135 if( ij > n + m_breve ) {
1146 q_D_hat = (n - n_R) - q_F_hat,
1151 q_hat = q_plus_hat + q_F_hat + q_C_hat,
1152 q_hat_max = n_X + n,
1157 ij_map_.resize(q_hat_max);
1158 constr_norm_.resize(q_hat_max);
1159 bnds_.resize(q_hat_max);
1160 d_hat_.resize(q_hat_max);
1162 if( num_act_change ) {
1164 for( size_type k = 1; k <= num_act_change; ++k ) {
1165 const int ij = ij_act_change[k-1];
1166 const EBounds bnd = bnds[k-1];
1170 constr_norm_[s] = 1.0;
1172 d_hat_[s] = - g(-ij);
1179 EBounds x_init_bnd = x_init(ij);
1180 if( x_init_bnd == FREE ) {
1183 constr_norm_[s] = 1.0;
1185 d_hat_[s] = constraints.
get_bnd(ij,bnd);
1192 constr_norm_[s] = 1.0;
1194 d_hat_[s] = - g(ij);
1198 constr_norm_[s] = 1.0;
1200 d_hat_[s] = constraints.
get_bnd(ij,bnd) - b_X(l_x_X_map(ij));
1207 constr_norm_[s] = 1.0;
1209 d_hat_[s] = constraints.
get_bnd(ij,bnd);
1218 P_XF_hat_row_.resize(q_F_hat_max);
1219 P_XF_hat_col_.resize(q_F_hat_max);
1220 P_FC_hat_row_.resize(q_F_hat_max);
1221 P_FC_hat_col_.resize(q_F_hat_max);
1222 P_plus_hat_row_.resize(q_plus_hat_max);
1223 P_plus_hat_col_.resize(q_plus_hat_max);
1229 ij_map_t::const_iterator
1230 ij_itr = ij_map_.begin(),
1231 ij_itr_end = ij_itr + q_hat;
1232 for( size_type s = 1; ij_itr != ij_itr_end; ++ij_itr, ++s ) {
1233 const int ij = *ij_itr;
1236 const size_type i = -ij;
1239 P_XF_hat_row_[k_XF_hat] = i;
1240 P_XF_hat_col_[k_XF_hat] = s;
1243 else if( !(ij <= n && (x_init_ij = x_init(ij)) != FREE ) ) {
1244 const size_type j = ij;
1247 P_plus_hat_row_[k_plus_hat] = j;
1248 P_plus_hat_col_[k_plus_hat] = s;
1255 P_XF_hat_.initialize_and_sort(
1256 n,q_hat,q_F_hat,0,0,GPMSTP::BY_ROW
1257 ,q_F_hat ? &P_XF_hat_row_[0] : NULL
1258 ,q_F_hat ? &P_XF_hat_col_[0] : NULL
1261 P_plus_hat_.initialize_and_sort(
1262 n+m_breve,q_hat,q_plus_hat,0,0,GPMSTP::BY_ROW
1263 ,q_plus_hat ? &P_plus_hat_row_[0] : NULL
1264 ,q_plus_hat ? &P_plus_hat_col_[0] : NULL
1270 throw std::logic_error(
1271 error_msg(__FILE__,__LINE__,
"QPSchur::ActiveSet::initialize(...) : "
1272 "Error, q_C_hat != 0, now supported yet!"));
1275 P_FC_hat_.initialize_and_sort(
1276 q_hat,q_hat,q_C_hat,0,0,GPMSTP::BY_ROW
1277 ,q_C_hat ? &P_FC_hat_row_[0] : NULL
1278 ,q_C_hat ? &P_FC_hat_col_[0] : NULL
1283 Q_XD_hat_row_.resize(q_D_hat_max);
1284 Q_XD_hat_col_.resize(q_D_hat_max);
1289 GenPermMatrixSlice::const_iterator
1290 Q_X_itr = qp.
Q_X().begin();
1291 P_row_t::const_iterator
1292 XF_search = P_XF_hat_row_.begin(),
1293 XF_search_end = XF_search + q_F_hat;
1294 for( size_type l = 1; l <= n_X; ++l, ++Q_X_itr ) {
1295 const size_type i = Q_X_itr->row_i();
1297 for( ; XF_search != XF_search_end && *XF_search < i; ++XF_search ) ;
1298 if( XF_search == XF_search_end || (XF_search < XF_search_end && *XF_search > i) ) {
1301 Q_XD_hat_row_[k_XD_hat] = i;
1302 Q_XD_hat_col_[k_XD_hat] = k_XD_hat + 1;
1308 Q_XD_hat_.initialize(
1309 n,q_D_hat,q_D_hat,0,0,GPMSTP::BY_ROW
1310 , q_D_hat ? &Q_XD_hat_row_[0] : NULL
1311 , q_D_hat ? &Q_XD_hat_col_[0] : NULL
1316 l_fxfx_.resize(q_D_hat_max);
1318 for( size_type k = 0; k < q_D_hat; ++k ) {
1319 l_fxfx_[k] = l_x_X_map(Q_XD_hat_row_[k]);
1340 U_hat_.initialize( &G, m ? &qp.
A() : NULL, &constraints.
A_bar()
1341 , &qp.
Q_R(), &P_XF_hat_, &P_plus_hat_ );
1351 q_plus_hat_ = q_plus_hat;
1356 z_hat_.resize(q_hat_max);
1357 p_z_hat_.resize(q_hat_max);
1358 mu_D_hat_.resize(n_X);
1359 p_mu_D_hat_.resize(n_X);
1361 initialized_ =
true;
1364 if( out && (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
1366 <<
"\nPrint definition of Active-Set before the Schur complement is formed...\n";
1373 DMatrix S_store(q_hat+1,q_hat+1);
1374 DMatrixSliceSym S_sym( S_store(2,q_hat+1,1,q_hat), BLAS_Cpp::lower );
1375 MatrixSymPosDefCholFactor S(&S_store());
1378 , MatrixSymNonsing::DUMMY_ARG );
1384 if( q_F_hat && q_plus_hat ) {
1387 qp_->constraints().A_bar()
1393 if( q_F_hat && q_C_hat ) {
1395 throw std::logic_error(
1396 error_msg(__FILE__,__LINE__,
"QPSchur::ActiveSet::initialize(...) : "
1397 "Error, q_C_hat != 0, now supported yet!"));
1401 if( out && (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
1403 <<
"\nIninitial Schur Complement before it is nonsingular:\n"
1404 <<
"\nS_hat =\nLower triangular part (ignore nonzeros above diagonal)\n"
1409 schur_comp().update_interface().initialize(
1410 S_sym,q_hat_max,
true,MSADU::Inertia( q_plus_hat + q_C_hat, 0, q_F_hat )
1413 catch(
const MSADU::WarnNearSingularUpdateException& excpt) {
1414 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1416 <<
"\nActiveSet::initialize(...) : " << excpt.what()
1420 catch(
const MSADU::SingularUpdateException& excpt) {
1421 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1423 <<
"\nActiveSet::initialize(...) : " << excpt.what()
1426 if(salvage_init_schur_comp) {
1427 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1429 <<
"\nsalvage_init_schur_comp == true\n"
1430 <<
"We will attempt to add as many rows/cols of the "
1431 <<
"initial Schur complement as possible ...\n";
1441 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1443 <<
"\nsalvage_init_schur_comp == false\n"
1444 <<
"We will just throw this singularity exception out of here ...\n";
1449 if( out && (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
1451 <<
"\nSchur Complement after factorization:\n"
1453 << schur_comp().op_interface();
1457 schur_comp().update_interface().set_uninitialized();
1461 initialized_ =
true;
1466 initialized_ =
false;
1478 size_type ja, EBounds bnd_ja,
bool update_steps
1480 ,
bool force_refactorization
1481 ,
bool allow_any_cond
1488 using LinAlgOpPack::V_StMtV;
1496 assert_initialized();
1498 bool wrote_output =
false;
1501 &constraints = qp_->constraints();
1503 if( is_init_fixed(ja) && (*x_init_)(ja) == bnd_ja ) {
1510 const size_type q_hat = this->q_hat();
1511 const size_type sd = s_map(-
int(ja));
1512 const size_type la = qp_->l_x_X_map()(ja);
1515 wrote_output = remove_augmented_element(
1516 sd,force_refactorization
1517 ,MatrixSymAddDelUpdateable::EIGEN_VAL_POS
1518 ,out,output_level,allow_any_cond);
1521 itr = std::lower_bound( P_XF_hat_row_.begin(), P_XF_hat_row_.begin()+q_F_hat_, ja );
1523 const size_type p = itr - P_XF_hat_row_.begin();
1524 std::copy( P_XF_hat_row_.begin() + p + 1, P_XF_hat_row_.begin()+q_F_hat_,
1525 P_XF_hat_row_.begin() + p );
1526 std::copy( P_XF_hat_col_.begin() + p + 1, P_XF_hat_col_.begin()+q_F_hat_,
1527 P_XF_hat_col_.begin() + p );
1530 deincrement_indices( sd, &P_XF_hat_col_, q_F_hat_-1 );
1532 deincrement_indices( sd, &P_FC_hat_col_, q_C_hat_ );
1533 if( q_plus_hat_ > 0 )
1534 deincrement_indices( sd, &P_plus_hat_col_, q_plus_hat_ );
1538 const size_type q_D_hat = this->q_D_hat();
1540 l_fxfx_[q_D_hat] = la;
1542 mu_D_hat_[q_D_hat] = 0.0;
1545 p_mu_D_hat_[q_D_hat] = 1.0;
1547 insert_pair_sorted(ja,q_D_hat+1,q_D_hat+1,&Q_XD_hat_row_,&Q_XD_hat_col_);
1560 value_type d_p = 0.0;
1561 const size_type q_hat = this->q_hat();
1562 Workspace<value_type> t_hat_ws(wss,q_hat);
1563 DVectorSlice t_hat(t_hat_ws.size()?&t_hat_ws[0]:NULL,t_hat_ws.size());
1564 value_type alpha_hat = 0.0;
1565 bool changed_bounds =
false;
1568 if( ja <= n_ && !is_init_fixed(ja) ) {
1576 la = qp_->Q_R().lookup_col_j(ja);
1578 const eta_t u_p = eta_t(la,n_R_+m_);
1584 V_StMtV( &t_hat(), -1.0, U_hat_,
trans, r() );
1586 alpha_hat = -
dot( u_p(), r() );
1588 d_p = constraints.
get_bnd( ja, bnd_ja );
1590 changed_bounds =
false;
1592 else if ( is_init_fixed(ja) ) {
1604 sd = s_map(-
int(ja));
1606 const size_type la = qp_->l_x_X_map()(ja);
1614 d_p = constraints.
get_bnd( ja, bnd_ja ) - qp_->b_X()(la);
1616 changed_bounds =
true;
1624 const eta_t e_ja = eta_t(ja,n_+m_breve_);
1625 const MatrixOp &A_bar = constraints.
A_bar();
1626 DVector u_p( n_R_ + m_ );
1629 u_p(n_R_+1,n_R_+m_) = 0.0;
1636 V_StMtV( &t_hat(), -1.0, U_hat_,
trans, r() );
1637 Vp_StPtMtV( &t_hat(), 1.0, P_XF_hat_,
trans, A_bar, no_trans, e_ja() );
1640 alpha_hat = -
dot( u_p(), r() );
1644 d_p = constraints.
get_bnd( ja, bnd_ja );
1647 r.resize( n_ - n_R_ );
1649 d_p += -
dot( qp_->b_X(), r() );
1652 changed_bounds =
false;
1660 schur_comp().update_interface().augment_update(
1661 &t_hat(), alpha_hat, force_refactorization
1662 ,MatrixSymAddDelUpdateable::EIGEN_VAL_NEG
1663 ,MSADU::PivotTolerances(
1665 ,allow_any_cond ? 0.0 :
pivot_tols().singular_tol
1666 ,allow_any_cond ? 0.0 :
pivot_tols().wrong_inertia_tol ) );
1668 catch(
const MSADU::WarnNearSingularUpdateException& excpt) {
1669 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1671 <<
"\nActiveSet::add_constraint(...) : " << excpt.what()
1673 wrote_output =
true;
1678 schur_comp().update_interface().initialize(
1679 alpha_hat, (n_-n_R_) + n_-m_ );
1682 if( changed_bounds )
1686 const size_type q_hat_new = q_F_hat_ + q_C_hat_ + q_plus_hat_;
1688 ij_map_[q_hat_new - 1] = ja;
1690 constr_norm_[q_hat_new - 1] = 1.0;
1692 bnds_[q_hat_new - 1] = bnd_ja;
1694 d_hat_(q_hat_new) = d_p;
1696 z_hat_(q_hat_new) = 0.0;
1697 if( update_steps ) {
1699 p_z_hat_(q_hat_new) = 1.0;
1701 if( !changed_bounds ) {
1703 insert_pair_sorted(ja,q_hat_new,q_plus_hat_,&P_plus_hat_row_,&P_plus_hat_col_);
1708 insert_pair_sorted(sd,q_hat_new,q_C_hat_,&P_FC_hat_row_,&P_FC_hat_col_);
1712 reinitialize_matrices(test_);
1713 return wrote_output;
1718 ,
bool force_refactorization,
bool allow_any_cond
1725 using LinAlgOpPack::V_StMtV;
1726 using LinAlgOpPack::V_MtV;
1727 using LinAlgOpPack::Vp_MtV;
1736 assert_initialized();
1738 bool wrote_output =
false;
1746 q_hat = this->q_hat(),
1747 q_F_hat = this->q_F_hat(),
1748 q_plus_hat = this->q_plus_hat(),
1749 q_D_hat = this->q_D_hat();
1751 const size_type
id = -jd;
1753 const size_type ld = qp_->l_x_X_map()(-jd);
1756 {
for( kd = 1; kd <= q_D_hat; ++kd ) {
1757 if( l_fxfx_[kd-1] == ld )
break;
1766 &A_bar = qp_->constraints().A_bar();
1767 const MatrixSymOpNonsing
1770 &U_hat = this->U_hat();
1771 const GenPermMatrixSlice
1774 &P_XF_hat = this->P_XF_hat(),
1775 &P_plus_hat = this->P_plus_hat();
1784 DVector u_p(n_R_+m_);
1787 V_MtV( &u_p(n_R_+1,n_R_+m_), qp_->A(),
trans, e_id() );
1789 nrm_u_p = DenseLinAlgPack::norm_inf( u_p() );
1794 DVector Q_X_G_e_id(Q_X.cols());
1797 d_p = -g(
id) -
dot( b_X, Q_X_G_e_id() );
1804 Workspace<value_type>
1805 t_hat_ws(wss,q_hat);
1807 t_hat(&t_hat_ws[0],q_hat);
1817 Vp_MtV( &t_hat(), U_hat,
trans, r() );
1821 alpha_hat = sigma - ( nrm_u_p > 0.0 ?
dot(u_p(),r()) : 0.0 );
1829 schur_comp().update_interface().augment_update(
1830 &t_hat(), alpha_hat, force_refactorization
1831 ,MatrixSymAddDelUpdateable::EIGEN_VAL_POS );
1833 catch(
const MSADU::WarnNearSingularUpdateException& excpt) {
1834 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1836 <<
"\nActiveSet::drop_constraint(...) : " << excpt.what()
1838 wrote_output =
true;
1843 schur_comp().update_interface().initialize(
1844 alpha_hat, (n_-n_R_) + n_-m_ );
1850 std::copy( l_fxfx_.begin() + kd, l_fxfx_.begin() + q_D_hat
1851 , l_fxfx_.begin() + (kd-1) );
1853 std::copy( mu_D_hat_.begin() + kd, mu_D_hat_.begin() + q_D_hat
1854 , mu_D_hat_.begin() + (kd-1) );
1857 itr = std::lower_bound( Q_XD_hat_row_.begin(), Q_XD_hat_row_.begin()+q_D_hat, id );
1859 const size_type p = itr - Q_XD_hat_row_.begin();
1860 std::copy( Q_XD_hat_row_.begin() + p + 1, Q_XD_hat_row_.begin()+q_D_hat,
1861 Q_XD_hat_row_.begin() + p );
1862 std::copy( Q_XD_hat_col_.begin() + p + 1, Q_XD_hat_col_.begin()+q_D_hat,
1863 Q_XD_hat_col_.begin() + p );
1865 deincrement_indices( kd, &Q_XD_hat_col_, q_D_hat-1 );
1870 q_hat = this->q_hat();
1875 ij_map_[q_hat-1] = -id;
1879 bnds_[q_hat-1] = FREE;
1881 d_hat_[q_hat-1] = d_p;
1883 z_hat_[q_hat-1] = 0.0;
1885 insert_pair_sorted(
id,q_hat,q_F_hat_,&P_XF_hat_row_,&P_XF_hat_col_);
1891 const size_type q_hat = this->q_hat();
1892 const size_type sd = s_map(jd);
1894 wrote_output = remove_augmented_element(
1895 sd,force_refactorization
1896 ,MatrixSymAddDelUpdateable::EIGEN_VAL_NEG
1897 ,out,output_level,allow_any_cond
1899 if( is_init_fixed(jd) ) {
1902 const size_type sd1 = s_map(-jd);
1906 itr = std::lower_bound( P_FC_hat_row_.begin(), P_FC_hat_row_.begin()+q_C_hat_, sd1 );
1908 const size_type p = itr - P_FC_hat_row_.begin();
1909 std::copy( P_FC_hat_row_.begin() + p + 1, P_FC_hat_row_.begin()+q_C_hat_,
1910 P_FC_hat_row_.begin() + p );
1911 std::copy( P_FC_hat_col_.begin() + p + 1, P_FC_hat_col_.begin()+q_C_hat_,
1912 P_FC_hat_col_.begin() + p );
1918 itr = std::lower_bound( P_plus_hat_row_.begin(), P_plus_hat_row_.begin()+q_plus_hat_, jd );
1920 const size_type p = itr - P_plus_hat_row_.begin();
1921 std::copy( P_plus_hat_row_.begin() + p + 1, P_plus_hat_row_.begin()+q_plus_hat_,
1922 P_plus_hat_row_.begin() + p );
1923 std::copy( P_plus_hat_col_.begin() + p + 1, P_plus_hat_col_.begin()+q_plus_hat_,
1924 P_plus_hat_col_.begin() + p );
1929 deincrement_indices( sd, &P_XF_hat_col_, q_F_hat_ );
1931 deincrement_indices( sd, &P_FC_hat_col_, q_C_hat_ );
1932 if( q_plus_hat_ > 0 )
1933 deincrement_indices( sd, &P_plus_hat_col_, q_plus_hat_ );
1936 reinitialize_matrices(test_);
1937 return wrote_output;
1941 int jd, size_type ja, EBounds bnd_ja,
bool update_steps
1945 bool wrote_output =
false;
1946 if( drop_constraint( jd, out, output_level,
false,
true ) )
1947 wrote_output =
true;
1948 if( add_constraint( ja, bnd_ja, update_steps, out, output_level,
true,
true ) )
1949 wrote_output =
true;
1950 return wrote_output;
1956 assert_initialized();
1963 assert_initialized();
1969 assert_initialized();
1970 return q_plus_hat_ + q_F_hat_ + q_C_hat_;
1975 assert_initialized();
1981 assert_initialized();
1987 assert_initialized();
1993 assert_initialized();
1994 return (n_ - n_R_) - q_F_hat_;
2000 return ij_map_[s-1];
2005 ij_map_t::const_iterator
2006 begin = ij_map_.begin(),
2007 end = begin + q_hat(),
2008 itr = std::find( begin, end, ij );
2009 return ( itr != end ? (itr - begin) + 1 : 0 );
2015 return constr_norm_(s);
2027 return l_fxfx_[k-1];
2032 assert_initialized();
2038 assert_initialized();
2039 return schur_comp().op_interface();
2044 assert_initialized();
2050 assert_initialized();
2056 assert_initialized();
2062 assert_initialized();
2068 assert_initialized();
2069 return d_hat_(1,q_hat());
2074 assert_initialized();
2075 return z_hat_(1,q_hat());
2080 assert_initialized();
2081 return z_hat_(1,q_hat());
2086 assert_initialized();
2087 return p_z_hat_(1,q_hat());
2092 assert_initialized();
2093 return p_z_hat_(1,q_hat());
2098 assert_initialized();
2099 return mu_D_hat_(1,q_D_hat());
2104 assert_initialized();
2105 return mu_D_hat_(1,q_D_hat());
2110 assert_initialized();
2111 return p_mu_D_hat_(1,q_D_hat());
2116 assert_initialized();
2117 return p_mu_D_hat_(1,q_D_hat());
2122 assert_initialized();
2123 return j <= n_ && (*x_init_)(j) != FREE;
2128 return n_ - m_ == (n_ - n_R_) - q_F_hat_ + q_C_hat_ + q_plus_hat_;
2133 void QPSchur::ActiveSet::assert_initialized()
const
2136 throw std::logic_error(
2137 error_msg(__FILE__,__LINE__,
"QPSchur::ActiveSet::assert_initialized() : Error, "
2138 "The active set has not been initialized yet!") );
2141 void QPSchur::ActiveSet::assert_s( size_type s)
const
2146 void QPSchur::ActiveSet::reinitialize_matrices(
bool test)
2148 namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
2150 const size_type q_hat = this->q_hat();
2151 const size_type q_D_hat = this->q_D_hat();
2153 P_XF_hat_.initialize(
2154 n_,q_hat,q_F_hat_,0,0,GPMSTP::BY_ROW
2155 ,q_F_hat_ ? &P_XF_hat_row_[0] : NULL
2156 ,q_F_hat_ ? &P_XF_hat_col_[0] : NULL
2159 P_FC_hat_.initialize(
2160 q_hat,q_hat,q_C_hat_,0,0,GPMSTP::BY_ROW
2161 ,q_C_hat_ ? &P_FC_hat_row_[0] : NULL
2162 ,q_C_hat_ ? &P_FC_hat_col_[0] : NULL
2165 P_plus_hat_.initialize(
2166 n_+m_breve_,q_hat,q_plus_hat_,0,0,GPMSTP::BY_ROW
2167 ,q_plus_hat_ ? &P_plus_hat_row_[0] : NULL
2168 ,q_plus_hat_ ? &P_plus_hat_col_[0] : NULL
2171 Q_XD_hat_.initialize(
2172 n_,q_D_hat,q_D_hat,0,0,GPMSTP::BY_ROW
2173 ,q_D_hat ? &Q_XD_hat_row_[0] : NULL
2174 ,q_D_hat ? &Q_XD_hat_col_[0] : NULL
2178 &qp_->G(), m_ ? &qp_->A() : NULL, &qp_->constraints().A_bar()
2179 ,&qp_->Q_R(), &P_XF_hat_, &P_plus_hat_);
2182 bool QPSchur::ActiveSet::remove_augmented_element(
2183 size_type sd,
bool force_refactorization
2184 ,MatrixSymAddDelUpdateable::EEigenValType eigen_val_drop
2186 ,
bool allow_any_cond
2189 bool wrote_output =
false;
2190 const size_type q_hat = this->q_hat();
2193 schur_comp().update_interface().delete_update(
2194 sd,force_refactorization,eigen_val_drop
2195 ,MSADU::PivotTolerances(
2197 ,allow_any_cond ? 0.0 :
pivot_tols().singular_tol
2198 ,allow_any_cond ? 0.0 :
pivot_tols().wrong_inertia_tol ));
2200 catch(
const MSADU::WarnNearSingularUpdateException& excpt) {
2201 if( out && (
int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
2203 <<
"\nActiveSet::drop_constraint(...) : " << excpt.what()
2205 wrote_output =
true;
2209 std::copy( ij_map_.begin() + sd, ij_map_.begin() + q_hat
2210 , ij_map_.begin() + (sd-1) );
2212 std::copy( constr_norm_.begin() + sd, constr_norm_.begin() + q_hat
2213 , constr_norm_.begin() + (sd-1) );
2215 std::copy( bnds_.begin() + sd, bnds_.begin() + q_hat
2216 , bnds_.begin() + (sd-1) );
2218 std::copy( d_hat_.begin() + sd, d_hat_.begin() + q_hat
2219 , d_hat_.begin() + (sd-1) );
2221 std::copy( z_hat_.begin() + sd, z_hat_.begin() + q_hat
2222 , z_hat_.begin() + (sd-1) );
2224 std::copy( p_z_hat_.begin() + sd, p_z_hat_.begin() + q_hat
2225 , p_z_hat_.begin() + (sd-1) );
2226 return wrote_output;
2235 act_set_.pivot_tols(pivot_tols);
2240 return act_set_.pivot_tols();
2244 const schur_comp_ptr_t& schur_comp
2246 ,value_type max_real_runtime
2247 ,value_type feas_tol
2248 ,value_type loose_feas_tol
2249 ,value_type dual_infeas_tol
2250 ,value_type huge_primal_step
2251 ,value_type huge_dual_step
2252 ,value_type warning_tol
2253 ,value_type error_tol
2254 ,size_type iter_refine_min_iter
2255 ,size_type iter_refine_max_iter
2256 ,value_type iter_refine_opt_tol
2257 ,value_type iter_refine_feas_tol
2258 ,
bool iter_refine_at_solution
2259 ,
bool salvage_init_schur_comp
2262 :schur_comp_(schur_comp)
2263 ,max_iter_(max_iter)
2264 ,max_real_runtime_(max_real_runtime)
2265 ,feas_tol_(feas_tol)
2266 ,loose_feas_tol_(loose_feas_tol)
2267 ,dual_infeas_tol_(dual_infeas_tol)
2268 ,huge_primal_step_(huge_primal_step)
2269 ,huge_dual_step_(huge_dual_step)
2270 ,warning_tol_(warning_tol)
2271 ,error_tol_(error_tol)
2272 ,iter_refine_min_iter_(iter_refine_min_iter)
2273 ,iter_refine_max_iter_(iter_refine_max_iter)
2274 ,iter_refine_opt_tol_(iter_refine_opt_tol)
2275 ,iter_refine_feas_tol_(iter_refine_feas_tol)
2276 ,iter_refine_at_solution_(iter_refine_at_solution)
2277 ,salvage_init_schur_comp_(salvage_init_schur_comp)
2278 ,act_set_(schur_comp,pivot_tols)
2283 ,size_type num_act_change,
const int ij_act_change[],
const EBounds bnds[]
2285 ,DVectorSlice* x, SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
2286 ,size_type* iter, size_type* num_adds, size_type* num_drops
2292 using DenseLinAlgPack::norm_inf;
2293 using AbstractLinAlgPack::norm_inf;
2299 const value_type inf = std::numeric_limits<value_type>::max();
2302 output_level = NO_OUTPUT;
2304 const int dbl_min_w = 20;
2305 const int dbl_w = (out ? my_max(dbl_min_w,
int(out->precision()+8)) : 20 );
2308 act_set_.set_schur_comp( schur_comp_ );
2311 solve_return = SUBOPTIMAL_POINT;
2314 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2316 <<
"\n*** Entering QPSchur::solve_qp(...)\n";
2325 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2327 <<
"\n*** Dump the definition of the QP to be solved ***\n";
2332 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2334 <<
"\n*** Warm start info\n"
2335 <<
"\nNumber of variables = " << qp.
n()
2336 <<
"\nNumber of initially fixed variables (not in Ko) = " << (qp.
n() - qp.
n_R())
2337 <<
"\nNumber of changes to the initial KKT system (num_act_change) = " << num_act_change << endl;
2341 num_var_fixed = 0, num_var_freed = 0, num_gen_equ = 0, num_gen_inequ = 0;
2342 for( size_type s = 1; s <= num_act_change; ++s ) {
2343 const int ij = ij_act_change[s-1];
2344 const EBounds bnd = bnds[s-1];
2349 else if( bnd == EQUALITY )
2355 <<
"\n Number of initially fixed variables freed from a bound = " << num_var_freed
2356 <<
"\n Number of initially free variables fixed to a bound = " << num_var_fixed
2357 <<
"\n Number of general equality constraints added = " << num_gen_equ
2358 <<
"\n Number of general inequality constraints added = " << num_gen_inequ << endl;
2361 if( num_act_change > 0 && (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
2364 << right << setw(5) <<
"s"
2365 << right << setw(20) <<
"ij_act_change"
2366 << right << setw(10) <<
"bnds" << endl
2367 << right << setw(5) <<
"---"
2368 << right << setw(20) <<
"-------------"
2369 << right << setw(10) <<
"--------" << endl;
2370 for( size_type s = 1; s <= num_act_change; ++s )
2372 << right << setw(5) << s
2373 << right << setw(20) << ij_act_change[s-1]
2374 << right << setw(10) << bnd_str(bnds[s-1]) << endl;
2379 act_set_.
initialize( qp, num_act_change, ij_act_change, bnds
2380 , test_what == RUN_TESTS, salvage_init_schur_comp(), out, output_level );
2384 catch(
const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
2385 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2387 <<
"\n*** Error in initializing schur complement\n"
2388 << excpt.what() << std::endl
2389 <<
"\nSetting num_act_change = 0 and proceeding with a cold start...\n";
2391 act_set_.
initialize( qp, num_act_change = 0, ij_act_change, bnds
2392 , test_what == RUN_TESTS, salvage_init_schur_comp(), out, output_level );
2396 Workspace<value_type> vo_ws(wss,qp.
n_R()+qp.
m());
2397 DVectorSlice vo(&vo_ws[0],vo_ws.size());
2400 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2402 <<
"\nSolution to the initial KKT system, vo = inv(Ko)*fo:\n\n||vo||inf = " << norm_inf(vo) << std::endl;
2404 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2406 <<
"\nvo =\n" << vo;
2429 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2432 <<
"\n*** Removing constriants until we are dual feasible"
2434 <<
"\n*** Start by removing constraints within the Schur complement first\n";
2437 if( (
int)OUTPUT_ITER_SUMMARY <= (
int)output_level
2438 && (
int)output_level < (
int)OUTPUT_ITER_QUANTITIES
2439 && num_act_change > 0 )
2442 <<
"\nIf max_viol > 0 and jd != 0 then constraint jd will be dropped from the active set\n\n"
2443 << right << setw(dbl_w) <<
"max_viol"
2444 << right << setw(5) <<
"sd"
2445 << right << setw(5) <<
"jd" << endl
2446 << right << setw(dbl_w) <<
"--------------"
2447 << right << setw(5) <<
"----"
2448 << right << setw(5) <<
"----" << endl;
2450 for(
int k = num_act_change; k > 0; --k, ++(*iter) ) {
2453 return MAX_RUNTIME_EXEEDED_FAIL;
2455 DVectorSlice z_hat = act_set_.
z_hat();
2459 value_type max_viol = 0.0;
2461 DVectorSlice::iterator
2462 z_itr = z_hat.begin();
2463 const size_type q_hat = act_set_.
q_hat();
2465 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2467 <<
"\nLooking for a constraint with the maximum dual infeasiblity to drop...\n\n"
2468 << right << setw(5) <<
"s"
2469 << right << setw(dbl_w) <<
"z_hat"
2470 << right << setw(20) <<
"bnd"
2471 << right << setw(dbl_w) <<
"viol"
2472 << right << setw(dbl_w) <<
"max_viol"
2473 << right << setw(5) <<
"jd" << endl
2474 << right << setw(5) <<
"----"
2475 << right << setw(dbl_w) <<
"--------------"
2476 << right << setw(20) <<
"--------------"
2477 << right << setw(dbl_w) <<
"--------------"
2478 << right << setw(dbl_w) <<
"--------------"
2479 << right << setw(5) <<
"----" << endl;
2481 for(
int s = 1; s <= q_hat; ++s, ++z_itr ) {
2482 int j = act_set_.
ij_map(s);
2485 EBounds bnd = act_set_.
bnd(s);
2487 const int dual_feas_status
2488 = correct_dual_infeas(
2490 ,out,output_level,false
2491 ,
"z_hat(s)", &(*z_itr), &viol );
2492 if( dual_feas_status < 0 && viol < max_viol ) {
2497 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2499 << right << setw(5) << s
2500 << right << setw(dbl_w) << *z_itr
2501 << right << setw(20) << bnd_str(bnd)
2502 << right << setw(dbl_w) << viol
2503 << right << setw(dbl_w) << max_viol
2504 << right << setw(5) << jd << endl;
2509 if( (
int)OUTPUT_ITER_SUMMARY <= (
int)output_level
2510 && (
int)output_level < (
int)OUTPUT_ITER_QUANTITIES )
2513 << right << setw(dbl_w) << max_viol
2514 << right << setw(5) << act_set_.
s_map(jd)
2515 << right << setw(5) << jd << endl;
2517 if( jd == 0 )
break;
2523 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2525 <<
"\nPrinting active set after dropping constraint jd = " << jd <<
" ...\n";
2531 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2533 <<
"\nThere where " << (*num_drops)
2534 <<
" constraints dropped from the schur complement from the initial guess of the active set.\n";
2538 Workspace<value_type> v_ws(wss,qp.
n_R()+qp.
m());
2539 DVectorSlice v(&v_ws[0],v_ws.size());
2540 if( act_set_.
q_hat() > 0 ) {
2541 calc_v( qp.
Ko(), &qp.
fo(), act_set_.
U_hat(), act_set_.
z_hat(), &v );
2542 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2544 <<
"\nSolution to the system; v = inv(Ko)*(fo - U_hat*z_hat):\n"
2545 <<
"\n||v||inf = " << norm_inf(v) << std::endl;
2547 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2557 set_x( act_set_, v, x );
2558 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2560 <<
"\nCurrent guess for unknowns x:\n\n||x||inf = " << norm_inf(*x) << std::endl;
2562 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2571 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2572 *out <<
"\n*** Second, free initially fixed variables not in Ko\n";
2577 if( (
int)OUTPUT_ITER_SUMMARY <= (
int)output_level
2578 && (
int)output_level < (
int)OUTPUT_ITER_QUANTITIES )
2581 <<
"\nIf max_viol > 0 and id != 0 then the variable x(id) will be freed from its initial bound\n\n"
2582 << right << setw(dbl_w) <<
"max_viol"
2583 << right << setw(5) <<
"kd"
2584 << right << setw(5) <<
"id" << endl
2585 << right << setw(dbl_w) <<
"--------------"
2586 << right << setw(5) <<
"----"
2587 << right << setw(5) <<
"----" << endl;
2589 size_type q_D_hat = act_set_.
q_D_hat();
2590 while( q_D_hat > 0 ) {
2593 return MAX_RUNTIME_EXEEDED_FAIL;
2595 DVectorSlice mu_D_hat = act_set_.
mu_D_hat();
2596 calc_mu_D( act_set_, *x, v, &mu_D_hat );
2598 value_type max_viol = 0.0;
2601 DVectorSlice::iterator
2602 mu_D_itr = mu_D_hat.begin();
2604 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2606 <<
"\nLooking for a variable bound with the max dual infeasibility to drop...\n\n"
2607 << right << setw(5) <<
"k"
2608 << right << setw(dbl_w) <<
"mu_D_hat"
2609 << right << setw(20) <<
"bnd"
2610 << right << setw(dbl_w) <<
"viol"
2611 << right << setw(dbl_w) <<
"max_viol"
2612 << right << setw(5) <<
"id" << endl
2613 << right << setw(5) <<
"----"
2614 << right << setw(dbl_w) <<
"--------------"
2615 << right << setw(20) <<
"--------------"
2616 << right << setw(dbl_w) <<
"--------------"
2617 << right << setw(dbl_w) <<
"--------------"
2618 << right << setw(5) <<
"----" << endl;
2620 for(
int k = 1; k <= q_D_hat; ++k, ++mu_D_itr ) {
2622 i = i_x_X_map(act_set_.
l_fxfx(k));
2626 const int dual_feas_status
2627 = correct_dual_infeas(
2629 ,out,output_level,
false
2630 ,
"mu_D_hat(k)", &(*mu_D_itr), &viol );
2631 if( dual_feas_status < 0 && viol < max_viol ) {
2637 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2639 << right << setw(5) << k
2640 << right << setw(dbl_w) << *mu_D_itr
2641 << right << setw(20) << bnd_str(bnd)
2642 << right << setw(dbl_w) << viol
2643 << right << setw(dbl_w) << max_viol
2644 << right << setw(5) <<
id << endl;
2648 if( (
int)OUTPUT_ITER_SUMMARY <= (
int)output_level
2649 && (
int)output_level < (
int)OUTPUT_ITER_QUANTITIES )
2652 << right << setw(dbl_w) << max_viol
2653 << right << setw(5) << kd
2654 << right << setw(5) <<
id << endl;
2656 if(
id == 0 )
break;
2663 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2665 <<
"\nPrinting active set after freeing initially fixed variable id = " <<
id <<
" ...\n";
2668 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
2670 <<
"\nSolution to the new KKT system; z_hat = inv(S_hat)*(d_hat - U_hat'*vo), v = inv(Ko)*(fo - U_hat*z_hat):\n";
2674 , &act_set_.
z_hat() );
2675 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2677 <<
"\n||z_hat||inf = " << norm_inf(act_set_.
z_hat()) << std::endl;
2679 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2681 <<
"\nz_hat =\n" << act_set_.
z_hat();
2684 calc_v( qp.
Ko(), &qp.
fo(), act_set_.
U_hat(), act_set_.
z_hat(), &v );
2685 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2687 <<
"\n||v||inf = " << norm_inf(v) << std::endl;
2689 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2694 set_x( act_set_, v, x );
2695 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
2697 <<
"\nCurrent guess for unknowns x:\n\n||x||inf = " << norm_inf(*x) << std::endl;
2699 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2707 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2709 <<
"\nThere where " << (*num_adds)
2710 <<
" initially fixed variables not in Ko that were freed and added to the schur complement.\n";
2714 size_type iter_refine_num_resid = 0, iter_refine_num_solves = 0;
2716 PICK_VIOLATED_CONSTRAINT
2717 ,out, output_level, test_what
2719 ,x, iter, num_adds, num_drops
2720 ,&iter_refine_num_resid, &iter_refine_num_solves
2724 if( solve_return != OPTIMAL_SOLUTION )
2725 set_x( act_set_, v, x );
2728 if( solve_return != SUBOPTIMAL_POINT && act_set_.
q_hat() ) {
2729 const size_type q_hat = act_set_.
q_hat();
2730 DVectorSlice z_hat = act_set_.
z_hat();
2731 for( size_type s = 1; s <= q_hat; ++s ) {
2732 const int j = act_set_.
ij_map(s);
2733 value_type viol = 0.0;
2734 const EBounds bnd = act_set_.
bnd(s);
2737 const int dual_feas_status
2738 = correct_dual_infeas(
2740 ,out,output_level,
true,
"z_hat(s)",&z_hat(s),&viol );
2741 if( dual_feas_status < 0 ) {
2742 solve_return = SUBOPTIMAL_POINT;
2747 if( solve_return != SUBOPTIMAL_POINT && act_set_.
q_D_hat() ) {
2748 const GenPermMatrixSlice& Q_XD_hat = act_set_.
Q_XD_hat();
2749 DVectorSlice mu_D_hat = act_set_.
mu_D_hat();
2751 for( GenPermMatrixSlice::const_iterator itr = Q_XD_hat.begin(); itr != Q_XD_hat.end(); ++itr ) {
2752 const int i = itr->row_i();
2753 value_type viol = 0.0;
2754 const EBounds bnd = x_init(i);
2756 const int dual_feas_status
2757 = correct_dual_infeas(
2759 ,out,output_level,
true,
"mu_D_hat(k)",&mu_D_hat(itr->col_j()),&viol );
2760 if( dual_feas_status < 0 ) {
2761 solve_return = SUBOPTIMAL_POINT;
2769 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2770 switch(solve_return) {
2771 case OPTIMAL_SOLUTION:
2773 <<
"\n*** Solution found!\n";
2775 case MAX_ITER_EXCEEDED:
2777 <<
"\n*** Maximum iterations exceeded!\n";
2779 case MAX_RUNTIME_EXEEDED_DUAL_FEAS:
2781 <<
"\n*** Maximum runtime exceeded!\n";
2783 case MAX_ALLOWED_STORAGE_EXCEEDED:
2785 <<
"\n*** The maxinum size of the schur complement has been exceeded!\n";
2787 case INFEASIBLE_CONSTRAINTS:
2789 <<
"\n*** The constraints are infeasible!\n";
2793 <<
"\n*** The QP appears to be nonconvex but will return current point anyway!\n";
2795 case DUAL_INFEASIBILITY:
2797 <<
"\n*** The dual variables are infeasible (numerical roundoff?)!\n";
2799 case SUBOPTIMAL_POINT:
2801 <<
"\n*** The current point is suboptimal but we will return it anyway!\n";
2806 *out <<
"\nNumber of QP iteratons = " << *iter;
2807 *out <<
"\nNumber of iterative refinement residual calculations = " << iter_refine_num_resid;
2808 *out <<
"\nNumber of iterative refinement solves = " << iter_refine_num_solves << endl;
2809 *out <<
"\n||x||inf = " << norm_inf(*x);
2810 *out <<
"\nmu.nz() = " << mu->nz();
2811 *out <<
"\nmax(|mu(i)|) = " << norm_inf((*mu)());
2812 *out <<
"\nmin(|mu(i)|) = " << min_abs((*mu)());
2815 <<
"\nmax(|lambda(i)|) = " << norm_inf(*lambda)
2816 <<
"\nmin(|lambda(i)|) = " << min_abs(*lambda);
2819 <<
"\nlambda_breve.nz() = " << lambda_breve->nz()
2820 <<
"\nmax(|lambda_breve(i)|) = " << norm_inf((*lambda_breve)())
2821 <<
"\nmin(|lambda_breve(i)|) = " << min_abs((*lambda_breve)());
2825 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2826 *out <<
"\nx =\n" << (*x)();
2827 *out <<
"\nmu =\n" << (*mu)();
2829 *out <<
"\nlambda =\n" << (*lambda)();
2831 *out <<
"\nlambda_breve =\n" << (*lambda_breve)();
2834 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2836 <<
"\n*** Leaving QPSchur::solve_qp(...)\n";
2839 return solve_return;
2852 ,
const DVectorSlice& vo,
ActiveSet* act_set, DVectorSlice* v
2853 , DVectorSlice* x, size_type* iter, size_type* num_adds, size_type* num_drops
2854 , size_type* iter_refine_num_resid, size_type* iter_refine_num_solves
2864 using DenseLinAlgPack::norm_inf;
2871 using LinAlgOpPack::V_StMtV;
2873 using AbstractLinAlgPack::norm_inf;
2874 using AbstractLinAlgPack::V_MtV;
2876 using LinAlgOpPack::V_MtV;
2884 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
2886 <<
"\n*** Starting Primal-Dual Iterations ***\n";
2889 const int dbl_min_w = 20;
2890 const int dbl_w = (out ? my_max(dbl_min_w,
int(out->precision()+8)): 20 );
2895 &qp = act_set->
qp();
2902 Workspace<value_type>
2903 v_plus_ws(wss,v->dim()),
2904 z_hat_plus_ws(wss,(n-m)+(n-n_R)),
2905 p_v_ws(wss,v->dim());
2907 v_plus(&v_plus_ws[0],v_plus_ws.size()),
2909 p_v(&p_v_ws[0],p_v_ws.size());
2912 inf = std::numeric_limits<value_type>::max();
2918 size_type last_ja = 0;
2919 value_type con_ja_val;
2921 value_type norm_2_constr;
2924 bool assume_lin_dep_ja;
2925 value_type gamma_plus;
2927 const int summary_lines_counter_max = 15;
2928 int summary_lines_counter = 0;
2931 long int last_jd = 0;
2937 bool warned_degeneracy =
false;
2939 value_type dual_infeas_scale = 1.0;
2941 bool return_to_init_fixed =
false;
2944 bool using_iter_refinement =
false;
2946 for( itr = 0; itr <= max_iter_; ++itr, ++(*iter) ) {
2947 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
2949 <<
"\n************************************************"
2950 <<
"\n*** qp_iter = " << itr
2951 <<
"\n*** q_hat = " << act_set->
q_hat() << std::endl;
2953 bool schur_comp_update_failed =
false;
2954 switch( next_step ) {
2955 case PICK_VIOLATED_CONSTRAINT: {
2956 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
2958 <<
"\n*** PICK_VIOLATED_CONSTRAINT\n";
2962 return MAX_RUNTIME_EXEEDED_DUAL_FEAS;
2969 set_x( *act_set, *v, x );
2970 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
2972 <<
"\n||x||inf = " << norm_inf(*x) << std::endl;
2974 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
2978 if( test_what == RUN_TESTS ) {
2979 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
2981 <<
"\nChecking current iteration quantities ...\n";
2984 sep_line[] =
"\n--------------------------------------------------------------------------------\n";
2997 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3000 <<
"\nComputing: mu_D_hat_calc = -Q_XD_hat'*g - Q_XD_hat'*G*x\n"
3001 <<
" - Q_XD_hat'*A*v(n_R+1:n_R+m) - Q_XD_hat'*A_bar*P_plus_hat*z_hat ...\n";
3003 Workspace<value_type> mu_D_hat_calc_ws( wss, act_set->
q_D_hat() );
3004 DVectorSlice mu_D_hat_calc( &mu_D_hat_calc_ws[0], mu_D_hat_calc_ws.size() );
3005 calc_mu_D( *act_set, *x, *v, &mu_D_hat_calc );
3006 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3008 <<
"\n||mu_D_hat_calc||inf = " << norm_inf(mu_D_hat_calc) << std::endl;
3010 if( (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
3012 <<
"\nmu_D_hat_calc =\n" << mu_D_hat_calc;
3014 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3016 <<
"\nChecking mu_D_hat_calc == mu_D_hat\n";
3018 DVector mu_D_hat_diff(mu_D_hat_calc.dim());
3021 mu_D_hat_err = norm_inf(mu_D_hat_diff()) / (1.0 + norm_inf(mu_D_hat_calc()));
3022 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3024 <<
"\n||mu_D_hat_calc-mu_D_hat||inf/(1.0+||mu_D_hat_calc||inf) = "
3025 << mu_D_hat_err << std::endl;
3029 ,
"QPSchur::qp_algo(...) : Error, "
3030 "||mu_D_hat_calc-mu_D_hat||inf/(1.0+||mu_D_hat_calc||inf) = "
3031 << mu_D_hat_err <<
" >= error_tol = " << error_tol()
3033 if( mu_D_hat_err >= warning_tol() && (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
3035 <<
"\nWarning! ||mu_D_hat_calc-mu_D_hat||inf/(1.0+||mu_D_hat_calc||inf) = "
3036 << mu_D_hat_err <<
" >= warning_tol = " << warning_tol() << std::endl;
3039 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3045 , &b_a, &norm_2_constr, &bnd_ja, &can_ignore_ja );
3046 assume_lin_dep_ja =
false;
3048 return_to_init_fixed =
true;
3050 return_to_init_fixed =
false;
3052 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3054 <<
"\nja = " << ja << endl;
3057 <<
"\ncon_ja_val = " << con_ja_val
3058 <<
"\nb_a = " << b_a
3059 <<
"\nnorm_2_constr = " << norm_2_constr
3060 <<
"\nbnd_ja = " << bnd_str(bnd_ja)
3061 <<
"\ncan_ignore_ja = " << bool_str(can_ignore_ja)
3062 <<
"\nreturn_to_init_fixed = " << bool_str(return_to_init_fixed)
3069 if( (
int)output_level == (
int)OUTPUT_ITER_SUMMARY ) {
3070 if( summary_lines_counter <= 0 ) {
3071 summary_lines_counter = summary_lines_counter_max;
3074 << right << setw(6) <<
"itr"
3075 << right << setw(6) <<
"qhat"
3076 << right << setw(6) <<
"q(+)"
3077 << right << setw(6) <<
"q_F"
3078 << right << setw(6) <<
"q_C"
3079 << right << setw(6) <<
"q_D"
3080 << right << setw(8) <<
"change"
3081 << right << setw(9) <<
"type"
3082 << right << setw(6) <<
"j"
3083 << right << setw(10) <<
"bnd"
3084 << right << setw(dbl_w) <<
"viol, p_z(jd)"
3085 << right << setw(6) <<
"rank" << endl
3086 << right << setw(6) <<
"----"
3087 << right << setw(6) <<
"----"
3088 << right << setw(6) <<
"----"
3089 << right << setw(6) <<
"----"
3090 << right << setw(6) <<
"----"
3091 << right << setw(6) <<
"----"
3092 << right << setw(8) <<
"------"
3093 << right << setw(9) <<
"-------"
3094 << right << setw(6) <<
"----"
3095 << right << setw(10) <<
"--------"
3096 << right << setw(dbl_w) <<
"--------------"
3097 << right << setw(6) <<
"----" << endl;
3101 if( (
int)output_level == (
int)OUTPUT_ITER_SUMMARY ) {
3103 << right << setw(6) << itr
3104 << right << setw(6) << act_set->
q_hat()
3106 << right << setw(6) << act_set->
q_F_hat()
3107 << right << setw(6) << act_set->
q_C_hat()
3108 << right << setw(6) << act_set->
q_D_hat()
3109 << right << setw(8) << ( ja ?
"ADD" :
"-" )
3110 << right << setw(9);
3115 if( bnd_ja == qp.
x_init()(ja) )
3120 else if( ja <= n ) {
3127 << right << setw(6) << ja
3128 << right << setw(10) << ( ja ? bnd_str(bnd_ja) :
"-" )
3129 << right << setw(dbl_w);
3131 *out << (con_ja_val - b_a);
3135 *out << right << setw(6) <<
"-" << endl;
3137 bool found_solution =
false;
3138 const size_type sa = act_set->
s_map(ja);
3139 value_type scaled_viol = 0.0;
3141 found_solution =
true;
3144 const bool is_most_violated
3146 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3148 <<
"\n\nWarning, we have picked the constriant a(" << ja <<
") with violation\n"
3149 <<
"(a(ja)'*x - b_a) = (" << con_ja_val <<
" - " << b_a <<
") = " << (con_ja_val - b_a)
3150 <<
"\nto add to the active set but it is already part of the active set.\n";
3152 const EBounds act_bnd = ( sa != 0 ? act_set->
bnd(sa) : qp.
x_init()(ja) );
3153 if( act_bnd != bnd_ja ) {
3154 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3156 <<
"However, this is not the same bound as the active bound!\n";
3160 if( act_bnd ==
LOWER && act_b_a > b_a || act_bnd ==
UPPER && act_b_a < b_a ) {
3161 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3163 <<
"\nError, c_L_bar(" << ja <<
") = " << (act_b_a > b_a ? act_b_a : b_a)
3164 <<
" > c_U_bar(" << ja <<
") = " << (act_b_a < b_a ? act_b_a : b_a) <<
".\n";
3170 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3172 <<
"The constraints are infeasible! Terminating the QP algorithm!\n";
3174 return INFEASIBLE_CONSTRAINTS;
3177 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3179 <<
"This is the active bound so this is an indication of instability\n"
3180 <<
"in the calculations.\n";
3183 summary_lines_counter = 0;
3184 if( !is_most_violated ) {
3185 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3187 <<
"\nThis is not the most violated constraint so set\n"
3188 <<
"pick_violated_policy = MOST_VIOLATED ...\n";
3190 summary_lines_counter = 0;
3192 next_step = PICK_VIOLATED_CONSTRAINT;
3195 else if( !using_iter_refinement ) {
3196 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3198 <<
"\nThis is the most violated constraint but we are currently not using\n"
3199 <<
"iterative refinement so let's switch it on ...\n";
3201 summary_lines_counter = 0;
3203 *act_set, out, output_level, -1.0, &qp.
fo(), -1.0, act_set->
q_hat() ? &act_set->
d_hat() : NULL
3204 ,v, act_set->
q_hat() ? &act_set->
z_hat() : NULL
3205 ,iter_refine_num_resid, iter_refine_num_solves
3207 if( status == ITER_REFINE_IMPROVED || status == ITER_REFINE_CONVERGED ) {
3208 using_iter_refinement =
true;
3209 next_step = PICK_VIOLATED_CONSTRAINT;
3213 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3215 <<
"\nError, iterative refinement was not allowed or failed to improve the residuals.\n"
3216 <<
"Terminating the QP algorithm ...\n";
3218 return SUBOPTIMAL_POINT;
3222 scaled_viol = std::fabs(con_ja_val - b_a)/(1.0 + std::fabs(con_ja_val));
3223 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3225 <<
"\n\nThis is the most violated constraint, we are using iterative refinement and\n"
3226 <<
"|a("<<ja<<
")'*x - b_a| / (1 + |a("<<ja<<
")'*x|) = |"<<con_ja_val<<
" - "<<b_a
3227 <<
"| / (1 + |"<<con_ja_val<<
"|) = "<<scaled_viol<< (scaled_viol<loose_feas_tol()?
" < ":
" > ")
3228 <<
"loose_feas_tol = "<<loose_feas_tol();
3230 if( scaled_viol < loose_feas_tol() ) {
3231 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3233 <<
"\nTerminating the algorithm with a near optimal solution!\n";
3235 found_solution =
true;
3238 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3240 <<
"\nError! The QP algorithm is terminated!\n";
3242 return SUBOPTIMAL_POINT;
3247 && (scaled_viol = std::fabs(con_ja_val - b_a)/(1.0 + std::fabs(con_ja_val))) < feas_tol() )
3249 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3251 <<
"\nWith all the dof used up, the violated inequality constriant "
3252 <<
"a("<< ja <<
")'*x statisfies:\n"
3253 <<
"|a("<<ja<<
")'*x - b_a| / (1 + |a("<<ja<<
")'*x|) = |" << con_ja_val <<
" - " << b_a
3254 <<
"| / (1 + |"<<con_ja_val<<
"|) = "<<scaled_viol<<
" < feas_tol = "<<feas_tol();
3257 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3259 <<
"\nThis is the most violated constraint so we will consider this\n"
3260 <<
"a near degenerate constraint so we are done!\n";
3262 found_solution =
true;
3265 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3267 <<
"\nThis is not the most violated constraint so set pick_violated_policy = "
3268 <<
"MOST_VIOLATED and pick another violated constraint ....\n";
3271 next_step = PICK_VIOLATED_CONSTRAINT;
3275 if(found_solution) {
3280 if( iter_refine_at_solution() && !using_iter_refinement ) {
3283 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3285 <<
"\nWe think we have found the solution and are not currently using iterative refinement\n"
3286 <<
"and iter_refine_at_solution==true so perform iterative refinement ...\n";
3288 using_iter_refinement =
true;
3290 *act_set, out, output_level, -1.0, &qp.
fo(), -1.0, act_set->
q_hat() ? &act_set->
d_hat() : NULL
3291 ,v, act_set->
q_hat() ? &act_set->
z_hat() : NULL
3292 ,iter_refine_num_resid, iter_refine_num_solves
3295 case ITER_REFINE_ONE_STEP:
3296 case ITER_REFINE_IMPROVED:
3297 case ITER_REFINE_CONVERGED:
3298 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3300 <<
"\nIterative refinement may have altered the unknowns so go back and look for another violated constraint ...\n";
3302 summary_lines_counter = 0;
3303 next_step = PICK_VIOLATED_CONSTRAINT;
3305 case ITER_REFINE_NOT_PERFORMED:
3306 case ITER_REFINE_NOT_NEEDED:
3307 case ITER_REFINE_NOT_IMPROVED:
3308 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3310 <<
"\nIterative refinement did not alter the unknowns so exit with this solution...\n";
3312 set_x( *act_set, *v, x );
3318 if( iter_refine_at_solution() || using_iter_refinement ) {
3324 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3326 <<
"\nRecomputing final mu_D_hat at the solution ...\n";
3328 calc_mu_D( *act_set, *x, *v, &act_set->
mu_D_hat() );
3329 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3331 <<
"\n||mu_D_hat||inf = " << norm_inf(act_set->
mu_D_hat()) << std::endl;
3333 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3335 <<
"\nmu_D_hat =\n" << act_set->
mu_D_hat();
3339 return OPTIMAL_SOLUTION;
3342 case UPDATE_ACTIVE_SET: {
3343 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3345 <<
"\n*** UPDATE_ACTIVE_SET\n";
3356 assume_lin_dep_ja =
true;
3357 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3360 <<
"\nAll of the degrees of freedom are used up so "
3361 <<
"the constraint ja must be linearly dependent.\n"
3362 <<
"The augmented KKT system is definitely singular!\n";
3366 <<
"\nThis is an initially fixed variable that was freed and "
3367 <<
"now is being fixed again.\n"
3368 <<
"The augmented KKT system could be singular or nonsingular, "
3369 <<
"we don't know at this point.\n";
3372 <<
"\nThe schur complement for the new KKT system will not be "
3373 <<
"updated yet in case it is singular!\n";
3377 assume_lin_dep_ja =
false;
3379 if(act_set->
add_constraint( ja, bnd_ja,
false, out, output_level,
true ))
3380 summary_lines_counter = 0;
3383 if( (
int)output_level == (int)OUTPUT_ITER_SUMMARY ) {
3384 *out << right << setw(6) <<
"LI" << endl;
3386 --summary_lines_counter;
3389 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3390 *out <<
"\nNew KKT system is nonsingular! (linearly independent (LI) constraints)\n";
3392 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
3394 <<
"\nPrinting active set after adding constraint ja = " << ja
3399 catch(
const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
3401 if( (
int)output_level >= (
int)OUTPUT_ITER_SUMMARY ) {
3403 <<
"\n\nSchur complement update failed:\n"
3404 <<
"(" << excpt.what() <<
")\n"
3405 <<
"\nConstraint ja = " << ja <<
" appears to be linearly dependent!\n\n";
3407 summary_lines_counter = 0;
3410 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3412 <<
"\nQPSchur::qp_algo(...) : "
3413 <<
"Error, constraint j = "<< ja <<
" is linearly dependent "
3414 <<
"and there are no other constraints to drop.\n"
3415 <<
"The QP must be infeasible\n";
3417 return INFEASIBLE_CONSTRAINTS;
3419 assume_lin_dep_ja =
true;
3420 schur_comp_update_failed =
true;
3422 catch(
const MatrixSymAddDelUpdateable::WrongInertiaUpdateException& excpt ) {
3424 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3426 <<
"\n\nSchur complement appears to have the wrong inertia:\n"
3427 <<
"(" << excpt.what() <<
")\n"
3428 <<
"\nThe QP appears to be nonconvex.\n"
3429 <<
"\nWe have no choice but to terminate the primal-dual QP algorithm!\n";
3431 return NONCONVEX_QP;
3434 if( assume_lin_dep_ja && can_ignore_ja ) {
3436 next_step = PICK_VIOLATED_CONSTRAINT;
3440 beta = ( con_ja_val > b_a ? +1.0 : -1.0 );
3441 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3443 <<
"\nbeta = " << beta << endl;
3446 case COMPUTE_SEARCH_DIRECTION: {
3447 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3449 <<
"\n*** COMPUTE_SEARCH_DIRECTION\n";
3451 const EtaVector e_ja( ja, n + m_breve );
3452 if( assume_lin_dep_ja ) {
3453 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3455 <<
"\nThe KKT system for the trial active set is assumed or known to be singular!\n"
3456 <<
"\nComputing the steps from:"
3457 <<
"\np_z_hat = inv(S_hat)*(-v_a+U_hat'*inv(Ko)*u_a), p_v = inv(Ko)*(-u_a-U_hat*p_z_hat) ...\n";
3493 sa = act_set->
s_map(-
int(ja)),
3498 Workspace<value_type> v_a_ws(wss,act_set->
q_hat());
3499 DVectorSlice v_a(&v_a_ws[0],v_a_ws.size());
3504 d_a = ( bnd_ja == qp.
x_init()(ja)
3506 : b_a - qp.
b_X()(la) );
3511 if(!all_dof_used_up) {
3512 calc_v( qp.
Ko(), NULL, act_set->
U_hat(), act_set->
p_z_hat(), &p_v() );
3518 if( using_iter_refinement ) {
3519 if( (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3521 <<
"\n\nPerforming iterative refinement on p_v, p_z_hat system ...\n";
3523 summary_lines_counter = 0;
3527 *act_set, out, output_level, 0.0, NULL, +1.0, &v_a
3529 ,iter_refine_num_resid, iter_refine_num_solves
3533 if(!all_dof_used_up)
3534 gamma_plus = ( ( d_a -
dot(v_a(),act_set->
z_hat()) )
3537 gamma_plus = beta * inf;
3566 la = act_set->
qp().
Q_R().lookup_col_j(ja);
3568 const EtaVector u_a = EtaVector(la,n_R+m);
3569 const value_type d_a = b_a;
3573 if( act_set->
q_hat() ) {
3580 V_StV( &t1, -1.0, u_a() );
3584 if(!all_dof_used_up)
3594 if( using_iter_refinement ) {
3595 if( (
int)output_level >= (
int)OUTPUT_ITER_SUMMARY ) {
3597 <<
"\n\nPerforming iterative refinement on p_v, p_z_hat system ...\n";
3599 summary_lines_counter = 0;
3602 Workspace<value_type> dense_u_a_ws(wss,u_a().dim());
3603 DVectorSlice dense_u_a(&dense_u_a_ws[0],dense_u_a_ws.size());
3605 dense_u_a(u_a().begin()->index()+u_a().offset()) = 1.0;
3607 *act_set, out, output_level, +1.0, &dense_u_a, 0.0, NULL
3609 ,iter_refine_num_resid, iter_refine_num_solves
3611 summary_lines_counter = 0;
3614 if(!all_dof_used_up)
3615 gamma_plus = ( d_a -
dot(u_a(),*v) ) /
dot(u_a(),p_v());
3617 gamma_plus = beta * inf;
3629 Workspace<value_type> u_a_ws( wss, n_R + m );
3630 DVectorSlice u_a( &u_a_ws[0], u_a_ws.size() );
3631 Workspace<value_type> v_a_ws( wss, act_set->
q_hat() );
3632 DVectorSlice v_a( &v_a_ws[0], v_a_ws.size() );
3638 u_a(n_R+1,n_R+m) = 0.0;
3640 Workspace<value_type> t0_ws( wss, n-n_R );
3641 DVectorSlice t0( &t0_ws[0], t0_ws.size() );
3647 d_a = b_a - ( n > n_R ?
dot( qp.
b_X(), t0() ) : 0.0 );
3649 Workspace<value_type> t1_ws( wss, n_R + m );
3650 DVectorSlice t1( &t1_ws[0], t1_ws.size() );
3652 if( act_set->
q_hat() ) {
3654 Workspace<value_type> t2_ws( wss, act_set->
q_hat() );
3655 DVectorSlice t2( &t2_ws[0], t2_ws.size() );
3661 Vp_StV( &t2(), -1.0, v_a() );
3664 if(!all_dof_used_up) {
3666 V_StV( &t1, -1.0, u_a() );
3676 if( using_iter_refinement ) {
3677 if( (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3679 <<
"\n\nPerforming iterative refinement on p_v, p_z_hat system ...\n";
3681 summary_lines_counter = 0;
3685 *act_set, out, output_level, +1.0, &u_a, +1.0, act_set->
q_hat() ? &v_a : NULL
3687 ,iter_refine_num_resid, iter_refine_num_solves
3689 summary_lines_counter = 0;
3692 if(!all_dof_used_up)
3693 gamma_plus = ( ( d_a -
dot(u_a,*v) -
dot(v_a(),act_set->
z_hat()) )
3696 gamma_plus = beta * inf;
3700 if(!all_dof_used_up)
3705 if(!all_dof_used_up)
3706 gamma_plus = ( d_a -
dot(u_a,*v) ) /
dot(u_a,p_v());
3708 gamma_plus = beta * inf;
3712 if( schur_comp_update_failed && gamma_plus * beta < 0 ) {
3713 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3715 <<
"\nThe schur complement update failed and gamma_plus = " << gamma_plus <<
" is the wrong sign"
3716 <<
"\nso we will assume the sign error for (...)/+-0 was due to arbitrary roundoff"
3717 <<
"\nand therefore we will set gamma_plus = -gamma_plus\n";
3719 gamma_plus = -gamma_plus;
3722 if(act_set->
q_hat()) {
3723 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3725 <<
"\n||p_z_hat||inf = " << norm_inf(act_set->
p_z_hat()) << endl;
3727 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3728 *out <<
"\np_z_hat =\n" << act_set->
p_z_hat();
3731 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3733 <<
"\n||p_v||inf = " << norm_inf(p_v()) << endl;
3735 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
3736 *out <<
"\np_v =\n" << p_v();
3738 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3740 <<
"\ngamma_plus = " << gamma_plus << endl;
3745 calc_p_mu_D( *act_set, p_v(), act_set->
p_z_hat(), &ja, &act_set->
p_mu_D_hat() );
3746 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3748 <<
"\n||p_mu_D_hat||inf = " << norm_inf(act_set->
p_mu_D_hat()) << std::endl;
3750 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3752 <<
"\np_mu_D_hat =\n" << act_set->
p_mu_D_hat();
3759 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3761 <<
"\nThe KKT system for the new active set is or known to be nonsingular!\n"
3762 <<
"\nComputing the steps from:"
3763 <<
"\nz_hat_plus = inv(S_hat)*(d_hat-U_hat'vo), v_plus = inv(Ko)*(fo-U_hat*z_hat_plus) ...\n";
3769 const size_type q_hat = act_set->
q_hat();
3770 z_hat_plus.bind(DVectorSlice(&z_hat_plus_ws[0],q_hat));
3773 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3775 <<
"\n||z_hat_plus||inf = " << norm_inf(z_hat_plus()) << std::endl;
3777 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
3779 <<
"\nz_hat_plus =\n" << z_hat_plus();
3782 calc_v( qp.
Ko(), &qp.
fo(), act_set->
U_hat(), z_hat_plus
3784 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3786 <<
"\n||v_plus||inf = " << norm_inf(v_plus()) << std::endl;
3788 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
3790 <<
"\nv_plus =\n" << v_plus();
3792 if( using_iter_refinement ) {
3793 if( (
int)output_level >= (
int)OUTPUT_ITER_SUMMARY ) {
3795 <<
"\n\nPerforming iterative refinement on v_plus, z_hat_plus system ...\n";
3797 summary_lines_counter = 0;
3801 *act_set, out, output_level, -1.0, &qp.
fo(), -1.0, &act_set->
d_hat()
3802 ,&v_plus, &z_hat_plus
3803 ,iter_refine_num_resid, iter_refine_num_solves
3807 DVectorSlice p_z_hat = act_set->
p_z_hat();
3809 V_VmV( &p_z_hat(), z_hat_plus(), act_set->
z_hat() );
3811 V_VmV( &p_v(), v_plus(), *v );
3814 calc_p_mu_D( *act_set, p_v(), p_z_hat(), NULL, &act_set->
p_mu_D_hat() );
3816 const size_type sa = act_set->
s_map(ja);
3824 gamma_plus = z_hat_plus(sa);
3834 Vt_S( &p_z_hat(), 1.0 / gamma_plus );
3836 Vt_S( &p_v(), 1.0 / gamma_plus );
3838 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3840 <<
"\ngamma_plus = " << gamma_plus << std::endl;
3842 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3844 <<
"\n||p_z_hat||inf = " << norm_inf(p_z_hat()) << std::endl;
3846 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
3848 <<
"\np_z_hat =\n" << p_z_hat();
3850 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3852 <<
"\n||p_v||inf = " << norm_inf(p_v()) << std::endl;
3854 if( (
int)output_level >= (
int)OUTPUT_ITER_QUANTITIES ) {
3856 <<
"\np_v =\n" << p_v();
3861 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3863 <<
"\n||p_mu_D_hat||inf =\n" << norm_inf(act_set->
p_mu_D_hat()) << std::endl;
3865 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3867 <<
"\np_mu_D_hat =\n" << act_set->
p_mu_D_hat();
3872 case COMPUTE_STEP_LENGTHS: {
3874 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3876 <<
"\n*** COMPUTE_STEP_LENGTHS\n";
3879 const size_type q_hat = act_set->
q_hat();
3880 dual_infeas_scale = 1.0;
3891 t_P = beta * gamma_plus;
3893 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3895 <<
"\nWarning, A near degenerate inequality constraint ja = " << ja
3896 <<
" is being added that has the wrong sign with:\n"
3897 <<
" t_P = " << t_P << std::endl
3898 <<
" dual_infeas_scale = " << dual_infeas_scale << std::endl
3899 <<
" norm_2_constr = " << norm_2_constr << std::endl
3900 <<
" |t_P/(norm_2_constr*dual_infeas_scale)| = "
3901 << std::fabs(t_P/(norm_2_constr*dual_infeas_scale))
3902 <<
" <= dual_infeas_tol = " << dual_infeas_tol() << std::endl;
3904 summary_lines_counter = 0;
3905 if( !using_iter_refinement ) {
3906 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3907 *out <<
"We are not using iterative refinement yet so turn it on"
3908 <<
"\nthen recompute the steps ...\n";
3910 using_iter_refinement =
true;
3911 next_step = COMPUTE_SEARCH_DIRECTION;
3915 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
3917 <<
"We are already using iterative refinement so the QP algorithm is terminated!\n";
3919 return DUAL_INFEASIBILITY;
3922 t_P = beta * gamma_plus;
3923 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
3925 <<
"\nt_P = " << t_P << endl;
3935 value_type max_feas_viol = 0.0;
3941 if( act_set->
q_hat() ) {
3942 DVectorSlice z_hat = act_set->
z_hat();
3943 DVectorSlice p_z_hat = act_set->
p_z_hat();
3944 DVectorSlice::iterator
3945 z_itr = z_hat.begin(),
3946 p_z_itr = p_z_hat.begin();
3948 qq = assume_lin_dep_ja || (!assume_lin_dep_ja && return_to_init_fixed)
3949 ? q_hat : q_hat - 1;
3951 if( qq > 0 && (
int)output_level >= (int)OUTPUT_ACT_SET ) {
3953 <<
"\nComputing the maximum step for multiplers for dual feasibility\n\n"
3954 << right << setw(5) <<
"s"
3955 << right << setw(5) <<
"j"
3956 << right << setw(dbl_w) <<
"z_hat"
3957 << right << setw(dbl_w) <<
"p_z_hat"
3958 << right << setw(20) <<
"bnd"
3959 << right << setw(dbl_w) <<
"t"
3960 << right << setw(dbl_w) <<
"t_D"
3961 << right << setw(5) <<
"jd" << endl
3962 << right << setw(5) <<
"----"
3963 << right << setw(5) <<
"----"
3964 << right << setw(dbl_w) <<
"--------------"
3965 << right << setw(dbl_w) <<
"--------------"
3966 << right << setw(20) <<
"--------------"
3967 << right << setw(dbl_w) <<
"--------------"
3968 << right << setw(dbl_w) <<
"--------------"
3969 << right << setw(5) <<
"----" << endl;
3971 for(
int s = 1; s <= qq; ++s, ++z_itr, ++p_z_itr) {
3972 int j = act_set->
ij_map(s);
3974 namespace ns = QPSchurPack;
3975 EBounds bnd = act_set->
bnd(s);
3977 if( (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
3979 << right << setw(5) << s
3980 << right << setw(5) << j
3981 << right << setw(dbl_w) << *z_itr
3982 << right << setw(dbl_w) << *p_z_itr
3983 << right << setw(20) << bnd_str(bnd);
3987 bool j_is_degen =
false;
3989 const int dual_feas_status
3990 = correct_dual_infeas(
3992 ,out,output_level,
true,
"z_hat(s)",&(*z_itr),&viol
3993 ,
"p_z_hat(s)",&(*p_z_itr),
"z_hat_plus(s)"
3994 , (assume_lin_dep_ja ? NULL: &z_hat_plus(s) ) );
3995 if( dual_feas_status < 0 ) {
3996 if( !using_iter_refinement ) {
3997 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO )
3998 *out <<
"We are not using iterative refinement yet so turn it on"
3999 <<
"\nthen recompute the steps ...\n";
4000 using_iter_refinement =
true;
4001 next_step = COMPUTE_SEARCH_DIRECTION;
4005 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO )
4006 *out <<
"We are already using iterative refinement so the QP algorithm is terminated!\n";
4007 return DUAL_INFEASIBILITY;
4010 else if( dual_feas_status == 0 ) {
4015 const value_type feas_viol = beta*(*p_z_itr);
4016 if( bnd == EQUALITY )
4018 else if( bnd ==
LOWER && feas_viol <= 0.0 )
4020 else if( bnd ==
UPPER && feas_viol >= 0.0 )
4024 t = -beta*(*z_itr)/(*p_z_itr);
4028 if(j_is_degen) j_degen = j;
4029 max_feas_viol = feas_viol;
4034 if( (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
4036 << right << setw(dbl_w) << t
4037 << right << setw(dbl_w) << t_D
4038 << right << setw(5) << jd << endl;
4048 const size_type q_D_hat = act_set->
q_D_hat();
4049 DVectorSlice mu_D_hat = act_set->
mu_D_hat();
4050 DVectorSlice p_mu_D_hat = act_set->
p_mu_D_hat();
4052 qD = assume_lin_dep_ja && return_to_init_fixed ? q_D_hat-1 : q_D_hat;
4054 if( qD > 0 && (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
4056 <<
"\nComputing the maximum step for multiplers for dual feasibility\n\n"
4057 << right << setw(5) <<
"k"
4058 << right << setw(5) <<
"i"
4059 << right << setw(dbl_w) <<
"mu_D_hat"
4060 << right << setw(dbl_w) <<
"p_mu_D_hat"
4061 << right << setw(20) <<
"x_init"
4062 << right << setw(dbl_w) <<
"t"
4063 << right << setw(dbl_w) <<
"t_D"
4064 << right << setw(5) <<
"jd" << endl
4065 << right << setw(5) <<
"----"
4066 << right << setw(5) <<
"----"
4067 << right << setw(dbl_w) <<
"--------------"
4068 << right << setw(dbl_w) <<
"--------------"
4069 << right << setw(20) <<
"--------------"
4070 << right << setw(dbl_w) <<
"--------------"
4071 << right << setw(dbl_w) <<
"--------------"
4072 << right << setw(5) <<
"----" << endl;
4074 GenPermMatrixSlice::const_iterator
4075 Q_XD_itr = act_set->
Q_XD_hat().begin(),
4076 Q_XD_end = Q_XD_itr + qD;
4077 for( ; Q_XD_itr != Q_XD_end; ++Q_XD_itr ) {
4078 const size_type k = Q_XD_itr->col_j();
4079 const size_type i = Q_XD_itr->row_i();
4080 DVectorSlice::iterator
4081 mu_D_itr = mu_D_hat.begin() + (k-1),
4082 p_mu_D_itr = p_mu_D_hat.begin() + (k-1);
4083 const size_type l = act_set->
l_fxfx(k);
4084 EBounds bnd = qp.
x_init()(i);
4086 if( (
int)output_level >= (int)OUTPUT_ACT_SET ) {
4088 << right << setw(5) << k
4089 << right << setw(5) << i
4090 << right << setw(dbl_w) << *mu_D_itr
4091 << right << setw(dbl_w) << *p_mu_D_itr
4092 << right << setw(20) << bnd_str(bnd);
4096 bool j_is_degen =
false;
4098 const int dual_feas_status
4099 = correct_dual_infeas(
4101 ,out,output_level,
true,
"mu_D_hat(k)",&(*mu_D_itr),&viol
4102 ,
"p_mu_D_hat(k)",&(*p_mu_D_itr) );
4103 if( dual_feas_status < 0 ) {
4104 if( !using_iter_refinement ) {
4105 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO )
4106 *out <<
"We are not using iterative refinement yet so turn it on"
4107 <<
"\nthen recompute the steps ...\n";
4108 using_iter_refinement =
true;
4109 next_step = COMPUTE_SEARCH_DIRECTION;
4113 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO )
4114 *out <<
"We are already using iterative refinement so the QP algorithm is terminated!\n";
4115 return DUAL_INFEASIBILITY;
4118 else if( dual_feas_status == 0 ) {
4123 const value_type feas_viol = beta*(*p_mu_D_itr);
4124 if( bnd == EQUALITY )
4126 else if( bnd ==
LOWER && feas_viol <= 0.0 )
4128 else if( bnd ==
UPPER && feas_viol >= 0.0 )
4132 t = -beta*(*mu_D_itr)/(*p_mu_D_itr);
4136 if(j_is_degen) j_degen = jd;
4137 max_feas_viol = feas_viol;
4142 if( (
int)output_level >= (
int)OUTPUT_ACT_SET ) {
4144 << right << setw(dbl_w) << t
4145 << right << setw(dbl_w) << t_D
4146 << right << setw(5) << jd << endl;
4151 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
4153 <<
"\nt_D = " << t_D << endl
4154 <<
"jd = " << jd << endl;
4156 if( jd == j_degen && jd != 0 && t_D < t_P ) {
4157 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4159 <<
"\nWarning, the near degenerate constraint j = "
4160 << jd <<
" which had the incorrect sign\nand was adjusted "
4161 <<
"was selected to be dropped from the active set.\n";
4165 if( assume_lin_dep_ja && !schur_comp_update_failed && (
int)output_level == (
int)OUTPUT_ITER_SUMMARY ) {
4166 if( t_P < huge_primal_step() )
4167 *out << right << setw(6) <<
"LI" << endl;
4169 *out << right << setw(6) <<
"LD" << endl;
4171 --summary_lines_counter;
4174 if( t_D < t_P && (
int)output_level == (
int)OUTPUT_ITER_SUMMARY ) {
4176 << right << setw(6) << itr
4177 << right << setw(6) << act_set->
q_hat()
4179 << right << setw(6) << act_set->
q_F_hat()
4180 << right << setw(6) << act_set->
q_C_hat()
4181 << right << setw(6) << act_set->
q_D_hat()
4182 << right << setw(8) <<
"DROP"
4183 << right << setw(9);
4187 else if( jd <= n ) {
4188 if( bnd_jd == qp.
x_init()(jd) )
4197 << right << setw(6) << jd
4198 << right << setw(10) << bnd_str(bnd_jd)
4199 << right << setw(dbl_w) << max_feas_viol
4200 << right << setw(6) <<
"LI" << endl;
4204 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
4206 <<
"\n*** TAKE_STEP\n";
4208 if( t_P >= huge_primal_step() && t_D >= huge_dual_step() ) {
4209 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4211 <<
"Error, QP is infeasible, inconsistent constraint a("<<ja<<
") detected\n";
4213 if( using_iter_refinement ) {
4214 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4216 <<
"We are already using iterative refinement so the QP algorithm is terminated!\n";
4218 return INFEASIBLE_CONSTRAINTS;
4221 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4222 *out <<
"We are not using iterative refinement yet so turn it on";
4223 if(assume_lin_dep_ja)
4224 *out <<
"\nthen pick another violated constriant to add ... \n";
4226 *out <<
"\nthen recompute the steps ...\n";
4228 summary_lines_counter = 0;
4231 using_iter_refinement =
true;
4232 if(assume_lin_dep_ja) {
4234 *act_set, out, output_level, -1.0, &qp.
fo(), -1.0, act_set->
q_hat() ? &act_set->
d_hat() : NULL
4235 ,v, act_set->
q_hat() ? &act_set->
z_hat() : NULL
4236 ,iter_refine_num_resid, iter_refine_num_solves
4238 next_step = PICK_VIOLATED_CONSTRAINT;
4242 next_step = COMPUTE_SEARCH_DIRECTION;
4247 else if( t_P > t_D ) {
4248 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS ) {
4249 if( t_P >= huge_primal_step() ) {
4251 <<
"\n*** (b) Dual Step (t_P = " << t_P <<
" >= huge_primal_step = "
4252 << huge_primal_step() <<
")\n";
4256 <<
"\n*** (b) Partial Primal-Dual Step\n";
4260 if( ja == last_jd && jd == last_ja ) {
4261 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4263 <<
"\n\nQPSchur::qp_algo(...) : Error, the constraint "
4264 <<
"a(" << ja <<
") with violation\n"
4265 <<
"(a(ja)'*x - b_a) = (" << con_ja_val
4266 <<
" - " << b_a <<
") = " << (con_ja_val - b_a) <<
"\n"
4267 <<
"we are adding to the active set and the constraint constriant\n"
4268 <<
"a(" << jd <<
") we are dropping were just dropped and added respectively!\n"
4269 <<
"The algorithm is cycling!\n";
4271 if( using_iter_refinement ) {
4272 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4274 <<
"We are already using iterative refinement so the QP algorithm is terminated!\n";
4276 return SUBOPTIMAL_POINT;
4279 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4280 *out <<
"We are not using iterative refinement yet so turn it on";
4281 if(assume_lin_dep_ja)
4282 *out <<
"\nthen pick another violated constriant to add ... \n";
4284 *out <<
"\nthen recompute the steps ...\n";
4286 summary_lines_counter = 0;
4289 using_iter_refinement =
true;
4290 if(assume_lin_dep_ja) {
4292 *act_set, out, output_level, -1.0, &qp.
fo(), -1.0, act_set->
q_hat() ? &act_set->
d_hat() : NULL
4293 ,v, act_set->
q_hat() ? &act_set->
z_hat() : NULL
4294 ,iter_refine_num_resid, iter_refine_num_solves
4296 next_step = PICK_VIOLATED_CONSTRAINT;
4300 next_step = COMPUTE_SEARCH_DIRECTION;
4307 if( assume_lin_dep_ja ) {
4309 summary_lines_counter = 0;
4313 summary_lines_counter = 0;
4316 catch(
const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
4317 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4319 <<
"\n\nSchur complement appears to be singular and should not be:\n"
4321 <<
"\nThe QP appears to be nonconvex and we therefore terminate the primal-dual QP algorithm!\n";
4323 return NONCONVEX_QP;
4325 catch(
const MatrixSymAddDelUpdateable::WrongInertiaUpdateException& excpt ) {
4326 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4328 <<
"\n\nSchur complement appears to have the wrong inertia:\n"
4330 <<
"\nThe QP appears to be nonconvex and we therefore terminate the primal-dual QP algorithm!\n";
4332 return NONCONVEX_QP;
4335 if(act_set->
q_hat())
4338 Vp_StV( v, beta * t_D, p_v() );
4345 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS )
4348 <<
"\nUpdated primal and dual variables:\n"
4349 <<
"\n||v||inf = " << norm_inf(*v) << endl;
4350 if(act_set->
q_hat()) {
4352 <<
"||z_hat||inf = " << norm_inf(act_set->
z_hat()) << endl;
4356 <<
"max(|mu_D_hat(i)|) = " << norm_inf(act_set->
mu_D_hat()) << endl
4357 <<
"min(|mu_D_hat(i)|) = " << min_abs(act_set->
mu_D_hat()) << endl;
4360 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES )
4362 *out <<
"\nv = \n" << *v << endl;
4363 if(assume_lin_dep_ja) {
4365 <<
"\nPrinting active set after dropping constraint jd = " << jd
4366 <<
" and adding constraint ja = " << ja <<
" ...\n";
4370 <<
"\nPrinting active set after dropping constraint jd = " << jd
4376 assume_lin_dep_ja =
false;
4377 next_step = COMPUTE_SEARCH_DIRECTION;
4381 if( (
int)output_level >= (
int)OUTPUT_ITER_STEPS ) {
4383 <<
"\n*** (c) Full Primal-Dual Step\n";
4386 if( !assume_lin_dep_ja ) {
4387 act_set->
z_hat() = z_hat_plus;
4391 bool threw_exception =
false;
4393 if(act_set->
add_constraint( ja, bnd_ja,
true, out, output_level,
true,
true ))
4394 summary_lines_counter = 0;
4396 catch(
const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
4397 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4399 <<
"\n\nSchur complement appears to be singular and should not be:\n"
4400 << excpt.what() << std::endl;
4402 threw_exception =
true;
4404 catch(
const MatrixSymAddDelUpdateable::WrongInertiaUpdateException& excpt ) {
4405 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4407 <<
"\n\nSchur complement appears to have the wrong inertia:\n"
4408 << excpt.what() << std::endl;
4410 threw_exception =
true;
4412 if( threw_exception ) {
4413 if( !using_iter_refinement ) {
4414 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4415 *out <<
"We are not using iterative refinement yet so turn it on and\n"
4416 <<
"go back and pick a new violated constraint to add to the active set ...\n";
4418 using_iter_refinement =
true;
4419 if( (
int)output_level >= (
int)OUTPUT_ITER_SUMMARY ) {
4421 <<
"\n\nPerforming iterative refinement on v, z_hat system ...\n";
4423 summary_lines_counter = 0;
4427 *act_set, out, output_level, -1.0, &qp.
fo(), -1.0, &act_set->
d_hat()
4428 ,v, &act_set->
z_hat()
4429 ,iter_refine_num_resid, iter_refine_num_solves
4431 next_step = PICK_VIOLATED_CONSTRAINT;
4435 if( (
int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4436 *out <<
"Darn, we are already using iterative refinement!"
4437 <<
"\nThe QP appears to be nonconvex and we therefore terminate the primal-dual QP algorithm!\n";
4439 return NONCONVEX_QP;
4443 if(act_set->
q_hat())
4446 Vp_StV( v, beta * t_P, p_v() );
4453 if( (
int)output_level >= (int)OUTPUT_ITER_STEPS )
4456 <<
"\nUpdated primal and dual variables:\n"
4457 <<
"\n||v||inf = " << norm_inf(*v) << endl;
4458 if(act_set->
q_hat()) {
4460 <<
"||z_hat||inf = " << norm_inf(act_set->
z_hat()) << endl;
4464 <<
"max(|mu_D_hat(i)|) = " << norm_inf(act_set->
mu_D_hat()) << endl
4465 <<
"min(|mu_D_hat(i)|) = " << min_abs(act_set->
mu_D_hat()) << endl;
4468 if( (
int)output_level >= (int)OUTPUT_ITER_QUANTITIES )
4470 *out <<
"\nv = \n" << *v << endl;
4471 if( assume_lin_dep_ja ) {
4473 <<
"\nPrinting active set after adding constraint ja = " << ja
4478 if(act_set->
q_hat())
4479 *out <<
"\nz_hat =\n" << act_set->
z_hat();
4481 *out <<
"\nmu_D_hat =\n" << act_set->
mu_D_hat();
4484 assume_lin_dep_ja =
false;
4485 next_step = PICK_VIOLATED_CONSTRAINT;
4495 catch( std::exception& excpt ) {
4496 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4498 <<
"\n\n*** Caught a standard exception :\n"
4500 <<
"\n*** Rethrowing the exception ...\n";
4505 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4507 <<
"\n\n*** Caught an unknown exception. Rethrowing the exception ...\n";
4512 return MAX_ITER_EXCEEDED;
4518 using LinAlgOpPack::V_MtV;
4519 using LinAlgOpPack::Vp_MtV;
4523 if( act_set.
qp().
n() > act_set.
qp().
n_R() )
4530 const ActiveSet& act_set,
const DVectorSlice& v
4531 ,SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
4535 using LinAlgOpPack::V_MtV;
4536 using LinAlgOpPack::Vp_MtV;
4538 using LinAlgOpPack::V_MtV;
4539 using LinAlgOpPack::Vp_MtV;
4540 using AbstractLinAlgPack::V_MtV;
4541 using AbstractLinAlgPack::Vp_MtV;
4542 namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
4547 n = act_set.
qp().
n(),
4548 n_R = act_set.
qp().
n_R(),
4549 m = act_set.
qp().
m(),
4551 q_hat = act_set.
q_hat();
4561 typedef SpVector::element_type ele_t;
4562 mu->resize( n, n-m );
4563 lambda_breve->resize( m_breve, n-m );
4570 z_hat = act_set.
z_hat();
4571 for( size_type s = 1; s <= q_hat; ++s ) {
4572 const int ij = act_set.
ij_map(s);
4574 const size_type j = ij;
4576 mu->add_element(ele_t(j,z_hat(s)));
4578 lambda_breve->add_element(ele_t(j-n,z_hat(s)));
4583 lambda_breve->sort();
4586 *lambda = v(n_R+1,n_R+m);
4592 const value_type minutes = timer->
read() / 60;
4593 if( minutes >= max_real_runtime() ) {
4594 if( (
int)output_level >= (
int)OUTPUT_BASIC_INFO ) {
4596 <<
"\n*** Runtime = " << minutes <<
" min >= max_real_runtime = " << max_real_runtime() <<
" min!\n"
4597 <<
"Must terminite the algorithm!\n";
4609 ,
const value_type ao
4610 ,
const DVectorSlice *bo
4611 ,
const value_type aa
4612 ,
const DVectorSlice *ba
4615 ,size_type *iter_refine_num_resid
4616 ,size_type *iter_refine_num_solves
4625 using DenseLinAlgPack::norm_inf;
4633 typedef DenseLinAlgPack::value_type extra_value_type;
4635 const value_type small_num = std::numeric_limits<value_type>::min();
4637 const int int_w = 8;
4638 const char int_ul[] =
"------";
4639 const int dbl_min_w = 20;
4640 const int dbl_w = ( out ? my_max(dbl_min_w,
int(out->precision()+8)): 20 );
4641 const char dbl_ul[] =
"------------------";
4645 const MatrixSymOpNonsing
4647 &S_hat = act_set.
S_hat();
4649 &U_hat = act_set.
U_hat();
4654 q_hat = act_set.
q_hat();
4657 d_hat = (q_hat ? act_set.
d_hat() : DVectorSlice());
4659 Workspace<extra_value_type>
4660 ext_ro_ws(wss,n_R+m),
4661 ext_ra_ws(wss,q_hat);
4663 ext_ro(&ext_ro_ws[0],ext_ro_ws.size()),
4664 ext_ra(ext_ra_ws.size()?&ext_ra_ws[0]:NULL,ext_ra_ws.size());
4665 Workspace<value_type>
4669 del_v_ws(wss,n_R+m),
4670 del_z_ws(wss,q_hat),
4671 v_itr_ws(wss,n_R+m),
4672 z_itr_ws(wss,q_hat);
4674 ro(&ro_ws[0],ro_ws.size()),
4675 ra(ra_ws.size()?&ra_ws[0]:NULL,ra_ws.size()),
4676 t1(&t1_ws[0],t1_ws.size()),
4677 del_v(&del_v_ws[0],del_v_ws.size()),
4678 del_z(del_z_ws.size()?&del_z_ws[0]:NULL,del_z_ws.size()),
4679 v_itr(&v_itr_ws[0],v_itr_ws.size()),
4680 z_itr(z_itr_ws.size()?&z_itr_ws[0]:NULL,z_itr_ws.size());
4688 if( out && (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4690 <<
"\nBeginning iterative refinement ..."
4691 <<
"\niter_refine_opt_tol = " << iter_refine_opt_tol()
4692 <<
", iter_refine_feas_tol = " << iter_refine_feas_tol()
4693 <<
"\niter_refine_min_iter = " << iter_refine_min_iter()
4694 <<
", iter_refine_max_iter = " << iter_refine_max_iter() <<
"\n\n";
4697 << right << setw(int_w) <<
"ir_itr"
4698 << right << setw(dbl_w) <<
"roR_scaling"
4699 << right << setw(dbl_w) <<
"||roR||s"
4700 << left << setw(1) <<
" ";
4703 << right << setw(dbl_w) <<
"rom_scaling"
4704 << right << setw(dbl_w) <<
"||rom||s"
4705 << left << setw(1) <<
" ";
4709 << right << setw(dbl_w) <<
"ra_scaling"
4710 << right << setw(dbl_w) <<
"||ra||s"
4711 << left << setw(1) <<
" ";
4714 << right << setw(dbl_w) <<
"||del_v||/||v||inf"
4715 << right << setw(dbl_w) <<
"||del_z||/||z||inf"
4719 << right << setw(int_w) << int_ul
4720 << right << setw(dbl_w) << dbl_ul
4721 << right << setw(dbl_w) << dbl_ul
4722 << left << setw(1) <<
" ";
4725 << right << setw(dbl_w) << dbl_ul
4726 << right << setw(dbl_w) << dbl_ul
4727 << left << setw(1) <<
" ";
4731 << right << setw(dbl_w) << dbl_ul
4732 << right << setw(dbl_w) << dbl_ul
4733 << left << setw(1) <<
" ";
4736 << right << setw(dbl_w) << dbl_ul
4737 << right << setw(dbl_w) << dbl_ul
4745 roR_nrm_o, rom_nrm_o, ra_nrm_o,
4746 roR_nrm, rom_nrm, ra_nrm;
4747 for( size_type iter_refine_k = 0;
4748 iter_refine_k < iter_refine_max_iter() && return_status != ITER_REFINE_CONVERGED;
4761 ++(*iter_refine_num_resid);
4769 ,m ? &rom_scaling : NULL
4772 ,q_hat ? &ext_ra : NULL
4773 ,q_hat ? &ra_scaling : NULL
4775 std::copy(ext_ro.begin(),ext_ro.end(),ro.
begin());
4776 if(q_hat) std::copy(ext_ra.begin(),ext_ra.end(),ra.
begin());
4780 roR_nrm = norm_inf(ro(1,n_R));
4781 rom_nrm = (m ? norm_inf(ro(n_R+1,n_R+m)) : 0.0);
4782 ra_nrm = (q_hat ? norm_inf(ra) : 0.0);
4783 if( iter_refine_k == 0 ) {
4784 roR_nrm_o = roR_nrm;
4785 rom_nrm_o = rom_nrm;
4789 is_roR_conv = roR_nrm / (1.0 + roR_scaling) < iter_refine_opt_tol(),
4790 is_rom_conv = (m ? rom_nrm / (1.0 + rom_scaling) < iter_refine_feas_tol() : true ),
4791 is_ra_conv = (q_hat ? ra_nrm / (1.0 + ra_scaling) < iter_refine_feas_tol() : true ),
4792 is_conv = is_roR_conv && is_rom_conv && is_ra_conv;
4796 if( out && (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4798 << right << setw(int_w) << iter_refine_k
4799 << right << setw(dbl_w) << roR_scaling
4800 << right << setw(dbl_w) << (roR_nrm / (1.0 + roR_scaling))
4801 << left << setw(1) << (is_roR_conv ?
"*" :
" ");
4804 << right << setw(dbl_w) << rom_scaling
4805 << right << setw(dbl_w) << (rom_nrm /(1.0 + rom_scaling))
4806 << left << setw(1) << (is_rom_conv ?
"*" :
" ");
4810 << right << setw(dbl_w) << ra_scaling
4811 << right << setw(dbl_w) << (ra_nrm /(1.0 + ra_scaling))
4812 << left << setw(1) << (is_ra_conv ?
"*" :
" ");
4818 if( iter_refine_k + 1 < iter_refine_min_iter() ) {
4821 else if( is_conv ) {
4822 if( out && (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4824 << right << setw(dbl_w) <<
"-"
4825 << right << setw(dbl_w) <<
"-"
4828 if( iter_refine_k == 0 )
4829 return_status = ITER_REFINE_NOT_NEEDED;
4831 return_status = ITER_REFINE_CONVERGED;
4835 if( roR_nrm_o < roR_nrm && rom_nrm_o < rom_nrm && ra_nrm_o < rom_nrm ) {
4836 return_status = ITER_REFINE_NOT_IMPROVED;
4845 ++(*iter_refine_num_solves);
4849 calc_z( act_set.
S_hat(), ra, act_set.
U_hat(), &t1, &del_z );
4851 calc_v( Ko, &ro, U_hat, del_z, &del_v );
4858 Vp_StV( &v_itr, -1.0, del_v );
4860 Vp_StV( &z_itr, -1.0, del_z );
4864 if( out && (
int)output_level >= (
int)OUTPUT_ITER_SUMMARY ) {
4866 << right << setw(dbl_w) << norm_inf(del_v) / (norm_inf(v_itr) + small_num)
4867 << right << setw(dbl_w) << norm_inf(del_z) / (norm_inf(z_itr) + small_num)
4871 if( iter_refine_max_iter() == 0 ) {
4872 if( (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4874 <<
"\nWarning, iter_refine_max_iter == 0. Iterative refinement was not performed."
4875 <<
"\nLeaving the original solution intact ...\n";
4877 return_status = ITER_REFINE_NOT_PERFORMED;
4880 if( iter_refine_max_iter() == 1 ) {
4881 if( (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4883 <<
"\nWarning, iter_refine_max_iter == 1. Only one step of iterative refinement"
4884 <<
"was performed and the step is taken which out checking the residual ...\n";
4889 return_status = ITER_REFINE_ONE_STEP;
4891 else if( roR_nrm_o < roR_nrm && rom_nrm_o < rom_nrm && ra_nrm_o < rom_nrm ) {
4892 if( (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4894 <<
"\nNo progress was made in reducing the residuals."
4895 <<
"\nLeaving the original solution intact ...\n";
4897 return_status = ITER_REFINE_NOT_IMPROVED;
4904 if( return_status != ITER_REFINE_CONVERGED && return_status != ITER_REFINE_NOT_NEEDED ) {
4905 if( (
int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4907 <<
"\nThe residuals were not converged but they were not increased either."
4908 <<
"\nTake the new point anyway ...\n";
4910 return_status = ITER_REFINE_IMPROVED;
4914 return return_status;
4920 const ActiveSet& act_set, std::ostream& out
4934 const int int_w = 10;
4935 const char int_ul[] =
"--------";
4936 const int dbl_min_w = 20;
4937 const int dbl_w = my_max(dbl_min_w,
int(out.precision()+8));
4938 const char dbl_ul[] =
"------------------";
4940 out <<
"\n*** Dumping the current active set ***\n"
4941 <<
"\nDimensions of the current active set:\n"
4942 <<
"\nn = " << right << setw(int_w) << qp.
n() <<
" (Number of unknowns)"
4943 <<
"\nn_R = " << right << setw(int_w) << qp.
n_R() <<
" (Number of initially free variables in Ko)"
4944 <<
"\nm = " << right << setw(int_w) << qp.
m() <<
" (Number of initially fixed variables not in Ko)"
4945 <<
"\nm_breve = " << right << setw(int_w) << constraints.
m_breve() <<
" (Number of extra general equality/inequality constriants)"
4946 <<
"\nq_hat = " << right << setw(int_w) << act_set.
q_hat() <<
" (Number of augmentations to the initial KKT system Ko)"
4947 <<
"\nq_plus_hat = " << right << setw(int_w) << act_set.
q_plus_hat() <<
" (Number of added variable bounds and general constraints)"
4948 <<
"\nq_F_hat = " << right << setw(int_w) << act_set.
q_F_hat() <<
" (Number of initially fixed variables not at their initial bound)"
4949 <<
"\nq_C_hat = " << right << setw(int_w) << act_set.
q_C_hat() <<
" (Number of initially fixed variables at the other bound)"
4950 <<
"\nq_D_hat = " << right << setw(int_w) << act_set.
q_D_hat() <<
" (Number of initially fixed variables still fixed at initial bound)"
4954 out <<
"\nQuantities for augmentations to the initial KKT system:\n";
4955 const size_type q_hat = act_set.
q_hat();
4957 << right << setw(int_w) <<
"s"
4958 << right << setw(int_w) <<
"ij_map(s)"
4959 << right << setw(int_w) <<
"bnd(s)"
4960 << right << setw(dbl_w) <<
"constr_norm(s)"
4961 << right << setw(dbl_w) <<
"d_hat(s)"
4962 << right << setw(dbl_w) <<
"z_hat(s)"
4963 << right << setw(dbl_w) <<
"p_z_hat(s)"
4965 out << right << setw(int_w) << int_ul
4966 << right << setw(int_w) << int_ul
4967 << right << setw(int_w) << int_ul
4968 << right << setw(dbl_w) << dbl_ul
4969 << right << setw(dbl_w) << dbl_ul
4970 << right << setw(dbl_w) << dbl_ul
4971 << right << setw(dbl_w) << dbl_ul
4973 {
for( size_type s = 1; s <= q_hat; ++s ) {
4974 out << right << setw(int_w) << s
4975 << right << setw(int_w) << act_set.
ij_map(s)
4976 << right << setw(int_w) << bnd_str(act_set.
bnd(s))
4978 << right << setw(dbl_w) << act_set.
d_hat()(s)
4979 << right << setw(dbl_w) << act_set.
z_hat()(s)
4980 << right << setw(dbl_w) << act_set.
p_z_hat()(s)
4985 out <<
"\nP_XF_hat =\n" << act_set.
P_XF_hat();
4986 out <<
"\nP_FC_hat =\n" << act_set.
P_FC_hat();
4987 out <<
"\nP_plus_hat =\n" << act_set.
P_plus_hat();
4988 out <<
"\nU_hat =\n" << act_set.
U_hat();
4990 out <<
"\nS_hat =\n" << act_set.
S_hat();
4993 out <<
"\nQuantities for initially fixed variables which are still fixed at their initial bound:\n";
4994 const size_type q_D_hat = act_set.
q_D_hat();
4996 << right << setw(int_w) <<
"k"
4997 << right << setw(int_w) <<
"l_fxfx(k)"
4998 << right << setw(dbl_w) <<
"mu_D_hat(k)"
4999 << right << setw(dbl_w) <<
"p_mu_D_hat(s)"
5001 out << right << setw(int_w) << int_ul
5002 << right << setw(int_w) << int_ul
5003 << right << setw(dbl_w) << dbl_ul
5004 << right << setw(dbl_w) << dbl_ul
5006 {
for( size_type k = 1; k <= q_D_hat; ++k ) {
5007 out << right << setw(int_w) << k
5008 << right << setw(int_w) << act_set.
l_fxfx(k)
5009 << right << setw(dbl_w) << act_set.
mu_D_hat()(k)
5010 << right << setw(dbl_w) << act_set.
p_mu_D_hat()(k)
5015 out <<
"\nQ_XD_hat =\n" << act_set.
Q_XD_hat();
5017 out <<
"\n*** End dump of current active set ***\n";
5038 m_breve = constraints.
m_breve();
5040 out <<
"\n*** Original QP ***\n"
5043 <<
"\nm_breve = " << m_breve
5045 out <<
"\ng =\n" <<
g();
5046 out <<
"\nG =\n" <<
G();
5048 out <<
"\nA =\n" <<
A();
5050 throw std::logic_error(
5051 error_msg(__FILE__,__LINE__,
"QPSchurPack::QP::dump_qp(...) : Error, "
5052 "m != not supported yet!"));
5055 out <<
"\nA_bar =\n" << constraints.
A_bar();
5057 DVector c_L_bar(n+m_breve), c_U_bar(n+m_breve);
5058 {
for( size_type j = 1; j <= n+m_breve; ++j ){
5062 out <<
"\nc_L_bar =\n" << c_L_bar();
5063 out <<
"\nc_U_bar =\n" << c_U_bar();
5065 out <<
"\n*** Initial KKT system (fixed and free variables) ***\n"
5066 <<
"\nn_R = " <<
n_R
5068 out <<
"\nb_X =\n" <<
b_X();
5069 out <<
"\nQ_R =\n" <<
Q_R();
5070 out <<
"\nQ_X =\n" <<
Q_X();
5071 out <<
"\nKo =\n" <<
Ko();
5072 out <<
"\nfo =\n" <<
fo();
virtual size_type n() const =0
bool drop_constraint(int jd, std::ostream *out, EOutputLevel output_level, bool force_refactorization=true, bool allow_any_cond=false)
Drop a constraint from the active set then refactorize the schur complement (if forced).
static value_type DEGENERATE_MULT
Value for near degenerate lagrange multipliers.
virtual const DVectorSlice fo() const =0
size_type s_map(int ij) const
Map from a constraint or initially fixed variable to a row and column in the schur complement S_bar...
virtual const MatrixOp & A_bar() const =0
int ij_map(size_type s) const
Returns -i for row & column of S_bar for an initially fixed variable left out of Ko that became free ...
QPSchur(const schur_comp_ptr_t &schur_comp=Teuchos::null, size_type max_iter=100, value_type max_real_runtime=1e+20, value_type feas_tol=1e-8, value_type loose_feas_tol=1e-6, value_type dual_infeas_tol=1e-12, value_type huge_primal_step=1e+20, value_type huge_dual_step=1e+20, value_type warning_tol=1e-10, value_type error_tol=1e-5, size_type iter_refine_min_iter=1, size_type iter_refine_max_iter=3, value_type iter_refine_opt_tol=1e-12, value_type iter_refine_feas_tol=1e-12, bool iter_refine_at_solution=true, bool salvage_init_schur_comp=true, MSADU::PivotTolerances pivot_tols=MSADU::PivotTolerances(1e-8, 1e-11, 1e-11))
MatrixSymAddDelUpdateable MSADU
virtual const x_init_t & x_init() const =0
Return the status of a variable initially.
void initialize(const MatrixSymOp *G, const MatrixOp *A, const MatrixOp *A_bar, const GenPermMatrixSlice *Q_R, const GenPermMatrixSlice *P_XF_hat, const GenPermMatrixSlice *P_plus_hat)
Initialize.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
EBounds bnd(size_type s) const
Return which bound is active for the active constraint.
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual value_type get_bnd(size_type j, EBounds bnd) const =0
Return the bound for a constraint.
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
const MatrixSymOpNonsing & S_hat() const
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
size_type q_F_hat() const
Return the number of variables that where initially fixed but are currently free or fixed to another ...
Represents the matrix U_hat.
virtual void set_multipliers(const ActiveSet &act_set, const DVectorSlice &v, SpVector *mu, DVectorSlice *lambda, SpVector *lambda_breve)
Map from the active set to the sparse multipliers for the inequality constraints. ...
MSADU::PivotTolerances pivot_tols() const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
size_type q_C_hat() const
Return the number of variables that where initially fixed but are currently fixed to another bound...
const GenPermMatrixSlice & P_FC_hat() const
bool is_init_fixed(size_type j) const
Determine if a constriant was an initially fixed variable.
void refactorize_schur_comp()
Reinitialize the schur complement factorization for the current active set.
static void dump_act_set_quantities(const ActiveSet &act_set, std::ostream &out, bool print_S_hat=true)
Dump all the active set quantities for debugging.
DVectorSlice p_mu_D_hat()
virtual void pick_violated_policy(EPickPolicy pick_policy)=0
Set the policy used to pick a violated constraint.
EOutputLevel
Output level.
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
virtual const l_x_X_map_t & l_x_X_map() const =0
Map from full x(i) to initially fixed x_X(l).
Represents and manages the active set for the QPSchur algorithm.
void syr2k(const MatrixOp &M, BLAS_Cpp::Transp M_trans, value_type alpha, const GenPermMatrixSlice &P1, BLAS_Cpp::Transp P1_trans, const GenPermMatrixSlice &P2, BLAS_Cpp::Transp P2_trans, value_type beta, MatrixSymOp *symwo_lhs)
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
virtual const DVectorSlice g() const =0
virtual const GenPermMatrixSlice & Q_X() const =0
(Q_X().ordered_by() == BY_ROW)
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
virtual size_type cols() const
U_hat_t()
Construct uninitialized.
bool add_constraint(size_type ja, EBounds bnd_ja, bool update_steps, std::ostream *out, EOutputLevel output_level, bool force_refactorization=true, bool allow_any_cond=false)
Add a constraint to the active set then refactorize the schur complemnt (if forced).
const U_hat_t & U_hat() const
virtual const MatrixSymOpNonsing & Ko() const =0
virtual size_type m_breve() const =0
const GenPermMatrixSlice & P_plus_hat() const
const_iterator begin() const
Represents the QP to be solved by QPSchur {abstract}.
virtual void set_x(const ActiveSet &act_set, const DVectorSlice &v, DVectorSlice *x)
Set the values in x for all the variables.
value_type constr_norm(size_type s) const
Returns ||a(j)||2 where j = ij_map(s).
virtual const MatrixOp & A() const =0
If m == 0 then don't call this, it may throw an exception or worse.
bool timeout_return(StopWatchPack::stopwatch *timer, std::ostream *out, EOutputLevel output_level) const
Determine if time has run out and if we should return.
size_type q_D_hat() const
Return the number of variables that where initially fixed and are still currently fixed to their inti...
virtual const GenPermMatrixSlice & Q_R() const =0
(Q_R().ordered_by() == BY_ROW)
virtual Constraints & constraints()=0
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
size_type q_plus_hat() const
Return the number of constraints from A_bar added to the active set.
virtual void ignore(size_type j)=0
Inform to ignore the jth constraint the next time pick_violated(...) is called.
virtual const i_x_X_map_t & i_x_X_map() const =0
Map from initially fixed x_X(l) to full x(i).
ESolveReturn
solve_qp return values
bool drop_add_constraints(int jd, size_type ja, EBounds bnd_ja, bool update_steps, std::ostream *out, EOutputLevel output_level)
Drop a constraint from, then add a constraint to the active set and refactorize the schur complement...
size_type q_hat() const
Return the total size of the schur complement.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
bool all_dof_used_up() const
Returns true if all the degrees of freedom of the QP are used up.
const ActiveSet & act_set() const
Return a reference to the active set object.
virtual size_type m() const =0
virtual ESolveReturn solve_qp(QP &qp, size_type num_act_change, const int ij_act_change[], const EBounds bnds[], std::ostream *out, EOutputLevel output_level, ERunTests test_what, DVectorSlice *x, SpVector *mu, DVectorSlice *lambda, SpVector *lambda_breve, size_type *iter, size_type *num_adds, size_type *num_drops)
Solve a QP.
Utility class for a ranged check vector.
size_type l_fxfx(size_type k) const
Returns the indice of x_X(l) of the initially fixed variables that are still fixed at their original ...
const DVectorSlice d_hat() const
virtual void dump_qp(std::ostream &out)
Dump the definition of the QP to a stream.
const GenPermMatrixSlice & Q_XD_hat() const
virtual ESolveReturn qp_algo(EPDSteps first_step, std::ostream *out, EOutputLevel output_level, ERunTests test_what, const DVectorSlice &vo, ActiveSet *act_set, DVectorSlice *v, DVectorSlice *x, size_type *iter, size_type *num_adds, size_type *num_drops, size_type *iter_refine_num_resid, size_type *iter_refine_num_solves, StopWatchPack::stopwatch *timer)
Run the algorithm from a dual feasible point.
ERunTests
Enumeration for if to run internal tests or not.
Represents the extra constraints in the QP to be satisfied by the schur complement QP solver QPSchur ...
const GenPermMatrixSlice & P_XF_hat() const
virtual const MatrixSymOp & G() const =0
EIterRefineReturn iter_refine(const ActiveSet &act_set, std::ostream *out, EOutputLevel output_level, const value_type ao, const DVectorSlice *bo, const value_type aa, const DVectorSlice *ba, DVectorSlice *v, DVectorSlice *z, size_type *iter_refine_num_resid, size_type *iter_refine_num_solves)
Perform iterative refinement on the augmented KKT system for the current active set.
void M_StMtInvMtM(MatrixSymOp *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, const MatrixSymNonsing &mswof, MatrixSymNonsing::EMatrixDummyArg mwo_rhs)
void initialize(QP &qp, size_type num_act_change, const int ij_act_change[], const EBounds bnds[], bool test, bool salvage_init_schur_comp, std::ostream *out, EOutputLevel output_level)
Initialize with an additional active set.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
virtual void pick_violated(const DVectorSlice &x, size_type *j_viol, value_type *constr_val, value_type *viol_bnd_val, value_type *norm_2_constr, EBounds *bnd, bool *can_ignore) const =0
Pick a violated constraint.
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
virtual const DVectorSlice b_X() const =0
virtual size_type n_R() const =0
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)