MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_TangentialStepWithInequStd_Step.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 <math.h>
00043 
00044 #include <ostream>
00045 #include <sstream>
00046 
00047 #include "MoochoPack_TangentialStepWithInequStd_Step.hpp"
00048 #include "MoochoPack_moocho_algo_conversion.hpp"
00049 #include "MoochoPack_Exceptions.hpp"
00050 #include "IterationPack_print_algorithm_step.hpp"
00051 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
00052 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00055 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00056 #include "Teuchos_dyn_cast.hpp"
00057 
00058 namespace {
00059 template< class T >
00060 inline
00061 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00062 template< class T >
00063 inline
00064 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00065 } // end namespace
00066 
00067 namespace MoochoPack {
00068 
00069 TangentialStepWithInequStd_Step::TangentialStepWithInequStd_Step(
00070   const qp_solver_ptr_t       &qp_solver
00071   ,const qp_tester_ptr_t      &qp_tester
00072   ,value_type                 warm_start_frac
00073   ,EQPTesting                 qp_testing
00074   ,bool                       primal_feasible_point_error
00075   ,bool                       dual_feasible_point_error
00076   )
00077   :qp_solver_(qp_solver)
00078   ,qp_tester_(qp_tester)
00079   ,warm_start_frac_(warm_start_frac)
00080   ,qp_testing_(qp_testing)
00081   ,primal_feasible_point_error_(primal_feasible_point_error)
00082   ,dual_feasible_point_error_(dual_feasible_point_error)
00083   ,dl_iq_(dl_name)
00084   ,du_iq_(du_name)
00085 {}
00086 
00087 bool TangentialStepWithInequStd_Step::do_step(
00088   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00089   ,poss_type assoc_step_poss
00090   )
00091 {
00092 
00093   using Teuchos::RCP;
00094   using Teuchos::dyn_cast;
00095   using ::fabs;
00096   using LinAlgOpPack::Vt_S;
00097   using LinAlgOpPack::V_VpV;
00098   using LinAlgOpPack::V_VmV;
00099   using LinAlgOpPack::Vp_StV;
00100   using LinAlgOpPack::Vp_V;
00101   using LinAlgOpPack::V_StV;
00102   using LinAlgOpPack::V_MtV;
00103 //  using ConstrainedOptPack::min_abs;
00104   using AbstractLinAlgPack::max_near_feas_step;
00105   typedef VectorMutable::vec_mut_ptr_t   vec_mut_ptr_t;
00106 
00107   NLPAlgo &algo = rsqp_algo(_algo);
00108   NLPAlgoState &s = algo.rsqp_state();
00109   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00110   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00111   std::ostream &out = algo.track().journal_out();
00112   //const bool check_results = algo.algo_cntr().check_results();
00113 
00114   // print step header.
00115   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00116     using IterationPack::print_algorithm_step;
00117     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00118   }
00119 
00120   // problem dimensions
00121   const size_type
00122     //n  = algo.nlp().n(),
00123     m  = algo.nlp().m(),
00124     r  = s.equ_decomp().size();
00125 
00126   // Get the iteration quantity container objects
00127   IterQuantityAccess<value_type>
00128     &alpha_iq = s.alpha(),
00129     &zeta_iq  = s.zeta(),
00130     &eta_iq   = s.eta();
00131   IterQuantityAccess<VectorMutable>
00132     &dl_iq      = dl_iq_(s),
00133     &du_iq      = du_iq_(s),
00134     &nu_iq      = s.nu(),
00135     *c_iq       = m > 0  ? &s.c()       : NULL,
00136     *lambda_iq  = m > 0  ? &s.lambda()  : NULL,
00137     &rGf_iq     = s.rGf(),
00138     &w_iq       = s.w(),
00139     &qp_grad_iq = s.qp_grad(),
00140     &py_iq      = s.py(),
00141     &pz_iq      = s.pz(),
00142     &Ypy_iq     = s.Ypy(),
00143     &Zpz_iq     = s.Zpz();
00144   IterQuantityAccess<MatrixOp>
00145     &Z_iq   = s.Z(),
00146     //*Uz_iq  = (m > r)  ? &s.Uz() : NULL,
00147     *Uy_iq  = (m > r)  ? &s.Uy() : NULL;
00148   IterQuantityAccess<MatrixSymOp>
00149     &rHL_iq = s.rHL();
00150   IterQuantityAccess<ActSetStats>
00151     &act_set_stats_iq = act_set_stats_(s);
00152   
00153   // Accessed/modified/updated (just some)
00154   VectorMutable  *Ypy_k = (m ? &Ypy_iq.get_k(0) : NULL);
00155   const MatrixOp  &Z_k   = Z_iq.get_k(0);
00156   VectorMutable  &pz_k  = pz_iq.set_k(0);
00157   VectorMutable  &Zpz_k = Zpz_iq.set_k(0);
00158 
00159   // Comupte qp_grad which is an approximation to rGf + Z'*HL*Y*py
00160 
00161   // qp_grad = rGf
00162   VectorMutable
00163     &qp_grad_k = ( qp_grad_iq.set_k(0) = rGf_iq.get_k(0) );
00164 
00165   // qp_grad += zeta * w
00166   if( w_iq.updated_k(0) ) {
00167     if(zeta_iq.updated_k(0))
00168       Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
00169     else
00170       Vp_V( &qp_grad_k, w_iq.get_k(0) );
00171   }
00172 
00173   //
00174   // Set the bounds for:
00175   //
00176   //   dl <= Z*pz + Y*py <= du  ->  dl - Ypy <= Z*pz <= du - Ypz
00177 
00178   vec_mut_ptr_t
00179     bl = s.space_x().create_member(),
00180     bu = s.space_x().create_member();
00181 
00182   if(m) {
00183     // bl = dl_k - Ypy_k
00184     V_VmV( bl.get(), dl_iq.get_k(0), *Ypy_k );
00185     // bu = du_k - Ypy_k
00186     V_VmV( bu.get(), du_iq.get_k(0), *Ypy_k );
00187   }
00188   else {
00189     *bl = dl_iq.get_k(0);
00190     *bu = du_iq.get_k(0);
00191   }
00192 
00193   // Print out the QP bounds for the constraints
00194   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00195     out << "\nqp_grad_k = \n" << qp_grad_k;
00196   }
00197   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00198     out << "\nbl = \n" << *bl;
00199     out << "\nbu = \n" << *bu;
00200   }
00201 
00202   //
00203   // Determine if we should perform a warm start or not.
00204   //
00205   bool do_warm_start = false;
00206   if( act_set_stats_iq.updated_k(-1) ) {
00207     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00208       out << "\nDetermining if the QP should use a warm start ...\n";
00209     }
00210     // We need to see if we should preform a warm start for the next iteration
00211     ActSetStats &stats = act_set_stats_iq.get_k(-1);
00212     const size_type
00213       num_active = stats.num_active(),
00214       num_adds   = stats.num_adds(),
00215       num_drops  = stats.num_drops();
00216     const value_type
00217       frac_same
00218       = ( num_adds == ActSetStats::NOT_KNOWN || num_active == 0
00219         ? 0.0
00220         : my_max(((double)(num_active)-num_adds-num_drops) / num_active, 0.0 ) );
00221     do_warm_start = ( num_active > 0 && frac_same >= warm_start_frac() );
00222     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00223       out << "\nnum_active = " << num_active;
00224       if( num_active ) {
00225         out << "\nmax(num_active-num_adds-num_drops,0)/(num_active) = "
00226           << "max("<<num_active<<"-"<<num_adds<<"-"<<num_drops<<",0)/("<<num_active<<") = "
00227           << frac_same;
00228         if( do_warm_start )
00229           out << " >= ";
00230         else
00231           out << " < ";
00232         out << "warm_start_frac = " << warm_start_frac();
00233       }
00234       if( do_warm_start )
00235         out << "\nUse a warm start this time!\n";
00236       else
00237         out << "\nDon't use a warm start this time!\n";
00238     }
00239   }
00240 
00241   // Use active set from last iteration as an estimate for current active set
00242   // if we are to use a warm start.
00243   // 
00244   // ToDo: If the selection of dependent and independent variables changes
00245   // then you will have to adjust this or not perform a warm start at all!
00246   if( do_warm_start ) {
00247     nu_iq.set_k(0,-1);
00248   }
00249   else {
00250     nu_iq.set_k(0) = 0.0; // No guess of the active set
00251   }
00252   VectorMutable  &nu_k  = nu_iq.get_k(0);
00253 
00254   //
00255   // Setup the reduced QP subproblem
00256   //
00257   // The call to the QP is setup for the more flexible call to the QPSolverRelaxed
00258   // interface to deal with the three independent variabilities: using simple
00259   // bounds for pz or not, general inequalities included or not, and extra equality
00260   // constraints included or not.
00261   // If this method of calling the QP solver were not used then 4 separate
00262   // calls to solve_qp(...) would have to be included to handle the four possible
00263   // QP formulations.
00264   //
00265 
00266   // The numeric arguments for the QP solver (in the nomenclatrue of QPSolverRelaxed)
00267 
00268   const value_type  qp_bnd_inf = NLP::infinite_bound();
00269 
00270   const Vector            &qp_g       = qp_grad_k;
00271   const MatrixSymOp       &qp_G       = rHL_iq.get_k(0);
00272   const value_type        qp_etaL     = 0.0;
00273   vec_mut_ptr_t           qp_dL       = Teuchos::null;
00274   vec_mut_ptr_t           qp_dU       = Teuchos::null;
00275   Teuchos::RCP<const MatrixOp>
00276                           qp_E        = Teuchos::null;
00277   BLAS_Cpp::Transp        qp_trans_E  = BLAS_Cpp::no_trans;
00278   vec_mut_ptr_t           qp_b        = Teuchos::null;
00279   vec_mut_ptr_t           qp_eL       = Teuchos::null;
00280   vec_mut_ptr_t           qp_eU       = Teuchos::null;
00281   Teuchos::RCP<const MatrixOp>
00282                           qp_F        = Teuchos::null;
00283   BLAS_Cpp::Transp        qp_trans_F  = BLAS_Cpp::no_trans;
00284   vec_mut_ptr_t           qp_f        = Teuchos::null;
00285   value_type              qp_eta      = 0.0;
00286   VectorMutable           &qp_d       = pz_k;  // pz_k will be updated directly!
00287   vec_mut_ptr_t           qp_nu       = Teuchos::null;
00288   vec_mut_ptr_t           qp_mu       = Teuchos::null;
00289   vec_mut_ptr_t           qp_Ed       = Teuchos::null;
00290   vec_mut_ptr_t           qp_lambda   = Teuchos::null;
00291 
00292   //
00293   // Determine if we can use simple bounds on pz.
00294   // 
00295   // If we have a variable-reduction null-space matrix
00296   // (with any choice for Y) then:
00297   // 
00298   // d = Z*pz + (1-eta) * Y*py
00299   // 
00300   // [ d(var_dep)   ]  = [ D ] * pz  + (1-eta) * [ Ypy(var_dep)   ]
00301   // [ d(var_indep) ]    [ I ]                   [ Ypy(var_indep) ]
00302   // 
00303   // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ypy(var_indep) ==
00304   // 0.0 and in this case the bounds on d(var_indep) become simple bounds on
00305   // pz even with the relaxation.  Also, if dl(var_dep) and du(var_dep) are
00306   // unbounded, then we can also use simple bounds since we don't need the
00307   // relaxation and we can set eta=0. In this case we just have to subtract
00308   // from the upper and lower bounds on pz!
00309   // 
00310   // Otherwise, we can not use simple variable bounds and implement the
00311   // relaxation properly.
00312   // 
00313 
00314   const MatrixIdentConcat
00315     *Zvr = dynamic_cast<const MatrixIdentConcat*>( &Z_k );
00316   const Range1D
00317     var_dep   = Zvr ? Zvr->D_rng() : Range1D::Invalid,
00318     var_indep = Zvr ? Zvr->I_rng() : Range1D();
00319 
00320   RCP<Vector> Ypy_indep;
00321   const value_type
00322     Ypy_indep_norm_inf
00323     = ( m ? (Ypy_indep=Ypy_k->sub_view(var_indep))->norm_inf() : 0.0);
00324 
00325   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00326     out
00327       << "\nDetermine if we can use simple bounds on pz ...\n"
00328       << "    m = " << m << std::endl
00329       << "    dynamic_cast<const MatrixIdentConcat*>(&Z_k) = " << Zvr << std::endl
00330       << "    ||Ypy_k(var_indep)||inf = " << Ypy_indep_norm_inf << std::endl;
00331 
00332   const bool
00333     bounded_var_dep
00334     = (
00335       m > 0
00336       &&
00337       num_bounded( *bl->sub_view(var_dep), *bu->sub_view(var_dep), qp_bnd_inf )
00338       );
00339 
00340   const bool
00341     use_simple_pz_bounds
00342     = (
00343       m == 0
00344       ||
00345       ( Zvr != NULL && ( Ypy_indep_norm_inf == 0.0 || bounded_var_dep == 0 ) )
00346       );
00347 
00348   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00349     out
00350       << (use_simple_pz_bounds
00351           ? "\nUsing simple bounds on pz ...\n"
00352           : "\nUsing bounds on full Z*pz ...\n")
00353       << (bounded_var_dep
00354           ? "\nThere are finite bounds on dependent variables.  Adding extra inequality constrints for D*pz ...\n" 
00355           : "\nThere are no finite bounds on dependent variables.  There will be no extra inequality constraints added on D*pz ...\n" ) ;
00356 
00357   if( use_simple_pz_bounds ) {
00358     // Set simple bound constraints on pz
00359     qp_dL = bl->sub_view(var_indep);
00360     qp_dU = bu->sub_view(var_indep);
00361     qp_nu = nu_k.sub_view(var_indep); // nu_k(var_indep) will be updated directly!
00362     if( m && bounded_var_dep ) {
00363       // Set general inequality constraints for D*pz
00364       qp_E   = Teuchos::rcp(&Zvr->D(),false);
00365       qp_b   = Ypy_k->sub_view(var_dep);
00366       qp_eL  = bl->sub_view(var_dep);
00367       qp_eU  = bu->sub_view(var_dep);
00368       qp_mu  = nu_k.sub_view(var_dep);  // nu_k(var_dep) will be updated directly!
00369       qp_Ed  = Zpz_k.sub_view(var_dep); // Zpz_k(var_dep) will be updated directly!
00370     }
00371     else {
00372       // Leave these as NULL since there is no extra general inequality constraints
00373     }
00374   }
00375   else if( !use_simple_pz_bounds ) { // ToDo: Leave out parts for unbounded dependent variables!
00376     // There are no simple bounds! (leave qp_dL, qp_dU and qp_nu as null)
00377     // Set general inequality constraints for Z*pz
00378     qp_E   = Teuchos::rcp(&Z_k,false);
00379     qp_b   = Teuchos::rcp(Ypy_k,false);
00380     qp_eL  = bl;
00381     qp_eU  = bu;
00382     qp_mu  = Teuchos::rcp(&nu_k,false);
00383     qp_Ed  = Teuchos::rcp(&Zpz_k,false); // Zpz_k will be updated directly!
00384   }
00385   else {
00386     TEUCHOS_TEST_FOR_EXCEPT(true);
00387   }
00388 
00389   // Set the general equality constriants (if they exist)
00390   Range1D equ_undecomp = s.equ_undecomp();
00391   if( m > r && m > 0 ) {
00392     // qp_f = Uy_k * py_k + c_k(equ_undecomp)
00393     qp_f = s.space_c().sub_space(equ_undecomp)->create_member();
00394     V_MtV( qp_f.get(), Uy_iq->get_k(0), BLAS_Cpp::no_trans, py_iq.get_k(0) );
00395     Vp_V( qp_f.get(), *c_iq->get_k(0).sub_view(equ_undecomp) );
00396     // Must resize for the undecomposed constriants if it has not already been
00397     qp_F       = Teuchos::rcp(&Uy_iq->get_k(0),false);
00398     qp_lambda  = lambda_iq->set_k(0).sub_view(equ_undecomp); // lambda_k(equ_undecomp), will be updated directly!
00399   }
00400 
00401   // Setup the rest of the arguments
00402 
00403   QPSolverRelaxed::EOutputLevel
00404     qp_olevel;
00405   switch( olevel ) {
00406     case PRINT_NOTHING:
00407       qp_olevel = QPSolverRelaxed::PRINT_NONE;
00408       break;
00409     case PRINT_BASIC_ALGORITHM_INFO:
00410       qp_olevel = QPSolverRelaxed::PRINT_NONE;
00411       break;
00412     case PRINT_ALGORITHM_STEPS:
00413       qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00414       break;
00415     case PRINT_ACTIVE_SET:
00416       qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
00417       break;
00418     case PRINT_VECTORS:
00419       qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
00420       break;
00421     case PRINT_ITERATION_QUANTITIES:
00422       qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
00423       break;
00424     default:
00425       TEUCHOS_TEST_FOR_EXCEPT(true);
00426   }
00427   // ToDo: Set print options so that only vectors matrices etc
00428   // are only printed in the null space
00429 
00430   //
00431   // Solve the QP
00432   // 
00433   qp_solver().infinite_bound(qp_bnd_inf);
00434   const QPSolverStats::ESolutionType
00435     solution_type =
00436     qp_solver().solve_qp(
00437       int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00438       ,qp_olevel
00439       ,( algo.algo_cntr().check_results()
00440          ? QPSolverRelaxed::RUN_TESTS :  QPSolverRelaxed::NO_TESTS )
00441       ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
00442       ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
00443       ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
00444       ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
00445       ,NULL // obj_d
00446       ,&qp_eta, &qp_d
00447       ,qp_nu.get()
00448       ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
00449       ,qp_F.get() ? qp_lambda.get() : NULL
00450       ,NULL // qp_Fd
00451       );
00452 
00453   //
00454   // Check the optimality conditions for the QP
00455   //
00456   std::ostringstream omsg;
00457   bool throw_qp_failure = false;
00458   if(   qp_testing() == QP_TEST
00459     || ( qp_testing() == QP_TEST_DEFAULT && algo.algo_cntr().check_results() )  )
00460   {
00461     if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
00462       out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
00463     }
00464     if(!qp_tester().check_optimality_conditions(
00465       solution_type,qp_solver().infinite_bound()
00466       ,int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00467       ,int(olevel) >= int(PRINT_VECTORS) ? true : false
00468       ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
00469       ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
00470       ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
00471       ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
00472       ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
00473       ,NULL // obj_d
00474       ,&qp_eta, &qp_d
00475       ,qp_nu.get()
00476       ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
00477       ,qp_F.get() ? qp_lambda.get() : NULL
00478       ,NULL // qp_Fd
00479       ))
00480     {
00481       omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
00482       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00483         out << omsg.str();
00484       }
00485       throw_qp_failure = true;
00486     }
00487   }
00488 
00489   //
00490   // Set the solution
00491   //
00492   if( !use_simple_pz_bounds ) {
00493     // Everything is already updated!
00494   }
00495   else if( use_simple_pz_bounds ) {
00496     // Just have to set Zpz_k(var_indep) = pz_k
00497     *Zpz_k.sub_view(var_indep) = pz_k;
00498     if( m && !bounded_var_dep ) {
00499       // Must compute Zpz_k(var_dep) = D*pz
00500       LinAlgOpPack::V_MtV( &*Zpz_k.sub_view(var_dep), Zvr->D(), BLAS_Cpp::no_trans, pz_k );
00501       // ToDo: Remove the compuation of Zpz here unless you must
00502     }
00503   }
00504   else {
00505     TEUCHOS_TEST_FOR_EXCEPT(true);
00506   }
00507 
00508   // Set the solution statistics
00509   qp_solver_stats_(s).set_k(0) = qp_solver().get_qp_stats();
00510 
00511   // Cut back Ypy_k = (1-eta) * Ypy_k
00512   const value_type eps = std::numeric_limits<value_type>::epsilon();
00513   if( fabs(qp_eta - 0.0) > eps ) {
00514     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00515       out
00516         << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<").  Cutting back Ypy_k = (1.0 - eta)*Ypy_k  ...\n";
00517     }
00518     Vt_S( Ypy_k, 1.0 - qp_eta );
00519   }
00520 
00521   // eta_k
00522   eta_iq.set_k(0) = qp_eta;
00523 
00524   //
00525   // Modify the solution if we have to!
00526   // 
00527   switch(solution_type) {
00528     case QPSolverStats::OPTIMAL_SOLUTION:
00529       break;  // we are good!
00530     case QPSolverStats::PRIMAL_FEASIBLE_POINT:
00531     {
00532       omsg
00533         << "\n*** Alert! the returned QP solution is PRIMAL_FEASIBLE_POINT but not optimal!\n";
00534       if( primal_feasible_point_error() )
00535         omsg
00536           << "\n*** primal_feasible_point_error == true, this is an error!\n";
00537       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00538         out << omsg.str();
00539       }
00540       throw_qp_failure = primal_feasible_point_error();
00541       break;
00542     } 
00543     case QPSolverStats::DUAL_FEASIBLE_POINT:
00544     {
00545       omsg
00546         << "\n*** Alert! the returned QP solution is DUAL_FEASIBLE_POINT"
00547         << "\n*** but not optimal so we cut back the step ...\n";
00548       if( dual_feasible_point_error() )
00549         omsg
00550           << "\n*** dual_feasible_point_error == true, this is an error!\n";
00551       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00552         out << omsg.str();
00553       }
00554       // Cut back the step to fit in the bounds
00555       // 
00556       // dl <= u*(Ypy_k+Zpz_k) <= du
00557       //
00558       vec_mut_ptr_t
00559         zero  = s.space_x().create_member(0.0),
00560         d_tmp = s.space_x().create_member();
00561       V_VpV( d_tmp.get(), *Ypy_k, Zpz_k );
00562       const std::pair<value_type,value_type>
00563         u_steps = max_near_feas_step( *zero, *d_tmp, dl_iq.get_k(0), du_iq.get_k(0), 0.0 );
00564       const value_type
00565         u = my_min( u_steps.first, 1.0 ); // largest positive step size
00566       alpha_iq.set_k(0) = u;
00567       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00568         out << "\nFinding u s.t. dl <= u*(Ypy_k+Zpz_k) <= du\n"
00569           << "max step length u = " << u << std::endl
00570           << "alpha_k = u = " << alpha_iq.get_k(0) << std::endl;
00571       }
00572       throw_qp_failure = dual_feasible_point_error();
00573       break;
00574     } 
00575     case QPSolverStats::SUBOPTIMAL_POINT:
00576     {
00577       omsg
00578         << "\n*** Alert!, the returned QP solution is SUBOPTIMAL_POINT!\n";
00579       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00580         out << omsg.str();
00581       }
00582       throw_qp_failure = true;
00583       break;
00584     }
00585     default:
00586       TEUCHOS_TEST_FOR_EXCEPT(true);  // should not happen!
00587   }
00588 
00589   //
00590   // Output the final solution!
00591   //
00592   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00593     out << "\n||pz_k||inf    = " << s.pz().get_k(0).norm_inf()
00594       << "\nnu_k.nz()      = " << s.nu().get_k(0).nz()
00595       << "\nmax(|nu_k(i)|) = " << s.nu().get_k(0).norm_inf()
00596 //      << "\nmin(|nu_k(i)|) = " << min_abs( s.nu().get_k(0)() )
00597       ;
00598     if( m > r ) out << "\n||lambda_k(undecomp)||inf = " << s.lambda().get_k(0).norm_inf();
00599     out << "\n||Zpz_k||2     = " << s.Zpz().get_k(0).norm_2()
00600       ;
00601     if(qp_eta > 0.0) out << "\n||Ypy||2 = " << s.Ypy().get_k(0).norm_2();
00602     out << std::endl;
00603   }
00604 
00605   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00606     out << "\npz_k = \n" << s.pz().get_k(0);
00607     if(var_indep.size())
00608       out << "\nnu_k(var_indep) = \n" << *s.nu().get_k(0).sub_view(var_indep);
00609   }
00610 
00611   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00612     if(var_indep.size())
00613       out << "\nZpz(var_indep)_k = \n" << *s.Zpz().get_k(0).sub_view(var_indep);
00614     out << std::endl;
00615   }
00616 
00617   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00618     if(var_dep.size())
00619       out << "\nZpz(var_dep)_k = \n" << *s.Zpz().get_k(0).sub_view(var_dep);
00620     out << "\nZpz_k = \n" << s.Zpz().get_k(0);
00621     out << std::endl;
00622   }
00623 
00624   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00625     out << "\nnu_k = \n" << s.nu().get_k(0);
00626     if(var_dep.size())
00627       out << "\nnu_k(var_dep) = \n" << *s.nu().get_k(0).sub_view(var_dep);
00628     if( m > r )
00629       out << "\nlambda_k(equ_undecomp) = \n" << *s.lambda().get_k(0).sub_view(equ_undecomp);
00630     if(qp_eta > 0.0) out << "\nYpy = \n" << s.Ypy().get_k(0);
00631   }
00632 
00633   if( qp_eta == 1.0 ) {
00634     omsg
00635       << "TangentialStepWithInequStd_Step::do_step(...) : Error, a QP relaxation parameter\n"
00636       << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
00637       << "that the NLP's constraints are infeasible\n"
00638       << "Throwing an InfeasibleConstraints exception!\n";
00639     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00640       out << omsg.str();
00641     }
00642     throw InfeasibleConstraints(omsg.str());
00643   }
00644 
00645   if( throw_qp_failure )
00646     throw QPFailure( omsg.str(), qp_solver().get_qp_stats() );
00647 
00648   return true;
00649 }
00650 
00651 void TangentialStepWithInequStd_Step::print_step(
00652   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00653   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00654   ) const
00655 {
00656   out
00657     << L << "*** Calculate the null-space step by solving a constrained QP\n"
00658     << L << "qp_grad_k = rGf_k\n"
00659     << L << "if w_k is updated then set qp_grad_k = qp_grad_k + zeta_k * w_k\n"
00660     << L << "bl = dl_k - Ypy_k\n"
00661     << L << "bu = du_k - Ypy_k\n"
00662     << L << "etaL = 0.0\n"
00663     << L << "*** Determine if we can use simple bounds on pz or not\n"
00664     << L << "if num_bounded(bl_k(var_dep),bu_k(var_dep)) > 0 then\n"
00665     << L << "  bounded_var_dep = true\n"
00666     << L << "else\n"
00667     << L << "  bounded_var_dep = false\n"
00668     << L << "end\n"
00669     << L << "if( m==0\n"
00670     << L << "    or\n"
00671     << L << "    ( Z_k is a variable reduction null space matrix\n"
00672     << L << "      and\n"
00673     << L << "      ( ||Ypy_k(var_indep)||inf == 0 or bounded_var_dep==false ) )\n"
00674     << L << "  ) then\n"
00675     << L << "  use_simple_pz_bounds = true\n"
00676     << L << "else\n"
00677     << L << "  use_simple_pz_bounds = false\n"
00678     << L << "end\n"
00679     << L << "*** Setup QP arguments\n"
00680     << L << "qp_g = qp_grad_k\n"
00681     << L << "qp_G = rHL_k\n"
00682     << L << "if (use_simple_pz_bounds == true) then\n"
00683     << L << "  qp_dL = bl(var_indep), qp_dU = bu(var_indep))\n"
00684     << L << "  if (m > 0) then\n"
00685     << L << "    qp_E  = Z_k.D,       qp_b  = Ypy_k(var_dep)\n"
00686     << L << "    qp_eL = bl(var_dep), qp_eU = bu(var_dep)\n"
00687     << L << "  end\n"
00688     << L << "elseif (use_simple_pz_bounds == false) then\n"
00689     << L << "  qp_dL = -inf,  qp_dU = +inf\n"
00690     << L << "  qp_E  = Z_k,   qp_b  = Ypy_k\n"
00691     << L << "  qp_eL = bl,    qp_eU = bu\n"
00692     << L << "end\n"
00693     << L << "if (m > r) then\n"
00694     << L << "  qp_F  = V_k,     qp_f  = Uy_k*py_k + c_k(equ_undecomp)\n"
00695     << L << "else\n"
00696     << L << "  qp_F  = empty,   qp_f  = empty\n"
00697     << L << "end\n"
00698     << L << "Given active-set statistics (act_set_stats_km1)\n"
00699     << L << "  frac_same = max(num_active-num_adds-num_drops,0)/(num_active)\n"
00700     << L << "Use a warm start when frac_same >= warm_start_frac\n"
00701     << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
00702     << L << ",qp_nu, qp_mu and qp_lambda (" << typeName(qp_solver()) << "):\n"
00703     << L << "  min  qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
00704     << L << "  qp_d <: R^(n-r)\n"
00705     << L << "  s.t.\n"
00706     << L << "       etaL  <=  eta\n"
00707     << L << "       qp_dL <= qp_d                         <= qp_dU   [qp_nu]\n"
00708     << L << "       qp_eL <= qp_E * qp_d + (1-eta)*qp_b   <= qp_eU   [qp_mu]\n"
00709     << L << "                qp_F * d_qp + (1-eta) * qp_f  = 0       [qp_lambda]\n"
00710     << L << "if (qp_testing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
00711     << L << "and check_results==true) then\n"
00712     << L << "  Check the optimality conditions of the above QP\n"
00713     << L << "  if the optimality conditions do not check out then\n"
00714     << L << "    set throw_qp_failure = true\n"
00715     << L << "  end\n"
00716     << L << "end\n"
00717     << L << "*** Set the solution to the QP subproblem\n"
00718     << L << "pz_k  = qp_d\n"
00719     << L << "eta_k = qp_eta\n"
00720     << L << "if (use_simple_pz_bounds == true) then\n"
00721     << L << "  nu_k(var_dep)  = qp_mu,  nu_k(var_indep)  = qp_nu\n"
00722     << L << "  Zpz_k(var_dep) = qp_Ed,  Zpz_k(var_indep) = pz_k\n"
00723     << L << "elseif (use_simple_pz_bounds == false) then\n"
00724     << L << "  nu_k  = qp_mu\n"
00725     << L << "  Zpz_k = qp_Ed\n"
00726     << L << "end\n"
00727     << L << "if m > r then\n"
00728     << L << "  lambda_k(equ_undecomp) = qp_lambda\n"
00729     << L << "end\n"
00730     << L << "if (eta_k > 0) then set Ypy_k = (1-eta_k) * Ypy_k\n"
00731     << L << "if QP solution is suboptimal then\n"
00732     << L << "  throw_qp_failure = true\n"
00733     << L << "elseif QP solution is primal feasible (not optimal) then\n"
00734     << L << "  throw_qp_failure = primal_feasible_point_error\n"
00735     << L << "elseif QP solution is dual feasible (not optimal) then\n"
00736     << L << "  find max u s.t.\n"
00737     << L << "    dl_k <= u*(Ypy_k+Zpz_k) <= du_k\n"
00738     << L << "  alpha_k = u\n"
00739     << L << "  throw_qp_failure = dual_feasible_point_error\n"
00740     << L << "end\n"
00741     << L << "if (eta_k == 1.0) then\n"
00742     << L << "  The constraints are infeasible!\n"
00743     << L << "  throw InfeasibleConstraints(...)\n"
00744     << L << "end\n"
00745     << L << "if (throw_qp_failure == true) then\n"
00746     << L << "  throw QPFailure(...)\n"
00747     << L << "end\n"
00748     ;
00749 }
00750 
00751 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends