external_packages/qpOASES-3.0beta/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-2011 by Hans Joachim Ferreau, Andreas Potschka,
6  * Christian Kirches et al. All rights reserved.
7  *
8  * qpOASES is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2.1 of the License, or (at your option) any later version.
12  *
13  * qpOASES is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  * See the GNU Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with qpOASES; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  *
22  */
23 
24 
36 #include <qpOASES/QProblem.hpp>
37 
38 
39 BEGIN_NAMESPACE_QPOASES
40 
41 
42 /*****************************************************************************
43  * P U B L I C *
44  *****************************************************************************/
45 
46 
47 /*
48  * Q P r o b l e m
49  */
50 QProblem::QProblem( ) : QProblemB( )
51 {
52  freeConstraintMatrix = BT_FALSE;
53  A = 0;
54 
55  lbA = 0;
56  ubA = 0;
57 
58  sizeT = 0;
59  T = 0;
60  Q = 0;
61 
62  Ax = 0;
63  Ax_l = 0;
64  Ax_u = 0;
65 
66  constraintProduct = 0;
67 
68  tempA = 0;
69  ZFR_delta_xFRz = 0;
70  delta_xFRy = 0;
71  delta_xFRz = 0;
72  tempB = 0;
73  delta_yAC_TMP = 0;
74 }
75 
76 
77 /*
78  * Q P r o b l e m
79  */
80 QProblem::QProblem( int _nV, int _nC, HessianType _hessianType ) : QProblemB( _nV,_hessianType )
81 {
82  /* consistency checks */
83  if ( _nV <= 0 )
84  _nV = 1;
85 
86  if ( _nC < 0 )
87  {
88  _nC = 0;
89  THROWERROR( RET_INVALID_ARGUMENTS );
90  }
91 
92  if ( _nC > 0 )
93  {
94  freeConstraintMatrix = BT_FALSE;
95  A = 0;
96 
97  lbA = new real_t[_nC];
98  ubA = new real_t[_nC];
99  }
100  else
101  {
102  /* prevent segmentation faults in case nC == 0
103  * (avoiding checks for A!=0 around all calls to A->... */
104  freeConstraintMatrix = BT_TRUE;
105  A = new DenseMatrix( );
106 
107  lbA = 0;
108  ubA = 0;
109  }
110 
111  constraints.init( _nC );
112 
113  delete[] y; /* y of no constraints version too short! */
114  y = new real_t[_nV+_nC];
115 
116  sizeT = (int)getMin( _nC,_nV );
117  T = new real_t[sizeT*sizeT];
118  Q = new real_t[_nV*_nV];
119 
120  if ( _nC > 0 )
121  {
122  Ax = new real_t[_nC];
123  Ax_l = new real_t[_nC];
124  Ax_u = new real_t[_nC];
125  }
126  else
127  {
128  Ax = 0;
129  Ax_l = 0;
130  Ax_u = 0;
131  }
132 
133  constraintProduct = 0;
134 
135  tempA = new real_t[_nV]; /* nFR */
136  ZFR_delta_xFRz = new real_t[_nV]; /* nFR */
137  delta_xFRz = new real_t[_nV]; /* nZ */
138 
139  if ( _nC > 0 )
140  {
141  tempB = new real_t[_nC]; /* nAC */
142  delta_xFRy = new real_t[_nC]; /* nAC */
143  delta_yAC_TMP = new real_t[_nC]; /* nAC */
144  }
145  else
146  {
147  tempB = 0;
148  delta_xFRy = 0;
149  delta_yAC_TMP = 0;
150  }
151 
152  flipper.init( _nV,_nC );
153 }
154 
155 
156 /*
157  * Q P r o b l e m
158  */
159 QProblem::QProblem( const QProblem& rhs ) : QProblemB( rhs )
160 {
161  freeConstraintMatrix = BT_FALSE;
162  A = 0;
163 
164  copy( rhs );
165 }
166 
167 
168 /*
169  * ~ Q P r o b l e m
170  */
171 QProblem::~QProblem( )
172 {
173  clear( );
174 }
175 
176 
177 /*
178  * o p e r a t o r =
179  */
180 QProblem& QProblem::operator=( const QProblem& rhs )
181 {
182  if ( this != &rhs )
183  {
184  clear( );
185  QProblemB::operator=( rhs );
186  copy( rhs );
187  }
188 
189  return *this;
190 }
191 
192 
193 /*
194  * r e s e t
195  */
196 returnValue QProblem::reset( )
197 {
198  int i;
199  int nV = getNV( );
200  int nC = getNC( );
201 
202  if ( nV == 0 )
203  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
204 
205 
206  /* 1) Reset bounds, Cholesky decomposition and status flags. */
207  if ( QProblemB::reset( ) != SUCCESSFUL_RETURN )
208  return THROWERROR( RET_RESET_FAILED );
209 
210  /* 2) Reset constraints. */
211  constraints.init( nC );
212 
213  /* 3) Reset TQ factorisation. */
214  for( i=0; i<sizeT*sizeT; ++i )
215  T[i] = 0.0;
216 
217  for( i=0; i<nV*nV; ++i )
218  Q[i] = 0.0;
219 
220  /* 4) Reset constraint product pointer. */
221  constraintProduct = 0;
222 
223  return SUCCESSFUL_RETURN;
224 }
225 
226 
227 /*
228  * i n i t
229  */
230 returnValue QProblem::init( SymmetricMatrix *_H, const real_t* const _g, Matrix *_A,
231  const real_t* const _lb, const real_t* const _ub,
232  const real_t* const _lbA, const real_t* const _ubA,
233  int& nWSR, real_t* const cputime
234  )
235 {
236  if ( getNV( ) == 0 )
237  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
238 
239  /* 1) Consistency check. */
240  if ( isInitialised( ) == BT_TRUE )
241  {
242  THROWWARNING( RET_QP_ALREADY_INITIALISED );
243  reset( );
244  }
245 
246  /* 2) Setup QP data. */
247  if ( setupQPdata( _H,_g,_A,_lb,_ub,_lbA,_ubA ) != SUCCESSFUL_RETURN )
248  return THROWERROR( RET_INVALID_ARGUMENTS );
249 
250  /* 3) Call to main initialisation routine (without any additional information). */
251  return solveInitialQP( 0,0,0,0, nWSR,cputime );
252 }
253 
254 
255 /*
256  * i n i t
257  */
258 returnValue QProblem::init( const real_t* const _H, const real_t* const _g, const real_t* const _A,
259  const real_t* const _lb, const real_t* const _ub,
260  const real_t* const _lbA, const real_t* const _ubA,
261  int& nWSR, real_t* const cputime
262  )
263 {
264  if ( getNV( ) == 0 )
265  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
266 
267  /* 1) Consistency check. */
268  if ( isInitialised( ) == BT_TRUE )
269  {
270  THROWWARNING( RET_QP_ALREADY_INITIALISED );
271  reset( );
272  }
273 
274  /* 2) Setup QP data. */
275  if ( setupQPdata( _H,_g,_A,_lb,_ub,_lbA,_ubA ) != SUCCESSFUL_RETURN )
276  return THROWERROR( RET_INVALID_ARGUMENTS );
277 
278  /* 3) Call to main initialisation routine (without any additional information). */
279  return solveInitialQP( 0,0,0,0, nWSR,cputime );
280 }
281 
282 
283 /*
284  * i n i t
285  */
286 returnValue QProblem::init( const char* const H_file, const char* const g_file, const char* const A_file,
287  const char* const lb_file, const char* const ub_file,
288  const char* const lbA_file, const char* const ubA_file,
289  int& nWSR, real_t* const cputime
290  )
291 {
292  if ( getNV( ) == 0 )
293  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
294 
295  /* 1) Consistency check. */
296  if ( isInitialised( ) == BT_TRUE )
297  {
298  THROWWARNING( RET_QP_ALREADY_INITIALISED );
299  reset( );
300  }
301 
302  /* 2) Setup QP data from files. */
303  if ( setupQPdataFromFile( H_file,g_file,A_file,lb_file,ub_file,lbA_file,ubA_file ) != SUCCESSFUL_RETURN )
304  return THROWERROR( RET_UNABLE_TO_READ_FILE );
305 
306  /* 3) Call to main initialisation routine (without any additional information). */
307  return solveInitialQP( 0,0,0,0, nWSR,cputime );
308 }
309 
310 
311 /*
312  * i n i t
313  */
314 returnValue QProblem::init( SymmetricMatrix *_H, const real_t* const _g, Matrix *_A,
315  const real_t* const _lb, const real_t* const _ub,
316  const real_t* const _lbA, const real_t* const _ubA,
317  int& nWSR, real_t* const cputime,
318  const real_t* const xOpt, const real_t* const yOpt,
319  const Bounds* const guessedBounds, const Constraints* const guessedConstraints
320  )
321 {
322  int i;
323  int nV = getNV( );
324  int nC = getNC( );
325 
326  if ( nV == 0 )
327  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
328 
329  /* 1) Consistency checks. */
330  if ( isInitialised( ) == BT_TRUE )
331  {
332  THROWWARNING( RET_QP_ALREADY_INITIALISED );
333  reset( );
334  }
335 
336  if ( guessedBounds != 0 )
337  {
338  for( i=0; i<nV; ++i )
339  {
340  if ( guessedBounds->getStatus( i ) == ST_UNDEFINED )
341  return THROWERROR( RET_INVALID_ARGUMENTS );
342  }
343  }
344 
345  if ( guessedConstraints != 0 )
346  {
347  for( i=0; i<nC; ++i )
348  if ( guessedConstraints->getStatus( i ) == ST_UNDEFINED )
349  return THROWERROR( RET_INVALID_ARGUMENTS );
350  }
351 
352  /* exclude these possibilities in order to avoid inconsistencies */
353  if ( ( xOpt == 0 ) && ( yOpt != 0 ) && ( ( guessedBounds != 0 ) || ( guessedConstraints != 0 ) ) )
354  return THROWERROR( RET_INVALID_ARGUMENTS );
355 
356  /* 2) Setup QP data. */
357  if ( setupQPdata( _H,_g,_A,_lb,_ub,_lbA,_ubA ) != SUCCESSFUL_RETURN )
358  return THROWERROR( RET_INVALID_ARGUMENTS );
359 
360  /* 3) Call to main initialisation routine. */
361  return solveInitialQP( xOpt,yOpt,guessedBounds,guessedConstraints, nWSR,cputime );
362 }
363 
364 
365 /*
366  * i n i t
367  */
368 returnValue QProblem::init( const real_t* const _H, const real_t* const _g, const real_t* const _A,
369  const real_t* const _lb, const real_t* const _ub,
370  const real_t* const _lbA, const real_t* const _ubA,
371  int& nWSR, real_t* const cputime,
372  const real_t* const xOpt, const real_t* const yOpt,
373  const Bounds* const guessedBounds, const Constraints* const guessedConstraints
374  )
375 {
376  int i;
377  int nV = getNV( );
378  int nC = getNC( );
379 
380  if ( nV == 0 )
381  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
382 
383  /* 1) Consistency checks. */
384  if ( isInitialised( ) == BT_TRUE )
385  {
386  THROWWARNING( RET_QP_ALREADY_INITIALISED );
387  reset( );
388  }
389 
390  if ( guessedBounds != 0 )
391  {
392  for( i=0; i<nV; ++i )
393  {
394  if ( guessedBounds->getStatus( i ) == ST_UNDEFINED )
395  return THROWERROR( RET_INVALID_ARGUMENTS );
396  }
397  }
398 
399  if ( guessedConstraints != 0 )
400  {
401  for( i=0; i<nC; ++i )
402  if ( guessedConstraints->getStatus( i ) == ST_UNDEFINED )
403  return THROWERROR( RET_INVALID_ARGUMENTS );
404  }
405 
406  /* exclude these possibilities in order to avoid inconsistencies */
407  if ( ( xOpt == 0 ) && ( yOpt != 0 ) && ( ( guessedBounds != 0 ) || ( guessedConstraints != 0 ) ) )
408  return THROWERROR( RET_INVALID_ARGUMENTS );
409 
410  /* 2) Setup QP data. */
411  if ( setupQPdata( _H,_g,_A,_lb,_ub,_lbA,_ubA ) != SUCCESSFUL_RETURN )
412  return THROWERROR( RET_INVALID_ARGUMENTS );
413 
414  /* 3) Call to main initialisation routine. */
415  return solveInitialQP( xOpt,yOpt,guessedBounds,guessedConstraints, nWSR,cputime );
416 }
417 
418 
419 /*
420  * i n i t
421  */
422 returnValue QProblem::init( const char* const H_file, const char* const g_file, const char* const A_file,
423  const char* const lb_file, const char* const ub_file,
424  const char* const lbA_file, const char* const ubA_file,
425  int& nWSR, real_t* const cputime,
426  const real_t* const xOpt, const real_t* const yOpt,
427  const Bounds* const guessedBounds, const Constraints* const guessedConstraints
428  )
429 {
430  int i;
431  int nV = getNV( );
432  int nC = getNC( );
433 
434  if ( nV == 0 )
435  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
436 
437  /* 1) Consistency checks. */
438  if ( isInitialised( ) == BT_TRUE )
439  {
440  THROWWARNING( RET_QP_ALREADY_INITIALISED );
441  reset( );
442  }
443 
444  for( i=0; i<nV; ++i )
445  {
446  if ( guessedBounds->getStatus( i ) == ST_UNDEFINED )
447  return THROWERROR( RET_INVALID_ARGUMENTS );
448  }
449 
450  for( i=0; i<nC; ++i )
451  if ( guessedConstraints->getStatus( i ) == ST_UNDEFINED )
452  return THROWERROR( RET_INVALID_ARGUMENTS );
453 
454  /* exclude these possibilities in order to avoid inconsistencies */
455  if ( ( xOpt == 0 ) && ( yOpt != 0 ) && ( ( guessedBounds != 0 ) || ( guessedConstraints != 0 ) ) )
456  return THROWERROR( RET_INVALID_ARGUMENTS );
457 
458  /* 2) Setup QP data from files. */
459  if ( setupQPdataFromFile( H_file,g_file,A_file,lb_file,ub_file,lbA_file,ubA_file ) != SUCCESSFUL_RETURN )
460  return THROWERROR( RET_UNABLE_TO_READ_FILE );
461 
462  /* 3) Call to main initialisation routine. */
463  return solveInitialQP( xOpt,yOpt,guessedBounds,guessedConstraints, nWSR,cputime );
464 }
465 
466 returnValue QProblem::computeInitialCholesky()
467 {
468  returnValue returnvalueCholesky;
469 
470  if ( ( getNAC( ) + getNFX( ) ) == 0 )
471  {
472  /* Factorise full Hessian if no bounds/constraints are active. */
473  returnvalueCholesky = setupCholeskyDecomposition( );
474 
475  /* If Hessian is not positive definite, regularise and try again. */
476  if ( returnvalueCholesky == RET_HESSIAN_NOT_SPD )
477  {
478  if ( regulariseHessian( ) != SUCCESSFUL_RETURN )
479  return THROWERROR( RET_INIT_FAILED_REGULARISATION );
480 
481  if ( setupCholeskyDecomposition( ) != SUCCESSFUL_RETURN )
482  return THROWERROR( RET_INIT_FAILED_CHOLESKY );
483  }
484 
485  if ( returnvalueCholesky == RET_INDEXLIST_CORRUPTED )
486  return THROWERROR( RET_INIT_FAILED_CHOLESKY );
487  }
488  else
489  {
490  /* Factorise projected Hessian if there active bounds/constraints. */
491  returnvalueCholesky = setupCholeskyDecompositionProjected( );
492 
493  /* If Hessian is not positive definite, regularise and try again. */
494  if ( returnvalueCholesky == RET_HESSIAN_NOT_SPD )
495  {
496  if ( regulariseHessian( ) != SUCCESSFUL_RETURN )
497  return THROWERROR( RET_INIT_FAILED_REGULARISATION );
498 
499  if ( setupCholeskyDecompositionProjected( ) != SUCCESSFUL_RETURN )
500  return THROWERROR( RET_INIT_FAILED_CHOLESKY );
501  }
502 
503  if ( returnvalueCholesky == RET_INDEXLIST_CORRUPTED )
504  return THROWERROR( RET_INIT_FAILED_CHOLESKY );
505  }
506 
507  haveCholesky = BT_TRUE;
508  return SUCCESSFUL_RETURN;
509 }
510 
511 
512 /*
513  * h o t s t a r t
514  */
515 returnValue QProblem::hotstart( const real_t* const g_new,
516  const real_t* const lb_new, const real_t* const ub_new,
517  const real_t* const lbA_new, const real_t* const ubA_new,
518  int& nWSR, real_t* const cputime
519  )
520 {
521  if ( getNV( ) == 0 )
522  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
523 
524  ++count;
525 
526  returnValue returnvalue = SUCCESSFUL_RETURN;
527  int nFar, i;
528  int nV = getNV ();
529  int nC = getNC ();
530 
531  int nWSR_max = nWSR;
532  int nWSR_performed = 0;
533 
534  real_t cputime_remaining = INFTY;
535  real_t cputime_needed = 0.0;
536 
537  real_t farbound = options.initialFarBounds;
538  real_t t, rampVal;
539 
540  if ( options.enableFarBounds == BT_FALSE )
541  {
542  if (haveCholesky == BT_FALSE)
543  {
544  returnvalue = computeInitialCholesky();
545  if (returnvalue != SUCCESSFUL_RETURN)
546  return THROWERROR(returnvalue);
547  }
548 
549  if ( usingRegularisation( ) == BT_TRUE )
550  returnvalue = solveRegularisedQP( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,cputime,0 );
551  else
552  returnvalue = solveQP( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,cputime,0 );
553  }
554  else
555  {
556  real_t *ub_new_far = new real_t[nV];
557  real_t *lb_new_far = new real_t[nV];
558  real_t *ubA_new_far = new real_t[nC];
559  real_t *lbA_new_far = new real_t[nC];
560 
561  /* possibly extend initial far bounds to largest bound/constraint data */
562  if (ub_new)
563  for (i = 0; i < nV; i++)
564  if (ub_new[i] < INFTY && ub_new[i] > farbound) farbound = ub_new[i];
565  if (lb_new)
566  for (i = 0; i < nV; i++)
567  if (lb_new[i] > -INFTY && lb_new[i] < -farbound) farbound = -lb_new[i];
568  if (ubA_new)
569  for (i = 0; i < nC; i++)
570  if (ubA_new[i] < INFTY && ubA_new[i] > farbound) farbound = ubA_new[i];
571  if (lbA_new)
572  for (i = 0; i < nC; i++)
573  if (lbA_new[i] > -INFTY && lbA_new[i] < -farbound) farbound = -lbA_new[i];
574 
575  if ( options.enableRamping == BT_TRUE )
576  {
577  /* TODO: encapsule this part to avoid code duplication! */
578  for ( i=0; i<nV; ++i )
579  {
580  t = static_cast<real_t>(i) / static_cast<real_t>(nV+nC-1);
581  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
582 
583  if ( ( lb_new == 0 ) || ( lb_new[i] <= -rampVal ) )
584  lb_new_far[i] = - rampVal;
585  else
586  lb_new_far[i] = lb_new[i];
587  if ( ( ub_new == 0 ) || ( ub_new[i] >= rampVal ) )
588  ub_new_far[i] = rampVal;
589  else
590  ub_new_far[i] = ub_new[i];
591  }
592  for ( i=0; i<nC; ++i )
593  {
594  t = static_cast<real_t>(nV+i) / static_cast<real_t>(nV+nC-1);
595  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
596 
597  if ( ( lbA_new == 0 ) || ( lbA_new[i] <= -rampVal ) )
598  lbA_new_far[i] = - rampVal;
599  else
600  lbA_new_far[i] = lbA_new[i];
601  if ( ( ubA_new == 0 ) || ( ubA_new[i] >= rampVal ) )
602  ubA_new_far[i] = rampVal;
603  else
604  ubA_new_far[i] = ubA_new[i];
605  }
606  }
607  else
608  {
609  for ( i=0; i<nV; ++i )
610  {
611  lb_new_far[i] = lb_new[i];
612  ub_new_far[i] = ub_new[i];
613  }
614 
615  for ( i=0; i<nC; ++i )
616  {
617  lbA_new_far[i] = lbA_new[i];
618  ubA_new_far[i] = ubA_new[i];
619  }
620  }
621 
622  if (haveCholesky == BT_FALSE)
623  {
624  returnvalue = computeInitialCholesky();
625  if (returnvalue != SUCCESSFUL_RETURN)
626  goto farewell;
627  }
628 
629  for ( ;; )
630  {
631  nWSR = nWSR_max;
632  if ( cputime != 0 )
633  cputime_remaining = *cputime - cputime_needed;
634 
635  if ( usingRegularisation( ) == BT_TRUE )
636  returnvalue = solveRegularisedQP( g_new,lb_new_far,ub_new_far,lbA_new_far,ubA_new_far, nWSR,&cputime_remaining,nWSR_performed );
637  else
638  returnvalue = solveQP( g_new,lb_new_far,ub_new_far,lbA_new_far,ubA_new_far, nWSR,&cputime_remaining,nWSR_performed );
639 
640  nWSR_performed = nWSR;
641  cputime_needed += cputime_remaining;
642 
643  /* Check for active far-bounds and move them away */
644  nFar = 0;
645  farbound *= options.growFarBounds;
646 
647  real_t maxFarbound = 1e20;
648  if ( infeasible == BT_TRUE )
649  {
650  if ( farbound > maxFarbound )
651  {
652  returnvalue = RET_HOTSTART_STOPPED_INFEASIBILITY;
653  goto farewell;
654  }
655 
656  if ( options.enableRamping == BT_TRUE )
657  {
658  /* TODO: encapsule this part to avoid code duplication! */
659  for ( i=0; i<nV; ++i )
660  {
661  t = static_cast<real_t>(i) / static_cast<real_t>(nV+nC-1);
662  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
663 
664  if ( lb_new == 0 )
665  lb_new_far[i] = - rampVal;
666  else
667  lb_new_far[i] = getMax (- rampVal, lb_new[i]);
668 
669  if ( ub_new == 0 )
670  ub_new_far[i] = rampVal;
671  else
672  ub_new_far[i] = getMin (rampVal, ub_new[i]);
673  }
674  for ( i=0; i<nC; ++i )
675  {
676  t = static_cast<real_t>(nV+i) / static_cast<real_t>(nV+nC-1);
677  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
678 
679  if ( lbA_new == 0 )
680  lbA_new_far[i] = - rampVal;
681  else
682  lbA_new_far[i] = getMax (- rampVal, lbA_new[i]);
683 
684  if ( ubA_new == 0 )
685  ubA_new_far[i] = rampVal;
686  else
687  ubA_new_far[i] = getMin (rampVal, ubA_new[i]);
688  }
689  }
690  else
691  {
692  for ( i=0; i<nV; ++i )
693  {
694  lb_new_far[i] = lb_new[i];
695  ub_new_far[i] = ub_new[i];
696  }
697 
698  for ( i=0; i<nC; ++i )
699  {
700  lbA_new_far[i] = lbA_new[i];
701  ubA_new_far[i] = ubA_new[i];
702  }
703  }
704  }
705  else if ( status == QPS_SOLVED )
706  {
707  real_t tol = farbound * options.boundTolerance;
708  status = QPS_HOMOTOPYQPSOLVED;
709 
710  nFar = 0;
711  for ( i=0; i<nV; ++i )
712  {
713  if ( ( ( lb_new == 0 ) || ( lb_new_far[i] > lb_new[i] ) ) && ( fabs ( lb_new_far[i] - x[i] ) < tol ) )
714  ++nFar;
715  if ( ( ( ub_new == 0 ) || ( ub_new_far[i] < ub_new[i] ) ) && ( fabs ( ub_new_far[i] - x[i] ) < tol ) )
716  ++nFar;
717  }
718  for ( i=0; i<nC; ++i )
719  {
720  if ( ( ( lbA_new == 0 ) || ( lbA_new_far[i] > lbA_new[i] ) ) && ( fabs ( lbA_new_far[i] - Ax[i] ) < tol ) )
721  ++nFar;
722  if ( ( ( ubA_new == 0 ) || ( ubA_new_far[i] < ubA_new[i] ) ) && ( fabs ( ubA_new_far[i] - Ax[i] ) < tol ) )
723  ++nFar;
724  }
725 
726  if ( nFar == 0 )
727  break;
728 
729  if ( farbound > maxFarbound )
730  {
731  unbounded = BT_TRUE;
732  returnvalue = RET_HOTSTART_STOPPED_UNBOUNDEDNESS;
733  goto farewell;
734  }
735 
736  if ( options.enableRamping == BT_TRUE )
737  {
738  /* TODO: encapsule this part to avoid code duplication! */
739  for ( i=0; i<nV; ++i )
740  {
741  t = static_cast<real_t>(i) / static_cast<real_t>(nV+nC-1);
742  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
743 
744  if ( lb_new == 0 )
745  lb_new_far[i] = - rampVal;
746  else
747  lb_new_far[i] = getMax (- rampVal, lb_new[i]);
748 
749  if ( ub_new == 0 )
750  ub_new_far[i] = rampVal;
751  else
752  ub_new_far[i] = getMin (rampVal, ub_new[i]);
753  }
754  for ( i=0; i<nC; ++i )
755  {
756  t = static_cast<real_t>(nV+i) / static_cast<real_t>(nV+nC-1);
757  rampVal = farbound + (1.0-t) * ramp0 + t * ramp1;
758 
759  if ( lbA_new == 0 )
760  lbA_new_far[i] = - rampVal;
761  else
762  lbA_new_far[i] = getMax (- rampVal, lbA_new[i]);
763 
764  if ( ubA_new == 0 )
765  ubA_new_far[i] = rampVal;
766  else
767  ubA_new_far[i] = getMin (rampVal, ubA_new[i]);
768  }
769  }
770  else
771  {
772  for ( i=0; i<nV; ++i )
773  {
774  lb_new_far[i] = lb_new[i];
775  ub_new_far[i] = ub_new[i];
776  }
777 
778  for ( i=0; i<nC; ++i )
779  {
780  lbA_new_far[i] = lbA_new[i];
781  ubA_new_far[i] = ubA_new[i];
782  }
783  }
784  }
785  else
786  {
787  /* some other error */
788  break;
789  }
790 
791  /* swap ramp0 and ramp1 */
792  t = ramp0; ramp0 = ramp1; ramp1 = t;
793  }
794 
795  farewell:
796  if ( cputime != 0 )
797  *cputime = cputime_needed;
798  delete[] lbA_new_far; delete[] ubA_new_far;
799  delete[] lb_new_far; delete[] ub_new_far;
800  }
801 
802  return ( returnvalue != SUCCESSFUL_RETURN ) ? THROWERROR( returnvalue ) : returnvalue;
803 }
804 
805 
806 /*
807  * h o t s t a r t
808  */
809 returnValue QProblem::hotstart( const char* const g_file,
810  const char* const lb_file, const char* const ub_file,
811  const char* const lbA_file, const char* const ubA_file,
812  int& nWSR, real_t* const cputime
813  )
814 {
815  int nV = getNV( );
816  int nC = getNC( );
817 
818  if ( nV == 0 )
819  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
820 
821  /* consistency check */
822  if ( g_file == 0 )
823  return THROWERROR( RET_INVALID_ARGUMENTS );
824 
825 
826  /* 1) Allocate memory (if bounds exist). */
827  real_t* g_new = new real_t[nV];
828  real_t* lb_new = 0;
829  real_t* ub_new = 0;
830  real_t* lbA_new = 0;
831  real_t* ubA_new = 0;
832 
833  if ( lb_file != 0 )
834  lb_new = new real_t[nV];
835  if ( ub_file != 0 )
836  ub_new = new real_t[nV];
837  if ( lbA_file != 0 )
838  lbA_new = new real_t[nC];
839  if ( ubA_file != 0 )
840  ubA_new = new real_t[nC];
841 
842  /* 2) Load new QP vectors from file. */
843  returnValue returnvalue;
844  returnvalue = loadQPvectorsFromFile( g_file,lb_file,ub_file,lbA_file,ubA_file,
845  g_new,lb_new,ub_new,lbA_new,ubA_new
846  );
847  if ( returnvalue != SUCCESSFUL_RETURN )
848  {
849  if ( ubA_file != 0 )
850  delete[] ubA_new;
851  if ( lbA_file != 0 )
852  delete[] lbA_new;
853  if ( ub_file != 0 )
854  delete[] ub_new;
855  if ( lb_file != 0 )
856  delete[] lb_new;
857  delete[] g_new;
858 
859  return THROWERROR( RET_UNABLE_TO_READ_FILE );
860  }
861 
862  /* 3) Actually perform hotstart. */
863  returnvalue = hotstart( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,cputime );
864 
865  /* 4) Free memory. */
866  if ( ubA_file != 0 )
867  delete[] ubA_new;
868  if ( lbA_file != 0 )
869  delete[] lbA_new;
870  if ( ub_file != 0 )
871  delete[] ub_new;
872  if ( lb_file != 0 )
873  delete[] lb_new;
874  delete[] g_new;
875 
876  return returnvalue;
877 }
878 
879 
880 /*
881  * h o t s t a r t
882  */
883 returnValue QProblem::hotstart( const real_t* const g_new,
884  const real_t* const lb_new, const real_t* const ub_new,
885  const real_t* const lbA_new, const real_t* const ubA_new,
886  int& nWSR, real_t* const cputime,
887  const Bounds* const guessedBounds, const Constraints* const guessedConstraints
888  )
889 {
890  int nV = getNV( );
891  int nC = getNC( );
892 
893  if ( nV == 0 )
894  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
895 
896 
897  /* start runtime measurement */
898  real_t starttime = 0.0;
899  if ( cputime != 0 )
900  starttime = getCPUtime( );
901 
902 
903  /* 1) Update working sets according to guesses for working sets of bounds and constraints. */
904  if ( ( guessedBounds != 0 ) && ( guessedConstraints != 0 ) )
905  {
906  if ( setupAuxiliaryQP( guessedBounds,guessedConstraints ) != SUCCESSFUL_RETURN )
907  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
908  }
909 
910  if ( ( guessedBounds == 0 ) && ( guessedConstraints != 0 ) )
911  {
912  /* create empty bounds for setting up auxiliary QP */
913  Bounds emptyBounds( nV );
914  if ( emptyBounds.setupAllFree( ) != SUCCESSFUL_RETURN )
915  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
916 
917  if ( setupAuxiliaryQP( &emptyBounds,guessedConstraints ) != SUCCESSFUL_RETURN )
918  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
919  }
920 
921  if ( ( guessedBounds != 0 ) && ( guessedConstraints == 0 ) )
922  {
923  /* create empty constraints for setting up auxiliary QP */
924  Constraints emptyConstraints( nC );
925  if ( emptyConstraints.setupAllInactive( ) != SUCCESSFUL_RETURN )
926  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
927 
928  if ( setupAuxiliaryQP( guessedBounds,&emptyConstraints ) != SUCCESSFUL_RETURN )
929  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
930  }
931 
932  if ( ( guessedBounds == 0 ) && ( guessedConstraints == 0 ) )
933  {
934  /* create empty bounds and constraints for setting up auxiliary QP */
935  Bounds emptyBounds( nV );
936  Constraints emptyConstraints( nC );
937  if ( emptyBounds.setupAllFree( ) != SUCCESSFUL_RETURN )
938  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
939  if ( emptyConstraints.setupAllInactive( ) != SUCCESSFUL_RETURN )
940  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
941 
942  if ( setupAuxiliaryQP( &emptyBounds,&emptyConstraints ) != SUCCESSFUL_RETURN )
943  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
944  }
945 
946  status = QPS_AUXILIARYQPSOLVED;
947 
948 
949  /* 2) Perform usual homotopy. */
950 
951  /* Allow only remaining CPU time for usual hotstart. */
952  if ( cputime != 0 )
953  *cputime -= getCPUtime( ) - starttime;
954 
955  returnValue returnvalue = hotstart( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,cputime );
956 
957  /* stop runtime measurement */
958  if ( cputime != 0 )
959  *cputime = getCPUtime( ) - starttime;
960 
961  return returnvalue;
962 }
963 
964 
965 /*
966  * h o t s t a r t
967  */
968 returnValue QProblem::hotstart( const char* const g_file,
969  const char* const lb_file, const char* const ub_file,
970  const char* const lbA_file, const char* const ubA_file,
971  int& nWSR, real_t* const cputime,
972  const Bounds* const guessedBounds, const Constraints* const guessedConstraints
973  )
974 {
975  int nV = getNV( );
976  int nC = getNC( );
977 
978  if ( nV == 0 )
979  return THROWERROR( RET_QPOBJECT_NOT_SETUP );
980 
981  /* consistency check */
982  if ( g_file == 0 )
983  return THROWERROR( RET_INVALID_ARGUMENTS );
984 
985 
986  /* 1) Allocate memory (if bounds exist). */
987  real_t* g_new = new real_t[nV];
988  real_t* lb_new = 0;
989  real_t* ub_new = 0;
990  real_t* lbA_new = 0;
991  real_t* ubA_new = 0;
992 
993  if ( lb_file != 0 )
994  lb_new = new real_t[nV];
995  if ( ub_file != 0 )
996  ub_new = new real_t[nV];
997  if ( lbA_file != 0 )
998  lbA_new = new real_t[nC];
999  if ( ubA_file != 0 )
1000  ubA_new = new real_t[nC];
1001 
1002  /* 2) Load new QP vectors from file. */
1003  returnValue returnvalue;
1004  returnvalue = loadQPvectorsFromFile( g_file,lb_file,ub_file,lbA_file,ubA_file,
1005  g_new,lb_new,ub_new,lbA_new,ubA_new
1006  );
1007  if ( returnvalue != SUCCESSFUL_RETURN )
1008  {
1009  if ( ubA_file != 0 )
1010  delete[] ubA_new;
1011  if ( lbA_file != 0 )
1012  delete[] lbA_new;
1013  if ( ub_file != 0 )
1014  delete[] ub_new;
1015  if ( lb_file != 0 )
1016  delete[] lb_new;
1017  delete[] g_new;
1018 
1019  return THROWERROR( RET_UNABLE_TO_READ_FILE );
1020  }
1021 
1022  /* 3) Actually perform hotstart using initialised homotopy. */
1023  returnvalue = hotstart( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,cputime,
1024  guessedBounds,guessedConstraints
1025  );
1026 
1027  /* 4) Free memory. */
1028  if ( ubA_file != 0 )
1029  delete[] ubA_new;
1030  if ( lbA_file != 0 )
1031  delete[] lbA_new;
1032  if ( ub_file != 0 )
1033  delete[] ub_new;
1034  if ( lb_file != 0 )
1035  delete[] lb_new;
1036  delete[] g_new;
1037 
1038  return returnvalue;
1039 }
1040 
1041 
1042 /*
1043  *
1044  */
1045 returnValue QProblem::solveCurrentEQP ( const int n_rhs,
1046  const real_t* g_in,
1047  const real_t* lb_in,
1048  const real_t* ub_in,
1049  const real_t* lbA_in,
1050  const real_t* ubA_in,
1051  real_t* x_out,
1052  real_t* y_out
1053  )
1054 {
1055  returnValue returnvalue = SUCCESSFUL_RETURN;
1056  int ii, jj;
1057  int nV = getNV( );
1058  int nC = getNC( );
1059  int nFR = getNFR( );
1060  int nFX = getNFX( );
1061  int nAC = getNAC( );
1062 
1063  real_t *delta_xFX = new real_t[nFX];
1064  real_t *delta_xFR = new real_t[nFR];
1065  real_t *delta_yAC = new real_t[nAC];
1066  real_t *delta_yFX = new real_t[nFX];
1067 
1068  /* 1) Determine index arrays. */
1069  int* FR_idx;
1070  int* FX_idx;
1071  int* AC_idx;
1072 
1073  bounds.getFree( )->getNumberArray( &FR_idx );
1074  bounds.getFixed( )->getNumberArray( &FX_idx );
1075  constraints.getActive( )->getNumberArray( &AC_idx );
1076 
1077  for ( ii = 0 ; ii < n_rhs; ++ii )
1078  {
1079  returnvalue = determineStepDirection(
1080  g_in, lbA_in, ubA_in, lb_in, ub_in, BT_FALSE, BT_FALSE,
1081  delta_xFX, delta_xFR, delta_yAC, delta_yFX );
1082 
1083  for ( jj = 0; jj < nFX; ++jj )
1084  x_out[FX_idx[jj]] = delta_xFX[jj];
1085  for ( jj = 0; jj < nFR; ++jj )
1086  x_out[FR_idx[jj]] = delta_xFR[jj];
1087  for ( jj = 0; jj < nFX; ++jj )
1088  y_out[FX_idx[jj]] = delta_yFX[jj];
1089  for ( jj = 0; jj < nAC; ++jj )
1090  y_out[nV+AC_idx[jj]] = delta_yAC[jj];
1091  for ( jj = 0; jj < nFR; ++jj )
1092  y_out[FR_idx[jj]] = 0.0;
1093 
1094  g_in += nV;
1095  lb_in += nV;
1096  ub_in += nV;
1097  lbA_in += nC;
1098  ubA_in += nC;
1099  x_out += nV;
1100  y_out += nV+nC;
1101  }
1102 
1103 
1104  delete[] delta_yFX;
1105  delete[] delta_yAC;
1106  delete[] delta_xFR;
1107  delete[] delta_xFX;
1108 
1109  return returnvalue;
1110 }
1111 
1112 
1113 
1114 /*
1115  * g e t N Z
1116  */
1117 int QProblem::getNZ( ) const
1118 {
1119  /* nZ = nFR - nAC */
1120  return getNFR( ) - getNAC( );
1121 }
1122 
1123 
1124 /*
1125  * g e t D u a l S o l u t i o n
1126  */
1127 returnValue QProblem::getDualSolution( real_t* const yOpt ) const
1128 {
1129  int i;
1130 
1131  /* return optimal dual solution vector
1132  * only if current QP has been solved */
1133  if ( ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
1134  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
1135  ( getStatus( ) == QPS_SOLVED ) )
1136  {
1137  for( i=0; i<getNV( )+getNC( ); ++i )
1138  yOpt[i] = y[i];
1139 
1140  return SUCCESSFUL_RETURN;
1141  }
1142  else
1143  {
1144  return RET_QP_NOT_SOLVED;
1145  }
1146 }
1147 
1148 
1149 
1150 /*
1151  * s e t C o n s t r a i n t P r o d u c t
1152  */
1153 returnValue QProblem::setConstraintProduct( ConstraintProduct* const _constraintProduct )
1154 {
1155  constraintProduct = _constraintProduct;
1156 
1157  return SUCCESSFUL_RETURN;
1158 }
1159 
1160 
1161 /*
1162  * p r i n t P r o p e r t i e s
1163  */
1164 returnValue QProblem::printProperties( )
1165 {
1166  #ifndef __XPCTARGET__
1167  /* Do not print properties if print level is set to none! */
1168  if ( options.printLevel == PL_NONE )
1169  return SUCCESSFUL_RETURN;
1170 
1171  char myPrintfString[80];
1172 
1173  myPrintf( "\n################# qpOASES -- QP PROPERTIES #################\n" );
1174  myPrintf( "\n" );
1175 
1176  /* 1) Variables properties. */
1177  snprintf( myPrintfString,80, "Number of Variables: %4.1d\n",getNV( ) );
1178  myPrintf( myPrintfString );
1179 
1180  if ( bounds.hasNoLower( ) == BT_TRUE )
1181  myPrintf( "Variables are not bounded from below.\n" );
1182  else
1183  myPrintf( "Variables are bounded from below.\n" );
1184 
1185  if ( bounds.hasNoUpper( ) == BT_TRUE )
1186  myPrintf( "Variables are not bounded from above.\n" );
1187  else
1188  myPrintf( "Variables are bounded from above.\n" );
1189 
1190  myPrintf( "\n" );
1191 
1192 
1193  /* 2) Constraints properties. */
1194  snprintf( myPrintfString,80, "Total number of Constraints: %4.1d\n",getNC( ) );
1195  myPrintf( myPrintfString );
1196 
1197  snprintf( myPrintfString,80, "Number of Equality Constraints: %4.1d\n",getNEC( ) );
1198  myPrintf( myPrintfString );
1199 
1200  snprintf( myPrintfString,80, "Number of Inequality Constraints: %4.1d\n",getNC( )-getNEC( ) );
1201  myPrintf( myPrintfString );
1202 
1203  if ( getNC( ) > 0 )
1204  {
1205  if ( constraints.hasNoLower( ) == BT_TRUE )
1206  myPrintf( "Constraints are not bounded from below.\n" );
1207  else
1208  myPrintf( "Constraints are bounded from below.\n" );
1209 
1210  if ( constraints.hasNoUpper( ) == BT_TRUE )
1211  myPrintf( "Constraints are not bounded from above.\n" );
1212  else
1213  myPrintf( "Constraints are bounded from above.\n" );
1214  }
1215 
1216  myPrintf( "\n" );
1217 
1218 
1219  /* 3) Further properties. */
1220  switch ( hessianType )
1221  {
1222  case HST_ZERO:
1223  myPrintf( "Hessian is zero matrix (i.e. actually an LP is solved).\n" );
1224  break;
1225 
1226  case HST_IDENTITY:
1227  myPrintf( "Hessian is identity matrix.\n" );
1228  break;
1229 
1230  case HST_POSDEF:
1231  myPrintf( "Hessian matrix is (strictly) positive definite.\n" );
1232  break;
1233 
1234  case HST_POSDEF_NULLSPACE:
1235  myPrintf( "Hessian matrix is positive definite on null space of active constraints.\n" );
1236  break;
1237 
1238  case HST_SEMIDEF:
1239  myPrintf( "Hessian matrix is positive semi-definite.\n" );
1240  break;
1241 
1242  default:
1243  myPrintf( "Hessian matrix has unknown type.\n" );
1244  break;
1245  }
1246 
1247  if ( infeasible == BT_TRUE )
1248  myPrintf( "QP was found to be infeasible.\n" );
1249  else
1250  myPrintf( "QP seems to be feasible.\n" );
1251 
1252  if ( unbounded == BT_TRUE )
1253  myPrintf( "QP was found to be unbounded from below.\n" );
1254  else
1255  myPrintf( "QP seems to be bounded from below.\n" );
1256 
1257  myPrintf( "\n" );
1258 
1259 
1260  /* 4) QP object properties. */
1261  switch ( status )
1262  {
1263  case QPS_NOTINITIALISED:
1264  myPrintf( "Status of QP object: freshly instantiated or reset.\n" );
1265  break;
1266 
1267  case QPS_PREPARINGAUXILIARYQP:
1268  myPrintf( "Status of QP object: an auxiliary QP is currently setup.\n" );
1269  break;
1270 
1271  case QPS_AUXILIARYQPSOLVED:
1272  myPrintf( "Status of QP object: an auxilary QP was solved.\n" );
1273  break;
1274 
1275  case QPS_PERFORMINGHOMOTOPY:
1276  myPrintf( "Status of QP object: a homotopy step is performed.\n" );
1277  break;
1278 
1279  case QPS_HOMOTOPYQPSOLVED:
1280  myPrintf( "Status of QP object: an intermediate QP along the homotopy path was solved.\n" );
1281  break;
1282 
1283  case QPS_SOLVED:
1284  myPrintf( "Status of QP object: solution of the actual QP was found.\n" );
1285  break;
1286  }
1287 
1288  switch ( options.printLevel )
1289  {
1290  case PL_LOW:
1291  myPrintf( "Print level of QP object is low, i.e. only error are printed.\n" );
1292  break;
1293 
1294  case PL_MEDIUM:
1295  myPrintf( "Print level of QP object is medium, i.e. error and warnings are printed.\n" );
1296  break;
1297 
1298  case PL_HIGH:
1299  myPrintf( "Print level of QP object is high, i.e. all available output is printed.\n" );
1300  break;
1301 
1302  default:
1303  break;
1304  }
1305 
1306  myPrintf( "\n" );
1307  #endif
1308 
1309  return SUCCESSFUL_RETURN;
1310 }
1311 
1312 
1313 
1314 /*****************************************************************************
1315  * P R O T E C T E D *
1316  *****************************************************************************/
1317 
1318 /*
1319  * c l e a r
1320  */
1321 returnValue QProblem::clear( )
1322 {
1323  if ( ( freeConstraintMatrix == BT_TRUE ) && ( A != 0 ) )
1324  {
1325  delete A;
1326  A = 0;
1327  }
1328 
1329  if ( lbA != 0 )
1330  {
1331  delete[] lbA;
1332  lbA = 0;
1333  }
1334 
1335  if ( ubA != 0 )
1336  {
1337  delete[] ubA;
1338  ubA = 0;
1339  }
1340 
1341  if ( T != 0 )
1342  {
1343  delete[] T;
1344  T = 0;
1345  }
1346 
1347  if ( Q != 0 )
1348  {
1349  delete[] Q;
1350  Q = 0;
1351  }
1352 
1353  if ( Ax != 0 )
1354  {
1355  delete[] Ax;
1356  Ax = 0;
1357  }
1358 
1359  if ( Ax_l != 0 )
1360  {
1361  delete[] Ax_l;
1362  Ax_l = 0;
1363  }
1364 
1365  if ( Ax_u != 0 )
1366  {
1367  delete[] Ax_u;
1368  Ax_u = 0;
1369  }
1370 
1371  if ( tempA != 0 )
1372  {
1373  delete[] tempA;
1374  tempA = 0;
1375  }
1376 
1377  if ( ZFR_delta_xFRz != 0 )
1378  {
1379  delete[] ZFR_delta_xFRz;
1380  ZFR_delta_xFRz = 0;
1381  }
1382 
1383  if ( delta_xFRy != 0 )
1384  {
1385  delete[] delta_xFRy;
1386  delta_xFRy = 0;
1387  }
1388 
1389  if ( delta_xFRz != 0 )
1390  {
1391  delete[] delta_xFRz;
1392  delta_xFRz = 0;
1393  }
1394 
1395  if ( tempB != 0 )
1396  {
1397  delete[] tempB;
1398  tempB = 0;
1399  }
1400 
1401  if ( delta_yAC_TMP != 0 )
1402  {
1403  delete[] delta_yAC_TMP;
1404  delta_yAC_TMP = 0;
1405  }
1406 
1407  return SUCCESSFUL_RETURN;
1408 }
1409 
1410 
1411 /*
1412  * c o p y
1413  */
1414 returnValue QProblem::copy( const QProblem& rhs
1415  )
1416 {
1417  int _nV = rhs.getNV( );
1418  int _nC = rhs.getNC( );
1419 
1420  constraints = rhs.constraints;
1421 
1422  if ( ( freeConstraintMatrix == BT_TRUE ) && ( A != 0 ) )
1423  {
1424  delete A;
1425  A = 0;
1426  }
1427 
1428  freeConstraintMatrix = rhs.freeConstraintMatrix;
1429 
1430  if ( freeConstraintMatrix == BT_TRUE )
1431  A = rhs.A->duplicate();
1432  else
1433  A = rhs.A;
1434 
1435  if ( rhs.lbA != 0 )
1436  {
1437  lbA = new real_t[_nC];
1438  setLBA( rhs.lbA );
1439  }
1440  else
1441  lbA = 0;
1442 
1443  if ( rhs.ubA != 0 )
1444  {
1445  ubA = new real_t[_nC];
1446  setUBA( rhs.ubA );
1447  }
1448  else
1449  ubA = 0;
1450 
1451  if ( rhs.y != 0 )
1452  {
1453  delete[] y; /* y of no constraints version too short! */
1454  y = new real_t[_nV+_nC];
1455  memcpy( y,rhs.y,(_nV+_nC)*sizeof(real_t) );
1456  }
1457  else
1458  y = 0;
1459 
1460  sizeT = rhs.sizeT;
1461 
1462  if ( rhs.T != 0 )
1463  {
1464  T = new real_t[sizeT*sizeT];
1465  memcpy( T,rhs.T,sizeT*sizeT*sizeof(real_t) );
1466  }
1467  else
1468  T = 0;
1469 
1470  if ( rhs.Q != 0 )
1471  {
1472  Q = new real_t[_nV*_nV];
1473  memcpy( Q,rhs.Q,_nV*_nV*sizeof(real_t) );
1474  }
1475  else
1476  Q = 0;
1477 
1478  if ( rhs.Ax != 0 )
1479  {
1480  Ax = new real_t[_nC];
1481  memcpy( Ax,rhs.Ax,_nC*sizeof(real_t) );
1482  }
1483  else
1484  Ax = 0;
1485 
1486  if ( rhs.Ax_l != 0 )
1487  {
1488  Ax_l = new real_t[_nC];
1489  memcpy( Ax_l,rhs.Ax_l,_nC*sizeof(real_t) );
1490  }
1491  else
1492  Ax_l = 0;
1493 
1494  if ( rhs.Ax_u != 0 )
1495  {
1496  Ax_u = new real_t[_nC];
1497  memcpy( Ax_u,rhs.Ax_u,_nC*sizeof(real_t) );
1498  }
1499  else
1500  Ax_u = 0;
1501 
1502  if ( rhs.constraintProduct != 0 )
1503  constraintProduct = rhs.constraintProduct;
1504  else
1505  constraintProduct = 0;
1506 
1507  tempA = new real_t[_nV]; /* nFR */
1508  ZFR_delta_xFRz = new real_t[_nV]; /* nFR */
1509  delta_xFRz = new real_t[_nV]; /* nZ */
1510 
1511  if ( _nC > 0 )
1512  {
1513  delta_xFRy = new real_t[_nC]; /* nAC */
1514  tempB = new real_t[_nC]; /* nAC */
1515  delta_yAC_TMP = new real_t[_nC]; /* nAC */
1516  }
1517  else
1518  {
1519  delta_xFRy = 0;
1520  tempB = 0;
1521  delta_yAC_TMP = 0;
1522  }
1523 
1524  return SUCCESSFUL_RETURN;
1525 }
1526 
1527 
1528 
1529 /*
1530  * s o l v e I n i t i a l Q P
1531  */
1532 returnValue QProblem::solveInitialQP( const real_t* const xOpt, const real_t* const yOpt,
1533  const Bounds* const guessedBounds, const Constraints* const guessedConstraints,
1534  int& nWSR, real_t* const cputime
1535  )
1536 {
1537  int i;
1538 
1539  /* some definitions */
1540  int nV = getNV( );
1541  int nC = getNC( );
1542 
1543 
1544  /* start runtime measurement */
1545  real_t starttime = 0.0;
1546  if ( cputime != 0 )
1547  starttime = getCPUtime( );
1548 
1549  status = QPS_NOTINITIALISED;
1550 
1551 
1552  /* I) ANALYSE QP DATA: */
1553  #ifdef __MANY_CONSTRAINTS__
1554  /* 0) Checks if l1-norm of each constraint is not greater than 1. */
1555  int j;
1556  real_t l1;
1557  for( i=0; i<nC; ++i )
1558  {
1559  l1 = 0.0;
1560  for( j=0; j<nV; ++j )
1561  l1 += fabs( AA(i,j) ); /* Andreas: do many constraints case later */
1562 
1563  if ( l1 > 1.0 + 10.0*EPS )
1564  return THROWERROR( RET_CONSTRAINTS_ARE_NOT_SCALED );
1565  }
1566  #endif
1567 
1568  /* 1) Check if Hessian happens to be the identity matrix. */
1569  if ( determineHessianType( ) != SUCCESSFUL_RETURN )
1570  return THROWERROR( RET_INIT_FAILED );
1571 
1572  /* 2) Setup type of bounds and constraints (i.e. unbounded, implicitly fixed etc.). */
1573  if ( setupSubjectToType( ) != SUCCESSFUL_RETURN )
1574  return THROWERROR( RET_INIT_FAILED );
1575 
1576  status = QPS_PREPARINGAUXILIARYQP;
1577 
1578 
1579  /* II) SETUP AUXILIARY QP WITH GIVEN OPTIMAL SOLUTION: */
1580  /* 1) Setup bounds and constraints data structure. */
1581  if ( bounds.setupAllFree( ) != SUCCESSFUL_RETURN )
1582  return THROWERROR( RET_INIT_FAILED );
1583 
1584  if ( constraints.setupAllInactive( ) != SUCCESSFUL_RETURN )
1585  return THROWERROR( RET_INIT_FAILED );
1586 
1587  /* 2) Setup optimal primal/dual solution for auxiliary QP. */
1588  if ( setupAuxiliaryQPsolution( xOpt,yOpt ) != SUCCESSFUL_RETURN )
1589  return THROWERROR( RET_INIT_FAILED );
1590 
1591  /* 3) Obtain linear independent working set for auxiliary QP. */
1592  Bounds auxiliaryBounds( nV );
1593  Constraints auxiliaryConstraints( nC );
1594 
1595  if ( obtainAuxiliaryWorkingSet( xOpt,yOpt,guessedBounds,guessedConstraints,
1596  &auxiliaryBounds,&auxiliaryConstraints ) != SUCCESSFUL_RETURN )
1597  return THROWERROR( RET_INIT_FAILED );
1598 
1599  /* 4) Setup working set of auxiliary QP and setup matrix factorisations. */
1600  /* a) Regularise Hessian if necessary. */
1601  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_SEMIDEF ) )
1602  {
1603  if ( regulariseHessian( ) != SUCCESSFUL_RETURN )
1604  return THROWERROR( RET_INIT_FAILED_REGULARISATION );
1605  }
1606 
1607  /* b) TQ factorisation. */
1608  if ( setupTQfactorisation( ) != SUCCESSFUL_RETURN )
1609  return THROWERROR( RET_INIT_FAILED_TQ );
1610 
1611  /* c) Working set of auxiliary QP. */
1612  if ( setupAuxiliaryWorkingSet( &auxiliaryBounds,&auxiliaryConstraints,BT_TRUE ) != SUCCESSFUL_RETURN )
1613  return THROWERROR( RET_INIT_FAILED );
1614 
1615  haveCholesky = BT_FALSE;
1616 
1617  /* 5) Store original QP formulation... */
1618  real_t* g_original = new real_t[nV];
1619  real_t* lb_original = new real_t[nV];
1620  real_t* ub_original = new real_t[nV];
1621  real_t* lbA_original = new real_t[nC];
1622  real_t* ubA_original = new real_t[nC];
1623 
1624  for( i=0; i<nV; ++i )
1625  {
1626  g_original[i] = g[i];
1627  lb_original[i] = lb[i];
1628  ub_original[i] = ub[i];
1629  }
1630 
1631  for( i=0; i<nC; ++i )
1632  {
1633  lbA_original[i] = lbA[i];
1634  ubA_original[i] = ubA[i];
1635  }
1636 
1637  /* ... and setup QP data of an auxiliary QP having an optimal solution
1638  * as specified by the user (or xOpt = yOpt = 0, by default). */
1639  if ( setupAuxiliaryQPgradient( ) != SUCCESSFUL_RETURN )
1640  {
1641  delete[] ubA_original; delete[] lbA_original; delete[] ub_original; delete[] lb_original; delete[] g_original;
1642  return THROWERROR( RET_INIT_FAILED );
1643  }
1644 
1645  if ( setupAuxiliaryQPbounds( &auxiliaryBounds,&auxiliaryConstraints,BT_TRUE ) != SUCCESSFUL_RETURN )
1646  {
1647  delete[] ubA_original; delete[] lbA_original; delete[] ub_original; delete[] lb_original; delete[] g_original;
1648  return THROWERROR( RET_INIT_FAILED );
1649  }
1650 
1651  status = QPS_AUXILIARYQPSOLVED;
1652 
1653 
1654  if ( options.enableRamping == BT_TRUE )
1655  performRamping( );
1656 
1657 
1658  /* III) SOLVE ACTUAL INITIAL QP: */
1659  /* Allow only remaining CPU time for usual hotstart. */
1660  if ( cputime != 0 )
1661  *cputime -= getCPUtime( ) - starttime;
1662 
1663  /* Use hotstart method to find the solution of the original initial QP,... */
1664  returnValue returnvalue = hotstart( g_original,lb_original,ub_original,lbA_original,ubA_original, nWSR,cputime );
1665 
1666  /* ... deallocate memory,... */
1667  delete[] ubA_original; delete[] lbA_original; delete[] ub_original; delete[] lb_original; delete[] g_original;
1668 
1669  /* ... check for infeasibility and unboundedness... */
1670  if ( isInfeasible( ) == BT_TRUE )
1671  return THROWERROR( RET_INIT_FAILED_INFEASIBILITY );
1672 
1673  if ( isUnbounded( ) == BT_TRUE )
1674  return THROWERROR( RET_INIT_FAILED_UNBOUNDEDNESS );
1675 
1676  /* ... and internal errors. */
1677  if ( ( returnvalue != SUCCESSFUL_RETURN ) && ( returnvalue != RET_MAX_NWSR_REACHED ) )
1678  return THROWERROR( RET_INIT_FAILED_HOTSTART );
1679 
1680 
1681  /* stop runtime measurement */
1682  if ( cputime != 0 )
1683  *cputime = getCPUtime( ) - starttime;
1684 
1685  THROWINFO( RET_INIT_SUCCESSFUL );
1686 
1687  return returnvalue;
1688 }
1689 
1690 
1691 /*
1692  * s o l v e Q P
1693  */
1694 returnValue QProblem::solveQP( const real_t* const g_new,
1695  const real_t* const lb_new, const real_t* const ub_new,
1696  const real_t* const lbA_new, const real_t* const ubA_new,
1697  int& nWSR, real_t* const cputime, int nWSRperformed
1698  )
1699 {
1700  int iter;
1701  int nV = getNV( );
1702  int nC = getNC( );
1703 
1704  /* consistency check */
1705  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
1706  ( getStatus( ) == QPS_PREPARINGAUXILIARYQP ) ||
1707  ( getStatus( ) == QPS_PERFORMINGHOMOTOPY ) )
1708  {
1709  return THROWERROR( RET_HOTSTART_FAILED_AS_QP_NOT_INITIALISED );
1710  }
1711 
1712  /* start runtime measurement */
1713  real_t starttime = 0.0;
1714  if ( cputime != 0 )
1715  starttime = getCPUtime( );
1716 
1717 
1718  /* I) PREPARATIONS */
1719  /* 1) Allocate delta vectors of gradient and (constraints') bounds,
1720  * index arrays and step direction arrays. */
1721  real_t* delta_xFR = new real_t[nV];
1722  real_t* delta_xFX = new real_t[nV];
1723  real_t* delta_yAC = new real_t[nC];
1724  real_t* delta_yFX = new real_t[nV];
1725 
1726  real_t* delta_g = new real_t[nV];
1727  real_t* delta_lb = new real_t[nV];
1728  real_t* delta_ub = new real_t[nV];
1729  real_t* delta_lbA = new real_t[nC];
1730  real_t* delta_ubA = new real_t[nC];
1731 
1732  returnValue returnvalue;
1733  BooleanType Delta_bC_isZero, Delta_bB_isZero;
1734 
1735  int BC_idx;
1736  SubjectToStatus BC_status;
1737  BooleanType BC_isBound;
1738 
1739  real_t homotopyLength;
1740 
1741  char messageString[80];
1742 
1743  /* 2) Update type of bounds and constraints, e.g.
1744  * a former equality constraint might have become a normal one etc. */
1745  if ( setupSubjectToType( lb_new,ub_new,lbA_new,ubA_new ) != SUCCESSFUL_RETURN )
1746  return THROWERROR( RET_HOTSTART_FAILED );
1747 
1748  /* 3) Reset status flags. */
1749  infeasible = BT_FALSE;
1750  unbounded = BT_FALSE;
1751 
1752 
1753  /* II) MAIN HOMOTOPY LOOP */
1754  for( iter=nWSRperformed; iter<nWSR; ++iter )
1755  {
1756  idxAddB = idxRemB = idxAddC = idxRemC = -1;
1757 
1758  if ( isCPUtimeLimitExceeded( cputime,starttime,iter-nWSRperformed ) == BT_TRUE )
1759  {
1760  /* If CPU time limit is exceeded, stop homotopy loop immediately!
1761  * Assign number of working set recalculations (runtime measurement is stopped later). */
1762  nWSR = iter;
1763  break;
1764  }
1765 
1766  status = QPS_PERFORMINGHOMOTOPY;
1767 
1768  #ifndef __XPCTARGET__
1769  snprintf( messageString,80,"%d ...",iter );
1770  getGlobalMessageHandler( )->throwInfo( RET_ITERATION_STARTED,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
1771  #endif
1772 
1773  /* 2) Detemination of shift direction of the gradient and the (constraints') bounds. */
1774  returnvalue = determineDataShift( g_new,lbA_new,ubA_new,lb_new,ub_new,
1775  delta_g,delta_lbA,delta_ubA,delta_lb,delta_ub,
1776  Delta_bC_isZero, Delta_bB_isZero
1777  );
1778  if ( returnvalue != SUCCESSFUL_RETURN )
1779  {
1780  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1781  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1782 
1783  /* Assign number of working set recalculations and stop runtime measurement. */
1784  nWSR = iter;
1785  if ( cputime != 0 )
1786  *cputime = getCPUtime( ) - starttime;
1787 
1788  THROWERROR( RET_SHIFT_DETERMINATION_FAILED );
1789  return returnvalue;
1790  }
1791 
1792  /* 3) Determination of step direction of X and Y. */
1793  returnvalue = determineStepDirection( delta_g,delta_lbA,delta_ubA,delta_lb,delta_ub,
1794  Delta_bC_isZero, Delta_bB_isZero,
1795  delta_xFX,delta_xFR,delta_yAC,delta_yFX
1796  );
1797  if ( returnvalue != SUCCESSFUL_RETURN )
1798  {
1799  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1800  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1801 
1802  /* Assign number of working set recalculations and stop runtime measurement. */
1803  nWSR = iter;
1804  if ( cputime != 0 )
1805  *cputime = getCPUtime( ) - starttime;
1806 
1807  THROWERROR( RET_STEPDIRECTION_DETERMINATION_FAILED );
1808  return returnvalue;
1809  }
1810 
1811  /* 4) Determination of step length TAU.
1812  * This step along the homotopy path is also taken (without changing working set). */
1813  returnvalue = performStep( delta_g, delta_lbA,delta_ubA,delta_lb,delta_ub,
1814  delta_xFX,delta_xFR,delta_yAC,delta_yFX,
1815  BC_idx,BC_status,BC_isBound
1816  );
1817  if ( returnvalue != SUCCESSFUL_RETURN )
1818  {
1819  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1820  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1821 
1822  /* Assign number of working set recalculations and stop runtime measurement. */
1823  nWSR = iter;
1824  if ( cputime != 0 )
1825  *cputime = getCPUtime( ) - starttime;
1826 
1827  THROWERROR( RET_STEPLENGTH_DETERMINATION_FAILED );
1828  return returnvalue;
1829  }
1830 
1831  /* 5) Termination criterion. */
1832  homotopyLength = relativeHomotopyLength(g_new, lb_new, ub_new, lbA_new, ubA_new);
1833  if ( homotopyLength <= options.terminationTolerance )
1834  {
1835  status = QPS_SOLVED;
1836 
1837  THROWINFO( RET_OPTIMAL_SOLUTION_FOUND );
1838 
1839  if ( printIteration( iter,BC_idx,BC_status,BC_isBound ) != SUCCESSFUL_RETURN )
1840  THROWERROR( RET_PRINT_ITERATION_FAILED ); /* do not pass this as return value! */
1841 
1842  nWSR = iter;
1843  if ( cputime != 0 )
1844  *cputime = getCPUtime( ) - starttime;
1845 
1846  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1847  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1848 
1849  return SUCCESSFUL_RETURN;
1850  }
1851 
1852 
1853  /* 6) Change active set. */
1854  returnvalue = changeActiveSet( BC_idx,BC_status,BC_isBound );
1855  if ( returnvalue != SUCCESSFUL_RETURN )
1856  {
1857  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1858  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1859 
1860  /* Assign number of working set recalculations and stop runtime measurement. */
1861  nWSR = iter;
1862  if ( cputime != 0 )
1863  *cputime = getCPUtime( ) - starttime;
1864 
1865  /* Checks for infeasibility... */
1866  if ( isInfeasible( ) == BT_TRUE )
1867  {
1868  status = QPS_HOMOTOPYQPSOLVED;
1869  return setInfeasibilityFlag( RET_HOTSTART_STOPPED_INFEASIBILITY );
1870  }
1871 
1872  /* ...unboundedness... */
1873  if ( unbounded == BT_TRUE ) /* not necessary since objective function convex! */
1874  return THROWERROR( RET_HOTSTART_STOPPED_UNBOUNDEDNESS );
1875 
1876  /* ... and throw unspecific error otherwise */
1877  THROWERROR( RET_HOMOTOPY_STEP_FAILED );
1878  return returnvalue;
1879  }
1880 
1881  /* 6.5) Possibly refactorise projected Hessian from scratch. */
1882  if (options.enableCholeskyRefactorisation > 0 && iter % options.enableCholeskyRefactorisation == 0)
1883  {
1884  returnvalue = setupCholeskyDecompositionProjected();
1885  if (returnvalue != SUCCESSFUL_RETURN)
1886  {
1887  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1888  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1889  return returnvalue;
1890  }
1891  }
1892 
1893 
1894 #ifdef __DEBUG_ITER__
1895  nFR = getNFR();
1896  nAC = getNAC();
1897  int i;
1898  real_t stat, bfeas, cfeas, bcmpl, ccmpl, errUnitary, errTQ, Tmaxomin;
1899  real_t *grad = new real_t[nV];
1900  real_t *AX = new real_t[nC];
1901  real_t Tmin, Tmax;
1902 
1903  stat = bfeas = cfeas = bcmpl = ccmpl = errUnitary = errTQ = Tmaxomin = 0.0;
1904 
1905  /* stationarity */
1906  for (i = 0; i < nV; i++) grad[i] = g[i] - y[i];
1907  H->times(1, 1.0, x, nV, 1.0, grad, nV);
1908  A->transTimes(1, -1.0, y+nV, nC, 1.0, grad, nV);
1909  for (i = 0; i < nV; i++) if (fabs(grad[i]) > stat) stat = fabs(grad[i]);
1910 
1911  /* feasibility */
1912  for (i = 0; i < nV; i++) if (lb[i] - x[i] > bfeas) bfeas = lb[i] - x[i];
1913  for (i = 0; i < nV; i++) if (x[i] - ub[i] > bfeas) bfeas = x[i] - ub[i];
1914  A->times(1, 1.0, x, nV, 0.0, AX, nC);
1915  for (i = 0; i < nC; i++) if (lbA[i] - AX[i] > cfeas) cfeas = lbA[i] - AX[i];
1916  for (i = 0; i < nC; i++) if (AX[i] - ubA[i] > cfeas) cfeas = AX[i] - ubA[i];
1917 
1918  /* complementarity */
1919  for (i = 0; i < nV; i++) if (y[i] > +EPS && fabs((lb[i] - x[i])*y[i]) > bcmpl) bcmpl = fabs((lb[i] - x[i])*y[i]);
1920  for (i = 0; i < nV; i++) if (y[i] < -EPS && fabs((ub[i] - x[i])*y[i]) > bcmpl) bcmpl = fabs((ub[i] - x[i])*y[i]);
1921  for (i = 0; i < nC; i++) if (y[nV+i] > +EPS && fabs((lbA[i]-AX[i])*y[nV+i]) > ccmpl) ccmpl = fabs((lbA[i]-AX[i])*y[nV+i]);
1922  for (i = 0; i < nC; i++) if (y[nV+i] < -EPS && fabs((ubA[i]-AX[i])*y[nV+i]) > ccmpl) ccmpl = fabs((ubA[i]-AX[i])*y[nV+i]);
1923 
1924  Tmin = 1.0e16; Tmax = 0.0;
1925  for (i = 0; i < nAC; i++)
1926  if (fabs(TT(i,sizeT-i-1)) < Tmin)
1927  Tmin = fabs(TT(i,sizeT-i-1));
1928  else if (fabs(TT(i,sizeT-i-1)) > Tmax)
1929  Tmax = fabs(TT(i,sizeT-i-1));
1930  Tmaxomin = Tmax/Tmin;
1931 
1932  if (iter % 10 == 0)
1933  fprintf(stderr, "\n%5s %4s %4s %4s %4s %9s %9s %9s %9s %9s %9s %9s %9s\n",
1934  "iter", "addB", "remB", "addC", "remC", "hom len", "tau", "stat",
1935  "bfeas", "cfeas", "bcmpl", "ccmpl", "Tmin");
1936  fprintf(stderr, "%5d ", iter);
1937  if (idxAddB >= 0) fprintf(stderr, "%4d ", idxAddB);
1938  else fprintf(stderr, "%4s ", " ");
1939  if (idxRemB >= 0) fprintf(stderr, "%4d ", idxRemB);
1940  else fprintf(stderr, "%4s ", " ");
1941  if (idxAddC >= 0) fprintf(stderr, "%4d ", idxAddC);
1942  else fprintf(stderr, "%4s ", " ");
1943  if (idxRemC >= 0) fprintf(stderr, "%4d ", idxRemC);
1944  else fprintf(stderr, "%4s ", " ");
1945  fprintf(stderr, "%9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
1946  homotopyLength, tau, stat, bfeas, cfeas, bcmpl, ccmpl, Tmin);
1947 
1948  delete[] AX;
1949  delete[] grad;
1950 #else
1951  if (options.printLevel == -1)
1952  {
1953  if (iter % 10 == 0)
1954  fprintf(stdout, "\n%5s %5s %5s %5s %5s %9s %9s\n",
1955  "iter", "addB", "remB", "addC", "remC", "hom len", "tau");
1956  fprintf(stdout, "%5d ", iter);
1957  if (idxAddB >= 0) fprintf(stdout, "%5d ", idxAddB);
1958  else fprintf(stdout, "%5s ", " ");
1959  if (idxRemB >= 0) fprintf(stdout, "%5d ", idxRemB);
1960  else fprintf(stdout, "%5s ", " ");
1961  if (idxAddC >= 0) fprintf(stdout, "%5d ", idxAddC);
1962  else fprintf(stdout, "%5s ", " ");
1963  if (idxRemC >= 0) fprintf(stdout, "%5d ", idxRemC);
1964  else fprintf(stdout, "%5s ", " ");
1965  fprintf(stdout, "%9.2e %9.2e\n", homotopyLength, tau);
1966  }
1967 #endif /* __DEBUG_ITER__ */
1968 
1969 
1970  /* 7) Output information of successful QP iteration. */
1971  status = QPS_HOMOTOPYQPSOLVED;
1972 
1973  if ( printIteration( iter,BC_idx,BC_status,BC_isBound ) != SUCCESSFUL_RETURN )
1974  THROWERROR( RET_PRINT_ITERATION_FAILED ); /* do not pass this as return value! */
1975 
1976  /* 8) Perform Ramping Strategy on zero homotopy step or drift correction (if desired). */
1977  if ( ( tau < EPS ) && ( options.enableRamping == BT_TRUE ) )
1978  performRamping( );
1979  else
1980  if ( (options.enableDriftCorrection > 0)
1981  && ((iter+1) % options.enableDriftCorrection == 0) )
1982  performDriftCorrection( ); /* always returns SUCCESSFUL_RETURN */
1983  }
1984 
1985  delete[] delta_yAC; delete[] delta_yFX; delete[] delta_xFX; delete[] delta_xFR;
1986  delete[] delta_ub; delete[] delta_lb; delete[] delta_ubA; delete[] delta_lbA; delete[] delta_g;
1987 
1988  /* stop runtime measurement */
1989  if ( cputime != 0 )
1990  *cputime = getCPUtime( ) - starttime;
1991 
1992 
1993  /* if programm gets to here, output information that QP could not be solved
1994  * within the given maximum numbers of working set changes */
1995  if ( options.printLevel == PL_HIGH )
1996  {
1997  #ifndef __XPCTARGET__
1998  snprintf( messageString,80,"(nWSR = %d)",iter );
1999  return getGlobalMessageHandler( )->throwWarning( RET_MAX_NWSR_REACHED,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
2000  #else
2001  return RET_MAX_NWSR_REACHED;
2002  #endif
2003  }
2004  else
2005  {
2006  return RET_MAX_NWSR_REACHED;
2007  }
2008 }
2009 
2010 
2011 /*
2012  * s o l v e R e g u l a r i s e d Q P
2013  */
2014 returnValue QProblem::solveRegularisedQP( const real_t* const g_new,
2015  const real_t* const lb_new, const real_t* const ub_new,
2016  const real_t* const lbA_new, const real_t* const ubA_new,
2017  int& nWSR, real_t* const cputime, int nWSRperformed
2018  )
2019 {
2020  int i, step;
2021  int nV = getNV( );
2022 
2023 
2024  /* Stop here if QP has not been regularised (i.e. normal QP solution). */
2025  if ( usingRegularisation( ) == BT_FALSE )
2026  return solveQP( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,cputime,nWSRperformed );
2027 
2028 
2029  /* I) SOLVE USUAL REGULARISED QP */
2030  returnValue returnvalue;
2031 
2032  int nWSR_max = nWSR;
2033  int nWSR_total = nWSRperformed;
2034 
2035  real_t cputime_total = 0.0;
2036  real_t cputime_cur = 0.0;
2037 
2038  if ( cputime == 0 )
2039  {
2040  returnvalue = solveQP( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,0,nWSRperformed );
2041  }
2042  else
2043  {
2044  cputime_cur = *cputime;
2045  returnvalue = solveQP( g_new,lb_new,ub_new,lbA_new,ubA_new, nWSR,&cputime_cur,nWSRperformed );
2046  }
2047  nWSR_total = nWSR;
2048  cputime_total += cputime_cur;
2049 
2050  /* Only continue if QP solution has been successful. */
2051  if ( returnvalue != SUCCESSFUL_RETURN )
2052  {
2053  if ( cputime != 0 )
2054  *cputime = cputime_total;
2055 
2056  if ( returnvalue == RET_MAX_NWSR_REACHED )
2057  THROWWARNING( RET_NO_REGSTEP_NWSR );
2058 
2059  return returnvalue;
2060  }
2061 
2062 
2063  /* II) PERFORM SUCCESSIVE REGULARISATION STEPS */
2064  real_t* gMod = new real_t[nV];
2065 
2066  for( step=0; step<options.numRegularisationSteps; ++step )
2067  {
2068  /* 1) Modify gradient: gMod = g - eps*xOpt
2069  * (assuming regularisation matrix to be eps*Id). */
2070  for( i=0; i<nV; ++i )
2071  gMod[i] = g_new[i] - options.epsRegularisation*x[i];
2072 
2073  /* 2) Solve regularised QP with modified gradient allowing
2074  * only as many working set recalculations and CPU time
2075  * as have been left from previous QP solutions. */
2076  if ( cputime == 0 )
2077  {
2078  nWSR = nWSR_max;
2079  returnvalue = solveQP( gMod,lb_new,ub_new,lbA_new,ubA_new, nWSR,0,nWSR_total );
2080  }
2081  else
2082  {
2083  nWSR = nWSR_max;
2084  cputime_cur = *cputime - cputime_total;
2085  returnvalue = solveQP( gMod,lb_new,ub_new,lbA_new,ubA_new, nWSR,&cputime_cur,nWSR_total );
2086  }
2087 
2088  nWSR_total = nWSR;
2089  cputime_total += cputime_cur;
2090 
2091  /* Only continue if QP solution has been successful. */
2092  if ( returnvalue != SUCCESSFUL_RETURN )
2093  {
2094  delete[] gMod;
2095 
2096  if ( cputime != 0 )
2097  *cputime = cputime_total;
2098 
2099  if ( returnvalue == RET_MAX_NWSR_REACHED )
2100  THROWWARNING( RET_FEWER_REGSTEPS_NWSR );
2101 
2102  return returnvalue;
2103  }
2104  }
2105 
2106  delete[] gMod;
2107 
2108  if ( cputime != 0 )
2109  *cputime = cputime_total;
2110 
2111  return SUCCESSFUL_RETURN;
2112 }
2113 
2114 
2115 /*
2116  * s e t u p S u b j e c t T o T y p e
2117  */
2118 returnValue QProblem::setupSubjectToType( )
2119 {
2120  return setupSubjectToType( lb,ub,lbA,ubA );
2121 }
2122 
2123 
2124 /*
2125  * s e t u p S u b j e c t T o T y p e
2126  */
2127 returnValue QProblem::setupSubjectToType( const real_t* const lb_new, const real_t* const ub_new,
2128  const real_t* const lbA_new, const real_t* const ubA_new
2129  )
2130 {
2131  int i;
2132  int nC = getNC( );
2133 
2134 
2135  /* I) SETUP SUBJECTTOTYPE FOR BOUNDS */
2136  if ( QProblemB::setupSubjectToType( lb_new,ub_new ) != SUCCESSFUL_RETURN )
2137  return THROWERROR( RET_SETUPSUBJECTTOTYPE_FAILED );
2138 
2139 
2140  /* II) SETUP SUBJECTTOTYPE FOR CONSTRAINTS */
2141  /* 1) Check if lower constraints' bounds are present. */
2142  constraints.setNoLower( BT_TRUE );
2143  if ( lbA_new != 0 )
2144  {
2145  for( i=0; i<nC; ++i )
2146  {
2147  if ( lbA_new[i] > -INFTY )
2148  {
2149  constraints.setNoLower( BT_FALSE );
2150  break;
2151  }
2152  }
2153  }
2154 
2155  /* 2) Check if upper constraints' bounds are present. */
2156  constraints.setNoUpper( BT_TRUE );
2157  if ( ubA_new != 0 )
2158  {
2159  for( i=0; i<nC; ++i )
2160  {
2161  if ( ubA_new[i] < INFTY )
2162  {
2163  constraints.setNoUpper( BT_FALSE );
2164  break;
2165  }
2166  }
2167  }
2168 
2169  /* 3) Determine implicit equality constraints and unbounded constraints. */
2170  if ( ( lbA_new != 0 ) && ( ubA_new != 0 ) )
2171  {
2172  for( i=0; i<nC; ++i )
2173  {
2174  if ( ( lbA_new[i] < -INFTY + options.boundTolerance ) && ( ubA_new[i] > INFTY - options.boundTolerance )
2175  && (options.enableFarBounds == BT_FALSE))
2176  {
2177  constraints.setType( i,ST_UNBOUNDED );
2178  }
2179  else
2180  {
2181  if ( options.enableEqualities && lbA_new[i] > ubA_new[i] - options.boundTolerance )
2182  constraints.setType( i,ST_EQUALITY );
2183  else
2184  constraints.setType( i,ST_BOUNDED );
2185  }
2186  }
2187  }
2188  else
2189  {
2190  if ( ( lbA_new == 0 ) && ( ubA_new == 0 ) )
2191  {
2192  for( i=0; i<nC; ++i )
2193  constraints.setType( i,ST_UNBOUNDED );
2194  }
2195  else
2196  {
2197  for( i=0; i<nC; ++i )
2198  constraints.setType( i,ST_BOUNDED );
2199  }
2200  }
2201 
2202  return SUCCESSFUL_RETURN;
2203 }
2204 
2205 
2206 /*
2207  * 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
2208  */
2209 returnValue QProblem::setupCholeskyDecompositionProjected( )
2210 {
2211  int i, j;
2212  int nV = getNV( );
2213  int nZ = getNZ( );
2214 
2215  /* 1) Initialises R with all zeros. */
2216  for( i=0; i<nV*nV; ++i )
2217  R[i] = 0.0;
2218 
2219  /* 2) Calculate Cholesky decomposition of projected Hessian Z'*H*Z. */
2220  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_IDENTITY ) )
2221  {
2222  if ( hessianType == HST_ZERO )
2223  {
2224  /* if Hessian is zero matrix, it is assumed that it has been
2225  * regularised and thus its Cholesky factor is the identity
2226  * matrix scaled by sqrt(eps). */
2227  if ( usingRegularisation( ) == BT_TRUE )
2228  {
2229  for( i=0; i<nV; ++i )
2230  RR(i,i) = sqrt( options.epsRegularisation );
2231  }
2232  else
2233  return THROWERROR( RET_CHOLESKY_OF_ZERO_HESSIAN );
2234  }
2235  else
2236  {
2237  /* if Hessian is identity, so is its Cholesky factor. */
2238  for( i=0; i<nV; ++i )
2239  RR(i,i) = 1.0;
2240  }
2241  }
2242  else
2243  {
2244  if ( nZ > 0 )
2245  {
2246  int* FR_idx;
2247  bounds.getFree( )->getNumberArray( &FR_idx );
2248 
2249  int* AC_idx;
2250  constraints.getActive( )->getNumberArray( &AC_idx );
2251 
2252  /* calculate Z'*H*Z */
2253  if (constraints.getActive ()->getLength () == 0) {
2254  /* make Z trivial */
2255  for ( j=0; j < nZ; ++j ) {
2256  for ( i=0; i < nV; ++i )
2257  QQ(i,j) = 0.0;
2258  QQ(FR_idx[j],j) = 1.0;
2259  }
2260  /* now Z is trivial, and so is Z'HZ */
2261  int nFR = getNFR ();
2262  for ( j=0; j < nFR; ++j )
2263  H->getCol (FR_idx[j], bounds.getFree (), 1.0, &R[j*nV]);
2264  } else {
2265  /* this is expensive if Z is large! */
2266  H->bilinear(bounds.getFree(), nZ, Q, nV, R, nV);
2267  }
2268 
2269  /* R'*R = Z'*H*Z */
2270  long info = 0;
2271  unsigned long _nZ = nZ, _nV = nV;
2272 
2273  POTRF ("U", &_nZ, R, &_nV, &info);
2274 
2275  /* <0 = invalid call, =0 ok, >0 not spd */
2276  if (info > 0) {
2277  hessianType = HST_SEMIDEF;
2278  return RET_HESSIAN_NOT_SPD;
2279  }
2280 
2281  /* zero first subdiagonal to make givens updates work */
2282  for (i=0;i<nZ-1;++i)
2283  RR(i+1,i) = 0.0;
2284  }
2285  }
2286 
2287  return SUCCESSFUL_RETURN;
2288 }
2289 
2290 
2291 /*
2292  * s e t u p T Q f a c t o r i s a t i o n
2293  */
2294 returnValue QProblem::setupTQfactorisation( )
2295 {
2296  int i, ii;
2297  int nV = getNV( );
2298  int nFR = getNFR( );
2299 
2300  int* FR_idx;
2301  bounds.getFree( )->getNumberArray( &FR_idx );
2302 
2303  /* 1) Set Q to unity matrix. */
2304  for( i=0; i<nV*nV; ++i )
2305  Q[i] = 0.0;
2306 
2307  for( i=0; i<nFR; ++i )
2308  {
2309  ii = FR_idx[i];
2310  QQ(ii,i) = 1.0;
2311  }
2312 
2313  /* 2) Set T to zero matrix. */
2314  for( i=0; i<sizeT*sizeT; ++i )
2315  T[i] = 0.0;
2316 
2317  return SUCCESSFUL_RETURN;
2318 }
2319 
2320 
2321 /*
2322  * 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
2323  */
2324 returnValue QProblem::obtainAuxiliaryWorkingSet( const real_t* const xOpt, const real_t* const yOpt,
2325  const Bounds* const guessedBounds, const Constraints* const guessedConstraints,
2326  Bounds* auxiliaryBounds, Constraints* auxiliaryConstraints
2327  ) const
2328 {
2329  int i = 0;
2330  int nV = getNV( );
2331  int nC = getNC( );
2332 
2333 
2334  /* 1) Ensure that desiredBounds is allocated (and different from guessedBounds). */
2335  if ( ( auxiliaryBounds == 0 ) || ( auxiliaryBounds == guessedBounds ) )
2336  return THROWERROR( RET_INVALID_ARGUMENTS );
2337 
2338  if ( ( auxiliaryConstraints == 0 ) || ( auxiliaryConstraints == guessedConstraints ) )
2339  return THROWERROR( RET_INVALID_ARGUMENTS );
2340 
2341 
2342  SubjectToStatus guessedStatus;
2343 
2344  /* 2) Setup working set of bounds for auxiliary initial QP. */
2345  if ( QProblemB::obtainAuxiliaryWorkingSet( xOpt,yOpt,guessedBounds, auxiliaryBounds ) != SUCCESSFUL_RETURN )
2346  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2347 
2348  /* 3) Setup working set of constraints for auxiliary initial QP. */
2349  if ( guessedConstraints != 0 )
2350  {
2351  /* If an initial working set is specific, use it!
2352  * Moreover, add all equality constraints if specified. */
2353  for( i=0; i<nC; ++i )
2354  {
2355  /* Add constraint only if it is not (goint to be) disabled! */
2356  guessedStatus = guessedConstraints->getStatus( i );
2357 
2358  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
2359  if ( constraints.getType( i ) == ST_EQUALITY )
2360  {
2361  if ( auxiliaryConstraints->setupConstraint( i,ST_LOWER ) != SUCCESSFUL_RETURN )
2362  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2363  }
2364  else
2365  #endif
2366  {
2367  if ( auxiliaryConstraints->setupConstraint( i,guessedStatus ) != SUCCESSFUL_RETURN )
2368  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2369  }
2370  }
2371  }
2372  else /* No initial working set specified. */
2373  {
2374  /* Obtain initial working set by "clipping". */
2375  if ( ( xOpt != 0 ) && ( yOpt == 0 ) )
2376  {
2377  for( i=0; i<nC; ++i )
2378  {
2379  if ( Ax[i] - lbA[i] <= options.boundTolerance )
2380  {
2381  if ( auxiliaryConstraints->setupConstraint( i,ST_LOWER ) != SUCCESSFUL_RETURN )
2382  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2383  continue;
2384  }
2385 
2386  if ( ubA[i] - Ax_u[i] <= options.boundTolerance )
2387  {
2388  if ( auxiliaryConstraints->setupConstraint( i,ST_UPPER ) != SUCCESSFUL_RETURN )
2389  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2390  continue;
2391  }
2392 
2393  /* Moreover, add all equality constraints if specified. */
2394  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
2395  if ( constraints.getType( i ) == ST_EQUALITY )
2396  {
2397  if ( auxiliaryConstraints->setupConstraint( i,ST_LOWER ) != SUCCESSFUL_RETURN )
2398  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2399  }
2400  else
2401  #endif
2402  {
2403  if ( auxiliaryConstraints->setupConstraint( i,ST_INACTIVE ) != SUCCESSFUL_RETURN )
2404  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2405  }
2406  }
2407  }
2408 
2409  /* Obtain initial working set in accordance to sign of dual solution vector. */
2410  if ( ( xOpt == 0 ) && ( yOpt != 0 ) )
2411  {
2412  for( i=0; i<nC; ++i )
2413  {
2414  if ( yOpt[nV+i] > EPS )
2415  {
2416  if ( auxiliaryConstraints->setupConstraint( i,ST_LOWER ) != SUCCESSFUL_RETURN )
2417  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2418  continue;
2419  }
2420 
2421  if ( yOpt[nV+i] < -EPS )
2422  {
2423  if ( auxiliaryConstraints->setupConstraint( i,ST_UPPER ) != SUCCESSFUL_RETURN )
2424  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2425  continue;
2426  }
2427 
2428  /* Moreover, add all equality constraints if specified. */
2429  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
2430  if ( constraints.getType( i ) == ST_EQUALITY )
2431  {
2432  if ( auxiliaryConstraints->setupConstraint( i,ST_LOWER ) != SUCCESSFUL_RETURN )
2433  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2434  }
2435  else
2436  #endif
2437  {
2438  if ( auxiliaryConstraints->setupConstraint( i,ST_INACTIVE ) != SUCCESSFUL_RETURN )
2439  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2440  }
2441  }
2442  }
2443 
2444  /* If xOpt and yOpt are null pointer and no initial working is specified,
2445  * start with empty working set (or implicitly fixed bounds and equality constraints only)
2446  * for auxiliary QP. */
2447  if ( ( xOpt == 0 ) && ( yOpt == 0 ) )
2448  {
2449  for( i=0; i<nC; ++i )
2450  {
2451  /* Only add all equality constraints if specified. */
2452  #ifdef __ALWAYS_INITIALISE_WITH_ALL_EQUALITIES__
2453  if ( constraints.getType( i ) == ST_EQUALITY )
2454  {
2455  if ( auxiliaryConstraints->setupConstraint( i,ST_LOWER ) != SUCCESSFUL_RETURN )
2456  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2457  }
2458  else
2459  #endif
2460  {
2461  if ( auxiliaryConstraints->setupConstraint( i,ST_INACTIVE ) != SUCCESSFUL_RETURN )
2462  return THROWERROR( RET_OBTAINING_WORKINGSET_FAILED );
2463  }
2464  }
2465  }
2466  }
2467 
2468  return SUCCESSFUL_RETURN;
2469 }
2470 
2471 
2472 
2473 /*
2474  * s e t u p A u x i l i a r y W o r k i n g S e t
2475  */
2476 returnValue QProblem::setupAuxiliaryWorkingSet( const Bounds* const auxiliaryBounds,
2477  const Constraints* const auxiliaryConstraints,
2478  BooleanType setupAfresh
2479  )
2480 {
2481  int i;
2482  int nV = getNV( );
2483  int nC = getNC( );
2484  BooleanType WSisTrivial = BT_TRUE;
2485 
2486  /* consistency checks */
2487  if ( auxiliaryBounds != 0 )
2488  {
2489  for( i=0; i<nV; ++i )
2490  if ( ( bounds.getStatus( i ) == ST_UNDEFINED ) || ( auxiliaryBounds->getStatus( i ) == ST_UNDEFINED ) )
2491  return THROWERROR( RET_UNKNOWN_BUG );
2492  }
2493  else
2494  {
2495  return THROWERROR( RET_INVALID_ARGUMENTS );
2496  }
2497 
2498  if ( auxiliaryConstraints != 0 )
2499  {
2500  for( i=0; i<nC; ++i )
2501  if ( ( constraints.getStatus( i ) == ST_UNDEFINED ) || ( auxiliaryConstraints->getStatus( i ) == ST_UNDEFINED ) )
2502  return THROWERROR( RET_UNKNOWN_BUG );
2503  }
2504  else
2505  {
2506  return THROWERROR( RET_INVALID_ARGUMENTS );
2507  }
2508 
2509  /* Check for trivial working set (all and only bounds active) */
2510  for (i = 0; i < nV; i++)
2511  if (auxiliaryBounds->getStatus(i) == ST_INACTIVE)
2512  {
2513  WSisTrivial = BT_FALSE;
2514  break;
2515  }
2516  for (i = 0; i < nC; i++)
2517  if (auxiliaryConstraints->getStatus(i) != ST_INACTIVE)
2518  {
2519  WSisTrivial = BT_FALSE;
2520  break;
2521  }
2522 
2523  if (WSisTrivial == BT_TRUE)
2524  {
2525  for (i = 0; i < nV; i++)
2526  bounds.moveFreeToFixed(i, auxiliaryBounds->getStatus(i));
2527 
2528  return SUCCESSFUL_RETURN;
2529  }
2530 
2531 
2532  /* I) SETUP CHOLESKY FLAG:
2533  * Cholesky decomposition shall only be updated if working set
2534  * shall be updated (i.e. NOT setup afresh!) */
2535  BooleanType updateCholesky;
2536  if ( setupAfresh == BT_TRUE )
2537  updateCholesky = BT_FALSE;
2538  else
2539  updateCholesky = BT_TRUE;
2540 
2541 
2542  BooleanType was_fulli = options.enableFullLITests;
2543  real_t backupEpsLITests = options.epsLITests;
2544 
2545  options.enableFullLITests = BT_FALSE;
2546  /* options.epsLITests = 1e-1; */
2547 
2548  /* II) REMOVE FORMERLY ACTIVE (CONSTRAINTS') BOUNDS (IF NECESSARY): */
2549  if ( setupAfresh == BT_FALSE )
2550  {
2551  /* 1) Remove all active constraints that shall be inactive or disabled AND
2552  * all active constraints that are active at the wrong bound. */
2553  for( i=0; i<nC; ++i )
2554  {
2555  if ( ( constraints.getStatus( i ) == ST_LOWER ) && ( auxiliaryConstraints->getStatus( i ) != ST_LOWER ) )
2556  if ( removeConstraint( i,updateCholesky,BT_FALSE,options.enableNZCTests ) != SUCCESSFUL_RETURN )
2557  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2558 
2559  if ( ( constraints.getStatus( i ) == ST_UPPER ) && ( auxiliaryConstraints->getStatus( i ) != ST_UPPER ) )
2560  if ( removeConstraint( i,updateCholesky,BT_FALSE,options.enableNZCTests ) != SUCCESSFUL_RETURN )
2561  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2562  }
2563 
2564  /* 2) Remove all active bounds that shall be inactive AND
2565  * all active bounds that are active at the wrong bound. */
2566  for( i=0; i<nV; ++i )
2567  {
2568  if ( ( bounds.getStatus( i ) == ST_LOWER ) && ( auxiliaryBounds->getStatus( i ) != ST_LOWER ) )
2569  if ( removeBound( i,updateCholesky,BT_FALSE,options.enableNZCTests ) != SUCCESSFUL_RETURN )
2570  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2571 
2572  if ( ( bounds.getStatus( i ) == ST_UPPER ) && ( auxiliaryBounds->getStatus( i ) != ST_UPPER ) )
2573  if ( removeBound( i,updateCholesky,BT_FALSE,options.enableNZCTests ) != SUCCESSFUL_RETURN )
2574  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2575  }
2576  }
2577 
2578 
2579  /* III) ADD NEWLY ACTIVE (CONSTRAINTS') BOUNDS: */
2580 
2581  /* 1) Add all equality bounds. */
2582  for( i=0; i<nV; ++i )
2583  {
2584  if ( ( bounds.getType( i ) == ST_EQUALITY ) && ( ( bounds.getStatus( i ) == ST_INACTIVE ) && ( auxiliaryBounds->getStatus( i ) != ST_INACTIVE ) ) )
2585  {
2586  /* No check for linear independence necessary. */
2587  if ( addBound( i,auxiliaryBounds->getStatus( i ),updateCholesky ) != SUCCESSFUL_RETURN )
2588  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2589  }
2590  }
2591 
2592  /* 2) Add all equality constraints. */
2593  for( i=0; i<nC; ++i )
2594  {
2595  if ( ( constraints.getType( i ) == ST_EQUALITY ) && ( ( constraints.getStatus( i ) == ST_INACTIVE ) && ( auxiliaryConstraints->getStatus( i ) != ST_INACTIVE ) ) )
2596  {
2597  /* Add constraint only if it is linearly independent from the current working set. */
2598  if ( addConstraint_checkLI( i ) == RET_LINEARLY_INDEPENDENT )
2599  {
2600  if ( addConstraint( i,auxiliaryConstraints->getStatus( i ),updateCholesky ) != SUCCESSFUL_RETURN )
2601  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2602  }
2603  else
2604  {
2605  /* Equalities are not linearly independent! */
2606  constraints.setType(i, ST_BOUNDED);
2607  }
2608  }
2609  }
2610 
2611 
2612  /* 3) Add all inactive bounds that shall be active AND
2613  * all formerly active bounds that have been active at the wrong bound. */
2614  for( i=0; i<nV; ++i )
2615  {
2616  if ( ( bounds.getType( i ) != ST_EQUALITY ) && ( ( bounds.getStatus( i ) == ST_INACTIVE ) && ( auxiliaryBounds->getStatus( i ) != ST_INACTIVE ) ) )
2617  {
2618  /* Add bound only if it is linearly independent from the current working set. */
2619  if ( addBound_checkLI( i ) == RET_LINEARLY_INDEPENDENT )
2620  {
2621  if ( addBound( i,auxiliaryBounds->getStatus( i ),updateCholesky ) != SUCCESSFUL_RETURN )
2622  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2623  }
2624  }
2625  }
2626 
2627  /* 4) Add all inactive constraints that shall be active AND
2628  * all formerly active constraints that have been active at the wrong bound. */
2629  for( i=0; i<nC; ++i )
2630  {
2631  if ( ( constraints.getType( i ) != ST_EQUALITY ) && ( auxiliaryConstraints->getStatus( i ) != ST_INACTIVE ) )
2632  {
2633  /* formerly inactive */
2634  if ( constraints.getStatus( i ) == ST_INACTIVE )
2635  {
2636  /* Add constraint only if it is linearly independent from the current working set. */
2637  if ( addConstraint_checkLI( i ) == RET_LINEARLY_INDEPENDENT )
2638  {
2639  if ( addConstraint( i,auxiliaryConstraints->getStatus( i ),updateCholesky ) != SUCCESSFUL_RETURN )
2640  return THROWERROR( RET_SETUP_WORKINGSET_FAILED );
2641  }
2642  }
2643  }
2644  }
2645 
2646  options.enableFullLITests = was_fulli;
2647  options.epsLITests = backupEpsLITests;
2648 
2649  return SUCCESSFUL_RETURN;
2650 }
2651 
2652 
2653 /*
2654  * s e t u p A u x i l i a r y Q P s o l u t i o n
2655  */
2656 returnValue QProblem::setupAuxiliaryQPsolution( const real_t* const xOpt, const real_t* const yOpt
2657  )
2658 {
2659  int i, j;
2660  int nV = getNV( );
2661  int nC = getNC( );
2662 
2663 
2664  /* Setup primal/dual solution vector for auxiliary initial QP:
2665  * if a null pointer is passed, a zero vector is assigned;
2666  * old solution vector is kept if pointer to internal solution vevtor is passed. */
2667  if ( xOpt != 0 )
2668  {
2669  if ( xOpt != x )
2670  for( i=0; i<nV; ++i )
2671  x[i] = xOpt[i];
2672 
2673  A->times(1, 1.0, x, nV, 0.0, Ax, nC);
2674 
2675  for ( j=0; j<nC; ++j )
2676  {
2677  Ax_l[j] = Ax[j];
2678  Ax_u[j] = Ax[j];
2679  }
2680  }
2681  else
2682  {
2683  for( i=0; i<nV; ++i )
2684  x[i] = 0.0;
2685 
2686  for ( j=0; j<nC; ++j )
2687  {
2688  Ax[j] = 0.0;
2689  Ax_l[j] = 0.0;
2690  Ax_u[j] = 0.0;
2691  }
2692  }
2693 
2694  if ( yOpt != 0 )
2695  {
2696  if ( yOpt != y )
2697  for( i=0; i<nV+nC; ++i )
2698  y[i] = yOpt[i];
2699  }
2700  else
2701  {
2702  for( i=0; i<nV+nC; ++i )
2703  y[i] = 0.0;
2704  }
2705 
2706  return SUCCESSFUL_RETURN;
2707 }
2708 
2709 
2710 /*
2711  * s e t u p A u x i l i a r y Q P g r a d i e n t
2712  */
2713 returnValue QProblem::setupAuxiliaryQPgradient( )
2714 {
2715  int i;
2716  int nV = getNV( );
2717  int nC = getNC( );
2718 
2719 
2720  /* Setup gradient vector: g = -H*x + [Id A]'*[yB yC]
2721  * = yB - H*x + A'*yC. */
2722  switch ( hessianType )
2723  {
2724  case HST_ZERO:
2725  if ( usingRegularisation( ) == BT_FALSE )
2726  for ( i=0; i<nV; ++i )
2727  g[i] = y[i];
2728  else
2729  for ( i=0; i<nV; ++i )
2730  g[i] = y[i] - options.epsRegularisation*x[i];
2731  break;
2732 
2733  case HST_IDENTITY:
2734  for ( i=0; i<nV; ++i )
2735  g[i] = y[i] - x[i];
2736  break;
2737 
2738  default:
2739  /* y'*Id */
2740  for ( i=0; i<nV; ++i )
2741  g[i] = y[i];
2742 
2743  /* - H*x */
2744  H->times(1, -1.0, x, nV, 1.0, g, nV);
2745  break;
2746  }
2747 
2748  /* + A'*yC */
2749  A->transTimes(1, 1.0, y + nV, nC, 1.0, g, nV);
2750 
2751  return SUCCESSFUL_RETURN;
2752 }
2753 
2754 
2755 /*
2756  * s e t u p A u x i l i a r y Q P b o u n d s
2757  */
2758 returnValue QProblem::setupAuxiliaryQPbounds( const Bounds* const auxiliaryBounds,
2759  const Constraints* const auxiliaryConstraints,
2760  BooleanType useRelaxation
2761  )
2762 {
2763  int i;
2764  int nV = getNV( );
2765  int nC = getNC( );
2766 
2767 
2768  /* 1) Setup bound vectors. */
2769  for ( i=0; i<nV; ++i )
2770  {
2771  switch ( bounds.getStatus( i ) )
2772  {
2773  case ST_INACTIVE:
2774  if ( useRelaxation == BT_TRUE )
2775  {
2776  if ( bounds.getType( i ) == ST_EQUALITY )
2777  {
2778  lb[i] = x[i];
2779  ub[i] = x[i];
2780  }
2781  else
2782  {
2783  /* If a bound is inactive although it was supposed to be
2784  * active by the auxiliaryBounds, it could not be added
2785  * due to linear dependence. Thus set it "strongly inactive". */
2786  if ( auxiliaryBounds->getStatus( i ) == ST_LOWER )
2787  lb[i] = x[i];
2788  else
2789  lb[i] = x[i] - options.boundRelaxation;
2790 
2791  if ( auxiliaryBounds->getStatus( i ) == ST_UPPER )
2792  ub[i] = x[i];
2793  else
2794  ub[i] = x[i] + options.boundRelaxation;
2795  }
2796  }
2797  break;
2798 
2799  case ST_LOWER:
2800  lb[i] = x[i];
2801  if ( bounds.getType( i ) == ST_EQUALITY )
2802  {
2803  ub[i] = x[i];
2804  }
2805  else
2806  {
2807  if ( useRelaxation == BT_TRUE )
2808  ub[i] = x[i] + options.boundRelaxation;
2809  }
2810  break;
2811 
2812  case ST_UPPER:
2813  ub[i] = x[i];
2814  if ( bounds.getType( i ) == ST_EQUALITY )
2815  {
2816  lb[i] = x[i];
2817  }
2818  else
2819  {
2820  if ( useRelaxation == BT_TRUE )
2821  lb[i] = x[i] - options.boundRelaxation;
2822  }
2823  break;
2824 
2825  default:
2826  return THROWERROR( RET_UNKNOWN_BUG );
2827  }
2828  }
2829 
2830  /* 2) Setup constraints vectors. */
2831  for ( i=0; i<nC; ++i )
2832  {
2833  switch ( constraints.getStatus( i ) )
2834  {
2835  case ST_INACTIVE:
2836  if ( useRelaxation == BT_TRUE )
2837  {
2838  if ( constraints.getType( i ) == ST_EQUALITY )
2839  {
2840  lbA[i] = Ax_l[i];
2841  ubA[i] = Ax_u[i];
2842  }
2843  else
2844  {
2845  /* If a constraint is inactive although it was supposed to be
2846  * active by the auxiliaryConstraints, it could not be added
2847  * due to linear dependence. Thus set it "strongly inactive". */
2848  if ( auxiliaryConstraints->getStatus( i ) == ST_LOWER )
2849  lbA[i] = Ax_l[i];
2850  else
2851  lbA[i] = Ax_l[i] - options.boundRelaxation;
2852 
2853  if ( auxiliaryConstraints->getStatus( i ) == ST_UPPER )
2854  ubA[i] = Ax_u[i];
2855  else
2856  ubA[i] = Ax_u[i] + options.boundRelaxation;
2857  }
2858  }
2859  break;
2860 
2861  case ST_LOWER:
2862  lbA[i] = Ax_l[i];
2863  if ( constraints.getType( i ) == ST_EQUALITY )
2864  {
2865  ubA[i] = Ax_l[i];
2866  }
2867  else
2868  {
2869  if ( useRelaxation == BT_TRUE )
2870  ubA[i] = Ax_l[i] + options.boundRelaxation;
2871  }
2872  break;
2873 
2874  case ST_UPPER:
2875  ubA[i] = Ax_u[i];
2876  if ( constraints.getType( i ) == ST_EQUALITY )
2877  {
2878  lbA[i] = Ax_u[i];
2879  }
2880  else
2881  {
2882  if ( useRelaxation == BT_TRUE )
2883  lbA[i] = Ax_u[i] - options.boundRelaxation;
2884  }
2885  break;
2886 
2887  default:
2888  return THROWERROR( RET_UNKNOWN_BUG );
2889  }
2890  Ax_l[i] = Ax_l[i] - lbA[i];
2891  Ax_u[i] = ubA[i] - Ax_u[i];
2892  }
2893 
2894  return SUCCESSFUL_RETURN;
2895 }
2896 
2897 
2898 /*
2899  * a d d C o n s t r a i n t
2900  */
2901 returnValue QProblem::addConstraint( int number, SubjectToStatus C_status,
2902  BooleanType updateCholesky,
2903  BooleanType ensureLI
2904  )
2905 {
2906  int i, j, ii;
2907 
2908  /* consistency checks */
2909  if ( constraints.getStatus( number ) != ST_INACTIVE )
2910  return THROWERROR( RET_CONSTRAINT_ALREADY_ACTIVE );
2911 
2912  if ( ( constraints.getNC( ) - getNAC( ) ) == constraints.getNUC( ) )
2913  return THROWERROR( RET_ALL_CONSTRAINTS_ACTIVE );
2914 
2915  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
2916  ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
2917  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
2918  ( getStatus( ) == QPS_SOLVED ) )
2919  {
2920  return THROWERROR( RET_UNKNOWN_BUG );
2921  }
2922 
2923 
2924  /* I) ENSURE LINEAR INDEPENDENCE OF THE WORKING SET,
2925  * i.e. remove a constraint or bound if linear dependence occurs. */
2926  /* check for LI only if Cholesky decomposition shall be updated! */
2927  if ( updateCholesky == BT_TRUE && ensureLI == BT_TRUE )
2928  {
2929  returnValue ensureLIreturnvalue = addConstraint_ensureLI( number,C_status );
2930 
2931  switch ( ensureLIreturnvalue )
2932  {
2933  case SUCCESSFUL_RETURN:
2934  break;
2935 
2936  case RET_LI_RESOLVED:
2937  break;
2938 
2939  case RET_ENSURELI_FAILED_NOINDEX:
2940  return RET_ADDCONSTRAINT_FAILED_INFEASIBILITY;
2941 
2942  case RET_ENSURELI_FAILED_CYCLING:
2943  return RET_ADDCONSTRAINT_FAILED_INFEASIBILITY;
2944 
2945  default:
2946  return THROWERROR( RET_ENSURELI_FAILED );
2947  }
2948  }
2949 
2950  /* some definitions */
2951  int nV = getNV( );
2952  int nFR = getNFR( );
2953  int nAC = getNAC( );
2954  int nZ = getNZ( );
2955 
2956  int tcol = sizeT - nAC;
2957 
2958 
2959  int* FR_idx;
2960  bounds.getFree( )->getNumberArray( &FR_idx );
2961 
2962  real_t* aFR = new real_t[nFR];
2963  real_t* wZ = new real_t[nZ];
2964  for( i=0; i<nZ; ++i )
2965  wZ[i] = 0.0;
2966 
2967 
2968  /* II) ADD NEW ACTIVE CONSTRAINT TO MATRIX T: */
2969  /* 1) Add row [wZ wY] = aFR'*[Z Y] to the end of T: assign aFR. */
2970  A->getRow(number, bounds.getFree(), 1.0, aFR);
2971 
2972  /* calculate wZ */
2973  for( i=0; i<nFR; ++i )
2974  {
2975  ii = FR_idx[i];
2976  for( j=0; j<nZ; ++j )
2977  wZ[j] += aFR[i] * QQ(ii,j);
2978  }
2979 
2980  /* 2) Calculate wY and store it directly into T. */
2981  if ( nAC > 0 )
2982  {
2983  for( j=0; j<nAC; ++j )
2984  TT(nAC,tcol+j) = 0.0;
2985  for( i=0; i<nFR; ++i )
2986  {
2987  ii = FR_idx[i];
2988  for( j=0; j<nAC; ++j )
2989  TT(nAC,tcol+j) += aFR[i] * QQ(ii,nZ+j);
2990  }
2991  }
2992 
2993  delete[] aFR;
2994 
2995 
2996  real_t c, s, nu;
2997 
2998  if ( nZ > 0 )
2999  {
3000  /* II) RESTORE TRIANGULAR FORM OF T: */
3001  /* Use column-wise Givens rotations to restore reverse triangular form
3002  * of T, simultanenous change of Q (i.e. Z) and R. */
3003  for( j=0; j<nZ-1; ++j )
3004  {
3005  computeGivens( wZ[j+1],wZ[j], wZ[j+1],wZ[j],c,s );
3006  nu = s/(1.0+c);
3007 
3008  for( i=0; i<nFR; ++i )
3009  {
3010  ii = FR_idx[i];
3011  applyGivens( c,s,nu,QQ(ii,1+j),QQ(ii,j), QQ(ii,1+j),QQ(ii,j) );
3012  }
3013 
3014  if ( ( updateCholesky == BT_TRUE ) &&
3015  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3016  {
3017  for( i=0; i<=j+1; ++i )
3018  applyGivens( c,s,nu,RR(i,1+j),RR(i,j), RR(i,1+j),RR(i,j) );
3019  }
3020  }
3021 
3022  TT(nAC,tcol-1) = wZ[nZ-1];
3023 
3024 
3025  if ( ( updateCholesky == BT_TRUE ) &&
3026  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3027  {
3028  /* III) RESTORE TRIANGULAR FORM OF R:
3029  * Use row-wise Givens rotations to restore upper triangular form of R. */
3030  for( i=0; i<nZ-1; ++i )
3031  {
3032  computeGivens( RR(i,i),RR(1+i,i), RR(i,i),RR(1+i,i),c,s );
3033  nu = s/(1.0+c);
3034 
3035  for( j=(1+i); j<(nZ-1); ++j ) /* last column of R is thrown away */
3036  applyGivens( c,s,nu,RR(i,j),RR(1+i,j), RR(i,j),RR(1+i,j) );
3037  }
3038  /* last column of R is thrown away */
3039  for( i=0; i<nZ; ++i )
3040  RR(i,nZ-1) = 0.0;
3041  }
3042  }
3043 
3044  delete[] wZ;
3045 
3046 
3047  /* IV) UPDATE INDICES */
3048  idxAddC = number;
3049  if ( constraints.moveInactiveToActive( number,C_status ) != SUCCESSFUL_RETURN )
3050  return THROWERROR( RET_ADDCONSTRAINT_FAILED );
3051 
3052 
3053  return SUCCESSFUL_RETURN;
3054 }
3055 
3056 
3057 
3058 /*
3059  * a d d C o n s t r a i n t _ c h e c k L I
3060  */
3061 returnValue QProblem::addConstraint_checkLI( int number )
3062 {
3063  returnValue returnvalue = RET_LINEARLY_DEPENDENT;
3064 
3065  int i, j, ii;
3066  int nV = getNV( );
3067  int nFR = getNFR( );
3068  int nZ = getNZ( );
3069  int nC = getNC( );
3070  int nAC = getNAC();
3071  int nFX = getNFX();
3072  int *FR_idx;
3073 
3074  bounds.getFree( )->getNumberArray( &FR_idx );
3075 
3076 
3077  if (options.enableFullLITests)
3078  {
3079  /*
3080  * expensive LI test. Backsolve with refinement using special right
3081  * hand side. This gives an estimate for what should be considered
3082  * "zero". We then check linear independence relative to this estimate.
3083  */
3084 
3085  int *FX_idx, *AC_idx, *IAC_idx;
3086 
3087  real_t *delta_g = new real_t[nV];
3088  real_t *delta_xFX = new real_t[nFX];
3089  real_t *delta_xFR = new real_t[nFR];
3090  real_t *delta_yAC = new real_t[nAC];
3091  real_t *delta_yFX = new real_t[nFX];
3092 
3093  bounds.getFixed( )->getNumberArray( &FX_idx );
3094  constraints.getActive( )->getNumberArray( &AC_idx );
3095  constraints.getInactive( )->getNumberArray( &IAC_idx );
3096 
3097  int dim = (nC>nV)?nC:nV;
3098  real_t *nul = new real_t[dim];
3099  for (ii = 0; ii < dim; ++ii)
3100  nul[ii]=0.0;
3101 
3102  A->getRow (number, 0, 1.0, delta_g);
3103 
3104  returnvalue = determineStepDirection ( delta_g,
3105  nul, nul, nul, nul,
3106  BT_FALSE, BT_FALSE,
3107  delta_xFX, delta_xFR, delta_yAC, delta_yFX);
3108  delete[] nul;
3109 
3110  /* compute the weight in inf-norm */
3111  real_t weight = 0.0;
3112  for (ii = 0; ii < nAC; ++ii)
3113  {
3114  real_t a = fabs (delta_yAC[ii]);
3115  if (weight < a) weight = a;
3116  }
3117  for (ii = 0; ii < nFX; ++ii)
3118  {
3119  real_t a = fabs (delta_yFX[ii]);
3120  if (weight < a) weight = a;
3121  }
3122 
3123  /* look at the "zero" in a relative inf-norm */
3124  real_t zero = 0.0;
3125  for (ii = 0; ii < nFX; ++ii)
3126  {
3127  real_t a = fabs (delta_xFX[ii]);
3128  if (zero < a) zero = a;
3129  }
3130  for (ii = 0; ii < nFR; ++ii)
3131  {
3132  real_t a = fabs (delta_xFR[ii]);
3133  if (zero < a) zero = a;
3134  }
3135 
3136  /* relative test against zero in inf-norm */
3137  if (zero > options.epsLITests * weight)
3138  returnvalue = RET_LINEARLY_INDEPENDENT;
3139 
3140  delete[] delta_yFX;
3141  delete[] delta_yAC;
3142  delete[] delta_xFR;
3143  delete[] delta_xFX;
3144  delete[] delta_g;
3145 
3146  }
3147  else
3148  {
3149  /*
3150  * cheap LI test for constraint. Check if constraint <number> is
3151  * linearly independent from the the active ones (<=> is element of null
3152  * space of Afr).
3153  */
3154 
3155  real_t *Arow = new real_t[nFR];
3156  A->getRow(number, bounds.getFree(), 1.0, Arow);
3157 
3158  real_t sum, l2;
3159 
3160  l2 = 0.0;
3161  for (i = 0; i < nFR; i++)
3162  l2 += Arow[i]*Arow[i];
3163 
3164  for( j=0; j<nZ; ++j )
3165  {
3166  sum = 0.0;
3167  for( i=0; i<nFR; ++i )
3168  {
3169  ii = FR_idx[i];
3170  sum += Arow[i] * QQ(ii,j);
3171  }
3172 
3173  if ( fabs( sum ) > options.epsLITests*l2 )
3174  {
3175  /*fprintf(stderr, "LI test: |sum| = %9.2e, l2 = %9.2e, var = %d\n", fabs(sum), l2, jj+1); */
3176  returnvalue = RET_LINEARLY_INDEPENDENT;
3177  break;
3178  }
3179  }
3180 
3181  delete[] Arow;
3182  }
3183 
3184  return THROWINFO( returnvalue );
3185 }
3186 
3187 
3188 /*
3189  * a d d C o n s t r a i n t _ e n s u r e L I
3190  */
3191 returnValue QProblem::addConstraint_ensureLI( int number, SubjectToStatus C_status )
3192 {
3193  int i, j, ii, jj;
3194  int nV = getNV( );
3195  int nFR = getNFR( );
3196  int nFX = getNFX( );
3197  int nAC = getNAC( );
3198  int nZ = getNZ( );
3199 
3200 
3201  /* I) Check if new constraint is linearly independent from the active ones. */
3202  returnValue returnvalueCheckLI = addConstraint_checkLI( number );
3203 
3204  if ( returnvalueCheckLI == RET_INDEXLIST_CORRUPTED )
3205  return THROWERROR( RET_ENSURELI_FAILED );
3206 
3207  if ( returnvalueCheckLI == RET_LINEARLY_INDEPENDENT )
3208  return SUCCESSFUL_RETURN;
3209 
3210 
3211  /* II) NEW CONSTRAINT IS LINEARLY DEPENDENT: */
3212  /* 1) Determine coefficients of linear combination,
3213  * cf. M.J. Best. Applied Mathematics and Parallel Computing, chapter:
3214  * An Algorithm for the Solution of the Parametric Quadratic Programming
3215  * Problem, pages 57-76. Physica-Verlag, Heidelberg, 1996. */
3216  int* FR_idx;
3217  bounds.getFree( )->getNumberArray( &FR_idx );
3218 
3219  int* FX_idx;
3220  bounds.getFixed( )->getNumberArray( &FX_idx );
3221 
3222  real_t* xiC = new real_t[nAC];
3223  real_t* xiC_TMP = new real_t[nAC];
3224  real_t* xiB = new real_t[nFX];
3225  real_t* Arow = new real_t[nFR];
3226 
3227  A->getRow(number, bounds.getFree(), C_status == ST_LOWER ? 1.0 : -1.0, Arow);
3228 
3229  /* 2) Calculate xiC */
3230  if ( nAC > 0 )
3231  {
3232  for( i=0; i<nAC; ++i )
3233  {
3234  xiC_TMP[i] = 0.0;
3235  for( j=0; j<nFR; ++j )
3236  {
3237  jj = FR_idx[j];
3238  xiC_TMP[i] += QQ(jj,nZ+i) * Arow[j];
3239  }
3240  }
3241 
3242  if ( backsolveT( xiC_TMP, BT_TRUE, xiC ) != SUCCESSFUL_RETURN )
3243  {
3244  delete[] Arow;
3245  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3246  return THROWERROR( RET_ENSURELI_FAILED_TQ );
3247  }
3248  }
3249 
3250  /* 3) Calculate xiB. */
3251  int* AC_idx;
3252  constraints.getActive( )->getNumberArray( &AC_idx );
3253 
3254  A->getRow(number, bounds.getFixed(), C_status == ST_LOWER ? 1.0 : -1.0, xiB);
3255  A->transTimes(constraints.getActive(), bounds.getFixed(), 1, -1.0, xiC, nAC, 1.0, xiB, nFX);
3256 
3257  /* III) DETERMINE CONSTRAINT/BOUND TO BE REMOVED. */
3258  real_t y_min = options.maxDualJump;
3259  int y_min_number = -1;
3260  int y_min_number_bound = -1;
3261  BooleanType y_min_isBound = BT_FALSE;
3262 
3263  real_t* num = new real_t[nV];
3264 
3265  /* 1) Constraints. */
3266  for( i=0; i<nAC; ++i )
3267  {
3268  ii = AC_idx[i];
3269  num[i] = y[nV+ii];
3270  }
3271 
3272  performRatioTest( nAC,AC_idx,&constraints, num,xiC, options.epsDen,options.epsDen, y_min,y_min_number );
3273 
3274 
3275  /* 2) Bounds. */
3276  for( i=0; i<nFX; ++i )
3277  {
3278  ii = FX_idx[i];
3279  num[i] = y[ii];
3280  }
3281 
3282  performRatioTest( nFX,FX_idx,&bounds, num,xiB, options.epsDen,options.epsDen, y_min,y_min_number_bound );
3283 
3284  if ( y_min_number_bound >= 0 )
3285  {
3286  y_min_number = y_min_number_bound;
3287  y_min_isBound = BT_TRUE;
3288  }
3289 
3290  delete[] num;
3291 
3292  /* setup output preferences */
3293  char messageString[80];
3294 
3295  /* IV) REMOVE CONSTRAINT/BOUND FOR RESOLVING LINEAR DEPENDENCE: */
3296  if ( y_min_number >= 0 )
3297  {
3298  /* Update Lagrange multiplier... */
3299  for( i=0; i<nAC; ++i )
3300  {
3301  ii = AC_idx[i];
3302  y[nV+ii] -= y_min * xiC[i];
3303  }
3304  for( i=0; i<nFX; ++i )
3305  {
3306  ii = FX_idx[i];
3307  y[ii] -= y_min * xiB[i];
3308  }
3309 
3310  /* ... also for newly active constraint... */
3311  if ( C_status == ST_LOWER )
3312  y[nV+number] = y_min;
3313  else
3314  y[nV+number] = -y_min;
3315 
3316  /* ... and for constraint to be removed. */
3317  if ( y_min_isBound == BT_TRUE )
3318  {
3319  #ifndef __XPCTARGET__
3320  snprintf( messageString,80,"bound no. %d.",y_min_number );
3321  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3322  #endif
3323 
3324  if ( removeBound( y_min_number,BT_TRUE,BT_FALSE,BT_FALSE ) != SUCCESSFUL_RETURN )
3325  {
3326  delete[] Arow;
3327  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3328 
3329  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
3330  }
3331  y[y_min_number] = 0.0;
3332  }
3333  else
3334  {
3335  #ifndef __XPCTARGET__
3336  snprintf( messageString,80,"constraint no. %d.",y_min_number );
3337  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3338  #endif
3339 
3340  if ( removeConstraint( y_min_number,BT_TRUE,BT_FALSE,BT_FALSE ) != SUCCESSFUL_RETURN )
3341  {
3342  delete[] Arow;
3343  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3344 
3345  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
3346  }
3347 
3348  y[nV+y_min_number] = 0.0;
3349  }
3350  }
3351  else
3352  {
3353  /* no constraint/bound can be removed => QP is infeasible! */
3354  delete[] Arow;
3355  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3356 
3357  return setInfeasibilityFlag( RET_ENSURELI_FAILED_NOINDEX );
3358  }
3359 
3360  delete[] Arow;
3361  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3362 
3363  return getGlobalMessageHandler( )->throwInfo( RET_LI_RESOLVED,0,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3364 }
3365 
3366 
3367 
3368 /*
3369  * a d d B o u n d
3370  */
3371 returnValue QProblem::addBound( int number, SubjectToStatus B_status,
3372  BooleanType updateCholesky,
3373  BooleanType ensureLI
3374  )
3375 {
3376  int i, j, ii;
3377 
3378  /* consistency checks */
3379  if ( bounds.getStatus( number ) != ST_INACTIVE )
3380  return THROWERROR( RET_BOUND_ALREADY_ACTIVE );
3381 
3382  if ( getNFR( ) == bounds.getNUV( ) )
3383  return THROWERROR( RET_ALL_BOUNDS_ACTIVE );
3384 
3385  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
3386  ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
3387  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
3388  ( getStatus( ) == QPS_SOLVED ) )
3389  {
3390  return THROWERROR( RET_UNKNOWN_BUG );
3391  }
3392 
3393 
3394  /* I) ENSURE LINEAR INDEPENDENCE OF THE WORKING SET,
3395  * i.e. remove a constraint or bound if linear dependence occurs. */
3396  /* check for LI only if Cholesky decomposition shall be updated! */
3397  if ( updateCholesky == BT_TRUE && ensureLI == BT_TRUE )
3398  {
3399  returnValue ensureLIreturnvalue = addBound_ensureLI( number,B_status );
3400 
3401  switch ( ensureLIreturnvalue )
3402  {
3403  case SUCCESSFUL_RETURN:
3404  break;
3405 
3406  case RET_LI_RESOLVED:
3407  break;
3408 
3409  case RET_ENSURELI_FAILED_NOINDEX:
3410  return RET_ADDBOUND_FAILED_INFEASIBILITY;
3411 
3412  case RET_ENSURELI_FAILED_CYCLING:
3413  return RET_ADDBOUND_FAILED_INFEASIBILITY;
3414 
3415  default:
3416  return THROWERROR( RET_ENSURELI_FAILED );
3417  }
3418  }
3419 
3420 
3421  /* some definitions */
3422  int nV = getNV( );
3423  int nFR = getNFR( );
3424  int nAC = getNAC( );
3425  int nZ = getNZ( );
3426 
3427  int tcol = sizeT - nAC;
3428 
3429 
3430  /* II) SWAP INDEXLIST OF FREE VARIABLES:
3431  * move the variable to be fixed to the end of the list of free variables. */
3432  int lastfreenumber = bounds.getFree( )->getLastNumber( );
3433  if ( lastfreenumber != number )
3434  if ( bounds.swapFree( number,lastfreenumber ) != SUCCESSFUL_RETURN )
3435  THROWERROR( RET_ADDBOUND_FAILED );
3436 
3437 
3438  int* FR_idx;
3439  bounds.getFree( )->getNumberArray( &FR_idx );
3440 
3441  real_t* w = new real_t[nFR];
3442 
3443 
3444  /* III) ADD NEW ACTIVE BOUND TO TOP OF MATRIX T: */
3445  /* 1) add row [wZ wY] = [Z Y](number) at the top of T: assign w */
3446  for( i=0; i<nFR; ++i )
3447  w[i] = QQ(FR_idx[nFR-1],i);
3448 
3449 
3450  /* 2) Use column-wise Givens rotations to restore reverse triangular form
3451  * of the first row of T, simultanenous change of Q (i.e. Z) and R. */
3452  real_t c, s, nu;
3453 
3454  for( j=0; j<nZ-1; ++j )
3455  {
3456  computeGivens( w[j+1],w[j], w[j+1],w[j],c,s );
3457  nu = s/(1.0+c);
3458 
3459  for( i=0; i<nFR; ++i )
3460  {
3461  ii = FR_idx[i];
3462  applyGivens( c,s,nu,QQ(ii,1+j),QQ(ii,j), QQ(ii,1+j),QQ(ii,j) );
3463  }
3464 
3465  if ( ( updateCholesky == BT_TRUE ) &&
3466  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3467  {
3468  for( i=0; i<=j+1; ++i )
3469  applyGivens( c,s,nu,RR(i,1+j),RR(i,j), RR(i,1+j),RR(i,j) );
3470  }
3471  }
3472 
3473 
3474  if ( nAC > 0 ) /* ( nAC == 0 ) <=> ( nZ == nFR ) <=> Y and T are empty => nothing to do */
3475  {
3476  /* store new column a in a temporary vector instead of shifting T one column to the left */
3477  real_t* tmp = new real_t[nAC];
3478  for( i=0; i<nAC; ++i )
3479  tmp[i] = 0.0;
3480 
3481  {
3482  j = nZ-1;
3483 
3484  computeGivens( w[j+1],w[j], w[j+1],w[j],c,s );
3485  nu = s/(1.0+c);
3486 
3487  for( i=0; i<nFR; ++i )
3488  {
3489  ii = FR_idx[i];
3490  applyGivens( c,s,nu,QQ(ii,1+j),QQ(ii,j), QQ(ii,1+j),QQ(ii,j) );
3491  }
3492 
3493  applyGivens( c,s,nu,TT(nAC-1,tcol),tmp[nAC-1], tmp[nAC-1],TT(nAC-1,tcol) );
3494  }
3495 
3496  for( j=nZ; j<nFR-1; ++j )
3497  {
3498  computeGivens( w[j+1],w[j], w[j+1],w[j],c,s );
3499  nu = s/(1.0+c);
3500 
3501  for( i=0; i<nFR; ++i )
3502  {
3503  ii = FR_idx[i];
3504  applyGivens( c,s,nu,QQ(ii,1+j),QQ(ii,j), QQ(ii,1+j),QQ(ii,j) );
3505  }
3506 
3507  for( i=(nFR-2-j); i<nAC; ++i )
3508  applyGivens( c,s,nu,TT(i,1+tcol-nZ+j),tmp[i], tmp[i],TT(i,1+tcol-nZ+j) );
3509  }
3510 
3511  delete[] tmp;
3512  }
3513 
3514  delete[] w;
3515 
3516 
3517  if ( ( updateCholesky == BT_TRUE ) &&
3518  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3519  {
3520  /* IV) RESTORE TRIANGULAR FORM OF R:
3521  * use row-wise Givens rotations to restore upper triangular form of R */
3522  for( i=0; i<nZ-1; ++i )
3523  {
3524  computeGivens( RR(i,i),RR(1+i,i), RR(i,i),RR(1+i,i),c,s );
3525  nu = s/(1.0+c);
3526 
3527  for( j=(1+i); j<nZ-1; ++j ) /* last column of R is thrown away */
3528  applyGivens( c,s,nu,RR(i,j),RR(1+i,j), RR(i,j),RR(1+i,j) );
3529  }
3530  /* last column of R is thrown away */
3531  for( i=0; i<nZ; ++i )
3532  RR(i,nZ-1) = 0.0;
3533  }
3534 
3535 
3536  /* V) UPDATE INDICES */
3537  idxAddB = number;
3538  if ( bounds.moveFreeToFixed( number,B_status ) != SUCCESSFUL_RETURN )
3539  return THROWERROR( RET_ADDBOUND_FAILED );
3540 
3541  return SUCCESSFUL_RETURN;
3542 }
3543 
3544 
3545 /*
3546  * a d d B o u n d _ c h e c k L I
3547  */
3548 returnValue QProblem::addBound_checkLI( int number )
3549 {
3550  int i, ii;
3551  int nV = getNV( ); /* for QQ() macro */
3552  int nFR = getNFR( );
3553  int nAC = getNAC();
3554  int nFX = getNFX();
3555  int nC = getNC( );
3556  returnValue returnvalue = RET_LINEARLY_DEPENDENT;
3557 
3558  if (options.enableFullLITests)
3559  {
3560  /*
3561  * expensive LI test. Backsolve with refinement using special right
3562  * hand side. This gives an estimate for what should be considered
3563  * "zero". We then check linear independence relative to this estimate.
3564  */
3565 
3566  /*
3567  * expensive LI test. Backsolve with refinement using special right
3568  * hand side. This gives an estimate for what should be considered
3569  * "zero". We then check linear independence relative to this estimate.
3570  */
3571 
3572  real_t *delta_g = new real_t[nV];
3573  real_t *delta_xFX = new real_t[nFX];
3574  real_t *delta_xFR = new real_t[nFR];
3575  real_t *delta_yAC = new real_t[nAC];
3576  real_t *delta_yFX = new real_t[nFX];
3577 
3578  for (ii = 0; ii < nV; ++ii)
3579  delta_g[ii] = 0.0;
3580  delta_g[number] = 1.0; /* sign doesn't matter here */
3581 
3582  int dim = (nC>nV)?nC:nV;
3583  real_t *nul = new real_t[dim];
3584  for (ii = 0; ii < dim; ++ii)
3585  nul[ii]=0.0;
3586 
3587  returnvalue = determineStepDirection (delta_g,
3588  nul, nul, nul, nul,
3589  BT_FALSE, BT_FALSE,
3590  delta_xFX, delta_xFR, delta_yAC, delta_yFX);
3591 
3592  delete[] nul;
3593 
3594  /* compute the weight in inf-norm */
3595  real_t weight = 0.0;
3596  for (ii = 0; ii < nAC; ++ii)
3597  {
3598  real_t a = fabs (delta_yAC[ii]);
3599  if (weight < a) weight = a;
3600  }
3601  for (ii = 0; ii < nFX; ++ii)
3602  {
3603  real_t a = fabs (delta_yFX[ii]);
3604  if (weight < a) weight = a;
3605  }
3606 
3607  /* look at the "zero" in a relative inf-norm */
3608  real_t zero = 0.0;
3609  for (ii = 0; ii < nFX; ++ii)
3610  {
3611  real_t a = fabs (delta_xFX[ii]);
3612  if (zero < a) zero = a;
3613  }
3614  for (ii = 0; ii < nFR; ++ii)
3615  {
3616  real_t a = fabs (delta_xFR[ii]);
3617  if (zero < a) zero = a;
3618  }
3619 
3620  /* relative test against zero in inf-norm */
3621  if (zero > options.epsLITests * weight)
3622  returnvalue = RET_LINEARLY_INDEPENDENT;
3623 
3624  delete[] delta_yFX;
3625  delete[] delta_yAC;
3626  delete[] delta_xFR;
3627  delete[] delta_xFX;
3628  delete[] delta_g;
3629 
3630  }
3631  else
3632  {
3633  /*
3634  * cheap LI test for simple bound. Check if constraint <number> is
3635  * linearly independent from the the active ones (<=> is element of null
3636  * space of Afr).
3637  */
3638 
3639  /* some definitions */
3640  int nZ = getNZ( );
3641 
3642  for( i=0; i<nZ; ++i )
3643  if ( fabs( QQ(number,i) ) > options.epsLITests )
3644  {
3645  returnvalue = RET_LINEARLY_INDEPENDENT;
3646  break;
3647  }
3648  }
3649 
3650  return THROWINFO( returnvalue );
3651 }
3652 
3653 
3654 /*
3655  * a d d B o u n d _ e n s u r e L I
3656  */
3657 returnValue QProblem::addBound_ensureLI( int number, SubjectToStatus B_status )
3658 {
3659  int i, ii;
3660  int nV = getNV( );
3661  int nFX = getNFX( );
3662  int nAC = getNAC( );
3663  int nZ = getNZ( );
3664 
3665 
3666  /* I) Check if new constraint is linearly independent from the active ones. */
3667  returnValue returnvalueCheckLI = addBound_checkLI( number );
3668 
3669  if ( returnvalueCheckLI == RET_INDEXLIST_CORRUPTED )
3670  return THROWERROR( RET_ENSURELI_FAILED );
3671 
3672  if ( returnvalueCheckLI == RET_LINEARLY_INDEPENDENT )
3673  return SUCCESSFUL_RETURN;
3674 
3675 
3676  /* II) NEW BOUND IS LINEARLY DEPENDENT: */
3677  /* 1) Determine coefficients of linear combination,
3678  * cf. M.J. Best. Applied Mathematics and Parallel Computing, chapter:
3679  * An Algorithm for the Solution of the Parametric Quadratic Programming
3680  * Problem, pages 57-76. Physica-Verlag, Heidelberg, 1996. */
3681  int* FR_idx;
3682  bounds.getFree( )->getNumberArray( &FR_idx );
3683 
3684  int* FX_idx;
3685  bounds.getFixed( )->getNumberArray( &FX_idx );
3686 
3687  int* AC_idx;
3688  constraints.getActive( )->getNumberArray( &AC_idx );
3689 
3690  real_t* xiC = new real_t[nAC];
3691  real_t* xiC_TMP = new real_t[nAC];
3692  real_t* xiB = new real_t[nFX];
3693 
3694  /* 2) Calculate xiC. */
3695  if ( nAC > 0 )
3696  {
3697  if ( B_status == ST_LOWER )
3698  {
3699  for( i=0; i<nAC; ++i )
3700  xiC_TMP[i] = QQ(number,nZ+i);
3701  }
3702  else
3703  {
3704  for( i=0; i<nAC; ++i )
3705  xiC_TMP[i] = -QQ(number,nZ+i);
3706  }
3707 
3708  if ( backsolveT( xiC_TMP, BT_TRUE, xiC ) != SUCCESSFUL_RETURN )
3709  {
3710  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3711  return THROWERROR( RET_ENSURELI_FAILED_TQ );
3712  }
3713  }
3714 
3715  /* 3) Calculate xiB. */
3716  A->transTimes(constraints.getActive(), bounds.getFixed(), 1, -1.0, xiC, nAC, 0.0, xiB, nFX);
3717 
3718 
3719  /* III) DETERMINE CONSTRAINT/BOUND TO BE REMOVED. */
3720  real_t y_min = options.maxDualJump;
3721  int y_min_number = -1;
3722  int y_min_number_bound = -1;
3723  BooleanType y_min_isBound = BT_FALSE;
3724 
3725  real_t* num = new real_t[nV];
3726 
3727  /* 1) Constraints. */
3728  for( i=0; i<nAC; ++i )
3729  {
3730  ii = AC_idx[i];
3731  num[i] = y[nV+ii];
3732  }
3733 
3734  performRatioTest( nAC,AC_idx,&constraints, num,xiC, options.epsDen,options.epsDen, y_min,y_min_number );
3735 
3736 
3737  /* 2) Bounds. */
3738  for( i=0; i<nFX; ++i )
3739  {
3740  ii = FX_idx[i];
3741  num[i] = y[ii];
3742  }
3743 
3744  performRatioTest( nFX,FX_idx,&bounds, num,xiB, options.epsDen,options.epsDen, y_min,y_min_number_bound );
3745 
3746  if ( y_min_number_bound >= 0 )
3747  {
3748  y_min_number = y_min_number_bound;
3749  y_min_isBound = BT_TRUE;
3750  }
3751 
3752  delete[] num;
3753 
3754 
3755  /* IV) REMOVE CONSTRAINT/BOUND FOR RESOLVING LINEAR DEPENDENCE: */
3756  char messageString[80];
3757 
3758  if ( y_min_number >= 0 )
3759  {
3760  /* Update Lagrange multiplier... */
3761  for( i=0; i<nAC; ++i )
3762  {
3763  ii = AC_idx[i];
3764  y[nV+ii] -= y_min * xiC[i];
3765  }
3766  for( i=0; i<nFX; ++i )
3767  {
3768  ii = FX_idx[i];
3769  y[ii] -= y_min * xiB[i];
3770  }
3771 
3772  /* ... also for newly active bound ... */
3773  if ( B_status == ST_LOWER )
3774  y[number] = y_min;
3775  else
3776  y[number] = -y_min;
3777 
3778  /* ... and for bound to be removed. */
3779  if ( y_min_isBound == BT_TRUE )
3780  {
3781  #ifndef __XPCTARGET__
3782  snprintf( messageString,80,"bound no. %d.",y_min_number );
3783  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3784  #endif
3785 
3786  if ( removeBound( y_min_number,BT_TRUE,BT_FALSE,BT_FALSE ) != SUCCESSFUL_RETURN )
3787  {
3788  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3789 
3790  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
3791  }
3792  y[y_min_number] = 0.0;
3793  }
3794  else
3795  {
3796  #ifndef __XPCTARGET__
3797  snprintf( messageString,80,"constraint no. %d.",y_min_number );
3798  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3799  #endif
3800 
3801  if ( removeConstraint( y_min_number,BT_TRUE,BT_FALSE,BT_FALSE ) != SUCCESSFUL_RETURN )
3802  {
3803  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3804 
3805  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
3806  }
3807 
3808  y[nV+y_min_number] = 0.0;
3809  }
3810  }
3811  else
3812  {
3813  /* no constraint/bound can be removed => QP is infeasible! */
3814  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3815 
3816  return setInfeasibilityFlag( RET_ENSURELI_FAILED_NOINDEX );
3817  }
3818 
3819  delete[] xiB; delete[] xiC_TMP; delete[] xiC;
3820 
3821  return getGlobalMessageHandler( )->throwInfo( RET_LI_RESOLVED,0,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
3822 }
3823 
3824 
3825 
3826 /*
3827  * r e m o v e C o n s t r a i n t
3828  */
3829 returnValue QProblem::removeConstraint( int number,
3830  BooleanType updateCholesky,
3831  BooleanType allowFlipping,
3832  BooleanType ensureNZC
3833  )
3834 {
3835  int i, j, ii, jj;
3836  returnValue returnvalue = SUCCESSFUL_RETURN;
3837  BooleanType hasFlipped = BT_FALSE;
3838 
3839  /* consistency check */
3840  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
3841  ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
3842  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
3843  ( getStatus( ) == QPS_SOLVED ) )
3844  {
3845  return THROWERROR( RET_UNKNOWN_BUG );
3846  }
3847 
3848  /* some definitions */
3849  int nV = getNV( );
3850  int nFR = getNFR( );
3851  int nAC = getNAC( );
3852  int nZ = getNZ( );
3853 
3854  int tcol = sizeT - nAC;
3855  int number_idx = constraints.getActive( )->getIndex( number );
3856 
3857  int addIdx;
3858  BooleanType addBoundNotConstraint;
3859  SubjectToStatus addStatus;
3860  BooleanType exchangeHappened = BT_FALSE;
3861 
3862 
3863  /* consistency checks */
3864  if ( constraints.getStatus( number ) == ST_INACTIVE )
3865  return THROWERROR( RET_CONSTRAINT_NOT_ACTIVE );
3866 
3867  if ( ( number_idx < 0 ) || ( number_idx >= nAC ) )
3868  return THROWERROR( RET_CONSTRAINT_NOT_ACTIVE );
3869 
3870 
3871  int* FR_idx;
3872  bounds.getFree( )->getNumberArray( &FR_idx );
3873 
3874  /* N) PERFORM ZERO CURVATURE TEST. */
3875  if (ensureNZC == BT_TRUE)
3876  {
3877  returnvalue = ensureNonzeroCurvature(BT_FALSE, number, exchangeHappened, addBoundNotConstraint, addIdx, addStatus);
3878 
3879  if (returnvalue != SUCCESSFUL_RETURN)
3880  return returnvalue;
3881  }
3882 
3883  /* save index sets and decompositions for flipping bounds strategy */
3884  if ( ( exchangeHappened == BT_FALSE ) && ( options.enableFlippingBounds == BT_TRUE ) && ( allowFlipping == BT_TRUE ) )
3885  flipper.set( &bounds,R,&constraints,Q,T );
3886 
3887  /* I) REMOVE <number>th ROW FROM T,
3888  * i.e. shift rows number+1 through nAC upwards (instead of the actual
3889  * constraint number its corresponding index within matrix A is used). */
3890  if ( number_idx < nAC-1 )
3891  {
3892  for( i=(number_idx+1); i<nAC; ++i )
3893  for( j=(nAC-i-1); j<nAC; ++j )
3894  TT(i-1,tcol+j) = TT(i,tcol+j);
3895  /* gimmick: write zeros into the last row of T */
3896  for( j=0; j<nAC; ++j )
3897  TT(nAC-1,tcol+j) = 0.0;
3898 
3899 
3900  /* II) RESTORE TRIANGULAR FORM OF T,
3901  * use column-wise Givens rotations to restore reverse triangular form
3902  * of T simultanenous change of Q (i.e. Y). */
3903  real_t c, s, nu;
3904 
3905  for( j=(nAC-2-number_idx); j>=0; --j )
3906  {
3907  computeGivens( TT(nAC-2-j,tcol+1+j),TT(nAC-2-j,tcol+j), TT(nAC-2-j,tcol+1+j),TT(nAC-2-j,tcol+j),c,s );
3908  nu = s/(1.0+c);
3909 
3910  for( i=(nAC-j-1); i<(nAC-1); ++i )
3911  applyGivens( c,s,nu,TT(i,tcol+1+j),TT(i,tcol+j), TT(i,tcol+1+j),TT(i,tcol+j) );
3912 
3913  for( i=0; i<nFR; ++i )
3914  {
3915  ii = FR_idx[i];
3916  applyGivens( c,s,nu,QQ(ii,nZ+1+j),QQ(ii,nZ+j), QQ(ii,nZ+1+j),QQ(ii,nZ+j) );
3917  }
3918  }
3919  }
3920  else
3921  {
3922  /* gimmick: write zeros into the last row of T */
3923  for( j=0; j<nAC; ++j )
3924  TT(nAC-1,tcol+j) = 0.0;
3925  }
3926 
3927 
3928  if ( ( updateCholesky == BT_TRUE ) &&
3929  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
3930  {
3931  /* III) UPDATE CHOLESKY DECOMPOSITION,
3932  * calculate new additional column (i.e. [r sqrt(rho2)]')
3933  * of the Cholesky factor R. */
3934  real_t* Hz = new real_t[nFR];
3935  real_t* z = new real_t[nFR];
3936  real_t rho2 = 0.0;
3937 
3938  /* 1) Calculate Hz = H*z, where z is the new rightmost column of Z
3939  * (i.e. the old leftmost column of Y). */
3940  for( j=0; j<nFR; ++j )
3941  z[j] = QQ(FR_idx[j],nZ);
3942  H->times(bounds.getFree(), bounds.getFree(), 1, 1.0, z, nFR, 0.0, Hz, nFR);
3943  delete[] z;
3944 
3945  if ( nZ > 0 )
3946  {
3947  real_t* ZHz = new real_t[nZ];
3948  for ( i=0; i<nZ; ++i )
3949  ZHz[i] = 0.0;
3950  real_t* r = new real_t[nZ];
3951 
3952  /* 2) Calculate ZHz = Z'*Hz (old Z). */
3953  for( j=0; j<nFR; ++j )
3954  {
3955  jj = FR_idx[j];
3956 
3957  for( i=0; i<nZ; ++i )
3958  ZHz[i] += QQ(jj,i) * Hz[j];
3959  }
3960 
3961  /* 3) Calculate r = R^-T * ZHz. */
3962  if ( backsolveR( ZHz,BT_TRUE,r ) != SUCCESSFUL_RETURN )
3963  {
3964  delete[] Hz; delete[] r; delete[] ZHz;
3965  return THROWERROR( RET_REMOVECONSTRAINT_FAILED );
3966  }
3967 
3968  /* 4) Calculate rho2 = rho^2 = z'*Hz - r'*r
3969  * and store r into R. */
3970  for( i=0; i<nZ; ++i )
3971  {
3972  rho2 -= r[i]*r[i];
3973  RR(i,nZ) = r[i];
3974  }
3975 
3976  delete[] r; delete[] ZHz;
3977  }
3978 
3979  /* 5) Store rho into R. */
3980  for( j=0; j<nFR; ++j )
3981  rho2 += QQ(FR_idx[j],nZ) * Hz[j];
3982 
3983  delete[] Hz;
3984 
3985  if ( ( options.enableFlippingBounds == BT_TRUE ) && ( allowFlipping == BT_TRUE ) && ( exchangeHappened == BT_FALSE ) )
3986  {
3987  if ( rho2 > options.epsFlipping )
3988  RR(nZ,nZ) = sqrt( rho2 );
3989  else
3990  {
3991  hessianType = HST_SEMIDEF;
3992 
3993  flipper.get( &bounds,R,&constraints,Q,T );
3994  constraints.flipFixed(number);
3995  idxAddC = number;
3996 
3997  switch (constraints.getStatus(number))
3998  {
3999  case ST_LOWER:
4000  lbA[number] = ubA[number]; Ax_l[number] = -Ax_u[number]; break;
4001  case ST_UPPER:
4002  ubA[number] = lbA[number]; Ax_u[number] = -Ax_l[number]; break;
4003  default:
4004  return THROWERROR( RET_MOVING_BOUND_FAILED );
4005  }
4006 
4007  hasFlipped = BT_TRUE;
4008  }
4009  }
4010  else if ( exchangeHappened == BT_FALSE )
4011  {
4012  if ( rho2 > ZERO )
4013  RR(nZ,nZ) = sqrt( rho2 );
4014  else
4015  {
4016  if ( allowFlipping == BT_FALSE )
4017  {
4018  RR(nZ,nZ) = 100.0*EPS;
4019  }
4020  else
4021  {
4022  hessianType = HST_SEMIDEF;
4023  return THROWERROR( RET_HESSIAN_NOT_SPD );
4024  }
4025  }
4026  }
4027  }
4028 
4029 
4030  /* IV) UPDATE INDICES */
4031  if ( hasFlipped == BT_FALSE )
4032  {
4033  idxRemC = number;
4034  if ( constraints.moveActiveToInactive( number ) != SUCCESSFUL_RETURN )
4035  return THROWERROR( RET_REMOVECONSTRAINT_FAILED );
4036  }
4037 
4038  if (exchangeHappened == BT_TRUE)
4039  {
4040  /* add bound or constraint */
4041 
4042  /* hessianType = HST_SEMIDEF; */
4043  RR(nZ,nZ) = 0.0;
4044 
4045  if ( addBoundNotConstraint )
4046  addBound(addIdx, addStatus, BT_TRUE, BT_FALSE);
4047  else
4048  addConstraint(addIdx, addStatus, BT_TRUE, BT_FALSE);
4049  }
4050 
4051  return SUCCESSFUL_RETURN;
4052 }
4053 
4054 
4055 /*
4056  * r e m o v e B o u n d
4057  */
4058 returnValue QProblem::removeBound( int number,
4059  BooleanType updateCholesky,
4060  BooleanType allowFlipping,
4061  BooleanType ensureNZC
4062  )
4063 {
4064  int i, j, ii, jj;
4065  returnValue returnvalue = SUCCESSFUL_RETURN;
4066  int addIdx;
4067  BooleanType addBoundNotConstraint;
4068  SubjectToStatus addStatus;
4069  BooleanType exchangeHappened = BT_FALSE;
4070 
4071 
4072  /* consistency checks */
4073  if ( bounds.getStatus( number ) == ST_INACTIVE )
4074  return THROWERROR( RET_BOUND_NOT_ACTIVE );
4075 
4076  if ( ( getStatus( ) == QPS_NOTINITIALISED ) ||
4077  ( getStatus( ) == QPS_AUXILIARYQPSOLVED ) ||
4078  ( getStatus( ) == QPS_HOMOTOPYQPSOLVED ) ||
4079  ( getStatus( ) == QPS_SOLVED ) )
4080  {
4081  return THROWERROR( RET_UNKNOWN_BUG );
4082  }
4083 
4084  /* some definitions */
4085  int nV = getNV( );
4086  int nFR = getNFR( );
4087  int nAC = getNAC( );
4088  int nZ = getNZ( );
4089 
4090  int tcol = sizeT - nAC;
4091 
4092  /* N) PERFORM ZERO CURVATURE TEST. */
4093  if (ensureNZC == BT_TRUE)
4094  {
4095  returnvalue = ensureNonzeroCurvature(BT_TRUE, number, exchangeHappened, addBoundNotConstraint, addIdx, addStatus);
4096 
4097  if (returnvalue != SUCCESSFUL_RETURN)
4098  return returnvalue;
4099  }
4100 
4101  /* save index sets and decompositions for flipping bounds strategy */
4102  if ( ( options.enableFlippingBounds == BT_TRUE ) && ( allowFlipping == BT_TRUE ) && ( exchangeHappened == BT_FALSE ) )
4103  flipper.set( &bounds,R,&constraints,Q,T );
4104 
4105  /* I) UPDATE INDICES */
4106  idxRemB = number;
4107  if ( bounds.moveFixedToFree( number ) != SUCCESSFUL_RETURN )
4108  return THROWERROR( RET_REMOVEBOUND_FAILED );
4109 
4110  int* FR_idx;
4111  bounds.getFree( )->getNumberArray( &FR_idx );
4112 
4113  /* I) APPEND <nFR+1>th UNITY VECTOR TO Q. */
4114  int nnFRp1 = FR_idx[nFR];
4115  for( i=0; i<nFR; ++i )
4116  {
4117  ii = FR_idx[i];
4118  QQ(ii,nFR) = 0.0;
4119  QQ(nnFRp1,i) = 0.0;
4120  }
4121  QQ(nnFRp1,nFR) = 1.0;
4122 
4123  if ( nAC > 0 )
4124  {
4125  /* store new column a in a temporary vector instead of shifting T one column to the left and appending a */
4126  int* AC_idx;
4127  constraints.getActive( )->getNumberArray( &AC_idx );
4128 
4129  real_t* tmp = new real_t[nAC];
4130  A->getCol(number, constraints.getActive(), 1.0, tmp);
4131 
4132 
4133  /* II) RESTORE TRIANGULAR FORM OF T,
4134  * use column-wise Givens rotations to restore reverse triangular form
4135  * of T = [T A(:,number)], simultanenous change of Q (i.e. Y and Z). */
4136  real_t c, s, nu;
4137 
4138  for( j=(nAC-1); j>=0; --j )
4139  {
4140  computeGivens( tmp[nAC-1-j],TT(nAC-1-j,tcol+j),TT(nAC-1-j,tcol+j),tmp[nAC-1-j],c,s );
4141  nu = s/(1.0+c);
4142 
4143  for( i=(nAC-j); i<nAC; ++i )
4144  applyGivens( c,s,nu,tmp[i],TT(i,tcol+j),TT(i,tcol+j),tmp[i] );
4145 
4146  for( i=0; i<=nFR; ++i )
4147  {
4148  ii = FR_idx[i];
4149  /* nZ+1+nAC = nFR+1 / nZ+(1) = nZ+1 */
4150  applyGivens( c,s,nu,QQ(ii,nZ+1+j),QQ(ii,nZ+j),QQ(ii,nZ+1+j),QQ(ii,nZ+j) );
4151  }
4152  }
4153 
4154  delete[] tmp;
4155  }
4156 
4157 
4158  if ( ( updateCholesky == BT_TRUE ) &&
4159  ( hessianType != HST_ZERO ) && ( hessianType != HST_IDENTITY ) )
4160  {
4161  /* III) UPDATE CHOLESKY DECOMPOSITION,
4162  * calculate new additional column (i.e. [r sqrt(rho2)]')
4163  * of the Cholesky factor R: */
4164  real_t z2 = QQ(nnFRp1,nZ);
4165  real_t rho2 = H->diag(nnFRp1)*z2*z2; /* rho2 = h2*z2*z2 */
4166 
4167  if ( nFR > 0 )
4168  {
4169  /* Attention: Index list of free variables has already grown by one! */
4170  real_t* Hz = new real_t[nFR+1];
4171  real_t* z = new real_t[nFR+1];
4172  /* 1) Calculate R'*r = Zfr'*Hfr*z1 + z2*Zfr'*h1 =: Zfr'*Hz + z2*Zfr'*h1 =: rhs and
4173  * rho2 = z1'*Hfr*z1 + 2*z2*h1'*z1 + h2*z2^2 - r'*r =: z1'*Hz + 2*z2*h1'*z1 + h2*z2^2 - r'r */
4174  for( j=0; j<nFR; ++j )
4175  z[j] = QQ(FR_idx[j],nZ);
4176  z[nFR] = 0.0;
4177  H->times(bounds.getFree(), bounds.getFree(), 1, 1.0, z, nFR+1, 0.0, Hz, nFR+1);
4178 
4179  H->getCol(nnFRp1, bounds.getFree(), 1.0, z);
4180 
4181  if ( nZ > 0 )
4182  {
4183  real_t* r = new real_t[nZ];
4184  real_t* rhs = new real_t[nZ];
4185  for( i=0; i<nZ; ++i )
4186  rhs[i] = 0.0;
4187 
4188  /* 2) Calculate rhs. */
4189  for( j=0; j<nFR; ++j )
4190  {
4191  jj = FR_idx[j];
4192  for( i=0; i<nZ; ++i )
4193  /* Zfr' * ( Hz + z2*h1 ) */
4194  rhs[i] += QQ(jj,i) * ( Hz[j] + z2 * z[j] );
4195  }
4196 
4197  /* 3) Calculate r = R^-T * rhs. */
4198  if ( backsolveR( rhs,BT_TRUE,BT_TRUE,r ) != SUCCESSFUL_RETURN )
4199  {
4200  delete[] z;
4201  delete[] Hz; delete[] r; delete[] rhs;
4202  return THROWERROR( RET_REMOVEBOUND_FAILED );
4203  }
4204 
4205 
4206  /* 4) Calculate rho2 = rho^2 = z'*Hz - r'*r
4207  * and store r into R. */
4208  for( i=0; i<nZ; ++i )
4209  {
4210  rho2 -= r[i]*r[i];
4211  RR(i,nZ) = r[i];
4212  }
4213 
4214  delete[] rhs; delete[] r;
4215  }
4216 
4217  for( j=0; j<nFR; ++j )
4218  {
4219  jj = FR_idx[j];
4220  /* z1' * ( Hz + 2*z2*h1 ) */
4221  rho2 += QQ(jj,nZ) * ( Hz[j] + 2.0*z2*z[j] );
4222  }
4223 
4224  delete[] z;
4225  delete[] Hz;
4226  }
4227 
4228  /* 5) Store rho into R. */
4229  if ( ( options.enableFlippingBounds == BT_TRUE ) && ( allowFlipping == BT_TRUE ) && ( exchangeHappened == BT_FALSE ) )
4230  {
4231  if ( rho2 > options.epsFlipping )
4232  RR(nZ,nZ) = sqrt( rho2 );
4233  else
4234  {
4235  hessianType = HST_SEMIDEF;
4236 
4237  flipper.get( &bounds,R,&constraints,Q,T );
4238  bounds.flipFixed(number);
4239  idxAddB = number;
4240 
4241  switch (bounds.getStatus(number))
4242  {
4243  case ST_LOWER: lb[number] = ub[number]; break;
4244  case ST_UPPER: ub[number] = lb[number]; break;
4245  default: return THROWERROR( RET_MOVING_BOUND_FAILED );
4246  }
4247 
4248  }
4249  }
4250  else if ( exchangeHappened == BT_FALSE )
4251  {
4252  if ( rho2 > ZERO )
4253  RR(nZ,nZ) = sqrt( rho2 );
4254  else
4255  {
4256  if ( allowFlipping == BT_FALSE )
4257  RR(nZ,nZ) = 100.0*EPS;
4258  else
4259  {
4260  hessianType = HST_SEMIDEF;
4261  return THROWERROR( RET_HESSIAN_NOT_SPD );
4262  }
4263  }
4264  }
4265  else
4266  {
4267  /* add bound or constraint */
4268 
4269  /* hessianType = HST_SEMIDEF; */
4270  RR(nZ,nZ) = 0.0;
4271 
4272  if ( addBoundNotConstraint )
4273  addBound(addIdx, addStatus, BT_TRUE, BT_FALSE);
4274  else
4275  addConstraint(addIdx, addStatus, BT_TRUE, BT_FALSE);
4276  }
4277  }
4278 
4279  return SUCCESSFUL_RETURN;
4280 }
4281 
4282 
4283 
4284 returnValue QProblem::performPlainRatioTest( int nIdx,
4285  const int* const idxList,
4286  const real_t* const num,
4287  const real_t* const den,
4288  real_t epsNum,
4289  real_t epsDen,
4290  real_t& t,
4291  int& BC_idx
4292  ) const
4293 {
4294  int i;
4295  for (i = 0; i < nIdx; i++)
4296  if ( (num[i] > epsNum) && (den[i] > epsDen) && (t * den[i] > num[i]) )
4297  {
4298  t = num[i] / den[i];
4299  BC_idx = idxList[i];
4300  }
4301 
4302  return SUCCESSFUL_RETURN;
4303 }
4304 
4305 
4306 returnValue QProblem::ensureNonzeroCurvature(
4307  BooleanType removeBoundNotConstraint,
4308  int remIdx,
4309  BooleanType &exchangeHappened,
4310  BooleanType &addBoundNotConstraint,
4311  int &addIdx,
4312  SubjectToStatus &addStatus
4313  )
4314 {
4315  int i, ii;
4316  int addLBndIdx = -1, addLCnstrIdx = -1, addUBndIdx = -1, addUCnstrIdx = -1; /* exchange indices */
4317  int *FX_idx, *AC_idx, *IAC_idx;
4318  returnValue returnvalue = SUCCESSFUL_RETURN;
4319 
4320  int nV = getNV( );
4321  int nFR = getNFR( );
4322  int nAC = getNAC( );
4323  int nC = getNC( );
4324  int nFX = getNFX( );
4325  int nIAC = getNIAC( );
4326 
4327  int* FR_idx;
4328  bounds.getFree( )->getNumberArray( &FR_idx );
4329 
4330  real_t *delta_g = new real_t[nV];
4331  real_t *delta_xFX = new real_t[nFX];
4332  real_t *delta_xFR = new real_t[nFR];
4333  real_t *delta_yAC = new real_t[nAC];
4334  real_t *delta_yFX = new real_t[nFX];
4335 
4336  bounds.getFixed( )->getNumberArray( &FX_idx );
4337  constraints.getActive( )->getNumberArray( &AC_idx );
4338  constraints.getInactive( )->getNumberArray( &IAC_idx );
4339 
4340  addBoundNotConstraint = BT_TRUE;
4341  addStatus = ST_INACTIVE;
4342  exchangeHappened = BT_FALSE;
4343 
4344  if (removeBoundNotConstraint)
4345  {
4346  int dim = nV < nC ? nC : nV;
4347  real_t *nul = new real_t[dim];
4348  real_t *ek = new real_t[nV]; /* minus e_k (bound k is removed) */
4349  for (ii = 0; ii < dim; ++ii)
4350  nul[ii]=0.0;
4351  for (ii = 0; ii < nV; ++ii)
4352  ek[ii]=0.0;
4353  ek[remIdx] = bounds.getStatus(remIdx) == ST_LOWER ? 1.0 : -1.0;
4354 
4355  returnvalue = determineStepDirection (nul, nul, nul, ek, ek,
4356  BT_FALSE, BT_FALSE,
4357  delta_xFX, delta_xFR, delta_yAC, delta_yFX);
4358  delete[] ek;
4359  delete[] nul;
4360  }
4361  else
4362  {
4363  real_t *nul = new real_t[nV];
4364  real_t *ek = new real_t[nC]; /* minus e_k (constraint k is removed) */
4365  for (ii = 0; ii < nV; ++ii)
4366  nul[ii]=0.0;
4367  for (ii = 0; ii < nC; ++ii)
4368  ek[ii]=0.0;
4369  ek[remIdx] = constraints.getStatus(remIdx) == ST_LOWER ? 1.0 : -1.0;
4370 
4371  returnvalue = determineStepDirection (nul,
4372  ek, ek, nul, nul,
4373  BT_FALSE, BT_TRUE,
4374  delta_xFX, delta_xFR, delta_yAC, delta_yFX);
4375  delete[] ek;
4376  delete[] nul;
4377  }
4378 
4379  /* compute the weight in inf-norm */
4380  real_t normXi = 0.0;
4381  for (ii = 0; ii < nAC; ++ii)
4382  {
4383  real_t a = fabs (delta_yAC[ii]);
4384  if (normXi < a) normXi = a;
4385  }
4386  for (ii = 0; ii < nFX; ++ii)
4387  {
4388  real_t a = fabs (delta_yFX[ii]);
4389  if (normXi < a) normXi = a;
4390  }
4391 
4392  /* look at the "zero" in a relative inf-norm */
4393  real_t normS = 0.0;
4394  for (ii = 0; ii < nFX; ++ii)
4395  {
4396  real_t a = fabs (delta_xFX[ii]);
4397  if (normS < a) normS = a;
4398  }
4399  for (ii = 0; ii < nFR; ++ii)
4400  {
4401  real_t a = fabs (delta_xFR[ii]);
4402  if (normS < a) normS = a;
4403  }
4404 
4405  /* relative test against zero in inf-norm */
4406  if (normXi < options.epsNZCTests * normS)
4407  {
4408  /* determine jump in x via ratio tests */
4409  real_t sigmaLBnd, sigmaLCnstr, sigmaUBnd, sigmaUCnstr, sigma;
4410 
4411  /* bounds */
4412 
4413  /* compress x-u */
4414  real_t *x_W = new real_t[getMax(1,nFR)];
4415  for (i = 0; i < nFR; i++)
4416  {
4417  ii = FR_idx[i];
4418  x_W[i] = ub[ii] - x[ii];
4419  }
4420  /* performRatioTest( nFR,FR_idx,&bounds, x_W,delta_xFR, options.epsDen,options.epsDen, sigmaUBnd,addUBndIdx ); */
4421  sigmaUBnd = options.maxPrimalJump;
4422  addUBndIdx = -1;
4423  performPlainRatioTest(nFR, FR_idx, x_W, delta_xFR, options.epsNum, options.epsDen, sigmaUBnd, addUBndIdx);
4424  if (removeBoundNotConstraint == BT_TRUE && bounds.getStatus(remIdx) == ST_LOWER)
4425  {
4426  /* also consider bound which is to be removed */
4427  real_t eins = 1.0;
4428  x_W[0] = ub[remIdx] - x[remIdx];
4429  performPlainRatioTest(1, &remIdx, x_W, &eins, options.epsNum, options.epsDen, sigmaUBnd, addUBndIdx);
4430  }
4431 
4432  /* compress x-l */
4433  for (i = 0; i < nFR; i++)
4434  {
4435  ii = FR_idx[i];
4436  x_W[i] = x[ii] - lb[ii];
4437  }
4438  for (i = 0; i < nFR; i++)
4439  delta_xFR[i] = -delta_xFR[i];
4440  /* performRatioTest( nFR,FR_idx,&bounds, x_W,delta_xFR, options.epsDen,options.epsDen, sigmaLBnd,addLBndIdx ); */
4441  sigmaLBnd = options.maxPrimalJump;
4442  addLBndIdx = -1;
4443  performPlainRatioTest(nFR, FR_idx, x_W, delta_xFR, options.epsNum, options.epsDen, sigmaLBnd, addLBndIdx);
4444  if (removeBoundNotConstraint == BT_TRUE && bounds.getStatus(remIdx) == ST_UPPER)
4445  {
4446  /* also consider bound which is to be removed */
4447  real_t eins = 1.0;
4448  x_W[0] = x[remIdx] - lb[remIdx];
4449  performPlainRatioTest(1, &remIdx, x_W, &eins, options.epsNum, options.epsDen, sigmaLBnd, addLBndIdx);
4450  }
4451  for (i = 0; i < nFR; i++)
4452  delta_xFR[i] = -delta_xFR[i];
4453 
4454  delete[] x_W;
4455 
4456  /* constraints */
4457 
4458  /* compute As (compressed to inactive constraints) */
4459  real_t *As = new real_t[nIAC];
4460  A->times(constraints.getInactive(), bounds.getFixed(), 1, 1.0, delta_xFX, nFX, 0.0, As, nIAC);
4461  A->times(constraints.getInactive(), bounds.getFree(), 1, 1.0, delta_xFR, nFR, 1.0, As, nIAC);
4462 
4463  /* compress Ax_u */
4464  real_t *Ax_W = new real_t[nIAC];
4465  for (i = 0; i < nIAC; i++)
4466  {
4467  ii = IAC_idx[i];
4468  Ax_W[i] = Ax_u[ii];
4469  }
4470  /* performRatioTest( nIAC,IAC_idx,&constraints, Ax_W,As, options.epsDen,options.epsDen, sigmaUCnstr,addUCnstrIdx ); */
4471  sigmaUCnstr = options.maxPrimalJump;
4472  addUCnstrIdx = -1;
4473  performPlainRatioTest(nIAC, IAC_idx, Ax_W, As, options.epsNum, options.epsDen, sigmaUCnstr, addUCnstrIdx);
4474  if (removeBoundNotConstraint == BT_FALSE && constraints.getStatus(remIdx) == ST_LOWER)
4475  {
4476  /* also consider constraint which is to be removed */
4477  real_t eins = 1.0;
4478  performPlainRatioTest(1, &remIdx, &Ax_u[remIdx], &eins, options.epsNum, options.epsDen, sigmaUCnstr, addUCnstrIdx);
4479  }
4480 
4481  /* compress Ax_l */
4482  for (i = 0; i < nIAC; i++)
4483  {
4484  ii = IAC_idx[i];
4485  Ax_W[i] = Ax_l[ii];
4486  }
4487  for (i = 0; i < nIAC; i++)
4488  As[i] = -As[i];
4489  /* performRatioTest( nIAC,IAC_idx,&constraints, Ax_W,As, options.epsDen,options.epsDen, sigmaLCnstr,addLCnstrIdx ); */
4490  sigmaLCnstr = options.maxPrimalJump;
4491  addLCnstrIdx = -1;
4492  performPlainRatioTest(nIAC, IAC_idx, Ax_W, As, options.epsNum, options.epsDen, sigmaLCnstr, addLCnstrIdx);
4493  if (removeBoundNotConstraint == BT_FALSE && constraints.getStatus(remIdx) == ST_UPPER)
4494  {
4495  /* also consider constraint which is to be removed */
4496  real_t eins = 1.0;
4497  performPlainRatioTest(1, &remIdx, &Ax_l[remIdx], &eins, options.epsNum, options.epsDen, sigmaLCnstr, addLCnstrIdx);
4498  }
4499 
4500  /* perform primal jump */
4501  sigma = options.maxPrimalJump;
4502  if (sigmaUCnstr < sigma) { sigma = sigmaUCnstr; addStatus = ST_UPPER; addBoundNotConstraint = BT_FALSE; addIdx = addUCnstrIdx; }
4503  if (sigmaLCnstr < sigma) { sigma = sigmaLCnstr; addStatus = ST_LOWER; addBoundNotConstraint = BT_FALSE; addIdx = addLCnstrIdx; }
4504  if (sigmaUBnd < sigma) { sigma = sigmaUBnd; addStatus = ST_UPPER; addBoundNotConstraint = BT_TRUE; addIdx = addUBndIdx; }
4505  if (sigmaLBnd < sigma) { sigma = sigmaLBnd; addStatus = ST_LOWER; addBoundNotConstraint = BT_TRUE; addIdx = addLBndIdx; }
4506 
4507  if (sigma >= options.maxPrimalJump)
4508  {
4509  unbounded = BT_TRUE;
4510  returnvalue = RET_HOTSTART_STOPPED_UNBOUNDEDNESS;
4511  }
4512  else
4513  {
4514  for (i = 0; i < nFR; i++)
4515  x[FR_idx[i]] += sigma * delta_xFR[i];
4516 
4517  for (i = 0; i < nFX; i++)
4518  x[FX_idx[i]] += sigma * delta_xFX[i];
4519 
4520  /* update Ax, Ax_u, and Ax_l */
4521  A->times(1, 1.0, x, nV, 0.0, Ax, nC);
4522  for (i = 0; i < nC; i++) Ax_u[i] = ubA[i] - Ax[i];
4523  for (i = 0; i < nC; i++) Ax_l[i] = Ax[i] - lbA[i];
4524 
4525  /* change working set later */
4526  exchangeHappened = BT_TRUE;
4527  }
4528 
4529  delete[] Ax_W;
4530  delete[] As;
4531  }
4532 
4533  delete[] delta_yFX;
4534  delete[] delta_yAC;
4535  delete[] delta_xFR;
4536  delete[] delta_xFX;
4537  delete[] delta_g;
4538 
4539  return returnvalue;
4540 }
4541 
4542 
4543 
4544 /*
4545  * b a c k s o l v e T
4546  */
4547 returnValue QProblem::backsolveT( const real_t* const b, BooleanType transposed, real_t* const a ) const
4548 {
4549  int i, j;
4550  int nT = getNAC( );
4551  int tcol = sizeT - nT;
4552 
4553  real_t sum;
4554 
4555  /* nothing to do */
4556  if ( nT <= 0 )
4557  return SUCCESSFUL_RETURN;
4558 
4559 
4560  /* Solve Ta = b, where T might be transposed. */
4561  if ( transposed == BT_FALSE )
4562  {
4563  /* solve Ta = b */
4564  for( i=0; i<nT; ++i )
4565  {
4566  sum = b[i];
4567  for( j=0; j<i; ++j )
4568  sum -= TT(i,sizeT-1-j) * a[nT-1-j];
4569 
4570  if ( fabs( TT(i,sizeT-1-i) ) > EPS )
4571  a[nT-1-i] = sum / TT(i,sizeT-1-i);
4572  else
4573  return THROWERROR( RET_DIV_BY_ZERO );
4574  }
4575  }
4576  else
4577  {
4578  /* solve T^T*a = b */
4579  for( i=0; i<nT; ++i )
4580  {
4581  sum = b[i];
4582  for( j=0; j<i; ++j )
4583  sum -= TT(nT-1-j,tcol+i) * a[nT-1-j];
4584 
4585  if ( fabs( TT(nT-1-i,tcol+i) ) > EPS )
4586  a[nT-1-i] = sum / TT(nT-1-i,tcol+i);
4587  else
4588  return THROWERROR( RET_DIV_BY_ZERO );
4589  }
4590  }
4591 
4592 
4593  return SUCCESSFUL_RETURN;
4594 }
4595 
4596 
4597 /*
4598  * d e t e r m i n e D a t a S h i f t
4599  */
4600 returnValue QProblem::determineDataShift( const real_t* const g_new, const real_t* const lbA_new, const real_t* const ubA_new,
4601  const real_t* const lb_new, const real_t* const ub_new,
4602  real_t* const delta_g, real_t* const delta_lbA, real_t* const delta_ubA,
4603  real_t* const delta_lb, real_t* const delta_ub,
4604  BooleanType& Delta_bC_isZero, BooleanType& Delta_bB_isZero
4605  )
4606 {
4607  int i, ii;
4608  int nC = getNC( );
4609  int nAC = getNAC( );
4610 
4611  int* FX_idx;
4612  int* AC_idx;
4613 
4614  bounds.getFixed( )->getNumberArray( &FX_idx );
4615  constraints.getActive( )->getNumberArray( &AC_idx );
4616 
4617 
4618 
4619  /* I) DETERMINE DATA SHIFT FOR BOUNDS */
4620  QProblemB::determineDataShift( g_new,lb_new,ub_new,
4621  delta_g,delta_lb,delta_ub,
4622  Delta_bB_isZero );
4623 
4624 
4625  /* II) DETERMINE DATA SHIFT FOR CONSTRAINTS */
4626  /* 1) Calculate shift directions. */
4627  for( i=0; i<nC; ++i )
4628  {
4629  /* if lower constraints' bounds are to be disabled or do not exist, shift them to -infinity */
4630  if ( lbA_new != 0 )
4631  delta_lbA[i] = lbA_new[i] - lbA[i];
4632  else
4633  delta_lbA[i] = -INFTY - lbA[i];
4634  }
4635 
4636  for( i=0; i<nC; ++i )
4637  {
4638  /* if upper constraints' bounds are to be disabled or do not exist, shift them to infinity */
4639  if ( ubA_new != 0 )
4640  delta_ubA[i] = ubA_new[i] - ubA[i];
4641  else
4642  delta_ubA[i] = INFTY - ubA[i];
4643  }
4644 
4645  /* 2) Determine if active constraints' bounds are to be shifted. */
4646  Delta_bC_isZero = BT_TRUE;
4647 
4648  for ( i=0; i<nAC; ++i )
4649  {
4650  ii = AC_idx[i];
4651 
4652  if ( ( fabs( delta_lbA[ii] ) > EPS ) || ( fabs( delta_ubA[ii] ) > EPS ) )
4653  {
4654  Delta_bC_isZero = BT_FALSE;
4655  break;
4656  }
4657  }
4658 
4659  return SUCCESSFUL_RETURN;
4660 }
4661 
4662 
4663 /*
4664  * d e t e r m i n e S t e p D i r e c t i o n
4665  */
4666 returnValue QProblem::determineStepDirection( const real_t* const delta_g, const real_t* const delta_lbA, const real_t* const delta_ubA,
4667  const real_t* const delta_lb, const real_t* const delta_ub,
4668  BooleanType Delta_bC_isZero, BooleanType Delta_bB_isZero,
4669  real_t* const delta_xFX, real_t* const delta_xFR,
4670  real_t* const delta_yAC, real_t* const delta_yFX
4671  )
4672 {
4673  int i, j, ii, jj, r;
4674  int nV = getNV( );
4675  int nFR = getNFR( );
4676  int nFX = getNFX( );
4677  int nAC = getNAC( );
4678  int nZ = getNZ( );
4679 
4680  int* FR_idx;
4681  int* FX_idx;
4682  int* AC_idx;
4683 
4684  bounds.getFree( )->getNumberArray( &FR_idx );
4685  bounds.getFixed( )->getNumberArray( &FX_idx );
4686  constraints.getActive( )->getNumberArray( &AC_idx );
4687 
4688 
4689  /* I) DETERMINE delta_xFX (this is exact, does not need refinement) */
4690  if ( Delta_bB_isZero == BT_FALSE )
4691  {
4692  for( i=0; i<nFX; ++i )
4693  {
4694  ii = FX_idx[i];
4695 
4696  if ( bounds.getStatus( ii ) == ST_LOWER )
4697  delta_xFX[i] = delta_lb[ii];
4698  else
4699  delta_xFX[i] = delta_ub[ii];
4700  }
4701  }
4702  else
4703  {
4704  for( i=0; i<nFX; ++i )
4705  delta_xFX[i] = 0.0;
4706  }
4707 
4708 
4709  /* tempA and tempB hold the residuals in gFR and bA (= lbA or ubA)
4710  * delta_xFR, delta_yAC hold the steps that get refined */
4711  for ( i=0; i<nFR; ++i )
4712  {
4713  ii = FR_idx[i];
4714  tempA[i] = delta_g[ii];
4715  delta_xFR[i] = 0.0;
4716  }
4717  for ( i=0; i<nAC; ++i )
4718  delta_yAC[i] = 0.0;
4719  if ( Delta_bC_isZero == BT_FALSE )
4720  {
4721  for ( i=0; i<nAC; ++i )
4722  {
4723  ii = AC_idx[i];
4724  if ( constraints.getStatus( ii ) == ST_LOWER )
4725  tempB[i] = delta_lbA[ii];
4726  else
4727  tempB[i] = delta_ubA[ii];
4728  }
4729  }
4730  else
4731  {
4732  for ( i=0; i<nAC; ++i )
4733  tempB[i] = 0.0;
4734  }
4735 
4736  /* Iterative refinement loop for delta_xFRz, delta_xFRy, delta_yAC */
4737  for ( r=0; r<=options.numRefinementSteps; ++r )
4738  {
4739  /* II) DETERMINE delta_xFR */
4740  if ( nFR > 0 )
4741  {
4742  for( i=0; i<nFR; ++i )
4743  delta_xFR_TMP[i] = 0.0;
4744 
4745  /* 1) Determine delta_xFRy. */
4746  if ( nAC > 0 )
4747  {
4748  if ( ( Delta_bC_isZero == BT_TRUE ) && ( Delta_bB_isZero == BT_TRUE ) )
4749  {
4750  for( i=0; i<nAC; ++i )
4751  delta_xFRy[i] = 0.0;
4752  }
4753  else
4754  {
4755  /* compute bA - A * delta_xFX. tempB already holds bA->
4756  * in refinements r>=1, delta_xFX is exactly zero */
4757  if ( ( Delta_bB_isZero == BT_FALSE ) && ( r == 0 ) )
4758  A->times(constraints.getActive(), bounds.getFixed(), 1, -1.0, delta_xFX, nFX, 1.0, tempB, nAC);
4759 
4760  if ( backsolveT( tempB, BT_FALSE, delta_xFRy ) != SUCCESSFUL_RETURN )
4761  return THROWERROR( RET_STEPDIRECTION_FAILED_TQ );
4762 
4763  for( i=0; i<nFR; ++i )
4764  {
4765  ii = FR_idx[i];
4766  for( j=0; j<nAC; ++j )
4767  delta_xFR_TMP[i] += QQ(ii,nZ+j) * delta_xFRy[j];
4768  }
4769  }
4770  }
4771 
4772 
4773  /* 2) Determine delta_xFRz. */
4774  for( i=0; i<nZ; ++i )
4775  delta_xFRz[i] = 0.0;
4776 
4777  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_IDENTITY ) )
4778  {
4779  /* compute Z*delta_gFR [/eps] (delta_gFR is stored in tempA) */
4780  for( j=0; j<nFR; ++j )
4781  {
4782  jj = FR_idx[j];
4783  for( i=0; i<nZ; ++i )
4784  delta_xFRz[i] -= QQ(jj,i) * tempA[j];
4785  }
4786 
4787  if ( hessianType == HST_ZERO )
4788  {
4789  for( i=0; i<nZ; ++i )
4790  delta_xFRz[i] /= options.epsRegularisation;
4791  }
4792  }
4793  else
4794  {
4795  /* compute HMX*delta_xFX. DESTROY delta_gFR that was in tempA */
4796  if ( ( Delta_bB_isZero == BT_FALSE ) && ( r == 0 ) )
4797  H->times(bounds.getFree(), bounds.getFixed(), 1, 1.0, delta_xFX, nFX, 1.0, tempA, nFR);
4798 
4799  /* compute HFR*delta_xFRy */
4800  if ( ( nAC > 0 ) && ( ( Delta_bC_isZero == BT_FALSE ) || ( Delta_bB_isZero == BT_FALSE ) ) )
4801  H->times(bounds.getFree(), bounds.getFree(), 1, 1.0, delta_xFR_TMP, nFR, 1.0, tempA, nFR);
4802 
4803  /* compute ZFR_delta_xFRz = (Z'*HFR*Z) \ Z * (HFR*delta_xFR + HMX*delta_xFX + delta_gFR) */
4804  if ( nZ > 0 )
4805  {
4806  for( j=0; j<nFR; ++j )
4807  {
4808  jj = FR_idx[j];
4809  for( i=0; i<nZ; ++i )
4810  delta_xFRz[i] -= QQ(jj,i) * tempA[j];
4811  }
4812 
4813  if ( backsolveR( delta_xFRz,BT_TRUE,delta_xFRz ) != SUCCESSFUL_RETURN )
4814  return THROWERROR( RET_STEPDIRECTION_FAILED_CHOLESKY );
4815 
4816  if ( backsolveR( delta_xFRz,BT_FALSE,delta_xFRz ) != SUCCESSFUL_RETURN )
4817  return THROWERROR( RET_STEPDIRECTION_FAILED_CHOLESKY );
4818  }
4819  }
4820 
4821  /* compute Z * ZFR_delta_xFRz */
4822  if ( nZ > 0 )
4823  {
4824  for( i=0; i<nFR; ++i )
4825  {
4826  ZFR_delta_xFRz[i] = 0.0;
4827 
4828  ii = FR_idx[i];
4829  for( j=0; j<nZ; ++j )
4830  ZFR_delta_xFRz[i] += QQ(ii,j) * delta_xFRz[j];
4831 
4832  delta_xFR_TMP[i] += ZFR_delta_xFRz[i];
4833  }
4834  }
4835  }
4836 
4837  /* III) DETERMINE delta_yAC */
4838  if ( nAC > 0 ) /* => ( nFR = nZ + nAC > 0 ) */
4839  {
4840  if ( ( hessianType == HST_ZERO ) || ( hessianType == HST_IDENTITY ) )
4841  {
4842  /* if zero: delta_yAC = (T')^-1 * ( Yfr*delta_gFR + eps*delta_xFRy ),
4843  * if identity: delta_yAC = (T')^-1 * ( Yfr*delta_gFR + delta_xFRy )
4844  *
4845  * DESTROY residual_bA that was stored in tempB
4846  * If we come here, residual_gFR in tempA is STILL VALID
4847  */
4848  if ( hessianType == HST_IDENTITY )
4849  {
4850  for( j=0; j<nAC; ++j )
4851  tempB[j] = delta_xFRy[j];
4852  }
4853  else /* hessianType == HST_ZERO */
4854  {
4855  if ( usingRegularisation( ) == BT_TRUE )
4856  {
4857  for( j=0; j<nAC; ++j )
4858  tempB[j] = options.epsRegularisation*delta_xFRy[j];
4859  }
4860  else
4861  {
4862  for( j=0; j<nAC; ++j )
4863  tempB[j] = 0.0;
4864  }
4865  }
4866 
4867  for( j=0; j<nAC; ++j )
4868  {
4869  for( i=0; i<nFR; ++i )
4870  {
4871  ii = FR_idx[i];
4872  tempB[j] += QQ(ii,nZ+j) * tempA[i];
4873  }
4874  }
4875  }
4876  else
4877  {
4878  /* Compute HFR * delta_xFR + HMX*delta_xFX
4879  * Here, tempA holds (HFR*delta_xFRy + HMX*delta_xFX) */
4880  if ( nZ > 0 )
4881  H->times(bounds.getFree(), bounds.getFree(), 1, 1.0, ZFR_delta_xFRz, nFR, 1.0, tempA, nFR);
4882 
4883  for( i=0; i<nAC; ++i)
4884  {
4885  tempB[i] = 0.0;
4886  for( j=0; j<nFR; ++j )
4887  {
4888  jj = FR_idx[j];
4889  tempB[i] += QQ(jj,nZ+i) * tempA[j];
4890  }
4891  }
4892  }
4893 
4894  if ( backsolveT( tempB,BT_TRUE,delta_yAC_TMP ) != SUCCESSFUL_RETURN )
4895  return THROWERROR( RET_STEPDIRECTION_FAILED_TQ );
4896  }
4897 
4898  /* refine the solution found so far */
4899  for ( i=0; i<nFR; ++i )
4900  delta_xFR[i] += delta_xFR_TMP[i];
4901  for ( i=0; i<nAC; ++i )
4902  delta_yAC[i] += delta_yAC_TMP[i];
4903 
4904  if ( options.numRefinementSteps > 0 )
4905  {
4906  /* compute residuals in tempA and tempB, and max-norm */
4907  for ( i=0; i<nFR; ++i )
4908  {
4909  ii = FR_idx[i];
4910  tempA[i] = delta_g[ii];
4911  }
4912 
4913  switch ( hessianType )
4914  {
4915  case HST_ZERO:
4916  break;
4917 
4918  case HST_IDENTITY:
4919  for ( i=0; i<nFR; ++i )
4920  tempA[i] += delta_xFR[i];
4921  break;
4922 
4923  default:
4924  H->times(bounds.getFree(), bounds.getFree(), 1, 1.0, delta_xFR, nFR, 1.0, tempA, nFR);
4925  H->times(bounds.getFree(), bounds.getFixed(), 1, 1.0, delta_xFX, nFX, 1.0, tempA, nFR);
4926  break;
4927  }
4928 
4929  A->transTimes(constraints.getActive(), bounds.getFree(), 1, -1.0, delta_yAC, nAC, 1.0, tempA, nFR);
4930  real_t rnrm = 0.0;
4931  for ( i=0; i<nFR; ++i )
4932  if (rnrm < fabs (tempA[i]))
4933  rnrm = fabs (tempA[i]);
4934 
4935  if (!Delta_bC_isZero)
4936  {
4937  for ( i=0; i<nAC; ++i )
4938  {
4939  ii = AC_idx[i];
4940  if ( constraints.getStatus( ii ) == ST_LOWER )
4941  tempB[i] = delta_lbA[ii];
4942  else
4943  tempB[i] = delta_ubA[ii];
4944  }
4945  }
4946  else
4947  {
4948  for ( i=0; i<nAC; ++i )
4949  tempB[i] = 0.0;
4950  }
4951  A->times(constraints.getActive(), bounds.getFree(), 1, -1.0, delta_xFR, nFR, 1.0, tempB, nAC);
4952  A->times(constraints.getActive(), bounds.getFixed(), 1, -1.0, delta_xFX, nFX, 1.0, tempB, nAC);
4953  for ( i=0; i<nAC; ++i )
4954  if (rnrm < fabs (tempB[i]))
4955  rnrm = fabs (tempB[i]);
4956 
4957  /* early termination of residual norm small enough */
4958  if ( rnrm < options.epsIterRef )
4959  break;
4960  }
4961  } /* end of refinement loop for delta_xFRz, delta_xFRy, delta_yAC */
4962 
4963 
4964  /* IV) DETERMINE delta_yFX */
4965  if ( nFX > 0 )
4966  {
4967  for( i=0; i<nFX; ++i )
4968  delta_yFX[i] = delta_g[FX_idx[i]];
4969 
4970  A->transTimes(constraints.getActive(), bounds.getFixed(), 1, -1.0, delta_yAC, nAC, 1.0, delta_yFX, nFX);
4971 
4972  if ( hessianType == HST_ZERO )
4973  {
4974  if ( usingRegularisation( ) == BT_TRUE )
4975  for( i=0; i<nFX; ++i )
4976  delta_yFX[i] += options.epsRegularisation*delta_xFX[i];
4977  }
4978  else if ( hessianType == HST_IDENTITY )
4979  {
4980  for( i=0; i<nFX; ++i )
4981  delta_yFX[i] += delta_xFX[i];
4982  }
4983  else
4984  {
4985  H->times(bounds.getFixed(), bounds.getFree(), 1, 1.0, delta_xFR, nFR, 1.0, delta_yFX, nFX);
4986  H->times(bounds.getFixed(), bounds.getFixed(), 1, 1.0, delta_xFX, nFX, 1.0, delta_yFX, nFX);
4987  }
4988  }
4989 
4990  return SUCCESSFUL_RETURN;
4991 }
4992 
4993 
4994 /*
4995  * p e r f o r m S t e p
4996  */
4997 returnValue QProblem::performStep( const real_t* const delta_g,
4998  const real_t* const delta_lbA, const real_t* const delta_ubA,
4999  const real_t* const delta_lb, const real_t* const delta_ub,
5000  const real_t* const delta_xFX, const real_t* const delta_xFR,
5001  const real_t* const delta_yAC, const real_t* const delta_yFX,
5002  int& BC_idx, SubjectToStatus& BC_status, BooleanType& BC_isBound
5003  )
5004 {
5005  int i, j, ii, jj;
5006  int nV = getNV( );
5007  int nC = getNC( );
5008  int nFR = getNFR( );
5009  int nFX = getNFX( );
5010  int nAC = getNAC( );
5011  int nIAC = getNIAC( );
5012 
5013  int* FR_idx;
5014  int* FX_idx;
5015  int* AC_idx;
5016  int* IAC_idx;
5017 
5018  bounds.getFree( )->getNumberArray( &FR_idx );
5019  bounds.getFixed( )->getNumberArray( &FX_idx );
5020  constraints.getActive( )->getNumberArray( &AC_idx );
5021  constraints.getInactive( )->getNumberArray( &IAC_idx );
5022 
5023 
5024  /* initialise maximum steplength array */
5025  tau = 1.0;
5026  BC_idx = -1;
5027  BC_status = ST_UNDEFINED;
5028 
5029  int BC_idx_tmp = -1;
5030 
5031  real_t* num = new real_t[ getMax( nV,nC ) ];
5032  real_t* den = new real_t[ getMax( nV,nC ) ];
5033 
5034  real_t* delta_Ax_l = new real_t[nC];
5035  real_t* delta_Ax_u = new real_t[nC];
5036  real_t* delta_Ax = new real_t[nC];
5037 
5038  real_t* delta_x = new real_t[nV];
5039  for( j=0; j<nFR; ++j )
5040  {
5041  jj = FR_idx[j];
5042  delta_x[jj] = delta_xFR[j];
5043  }
5044  for( j=0; j<nFX; ++j )
5045  {
5046  jj = FX_idx[j];
5047  delta_x[jj] = delta_xFX[j];
5048  }
5049 
5050 
5051  /* I) DETERMINE MAXIMUM DUAL STEPLENGTH: */
5052  /* 1) Ensure that active dual constraints' bounds remain valid
5053  * (ignoring inequality constraints). */
5054  for( i=0; i<nAC; ++i )
5055  {
5056  ii = AC_idx[i];
5057 
5058  num[i] = y[nV+ii];
5059  den[i] = -delta_yAC[i];
5060  }
5061 
5062  performRatioTest( nAC,AC_idx,&constraints, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
5063 
5064  if ( BC_idx_tmp >= 0 )
5065  {
5066  BC_idx = BC_idx_tmp;
5067  BC_status = ST_INACTIVE;
5068  BC_isBound = BT_FALSE;
5069  }
5070 
5071 
5072  /* 2) Ensure that active dual bounds remain valid
5073  * (ignoring implicitly fixed variables). */
5074  for( i=0; i<nFX; ++i )
5075  {
5076  ii = FX_idx[i];
5077  num[i] = y[ii];
5078  den[i] = -delta_yFX[i];
5079  }
5080 
5081  performRatioTest( nFX,FX_idx,&bounds, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
5082 
5083  if ( BC_idx_tmp >= 0 )
5084  {
5085  BC_idx = BC_idx_tmp;
5086  BC_status = ST_INACTIVE;
5087  BC_isBound = BT_TRUE;
5088  }
5089 
5090 
5091  /* II) DeETERMINE MAXIMUM PRIMAL STEPLENGTH */
5092  /* 1) Ensure that inactive constraints' bounds remain valid
5093  * (ignoring unbounded constraints). */
5094 
5095  #ifdef __MANY_CONSTRAINTS__
5096  real_t delta_x_max = 0.0;
5097  for( i=0; i<nV; ++i )
5098  {
5099  if ( fabs( delta_x[i] ) > delta_x_max )
5100  delta_x_max = fabs( delta_x[i] );
5101  }
5102 
5103  for( i=0; i<nIAC; ++i )
5104  {
5105  ii = IAC_idx[i];
5106 
5107  if ( constraints.getType( ii ) != ST_UNBOUNDED )
5108  {
5109  delta_Ax_l[ii] = -delta_x_max;
5110  if ( delta_lbA[ii] > 0.0 )
5111  delta_Ax_l[ii] -= delta_lbA[ii];
5112 
5113  delta_Ax_u[ii] = -delta_x_max;
5114  if ( delta_ubA[ii] > 0.0 )
5115  delta_Ax_u[ii] -= delta_ubA[ii];
5116 
5117  if ( ( -delta_Ax_l[ii] >= Ax_l[ii] ) || ( -delta_Ax_u[ii] >= Ax_u[ii] ) )
5118  {
5119  /* calculate products A*x and A*deltaX */
5120  if ( constraintProduct == 0 )
5121  {
5122  Ax[ii] = 0.0;
5123  delta_Ax[ii] = 0.0;
5124  for( j=0; j<nV; ++j )
5125  {
5126  Ax[ii] += AA(ii,j) * x[j]; /* Andreas: Do not touch many constraints case yet */
5127  delta_Ax[ii] += AA(ii,j) * delta_x[j];
5128  }
5129  }
5130  else
5131  {
5132  if ( (*constraintProduct)( ii,x, &(Ax[ii]) ) != 0 )
5133  {
5134  delete[] den; delete[] num;
5135  delete[] delta_Ax; delete[] delta_Ax_u; delete[] delta_Ax; delete[] delta_x;
5136  return THROWERROR( RET_ERROR_IN_CONSTRAINTPRODUCT );
5137  }
5138 
5139  if ( (*constraintProduct)( ii,delta_x, &(delta_Ax[ii]) ) != 0 )
5140  {
5141  delete[] den; delete[] num;
5142  delete[] delta_Ax; delete[] delta_Ax_u; delete[] delta_Ax_l; delete[] delta_x;
5143  return THROWERROR( RET_ERROR_IN_CONSTRAINTPRODUCT );
5144  }
5145  }
5146 
5147  Ax_u[ii] = ubA[ii] - Ax[ii];
5148  Ax_l[ii] = Ax[ii] - lbA[ii];
5149 
5150  /* inactive lower constraints' bounds */
5151  if ( constraints.hasNoLower( ) == BT_FALSE )
5152  {
5153  num[i] = getMax( Ax_l[ii],0.0 );
5154  den[i] = delta_lbA[ii] - delta_Ax[ii];
5155 
5156  if ( isBlocking( ii,num[i],den[i], options.epsNum,options.epsDen, tau ) == BT_TRUE )
5157  {
5158  tau = num[i]/den[i];
5159  BC_idx = ii;
5160  BC_status = ST_LOWER;
5161  BC_isBound = BT_FALSE;
5162  }
5163  }
5164  delta_Ax_l[ii] = delta_Ax[ii] - delta_lbA[ii];
5165 
5166  /* inactive upper constraints' bounds */
5167  if ( constraints.hasNoUpper( ) == BT_FALSE )
5168  {
5169  num[i] = getMax( Ax_u[ii],0.0 );
5170  den[i] = delta_Ax[ii] - delta_ubA[ii];
5171 
5172  if ( isBlocking( ii,num[i],den[i], options.epsNum,options.epsDen, tau ) == BT_TRUE )
5173  {
5174  tau = num[i]/den[i];
5175  BC_status = ST_UPPER;
5176  BC_isBound = BT_FALSE;
5177  }
5178  }
5179  delta_Ax_u[ii] = delta_ubA[ii] - delta_Ax[ii];
5180  }
5181  }
5182  }
5183  #else /* NOT __MANY_CONSTRAINTS__ */
5184  /* calculate product A*x */
5185  if ( constraintProduct == 0 )
5186  {
5187  A->times(constraints.getInactive(), 0, 1, 1.0, delta_x, nV, 0.0, delta_Ax, nC, BT_FALSE);
5188  }
5189  else
5190  {
5191  for( i=0; i<nIAC; ++i )
5192  {
5193  ii = IAC_idx[i];
5194 
5195  if ( constraints.getType( ii ) != ST_UNBOUNDED )
5196  {
5197  if ( (*constraintProduct)( ii,delta_x, &(delta_Ax[ii]) ) != 0 )
5198  {
5199  delete[] den; delete[] num;
5200  delete[] delta_Ax; delete[] delta_Ax_u; delete[] delta_Ax_l; delete[] delta_x;
5201  return THROWERROR( RET_ERROR_IN_CONSTRAINTPRODUCT );
5202  }
5203  }
5204  }
5205  }
5206 
5207  if ( constraints.hasNoLower( ) == BT_FALSE )
5208  {
5209  for( i=0; i<nIAC; ++i )
5210  {
5211  ii = IAC_idx[i];
5212  num[i] = getMax( Ax_l[ii],0.0 );
5213  den[i] = delta_lbA[ii] - delta_Ax[ii];
5214  }
5215 
5216  performRatioTest( nIAC,IAC_idx,&constraints, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
5217 
5218  if ( BC_idx_tmp >= 0 )
5219  {
5220  BC_idx = BC_idx_tmp;
5221  BC_status = ST_LOWER;
5222  BC_isBound = BT_FALSE;
5223  }
5224  }
5225 
5226  if ( constraints.hasNoUpper( ) == BT_FALSE )
5227  {
5228  for( i=0; i<nIAC; ++i )
5229  {
5230  ii = IAC_idx[i];
5231  num[i] = getMax( Ax_u[ii],0.0 );
5232  den[i] = delta_Ax[ii] - delta_ubA[ii];
5233  }
5234 
5235  performRatioTest( nIAC,IAC_idx,&constraints, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
5236 
5237  if ( BC_idx_tmp >= 0 )
5238  {
5239  BC_idx = BC_idx_tmp;
5240  BC_status = ST_UPPER;
5241  BC_isBound = BT_FALSE;
5242  }
5243  }
5244 
5245 
5246  for( i=0; i<nIAC; ++i )
5247  {
5248  ii = IAC_idx[i];
5249 
5250  if ( constraints.getType( ii ) != ST_UNBOUNDED )
5251  {
5252  delta_Ax_l[ii] = delta_Ax[ii] - delta_lbA[ii];
5253  delta_Ax_u[ii] = delta_ubA[ii] - delta_Ax[ii];
5254  }
5255  }
5256  #endif /* __MANY_CONSTRAINTS__ */
5257 
5258 
5259  /* 2) Ensure that inactive bounds remain valid
5260  * (ignoring unbounded variables). */
5261  /* inactive lower bounds */
5262  if ( bounds.hasNoLower( ) == BT_FALSE )
5263  {
5264  for( i=0; i<nFR; ++i )
5265  {
5266  ii = FR_idx[i];
5267  num[i] = getMax( x[ii] - lb[ii],0.0 );
5268  den[i] = delta_lb[ii] - delta_xFR[i];
5269  }
5270 
5271  performRatioTest( nFR,FR_idx,&bounds, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
5272 
5273  if ( BC_idx_tmp >= 0 )
5274  {
5275  BC_idx = BC_idx_tmp;
5276  BC_status = ST_LOWER;
5277  BC_isBound = BT_TRUE;
5278  }
5279  }
5280 
5281  /* inactive upper bounds */
5282  if ( bounds.hasNoUpper( ) == BT_FALSE )
5283  {
5284  for( i=0; i<nFR; ++i )
5285  {
5286  ii = FR_idx[i];
5287  num[i] = getMax( ub[ii] - x[ii],0.0 );
5288  den[i] = delta_xFR[i] - delta_ub[ii];
5289  }
5290 
5291  performRatioTest( nFR,FR_idx,&bounds, num,den, options.epsNum,options.epsDen, tau,BC_idx_tmp );
5292 
5293  if ( BC_idx_tmp >= 0 )
5294  {
5295  BC_idx = BC_idx_tmp;
5296  BC_status = ST_UPPER;
5297  BC_isBound = BT_TRUE;
5298  }
5299  }
5300 
5301  delete[] den;
5302  delete[] num;
5303  delete[] delta_x;
5304 
5305 
5306  #ifndef __XPCTARGET__
5307  char messageString[80];
5308 
5309  if ( BC_status == ST_UNDEFINED )
5310  snprintf( messageString,80,"Stepsize is %.16e!",tau );
5311  else
5312  snprintf( messageString,80,"Stepsize is %.16e! (BC_idx = %d, BC_isBound = %d, BC_status = %d)",tau,BC_idx,BC_isBound,BC_status );
5313 
5314  getGlobalMessageHandler( )->throwInfo( RET_STEPSIZE_NONPOSITIVE,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
5315  #endif
5316 
5317 
5318  /* III) PERFORM STEP ALONG HOMOTOPY PATH */
5319  if ( tau > ZERO )
5320  {
5321  /* 1) Perform step in primal und dual space... */
5322  for( i=0; i<nFR; ++i )
5323  {
5324  ii = FR_idx[i];
5325  x[ii] += tau*delta_xFR[i];
5326  }
5327 
5328  for( i=0; i<nFX; ++i )
5329  {
5330  ii = FX_idx[i];
5331  x[ii] += tau*delta_xFX[i];
5332  y[ii] += tau*delta_yFX[i];
5333  }
5334 
5335  for( i=0; i<nAC; ++i )
5336  {
5337  ii = AC_idx[i];
5338  y[nV+ii] += tau*delta_yAC[i];
5339  }
5340 
5341  /* 2) Shift QP data. */
5342  for( i=0; i<nV; ++i )
5343  {
5344  g[i] += tau*delta_g[i];
5345  lb[i] += tau*delta_lb[i];
5346  ub[i] += tau*delta_ub[i];
5347  }
5348 
5349  for( i=0; i<nC; ++i )
5350  {
5351  lbA[i] += tau*delta_lbA[i];
5352  ubA[i] += tau*delta_ubA[i];
5353  }
5354 
5355  /* recompute Ax. */
5356  /*A->times(1, 1.0, x, nV, 0.0, Ax, nC);
5357  for( i=0; i<nC; ++i )
5358  {
5359  Ax_u[i] = ubA[i] - Ax[i];
5360  Ax_l[i] = Ax[i] - lbA[i];
5361  }*/
5362 
5363  A->times( constraints.getActive(),0, 1, 1.0, x, nV, 0.0, Ax, nC, BT_FALSE );
5364  for( i=0; i<nAC; ++i )
5365  {
5366  ii = AC_idx[i];
5367  Ax_u[ii] = ubA[ii] - Ax[ii];
5368  Ax_l[ii] = Ax[ii] - lbA[ii];
5369  }
5370  for( i=0; i<nIAC; ++i )
5371  {
5372  ii = IAC_idx[i];
5373  if ( constraints.getType( ii ) != ST_UNBOUNDED )
5374  {
5375  Ax[ii] += tau*delta_Ax[ii];
5376  Ax_l[ii] += tau*delta_Ax_l[ii];
5377  Ax_u[ii] += tau*delta_Ax_u[ii];
5378  }
5379  }
5380  }
5381  else
5382  {
5383  /* print a stepsize warning if stepsize is zero */
5384  #ifndef __XPCTARGET__
5385  snprintf( messageString,80,"Stepsize is %.16e",tau );
5386  getGlobalMessageHandler( )->throwWarning( RET_STEPSIZE,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
5387  #endif
5388  }
5389 
5390  delete[] delta_Ax; delete[] delta_Ax_u; delete[] delta_Ax_l;
5391 
5392  return SUCCESSFUL_RETURN;
5393 }
5394 
5395 
5396 /*
5397  * c h a n g e A c t i v e S e t
5398  */
5399 returnValue QProblem::changeActiveSet( int BC_idx, SubjectToStatus BC_status, BooleanType BC_isBound )
5400 {
5401  int nV = getNV( );
5402 
5403  char messageString[80];
5404 
5405  switch ( BC_status )
5406  {
5407  /* Optimal solution found as no working set change detected. */
5408  case ST_UNDEFINED:
5409  return RET_OPTIMAL_SOLUTION_FOUND;
5410 
5411  /* Remove one variable from active set. */
5412  case ST_INACTIVE:
5413  if ( BC_isBound == BT_TRUE )
5414  {
5415  #ifndef __XPCTARGET__
5416  snprintf( messageString,80,"bound no. %d.", BC_idx );
5417  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
5418  #endif
5419 
5420  if ( removeBound( BC_idx,BT_TRUE,BT_TRUE,options.enableNZCTests ) != SUCCESSFUL_RETURN )
5421  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
5422 
5423  y[BC_idx] = 0.0;
5424  }
5425  else
5426  {
5427  #ifndef __XPCTARGET__
5428  snprintf( messageString,80,"constraint no. %d.", BC_idx );
5429  getGlobalMessageHandler( )->throwInfo( RET_REMOVE_FROM_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
5430  #endif
5431 
5432  if ( removeConstraint( BC_idx,BT_TRUE,BT_TRUE,options.enableNZCTests ) != SUCCESSFUL_RETURN )
5433  return THROWERROR( RET_REMOVE_FROM_ACTIVESET_FAILED );
5434 
5435  y[nV+BC_idx] = 0.0;
5436  }
5437  break;
5438 
5439 
5440  /* Add one variable to active set. */
5441  default:
5442  returnValue returnvalue;
5443  if ( BC_isBound == BT_TRUE )
5444  {
5445  #ifndef __XPCTARGET__
5446  if ( BC_status == ST_LOWER )
5447  snprintf( messageString,80,"lower bound no. %d.", BC_idx );
5448  else
5449  snprintf( messageString,80,"upper bound no. %d.", BC_idx );
5450  getGlobalMessageHandler( )->throwInfo( RET_ADD_TO_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
5451  #endif
5452 
5453  returnvalue = addBound( BC_idx,BC_status,BT_TRUE );
5454  if ( returnvalue == RET_ADDBOUND_FAILED_INFEASIBILITY )
5455  return returnvalue;
5456  if ( returnvalue != SUCCESSFUL_RETURN )
5457  return THROWERROR( RET_ADD_TO_ACTIVESET_FAILED );
5458  }
5459  else
5460  {
5461  #ifndef __XPCTARGET__
5462  if ( BC_status == ST_LOWER )
5463  snprintf( messageString,80,"lower constraint's bound no. %d.", BC_idx );
5464  else
5465  snprintf( messageString,80,"upper constraint's bound no. %d.", BC_idx );
5466  getGlobalMessageHandler( )->throwInfo( RET_ADD_TO_ACTIVESET,messageString,__FUNCTION__,__FILE__,__LINE__,VS_VISIBLE );
5467  #endif
5468 
5469  returnvalue = addConstraint( BC_idx,BC_status,BT_TRUE );
5470  if ( returnvalue == RET_ADDCONSTRAINT_FAILED_INFEASIBILITY )
5471  return returnvalue;
5472  if ( returnvalue != SUCCESSFUL_RETURN )
5473  return THROWERROR( RET_ADD_TO_ACTIVESET_FAILED );
5474  }
5475  }
5476 
5477  return SUCCESSFUL_RETURN;
5478 }
5479 
5480 
5481 
5482 /*
5483  * r e l a t i v e H o m o t o p y L e n g t h
5484  */
5485 real_t QProblem::relativeHomotopyLength( const real_t* const g_new, const real_t* const lb_new, const real_t* const ub_new,
5486  const real_t* const lbA_new, const real_t* const ubA_new)
5487 {
5488  int nC = getNC( ), i;
5489  real_t len = QProblemB::relativeHomotopyLength(g_new, lb_new, ub_new);
5490  real_t d, s;
5491 
5492  /* lower constraint bounds */
5493  for (i = 0; i < nC && lbA_new; i++)
5494  {
5495  s = fabs(lbA_new[i]);
5496  if (s < 1.0) s = 1.0;
5497  d = fabs(lbA_new[i] - lbA[i]) / s;
5498  if (d > len) len = d;
5499  }
5500 
5501  /* upper constraint bounds */
5502  for (i = 0; i < nC && ubA_new; i++)
5503  {
5504  s = fabs(ubA_new[i]);
5505  if (s < 1.0) s = 1.0;
5506  d = fabs(ubA_new[i] - ubA[i]) / s;
5507  if (d > len) len = d;
5508  }
5509 
5510  return len;
5511 }
5512 
5513 
5514 /*
5515  * p e r f o r m R a m p i n g
5516  */
5517 returnValue QProblem::performRamping( )
5518 {
5519  int nV = getNV( ), nC = getNC( ), bstat, cstat, i;
5520  real_t t, rampVal;
5521 
5522  /* ramp inactive variable bounds and active dual bound variables */
5523  for (i = 0; i < nV; i++)
5524  {
5525  switch (bounds.getType(i))
5526  {
5527  case ST_EQUALITY: lb[i] = x[i]; ub[i] = x[i]; continue; /* reestablish exact feasibility */
5528  case ST_UNBOUNDED: continue;
5529  default: break;
5530  }
5531 
5532  t = static_cast<real_t>(i) / static_cast<real_t>(nV+nC-1);
5533  rampVal = (1.0-t) * ramp0 + t * ramp1;
5534  bstat = bounds.getStatus(i);
5535  if (bstat != ST_LOWER) { lb[i] = x[i] - rampVal; }
5536  if (bstat != ST_UPPER) { ub[i] = x[i] + rampVal; }
5537  if (bstat == ST_LOWER) { lb[i] = x[i]; y[i] = +rampVal; }
5538  if (bstat == ST_UPPER) { ub[i] = x[i]; y[i] = -rampVal; }
5539  if (bstat == ST_INACTIVE) y[i] = 0.0; /* reestablish exact complementarity */
5540  }
5541 
5542  /* ramp inactive constraints and active dual constraint variables */
5543  for (i = 0; i < nC; i++)
5544  {
5545  switch (constraints.getType(i))
5546  {
5547  case ST_EQUALITY: lbA[i] = Ax[i]; ubA[i] = Ax[i]; continue; /* reestablish exact feasibility */
5548  case ST_BOUNDED: continue;
5549  default: break;
5550  }
5551 
5552  t = static_cast<real_t>(nV+i) / static_cast<real_t>(nV+nC-1);
5553  rampVal = (1.0-t) * ramp0 + t * ramp1;
5554  cstat = bounds.getStatus(i);
5555  if (cstat != ST_LOWER) { lbA[i] = Ax[i] - rampVal; }
5556  if (cstat != ST_UPPER) { ubA[i] = Ax[i] + rampVal; }
5557  if (cstat == ST_LOWER) { lbA[i] = Ax[i]; y[nV+i] = +rampVal; }
5558  if (cstat == ST_UPPER) { ubA[i] = Ax[i]; y[nV+i] = -rampVal; }
5559  if (cstat == ST_INACTIVE) y[nV+i] = 0.0; /* reestablish exact complementarity */
5560  }
5561 
5562  /* reestablish exact stationarity */
5563  setupAuxiliaryQPgradient( );
5564 
5565  /* swap ramp direction */
5566  t = ramp0; ramp0 = ramp1; ramp1 = t;
5567 
5568  return SUCCESSFUL_RETURN;
5569 }
5570 
5571 
5572 /*
5573  * p e r f o r m D r i f t C o r r e c t i o n
5574  */
5575 returnValue QProblem::performDriftCorrection( )
5576 {
5577  int i;
5578  int nV = getNV ();
5579  int nC = getNC ();
5580 
5581  for ( i=0; i<nV; ++i )
5582  {
5583  switch ( bounds.getType ( i ) )
5584  {
5585  case ST_BOUNDED:
5586  switch ( bounds.getStatus ( i ) )
5587  {
5588  case ST_LOWER:
5589  lb[i] = x[i];
5590  ub[i] = getMax (ub[i], x[i]);
5591  y[i] = getMax (y[i], 0.0);
5592  break;
5593  case ST_UPPER:
5594  lb[i] = getMin (lb[i], x[i]);
5595  ub[i] = x[i];
5596  y[i] = getMin (y[i], 0.0);
5597  break;
5598  case ST_INACTIVE:
5599  lb[i] = getMin (lb[i], x[i]);
5600  ub[i] = getMax (ub[i], x[i]);
5601  y[i] = 0.0;
5602  break;
5603  case ST_UNDEFINED:
5604  break;
5605  }
5606  break;
5607  case ST_EQUALITY:
5608  lb[i] = x[i];
5609  ub[i] = x[i];
5610  break;
5611  case ST_UNBOUNDED:
5612  case ST_UNKNOWN:
5613  break;
5614  }
5615  }
5616 
5617  for ( i=0; i<nC; ++i )
5618  {
5619  switch ( constraints.getType ( i ) )
5620  {
5621  case ST_BOUNDED:
5622  switch ( constraints.getStatus ( i ) )
5623  {
5624  case ST_LOWER:
5625  lbA[i] = Ax[i];
5626  Ax_l[i] = 0.0;
5627  ubA[i] = getMax (ubA[i], Ax[i]);
5628  Ax_u[i] = ubA[i] - Ax[i];
5629  y[i+nV] = getMax (y[i+nV], 0.0);
5630  break;
5631  case ST_UPPER:
5632  lbA[i] = getMin (lbA[i], Ax[i]);
5633  Ax_l[i] = Ax[i] - lbA[i];
5634  ubA[i] = Ax[i];
5635  Ax_u[i] = 0.0;
5636  y[i+nV] = getMin (y[i+nV], 0.0);
5637  break;
5638  case ST_INACTIVE:
5639  lbA[i] = getMin (lbA[i], Ax[i]);
5640  Ax_l[i] = Ax[i] - lbA[i];
5641  ubA[i] = getMax (ubA[i], Ax[i]);
5642  Ax_u[i] = ubA[i] - Ax[i];
5643  y[i+nV] = 0.0;
5644  break;
5645  case ST_UNDEFINED:
5646  break;
5647  }
5648  break;
5649  case ST_EQUALITY:
5650  lbA[i] = Ax[i];
5651  Ax_l[i] = 0.0;
5652  ubA[i] = Ax[i];
5653  Ax_u[i] = 0.0;
5654  break;
5655  case ST_UNBOUNDED:
5656  case ST_UNKNOWN:
5657  break;
5658  }
5659  }
5660 
5661  setupAuxiliaryQPgradient( );
5662 
5663  return SUCCESSFUL_RETURN;
5664 }
5665 
5666 /*
5667  * s e t u p A u x i l i a r y Q P
5668  */
5669 returnValue QProblem::setupAuxiliaryQP( const Bounds* const guessedBounds, const Constraints* const guessedConstraints )
5670 {
5671  int i, j;
5672  int nV = getNV( );
5673  int nC = getNC( );
5674 
5675  /* consistency check */
5676  if ( ( guessedBounds == 0 ) || ( guessedConstraints == 0 ) )
5677  return THROWERROR( RET_INVALID_ARGUMENTS );
5678 
5679  /* nothing to do */
5680  if ( ( guessedBounds == &bounds ) && ( guessedConstraints == &constraints ) )
5681  return SUCCESSFUL_RETURN;
5682 
5683  status = QPS_PREPARINGAUXILIARYQP;
5684 
5685 
5686  /* I) SETUP WORKING SET ... */
5687  if ( shallRefactorise( guessedBounds,guessedConstraints ) == BT_TRUE )
5688  {
5689  /* ... WITH REFACTORISATION: */
5690  /* 1) Reset bounds/constraints ... */
5691  bounds.init( nV );
5692  constraints.init( nC );
5693 
5694  /* ... and set them up afresh. */
5695  if ( setupSubjectToType( ) != SUCCESSFUL_RETURN )
5696  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5697 
5698  if ( bounds.setupAllFree( ) != SUCCESSFUL_RETURN )
5699  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5700 
5701  if ( constraints.setupAllInactive( ) != SUCCESSFUL_RETURN )
5702  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5703 
5704  /* 2) Setup TQ factorisation. */
5705  if ( setupTQfactorisation( ) != SUCCESSFUL_RETURN )
5706  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5707 
5708  /* 3) Setup guessed working sets afresh. */
5709  if ( setupAuxiliaryWorkingSet( guessedBounds,guessedConstraints,BT_TRUE ) != SUCCESSFUL_RETURN )
5710  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5711 
5712  /* 4) Calculate Cholesky decomposition. */
5713  if ( ( getNAC( ) + getNFX( ) ) == 0 )
5714  {
5715  /* Factorise full Hessian if no bounds/constraints are active. */
5716  if ( setupCholeskyDecomposition( ) != SUCCESSFUL_RETURN )
5717  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5718  }
5719  else
5720  {
5721  /* Factorise projected Hessian if there active bounds/constraints. */
5722  if ( setupCholeskyDecompositionProjected( ) != SUCCESSFUL_RETURN )
5723  return THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5724  }
5725  }
5726  else
5727  {
5728  /* ... WITHOUT REFACTORISATION: */
5729  if ( setupAuxiliaryWorkingSet( guessedBounds,guessedConstraints,BT_FALSE ) != SUCCESSFUL_RETURN )
5730  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5731  }
5732 
5733 
5734  /* II) SETUP AUXILIARY QP DATA: */
5735  /* 1) Ensure that dual variable is zero for fixed bounds and active constraints. */
5736  for ( i=0; i<nV; ++i )
5737  if ( bounds.getStatus( i ) != ST_INACTIVE )
5738  y[i] = 0.0;
5739 
5740  for ( i=0; i<nC; ++i )
5741  if ( constraints.getStatus( i ) != ST_INACTIVE )
5742  y[nV+i] = 0.0;
5743 
5744  /* 2) Setup gradient and (constraints') bound vectors. */
5745  if ( setupAuxiliaryQPgradient( ) != SUCCESSFUL_RETURN )
5746  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5747 
5748  A->times(1, 1.0, x, nV, 0.0, Ax, nC);
5749  for ( j=0; j<nC; ++j )
5750  {
5751  Ax_l[j] = Ax[j];
5752  Ax_u[j] = Ax[j];
5753  }
5754 
5755  /* (also sets Ax_l and Ax_u) */
5756  if ( setupAuxiliaryQPbounds( 0,0,BT_FALSE ) != SUCCESSFUL_RETURN )
5757  THROWERROR( RET_SETUP_AUXILIARYQP_FAILED );
5758 
5759  return SUCCESSFUL_RETURN;
5760 }
5761 
5762 
5763 /*
5764  * s h a l l R e f a c t o r i s e
5765  */
5766 
5767 BooleanType QProblem::shallRefactorise( const Bounds* const guessedBounds,
5768  const Constraints* const guessedConstraints
5769  ) const
5770 {
5771  int i;
5772  int nV = getNV( );
5773  int nC = getNC( );
5774 
5775  /* always refactorise if Hessian is not known to be positive definite */
5776  if ( getHessianType( ) == HST_SEMIDEF )
5777  return BT_TRUE;
5778 
5779  /* 1) Determine number of bounds that have same status
5780  * in guessed AND current bounds.*/
5781  int differenceNumberBounds = 0;
5782 
5783  for( i=0; i<nV; ++i )
5784  if ( guessedBounds->getStatus( i ) != bounds.getStatus( i ) )
5785  ++differenceNumberBounds;
5786 
5787  /* 2) Determine number of constraints that have same status
5788  * in guessed AND current constraints.*/
5789  int differenceNumberConstraints = 0;
5790 
5791  for( i=0; i<nC; ++i )
5792  if ( guessedConstraints->getStatus( i ) != constraints.getStatus( i ) )
5793  ++differenceNumberConstraints;
5794 
5795  /* 3) Decide wheter to refactorise or not. */
5796  if ( 2*(differenceNumberBounds+differenceNumberConstraints) > guessedConstraints->getNAC( )+guessedBounds->getNFX( ) )
5797  return BT_TRUE;
5798  else
5799  return BT_FALSE;
5800 }
5801 
5802 
5803 /*
5804  * s e t u p Q P d a t a
5805  */
5806 returnValue QProblem::setupQPdata( SymmetricMatrix *_H, const real_t* const _g, Matrix *_A,
5807  const real_t* const _lb, const real_t* const _ub,
5808  const real_t* const _lbA, const real_t* const _ubA
5809  )
5810 {
5811  int i;
5812  int nC = getNC( );
5813 
5814 
5815  /* 1) Load Hessian matrix as well as lower and upper bounds vectors. */
5816  if ( QProblemB::setupQPdata( _H,_g,_lb,_ub ) != SUCCESSFUL_RETURN )
5817  return THROWERROR( RET_INVALID_ARGUMENTS );
5818 
5819  /* 2) Load constraint matrix. */
5820  if ( ( nC > 0 ) && ( _A == 0 ) )
5821  return THROWERROR( RET_INVALID_ARGUMENTS );
5822 
5823  if ( nC > 0 )
5824  {
5825  setA( _A );
5826 
5827  /* 3) Setup lower constraints' bounds vector. */
5828  if ( _lbA != 0 )
5829  {
5830  setLBA( _lbA );
5831  }
5832  else
5833  {
5834  /* if no lower constraints' bounds are specified, set them to -infinity */
5835  for( i=0; i<nC; ++i )
5836  lbA[i] = -INFTY;
5837  }
5838 
5839  /* 4) Setup upper constraints' bounds vector. */
5840  if ( _ubA != 0 )
5841  {
5842  setUBA( _ubA );
5843  }
5844  else
5845  {
5846  /* if no upper constraints' bounds are specified, set them to infinity */
5847  for( i=0; i<nC; ++i )
5848  ubA[i] = INFTY;
5849  }
5850  }
5851 
5852  return SUCCESSFUL_RETURN;
5853 }
5854 
5855 
5856 /*
5857  * s e t u p Q P d a t a
5858  */
5859 returnValue QProblem::setupQPdata( const real_t* const _H, const real_t* const _g, const real_t* const _A,
5860  const real_t* const _lb, const real_t* const _ub,
5861  const real_t* const _lbA, const real_t* const _ubA
5862  )
5863 {
5864  int i;
5865  int nC = getNC( );
5866 
5867 
5868  /* 1) Load Hessian matrix as well as lower and upper bounds vectors. */
5869  if ( QProblemB::setupQPdata( _H,_g,_lb,_ub ) != SUCCESSFUL_RETURN )
5870  return THROWERROR( RET_INVALID_ARGUMENTS );
5871 
5872  /* 2) Load constraint matrix. */
5873  if ( ( nC > 0 ) && ( _A == 0 ) )
5874  return THROWERROR( RET_INVALID_ARGUMENTS );
5875 
5876  if ( nC > 0 )
5877  {
5878  setA( _A );
5879 
5880  /* 3) Setup lower constraints' bounds vector. */
5881  if ( _lbA != 0 )
5882  {
5883  setLBA( _lbA );
5884  }
5885  else
5886  {
5887  /* if no lower constraints' bounds are specified, set them to -infinity */
5888  for( i=0; i<nC; ++i )
5889  lbA[i] = -INFTY;
5890  }
5891 
5892  /* 4) Setup upper constraints' bounds vector. */
5893  if ( _ubA != 0 )
5894  {
5895  setUBA( _ubA );
5896  }
5897  else
5898  {
5899  /* if no upper constraints' bounds are specified, set them to infinity */
5900  for( i=0; i<nC; ++i )
5901  ubA[i] = INFTY;
5902  }
5903  }
5904 
5905  return SUCCESSFUL_RETURN;
5906 }
5907 
5908 
5909 /*
5910  * s e t u p Q P d a t a F r o m F i l e
5911  */
5912 returnValue QProblem::setupQPdataFromFile( const char* const H_file, const char* const g_file, const char* const A_file,
5913  const char* const lb_file, const char* const ub_file,
5914  const char* const lbA_file, const char* const ubA_file
5915  )
5916 {
5917  int i;
5918  int nV = getNV( );
5919  int nC = getNC( );
5920 
5921  returnValue returnvalue;
5922 
5923 
5924  /* 1) Load Hessian matrix as well as lower and upper bounds vectors from files. */
5925  returnvalue = QProblemB::setupQPdataFromFile( H_file,g_file,lb_file,ub_file );
5926  if ( returnvalue != SUCCESSFUL_RETURN )
5927  return THROWERROR( returnvalue );
5928 
5929  /* 2) Load constraint matrix from file. */
5930  if ( ( nC > 0 ) && ( A_file == 0 ) )
5931  return THROWERROR( RET_INVALID_ARGUMENTS );
5932 
5933  if ( nC > 0 )
5934  {
5935  real_t* _A = new real_t[nC * nV];
5936  returnvalue = readFromFile( _A, nC,nV, A_file );
5937  if ( returnvalue != SUCCESSFUL_RETURN )
5938  {
5939  delete[] _A;
5940  return THROWERROR( returnvalue );
5941  }
5942  setA( _A );
5943  A->doFreeMemory( );
5944 
5945  /* 3) Load lower constraints' bounds vector from file. */
5946  if ( lbA_file != 0 )
5947  {
5948  returnvalue = readFromFile( lbA, nC, lbA_file );
5949  if ( returnvalue != SUCCESSFUL_RETURN )
5950  return THROWERROR( returnvalue );
5951  }
5952  else
5953  {
5954  /* if no lower constraints' bounds are specified, set them to -infinity */
5955  for( i=0; i<nC; ++i )
5956  lbA[i] = -INFTY;
5957  }
5958 
5959  /* 4) Load upper constraints' bounds vector from file. */
5960  if ( ubA_file != 0 )
5961  {
5962  returnvalue = readFromFile( ubA, nC, ubA_file );
5963  if ( returnvalue != SUCCESSFUL_RETURN )
5964  return THROWERROR( returnvalue );
5965  }
5966  else
5967  {
5968  /* if no upper constraints' bounds are specified, set them to infinity */
5969  for( i=0; i<nC; ++i )
5970  ubA[i] = INFTY;
5971  }
5972  }
5973 
5974  return SUCCESSFUL_RETURN;
5975 }
5976 
5977 
5978 /*
5979  * l o a d Q P v e c t o r s F r o m F i l e
5980  */
5981 returnValue QProblem::loadQPvectorsFromFile( const char* const g_file, const char* const lb_file, const char* const ub_file,
5982  const char* const lbA_file, const char* const ubA_file,
5983  real_t* const g_new, real_t* const lb_new, real_t* const ub_new,
5984  real_t* const lbA_new, real_t* const ubA_new
5985  ) const
5986 {
5987  int nC = getNC( );
5988 
5989  returnValue returnvalue;
5990 
5991 
5992  /* 1) Load gradient vector as well as lower and upper bounds vectors from files. */
5993  returnvalue = QProblemB::loadQPvectorsFromFile( g_file,lb_file,ub_file, g_new,lb_new,ub_new );
5994  if ( returnvalue != SUCCESSFUL_RETURN )
5995  return THROWERROR( returnvalue );
5996 
5997  if ( nC > 0 )
5998  {
5999  /* 2) Load lower constraints' bounds vector from file. */
6000  if ( lbA_file != 0 )
6001  {
6002  if ( lbA_new != 0 )
6003  {
6004  returnvalue = readFromFile( lbA_new, nC, lbA_file );
6005  if ( returnvalue != SUCCESSFUL_RETURN )
6006  return THROWERROR( returnvalue );
6007  }
6008  else
6009  {
6010  /* If filename is given, storage must be provided! */
6011  return THROWERROR( RET_INVALID_ARGUMENTS );
6012  }
6013  }
6014 
6015  /* 3) Load upper constraints' bounds vector from file. */
6016  if ( ubA_file != 0 )
6017  {
6018  if ( ubA_new != 0 )
6019  {
6020  returnvalue = readFromFile( ubA_new, nC, ubA_file );
6021  if ( returnvalue != SUCCESSFUL_RETURN )
6022  return THROWERROR( returnvalue );
6023  }
6024  else
6025  {
6026  /* If filename is given, storage must be provided! */
6027  return THROWERROR( RET_INVALID_ARGUMENTS );
6028  }
6029  }
6030  }
6031 
6032  return SUCCESSFUL_RETURN;
6033 }
6034 
6035 
6036 /*
6037  * p r i n t I t e r a t i o n
6038  */
6039 returnValue QProblem::printIteration( int iteration,
6040  int BC_idx, SubjectToStatus BC_status, BooleanType BC_isBound
6041  )
6042 {
6043  #ifndef __XPCTARGET__
6044  char myPrintfString[80];
6045  char info[8];
6046 
6047  /* consistency check */
6048  if ( iteration < 0 )
6049  return THROWERROR( RET_INVALID_ARGUMENTS );
6050 
6051  /* nothing to do */
6052  if ( options.printLevel != PL_MEDIUM )
6053  return SUCCESSFUL_RETURN;
6054 
6055 
6056  /* 1) Print header at first iteration. */
6057  if ( iteration == 0 )
6058  {
6059  snprintf( myPrintfString,80,"\n\n#################### qpOASES -- QP NO. %3.0d #####################\n\n", count );
6060  myPrintf( myPrintfString );
6061 
6062  myPrintf( " Iter | StepLength | Info | nFX | nAC \n" );
6063  myPrintf( " ----------+------------------+------------------+---------+--------- \n" );
6064  }
6065 
6066  /* 2) Print iteration line. */
6067  if ( BC_status == ST_UNDEFINED )
6068  {
6069  if ( hessianType == HST_ZERO )
6070  snprintf( info,3,"LP" );
6071  else
6072  snprintf( info,3,"QP" );
6073 
6074  snprintf( myPrintfString,80," %5.1d | %1.6e | %s SOLVED | %4.1d | %4.1d \n", iteration,tau,info,getNFX( ),getNAC( ) );
6075  myPrintf( myPrintfString );
6076  }
6077  else
6078  {
6079  if ( BC_status == ST_INACTIVE )
6080  snprintf( info,5,"REM " );
6081  else
6082  snprintf( info,5,"ADD " );
6083 
6084  if ( BC_isBound == BT_TRUE )
6085  snprintf( &(info[4]),4,"BND" );
6086  else
6087  snprintf( &(info[4]),4,"CON" );
6088 
6089  snprintf( myPrintfString,80," %5.1d | %1.6e | %s %4.1d | %4.1d | %4.1d \n", iteration,tau,info,BC_idx,getNFX( ),getNAC( ) );
6090  myPrintf( myPrintfString );
6091  }
6092  #endif
6093 
6094  return SUCCESSFUL_RETURN;
6095 }
6096 
6097 
6098 END_NAMESPACE_QPOASES
6099 
6100 
6101 /*
6102  * end of file
6103  */
RET_REMOVEBOUND_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:146
getMax
int getMax(int x, int y)
RET_STEPSIZE_NONPOSITIVE
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:139
RET_LINEARLY_DEPENDENT
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:153
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
QProblemB::getHessianType
HessianType getHessianType() const
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
QProblem::setupQPdataFromFile
returnValue setupQPdataFromFile(const char *const H_file, const char *const g_file, const char *const A_file, const char *const lb_file, const char *const ub_file, const char *const lbA_file, const char *const ubA_file)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:5912
QProblem::solveCurrentEQP
returnValue solveCurrentEQP(const int n_rhs, const real_t *g_in, const real_t *lb_in, const real_t *ub_in, const real_t *lbA_in, const real_t *ubA_in, real_t *x_out, real_t *y_out)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:1045
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
RET_UNABLE_TO_READ_FILE
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:172
ConstraintProduct
Interface for specifying user-defined evaluations of constraint products.
Definition: qpOASES-3.0beta/include/qpOASES/ConstraintProduct.hpp:57
QProblemB::isInfeasible
BooleanType isInfeasible() const
RET_UNKNOWN_BUG
Definition: acado_types.hpp:845
ST_LOWER
#define ST_LOWER
Definition: qpOASES_wrapper.h:73
RET_STEPLENGTH_DETERMINATION_FAILED
Definition: acado_types.hpp:1096
Bounds::getFree
Indexlist * getFree()
RET_HESSIAN_NOT_SPD
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:164
QProblem::shallRefactorise
BooleanType shallRefactorise(const Bounds *const guessedBounds, const Constraints *const guessedConstraints) const
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:5767
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
__LINE__
#define __LINE__
Definition: acado_message_handling.hpp:98
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
QProblem::flipper
Flipper flipper
Definition: QProblem.h:101
SUCCESSFUL_RETURN
Definition: acado_types.hpp:836
QProblem::solveRegularisedQP
returnValue solveRegularisedQP(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, int nWSRperformed=0)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:2014
RET_QPOBJECT_NOT_SETUP
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:92
RET_ITERATION_STARTED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:120
QProblem::relativeHomotopyLength
real_t relativeHomotopyLength(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)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:5485
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
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
DenseMatrix::getCol
virtual returnValue getCol(int cNum, const Indexlist *const irows, real_t alpha, real_t *col) const
Definition: qpOASES-3.0beta/src/Matrices.cpp:134
RET_MOVING_BOUND_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:94
RET_HOTSTART_FAILED_AS_QP_NOT_INITIALISED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:119
QProblemB::getNFX
int getNFX()
RET_INVALID_ARGUMENTS
Definition: acado_types.hpp:839
QProblemB::computeGivens
void computeGivens(real_t xold, real_t yold, real_t &xnew, real_t &ynew, real_t &c, real_t &s) const
EPS
BEGIN_NAMESPACE_ACADO const double EPS
Definition: acado_constants.hpp:43
QQ
#define QQ(I, J)
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:135
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
TT
#define TT(I, J)
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:138
MessageHandling::throwInfo
returnValue throwInfo(returnValue Inumber, const char *additionaltext, const char *functionname, const char *filename, const unsigned long linenumber, VisibilityStatus localVisibilityStatus)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/MessageHandling.cpp:315
returnValue
Allows to pass back messages to the calling function.
Definition: acado_types.hpp:1139
RET_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
HST_POSDEF
#define HST_POSDEF
Definition: qpOASES_wrapper.h:66
getCPUtime
real_t getCPUtime()
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/Utils.cpp:433
QProblem::setLBA
returnValue setLBA(const real_t *const lbA_new)
ZERO
const double ZERO
Definition: acado_constants.hpp:62
QProblem::hessianType
HessianType hessianType
Definition: QProblem.h:87
RET_INIT_FAILED_REGULARISATION
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:102
RR
#define RR(I, J)
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:132
QProblem::ensureNonzeroCurvature
returnValue ensureNonzeroCurvature(BooleanType removeBoundNotConstraint, int remIdx, BooleanType &exchangeHappened, BooleanType &addBoundNotConstraint, int &addIdx, SubjectToStatus &addStatus)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:4306
PL_MEDIUM
#define PL_MEDIUM
Definition: qpOASES_wrapper.h:85
QProblemB::getNV
int getNV() const
getMin
real_t getMin(real_t x, real_t y)
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
Constraints::getInactive
Indexlist * getInactive()
QProblemB::setInfeasibilityFlag
returnValue setInfeasibilityFlag(returnValue returnvalue)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2006
RET_INIT_FAILED_UNBOUNDEDNESS
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:105
HessianType
HessianType
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:116
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
Flipper::set
returnValue set(const Bounds *const _bounds, const double *const _R, const Constraints *const _constraints=0, const double *const _Q=0, const double *const _T=0)
Definition: qpOASES-3.0beta/src/Flipper.cpp:176
__FUNCTION__
#define __FUNCTION__
Definition: acado_message_handling.hpp:90
MessageHandling::throwWarning
returnValue throwWarning(returnValue Wnumber, const char *additionaltext, const char *functionname, const char *filename, const unsigned long linenumber, VisibilityStatus localVisibilityStatus)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/MessageHandling.cpp:291
HST_IDENTITY
#define HST_IDENTITY
Definition: qpOASES_wrapper.h:65
RET_ADD_TO_ACTIVESET_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:150
PL_NONE
#define PL_NONE
Definition: qpOASES_wrapper.h:83
RET_REMOVECONSTRAINT_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:145
readFromFile
returnValue readFromFile(real_t *data, int nrow, int ncol, const char *datafilename)
Definition: external_packages/qpOASES-3.0beta/src/Utils.cpp:283
myPrintf
returnValue myPrintf(const char *s)
Definition: external_packages/qpOASES-3.0beta/src/Utils.cpp:242
QProblem::performPlainRatioTest
returnValue performPlainRatioTest(int nIdx, const int *const idxList, const real_t *const num, const real_t *const den, real_t epsNum, real_t epsDen, real_t &t, int &BC_idx) const
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:4284
RET_LI_RESOLVED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:155
QProblem::ramp1
real_t ramp1
Definition: QProblem.h:95
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
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_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
DenseMatrix::diag
virtual real_t diag(int i) const
Definition: qpOASES-3.0beta/src/Matrices.cpp:83
QProblem::changeActiveSet
returnValue changeActiveSet(int BC_idx, SubjectToStatus BC_status, BooleanType BC_isBound)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:5399
PL_HIGH
#define PL_HIGH
Definition: qpOASES_wrapper.h:86
RET_HOTSTART_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:117
ST_UNBOUNDED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:89
QProblemB::usingRegularisation
BooleanType usingRegularisation() const
QProblemB::relativeHomotopyLength
real_t relativeHomotopyLength(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2129
RET_CHOLESKY_OF_ZERO_HESSIAN
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:180
ST_BOUNDED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:90
QProblem::ramp0
real_t ramp0
Definition: QProblem.h:94
QProblemB::isCPUtimeLimitExceeded
BooleanType isCPUtimeLimitExceeded(const real_t *const cputime, real_t starttime, int nWSR) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2021
SubjectTo::getStatus
SubjectToStatus getStatus(int i) const
DenseMatrix
Interfaces matrix-vector operations tailored to general dense matrices.
Definition: qpOASES-3.0beta/include/qpOASES/Matrices.hpp:273
QProblem::setConstraintProduct
returnValue setConstraintProduct(ConstraintProduct *const _constraintProduct)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:1153
QProblem::bounds
Bounds bounds
Definition: QProblem.h:72
QProblemB::performRatioTest
returnValue performRatioTest(int nIdx, const int *const idxList, const SubjectTo *const subjectTo, const real_t *const num, const real_t *const den, real_t epsNum, real_t epsDen, real_t &t, int &BC_idx) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:2081
QProblemB::isBlocking
BooleanType isBlocking(real_t num, real_t den, real_t epsNum, real_t epsDen, real_t &t) const
QProblemB::determineDataShift
returnValue determineDataShift(const real_t *const g_new, const real_t *const lb_new, const real_t *const ub_new, real_t *const delta_g, real_t *const delta_lb, real_t *const delta_ub, BooleanType &Delta_bB_isZero)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1724
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
nu
const int nu
Definition: rocket_with_templates_ocp.hpp:38
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
RET_CONSTRAINTS_ARE_NOT_SCALED
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:181
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()
QProblem::AA
DenseMatrix AA
Definition: QProblem.h:104
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
QProblem::getNEC
int getNEC() const
QProblemB::getNFR
int getNFR()
QPS_PERFORMINGHOMOTOPY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:79
QProblem::determineStepDirection
returnValue determineStepDirection(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: external_packages/qpOASES-3.0beta/src/QProblem.cpp:4666
Matrix
Abstract base class for interfacing tailored matrix-vector operations.
Definition: qpOASES-3.0beta/include/qpOASES/Matrices.hpp:112
QProblemB::setupQPdataFromFile
returnValue setupQPdataFromFile(const char *const H_file, const char *const g_file, const char *const lb_file, const char *const ub_file)
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1873
RET_ADDBOUND_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:143
QProblem::options
Options options
Definition: QProblem.h:98
RET_NO_REGSTEP_NWSR
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:178
SubjectToStatus
SubjectToStatus
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:97
RET_INIT_FAILED_CHOLESKY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:102
QProblem::loadQPvectorsFromFile
returnValue loadQPvectorsFromFile(const char *const g_file, const char *const lb_file, const char *const ub_file, const char *const lbA_file, const char *const ubA_file, real_t *const g_new, real_t *const lb_new, real_t *const ub_new, real_t *const lbA_new, real_t *const ubA_new) const
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:5981
PL_LOW
#define PL_LOW
Definition: qpOASES_wrapper.h:84
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
END_NAMESPACE_QPOASES
#define END_NAMESPACE_QPOASES
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:120
QProblem::solveQP
returnValue solveQP(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, int nWSRperformed=0)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:1694
QProblem::getNC
int getNC() const
RET_QP_ALREADY_INITIALISED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:97
RET_ERROR_IN_CONSTRAINTPRODUCT
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:182
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
QProblem::copy
returnValue copy(const QProblem &rhs)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:1414
QProblem::performStep
returnValue performStep(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, int &BC_idx, SubjectToStatus &BC_status, BooleanType &BC_isBound)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:4997
HST_POSDEF_NULLSPACE
#define HST_POSDEF_NULLSPACE
Definition: qpOASES_wrapper.h:67
RET_ADD_TO_ACTIVESET
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:148
QProblem::setupAuxiliaryQP
virtual returnValue setupAuxiliaryQP(const Bounds *const guessedBounds, const Constraints *const guessedConstraints)
Definition: external_packages/qpOASES-3.0beta/src/QProblem.cpp:5669
QProblemB::isInitialised
BooleanType isInitialised() const
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
QProblem::determineDataShift
returnValue determineDataShift(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: external_packages/qpOASES-3.0beta/src/QProblem.cpp:4600
RET_SETUP_AUXILIARYQP_FAILED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:109
RET_REMOVE_FROM_ACTIVESET
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:147
RET_INIT_SUCCESSFUL
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:106
QProblem::delta_xFR_TMP
real_t delta_xFR_TMP[NVMAX]
Definition: QProblem.h:92
ST_UNKNOWN
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:92
HST_ZERO
#define HST_ZERO
Definition: qpOASES_wrapper.h:64
Bounds::getNUV
int getNUV() const
DenseMatrix::times
returnValue times(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
Definition: qpOASES-3.0beta/src/Matrices.cpp:151
THROWWARNING
#define THROWWARNING(retval)
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:394
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
Bounds::getNFX
int getNFX()
RET_STEPSIZE
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:138
QPS_NOTINITIALISED
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Types.hpp:74
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
Constraints::getActive
Indexlist * getActive()
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
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
SubjectTo::hasNoUpper
BooleanType hasNoUpper() const
Bounds
Manages working sets of bounds (= box constraints).
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/Bounds.hpp:50
SubjectTo::hasNoLower
BooleanType hasNoLower() const
ST_UPPER
#define ST_UPPER
Definition: qpOASES_wrapper.h:75
BEGIN_NAMESPACE_QPOASES
#define BEGIN_NAMESPACE_QPOASES
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/Types.hpp:117
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
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
QProblem::haveCholesky
BooleanType haveCholesky
Definition: QProblem.h:75
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
Flipper::init
returnValue init(int _nV=0, int _nC=0)
Definition: qpOASES-3.0beta/src/Flipper.cpp:116
QProblemB::loadQPvectorsFromFile
returnValue loadQPvectorsFromFile(const char *const g_file, const char *const lb_file, const char *const ub_file, real_t *const g_new, real_t *const lb_new, real_t *const ub_new) const
Definition: external_packages/qpOASES-3.0beta/src/QProblemB.cpp:1945
getGlobalMessageHandler
MessageHandling * getGlobalMessageHandler()
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/SRC/MessageHandling.cpp:515
QProblem::setA
returnValue setA(const real_t *const A_new)
RET_DIV_BY_ZERO
Definition: acado_types.hpp:837
QProblem::setUBA
returnValue setUBA(const real_t *const ubA_new)
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
Flipper::get
returnValue get(Bounds *const _bounds, double *const R, Constraints *const _constraints=0, double *const _Q=0, double *const _T=0) const
Definition: qpOASES-3.0beta/src/Flipper.cpp:147
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_FEWER_REGSTEPS_NWSR
Definition: external_packages/qpOASES-3.0beta/include/qpOASES/MessageHandling.hpp:179
RET_HOTSTART_STOPPED_INFEASIBILITY
Definition: examples/code_generation/mpc_mhe/getting_started_export/qpoases/INCLUDE/MessageHandling.hpp:126
SymmetricMatrix
Abstract base class for interfacing matrix-vector operations tailored to symmetric matrices...
Definition: qpOASES-3.0beta/include/qpOASES/Matrices.hpp:246
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