ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_QPSolverRelaxedQPSchur.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <assert.h>
00043 
00044 #include <vector>
00045 
00046 #include "ConstrainedOptPack_QPSolverRelaxedQPSchur.hpp"
00047 #include "AbstractLinAlgPack_MatrixOp.hpp"
00048 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00049 #include "AbstractLinAlgPack_VectorSpaceFactory.hpp"
00050 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
00051 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00052 #include "AbstractLinAlgPack_VectorSpaceSerial.hpp"
00053 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00054 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00055 #include "Teuchos_dyn_cast.hpp"
00056 #include "ProfileHackPack_profile_hack.hpp"
00057 
00058 namespace ConstrainedOptPack {
00059 
00060 QPSolverRelaxedQPSchur::QPSolverRelaxedQPSchur(
00061   const init_kkt_sys_ptr_t&  init_kkt_sys
00062   ,const constraints_ptr_t&  constraints
00063   ,value_type     max_qp_iter_frac
00064   ,value_type         max_real_runtime
00065   ,QPSchurPack::ConstraintsRelaxedStd::EInequalityPickPolicy
00066                       inequality_pick_policy
00067   ,ELocalOutputLevel  print_level
00068   ,value_type     bounds_tol
00069   ,value_type     inequality_tol
00070   ,value_type     equality_tol
00071   ,value_type     loose_feas_tol
00072   ,value_type     dual_infeas_tol
00073   ,value_type     huge_primal_step
00074   ,value_type     huge_dual_step
00075   ,value_type     bigM
00076   ,value_type     warning_tol
00077   ,value_type     error_tol
00078   ,size_type          iter_refine_min_iter
00079   ,size_type          iter_refine_max_iter
00080   ,value_type         iter_refine_opt_tol
00081   ,value_type         iter_refine_feas_tol
00082   ,bool               iter_refine_at_solution
00083   ,value_type         pivot_warning_tol
00084   ,value_type         pivot_singular_tol
00085   ,value_type         pivot_wrong_inertia_tol
00086   ,bool               add_equalities_initially
00087   )
00088   :init_kkt_sys_(init_kkt_sys)
00089   ,constraints_(constraints)
00090   ,max_qp_iter_frac_(max_qp_iter_frac)
00091   ,max_real_runtime_(max_real_runtime)
00092   ,inequality_pick_policy_(inequality_pick_policy)
00093   ,print_level_(print_level)
00094   ,bounds_tol_(bounds_tol)
00095   ,inequality_tol_(inequality_tol)
00096   ,equality_tol_(equality_tol)
00097   ,loose_feas_tol_(loose_feas_tol)
00098   ,dual_infeas_tol_(dual_infeas_tol)
00099   ,huge_primal_step_(huge_primal_step)
00100   ,huge_dual_step_(huge_dual_step)
00101   ,bigM_(bigM)
00102   ,warning_tol_(warning_tol)
00103   ,error_tol_(error_tol)
00104   ,iter_refine_min_iter_(iter_refine_min_iter)
00105   ,iter_refine_max_iter_(iter_refine_max_iter)
00106   ,iter_refine_opt_tol_(iter_refine_opt_tol)
00107   ,iter_refine_feas_tol_(iter_refine_feas_tol)
00108   ,iter_refine_at_solution_(iter_refine_at_solution)
00109   ,pivot_warning_tol_(pivot_warning_tol)
00110   ,pivot_singular_tol_(pivot_singular_tol)
00111   ,pivot_wrong_inertia_tol_(pivot_wrong_inertia_tol)
00112   ,add_equalities_initially_(add_equalities_initially)
00113 {}
00114 
00115 QPSolverRelaxedQPSchur::~QPSolverRelaxedQPSchur()
00116 {
00117   this->release_memory();
00118 }
00119 
00120 // Overridden from QPSolverRelaxed
00121 
00122 QPSolverStats
00123 QPSolverRelaxedQPSchur::get_qp_stats() const
00124 {
00125   return qp_stats_;
00126 }
00127 
00128 void QPSolverRelaxedQPSchur::release_memory()
00129 {
00130   // Nothing to release!
00131 }
00132 
00133 QPSolverStats::ESolutionType
00134 QPSolverRelaxedQPSchur::imp_solve_qp(
00135   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00136   ,const Vector& g, const MatrixSymOp& G
00137   ,value_type etaL
00138   ,const Vector* dL, const Vector* dU
00139   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00140   ,const Vector* eL, const Vector* eU
00141   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00142   ,value_type* obj_d
00143   ,value_type* eta, VectorMutable* d
00144   ,VectorMutable* nu
00145   ,VectorMutable* mu, VectorMutable* Ed
00146   ,VectorMutable* lambda, VectorMutable* Fd
00147   )
00148 {
00149   namespace mmp = MemMngPack;
00150   using Teuchos::dyn_cast;
00151   using LinAlgOpPack::V_mV;
00152   typedef QPSchurPack::ConstraintsRelaxedStd constr_t;
00153 
00154 #ifdef PROFILE_HACK_ENABLED
00155   ProfileHackPack::ProfileTiming profile_timing( "QPSolverRelaxedQPSchur::imp_solve_qp(...)" );
00156 #endif
00157 
00158   const size_type
00159     nd   = g.dim(),
00160     m_in = E ? BLAS_Cpp::rows(E->rows(),E->cols(),trans_E) : 0,
00161     m_eq = F ? BLAS_Cpp::rows(F->rows(),F->cols(),trans_F) : 0;
00162 
00163   VectorDenseEncap  g_de(g);
00164 
00165   VectorSpace::space_ptr_t
00166     space_d_eta = d->space().small_vec_spc_fcty()->create_vec_spc(nd+1); // ToDo: Generalize!
00167 
00168   // ///////////////////////////////
00169   // Setup the initial KKT system
00170 
00171   InitKKTSystem::i_x_free_t     i_x_free;
00172   InitKKTSystem::i_x_fixed_t    i_x_fixed;
00173   InitKKTSystem::bnd_fixed_t    bnd_fixed;
00174   InitKKTSystem::j_f_decomp_t   j_f_decomp;
00175   size_type n_R_tmp;
00176   init_kkt_sys().initialize_kkt_system(
00177     g,G,etaL,dL,dU,F,trans_F,f,d,nu
00178     ,&n_R_tmp,&i_x_free,&i_x_fixed,&bnd_fixed,&j_f_decomp
00179     ,&b_X_,&Ko_,&fo_ );
00180   const size_type
00181     n_R = n_R_tmp,
00182     n_X = nd + 1 - n_R; // fixed variables in d and eta
00183   TEUCHOS_TEST_FOR_EXCEPT( !(  i_x_free.size() == 0 || i_x_free.size()  >= n_R  ) );  // Todo: Make an exception!
00184   TEUCHOS_TEST_FOR_EXCEPT( !(  i_x_fixed.size() >= n_X  ) );  // Todo: Make an exception!
00185   TEUCHOS_TEST_FOR_EXCEPT( !(  bnd_fixed.size() >= n_X  ) ); // Todo: Make and exception!
00186 
00187   // //////////////////////////////
00188   // Initialize constraints object
00189 
00190   // Setup j_f_undecomp
00191   const bool all_f_undecomp = F ? j_f_decomp.size() == 0 : true;
00192   const size_type
00193     m_undecomp = F ? f->dim() - j_f_decomp.size() : 0;
00194   typedef std::vector<size_type> j_f_undecomp_t;
00195   j_f_undecomp_t j_f_undecomp;
00196   if( m_undecomp && !all_f_undecomp ) {
00197     j_f_undecomp.resize(m_undecomp);
00198     // Create a full lookup array to determine if a constraint
00199     // is decomposed or not.  We need this to fill the array
00200     // j_f_undecomp[] (which is sorted).
00201     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
00202   }
00203 
00204   // initialize constraints object
00205   constraints_->initialize(
00206     space_d_eta
00207     ,etaL,dL,dU,E,trans_E,b,eL,eU,F,trans_F,f
00208     ,m_undecomp, m_undecomp && !all_f_undecomp ? &j_f_undecomp[0] : NULL
00209     ,Ed
00210     ,!add_equalities_initially()  // If we add equalities the the schur complement intially
00211                                   // then we don't need to check if they are violated.
00212     );
00213   // ToDo: Add j_f_decomp to the above constraints class!
00214 
00215   // ///////////////////////////
00216   // Initialize the QP object
00217 
00218   // g_relaxed_ = [ g; bigM ]
00219   g_relaxed_.resize(nd+1);
00220   g_relaxed_(1,nd) = g_de();
00221   g_relaxed_(nd+1) = bigM();
00222 
00223   // G_relaxed_ = [ G, zeros(...); zeros(...), bigM ]
00224   bigM_vec_.initialize(1); // dim == 1
00225   bigM_vec_ = bigM(); 
00226   G_relaxed_.initialize(
00227     Teuchos::rcp(&dyn_cast<const MatrixSymOpNonsing>(G),false)
00228     ,Teuchos::rcp(&bigM_vec_,false)
00229     ,space_d_eta
00230     );
00231   
00232   qp_.initialize(
00233     g_relaxed_(),G_relaxed_,NULL
00234     ,n_R, i_x_free.size() ? &i_x_free[0] : NULL
00235     ,&i_x_fixed[0],&bnd_fixed[0]
00236     ,b_X_(),*Ko_,fo_(),constraints_.get()
00237     ,out,test_what==RUN_TESTS,warning_tol(),error_tol()
00238     ,int(olevel)>=int(PRINT_ITER_VECTORS)
00239     );
00240 
00241   // ///////////////////////////////////////////////////////
00242   // Setup for a warm start (changes to initial KKT system)
00243 
00244   typedef std::vector<int>          ij_act_change_t;
00245   typedef std::vector<EBounds>        bnds_t;
00246   size_type     num_act_change = 0; // The default is a cold start
00247   const size_type     max_num_act_change = 2*nd;
00248   ij_act_change_t   ij_act_change(max_num_act_change);
00249   bnds_t        bnds(max_num_act_change);
00250   // Go ahead and add the equality constraints.  If these are linearly
00251   // dependent let's hope that QPSchur can handle this and still do a
00252   // good job of things.  This is a scary thing to do!
00253   if( m_eq && add_equalities_initially() ) {
00254     for( size_type j = 1; j <= m_eq; ++j ) {
00255       ij_act_change[num_act_change] = (nd + 1) + m_in + j;
00256       bnds[num_act_change]          = EQUALITY;
00257       ++num_act_change;
00258     }
00259   }
00260   // We will add fixed (EQUALITY) variable bounds to the initial active set
00261   // (if it is not already an intially fixed variable).  If fixing a variable
00262   // causes the KKT system to become singular then we are in real trouble!
00263   // We should add these eairly on since they will not be freed.
00264   if( dL ) {
00265     const QPSchurPack::QP::x_init_t &x_init = qp_.x_init();
00266     const value_type inf_bnd = this->infinite_bound();
00267     VectorDenseEncap dL_de(*dL);
00268     VectorDenseEncap dU_de(*dU);
00269     // read iterators
00270     AbstractLinAlgPack::sparse_bounds_itr
00271       dLU_itr( dL_de().begin(), dL_de().end()
00272           ,dU_de().begin(), dU_de().end()
00273           ,inf_bnd );
00274     for( ; !dLU_itr.at_end(); ++dLU_itr ) {
00275       if( dLU_itr.lbound() == dLU_itr.ubound() && x_init(dLU_itr.index()) == FREE ) {
00276         ij_act_change[num_act_change] = dLU_itr.index();
00277         bnds[num_act_change]          = EQUALITY;
00278         ++num_act_change;
00279       }
00280     }
00281   }
00282   // Add inequality constriants to the list from nu and mu
00283   if( ( nu && nu->nz() ) || ( m_in && mu->nz() ) ) {
00284     //
00285     // Setup num_act_change, ij_act_change, and bnds for a warm start!
00286     //
00287     const size_type
00288       nu_nz = nu ? nu->nz() : 0,
00289       mu_nz = mu ? mu->nz() : 0;
00290     // Combine all the multipliers for the bound and general inequality
00291     // constraints and sort them from the largest to the smallest.  Hopefully
00292     // the constraints with the larger multiplier values will not be dropped
00293     // from the active set.
00294     SpVector gamma( nd + 1 + m_in , nu_nz + mu_nz );
00295     typedef SpVector::element_type ele_t;
00296     if(nu && nu_nz) {
00297       VectorDenseEncap nu_de(*nu);
00298       DVectorSlice::const_iterator
00299         nu_itr = nu_de().begin(),
00300         nu_end = nu_de().end();
00301       index_type i = 1;
00302       while( nu_itr != nu_end ) {
00303         for( ; *nu_itr == 0.0; ++nu_itr, ++i );
00304         gamma.add_element(ele_t(i,*nu_itr));
00305         ++nu_itr; ++i;
00306       }
00307     }
00308     if(mu && mu_nz) {
00309       VectorDenseEncap mu_de(*mu);
00310       DVectorSlice::const_iterator
00311         mu_itr = mu_de().begin(),
00312         mu_end = mu_de().end();
00313       index_type i = 1;
00314       while( mu_itr != mu_end ) {
00315         for( ; *mu_itr == 0.0; ++mu_itr, ++i );
00316         gamma.add_element(ele_t(i+nd,*mu_itr));
00317         ++mu_itr; ++i;
00318       }
00319     }
00320     std::sort( gamma.begin(), gamma.end()
00321       , AbstractLinAlgPack::SortByDescendingAbsValue() );
00322     // Now add the inequality constraints in decreasing order (if they are
00323     // not already initially fixed variables)
00324     const QPSchurPack::QP::x_init_t &x_init = qp_.x_init();
00325     if(gamma.nz()) {
00326       const SpVector::difference_type o = gamma.offset();
00327       for( SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
00328         const size_type i =  itr->index() + o;
00329         if( i <= nd && x_init(i) != FREE )
00330           continue; // This variable is already initially fixed
00331         // This is not an initially fixed variable so add it
00332         ij_act_change[num_act_change] = i;
00333         bnds[num_act_change]
00334           = itr->value() > 0.0 ? UPPER : LOWER;
00335         ++num_act_change;
00336       }
00337     }
00338   }
00339   // We need to loop through x_init() and nu() in order and look for variables
00340   // that are initially fixed in x_init() but are not present in nu().  For these
00341   // variables we need to free them in ij_act_change[].
00342   if( dL && nu->nz() ) {
00343     QPSchurPack::QP::x_init_t::const_iterator
00344       x_init_itr = qp_.x_init().begin();
00345     VectorDenseEncap nu_de(*nu);
00346     DVectorSlice::const_iterator
00347       nu_itr = nu_de().begin();
00348     for( size_type i = 1; i <= nd; ++i, ++x_init_itr, ++nu_itr ) {
00349       if( *x_init_itr != FREE && *x_init_itr != EQUALITY ) {
00350         // This is an initially fixed upper or lower bound.
00351         // Look for lagrange multiplier stating that it is
00352         // still fixed.
00353         if( *nu_itr != 0.0 ) {
00354           // This active bound is present but lets make sure
00355           // that it is still the same bound
00356           if( ( *x_init_itr == LOWER && *nu_itr > 0 )
00357             || ( *x_init_itr == UPPER && *nu_itr < 0 ) )
00358           {
00359             // The bound has changed from upper to lower or visa-versa!
00360             ij_act_change[num_act_change] = i;
00361             bnds[num_act_change]
00362               = *nu_itr > 0.0 ? UPPER : LOWER;
00363             ++num_act_change;
00364           }
00365         }
00366         else {
00367           // This initially fixed variable is not fixed in nu so lets free it!
00368           ij_act_change[num_act_change] = -i;
00369           bnds[num_act_change]          = FREE;
00370           ++num_act_change;
00371         }
00372       }
00373     }
00374   }
00375   // Consider the relaxation variable!
00376   if( *eta > etaL) {
00377     ij_act_change[num_act_change] = -int(nd+1);
00378     bnds[num_act_change]          = FREE;
00379     ++num_act_change;
00380   }   
00381 
00382   // Set the output level
00383   QPSchur::EOutputLevel qpschur_olevel;
00384   switch( print_level() ) {
00385     case USE_INPUT_ARG: {
00386       // Use the input print level
00387       switch( olevel ) {
00388         case PRINT_NONE:
00389           qpschur_olevel = QPSchur::NO_OUTPUT;
00390           break;
00391         case PRINT_BASIC_INFO:
00392           qpschur_olevel = QPSchur::OUTPUT_BASIC_INFO;
00393           break;
00394         case PRINT_ITER_SUMMARY:
00395           qpschur_olevel = QPSchur::OUTPUT_ITER_SUMMARY;
00396           break;
00397         case PRINT_ITER_STEPS:
00398           qpschur_olevel = QPSchur::OUTPUT_ITER_STEPS;
00399           break;
00400         case PRINT_ITER_ACT_SET:
00401         case PRINT_ITER_VECTORS:
00402           qpschur_olevel = QPSchur::OUTPUT_ACT_SET;
00403           break;
00404         case PRINT_EVERY_THING:
00405           qpschur_olevel = QPSchur::OUTPUT_ITER_QUANTITIES;
00406           break;
00407         default:
00408           TEUCHOS_TEST_FOR_EXCEPT(true);
00409       }
00410       break;
00411     }
00412     case NO_OUTPUT:
00413       qpschur_olevel = QPSchur::NO_OUTPUT;
00414       break;
00415     case OUTPUT_BASIC_INFO:
00416       qpschur_olevel = QPSchur::OUTPUT_BASIC_INFO;
00417       break;
00418     case OUTPUT_ITER_SUMMARY:
00419       qpschur_olevel = QPSchur::OUTPUT_ITER_SUMMARY;
00420       break;
00421     case OUTPUT_ITER_STEPS:
00422       qpschur_olevel = QPSchur::OUTPUT_ITER_STEPS;
00423       break;
00424     case OUTPUT_ACT_SET:
00425       qpschur_olevel = QPSchur::OUTPUT_ACT_SET;
00426       break;
00427     case OUTPUT_ITER_QUANTITIES:
00428       qpschur_olevel = QPSchur::OUTPUT_ITER_QUANTITIES;
00429       break;
00430     default:
00431       TEUCHOS_TEST_FOR_EXCEPT(true);
00432   }
00433 
00434   //
00435   // Set options for ConstraintsRelaxedStd.
00436   // 
00437   if( bounds_tol() > 0.0 )
00438     constraints_->bounds_tol(bounds_tol());
00439   if( inequality_tol() > 0.0 )
00440     constraints_->inequality_tol(inequality_tol());
00441   if( equality_tol() > 0.0 )
00442     constraints_->equality_tol(equality_tol());
00443   constraints_->inequality_pick_policy(inequality_pick_policy());
00444 
00445   //
00446   // Set options for QPSchur.
00447   // 
00448   qp_solver_.max_iter(static_cast<index_type>(max_qp_iter_frac()*nd) );
00449   qp_solver_.max_real_runtime( max_real_runtime() );
00450   qp_solver_.feas_tol( constraints_->bounds_tol() );  // Let's assume the bound tolerance is the tightest
00451   if(loose_feas_tol() > 0.0)
00452     qp_solver_.loose_feas_tol( loose_feas_tol() );
00453   else
00454     qp_solver_.loose_feas_tol( 10.0 * qp_solver_.feas_tol() );
00455   if(dual_infeas_tol() > 0.0)
00456     qp_solver_.dual_infeas_tol( dual_infeas_tol() );
00457   if(huge_primal_step() > 0.0)
00458     qp_solver_.huge_primal_step( huge_primal_step() );
00459   if(huge_dual_step() > 0.0)
00460     qp_solver_.huge_dual_step( huge_dual_step() );
00461   qp_solver_.set_schur_comp( QPSchur::schur_comp_ptr_t( &schur_comp_, false ) );
00462   qp_solver_.warning_tol( warning_tol() );
00463   qp_solver_.error_tol( error_tol() );
00464   qp_solver_.iter_refine_min_iter( iter_refine_min_iter() );
00465   qp_solver_.iter_refine_max_iter( iter_refine_max_iter() );
00466   qp_solver_.iter_refine_opt_tol( iter_refine_opt_tol() );
00467   qp_solver_.iter_refine_feas_tol( iter_refine_feas_tol() );
00468   qp_solver_.iter_refine_at_solution( iter_refine_at_solution() );
00469   qp_solver_.pivot_tols(
00470     MatrixSymAddDelUpdateable::PivotTolerances(
00471       pivot_warning_tol(), pivot_singular_tol(), pivot_wrong_inertia_tol()
00472       ));
00473   
00474   //
00475   // Solve the QP with QPSchur
00476   // 
00477   DVector _x(nd+1);   // solution vector [ d; eta ]
00478   SpVector _mu;     // lagrange multipliers for variable bounds [ nu; kappa ]
00479   SpVector _lambda_breve; // solution for extra general constraints [ mu; lambda ]
00480   size_type qp_iter = 0, num_adds = 0, num_drops = 0;
00481   QPSchur::ESolveReturn
00482     solve_returned
00483       = qp_solver_.solve_qp(
00484         qp_
00485         ,num_act_change, num_act_change ? &ij_act_change[0] : NULL
00486         ,num_act_change ? &bnds[0] : NULL
00487         ,out, qpschur_olevel
00488         ,test_what==RUN_TESTS ? QPSchur::RUN_TESTS : QPSchur::NO_TESTS
00489         ,&_x(), &_mu, NULL, &_lambda_breve
00490         ,&qp_iter, &num_adds, &num_drops
00491         );
00492 
00493   // Set the solution
00494 
00495   // d
00496   (VectorDenseMutableEncap(*d))() = _x(1,nd);
00497   // nu
00498   if( nu ) {
00499     *nu = 0.0;
00500     const SpVector::difference_type o = _mu.offset();
00501     if( _mu.nz() ) {
00502       for(SpVector::const_iterator _mu_itr = _mu.begin(); _mu_itr != _mu.end(); ++_mu_itr) {
00503         typedef SpVector::element_type ele_t;
00504         if( _mu_itr->index() + o <= nd ) // don't add multiplier for eta <= etaL
00505           nu->set_ele( _mu_itr->index() + o, _mu_itr->value() );
00506       }
00507     }
00508   }
00509   // mu, lambda 
00510   if( m_in || m_eq ) {
00511     *eta = _x(nd+1);  // must be non-null
00512     *mu = 0.0;
00513     const SpVector::difference_type o = _lambda_breve.offset();
00514     if(_lambda_breve.nz()) {
00515       for(SpVector::const_iterator itr = _lambda_breve.begin();
00516         itr != _lambda_breve.end();
00517         ++itr)
00518       {
00519         typedef SpVector::element_type ele_t;
00520         if( itr->index() + o <= m_in ) {
00521           mu->set_ele( itr->index() + o, itr->value() );
00522         }
00523         else {
00524           lambda->set_ele( itr->index() + o - m_in, itr->value() );
00525         }
00526       }
00527     }
00528   }
00529   // obj_d (This could be updated within QPSchur in the future)
00530   if(obj_d) {
00531     // obj_d = g'*d + 1/2 * d' * G * g
00532     *obj_d = AbstractLinAlgPack::dot(g,*d)
00533       + 0.5 * AbstractLinAlgPack::transVtMtV(*d,G,BLAS_Cpp::no_trans,*d);
00534   }
00535   // Ed
00536   if(Ed && E) {
00537     switch(constraints_->inequality_pick_policy()) {
00538       case constr_t::ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY:
00539         if(solve_returned == QPSchur::OPTIMAL_SOLUTION)
00540           break; // Ed already computed (see ConstraintsRelaxedStd::pick_violated())
00541       case constr_t::ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY:
00542         break; // Ed already computed (see ConstraintsRelaxedStd::pick_violated())
00543       default:
00544         // We need to compute Ed
00545         LinAlgOpPack::V_MtV( Ed, *E, trans_E, *d );
00546     }
00547   }
00548   // Fd (This could be updated within ConstraintsRelaxedStd in the future)
00549   if(Fd) {
00550     LinAlgOpPack::V_MtV( Fd, *F, trans_F, *d );
00551   }
00552   // Set the QP statistics
00553   QPSolverStats::ESolutionType solution_type;
00554   QPSolverStats::EConvexity convexity = QPSolverStats::CONVEX;
00555   switch( solve_returned ) {
00556     case QPSchur::OPTIMAL_SOLUTION:
00557       solution_type = QPSolverStats::OPTIMAL_SOLUTION;
00558       break;
00559     case QPSchur::MAX_ITER_EXCEEDED:
00560       solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
00561       break;
00562     case QPSchur::MAX_RUNTIME_EXEEDED_FAIL:
00563       solution_type = QPSolverStats::SUBOPTIMAL_POINT;
00564       break;
00565     case QPSchur::MAX_RUNTIME_EXEEDED_DUAL_FEAS:
00566       solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
00567       break;
00568     case QPSchur::MAX_ALLOWED_STORAGE_EXCEEDED:
00569       solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
00570       break;
00571     case QPSchur::INFEASIBLE_CONSTRAINTS:
00572     case QPSchur::NONCONVEX_QP:
00573       convexity = QPSolverStats::NONCONVEX;
00574     case QPSchur::DUAL_INFEASIBILITY:
00575     case QPSchur::SUBOPTIMAL_POINT:
00576       solution_type = QPSolverStats::SUBOPTIMAL_POINT;
00577       break;
00578     default:
00579       TEUCHOS_TEST_FOR_EXCEPT(true);
00580   }
00581   qp_stats_.set_stats(
00582     solution_type,convexity,qp_iter,num_adds,num_drops
00583     , num_act_change > 0 || n_X > 1, *eta > 0.0 );
00584 
00585   return qp_stats_.solution_type();
00586 }
00587 
00588 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends