Scippy

SCIP

Solving Constraint Integer Programs

nlpi_filtersqp.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file nlpi_filtersqp.c
17  * @ingroup DEFPLUGINS_NLPI
18  * @brief filterSQP NLP interface
19  * @author Stefan Vigerske
20  *
21  * @todo scaling
22  */
23 
24 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
25 
26 #include <string.h>
27 #if defined(_WIN32)
28 #include <windows.h>
29 #else
30 #include <sys/time.h>
31 #endif
32 
33 /* fallback to non-thread version for windows, because pthread does not exist */
34 #if defined(_MSC_VER) && defined(SCIP_THREADSAFE)
35 #undef SCIP_THREADSAFE
36 #endif
37 
38 #ifdef SCIP_THREADSAFE
39 #include <pthread.h>
40 #endif
41 
42 #include "scip/nlpi_filtersqp.h"
43 #include "scip/nlpioracle.h"
44 #include "scip/scip_general.h"
45 #include "scip/scip_message.h"
46 #include "scip/scip_mem.h"
47 #include "scip/scip_numerics.h"
48 #include "scip/scip_nlpi.h"
49 #include "scip/scip_randnumgen.h"
50 #include "scip/scip_solve.h"
51 #include "scip/pub_misc.h"
52 
53 #define NLPI_NAME "filtersqp" /**< short concise name of solver */
54 #define NLPI_DESC "Sequential Quadratic Programming trust region solver by R. Fletcher and S. Leyffer" /**< description of solver */
55 #define NLPI_PRIORITY -1000 /**< priority of NLP solver */
56 
57 #define RANDSEED 26051979 /**< initial random seed */
58 #define MAXPERTURB 0.01 /**< maximal perturbation of bounds in starting point heuristic */
59 #define MAXNRUNS 3 /**< maximal number of FilterSQP calls per NLP solve (several calls if increasing workspace or decreasing eps) */
60 #define WORKSPACEGROWTHFACTOR 2 /**< factor by which to increase workspace */
61 #define MINEPS 1e-14 /**< minimal FilterSQP epsilon */
62 #define OPTTOLFACTOR 0.5 /**< factor to apply to optimality tolerance, because FilterSQP do scaling */
63 
64 /*
65  * Data structures
66  */
67 
68 typedef int fint;
69 typedef double real;
70 typedef long ftnlen;
71 
72 struct SCIP_Time
73 {
74  time_t sec; /**< seconds part of time since epoch */
75  long usec; /**< micro-seconds part of time */
76 };
77 typedef struct SCIP_Time SCIP_TIME;
78 
79 struct SCIP_NlpiData
80 {
81  SCIP_RANDNUMGEN* randnumgen; /**< random number generator, if we have to make up a starting point */
82  SCIP_TIME starttime; /**< time at start of solve */
83 };
84 
85 struct SCIP_NlpiProblem
86 {
87  SCIP* scip; /**< SCIP data structure */
88  SCIP_NLPIORACLE* oracle; /**< Oracle-helper to store and evaluate NLP */
89 
90  int varssize; /**< length of variables-related arrays, if allocated */
91  int conssize; /**< length of constraints-related arrays, if allocated */
92 
93  SCIP_Real* initguess; /**< initial values for primal variables, or NULL if not known, size varssize */
94  SCIP_Bool warmstart; /**< whether we could warmstart the next solve */
95 
96  SCIP_Real* primalvalues; /**< primal values of variables in solution, size varssize */
97  SCIP_Real* consdualvalues; /**< dual values of constraints in solution, size conssize */
98  SCIP_Real* varlbdualvalues; /**< dual values of variable lower bounds in solution, size varssize */
99  SCIP_Real* varubdualvalues; /**< dual values of variable upper bounds in solution, size varssize */
100 
101  SCIP_NLPSOLSTAT solstat; /**< solution status from last NLP solve */
102  SCIP_NLPTERMSTAT termstat; /**< termination status from last NLP solve */
103  SCIP_Real solvetime; /**< time spend for last NLP solve */
104  int niterations; /**< number of iterations for last NLP solve */
105 
106  fint istat[14]; /**< integer solution statistics from last FilterSQP call */
107  real rstat[7]; /**< real solution statistics from last FilterSQP call */
108  real fmin; /**< lower bound on objective value */
109  SCIP_Real maxtime; /**< time limit */
110 
111  /* cached FilterSQP data */
112  real* x; /**< variable values, size varssize */
113  real* c; /**< constraint value, size conssize */
114  real* lam; /**< duals, size varssize + conssize */
115  real* bl; /**< variable lower bounds and constraint lhs, size varssize + conssize */
116  real* bu; /**< variable upper bounds and constraint rhs, size varssize + conssize */
117  real* s; /**< scaling factors, size varssize + conssize */
118  char* cstype; /**< constraint linearity, size conssize */
119  real* a; /**< gradients values, size la[0]-1 */
120  fint* la; /**< gradients indices, size lasize */
121  int lasize; /**< length of la array */
122  fint* hessiannz; /**< nonzero information about Hessian, size hessiannzsize */
123  int hessiannzsize;/**< length of hessiannz array */
124  real* ws; /**< real workspace, size mxwk */
125  fint* lws; /**< integer workspace, size mxiwk */
126  fint mxwk; /**< size of real workspace */
127  fint mxiwk; /**< size of integer workspace */
128  real* evalbuffer; /**< buffer to cache evaluation results before passing it to FilterSQP, size evalbufsize */
129  int evalbufsize; /**< size of evaluation buffer */
130 };
131 
132 /*
133  * Local methods
134  */
135 
136 #ifdef FNAME_LCASE_DECOR
137 # define F77_FUNC(name,NAME) name ## _
138 #endif
139 #ifdef FNAME_UCASE_DECOR
140 # define F77_FUNC(name,NAME) NAME ## _
141 #endif
142 #ifdef FNAME_LCASE_NODECOR
143 # define F77_FUNC(name,NAME) name
144 #endif
145 #ifdef FNAME_UCASE_NODECOR
146 # define F77_FUNC(name,NAME) NAME
147 #endif
148 
149 /** FilterSQP main routine.
150  *
151  * Array a has length nnza, which is the number of nonzeros in the gradient of the objective and the Jacobian.
152  * The first entries of a is the objective gradient, next are the gradients of the constraints.
153  *
154  * Array la has length lamax, which is at least nnza+m+2.
155  * la contains the index information of a row-oriented sparse matrix storage. It stores the number of nonzeros, the column indices, and the row starts:
156  * la[0] must be set to nnza+1.
157  * la[1]..la[nnza] are indices of the variables corresponding to the entries in a (colidx).
158  * la[nnza+1]..la[nnza+1+m] contain the index where each row starts in a and la (rowstart).
159  */
160 void F77_FUNC(filtersqp,FILTERSQP)(
161  fint* n, /**< number of variables */
162  fint* m, /**< number of constraints (excluding simple bounds) */
163  fint* kmax, /**< maximum size of null-space (at most n) */
164  fint* maxa, /**< maximum number of nonzero entries allowed in Jacobian */
165  fint* maxf, /**< maximum size of the filter - typically 100 */
166  fint* mlp, /**< maximum level parameter for resolving degeneracy in bqpd - typically 100 */
167  fint* mxwk, /**< maximum size of real workspace ws */
168  fint* mxiwk, /**< maximum size of integer workspace lws */
169  fint* iprint, /**< print flag: 0 is quiet, higher is more */
170  fint* nout, /**< output channel - 6 for screen */
171  fint* ifail, /**< fail flag and warmstart indicator */
172  real* rho, /**< initial trust-region radius - default 10 */
173  real* x, /**< starting point and final solution (array of length n) */
174  real* c, /**< final values of general constraints (array of length m) */
175  real* f, /**< final objective value */
176  real* fmin, /**< lower bound on objective value (as termination criteria) */
177  real* bl, /**< lower bounds of variables and constraints (array of length n+m) */
178  real* bu, /**< upper bounds of variables and constraints (array of length n+m) */
179  real* s, /**< scaling factors (array of length n+m) */
180  real* a, /**< objective gradient (always dense) and Jacobian (sparse or dense) entries */
181  fint* la, /**< column indices and length of rows of entries in a (if sparse) or leading dimension of a (if dense) */
182  real* ws, /**< real workspace (array of length mxwk) */
183  fint* lws, /**< integer workspace (array of length mxiwk) */
184  real* lam, /**< Lagrangian multipliers of simple bounds and general constraints (array of length n+m) */
185  char* cstype, /**< indicator whether constraint is linear ('L') or nonlinear ('N') (array of length m) */
186  real* user, /**< real workspace passed through to user routines */
187  fint* iuser, /**< integer workspace passed through to user routines */
188  fint* maxiter, /**< iteration limit for SQP solver */
189  fint* istat, /**< integer space for solution statistics (array of length 14) */
190  real* rstat, /**< real space for solution statitics (array of length 7) */
191  ftnlen cstype_len /**< 1 ??? */
192  );
193 
194 void F77_FUNC(objfun,OBJFUN)(real *x, fint *n, real *f, real *user, fint *iuser,
195  fint *errflag);
196 
197 void F77_FUNC(confun,CONFUN)(real *x, fint *n, fint *m, real *c, real *a, fint *la,
198  real *user, fint *iuser, fint *errflag);
199 
200 /* TODO what are the arguments of this function and does it need to be implemented?
201  * it's not in the filterSQP manual, but its an undefined symbol in the filterSQP library
202 void F77_FUNC(objgrad,OBJGRAD)(fint *, fint *, fint *, real *, real *, fint *, fint
203  *, real *, fint *, fint *);
204 */
205 void F77_FUNC(objgrad,OBJGRAD)(void);
206 
207 void F77_FUNC(gradient,GRADIENT)(fint *n, fint *m, fint *mxa, real *x, real *a, fint *la,
208  fint *maxa, real *user, fint *iuser, fint *errflag);
209 
210 void F77_FUNC(hessian,HESSIAN)(real *x, fint *n, fint *m, fint *phase, real *lam,
211  real *ws, fint *lws, real *user, fint *iuser,
212  fint *l_hess, fint *li_hess, fint *errflag);
213 
214 /** common block for problemname */
215 /*lint -esym(754,char_l,pname,*::char_l,*::pname) */
216 extern struct
217 {
218  fint char_l;
219  char pname[10];
220 } F77_FUNC(cpname,CPNAME);
221 /*lint -esym(752,cpname_) */
222 
223 /** common block for Hessian storage set to 0, i.e. NO Hessian */
224 /*lint -esym(754,*::phr,*::phc) */
225 extern struct
226 {
228 } F77_FUNC(hessc,HESSC);
229 /*lint -esym(754,phr,phc) */
230 
231 /** common block for upper bound on filter */
232 extern struct
233 {
235 } F77_FUNC(ubdc,UBDC);
236 
237 /** common block for infinity & epsilon */
238 extern struct
239 {
241 } F77_FUNC(nlp_eps_inf,NLP_EPS_INF);
242 
243 /** common block for printing from QP solver */
244 /*lint -esym(754,*::n_bqpd_calls,*::n_bqpd_prfint) */
245 extern struct
246 {
248 } F77_FUNC(bqpd_count,BQPD_COUNT);
249 /*lint -esym(752,bqpd_count_) */
250 /*lint -esym(754,n_bqpd_calls,n_bqpd_prfint) */
251 
252 /** common for scaling: scale_mode = 0 (none), 1 (variables), 2 (vars+cons) */
253 /*lint -esym(754,*::phe) */
254 extern struct
255 {
257 } F77_FUNC(scalec,SCALEC);
258 /*lint -esym(754,phe) */
259 
260 #ifdef SCIP_THREADSAFE
261 static pthread_mutex_t filtersqpmutex = PTHREAD_MUTEX_INITIALIZER; /*lint !e708*/
262 #endif
263 
264 static
266 {
267  SCIP_TIME t;
268 #ifndef _WIN32
269  struct timeval tp; /*lint !e86*/
270 #endif
271 
272 #ifdef _WIN32
273  t.sec = time(NULL);
274  t.usec = 0;
275 
276 #else
277  (void) gettimeofday(&tp, NULL);
278  t.sec = tp.tv_sec;
279  t.usec = tp.tv_usec;
280 #endif
281 
282  return t;
283 }
284 
285 /* gives time since starttime (in problem) */
286 static
288  SCIP_NLPIDATA* nlpidata /**< NLPI data */
289  )
290 {
291  SCIP_TIME now;
292 
293  assert(nlpidata != NULL);
294 
295  now = gettime();
296 
297  /* now should be after startime */
298  assert(now.sec >= nlpidata->starttime.sec);
299  assert(now.sec > nlpidata->starttime.sec || now.usec >= nlpidata->starttime.usec);
300 
301  return (SCIP_Real)(now.sec - nlpidata->starttime.sec) + 1e-6 * (SCIP_Real)(now.usec - nlpidata->starttime.usec);
302 }
303 
304 static
306  SCIP_NLPIDATA* nlpidata, /**< NLPI data */
307  SCIP_NLPIPROBLEM* nlpiproblem /**< NLPI problem */
308  )
309 {
310  if( nlpiproblem->maxtime == SCIP_REAL_MAX ) /*lint !e777 */
311  return FALSE;
312 
313  return timeelapsed(nlpidata) >= nlpiproblem->maxtime;
314 }
315 
316 /** Objective function evaluation */ /*lint -e{715} */
317 void F77_FUNC(objfun,OBJFUN)(
318  real* x, /**< value of current variables (array of length n) */
319  fint* n, /**< number of variables */
320  real* f, /**< buffer to store value of objective function */
321  real* user, /**< user real workspace */
322  fint* iuser, /**< user integer workspace */
323  fint* errflag /**< set to 1 if arithmetic exception occurs, otherwise 0 */
324  )
325 { /*lint --e{715} */
326  SCIP_NLPIPROBLEM* problem;
327  real val;
328 
329  problem = (SCIP_NLPIPROBLEM*)(void*)iuser;
330  assert(problem != NULL);
331 
332  if( SCIPisSolveInterrupted(problem->scip) || timelimitreached((SCIP_NLPIDATA*)(void*)user, problem) )
333  {
334  SCIPdebugMsg(problem->scip, "interrupted or timelimit reached, issuing arithmetic exception in objfun\n");
335  *errflag = 1;
336  return;
337  }
338 
339  *errflag = 1;
340  if( SCIPnlpiOracleEvalObjectiveValue(problem->scip, problem->oracle, x, &val) == SCIP_OKAY && SCIPisFinite(val) )
341  {
342  *errflag = 0;
343  *f = val;
344  }
345  else
346  {
347  SCIPdebugMsg(problem->scip, "arithmetic exception in objfun\n");
348  }
349 }
350 
351 /** Constraint functions evaluation */ /*lint -e{715} */
352 void F77_FUNC(confun,CONFUN)(
353  real* x, /**< value of current variables (array of length n) */
354  fint* n, /**< number of variables */
355  fint* m, /**< number of constraints */
356  real* c, /**< buffer to store values of constraints (array of length m) */
357  real* a, /**< Jacobian matrix entries */
358  fint* la, /**< Jacobian index information */
359  real* user, /**< user real workspace */
360  fint* iuser, /**< user integer workspace */
361  fint* errflag /**< set to 1 if arithmetic exception occurs, otherwise 0 */
362  )
363 { /*lint --e{715} */
364  SCIP_NLPIPROBLEM* problem;
365  real val;
366  int j;
367 
368  problem = (SCIP_NLPIPROBLEM*)(void*)iuser;
369  assert(problem != NULL);
370 
371  *errflag = 0;
372  for( j = 0; j < *m; ++j )
373  {
374  if( SCIPnlpiOracleEvalConstraintValue(problem->scip, problem->oracle, j, x, &val) != SCIP_OKAY || !SCIPisFinite(val) )
375  {
376  *errflag = 1;
377  SCIPdebugMsg(problem->scip, "arithmetic exception in confun for constraint %d\n", j);
378  break;
379  }
380  c[j] = val;
381  }
382 }
383 
384 /** Objective gradient and Jacobian evaluation
385  *
386  * \note If an arithmetic exception occurred, then the gradients must not be modified.
387  */ /*lint -e{715} */
388 void
389 F77_FUNC(gradient,GRADIENT)(
390  fint* n, /**< number of variables */
391  fint* m, /**< number of constraints */
392  fint* mxa, /**< actual number of entries in a */
393  real* x, /**< value of current variables (array of length n) */
394  real* a, /**< Jacobian matrix entries */
395  fint* la, /**< Jacobian index information: column indices and pointers to start of each row */
396  fint* maxa, /**< maximal size of a */
397  real* user, /**< user real workspace */
398  fint* iuser, /**< user integer workspace */
399  fint* errflag /**< set to 1 if arithmetic exception occurs, otherwise 0 */
400  )
401 { /*lint --e{715} */
402  SCIP_NLPIPROBLEM* problem;
403  SCIP_Real dummy;
404 
405  problem = (SCIP_NLPIPROBLEM*)(void*)iuser;
406  assert(problem != NULL);
407  assert(problem->evalbuffer != NULL);
408  assert(problem->evalbufsize >= *maxa);
409 
410  *errflag = 1;
411 
412  if( SCIPnlpiOracleEvalObjectiveGradient(problem->scip, problem->oracle, x, TRUE, &dummy, problem->evalbuffer) == SCIP_OKAY )
413  {
414  if( SCIPnlpiOracleEvalJacobian(problem->scip, problem->oracle, x, TRUE, NULL, problem->evalbuffer+*n) == SCIP_OKAY )
415  {
416  BMScopyMemoryArray(a, problem->evalbuffer, *maxa);
417  *errflag = 0;
418  }
419  else
420  {
421  SCIPdebugMsg(problem->scip, "arithmetic exception in gradient for constraints\n");
422  }
423  }
424  else
425  {
426  SCIPdebugMsg(problem->scip, "arithmetic exception in gradient for objective\n");
427  }
428 }
429 
430 /* Objective gradient evaluation */
431 /*
432 void F77_FUNC(objgrad,OBJGRAD)(
433  fint*,
434  fint*,
435  fint*,
436  real*,
437  real*,
438  fint*,
439  fint*,
440  real*,
441  fint*,
442  fint*
443  )
444 */
445 void F77_FUNC(objgrad,OBJGRAD)(void)
446 {
447  SCIPerrorMessage("Objgrad not implemented. Should not be called.\n");
448 }
449 
450 /** Hessian of the Lagrangian evaluation
451  *
452  * phase = 1 : Hessian of the Lagrangian without objective Hessian
453  *
454  * phase = 2 : Hessian of the Lagrangian (including objective Hessian)
455  *
456  * \note If an arithmetic exception occurred, then the Hessian must not be modified.
457  */ /*lint -e{715} */
458 void
459 F77_FUNC(hessian,HESSIAN)(
460  real* x, /**< value of current variables (array of length n) */
461  fint* n, /**< number of variables */
462  fint* m, /**< number of constraints */
463  fint* phase, /**< indicates what kind of Hessian matrix is required */
464  real* lam, /**< Lagrangian multipliers (array of length n+m) */
465  real* ws, /**< real workspace for Hessian, passed to Wdotd */
466  fint* lws, /**< integer workspace for Hessian, passed to Wdotd */
467  real* user, /**< user real workspace */
468  fint* iuser, /**< user integer workspace */
469  fint* l_hess, /**< space of Hessian real storage ws. On entry: maximal space allowed, on exit: actual amount used */
470  fint* li_hess, /**< space of Hessian integer storage lws. On entry: maximal space allowed, on exit: actual amount used */
471  fint* errflag /**< set to 1 if arithmetic exception occurs, otherwise 0 */
472  )
473 { /*lint --e{715} */
474  SCIP_NLPIPROBLEM* problem;
475  SCIP_Real* lambda;
476  int nnz;
477  int i;
478 
479  problem = (SCIP_NLPIPROBLEM*)(void*)iuser;
480  assert(problem != NULL);
481  assert(problem->evalbuffer != NULL);
482 
483  nnz = problem->hessiannz[0]-1;
484  assert(problem->evalbufsize >= nnz);
485 
486  *errflag = 1;
487 
488  /* initialize lambda to -lam */
489  BMSallocMemoryArray(&lambda, *m);
490  for( i = 0; i < *m; ++i )
491  lambda[i] = -lam[*n+i];
492 
493  if( SCIPnlpiOracleEvalHessianLag(problem->scip, problem->oracle, x, TRUE, TRUE, (*phase == 1) ? 0.0 : 1.0, lambda, problem->evalbuffer) == SCIP_OKAY )
494  {
495  *l_hess = nnz;
496 
497  BMScopyMemoryArray(ws, problem->evalbuffer, nnz);
498 
499  *errflag = 0;
500 
501  /* copy the complete problem->hessiannz into lws */
502  for( i = 0; i < nnz + *n + 2; ++i )
503  lws[i] = problem->hessiannz[i];
504  *li_hess = nnz + *n + 2;
505  }
506  else
507  {
508  SCIPdebugMsg(problem->scip, "arithmetic exception in hessian\n");
509  }
510 
511  BMSfreeMemoryArray(&lambda);
512 }
513 
514 
515 
516 static
518  SCIP* scip, /**< SCIP data structure */
519  SCIP_NLPIORACLE* oracle, /**< NLPI oracle */
520  fint** la, /**< buffer to store pointer to sparsity structure */
521  int* lasize, /**< buffer to store length of *la array */
522  real** a /**< buffer to store pointer to value buffer */
523  )
524 {
525  const int* offset;
526  const int* col;
527  int nnz; /* number of nonzeros in Jacobian */
528  int nvars;
529  int ncons;
530  int i;
531  int c;
532 
533  assert(la != NULL);
534  assert(lasize != NULL);
535  assert(a != NULL);
536  assert(*la == NULL);
537  assert(*a == NULL);
538 
539  nvars = SCIPnlpiOracleGetNVars(oracle);
540  ncons = SCIPnlpiOracleGetNConstraints(oracle);
541 
542  /* get jacobian sparsity in oracle format: offset are rowstarts in col and col are column indices */
543  SCIP_CALL( SCIPnlpiOracleGetJacobianSparsity(scip, oracle, &offset, &col) );
544  nnz = offset[ncons];
545 
546  /* la stores both column indices (first) and rowstarts (at the end) of the objective gradient and Jacobian
547  * For the objective gradient, always n entries are taken, the Jacobian is stored sparse
548  * la(0) = n+nnz+1 position where rowstarts start in la
549  * la(j) column index of objective gradient or Jacobian row, rowwise
550  * la(la(0)) position of first nonzero element for objective gradient in a()
551  * la(la(0)+i) position of first nonzero element for constraint i gradient in a(), i=1..m
552  * la(la(0)+m+1) = n+nnz first unused position in a
553  * where n = nvars and m = ncons
554  */
555 
556  /* need space for la(0) and column indices and rowstarts (1+ncons+1 for objective, constraints, and end (offsets[ncons])) */
557  *lasize = 1 + (nvars+nnz) + 1+ncons + 1;
558  SCIP_CALL( SCIPallocBlockMemoryArray(scip, la, *lasize) );
559 
560  (*la)[0] = nvars+nnz+1;
561 
562  /* the objective gradient is stored in sparse form */
563  for( i = 0; i < nvars; ++i )
564  (*la)[1+i] = 1+i; /* shift by 1 for Fortran */
565  (*la)[(*la)[0]] = 1; /* objective entries start at the beginning in a, shift by 1 for Fortran */
566 
567  /* column indicies are as given by col */
568  for( i = 0; i < nnz; ++i )
569  (*la)[1+nvars+i] = col[i] + 1; /* shift by 1 for Fortran */
570 
571  /* rowstarts are as given by offset, plus extra for objective gradient */
572  for( c = 0; c <= ncons; ++c )
573  (*la)[(*la)[0]+1+c] = offset[c] + nvars + 1; /* shift by nvars for objective, shift by 1 for Fortran */
574 
575  SCIP_CALL( SCIPallocBlockMemoryArray(scip, a, nvars + nnz) );
576 
577 #ifdef SCIP_MORE_DEBUG /* enable for debugging Jacobian */
578  for( i = 0; i < 1 + (nvars+nnz) + 1+ncons + 1; ++i )
579  printf("la[%2d] = %2d\n", i, (*la)[i]);
580 #endif
581 
582  return SCIP_OKAY;
583 }
584 
585 static
587  SCIP* scip, /**< SCIP data structure */
588  SCIP_NLPIORACLE* oracle, /**< NLPI oracle */
589  fint** la, /**< buffer to store pointer to Hessian sparsity structure */
590  int* lasize /**< buffer to store length of *la array */
591  )
592 {
593  const int* offset;
594  const int* col;
595  int nnz; /* number of nonzeros in Jacobian */
596  int nvars;
597  int i;
598  int v;
599 
600  assert(la != NULL);
601  assert(lasize != NULL);
602  assert(*la == NULL);
603 
604  nvars = SCIPnlpiOracleGetNVars(oracle);
605 
606  /* get Hessian sparsity in oracle format: offset are rowstarts in col and col are column indices */
607  SCIP_CALL( SCIPnlpiOracleGetHessianLagSparsity(scip, oracle, &offset, &col) );
608  nnz = offset[nvars];
609 
610  /* la stores both column indices (first) and rowstarts (at the end) of the (sparse) Hessian
611  * la(0) = nnz+1 position where rowstarts start in la
612  * la(j) column index of Hessian row, rowwise
613  * la(la(0)+i) position of first nonzero element for row i, i = 0..n-1
614  * la(la(0)+n) = nnz first unused position in Hessian
615  * where n = nvars
616  */
617 
618  /* need space for la(0) and column indices and rowstarts
619  * 1 for la(0)
620  * nnz for column indices
621  * nvars for rowstarts
622  * 1 for first unused position
623  */
624  *lasize = 1 + nnz + nvars + 1;
625  SCIP_CALL( SCIPallocBlockMemoryArray(scip, la, *lasize) );
626 
627  (*la)[0] = nnz+1;
628 
629  /* column indicies are as given by col */
630  for( i = 0; i < nnz; ++i )
631  (*la)[1+i] = col[i] + 1; /* shift by 1 for Fortran */
632 
633  /* rowstarts are as given by offset */
634  for( v = 0; v <= nvars; ++v )
635  (*la)[(*la)[0]+v] = offset[v] + 1; /* shift by 1 for Fortran */
636 
637  F77_FUNC(hessc,HESSC).phl = 1;
638 
639 #ifdef SCIP_MORE_DEBUG /* enable for debugging Hessian */
640  for( i = 0; i < 1 + nnz + nvars + 1; ++i )
641  printf("lw[%2d] = %2d\n", i, (*la)[i]);
642 #endif
643 
644  return SCIP_OKAY;
645 }
646 
647 /** setup starting point for FilterSQP */
648 static
650  SCIP_NLPIDATA* data, /**< NLPI data */
651  SCIP_NLPIPROBLEM* problem, /**< NLPI problem */
652  real* x, /**< array to store initial values */
653  SCIP_Bool* success /**< whether we could setup a point in which functions could be evaluated */
654  )
655 {
656  int i;
657  int n;
658  SCIP_Real val;
659 
660  assert(data != NULL);
661  assert(problem != NULL);
662  assert(x != NULL);
663  assert(success != NULL);
664 
665  n = SCIPnlpiOracleGetNVars(problem->oracle);
666 
667  /* setup starting point */
668  if( problem->initguess != NULL )
669  {
670  for( i = 0; i < n; ++i )
671  x[i] = problem->initguess[i];
672  }
673  else
674  {
675  SCIP_Real lb;
676  SCIP_Real ub;
677 
678  SCIPdebugMsg(problem->scip, "FilterSQP started without initial primal values; make up something by projecting 0 onto variable bounds and perturb\n");
679 
680  if( data->randnumgen == NULL )
681  {
682  SCIP_CALL( SCIPcreateRandom(problem->scip, &data->randnumgen, RANDSEED, TRUE) );
683  }
684 
685  for( i = 0; i < n; ++i )
686  {
687  lb = SCIPnlpiOracleGetVarLbs(problem->oracle)[i];
688  ub = SCIPnlpiOracleGetVarUbs(problem->oracle)[i];
689  if( lb > 0.0 )
690  x[i] = SCIPrandomGetReal(data->randnumgen, lb, lb + MAXPERTURB*MIN(1.0, ub-lb));
691  else if( ub < 0.0 )
692  x[i] = SCIPrandomGetReal(data->randnumgen, ub - MAXPERTURB*MIN(1.0, ub-lb), ub);
693  else
694  x[i] = SCIPrandomGetReal(data->randnumgen, MAX(lb, -MAXPERTURB*MIN(1.0, ub-lb)), MIN(ub, MAXPERTURB*MIN(1.0, ub-lb)));
695  }
696  }
697 
698  /* check whether objective and constraints can be evaluated and differentiated once in starting point
699  * NOTE: this does not check whether hessian can be computed!
700  */
701  *success = SCIPnlpiOracleEvalObjectiveValue(problem->scip, problem->oracle, x, &val) == SCIP_OKAY && SCIPisFinite(val);
702  *success &= SCIPnlpiOracleEvalObjectiveGradient(problem->scip, problem->oracle, x, FALSE, &val, problem->evalbuffer) == SCIP_OKAY; /*lint !e514*/
703  i = 0;
704  for( ; *success && i < SCIPnlpiOracleGetNConstraints(problem->oracle); ++i )
705  *success = SCIPnlpiOracleEvalConstraintValue(problem->scip, problem->oracle, i, x, &val) == SCIP_OKAY && SCIPisFinite(val);
706  *success &= SCIPnlpiOracleEvalJacobian(problem->scip, problem->oracle, x, FALSE, NULL, problem->evalbuffer) == SCIP_OKAY; /*lint !e514*/
707 
708  if( !*success )
709  {
710  SCIPdebugMsg(problem->scip, "could not evaluate or constraint %d in %s starting point or Jacobian not available\n", i-1, problem->initguess != NULL ? "provided" : "made up");
711 
712  if( problem->initguess != NULL )
713  {
714  /* forget given starting point and try to make up our own */
715  SCIPfreeBlockMemoryArray(problem->scip, &problem->initguess, problem->varssize);
716  SCIP_CALL( setupStart(data, problem, x, success) );
717  }
718  }
719 
720  return SCIP_OKAY;
721 }
722 
723 /** sets the solstat and termstat to unknown and other, resp. */
724 static
726  SCIP_NLPIPROBLEM* problem /**< data structure of problem */
727  )
728 {
729  assert(problem != NULL);
730 
731  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
732  problem->termstat = SCIP_NLPTERMSTAT_OTHER;
733 }
734 
735 /** store NLP solve parameters in nlpiproblem */
736 static
738  SCIP* scip, /**< SCIP data structure */
739  SCIP_NLPIPROBLEM* nlpiproblem, /**< NLP */
740  const SCIP_NLPPARAM param /**< solve parameters */
741  )
742 {
743  assert(scip != NULL);
744  assert(nlpiproblem != NULL);
745 
746  if( param.fastfail )
747  {
748  SCIPdebugMsg(scip, "fast fail parameter not supported by FilterSQP interface yet. Ignored.\n");
749  }
750 
751  nlpiproblem->fmin = param.lobjlimit;
752 
753  nlpiproblem->maxtime = param.timelimit;
754 
755  return SCIP_OKAY;
756 }
757 
758 /** processes results from FilterSQP call */
759 static
761  SCIP_NLPIDATA* nlpidata, /**< NLPI data */
762  SCIP_NLPIPROBLEM* problem, /**< NLPI problem */
763  fint ifail, /**< fail flag from FilterSQP call */
764  SCIP_Real feastol, /**< feasibility tolerance */
765  SCIP_Real opttol, /**< optimality tolerance */
766  real* x, /**< primal solution values from FilterSQP call, or NULL if stopped before filtersqp got called */
767  real* lam /**< dual solution values from FilterSQP call, or NULL if stopped before filtersqp got called */
768  )
769 {
770  int i;
771  int nvars;
772  int ncons;
773 
774  assert(problem != NULL);
775  assert(ifail >= 0);
776  assert((x != NULL) == (lam != NULL));
777 
778  problem->solvetime = timeelapsed(nlpidata);
779 
780  nvars = SCIPnlpiOracleGetNVars(problem->oracle);
781  ncons = SCIPnlpiOracleGetNConstraints(problem->oracle);
782 
783  if( ifail <= 8 && x != NULL )
784  {
785  /* FilterSQP terminated somewhat normally -> store away solution */
786 
787  /* make sure we have memory for solution */
788  if( problem->primalvalues == NULL )
789  {
790  assert(problem->varssize >= nvars); /* ensured in nlpiAddVariables */
791  assert(problem->conssize >= ncons); /* ensured in nlpiAddConstraints */
792  assert(problem->consdualvalues == NULL); /* if primalvalues == NULL, then also consdualvalues should be NULL */
793  assert(problem->varlbdualvalues == NULL); /* if primalvalues == NULL, then also varlbdualvalues should be NULL */
794  assert(problem->varubdualvalues == NULL); /* if primalvalues == NULL, then also varubdualvalues should be NULL */
795 
796  SCIP_CALL( SCIPallocBlockMemoryArray(problem->scip, &problem->primalvalues, problem->varssize) );
797  SCIP_CALL( SCIPallocBlockMemoryArray(problem->scip, &problem->consdualvalues, problem->conssize) );
798  SCIP_CALL( SCIPallocBlockMemoryArray(problem->scip, &problem->varlbdualvalues, problem->varssize) );
799  SCIP_CALL( SCIPallocBlockMemoryArray(problem->scip, &problem->varubdualvalues, problem->varssize) );
800  }
801  else
802  {
803  assert(problem->consdualvalues != NULL);
804  assert(problem->varlbdualvalues != NULL);
805  assert(problem->varubdualvalues != NULL);
806  }
807 
808  for( i = 0; i < nvars; ++i )
809  {
810  problem->primalvalues[i] = x[i];
811 
812  /* if dual from FilterSQP is positive, then it belongs to the lower bound, otherwise to the upper bound */
813  problem->varlbdualvalues[i] = MAX(0.0, lam[i]);
814  problem->varubdualvalues[i] = MAX(0.0, -lam[i]);
815  }
816 
817  /* duals from FilterSQP are negated */
818  for( i = 0; i < ncons; ++i )
819  problem->consdualvalues[i] = -lam[nvars + i];
820  }
821 
822  /* translate ifail to solution and termination status and decide whether we could warmstart next */
823  problem->warmstart = FALSE;
824  switch( ifail )
825  {
826  case 0: /* successful run, solution found */
827  assert(problem->rstat[4] <= feastol); /* should be feasible */
828  problem->solstat = (problem->rstat[0] <= opttol ? SCIP_NLPSOLSTAT_LOCOPT : SCIP_NLPSOLSTAT_FEASIBLE);
829  problem->termstat = SCIP_NLPTERMSTAT_OKAY;
830  problem->warmstart = TRUE;
831  break;
832  case 1: /* unbounded, feasible point with f(x) <= fmin */
833  assert(problem->rstat[4] <= feastol); /* should be feasible */
835  if( problem->fmin == SCIP_REAL_MIN ) /*lint !e777*/
836  problem->termstat = SCIP_NLPTERMSTAT_OKAY; /* fmin was not set */
837  else
839  break;
840  case 2: /* linear constraints are inconsistent */
842  problem->termstat = SCIP_NLPTERMSTAT_OKAY;
843  break;
844  case 3: /* (locally) nonlinear infeasible, minimal-infeasible solution found */
845  /* problem->solstat = (problem->rstat[0] <= opttol ? SCIP_NLPSOLSTAT_LOCINFEASIBLE : SCIP_NLPSOLSTAT_UNKNOWN); */
846  problem->solstat = SCIP_NLPSOLSTAT_LOCINFEASIBLE; /* TODO FilterSQP does not set rstat[0] in this case, assuming local infeasibility is valid */
847  problem->termstat = SCIP_NLPTERMSTAT_OKAY;
848  problem->warmstart = TRUE;
849  break;
850  case 4: /* terminate at point with h(x) <= eps (constraint violation below epsilon) but QP infeasible */
851  assert(problem->rstat[4] <= feastol); /* should be feasible */
854  problem->warmstart = TRUE;
855  break;
856  case 5: /* termination with rho < eps (trust region radius below epsilon) */
857  if( problem->rstat[4] <= feastol )
859  else
860  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
862  problem->warmstart = TRUE;
863  break;
864  case 6: /* termination with iter > max_iter */
865  if( problem->rstat[4] <= feastol )
867  else
868  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
870  problem->warmstart = TRUE;
871  break;
872  case 7: /* crash in user routine (IEEE error) could not be resolved, or timelimit reached, or interrupted */
873  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
874  if( problem->solvetime >= problem->maxtime )
875  {
877  problem->warmstart = TRUE;
878  }
879  else if( SCIPisSolveInterrupted(problem->scip) )
881  else
883  break;
884  case 8: /* unexpect ifail from QP solver */
885  if( problem->rstat[4] <= feastol )
887  else
888  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
889  problem->termstat = SCIP_NLPTERMSTAT_OTHER;
890  break;
891  case 9: /* not enough REAL workspace */
892  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
894  break;
895  case 10: /* not enough INTEGER workspace */
896  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
898  break;
899  default:
900  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
901  problem->termstat = SCIP_NLPTERMSTAT_OTHER;
902  break;
903  }
904 
905  return SCIP_OKAY;
906 }
907 
908 
909 /*
910  * Callback methods of NLP solver interface
911  */
912 
913 /** copy method of NLP interface (called when SCIP copies plugins) */
914 static
915 SCIP_DECL_NLPICOPY(nlpiCopyFilterSQP)
916 {
918 
919  return SCIP_OKAY; /*lint !e527*/
920 } /*lint !e715*/
921 
922 /** destructor of NLP interface to free nlpi data */
923 static
924 SCIP_DECL_NLPIFREE(nlpiFreeFilterSQP)
925 {
926  assert(nlpi != NULL);
927  assert(nlpidata != NULL);
928  assert(*nlpidata != NULL);
929 
930  if( (*nlpidata)->randnumgen != NULL )
931  {
932  SCIPfreeRandom(scip, &(*nlpidata)->randnumgen);
933  }
934 
935  SCIPfreeBlockMemory(scip, nlpidata);
936  assert(*nlpidata == NULL);
937 
938  return SCIP_OKAY;
939 }
940 
941 /** creates a problem instance */
942 static
943 SCIP_DECL_NLPICREATEPROBLEM(nlpiCreateProblemFilterSQP)
944 {
945  assert(nlpi != NULL);
946  assert(problem != NULL);
947 
949  (*problem)->scip = scip;
950 
951  SCIP_CALL( SCIPnlpiOracleCreate(scip, &(*problem)->oracle) );
952  SCIP_CALL( SCIPnlpiOracleSetProblemName(scip, (*problem)->oracle, name) );
953 
954  (*problem)->fmin = SCIP_REAL_MIN;
955  (*problem)->maxtime = SCIP_REAL_MAX;
956 
957  invalidateSolution(*problem);
958 
959  return SCIP_OKAY; /*lint !e527*/
960 } /*lint !e715*/
961 
962 /** free a problem instance */
963 static
964 SCIP_DECL_NLPIFREEPROBLEM(nlpiFreeProblemFilterSQP)
965 {
966  assert(nlpi != NULL);
967  assert(problem != NULL);
968  assert(*problem != NULL);
969 
970  if( (*problem)->oracle != NULL )
971  {
972  SCIP_CALL( SCIPnlpiOracleFree(scip, &(*problem)->oracle) );
973  }
974 
975  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->initguess, (*problem)->varssize);
976  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->primalvalues, (*problem)->varssize);
977  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->consdualvalues, (*problem)->conssize);
978  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->varlbdualvalues, (*problem)->varssize);
979  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->varubdualvalues, (*problem)->varssize);
980  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->cstype, (*problem)->conssize);
981  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->s, (*problem)->varssize + (*problem)->conssize); /*lint !e776 */
982  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->bu, (*problem)->varssize + (*problem)->conssize); /*lint !e776 */
983  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->bl, (*problem)->varssize + (*problem)->conssize); /*lint !e776 */
984  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->x, (*problem)->varssize);
985  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->c, (*problem)->conssize);
986  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->lam, (*problem)->varssize + (*problem)->conssize); /*lint !e776 */
987  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->a, (*problem)->la != NULL ? (*problem)->la[0]-1 : 0);
988  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->la, (*problem)->lasize);
989  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->hessiannz, (*problem)->hessiannzsize);
990  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->lws, (*problem)->mxiwk);
991  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->ws, (*problem)->mxwk);
992  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->evalbuffer, (*problem)->evalbufsize);
993 
994  SCIPfreeBlockMemory(scip, problem);
995  assert(*problem == NULL);
996 
997  return SCIP_OKAY;
998 } /*lint !e715*/
999 
1000 /** add variables */
1001 static
1002 SCIP_DECL_NLPIADDVARS(nlpiAddVarsFilterSQP)
1003 {
1004  int oldnvars;
1005 
1006  assert(nlpi != NULL);
1007  assert(problem != NULL);
1008  assert(problem->oracle != NULL);
1009 
1010  oldnvars = SCIPnlpiOracleGetNVars(problem->oracle);
1011 
1012  SCIP_CALL( SCIPnlpiOracleAddVars(scip, problem->oracle, nvars, lbs, ubs, varnames) );
1013 
1014  SCIPfreeBlockMemoryArrayNull(scip, &problem->initguess, problem->varssize);
1015  invalidateSolution(problem);
1016  problem->warmstart = FALSE;
1017 
1018  /* increase variables-related arrays in problem, if necessary */
1019  if( problem->varssize < SCIPnlpiOracleGetNVars(problem->oracle) )
1020  {
1021  int newsize = SCIPcalcMemGrowSize(scip, SCIPnlpiOracleGetNVars(problem->oracle));
1022 
1023  if( problem->primalvalues != NULL )
1024  {
1025  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->primalvalues, problem->varssize, newsize) );
1026  }
1027 
1028  if( problem->varlbdualvalues != NULL )
1029  {
1030  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->varlbdualvalues, problem->varssize, newsize) );
1031  }
1032 
1033  if( problem->varubdualvalues != NULL )
1034  {
1035  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->varubdualvalues, problem->varssize, newsize) );
1036  }
1037 
1038  if( problem->x != NULL )
1039  {
1040  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->x, problem->varssize, newsize) );
1041  }
1042 
1043  if( problem->lam != NULL )
1044  {
1045  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->lam, problem->varssize + problem->conssize, newsize + problem->conssize) ); /*lint !e776*/
1046  }
1047 
1048  if( problem->bl != NULL )
1049  {
1050  assert(problem->bu != NULL);
1051  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->bl, problem->varssize + problem->conssize, newsize + problem->conssize) ); /*lint !e776*/
1052  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->bu, problem->varssize + problem->conssize, newsize + problem->conssize) ); /*lint !e776*/
1053  }
1054 
1055  if( problem->s != NULL )
1056  {
1057  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->s, problem->varssize + problem->conssize, newsize + problem->conssize) ); /*lint !e776*/
1058  }
1059 
1060  problem->varssize = newsize;
1061  }
1062 
1063  /* update variable bounds in FilterSQP data */
1064  if( problem->bl != NULL )
1065  {
1066  int i;
1067  int nconss;
1068 
1069  nconss = SCIPnlpiOracleGetNConstraints(problem->oracle);
1070 
1071  /* bl and bu have first variable bounds and then constraint sides
1072  * copy the constraint sides to their new position before putting in the new variable bounds
1073  */
1074  for( i = nconss-1; i >= 0; --i )
1075  {
1076  problem->bl[oldnvars+nvars+i] = problem->bl[oldnvars+i];
1077  problem->bu[oldnvars+nvars+i] = problem->bu[oldnvars+i];
1078  }
1079 
1080  /* set bounds of new variable */
1081  for( i = 0; i < nvars; ++i )
1082  {
1083  problem->bl[oldnvars+i] = lbs[i];
1084  problem->bu[oldnvars+i] = ubs[i];
1085  }
1086  }
1087 
1088  /* gradients information is out of date now (objective gradient is stored in dense form) */
1089  SCIPfreeBlockMemoryArrayNull(scip, &problem->a, problem->la != NULL ? problem->la[0]-1 : 0);
1090  SCIPfreeBlockMemoryArrayNull(scip, &problem->la, problem->lasize);
1091 
1092  /* Hessian information is out of date now (no new entries in Hessian, but also empty cols shows up in sparsity info) */
1093  SCIPfreeBlockMemoryArrayNull(scip, &problem->hessiannz, problem->hessiannzsize);
1094 
1095  return SCIP_OKAY;
1096 } /*lint !e715*/
1097 
1098 
1099 /** add constraints */
1100 static
1101 SCIP_DECL_NLPIADDCONSTRAINTS(nlpiAddConstraintsFilterSQP)
1102 {
1103  int oldnconss;
1104 
1105  assert(nlpi != NULL);
1106  assert(problem != NULL);
1107  assert(problem->oracle != NULL);
1108 
1109  oldnconss = SCIPnlpiOracleGetNConstraints(problem->oracle);
1110 
1111  SCIP_CALL( SCIPnlpiOracleAddConstraints(scip, problem->oracle, nconss, lhss, rhss, nlininds, lininds, linvals, exprs, names) );
1112 
1113  invalidateSolution(problem);
1114  problem->warmstart = FALSE;
1115 
1116  /* increase constraints-related arrays in problem, if necessary */
1117  if( SCIPnlpiOracleGetNConstraints(problem->oracle) > problem->conssize )
1118  {
1119  int newsize = SCIPcalcMemGrowSize(scip, SCIPnlpiOracleGetNConstraints(problem->oracle));
1120 
1121  if( problem->consdualvalues != NULL )
1122  {
1123  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->consdualvalues, problem->conssize, newsize) );
1124  }
1125 
1126  if( problem->c != NULL )
1127  {
1128  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->c, problem->conssize, newsize) );
1129  }
1130 
1131  if( problem->lam != NULL )
1132  {
1133  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->lam, problem->varssize + problem->conssize, problem->varssize + newsize) ); /*lint !e776*/
1134  }
1135 
1136  if( problem->bl != NULL )
1137  {
1138  assert(problem->bu != NULL);
1139  assert(problem->cstype != NULL);
1140  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->bl, problem->varssize + problem->conssize, problem->varssize + newsize) ); /*lint !e776*/
1141  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->bu, problem->varssize + problem->conssize, problem->varssize + newsize) ); /*lint !e776*/
1142  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->cstype, problem->conssize, newsize) );
1143  }
1144 
1145  if( problem->s != NULL )
1146  {
1147  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->s, problem->varssize + problem->conssize, problem->varssize + newsize) ); /*lint !e776*/
1148  }
1149 
1150  problem->conssize = newsize;
1151  }
1152 
1153  /* update constraint sides and type in FilterSQP data */
1154  if( problem->bl != NULL )
1155  {
1156  int i;
1157  int nvars;
1158 
1159  nvars = SCIPnlpiOracleGetNVars(problem->oracle);
1160 
1161  for( i = 0; i < nconss; ++i )
1162  {
1163  problem->bl[nvars+oldnconss+i] = lhss[i];
1164  problem->bu[nvars+oldnconss+i] = rhss[i];
1165  problem->cstype[oldnconss+i] = SCIPnlpiOracleIsConstraintNonlinear(problem->oracle, oldnconss+i) ? 'N' : 'L';
1166  }
1167  }
1168 
1169  /* gradients information is out of date now */
1170  SCIPfreeBlockMemoryArrayNull(scip, &problem->a, problem->la != NULL ? problem->la[0]-1 : 0);
1171  SCIPfreeBlockMemoryArrayNull(scip, &problem->la, problem->lasize);
1172 
1173  /* Hessian information is out of date now */
1174  SCIPfreeBlockMemoryArrayNull(scip, &problem->hessiannz, problem->hessiannzsize);
1175 
1176  return SCIP_OKAY;
1177 } /*lint !e715*/
1178 
1179 /** sets or overwrites objective, a minimization problem is expected */
1180 static
1181 SCIP_DECL_NLPISETOBJECTIVE( nlpiSetObjectiveFilterSQP )
1182 {
1183  assert(nlpi != NULL);
1184  assert(problem != NULL);
1185  assert(problem->oracle != NULL);
1186 
1187  SCIP_CALL( SCIPnlpiOracleSetObjective(scip, problem->oracle,
1188  constant, nlins, lininds, linvals, expr) );
1189 
1190  invalidateSolution(problem);
1191 
1192  /* gradients info (la,a) should still be ok, as objective gradient is stored in dense form */
1193 
1194  /* Hessian information is out of date now */
1195  SCIPfreeBlockMemoryArrayNull(scip, &problem->hessiannz, problem->hessiannzsize);
1196 
1197  return SCIP_OKAY;
1198 } /*lint !e715*/
1199 
1200 /** change variable bounds */
1201 static
1202 SCIP_DECL_NLPICHGVARBOUNDS(nlpiChgVarBoundsFilterSQP)
1203 {
1204  assert(nlpi != NULL);
1205  assert(problem != NULL);
1206  assert(problem->oracle != NULL);
1207 
1208  SCIP_CALL( SCIPnlpiOracleChgVarBounds(scip, problem->oracle, nvars, indices, lbs, ubs) );
1209 
1210  invalidateSolution(problem);
1211 
1212  /* update bounds in FilterSQP data */
1213  if( problem->bl != NULL )
1214  {
1215  int i;
1216 
1217  for( i = 0; i < nvars; ++i )
1218  {
1219  problem->bl[indices[i]] = lbs[i];
1220  problem->bu[indices[i]] = ubs[i];
1221  }
1222  }
1223 
1224  return SCIP_OKAY;
1225 } /*lint !e715*/
1226 
1227 /** change constraint bounds */
1228 static
1229 SCIP_DECL_NLPICHGCONSSIDES(nlpiChgConsSidesFilterSQP)
1230 {
1231  assert(nlpi != NULL);
1232  assert(problem != NULL);
1233  assert(problem->oracle != NULL);
1234 
1235  SCIP_CALL( SCIPnlpiOracleChgConsSides(scip, problem->oracle, nconss, indices, lhss, rhss) );
1236 
1237  invalidateSolution(problem);
1238 
1239  /* update constraint sense in FilterSQP data */
1240  if( problem->bl != NULL )
1241  {
1242  int i;
1243  int nvars;
1244 
1245  nvars = SCIPnlpiOracleGetNVars(problem->oracle);
1246 
1247  for( i = 0; i < nconss; ++i )
1248  {
1249  problem->bl[nvars+indices[i]] = lhss[i];
1250  problem->bu[nvars+indices[i]] = rhss[i];
1251  }
1252  }
1253 
1254  return SCIP_OKAY;
1255 } /*lint !e715*/
1256 
1257 /** delete a set of variables */
1258 static
1259 SCIP_DECL_NLPIDELVARSET(nlpiDelVarSetFilterSQP)
1260 {
1261  assert(nlpi != NULL);
1262  assert(problem != NULL);
1263  assert(problem->oracle != NULL);
1264 
1265  SCIP_CALL( SCIPnlpiOracleDelVarSet(scip, problem->oracle, dstats) );
1266 
1267  /* @TODO keep initguess and bl, bu for remaining variables? */
1268 
1269  SCIPfreeBlockMemoryArrayNull(scip, &problem->initguess, problem->varssize);
1270  invalidateSolution(problem);
1271  problem->warmstart = FALSE;
1272 
1273  SCIPfreeBlockMemoryArrayNull(scip, &problem->bl, problem->varssize + problem->conssize); /*lint !e776 */
1274  SCIPfreeBlockMemoryArrayNull(scip, &problem->bu, problem->varssize + problem->conssize); /*lint !e776 */
1275  SCIPfreeBlockMemoryArrayNull(scip, &problem->cstype, problem->conssize); /* because we assume that cstype is allocated iff bl is allocated */
1276 
1277  /* gradients information is out of date now (objective gradient is stored in dense form) */
1278  SCIPfreeBlockMemoryArrayNull(scip, &problem->a, problem->la != NULL ? problem->la[0]-1 : 0);
1279  SCIPfreeBlockMemoryArrayNull(scip, &problem->la, problem->lasize);
1280 
1281  /* Hessian information is out of date now */
1282  SCIPfreeBlockMemoryArrayNull(scip, &problem->hessiannz, problem->hessiannzsize);
1283 
1284  return SCIP_OKAY;
1285 } /*lint !e715*/
1286 
1287 /** delete a set of constraints */
1288 static
1289 SCIP_DECL_NLPIDELCONSSET(nlpiDelConstraintSetFilterSQP)
1290 {
1291  assert(nlpi != NULL);
1292  assert(problem != NULL);
1293  assert(problem->oracle != NULL);
1294 
1295  SCIP_CALL( SCIPnlpiOracleDelConsSet(scip, problem->oracle, dstats) );
1296 
1297  invalidateSolution(problem);
1298  problem->warmstart = FALSE;
1299 
1300  SCIPfreeBlockMemoryArrayNull(scip, &problem->bl, problem->varssize + problem->conssize); /*lint !e776 */
1301  SCIPfreeBlockMemoryArrayNull(scip, &problem->bu, problem->varssize + problem->conssize); /*lint !e776 */
1302  SCIPfreeBlockMemoryArrayNull(scip, &problem->cstype, problem->conssize); /* because we assume that cstype is allocated iff bl is allocated */
1303 
1304  /* gradients information is out of date now */
1305  SCIPfreeBlockMemoryArrayNull(scip, &problem->a, problem->la != NULL ? problem->la[0]-1 : 0);
1306  SCIPfreeBlockMemoryArrayNull(scip, &problem->la, problem->lasize);
1307 
1308  /* Hessian information is out of date now */
1309  SCIPfreeBlockMemoryArrayNull(scip, &problem->hessiannz, problem->hessiannzsize);
1310 
1311  return SCIP_OKAY;
1312 } /*lint !e715*/
1313 
1314 /** changes (or adds) linear coefficients in a constraint or objective */
1315 static
1316 SCIP_DECL_NLPICHGLINEARCOEFS(nlpiChgLinearCoefsFilterSQP)
1317 {
1318  assert(nlpi != NULL);
1319  assert(problem != NULL);
1320  assert(problem->oracle != NULL);
1321 
1322  SCIP_CALL( SCIPnlpiOracleChgLinearCoefs(scip, problem->oracle, idx, nvals, varidxs, vals) );
1323 
1324  invalidateSolution(problem);
1325 
1326  /* gradients information (la,a) may have changed if elements were added or removed
1327  * (we only care that sparsity doesn't change, not about actual values in a)
1328  * TODO free only if coefficients were added or removed (SCIPnlpiOracleChgLinearCoefs() could give feedback)
1329  */
1330  SCIPfreeBlockMemoryArrayNull(scip, &problem->a, problem->la != NULL ? problem->la[0]-1 : 0);
1331  SCIPfreeBlockMemoryArrayNull(scip, &problem->la, problem->lasize);
1332 
1333  /* Hessian information should still be ok */
1334 
1335  return SCIP_OKAY;
1336 } /*lint !e715*/
1337 
1338 /** replaces the expression of a constraint or objective */
1339 static
1340 SCIP_DECL_NLPICHGEXPR(nlpiChgExprFilterSQP)
1341 {
1342  assert(nlpi != NULL);
1343  assert(problem != NULL);
1344  assert(problem->oracle != NULL);
1345 
1346  SCIP_CALL( SCIPnlpiOracleChgExpr(scip, problem->oracle, idxcons, expr) );
1347 
1348  invalidateSolution(problem);
1349 
1350  /* update constraint linearity in FilterSQP data, as we might have changed from linear to nonlinear now */
1351  if( problem->cstype != NULL && idxcons >= 0 )
1352  problem->cstype[idxcons] = expr != NULL ? 'N' : 'L';
1353 
1354  /* gradients information (la,a) may have changed */
1355  SCIPfreeBlockMemoryArrayNull(scip, &problem->a, problem->la != NULL ? problem->la[0]-1 : 0);
1356  SCIPfreeBlockMemoryArrayNull(scip, &problem->la, problem->lasize);
1357 
1358  /* Hessian information may have changed */
1359  SCIPfreeBlockMemoryArrayNull(scip, &problem->hessiannz, problem->hessiannzsize);
1360 
1361  return SCIP_OKAY;
1362 } /*lint !e715*/
1363 
1364 /** change the constant offset in the objective */
1365 static
1366 SCIP_DECL_NLPICHGOBJCONSTANT(nlpiChgObjConstantFilterSQP)
1367 {
1368  assert(nlpi != NULL);
1369  assert(problem != NULL);
1370  assert(problem->oracle != NULL);
1371 
1372  SCIP_CALL( SCIPnlpiOracleChgObjConstant(scip, problem->oracle, objconstant) );
1373 
1374  invalidateSolution(problem);
1375 
1376  return SCIP_OKAY;
1377 } /*lint !e715*/
1378 
1379 /** sets initial guess for primal variables */
1380 static
1381 SCIP_DECL_NLPISETINITIALGUESS(nlpiSetInitialGuessFilterSQP)
1382 {
1383  assert(nlpi != NULL);
1384  assert(problem != NULL);
1385  assert(problem->oracle != NULL);
1386 
1387  if( primalvalues != NULL )
1388  {
1389  if( problem->initguess == NULL )
1390  {
1391  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->initguess, problem->varssize) );
1392  }
1393  assert(SCIPnlpiOracleGetNVars(problem->oracle) <= problem->varssize);
1394  BMScopyMemoryArray(problem->initguess, primalvalues, SCIPnlpiOracleGetNVars(problem->oracle));
1395  }
1396  else
1397  {
1398  SCIPfreeBlockMemoryArrayNull(scip, &problem->initguess, problem->varssize);
1399  }
1400 
1401  return SCIP_OKAY;
1402 } /*lint !e715*/
1403 
1404 /** tries to solve NLP */
1405 static
1406 SCIP_DECL_NLPISOLVE(nlpiSolveFilterSQP)
1407 {
1408  SCIP_NLPIDATA* data;
1409  SCIP_Bool success;
1410  fint n;
1411  fint m;
1412  fint kmax;
1413  fint maxa;
1414  fint maxf;
1415  fint mlp;
1416  fint lh1;
1417  fint nout;
1418  fint ifail;
1419  fint maxiter;
1420  fint iprint;
1421  real rho;
1422  real f;
1423  real* user;
1424  fint* iuser;
1425  ftnlen cstype_len = 1;
1426  fint minmxwk;
1427  fint minmxiwk;
1428  int nruns;
1429  int i;
1430 
1431  SCIPdebugMsg(scip, "solve with parameters " SCIP_NLPPARAM_PRINT(param));
1432 
1433  data = SCIPnlpiGetData(nlpi);
1434  assert(data != NULL);
1435 
1436  SCIP_CALL( SCIPnlpiOracleResetEvalTime(scip, problem->oracle) );
1437 
1438  if( param.timelimit == 0.0 )
1439  {
1440  /* there is nothing we can do if we are not given any time */
1441  problem->niterations = 0;
1442  problem->solvetime = 0.0;
1443  problem->termstat = SCIP_NLPTERMSTAT_TIMELIMIT;
1444  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
1445 
1446  return SCIP_OKAY;
1447  }
1448 
1449  /* start measuring time */
1450  data->starttime = gettime();
1451 
1452  SCIP_CALL( handleNlpParam(scip, problem, param) );
1453 
1454  iprint = param.verblevel;
1455 
1456  /* if warmstart parameter is disabled, then we will not warmstart */
1457  if( !param.warmstart )
1458  problem->warmstart = FALSE;
1459 
1460  n = SCIPnlpiOracleGetNVars(problem->oracle);
1461  m = SCIPnlpiOracleGetNConstraints(problem->oracle);
1462  kmax = n; /* maximal nullspace dimension */
1463  maxf = 100; /* maximal filter length */
1464  mlp = 100; /* maximum level of degeneracy */
1465 
1466  /* TODO eventually, the output should be redirected to the message handler,
1467  * but even to just redirect to some other file, we would have to open the output-unit in Fortran
1468  */
1469  nout = 6; /* output to stdout for now */
1470  ifail = problem->warmstart ? -1 : 0; /* -1 for warmstart, otherwise 0 */
1471  rho = 10.0; /* initial trust-region radius */
1472 
1473  user = (real*)data;
1474  iuser = (fint*)problem;
1475  if( problem->warmstart ) /* if warmstart, then need to keep istat[0] */
1476  memset(problem->istat+1, 0, sizeof(problem->istat)-sizeof(*problem->istat));
1477  else
1478  memset(problem->istat, 0, sizeof(problem->istat));
1479  memset(problem->rstat, 0, sizeof(problem->rstat));
1480  problem->niterations = 0;
1481 
1482  if( problem->x == NULL )
1483  {
1484  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->x, problem->varssize) );
1485  }
1486  if( problem->c == NULL )
1487  {
1488  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->c, problem->conssize) );
1489  }
1490  if( problem->lam == NULL )
1491  {
1492  SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &problem->lam, problem->varssize + problem->conssize) ); /*lint !e776 */
1493  }
1494  else
1495  {
1496  BMSclearMemoryArray(problem->lam, problem->varssize + problem->conssize);
1497  }
1498  if( problem->s == NULL )
1499  {
1500  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->s, problem->varssize + problem->conssize) );
1501  }
1502 
1503  if( problem->la == NULL )
1504  {
1505  /* allocate la, a and initialize la for Objective Gradient and Jacobian */
1506  SCIP_CALL( setupGradients(scip, problem->oracle, &problem->la, &problem->lasize, &problem->a) );
1507  }
1508  /* maximal number entries in a = nvars+nnz */
1509  maxa = problem->la[0]-1;
1510 
1511  if( problem->hessiannz == NULL )
1512  {
1513  /* allocate and initialize problem->hessiannz for Hessian */
1514  SCIP_CALL( setupHessian(scip, problem->oracle, &problem->hessiannz, &problem->hessiannzsize) );
1515  }
1516 
1517  /* setup variable bounds, constraint sides, and constraint types */
1518  if( problem->bl == NULL )
1519  {
1520  assert(problem->bu == NULL);
1521  assert(problem->cstype == NULL);
1522 
1523  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->bl, problem->varssize + problem->conssize) );
1524  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->bu, problem->varssize + problem->conssize) );
1525  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->cstype, problem->conssize) );
1526 
1527  BMScopyMemoryArray(problem->bl, SCIPnlpiOracleGetVarLbs(problem->oracle), n);
1528  BMScopyMemoryArray(problem->bu, SCIPnlpiOracleGetVarUbs(problem->oracle), n);
1529  for( i = 0; i < m; ++i )
1530  {
1531  problem->bl[n+i] = SCIPnlpiOracleGetConstraintLhs(problem->oracle, i);
1532  problem->bu[n+i] = SCIPnlpiOracleGetConstraintRhs(problem->oracle, i);
1533  problem->cstype[i] = SCIPnlpiOracleIsConstraintNonlinear(problem->oracle, i) ? 'N' : 'L';
1534  }
1535  }
1536 
1537  /* buffer for evaluation results (used in setupStart already) */
1538  SCIP_CALL( SCIPensureBlockMemoryArray(scip, &problem->evalbuffer, &problem->evalbufsize, MAX3(n, problem->hessiannz[0], maxa)) );
1539 
1540  /* setup starting point */
1541  SCIP_CALL( setupStart(data, problem, problem->x, &success) );
1542  if( !success )
1543  {
1544  /* FilterSQP would crash if starting point cannot be evaluated, so give up */
1545  SCIP_CALL( processSolveOutcome(data, problem, 7, param.feastol, param.opttol, NULL, NULL) );
1546  return SCIP_OKAY;
1547  }
1548 
1549  /* setup workspace */
1550  /* initial guess of real workspace size */
1551  /* FilterSQP manual: mxwk = 21*n + 8*m + mlp + 8*maxf + kmax*(kmax+9)/2 + nprof, with nprof = 20*n as a good guess */
1552  /* Bonmin: mxwk = 21*n + 8*m + mlp + 8*maxf + lh1 + kmax*(kmax+9)/2 + mxwk0,
1553  * with lh1 = nnz_h+8+2*n+m and mxwk0 = 2000000 (parameter) */
1554  lh1 = problem->hessiannz[0]-1 + 8 + 2*n + m;
1555  minmxwk = 21*n + 8*m + mlp + 8*maxf + lh1 + kmax*(kmax+9)/2 + MAX(20*n,2000);
1556  SCIP_CALL( SCIPensureBlockMemoryArray(scip, &problem->ws, &problem->mxwk, minmxwk) );
1557 
1558  /* initial guess of integer workspace size */
1559  /* FilterSQP manual: mxiwk = 13*n + 4*m + mlp + 100 + kmax */
1560  /* Bonmin: mxiwk = 13*n + 4*m + mlp + lh1 + kmax + 113 + mxiwk0, with mxiwk0 = 500000 (parameter) */
1561  minmxiwk = 13*n + 4*m + mlp + lh1 + 100 + kmax + 113 + MAX(5*n,5000);
1562  if( !problem->warmstart ) /* if warmstart, then lws should remain untouched (n and m didn't change anyway) */
1563  {
1564  SCIP_CALL( SCIPensureBlockMemoryArray(scip, &problem->lws, &problem->mxiwk, minmxiwk) );
1565  }
1566  assert(problem->lws != NULL);
1567  /* in case of some evalerrors, not clearing ws could lead to valgrind warnings about use of uninitialized memory */
1568  BMSclearMemoryArray(problem->ws, problem->mxwk);
1569 
1570  /* from here on we are not thread-safe: if intended for multithread use, then protect filtersqp call with mutex
1571  * NOTE: we need to make sure that we do not return from nlpiSolve before unlocking the mutex
1572  */
1573 #ifdef SCIP_THREADSAFE
1574  (void) pthread_mutex_lock(&filtersqpmutex);
1575 #endif
1576 
1577  /* initialize global variables from filtersqp */
1578  /* FilterSQP eps is tolerance for both feasibility and optimality, and also for trust-region radius, etc. */
1579  F77_FUNC(nlp_eps_inf,NLP_EPS_INF).eps = MIN(param.feastol, param.opttol * OPTTOLFACTOR);
1580  F77_FUNC(nlp_eps_inf,NLP_EPS_INF).infty = SCIPinfinity(scip);
1581  F77_FUNC(ubdc,UBDC).ubd = 100.0;
1582  F77_FUNC(ubdc,UBDC).tt = 1.25;
1583  F77_FUNC(scalec,SCALEC).scale_mode = 0;
1584 
1585  for( nruns = 1; ; ++nruns )
1586  {
1587  maxiter = param.iterlimit - problem->niterations;
1588 
1589  F77_FUNC(filtersqp,FILTERSQP)(
1590  &n, &m, &kmax, &maxa,
1591  &maxf, &mlp, &problem->mxwk, &problem->mxiwk,
1592  &iprint, &nout, &ifail, &rho,
1593  problem->x, problem->c, &f, &problem->fmin, problem->bl,
1594  problem->bu, problem->s, problem->a, problem->la, problem->ws,
1595  problem->lws, problem->lam, problem->cstype, user,
1596  iuser, &maxiter, problem->istat,
1597  problem->rstat, cstype_len);
1598 
1599  problem->niterations += problem->istat[1];
1600 
1601  assert(ifail <= 10);
1602  /* if ifail >= 8 (probably the workspace was too small), then retry with larger workspace
1603  * if ifail == 0 (local optimal), but absolute violation of KKT too large, then retry with small eps
1604  */
1605  if( ifail < 8 && (ifail != 0 || problem->rstat[0] <= param.opttol) )
1606  break;
1607 
1608  if( param.verblevel > 0 )
1609  {
1610  SCIPinfoMessage(scip, NULL, "FilterSQP terminated with status %d in run %d, absolute KKT violation is %g\n", ifail, nruns, problem->rstat[0]);
1611  }
1612 
1613  /* if iteration or time limit exceeded or solve is interrupted, then don't retry */
1614  if( problem->niterations >= param.iterlimit || SCIPisSolveInterrupted(scip) || timelimitreached(data, problem) )
1615  {
1616  if( param.verblevel > 0 )
1617  {
1618  SCIPinfoMessage(scip, NULL, "Time or iteration limit reached or interrupted, not retrying\n");
1619  }
1620  break;
1621  }
1622 
1623  /* if maximal number of runs reached, then stop */
1624  if( nruns >= MAXNRUNS )
1625  {
1626  if( param.verblevel > 0 )
1627  {
1628  SCIPinfoMessage(scip, NULL, "Run limit reached, not retrying\n");
1629  }
1630  break;
1631  }
1632 
1633  if( ifail == 0 )
1634  {
1635  SCIP_Real epsfactor;
1636 
1637  if( F77_FUNC(nlp_eps_inf,NLP_EPS_INF).eps <= MINEPS )
1638  {
1639  if( param.verblevel > 0 )
1640  {
1641  SCIPinfoMessage(scip, NULL, "Already reached minimal epsilon, not retrying\n");
1642  }
1643  break;
1644  }
1645 
1646  epsfactor = param.opttol / problem->rstat[0];
1647  assert(epsfactor < 1.0); /* because of the if's above */
1648  epsfactor *= OPTTOLFACTOR;
1649 
1650  F77_FUNC(nlp_eps_inf,NLP_EPS_INF).eps = MAX(MINEPS, F77_FUNC(nlp_eps_inf,NLP_EPS_INF).eps * epsfactor);
1651  if( param.verblevel > 0 )
1652  {
1653  SCIPinfoMessage(scip, NULL, "Continue with eps = %g\n", F77_FUNC(nlp_eps_inf,NLP_EPS_INF).eps);
1654  }
1655  ifail = -1; /* do warmstart */
1656 
1657  continue;
1658  }
1659 
1660  /* increase real workspace, if ifail = 9 (real workspace too small) or ifail = 8 (unexpected ifail from QP solver, often also when workspace too small) */
1661  if( ifail == 8 || ifail == 9 )
1662  {
1663  int newsize = SCIPcalcMemGrowSize(scip, WORKSPACEGROWTHFACTOR*problem->mxwk);
1664  if( BMSreallocBlockMemoryArray(SCIPblkmem(scip), &problem->ws, problem->mxwk, newsize) == NULL )
1665  {
1666  /* realloc failed: give up NLP solve */
1667  problem->mxwk = 0;
1668  break;
1669  }
1670  problem->mxwk = newsize;
1671  }
1672 
1673  /* increase integer workspace, if ifail = 10 (integer workspace too small) or ifail = 8 (unexpected ifail from QP solver, often also when workspace too small) */
1674  if( ifail == 8 || ifail == 10 )
1675  {
1676  int newsize = SCIPcalcMemGrowSize(scip, WORKSPACEGROWTHFACTOR*problem->mxiwk);
1677  if( BMSreallocBlockMemoryArray(SCIPblkmem(scip), &problem->lws, problem->mxiwk, newsize) == NULL )
1678  {
1679  /* realloc failed: give up NLP solve */
1680  problem->mxiwk = 0;
1681  break;
1682  }
1683  problem->mxiwk = newsize;
1684 
1685  /* better don't try warmstart for the next trial; warmstart requires that lws is untouched, does extending count as touching? */
1686  ifail = 0;
1687  }
1688 
1689  /* reset starting point, in case it was overwritten by failed solve (return can only be SCIP_OKAY, because randnumgen must exist already)
1690  * NOTE: If initguess is NULL (no user-given starting point), then this will result in a slightly different starting point as in the previous setupStart() call (random numbers)
1691  * However, as no point was given, it shouldn't matter which point we actually start from.
1692  */
1693  (void) setupStart(data, problem, problem->x, &success);
1694  assert(success);
1695  }
1696 
1697 #ifdef SCIP_THREADSAFE
1698  (void) pthread_mutex_unlock(&filtersqpmutex);
1699 #endif
1700 
1701  SCIP_CALL( processSolveOutcome(data, problem, ifail, param.feastol, param.opttol, problem->x, problem->lam) );
1702 
1703  return SCIP_OKAY; /*lint !e527*/
1704 } /*lint !e715*/
1705 
1706 /** gives solution status */
1707 static
1708 SCIP_DECL_NLPIGETSOLSTAT(nlpiGetSolstatFilterSQP)
1709 {
1710  assert(problem != NULL);
1711 
1712  return problem->solstat;
1713 } /*lint !e715*/
1714 
1715 /** gives termination reason */
1716 static
1717 SCIP_DECL_NLPIGETTERMSTAT(nlpiGetTermstatFilterSQP)
1718 {
1719  assert(problem != NULL);
1720 
1721  return problem->termstat;
1722 } /*lint !e715*/
1723 
1724 /** gives primal and dual solution values */
1725 static
1726 SCIP_DECL_NLPIGETSOLUTION(nlpiGetSolutionFilterSQP)
1727 {
1728  assert(problem != NULL);
1729 
1730  if( primalvalues != NULL )
1731  {
1732  assert(problem->primalvalues != NULL);
1733 
1734  *primalvalues = problem->primalvalues;
1735  }
1736 
1737  if( consdualvalues != NULL )
1738  {
1739  assert(problem->consdualvalues != NULL);
1740 
1741  *consdualvalues = problem->consdualvalues;
1742  }
1743 
1744  if( varlbdualvalues != NULL )
1745  {
1746  assert(problem->varlbdualvalues != NULL);
1747 
1748  *varlbdualvalues = problem->varlbdualvalues;
1749  }
1750 
1751  if( varubdualvalues != NULL )
1752  {
1753  assert(problem->varubdualvalues != NULL);
1754 
1755  *varubdualvalues = problem->varubdualvalues;
1756  }
1757 
1758  if( objval != NULL )
1759  {
1760  if( problem->primalvalues != NULL )
1761  {
1762  /* TODO store last solution value instead of reevaluating the objective function */
1763  SCIP_CALL( SCIPnlpiOracleEvalObjectiveValue(scip, problem->oracle, problem->primalvalues, objval) );
1764  }
1765  else
1766  *objval = SCIP_INVALID;
1767  }
1768 
1769  return SCIP_OKAY; /*lint !e527*/
1770 } /*lint !e715*/
1771 
1772 /** gives solve statistics */
1773 static
1774 SCIP_DECL_NLPIGETSTATISTICS(nlpiGetStatisticsFilterSQP)
1775 {
1776  assert(problem != NULL);
1777  assert(statistics != NULL);
1778 
1779  statistics->niterations = problem->niterations;
1780  statistics->totaltime = problem->solvetime;
1781  statistics->evaltime = SCIPnlpiOracleGetEvalTime(scip, problem->oracle);
1782  statistics->consviol = problem->rstat[4];
1783  statistics->boundviol = 0.0;
1784 
1785  return SCIP_OKAY; /*lint !e527*/
1786 } /*lint !e715*/
1787 
1788 /*
1789  * NLP solver interface specific interface methods
1790  */
1791 
1792 /** create solver interface for filterSQP solver and include it into SCIP, if filterSQP is available */
1794  SCIP* scip /**< SCIP data structure */
1795  )
1796 {
1797  SCIP_NLPIDATA* nlpidata;
1798 
1799  /* create filterSQP solver interface data */
1800  SCIP_CALL( SCIPallocClearBlockMemory(scip, &nlpidata) );
1801 
1802  /* create solver interface */
1803  SCIP_CALL( SCIPincludeNlpi(scip,
1805  nlpiCopyFilterSQP, nlpiFreeFilterSQP, NULL,
1806  nlpiCreateProblemFilterSQP, nlpiFreeProblemFilterSQP, NULL,
1807  nlpiAddVarsFilterSQP, nlpiAddConstraintsFilterSQP, nlpiSetObjectiveFilterSQP,
1808  nlpiChgVarBoundsFilterSQP, nlpiChgConsSidesFilterSQP, nlpiDelVarSetFilterSQP, nlpiDelConstraintSetFilterSQP,
1809  nlpiChgLinearCoefsFilterSQP, nlpiChgExprFilterSQP,
1810  nlpiChgObjConstantFilterSQP, nlpiSetInitialGuessFilterSQP, nlpiSolveFilterSQP, nlpiGetSolstatFilterSQP, nlpiGetTermstatFilterSQP,
1811  nlpiGetSolutionFilterSQP, nlpiGetStatisticsFilterSQP,
1812  nlpidata) );
1813 
1815 
1816  return SCIP_OKAY;
1817 }
1818 
1819 /** gets string that identifies filterSQP */
1821  void
1822  )
1823 {
1824  return "filterSQP"; /* TODO version number? */
1825 }
1826 
1827 /** gets string that describes filterSQP */
1829  void
1830  )
1831 {
1832  return NLPI_DESC;
1833 }
1834 
1835 /** returns whether filterSQP is available, i.e., whether it has been linked in */
1837  void
1838  )
1839 {
1840  return TRUE;
1841 }
static SCIP_DECL_NLPIFREEPROBLEM(nlpiFreeProblemFilterSQP)
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:101
static SCIP_RETCODE processSolveOutcome(SCIP_NLPIDATA *nlpidata, SCIP_NLPIPROBLEM *problem, fint ifail, SCIP_Real feastol, SCIP_Real opttol, real *x, real *lam)
SCIP_RETCODE SCIPnlpiOracleEvalObjectiveGradient(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Bool isnewx, SCIP_Real *objval, SCIP_Real *objgrad)
Definition: nlpioracle.c:1959
#define SCIPreallocBlockMemoryArray(scip, ptr, oldnum, newnum)
Definition: scip_mem.h:90
enum SCIP_NlpTermStat SCIP_NLPTERMSTAT
Definition: type_nlpi.h:185
static SCIP_DECL_NLPIGETSOLSTAT(nlpiGetSolstatFilterSQP)
SCIP_NLPIDATA * SCIPnlpiGetData(SCIP_NLPI *nlpi)
Definition: nlpi.c:683
SCIP_RETCODE SCIPnlpiOracleChgLinearCoefs(SCIP *scip, SCIP_NLPIORACLE *oracle, int considx, int nentries, const int *varidxs, const SCIP_Real *newcoefs)
Definition: nlpioracle.c:1548
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:84
SCIP_RETCODE SCIPincludeNlpSolverFilterSQP(SCIP *scip)
public methods for memory management
SCIP_RETCODE SCIPnlpiOracleEvalConstraintValue(SCIP *scip, SCIP_NLPIORACLE *oracle, int considx, const SCIP_Real *x, SCIP_Real *conval)
Definition: nlpioracle.c:1903
static SCIP_DECL_NLPIADDCONSTRAINTS(nlpiAddConstraintsFilterSQP)
static SCIP_DECL_NLPICHGVARBOUNDS(nlpiChgVarBoundsFilterSQP)
#define SCIPallocClearBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:88
int SCIPcalcMemGrowSize(SCIP *scip, int num)
Definition: scip_mem.c:130
#define NLPI_PRIORITY
fint n_bqpd_calls
public solving methods
SCIP_Real * varubdualvalues
static SCIP_DECL_NLPISOLVE(nlpiSolveFilterSQP)
methods to store an NLP and request function, gradient, and Hessian values
SCIP_RETCODE SCIPnlpiOracleEvalHessianLag(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Bool isnewx_obj, SCIP_Bool isnewx_cons, SCIP_Real objfactor, const SCIP_Real *lambda, SCIP_Real *hessian)
Definition: nlpioracle.c:2371
const char * SCIPgetSolverDescFilterSQP(void)
SCIP_NLPIORACLE * oracle
#define FALSE
Definition: def.h:87
SCIP_RETCODE SCIPnlpiOracleAddConstraints(SCIP *scip, SCIP_NLPIORACLE *oracle, int nconss, const SCIP_Real *lhss, const SCIP_Real *rhss, const int *nlininds, int *const *lininds, SCIP_Real *const *linvals, SCIP_EXPR **exprs, const char **consnames)
Definition: nlpioracle.c:1158
static SCIP_DECL_NLPICHGCONSSIDES(nlpiChgConsSidesFilterSQP)
SCIP_Real SCIPinfinity(SCIP *scip)
#define TRUE
Definition: def.h:86
#define NLPI_NAME
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
void F77_FUNC(filtersqp, FILTERSQP)
static SCIP_DECL_NLPIADDVARS(nlpiAddVarsFilterSQP)
SCIP_Real * primalvalues
SCIP_RETCODE SCIPnlpiOracleSetProblemName(SCIP *scip, SCIP_NLPIORACLE *oracle, const char *name)
Definition: nlpioracle.c:1036
#define WORKSPACEGROWTHFACTOR
SCIP_Real SCIPnlpiOracleGetEvalTime(SCIP *scip, SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:2434
static SCIP_TIME gettime(void)
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:116
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:99
static SCIP_DECL_NLPIGETSOLUTION(nlpiGetSolutionFilterSQP)
SCIP_RETCODE SCIPnlpiOracleSetObjective(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real constant, int nlin, const int *lininds, const SCIP_Real *linvals, SCIP_EXPR *expr)
Definition: nlpioracle.c:1219
#define OPTTOLFACTOR
SCIP_Real SCIPnlpiOracleGetConstraintRhs(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1812
#define SCIPdebugMsg
Definition: scip_message.h:69
SCIP_VAR ** x
Definition: circlepacking.c:54
filterSQP NLP interface
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_Real timelimit
Definition: type_nlpi.h:63
real eps
SCIP_Real * varlbdualvalues
static SCIP_DECL_NLPIDELVARSET(nlpiDelVarSetFilterSQP)
public methods for numerical tolerances
static SCIP_DECL_NLPIGETTERMSTAT(nlpiGetTermstatFilterSQP)
SCIP_Bool SCIPisFilterSQPAvailableFilterSQP(void)
public methods for NLPI solver interfaces
long ftnlen
static SCIP_DECL_NLPISETINITIALGUESS(nlpiSetInitialGuessFilterSQP)
#define BMSfreeMemoryArray(ptr)
Definition: memory.h:140
#define SCIPerrorMessage
Definition: pub_message.h:55
int SCIPnlpiOracleGetNConstraints(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1717
struct SCIP_NlpiData SCIP_NLPIDATA
Definition: type_nlpi.h:43
enum SCIP_NlpSolStat SCIP_NLPSOLSTAT
Definition: type_nlpi.h:159
#define SCIP_NLPPARAM_PRINT(param)
Definition: type_nlpi.h:133
int SCIPnlpiOracleGetNVars(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1707
SCIP_RETCODE SCIPnlpiOracleResetEvalTime(SCIP *scip, SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:2418
SCIP_RETCODE SCIPnlpiOracleChgVarBounds(SCIP *scip, SCIP_NLPIORACLE *oracle, int nvars, const int *indices, const SCIP_Real *lbs, const SCIP_Real *ubs)
Definition: nlpioracle.c:1248
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition: scip_mem.c:48
SCIP_Bool warmstart
fint scale_mode
#define NULL
Definition: lpi_spx1.cpp:155
int fint
fint n_bqpd_prfint
SCIP_RETCODE SCIPnlpiOracleGetJacobianSparsity(SCIP *scip, SCIP_NLPIORACLE *oracle, const int **offset, const int **col)
Definition: nlpioracle.c:2018
SCIP_RETCODE SCIPnlpiOracleAddVars(SCIP *scip, SCIP_NLPIORACLE *oracle, int nvars, const SCIP_Real *lbs, const SCIP_Real *ubs, const char **varnames)
Definition: nlpioracle.c:1072
#define SCIP_CALL(x)
Definition: def.h:384
static SCIP_DECL_NLPICOPY(nlpiCopyFilterSQP)
#define SCIPensureBlockMemoryArray(scip, ptr, arraysizeptr, minsize)
Definition: scip_mem.h:98
SCIP_RETCODE SCIPnlpiOracleChgExpr(SCIP *scip, SCIP_NLPIORACLE *oracle, int considx, SCIP_EXPR *expr)
Definition: nlpioracle.c:1644
#define MAXNRUNS
static SCIP_DECL_NLPIFREE(nlpiFreeFilterSQP)
SCIP_Real SCIPnlpiOracleGetConstraintLhs(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1799
static SCIP_RETCODE handleNlpParam(SCIP *scip, SCIP_NLPIPROBLEM *nlpiproblem, const SCIP_NLPPARAM param)
static SCIP_Real timeelapsed(SCIP_NLPIDATA *nlpidata)
real infty
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
fint phr
public data structures and miscellaneous methods
SCIP_RETCODE SCIPnlpiOracleEvalJacobian(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Bool isnewx, SCIP_Real *convals, SCIP_Real *jacobi)
Definition: nlpioracle.c:2150
#define MINEPS
SCIP_RETCODE SCIPnlpiOracleDelVarSet(SCIP *scip, SCIP_NLPIORACLE *oracle, int *delstats)
Definition: nlpioracle.c:1320
static SCIP_Bool timelimitreached(SCIP_NLPIDATA *nlpidata, SCIP_NLPIPROBLEM *nlpiproblem)
#define SCIP_Bool
Definition: def.h:84
fint phl
#define RANDSEED
real ubd
#define MAX(x, y)
Definition: tclique_def.h:83
SCIP_RETCODE SCIPnlpiOracleCreate(SCIP *scip, SCIP_NLPIORACLE **oracle)
Definition: nlpioracle.c:974
SCIP_Real lobjlimit
Definition: type_nlpi.h:59
static SCIP_DECL_NLPICHGEXPR(nlpiChgExprFilterSQP)
const SCIP_Real * SCIPnlpiOracleGetVarLbs(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1727
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:127
#define NLPI_DESC
SCIP_Real * initguess
SCIP_Bool SCIPnlpiOracleIsConstraintNonlinear(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1838
static SCIP_DECL_NLPIGETSTATISTICS(nlpiGetStatisticsFilterSQP)
const SCIP_Real * SCIPnlpiOracleGetVarUbs(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1737
SCIP_RETCODE SCIPnlpiOracleFree(SCIP *scip, SCIP_NLPIORACLE **oracle)
Definition: nlpioracle.c:1004
SCIP_Real SCIPrandomGetReal(SCIP_RANDNUMGEN *randnumgen, SCIP_Real minrandval, SCIP_Real maxrandval)
Definition: misc.c:10025
#define SCIP_REAL_MAX
Definition: def.h:178
#define SCIP_REAL_MIN
Definition: def.h:179
general public methods
SCIP_RETCODE SCIPnlpiOracleEvalObjectiveValue(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Real *objval)
Definition: nlpioracle.c:1878
static SCIP_DECL_NLPIDELCONSSET(nlpiDelConstraintSetFilterSQP)
public methods for random numbers
time_t sec
static SCIP_DECL_NLPICHGOBJCONSTANT(nlpiChgObjConstantFilterSQP)
#define MAXPERTURB
real tt
SCIP_VAR * a
Definition: circlepacking.c:57
#define SCIP_Real
Definition: def.h:177
SCIP_RETCODE SCIPnlpiOracleGetHessianLagSparsity(SCIP *scip, SCIP_NLPIORACLE *oracle, const int **offset, const int **col)
Definition: nlpioracle.c:2277
static SCIP_DECL_NLPICREATEPROBLEM(nlpiCreateProblemFilterSQP)
static void invalidateSolution(SCIP_NLPIPROBLEM *problem)
public methods for message handling
#define SCIP_INVALID
Definition: def.h:197
static SCIP_RETCODE setupStart(SCIP_NLPIDATA *data, SCIP_NLPIPROBLEM *problem, real *x, SCIP_Bool *success)
#define SCIPisFinite(x)
Definition: pub_misc.h:1892
SCIP_RETCODE SCIPincludeExternalCodeInformation(SCIP *scip, const char *name, const char *description)
Definition: scip_general.c:704
SCIP_Real * consdualvalues
SCIP_NLPTERMSTAT termstat
SCIP_NLPSOLSTAT solstat
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:102
SCIP_RETCODE SCIPincludeNlpi(SCIP *scip, const char *name, const char *description, int priority, SCIP_DECL_NLPICOPY((*nlpicopy)), SCIP_DECL_NLPIFREE((*nlpifree)), SCIP_DECL_NLPIGETSOLVERPOINTER((*nlpigetsolverpointer)), SCIP_DECL_NLPICREATEPROBLEM((*nlpicreateproblem)), SCIP_DECL_NLPIFREEPROBLEM((*nlpifreeproblem)), SCIP_DECL_NLPIGETPROBLEMPOINTER((*nlpigetproblempointer)), SCIP_DECL_NLPIADDVARS((*nlpiaddvars)), SCIP_DECL_NLPIADDCONSTRAINTS((*nlpiaddconstraints)), SCIP_DECL_NLPISETOBJECTIVE((*nlpisetobjective)), SCIP_DECL_NLPICHGVARBOUNDS((*nlpichgvarbounds)), SCIP_DECL_NLPICHGCONSSIDES((*nlpichgconssides)), SCIP_DECL_NLPIDELVARSET((*nlpidelvarset)), SCIP_DECL_NLPIDELCONSSET((*nlpidelconsset)), SCIP_DECL_NLPICHGLINEARCOEFS((*nlpichglinearcoefs)), SCIP_DECL_NLPICHGEXPR((*nlpichgexpr)), SCIP_DECL_NLPICHGOBJCONSTANT((*nlpichgobjconstant)), SCIP_DECL_NLPISETINITIALGUESS((*nlpisetinitialguess)), SCIP_DECL_NLPISOLVE((*nlpisolve)), SCIP_DECL_NLPIGETSOLSTAT((*nlpigetsolstat)), SCIP_DECL_NLPIGETTERMSTAT((*nlpigettermstat)), SCIP_DECL_NLPIGETSOLUTION((*nlpigetsolution)), SCIP_DECL_NLPIGETSTATISTICS((*nlpigetstatistics)), SCIP_NLPIDATA *nlpidata)
Definition: scip_nlpi.c:98
static SCIP_RETCODE setupGradients(SCIP *scip, SCIP_NLPIORACLE *oracle, fint **la, int *lasize, real **a)
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:123
#define SCIPallocClearBlockMemory(scip, ptr)
Definition: scip_mem.h:82
static SCIP_RETCODE setupHessian(SCIP *scip, SCIP_NLPIORACLE *oracle, fint **la, int *lasize)
static SCIP_DECL_NLPISETOBJECTIVE(nlpiSetObjectiveFilterSQP)
SCIP_RETCODE SCIPnlpiOracleChgConsSides(SCIP *scip, SCIP_NLPIORACLE *oracle, int nconss, const int *indices, const SCIP_Real *lhss, const SCIP_Real *rhss)
Definition: nlpioracle.c:1285
fint phc
fint phe
double real
SCIP_RETCODE SCIPnlpiOracleDelConsSet(SCIP *scip, SCIP_NLPIORACLE *oracle, int *delstats)
Definition: nlpioracle.c:1462
#define BMSreallocBlockMemoryArray(mem, ptr, oldnum, newnum)
Definition: memory.h:451
SCIP_NLPPARAM_FASTFAIL fastfail
Definition: type_nlpi.h:66
static SCIP_DECL_NLPICHGLINEARCOEFS(nlpiChgLinearCoefsFilterSQP)
const char * SCIPgetSolverNameFilterSQP(void)
SCIP_Bool SCIPisSolveInterrupted(SCIP *scip)
Definition: scip_solve.c:3548
SCIP_RETCODE SCIPnlpiOracleChgObjConstant(SCIP *scip, SCIP_NLPIORACLE *oracle, SCIP_Real objconstant)
Definition: nlpioracle.c:1690