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

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