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