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   assert( 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   }
00366 
00367   if(E) {
00368     V_MtV( t_d.get(), *E, trans_not(trans_E), *mu );
00369     Vp_V( u_d.get(), *t_d );
00370     opt_scale += t_d->norm_inf();
00371     if(out) {
00372       *out << "||op(E)'*mu||inf = " << t_d->norm_inf() << endl;
00373       if(print_vectors)
00374         *out << "op(E)'*mu =\n" << *t_d;
00375     }
00376   }     
00377 
00378   if(F) {
00379     V_MtV( t_d.get(), *F, trans_not(trans_F), *lambda );
00380     Vp_V( u_d.get(), *t_d );
00381     opt_scale += t_d->norm_inf();
00382     if(out) {
00383       *out << "||op(F)'*lambda||inf = " << t_d->norm_inf() << endl;
00384       if(print_vectors)
00385         *out << "op(F)'*lambda =\n" << *t_d;
00386     }
00387   }
00388 
00389   if( *eta > etaL ) { // opt_scale + |(eta - etaL) * (b'*mu + f'*lambda)|
00390     const value_type
00391       term = ::fabs( (*eta - etaL) * (E ? dot(*b,*mu) : 0.0) + (F ? dot(*f,*lambda) : 0.0 ) );
00392     if(out) {
00393       *out << "|(eta - etaL) * (b'*mu + f'*lambda)| = " << term << endl;
00394     }
00395     opt_scale += term;
00396   }
00397 
00398   if(out && print_vectors)
00399     *out
00400       << "g + G*d + nu + op(E)'*mu - op(F)'*lambda =\n" << *u_d;
00401 
00402   Vt_S( u_d.get(), 1.0/(1.0+opt_scale) );
00403 
00404   err = sum( *u_d );
00405 
00406   if(out)
00407     *out
00408       << "\nopt_scale = " << opt_scale << endl
00409       << "opt_err = sum( | g + G*d + nu + op(E)'*mu - op(F)'*lambda | / (1+opt_scale) ) / nd\n"
00410       << "        = " << err << " / " << nd << " = " << (err/nd) << endl;
00411 
00412   err *= nd;
00413 
00414   if( force_opt_error_check ) {
00415     if( err >= opt_error_tol() ) {
00416       if(out)
00417         *out << "\nopt_err = " << err << " >= opt_error_tol = " << opt_error_tol() << endl;
00418       test_failed = true;
00419     }
00420     else if( err >= opt_warning_tol() ) {
00421       if(out)
00422         *out << "\nopt_err = " << err << " >= opt_error_tol = " << opt_error_tol() << endl;
00423     }
00424   }
00425 
00426   if(out) {
00427     *out
00428       << sep_line
00429       << "\nTesting feasibility of the constraints and the complementarity conditions ...\n";
00430     if(!force_inequality_error_check)
00431       *out
00432         << "The inequality feasibility error tolerance will not be enforced ...\n";
00433     if(!force_equality_error_check)
00434       *out
00435         << "The equality feasibility error tolerance will not be enforced ...\n";
00436     if(!force_complementarity_error_check)
00437       *out
00438         << "The complementarity conditions error tolerance will not be enforced ...\n";
00439   }
00440 
00442   // etaL - eta
00443   if(out)
00444     *out
00445       << sep_line
00446       << "\nChecking etaL - eta <= 0 ...\n";
00447   if(out)
00448     *out
00449       << "etaL - eta = " << (etaL - (*eta)) << endl;
00450   if( etaL - (*eta) > feas_warning_tol() ) {
00451     if(out)
00452       *out
00453         << "Warning, etaL - eta = " << etaL << " - " << (*eta)
00454         << " = " << (etaL - (*eta)) << " >  feas_warning_tol = "
00455         <<  feas_warning_tol() << endl;
00456   }
00457   if( force_inequality_error_check && etaL - (*eta) > feas_error_tol() ) {
00458     if(out)
00459       *out
00460         << "Error, etaL - eta = " << etaL << " - " << (*eta)
00461         << " = " << (etaL - (*eta)) << " >  feas_error_tol = "
00462         <<  feas_error_tol() << endl;
00463     test_failed = true;
00464   } 
00465 
00466   d_norm_inf = d->norm_inf();
00467 
00468   if(dL) {
00469 
00471     // dL - d <= 0
00472     if(out)
00473       *out
00474         << sep_line
00475         << "\nChecking dL - d <= 0 ...\n";
00476     V_VmV( u_d.get(), *dL, *d );
00477     if(out && print_vectors)
00478       *out
00479         << "dL - d =\n" << *u_d;
00480     Vt_S( u_d.get(), 1.0/(1.0+d_norm_inf) );
00481 
00482     err = max_element(*u_d);
00483     if(out)
00484       *out
00485         << "\nmax(dL-d) = " << err << endl;
00486     if( force_inequality_error_check )
00487       handle_error(
00488         out,err,"max(dU-d)",feas_error_tol(),"feas_error_tol"
00489         ,feas_warning_tol(),"feas_error_tol",&test_failed
00490         );
00491 
00492     // ToDo: Update below code!
00493 /*
00494     if(out)
00495       *out
00496         << sep_line
00497         << "\nChecking nuL(i) * (dL - d)(i) = 0  ...\n";
00498     set_complementarity( *nu, u(), *d, opt_scale, lower, &c );
00499     if(out && print_vectors)
00500       *out
00501         << "nuL(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ) =\n" << c();
00502     if(out) {
00503       *out
00504         << "Comparing:\n"
00505         << "u(i) = nuL(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ), v = 0 ...\n";
00506     }
00507     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00508       , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00509       , print_all_warnings, out )) test_failed = true;
00510 */
00511 
00513     // d - dU <= 0
00514     if(out)
00515       *out
00516         << sep_line
00517         << "\nChecking d - dU <= 0 ...\n";
00518     V_VmV( u_d.get(), *d, *dU );
00519     if(out && print_vectors)
00520       *out
00521         << "d - dU =\n" << *u_d;
00522     Vt_S( u_d.get(), 1.0/(1.0+d_norm_inf) );
00523 
00524     err = max_element(*u_d);
00525     if(out)
00526       *out
00527         << "\nmax(d-dU) = " << err << endl;
00528     if( force_inequality_error_check )
00529       handle_error(
00530         out,err,"max(d-dU)",feas_error_tol(),"feas_error_tol"
00531         ,feas_warning_tol(),"feas_error_tol",&test_failed
00532         );
00533     // ToDo: Update below code!
00534 /*
00535     if(out)
00536       *out
00537         << sep_line
00538         << "\nChecking nuU(i) * (d - dU)(i) = 0  ...\n";
00539     set_complementarity( *nu, u(), *d, opt_scale, upper, &c );
00540     if(out && print_vectors)
00541       *out
00542         << "nuU(i) * (d - dU)(i) / ( 1 + |d(i)| + opt_scale ) =\n" << c();
00543     if(out) {
00544       *out
00545         << "Comparing:\n"
00546         << "u(i) = nuU(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ), v = 0 ...\n";
00547     }
00548     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00549              , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00550              , print_all_warnings, out )) test_failed = true;
00551 */
00552   }
00553 
00554   if( E ) {
00555 
00557     // e = op(E)*d + b*eta
00558     if(out)
00559       *out
00560         << sep_line
00561         << "\nComputing e = op(E)*d - b*eta ...\n";
00562     VectorSpace::vec_mut_ptr_t
00563       e   = ( trans_E == no_trans ? E->space_cols() : E->space_rows() ).create_member(),
00564       t_e = e->space().create_member();
00565     V_MtV( e.get(), *E, trans_E, *d );
00566     Vp_StV( e.get(), -(*eta), *b );
00567     e_norm_inf = e->norm_inf();
00568     if(out && print_vectors)
00569       *out
00570         << "e = op(E)*d - b*eta  =\n" << *e;
00571 
00573     // eL - e <= 0
00574     if(out)
00575       *out
00576         << sep_line
00577         << "\nChecking eL - e <= 0 ...\n";
00578     V_VmV( t_e.get(), *eL, *e );
00579     if(out && print_vectors)
00580       *out
00581         << "eL - e =\n" << *t_e;
00582     Vt_S( t_e.get(), 1.0/(1.0+e_norm_inf) );
00583 
00584     err = max_element(*t_e);
00585     if(out)
00586       *out
00587         << "\nmax(eL-e) = " << err << endl;
00588     if( force_inequality_error_check )
00589       handle_error(
00590         out,err,"max(eL-e)",feas_error_tol(),"feas_error_tol"
00591         ,feas_warning_tol(),"feas_error_tol",&test_failed
00592         );
00593     // ToDo: Update below code!
00594 /*
00595     if(out)
00596       *out
00597         << sep_line
00598         << "\nChecking muL(i) * (eL - e)(i) = 0  ...\n";
00599     set_complementarity( *mu, u(), e(), opt_scale, lower, &c );
00600     if(out && print_vectors)
00601       *out
00602         << "muL(i) * (eL - e)(i) / ( 1 + |e(i)| + opt_scale ) =\n" << c();
00603     if(out) {
00604       *out
00605         << "Comparing:\n"
00606         << "u(i) = muL(i) * (eL - e)(i) / ( 1 + |e(i)| + opt_scale ), v = 0 ...\n";
00607     }
00608     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00609              , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00610              , print_all_warnings, out )) test_failed = true;
00611 */
00612     
00614     // e - eU <= 0
00615     if(out)
00616       *out
00617         << sep_line
00618         << "\nChecking e - eU <= 0 ...\n";
00619     V_VmV( t_e.get(), *e, *eU );
00620     if(out && print_vectors)
00621       *out
00622         << "\ne - eU =\n" << *t_e;
00623     Vt_S( t_e.get(), 1.0/(1.0+e_norm_inf) );
00624 
00625     err = max_element(*t_e);
00626     if(out)
00627       *out
00628         << "\nmax(e-eU) = " << err << endl;
00629     if( force_inequality_error_check )
00630       handle_error(
00631         out,err,"max(e-eU)",feas_error_tol(),"feas_error_tol"
00632         ,feas_warning_tol(),"feas_error_tol",&test_failed
00633         );
00634     // ToDo: Update below code!
00635 /*
00636     if(out)
00637       *out
00638         << sep_line
00639         << "\nChecking muU(i) * (e - eU)(i) = 0  ...\n";
00640     set_complementarity( *mu, u(), e(), opt_scale, upper, &c );
00641     if(out && print_vectors)
00642       *out
00643         << "\nmuU(i) * (e - eU)(i) / ( 1 + |e(i)| + opt_scale ) =\n" << c();
00644     if(out) {
00645       *out
00646         << "\nComparing:\n"
00647         << "u(i) = muU(i) * (e - eU)(i) / ( 1 + |e(i)| + opt_scale )\n"
00648         << "v = 0 ...\n";
00649     }
00650     if(!comp_v.comp( c(), 0.0, opt_warning_tol()
00651              , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
00652              , print_all_warnings, out )) test_failed = true;
00653 */
00654     
00655   }
00656 
00657   if( F ) {
00658     TEST_FOR_EXCEPT(true); // ToDo: Update below code!
00659 /*
00660 
00662     // r = - op(F)*d + eta * f 
00663     if(out)
00664       *out
00665         << sep_line
00666         << "\nComputing r = - op(F)*d + eta * f ...\n";
00667     V_StMtV( &r, -1.0, *F, trans_F, *d );
00668     Vp_StV( &r(), *eta, *f );
00669     if(out && print_vectors)
00670       *out
00671         << "\nr = - op(F)*d + eta * f =\n" << r();
00672 
00673     if(out) {
00674       *out
00675         << sep_line
00676         << "\nChecking r == f:\n"
00677         << "u = r, v = f ...\n";
00678     }
00679     if(!comp_v.comp( r(), *f, opt_warning_tol()
00680       , force_equality_error_check ? feas_error_tol() : really_big_error_tol
00681       , print_all_warnings, out )) test_failed = true;
00682 
00683 */
00684   }
00685 
00686   if(out) {
00687     *out
00688       << sep_line;
00689     if(solution_type != qps_t::SUBOPTIMAL_POINT) {
00690       if(test_failed) {
00691         *out
00692           << "\nDarn it!  At least one of the enforced QP optimality conditions were "
00693           << "not within the specified error tolerances!\n";
00694       }
00695       else {
00696         *out
00697           << "\nCongradulations!  All of the enforced QP optimality conditions were "
00698           << "within the specified error tolerances!\n";
00699       }
00700     }
00701     *out
00702       << "\n*** End checking QP optimality conditions ***\n";
00703   }
00704 
00705   return !test_failed;
00706 
00707 }
00708 
00709 } // end namespace ConstrainedOptPack

Generated on Thu Sep 18 12:35:15 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1