ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_QPSolverRelaxedQPKWIK.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 
43 #include "Moocho_ConfigDefs.hpp"
44 
45 
46 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_QPKWIK
47 
48 
49 #include <assert.h>
50 
51 #include <vector>
52 
53 #include "ConstrainedOptPack_QPSolverRelaxedQPKWIK.hpp"
54 #include "AbstractLinAlgPack_SpVectorClass.hpp"
55 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
56 #include "AbstractLinAlgPack_EtaVector.hpp"
57 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
58 #include "AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
59 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
60 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
61 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
62 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
63 #include "AbstractLinAlgPack_sparse_bounds.hpp"
64 #include "AbstractLinAlgPack_SpVectorOp.hpp"
65 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
66 #include "Teuchos_dyn_cast.hpp"
67 #include "Teuchos_Assert.hpp"
68 
69 namespace QPKWIKNEW_CppDecl {
70 
71 // Declarations that will link to the fortran object file.
72 // These may change for different platforms
73 
74 using FortranTypes::f_int; // INTEGER
75 using FortranTypes::f_real; // REAL
76 using FortranTypes::f_dbl_prec; // DOUBLE PRECISION
77 using FortranTypes::f_logical; // LOGICAL
78 
79 // ////////////////////////////////////////////////////
80 // Declarations to link with Fortran QPKWIK procedures
81 
82 namespace Fortran {
83 extern "C" {
84 
85 FORTRAN_FUNC_DECL_UL(void,QPKWIKNEW,qpkwiknew) (
86  const f_int& N, const f_int& M1, const f_int& M2, const f_int& M3
87  ,const f_dbl_prec GRAD[], f_dbl_prec UINV[], const f_int& LDUINV
88  ,const f_int IBND[], const f_dbl_prec BL[], const f_dbl_prec BU[]
89  ,const f_dbl_prec A[], const f_int& LDA, const f_dbl_prec YPY[]
90  ,const f_int& IYPY, const f_int& WARM, f_dbl_prec NUMPARAM[], const f_int& MAX_ITER
91  ,f_dbl_prec X[], f_int* NACTSTORE, f_int IACTSTORE[], f_int* INF
92  ,f_int* NACT, f_int IACT[], f_dbl_prec UR[], f_dbl_prec* EXTRA
93  ,f_int* ITER, f_int* NUM_ADDS, f_int* NUM_DROPS
94  ,f_int ISTATE[], const f_int& LRW, f_dbl_prec RW[]
95  );
96 
97 FORTRAN_FUNC_DECL_UL_(f_int,QPKWIKNEW_LISTATE,qpkwiknew_listate) (
98  const f_int& n, const f_int& m1, const f_int& m2, const f_int& m3);
99 
100 FORTRAN_FUNC_DECL_UL_(f_int,QPKWIKNEW_LRW,qpkwiknew_lrw) (
101  const f_int& n, const f_int& m1, const f_int& m2, const f_int& m3);
102 
103 } // end extern "C"
104 } // end namespace Fortran
105 
106 // //////////////////////////////////
107 // QPKWIK interface functions
108 
109 // Solve a QP using QPKWIK.
110 //
111 // See the Fortran file for documentation. C++ programs should use this interface.
112 inline
113 void qpkwiknew (
114  const f_int& n, const f_int& m1, const f_int& m2, const f_int& m3
115  ,const f_dbl_prec grad[], f_dbl_prec uinv[], const f_int& lduinv
116  ,const f_int ibnd[], const f_dbl_prec bl[], const f_dbl_prec bu[]
117  ,const f_dbl_prec a[], const f_int& lda, const f_dbl_prec ypy[]
118  ,const f_int& iypy, const f_int& warm, f_dbl_prec numparam[], const f_int& max_iter
119  ,f_dbl_prec x[], f_int* nactstore, f_int iactstore[], f_int* inf
120  ,f_int* nact, f_int iact[], f_dbl_prec ur[], f_dbl_prec* extra
121  ,f_int* iter, f_int* num_adds, f_int* num_drops
122  ,f_int istate[], const f_int& lrw, f_dbl_prec rw[]
123  )
124 {
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
129  , lrw, rw
130  );
131 }
132 
133 // Get the length of the integer state variables
134 inline
135 f_int qpkwiknew_listate(const f_int& n, const f_int& m1, const f_int& m2
136  , const f_int& m3)
137 {
138  return Fortran::FORTRAN_FUNC_CALL_UL_(QPKWIKNEW_LISTATE,qpkwiknew_listate) (n, m1, m2, m3);
139 }
140 
141 // Get the length of the real (double precision) workspace
142 inline
143 f_int qpkwiknew_lrw(const f_int& n, const f_int& m1, const f_int& m2
144  , const f_int& m3)
145 {
146  return Fortran::FORTRAN_FUNC_CALL_UL_(QPKWIKNEW_LRW,qpkwiknew_lrw) (n, m1, m2, m3);
147 }
148 
149 } // end namespace QPKWIKNEW_CppDecl
150 
151 // /////////////////////////////////////
152 // Local helpers
153 
154 namespace {
155 
156 template< class T >
157 inline
158 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
159 
160 using FortranTypes::f_int;
161 typedef DenseLinAlgPack::value_type value_type;
162 
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" };
165 
166 EConstraintType constraint_type( const f_int m1, const f_int m2, const f_int m3, const f_int j )
167 {
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;
175  return NU_L; // should never be exectuted
176 }
177 
178 f_int constraint_index( const f_int m1, const f_int m2, const f_int m3, const f_int ibnd[]
179  , const EConstraintType type, const f_int j )
180 {
181  switch(type) {
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;
188  }
190  return 0; // should never be exectuted
191 }
192 
193 } // end namespace
194 
195 // ///////////////////////////////////////
196 // Members for QPSolverRelaxedQPKWIK
197 
198 namespace ConstrainedOptPack {
199 
201  value_type max_qp_iter_frac
202  ,value_type infinite_bound
203  )
204  :max_qp_iter_frac_(max_qp_iter_frac)
205  ,infinite_bound_(infinite_bound)
206  ,N_(0)
207  ,M1_(0)
208  ,M2_(0)
209  ,M3_(0)
210 {
211  NUMPARAM_[0] = 1e-10; // SMALL
212  NUMPARAM_[1] = 1e-20; // VSMALL
213  NUMPARAM_[2] = 1e+20; // VLARGE
214 }
215 
217 {
218  this->release_memory();
219 }
220 
221 // Overridden from QPSolverRelaxed
222 
223 QPSolverStats
225 {
226  return qp_stats_;
227 }
228 
230 {
231  // Todo: resize to zero all the workspace!
232 }
233 
236  std::ostream* out, EOutputLevel olevel, ERunTests test_what
237  ,const Vector& g, const MatrixSymOp& G
238  ,value_type etaL
239  ,const Vector* dL, const Vector* dU
240  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
241  ,const Vector* eL, const Vector* eU
242  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
243  ,value_type* obj_d
244  ,value_type* eta, VectorMutable* d
245  ,VectorMutable* nu
246  ,VectorMutable* mu, VectorMutable* Ed
247  ,VectorMutable* lambda, VectorMutable* Fd
248  )
249 {
250  using Teuchos::dyn_cast;
251  using DenseLinAlgPack::nonconst_tri_ele;
252  using LinAlgOpPack::dot;
253  using LinAlgOpPack::V_StV;
254  using LinAlgOpPack::assign;
255  using LinAlgOpPack::V_StV;
256  using LinAlgOpPack::V_MtV;
260  using ConstrainedOptPack::MatrixExtractInvCholFactor;
261 
262  // /////////////////////////
263  // Map to QPKWIK input
264 
265  // Validate that rHL is of the proper type.
266  const MatrixExtractInvCholFactor &cG
267  = dyn_cast<const MatrixExtractInvCholFactor>(G);
268 
269  // Determine the number of sparse bounds on variables and inequalities.
270  // By default set for the dense case
271  const value_type inf = this->infinite_bound();
272  const size_type
273  nd = d->dim(),
274  m_in = E ? b->dim() : 0,
275  m_eq = F ? f->dim() : 0,
276  nvarbounds = dL ? num_bounded(*dL,*dU,inf) : 0,
277  ninequbounds = E ? num_bounded(*eL,*eU,inf) : 0,
278  nequalities = F ? f->dim() : 0;
279 
280  // Determine if this is a QP with a structure different from the
281  // one just solved.
282 
283  const bool same_qp_struct = ( N_ == nd && M1_ == nvarbounds && M2_ == ninequbounds && M3_ == nequalities );
284 
286  // Set the input parameters to be sent to QPKWIKNEW
287 
288  // N
289  N_ = nd;
290 
291  // M1
292  M1_ = nvarbounds;
293 
294  // M2
295  M2_ = ninequbounds;
296 
297  // M3
298  M3_ = nequalities;
299 
300  // GRAD
301  GRAD_ = VectorDenseEncap(g)();
302 
303  // UINV_AUG
304  //
305  // UINV_AUG = [ sqrt(bigM) 0 ]
306  // [ 0 L' ]
307  //
308  UINV_AUG_.resize(N_+1,N_+1);
309  cG.extract_inv_chol( &nonconst_tri_ele( UINV_AUG_(2,N_+1,2,N_+1), BLAS_Cpp::upper ) );
310  UINV_AUG_(1,1) = 1.0 / ::sqrt( NUMPARAM_[2] );
311  UINV_AUG_.col(1)(2,N_+1) = 0.0;
312  UINV_AUG_.row(1)(2,N_+1) = 0.0;
313 
314  // LDUINV_AUG
315  LDUINV_AUG_ = UINV_AUG_.rows();
316 
317  // IBND, BL , BU, A, LDA, YPY
318 
319  IBND_INV_.resize( nd + m_in);
320  std::fill( IBND_INV_.begin(), IBND_INV_.end(), 0 ); // Initialize the zero
321  IBND_.resize( my_max( 1, M1_ + M2_ ) );
322  BL_.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 ) );
326  YPY_.resize( my_max( 1, M1_ + M2_ ) );
327  if(M1_)
328  YPY_(1,M1_) = 0.0; // Must be for this QP interface
329 
330  // Initialize variable bound constraints
331  if( dL ) {
332  VectorDenseEncap dL_de(*dL);
333  VectorDenseEncap dU_de(*dU);
334  // read iterators
336  dLU_itr( dL_de().begin(), dL_de().end()
337  ,dU_de().begin(), dU_de().end()
338  ,inf );
339  // written iterators
340  IBND_t::iterator
341  IBND_itr = IBND_.begin(),
342  IBND_end = IBND_.begin() + M1_;
343  DVector::iterator
344  BL_itr = BL_.begin(),
345  BU_itr = BU_.begin(),
346  YPY_itr = YPY_.begin();
347  // Loop
348  for( size_type ibnd_i = 1; IBND_itr != IBND_end; ++ibnd_i, ++dLU_itr ) {
349  IBND_INV_[dLU_itr.index()-1] = ibnd_i;
350  *IBND_itr++ = dLU_itr.index();
351  *BL_itr++ = dLU_itr.lbound();
352  *BU_itr++ = dLU_itr.ubound();
353  *YPY_itr++ = 0.0; // Must be zero with this QP interface
354  }
355  }
356 
357  // Initialize inequality constraints
358 
359  if(M2_) {
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()
366  ,inf );
367  if( M2_ < m_in ) {
368  // Initialize BL, BU, YPY and A for sparse bounds on general inequalities
369  // written iterators
370  DVector::iterator
371  BL_itr = BL_.begin() + M1_,
372  BU_itr = BU_.begin() + M1_,
373  YPY_itr = YPY_.begin() + M1_;
374  IBND_t::iterator
375  ibnds_itr = IBND_.begin() + M1_;
376  // loop
377  for(size_type i = 1; i <= M2_; ++i, ++eLU_itr, ++ibnds_itr ) {
378  TEUCHOS_TEST_FOR_EXCEPT( !( !eLU_itr.at_end() ) );
379  const size_type k = eLU_itr.index();
380  *BL_itr++ = eLU_itr.lbound();
381  *BU_itr++ = eLU_itr.ubound();
382  *YPY_itr++ = b_de()(k);
383  *ibnds_itr = k; // Only for my record, not used by QPKWIK
384  IBND_INV_[nd+k-1] = M1_ + i;
385  // Add the corresponding row of op(E) to A
386  // y == A.row(i)'
387  // y' = e_k' * op(E) => y = op(E')*e_k
388  DVectorSlice y = A_.row(i);
389  EtaVector e_k(k,eL_de().dim());
390  V_MtV( &y( 1, N_ ), *E, BLAS_Cpp::trans_not(trans_E), e_k() ); // op(E')*e_k
391  }
392  }
393  else {
394  //
395  // Initialize BL, BU, YPY and A for dense bounds on general inequalities
396  //
397  // Initialize BL(M1+1:M1+M2), BU(M1+1:M1+M2)
398  // and IBND(M1+1:M1+M2) = identity (only for my record, not used by QPKWIK)
399  DVector::iterator
400  BL_itr = BL_.begin() + M1_,
401  BU_itr = BU_.begin() + M1_;
402  IBND_t::iterator
403  ibnds_itr = IBND_.begin() + M1_;
404  for(size_type i = 1; i <= m_in; ++i ) {
405  if( !eLU_itr.at_end() && eLU_itr.index() == i ) {
406  *BL_itr++ = eLU_itr.lbound();
407  *BU_itr++ = eLU_itr.ubound();
408  ++eLU_itr;
409  }
410  else {
411  *BL_itr++ = -inf;
412  *BU_itr++ = +inf;
413  }
414  *ibnds_itr++ = i;
415  IBND_INV_[nd+i-1]= M1_ + i;
416  }
417  // A(1:M2,1:N) = op(E)
418  assign( &A_(1,M2_,1,N_), *E, trans_E );
419  // YPY
420  YPY_(M1_+1,M1_+M2_) = b_de();
421  }
422  }
423 
424  // Initialize equalities
425 
426  if(M3_) {
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 );
429  }
430 
431  // IYPY
432  IYPY_ = 1; // ???
433 
434  // WARM
435  WARM_ = 0; // Cold start by default
436 
437  // MAX_ITER
438  MAX_ITER_ = static_cast<f_int>(max_qp_iter_frac() * N_);
439 
440  // INF
441  INF_ = ( same_qp_struct ? 1 : 0 );
442 
443  // Initilize output, internal state and workspace quantities.
444  if(!same_qp_struct) {
445  X_.resize(N_);
446  NACTSTORE_ = 0;
447  IACTSTORE_.resize(N_+1);
448  IACT_.resize(N_+1);
449  UR_.resize(N_+1);
450  ISTATE_.resize( QPKWIKNEW_CppDecl::qpkwiknew_listate(N_,M1_,M2_,M3_) );
451  LRW_ = QPKWIKNEW_CppDecl::qpkwiknew_lrw(N_,M1_,M2_,M3_);
452  RW_.resize(LRW_);
453  }
454 
455  // /////////////////////////////////////////////
456  // Setup a warm start form the input arguments
457  //
458  // Interestingly enough, QPKWIK sorts all of the
459  // constraints according to scaled multiplier values
460  // and mixes equality with inequality constriants.
461  // It seems to me that you should start with equality
462  // constraints first.
463 
464  WARM_ = 0;
465  NACTSTORE_ = 0;
466 
467  if( m_eq ) {
468  // Add equality constraints first since we know these will
469  // be part of the active set.
470  for( size_type j = 1; j <= m_eq; ++j ) {
471  IACTSTORE_[NACTSTORE_] = 2*M1_ + 2*M2_ + j;
472  ++NACTSTORE_;
473  }
474  }
475  if( ( nu && nu->nz() ) || ( mu && mu->nz() ) ) {
476  // Add inequality constraints
477  const size_type
478  nu_nz = nu ? nu->nz() : 0,
479  mu_nz = mu ? mu->nz() : 0;
480  // Combine all the multipliers for the bound and general inequality
481  // constraints and sort them from the largest to the smallest. Hopefully
482  // the constraints with the larger multiplier values will not be dropped
483  // from the active set.
484  SpVector gamma( nd + 1 + m_in , nu_nz + mu_nz );
485  typedef SpVector::element_type ele_t;
486  if(nu && nu_nz) {
487  VectorDenseEncap nu_de(*nu);
488  DVectorSlice::const_iterator
489  nu_itr = nu_de().begin(),
490  nu_end = nu_de().end();
491  index_type i = 1;
492  while( nu_itr != nu_end ) {
493  for( ; *nu_itr == 0.0; ++nu_itr, ++i );
494  gamma.add_element(ele_t(i,*nu_itr));
495  ++nu_itr; ++i;
496  }
497  }
498  if(mu && mu_nz) {
499  VectorDenseEncap mu_de(*mu);
500  DVectorSlice::const_iterator
501  mu_itr = mu_de().begin(),
502  mu_end = mu_de().end();
503  index_type i = 1;
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));
507  ++mu_itr; ++i;
508  }
509  }
510  std::sort( gamma.begin(), gamma.end()
512  // Now add the inequality constraints in decreasing order
513  const SpVector::difference_type o = gamma.offset();
514  for( SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
515  const size_type j = itr->index() + o;
516  const value_type val = itr->value();
517  if( j <= nd ) { // Variable bound
518  const size_type ibnd_i = IBND_INV_[j-1];
519  TEUCHOS_TEST_FOR_EXCEPT( !( ibnd_i ) );
520  IACTSTORE_[NACTSTORE_]
521  = (val < 0.0
522  ? ibnd_i // lower bound (see IACT(*))
523  : M1_ + M2_ + ibnd_i // upper bound (see IACT(*))
524  );
525  ++NACTSTORE_;
526  }
527  else if( j <= nd + m_in ) { // General inequality constraint
528  const size_type ibnd_i = IBND_INV_[j-1]; // offset into M1_ + ibnd_j
529  TEUCHOS_TEST_FOR_EXCEPT( !( ibnd_i ) );
530  IACTSTORE_[NACTSTORE_]
531  = (val < 0.0
532  ? ibnd_i // lower bound (see IACT(*))
533  : M1_ + M2_ + ibnd_i // upper bound (see IACT(*))
534  );
535  ++NACTSTORE_;
536  }
537  }
538  }
539  if( NACTSTORE_ > 0 )
540  WARM_ = 1;
541 
542  // /////////////////////////
543  // Call QPKWIK
544 
545  if( out && olevel > PRINT_NONE ) {
546  *out
547  << "\nCalling QPKWIK to solve QP problem ...\n";
548  }
549 
550  QPKWIKNEW_CppDecl::qpkwiknew(
551  N_, M1_, M2_, M3_, &GRAD_(1), &UINV_AUG_(1,1), LDUINV_AUG_, &IBND_[0]
552  ,&BL_(1), &BU_(1), &A_(1,1), LDA_, &YPY_(1), IYPY_, WARM_, NUMPARAM_, MAX_ITER_, &X_(1)
553  ,&NACTSTORE_, &IACTSTORE_[0], &INF_, &NACT_, &IACT_[0], &UR_[0], &EXTRA_
554  ,&ITER_, &NUM_ADDS_, &NUM_DROPS_, &ISTATE_[0], LRW_, &RW_[0]
555  );
556 
557  // ////////////////////////
558  // Map from QPKWIK output
559 
560  // eta
561  *eta = EXTRA_;
562  // d
563  (VectorDenseMutableEncap(*d))() = X_();
564  // nu (simple variable bounds) and mu (general inequalities)
565  if(nu) *nu = 0.0;
566  if(mu) *mu = 0.0;
567  // ToDo: Create VectorDenseMutableEncap views for faster access!
568  {for(size_type i = 1; i <= NACT_; ++i) {
569  size_type j = IACT_[i-1];
570  EConstraintType type = constraint_type(M1_,M2_,M3_,j);
571  FortranTypes::f_int idc = constraint_index(M1_,M2_,M3_,&IBND_[0],type,j);
572  switch(type) {
573  case NU_L:
574  nu->set_ele( idc , -UR_(i) );
575  break;
576  case GAMA_L:
577  mu->set_ele( IBND_[ M1_ + idc - 1 ], -UR_(i) );
578  break;
579  case NU_U:
580  nu->set_ele( idc, UR_(i)) ;
581  break;
582  case GAMA_U:
583  mu->set_ele( IBND_[ M1_ + idc - 1 ], UR_(i) );
584  break;
585  case LAMBDA:
586  lambda->set_ele( idc, UR_(i) );
587  break;
588  }
589  }}
590  // obj_d (This could be updated within QPKWIK in the future)
591  if(obj_d) {
592  // obj_d = g'*d + 1/2 * d' * G * g
593  *obj_d = dot(g,*d) + 0.5 * transVtMtV(*d,G,BLAS_Cpp::no_trans,*d);
594  }
595  // Ed (This could be updated within QPKWIK in the future)
596  if(Ed) {
597  V_MtV( Ed, *E, trans_E, *d );
598  }
599  // Fd (This could be updated within QPKWIK in the future)
600  if(Fd) {
601  V_MtV( Fd, *F, trans_F, *d );
602  }
603  // Set the QP statistics
604  QPSolverStats::ESolutionType solution_type;
605  if( INF_ >= 0 ) {
606  solution_type = QPSolverStats::OPTIMAL_SOLUTION;
607  }
608  else if( INF_ == -1 ) { // Infeasible constraints
610  true, QPSolverRelaxed::Infeasible
611  ,"QPSolverRelaxedQPKWIK::solve_qp(...) : Error, QP is infeasible" );
612  }
613  else if( INF_ == -2 ) { // LRW too small
614  TEUCHOS_TEST_FOR_EXCEPT( !( INF_ != -2 ) ); // Local programming error?
615  }
616  else if( INF_ == -3 ) { // Max iterations exceeded
617  solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
618  }
619  else {
620  TEUCHOS_TEST_FOR_EXCEPT(true); // Unknown return value!
621  }
622  qp_stats_.set_stats(
623  solution_type, QPSolverStats::CONVEX
624  ,ITER_, NUM_ADDS_, NUM_DROPS_
625  ,WARM_==1, *eta > 0.0 );
626 
627  return qp_stats_.solution_type();
628 }
629 
630 
631 } // end namespace ConstrainedOptPack
632 
633 
634 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_QPKWIK
635 
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T_To & dyn_cast(T_From &from)
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &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)
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
void assign(MatrixOp *M_lhs, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
size_t size_type
size_type num_bounded(const Vector &xl, const Vector &xu, value_type inf_bound)
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)
Transp trans_not(Transp _trans)
Transp
ESolutionType
Enumeration for the type of point returned from solve_qp(...).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)