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