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