MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_FeasibilityStepReducedStd_Strategy.cpp
Go to the documentation of this file.
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 "MoochoPack_FeasibilityStepReducedStd_Strategy.hpp"
00043 #include "MoochoPack_NLPAlgo.hpp"
00044 #include "MoochoPack_NLPAlgoState.hpp"
00045 #include "MoochoPack_Exceptions.hpp"
00046 #include "Teuchos_dyn_cast.hpp"
00047 
00048 namespace MoochoPack {
00049 
00050 FeasibilityStepReducedStd_Strategy::FeasibilityStepReducedStd_Strategy(
00051   const quasi_range_space_step_ptr_t   &quasi_range_space_step
00052   ,const qp_solver_ptr_t               &qp_solver
00053   ,const qp_tester_ptr_t               &qp_tester
00054   ,EQPObjective                        qp_objective
00055   ,EQPTesting                          qp_testing
00056   )
00057   :quasi_range_space_step_(quasi_range_space_step)
00058   ,qp_solver_(qp_solver)
00059   ,qp_tester_(qp_tester)
00060   ,qp_objective_(qp_objective)
00061   ,qp_testing_(qp_testing)
00062   ,dl_iq_(dl_name)
00063   ,du_iq_(du_name)
00064   ,current_k_(-1)
00065 {}
00066 
00067 bool FeasibilityStepReducedStd_Strategy::compute_feasibility_step(
00068   std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
00069   ,const Vector& xo, const Vector& c_xo, VectorMutable* w
00070     )
00071 {
00072   using Teuchos::dyn_cast;
00073 
00074 /* Todo: UPdate below code!
00075 
00076   // problem dimensions
00077   const size_type
00078     n = algo->nlp().n(),
00079     m = algo->nlp().m(),
00080     r  = s->equ_decomp().size();
00081 
00082   // Compute the quasi-range space step Ywy
00083   Workspace<value_type> Ywy_ws(wss,xo.size());
00084   DVectorSlice                Ywy(&Ywy_ws[0],Ywy_ws.size());
00085   if(!quasi_range_space_step().solve_quasi_range_space_step(
00086     out,olevel,algo,s,xo,c_xo,&Ywy ))
00087     return false;
00088 
00089   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00090     out << "\n||Ywy||2     = " << DenseLinAlgPack::norm_2(Ywy);
00091     out << std::endl;
00092   }
00093   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00094     out << "\nYwy = \n" << Ywy;
00095   }
00096 
00097   //
00098   // Set the bounds on the null space QP subproblem:
00099   //
00100   // d_bounds_k.l <= (xo - x_k) + (1-eta) * Ywy + Z*wz <= d_bounds_k.u
00101   // =>
00102   // bl <= Z*wz - eta * Ywy <= bu
00103   //
00104   // where:   bl = d_bounds_k.l - (xo - x_k) - Ywy
00105   //          bu = d_bounds_k.u - (xo - x_k) - Ywy
00106   //
00107   // Above, fix the variables that are at an active bound as equalities
00108   // to maintain the same active-set.
00109   //
00110   const SparseBounds
00111     &d_bounds = d_bounds_(*s).get_k(0);
00112   const SpVectorSlice
00113     &dl = d_bounds.l,
00114     &du = d_bounds.u;
00115   const DVector
00116     &x_k = s->x().get_k(0).v();
00117   const SpVector
00118     &nu_k = s->nu().get_k(0);
00119   TEUCHOS_TEST_FOR_EXCEPT( !( nu_k.is_sorted() ) );
00120   SpVector bl(n,n), bu(n,n);
00121   sparse_bounds_itr
00122     d_bnds_itr(dl.begin(),dl.end(),dl.offset(),du.begin(),du.end(),du.offset());
00123   SpVector::const_iterator
00124     nu_itr     = nu_k.begin(),
00125     nu_end     = nu_k.end();
00126   for( ; !d_bnds_itr.at_end(); ++d_bnds_itr ) {
00127     typedef SpVectorSlice::element_type ele_t;
00128     const size_type i = d_bnds_itr.indice();
00129     while( nu_itr != nu_end && nu_itr->indice() + nu_k.offset() < i )
00130       ++nu_itr;
00131     if( nu_itr != nu_end && nu_itr->indice() + nu_k.offset() == i ) {
00132       const value_type
00133         act_bnd = nu_itr->value() > 0.0 ? d_bnds_itr.ubound() : d_bnds_itr.lbound();
00134       bl.add_element(ele_t( i, act_bnd - xo(i) + x_k(i) - Ywy(i) ));
00135       bu.add_element(ele_t( i, act_bnd - xo(i) + x_k(i) - Ywy(i) ));
00136     }
00137     else {
00138       if( d_bnds_itr.lbound() != -d_bnds_itr.big_bnd() )
00139         bl.add_element(ele_t(i,d_bnds_itr.lbound()  - xo(i) + x_k(i) - Ywy(i) ));
00140       if( d_bnds_itr.ubound() != +d_bnds_itr.big_bnd() )
00141         bu.add_element(ele_t(i, d_bnds_itr.ubound() - xo(i) + x_k(i) - Ywy(i) ));
00142     }
00143   }
00144   bl.assume_sorted(true);
00145   bu.assume_sorted(true);
00146   //
00147   // Setup the objective function for the null space QP subproblem
00148   //
00149   // 
00150   // OBJ_MIN_FULL_STEP
00151   //    min 1/2 * (Y*wy + Z*wz)'*(Y*wy + Z*wz) = 1/2*wy'*Y'*Y*wy + (Z'*Y*wy)'*wz + 1/2*wz'*(Z'*Z)*wz
00152   //    => grad = (Z'*Y*wy), Hess = Z'*Z
00153   //
00154   // OBJ_MIN_WZ
00155   //    min 1/2 * wz'*wz => grad = 0, Hess = I
00156   //
00157   // OBJ_RSQP
00158   //    min qp_grad_k'*wz + 1/2 * wz'*rHL_k*wz
00159   //    => grad = qp_grad, Hess = rHL_k
00160   //
00161   const MatrixOp
00162     &Z_k = s->Z().get_k(0);
00163   if( current_k_ != s->k() ) {
00164     if( qp_objective() != OBJ_RSQP )
00165       grad_store_.resize(n-r);
00166     if( qp_objective() == OBJ_MIN_FULL_STEP )
00167       Hess_store_.resize(n-r+1,n-r+1);
00168   }
00169   DVectorSlice grad;
00170   switch(qp_objective())
00171   {
00172       case OBJ_MIN_FULL_STEP: // grad = (Z'*Ywy), Hess = Z'*Z
00173     {
00174       grad.bind( grad_store_() );
00175       if( current_k_ != s->k() ) {
00176         // grad = (Z'*Ywy)
00177         LinAlgOpPack::V_MtV( &grad, Z_k, BLAS_Cpp::trans, Ywy );
00178         // Hess = Z'*Z
00179         DMatrixSliceSym S(Hess_store_(2,n-r+1,1,n-r),BLAS_Cpp::lower); // Must be strictly lower triangular here!
00180         Z_k.syrk( BLAS_Cpp::trans, 1.0, 0.0, &S ); // S = 1.0*Z'*Z + 0.0*S
00181         MatrixSymPosDefCholFactor
00182           *H_ptr = NULL;
00183         if( Hess_ptr_.get() == NULL || dynamic_cast<const MatrixSymPosDefCholFactor*>(Hess_ptr_.get()) == NULL )
00184           Hess_ptr_ = new MatrixSymPosDefCholFactor;
00185         H_ptr = const_cast<MatrixSymPosDefCholFactor*>(dynamic_cast<const MatrixSymPosDefCholFactor*>(Hess_ptr_.get()));
00186         TEUCHOS_TEST_FOR_EXCEPT( !( H_ptr ) ); // Should not be null!
00187         H_ptr->init_setup(
00188           &Hess_store_()  // The original matrix is stored in the lower triangular part (below diagonal)!
00189           ,NULL           // Nothing to deallocate
00190           ,n-r
00191           ,true           // maintain the original factor
00192           ,false          // don't maintain the factor
00193           ,true           // allow the factor to be computed if needed
00194           ,true           // Set the view
00195           ,1.0            // Scale the matrix by one
00196           );
00197       }
00198       break;
00199     }
00200       case OBJ_MIN_NULL_SPACE_STEP: // grad = 0, Hess = I
00201     {
00202       grad.bind( grad_store_() );
00203       MatrixSymIdent
00204         *H_ptr = NULL;
00205       if( Hess_ptr_.get() == NULL || dynamic_cast<const MatrixSymIdent*>(Hess_ptr_.get()) == NULL )
00206         Hess_ptr_ = new MatrixSymIdent;
00207       if( current_k_ != s->k() ) {
00208         H_ptr = const_cast<MatrixSymIdent*>(dynamic_cast<const MatrixSymIdent*>(Hess_ptr_.get()));
00209         TEUCHOS_TEST_FOR_EXCEPT( !( H_ptr ) ); // Should not be null!
00210         H_ptr->init_setup(n-r,1.0);
00211         grad = 0.0;
00212       }
00213       break;
00214     }
00215       case OBJ_RSQP: // grad = qp_grad, Hess = rHL_k
00216     {
00217       grad.bind( s->qp_grad().get_k(0)() );
00218       Hess_ptr_ = Hess_ptr_t( &s->rHL().get_k(0), false ); // don't delete memory!
00219       break;
00220     }
00221       defaut:
00222       TEUCHOS_TEST_FOR_EXCEPT(true); // Not a valid option
00223   }
00224 
00225   //
00226   // Solve the null space subproblem
00227   //
00228 
00229   Workspace<value_type>  wz_ws(wss,n-r),Zwz_ws(wss,n);
00230   DVectorSlice                 wz(&wz_ws[0],wz_ws.size());
00231   DVectorSlice                 Zwz(&Zwz_ws[0],Zwz_ws.size());
00232   value_type                  qp_eta      = 0;
00233 
00234   bool throw_qp_failure = false;
00235 
00236   if( bl.nz() == 0 && bu.nz() == 0 && m-r == 0 ) {
00237     //
00238     // Just solve the unconstrainted QP
00239     //
00240     // wz = - inv(Hess)*grad
00241 #ifdef _WINDOWS
00242     const MatrixFactorized &Hess = dynamic_cast<const MatrixFactorized&>(*Hess_ptr_);
00243 #else
00244     const MatrixFactorized &Hess = dyn_cast<const MatrixFactorized>(*Hess_ptr_);
00245 #endif
00246     AbstractLinAlgPack::V_InvMtV( &wz, Hess, BLAS_Cpp::no_trans, grad );
00247     DenseLinAlgPack::Vt_S(&wz,-1.0);
00248     // Zwz = Z*wz
00249     LinAlgOpPack::V_MtV( &Zwz, Z_k, BLAS_Cpp::no_trans, wz );
00250   }
00251   else {
00252 
00253     //
00254     // Set the arguments to the QP subproblem
00255     //
00256 
00257     DVectorSlice      qp_g    = grad;
00258     const MatrixOp&   qp_G    = *Hess_ptr_;
00259     const value_type    qp_etaL   = 0.0;
00260     SpVectorSlice     qp_dL(NULL,0,0,n-r);  // If nz() == 0 then no simple bounds
00261     SpVectorSlice     qp_dU(NULL,0,0,n-r);
00262     const MatrixOp    *qp_E   = NULL;
00263     BLAS_Cpp::Transp    qp_trans_E  = BLAS_Cpp::no_trans;
00264     DVectorSlice             qp_b;
00265     SpVectorSlice     qp_eL(NULL,0,0,n);
00266     SpVectorSlice     qp_eU(NULL,0,0,n);
00267     const MatrixOp    *qp_F   = NULL;
00268     BLAS_Cpp::Transp    qp_trans_F  = BLAS_Cpp::no_trans;
00269     DVectorSlice        qp_f;
00270     DVectorSlice        qp_d    = wz;
00271     SpVector        *qp_nu    = NULL;
00272     SpVector        *qp_mu    = NULL;
00273     DVectorSlice        qp_Ed;
00274     DVectorSlice        qp_lambda;
00275 
00276     SpVector _nu_wz, _nu_Dwz,   // Possible storage for multiplers for separate inequality
00277       _nu;        // constriants for wz.
00278     DVector _Dwz;       // Possible storage for D*wz computed by QP solver?
00279 
00280     //
00281     // Determine if we can use simple bounds on wz.
00282     // 
00283     // If we have a variable reduction null space matrix
00284     // (with any choice for Y) then:
00285     // 
00286     // w = Z*wz + (1-eta) * Y*wy
00287     // 
00288     // [ w(dep)   ]  = [ D ] * wz  + (1-eta) * [ Ywy(dep)   ]
00289     // [ w(indep) ]    [ I ]                   [ Ywy(indep) ]
00290     // 
00291     // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ywy(indep) = 0 and
00292     // in this case the bounds on d(indep) become simple bounds on pz even
00293     // with the relaxation.
00294     // 
00295     // Otherwise, we can not use simple variable bounds and implement the
00296     // relaxation properly.
00297     // 
00298 
00299     const ZVarReductMatrix
00300       *Zvr = dynamic_cast<const ZVarReductMatrix*>( &Z_k );
00301     Range1D
00302       indep   = Zvr ? Zvr->indep() : Range1D(),
00303       dep   = Zvr ? Zvr->dep()   : Range1D();
00304 
00305     const bool
00306       use_simple_wz_bounds = ( Zvr!=NULL && DenseLinAlgPack::norm_inf(Ywy(indep))==0.0 );
00307 
00308     if( use_simple_wz_bounds ) {
00309 
00310       // Set simple bound constraints on pz
00311       qp_dL.bind( bl(indep) );
00312       qp_dU.bind( bu(indep) );
00313       qp_nu = &( _nu_wz = s->nu().get_k(0)(indep) );  // warm start?
00314     
00315       // Set general inequality constraints for D*pz
00316       qp_E = &Zvr->D();
00317       qp_b.bind( Ywy(dep) );
00318       qp_eL.bind( bl(dep) );
00319       qp_eU.bind( bu(dep) );
00320       qp_mu = &( _nu_Dwz = s->nu().get_k(0)(dep) ); // warm start?
00321       _Dwz.resize(r);
00322       qp_Ed.bind(_Dwz()); // Part of Zwz will be updated directly!
00323 
00324     }
00325     else {
00326 
00327       // Set general inequality constraints for Z*pz
00328       qp_E = &Z_k;
00329       qp_b.bind( Ywy() );
00330       qp_eL.bind( bl() );
00331       qp_eU.bind( bu() );
00332       qp_mu = &(_nu = s->nu().get_k(0));  // warm start??
00333       qp_Ed.bind(Zwz);  // Zwz
00334     }
00335 
00336     // Set the general equality constriants (if they exist)
00337     DVector q(m-r);
00338     Range1D undecomp = s->con_undecomp();
00339     if( m > r ) {
00340       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
00341     }
00342 
00343     // Setup the rest of the arguments
00344 
00345     QPSolverRelaxed::EOutputLevel
00346       qp_olevel;
00347     switch( olevel ) {
00348         case PRINT_NOTHING:
00349         qp_olevel = QPSolverRelaxed::PRINT_NONE;
00350         break;
00351         case PRINT_BASIC_ALGORITHM_INFO:
00352         qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00353         break;
00354         case PRINT_ALGORITHM_STEPS:
00355         qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00356         break;
00357         case PRINT_ACTIVE_SET:
00358         qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
00359         break;
00360         case PRINT_VECTORS:
00361         qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
00362         break;
00363         case PRINT_ITERATION_QUANTITIES:
00364         qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
00365         break;
00366         default:
00367         TEUCHOS_TEST_FOR_EXCEPT(true);
00368     }
00369 
00370     //
00371     // Solve the QP
00372     // 
00373     const QPSolverStats::ESolutionType
00374       solution_type =
00375       qp_solver().solve_qp(
00376         int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00377         , qp_olevel
00378         , algo->algo_cntr().check_results()
00379         ? QPSolverRelaxed::RUN_TESTS :  QPSolverRelaxed::NO_TESTS
00380         , qp_g, qp_G, qp_etaL, qp_dL, qp_dU
00381         , qp_E, qp_trans_E, qp_E ? &qp_b : NULL
00382         , qp_E ? &qp_eL : NULL, qp_E ? &qp_eU : NULL
00383         , qp_F, qp_trans_F, qp_F ? &qp_f : NULL
00384         , NULL
00385         , &qp_eta, &qp_d
00386         , qp_nu
00387         , qp_mu, qp_E ? &qp_Ed : NULL
00388         , qp_F ? &qp_lambda : NULL, NULL
00389         );
00390 
00391     //
00392     // Check the optimality conditions for the QP
00393     //
00394     std::ostringstream omsg;
00395     if(   qp_testing() == QP_TEST
00396         || ( qp_testing() == QP_TEST_DEFAULT && algo->algo_cntr().check_results() )  )
00397     {
00398       if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
00399         out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
00400       }
00401       if(!qp_tester().check_optimality_conditions(
00402         solution_type
00403         , int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00404         , int(olevel) >= int(PRINT_VECTORS) ? true : false
00405         , int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
00406         , qp_g, qp_G, qp_etaL, qp_dL, qp_dU
00407         , qp_E, qp_trans_E, qp_E ? &qp_b : NULL
00408         , qp_E ? &qp_eL : NULL, qp_E ? &qp_eU : NULL
00409         , qp_F, qp_trans_F, qp_F ? &qp_f : NULL
00410         , NULL
00411         , &qp_eta, &qp_d
00412         , qp_nu
00413         , qp_mu, qp_E ? &qp_Ed : NULL
00414         , qp_F ? &qp_lambda : NULL, NULL
00415         ))
00416       {
00417         omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
00418         if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00419           out << omsg.str();
00420         }
00421         throw_qp_failure = true;
00422       }
00423     }
00424 
00425     if( solution_type !=  QPSolverStats::OPTIMAL_SOLUTION ) {
00426       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00427         out << "\nCould not solve the QP!\n";
00428       }
00429       return false;
00430     }
00431 
00432     //
00433     // Set the solution
00434     //
00435     if( use_simple_wz_bounds ) {
00436       // Set Zwz
00437       Zwz(dep)   = _Dwz;
00438       Zwz(indep) = wz;
00439     }
00440     else {
00441       // Everything should already be set!
00442     }
00443 
00444     // Cut back Ywy = (1-eta) * Ywy
00445     const value_type eps = std::numeric_limits<value_type>::epsilon();
00446     if( fabs(qp_eta - 0.0) > eps ) {
00447       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00448         out
00449           << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<").  Cutting back Ywy_k = (1.0 - eta)*Ywy  ...\n";
00450       }
00451       DenseLinAlgPack::Vt_S( &Ywy , 1.0 - qp_eta );
00452     }
00453   }
00454 
00455   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00456     out << "\n||wz||inf    = " << DenseLinAlgPack::norm_inf(wz);
00457     out << "\n||Zwz||2     = " << DenseLinAlgPack::norm_2(Zwz);
00458     if(qp_eta > 0.0) out << "\n||Ypy||2 = " << DenseLinAlgPack::norm_2(Ywy);
00459     out << std::endl;
00460   }
00461   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00462     out << "\nwz = \n" << wz;
00463     out << "\nZwz = \n" << Zwz;
00464     if(qp_eta > 0.0) out << "\nYwy = \n" << Ywy;
00465   }
00466   if( qp_eta == 1.0 ) {
00467     std::ostringstream omsg;
00468     omsg
00469       << "FeasibilityStepReducedStd_Strategy::compute_feasibility_step(...) : "
00470       << "Error, a QP relaxation parameter\n"
00471       << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
00472       << "that the NLP's constraints are infeasible\n"
00473       << "Throwing an InfeasibleConstraints exception!\n";
00474     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00475       out << omsg.str();
00476     }
00477     throw_qp_failure = true;
00478 //    throw InfeasibleConstraints(omsg.str());
00479   }
00480 
00481   //
00482   // Set the full step
00483   //
00484   // w = Ywy + Zwz
00485   //
00486   DenseLinAlgPack::V_VpV( w, Ywy, Zwz );
00487 
00488   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00489     out << "\n||w||inf    = " << DenseLinAlgPack::norm_inf(*w);
00490     out << std::endl;
00491   }
00492   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00493     out << "\nw = \n" << *w;
00494   }
00495 
00496   current_k_ = s->k();
00497 
00498   if( throw_qp_failure )
00499     return false;
00500 
00501 */
00502   TEUCHOS_TEST_FOR_EXCEPT(true);
00503 
00504   return true;
00505 }
00506 
00507 void FeasibilityStepReducedStd_Strategy::print_step( std::ostream& out, const std::string& L ) const
00508 {
00509   out << L << "*** Computes feasibility steps by solving a constrained QP using the range and null\n"
00510     << L << "*** space decomposition\n"
00511     << L << "begin quais-range space step: \"" << typeName(quasi_range_space_step()) << "\"\n";
00512 
00513    quasi_range_space_step().print_step( out, L + "  " );
00514 
00515   out << L << "end quasi-range space step\n"
00516     << L << "if quasi-range space step failed then\n"
00517     << L << "  this feasibility step fails also!\n"
00518     << L << "end\n"
00519     << L << "Ywy = v\n"
00520     << L << "*** Set the bounds for bl <= Z*wz <= bu\n"
00521     << L << "bl = d_bounds_k.l - (xo - x_k) - Ywy\n"
00522     << L << "bu = d_bounds_k.u - (xo - x_k) - Ywy\n"
00523     << L << "Set bl(i) = bu(i) for those nu_k(i) != 0.0\n"
00524     << L << "if (qp_objective == OBJ_MIN_FULL_STEP) and (current_k != k) then\n"
00525     << L << "  grad = Z_k'*Ywy\n"
00526     << L << "  Hess = Z_k'*Z_k\n"
00527     << L << "elseif (qp_objective == OBJ_MIN_NULL_SPACE_STEP) and (current_k != k) then\n"
00528     << L << "  grad = 0\n"
00529     << L << "  Hess = I\n"
00530     << L << "elseif (qp_objective == OBJ_RSQP) and (current_k != k) then\n"
00531     << L << "  grad = qp_grad_k\n"
00532     << L << "  Hess = rHL_k\n"
00533     << L << "end\n"
00534     << L << "if check_results == true then\n"
00535     << L << "  assert that bl and bu are valid and sorted\n"
00536     << L << "end\n"
00537     << L << "etaL = 0.0\n"
00538     << L << "*** Determine if we can use simple bounds on pz or not\n"
00539     << L << "if Z_k is a variable reduction null space matrix and norm(Ypy_k(indep),0) == 0 then\n"
00540     << L << "  use_simple_wz_bounds = true\n"
00541     << L << "else\n"
00542     << L << "  use_simple_wz_bounds = false\n"
00543     << L << "end\n"
00544     << L << "*** Setup QP arguments\n"
00545     << L << "qp_g = qp_grad_k\n"
00546     << L << "qp_G = rHL_k\n"
00547     << L << "if use_simple_wz_bounds == true then\n"
00548     << L << "  qp_dL = bl(indep),  qp_dU = bu(indep)\n"
00549     << L << "  qp_E  = Z_k.D,      qp_b  = Ywy(dep)\n"
00550     << L << "  qp_eL = bl(dep),    qp_eU = bu(dep)\n"
00551     << L << "else\n"
00552     << L << "  qp_dL = -inf,       qp_dU = +inf\n"
00553     << L << "  qp_E  = Z_k,        qp_b  = Ywy\n"
00554     << L << "  qp_eL = bl,         qp_eU = bu\n"
00555     << L << "end\n"
00556     << L << "if m > r then\n"
00557     << L << "  qp_F  = V_k,        qp_f  = c_k(undecomp) + Gc_k(undecomp)'*Ywy\n"
00558     << L << "else\n"
00559     << L << "  qp_F  = empty,      qp_f  = empty\n"
00560     << L << "end\n"
00561     << L << "Use a warm start given the active-set in nu_k\n"
00562     << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
00563     << L << ",qp_nu, qp_mu and qp_lambda (" << typeName(qp_solver()) << "):\n"
00564     << L << "  min    qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
00565     << L << "  qp_d <: R^(n-r)\n"
00566     << L << "  s.t.\n"
00567     << L << "         etaL  <=  qp_eta\n"
00568     << L << "         qp_dL <= qp_d <= qp_dU                          [qp_nu]\n"
00569     << L << "         qp_eL <= qp_E * qp_d + (1-eta)*qp_b  <= qp_eU   [qp_mu]\n"
00570     << L << "         qp_F * d_qp + (1-eta) * qp_f = 0                [qp_lambda]\n"
00571     << L << "if (qp_teing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
00572     << L << "and check_results==true) then\n"
00573     << L << "  Check the optimality conditions of the above QP\n"
00574     << L << "  if the optimality conditions do not check out then\n"
00575     << L << "    set throw_qp_failure = true\n"
00576     << L << "  end\n"
00577     << L << "end\n"
00578     << L << "*** Set the solution to the QP subproblem\n"
00579     << L << "wz  = qp_d\n"
00580     << L << "eta = qp_eta\n"
00581     << L << "if use_simple_wz_bounds == true then\n"
00582     << L << "  Zwz(dep)   = qp_Ed,  Zwz(indep) = wz\n"
00583     << L << "else\n"
00584     << L << "  Zwz = qp_Ed\n"
00585     << L << "end\n"
00586     << L << "if eta > 0 then set Ywy = (1-eta) * Ywy\n"
00587     << L << "if QP solution is suboptimal then\n"
00588     << L << "  throw_qp_failure = true\n"
00589     << L << "elseif QP solution is primal feasible (not optimal) then\n"
00590     << L << "  throw_qp_failure = primal_feasible_point_error\n"
00591     << L << "elseif QP solution is dual feasible (not optimal) then\n"
00592     << L << "  find max u s.t.\n"
00593     << L << "    d_bounds_k.l <= (xo - x) + u*(Ywy+Zwz) <= d_bounds_k.u\n"
00594     << L << "  alpha_k = u\n"
00595     << L << "  throw_qp_failure = true\n"
00596     << L << "end\n"
00597     << L << "if eta == 1.0 then\n"
00598     << L << "  The constraints are infeasible!\n"
00599     << L << "  throw_qp_failure = true\n"
00600     << L << "end\n"
00601     << L << "current_k = k\n"
00602     << L << "w = Zwz + Ywy\n"
00603     << L << "if (throw_qp_failure == true) then\n"
00604     << L << "  The feasibility step computation has failed!\n"
00605     << L << "end\n"
00606     ;
00607 }
00608 
00609 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines