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