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