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

Generated on Thu Sep 18 12:34:16 2008 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by doxygen 1.3.9.1