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_QPSolverRelaxed.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 <assert.h>
43 
44 #include <limits>
45 
46 #include "ConstrainedOptPack_QPSolverRelaxed.hpp"
47 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
48 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
49 #include "AbstractLinAlgPack_VectorMutable.hpp"
50 #include "AbstractLinAlgPack_VectorOut.hpp"
51 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
52 #include "ProfileHackPack_profile_hack.hpp"
53 #include "Teuchos_Assert.hpp"
54 
55 namespace ConstrainedOptPack {
56 
57 // public
58 
60  :infinite_bound_(std::numeric_limits<value_type>::max())
61 {}
62 
65  std::ostream* out, EOutputLevel olevel, ERunTests test_what
66  ,const Vector& g, const MatrixSymOp& G
67  ,value_type etaL
68  ,const Vector& dL, const Vector& dU
69  ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
70  ,const Vector& eL, const Vector& eU
71  ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
72  ,value_type* obj_d
73  ,value_type* eta, VectorMutable* d
74  ,VectorMutable* nu
75  ,VectorMutable* mu, VectorMutable* Ed
76  ,VectorMutable* lambda, VectorMutable* Fd
77  )
78 {
79  return solve_qp(out,olevel,test_what,g,G,etaL,&dL,&dU
80  ,&E,trans_E,&b,&eL,&eU,&F,trans_F,&f
81  ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
82 }
83 
86  std::ostream* out, EOutputLevel olevel, ERunTests test_what
87  ,const Vector& g, const MatrixSymOp& G
88  ,value_type etaL
89  ,const Vector& dL, const Vector& dU
90  ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
91  ,const Vector& eL, const Vector& eU
92  ,value_type* obj_d
93  ,value_type* eta, VectorMutable* d
94  ,VectorMutable* nu
95  ,VectorMutable* mu, VectorMutable* Ed
96  )
97 {
98  return solve_qp(out,olevel,test_what,g,G,etaL,&dL,&dU
99  ,&E,trans_E,&b,&eL,&eU,NULL,BLAS_Cpp::no_trans,NULL
100  ,obj_d,eta,d,nu,mu,Ed,NULL,NULL);
101 }
102 
105  std::ostream* out, EOutputLevel olevel, ERunTests test_what
106  ,const Vector& g, const MatrixSymOp& G
107  ,value_type etaL
108  ,const Vector& dL, const Vector& dU
109  ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
110  ,value_type* obj_d
111  ,value_type* eta, VectorMutable* d
112  ,VectorMutable* nu
113  ,VectorMutable* lambda, VectorMutable* Fd
114  )
115 {
116  return solve_qp(out,olevel,test_what,g,G,etaL,&dL,&dU
117  ,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL,&F,trans_F,&f
118  ,obj_d,eta,d,nu,NULL,NULL,lambda,Fd);
119 }
120 
121 
124  std::ostream* out, EOutputLevel olevel, ERunTests test_what
125  ,const Vector& g, const MatrixSymOp& G
126  ,const Vector& dL, const Vector& dU
127  ,value_type* obj_d
128  ,VectorMutable* d
129  ,VectorMutable* nu
130  )
131 {
132  value_type eta = 0, etaL = 0;
133  return solve_qp(
134  out,olevel,test_what,g,G,etaL,&dL,&dU
135  ,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL
136  ,NULL,BLAS_Cpp::no_trans,NULL
137  ,obj_d,&eta,d,nu,NULL,NULL,NULL,NULL);
138 }
139 
142  std::ostream* out, EOutputLevel olevel, ERunTests test_what
143  ,const Vector& g, const MatrixSymOp& G
144  ,value_type etaL
145  ,const Vector* dL, const Vector* dU
146  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
147  ,const Vector* eL, const Vector* eU
148  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
149  ,value_type* obj_d
150  ,value_type* eta, VectorMutable* d
151  ,VectorMutable* nu
152  ,VectorMutable* mu, VectorMutable* Ed
153  ,VectorMutable* lambda, VectorMutable* Fd
154  )
155 {
156 #ifdef PROFILE_HACK_ENABLED
157  ProfileHackPack::ProfileTiming profile_timing( "QPSolverRelaxed::solve_qp(...)" );
158 #endif
160  infinite_bound(),g,G,etaL,dL,dU
161  ,E,trans_E,b,eL,eU,F,trans_F,f
162  ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
164  infinite_bound(),out,olevel,g,G,etaL,dL,dU,E,trans_E,b,eL,eU
165  ,F,trans_F,f,eta,d,nu,mu,lambda );
167  solve_return = imp_solve_qp(
168  out,olevel,test_what,g,G,etaL,dL,dU
169  ,E,trans_E,b,eL,eU,F,trans_F,f
170  ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
172  infinite_bound(),out,olevel,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
173  return solve_return;
174 }
175 
177  const value_type infinite_bound
178  ,const Vector& g, const MatrixSymOp& G
179  ,value_type etaL
180  ,const Vector* dL, const Vector* dU
181  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
182  ,const Vector* eL, const Vector* eU
183  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
184  ,const value_type* obj_d
185  ,const value_type* eta, const Vector* d
186  ,const Vector* nu
187  ,const Vector* mu, const Vector* Ed
188  ,const Vector* lambda, const Vector* Fd
189  )
190 {
191  // Validate output arguments
193  !d, std::invalid_argument
194  ,"QPSolverRelaxed::validate_input(...) : Error, "
195  "If d!=NULL is not allowed." );
197  ( E || F ) && !eta, std::invalid_argument
198  ,"QPSolverRelaxed::validate_input(...) : Error, "
199  "If eta!=NULL is not allowed if E!=NULL or F!=NULL." );
200 
201  // Validate the sets of constraints arguments
203  dL && ( !dU || !nu ), std::invalid_argument
204  ,"QPSolverRelaxed::validate_input(...) : Error, "
205  "If dL!=NULL then dU!=NULL and nu!=NULL must also be true." );
207  E && ( !b || !eL || !eU || !mu ), std::invalid_argument
208  ,"QPSolverRelaxed::validate_input(...) : Error, "
209  "If E!=NULL then b!=NULL, eL!=NULL, eU!=NULL and mu!=NULL must also "
210  "be true." );
212  F && ( !f || !lambda ), std::invalid_argument
213  ,"QPSolverRelaxed::validate_input(...) : Error, "
214  "If F!=NULL then f!=NULL and lambda!=NULL must also "
215  "be true." );
216 
217  // ToDo: Replace the below code with checks of compatibility of the vector spaces!
218 
219  // Validate the sizes of the arguments
220  const size_type
221  nd = d->dim();
223  g.dim() != nd, std::invalid_argument
224  ,"QPSolverRelaxed::validate_input(...) : Error, "
225  "g.dim() != d->dim()." );
227  G.rows() != nd || G.cols() != nd, std::invalid_argument
228  ,"QPSolverRelaxed::validate_input(...) : Error, "
229  "G.rows() != d->dim() or G.cols() != d->dim()." );
231  dL && dL->dim() != nd, std::invalid_argument
232  ,"QPSolverRelaxed::validate_input(...) : Error, "
233  "dL->dim() = " << dL->dim() << " != d->dim() = " << nd << "." );
235  dU && dU->dim() != nd, std::invalid_argument
236  ,"QPSolverRelaxed::validate_input(...) : Error, "
237  "dU->dim() = " << dU->dim() << " != d->dim() = " << nd << "." );
238  if( E ) {
239  const size_type
240  m_in = BLAS_Cpp::rows( E->rows(), E->cols(), trans_E );
242  BLAS_Cpp::cols( E->rows(), E->cols(), trans_E ) != nd, std::invalid_argument
243  ,"QPSolverRelaxed::validate_input(...) : Error, op(E).cols() != d->dim()." );
244  TEUCHOS_TEST_FOR_EXCEPTION(
245  b->dim() != m_in, std::invalid_argument
246  ,"QPSolverRelaxed::validate_input(...) : Error, b->dim() != op(E).rows()." );
247  TEUCHOS_TEST_FOR_EXCEPTION(
248  eL->dim() != m_in, std::invalid_argument
249  ,"QPSolverRelaxed::validate_input(...) : Error, eL->dim() != op(E).rows()." );
250  TEUCHOS_TEST_FOR_EXCEPTION(
251  eU->dim() != m_in, std::invalid_argument
252  ,"QPSolverRelaxed::validate_input(...) : Error, eU->dim() != op(E).rows()." );
253  TEUCHOS_TEST_FOR_EXCEPTION(
254  Ed && Ed->dim() != m_in, std::invalid_argument
255  ,"QPSolverRelaxed::validate_input(...) : Error, Ed->dim() != op(E).rows()." );
256  }
257  if( F ) {
258  const size_type
259  m_eq = BLAS_Cpp::rows( F->rows(), F->cols(), trans_F );
261  BLAS_Cpp::cols( F->rows(), F->cols(), trans_F ) != nd, std::invalid_argument
262  ,"QPSolverRelaxed::validate_input(...) : Error, op(F).cols() != d->dim()." );
263  TEUCHOS_TEST_FOR_EXCEPTION(
264  f->dim() != m_eq, std::invalid_argument
265  ,"QPSolverRelaxed::validate_input(...) : Error, f->dim() != op(F).rows()." );
266  TEUCHOS_TEST_FOR_EXCEPTION(
267  lambda->dim() != m_eq, std::invalid_argument
268  ,"QPSolverRelaxed::validate_input(...) : Error, lambda->dim() != op(F).rows()." );
269  TEUCHOS_TEST_FOR_EXCEPTION(
270  Fd && Fd->dim() != m_eq, std::invalid_argument
271  ,"QPSolverRelaxed::validate_input(...) : Error, Fd->dim() != op(F).rows()." );
272  }
273 }
274 
276  const value_type infinite_bound
277  ,std::ostream* out, EOutputLevel olevel
278  ,const Vector& g, const MatrixSymOp& G
279  ,value_type etaL
280  ,const Vector* dL, const Vector* dU
281  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
282  ,const Vector* eL, const Vector* eU
283  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
284  ,value_type* eta, VectorMutable* d
285  ,VectorMutable* nu
286  ,VectorMutable* mu
287  ,VectorMutable* lambda
288  )
289 {
291  if( out && (int)olevel >= (int)PRINT_ITER_STEPS ) {
292  *out<< "\n*** Printing input to QPSolverRelaxed::solve_qp(...) ...\n";
293  // g
294  *out << "\n||g||inf = " << g.norm_inf() << std::endl;
295  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
296  *out<< "g =\n" << g;
297  // G
298  if( (int)olevel >= (int)PRINT_EVERY_THING )
299  *out<< "\nG =\n" << G;
300  // etaL
301  *out << "\netaL = " << etaL << std::endl;
302  // eta
303  *out << "\neta = " << *eta << std::endl;
304  if(dL) {
305  // dL, dU
306  *out << "\ndL.dim() = " << dL->dim();
307  *out << "\ndU.dim() = " << dU->dim();
308  *out << "\nnum_bounded(dL,dU) = " << num_bounded(*dL,*dU,infinite_bound) << std::endl;
309  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
310  *out << "dL =\n" << *dL;
311  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
312  *out << "dU =\n" << *dU;
313  }
314  else {
315  *out << "\ndL = -inf";
316  *out << "\ndU = +inf";
317  }
318  // d
319  *out << "\n||d||inf = " << d->norm_inf() << std::endl;
320  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
321  *out<< "d =\n" << *d;
322  // nu
323  if(nu) {
324  *out<< "\nnu->nz() = " << nu->nz() << std::endl
325  << "||nu||inf = " << nu->norm_inf() << std::endl;
326  if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
327  *out<< "nu =\n" << *nu;
328  }
329  if(E) {
330  // op(E)
331  if( (int)olevel >= (int)PRINT_EVERY_THING )
332  *out<< "\nE" << std::endl << *E
333  << "trans_E = " << BLAS_Cpp::trans_to_string(trans_E) << std::endl;
334  // b
335  *out << "\n||b||inf = " << b->norm_inf() << std::endl;
336  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
337  *out<< "b =\n" << *b;
338  // eL, eU
339  *out<< "\neL.dim() = " << eL->dim();
340  *out<< "\neU.dim() = " << eU->dim();
341  *out << "\nnum_bounded(eL,eU) = " << num_bounded(*eL,*eU,infinite_bound) << std::endl;
342  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
343  *out<< "eL =\n" << *eL;
344  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
345  *out<< "eU =\n" << *eU;
346  // mu
347  *out<< "\nmu.nz() = " << mu->nz() << std::endl
348  << "||mu||inf = " << mu->norm_inf() << std::endl;
349  if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
350  *out<< "mu =\n" << *mu;
351  }
352  if(F) {
353  // op(F)
354  if( (int)olevel >= (int)PRINT_EVERY_THING )
355  *out<< "\nF" << std::endl << *F
356  << "trans_F = " << BLAS_Cpp::trans_to_string(trans_F) << std::endl;
357  // f
358  *out<< "\n||f||inf = " << f->norm_inf() << std::endl;
359  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
360  *out<< "f =\n" << *f;
361  // lambda
362  *out<< "\n||lambda||inf = " << lambda->norm_inf() << std::endl;
363  if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
364  *out<< "lambda =\n" << *lambda;
365  }
366  *out<< "\nEnd input to QPSolverRelaxed::solve_qp(...)\n";
367  }
368 }
369 
371  const value_type infinite_bound
372  ,std::ostream* out, EOutputLevel olevel
373  ,const value_type* obj_d
374  ,const value_type* eta, const Vector* d
375  ,const Vector* nu
376  ,const Vector* mu, const Vector* Ed
377  ,const Vector* lambda, const Vector* Fd
378  )
379 {
380  if( out && (int)olevel > (int)PRINT_ITER_STEPS ) {
381  *out<< "\n*** Printing output from QPSolverRelaxed::solve_qp(...) ...\n";
382  // obj_d
383  if(obj_d)
384  *out << "\nobj_d = " << *obj_d << std::endl;
385  // eta
386  *out << "\neta = " << *eta << std::endl;
387  // d
388  *out << "\n||d||inf = " << d->norm_inf() << std::endl;
389  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
390  *out<< "d =\n" << *d;
391  // nu
392  if(nu) {
393  *out<< "\nnu.nz() = " << nu->nz() << std::endl
394  << "||nu||inf = " << nu->norm_inf() << std::endl;
395  if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
396  *out<< "nu =\n" << *nu;
397  }
398  // Ed
399  if(Ed) {
400  *out << "\n||Ed||inf = " << Ed->norm_inf() << std::endl;
401  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
402  *out<< "Ed =\n" << *Ed;
403  }
404  // mu
405  if(mu) {
406  *out<< "\nmu.nz() = " << mu->nz() << std::endl
407  << "||mu||inf = " << mu->norm_inf() << std::endl;
408  if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
409  *out<< "mu =\n" << *mu;
410  }
411  // lambda
412  if(lambda) {
413  *out<< "\n||lambda||inf = " << lambda->norm_inf() << std::endl;
414  if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
415  *out<< "lambda =\n" << *lambda;
416  }
417  // Fd
418  if(Fd) {
419  *out << "\n||Fd||inf = " << Fd->norm_inf() << std::endl;
420  if( (int)olevel >= (int)PRINT_ITER_VECTORS )
421  *out<< "Fd =\n" << *Fd;
422  }
423  *out<< "\nEnd output from QPSolverRelaxed::solve_qp(...)\n";
424  }
425 }
426 
427 } // end namespace ConstrainedOptPack
ERunTests
Enumeration for if to run internal tests or not.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual QPSolverStats::ESolutionType imp_solve_qp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, 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, value_type *obj_d, value_type *eta, VectorMutable *d, VectorMutable *nu, VectorMutable *mu, VectorMutable *Ed, VectorMutable *lambda, VectorMutable *Fd)=0
Subclasses are to override this to implement the QP algorithm.
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
static void print_qp_output(const value_type infinite_bound, std::ostream *out, EOutputLevel olevel, 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)
Utility (static) function for printing the output input/output arguments after the QP solver is run...
static void print_qp_input(const value_type infinite_bound, std::ostream *out, EOutputLevel olevel, 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, value_type *eta, VectorMutable *d, VectorMutable *nu, VectorMutable *mu, VectorMutable *lambda)
Utility (static) function for printing the input input/output arguments before the QP solver is run...
size_t size_type
size_type num_bounded(const Vector &xl, const Vector &xu, value_type inf_bound)
virtual QPSolverStats::ESolutionType solve_qp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, 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, value_type *obj_d, value_type *eta, VectorMutable *d, VectorMutable *nu, VectorMutable *mu, VectorMutable *Ed, VectorMutable *lambda, VectorMutable *Fd)
Solve the QP.
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()...
const char * trans_to_string(Transp _trans)
EOutputLevel
Enumeration for the amount of output to create from solve_qp().
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
ESolutionType
Enumeration for the type of point returned from solve_qp(...).