46 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_QPKWIK
69 namespace QPKWIKNEW_CppDecl {
125 Fortran::FORTRAN_FUNC_CALL_UL(QPKWIKNEW,qpkwiknew) (
126 n, m1, m2, m3, grad, uinv, lduinv
127 , ibnd, bl, bu,
a, lda, ypy, iypy, warm, numparam, max_iter, x, nactstore
128 , iactstore,
inf, nact, iact, ur, extra, iter, num_adds, num_drops, istate
138 return Fortran::FORTRAN_FUNC_CALL_UL_(QPKWIKNEW_LISTATE,qpkwiknew_listate) (n, m1, m2, m3);
146 return Fortran::FORTRAN_FUNC_CALL_UL_(QPKWIKNEW_LRW,qpkwiknew_lrw) (n, m1, m2, m3);
158 T my_max(
const T& v1,
const T& v2 ) {
return v1 > v2 ? v1 : v2; }
163 enum EConstraintType { NU_L, NU_U, GAMA_L, GAMA_U, LAMBDA, RELAXATION };
164 char constraint_type_name[6][15] = {
"NU_L",
"NU_U",
"GAMA_L",
"GAMA_U",
"LAMBDA",
"RELAXATION" };
166 EConstraintType constraint_type(
const f_int m1,
const f_int m2,
const f_int m3,
const f_int j )
168 if (1 <= j && j <= m1 )
return NU_L;
169 else if(m1+1 <= j && j <= m1+m2 )
return GAMA_L;
170 else if(m1+m2+1 <= j && j <= 2*m1+m2 )
return NU_U;
171 else if(2*m1+m2+1 <= j && j <= 2*m1+2*m2 )
return GAMA_U;
172 else if(2*m1+2*m2+1 <= j && j <= 2*m1+2*m2+m3)
return LAMBDA;
173 else if( j == 2*m1+2*m2+m3 + 1 )
return RELAXATION;
179 ,
const EConstraintType type,
const f_int j )
182 case NU_L :
return ibnd[j-1];
183 case GAMA_L :
return j-m1;
184 case NU_U :
return ibnd[j-m1-m2-1];
185 case GAMA_U :
return j-2*m1-m2;
186 case LAMBDA :
return j-2*m1-2*m2;
187 case RELAXATION :
return 0;
198 namespace ConstrainedOptPack {
201 value_type max_qp_iter_frac
202 ,value_type infinite_bound
204 :max_qp_iter_frac_(max_qp_iter_frac)
205 ,infinite_bound_(infinite_bound)
211 NUMPARAM_[0] = 1e-10;
212 NUMPARAM_[1] = 1e-20;
213 NUMPARAM_[2] = 1e+20;
236 std::ostream*
out, EOutputLevel olevel, ERunTests test_what
237 ,
const Vector& g,
const MatrixSymOp& G
239 ,
const Vector* dL,
const Vector* dU
241 ,
const Vector* eL,
const Vector* eU
244 ,value_type* eta, VectorMutable* d
246 ,VectorMutable* mu, VectorMutable* Ed
247 ,VectorMutable* lambda, VectorMutable* Fd
260 using ConstrainedOptPack::MatrixExtractInvCholFactor;
266 const MatrixExtractInvCholFactor &cG
267 =
dyn_cast<
const MatrixExtractInvCholFactor>(G);
271 const value_type
inf = this->infinite_bound();
274 m_in = E ? b->dim() : 0,
275 m_eq = F ? f->dim() : 0,
277 ninequbounds = E ? num_bounded(*eL,*eU,inf) : 0,
278 nequalities = F ? f->dim() : 0;
283 const bool same_qp_struct = (
N_ == nd &&
M1_ == nvarbounds &&
M2_ == ninequbounds &&
M3_ == nequalities );
301 GRAD_ = VectorDenseEncap(g)();
321 IBND_.resize( my_max( 1, M1_ + M2_ ) );
323 BU_.
resize( my_max( 1, M1_ + M2_ + M3_ ) );
324 LDA_ = my_max( 1, M2_ + M3_ );
325 A_.resize(
LDA_, ( M2_ + M3_ > 0 ?
N_ : 1 ) );
332 VectorDenseEncap dL_de(*dL);
333 VectorDenseEncap dU_de(*dU);
336 dLU_itr( dL_de().begin(), dL_de().end()
337 ,dU_de().begin(), dU_de().end()
341 IBND_itr =
IBND_.begin(),
344 BL_itr =
BL_.begin(),
345 BU_itr =
BU_.begin(),
346 YPY_itr =
YPY_.begin();
348 for(
size_type ibnd_i = 1; IBND_itr != IBND_end; ++ibnd_i, ++dLU_itr ) {
350 *IBND_itr++ = dLU_itr.index();
351 *BL_itr++ = dLU_itr.lbound();
352 *BU_itr++ = dLU_itr.ubound();
360 VectorDenseEncap eL_de(*eL);
361 VectorDenseEncap eU_de(*eU);
362 VectorDenseEncap b_de(*b);
364 eLU_itr( eL_de().begin(), eL_de().end()
365 ,eU_de().begin(), eU_de().end()
371 BL_itr =
BL_.begin() +
M1_,
372 BU_itr =
BU_.begin() +
M1_,
377 for(
size_type i = 1; i <=
M2_; ++i, ++eLU_itr, ++ibnds_itr ) {
380 *BL_itr++ = eLU_itr.lbound();
381 *BU_itr++ = eLU_itr.ubound();
382 *YPY_itr++ = b_de()(k);
389 EtaVector e_k(k,eL_de().dim());
400 BL_itr =
BL_.begin() +
M1_,
401 BU_itr =
BU_.begin() +
M1_;
405 if( !eLU_itr.at_end() && eLU_itr.index() == i ) {
406 *BL_itr++ = eLU_itr.lbound();
407 *BU_itr++ = eLU_itr.ubound();
420 YPY_(M1_+1,M1_+M2_) = b_de();
427 V_StV( &
BU_( M1_ + M2_ + 1, M1_ + M2_ + M3_ ), -1.0, VectorDenseEncap(*f)() );
428 assign( &
A_( M2_ + 1, M2_ + M3_, 1,
N_ ), *F, trans_F );
441 INF_ = ( same_qp_struct ? 1 : 0 );
444 if(!same_qp_struct) {
450 ISTATE_.resize( QPKWIKNEW_CppDecl::qpkwiknew_listate(
N_,M1_,M2_,M3_) );
451 LRW_ = QPKWIKNEW_CppDecl::qpkwiknew_lrw(
N_,M1_,M2_,M3_);
475 if( ( nu && nu->nz() ) || ( mu && mu->nz() ) ) {
478 nu_nz = nu ? nu->nz() : 0,
479 mu_nz = mu ? mu->nz() : 0;
484 SpVector gamma( nd + 1 + m_in , nu_nz + mu_nz );
485 typedef SpVector::element_type ele_t;
487 VectorDenseEncap nu_de(*nu);
488 DVectorSlice::const_iterator
489 nu_itr = nu_de().begin(),
490 nu_end = nu_de().end();
492 while( nu_itr != nu_end ) {
493 for( ; *nu_itr == 0.0; ++nu_itr, ++i );
494 gamma.add_element(ele_t(i,*nu_itr));
499 VectorDenseEncap mu_de(*mu);
500 DVectorSlice::const_iterator
501 mu_itr = mu_de().begin(),
502 mu_end = mu_de().end();
504 while( mu_itr != mu_end ) {
505 for( ; *mu_itr == 0.0; ++mu_itr, ++i );
506 gamma.add_element(ele_t(i+nd,*mu_itr));
510 std::sort( gamma.begin(), gamma.end()
513 const SpVector::difference_type o = gamma.offset();
514 for( SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
516 const value_type val = itr->value();
527 else if( j <= nd + m_in ) {
547 <<
"\nCalling QPKWIK to solve QP problem ...\n";
550 QPKWIKNEW_CppDecl::qpkwiknew(
552 ,&
BL_(1), &
BU_(1), &
A_(1,1),
LDA_, &
YPY_(1),
IYPY_,
WARM_,
NUMPARAM_,
MAX_ITER_, &
X_(1)
563 (VectorDenseMutableEncap(*d))() =
X_();
570 EConstraintType type = constraint_type(M1_,M2_,M3_,j);
574 nu->set_ele( idc , -
UR_(i) );
577 mu->set_ele(
IBND_[ M1_ + idc - 1 ], -
UR_(i) );
580 nu->set_ele( idc,
UR_(i)) ;
583 mu->set_ele(
IBND_[ M1_ + idc - 1 ],
UR_(i) );
586 lambda->set_ele( idc,
UR_(i) );
597 V_MtV( Ed, *E, trans_E, *d );
601 V_MtV( Fd, *F, trans_F, *d );
608 else if(
INF_ == -1 ) {
610 true, QPSolverRelaxed::Infeasible
611 ,
"QPSolverRelaxedQPKWIK::solve_qp(...) : Error, QP is infeasible" );
613 else if(
INF_ == -2 ) {
616 else if(
INF_ == -3 ) {
625 ,
WARM_==1, *eta > 0.0 );
634 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_QPKWIK
Iterate through a set of sparse bounds.
void set_stats(ESolutionType solution_type, EConvexity convexity, int num_qp_iter, int num_adds, int num_drops, bool warm_start, bool infeasible_qp)
Initialize the statistics.
AbstractLinAlgPack::size_type size_type
const f_int const f_int const f_int const f_dbl_prec const f_dbl_prec * BL
const LAPACK_C_Decl::f_int LAPACK_C_Decl::f_dbl_prec const LAPACK_C_Decl::f_int & LDA
RTOp_index_type index_type
FortranTypes::f_int f_int
SparseVector< SparseElement< index_type, value_type >, std::allocator< SparseElement< index_type, value_type > > > SpVector
const f_int const f_int const f_dbl_prec const f_dbl_prec * X
FORTRAN_FUNC_DECL_UL_(void, QPOPT_SET_DEFAULTS, qpopt_set_defaults)()
FortranTypes::f_logical f_logical
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
int resize(OrdinalType length_in)
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
QPSolverStats::ESolutionType imp_solve_qp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const Vector &g, const MatrixSymOp &G, 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, value_type *obj_d, value_type *eta, VectorMutable *d, VectorMutable *nu, VectorMutable *mu, VectorMutable *Ed, VectorMutable *lambda, VectorMutable *Fd)
FortranTypes::f_int f_int
const f_int const f_int const f_int const f_dbl_prec const f_dbl_prec const f_dbl_prec const f_dbl_prec const f_dbl_prec f_int * ISTATE
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * op(M_rhs2) * v_rhs3
T_To & dyn_cast(T_From &from)
const f_int const f_int const f_int const f_dbl_prec const f_dbl_prec const f_dbl_prec * BU
Create an eta vector (scaled by alpha = default 1).
QPSolverStats get_qp_stats() const
size_type num_bounded(const Vector &xl, const Vector &xu, value_type inf_bound)
Count the number of finitly bounded elements in xl <= x <= xu.
DMatrixSliceTriEle nonconst_tri_ele(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a triangular element-wise matrix.
const f_int f_dbl_prec a[]
FORTRAN_FUNC_DECL_UL(void, MA28AD, ma28ad)(const f_int &n
QPSolverRelaxedQPKWIK(value_type max_qp_iter_frac=10.0, value_type infinite_bound=1e+20)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
void sqrt(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = sqrt(vs_rhs)
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
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
FortranTypes::f_dbl_prec value_type
Typedef for the value type of elements that is used for the library.
const f_int const f_int & N
Function object class for sorting a sparse vectors in descending order by abs(v(i)).
FortranTypes::f_dbl_prec f_dbl_prec
ESolutionType
Enumeration for the type of point returned from solve_qp(...).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
FortranTypes::f_real f_real
ESolutionType solution_type() const
const f_int const f_int const f_int const f_dbl_prec const f_dbl_prec const f_dbl_prec const f_dbl_prec const f_dbl_prec f_int f_dbl_prec f_int f_int & ITER