69 convert_bnd_type(
int bnd_type )
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
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;
262 switch(pick_policy) {
277 switch(inequality_pick_policy_) {
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;
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;
442 int max_inequality_viol_type = -2;
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" );
456 VectorSpace::vec_mut_ptr_t e =
eL_->space().create_member();
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 ) {
494 if(
Ed_ && !Ed_computed ) {
501 *norm_2_constr = 0.0;
509 true, std::logic_error
510 ,
"ConstraintsRelaxedStd::ignore(...) : Error, "
511 "This operation is not supported yet!" );
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() <<
"]" );
530 const SpVectorSlice::element_type *ele_ptr = NULL;
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 ) {
562 return eL_->get_ele(j_local);
564 return eU_->get_ele(j_local);
574 return -
A_bar_.
f()->get_ele(j_local);
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 ) {
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;
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;
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);
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
869 x1 = x.sub_view(d_rng);
871 x2 = x.get_ele(nd()+1);
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() );
element_type * lookup_element(size_type i)
const MatrixConstraints & A_bar_relaxed() const
AbstractLinAlgPack::size_type size_type
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta=1.0)
vs_lhs = alpha * op(mwo_rhs1) * vs_rhs2 + beta * vs_lhs.
size_type m_breve() const
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= 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)
v_lhs = alpha * op(P_rhs1) * op(M_rhs2) * v_rhs3 + beta * v_rhs
RTOp_value_type value_type
void cache_last_added(size_type last_added_j, value_type last_added_bound, EBounds last_added_bound_type) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
EBounds last_added_bound_type_
size_type next_undecomp_f_k_
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.
RTOp_index_type size_type
const VectorSpace & space_cols() const
const MatrixOp & A_bar() const
Represents the constraints matrix.
const LAPACK_C_Decl::f_int & M
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)
Compute the maximum violation from a set of inequality constraints vL <= v <= vU. ...
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)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
EPickPolicy pick_violated_policy() const
const MatrixOp * E() const
const f_int f_dbl_prec a[]
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)
result = v_rhs1' * v_rhs2
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
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)
Return the opposite of the transpose argument.
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
value_type last_added_bound_
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
Sparse storage element type.
const MatrixOp * F() const
MatrixConstraints()
Construct to unitinitialized.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return columns of a possible transposed matrix.
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
v_lhs = V_rhs1 - V_rhs2.
RangePack::Range1D Range1D
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)