external_packages/qpOASES-3.0beta/src/QProblemB.cpp
Go to the documentation of this file.
1 /*
2  * This file is part of qpOASES.
3  *
4  * qpOASES -- An Implementation of the Online Active Set Strategy.
5  * Copyright (C) 2007-2011 by Hans Joachim Ferreau, Andreas Potschka,
6  * Christian Kirches et al. All rights reserved.
7  *
8  * qpOASES is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2.1 of the License, or (at your option) any later version.
12  *
13  * qpOASES is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  * See the GNU Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with qpOASES; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  *
22  */
23 
24 
36 #include <qpOASES/QProblemB.hpp>
37 
38 BEGIN_NAMESPACE_QPOASES
39 
40 
41 /*****************************************************************************
42  * P U B L I C *
43  *****************************************************************************/
44 
45 
46 /*
47  * Q P r o b l e m B
48  */
49 QProblemB::QProblemB( )
50 {
51  /* print copyright notice */
52  if (options.printLevel > PL_NONE)
53  printCopyrightNotice( );
54 
55  /* reset global message handler */
56  getGlobalMessageHandler( )->reset( );
57 
58  freeHessian = BT_FALSE;
59  H = 0;
60 
61  g = 0;
62  lb = 0;
63  ub = 0;
64 
65  R = 0;
66  haveCholesky = BT_FALSE;
67 
68  x = 0;
69  y = 0;
70 
71  tau = 0.0;
72 
73  hessianType = HST_UNKNOWN;
74  isRegularised = BT_FALSE;
75 
76  infeasible = BT_FALSE;
77  unbounded = BT_FALSE;
78 
79  status = QPS_NOTINITIALISED;
80 
81  count = 0;
82 
83  ramp0 = options.initialRamping;
84  ramp1 = options.finalRamping;
85 
86  delta_xFR_TMP = 0;
87 
88  setPrintLevel( options.printLevel );
89 }
90 
91 
92 /*
93  * Q P r o b l e m B
94  */
95 QProblemB::QProblemB( int _nV, HessianType _hessianType )
96 {
97  /* print copyright notice */
98  if (options.printLevel > PL_NONE)
99  printCopyrightNotice( );
100 
101  /* consistency check */
102  if ( _nV <= 0 )
103  {
104  _nV = 1;
105  THROWERROR( RET_INVALID_ARGUMENTS );
106  }
107 
108  /* reset global message handler */
109  getGlobalMessageHandler( )->reset( );
110 
111  freeHessian = BT_FALSE;
112  H = 0;
113 
114  g = new real_t[_nV];
115  lb = new real_t[_nV];
116  ub = new real_t[_nV];
117 
118  bounds.init( _nV );
119 
120  R = new real_t[_nV*_nV];
121 
122  x = new real_t[_nV];
123  y = new real_t[_nV];
124 
125  tau = 0.0;
126 
127  hessianType = _hessianType;
128  isRegularised = BT_FALSE;
129 
130  infeasible = BT_FALSE;
131  unbounded = BT_FALSE;
132 
133  status = QPS_NOTINITIALISED;
134 
135  count = 0;
136 
137  ramp0 = options.initialRamping;
138  ramp1 = options.finalRamping;
139 
140  delta_xFR_TMP = new real_t[_nV];
141 
142  setPrintLevel( options.printLevel );
143 
144  flipper.init( _nV );
145 }
146 
147 
148 /*
149  * Q P r o b l e m B
150  */
151 QProblemB::QProblemB( const QProblemB& rhs )
152 {
153  freeHessian = BT_FALSE;
154  H = 0;
155 
156  copy( rhs );
157 }
158 
159 
160 /*
161  * ~ Q P r o b l e m B
162  */
163 QProblemB::~QProblemB( )
164 {
165  clear( );
166 
167  /* reset global message handler */
168  getGlobalMessageHandler( )->reset( );
169 }
170 
171 
172 /*
173  * o p e r a t o r =
174  */
175 QProblemB& QProblemB::operator=( const QProblemB& rhs )
176 {
177  if ( this != &rhs )
178  {
179  clear( );
180  copy( rhs );
181  }
182 
183  return *this;
184 }
185 
186 
187 /*
188  * r e s e t
189  */
190 returnValue QProblemB::reset( )
191 {
192  int i;
193  int nV = getNV( );
194 
195  if ( nV == 0 )
196  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
197 
198  /* 1) Reset bounds. */
199  bounds.init( nV );
200 
201  /* 2) Reset Cholesky decomposition. */
202  for( i=0; i<nV*nV; ++i )
203  R[i] = 0.0;
204 
205  haveCholesky = BT_FALSE;
206 
207  /* 3) Reset steplength and status flags. */
208  tau = 0.0;
209 
210  hessianType = HST_UNKNOWN;
211  isRegularised = BT_FALSE;
212 
213  infeasible = BT_FALSE;
214  unbounded = BT_FALSE;
215 
216  status = QPS_NOTINITIALISED;
217 
218  ramp0 = options.initialRamping;
219  ramp1 = options.finalRamping;
220 
221  return SUCCESSFUL_RETURN;
222 }
223 
224 
225 /*
226  * i n i t
227  */
228 returnValue QProblemB::init( SymmetricMatrix *_H, const real_t* const _g,
229  const real_t* const _lb, const real_t* const _ub,
230  int& nWSR, real_t* const cputime
231  )
232 {
233  if ( getNV( ) == 0 )
234  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
235 
236  /* 1) Consistency check. */
237  if ( isInitialised( ) == BT_TRUE )
238  {
239  THROWWARNING( RET_QP_ALREADY_INITIALISED );
240  reset( );
241  }
242 
243  /* 2) Setup QP data. */
244  if ( setupQPdata( _H,_g,_lb,_ub ) != SUCCESSFUL_RETURN )
245  return THROWERROR( RET_INVALID_ARGUMENTS );
246 
247  /* 3) Call to main initialisation routine (without any additional information). */
248  return solveInitialQP( 0,0,0, nWSR,cputime );
249 }
250 
251 
252 /*
253  * i n i t
254  */
255 returnValue QProblemB::init( const real_t* const _H, const real_t* const _g,
256  const real_t* const _lb, const real_t* const _ub,
257  int& nWSR, real_t* const cputime
258  )
259 {
260  if ( getNV( ) == 0 )
261  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
262 
263  /* 1) Consistency check. */
264  if ( isInitialised( ) == BT_TRUE )
265  {
266  THROWWARNING( RET_QP_ALREADY_INITIALISED );
267  reset( );
268  }
269 
270  /* 2) Setup QP data. */
271  if ( setupQPdata( _H,_g,_lb,_ub ) != SUCCESSFUL_RETURN )
272  return THROWERROR( RET_INVALID_ARGUMENTS );
273 
274  /* 3) Call to main initialisation routine (without any additional information). */
275  return solveInitialQP( 0,0,0, nWSR,cputime );
276 }
277 
278 
279 /*
280  * i n i t
281  */
282 returnValue QProblemB::init( const char* const H_file, const char* const g_file,
283  const char* const lb_file, const char* const ub_file,
284  int& nWSR, real_t* const cputime
285  )
286 {
287  if ( getNV( ) == 0 )
288  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
289 
290  /* 1) Consistency check. */
291  if ( isInitialised( ) == BT_TRUE )
292  {
293  THROWWARNING( RET_QP_ALREADY_INITIALISED );
294  reset( );
295  }
296 
297  /* 2) Setup QP data from files. */
298  if ( setupQPdataFromFile( H_file,g_file,lb_file,ub_file ) != SUCCESSFUL_RETURN )
299  return THROWERROR( RET_UNABLE_TO_READ_FILE );
300 
301  /* 3) Call to main initialisation routine (without any additional information). */
302  return solveInitialQP( 0,0,0, nWSR,cputime );
303 }
304 
305 
306 /*
307  * i n i t
308  */
309 returnValue QProblemB::init( SymmetricMatrix *_H, const real_t* const _g,
310  const real_t* const _lb, const real_t* const _ub,
311  int& nWSR, real_t* const cputime,
312  const real_t* const xOpt, const real_t* const yOpt,
313  const Bounds* const guessedBounds
314  )
315 {
316  int i;
317  int nV = getNV( );
318 
319  if ( nV == 0 )
320  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
321 
322  /* 1) Consistency checks. */
323  if ( isInitialised( ) == BT_TRUE )
324  {
325  THROWWARNING( RET_QP_ALREADY_INITIALISED );
326  reset( );
327  }
328 
329  if ( guessedBounds != 0 )
330  {
331  for( i=0; i<nV; ++i )
332  {
333  if ( guessedBounds->getStatus( i ) == ST_UNDEFINED )
334  return THROWERROR( RET_INVALID_ARGUMENTS );
335  }
336  }
337 
338  /* exclude this possibility in order to avoid inconsistencies */
339  if ( ( xOpt == 0 ) && ( yOpt != 0 ) && ( guessedBounds != 0 ) )
340  return THROWERROR( RET_INVALID_ARGUMENTS );
341 
342  /* 2) Setup QP data. */
343  if ( setupQPdata( _H,_g,_lb,_ub ) != SUCCESSFUL_RETURN )
344  return THROWERROR( RET_INVALID_ARGUMENTS );
345 
346  /* 3) Call to main initialisation routine. */
347  return solveInitialQP( xOpt,yOpt,guessedBounds, nWSR,cputime );
348 }
349 
350 
351 /*
352  * i n i t
353  */
354 returnValue QProblemB::init( const real_t* const _H, const real_t* const _g,
355  const real_t* const _lb, const real_t* const _ub,
356  int& nWSR, real_t* const cputime,
357  const real_t* const xOpt, const real_t* const yOpt,
358  const Bounds* const guessedBounds
359  )
360 {
361  int i;
362  int nV = getNV( );
363 
364  if ( nV == 0 )
365  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
366 
367  /* 1) Consistency checks. */
368  if ( isInitialised( ) == BT_TRUE )
369  {
370  THROWWARNING( RET_QP_ALREADY_INITIALISED );
371  reset( );
372  }
373 
374  if ( guessedBounds != 0 )
375  {
376  for( i=0; i<nV; ++i )
377  {
378  if ( guessedBounds->getStatus( i ) == ST_UNDEFINED )
379  return THROWERROR( RET_INVALID_ARGUMENTS );
380  }
381  }
382 
383  /* exclude this possibility in order to avoid inconsistencies */
384  if ( ( xOpt == 0 ) && ( yOpt != 0 ) && ( guessedBounds != 0 ) )
385  return THROWERROR( RET_INVALID_ARGUMENTS );
386 
387  /* 2) Setup QP data. */
388  if ( setupQPdata( _H,_g,_lb,_ub ) != SUCCESSFUL_RETURN )
389  return THROWERROR( RET_INVALID_ARGUMENTS );
390 
391  /* 3) Call to main initialisation routine. */
392  return solveInitialQP( xOpt,yOpt,guessedBounds, nWSR,cputime );
393 }
394 
395 
396 /*
397  * i n i t
398  */
399 returnValue QProblemB::init( const char* const H_file, const char* const g_file,
400  const char* const lb_file, const char* const ub_file,
401  int& nWSR, real_t* const cputime,
402  const real_t* const xOpt, const real_t* const yOpt,
403  const Bounds* const guessedBounds
404  )
405 {
406  int i;
407  int nV = getNV( );
408 
409  if ( nV == 0 )
410  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
411 
412  /* 1) Consistency checks. */
413  if ( isInitialised( ) == BT_TRUE )
414  {
415  THROWWARNING( RET_QP_ALREADY_INITIALISED );
416  reset( );
417  }
418 
419  for( i=0; i<nV; ++i )
420  {
421  if ( guessedBounds->getStatus( i ) == ST_UNDEFINED )
422  return THROWERROR( RET_INVALID_ARGUMENTS );
423  }
424 
425  /* exclude this possibility in order to avoid inconsistencies */
426  if ( ( xOpt == 0 ) && ( yOpt != 0 ) && ( guessedBounds != 0 ) )
427  return THROWERROR( RET_INVALID_ARGUMENTS );
428 
429  /* 2) Setup QP data from files. */
430  if ( setupQPdataFromFile( H_file,g_file,lb_file,ub_file ) != SUCCESSFUL_RETURN )
431  return THROWERROR( RET_UNABLE_TO_READ_FILE );
432 
433  /* 3) Call to main initialisation routine. */
434  return solveInitialQP( xOpt,yOpt,guessedBounds, nWSR,cputime );
435 }
436 
437 
438 /*
439  * h o t s t a r t
440  */
441 returnValue QProblemB::hotstart( const real_t* const g_new,
442  const real_t* const lb_new, const real_t* const ub_new,
443  int& nWSR, real_t* const cputime
444  )
445 {
446  if ( getNV( ) == 0 )
447  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
448 
449  ++count;
450 
451  returnValue returnvalue = SUCCESSFUL_RETURN;
452  int nFar, i;
453  int nV = getNV ();
454 
455  int nWSR_max = nWSR;
456  int nWSR_performed = 0;
457 
458  real_t cputime_remaining = INFTY;
459  real_t cputime_needed = 0.0;
460 
461  real_t farbound = options.initialFarBounds;
462  real_t t, rampVal;
463 
464  if ( options.enableFarBounds == BT_FALSE )
465  {
466  /*if (haveCholesky == BT_FALSE)
467  {
468  returnvalue = computeInitialCholesky();
469  if (returnvalue != SUCCESSFUL_RETURN)
470  return THROWERROR(returnvalue);
471  }*/
472 
473  if ( usingRegularisation( ) == BT_TRUE )
474  returnvalue = solveRegularisedQP( g_new,lb_new,ub_new, nWSR,cputime,0 );
475  else
476  returnvalue = solveQP( g_new,lb_new,ub_new, nWSR,cputime,0 );
477  }
478  else
479  {
480  real_t *ub_new_far = new real_t[nV];
481  real_t *lb_new_far = new real_t[nV];
482 
483  /* possibly extend initial far bounds to largest bound/constraint data */
484  if (ub_new)
485  for (i = 0; i < nV; i++)
486  if (ub_new[i] < INFTY && ub_new[i] > farbound) farbound = ub_new[i];
487  if (lb_new)
488  for (i = 0; i < nV; i++)
489  if (lb_new[i] > -INFTY && lb_new[i] < -farbound) farbound = -lb_new[i];
490 
491  if ( options.enableRamping == BT_TRUE )
492  {
493  /* TODO: encapsule this part to avoid code duplication! */
494  for ( i=0; i<nV; ++i )
495  {
496  t = static_cast<real_t>(i) / static_cast<real_t>(nV-1);
497  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
498 
499  if ( ( lb_new == 0 ) || ( lb_new[i] <= -rampVal ) )
500  lb_new_far[i] = - rampVal;
501  else
502  lb_new_far[i] = lb_new[i];
503  if ( ( ub_new == 0 ) || ( ub_new[i] >= rampVal ) )
504  ub_new_far[i] = rampVal;
505  else
506  ub_new_far[i] = ub_new[i];
507  }
508  }
509  else
510  {
511  for ( i=0; i<nV; ++i )
512  {
513  lb_new_far[i] = lb_new[i];
514  ub_new_far[i] = ub_new[i];
515  }
516  }
517 
518  /*if (haveCholesky == BT_FALSE)
519  {
520  returnvalue = computeInitialCholesky();
521  if (returnvalue != SUCCESSFUL_RETURN)
522  goto farewell;
523  }*/
524 
525  for ( ;; )
526  {
527  nWSR = nWSR_max;
528  if ( cputime != 0 )
529  cputime_remaining = *cputime - cputime_needed;
530 
531  if ( usingRegularisation( ) == BT_TRUE )
532  returnvalue = solveRegularisedQP( g_new,lb_new_far,ub_new_far, nWSR,&cputime_remaining,nWSR_performed );
533  else
534  returnvalue = solveQP( g_new,lb_new_far,ub_new_far, nWSR,&cputime_remaining,nWSR_performed );
535 
536  nWSR_performed = nWSR;
537  cputime_needed += cputime_remaining;
538 
539  /* Check for active far-bounds and move them away */
540  nFar = 0;
541  farbound *= options.growFarBounds;
542 
543  real_t maxFarbound = 1e20;
544  if ( infeasible == BT_TRUE )
545  {
546  if ( farbound > maxFarbound )
547  {
548  returnvalue = RET_HOTSTART_STOPPED_INFEASIBILITY;
549  goto farewell;
550  }
551 
552  if ( options.enableRamping == BT_TRUE )
553  {
554  /* TODO: encapsule this part to avoid code duplication! */
555  for ( i=0; i<nV; ++i )
556  {
557  t = static_cast<real_t>(i) / static_cast<real_t>(nV-1);
558  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
559 
560  if ( lb_new == 0 )
561  lb_new_far[i] = - rampVal;
562  else
563  lb_new_far[i] = getMax (- rampVal, lb_new[i]);
564 
565  if ( ub_new == 0 )
566  ub_new_far[i] = rampVal;
567  else
568  ub_new_far[i] = getMin (rampVal, ub_new[i]);
569  }
570  }
571  else
572  {
573  for ( i=0; i<nV; ++i )
574  {
575  lb_new_far[i] = lb_new[i];
576  ub_new_far[i] = ub_new[i];
577  }
578  }
579  }
580  else if ( status == QPS_SOLVED )
581  {
582  real_t tol = farbound * options.boundTolerance;
583  status = QPS_HOMOTOPYQPSOLVED;
584 
585  nFar = 0;
586  for ( i=0; i<nV; ++i )
587  {
588  if ( ( ( lb_new == 0 ) || ( lb_new_far[i] > lb_new[i] ) ) && ( fabs ( lb_new_far[i] - x[i] ) < tol ) )
589  ++nFar;
590  if ( ( ( ub_new == 0 ) || ( ub_new_far[i] < ub_new[i] ) ) && ( fabs ( ub_new_far[i] - x[i] ) < tol ) )
591  ++nFar;
592  }
593 
594  if ( nFar == 0 )
595  break;
596 
597  if ( farbound > maxFarbound )
598  {
599  unbounded = BT_TRUE;
600  returnvalue = RET_HOTSTART_STOPPED_UNBOUNDEDNESS;
601  goto farewell;
602  }
603 
604  if ( options.enableRamping == BT_TRUE )
605  {
606  /* TODO: encapsule this part to avoid code duplication! */
607  for ( i=0; i<nV; ++i )
608  {
609  t = static_cast<real_t>(i) / static_cast<real_t>(nV-1);
610  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
611 
612  if ( lb_new == 0 )
613  lb_new_far[i] = - rampVal;
614  else
615  lb_new_far[i] = getMax (- rampVal, lb_new[i]);
616 
617  if ( ub_new == 0 )
618  ub_new_far[i] = rampVal;
619  else
620  ub_new_far[i] = getMin (rampVal, ub_new[i]);
621  }
622  }
623  else
624  {
625  for ( i=0; i<nV; ++i )
626  {
627  lb_new_far[i] = lb_new[i];
628  ub_new_far[i] = ub_new[i];
629  }
630  }
631  }
632  else
633  {
634  /* some other error */
635  break;
636  }
637 
638  /* swap ramp0 and ramp1 */
639  t = ramp0; ramp0 = ramp1; ramp1 = t;
640  }
641 
642  farewell:
643  if ( cputime != 0 )
644  *cputime = cputime_needed;
645  delete[] lb_new_far; delete[] ub_new_far;
646  }
647 
648  return ( returnvalue != SUCCESSFUL_RETURN ) ? THROWERROR( returnvalue ) : returnvalue;
649 }
650 
651 
652 /*
653  * h o t s t a r t
654  */
655 returnValue QProblemB::hotstart( const char* const g_file,
656  const char* const lb_file, const char* const ub_file,
657  int& nWSR, real_t* const cputime
658  )
659 {
660  int nV = getNV( );
661 
662  if ( nV == 0 )
663  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
664 
665  /* consistency check */
666  if ( g_file == 0 )
667  return THROWERROR( RET_INVALID_ARGUMENTS );
668 
669 
670  /* 1) Allocate memory (if bounds exist). */
671  real_t* g_new = new real_t[nV];
672  real_t* lb_new = 0;
673  real_t* ub_new = 0;
674 
675  if ( lb_file != 0 )
676  lb_new = new real_t[nV];
677  if ( ub_file != 0 )
678  ub_new = new real_t[nV];
679 
680  /* 2) Load new QP vectors from file. */
681  returnValue returnvalue;
682  returnvalue = loadQPvectorsFromFile( g_file,lb_file,ub_file,
683  g_new,lb_new,ub_new
684  );
685  if ( returnvalue != SUCCESSFUL_RETURN )
686  {
687  if ( ub_file != 0 )
688  delete[] ub_new;
689  if ( lb_file != 0 )
690  delete[] lb_new;
691  delete[] g_new;
692 
693  return THROWERROR( RET_UNABLE_TO_READ_FILE );
694  }
695 
696  /* 3) Actually perform hotstart. */
697  returnvalue = hotstart( g_new,lb_new,ub_new, nWSR,cputime );
698 
699  /* 4) Free memory. */
700  if ( ub_file != 0 )
701  delete[] ub_new;
702  if ( lb_file != 0 )
703  delete[] lb_new;
704  delete[] g_new;
705 
706  return returnvalue;
707 }
708 
709 
710 /*
711  * h o t s t a r t
712  */
713 returnValue QProblemB::hotstart( const real_t* const g_new,
714  const real_t* const lb_new, const real_t* const ub_new,
715  int& nWSR, real_t* const cputime,
716  const Bounds* const guessedBounds
717  )
718 {
719  int nV = getNV( );
720 
721  if ( nV == 0 )
722  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
723 
724 
725  /* start runtime measurement */
726  real_t starttime = 0.0;
727  if ( cputime != 0 )
728  starttime = getCPUtime( );
729 
730 
731  /* 1) Update working set according to guess for working set of bounds. */
732  if ( guessedBounds != 0 )
733  {
734  if ( setupAuxiliaryQP( guessedBounds ) != SUCCESSFUL_RETURN )
735  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
736  }
737  else
738  {
739  /* create empty bounds for setting up auxiliary QP */
740  Bounds emptyBounds( nV );
741  if ( emptyBounds.setupAllFree( ) != SUCCESSFUL_RETURN )
742  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
743 
744  if ( setupAuxiliaryQP( &emptyBounds ) != SUCCESSFUL_RETURN )
745  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
746  }
747 
748  /* 2) Perform usual homotopy. */
749 
750  /* Allow only remaining CPU time for usual hotstart. */
751  if ( cputime != 0 )
752  *cputime -= getCPUtime( ) - starttime;
753 
754  returnValue returnvalue = hotstart( g_new,lb_new,ub_new, nWSR,cputime );
755 
756  /* stop runtime measurement */
757  if ( cputime != 0 )
758  *cputime = getCPUtime( ) - starttime;
759 
760  return returnvalue;
761 }
762 
763 
764 /*
765  * h o t s t a r t
766  */
767 returnValue QProblemB::hotstart( const char* const g_file,
768  const char* const lb_file, const char* const ub_file,
769  int& nWSR, real_t* const cputime,
770  const Bounds* const guessedBounds
771  )
772 {
773  int nV = getNV( );
774 
775  if ( nV == 0 )
776  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
777 
778  /* consistency check */
779  if ( g_file == 0 )
780  return THROWERROR( RET_INVALID_ARGUMENTS );
781 
782 
783  /* 1) Allocate memory (if bounds exist). */
784  real_t* g_new = new real_t[nV];
785  real_t* lb_new = 0;
786  real_t* ub_new = 0;
787 
788  if ( lb_file != 0 )
789  lb_new = new real_t[nV];
790  if ( ub_file != 0 )
791  ub_new = new real_t[nV];
792 
793  /* 2) Load new QP vectors from file. */
794  returnValue returnvalue;
795  returnvalue = loadQPvectorsFromFile( g_file,lb_file,ub_file,
796  g_new,lb_new,ub_new
797  );
798  if ( returnvalue != SUCCESSFUL_RETURN )
799  {
800  if ( ub_file != 0 )
801  delete[] ub_new;
802  if ( lb_file != 0 )
803  delete[] lb_new;
804  delete[] g_new;
805 
806  return THROWERROR( RET_UNABLE_TO_READ_FILE );
807  }
808 
809  /* 3) Actually perform hotstart using initialised homotopy. */
810  returnvalue = hotstart( g_new,lb_new,ub_new, nWSR,cputime,
811  guessedBounds
812  );
813 
814  /* 4) Free memory. */
815  if ( ub_file != 0 )
816  delete[] ub_new;
817  if ( lb_file != 0 )
818  delete[] lb_new;
819  delete[] g_new;
820 
821  return returnvalue;
822 }
823 
824 
825 /*
826  * g e t N Z
827  */
828 int QProblemB::getNZ( ) const
829 {
830  /* if no constraints are present: nZ=nFR */
831  return getNFR( );
832 }
833 
834 
835 /*
836  * g e t O b j V a l
837  */
838 real_t QProblemB::getObjVal( ) const
839 {
840  real_t objVal;
841 
842  /* calculated optimal objective function value
843  * only if current QP has been solved */
844  if ( ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
845  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
846  ( getStatus( ) == QPS_SOLVED ) )
847  {
848  objVal = getObjVal( x );
849  }
850  else
851  {
852  objVal = INFTY;
853  }
854 
855  return objVal;
856 }
857 
858 
859 /*
860  * g e t O b j V a l
861  */
862 real_t QProblemB::getObjVal( const real_t* const _x ) const
863 {
864  int i;
865  int nV = getNV( );
866 
867  if ( nV == 0 )
868  return 0.0;
869 
870  real_t objVal = 0.0;
871 
872  for( i=0; i<nV; ++i )
873  objVal += _x[i]*g[i];
874 
875  switch ( hessianType )
876  {
877  case HST_ZERO:
878  break;
879 
880  case HST_IDENTITY:
881  for( i=0; i<nV; ++i )
882  objVal += 0.5*_x[i]*_x[i];
883  break;
884 
885  default:
886  real_t *Hx = new real_t[nV];
887  H->times(1, 1.0, _x, nV, 0.0, Hx, nV);
888  for( i=0; i<nV; ++i )
889  objVal += 0.5*_x[i]*Hx[i];
890  delete[] Hx;
891  break;
892  }
893 
894  /* When using regularisation, the objective function value
895  * needs to be modified as follows:
896  * objVal = objVal - 0.5*_x*(Hmod-H)*_x - _x'*(gMod-g)
897  * = objVal - 0.5*_x*eps*_x * - _x'*(-eps*_x)
898  * = objVal + 0.5*_x*eps*_x */
899  if ( usingRegularisation( ) == BT_TRUE )
900  {
901  for( i=0; i<nV; ++i )
902  objVal += 0.5*_x[i]*options.epsRegularisation*_x[i];
903  }
904 
905  return objVal;
906 }
907 
908 
909 /*
910  * g e t P r i m a l S o l u t i o n
911  */
912 returnValue QProblemB::getPrimalSolution( real_t* const xOpt ) const
913 {
914  int i;
915 
916  /* return optimal primal solution vector
917  * only if current QP has been solved */
918  if ( ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
919  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
920  ( getStatus( ) == QPS_SOLVED ) )
921  {
922  for( i=0; i<getNV( ); ++i )
923  xOpt[i] = x[i];
924 
925  return SUCCESSFUL_RETURN;
926  }
927  else
928  {
929  return RET_QP_NOT_SOLVED;
930  }
931 }
932 
933 
934 /*
935  * g e t D u a l S o l u t i o n
936  */
937 returnValue QProblemB::getDualSolution( real_t* const yOpt ) const
938 {
939  int i;
940 
941  /* return optimal dual solution vector
942  * only if current QP has been solved */
943  if ( ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
944  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
945  ( getStatus( ) == QPS_SOLVED ) )
946  {
947  for( i=0; i<getNV( ); ++i )
948  yOpt[i] = y[i];
949 
950  return SUCCESSFUL_RETURN;
951  }
952  else
953  {
954  return RET_QP_NOT_SOLVED;
955  }
956 }
957 
958 
959 /*
960  * s e t P r i n t L e v e l
961  */
962 returnValue QProblemB::setPrintLevel( PrintLevel _printLevel )
963 {
964  #ifndef __MATLAB__
965  if ( ( options.printLevel == PL_HIGH ) && ( options.printLevel != _printLevel ) )
966  THROWINFO( RET_PRINTLEVEL_CHANGED );
967  #endif
968 
969  options.printLevel = _printLevel;
970 
971  /* update message handler preferences */
972  switch ( options.printLevel )
973  {
974  case PL_TABULAR:
975  case PL_NONE:
976  getGlobalMessageHandler( )->setErrorVisibilityStatus( VS_HIDDEN );
977  getGlobalMessageHandler( )->setWarningVisibilityStatus( VS_HIDDEN );
978  getGlobalMessageHandler( )->setInfoVisibilityStatus( VS_HIDDEN );
979  break;
980 
981  case PL_LOW:
982  #ifndef __XPCTARGET__
983  getGlobalMessageHandler( )->setErrorVisibilityStatus( VS_VISIBLE );
984  getGlobalMessageHandler( )->setWarningVisibilityStatus( VS_HIDDEN );
985  getGlobalMessageHandler( )->setInfoVisibilityStatus( VS_HIDDEN );
986  #endif
987  break;
988 
989  case PL_MEDIUM:
990  #ifndef __XPCTARGET__
991  getGlobalMessageHandler( )->setErrorVisibilityStatus( VS_VISIBLE );
992  getGlobalMessageHandler( )->setWarningVisibilityStatus( VS_VISIBLE );
993  getGlobalMessageHandler( )->setInfoVisibilityStatus( VS_HIDDEN );
994  #endif
995  break;
996 
997  default: /* PL_HIGH */
998  #ifndef __XPCTARGET__
999  getGlobalMessageHandler( )->setErrorVisibilityStatus( VS_VISIBLE );
1000  getGlobalMessageHandler( )->setWarningVisibilityStatus( VS_VISIBLE );
1001  getGlobalMessageHandler( )->setInfoVisibilityStatus( VS_VISIBLE );
1002  #endif
1003  break;
1004  }
1005 
1006  return SUCCESSFUL_RETURN;
1007 }
1008 
1009 
1010 /*
1011  * p r i n t P r o p e r t i e s
1012  */
1013 returnValue QProblemB::printProperties( )
1014 {
1015  #ifndef __XPCTARGET__
1016  #ifndef __DSPACE__
1017  /* Do not print properties if print level is set to none! */
1018  if ( options.printLevel == PL_NONE )
1019  return SUCCESSFUL_RETURN;
1020 
1021  char myPrintfString[80];
1022 
1023  myPrintf( "\n################# qpOASES -- QP PROPERTIES #################\n" );
1024  myPrintf( "\n" );
1025 
1026  /* 1) Variables properties. */
1027  snprintf( myPrintfString,80, "Number of Variables: %4.1d\n",getNV( ) );
1028  myPrintf( myPrintfString );
1029 
1030  if ( bounds.hasNoLower( ) == BT_TRUE )
1031  myPrintf( "Variables are not bounded from below.\n" );
1032  else
1033  myPrintf( "Variables are bounded from below.\n" );
1034 
1035  if ( bounds.hasNoUpper( ) == BT_TRUE )
1036  myPrintf( "Variables are not bounded from above.\n" );
1037  else
1038  myPrintf( "Variables are bounded from above.\n" );
1039 
1040  myPrintf( "\n" );
1041 
1042 
1043  /* 2) Further properties. */
1044  switch ( hessianType )
1045  {
1046  case HST_ZERO:
1047  myPrintf( "Hessian is zero matrix (i.e. actually an LP is solved).\n" );
1048  break;
1049 
1050  case HST_IDENTITY:
1051  myPrintf( "Hessian is identity matrix.\n" );
1052  break;
1053 
1054  case HST_POSDEF:
1055  myPrintf( "Hessian matrix is (strictly) positive definite.\n" );
1056  break;
1057 
1058  case HST_POSDEF_NULLSPACE:
1059  myPrintf( "Hessian matrix is positive definite on null space of active constraints.\n" );
1060  break;
1061 
1062  case HST_SEMIDEF:
1063  myPrintf( "Hessian matrix is positive semi-definite.\n" );
1064  break;
1065 
1066  default:
1067  myPrintf( "Hessian matrix has unknown type.\n" );
1068  break;
1069  }
1070 
1071  if ( infeasible == BT_TRUE )
1072  myPrintf( "QP was found to be infeasible.\n" );
1073  else
1074  myPrintf( "QP seems to be feasible.\n" );
1075 
1076  if ( unbounded == BT_TRUE )
1077  myPrintf( "QP was found to be unbounded from below.\n" );
1078  else
1079  myPrintf( "QP seems to be bounded from below.\n" );
1080 
1081  myPrintf( "\n" );
1082 
1083 
1084  /* 3) QP object properties. */
1085  switch ( status )
1086  {
1087  case QPS_NOTINITIALISED:
1088  myPrintf( "Status of QP object: freshly instantiated or reset.\n" );
1089  break;
1090 
1091  case QPS_PREPARINGAUXILIARYQP:
1092  myPrintf( "Status of QP object: an auxiliary QP is currently setup.\n" );
1093  break;
1094 
1095  case QPS_AUXILIARYQPSOLVED:
1096  myPrintf( "Status of QP object: an auxilary QP was solved.\n" );
1097  break;
1098 
1099  case QPS_PERFORMINGHOMOTOPY:
1100  myPrintf( "Status of QP object: a homotopy step is performed.\n" );
1101  break;
1102 
1103  case QPS_HOMOTOPYQPSOLVED:
1104  myPrintf( "Status of QP object: an intermediate QP along the homotopy path was solved.\n" );
1105  break;
1106 
1107  case QPS_SOLVED:
1108  myPrintf( "Status of QP object: solution of the actual QP was found.\n" );
1109  break;
1110  }
1111 
1112  switch ( options.printLevel )
1113  {
1114  case PL_LOW:
1115  myPrintf( "Print level of QP object is low, i.e. only error are printed.\n" );
1116  break;
1117 
1118  case PL_MEDIUM:
1119  myPrintf( "Print level of QP object is medium, i.e. error and warnings are printed.\n" );
1120  break;
1121 
1122  case PL_HIGH:
1123  myPrintf( "Print level of QP object is high, i.e. all available output is printed.\n" );
1124  break;
1125 
1126  default:
1127  break;
1128  }
1129 
1130  myPrintf( "\n" );
1131  #endif
1132  #endif
1133 
1134  return SUCCESSFUL_RETURN;
1135 }
1136 
1137 
1138 returnValue QProblemB::printOptions( ) const
1139 {
1140  return options.print( );
1141 }
1142 
1143 
1144 
1145 /*****************************************************************************
1146  * P R O T E C T E D *
1147  *****************************************************************************/
1148 
1149 /*
1150  * c l e a r
1151  */
1152 returnValue QProblemB::clear( )
1153 {
1154  if ( ( freeHessian == BT_TRUE ) && ( H != 0 ) )
1155  {
1156  delete H;
1157  H = 0;
1158  }
1159 
1160  if ( g != 0 )
1161  {
1162  delete[] g;
1163  g = 0;
1164  }
1165 
1166  if ( lb != 0 )
1167  {
1168  delete[] lb;
1169  lb = 0;
1170  }
1171 
1172  if ( ub != 0 )
1173  {
1174  delete[] ub;
1175  ub = 0;
1176  }
1177 
1178  if ( R != 0 )
1179  {
1180  delete[] R;
1181  R = 0;
1182  }
1183 
1184  if ( x != 0 )
1185  {
1186  delete[] x;
1187  x = 0;
1188  }
1189 
1190  if ( y != 0 )
1191  {
1192  delete[] y;
1193  y = 0;
1194  }
1195 
1196  if ( delta_xFR_TMP != 0 )
1197  {
1198  delete[] delta_xFR_TMP;
1199  delta_xFR_TMP = 0;
1200  }
1201 
1202  return SUCCESSFUL_RETURN;
1203 }
1204 
1205 
1206 /*
1207  * c o p y
1208  */
1209 returnValue QProblemB::copy( const QProblemB& rhs
1210  )
1211 {
1212  int _nV = rhs.getNV( );
1213 
1214  bounds = rhs.bounds;
1215 
1216  freeHessian = rhs.freeHessian;
1217 
1218  if ( freeHessian == BT_TRUE )
1219  H = dynamic_cast<SymmetricMatrix *>(rhs.H->duplicate());
1220  else
1221  H = rhs.H;
1222 
1223  if ( rhs.g != 0 )
1224  {
1225  g = new real_t[_nV];
1226  setG( rhs.g );
1227  }
1228  else
1229  g = 0;
1230 
1231  if ( rhs.lb != 0 )
1232  {
1233  lb = new real_t[_nV];
1234  setLB( rhs.lb );
1235  }
1236  else
1237  lb = 0;
1238 
1239  if ( rhs.ub != 0 )
1240  {
1241  ub = new real_t[_nV];
1242  setUB( rhs.ub );
1243  }
1244  else
1245  ub = 0;
1246 
1247  if ( rhs.R != 0 )
1248  {
1249  R = new real_t[_nV*_nV];
1250  memcpy( R,rhs.R,_nV*_nV*sizeof(real_t) );
1251  }
1252  else
1253  R = 0;
1254 
1255  haveCholesky = rhs.haveCholesky;
1256 
1257  if ( rhs.x != 0 )
1258  {
1259  x = new real_t[_nV];
1260  memcpy( x,rhs.x,_nV*sizeof(real_t) );
1261  }
1262  else
1263  x = 0;
1264 
1265  if ( rhs.y != 0 )
1266  {
1267  y = new real_t[_nV];
1268  memcpy( y,rhs.y,_nV*sizeof(real_t) );
1269  }
1270  else
1271  y = 0;
1272 
1273  tau = rhs.tau;
1274 
1275  hessianType = rhs.hessianType;
1276  isRegularised = rhs.isRegularised;
1277 
1278  infeasible = rhs.infeasible;
1279  unbounded = rhs.unbounded;
1280 
1281  status = rhs.status;
1282 
1283  count = rhs.count;
1284 
1285  ramp0 = rhs.ramp0;
1286  ramp1 = rhs.ramp1;
1287 
1288  delta_xFR_TMP = new real_t[_nV]; /* nFR */
1289 
1290  options = rhs.options;
1291  setPrintLevel( options.printLevel );
1292 
1293  flipper = rhs.flipper;
1294 
1295  return SUCCESSFUL_RETURN;
1296 }
1297 
1298 
1299 /*
1300  * d e t e r m i n e H e s s i a n T y p e
1301  */
1302 returnValue QProblemB::determineHessianType( )
1303 {
1304  int i;
1305  int nV = getNV( );
1306 
1307  /* if Hessian type has been set by user, do NOT change it! */
1308  if ( hessianType != HST_UNKNOWN )
1309  return SUCCESSFUL_RETURN;
1310 
1311  /* if Hessian has not been allocated, assume it to be all zeros! */
1312  if ( H == 0 )
1313  {
1314  hessianType = HST_ZERO;
1315  return SUCCESSFUL_RETURN;
1316  }
1317 
1318 
1319  /* 1) If Hessian has outer-diagonal elements,
1320  * Hessian is assumed to be positive definite. */
1321  hessianType = HST_POSDEF;
1322  if (H->isDiag() == BT_FALSE)
1323  return SUCCESSFUL_RETURN;
1324 
1325  /* 2) Otherwise it is diagonal and test for identity or zero matrix is performed. */
1326  /* hessianType = HST_DIAGONAL; */
1327 
1328  BooleanType isIdentity = BT_TRUE;
1329  BooleanType isZero = BT_TRUE;
1330 
1331  for ( i=0; i<nV; ++i )
1332  {
1333  if ( H->diag(i) < -ZERO ) /* TODO: Should work with flipping bounds enabled! No error message needed. */
1334  return THROWERROR( RET_HESSIAN_INDEFINITE );
1335 
1336  if ( fabs( H->diag(i) - 1.0 ) > EPS )
1337  isIdentity = BT_FALSE;
1338 
1339  if ( fabs( H->diag(i) ) > EPS )
1340  isZero = BT_FALSE;
1341  }
1342 
1343  if ( isIdentity == BT_TRUE )
1344  hessianType = HST_IDENTITY;
1345 
1346  if ( isZero == BT_TRUE )
1347  hessianType = HST_ZERO;
1348 
1349  return SUCCESSFUL_RETURN;
1350 }
1351 
1352 
1353 /*
1354  * s e t u p S u b j e c t T o T y p e
1355  */
1356 returnValue QProblemB::setupSubjectToType( )
1357 {
1358  return setupSubjectToType( lb,ub );
1359 }
1360 
1361 
1362 /*
1363  * s e t u p S u b j e c t T o T y p e
1364  */
1365 returnValue QProblemB::setupSubjectToType( const real_t* const lb_new, const real_t* const ub_new )
1366 {
1367  int i;
1368  int nV = getNV( );
1369 
1370 
1371  /* 1) Check if lower bounds are present. */
1372  bounds.setNoLower( BT_TRUE );
1373  if ( lb_new != 0 )
1374  {
1375  for( i=0; i<nV; ++i )
1376  {
1377  if ( lb_new[i] > -INFTY )
1378  {
1379  bounds.setNoLower( BT_FALSE );
1380  break;
1381  }
1382  }
1383  }
1384 
1385  /* 2) Check if upper bounds are present. */
1386  bounds.setNoUpper( BT_TRUE );
1387  if ( ub_new != 0 )
1388  {
1389  for( i=0; i<nV; ++i )
1390  {
1391  if ( ub_new[i] < INFTY )
1392  {
1393  bounds.setNoUpper( BT_FALSE );
1394  break;
1395  }
1396  }
1397  }
1398 
1399  /* 3) Determine implicitly fixed and unbounded variables. */
1400  if ( ( lb_new != 0 ) && ( ub_new != 0 ) )
1401  {
1402  for( i=0; i<nV; ++i )
1403  {
1404  if ( ( lb_new[i] < -INFTY + options.boundTolerance ) && ( ub_new[i] > INFTY - options.boundTolerance )
1405  && (options.enableFarBounds == BT_FALSE))
1406  {
1407  bounds.setType( i,ST_UNBOUNDED );
1408  }
1409  else
1410  {
1411  if ( options.enableEqualities && lb_new[i] > ub_new[i] - options.boundTolerance )
1412  bounds.setType( i,ST_EQUALITY );
1413  else
1414  bounds.setType( i,ST_BOUNDED );
1415  }
1416  }
1417  }
1418  else
1419  {
1420  if ( ( lb_new == 0 ) && ( ub_new == 0 ) )
1421  {
1422  for( i=0; i<nV; ++i )
1423  bounds.setType( i,ST_UNBOUNDED );
1424  }
1425  else
1426  {
1427  for( i=0; i<nV; ++i )
1428  bounds.setType( i,ST_BOUNDED );
1429  }
1430  }
1431 
1432  return SUCCESSFUL_RETURN;
1433 }
1434 
1435 
1436 /*
1437  * s e t u p C h o l e s k y D e c o m p o s i t i o n
1438  */
1439 returnValue QProblemB::setupCholeskyDecomposition( )
1440 {
1441  int i, j;
1442  int nV = getNV( );
1443  int nFR = getNFR( );
1444 
1445  /* 1) Initialises R with all zeros. */
1446  for( i=0; i<nV*nV; ++i )
1447  R[i] = 0.0;
1448 
1449  /* 2) Calculate Cholesky decomposition of H (projected to free variables). */
1450  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_IDENTITY ) )
1451  {
1452  if ( hessianType == HST_ZERO )
1453  {
1454  /* if Hessian is zero matrix, it is assumed that it has been
1455  * regularised and thus its Cholesky factor is the identity
1456  * matrix scaled by sqrt(eps). */
1457  if ( usingRegularisation( ) == BT_TRUE )
1458  {
1459  for( i=0; i<nV; ++i )
1460  RR(i,i) = sqrt( options.epsRegularisation );
1461  }
1462  else
1463  return THROWERROR( RET_CHOLESKY_OF_ZERO_HESSIAN );
1464  }
1465  else
1466  {
1467  /* if Hessian is identity, so is its Cholesky factor. */
1468  for( i=0; i<nV; ++i )
1469  RR(i,i) = 1.0;
1470  }
1471  }
1472  else
1473  {
1474  if ( nFR > 0 )
1475  {
1476  int* FR_idx;
1477  bounds.getFree( )->getNumberArray( &FR_idx );
1478 
1479  /* get H */
1480  for ( j=0; j < nFR; ++j )
1481  H->getCol (FR_idx[j], bounds.getFree (), 1.0, &R[j*nV]);
1482 
1483  /* R'*R = H */
1484  long info = 0;
1485  unsigned long _nFR = nFR, _nV = nV;
1486 
1487  POTRF ("U", &_nFR, R, &_nV, &info);
1488 
1489  /* <0 = invalid call, =0 ok, >0 not spd */
1490  if (info > 0) {
1491  hessianType = HST_SEMIDEF;
1492  return RET_HESSIAN_NOT_SPD;
1493  }
1494 
1495  /* zero first subdiagonal to make givens updates work */
1496  for (i=0;i<nFR-1;++i)
1497  RR(i+1,i) = 0.0;
1498 
1499  }
1500  }
1501 
1502  return SUCCESSFUL_RETURN;
1503 }
1504 
1505 
1506 /*
1507  * o b t a i n A u x i l i a r y W o r k i n g S e t
1508  */
1509 returnValue QProblemB::obtainAuxiliaryWorkingSet( const real_t* const xOpt, const real_t* const yOpt,
1510  const Bounds* const guessedBounds, Bounds* auxiliaryBounds
1511  ) const
1512 {
1513  int i = 0;
1514  int nV = getNV( );
1515 
1516 
1517  /* 1) Ensure that desiredBounds is allocated (and different from guessedBounds). */
1518  if ( ( auxiliaryBounds == 0 ) || ( auxiliaryBounds == guessedBounds ) )
1519  return THROWERROR( RET_INVALID_ARGUMENTS );
1520 
1521 
1522  /* 2) Setup working set for auxiliary initial QP. */
1523  if ( guessedBounds != 0 )
1524  {
1525  /* If an initial working set is specific, use it!
1526  * Moreover, add all implictly fixed variables if specified. */
1527  for( i=0; i<nV; ++i )
1528  {
1529  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
1530  if ( bounds.getType( i ) == ST_EQUALITY )
1531  {
1532  if ( auxiliaryBounds->setupBound( i,ST_LOWER ) != SUCCESSFUL_RETURN )
1533  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1534  }
1535  else
1536  #endif
1537  {
1538  if ( auxiliaryBounds->setupBound( i,guessedBounds->getStatus( i ) ) != SUCCESSFUL_RETURN )
1539  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1540  }
1541  }
1542  }
1543  else /* No initial working set specified. */
1544  {
1545  if ( ( xOpt != 0 ) && ( yOpt == 0 ) )
1546  {
1547  /* Obtain initial working set by "clipping". */
1548  for( i=0; i<nV; ++i )
1549  {
1550  if ( xOpt[i] <= lb[i] + options.boundTolerance )
1551  {
1552  if ( auxiliaryBounds->setupBound( i,ST_LOWER ) != SUCCESSFUL_RETURN )
1553  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1554  continue;
1555  }
1556 
1557  if ( xOpt[i] >= ub[i] - options.boundTolerance )
1558  {
1559  if ( auxiliaryBounds->setupBound( i,ST_UPPER ) != SUCCESSFUL_RETURN )
1560  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1561  continue;
1562  }
1563 
1564  /* Moreover, add all implictly fixed variables if specified. */
1565  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
1566  if ( bounds.getType( i ) == ST_EQUALITY )
1567  {
1568  if ( auxiliaryBounds->setupBound( i,ST_LOWER ) != SUCCESSFUL_RETURN )
1569  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1570  }
1571  else
1572  #endif
1573  {
1574  if ( auxiliaryBounds->setupBound( i,ST_INACTIVE ) != SUCCESSFUL_RETURN )
1575  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1576  }
1577  }
1578  }
1579 
1580  if ( ( xOpt == 0 ) && ( yOpt != 0 ) )
1581  {
1582  /* Obtain initial working set in accordance to sign of dual solution vector. */
1583  for( i=0; i<nV; ++i )
1584  {
1585  if ( yOpt[i] > EPS )
1586  {
1587  if ( auxiliaryBounds->setupBound( i,ST_LOWER ) != SUCCESSFUL_RETURN )
1588  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1589  continue;
1590  }
1591 
1592  if ( yOpt[i] < -EPS )
1593  {
1594  if ( auxiliaryBounds->setupBound( i,ST_UPPER ) != SUCCESSFUL_RETURN )
1595  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1596  continue;
1597  }
1598 
1599  /* Moreover, add all implictly fixed variables if specified. */
1600  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
1601  if ( bounds.getType( i ) == ST_EQUALITY )
1602  {
1603  if ( auxiliaryBounds->setupBound( i,ST_LOWER ) != SUCCESSFUL_RETURN )
1604  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1605  }
1606  else
1607  #endif
1608  {
1609  if ( auxiliaryBounds->setupBound( i,ST_INACTIVE ) != SUCCESSFUL_RETURN )
1610  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1611  }
1612  }
1613  }
1614 
1615  /* If xOpt and yOpt are null pointer and no initial working is specified,
1616  * start with empty working set (or implicitly fixed bounds only)
1617  * for auxiliary QP. */
1618  if ( ( xOpt == 0 ) && ( yOpt == 0 ) )
1619  {
1620  for( i=0; i<nV; ++i )
1621  {
1622  switch( bounds.getType( i ) )
1623  {
1624  case ST_UNBOUNDED:
1625  if ( auxiliaryBounds->setupBound( i,ST_INACTIVE ) != SUCCESSFUL_RETURN )
1626  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1627  break;
1628 
1629  /* Only add all implictly fixed variables if specified. */
1630  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
1631  case ST_EQUALITY:
1632  if ( auxiliaryBounds->setupBound( i,ST_LOWER ) != SUCCESSFUL_RETURN )
1633  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1634  break;
1635  #endif
1636 
1637  default:
1638  if ( auxiliaryBounds->setupBound( i,options.initialStatusBounds ) != SUCCESSFUL_RETURN )
1639  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
1640  break;
1641  }
1642  }
1643  }
1644  }
1645 
1646  return SUCCESSFUL_RETURN;
1647 }
1648 
1649 
1650 /*
1651  * b a c k s o l v e R
1652  */
1653 returnValue QProblemB::backsolveR( const real_t* const b, BooleanType transposed,
1654  real_t* const a
1655  ) const
1656 {
1657  /* Call standard backsolve procedure (i.e. removingBound == BT_FALSE). */
1658  return backsolveR( b,transposed,BT_FALSE,a );
1659 }
1660 
1661 
1662 /*
1663  * b a c k s o l v e R
1664  */
1665 returnValue QProblemB::backsolveR( const real_t* const b, BooleanType transposed,
1666  BooleanType removingBound,
1667  real_t* const a
1668  ) const
1669 {
1670  int i, j;
1671  int nV = getNV( );
1672  int nR = getNZ( );
1673 
1674  real_t sum;
1675 
1676  /* if backsolve is called while removing a bound, reduce nZ by one. */
1677  if ( removingBound == BT_TRUE )
1678  --nR;
1679 
1680  /* nothing to do */
1681  if ( nR <= 0 )
1682  return SUCCESSFUL_RETURN;
1683 
1684 
1685  /* Solve Ra = b, where R might be transposed. */
1686  if ( transposed == BT_FALSE )
1687  {
1688  /* solve Ra = b */
1689  for( i=(nR-1); i>=0; --i )
1690  {
1691  sum = b[i];
1692  for( j=(i+1); j<nR; ++j )
1693  sum -= RR(i,j) * a[j];
1694 
1695  if ( fabs( RR(i,i) ) >= ZERO*fabs( sum ) )
1696  a[i] = sum / RR(i,i);
1697  else
1698  return THROWERROR( RET_DIV_BY_ZERO );
1699  }
1700  }
1701  else
1702  {
1703  /* solve R^T*a = b */
1704  for( i=0; i<nR; ++i )
1705  {
1706  sum = b[i];
1707  for( j=0; j<i; ++j )
1708  sum -= RR(j,i) * a[j];
1709 
1710  if ( fabs( RR(i,i) ) >= ZERO*fabs( sum ) )
1711  a[i] = sum / RR(i,i);
1712  else
1713  return THROWERROR( RET_DIV_BY_ZERO );
1714  }
1715  }
1716 
1717  return SUCCESSFUL_RETURN;
1718 }
1719 
1720 
1721 /*
1722  * d e t e r m i n e D a t a S h i f t
1723  */
1724 returnValue QProblemB::determineDataShift( const real_t* const g_new, const real_t* const lb_new, const real_t* const ub_new,
1725  real_t* const delta_g, real_t* const delta_lb, real_t* const delta_ub,
1726  BooleanType& Delta_bB_isZero
1727  )
1728 {
1729  int i, ii;
1730  int nV = getNV( );
1731  int nFX = getNFX( );
1732 
1733  int* FX_idx;
1734  bounds.getFixed( )->getNumberArray( &FX_idx );
1735 
1736 
1737  /* 1) Calculate shift directions. */
1738  for( i=0; i<nV; ++i )
1739  delta_g[i] = g_new[i] - g[i];
1740 
1741  if ( lb_new != 0 )
1742  {
1743  for( i=0; i<nV; ++i )
1744  delta_lb[i] = lb_new[i] - lb[i];
1745  }
1746  else
1747  {
1748  /* if no lower bounds exist, assume the new lower bounds to be -infinity */
1749  for( i=0; i<nV; ++i )
1750  delta_lb[i] = -INFTY - lb[i];
1751  }
1752 
1753  if ( ub_new != 0 )
1754  {
1755  for( i=0; i<nV; ++i )
1756  delta_ub[i] = ub_new[i] - ub[i];
1757  }
1758  else
1759  {
1760  /* if no upper bounds exist, assume the new upper bounds to be infinity */
1761  for( i=0; i<nV; ++i )
1762  delta_ub[i] = INFTY - ub[i];
1763  }
1764 
1765  /* 2) Determine if active bounds are to be shifted. */
1766  Delta_bB_isZero = BT_TRUE;
1767 
1768  for ( i=0; i<nFX; ++i )
1769  {
1770  ii = FX_idx[i];
1771 
1772  if ( ( fabs( delta_lb[ii] ) > EPS ) || ( fabs( delta_ub[ii] ) > EPS ) )
1773  {
1774  Delta_bB_isZero = BT_FALSE;
1775  break;
1776  }
1777  }
1778 
1779  return SUCCESSFUL_RETURN;
1780 }
1781 
1782 
1783 
1784 /*
1785  * s e t u p Q P d a t a
1786  */
1787 returnValue QProblemB::setupQPdata( SymmetricMatrix *_H, const real_t* const _g,
1788  const real_t* const _lb, const real_t* const _ub
1789  )
1790 {
1791  int i;
1792  int nV = getNV( );
1793 
1794  /* 1) Setup Hessian matrix. */
1795  setH( _H );
1796 
1797  /* 2) Setup gradient vector. */
1798  if ( _g == 0 )
1799  return THROWERROR( RET_INVALID_ARGUMENTS );
1800  else
1801  setG( _g );
1802 
1803  /* 3) Setup lower bounds vector. */
1804  if ( _lb != 0 )
1805  setLB( _lb );
1806  else
1807  /* if no lower bounds are specified, set them to -infinity */
1808  for( i=0; i<nV; ++i )
1809  lb[i] = -INFTY;
1810 
1811  /* 4) Setup upper bounds vector. */
1812  if ( _ub != 0 )
1813  setUB( _ub );
1814  else
1815  /* if no upper bounds are specified, set them to infinity */
1816  for( i=0; i<nV; ++i )
1817  ub[i] = INFTY;
1818 
1819  return SUCCESSFUL_RETURN;
1820 }
1821 
1822 
1823 /*
1824  * s e t u p Q P d a t a
1825  */
1826 returnValue QProblemB::setupQPdata( const real_t* const _H, const real_t* const _g,
1827  const real_t* const _lb, const real_t* const _ub
1828  )
1829 {
1830  int i;
1831  int nV = getNV( );
1832 
1833  /* 1) Setup Hessian matrix. */
1834  setH( _H );
1835 
1836  /* 2) Setup gradient vector. */
1837  if ( _g == 0 )
1838  return THROWERROR( RET_INVALID_ARGUMENTS );
1839  else
1840  setG( _g );
1841 
1842  /* 3) Setup lower bounds vector. */
1843  if ( _lb != 0 )
1844  {
1845  setLB( _lb );
1846  }
1847  else
1848  {
1849  /* if no lower bounds are specified, set them to -infinity */
1850  for( i=0; i<nV; ++i )
1851  lb[i] = -INFTY;
1852  }
1853 
1854  /* 4) Setup upper bounds vector. */
1855  if ( _ub != 0 )
1856  {
1857  setUB( _ub );
1858  }
1859  else
1860  {
1861  /* if no upper bounds are specified, set them to infinity */
1862  for( i=0; i<nV; ++i )
1863  ub[i] = INFTY;
1864  }
1865 
1866  return SUCCESSFUL_RETURN;
1867 }
1868 
1869 
1870 /*
1871  * s e t u p Q P d a t a F r o m F i l e
1872  */
1873 returnValue QProblemB::setupQPdataFromFile( const char* const H_file, const char* const g_file,
1874  const char* const lb_file, const char* const ub_file
1875  )
1876 {
1877  int i;
1878  int nV = getNV( );
1879 
1880  returnValue returnvalue;
1881 
1882 
1883  /* 1) Load Hessian matrix from file. */
1884  if ( H_file != 0 )
1885  {
1886  real_t* _H = new real_t[nV * nV];
1887  returnvalue = readFromFile( _H, nV,nV, H_file );
1888  if ( returnvalue != SUCCESSFUL_RETURN )
1889  {
1890  delete[] _H;
1891  return THROWERROR( returnvalue );
1892  }
1893  setH( _H );
1894  H->doFreeMemory( );
1895  }
1896  else
1897  {
1898  real_t* _H = 0;
1899  setH( _H );
1900  }
1901 
1902  /* 2) Load gradient vector from file. */
1903  if ( g_file == 0 )
1904  return THROWERROR( RET_INVALID_ARGUMENTS );
1905 
1906  returnvalue = readFromFile( g, nV, g_file );
1907  if ( returnvalue != SUCCESSFUL_RETURN )
1908  return THROWERROR( returnvalue );
1909 
1910  /* 3) Load lower bounds vector from file. */
1911  if ( lb_file != 0 )
1912  {
1913  returnvalue = readFromFile( lb, nV, lb_file );
1914  if ( returnvalue != SUCCESSFUL_RETURN )
1915  return THROWERROR( returnvalue );
1916  }
1917  else
1918  {
1919  /* if no lower bounds are specified, set them to -infinity */
1920  for( i=0; i<nV; ++i )
1921  lb[i] = -INFTY;
1922  }
1923 
1924  /* 4) Load upper bounds vector from file. */
1925  if ( ub_file != 0 )
1926  {
1927  returnvalue = readFromFile( ub, nV, ub_file );
1928  if ( returnvalue != SUCCESSFUL_RETURN )
1929  return THROWERROR( returnvalue );
1930  }
1931  else
1932  {
1933  /* if no upper bounds are specified, set them to infinity */
1934  for( i=0; i<nV; ++i )
1935  ub[i] = INFTY;
1936  }
1937 
1938  return SUCCESSFUL_RETURN;
1939 }
1940 
1941 
1942 /*
1943  * l o a d Q P v e c t o r s F r o m F i l e
1944  */
1945 returnValue QProblemB::loadQPvectorsFromFile( const char* const g_file, const char* const lb_file, const char* const ub_file,
1946  real_t* const g_new, real_t* const lb_new, real_t* const ub_new
1947  ) const
1948 {
1949  int nV = getNV( );
1950 
1951  returnValue returnvalue;
1952 
1953 
1954  /* 1) Load gradient vector from file. */
1955  if ( ( g_file != 0 ) && ( g_new != 0 ) )
1956  {
1957  returnvalue = readFromFile( g_new, nV, g_file );
1958  if ( returnvalue != SUCCESSFUL_RETURN )
1959  return THROWERROR( returnvalue );
1960  }
1961  else
1962  {
1963  /* At least gradient vector needs to be specified! */
1964  return THROWERROR( RET_INVALID_ARGUMENTS );
1965  }
1966 
1967  /* 2) Load lower bounds vector from file. */
1968  if ( lb_file != 0 )
1969  {
1970  if ( lb_new != 0 )
1971  {
1972  returnvalue = readFromFile( lb_new, nV, lb_file );
1973  if ( returnvalue != SUCCESSFUL_RETURN )
1974  return THROWERROR( returnvalue );
1975  }
1976  else
1977  {
1978  /* If filename is given, storage must be provided! */
1979  return THROWERROR( RET_INVALID_ARGUMENTS );
1980  }
1981  }
1982 
1983  /* 3) Load upper bounds vector from file. */
1984  if ( ub_file != 0 )
1985  {
1986  if ( ub_new != 0 )
1987  {
1988  returnvalue = readFromFile( ub_new, nV, ub_file );
1989  if ( returnvalue != SUCCESSFUL_RETURN )
1990  return THROWERROR( returnvalue );
1991  }
1992  else
1993  {
1994  /* If filename is given, storage must be provided! */
1995  return THROWERROR( RET_INVALID_ARGUMENTS );
1996  }
1997  }
1998 
1999  return SUCCESSFUL_RETURN;
2000 }
2001 
2002 
2003 /*
2004  * s e t I n f e a s i b i l i t y F l a g
2005  */
2006 returnValue QProblemB::setInfeasibilityFlag( returnValue returnvalue
2007  )
2008 {
2009  infeasible = BT_TRUE;
2010 
2011  if ( options.enableFarBounds == BT_FALSE )
2012  THROWERROR( returnvalue );
2013 
2014  return returnvalue;
2015 }
2016 
2017 
2018 /*
2019  * i s C P U t i m e L i m i t E x c e e d e d
2020  */
2021 BooleanType QProblemB::isCPUtimeLimitExceeded( const real_t* const cputime,
2022  real_t starttime,
2023  int nWSR
2024  ) const
2025 {
2026  /* Always perform next QP iteration if no CPU time limit is given. */
2027  if ( cputime == 0 )
2028  return BT_FALSE;
2029 
2030  /* Always perform first QP iteration. */
2031  if ( nWSR <= 0 )
2032  return BT_FALSE;
2033 
2034  real_t elapsedTime = getCPUtime( ) - starttime;
2035  real_t timePerIteration = elapsedTime / ((real_t) nWSR);
2036 
2037  /* Determine if next QP iteration exceed CPU time limit
2038  * considering the (current) average CPU time per iteration. */
2039  if ( ( elapsedTime + timePerIteration*1.25 ) <= ( *cputime ) )
2040  return BT_FALSE;
2041  else
2042  return BT_TRUE;
2043 }
2044 
2045 
2046 /*
2047  * r e g u l a r i s e H e s s i a n
2048  */
2049 returnValue QProblemB::regulariseHessian( )
2050 {
2051  /* Do nothing if Hessian regularisation is disbaled! */
2052  if ( options.enableRegularisation == BT_FALSE )
2053  return SUCCESSFUL_RETURN;
2054 
2055  /* Regularisation of identity Hessian not possible. */
2056  if ( hessianType == HST_IDENTITY )
2057  return THROWERROR( RET_CANNOT_REGULARISE_IDENTITY );
2058 
2059  /* Determine regularisation parameter. */
2060  if ( usingRegularisation( ) == BT_TRUE )
2061  return THROWERROR( RET_HESSIAN_ALREADY_REGULARISED );
2062  else
2063  {
2064  /* Regularisation of zero Hessian is done implicitly. */
2065  if ( hessianType != HST_ZERO )
2066  if ( H->addToDiag( options.epsRegularisation ) == RET_NO_DIAGONAL_AVAILABLE )
2067  return THROWERROR( RET_CANNOT_REGULARISE_SPARSE );
2068 
2069  isRegularised = BT_TRUE;
2070  THROWINFO( RET_USING_REGULARISATION );
2071  }
2072 
2073  return SUCCESSFUL_RETURN;
2074 }
2075 
2076 
2077 
2078 /*
2079  * p e r f o r m R a t i o T e s t
2080  */
2081 returnValue QProblemB::performRatioTest( const int nIdx,
2082  const int* const idxList,
2083  const SubjectTo* const subjectTo,
2084  const real_t* const num,
2085  const real_t* const den,
2086  real_t epsNum,
2087  real_t epsDen,
2088  real_t& t,
2089  int& BC_idx
2090  ) const
2091 {
2092  int i, ii;
2093 
2094  BC_idx = -1;
2095 
2096  for( i=0; i<nIdx; ++i )
2097  {
2098  ii = idxList[i];
2099 
2100  if ( subjectTo->getType( ii ) != ST_EQUALITY )
2101  {
2102  if ( ( subjectTo->getStatus( ii ) == ST_LOWER ) || ( subjectTo->getStatus( ii ) == ST_INACTIVE ) )
2103  {
2104  if ( isBlocking( num[i],den[i],epsNum,epsDen,t ) == BT_TRUE )
2105  {
2106  t = num[i] / den[i];
2107  BC_idx = ii;
2108  }
2109  }
2110  else
2111  {
2112  if ( isBlocking( -num[i],-den[i],epsNum,epsDen,t ) == BT_TRUE )
2113  {
2114  t = num[i] / den[i];
2115  BC_idx = ii;
2116  }
2117  }
2118  }
2119  }
2120 
2121  return SUCCESSFUL_RETURN;
2122 }
2123 
2124 
2125 
2126 /*
2127  * r e l a t i v e H o m o t o p y L e n g t h
2128  */
2129 real_t QProblemB::relativeHomotopyLength(const real_t* const g_new, const real_t* const lb_new, const real_t* const ub_new)
2130 {
2131  int nV = getNV( ), i;
2132  real_t d, s, len = 0.0;
2133 
2134  /* gradient */
2135  for (i = 0; i < nV; i++)
2136  {
2137  s = fabs(g_new[i]);
2138  if (s < 1.0) s = 1.0;
2139  d = fabs(g_new[i] - g[i]) / s;
2140  if (d > len) len = d;
2141  }
2142 
2143  /* lower bounds */
2144  for (i = 0; i < nV && lb_new; i++)
2145  {
2146  s = fabs(lb_new[i]);
2147  if (s < 1.0) s = 1.0;
2148  d = fabs(lb_new[i] - lb[i]) / s;
2149  if (d > len) len = d;
2150  }
2151 
2152  /* upper bounds */
2153  for (i = 0; i < nV && ub_new; i++)
2154  {
2155  s = fabs(ub_new[i]);
2156  if (s < 1.0) s = 1.0;
2157  d = fabs(ub_new[i] - ub[i]) / s;
2158  if (d > len) len = d;
2159  }
2160 
2161  return len;
2162 }
2163 
2164 
2165 /*
2166  * p e r f o r m R a m p i n g
2167  */
2168 returnValue QProblemB::performRamping( )
2169 {
2170  int nV = getNV( ), bstat, i;
2171  real_t t, rampVal;
2172 
2173  /* ramp inactive bounds and active dual variables */
2174  for (i = 0; i < nV; i++)
2175  {
2176  switch (bounds.getType(i))
2177  {
2178  case ST_EQUALITY: lb[i] = x[i]; ub[i] = x[i]; continue; /* reestablish exact feasibility */
2179  case ST_BOUNDED: continue;
2180  default: break;
2181  }
2182 
2183  t = static_cast<real_t>(i) / static_cast<real_t>(nV-1);
2184  rampVal = (1.0-t) * ramp0 + t * ramp1;
2185  bstat = bounds.getStatus(i);
2186  if (bstat != ST_LOWER) { lb[i] = x[i] - rampVal; }
2187  if (bstat != ST_UPPER) { ub[i] = x[i] + rampVal; }
2188  if (bstat == ST_LOWER) { lb[i] = x[i]; y[i] = +rampVal; }
2189  if (bstat == ST_UPPER) { ub[i] = x[i]; y[i] = -rampVal; }
2190  if (bstat == ST_INACTIVE) y[i] = 0.0; /* reestablish exact complementarity */
2191  }
2192 
2193  /* reestablish exact stationarity */
2194  setupAuxiliaryQPgradient( );
2195 
2196  /* swap ramp direction */
2197  t = ramp0; ramp0 = ramp1; ramp1 = t;
2198 
2199  return SUCCESSFUL_RETURN;
2200 }
2201 
2202 
2203 
2204 /*****************************************************************************
2205  * P R I V A T E *
2206  *****************************************************************************/
2207 
2208 /*
2209  * s o l v e I n i t i a l Q P
2210  */
2211 returnValue QProblemB::solveInitialQP( const real_t* const xOpt, const real_t* const yOpt,
2212  const Bounds* const guessedBounds,
2213  int& nWSR, real_t* const cputime
2214  )
2215 {
2216  int i;
2217  int nV = getNV( );
2218 
2219 
2220  /* start runtime measurement */
2221  real_t starttime = 0.0;
2222  if ( cputime != 0 )
2223  starttime = getCPUtime( );
2224 
2225 
2226  status = QPS_NOTINITIALISED;
2227 
2228  /* I) ANALYSE QP DATA: */
2229  /* 1) Check if Hessian happens to be the identity matrix. */
2230  if ( determineHessianType( ) != SUCCESSFUL_RETURN )
2231  return THROWERROR( RET_INIT_FAILED );
2232 
2233  /* 2) Setup type of bounds (i.e. unbounded, implicitly fixed etc.). */
2234  if ( setupSubjectToType( ) != SUCCESSFUL_RETURN )
2235  return THROWERROR( RET_INIT_FAILED );
2236 
2237  status = QPS_PREPARINGAUXILIARYQP;
2238 
2239 
2240  /* II) SETUP AUXILIARY QP WITH GIVEN OPTIMAL SOLUTION: */
2241  /* 1) Setup bounds data structure. */
2242  if ( bounds.setupAllFree( ) != SUCCESSFUL_RETURN )
2243  return THROWERROR( RET_INIT_FAILED );
2244 
2245  /* 2) Setup optimal primal/dual solution for auxiliary QP. */
2246  if ( setupAuxiliaryQPsolution( xOpt,yOpt ) != SUCCESSFUL_RETURN )
2247  return THROWERROR( RET_INIT_FAILED );
2248 
2249  /* 3) Obtain linear independent working set for auxiliary QP. */
2250  Bounds auxiliaryBounds( nV );
2251  if ( obtainAuxiliaryWorkingSet( xOpt,yOpt,guessedBounds, &auxiliaryBounds ) != SUCCESSFUL_RETURN )
2252  return THROWERROR( RET_INIT_FAILED );
2253 
2254  /* 4) Setup working set of auxiliary QP and setup cholesky decomposition. */
2255  /* a) Working set of auxiliary QP. */
2256  if ( setupAuxiliaryWorkingSet( &auxiliaryBounds,BT_TRUE ) != SUCCESSFUL_RETURN )
2257  return THROWERROR( RET_INIT_FAILED );
2258 
2259  /* b) Regularise Hessian if necessary. */
2260  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_SEMIDEF ) )
2261  {
2262  if ( regulariseHessian( ) != SUCCESSFUL_RETURN )
2263  return THROWERROR( RET_INIT_FAILED_REGULARISATION );
2264  }
2265 
2266  haveCholesky = BT_FALSE;
2267 
2268  /* c) Cholesky decomposition. */
2269  returnValue returnvalueCholesky = setupCholeskyDecomposition( );
2270 
2271  /* If Hessian is not positive definite, regularise and try again. */
2272  if ( returnvalueCholesky == RET_HESSIAN_NOT_SPD )
2273  {
2274  if ( regulariseHessian( ) != SUCCESSFUL_RETURN )
2275  return THROWERROR( RET_INIT_FAILED );
2276 
2277  if ( setupCholeskyDecomposition( ) != SUCCESSFUL_RETURN )
2278  return THROWERROR( RET_INIT_FAILED_CHOLESKY );
2279  }
2280 
2281  if ( returnvalueCholesky == RET_INDEXLIST_CORRUPTED )
2282  return THROWERROR( RET_INIT_FAILED_CHOLESKY );
2283 
2284  haveCholesky = BT_TRUE;
2285 
2286  /* 5) Store original QP formulation... */
2287  real_t* g_original = new real_t[nV];
2288  real_t* lb_original = new real_t[nV];
2289  real_t* ub_original = new real_t[nV];
2290 
2291  for( i=0; i<nV; ++i )
2292  {
2293  g_original[i] = g[i];
2294  lb_original[i] = lb[i];
2295  ub_original[i] = ub[i];
2296  }
2297 
2298  /* ... and setup QP data of an auxiliary QP having an optimal solution
2299  * as specified by the user (or xOpt = yOpt = 0, by default). */
2300  if ( setupAuxiliaryQPgradient( ) != SUCCESSFUL_RETURN )
2301  {
2302  delete[] ub_original; delete[] lb_original; delete[] g_original;
2303  return THROWERROR( RET_INIT_FAILED );
2304  }
2305 
2306  if ( setupAuxiliaryQPbounds( BT_TRUE ) != SUCCESSFUL_RETURN )
2307  {
2308  delete[] ub_original; delete[] lb_original; delete[] g_original;
2309  return THROWERROR( RET_INIT_FAILED );
2310  }
2311 
2312  status = QPS_AUXILIARYQPSOLVED;
2313 
2314 
2315  /* III) SOLVE ACTUAL INITIAL QP: */
2316 
2317  /* Allow only remaining CPU time for usual hotstart. */
2318  if ( cputime != 0 )
2319  *cputime -= getCPUtime( ) - starttime;
2320 
2321  /* Use hotstart method to find the solution of the original initial QP,... */
2322  returnValue returnvalue = hotstart( g_original,lb_original,ub_original, nWSR,cputime );
2323 
2324  /* ... deallocate memory,... */
2325  delete[] ub_original; delete[] lb_original; delete[] g_original;
2326 
2327  /* ... check for infeasibility and unboundedness... */
2328  if ( isInfeasible( ) == BT_TRUE )
2329  return THROWERROR( RET_INIT_FAILED_INFEASIBILITY );
2330 
2331  if ( isUnbounded( ) == BT_TRUE )
2332  return THROWERROR( RET_INIT_FAILED_UNBOUNDEDNESS );
2333 
2334  /* ... and internal errors. */
2335  if ( ( returnvalue != SUCCESSFUL_RETURN ) && ( returnvalue != RET_MAX_NWSR_REACHED ) )
2336  return THROWERROR( RET_INIT_FAILED_HOTSTART );
2337 
2338 
2339  /* stop runtime measurement */
2340  if ( cputime != 0 )
2341  *cputime = getCPUtime( ) - starttime;
2342 
2343  THROWINFO( RET_INIT_SUCCESSFUL );
2344 
2345  return returnvalue;
2346 }
2347 
2348 
2349 /*
2350  * s o l v e Q P
2351  */
2352 returnValue QProblemB::solveQP( const real_t* const g_new,
2353  const real_t* const lb_new, const real_t* const ub_new,
2354  int& nWSR, real_t* const cputime, int nWSRperformed
2355  )
2356 {
2357  int iter;
2358  int nV = getNV( );
2359 
2360  /* consistency check */
2361  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
2362  ( getStatus( ) == QPS_PREPARINGAUXILIARYQP ) ||
2363  ( getStatus( ) == QPS_PERFORMINGHOMOTOPY ) )
2364  {
2365  return THROWERROR( RET_HOTSTART_FAILED_AS_QP_NOT_INITIALISED );
2366  }
2367 
2368  /* start runtime measurement */
2369  real_t starttime = 0.0;
2370  if ( cputime != 0 )
2371  starttime = getCPUtime( );
2372 
2373 
2374  /* I) PREPARATIONS */
2375  /* 1) Allocate delta vectors of gradient and bounds,
2376  * index arrays and step direction arrays. */
2377  real_t* delta_xFR = new real_t[nV];
2378  real_t* delta_xFX = new real_t[nV];
2379  real_t* delta_yFX = new real_t[nV];
2380 
2381  real_t* delta_g = new real_t[nV];
2382  real_t* delta_lb = new real_t[nV];
2383  real_t* delta_ub = new real_t[nV];
2384 
2385  returnValue returnvalue;
2386  BooleanType Delta_bB_isZero;
2387 
2388  int BC_idx;
2389  SubjectToStatus BC_status;
2390 
2391  real_t homotopyLength;
2392 
2393  char messageString[80];
2394 
2395  /* 2) Update type of bounds, e.g. a formerly implicitly fixed
2396  * variable might have become a normal one etc. */
2397  if ( setupSubjectToType( lb_new,ub_new ) != SUCCESSFUL_RETURN )
2398  return THROWERROR( RET_HOTSTART_FAILED );
2399 
2400  /* 3) Reset status flags. */
2401  infeasible = BT_FALSE;
2402  unbounded = BT_FALSE;
2403 
2404 
2405  /* II) MAIN HOMOTOPY LOOP */
2406  for( iter=nWSRperformed; iter<nWSR; ++iter )
2407  {
2408  if ( isCPUtimeLimitExceeded( cputime,starttime,iter-nWSRperformed ) == BT_TRUE )
2409  {
2410  /* Assign number of working set recalculations and stop runtime measurement. */
2411  nWSR = iter;
2412  if ( cputime != 0 )
2413  *cputime = getCPUtime( ) - starttime;
2414 
2415  break;
2416  }
2417 
2418  status = QPS_PERFORMINGHOMOTOPY;
2419 
2420  #ifndef __XPCTARGET__
2421  snprintf( messageString,80,"%d ...",iter );
2422  getGlobalMessageHandler( )->throwInfo( RET_ITERATION_STARTED,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
2423  #endif
2424 
2425  /* 2) Initialise shift direction of the gradient and the bounds. */
2426  returnvalue = determineDataShift( g_new,lb_new,ub_new,
2427  delta_g,delta_lb,delta_ub,
2428  Delta_bB_isZero
2429  );
2430  if ( returnvalue != SUCCESSFUL_RETURN )
2431  {
2432  delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
2433  delete[] delta_ub; delete[] delta_lb; delete[] delta_g;
2434 
2435  /* Assign number of working set recalculations and stop runtime measurement. */
2436  nWSR = iter;
2437  if ( cputime != 0 )
2438  *cputime = getCPUtime( ) - starttime;
2439 
2440  THROWERROR( RET_SHIFT_DETERMINATION_FAILED );
2441  return returnvalue;
2442  }
2443 
2444  /* 3) Determination of step direction of X and Y. */
2445  returnvalue = determineStepDirection( delta_g,delta_lb,delta_ub,
2446  Delta_bB_isZero,
2447  delta_xFX,delta_xFR,delta_yFX
2448  );
2449  if ( returnvalue != SUCCESSFUL_RETURN )
2450  {
2451  delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
2452  delete[] delta_ub; delete[] delta_lb; delete[] delta_g;
2453 
2454  /* Assign number of working set recalculations and stop runtime measurement. */
2455  nWSR = iter;
2456  if ( cputime != 0 )
2457  *cputime = getCPUtime( ) - starttime;
2458 
2459  THROWERROR( RET_STEPDIRECTION_DETERMINATION_FAILED );
2460  return returnvalue;
2461  }
2462 
2463 
2464  /* 4) Determination of step length TAU.
2465  * This step along the homotopy path is also taken (without changing working set). */
2466  returnvalue = performStep( delta_g,delta_lb,delta_ub,
2467  delta_xFX,delta_xFR,delta_yFX,
2468  BC_idx,BC_status
2469  );
2470  if ( returnvalue != SUCCESSFUL_RETURN )
2471  {
2472  delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
2473  delete[] delta_ub; delete[] delta_lb; delete[] delta_g;
2474 
2475  /* Assign number of working set recalculations and stop runtime measurement. */
2476  nWSR = iter;
2477  if ( cputime != 0 )
2478  *cputime = getCPUtime( ) - starttime;
2479 
2480  THROWERROR( RET_STEPLENGTH_DETERMINATION_FAILED );
2481  return returnvalue;
2482  }
2483 
2484  /* 5) Termination criterion. */
2485  homotopyLength = relativeHomotopyLength(g_new, lb_new, ub_new);
2486  if ( homotopyLength <= options.terminationTolerance )
2487  {
2488  status = QPS_SOLVED;
2489 
2490  THROWINFO( RET_OPTIMAL_SOLUTION_FOUND );
2491 
2492  if ( printIteration( iter,BC_idx,BC_status ) != SUCCESSFUL_RETURN )
2493  THROWERROR( RET_PRINT_ITERATION_FAILED ); /* do not pass this as return value! */
2494 
2495  nWSR = iter;
2496  if ( cputime != 0 )
2497  *cputime = getCPUtime( ) - starttime;
2498 
2499  delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
2500  delete[] delta_ub; delete[] delta_lb; delete[] delta_g;
2501 
2502  return SUCCESSFUL_RETURN;
2503  }
2504 
2505 
2506  /* 6) Change active set. */
2507  returnvalue = changeActiveSet( BC_idx,BC_status );
2508 
2509  if ( returnvalue != SUCCESSFUL_RETURN )
2510  {
2511  delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
2512  delete[] delta_ub; delete[] delta_lb; delete[] delta_g;
2513 
2514  /* Assign number of working set recalculations and stop runtime measurement. */
2515  nWSR = iter;
2516  if ( cputime != 0 )
2517  *cputime = getCPUtime( ) - starttime;
2518 
2519  /* checks for infeasibility... */
2520  if ( infeasible == BT_TRUE )
2521  {
2522  status = QPS_HOMOTOPYQPSOLVED;
2523  return setInfeasibilityFlag( RET_HOTSTART_STOPPED_INFEASIBILITY );
2524  }
2525 
2526  /* ...unboundedness... */
2527  if ( unbounded == BT_TRUE ) /* not necessary since objective function convex! */
2528  return THROWERROR( RET_HOTSTART_STOPPED_UNBOUNDEDNESS );
2529 
2530  /* ... and throw unspecific error otherwise */
2531  THROWERROR( RET_HOMOTOPY_STEP_FAILED );
2532  return returnvalue;
2533  }
2534 
2535  /* 7) Perform Ramping Strategy on zero homotopy step or drift correction (if desired). */
2536  if ( ( tau <= EPS ) && ( options.enableRamping == BT_TRUE ) )
2537  performRamping( );
2538  else
2539  if ( (options.enableDriftCorrection > 0)
2540  && ((iter+1) % options.enableDriftCorrection == 0) )
2541  performDriftCorrection( ); /* always returns SUCCESSFUL_RETURN */
2542 
2543  /* 8) Output information of successful QP iteration. */
2544  status = QPS_HOMOTOPYQPSOLVED;
2545 
2546  if ( printIteration( iter,BC_idx,BC_status ) != SUCCESSFUL_RETURN )
2547  THROWERROR( RET_PRINT_ITERATION_FAILED ); /* do not pass this as return value! */
2548  }
2549 
2550  delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
2551  delete[] delta_ub; delete[] delta_lb; delete[] delta_g;
2552 
2553  /* stop runtime measurement */
2554  if ( cputime != 0 )
2555  *cputime = getCPUtime( ) - starttime;
2556 
2557 
2558  /* if programm gets to here, output information that QP could not be solved
2559  * within the given maximum numbers of working set changes */
2560  if ( options.printLevel == PL_HIGH )
2561  {
2562  #ifndef __XPCTARGET__
2563  snprintf( messageString,80,"(nWSR = %d)",iter );
2564  return getGlobalMessageHandler( )->throwWarning( RET_MAX_NWSR_REACHED,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
2565  #else
2566  return RET_MAX_NWSR_REACHED;
2567  #endif
2568  }
2569  else
2570  {
2571  return RET_MAX_NWSR_REACHED;
2572  }
2573 }
2574 
2575 
2576 /*
2577  * s o l v e R e g u l a r i s e d Q P
2578  */
2579 returnValue QProblemB::solveRegularisedQP( const real_t* const g_new,
2580  const real_t* const lb_new, const real_t* const ub_new,
2581  int& nWSR, real_t* const cputime, int nWSRperformed
2582  )
2583 {
2584  int i, step;
2585  int nV = getNV( );
2586 
2587 
2588  /* Stop here if QP has not been regularised (i.e. normal QP solution). */
2589  if ( usingRegularisation( ) == BT_FALSE )
2590  return solveQP( g_new,lb_new,ub_new, nWSR,cputime,nWSRperformed );
2591 
2592 
2593  /* I) SOLVE USUAL REGULARISED QP */
2594  returnValue returnvalue;
2595 
2596  int nWSR_max = nWSR;
2597  int nWSR_total = nWSRperformed;
2598 
2599  real_t cputime_total = 0.0;
2600  real_t cputime_cur = 0.0;
2601 
2602  if ( cputime == 0 )
2603  {
2604  returnvalue = solveQP( g_new,lb_new,ub_new, nWSR,0,nWSRperformed );
2605  }
2606  else
2607  {
2608  cputime_cur = *cputime;
2609  returnvalue = solveQP( g_new,lb_new,ub_new, nWSR,&cputime_cur,nWSRperformed );
2610  }
2611  nWSR_total = nWSR;
2612  cputime_total += cputime_cur;
2613 
2614 
2615  /* Only continue if QP solution has been successful. */
2616  if ( returnvalue != SUCCESSFUL_RETURN )
2617  {
2618  if ( cputime != 0 )
2619  *cputime = cputime_total;
2620 
2621  if ( returnvalue == RET_MAX_NWSR_REACHED )
2622  THROWWARNING( RET_NO_REGSTEP_NWSR );
2623 
2624  return returnvalue;
2625  }
2626 
2627 
2628  /* II) PERFORM SUCCESSIVE REGULARISATION STEPS */
2629  real_t* gMod = new real_t[nV];
2630 
2631  for( step=0; step<options.numRegularisationSteps; ++step )
2632  {
2633  /* 1) Modify gradient: gMod = g - eps*xOpt
2634  * (assuming regularisation matrix to be eps*Id). */
2635  for( i=0; i<nV; ++i )
2636  gMod[i] = g_new[i] - options.epsRegularisation*x[i];
2637 
2638  /* 2) Solve regularised QP with modified gradient allowing
2639  * only as many working set recalculations and CPU time
2640  * as have been left from previous QP solutions. */
2641  if ( cputime == 0 )
2642  {
2643  nWSR = nWSR_max;
2644  returnvalue = solveQP( gMod,lb_new,ub_new, nWSR,0,nWSR_total );
2645  }
2646  else
2647  {
2648  nWSR = nWSR_max;
2649  cputime_cur = *cputime - cputime_total;
2650  returnvalue = solveQP( gMod,lb_new,ub_new, nWSR,&cputime_cur,nWSR_total );
2651  }
2652 
2653  nWSR_total = nWSR;
2654  cputime_total += cputime_cur;
2655 
2656  /* Only continue if QP solution has been successful. */
2657  if ( returnvalue != SUCCESSFUL_RETURN )
2658  {
2659  delete[] gMod;
2660 
2661  if ( cputime != 0 )
2662  *cputime = cputime_total;
2663 
2664  if ( returnvalue == RET_MAX_NWSR_REACHED )
2665  THROWWARNING( RET_FEWER_REGSTEPS_NWSR );
2666 
2667  return returnvalue;
2668  }
2669  }
2670 
2671  delete[] gMod;
2672 
2673  if ( cputime != 0 )
2674  *cputime = cputime_total;
2675 
2676  return SUCCESSFUL_RETURN;
2677 }
2678 
2679 
2680 /*
2681  * s e t u p A u x i l i a r y W o r k i n g S e t
2682  */
2683 returnValue QProblemB::setupAuxiliaryWorkingSet( const Bounds* const auxiliaryBounds,
2684  BooleanType setupAfresh
2685  )
2686 {
2687  int i;
2688  int nV = getNV( );
2689 
2690  /* consistency checks */
2691  if ( auxiliaryBounds != 0 )
2692  {
2693  for( i=0; i<nV; ++i )
2694  if ( ( bounds.getStatus( i ) == ST_UNDEFINED ) || ( auxiliaryBounds->getStatus( i ) == ST_UNDEFINED ) )
2695  return THROWERROR( RET_UNKNOWN_BUG );
2696  }
2697  else
2698  {
2699  return THROWERROR( RET_INVALID_ARGUMENTS );
2700  }
2701 
2702 
2703  /* I) SETUP CHOLESKY FLAG:
2704  * Cholesky decomposition shall only be updated if working set
2705  * shall be updated (i.e. NOT setup afresh!) */
2706  BooleanType updateCholesky;
2707  if ( setupAfresh == BT_TRUE )
2708  updateCholesky = BT_FALSE;
2709  else
2710  updateCholesky = BT_TRUE;
2711 
2712 
2713  /* II) REMOVE FORMERLY ACTIVE BOUNDS (IF NECESSARY): */
2714  if ( setupAfresh == BT_FALSE )
2715  {
2716  /* Remove all active bounds that shall be inactive AND
2717  * all active bounds that are active at the wrong bound. */
2718  for( i=0; i<nV; ++i )
2719  {
2720  if ( ( bounds.getStatus( i ) == ST_LOWER ) && ( auxiliaryBounds->getStatus( i ) != ST_LOWER ) )
2721  if ( removeBound( i,updateCholesky ) != SUCCESSFUL_RETURN )
2722  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2723 
2724  if ( ( bounds.getStatus( i ) == ST_UPPER ) && ( auxiliaryBounds->getStatus( i ) != ST_UPPER ) )
2725  if ( removeBound( i,updateCholesky ) != SUCCESSFUL_RETURN )
2726  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2727  }
2728  }
2729 
2730 
2731  /* III) ADD NEWLY ACTIVE BOUNDS: */
2732  /* Add all inactive bounds that shall be active AND
2733  * all formerly active bounds that have been active at the wrong bound. */
2734  for( i=0; i<nV; ++i )
2735  {
2736  if ( ( bounds.getStatus( i ) == ST_INACTIVE ) && ( auxiliaryBounds->getStatus( i ) != ST_INACTIVE ) )
2737  {
2738  if ( addBound( i,auxiliaryBounds->getStatus( i ),updateCholesky ) != SUCCESSFUL_RETURN )
2739  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2740  }
2741  }
2742 
2743  return SUCCESSFUL_RETURN;
2744 }
2745 
2746 
2747 /*
2748  * s e t u p A u x i l i a r y Q P s o l u t i o n
2749  */
2750 returnValue QProblemB::setupAuxiliaryQPsolution( const real_t* const xOpt, const real_t* const yOpt
2751  )
2752 {
2753  int i;
2754  int nV = getNV( );
2755 
2756 
2757  /* Setup primal/dual solution vectors for auxiliary initial QP:
2758  * if a null pointer is passed, a zero vector is assigned;
2759  * old solution vector is kept if pointer to internal solution vector is passed. */
2760  if ( xOpt != 0 )
2761  {
2762  if ( xOpt != x )
2763  for( i=0; i<nV; ++i )
2764  x[i] = xOpt[i];
2765  }
2766  else
2767  {
2768  for( i=0; i<nV; ++i )
2769  x[i] = 0.0;
2770  }
2771 
2772  if ( yOpt != 0 )
2773  {
2774  if ( yOpt != y )
2775  for( i=0; i<nV; ++i )
2776  y[i] = yOpt[i];
2777  }
2778  else
2779  {
2780  for( i=0; i<nV; ++i )
2781  y[i] = 0.0;
2782  }
2783 
2784  return SUCCESSFUL_RETURN;
2785 }
2786 
2787 
2788 /*
2789  * s e t u p A u x i l i a r y Q P g r a d i e n t
2790  */
2791 returnValue QProblemB::setupAuxiliaryQPgradient( )
2792 {
2793  int i;
2794  int nV = getNV( );
2795 
2796  /* Setup gradient vector: g = -H*x + y'*Id. */
2797  switch ( hessianType )
2798  {
2799  case HST_ZERO:
2800  if ( usingRegularisation( ) == BT_FALSE )
2801  for ( i=0; i<nV; ++i )
2802  g[i] = y[i];
2803  else
2804  for ( i=0; i<nV; ++i )
2805  g[i] = y[i] - options.epsRegularisation*x[i];
2806  break;
2807 
2808  case HST_IDENTITY:
2809  for ( i=0; i<nV; ++i )
2810  g[i] = y[i] - x[i];
2811  break;
2812 
2813  default:
2814  /* y'*Id */
2815  for ( i=0; i<nV; ++i )
2816  g[i] = y[i];
2817 
2818  /* -H*x */
2819  H->times(1, -1.0, x, nV, 1.0, g, nV);
2820 
2821  break;
2822  }
2823 
2824  return SUCCESSFUL_RETURN;
2825 }
2826 
2827 
2828 /*
2829  * s e t u p A u x i l i a r y Q P b o u n d s
2830  */
2831 returnValue QProblemB::setupAuxiliaryQPbounds( BooleanType useRelaxation )
2832 {
2833  int i;
2834  int nV = getNV( );
2835 
2836 
2837  /* Setup bound vectors. */
2838  for ( i=0; i<nV; ++i )
2839  {
2840  switch ( bounds.getStatus( i ) )
2841  {
2842  case ST_INACTIVE:
2843  if ( useRelaxation == BT_TRUE )
2844  {
2845  if ( bounds.getType( i ) == ST_EQUALITY )
2846  {
2847  lb[i] = x[i];
2848  ub[i] = x[i];
2849  }
2850  else
2851  {
2852  lb[i] = x[i] - options.boundRelaxation;
2853  ub[i] = x[i] + options.boundRelaxation;
2854  }
2855  }
2856  break;
2857 
2858  case ST_LOWER:
2859  lb[i] = x[i];
2860  if ( bounds.getType( i ) == ST_EQUALITY )
2861  {
2862  ub[i] = x[i];
2863  }
2864  else
2865  {
2866  if ( useRelaxation == BT_TRUE )
2867  ub[i] = x[i] + options.boundRelaxation;
2868  }
2869  break;
2870 
2871  case ST_UPPER:
2872  ub[i] = x[i];
2873  if ( bounds.getType( i ) == ST_EQUALITY )
2874  {
2875  lb[i] = x[i];
2876  }
2877  else
2878  {
2879  if ( useRelaxation == BT_TRUE )
2880  lb[i] = x[i] - options.boundRelaxation;
2881  }
2882  break;
2883 
2884  default:
2885  return THROWERROR( RET_UNKNOWN_BUG );
2886  }
2887  }
2888 
2889  return SUCCESSFUL_RETURN;
2890 }
2891 
2892 
2893 /*
2894  * s e t u p A u x i l i a r y Q P
2895  */
2896 returnValue QProblemB::setupAuxiliaryQP( const Bounds* const guessedBounds )
2897 {
2898  int i;
2899  int nV = getNV( );
2900 
2901  /* nothing to do */
2902  if ( guessedBounds == &bounds )
2903  return SUCCESSFUL_RETURN;
2904 
2905  status = QPS_PREPARINGAUXILIARYQP;
2906 
2907 
2908  /* I) SETUP WORKING SET ... */
2909  if ( shallRefactorise( guessedBounds ) == BT_TRUE )
2910  {
2911  /* ... WITH REFACTORISATION: */
2912  /* 1) Reset bounds ... */
2913  bounds.init( nV );
2914 
2915  /* ... and set them up afresh. */
2916  if ( setupSubjectToType( ) != SUCCESSFUL_RETURN )
2917  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2918 
2919  if ( bounds.setupAllFree( ) != SUCCESSFUL_RETURN )
2920  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2921 
2922  /* 2) Setup guessed working set afresh. */
2923  if ( setupAuxiliaryWorkingSet( guessedBounds,BT_TRUE ) != SUCCESSFUL_RETURN )
2924  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2925 
2926  /* 3) Calculate Cholesky decomposition. */
2927  if ( setupCholeskyDecomposition( ) != SUCCESSFUL_RETURN )
2928  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2929  }
2930  else
2931  {
2932  /* ... WITHOUT REFACTORISATION: */
2933  if ( setupAuxiliaryWorkingSet( guessedBounds,BT_FALSE ) != SUCCESSFUL_RETURN )
2934  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2935  }
2936 
2937 
2938  /* II) SETUP AUXILIARY QP DATA: */
2939  /* 1) Ensure that dual variable is zero for fixed bounds. */
2940  for ( i=0; i<nV; ++i )
2941  if ( bounds.getStatus( i ) != ST_INACTIVE )
2942  y[i] = 0.0;
2943 
2944  /* 2) Setup gradient and bound vectors. */
2945  if ( setupAuxiliaryQPgradient( ) != SUCCESSFUL_RETURN )
2946  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2947 
2948  if ( setupAuxiliaryQPbounds( BT_FALSE ) != SUCCESSFUL_RETURN )
2949  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
2950 
2951  return SUCCESSFUL_RETURN;
2952 }
2953 
2954 
2955 /*
2956  * d e t e r m i n e S t e p D i r e c t i o n
2957  */
2958 returnValue QProblemB::determineStepDirection( const real_t* const delta_g, const real_t* const delta_lb, const real_t* const delta_ub,
2959  BooleanType Delta_bB_isZero,
2960  real_t* const delta_xFX, real_t* const delta_xFR,
2961  real_t* const delta_yFX
2962  )
2963 {
2964  int i, ii;
2965  int r;
2966  int nFR = getNFR( );
2967  int nFX = getNFX( );
2968 
2969  int* FR_idx;
2970  int* FX_idx;
2971 
2972  bounds.getFree( )->getNumberArray( &FR_idx );
2973  bounds.getFixed( )->getNumberArray( &FX_idx );
2974 
2975 
2976  /* This routine computes
2977  * delta_xFX := delta_b
2978  * delta_xFR := R \ R' \ -( delta_g + HMX*delta_xFX )
2979  * delta_yFX := HMX'*delta_xFR + HFX*delta_xFX { + eps*delta_xFX }
2980  */
2981 
2982  /* I) DETERMINE delta_xFX := delta_{l|u}b */
2983  if ( Delta_bB_isZero == BT_FALSE )
2984  {
2985  for( i=0; i<nFX; ++i )
2986  {
2987  ii = FX_idx[i];
2988 
2989  if ( bounds.getStatus( ii ) == ST_LOWER )
2990  delta_xFX[i] = delta_lb[ii];
2991  else
2992  delta_xFX[i] = delta_ub[ii];
2993  }
2994  }
2995  else
2996  {
2997  for( i=0; i<nFX; ++i )
2998  delta_xFX[i] = 0.0;
2999  }
3000 
3001 
3002  /* delta_xFR_TMP holds the residual, initialized with right hand side
3003  * delta_xFR holds the step that gets refined incrementally */
3004  for ( i=0; i<nFR; ++i )
3005  {
3006  ii = FR_idx[i];
3007  delta_xFR_TMP[i] = - delta_g[ii];
3008  delta_xFR[i] = 0.0;
3009  }
3010 
3011 
3012  /* Iterative refinement loop for delta_xFR */
3013  for ( r=0; r<=options.numRefinementSteps; ++r )
3014  {
3015  /* II) DETERMINE delta_xFR */
3016  if ( nFR > 0 )
3017  {
3018  /* Add - HMX*delta_xFX
3019  * This is skipped if delta_b=0 or mixed part HM=0 (H=0 or H=Id) */
3020  if ( ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) && ( Delta_bB_isZero == BT_FALSE ) && ( r == 0 ) )
3021  H->times(bounds.getFree(), bounds.getFixed(), 1, -1.0, delta_xFX, nFX, 1.0, delta_xFR_TMP, nFR);
3022 
3023  /* Determine R' \ ( - HMX*delta_xFX - delta_gFR ) where R'R = HFR */
3024  if ( backsolveR( delta_xFR_TMP,BT_TRUE,delta_xFR_TMP ) != SUCCESSFUL_RETURN )
3025  return THROWERROR( RET_STEPDIRECTION_FAILED_CHOLESKY );
3026 
3027  /* Determine HFR \ ( - HMX*delta_xFX - delta_gFR ) */
3028  if ( backsolveR( delta_xFR_TMP,BT_FALSE,delta_xFR_TMP ) != SUCCESSFUL_RETURN )
3029  return THROWERROR( RET_STEPDIRECTION_FAILED_CHOLESKY );
3030  }
3031 
3032  /* refine solution found for delta_xFR so far */
3033  for ( i=0; i<nFR; ++i )
3034  delta_xFR[i] += delta_xFR_TMP[i];
3035 
3036  if ( options.numRefinementSteps > 0 )
3037  {
3038  real_t rnrm = 0.0;
3039  /* compute new residual in delta_xFR_TMP:
3040  * residual := - HFR*delta_xFR - HMX*delta_xFX - delta_gFR
3041  * set to -delta_gFR */
3042  for ( i=0; i<nFR; ++i )
3043  {
3044  ii = FR_idx[i];
3045  delta_xFR_TMP[i] = -delta_g[ii];
3046  }
3047  /* add - HFR*delta_xFR */
3048  switch ( hessianType )
3049  {
3050  case HST_ZERO:
3051  break;
3052 
3053  case HST_IDENTITY:
3054  for ( i=0; i<nFR; ++i )
3055  {
3056  delta_xFR_TMP[i] -= delta_xFR[i];
3057 
3058  /* compute max norm */
3059  if (rnrm < fabs (delta_xFR_TMP[i]))
3060  rnrm = fabs (delta_xFR_TMP[i]);
3061  }
3062  break;
3063 
3064  default:
3065  H->times(bounds.getFree(), bounds.getFree(), 1, -1.0, delta_xFR, nFR, 1.0, delta_xFR_TMP, nFR);
3066  H->times(bounds.getFree(), bounds.getFixed(), 1, -1.0, delta_xFX, nFX, 1.0, delta_xFR_TMP, nFR);
3067 
3068  /* compute max norm */
3069  for ( i=0; i<nFR; ++i )
3070  if (rnrm < fabs (delta_xFR_TMP[i]))
3071  rnrm = fabs (delta_xFR_TMP[i]);
3072 
3073  break;
3074  }
3075 
3076  /* early termination of residual norm small enough */
3077  if ( rnrm < options.epsIterRef )
3078  break;
3079  }
3080 
3081  } /* end of refinement loop for delta_xFR */
3082 
3083  /* III) DETERMINE delta_yFX */
3084  if ( nFX > 0 )
3085  {
3086  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_IDENTITY ) )
3087  {
3088  for( i=0; i<nFX; ++i )
3089  {
3090  /* set to delta_g */
3091  ii = FX_idx[i];
3092  delta_yFX[i] = delta_g[ii];
3093 
3094  /* add HFX*delta_xFX = {0|I}*delta_xFX */
3095  if ( hessianType == HST_ZERO )
3096  {
3097  if ( usingRegularisation( ) == BT_TRUE )
3098  delta_yFX[i] += options.epsRegularisation*delta_xFX[i];
3099  }
3100  else
3101  delta_yFX[i] += delta_xFX[i];
3102  }
3103  }
3104  else
3105  {
3106  for( i=0; i<nFX; ++i )
3107  {
3108  /* set to delta_g */
3109  ii = FX_idx[i];
3110  delta_yFX[i] = delta_g[ii];
3111  }
3112  H->times(bounds.getFixed(), bounds.getFree(), 1, 1.0, delta_xFR, nFR, 1.0, delta_yFX, nFX);
3113  if (Delta_bB_isZero == BT_FALSE)
3114  H->times(bounds.getFixed(), bounds.getFixed(), 1, 1.0, delta_xFX, nFX, 1.0, delta_yFX, nFX);
3115  }
3116  }
3117 
3118  return SUCCESSFUL_RETURN;
3119 }
3120 
3121 
3122 /*
3123  * p e r f o r m S t e p
3124  */
3125 returnValue QProblemB::performStep( const real_t* const delta_g,
3126  const real_t* const delta_lb, const real_t* const delta_ub,
3127  const real_t* const delta_xFX,
3128  const real_t* const delta_xFR,
3129  const real_t* const delta_yFX,
3130  int& BC_idx, SubjectToStatus& BC_status
3131  )
3132 {
3133  int i, ii;
3134  int nV = getNV( );
3135  int nFR = getNFR( );
3136  int nFX = getNFX( );
3137 
3138  int* FR_idx;
3139  int* FX_idx;
3140 
3141  bounds.getFree( )->getNumberArray( &FR_idx );
3142  bounds.getFixed( )->getNumberArray( &FX_idx );
3143 
3144  tau = 1.0;
3145  BC_idx = -1;
3146  BC_status = ST_UNDEFINED;
3147 
3148  int BC_idx_tmp = -1;
3149 
3150  real_t* num = new real_t[nV];
3151  real_t* den = new real_t[nV];
3152 
3153 
3154  /* I) DETERMINE MAXIMUM DUAL STEPLENGTH, i.e. ensure that
3155  * active dual bounds remain valid (ignoring implicitly fixed variables): */
3156  for( i=0; i<nFX; ++i )
3157  {
3158  ii = FX_idx[i];
3159  num[i] = y[ii];
3160  den[i] = -delta_yFX[i];
3161  }
3162 
3163  performRatioTest( nFX,FX_idx,&bounds, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
3164 
3165  if ( BC_idx_tmp >= 0 )
3166  {
3167  BC_idx = BC_idx_tmp;
3168  BC_status = ST_INACTIVE;
3169  }
3170 
3171 
3172  /* II) DETERMINE MAXIMUM PRIMAL STEPLENGTH, i.e. ensure that
3173  * inactive bounds remain valid (ignoring unbounded variables). */
3174  /* 1) Inactive lower bounds. */
3175  if ( bounds.hasNoLower( ) == BT_FALSE )
3176  {
3177  for( i=0; i<nFR; ++i )
3178  {
3179  ii = FR_idx[i];
3180  num[i] = getMax( x[ii] - lb[ii],0.0 );
3181  den[i] = delta_lb[ii] - delta_xFR[i];
3182  }
3183 
3184  performRatioTest( nFR,FR_idx,&bounds, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
3185 
3186  if ( BC_idx_tmp >= 0 )
3187  {
3188  BC_idx = BC_idx_tmp;
3189  BC_status = ST_LOWER;
3190  }
3191  }
3192 
3193  /* 2) Inactive upper bounds. */
3194  if ( bounds.hasNoUpper( ) == BT_FALSE )
3195  {
3196  for( i=0; i<nFR; ++i )
3197  {
3198  ii = FR_idx[i];
3199  num[i] = getMax( ub[ii] - x[ii],0.0 );
3200  den[i] = delta_xFR[i] - delta_ub[ii];
3201  }
3202 
3203  performRatioTest( nFR,FR_idx,&bounds, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
3204 
3205  if ( BC_idx_tmp >= 0 )
3206  {
3207  BC_idx = BC_idx_tmp;
3208  BC_status = ST_UPPER;
3209  }
3210  }
3211 
3212  delete[] den;
3213  delete[] num;
3214 
3215 
3216  #ifndef __XPCTARGET__
3217  char messageString[80];
3218 
3219  if ( BC_status == ST_UNDEFINED )
3220  snprintf( messageString,80,"Stepsize is %.10e!",tau );
3221  else
3222  snprintf( messageString,80,"Stepsize is %.10e! (BC_idx = %d, BC_status = %d)",tau,BC_idx,BC_status );
3223 
3224  getGlobalMessageHandler( )->throwInfo( RET_STEPSIZE_NONPOSITIVE,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3225  #endif
3226 
3227 
3228  /* III) PERFORM STEP ALONG HOMOTOPY PATH */
3229  if ( tau > ZERO )
3230  {
3231  /* 1) Perform step in primal und dual space. */
3232  for( i=0; i<nFR; ++i )
3233  {
3234  ii = FR_idx[i];
3235  x[ii] += tau*delta_xFR[i];
3236  }
3237 
3238  for( i=0; i<nFX; ++i )
3239  {
3240  ii = FX_idx[i];
3241  x[ii] += tau*delta_xFX[i];
3242  y[ii] += tau*delta_yFX[i];
3243  }
3244 
3245  /* 2) Shift QP data. */
3246  for( i=0; i<nV; ++i )
3247  {
3248  g[i] += tau*delta_g[i];
3249  lb[i] += tau*delta_lb[i];
3250  ub[i] += tau*delta_ub[i];
3251  }
3252  }
3253  else
3254  {
3255  /* print a warning if stepsize is zero */
3256  #ifndef __XPCTARGET__
3257  snprintf( messageString,80,"Stepsize is %.6e",tau );
3258  getGlobalMessageHandler( )->throwWarning( RET_STEPSIZE,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3259  #endif
3260  }
3261 
3262 
3263  return SUCCESSFUL_RETURN;
3264 }
3265 
3266 
3267 /*
3268  * c h a n g e A c t i v e S e t
3269  */
3270 returnValue QProblemB::changeActiveSet( int BC_idx, SubjectToStatus BC_status )
3271 {
3272  char messageString[80];
3273 
3274  /* IV) UPDATE ACTIVE SET */
3275  switch ( BC_status )
3276  {
3277  /* Optimal solution found as no working set change detected. */
3278  case ST_UNDEFINED:
3279  return RET_OPTIMAL_SOLUTION_FOUND;
3280 
3281 
3282  /* Remove one variable from active set. */
3283  case ST_INACTIVE:
3284  #ifndef __XPCTARGET__
3285  snprintf( messageString,80,"bound no. %d.", BC_idx );
3286  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3287  #endif
3288 
3289  if ( removeBound( BC_idx,BT_TRUE ) != SUCCESSFUL_RETURN )
3290  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
3291 
3292  y[BC_idx] = 0.0;
3293  break;
3294 
3295 
3296  /* Add one variable to active set. */
3297  default:
3298  #ifndef __XPCTARGET__
3299  if ( BC_status == ST_LOWER )
3300  snprintf( messageString,80,"lower bound no. %d.", BC_idx );
3301  else
3302  snprintf( messageString,80,"upper bound no. %d.", BC_idx );
3303  getGlobalMessageHandler( )->throwInfo( RET_ADD_TO_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3304  #endif
3305 
3306  if ( addBound( BC_idx,BC_status,BT_TRUE ) != SUCCESSFUL_RETURN )
3307  return THROWERROR( RET_ADD_TO_ACTIVESET_FAILED );
3308  break;
3309  }
3310 
3311  return SUCCESSFUL_RETURN;
3312 }
3313 
3314 
3315 
3316 /*
3317  * p e r f o r m D r i f t C o r r e c t i o n
3318  */
3319 returnValue QProblemB::performDriftCorrection( )
3320 {
3321  int i;
3322  int nV = getNV ();
3323 
3324  for ( i=0; i<nV; ++i )
3325  {
3326  switch ( bounds.getType ( i ) )
3327  {
3328  case ST_BOUNDED:
3329  switch ( bounds.getStatus ( i ) )
3330  {
3331  case ST_LOWER:
3332  lb[i] = x[i];
3333  ub[i] = getMax (ub[i], x[i]);
3334  y[i] = getMax (y[i], 0.0);
3335  break;
3336  case ST_UPPER:
3337  lb[i] = getMin (lb[i], x[i]);
3338  ub[i] = x[i];
3339  y[i] = getMin (y[i], 0.0);
3340  break;
3341  case ST_INACTIVE:
3342  lb[i] = getMin (lb[i], x[i]);
3343  ub[i] = getMax (ub[i], x[i]);
3344  y[i] = 0.0;
3345  break;
3346  case ST_UNDEFINED:
3347  break;
3348  }
3349  break;
3350  case ST_EQUALITY:
3351  lb[i] = x[i];
3352  ub[i] = x[i];
3353  break;
3354  case ST_UNBOUNDED:
3355  case ST_UNKNOWN:
3356  break;
3357  }
3358  }
3359 
3360  setupAuxiliaryQPgradient( );
3361 
3362  return SUCCESSFUL_RETURN;
3363 }
3364 
3365 
3366 /*
3367  * s h a l l R e f a c t o r i s e
3368  */
3369 BooleanType QProblemB::shallRefactorise( const Bounds* const guessedBounds ) const
3370 {
3371  int i;
3372  int nV = getNV( );
3373 
3374  /* always refactorise if Hessian is not known to be positive definite */
3375  if ( getHessianType( ) == HST_SEMIDEF )
3376  return BT_TRUE;
3377 
3378  /* 1) Determine number of bounds that have same status
3379  * in guessed AND current bounds.*/
3380  int differenceNumber = 0;
3381 
3382  for( i=0; i<nV; ++i )
3383  if ( guessedBounds->getStatus( i ) != bounds.getStatus( i ) )
3384  ++differenceNumber;
3385 
3386  /* 2) Decide wheter to refactorise or not. */
3387  if ( 2*differenceNumber > guessedBounds->getNFX( ) )
3388  return BT_TRUE;
3389  else
3390  return BT_FALSE;
3391 }
3392 
3393 
3394 /*
3395  * a d d B o u n d
3396  */
3397 returnValue QProblemB::addBound( int number, SubjectToStatus B_status,
3398  BooleanType updateCholesky
3399  )
3400 {
3401  int i, j;
3402  int nV = getNV( );
3403  int nFR = getNFR( );
3404 
3405 
3406  /* consistency check */
3407  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
3408  ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
3409  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
3410  ( getStatus( ) == QPS_SOLVED ) )
3411  {
3412  return THROWERROR( RET_UNKNOWN_BUG );
3413  }
3414 
3415  /* Perform cholesky updates only if QProblemB has been initialised! */
3416  if ( getStatus( ) == QPS_PREPARINGAUXILIARYQP )
3417  {
3418  /* UPDATE INDICES */
3419  if ( bounds.moveFreeToFixed( number,B_status ) != SUCCESSFUL_RETURN )
3420  return THROWERROR( RET_ADDBOUND_FAILED );
3421 
3422  return SUCCESSFUL_RETURN;
3423  }
3424 
3425 
3426  /* I) PERFORM CHOLESKY UPDATE: */
3427  if ( ( updateCholesky == BT_TRUE ) &&
3428  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3429  {
3430  /* 1) Index of variable to be added within the list of free variables. */
3431  int number_idx = bounds.getFree( )->getIndex( number );
3432 
3433  real_t c, s, nu;
3434 
3435  /* 2) Use row-wise Givens rotations to restore upper triangular form of R. */
3436  for( i=number_idx+1; i<nFR; ++i )
3437  {
3438  computeGivens( RR(i-1,i),RR(i,i), RR(i-1,i),RR(i,i),c,s );
3439  nu = s/(1.0+c);
3440 
3441  for( j=(1+i); j<nFR; ++j ) /* last column of R is thrown away */
3442  applyGivens( c,s,nu,RR(i-1,j),RR(i,j), RR(i-1,j),RR(i,j) );
3443  }
3444 
3445  /* 3) Delete <number_idx>th column and ... */
3446  for( i=0; i<nFR-1; ++i )
3447  for( j=number_idx+1; j<nFR; ++j )
3448  RR(i,j-1) = RR(i,j);
3449  /* ... last column of R. */
3450  for( i=0; i<nFR; ++i )
3451  RR(i,nFR-1) = 0.0;
3452  }
3453 
3454  /* II) UPDATE INDICES */
3455  if ( bounds.moveFreeToFixed( number,B_status ) != SUCCESSFUL_RETURN )
3456  return THROWERROR( RET_ADDBOUND_FAILED );
3457 
3458 
3459  return SUCCESSFUL_RETURN;
3460 }
3461 
3462 
3463 /*
3464  * r e m o v e B o u n d
3465  */
3466 returnValue QProblemB::removeBound( int number,
3467  BooleanType updateCholesky
3468  )
3469 {
3470  int i;
3471  int nV = getNV( );
3472  int nFR = getNFR( );
3473 
3474 
3475  /* consistency check */
3476  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
3477  ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
3478  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
3479  ( getStatus( ) == QPS_SOLVED ) )
3480  {
3481  return THROWERROR( RET_UNKNOWN_BUG );
3482  }
3483 
3484  /* save index sets and decompositions for flipping bounds strategy */
3485  if ( options.enableFlippingBounds == BT_TRUE )
3486  flipper.set( &bounds,R );
3487 
3488  /* I) UPDATE INDICES */
3489  if ( bounds.moveFixedToFree( number ) != SUCCESSFUL_RETURN )
3490  return THROWERROR( RET_REMOVEBOUND_FAILED );
3491 
3492  /* Perform cholesky updates only if QProblemB has been initialised! */
3493  if ( getStatus( ) == QPS_PREPARINGAUXILIARYQP )
3494  return SUCCESSFUL_RETURN;
3495 
3496 
3497  /* II) PERFORM CHOLESKY UPDATE */
3498  if ( ( updateCholesky == BT_TRUE ) &&
3499  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3500  {
3501  int* FR_idx;
3502  bounds.getFree( )->getNumberArray( &FR_idx );
3503 
3504  /* 1) Calculate new column of cholesky decomposition. */
3505  real_t* rhs = new real_t[nFR+1];
3506  real_t* r = new real_t[nFR];
3507 
3508  real_t r0;
3509  switch ( hessianType )
3510  {
3511  case HST_ZERO:
3512  if ( usingRegularisation( ) == BT_FALSE )
3513  r0 = 0.0;
3514  else
3515  r0 = options.epsRegularisation;
3516  for( i=0; i<nFR; ++i )
3517  rhs[i] = 0.0;
3518  break;
3519 
3520  case HST_IDENTITY:
3521  r0 = 1.0;
3522  for( i=0; i<nFR; ++i )
3523  rhs[i] = 0.0;
3524  break;
3525 
3526  default:
3527  H->getRow(number, bounds.getFree(), 1.0, rhs);
3528  r0 = H->diag(number);
3529  break;
3530  }
3531 
3532  if ( backsolveR( rhs,BT_TRUE,BT_TRUE,r ) != SUCCESSFUL_RETURN )
3533  {
3534  delete[] rhs; delete[] r;
3535  return THROWERROR( RET_REMOVEBOUND_FAILED );
3536  }
3537 
3538  for( i=0; i<nFR; ++i )
3539  r0 -= r[i]*r[i];
3540 
3541  /* 2) Store new column into R. */
3542  for( i=0; i<nFR; ++i )
3543  RR(i,nFR) = r[i];
3544 
3545  if ( options.enableFlippingBounds == BT_TRUE )
3546  {
3547  if ( r0 > options.epsFlipping )
3548  RR(nFR,nFR) = sqrt( r0 );
3549  else
3550  {
3551  hessianType = HST_SEMIDEF;
3552 
3553  flipper.get( &bounds,R );
3554  bounds.flipFixed(number);
3555 
3556  switch (bounds.getStatus(number))
3557  {
3558  case ST_LOWER: lb[number] = ub[number]; break;
3559  case ST_UPPER: ub[number] = lb[number]; break;
3560  default: delete[] rhs; delete[] r; return THROWERROR( RET_MOVING_BOUND_FAILED );
3561  }
3562 
3563  }
3564  }
3565  else
3566  {
3567  if ( r0 > ZERO )
3568  RR(nFR,nFR) = sqrt( r0 );
3569  else
3570  {
3571  delete[] rhs; delete[] r;
3572 
3573  hessianType = HST_SEMIDEF;
3574  return THROWERROR( RET_HESSIAN_NOT_SPD );
3575  }
3576  }
3577 
3578  delete[] rhs; delete[] r;
3579  }
3580 
3581  return SUCCESSFUL_RETURN;
3582 }
3583 
3584 
3585 /*
3586  * p r i n t I t e r a t i o n
3587  */
3588 returnValue QProblemB::printIteration( int iteration,
3589  int BC_idx, SubjectToStatus BC_status
3590  )
3591 {
3592  #ifndef __XPCTARGET__
3593  char myPrintfString[160];
3594 
3595  /* consistency check */
3596  if ( iteration < 0 )
3597  return THROWERROR( RET_INVALID_ARGUMENTS );
3598 
3599  /* nothing to do */
3600  if ( options.printLevel != PL_MEDIUM )
3601  return SUCCESSFUL_RETURN;
3602 
3603 
3604  /* 1) Print header at first iteration. */
3605  if ( iteration == 0 )
3606  {
3607  snprintf( myPrintfString,160,"\n\n################# qpOASES -- QP NO. %3.0d ##################\n\n", count );
3608  myPrintf( myPrintfString );
3609 
3610  myPrintf( " Iter | StepLength | Info | nFX \n" );
3611  myPrintf( " ----------+------------------+------------------+--------- \n" );
3612  }
3613 
3614  /* 2) Print iteration line. */
3615  if ( BC_status == ST_UNDEFINED )
3616  {
3617  if ( hessianType == HST_ZERO )
3618  snprintf( myPrintfString,80," %5.1d | %1.6e | LP SOLVED | %4.1d \n", iteration,tau,getNFX( ) );
3619  else
3620  snprintf( myPrintfString,80," %5.1d | %1.6e | QP SOLVED | %4.1d \n", iteration,tau,getNFX( ) );
3621  myPrintf( myPrintfString );
3622  }
3623  else
3624  {
3625  char info[8];
3626 
3627  if ( BC_status == ST_INACTIVE )
3628  snprintf( info,8,"REM BND" );
3629  else
3630  snprintf( info,8,"ADD BND" );
3631 
3632  snprintf( myPrintfString,80," %5.1d | %1.6e | %s %4.1d | %4.1d \n", iteration,tau,info,BC_idx,getNFX( ) );
3633  myPrintf( myPrintfString );
3634  }
3635  #endif
3636 
3637  return SUCCESSFUL_RETURN;
3638 }
3639 
3640 
3641 
3642 END_NAMESPACE_QPOASES
3643 
3644 
3645 /*
3646  * end of file
3647  */
RET_REMOVEBOUND_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:146
getMax
int getMax(int x, int y)
RET_STEPSIZE_NONPOSITIVE
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:139
QProblemB::getHessianType
HessianType getHessianType() const
RET_INIT_FAILED_HOTSTART
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:103
SubjectTo::setNoLower
void setNoLower(BooleanType _status)
SubjectTo::setType
returnValue setType(int i, SubjectToType value)
sqrt
IntermediateState sqrt(const Expression &arg)
Definition: acado_syntax.cpp:50
RET_UNABLE_TO_READ_FILE
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:172
QProblemB::isInfeasible
BooleanType isInfeasible() const
RET_UNKNOWN_BUG
Definition: acado_types.hpp:845
ST_LOWER
#define ST_LOWER
Definition: qpOASES_wrapper.h:73
RET_STEPLENGTH_DETERMINATION_FAILED
Definition: acado_types.hpp:1096
Bounds::getFree
Indexlist * getFree()
QProblemB::copy
returnValue copy(const QProblemB &rhs)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1209
RET_HESSIAN_NOT_SPD
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:164
RET_SHIFT_DETERMINATION_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:121
__LINE__
#define __LINE__
Definition: acado_message_handling.hpp:98
isZero
BooleanType isZero(real_t x, real_t TOL=ZERO)
QProblemB
Implements the online active set strategy for box-constrained QPs.
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/QProblemB.hpp:55
SUCCESSFUL_RETURN
Definition: acado_types.hpp:836
RET_QPOBJECT_NOT_SETUP
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:92
RET_ITERATION_STARTED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:120
QProblemB::isUnbounded
BooleanType isUnbounded() const
INFTY
const double INFTY
Definition: acado_constants.hpp:65
RET_HOMOTOPY_STEP_FAILED
Definition: acado_types.hpp:1093
RET_PRINT_ITERATION_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:167
RET_MOVING_BOUND_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:94
RET_HOTSTART_FAILED_AS_QP_NOT_INITIALISED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:119
QProblemB::getNFX
int getNFX()
RET_INVALID_ARGUMENTS
Definition: acado_types.hpp:839
QProblemB::computeGivens
void computeGivens(real_t xold, real_t yold, real_t &xnew, real_t &ynew, real_t &c, real_t &s) const
EPS
BEGIN_NAMESPACE_ACADO const double EPS
Definition: acado_constants.hpp:43
RET_NO_DIAGONAL_AVAILABLE
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:201
RET_MAX_NWSR_REACHED
Definition: acado_types.hpp:1098
THROWINFO
#define THROWINFO(retval)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:397
QProblemB::obtainAuxiliaryWorkingSet
returnValue obtainAuxiliaryWorkingSet(const real_t *const xOpt, const real_t *const yOpt, const Bounds *const guessedBounds, Bounds *auxiliaryBounds) const
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:900
QProblemB::shallRefactorise
BooleanType shallRefactorise(const Bounds *const guessedBounds) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:3369
MessageHandling::throwInfo
returnValue throwInfo(returnValue Inumber, const char *additionaltext, const char *functionname, const char *filename, const unsigned long linenumber, VisibilityStatus localVisibilityStatus)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/MessageHandling.cpp:315
returnValue
Allows to pass back messages to the calling function.
Definition: acado_types.hpp:1139
RET_INDEXLIST_CORRUPTED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:83
HST_POSDEF
#define HST_POSDEF
Definition: qpOASES_wrapper.h:66
getCPUtime
real_t getCPUtime()
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/Utils.cpp:433
QProblemB::solveRegularisedQP
returnValue solveRegularisedQP(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new, int &nWSR, real_t *const cputime, int nWSRperformed=0)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2579
SubjectTo
Base class for managing working sets of bounds and constraints.
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/SubjectTo.hpp:50
QProblemB::determineStepDirection
returnValue determineStepDirection(const real_t *const delta_g, const real_t *const delta_lb, const real_t *const delta_ub, BooleanType Delta_bB_isZero, real_t *const delta_xFX, real_t *const delta_xFR, real_t *const delta_yFX)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2958
ZERO
const double ZERO
Definition: acado_constants.hpp:62
Options::print
returnValue print() const
Definition: qpOASES-3.0beta/src/Options.cpp:256
QProblemB::setupAuxiliaryWorkingSet
returnValue setupAuxiliaryWorkingSet(const Bounds *const auxiliaryBounds, BooleanType setupAfresh)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:1029
RET_INIT_FAILED_REGULARISATION
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:102
RR
#define RR(I, J)
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:132
PL_MEDIUM
#define PL_MEDIUM
Definition: qpOASES_wrapper.h:85
QProblemB::getNV
int getNV() const
QProblemB::setH
returnValue setH(const real_t *const H_new)
getMin
real_t getMin(real_t x, real_t y)
THROWERROR
#define THROWERROR(retval)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:391
VS_HIDDEN
Definition: acado_message_handling.hpp:54
QProblemB::setInfeasibilityFlag
returnValue setInfeasibilityFlag(returnValue returnvalue)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2006
RET_INIT_FAILED_UNBOUNDEDNESS
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:105
QProblemB::setupAuxiliaryQPsolution
returnValue setupAuxiliaryQPsolution(const real_t *const xOpt, const real_t *const yOpt)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:1096
QProblemB::backsolveR
returnValue backsolveR(const real_t *const b, BooleanType transposed, real_t *const a)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:1356
QProblemB::init
returnValue init(const real_t *const _H, const real_t *const _g, const real_t *const _lb, const real_t *const _ub, int &nWSR, const real_t *const yOpt=0, real_t *const cputime=0)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:268
HessianType
HessianType
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:116
QPS_SOLVED
Definition: acado_types.hpp:674
Flipper::set
returnValue set(const Bounds *const _bounds, const double *const _R, const Constraints *const _constraints=0, const double *const _Q=0, const double *const _T=0)
Definition: qpOASES-3.0beta/src/Flipper.cpp:176
__FUNCTION__
#define __FUNCTION__
Definition: acado_message_handling.hpp:90
MessageHandling::throwWarning
returnValue throwWarning(returnValue Wnumber, const char *additionaltext, const char *functionname, const char *filename, const unsigned long linenumber, VisibilityStatus localVisibilityStatus)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/MessageHandling.cpp:291
HST_IDENTITY
#define HST_IDENTITY
Definition: qpOASES_wrapper.h:65
RET_ADD_TO_ACTIVESET_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:150
PL_NONE
#define PL_NONE
Definition: qpOASES_wrapper.h:83
readFromFile
returnValue readFromFile(real_t *data, int nrow, int ncol, const char *datafilename)
Definition: external_packages/qpOASES-3.0beta/src/Utils.cpp:283
myPrintf
returnValue myPrintf(const char *s)
Definition: external_packages/qpOASES-3.0beta/src/Utils.cpp:242
ST_UNDEFINED
#define ST_UNDEFINED
Definition: qpOASES_wrapper.h:78
RET_CANNOT_REGULARISE_IDENTITY
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:176
RET_HESSIAN_INDEFINITE
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:162
RET_HESSIAN_ALREADY_REGULARISED
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:175
ST_INACTIVE
#define ST_INACTIVE
Definition: qpOASES_wrapper.h:74
RET_STEPDIRECTION_DETERMINATION_FAILED
Definition: acado_types.hpp:1094
Bounds::setupBound
returnValue setupBound(int _number, SubjectToStatus _status)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/Bounds.cpp:120
RET_SETUP_WORKINGSET_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:108
QProblemB::changeActiveSet
returnValue changeActiveSet(int BC_idx, SubjectToStatus BC_status)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:3270
RET_QP_NOT_SOLVED
Definition: acado_types.hpp:1083
PL_HIGH
#define PL_HIGH
Definition: qpOASES_wrapper.h:86
RET_HOTSTART_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:117
ST_UNBOUNDED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:89
QProblemB::usingRegularisation
BooleanType usingRegularisation() const
QProblemB::relativeHomotopyLength
real_t relativeHomotopyLength(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2129
RET_CHOLESKY_OF_ZERO_HESSIAN
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:180
ST_BOUNDED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:90
QProblemB::isCPUtimeLimitExceeded
BooleanType isCPUtimeLimitExceeded(const real_t *const cputime, real_t starttime, int nWSR) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2021
SubjectTo::getStatus
SubjectToStatus getStatus(int i) const
QProblemB::hotstart
returnValue hotstart(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new, int &nWSR, real_t *const cputime)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:285
QProblemB::performRatioTest
returnValue performRatioTest(int nIdx, const int *const idxList, const SubjectTo *const subjectTo, const real_t *const num, const real_t *const den, real_t epsNum, real_t epsDen, real_t &t, int &BC_idx) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2081
QProblemB::isBlocking
BooleanType isBlocking(real_t num, real_t den, real_t epsNum, real_t epsDen, real_t &t) const
QProblemB::determineDataShift
returnValue determineDataShift(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new, real_t *const delta_g, real_t *const delta_lb, real_t *const delta_ub, BooleanType &Delta_bB_isZero)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1724
QProblemB::addBound
returnValue addBound(int number, SubjectToStatus B_status, BooleanType updateCholesky)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:1226
BooleanType
BooleanType
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:44
QPS_AUXILIARYQPSOLVED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:77
RET_CANNOT_REGULARISE_SPARSE
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:177
QProblemB::setupQPdata
returnValue setupQPdata(const real_t *const _H, const real_t *const _g, const real_t *const _lb, const real_t *const _ub)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:1505
nu
const int nu
Definition: rocket_with_templates_ocp.hpp:38
MessageHandling::setInfoVisibilityStatus
void setInfoVisibilityStatus(VisibilityStatus _infoVisibility)
QProblemB::printIteration
returnValue printIteration(int iteration, int BC_idx, SubjectToStatus BC_status)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:3588
QProblemB::performStep
returnValue performStep(const real_t *const delta_g, const real_t *const delta_lb, const real_t *const delta_ub, const real_t *const delta_xFX, const real_t *const delta_xFR, const real_t *const delta_yFX, int &BC_idx, SubjectToStatus &BC_status)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:3125
RET_HOTSTART_STOPPED_UNBOUNDEDNESS
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:127
RET_INIT_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:100
RET_PRINTLEVEL_CHANGED
Definition: acado_types.hpp:846
printCopyrightNotice
returnValue printCopyrightNotice()
Definition: external_packages/qpOASES-3.0beta/src/Utils.cpp:265
QProblemB::getNFR
int getNFR()
QPS_PERFORMINGHOMOTOPY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:79
QProblemB::setupQPdataFromFile
returnValue setupQPdataFromFile(const char *const H_file, const char *const g_file, const char *const lb_file, const char *const ub_file)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1873
RET_ADDBOUND_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:143
RET_NO_REGSTEP_NWSR
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:178
SubjectToStatus
SubjectToStatus
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:97
RET_INIT_FAILED_CHOLESKY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:102
PL_LOW
#define PL_LOW
Definition: qpOASES_wrapper.h:84
VS_VISIBLE
Definition: acado_message_handling.hpp:53
rhs
void rhs(const real_t *x, real_t *f)
Definition: examples/code_generation/simulation/externModel_export/model.c:5
END_NAMESPACE_QPOASES
#define END_NAMESPACE_QPOASES
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:120
PrintLevel
PrintLevel
Definition: acado_types.hpp:292
RET_QP_ALREADY_INITIALISED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:97
RET_USING_REGULARISATION
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:172
SubjectTo::getType
SubjectToType getType(int i) const
HST_POSDEF_NULLSPACE
#define HST_POSDEF_NULLSPACE
Definition: qpOASES_wrapper.h:67
RET_ADD_TO_ACTIVESET
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:148
QProblemB::setupAuxiliaryQP
returnValue setupAuxiliaryQP(const Bounds *const guessedBounds)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2896
QProblemB::isInitialised
BooleanType isInitialised() const
QProblemB::getStatus
QProblemStatus getStatus() const
BT_TRUE
#define BT_TRUE
Definition: acado_types.hpp:47
RET_OPTIMAL_SOLUTION_FOUND
Definition: acado_types.hpp:1097
RET_SETUP_AUXILIARYQP_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:109
RET_REMOVE_FROM_ACTIVESET
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:147
RET_INIT_SUCCESSFUL
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:106
ST_UNKNOWN
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:92
HST_ZERO
#define HST_ZERO
Definition: qpOASES_wrapper.h:64
QProblemB::setUB
returnValue setUB(const real_t *const ub_new)
THROWWARNING
#define THROWWARNING(retval)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:394
MessageHandling::setErrorVisibilityStatus
void setErrorVisibilityStatus(VisibilityStatus _errorVisibility)
RET_OBTAINING_WORKINGSET_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:107
HST_SEMIDEF
#define HST_SEMIDEF
Definition: qpOASES_wrapper.h:68
Bounds::getNFX
int getNFX()
RET_STEPSIZE
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:138
QPS_NOTINITIALISED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:74
QProblemB::setG
returnValue setG(const real_t *const g_new)
QPS_HOMOTOPYQPSOLVED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:81
HST_UNKNOWN
#define HST_UNKNOWN
Definition: qpOASES_wrapper.h:70
QProblemB::removeBound
returnValue removeBound(int number, BooleanType updateCholesky)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:1287
SubjectTo::setNoUpper
void setNoUpper(BooleanType _status)
QPS_PREPARINGAUXILIARYQP
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:75
QProblemB::solveQP
returnValue solveQP(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new, int &nWSR, real_t *const cputime, int nWSRperformed=0)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2352
RET_INIT_FAILED_INFEASIBILITY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:104
__FILE__
#define __FILE__
Definition: acado_message_handling.hpp:94
BT_FALSE
#define BT_FALSE
Definition: acado_types.hpp:49
PL_TABULAR
#define PL_TABULAR
Definition: qpOASES_wrapper.h:82
ST_EQUALITY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:91
MessageHandling::setWarningVisibilityStatus
void setWarningVisibilityStatus(VisibilityStatus _warningVisibility)
SubjectTo::hasNoUpper
BooleanType hasNoUpper() const
Bounds
Manages working sets of bounds (= box constraints).
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Bounds.hpp:50
SubjectTo::hasNoLower
BooleanType hasNoLower() const
ST_UPPER
#define ST_UPPER
Definition: qpOASES_wrapper.h:75
BEGIN_NAMESPACE_QPOASES
#define BEGIN_NAMESPACE_QPOASES
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:117
real_t
double real_t
Definition: AD_test.c:10
RET_STEPDIRECTION_FAILED_CHOLESKY
Definition: acado_types.hpp:1095
QProblemB::applyGivens
void applyGivens(real_t c, real_t s, real_t xold, real_t yold, real_t &xnew, real_t &ynew) const
Flipper::init
returnValue init(int _nV=0, int _nC=0)
Definition: qpOASES-3.0beta/src/Flipper.cpp:116
QProblemB::loadQPvectorsFromFile
returnValue loadQPvectorsFromFile(const char *const g_file, const char *const lb_file, const char *const ub_file, real_t *const g_new, real_t *const lb_new, real_t *const ub_new) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1945
getGlobalMessageHandler
MessageHandling * getGlobalMessageHandler()
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/MessageHandling.cpp:515
QProblemB::solveInitialQP
returnValue solveInitialQP(const real_t *const xOpt, const real_t *const yOpt, const Bounds *const guessedBounds, int &nWSR, real_t *const cputime)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/QProblemB.cpp:791
RET_DIV_BY_ZERO
Definition: acado_types.hpp:837
RET_REMOVE_FROM_ACTIVESET_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:149
Bounds::getFixed
Indexlist * getFixed()
QProblemB::setLB
returnValue setLB(const real_t *const lb_new)
Flipper::get
returnValue get(Bounds *const _bounds, double *const R, Constraints *const _constraints=0, double *const _Q=0, double *const _T=0) const
Definition: qpOASES-3.0beta/src/Flipper.cpp:147
Bounds::moveFreeToFixed
returnValue moveFreeToFixed(int _number, SubjectToStatus _status)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/Bounds.cpp:217
RET_FEWER_REGSTEPS_NWSR
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:179
RET_HOTSTART_STOPPED_INFEASIBILITY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:126
SymmetricMatrix
Abstract base class for interfacing matrix-vector operations tailored to symmetric matrices...
Definition: qpOASES-3.0beta/include/qpOASES/Matrices.hpp:246

acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:02