80 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_LOQO
87 #include "ConstrainedOptPack/src/AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
89 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
92 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_EtaVector.hpp"
95 #include "Midynamic_cast_verbose.h"
96 #include "MiWorkspacePack.h"
103 namespace LinAlgOpPack {
109 namespace ConstrainedOptPack {
114 void QPSolverRelaxedLOQO::InitLOQOHessianJacobian::init_hess_jacob(
117 ,
const int loqo_b_stat[],
const size_type num_inequal
123 LOQO* loqo_lp = (LOQO*)_loqo_lp;
127 m_in = E ? b->size() : 0,
128 m_eq = F ? f->size() : 0;
139 loqo_lp->qnz = nd*nd + 1;
140 MALLOC( loqo_lp->Q, loqo_lp->qnz,
double );
141 MALLOC( loqo_lp->iQ, loqo_lp->qnz,
int );
142 MALLOC( loqo_lp->kQ, nd+2,
int );
145 loqo_lp->kQ[j-1] = nd*(j-1);
147 loqo_lp->iQ[ loqo_lp->kQ[j-1] + (i-1) ] = i-1;
149 loqo_lp->kQ[nd] = nd*nd;
150 loqo_lp->iQ[loqo_lp->kQ[nd]] = nd;
151 loqo_lp->kQ[nd+1] = nd*nd + 1;
156 loqo_lp->Q[nd*nd] = bigM;
163 loqo_lp->nz = (num_inequal+m_eq) * (nd+1);
164 MALLOC( loqo_lp->A, loqo_lp->nz,
double );
165 MALLOC( loqo_lp->iA, loqo_lp->nz,
int );
166 MALLOC( loqo_lp->kA, nd+2,
int );
168 if( num_inequal == m_in ) {
171 {
for(
size_type j = 1; j <= nd+1; ++j ) {
172 loqo_lp->kA[j-1] = (m_in+m_eq)*(j-1);
173 for(
size_type i = 1; i <= m_in+m_eq; ++i )
174 loqo_lp->iA[ loqo_lp->kA[j-1] + (i-1) ] = i-1;
176 loqo_lp->kA[nd+1] = (m_in+m_eq)*(nd+1);
178 DMatrixSlice A( loqo_lp->A, loqo_lp->nz, loqo_lp->m, loqo_lp->m, nd+1 );
196 DMatrixSlice A( loqo_lp->A, loqo_lp->nz, loqo_lp->m, loqo_lp->m, nd+1 );
197 for(
size_type k = 1; k <= num_inequal; ++k ) {
198 const int j = loqo_b_stat[k-1];
209 QPSolverRelaxedLOQO::QPSolverRelaxedLOQO(
210 const init_hess_jacob_ptr_t init_hess_jacob
214 :init_hess_jacob_(init_hess_jacob)
216 ,nonbinding_lag_mult_(nonbinding_lag_mult)
219 nonbinding_lag_mult_ = 1e-6;
222 QPSolverRelaxedLOQO::~QPSolverRelaxedLOQO()
224 this->release_memory();
230 QPSolverRelaxedLOQO::get_qp_stats()
const
235 void QPSolverRelaxedLOQO::release_memory()
241 QPSolverRelaxedLOQO::imp_solve_qp(
242 std::ostream*
out, EOutputLevel olevel, ERunTests test_what
265 m_in = E ? b->size() : 0,
266 m_eq = F ? f->size() : 0;
272 LOQO *loqo_lp = openlp();
281 MALLOC( loqo_lp->b, m_in+m_eq,
double );
284 MALLOC( loqo_lp->r, m_in+m_eq,
double );
293 Workspace<int> loqo_b_stat_ws(wss,m_in);
295 std::fill( loqo_b_stat.begin(), loqo_b_stat.end(), 0 );
302 eLU_itr( eL->begin(), eL->end(), eL->offset()
303 , eU->begin(), eU->end(), eU->offset(), inf_bnd );
305 DVectorSlice::iterator
306 b_itr = loqo_b.begin(),
307 r_itr = loqo_r.begin();
309 b_stat_itr = loqo_b_stat.begin();
311 for(
int k = 1; !eLU_itr.at_end(); ++k, ++eLU_itr, ++b_itr, ++r_itr, ++b_stat_itr, ++num_inequal )
314 if(eLU_itr.lbound() > -inf_bnd) {
315 *b_itr = eLU_itr.lbound();
316 *r_itr = eLU_itr.ubound() >= inf_bnd ? real_big : eLU_itr.ubound() - eLU_itr.lbound();
321 *b_itr = -eLU_itr.ubound();
322 *r_itr = eLU_itr.lbound() <= -inf_bnd ? real_big : - eLU_itr.lbound() + eLU_itr.ubound();
329 loqo_r(num_inequal+1,num_inequal+m_eq) = 0.0;
337 loqo_lp->m = num_inequal + m_eq;
344 MALLOC( loqo_lp->c, nd+1,
double );
347 loqo_c(nd+1) = bigM();
350 MALLOC( loqo_lp->l, nd+1,
double );
352 std::fill( loqo_l.begin(), loqo_l.end(), -real_big );
354 SpVectorSlice::const_iterator
357 for( ; dL_itr != dL_end; ++dL_itr )
358 loqo_l( dL_itr->indice() + dL.offset() ) = dL_itr->value();
363 MALLOC( loqo_lp->u, nd+1,
double );
365 std::fill( loqo_u.begin(), loqo_u.end(), +real_big );
367 SpVectorSlice::const_iterator
370 for( ; dU_itr != dU_end; ++dU_itr )
371 loqo_u( dU_itr->indice() + dU.offset() ) = dU_itr->value();
373 loqo_u(nd+1) = +real_big;
379 init_hess_jacob().init_hess_jacob(
380 G,bigM(),E,trans_E,b,&loqo_b_stat[0],num_inequal,F,trans_F,f
387 MALLOC( loqo_lp->x, nd+1,
double );
397 loqo_lp->quadratic = 1;
401 loqo_lp->verbose = 0;
403 case PRINT_BASIC_INFO:
404 loqo_lp->verbose = 1;
406 case PRINT_ITER_SUMMARY:
407 loqo_lp->verbose = 2;
409 case PRINT_ITER_STEPS:
410 loqo_lp->verbose = 3;
412 case PRINT_ITER_ACT_SET:
413 loqo_lp->verbose = 4;
415 case PRINT_ITER_VECTORS:
416 loqo_lp->verbose = 5;
418 case PRINT_EVERY_THING:
419 loqo_lp->verbose = 6;
429 if( out && olevel >= PRINT_BASIC_INFO ) {
430 *out <<
"\nSolving QP using LOQO ...\n";
434 const int loqo_status = solvelp(loqo_lp);
436 if( out && olevel >= PRINT_BASIC_INFO ) {
437 *out <<
"\nLOQO returned status = " << loqo_status <<
"\n";
448 *d = loqo_x_sol(1,nd);
451 *eta = loqo_x_sol(nd+1);
455 *obj_d = loqo_lp->primal_obj - (*eta + 0.5 * (*eta)*(*eta)) * bigM();
463 loqo_z(loqo_lp->z,loqo_lp->n),
464 loqo_s(loqo_lp->s,loqo_lp->n);
465 DVectorSlice::const_iterator
466 z_itr = loqo_z.begin(),
467 s_itr = loqo_s.begin();
468 typedef SpVector::element_type ele_t;
469 for(
size_type i = 1; i <= nd; ++i, ++z_itr, ++s_itr ) {
470 if( *z_itr > *s_itr && *z_itr >= nonbinding_lag_mult() ) {
472 nu->add_element(ele_t(i,-(*z_itr)));
474 else if( *s_itr > *z_itr && *s_itr >= nonbinding_lag_mult() ) {
476 nu->add_element(ele_t(i,+(*s_itr)));
480 nu->assume_sorted(
true);
485 mu->resize(m_in,num_inequal);
487 b_stat_itr = loqo_b_stat.begin();
491 loqo_v(loqo_lp->v,loqo_lp->m),
492 loqo_q(loqo_lp->q,loqo_lp->m);
493 DVectorSlice::const_iterator
494 v_itr = loqo_v.begin(),
495 q_itr = loqo_q.begin();
497 typedef SpVector::element_type ele_t;
498 for(
size_type k = 1; k <= num_inequal; ++k, ++b_stat_itr, ++v_itr, ++q_itr ) {
499 const int j = *b_stat_itr;
500 if( *v_itr > *q_itr && *v_itr >= nonbinding_lag_mult() ) {
503 mu->add_element(ele_t(-j,+(*v_itr)));
505 mu->add_element(ele_t(+j,-(*v_itr)));
507 else if( *q_itr > *v_itr && *q_itr >= nonbinding_lag_mult() ) {
509 mu->add_element(ele_t(+j,+(*q_itr)));
523 loqo_y(loqo_lp->y,loqo_lp->m);
524 DVectorSlice::const_iterator
525 y_itr = loqo_y.begin() + num_inequal;
526 DVectorSlice::iterator
527 lambda_itr = lambda->
begin();
529 for(
size_type k = 1; k <= m_eq; ++k, ++y_itr, ++lambda_itr ) {
530 *lambda_itr = -(*y_itr);
544 switch( loqo_status ) {
558 ,
false, *eta > 0.0 );
567 return qp_stats_.solution_type();
573 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_LOQO
Iterate through a set of sparse bounds.
AbstractLinAlgPack::size_type size_type
SparseVector< SparseElement< index_type, value_type >, std::allocator< SparseElement< index_type, value_type > > > SpVector
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
const LAPACK_C_Decl::f_int LAPACK_C_Decl::f_dbl_prec * A
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
C++ Standard Library compatable iterator class for accesing nonunit stride arrays of data...
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
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)
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
void Vt_S(DVectorSlice *vs_lhs, value_type alpha)
vs_lhs *= alpha (BLAS xSCAL) (*** Note that alpha == 0.0 is handeled as vs_lhs = 0.0)
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
DenseLinAlgPack::DMatrixSlice DMatrixSlice
ESolutionType
Enumeration for the type of point returned from solve_qp(...).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)