ConstrainedOptPack_QPSolverRelaxedTester.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 // 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 <math.h>
00030 
00031 #include <limits>
00032 
00033 #include "ConstrainedOptPack_QPSolverRelaxedTester.hpp"
00034 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00035 #include "AbstractLinAlgPack_VectorSpace.hpp"
00036 #include "AbstractLinAlgPack_VectorMutable.hpp"
00037 #include "AbstractLinAlgPack_VectorOut.hpp"
00038 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00039 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00040 
00041 namespace {
00042 
00043 //
00044 const char* solution_type_str( ConstrainedOptPack::QPSolverStats::ESolutionType solution_type )
00045 {
00046   typedef ConstrainedOptPack::QPSolverStats qpst;
00047   switch( solution_type ) {
00048   
00049   case qpst::OPTIMAL_SOLUTION:
00050     return "OPTIMAL_SOLUTION";
00051   case qpst::PRIMAL_FEASIBLE_POINT:
00052     return "PRIMAL_FEASIBLE_POINT";
00053   case qpst::DUAL_FEASIBLE_POINT:
00054     return "DUAL_FEASIBLE_POINT";
00055   case qpst::SUBOPTIMAL_POINT:
00056     return "SUBOPTIMAL_POINT";
00057   default:
00058     TEST_FOR_EXCEPT(true);
00059   }
00060   return "";  // will never be executed.
00061 }
00062 
00063 /* ToDo: Update this code!
00064 
00065 // Compute the scaled complementarity conditions.
00066 // 
00067 // If uplo == upper then:
00068 // 
00069 //                / gamma(i) * constr_resid(i) / ( 1 + |constr(i)| + opt_scale ), for gamma(i) > 0
00070 // comp_err(i) = |
00071 //                \ 0 otherwise
00072 // 
00073 // If uplo == lower then:
00074 // 
00075 //                / gamma(i) * constr_resid(i) /  ( 1 + |constr(i)| + opt_scale ), for gamma(i) < 0
00076 // comp_err(i) = |
00077 //                \ 0 otherwise
00078 // 
00079 // 
00080 void set_complementarity(
00081   const AbstractLinAlgPack::SpVector  &gamma
00082   ,const DenseLinAlgPack::DVectorSlice    &constr_resid
00083   ,const DenseLinAlgPack::DVectorSlice     &constr
00084   ,const DenseLinAlgPack::value_type      opt_scale
00085   ,BLAS_Cpp::Uplo         uplo
00086   ,DenseLinAlgPack::DVector       *comp_err
00087   )
00088 {
00089   TEST_FOR_EXCEPT( !(  gamma.size() == constr_resid.size() && gamma.size() == constr.size()  ) );
00090   comp_err->resize( gamma.size() );
00091   *comp_err = 0.0;
00092   const AbstractLinAlgPack::SpVector::difference_type o = gamma.offset();
00093   if( gamma.nz() ) {
00094     for( AbstractLinAlgPack::SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
00095       const DenseLinAlgPack::size_type i = itr->indice() + o;
00096       if( itr->value() > 0 && uplo == BLAS_Cpp::upper )
00097         (*comp_err)(i) = itr->value() * constr_resid(i) / ( 1.0 + ::fabs(constr(i)) + opt_scale );
00098       else if( itr->value() < 0 && uplo == BLAS_Cpp::lower )
00099         (*comp_err)(i) = itr->value() * constr_resid(i) / ( 1.0 + ::fabs(constr(i)) + opt_scale );
00100     }
00101   }
00102 }
00103 
00104 */
00105 
00106 // Handle the error reporting
00107 void handle_error(
00108   std::ostream                            *out
00109   ,const AbstractLinAlgPack::value_type   err
00110   ,const char                             err_name[]
00111   ,const AbstractLinAlgPack::value_type   error_tol
00112   ,const char                             error_tol_name[]
00113   ,const AbstractLinAlgPack::value_type   warning_tol
00114   ,const char                             warning_tol_name[]
00115   ,bool                                   *test_failed
00116   )
00117 {
00118   if( err >= error_tol ) {
00119     if(out)
00120       *out << "\n" << err_name << " = " << err << " >= " << error_tol_name << " = " << error_tol << std::endl;
00121     *test_failed = true;
00122   }
00123   else if( err >= warning_tol ) {
00124     if(out)
00125       *out << "\n" << err_name << " = " << err << " >= " << warning_tol_name << " = " << warning_tol << std::endl;
00126   }
00127 }
00128 
00129 } // end namespace
00130 
00131 namespace ConstrainedOptPack {
00132 
00133 // public
00134 
00135 QPSolverRelaxedTester::QPSolverRelaxedTester(
00136   value_type   opt_warning_tol
00137   ,value_type  opt_error_tol
00138   ,value_type  feas_warning_tol
00139   ,value_type  feas_error_tol
00140   ,value_type  comp_warning_tol
00141   ,value_type  comp_error_tol
00142   )
00143   :opt_warning_tol_(opt_warning_tol)
00144   ,opt_error_tol_(opt_error_tol)
00145   ,feas_warning_tol_(feas_warning_tol)
00146   ,feas_error_tol_(feas_error_tol)
00147   ,comp_warning_tol_(comp_warning_tol)
00148   ,comp_error_tol_(comp_error_tol)
00149 {}
00150 
00151 bool QPSolverRelaxedTester::check_optimality_conditions(
00152   QPSolverStats::ESolutionType solution_type
00153   ,const value_type infinite_bound
00154   ,std::ostream* out, bool print_all_warnings, bool print_vectors
00155   ,const Vector& g, const MatrixSymOp& G
00156   ,value_type etaL
00157   ,const Vector& dL, const Vector& dU
00158   ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
00159   ,const Vector& eL, const Vector& eU
00160   ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
00161   ,const value_type* obj_d
00162   ,const value_type* eta, const Vector* d
00163   ,const Vector* nu
00164   ,const Vector* mu, const Vector* Ed
00165   ,const Vector* lambda, const Vector* Fd
00166   )
00167 {
00168   return check_optimality_conditions(
00169     solution_type,infinite_bound,out,print_all_warnings,print_vectors
00170     ,g,G,etaL,&dL,&dU,&E,trans_E,&b,&eL,&eU,&F,trans_F,&f
00171     ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00172 }
00173 
00174 bool QPSolverRelaxedTester::check_optimality_conditions(
00175   QPSolverStats::ESolutionType solution_type
00176   ,const value_type infinite_bound
00177   ,std::ostream* out, bool print_all_warnings, bool print_vectors
00178   ,const Vector& g, const MatrixSymOp& G
00179   ,value_type etaL
00180   ,const Vector& dL, const Vector& dU
00181   ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
00182   ,const Vector& eL, const Vector& eU
00183   ,const value_type* obj_d
00184   ,const value_type* eta, const Vector* d
00185   ,const Vector* nu
00186   ,const Vector* mu, const Vector* Ed
00187   )
00188 {
00189   return check_optimality_conditions(
00190     solution_type,infinite_bound,out,print_all_warnings,print_vectors
00191     ,g,G,etaL,&dL,&dU,&E,trans_E,&b,&eL,&eU,NULL,BLAS_Cpp::no_trans,NULL
00192     ,obj_d,eta,d,nu,mu,Ed,NULL,NULL);
00193 }
00194 
00195 bool QPSolverRelaxedTester::check_optimality_conditions(
00196   QPSolverStats::ESolutionType solution_type
00197   ,const value_type infinite_bound
00198   ,std::ostream* out, bool print_all_warnings, bool print_vectors
00199   ,const Vector& g, const MatrixSymOp& G
00200   ,value_type etaL
00201   ,const Vector& dL, const Vector& dU
00202   ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
00203   ,const value_type* obj_d
00204   ,const value_type* eta, const Vector* d
00205   ,const Vector* nu
00206   ,const Vector* lambda, const Vector* Fd
00207   )
00208 {
00209   return check_optimality_conditions(
00210     solution_type,infinite_bound,out,print_all_warnings,print_vectors
00211     ,g,G,etaL,&dL,&dU,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL,&F,trans_F,&f
00212     ,obj_d,eta,d,nu,NULL,NULL,lambda,Fd );
00213 }
00214 
00215 bool QPSolverRelaxedTester::check_optimality_conditions(
00216   QPSolverStats::ESolutionType solution_type
00217   ,const value_type infinite_bound
00218   ,std::ostream* out, bool print_all_warnings, bool print_vectors
00219   ,const Vector& g, const MatrixSymOp& G
00220   ,const Vector& dL, const Vector& dU
00221   ,const value_type* obj_d
00222   ,const Vector* d
00223   ,const Vector* nu
00224   )
00225 {
00226   return check_optimality_conditions(
00227     solution_type,infinite_bound,out,print_all_warnings,print_vectors
00228     ,g,G,0.0,&dL,&dU,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL,NULL,BLAS_Cpp::no_trans,NULL
00229     ,obj_d,NULL,d,nu,NULL,NULL,NULL,NULL);
00230 }
00231 
00232 bool QPSolverRelaxedTester::check_optimality_conditions(
00233   QPSolverStats::ESolutionType solution_type
00234   ,const value_type infinite_bound
00235   ,std::ostream* out, bool print_all_warnings, bool print_vectors
00236   ,const Vector& g, const MatrixSymOp& G
00237   ,value_type etaL
00238   ,const Vector* dL, const Vector* dU
00239   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00240   ,const Vector* eL, const Vector* eU
00241   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00242   ,const value_type* obj_d
00243   ,const value_type* eta, const Vector* d
00244   ,const Vector* nu
00245   ,const Vector* mu, const Vector* Ed
00246   ,const Vector* lambda, const Vector* Fd
00247   )
00248 {
00249   QPSolverRelaxed::validate_input(
00250     infinite_bound,g,G,etaL,dL,dU
00251     ,E,trans_E,b,eL,eU,F,trans_F,f
00252     ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00253 
00254   return imp_check_optimality_conditions(
00255     solution_type,infinite_bound
00256     ,out,print_all_warnings,print_vectors,g,G,etaL,dL,dU
00257     ,E,trans_E,b,eL,eU,F,trans_F,f
00258     ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00259 }
00260 
00261 // protected
00262 
00263 bool QPSolverRelaxedTester::imp_check_optimality_conditions(
00264   QPSolverStats::ESolutionType solution_type
00265   ,const value_type infinite_bound
00266   ,std::ostream* out, bool print_all_warnings, bool print_vectors
00267   ,const Vector& g, const MatrixSymOp& G
00268   ,value_type etaL
00269   ,const Vector* dL, const Vector* dU
00270   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00271   ,const Vector* eL, const Vector* eU
00272   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00273   ,const value_type* obj_d
00274   ,const value_type* eta, const Vector* d
00275   ,const Vector* nu
00276   ,const Vector* mu, const Vector* Ed
00277   ,const Vector* lambda, const Vector* Fd
00278   )
00279 {
00280   using std::endl;
00281   using BLAS_Cpp::trans_not;
00282   using BLAS_Cpp::no_trans;
00283   using BLAS_Cpp::trans;
00284   using BLAS_Cpp::upper;
00285   using BLAS_Cpp::lower;
00286   using LinAlgOpPack::sum;
00287   using LinAlgOpPack::dot;
00288   using LinAlgOpPack::Vt_S;
00289   using LinAlgOpPack::V_VmV;
00290   using LinAlgOpPack::Vp_StV;
00291   using LinAlgOpPack::V_MtV;
00292   using LinAlgOpPack::Vp_MtV;
00293   using LinAlgOpPack::V_StMtV;
00294   using LinAlgOpPack::Vp_V;
00295   using AbstractLinAlgPack::max_element;
00296   typedef QPSolverStats qps_t;
00297 
00298   bool test_failed = false;
00299 
00300   const size_type
00301     nd = d->dim();
00302 
00303   const value_type
00304     really_big_error_tol = std::numeric_limits<value_type>::max();
00305 
00306   value_type opt_scale = 0.0;
00307   VectorSpace::vec_mut_ptr_t
00308     t_d = d->space().create_member(),
00309     u_d = d->space().create_member();
00310 
00311   value_type
00312     err = 0,
00313     d_norm_inf, // holds ||d||inf
00314     e_norm_inf; // holds ||e||inf
00315 
00316   if(out)
00317     *out
00318       << "\n*** Begin checking QP optimality conditions ***\n"
00319       << "\nThe solution type is " << solution_type_str(solution_type) << endl;
00320 
00321   bool force_opt_error_check
00322     = solution_type==qps_t::OPTIMAL_SOLUTION || solution_type==qps_t::DUAL_FEASIBLE_POINT;
00323   const bool force_inequality_error_check
00324     = solution_type==qps_t::OPTIMAL_SOLUTION || solution_type==qps_t::PRIMAL_FEASIBLE_POINT;
00325   const bool force_equality_error_check
00326     = solution_type!=qps_t::SUBOPTIMAL_POINT;
00327   const bool force_complementarity_error_check
00328     = solution_type!=qps_t::SUBOPTIMAL_POINT;
00329 
00330   const char sep_line[] = "\n--------------------------------------------------------------------------------\n";
00331 
00333   // Checking d(L)/d(d) = 0
00334   if(out)
00335     *out
00336       << sep_line
00337       << "\nChecking d(L)/d(d) = g + G*d + nu + op(E)'*mu - op(F)'*lambda == 0 ...\n";
00338 
00339   if(out && !force_opt_error_check)
00340     *out
00341       << "The optimality error tolerance will not be enforced ...\n";
00342 
00343   opt_scale = 0.0;
00344 
00345   *u_d = g;
00346   opt_scale += g.norm_inf();
00347 
00348   if(out) {
00349     *out << "||g||inf = " << g.norm_inf() << endl;
00350   }
00351   
00352   V_MtV( t_d.get(), G, no_trans, *d );
00353   Vp_V( u_d.get(), *t_d );
00354   opt_scale += t_d->norm_inf();
00355 
00356   if(out) {
00357     *out << "||G*d||inf = " << t_d->norm_inf() << endl;
00358     if(print_vectors)
00359       *out << "g + G*d =\n" << *u_d;
00360   }
00361 
00362   if( nu ) {
00363     Vp_V( u_d.get(), *nu );
00364     opt_scale += nu->norm_inf();
00365     if(out)
00366       *out << "||nu||inf = " << nu->norm_inf() << endl;
00367   }
00368 
00369   if(E) {
00370     V_MtV( t_d.get(), *E, trans_not(trans_E), *mu );
00371     Vp_V( u_d.get(), *t_d );
00372     opt_scale += t_d->norm_inf();
00373     if(out) {
00374       *out << "||op(E)'*mu||inf = " << t_d->norm_inf() << endl;
00375       if(print_vectors)
00376         *out << "op(E)'*mu =\n" << *t_d;
00377     }
00378   }     
00379 
00380   if(F) {
00381     V_MtV( t_d.get(), *F, trans_not(trans_F), *lambda );
00382     Vp_V( u_d.get(), *t_d );
00383     opt_scale += t_d->norm_inf();
00384     if(out) {
00385       *out << "||op(F)'*lambda||inf = " << t_d->norm_inf() << endl;
00386       if(print_vectors)
00387         *out << "op(F)'*lambda =\n" << *t_d;
00388     }
00389   }
00390 
00391   if( *eta > etaL ) { // opt_scale + |(eta - etaL) * (b'*mu + f'*lambda)|
00392     const value_type
00393       term = ::fabs( (*eta - etaL) * (E ? dot(*b,*mu) : 0.0) + (F ? dot(*f,*lambda) : 0.0 ) );
00394     if(out) {
00395       *out << "|(eta - etaL) * (b'*mu + f'*lambda)| = " << term << endl;
00396     }
00397     opt_scale += term;
00398   }
00399 
00400   if(out && print_vectors)
00401     *out
00402       << "g + G*d + nu + op(E)'*mu - op(F)'*lambda =\n" << *u_d;
00403 
00404   Vt_S( u_d.get(), 1.0/(1.0+opt_scale) );
00405 
00406   err = sum( *u_d );
00407 
00408   if(out)
00409     *out
00410       << "\nopt_scale = " << opt_scale << endl
00411       << "opt_err = sum( | g + G*d + nu + op(E)'*mu - op(F)'*lambda | / (1+opt_scale) ) / nd\n"
00412       << "        = " << err << " / " << nd << " = " << (err/nd) << endl;
00413 
00414   err *= nd;
00415 
00416   if( force_opt_error_check ) {
00417     if( err >= opt_error_tol() ) {
00418       if(out)
00419         *out << "\nopt_err = " << err << " >= opt_error_tol = " << opt_error_tol() << endl;
00420       test_failed = true;
00421     }
00422     else if( err >= opt_warning_tol() ) {
00423       if(out)
00424         *out << "\nopt_err = " << err << " >= opt_error_tol = " << opt_error_tol() << endl;
00425     }
00426   }
00427 
00428   if(out) {
00429     *out
00430       << sep_line
00431       << "\nTesting feasibility of the constraints and the complementarity conditions ...\n";
00432     if(!force_inequality_error_check)
00433       *out
00434         << "The inequality feasibility error tolerance will not be enforced ...\n";
00435     if(!force_equality_error_check)
00436       *out
00437         << "The equality feasibility error tolerance will not be enforced ...\n";
00438     if(!force_complementarity_error_check)
00439       *out
00440         << "The complementarity conditions error tolerance will not be enforced ...\n";
00441   }
00442 
00444   // etaL - eta
00445   if(out)
00446     *out
00447       << sep_line
00448       << "\nChecking etaL - eta <= 0 ...\n";
00449   if(out)
00450     *out
00451       << "etaL - eta = " << (etaL - (*eta)) << endl;
00452   if( etaL - (*eta) > feas_warning_tol() ) {
00453     if(out)
00454       *out
00455         << "Warning, etaL - eta = " << etaL << " - " << (*eta)
00456         << " = " << (etaL - (*eta)) << " >  feas_warning_tol = "
00457         <<  feas_warning_tol() << endl;
00458   }
00459   if( force_inequality_error_check && etaL - (*eta) > feas_error_tol() ) {
00460     if(out)
00461       *out
00462         << "Error, etaL - eta = " << etaL << " - " << (*eta)
00463         << " = " << (etaL - (*eta)) << " >  feas_error_tol = "
00464         <<  feas_error_tol() << endl;
00465     test_failed = true;
00466   } 
00467 
00468   d_norm_inf = d->norm_inf();
00469 
00470   if(dL) {
00471 
00473     // dL - d <= 0
00474     if(out)
00475       *out
00476         << sep_line
00477         << "\nChecking dL - d <= 0 ...\n";
00478     V_VmV( u_d.get(), *dL, *d );
00479     if(out && print_vectors)
00480       *out
00481         << "dL - d =\n" << *u_d;
00482     Vt_S( u_d.get(), 1.0/(1.0+d_norm_inf) );
00483 
00484     err = max_element(*u_d);
00485     if(out)
00486       *out
00487         << "\nmax(dL-d) = " << err << endl;
00488     if( force_inequality_error_check )
00489       handle_error(
00490         out,err,"max(dU-d)",feas_error_tol(),"feas_error_tol"
00491         ,feas_warning_tol(),"feas_error_tol",&test_failed
00492         );
00493 
00494     // ToDo: Update below code!
00495 /*
00496     if(out)
00497       *out
00498         << sep_line
00499         << "\nChecking nuL(i) * (dL - d)(i) = 0  ...\n";
00500     set_complementarity( *nu, u(), *d, opt_scale, lower, &c );
00501     if(out && print_vectors)
00502       *out
00503         << "nuL(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ) =\n" << c();
00504     if(out) {
00505       *out
00506         << "Comparing:\n"
00507         << "u(i) = nuL(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ), v = 0 ...\n";
00508     }
00509     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00510       , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00511       , print_all_warnings, out )) test_failed = true;
00512 */
00513 
00515     // d - dU <= 0
00516     if(out)
00517       *out
00518         << sep_line
00519         << "\nChecking d - dU <= 0 ...\n";
00520     V_VmV( u_d.get(), *d, *dU );
00521     if(out && print_vectors)
00522       *out
00523         << "d - dU =\n" << *u_d;
00524     Vt_S( u_d.get(), 1.0/(1.0+d_norm_inf) );
00525 
00526     err = max_element(*u_d);
00527     if(out)
00528       *out
00529         << "\nmax(d-dU) = " << err << endl;
00530     if( force_inequality_error_check )
00531       handle_error(
00532         out,err,"max(d-dU)",feas_error_tol(),"feas_error_tol"
00533         ,feas_warning_tol(),"feas_error_tol",&test_failed
00534         );
00535     // ToDo: Update below code!
00536 /*
00537     if(out)
00538       *out
00539         << sep_line
00540         << "\nChecking nuU(i) * (d - dU)(i) = 0  ...\n";
00541     set_complementarity( *nu, u(), *d, opt_scale, upper, &c );
00542     if(out && print_vectors)
00543       *out
00544         << "nuU(i) * (d - dU)(i) / ( 1 + |d(i)| + opt_scale ) =\n" << c();
00545     if(out) {
00546       *out
00547         << "Comparing:\n"
00548         << "u(i) = nuU(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ), v = 0 ...\n";
00549     }
00550     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00551              , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00552              , print_all_warnings, out )) test_failed = true;
00553 */
00554   }
00555 
00556   if( E ) {
00557 
00559     // e = op(E)*d + b*eta
00560     if(out)
00561       *out
00562         << sep_line
00563         << "\nComputing e = op(E)*d - b*eta ...\n";
00564     VectorSpace::vec_mut_ptr_t
00565       e   = ( trans_E == no_trans ? E->space_cols() : E->space_rows() ).create_member(),
00566       t_e = e->space().create_member();
00567     V_MtV( e.get(), *E, trans_E, *d );
00568     Vp_StV( e.get(), -(*eta), *b );
00569     e_norm_inf = e->norm_inf();
00570     if(out && print_vectors)
00571       *out
00572         << "e = op(E)*d - b*eta  =\n" << *e;
00573 
00575     // eL - e <= 0
00576     if(out)
00577       *out
00578         << sep_line
00579         << "\nChecking eL - e <= 0 ...\n";
00580     V_VmV( t_e.get(), *eL, *e );
00581     if(out && print_vectors)
00582       *out
00583         << "eL - e =\n" << *t_e;
00584     Vt_S( t_e.get(), 1.0/(1.0+e_norm_inf) );
00585 
00586     err = max_element(*t_e);
00587     if(out)
00588       *out
00589         << "\nmax(eL-e) = " << err << endl;
00590     if( force_inequality_error_check )
00591       handle_error(
00592         out,err,"max(eL-e)",feas_error_tol(),"feas_error_tol"
00593         ,feas_warning_tol(),"feas_error_tol",&test_failed
00594         );
00595     // ToDo: Update below code!
00596 /*
00597     if(out)
00598       *out
00599         << sep_line
00600         << "\nChecking muL(i) * (eL - e)(i) = 0  ...\n";
00601     set_complementarity( *mu, u(), e(), opt_scale, lower, &c );
00602     if(out && print_vectors)
00603       *out
00604         << "muL(i) * (eL - e)(i) / ( 1 + |e(i)| + opt_scale ) =\n" << c();
00605     if(out) {
00606       *out
00607         << "Comparing:\n"
00608         << "u(i) = muL(i) * (eL - e)(i) / ( 1 + |e(i)| + opt_scale ), v = 0 ...\n";
00609     }
00610     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00611              , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00612              , print_all_warnings, out )) test_failed = true;
00613 */
00614     
00616     // e - eU <= 0
00617     if(out)
00618       *out
00619         << sep_line
00620         << "\nChecking e - eU <= 0 ...\n";
00621     V_VmV( t_e.get(), *e, *eU );
00622     if(out && print_vectors)
00623       *out
00624         << "\ne - eU =\n" << *t_e;
00625     Vt_S( t_e.get(), 1.0/(1.0+e_norm_inf) );
00626 
00627     err = max_element(*t_e);
00628     if(out)
00629       *out
00630         << "\nmax(e-eU) = " << err << endl;
00631     if( force_inequality_error_check )
00632       handle_error(
00633         out,err,"max(e-eU)",feas_error_tol(),"feas_error_tol"
00634         ,feas_warning_tol(),"feas_error_tol",&test_failed
00635         );
00636     // ToDo: Update below code!
00637 /*
00638     if(out)
00639       *out
00640         << sep_line
00641         << "\nChecking muU(i) * (e - eU)(i) = 0  ...\n";
00642     set_complementarity( *mu, u(), e(), opt_scale, upper, &c );
00643     if(out && print_vectors)
00644       *out
00645         << "\nmuU(i) * (e - eU)(i) / ( 1 + |e(i)| + opt_scale ) =\n" << c();
00646     if(out) {
00647       *out
00648         << "\nComparing:\n"
00649         << "u(i) = muU(i) * (e - eU)(i) / ( 1 + |e(i)| + opt_scale )\n"
00650         << "v = 0 ...\n";
00651     }
00652     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00653              , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00654              , print_all_warnings, out )) test_failed = true;
00655 */
00656     
00657   }
00658 
00659   if( F ) {
00660     TEST_FOR_EXCEPT(true); // ToDo: Update below code!
00661 /*
00662 
00664     // r = - op(F)*d + eta * f 
00665     if(out)
00666       *out
00667         << sep_line
00668         << "\nComputing r = - op(F)*d + eta * f ...\n";
00669     V_StMtV( &r, -1.0, *F, trans_F, *d );
00670     Vp_StV( &r(), *eta, *f );
00671     if(out && print_vectors)
00672       *out
00673         << "\nr = - op(F)*d + eta * f =\n" << r();
00674 
00675     if(out) {
00676       *out
00677         << sep_line
00678         << "\nChecking r == f:\n"
00679         << "u = r, v = f ...\n";
00680     }
00681     if(!comp_v.comp( r(), *f, opt_warning_tol()
00682       , force_equality_error_check ? feas_error_tol() : really_big_error_tol
00683       , print_all_warnings, out )) test_failed = true;
00684 
00685 */
00686   }
00687 
00688   if(out) {
00689     *out
00690       << sep_line;
00691     if(solution_type != qps_t::SUBOPTIMAL_POINT) {
00692       if(test_failed) {
00693         *out
00694           << "\nDarn it!  At least one of the enforced QP optimality conditions were "
00695           << "not within the specified error tolerances!\n";
00696       }
00697       else {
00698         *out
00699           << "\nCongradulations!  All of the enforced QP optimality conditions were "
00700           << "within the specified error tolerances!\n";
00701       }
00702     }
00703     *out
00704       << "\n*** End checking QP optimality conditions ***\n";
00705   }
00706 
00707   return !test_failed;
00708 
00709 }
00710 
00711 } // end namespace ConstrainedOptPack

Generated on Wed May 12 21:52:29 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7