53 #include "ConstrainedOptPack_ConstraintsRelaxedStd.hpp"
54 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
55 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
56 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
57 #include "AbstractLinAlgPack_MatrixOp.hpp"
58 #include "AbstractLinAlgPack_VectorMutable.hpp"
59 #include "AbstractLinAlgPack_SpVectorClass.hpp"
60 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
61 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
62 #include "AbstractLinAlgPack_SpVectorOp.hpp"
63 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
64 #include "Teuchos_Assert.hpp"
68 ConstrainedOptPack::EBounds
69 convert_bnd_type(
int bnd_type )
73 return ConstrainedOptPack::LOWER;
75 return ConstrainedOptPack::EQUALITY;
77 return ConstrainedOptPack::UPPER;
81 return ConstrainedOptPack::LOWER;
85 AbstractLinAlgPack::value_type get_sparse_element(
86 const AbstractLinAlgPack::SpVectorSlice& v
90 const AbstractLinAlgPack::SpVectorSlice::element_type
91 *ele_ptr = v.lookup_element(i);
92 return ele_ptr ? ele_ptr->value() : 0.0;
97 namespace ConstrainedOptPack {
98 namespace QPSchurPack {
103 :inequality_pick_policy_(ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY)
111 ,last_added_bound_(0.0)
112 ,last_added_bound_type_(FREE)
113 ,next_undecomp_f_k_(0)
117 const VectorSpace::space_ptr_t &space_d_eta
133 ,value_type bounds_tol
134 ,value_type inequality_tol
135 ,value_type equality_tol
139 nd = space_d_eta->dim() - 1,
147 dL && !dU, std::invalid_argument
148 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
149 "If dL!=NULL then dU!=NULL must also be true." );
150 TEUCHOS_TEST_FOR_EXCEPTION(
151 E && ( !b || !eL || !eU ), std::invalid_argument
152 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
153 "If E!=NULL then b!=NULL, eL!=NULL and eU!=NULL must also be true." );
154 TEUCHOS_TEST_FOR_EXCEPTION(
155 F && !f, std::invalid_argument
156 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
157 "If F!=NULL then f!=NULL must also be true." );
161 const size_type dL_dim = dL->dim(), dU_dim = dU->dim();
162 TEUCHOS_TEST_FOR_EXCEPTION(
163 dL_dim != nd, std::invalid_argument
164 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
165 "dL.dim() != d->dim()." );
166 TEUCHOS_TEST_FOR_EXCEPTION(
167 dU_dim != nd, std::invalid_argument
168 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
169 "dU.dim() != d->dim()." );
179 Ed_dim = Ed ? Ed->dim() : 0;
181 TEUCHOS_TEST_FOR_EXCEPTION(
182 opE_cols != nd, std::invalid_argument
183 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
184 "op(E).cols() != nd." );
185 TEUCHOS_TEST_FOR_EXCEPTION(
186 b_dim != m_in, std::invalid_argument
187 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
188 "b->dim() != op(E).rows()." );
189 TEUCHOS_TEST_FOR_EXCEPTION(
190 eL_dim != m_in, std::invalid_argument
191 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
192 "eL->dim() != op(E).rows()." );
193 TEUCHOS_TEST_FOR_EXCEPTION(
194 eU_dim != m_in, std::invalid_argument
195 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
196 "eU->dim() != op(E).rows()." );
197 TEUCHOS_TEST_FOR_EXCEPTION(
198 Ed && Ed_dim != m_in, std::invalid_argument
199 ,
"ConstraintsRelaxedStd::initialize(...) : Error, "
200 "Ed->dim() != op(E).rows()." );
209 TEUCHOS_TEST_FOR_EXCEPTION(
210 opF_cols != nd, std::invalid_argument
211 ,
"QPSolverRelaxed::solve_qp(...) : Error, "
212 "op(F).cols() != nd." );
213 TEUCHOS_TEST_FOR_EXCEPTION(
214 f_dim != m_eq, std::invalid_argument
215 ,
"QPSolverRelaxed::solve_qp(...) : Error, "
216 "f->dim() != op(F).rows()." );
221 space_d_eta,m_in,m_eq,E,trans_E,b,F,trans_F,f,m_undecomp,j_f_undecomp);
229 bounds_tol_ = bounds_tol;
230 inequality_tol_ = inequality_tol;
231 equality_tol_ = equality_tol;
233 next_undecomp_f_k_ = m_undecomp ? 1 : 0;
247 return A_bar_.
nd() + 1;
252 return A_bar_.
m_in() + A_bar_.
m_eq();
262 switch(pick_policy) {
264 inequality_pick_policy_ = ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY;
267 inequality_pick_policy_ = ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY;
277 switch(inequality_pick_policy_) {
278 case ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY:
280 case ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY:
282 case ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY:
283 return MOST_VIOLATED;
291 const DVectorSlice& x_in,
size_type* j_viol, value_type* constr_val
292 ,value_type* viol_bnd_val, value_type* norm_2_constr, EBounds* bnd,
bool* can_ignore
295 namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
300 x_in.dim() != A_bar_.
nd()+1, std::length_error
301 ,
"ConstraintsRelaxedStd::pick_violated(...) : Error, "
302 "x is the wrong length" );
308 VectorSpace::vec_mut_ptr_t
310 (VectorDenseMutableEncap(*x)()) = x_in;
311 VectorSpace::vec_mut_ptr_t
312 d = x->sub_view(1,nd);
314 eta = x->get_ele(nd+1);
316 bool Ed_computed =
false;
320 if( check_F_ && A_bar_.
F() ) {
398 value_type max_bound_viol = 0.0;
399 value_type max_bound_d_viol = 0.0;
400 value_type max_bound_dLU_viol = 0.0;
401 int max_bound_viol_type = -2;
402 if( dL_ && ( dL_->nz() || dU_->nz() ) ) {
406 ,&max_bound_viol_j, &max_bound_viol
407 ,&max_bound_d_viol, &max_bound_viol_type, &max_bound_dLU_viol
409 if( max_bound_viol > bounds_tol_ ) {
411 *j_viol = max_bound_viol_j;
412 *constr_val = max_bound_d_viol;
413 *norm_2_constr = 1.0;
414 *bnd = convert_bnd_type(max_bound_viol_type);
415 *viol_bnd_val = max_bound_dLU_viol;
419 max_bound_viol_j = 0;
423 if( ( inequality_pick_policy_ == ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY
424 || inequality_pick_policy_ == ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY )
429 last_added_j_ = *j_viol;
430 last_added_bound_type_ = *bnd;
431 last_added_bound_ = *viol_bnd_val;
439 value_type max_inequality_viol = 0.0;
440 value_type max_inequality_e_viol = 0.0;
441 value_type max_inequality_eLU_viol = 0.0;
442 int max_inequality_viol_type = -2;
444 if( inequality_pick_policy_ == ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY ) {
448 true, std::logic_error
449 ,
"ConstraintsRelaxedStd::pick_violated(...) : Error,\n"
450 "The option ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY has not been implemented yet\n" );
454 if( A_bar_.
m_in() && ( eL_->nz() || eU_->nz() ) ) {
456 VectorSpace::vec_mut_ptr_t e = eL_->space().create_member();
457 LinAlgOpPack::V_MtV( e.get(), *A_bar_.
E(), A_bar_.
trans_E(), *d );
466 ,&max_inequality_viol_j, &max_inequality_viol
467 ,&max_inequality_e_viol, &max_inequality_viol_type, &max_inequality_eLU_viol
469 if( max_inequality_viol > max_bound_viol
470 && max_inequality_viol > inequality_tol_ )
472 *j_viol = max_inequality_viol_j + nd + 1;
473 *constr_val = max_inequality_e_viol;
474 *norm_2_constr = 1.0;
475 *bnd = convert_bnd_type(max_inequality_viol_type);
476 *viol_bnd_val = max_inequality_eLU_viol;
480 max_inequality_viol_j = 0;
485 if( max_bound_viol_j || max_inequality_viol_j ) {
487 last_added_j_ = *j_viol;
488 last_added_bound_type_ = *bnd;
489 last_added_bound_ = *viol_bnd_val;
494 if( Ed_ && !Ed_computed ) {
496 LinAlgOpPack::V_MtV( Ed_, *A_bar_.
E(), A_bar_.
trans_E(), *d );
501 *norm_2_constr = 0.0;
509 true, std::logic_error
510 ,
"ConstraintsRelaxedStd::ignore(...) : Error, "
511 "This operation is not supported yet!" );
516 const value_type inf = std::numeric_limits<value_type>::max();
519 j > A_bar_.cols(), std::range_error
520 ,
"ConstraintsRelaxedStd::get_bnd(j,bnd) : Error, "
521 "j = " << j <<
" is not in range [1," << A_bar_.cols() <<
"]" );
524 if( j == last_added_j_ && bnd == last_added_bound_type_ ) {
525 return last_added_bound_;
530 const SpVectorSlice::element_type *ele_ptr = NULL;
531 if( j_local <= A_bar_.
nd() ) {
536 return dL_->get_ele(j_local);
538 return dU_->get_ele(j_local);
544 return ( bnd ==
LOWER ? -1.0 : +1.0 ) * inf;
547 else if( (j_local -= A_bar_.
nd()) <= 1 ) {
558 else if( (j_local -= 1) <= A_bar_.
m_in() ) {
562 return eL_->get_ele(j_local);
564 return eU_->get_ele(j_local);
569 else if( (j_local -= A_bar_.
m_in()) <= A_bar_.
m_eq() ) {
574 return -A_bar_.
f()->get_ele(j_local);
582 void ConstraintsRelaxedStd::cache_last_added(
583 size_type last_added_j, value_type last_added_bound
584 ,EBounds last_added_bound_type
587 last_added_j_ = last_added_j;
588 last_added_bound_ = last_added_bound;
589 last_added_bound_type_ = last_added_bound_type;
604 ,space_cols_(Teuchos::null)
609 const VectorSpace::space_ptr_t &space_d_eta
622 namespace mmp = MemMngPack;
623 namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
625 const size_type nd = space_d_eta->dim() - 1;
628 const bool test_setup =
true;
629 if( m_undecomp > 0 && f->dim() > m_undecomp ) {
630 P_u_row_i_.resize(m_undecomp);
631 P_u_col_j_.resize(m_undecomp);
633 *j_f_u = j_f_undecomp;
635 row_i_itr = P_u_row_i_.begin();
637 col_j_itr = P_u_col_j_.begin();
638 for(
size_type i = 1; i <= m_undecomp; ++i, ++j_f_u, ++row_i_itr, ++col_j_itr ) {
642 P_u_.initialize(nd,m_undecomp,m_undecomp,0,0,GPMSIP::BY_ROW_AND_COL
643 ,&P_u_row_i_[0],&P_u_col_j_[0],test_setup);
645 else if( m_undecomp > 0) {
647 P_u_.initialize(m_undecomp,m_undecomp,GenPermMatrixSlice::IDENTITY_MATRIX);
651 space_cols_ = space_d_eta;
654 VectorSpace::space_ptr_t row_spaces[3];
655 int num_row_spaces = 1;
656 row_spaces[0] = space_d_eta;
663 VectorSpace::space_ptr_t
670 row_spaces[num_row_spaces++] = vs;
672 space_rows_.initialize( row_spaces, num_row_spaces, space_d_eta->small_vec_spc_fcty() );
675 nd_ = space_d_eta->dim() - 1;
777 ,
const Vector& x, value_type b
782 namespace mmp = MemMngPack;
798 E_start = nd() + 1 + 1,
799 F_start = E_start + m_in(),
800 F_end = F_start + m_eq() - 1;
802 d_rng = Range1D(1,nd()),
803 E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
804 F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
821 VectorMutable::vec_mut_ptr_t
822 y1 = y->sub_view(d_rng);
824 y2 = y->get_ele(nd()+1);
826 x1 = x.sub_view(d_rng);
828 x2 = x.get_ele(nd()+1);
830 x3 = m_in() ? x.sub_view(E_rng) : Teuchos::null,
831 x4 = m_eq() ? x.sub_view(F_rng) : Teuchos::null;
834 Vp_StV( y1.get(), a, *x1 );
842 y2 += - a *
dot( *this->b(), *x3 );
844 y2 += - a *
dot( *f(), *x4 );
845 y->set_ele(nd()+1,y2);
861 VectorMutable::vec_mut_ptr_t
862 y1 = y->sub_view(d_rng);
864 y2 = y->get_ele(nd()+1);
865 VectorMutable::vec_mut_ptr_t
866 y3 = m_in() ? y->sub_view(E_rng) : Teuchos::null,
867 y4 = m_eq() ? y->sub_view(F_rng) : Teuchos::null;
869 x1 = x.sub_view(d_rng);
871 x2 = x.get_ele(nd()+1);
873 Vp_StV( y1.get(), a, *x1 );
876 y->set_ele(nd()+1,y2);
879 Vp_StMtV( y3.get(), a, *E(), trans_E(), *x1 );
880 Vp_StV( y3.get(), - a * x2, *this->b() );
884 Vp_StMtV( y4.get(), a, *F(), trans_F(), *x1 );
885 Vp_StV( y4.get(), - a * x2, *f() );
894 VectorMutable* y_out, value_type a
897 ,
const SpVectorSlice& x, value_type beta
const MatrixConstraints & A_bar_relaxed() const
size_type m_breve() const
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
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)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
void Vp_StPtMtV(VectorMutable *vs_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const SpVectorSlice &sv_rhs3, value_type beta) const
value_type get_bnd(size_type j, EBounds bnd) const
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
Here the next violated constraint to add to the active set is selected.
const VectorSpace & space_rows() const
BLAS_Cpp::Transp trans_E() const
MatrixOp & operator=(const MatrixOp &m)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void initialize(const VectorSpace::space_ptr_t &space_d_eta, const size_type m_in, const size_type m_eq, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, size_type m_undecomp, const size_type j_f_undecomp[])
Initialize.
const VectorSpace & space_cols() const
const MatrixOp & A_bar() const
Represents the constraints matrix.
bool max_inequ_viol(const AbstractLinAlgPack::Vector &v, const AbstractLinAlgPack::Vector &vL, const AbstractLinAlgPack::Vector &vU, AbstractLinAlgPack::size_type *max_viol_i, AbstractLinAlgPack::value_type *max_viol, AbstractLinAlgPack::value_type *v_i, int *bnd_type, AbstractLinAlgPack::value_type *vLU_i)
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)
EPickPolicy pick_violated_policy() const
const MatrixOp * E() const
void Vp_StMtV(VectorMutable *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &vs_rhs2, value_type beta) const
ConstraintsRelaxedStd()
Constructs to uninitialized.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
void initialize(const VectorSpace::space_ptr_t &space_d_eta, value_type etaL, const Vector *dL, const Vector *dU, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const Vector *eL, const Vector *eU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, size_type m_undecomp, const size_type j_f_undecomp[], VectorMutable *Ed, bool check_F=true, value_type bounds_tol=1e-10, value_type inequality_tol=1e-10, value_type equality_tol=1e-10)
Initialize constriants.
Transp trans_not(Transp _trans)
const MatrixOp * F() const
MatrixConstraints()
Construct to unitinitialized.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)