ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_QPSolverRelaxedTester.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <math.h>
43 
44 #include <limits>
45 
46 #include "ConstrainedOptPack_QPSolverRelaxedTester.hpp"
47 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
48 #include "AbstractLinAlgPack_VectorSpace.hpp"
49 #include "AbstractLinAlgPack_VectorMutable.hpp"
50 #include "AbstractLinAlgPack_VectorOut.hpp"
51 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
52 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
53 
54 namespace {
55 
56 //
57 const char* solution_type_str( ConstrainedOptPack::QPSolverStats::ESolutionType solution_type )
58 {
60  switch( solution_type ) {
61 
62  case qpst::OPTIMAL_SOLUTION:
63  return "OPTIMAL_SOLUTION";
64  case qpst::PRIMAL_FEASIBLE_POINT:
65  return "PRIMAL_FEASIBLE_POINT";
66  case qpst::DUAL_FEASIBLE_POINT:
67  return "DUAL_FEASIBLE_POINT";
68  case qpst::SUBOPTIMAL_POINT:
69  return "SUBOPTIMAL_POINT";
70  default:
72  }
73  return ""; // will never be executed.
74 }
75 
76 /* ToDo: Update this code!
77 
78 // Compute the scaled complementarity conditions.
79 //
80 // If uplo == upper then:
81 //
82 // / gamma(i) * constr_resid(i) / ( 1 + |constr(i)| + opt_scale ), for gamma(i) > 0
83 // comp_err(i) = |
84 // \ 0 otherwise
85 //
86 // If uplo == lower then:
87 //
88 // / gamma(i) * constr_resid(i) / ( 1 + |constr(i)| + opt_scale ), for gamma(i) < 0
89 // comp_err(i) = |
90 // \ 0 otherwise
91 //
92 //
93 void set_complementarity(
94  const AbstractLinAlgPack::SpVector &gamma
95  ,const DenseLinAlgPack::DVectorSlice &constr_resid
96  ,const DenseLinAlgPack::DVectorSlice &constr
97  ,const DenseLinAlgPack::value_type opt_scale
98  ,BLAS_Cpp::Uplo uplo
99  ,DenseLinAlgPack::DVector *comp_err
100  )
101 {
102  TEUCHOS_TEST_FOR_EXCEPT( !( gamma.size() == constr_resid.size() && gamma.size() == constr.size() ) );
103  comp_err->resize( gamma.size() );
104  *comp_err = 0.0;
105  const AbstractLinAlgPack::SpVector::difference_type o = gamma.offset();
106  if( gamma.nz() ) {
107  for( AbstractLinAlgPack::SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
108  const DenseLinAlgPack::size_type i = itr->indice() + o;
109  if( itr->value() > 0 && uplo == BLAS_Cpp::upper )
110  (*comp_err)(i) = itr->value() * constr_resid(i) / ( 1.0 + ::fabs(constr(i)) + opt_scale );
111  else if( itr->value() < 0 && uplo == BLAS_Cpp::lower )
112  (*comp_err)(i) = itr->value() * constr_resid(i) / ( 1.0 + ::fabs(constr(i)) + opt_scale );
113  }
114  }
115 }
116 
117 */
118 
119 // Handle the error reporting
120 void handle_error(
121  std::ostream *out
122  ,const AbstractLinAlgPack::value_type err
123  ,const char err_name[]
124  ,const AbstractLinAlgPack::value_type error_tol
125  ,const char error_tol_name[]
126  ,const AbstractLinAlgPack::value_type warning_tol
127  ,const char warning_tol_name[]
128  ,bool *test_failed
129  )
130 {
131  if( err >= error_tol ) {
132  if(out)
133  *out << "\n" << err_name << " = " << err << " >= " << error_tol_name << " = " << error_tol << std::endl;
134  *test_failed = true;
135  }
136  else if( err >= warning_tol ) {
137  if(out)
138  *out << "\n" << err_name << " = " << err << " >= " << warning_tol_name << " = " << warning_tol << std::endl;
139  }
140 }
141 
142 } // end namespace
143 
144 namespace ConstrainedOptPack {
145 
146 // public
147 
149  value_type opt_warning_tol
150  ,value_type opt_error_tol
151  ,value_type feas_warning_tol
152  ,value_type feas_error_tol
153  ,value_type comp_warning_tol
154  ,value_type comp_error_tol
155  )
156  :opt_warning_tol_(opt_warning_tol)
157  ,opt_error_tol_(opt_error_tol)
158  ,feas_warning_tol_(feas_warning_tol)
159  ,feas_error_tol_(feas_error_tol)
160  ,comp_warning_tol_(comp_warning_tol)
161  ,comp_error_tol_(comp_error_tol)
162 {}
163 
165  QPSolverStats::ESolutionType solution_type
166  ,const value_type infinite_bound
167  ,std::ostream* out, bool print_all_warnings, bool print_vectors
168  ,const Vector& g, const MatrixSymOp& G
169  ,value_type etaL
170  ,const Vector& dL, const Vector& dU
171  ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
172  ,const Vector& eL, const Vector& eU
173  ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
174  ,const value_type* obj_d
175  ,const value_type* eta, const Vector* d
176  ,const Vector* nu
177  ,const Vector* mu, const Vector* Ed
178  ,const Vector* lambda, const Vector* Fd
179  )
180 {
182  solution_type,infinite_bound,out,print_all_warnings,print_vectors
183  ,g,G,etaL,&dL,&dU,&E,trans_E,&b,&eL,&eU,&F,trans_F,&f
184  ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
185 }
186 
188  QPSolverStats::ESolutionType solution_type
189  ,const value_type infinite_bound
190  ,std::ostream* out, bool print_all_warnings, bool print_vectors
191  ,const Vector& g, const MatrixSymOp& G
192  ,value_type etaL
193  ,const Vector& dL, const Vector& dU
194  ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
195  ,const Vector& eL, const Vector& eU
196  ,const value_type* obj_d
197  ,const value_type* eta, const Vector* d
198  ,const Vector* nu
199  ,const Vector* mu, const Vector* Ed
200  )
201 {
203  solution_type,infinite_bound,out,print_all_warnings,print_vectors
204  ,g,G,etaL,&dL,&dU,&E,trans_E,&b,&eL,&eU,NULL,BLAS_Cpp::no_trans,NULL
205  ,obj_d,eta,d,nu,mu,Ed,NULL,NULL);
206 }
207 
209  QPSolverStats::ESolutionType solution_type
210  ,const value_type infinite_bound
211  ,std::ostream* out, bool print_all_warnings, bool print_vectors
212  ,const Vector& g, const MatrixSymOp& G
213  ,value_type etaL
214  ,const Vector& dL, const Vector& dU
215  ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
216  ,const value_type* obj_d
217  ,const value_type* eta, const Vector* d
218  ,const Vector* nu
219  ,const Vector* lambda, const Vector* Fd
220  )
221 {
223  solution_type,infinite_bound,out,print_all_warnings,print_vectors
224  ,g,G,etaL,&dL,&dU,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL,&F,trans_F,&f
225  ,obj_d,eta,d,nu,NULL,NULL,lambda,Fd );
226 }
227 
229  QPSolverStats::ESolutionType solution_type
230  ,const value_type infinite_bound
231  ,std::ostream* out, bool print_all_warnings, bool print_vectors
232  ,const Vector& g, const MatrixSymOp& G
233  ,const Vector& dL, const Vector& dU
234  ,const value_type* obj_d
235  ,const Vector* d
236  ,const Vector* nu
237  )
238 {
240  solution_type,infinite_bound,out,print_all_warnings,print_vectors
241  ,g,G,0.0,&dL,&dU,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL,NULL,BLAS_Cpp::no_trans,NULL
242  ,obj_d,NULL,d,nu,NULL,NULL,NULL,NULL);
243 }
244 
246  QPSolverStats::ESolutionType solution_type
247  ,const value_type infinite_bound
248  ,std::ostream* out, bool print_all_warnings, bool print_vectors
249  ,const Vector& g, const MatrixSymOp& G
250  ,value_type etaL
251  ,const Vector* dL, const Vector* dU
252  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
253  ,const Vector* eL, const Vector* eU
254  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
255  ,const value_type* obj_d
256  ,const value_type* eta, const Vector* d
257  ,const Vector* nu
258  ,const Vector* mu, const Vector* Ed
259  ,const Vector* lambda, const Vector* Fd
260  )
261 {
263  infinite_bound,g,G,etaL,dL,dU
264  ,E,trans_E,b,eL,eU,F,trans_F,f
265  ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
266 
268  solution_type,infinite_bound
269  ,out,print_all_warnings,print_vectors,g,G,etaL,dL,dU
270  ,E,trans_E,b,eL,eU,F,trans_F,f
271  ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
272 }
273 
274 // protected
275 
277  QPSolverStats::ESolutionType solution_type
278  ,const value_type infinite_bound
279  ,std::ostream* out, bool print_all_warnings, bool print_vectors
280  ,const Vector& g, const MatrixSymOp& G
281  ,value_type etaL
282  ,const Vector* dL, const Vector* dU
283  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
284  ,const Vector* eL, const Vector* eU
285  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
286  ,const value_type* obj_d
287  ,const value_type* eta, const Vector* d
288  ,const Vector* nu
289  ,const Vector* mu, const Vector* Ed
290  ,const Vector* lambda, const Vector* Fd
291  )
292 {
293  using std::endl;
294  using BLAS_Cpp::trans_not;
295  using BLAS_Cpp::no_trans;
296  using BLAS_Cpp::trans;
297  using BLAS_Cpp::upper;
298  using BLAS_Cpp::lower;
299  using LinAlgOpPack::sum;
300  using LinAlgOpPack::dot;
301  using LinAlgOpPack::Vt_S;
302  using LinAlgOpPack::V_VmV;
303  using LinAlgOpPack::Vp_StV;
304  using LinAlgOpPack::V_MtV;
305  using LinAlgOpPack::Vp_MtV;
306  using LinAlgOpPack::V_StMtV;
307  using LinAlgOpPack::Vp_V;
309  typedef QPSolverStats qps_t;
310 
311  bool test_failed = false;
312 
313  const size_type
314  nd = d->dim();
315 
316  const value_type
317  really_big_error_tol = std::numeric_limits<value_type>::max();
318 
319  value_type opt_scale = 0.0;
320  VectorSpace::vec_mut_ptr_t
321  t_d = d->space().create_member(),
322  u_d = d->space().create_member();
323 
324  value_type
325  err = 0,
326  d_norm_inf, // holds ||d||inf
327  e_norm_inf; // holds ||e||inf
328 
329  if(out)
330  *out
331  << "\n*** Begin checking QP optimality conditions ***\n"
332  << "\nThe solution type is " << solution_type_str(solution_type) << endl;
333 
334  bool force_opt_error_check
335  = solution_type==qps_t::OPTIMAL_SOLUTION || solution_type==qps_t::DUAL_FEASIBLE_POINT;
336  const bool force_inequality_error_check
337  = solution_type==qps_t::OPTIMAL_SOLUTION || solution_type==qps_t::PRIMAL_FEASIBLE_POINT;
338  const bool force_equality_error_check
339  = solution_type!=qps_t::SUBOPTIMAL_POINT;
340  const bool force_complementarity_error_check
341  = solution_type!=qps_t::SUBOPTIMAL_POINT;
342 
343  const char sep_line[] = "\n--------------------------------------------------------------------------------\n";
344 
346  // Checking d(L)/d(d) = 0
347  if(out)
348  *out
349  << sep_line
350  << "\nChecking d(L)/d(d) = g + G*d + nu + op(E)'*mu - op(F)'*lambda == 0 ...\n";
351 
352  if(out && !force_opt_error_check)
353  *out
354  << "The optimality error tolerance will not be enforced ...\n";
355 
356  opt_scale = 0.0;
357 
358  *u_d = g;
359  opt_scale += g.norm_inf();
360 
361  if(out) {
362  *out << "||g||inf = " << g.norm_inf() << endl;
363  }
364 
365  V_MtV( t_d.get(), G, no_trans, *d );
366  Vp_V( u_d.get(), *t_d );
367  opt_scale += t_d->norm_inf();
368 
369  if(out) {
370  *out << "||G*d||inf = " << t_d->norm_inf() << endl;
371  if(print_vectors)
372  *out << "g + G*d =\n" << *u_d;
373  }
374 
375  if( nu ) {
376  Vp_V( u_d.get(), *nu );
377  opt_scale += nu->norm_inf();
378  if(out)
379  *out << "||nu||inf = " << nu->norm_inf() << endl;
380  }
381 
382  if(E) {
383  V_MtV( t_d.get(), *E, trans_not(trans_E), *mu );
384  Vp_V( u_d.get(), *t_d );
385  opt_scale += t_d->norm_inf();
386  if(out) {
387  *out << "||op(E)'*mu||inf = " << t_d->norm_inf() << endl;
388  if(print_vectors)
389  *out << "op(E)'*mu =\n" << *t_d;
390  }
391  }
392 
393  if(F) {
394  V_MtV( t_d.get(), *F, trans_not(trans_F), *lambda );
395  Vp_V( u_d.get(), *t_d );
396  opt_scale += t_d->norm_inf();
397  if(out) {
398  *out << "||op(F)'*lambda||inf = " << t_d->norm_inf() << endl;
399  if(print_vectors)
400  *out << "op(F)'*lambda =\n" << *t_d;
401  }
402  }
403 
404  if( *eta > etaL ) { // opt_scale + |(eta - etaL) * (b'*mu + f'*lambda)|
405  const value_type
406  term = ::fabs( (*eta - etaL) * (E ? dot(*b,*mu) : 0.0) + (F ? dot(*f,*lambda) : 0.0 ) );
407  if(out) {
408  *out << "|(eta - etaL) * (b'*mu + f'*lambda)| = " << term << endl;
409  }
410  opt_scale += term;
411  }
412 
413  if(out && print_vectors)
414  *out
415  << "g + G*d + nu + op(E)'*mu - op(F)'*lambda =\n" << *u_d;
416 
417  Vt_S( u_d.get(), 1.0/(1.0+opt_scale) );
418 
419  err = sum( *u_d );
420 
421  if(out)
422  *out
423  << "\nopt_scale = " << opt_scale << endl
424  << "opt_err = sum( | g + G*d + nu + op(E)'*mu - op(F)'*lambda | / (1+opt_scale) ) / nd\n"
425  << " = " << err << " / " << nd << " = " << (err/nd) << endl;
426 
427  err *= nd;
428 
429  if( force_opt_error_check ) {
430  if( err >= opt_error_tol() ) {
431  if(out)
432  *out << "\nopt_err = " << err << " >= opt_error_tol = " << opt_error_tol() << endl;
433  test_failed = true;
434  }
435  else if( err >= opt_warning_tol() ) {
436  if(out)
437  *out << "\nopt_err = " << err << " >= opt_error_tol = " << opt_error_tol() << endl;
438  }
439  }
440 
441  if(out) {
442  *out
443  << sep_line
444  << "\nTesting feasibility of the constraints and the complementarity conditions ...\n";
445  if(!force_inequality_error_check)
446  *out
447  << "The inequality feasibility error tolerance will not be enforced ...\n";
448  if(!force_equality_error_check)
449  *out
450  << "The equality feasibility error tolerance will not be enforced ...\n";
451  if(!force_complementarity_error_check)
452  *out
453  << "The complementarity conditions error tolerance will not be enforced ...\n";
454  }
455 
457  // etaL - eta
458  if(out)
459  *out
460  << sep_line
461  << "\nChecking etaL - eta <= 0 ...\n";
462  if(out)
463  *out
464  << "etaL - eta = " << (etaL - (*eta)) << endl;
465  if( etaL - (*eta) > feas_warning_tol() ) {
466  if(out)
467  *out
468  << "Warning, etaL - eta = " << etaL << " - " << (*eta)
469  << " = " << (etaL - (*eta)) << " > feas_warning_tol = "
470  << feas_warning_tol() << endl;
471  }
472  if( force_inequality_error_check && etaL - (*eta) > feas_error_tol() ) {
473  if(out)
474  *out
475  << "Error, etaL - eta = " << etaL << " - " << (*eta)
476  << " = " << (etaL - (*eta)) << " > feas_error_tol = "
477  << feas_error_tol() << endl;
478  test_failed = true;
479  }
480 
481  d_norm_inf = d->norm_inf();
482 
483  if(dL) {
484 
486  // dL - d <= 0
487  if(out)
488  *out
489  << sep_line
490  << "\nChecking dL - d <= 0 ...\n";
491  V_VmV( u_d.get(), *dL, *d );
492  if(out && print_vectors)
493  *out
494  << "dL - d =\n" << *u_d;
495  Vt_S( u_d.get(), 1.0/(1.0+d_norm_inf) );
496 
497  err = max_element(*u_d);
498  if(out)
499  *out
500  << "\nmax(dL-d) = " << err << endl;
501  if( force_inequality_error_check )
502  handle_error(
503  out,err,"max(dU-d)",feas_error_tol(),"feas_error_tol"
504  ,feas_warning_tol(),"feas_error_tol",&test_failed
505  );
506 
507  // ToDo: Update below code!
508 /*
509  if(out)
510  *out
511  << sep_line
512  << "\nChecking nuL(i) * (dL - d)(i) = 0 ...\n";
513  set_complementarity( *nu, u(), *d, opt_scale, lower, &c );
514  if(out && print_vectors)
515  *out
516  << "nuL(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ) =\n" << c();
517  if(out) {
518  *out
519  << "Comparing:\n"
520  << "u(i) = nuL(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ), v = 0 ...\n";
521  }
522  if(!comp_v.comp( c(), 0.0, opt_warning_tol()
523  , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
524  , print_all_warnings, out )) test_failed = true;
525 */
526 
528  // d - dU <= 0
529  if(out)
530  *out
531  << sep_line
532  << "\nChecking d - dU <= 0 ...\n";
533  V_VmV( u_d.get(), *d, *dU );
534  if(out && print_vectors)
535  *out
536  << "d - dU =\n" << *u_d;
537  Vt_S( u_d.get(), 1.0/(1.0+d_norm_inf) );
538 
539  err = max_element(*u_d);
540  if(out)
541  *out
542  << "\nmax(d-dU) = " << err << endl;
543  if( force_inequality_error_check )
544  handle_error(
545  out,err,"max(d-dU)",feas_error_tol(),"feas_error_tol"
546  ,feas_warning_tol(),"feas_error_tol",&test_failed
547  );
548  // ToDo: Update below code!
549 /*
550  if(out)
551  *out
552  << sep_line
553  << "\nChecking nuU(i) * (d - dU)(i) = 0 ...\n";
554  set_complementarity( *nu, u(), *d, opt_scale, upper, &c );
555  if(out && print_vectors)
556  *out
557  << "nuU(i) * (d - dU)(i) / ( 1 + |d(i)| + opt_scale ) =\n" << c();
558  if(out) {
559  *out
560  << "Comparing:\n"
561  << "u(i) = nuU(i) * (dL - d)(i) / ( 1 + |d(i)| + opt_scale ), v = 0 ...\n";
562  }
563  if(!comp_v.comp( c(), 0.0, opt_warning_tol()
564  , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
565  , print_all_warnings, out )) test_failed = true;
566 */
567  }
568 
569  if( E ) {
570 
572  // e = op(E)*d + b*eta
573  if(out)
574  *out
575  << sep_line
576  << "\nComputing e = op(E)*d - b*eta ...\n";
577  VectorSpace::vec_mut_ptr_t
578  e = ( trans_E == no_trans ? E->space_cols() : E->space_rows() ).create_member(),
579  t_e = e->space().create_member();
580  V_MtV( e.get(), *E, trans_E, *d );
581  Vp_StV( e.get(), -(*eta), *b );
582  e_norm_inf = e->norm_inf();
583  if(out && print_vectors)
584  *out
585  << "e = op(E)*d - b*eta =\n" << *e;
586 
588  // eL - e <= 0
589  if(out)
590  *out
591  << sep_line
592  << "\nChecking eL - e <= 0 ...\n";
593  V_VmV( t_e.get(), *eL, *e );
594  if(out && print_vectors)
595  *out
596  << "eL - e =\n" << *t_e;
597  Vt_S( t_e.get(), 1.0/(1.0+e_norm_inf) );
598 
599  err = max_element(*t_e);
600  if(out)
601  *out
602  << "\nmax(eL-e) = " << err << endl;
603  if( force_inequality_error_check )
604  handle_error(
605  out,err,"max(eL-e)",feas_error_tol(),"feas_error_tol"
606  ,feas_warning_tol(),"feas_error_tol",&test_failed
607  );
608  // ToDo: Update below code!
609 /*
610  if(out)
611  *out
612  << sep_line
613  << "\nChecking muL(i) * (eL - e)(i) = 0 ...\n";
614  set_complementarity( *mu, u(), e(), opt_scale, lower, &c );
615  if(out && print_vectors)
616  *out
617  << "muL(i) * (eL - e)(i) / ( 1 + |e(i)| + opt_scale ) =\n" << c();
618  if(out) {
619  *out
620  << "Comparing:\n"
621  << "u(i) = muL(i) * (eL - e)(i) / ( 1 + |e(i)| + opt_scale ), v = 0 ...\n";
622  }
623  if(!comp_v.comp( c(), 0.0, opt_warning_tol()
624  , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
625  , print_all_warnings, out )) test_failed = true;
626 */
627 
629  // e - eU <= 0
630  if(out)
631  *out
632  << sep_line
633  << "\nChecking e - eU <= 0 ...\n";
634  V_VmV( t_e.get(), *e, *eU );
635  if(out && print_vectors)
636  *out
637  << "\ne - eU =\n" << *t_e;
638  Vt_S( t_e.get(), 1.0/(1.0+e_norm_inf) );
639 
640  err = max_element(*t_e);
641  if(out)
642  *out
643  << "\nmax(e-eU) = " << err << endl;
644  if( force_inequality_error_check )
645  handle_error(
646  out,err,"max(e-eU)",feas_error_tol(),"feas_error_tol"
647  ,feas_warning_tol(),"feas_error_tol",&test_failed
648  );
649  // ToDo: Update below code!
650 /*
651  if(out)
652  *out
653  << sep_line
654  << "\nChecking muU(i) * (e - eU)(i) = 0 ...\n";
655  set_complementarity( *mu, u(), e(), opt_scale, upper, &c );
656  if(out && print_vectors)
657  *out
658  << "\nmuU(i) * (e - eU)(i) / ( 1 + |e(i)| + opt_scale ) =\n" << c();
659  if(out) {
660  *out
661  << "\nComparing:\n"
662  << "u(i) = muU(i) * (e - eU)(i) / ( 1 + |e(i)| + opt_scale )\n"
663  << "v = 0 ...\n";
664  }
665  if(!comp_v.comp( c(), 0.0, opt_warning_tol()
666  , force_complementarity_error_check ? comp_error_tol() : really_big_error_tol
667  , print_all_warnings, out )) test_failed = true;
668 */
669 
670  }
671 
672  if( F ) {
673  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Update below code!
674 /*
675 
677  // r = - op(F)*d + eta * f
678  if(out)
679  *out
680  << sep_line
681  << "\nComputing r = - op(F)*d + eta * f ...\n";
682  V_StMtV( &r, -1.0, *F, trans_F, *d );
683  Vp_StV( &r(), *eta, *f );
684  if(out && print_vectors)
685  *out
686  << "\nr = - op(F)*d + eta * f =\n" << r();
687 
688  if(out) {
689  *out
690  << sep_line
691  << "\nChecking r == f:\n"
692  << "u = r, v = f ...\n";
693  }
694  if(!comp_v.comp( r(), *f, opt_warning_tol()
695  , force_equality_error_check ? feas_error_tol() : really_big_error_tol
696  , print_all_warnings, out )) test_failed = true;
697 
698 */
699  }
700 
701  if(out) {
702  *out
703  << sep_line;
704  if(solution_type != qps_t::SUBOPTIMAL_POINT) {
705  if(test_failed) {
706  *out
707  << "\nDarn it! At least one of the enforced QP optimality conditions were "
708  << "not within the specified error tolerances!\n";
709  }
710  else {
711  *out
712  << "\nCongradulations! All of the enforced QP optimality conditions were "
713  << "within the specified error tolerances!\n";
714  }
715  }
716  *out
717  << "\n*** End checking QP optimality conditions ***\n";
718  }
719 
720  return !test_failed;
721 
722 }
723 
724 } // end namespace ConstrainedOptPack
virtual bool imp_check_optimality_conditions(QPSolverStats::ESolutionType solution_type, const value_type infinite_bound, std::ostream *out, bool print_all_warnings, bool print_vectors, const Vector &g, const MatrixSymOp &G, value_type etaL, const Vector *dL, const Vector *dU, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const Vector *eL, const Vector *eU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, const value_type *obj_d, const value_type *eta, const Vector *d, const Vector *nu, const Vector *mu, const Vector *Ed, const Vector *lambda, const Vector *Fd)
Subclasses are to override this to implement the testing code.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
value_type max_element(const Vector &v)
virtual bool check_optimality_conditions(QPSolverStats::ESolutionType solution_type, const value_type infinite_bound, std::ostream *out, bool print_all_warnings, bool print_vectors, const Vector &g, const MatrixSymOp &G, value_type etaL, const Vector &dL, const Vector &dU, const MatrixOp &E, BLAS_Cpp::Transp trans_E, const Vector &b, const Vector &eL, const Vector &eU, const MatrixOp &F, BLAS_Cpp::Transp trans_F, const Vector &f, const value_type *obj_d, const value_type *eta, const Vector *d, const Vector *nu, const Vector *mu, const Vector *Ed, const Vector *lambda, const Vector *Fd)
Check the optimality conditions for the solved (or partially solved) QP.
size_t size_type
Class for storing statistics about a run of a (active set?) QP solver.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
Transp trans_not(Transp _trans)
static void validate_input(const value_type infinite_bound, const Vector &g, const MatrixSymOp &G, value_type etaL, const Vector *dL, const Vector *dU, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const Vector *eL, const Vector *eU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, const value_type *obj_d, const value_type *eta, const Vector *d, const Vector *nu, const Vector *mu, const Vector *Ed, const Vector *lambda, const Vector *Fd)
This is a (static) function that is used as a utility to validate the input arguments to solve_qp()...
value_type sum(const Vector &v_rhs)
QPSolverRelaxedTester(value_type opt_warning_tol=1e-10, value_type opt_error_tol=1e-5, value_type feas_warning_tol=1e-10, value_type feas_error_tol=1e-5, value_type comp_warning_tol=1e-10, value_type comp_error_tol=1e-5)
Transp
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
ESolutionType
Enumeration for the type of point returned from solve_qp(...).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)