Scippy

SCIP

Solving Constraint Integer Programs

cons_quadratic.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-2019 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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_quadratic.c
17  * @brief constraint handler for quadratic constraints \f$\textrm{lhs} \leq \sum_{i,j=1}^n a_{i,j} x_i x_j + \sum_{i=1}^n b_i x_i \leq \textrm{rhs}\f$
18  * @author Stefan Vigerske
19  * @author Benjamin Mueller
20  * @author Felipe Serrano
21  *
22  * @todo SCIP might fix linear variables on +/- infinity; remove them in presolve and take care later
23  * @todo round constraint sides to integers if all coefficients and variables are (impl.) integer
24  * @todo constraints in one variable should be replaced by linear or bounddisjunction constraint
25  * @todo check if some quadratic terms appear in several constraints and try to simplify (e.g., nous1)
26  * @todo skip separation in enfolp if for current LP (check LP id) was already separated
27  * @todo watch unbounded variables to enable/disable propagation
28  * @todo sort order in bilinvar1/bilinvar2 such that the var which is involved in more terms is in bilinvar1, and use this info propagate and AddLinearReform
29  * @todo underestimate for multivariate concave quadratic terms as in cons_nonlinear
30  */
31 
32 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
33 
34 #define SCIP_PRIVATE_ROWPREP
35 
36 #include "blockmemshell/memory.h"
37 #include <ctype.h>
38 #include "nlpi/nlpi.h"
39 #include "nlpi/nlpi_ipopt.h"
40 #include "nlpi/pub_expr.h"
41 #include "nlpi/type_expr.h"
42 #include "scip/cons_and.h"
43 #include "scip/cons_bounddisjunction.h"
44 #include "scip/cons_linear.h"
45 #include "scip/cons_nonlinear.h"
46 #include "scip/cons_quadratic.h"
47 #include "scip/cons_varbound.h"
48 #include "scip/debug.h"
49 #include "scip/heur_subnlp.h"
50 #include "scip/heur_trysol.h"
51 #include "scip/intervalarith.h"
52 #include "scip/pub_cons.h"
53 #include "scip/pub_event.h"
54 #include "scip/pub_heur.h"
55 #include "scip/pub_lp.h"
56 #include "scip/pub_message.h"
57 #include "scip/pub_misc.h"
58 #include "scip/pub_misc_sort.h"
59 #include "scip/pub_nlp.h"
60 #include "scip/pub_sol.h"
61 #include "scip/pub_tree.h"
62 #include "scip/pub_var.h"
63 #include "scip/scip_branch.h"
64 #include "scip/scip_cons.h"
65 #include "scip/scip_copy.h"
66 #include "scip/scip_cut.h"
67 #include "scip/scip_event.h"
68 #include "scip/scip_general.h"
69 #include "scip/scip_heur.h"
70 #include "scip/scip_lp.h"
71 #include "scip/scip_mem.h"
72 #include "scip/scip_message.h"
73 #include "scip/scip_nlp.h"
74 #include "scip/scip_nonlinear.h"
75 #include "scip/scip_numerics.h"
76 #include "scip/scip_param.h"
77 #include "scip/scip_prob.h"
78 #include "scip/scip_probing.h"
79 #include "scip/scip_sepa.h"
80 #include "scip/scip_sol.h"
81 #include "scip/scip_solve.h"
82 #include "scip/scip_solvingstats.h"
83 #include "scip/scip_tree.h"
84 #include "scip/scip_var.h"
85 #include <string.h>
86 
87 /* constraint handler properties */
88 #define CONSHDLR_NAME "quadratic"
89 #define CONSHDLR_DESC "quadratic constraints of the form lhs <= b' x + x' A x <= rhs"
90 #define CONSHDLR_SEPAPRIORITY 10 /**< priority of the constraint handler for separation */
91 #define CONSHDLR_ENFOPRIORITY -50 /**< priority of the constraint handler for constraint enforcing */
92 #define CONSHDLR_CHECKPRIORITY -4000000 /**< priority of the constraint handler for checking feasibility */
93 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
94 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
95 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
96  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
97 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
98 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
99 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
100 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
101 
102 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP /**< propagation timing mask of the constraint handler */
103 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_ALWAYS /**< presolving timing of the constraint handler (fast, medium, or exhaustive) */
104 
105 #define MAXDNOM 10000LL /**< maximal denominator for simple rational fixed values */
106 #define NONLINCONSUPGD_PRIORITY 40000 /**< priority of upgrading nonlinear constraints */
107 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
108 
109 /* Activating this define enables reformulation of bilinear terms x*y with implications from x to y into linear terms.
110  * However, implications are not enforced by SCIP. Thus, if, e.g., the used implication was derived from this constraint and we then reformulate the constraint,
111  * then the implication may not be enforced in a solution.
112  * This issue need to be fixed before this feature can be enabled.
113  */
114 /* #define CHECKIMPLINBILINEAR */
115 
116 /* enable new propagation for bivariate quadratic terms */
117 #define PROPBILINNEW
118 
119 /* epsilon for differentiating between a boundary and interior point */
120 #define INTERIOR_EPS 1e-1
121 
122 /* scaling factor for gauge function */
123 #define GAUGESCALE 0.99999
124 
125 #define ROWPREP_SCALEUP_VIOLNONZERO (10.0*SCIPepsilon(scip)) /**< minimal violation for considering up-scaling of rowprep (we want to avoid upscaling very small violations) */
126 #define ROWPREP_SCALEUP_MINVIOLFACTOR 2.0 /**< scale up will target a violation of ~MINVIOLFACTOR*minviol, where minviol is given by caller */
127 #define ROWPREP_SCALEUP_MAXMINCOEF (1.0 / SCIPfeastol(scip)) /**< scale up only if min. coef is below this number (before scaling) */
128 #define ROWPREP_SCALEUP_MAXMAXCOEF SCIPgetHugeValue(scip) /**< scale up only if max. coef will not exceed this number by scaling */
129 #define ROWPREP_SCALEUP_MAXSIDE SCIPgetHugeValue(scip) /**< scale up only if side will not exceed this number by scaling */
130 #define ROWPREP_SCALEDOWN_MINMAXCOEF (1.0 / SCIPfeastol(scip)) /**< scale down if max. coef is at least this number (before scaling) */
131 #define ROWPREP_SCALEDOWN_MINCOEF SCIPfeastol(scip) /**< scale down only if min. coef does not drop below this number by scaling */
132 
133 /*
134  * Data structures
135  */
136 
137 /** eventdata for variable bound change events in quadratic constraints */
138 struct SCIP_QuadVarEventData
139 {
140  SCIP_CONS* cons; /**< the constraint */
141  int varidx; /**< the index of the variable which bound change is caught, positive for linear variables, negative for quadratic variables */
142  int filterpos; /**< position of eventdata in SCIP's event filter */
143 };
144 
145 /** Data of a quadratic constraint. */
146 struct SCIP_ConsData
147 {
148  SCIP_Real lhs; /**< left hand side of constraint */
149  SCIP_Real rhs; /**< right hand side of constraint */
150 
151  int nlinvars; /**< number of linear variables */
152  int linvarssize; /**< length of linear variable arrays */
153  SCIP_VAR** linvars; /**< linear variables */
154  SCIP_Real* lincoefs; /**< coefficients of linear variables */
155  SCIP_QUADVAREVENTDATA** lineventdata; /**< eventdata for bound change of linear variable */
156 
157  int nquadvars; /**< number of variables in quadratic terms */
158  int quadvarssize; /**< length of quadratic variable terms arrays */
159  SCIP_QUADVARTERM* quadvarterms; /**< array with quadratic variable terms */
160 
161  int nbilinterms; /**< number of bilinear terms */
162  int bilintermssize; /**< length of bilinear term arrays */
163  SCIP_BILINTERM* bilinterms; /**< bilinear terms array */
164  int* bilintermsidx; /**< unique index of each bilinear term xy in the bilinestimators array of the constraint handler data */
165 
166  SCIP_NLROW* nlrow; /**< a nonlinear row representation of this constraint */
167 
168  unsigned int linvarssorted:1; /**< are the linear variables already sorted? */
169  unsigned int linvarsmerged:1; /**< are equal linear variables already merged? */
170  unsigned int quadvarssorted:1; /**< are the quadratic variables already sorted? */
171  unsigned int quadvarsmerged:1; /**< are equal quadratic variables already merged? */
172  unsigned int bilinsorted:1; /**< are the bilinear terms already sorted? */
173  unsigned int bilinmerged:1; /**< are equal bilinear terms (and bilinear terms with zero coefficient) already merged? */
174 
175  unsigned int isconvex:1; /**< is quadratic function is convex ? */
176  unsigned int isconcave:1; /**< is quadratic function is concave ? */
177  unsigned int iscurvchecked:1; /**< is quadratic function checked on convexity or concavity ? */
178  unsigned int isremovedfixings:1; /**< did we removed fixed/aggr/multiaggr variables ? */
179  unsigned int ispropagated:1; /**< was the constraint propagated with respect to the current bounds ? */
180  unsigned int ispresolved:1; /**< did we checked for possibilities of upgrading or implicit integer variables ? */
181  unsigned int initialmerge:1; /**< did we perform an initial merge and clean in presolving yet ? */
182 #ifdef CHECKIMPLINBILINEAR
183  unsigned int isimpladded:1; /**< has there been an implication added for a binary variable in a bilinear term? */
184 #endif
185  unsigned int isgaugeavailable:1; /**< is the gauge function computed? */
186  unsigned int isedavailable:1; /**< is the eigen decomposition of A available? */
187 
188  SCIP_Real minlinactivity; /**< sum of minimal activities of all linear terms with finite minimal activity */
189  SCIP_Real maxlinactivity; /**< sum of maximal activities of all linear terms with finite maximal activity */
190  int minlinactivityinf; /**< number of linear terms with infinite minimal activity */
191  int maxlinactivityinf; /**< number of linear terms with infinity maximal activity */
192  SCIP_INTERVAL quadactivitybounds; /**< bounds on the activity of the quadratic term, if up to date, otherwise empty interval */
193  SCIP_Real activity; /**< activity of quadratic function w.r.t. current solution */
194  SCIP_Real lhsviol; /**< violation of lower bound by current solution (used temporarily inside constraint handler) */
195  SCIP_Real rhsviol; /**< violation of lower bound by current solution (used temporarily inside constraint handler) */
196 
197  int linvar_maydecrease; /**< index of a variable in linvars that may be decreased without making any other constraint infeasible, or -1 if none */
198  int linvar_mayincrease; /**< index of a variable in linvars that may be increased without making any other constraint infeasible, or -1 if none */
199 
200  SCIP_VAR** sepaquadvars; /**< variables corresponding to quadvarterms to use in separation, only available in solving stage */
201  int* sepabilinvar2pos; /**< position of second variable in bilinear terms to use in separation, only available in solving stage */
202  SCIP_Real lincoefsmin; /**< minimal absolute value of coefficients in linear part, only available in solving stage */
203  SCIP_Real lincoefsmax; /**< maximal absolute value of coefficients in linear part, only available in solving stage */
204 
205  SCIP_Real* factorleft; /**< coefficients of left factor if constraint function is factorable */
206  SCIP_Real* factorright; /**< coefficients of right factor if constraint function is factorable */
207 
208  SCIP_Real* gaugecoefs; /**< coefficients of the gauge function */
209  SCIP_Real gaugeconst; /**< constant of the gauge function */
210  SCIP_Real* interiorpoint; /**< interior point of the region defined by the convex function */
211  SCIP_Real interiorpointval; /**< function value at interior point */
212 
213  SCIP_Real* eigenvalues; /**< eigenvalues of A */
214  SCIP_Real* eigenvectors; /**< orthonormal eigenvectors of A; if A = P D P^T, then eigenvectors is P^T */
215  SCIP_Real* bp; /**< stores b * P where b are the linear coefficients of the quadratic vars */
216  SCIP_Real maxnonconvexity; /**< nonconvexity measure: estimate on largest absolute value of nonconvex eigenvalues */
217 
218  SCIP_Bool isdisaggregated; /**< has the constraint already been disaggregated? if might happen that more disaggreation would be potentially
219  possible, but we reached the maximum number of sparsity components during presolveDisaggregate() */
220 };
221 
222 /** quadratic constraint update method */
223 struct SCIP_QuadConsUpgrade
224 {
225  SCIP_DECL_QUADCONSUPGD((*quadconsupgd)); /**< method to call for upgrading quadratic constraint */
226  int priority; /**< priority of upgrading method */
227  SCIP_Bool active; /**< is upgrading enabled */
228 };
229 typedef struct SCIP_QuadConsUpgrade SCIP_QUADCONSUPGRADE; /**< quadratic constraint update method */
230 
231 /** structure to store everything needed for using linear inequalities to improve upon the McCormick relaxation */
232 struct BilinearEstimator
233 {
234  SCIP_VAR* x; /**< first variable */
235  SCIP_VAR* y; /**< second variable */
236  SCIP_Real inequnderest[6]; /**< at most two inequalities that can be used to underestimate xy; stored as (xcoef,ycoef,constant) with xcoef x <= ycoef y + constant */
237  SCIP_Real ineqoverest[6]; /**< at most two inequalities that can be used to overestimate xy; stored as (xcoef,ycoef,constant) with xcoef x <= ycoef y + constant */
238  SCIP_Real maxnonconvexity; /**< estimate on largest absolute value of nonconvex eigenvalues of all quadratic constraint containing xy */
239  int ninequnderest; /**< total number of inequalities for underestimating xy */
240  int nineqoverest; /**< total number of inequalities for overestimating xy */
241  int nunderest; /**< number of constraints that require to underestimate xy */
242  int noverest; /**< number of constraints that require to overestimate xy */
243 
244  SCIP_Real lastimprfac; /**< last achieved improvement factor */
245 };
246 typedef struct BilinearEstimator BILINESTIMATOR;
247 
248 /** constraint handler data */
249 struct SCIP_ConshdlrData
250 {
251  int replacebinaryprodlength; /**< length of linear term which when multiplied with a binary variable is replaced by an auxiliary variable and an equivalent linear formulation */
252  int empathy4and; /**< how much empathy we have for using the AND constraint handler: 0 avoid always; 1 use sometimes; 2 use as often as possible */
253  SCIP_Bool binreforminitial; /**< whether to make constraints added due to replacing products with binary variables initial */
254  SCIP_Bool binreformbinaryonly;/**< whether to consider only binary variables when reformulating products with binary variables */
255  SCIP_Real binreformmaxcoef; /**< factor on 1/feastol to limit coefficients and coef range in linear constraints created by binary reformulation */
256  SCIP_Real cutmaxrange; /**< maximal range (maximal coef / minimal coef) of a cut in order to be added to LP */
257  SCIP_Bool linearizeheursol; /**< whether linearizations of convex quadratic constraints should be added to cutpool when some heuristics finds a new solution */
258  SCIP_Bool checkcurvature; /**< whether functions should be checked for convexity/concavity */
259  SCIP_Bool checkfactorable; /**< whether functions should be checked to be factorable */
260  char checkquadvarlocks; /**< whether quadratic variables contained in a single constraint should be forced to be at their lower or upper bounds ('d'isable, change 't'ype, add 'b'ound disjunction) */
261  SCIP_Bool linfeasshift; /**< whether to make solutions in check feasible if possible */
262  int maxdisaggrsize; /**< maximum number of components when disaggregating a quadratic constraint (<= 1: off) */
263  char disaggrmergemethod; /**< method on merging blocks in disaggregation */
264  int maxproprounds; /**< limit on number of propagation rounds for a single constraint within one round of SCIP propagation during solve */
265  int maxproproundspresolve; /**< limit on number of propagation rounds for a single constraint within one presolving round */
266  SCIP_Real sepanlpmincont; /**< minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation */
267  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
268  SCIP_Bool gaugecuts; /**< should convex quadratics generated strong cuts via gauge function? */
269  SCIP_Bool projectedcuts; /**< should convex quadratics generated strong cuts via projections? */
270  char interiorcomputation;/**< how the interior point should be computed: 'a'ny point per constraint, 'm'ost interior per constraint */
271  char branchscoring; /**< method to use to compute score of branching candidates */
272  int enfolplimit; /**< maximum number of enforcement round before declaring the LP relaxation
273  * infeasible (-1: no limit); WARNING: if this parameter is not set to -1,
274  * SCIP might declare sub-optimal solutions optimal or feasible instances
275  * infeasible; thus, the result returned by SCIP might be incorrect!
276  */
277  SCIP_HEUR* subnlpheur; /**< a pointer to the subnlp heuristic, if available */
278  SCIP_HEUR* trysolheur; /**< a pointer to the trysol heuristic, if available */
279  SCIP_EVENTHDLR* eventhdlr; /**< our handler for variable bound change events */
280  int newsoleventfilterpos; /**< filter position of new solution event handler, if caught */
281  SCIP_Bool sepanlp; /**< where linearization of the NLP relaxation solution added? */
282  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
283  int nenforounds; /**< counter on number of enforcement rounds for the current node */
284  SCIP_QUADCONSUPGRADE** quadconsupgrades; /**< quadratic constraint upgrade methods for specializing quadratic constraints */
285  int quadconsupgradessize; /**< size of quadconsupgrade array */
286  int nquadconsupgrades; /**< number of quadratic constraint upgrade methods */
287 
288  BILINESTIMATOR* bilinestimators; /**< array containing all required information for using stronger estimators for each bilinear term in all quadratic constraints */
289  int nbilinterms; /**< number of bilinear terms in all quadratic constraints */
290 
291  SCIP_Bool usebilinineqbranch; /**< should linear inequalities be considered when computing the branching scores for bilinear terms? */
292  SCIP_Bool storedbilinearterms; /**< did we already try to store all bilinear terms? */
293 
294  SCIP_Real minscorebilinterms; /**< minimal required score in order to use linear inequalities for tighter bilinear relaxations */
295  SCIP_Real mincurvcollectbilinterms;/**< minimal curvature of constraints to be considered when returning bilinear terms to other plugins */
296  int bilinineqmaxseparounds; /**< maximum number of separation rounds to use linear inequalities for the bilinear term relaxation in a local node */
297 };
298 
299 
300 /*
301  * local methods for managing quadratic constraint update methods
302  */
303 
304 
305 /** checks whether a quadratic constraint upgrade method has already be registered */
306 static
307 SCIP_Bool conshdlrdataHasUpgrade(
308  SCIP* scip, /**< SCIP data structure */
309  SCIP_CONSHDLRDATA* conshdlrdata, /**< constraint handler data */
310  SCIP_DECL_QUADCONSUPGD((*quadconsupgd)), /**< method to call for upgrading quadratic constraint */
311  const char* conshdlrname /**< name of the constraint handler */
312  )
313 {
314  int i;
315 
316  assert(scip != NULL);
317  assert(conshdlrdata != NULL);
318  assert(quadconsupgd != NULL);
319  assert(conshdlrname != NULL);
320 
321  for( i = conshdlrdata->nquadconsupgrades - 1; i >= 0; --i )
322  {
323  if( conshdlrdata->quadconsupgrades[i]->quadconsupgd == quadconsupgd )
324  {
325  SCIPwarningMessage(scip, "Try to add already known upgrade message for constraint handler <%s>.\n", conshdlrname);
326  return TRUE;
327  }
328  }
329 
330  return FALSE;
331 }
332 
333 /*
334  * Local methods
335  */
336 
337 /** translate from one value of infinity to another
338  *
339  * if val is >= infty1, then give infty2, else give val
340  */
341 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
342 
343 /** catches variable bound change events on a linear variable in a quadratic constraint */
344 static
345 SCIP_RETCODE catchLinearVarEvents(
346  SCIP* scip, /**< SCIP data structure */
347  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
348  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
349  int linvarpos /**< position of variable in linear variables array */
350  )
351 {
352  SCIP_CONSDATA* consdata;
353  SCIP_QUADVAREVENTDATA* eventdata;
354  SCIP_EVENTTYPE eventtype;
355 
356  assert(scip != NULL);
357  assert(eventhdlr != NULL);
358  assert(cons != NULL);
359 
360  consdata = SCIPconsGetData(cons);
361  assert(consdata != NULL);
362 
363  assert(linvarpos >= 0);
364  assert(linvarpos < consdata->nlinvars);
365  assert(consdata->lineventdata != NULL);
366 
367  SCIP_CALL( SCIPallocBlockMemory(scip, &eventdata) );
368 
369  eventdata->cons = cons;
370  eventdata->varidx = linvarpos;
371 
372  eventtype = SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
373  if( !SCIPisInfinity(scip, consdata->rhs) )
374  {
375  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest
376  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
377  if( consdata->lincoefs[linvarpos] > 0.0 )
378  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
379  else
380  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
381  }
382  if( !SCIPisInfinity(scip, -consdata->lhs) )
383  {
384  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest
385  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
386  if( consdata->lincoefs[linvarpos] > 0.0 )
387  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
388  else
389  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
390  }
391 
392  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->linvars[linvarpos], eventtype, eventhdlr, (SCIP_EVENTDATA*)eventdata, &eventdata->filterpos) );
393 
394  consdata->lineventdata[linvarpos] = eventdata;
395 
396  /* invalidate activity information
397  * NOTE: It could happen that a constraint gets temporary deactivated and some variable bounds change. In this case
398  * we do not recognize those bound changes with the variable events and thus we have to recompute the activities.
399  */
400  consdata->minlinactivity = SCIP_INVALID;
401  consdata->maxlinactivity = SCIP_INVALID;
402  consdata->minlinactivityinf = -1;
403  consdata->maxlinactivityinf = -1;
404 
405  return SCIP_OKAY;
406 }
407 
408 /** drops variable bound change events on a linear variable in a quadratic constraint */
409 static
410 SCIP_RETCODE dropLinearVarEvents(
411  SCIP* scip, /**< SCIP data structure */
412  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
413  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
414  int linvarpos /**< position of variable in linear variables array */
415  )
416 {
417  SCIP_CONSDATA* consdata;
418  SCIP_EVENTTYPE eventtype;
419 
420  assert(scip != NULL);
421  assert(eventhdlr != NULL);
422  assert(cons != NULL);
423 
424  consdata = SCIPconsGetData(cons);
425  assert(consdata != NULL);
426 
427  assert(linvarpos >= 0);
428  assert(linvarpos < consdata->nlinvars);
429  assert(consdata->lineventdata != NULL);
430  assert(consdata->lineventdata[linvarpos] != NULL);
431  assert(consdata->lineventdata[linvarpos]->cons == cons);
432  assert(consdata->lineventdata[linvarpos]->varidx == linvarpos);
433  assert(consdata->lineventdata[linvarpos]->filterpos >= 0);
434 
435  eventtype = SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
436  if( !SCIPisInfinity(scip, consdata->rhs) )
437  {
438  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest
439  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
440  if( consdata->lincoefs[linvarpos] > 0.0 )
441  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
442  else
443  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
444  }
445  if( !SCIPisInfinity(scip, -consdata->lhs) )
446  {
447  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest
448  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
449  if( consdata->lincoefs[linvarpos] > 0.0 )
450  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
451  else
452  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
453  }
454 
455  SCIP_CALL( SCIPdropVarEvent(scip, consdata->linvars[linvarpos], eventtype, eventhdlr, (SCIP_EVENTDATA*)consdata->lineventdata[linvarpos], consdata->lineventdata[linvarpos]->filterpos) );
456 
457  SCIPfreeBlockMemory(scip, &consdata->lineventdata[linvarpos]); /*lint !e866 */
458 
459  return SCIP_OKAY;
460 }
461 
462 /** catches variable bound change events on a quadratic variable in a quadratic constraint */
463 static
464 SCIP_RETCODE catchQuadVarEvents(
465  SCIP* scip, /**< SCIP data structure */
466  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
467  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
468  int quadvarpos /**< position of variable in quadratic variables array */
469  )
470 {
471  SCIP_CONSDATA* consdata;
472  SCIP_QUADVAREVENTDATA* eventdata;
473  SCIP_EVENTTYPE eventtype;
474 
475  assert(scip != NULL);
476  assert(eventhdlr != NULL);
477  assert(cons != NULL);
478 
479  consdata = SCIPconsGetData(cons);
480  assert(consdata != NULL);
481 
482  assert(quadvarpos >= 0);
483  assert(quadvarpos < consdata->nquadvars);
484  assert(consdata->quadvarterms[quadvarpos].eventdata == NULL);
485 
486  SCIP_CALL( SCIPallocBlockMemory(scip, &eventdata) );
487 
488  eventtype = SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
489 #ifdef CHECKIMPLINBILINEAR
490  eventtype |= SCIP_EVENTTYPE_IMPLADDED;
491 #endif
492  eventdata->cons = cons;
493  eventdata->varidx = -quadvarpos-1;
494  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->quadvarterms[quadvarpos].var, eventtype, eventhdlr, (SCIP_EVENTDATA*)eventdata, &eventdata->filterpos) );
495 
496  consdata->quadvarterms[quadvarpos].eventdata = eventdata;
497 
498  /* invalidate activity information
499  * NOTE: It could happen that a constraint gets temporary deactivated and some variable bounds change. In this case
500  * we do not recognize those bound changes with the variable events and thus we have to recompute the activities.
501  */
502  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
503 
504  return SCIP_OKAY;
505 }
506 
507 /** catches variable bound change events on a quadratic variable in a quadratic constraint */
508 static
509 SCIP_RETCODE dropQuadVarEvents(
510  SCIP* scip, /**< SCIP data structure */
511  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
512  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
513  int quadvarpos /**< position of variable in quadratic variables array */
514  )
515 {
516  SCIP_CONSDATA* consdata;
517  SCIP_EVENTTYPE eventtype;
518 
519  assert(scip != NULL);
520  assert(eventhdlr != NULL);
521  assert(cons != NULL);
522 
523  consdata = SCIPconsGetData(cons);
524  assert(consdata != NULL);
525 
526  assert(quadvarpos >= 0);
527  assert(quadvarpos < consdata->nquadvars);
528  assert(consdata->quadvarterms[quadvarpos].eventdata != NULL);
529  assert(consdata->quadvarterms[quadvarpos].eventdata->cons == cons);
530  assert(consdata->quadvarterms[quadvarpos].eventdata->varidx == -quadvarpos-1);
531  assert(consdata->quadvarterms[quadvarpos].eventdata->filterpos >= 0);
532 
533  eventtype = SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
534 #ifdef CHECKIMPLINBILINEAR
535  eventtype |= SCIP_EVENTTYPE_IMPLADDED;
536 #endif
537 
538  SCIP_CALL( SCIPdropVarEvent(scip, consdata->quadvarterms[quadvarpos].var, eventtype, eventhdlr, (SCIP_EVENTDATA*)consdata->quadvarterms[quadvarpos].eventdata, consdata->quadvarterms[quadvarpos].eventdata->filterpos) );
539 
540  SCIPfreeBlockMemory(scip, &consdata->quadvarterms[quadvarpos].eventdata);
541 
542  return SCIP_OKAY;
543 }
544 
545 /** catch variable events */
546 static
547 SCIP_RETCODE catchVarEvents(
548  SCIP* scip, /**< SCIP data structure */
549  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
550  SCIP_CONS* cons /**< constraint for which to catch bound change events */
551  )
552 {
553  SCIP_CONSDATA* consdata;
554  SCIP_VAR* var;
555  int i;
556 
557  assert(scip != NULL);
558  assert(cons != NULL);
559  assert(eventhdlr != NULL);
560 
561  consdata = SCIPconsGetData(cons);
562  assert(consdata != NULL);
563  assert(consdata->lineventdata == NULL);
564 
565  /* we will update isremovedfixings, so reset it to TRUE first */
566  consdata->isremovedfixings = TRUE;
567 
568  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize) );
569  for( i = 0; i < consdata->nlinvars; ++i )
570  {
571  SCIP_CALL( catchLinearVarEvents(scip, eventhdlr, cons, i) );
572 
573  var = consdata->linvars[i];
574  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
575  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
576  }
577 
578  for( i = 0; i < consdata->nquadvars; ++i )
579  {
580  assert(consdata->quadvarterms[i].eventdata == NULL);
581 
582  SCIP_CALL( catchQuadVarEvents(scip, eventhdlr, cons, i) );
583 
584  var = consdata->quadvarterms[i].var;
585  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
586  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
587  }
588 
589  consdata->ispropagated = FALSE;
590 
591  return SCIP_OKAY;
592 }
593 
594 /** drop variable events */
595 static
596 SCIP_RETCODE dropVarEvents(
597  SCIP* scip, /**< SCIP data structure */
598  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
599  SCIP_CONS* cons /**< constraint for which to drop bound change events */
600  )
601 {
602  SCIP_CONSDATA* consdata;
603  int i;
604 
605  assert(scip != NULL);
606  assert(eventhdlr != NULL);
607  assert(cons != NULL);
608 
609  consdata = SCIPconsGetData(cons);
610  assert(consdata != NULL);
611 
612  if( consdata->lineventdata != NULL )
613  {
614  for( i = 0; i < consdata->nlinvars; ++i )
615  {
616  if( consdata->lineventdata[i] != NULL )
617  {
618  SCIP_CALL( dropLinearVarEvents(scip, eventhdlr, cons, i) );
619  }
620  }
621  SCIPfreeBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize);
622  }
623 
624  for( i = 0; i < consdata->nquadvars; ++i )
625  {
626  if( consdata->quadvarterms[i].eventdata != NULL )
627  {
628  SCIP_CALL( dropQuadVarEvents(scip, eventhdlr, cons, i) );
629  }
630  }
631 
632  return SCIP_OKAY;
633 }
634 
635 /** locks a linear variable in a constraint */
636 static
637 SCIP_RETCODE lockLinearVariable(
638  SCIP* scip, /**< SCIP data structure */
639  SCIP_CONS* cons, /**< constraint where to lock a variable */
640  SCIP_VAR* var, /**< variable to lock */
641  SCIP_Real coef /**< coefficient of variable in constraint */
642  )
643 {
644  SCIP_CONSDATA* consdata;
645 
646  assert(scip != NULL);
647  assert(cons != NULL);
648  assert(var != NULL);
649  assert(coef != 0.0);
650 
651  consdata = SCIPconsGetData(cons);
652  assert(consdata != NULL);
653 
654  if( coef > 0.0 )
655  {
656  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
657  }
658  else
659  {
660  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
661  }
662 
663  return SCIP_OKAY;
664 }
665 
666 /** unlocks a linear variable in a constraint */
667 static
668 SCIP_RETCODE unlockLinearVariable(
669  SCIP* scip, /**< SCIP data structure */
670  SCIP_CONS* cons, /**< constraint where to unlock a variable */
671  SCIP_VAR* var, /**< variable to unlock */
672  SCIP_Real coef /**< coefficient of variable in constraint */
673  )
674 {
675  SCIP_CONSDATA* consdata;
676 
677  assert(scip != NULL);
678  assert(cons != NULL);
679  assert(var != NULL);
680  assert(coef != 0.0);
681 
682  consdata = SCIPconsGetData(cons);
683  assert(consdata != NULL);
684 
685  if( coef > 0.0 )
686  {
687  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
688  }
689  else
690  {
691  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
692  }
693 
694  return SCIP_OKAY;
695 }
696 
697 /** locks a quadratic variable in a constraint */
698 static
699 SCIP_RETCODE lockQuadraticVariable(
700  SCIP* scip, /**< SCIP data structure */
701  SCIP_CONS* cons, /**< constraint where to lock a variable */
702  SCIP_VAR* var /**< variable to lock */
703  )
704 {
705  SCIP_CALL( SCIPlockVarCons(scip, var, cons, TRUE, TRUE) );
706 
707  return SCIP_OKAY;
708 }
709 
710 /** unlocks a quadratic variable in a constraint */
711 static
712 SCIP_RETCODE unlockQuadraticVariable(
713  SCIP* scip, /**< SCIP data structure */
714  SCIP_CONS* cons, /**< constraint where to unlock a variable */
715  SCIP_VAR* var /**< variable to unlock */
716  )
717 {
718  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, TRUE, TRUE) );
719 
720  return SCIP_OKAY;
721 }
722 
723 /** computes the minimal and maximal activity for the linear part in a constraint data
724  *
725  * Only sums up terms that contribute finite values.
726  * Gives the number of terms that contribute infinite values.
727  * Only computes those activities where the corresponding side of the constraint is finite.
728  */
729 static
730 void consdataUpdateLinearActivity(
731  SCIP* scip, /**< SCIP data structure */
732  SCIP_CONSDATA* consdata, /**< constraint data */
733  SCIP_Real intervalinfty /**< infinity value used in interval operations */
734  )
735 { /*lint --e{666}*/
736  SCIP_ROUNDMODE prevroundmode;
737  int i;
738  SCIP_Real bnd;
739 
740  assert(scip != NULL);
741  assert(consdata != NULL);
742 
743  /* if variable bounds are not strictly consistent, then the activity update methods may yield inconsistent activities
744  * in this case, we also recompute the activities
745  */
746  if( consdata->minlinactivity != SCIP_INVALID && consdata->maxlinactivity != SCIP_INVALID && /*lint !e777 */
747  (consdata->minlinactivityinf > 0 || consdata->maxlinactivityinf > 0 || consdata->minlinactivity <= consdata->maxlinactivity) )
748  {
749  /* activities should be up-to-date */
750  assert(consdata->minlinactivityinf >= 0);
751  assert(consdata->maxlinactivityinf >= 0);
752  return;
753  }
754 
755  consdata->minlinactivityinf = 0;
756  consdata->maxlinactivityinf = 0;
757 
758  /* if lhs is -infinite, then we do not compute a maximal activity, so we set it to infinity
759  * if rhs is infinite, then we do not compute a minimal activity, so we set it to -infinity
760  */
761  consdata->minlinactivity = SCIPisInfinity(scip, consdata->rhs) ? -intervalinfty : 0.0;
762  consdata->maxlinactivity = SCIPisInfinity(scip, -consdata->lhs) ? intervalinfty : 0.0;
763 
764  if( consdata->nlinvars == 0 )
765  return;
766 
767  /* if the activities computed here should be still up-to-date after bound changes,
768  * variable events need to be caught */
769  assert(consdata->lineventdata != NULL);
770 
771  prevroundmode = SCIPintervalGetRoundingMode();
772 
773  if( !SCIPisInfinity(scip, consdata->rhs) )
774  {
775  /* compute minimal activity only if there is a finite right hand side */
776  SCIPintervalSetRoundingModeDownwards();
777 
778  for( i = 0; i < consdata->nlinvars; ++i )
779  {
780  assert(consdata->lineventdata[i] != NULL);
781  if( consdata->lincoefs[i] >= 0.0 )
782  {
783  bnd = MIN(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
784  if( SCIPisInfinity(scip, -bnd) )
785  {
786  ++consdata->minlinactivityinf;
787  continue;
788  }
789  assert(!SCIPisInfinity(scip, bnd)); /* do not like variables that are fixed at +infinity */
790  }
791  else
792  {
793  bnd = MAX(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
794  if( SCIPisInfinity(scip, bnd) )
795  {
796  ++consdata->minlinactivityinf;
797  continue;
798  }
799  assert(!SCIPisInfinity(scip, -bnd)); /* do not like variables that are fixed at -infinity */
800  }
801  consdata->minlinactivity += consdata->lincoefs[i] * bnd;
802  }
803  }
804 
805  if( !SCIPisInfinity(scip, -consdata->lhs) )
806  {
807  /* compute maximal activity only if there is a finite left hand side */
808  SCIPintervalSetRoundingModeUpwards();
809 
810  for( i = 0; i < consdata->nlinvars; ++i )
811  {
812  assert(consdata->lineventdata[i] != NULL);
813  if( consdata->lincoefs[i] >= 0.0 )
814  {
815  bnd = MAX(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
816  if( SCIPisInfinity(scip, bnd) )
817  {
818  ++consdata->maxlinactivityinf;
819  continue;
820  }
821  assert(!SCIPisInfinity(scip, -bnd)); /* do not like variables that are fixed at -infinity */
822  }
823  else
824  {
825  bnd = MIN(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
826  if( SCIPisInfinity(scip, -bnd) )
827  {
828  ++consdata->maxlinactivityinf;
829  continue;
830  }
831  assert(!SCIPisInfinity(scip, bnd)); /* do not like variables that are fixed at +infinity */
832  }
833  consdata->maxlinactivity += consdata->lincoefs[i] * bnd;
834  }
835  }
836 
837  SCIPintervalSetRoundingMode(prevroundmode);
838 
839  assert(consdata->minlinactivityinf > 0 || consdata->maxlinactivityinf > 0 || consdata->minlinactivity <= consdata->maxlinactivity);
840 }
841 
842 /** update the linear activities after a change in the lower bound of a variable */
843 static
844 void consdataUpdateLinearActivityLbChange(
845  SCIP* scip, /**< SCIP data structure */
846  SCIP_CONSDATA* consdata, /**< constraint data */
847  SCIP_Real coef, /**< coefficient of variable in constraint */
848  SCIP_Real oldbnd, /**< previous lower bound of variable */
849  SCIP_Real newbnd /**< new lower bound of variable */
850  )
851 {
852  SCIP_ROUNDMODE prevroundmode;
853 
854  assert(scip != NULL);
855  assert(consdata != NULL);
856  /* we can't deal with lower bounds at infinity */
857  assert(!SCIPisInfinity(scip, oldbnd));
858  assert(!SCIPisInfinity(scip, newbnd));
859 
860  /* @todo since we check the linear activity for consistency later anyway, we may skip changing the rounding mode here */
861 
862  /* assume lhs <= a*x + y <= rhs, then the following bound changes can be deduced:
863  * a > 0: y <= rhs - a*lb(x), y >= lhs - a*ub(x)
864  * a < 0: y <= rhs - a*ub(x), y >= lhs - a*lb(x)
865  */
866 
867  if( coef > 0.0 )
868  {
869  /* we should only be called if rhs is finite */
870  assert(!SCIPisInfinity(scip, consdata->rhs));
871 
872  /* we have no min activities computed so far, so cannot update */
873  if( consdata->minlinactivity == SCIP_INVALID ) /*lint !e777 */
874  return;
875 
876  prevroundmode = SCIPintervalGetRoundingMode();
877  SCIPintervalSetRoundingModeDownwards();
878 
879  /* update min activity */
880  if( SCIPisInfinity(scip, -oldbnd) )
881  {
882  --consdata->minlinactivityinf;
883  assert(consdata->minlinactivityinf >= 0);
884  }
885  else
886  {
887  SCIP_Real minuscoef;
888  minuscoef = -coef;
889  consdata->minlinactivity += minuscoef * oldbnd;
890  }
891 
892  if( SCIPisInfinity(scip, -newbnd) )
893  {
894  ++consdata->minlinactivityinf;
895  }
896  else
897  {
898  consdata->minlinactivity += coef * newbnd;
899  }
900 
901  SCIPintervalSetRoundingMode(prevroundmode);
902  }
903  else
904  {
905  /* we should only be called if lhs is finite */
906  assert(!SCIPisInfinity(scip, -consdata->lhs));
907 
908  /* we have no max activities computed so far, so cannot update */
909  if( consdata->maxlinactivity == SCIP_INVALID ) /*lint !e777 */
910  return;
911 
912  prevroundmode = SCIPintervalGetRoundingMode();
913  SCIPintervalSetRoundingModeUpwards();
914 
915  /* update max activity */
916  if( SCIPisInfinity(scip, -oldbnd) )
917  {
918  --consdata->maxlinactivityinf;
919  assert(consdata->maxlinactivityinf >= 0);
920  }
921  else
922  {
923  SCIP_Real minuscoef;
924  minuscoef = -coef;
925  consdata->maxlinactivity += minuscoef * oldbnd;
926  }
927 
928  if( SCIPisInfinity(scip, -newbnd) )
929  {
930  ++consdata->maxlinactivityinf;
931  }
932  else
933  {
934  consdata->maxlinactivity += coef * newbnd;
935  }
936 
937  SCIPintervalSetRoundingMode(prevroundmode);
938  }
939 }
940 
941 /** update the linear activities after a change in the upper bound of a variable */
942 static
943 void consdataUpdateLinearActivityUbChange(
944  SCIP* scip, /**< SCIP data structure */
945  SCIP_CONSDATA* consdata, /**< constraint data */
946  SCIP_Real coef, /**< coefficient of variable in constraint */
947  SCIP_Real oldbnd, /**< previous lower bound of variable */
948  SCIP_Real newbnd /**< new lower bound of variable */
949  )
950 {
951  SCIP_ROUNDMODE prevroundmode;
952 
953  assert(scip != NULL);
954  assert(consdata != NULL);
955  /* we can't deal with upper bounds at -infinity */
956  assert(!SCIPisInfinity(scip, -oldbnd));
957  assert(!SCIPisInfinity(scip, -newbnd));
958 
959  /* @todo since we check the linear activity for consistency later anyway, we may skip changing the rounding mode here */
960 
961  /* assume lhs <= a*x + y <= rhs, then the following bound changes can be deduced:
962  * a > 0: y <= rhs - a*lb(x), y >= lhs - a*ub(x)
963  * a < 0: y <= rhs - a*ub(x), y >= lhs - a*lb(x)
964  */
965 
966  if( coef > 0.0 )
967  {
968  /* we should only be called if lhs is finite */
969  assert(!SCIPisInfinity(scip, -consdata->lhs));
970 
971  /* we have no max activities computed so far, so cannot update */
972  if( consdata->maxlinactivity == SCIP_INVALID ) /*lint !e777 */
973  return;
974 
975  prevroundmode = SCIPintervalGetRoundingMode();
976  SCIPintervalSetRoundingModeUpwards();
977 
978  /* update max activity */
979  if( SCIPisInfinity(scip, oldbnd) )
980  {
981  --consdata->maxlinactivityinf;
982  assert(consdata->maxlinactivityinf >= 0);
983  }
984  else
985  {
986  SCIP_Real minuscoef;
987  minuscoef = -coef;
988  consdata->maxlinactivity += minuscoef * oldbnd;
989  }
990 
991  if( SCIPisInfinity(scip, newbnd) )
992  {
993  ++consdata->maxlinactivityinf;
994  }
995  else
996  {
997  consdata->maxlinactivity += coef * newbnd;
998  }
999 
1000  SCIPintervalSetRoundingMode(prevroundmode);
1001  }
1002  else
1003  {
1004  /* we should only be called if rhs is finite */
1005  assert(!SCIPisInfinity(scip, consdata->rhs));
1006 
1007  /* we have no min activities computed so far, so cannot update */
1008  if( consdata->minlinactivity == SCIP_INVALID ) /*lint !e777 */
1009  return;
1010 
1011  prevroundmode = SCIPintervalGetRoundingMode();
1012  SCIPintervalSetRoundingModeDownwards();
1013 
1014  /* update min activity */
1015  if( SCIPisInfinity(scip, oldbnd) )
1016  {
1017  --consdata->minlinactivityinf;
1018  assert(consdata->minlinactivityinf >= 0);
1019  }
1020  else
1021  {
1022  SCIP_Real minuscoef;
1023  minuscoef = -coef;
1024  consdata->minlinactivity += minuscoef * oldbnd;
1025  }
1026 
1027  if( SCIPisInfinity(scip, newbnd) )
1028  {
1029  ++consdata->minlinactivityinf;
1030  }
1031  else
1032  {
1033  consdata->minlinactivity += coef * newbnd;
1034  }
1035 
1036  SCIPintervalSetRoundingMode(prevroundmode);
1037  }
1038 }
1039 
1040 /** returns whether a quadratic variable domain can be reduced to its lower or upper bound; this is the case if the
1041  * quadratic variable is in just one single quadratic constraint and (sqrcoef > 0 and LHS = -infinity), or
1042  * (sqrcoef < 0 and RHS = +infinity) hold
1043  */
1044 static
1045 SCIP_Bool hasQuadvarHpProperty(
1046  SCIP* scip, /**< SCIP data structure */
1047  SCIP_CONSDATA* consdata, /**< constraint data */
1048  int idx /**< index of quadratic variable */
1049  )
1050 {
1051  SCIP_VAR* var;
1052  SCIP_Real quadcoef;
1053  SCIP_Bool haslhs;
1054  SCIP_Bool hasrhs;
1055 
1056  assert(scip != NULL);
1057  assert(consdata != NULL);
1058  assert(idx >= 0 && idx < consdata->nquadvars);
1059 
1060  var = consdata->quadvarterms[idx].var;
1061  assert(var != NULL);
1062 
1063  quadcoef = consdata->quadvarterms[idx].sqrcoef;
1064  haslhs = !SCIPisInfinity(scip, -consdata->lhs);
1065  hasrhs = !SCIPisInfinity(scip, consdata->rhs);
1066 
1067  return SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 1
1068  && SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 1 && SCIPisZero(scip, SCIPvarGetObj(var))
1069  && SCIPvarGetType(var) != SCIP_VARTYPE_BINARY && ((quadcoef < 0.0 && !haslhs) || (quadcoef > 0.0 && !hasrhs));
1070 }
1071 
1072 /** processes variable fixing or bound change event */
1073 static
1074 SCIP_DECL_EVENTEXEC(processVarEvent)
1075 {
1076  SCIP_CONS* cons;
1077  SCIP_CONSDATA* consdata;
1078  SCIP_EVENTTYPE eventtype;
1079  int varidx;
1080 
1081  assert(scip != NULL);
1082  assert(event != NULL);
1083  assert(eventdata != NULL);
1084  assert(eventhdlr != NULL);
1085 
1086  cons = ((SCIP_QUADVAREVENTDATA*)eventdata)->cons;
1087  assert(cons != NULL);
1088  consdata = SCIPconsGetData(cons);
1089  assert(consdata != NULL);
1090 
1091  varidx = ((SCIP_QUADVAREVENTDATA*)eventdata)->varidx;
1092  assert(varidx < 0 || varidx < consdata->nlinvars);
1093  assert(varidx >= 0 || -varidx-1 < consdata->nquadvars);
1094 
1095  eventtype = SCIPeventGetType(event);
1096 
1097  /* process local bound changes */
1098  if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
1099  {
1100  if( varidx < 0 )
1101  {
1102  /* mark activity bounds for quad term as not up to date anymore */
1103  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
1104  }
1105  else
1106  {
1107  /* update activity bounds for linear terms */
1108  if( eventtype & SCIP_EVENTTYPE_LBCHANGED )
1109  consdataUpdateLinearActivityLbChange(scip, consdata, consdata->lincoefs[varidx], SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
1110  else
1111  consdataUpdateLinearActivityUbChange(scip, consdata, consdata->lincoefs[varidx], SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
1112  }
1113 
1114  if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
1115  {
1116  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
1117  consdata->ispropagated = FALSE;
1118  }
1119  }
1120 
1121  /* process global bound changes */
1122  if( eventtype & SCIP_EVENTTYPE_GBDCHANGED )
1123  {
1124  SCIP_VAR* var;
1125 
1126  var = varidx < 0 ? consdata->quadvarterms[-varidx-1].var : consdata->linvars[varidx];
1127  assert(var != NULL);
1128 
1129  if( varidx < 0 )
1130  {
1131  SCIP_QUADVARTERM* quadvarterm;
1132 
1133  quadvarterm = &consdata->quadvarterms[-varidx-1];
1134 
1135  /* if an integer variable x with a x^2 is tightened to [0,1], then we can replace the x^2 by x, which is done in mergeAndCleanQuadVarTerms()
1136  * we currently do this only if the binary variable does not show up in any bilinear terms
1137  */
1138  if( SCIPgetStage(scip) < SCIP_STAGE_SOLVING && SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER && SCIPvarIsBinary(var) &&
1139  quadvarterm->sqrcoef != 0.0 && quadvarterm->nadjbilin == 0 )
1140  {
1141  consdata->quadvarsmerged = FALSE;
1142  consdata->initialmerge = FALSE;
1143  }
1144  }
1145 
1146  if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
1147  consdata->isremovedfixings = FALSE;
1148  }
1149 
1150  /* process variable fixing event */
1151  if( eventtype & SCIP_EVENTTYPE_VARFIXED )
1152  {
1153  consdata->isremovedfixings = FALSE;
1154  }
1155 
1156 #ifdef CHECKIMPLINBILINEAR
1157  if( eventtype & SCIP_EVENTTYPE_IMPLADDED )
1158  {
1159  assert(varidx < 0); /* we catch impladded events only for quadratic variables */
1160  /* if variable is binary (quite likely if an implication has been added) and occurs in a bilinear term, then mark that we should check implications */
1161  if( SCIPvarIsBinary(SCIPeventGetVar(event)) && consdata->quadvarterms[-varidx-1].nadjbilin > 0 )
1162  consdata->isimpladded = TRUE;
1163  }
1164 #endif
1165 
1166  return SCIP_OKAY;
1167 }
1168 
1169 /** ensures, that linear vars and coefs arrays can store at least num entries */
1170 static
1171 SCIP_RETCODE consdataEnsureLinearVarsSize(
1172  SCIP* scip, /**< SCIP data structure */
1173  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1174  int num /**< minimum number of entries to store */
1175  )
1176 {
1177  assert(scip != NULL);
1178  assert(consdata != NULL);
1179  assert(consdata->nlinvars <= consdata->linvarssize);
1180 
1181  if( num > consdata->linvarssize )
1182  {
1183  int newsize;
1184 
1185  newsize = SCIPcalcMemGrowSize(scip, num);
1186  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->linvars, consdata->linvarssize, newsize) );
1187  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->lincoefs, consdata->linvarssize, newsize) );
1188  if( consdata->lineventdata != NULL )
1189  {
1190  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize, newsize) );
1191  }
1192  consdata->linvarssize = newsize;
1193  }
1194  assert(num <= consdata->linvarssize);
1195 
1196  return SCIP_OKAY;
1197 }
1198 
1199 /** ensures, that quadratic variable terms array can store at least num entries */
1200 static
1201 SCIP_RETCODE consdataEnsureQuadVarTermsSize(
1202  SCIP* scip, /**< SCIP data structure */
1203  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1204  int num /**< minimum number of entries to store */
1205  )
1206 {
1207  assert(scip != NULL);
1208  assert(consdata != NULL);
1209  assert(consdata->nquadvars <= consdata->quadvarssize);
1210 
1211  if( num > consdata->quadvarssize )
1212  {
1213  int newsize;
1214 
1215  newsize = SCIPcalcMemGrowSize(scip, num);
1216  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->quadvarterms, consdata->quadvarssize, newsize) );
1217  consdata->quadvarssize = newsize;
1218  }
1219  assert(num <= consdata->quadvarssize);
1220 
1221  return SCIP_OKAY;
1222 }
1223 
1224 /** ensures, that adjacency array can store at least num entries */
1225 static
1226 SCIP_RETCODE consdataEnsureAdjBilinSize(
1227  SCIP* scip, /**< SCIP data structure */
1228  SCIP_QUADVARTERM* quadvarterm, /**< quadratic variable term */
1229  int num /**< minimum number of entries to store */
1230  )
1231 {
1232  assert(scip != NULL);
1233  assert(quadvarterm != NULL);
1234  assert(quadvarterm->nadjbilin <= quadvarterm->adjbilinsize);
1235 
1236  if( num > quadvarterm->adjbilinsize )
1237  {
1238  int newsize;
1239 
1240  newsize = SCIPcalcMemGrowSize(scip, num);
1241  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &quadvarterm->adjbilin, quadvarterm->adjbilinsize, newsize) );
1242  quadvarterm->adjbilinsize = newsize;
1243  }
1244  assert(num <= quadvarterm->adjbilinsize);
1245 
1246  return SCIP_OKAY;
1247 }
1248 
1249 /** ensures, that bilinear term arrays can store at least num entries */
1250 static
1251 SCIP_RETCODE consdataEnsureBilinSize(
1252  SCIP* scip, /**< SCIP data structure */
1253  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1254  int num /**< minimum number of entries to store */
1255  )
1256 {
1257  assert(scip != NULL);
1258  assert(consdata != NULL);
1259  assert(consdata->nbilinterms <= consdata->bilintermssize);
1260 
1261  if( num > consdata->bilintermssize )
1262  {
1263  int newsize;
1264 
1265  newsize = SCIPcalcMemGrowSize(scip, num);
1266  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->bilinterms, consdata->bilintermssize, newsize) );
1267  consdata->bilintermssize = newsize;
1268  }
1269  assert(num <= consdata->bilintermssize);
1270 
1271  return SCIP_OKAY;
1272 }
1273 
1274 /** creates empty constraint data structure */
1275 static
1276 SCIP_RETCODE consdataCreateEmpty(
1277  SCIP* scip, /**< SCIP data structure */
1278  SCIP_CONSDATA** consdata /**< a buffer to store pointer to new constraint data */
1279  )
1280 {
1281  assert(scip != NULL);
1282  assert(consdata != NULL);
1283 
1284  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
1285  BMSclearMemory(*consdata);
1286 
1287  (*consdata)->lhs = -SCIPinfinity(scip);
1288  (*consdata)->rhs = SCIPinfinity(scip);
1289 
1290  (*consdata)->linvarssorted = TRUE;
1291  (*consdata)->linvarsmerged = TRUE;
1292  (*consdata)->quadvarssorted = TRUE;
1293  (*consdata)->quadvarsmerged = TRUE;
1294  (*consdata)->bilinsorted = TRUE;
1295  (*consdata)->bilinmerged = TRUE;
1296 
1297  (*consdata)->isremovedfixings = TRUE;
1298  (*consdata)->ispropagated = TRUE;
1299  (*consdata)->initialmerge = FALSE;
1300 
1301  (*consdata)->linvar_maydecrease = -1;
1302  (*consdata)->linvar_mayincrease = -1;
1303 
1304  (*consdata)->minlinactivity = SCIP_INVALID;
1305  (*consdata)->maxlinactivity = SCIP_INVALID;
1306  (*consdata)->minlinactivityinf = -1;
1307  (*consdata)->maxlinactivityinf = -1;
1308 
1309  (*consdata)->isgaugeavailable = FALSE;
1310  (*consdata)->isedavailable = FALSE;
1311 
1312  return SCIP_OKAY;
1313 }
1314 
1315 /** creates constraint data structure */
1316 static
1317 SCIP_RETCODE consdataCreate(
1318  SCIP* scip, /**< SCIP data structure */
1319  SCIP_CONSDATA** consdata, /**< a buffer to store pointer to new constraint data */
1320  SCIP_Real lhs, /**< left hand side of constraint */
1321  SCIP_Real rhs, /**< right hand side of constraint */
1322  int nlinvars, /**< number of linear variables */
1323  SCIP_VAR** linvars, /**< array of linear variables */
1324  SCIP_Real* lincoefs, /**< array of coefficients of linear variables */
1325  int nquadvars, /**< number of quadratic variables */
1326  SCIP_QUADVARTERM* quadvarterms, /**< array of quadratic variable terms */
1327  int nbilinterms, /**< number of bilinear terms */
1328  SCIP_BILINTERM* bilinterms, /**< array of bilinear terms */
1329  SCIP_Bool capturevars /**< whether we should capture variables */
1330  )
1331 {
1332  int i;
1333 
1334  assert(scip != NULL);
1335  assert(consdata != NULL);
1336 
1337  assert(nlinvars == 0 || linvars != NULL);
1338  assert(nlinvars == 0 || lincoefs != NULL);
1339  assert(nquadvars == 0 || quadvarterms != NULL);
1340  assert(nbilinterms == 0 || bilinterms != NULL);
1341 
1342  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
1343  BMSclearMemory(*consdata);
1344 
1345  (*consdata)->minlinactivity = SCIP_INVALID;
1346  (*consdata)->maxlinactivity = SCIP_INVALID;
1347  (*consdata)->minlinactivityinf = -1;
1348  (*consdata)->maxlinactivityinf = -1;
1349 
1350  (*consdata)->lhs = lhs;
1351  (*consdata)->rhs = rhs;
1352 
1353  if( nlinvars > 0 )
1354  {
1355  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->linvars, linvars, nlinvars) );
1356  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->lincoefs, lincoefs, nlinvars) );
1357  (*consdata)->nlinvars = nlinvars;
1358  (*consdata)->linvarssize = nlinvars;
1359 
1360  if( capturevars )
1361  for( i = 0; i < nlinvars; ++i )
1362  {
1363  SCIP_CALL( SCIPcaptureVar(scip, linvars[i]) );
1364  }
1365  }
1366  else
1367  {
1368  (*consdata)->linvarssorted = TRUE;
1369  (*consdata)->linvarsmerged = TRUE;
1370  (*consdata)->minlinactivity = 0.0;
1371  (*consdata)->maxlinactivity = 0.0;
1372  (*consdata)->minlinactivityinf = 0;
1373  (*consdata)->maxlinactivityinf = 0;
1374  }
1375 
1376  if( nquadvars > 0 )
1377  {
1378  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->quadvarterms, quadvarterms, nquadvars) );
1379 
1380  for( i = 0; i < nquadvars; ++i )
1381  {
1382  (*consdata)->quadvarterms[i].eventdata = NULL;
1383  if( quadvarterms[i].nadjbilin )
1384  {
1385  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->quadvarterms[i].adjbilin, quadvarterms[i].adjbilin, quadvarterms[i].nadjbilin) );
1386  (*consdata)->quadvarterms[i].adjbilinsize = quadvarterms[i].nadjbilin;
1387  }
1388  else
1389  {
1390  assert((*consdata)->quadvarterms[i].nadjbilin == 0);
1391  (*consdata)->quadvarterms[i].adjbilin = NULL;
1392  (*consdata)->quadvarterms[i].adjbilinsize = 0;
1393  }
1394  if( capturevars )
1395  {
1396  SCIP_CALL( SCIPcaptureVar(scip, quadvarterms[i].var) );
1397  }
1398  }
1399 
1400  (*consdata)->nquadvars = nquadvars;
1401  (*consdata)->quadvarssize = nquadvars;
1402  SCIPintervalSetEmpty(&(*consdata)->quadactivitybounds);
1403  }
1404  else
1405  {
1406  (*consdata)->quadvarssorted = TRUE;
1407  (*consdata)->quadvarsmerged = TRUE;
1408  SCIPintervalSet(&(*consdata)->quadactivitybounds, 0.0);
1409  }
1410 
1411  if( nbilinterms > 0 )
1412  {
1413  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->bilinterms, bilinterms, nbilinterms) );
1414  (*consdata)->nbilinterms = nbilinterms;
1415  (*consdata)->bilintermssize = nbilinterms;
1416  }
1417  else
1418  {
1419  (*consdata)->bilinsorted = TRUE;
1420  (*consdata)->bilinmerged = TRUE;
1421  }
1422 
1423  (*consdata)->linvar_maydecrease = -1;
1424  (*consdata)->linvar_mayincrease = -1;
1425 
1426  (*consdata)->activity = SCIP_INVALID;
1427  (*consdata)->lhsviol = SCIPisInfinity(scip, -lhs) ? 0.0 : SCIP_INVALID;
1428  (*consdata)->rhsviol = SCIPisInfinity(scip, rhs) ? 0.0 : SCIP_INVALID;
1429 
1430  (*consdata)->isgaugeavailable = FALSE;
1431 
1432  return SCIP_OKAY;
1433 }
1434 
1435 /** frees constraint data structure */
1436 static
1437 SCIP_RETCODE consdataFree(
1438  SCIP* scip, /**< SCIP data structure */
1439  SCIP_CONSDATA** consdata /**< pointer to constraint data to free */
1440  )
1441 {
1442  int i;
1443 
1444  assert(scip != NULL);
1445  assert(consdata != NULL);
1446  assert(*consdata != NULL);
1447 
1448  /* free sepa arrays, may exists if constraint is deleted in solving stage */
1449  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->sepaquadvars, (*consdata)->nquadvars);
1450  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->sepabilinvar2pos, (*consdata)->nbilinterms);
1451 
1452  /* release linear variables and free linear part */
1453  if( (*consdata)->linvarssize > 0 )
1454  {
1455  for( i = 0; i < (*consdata)->nlinvars; ++i )
1456  {
1457  assert((*consdata)->lineventdata == NULL || (*consdata)->lineventdata[i] == NULL); /* variable events should have been dropped earlier */
1458  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->linvars[i]) );
1459  }
1460  SCIPfreeBlockMemoryArray(scip, &(*consdata)->linvars, (*consdata)->linvarssize);
1461  SCIPfreeBlockMemoryArray(scip, &(*consdata)->lincoefs, (*consdata)->linvarssize);
1462  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->lineventdata, (*consdata)->linvarssize);
1463  }
1464  assert((*consdata)->linvars == NULL);
1465  assert((*consdata)->lincoefs == NULL);
1466  assert((*consdata)->lineventdata == NULL);
1467 
1468  /* release quadratic variables and free quadratic variable term part */
1469  for( i = 0; i < (*consdata)->nquadvars; ++i )
1470  {
1471  assert((*consdata)->quadvarterms[i].eventdata == NULL); /* variable events should have been dropped earlier */
1472  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->quadvarterms[i].adjbilin, (*consdata)->quadvarterms[i].adjbilinsize);
1473  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->quadvarterms[i].var) );
1474  }
1475  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->quadvarterms, (*consdata)->quadvarssize);
1476 
1477  /* free bilinear terms */
1478  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->bilinterms, (*consdata)->bilintermssize);
1479 
1480  /* free nonlinear row representation */
1481  if( (*consdata)->nlrow != NULL )
1482  {
1483  SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
1484  }
1485 
1486  /* free interior point information, may exists if constraint is deleted in solving stage */
1487  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->interiorpoint, (*consdata)->nquadvars);
1488  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->gaugecoefs, (*consdata)->nquadvars);
1489 
1490  /* free eigen decomposition information */
1491  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->eigenvalues, (*consdata)->nquadvars);
1492  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->eigenvectors, (int)((*consdata)->nquadvars*(*consdata)->nquadvars));
1493  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->bp, (*consdata)->nquadvars);
1494 
1495  /* free unique indices of bilinear terms array */
1496  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->bilintermsidx, (*consdata)->nbilinterms);
1497 
1498  SCIPfreeBlockMemory(scip, consdata);
1499  *consdata = NULL;
1500 
1501  return SCIP_OKAY;
1502 }
1503 
1504 /** sorts linear part of constraint data */
1505 static
1506 void consdataSortLinearVars(
1507  SCIP_CONSDATA* consdata /**< quadratic constraint data */
1508  )
1509 {
1510  assert(consdata != NULL);
1511 
1512  if( consdata->linvarssorted )
1513  return;
1514 
1515  if( consdata->nlinvars <= 1 )
1516  {
1517  consdata->linvarssorted = TRUE;
1518  return;
1519  }
1520 
1521  if( consdata->lineventdata == NULL )
1522  {
1523  SCIPsortPtrReal((void**)consdata->linvars, consdata->lincoefs, SCIPvarComp, consdata->nlinvars);
1524  }
1525  else
1526  {
1527  int i;
1528 
1529  SCIPsortPtrPtrReal((void**)consdata->linvars, (void**)consdata->lineventdata, consdata->lincoefs, SCIPvarComp, consdata->nlinvars);
1530 
1531  /* update variable indices in event data */
1532  for( i = 0; i < consdata->nlinvars; ++i )
1533  if( consdata->lineventdata[i] != NULL )
1534  consdata->lineventdata[i]->varidx = i;
1535  }
1536 
1537  consdata->linvarssorted = TRUE;
1538 }
1539 
1540 #ifdef SCIP_DISABLED_CODE /* no-one needs this routine currently */
1541 /** returns the position of variable in the linear coefficients array of a constraint, or -1 if not found */
1542 static
1543 int consdataFindLinearVar(
1544  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1545  SCIP_VAR* var /**< variable to search for */
1546  )
1547 {
1548  int pos;
1549 
1550  assert(consdata != NULL);
1551  assert(var != NULL);
1552 
1553  if( consdata->nlinvars == 0 )
1554  return -1;
1555 
1556  consdataSortLinearVars(consdata);
1557 
1558  if( !SCIPsortedvecFindPtr((void**)consdata->linvars, SCIPvarComp, (void*)var, consdata->nlinvars, &pos) )
1559  pos = -1;
1560 
1561  return pos;
1562 }
1563 #endif
1564 
1565 /** index comparison method for quadratic variable terms: compares two indices of the quadratic variable set in the quadratic constraint */
1566 static
1567 SCIP_DECL_SORTINDCOMP(quadVarTermComp)
1568 { /*lint --e{715}*/
1569  SCIP_CONSDATA* consdata = (SCIP_CONSDATA*)dataptr;
1570 
1571  assert(consdata != NULL);
1572  assert(0 <= ind1 && ind1 < consdata->nquadvars);
1573  assert(0 <= ind2 && ind2 < consdata->nquadvars);
1574 
1575  return SCIPvarCompare(consdata->quadvarterms[ind1].var, consdata->quadvarterms[ind2].var);
1576 }
1577 
1578 /** sorting of quadratic variable terms */
1579 static
1580 SCIP_RETCODE consdataSortQuadVarTerms(
1581  SCIP* scip, /**< SCIP data structure */
1582  SCIP_CONSDATA* consdata /**< quadratic constraint data */
1583  )
1584 {
1585  int* perm;
1586  int i;
1587  int nexti;
1588  int v;
1589  SCIP_QUADVARTERM quadterm;
1590 
1591  assert(scip != NULL);
1592  assert(consdata != NULL);
1593 
1594  if( consdata->quadvarssorted )
1595  return SCIP_OKAY;
1596 
1597  if( consdata->nquadvars == 0 )
1598  {
1599  consdata->quadvarssorted = TRUE;
1600  return SCIP_OKAY;
1601  }
1602 
1603  /* get temporary memory to store the sorted permutation */
1604  SCIP_CALL( SCIPallocBufferArray(scip, &perm, consdata->nquadvars) );
1605 
1606  /* call bubble sort */
1607  SCIPsort(perm, quadVarTermComp, (void*)consdata, consdata->nquadvars);
1608 
1609  /* permute the quadratic variable terms according to the resulting permutation */
1610  for( v = 0; v < consdata->nquadvars; ++v )
1611  {
1612  if( perm[v] != v )
1613  {
1614  quadterm = consdata->quadvarterms[v];
1615 
1616  i = v;
1617  do
1618  {
1619  assert(0 <= perm[i] && perm[i] < consdata->nquadvars);
1620  assert(perm[i] != i);
1621  consdata->quadvarterms[i] = consdata->quadvarterms[perm[i]];
1622  if( consdata->quadvarterms[i].eventdata != NULL )
1623  {
1624  consdata->quadvarterms[i].eventdata->varidx = -i-1;
1625  }
1626  nexti = perm[i];
1627  perm[i] = i;
1628  i = nexti;
1629  }
1630  while( perm[i] != v );
1631  consdata->quadvarterms[i] = quadterm;
1632  if( consdata->quadvarterms[i].eventdata != NULL )
1633  {
1634  consdata->quadvarterms[i].eventdata->varidx = -i-1;
1635  }
1636  perm[i] = i;
1637  }
1638  }
1639  consdata->quadvarssorted = TRUE;
1640 
1641  /* free temporary memory */
1642  SCIPfreeBufferArray(scip, &perm);
1643 
1644  return SCIP_OKAY;
1645 }
1646 
1647 /** returns the position of variable in the quadratic variable terms array of a constraint, or -1 if not found */
1648 static
1649 SCIP_RETCODE consdataFindQuadVarTerm(
1650  SCIP* scip, /**< SCIP data structure */
1651  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1652  SCIP_VAR* var, /**< variable to search for */
1653  int* pos /**< buffer where to store position of var in quadvarterms array, or -1 if not found */
1654  )
1655 {
1656  int left;
1657  int right;
1658  int cmpres;
1659 
1660  assert(consdata != NULL);
1661  assert(var != NULL);
1662  assert(pos != NULL);
1663 
1664  if( consdata->nquadvars == 0 )
1665  {
1666  *pos = -1;
1667  return SCIP_OKAY;
1668  }
1669 
1670  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
1671 
1672  left = 0;
1673  right = consdata->nquadvars - 1;
1674  while( left <= right )
1675  {
1676  int middle;
1677 
1678  middle = (left+right)/2;
1679  assert(0 <= middle && middle < consdata->nquadvars);
1680 
1681  cmpres = SCIPvarCompare(var, consdata->quadvarterms[middle].var);
1682 
1683  if( cmpres < 0 )
1684  right = middle - 1;
1685  else if( cmpres > 0 )
1686  left = middle + 1;
1687  else
1688  {
1689  *pos = middle;
1690  return SCIP_OKAY;
1691  }
1692  }
1693  assert(left == right+1);
1694 
1695  *pos = -1;
1696 
1697  return SCIP_OKAY;
1698 }
1699 
1700 /** index comparison method for bilinear terms: compares two index pairs of the bilinear term set in the quadratic constraint */
1701 static
1702 SCIP_DECL_SORTINDCOMP(bilinTermComp)
1703 { /*lint --e{715}*/
1704  SCIP_CONSDATA* consdata = (SCIP_CONSDATA*)dataptr;
1705  int var1cmp;
1706 
1707  assert(consdata != NULL);
1708  assert(0 <= ind1 && ind1 < consdata->nbilinterms);
1709  assert(0 <= ind2 && ind2 < consdata->nbilinterms);
1710 
1711  var1cmp = SCIPvarCompare(consdata->bilinterms[ind1].var1, consdata->bilinterms[ind2].var1);
1712  if( var1cmp != 0 )
1713  return var1cmp;
1714 
1715  return SCIPvarCompare(consdata->bilinterms[ind1].var2, consdata->bilinterms[ind2].var2);
1716 }
1717 
1718 #ifndef NDEBUG
1719 /** checks if all bilinear terms are sorted correctly */
1720 static
1721 SCIP_Bool consdataCheckBilinTermsSort(
1722  SCIP_CONSDATA* consdata
1723  )
1724 {
1725  int i;
1726 
1727  assert(consdata != NULL);
1728 
1729  /* nothing to check if the bilinear terms have not been sorted yet */
1730  if( !consdata->bilinsorted )
1731  return TRUE;
1732 
1733  for( i = 0; i < consdata->nbilinterms - 1; ++i )
1734  {
1735  if( bilinTermComp(consdata, i, i+1) > 0 )
1736  return FALSE;
1737  }
1738  return TRUE;
1739 }
1740 #endif
1741 
1742 /** sorting of bilinear terms */
1743 static
1744 SCIP_RETCODE consdataSortBilinTerms(
1745  SCIP* scip, /**< SCIP data structure */
1746  SCIP_CONSDATA* consdata /**< quadratic constraint data */
1747  )
1748 {
1749  int* perm;
1750  int* invperm;
1751  int i;
1752  int nexti;
1753  int v;
1754  SCIP_BILINTERM bilinterm;
1755 
1756  assert(scip != NULL);
1757  assert(consdata != NULL);
1758 
1759  if( consdata->bilinsorted )
1760  return SCIP_OKAY;
1761 
1762  if( consdata->nbilinterms == 0 )
1763  {
1764  consdata->bilinsorted = TRUE;
1765  return SCIP_OKAY;
1766  }
1767 
1768  /* get temporary memory to store the sorted permutation and the inverse permutation */
1769  SCIP_CALL( SCIPallocBufferArray(scip, &perm, consdata->nbilinterms) );
1770  SCIP_CALL( SCIPallocBufferArray(scip, &invperm, consdata->nbilinterms) );
1771 
1772  /* call bubble sort */
1773  SCIPsort(perm, bilinTermComp, (void*)consdata, consdata->nbilinterms);
1774 
1775  /* compute inverted permutation */
1776  for( v = 0; v < consdata->nbilinterms; ++v )
1777  {
1778  assert(0 <= perm[v] && perm[v] < consdata->nbilinterms);
1779  invperm[perm[v]] = v;
1780  }
1781 
1782  /* permute the bilinear terms according to the resulting permutation */
1783  for( v = 0; v < consdata->nbilinterms; ++v )
1784  {
1785  if( perm[v] != v )
1786  {
1787  bilinterm = consdata->bilinterms[v];
1788 
1789  i = v;
1790  do
1791  {
1792  assert(0 <= perm[i] && perm[i] < consdata->nbilinterms);
1793  assert(perm[i] != i);
1794  consdata->bilinterms[i] = consdata->bilinterms[perm[i]];
1795  nexti = perm[i];
1796  perm[i] = i;
1797  i = nexti;
1798  }
1799  while( perm[i] != v );
1800  consdata->bilinterms[i] = bilinterm;
1801  perm[i] = i;
1802  }
1803  }
1804 
1805  /* update the adjacency information in the quadratic variable terms */
1806  for( v = 0; v < consdata->nquadvars; ++v )
1807  for( i = 0; i < consdata->quadvarterms[v].nadjbilin; ++i )
1808  consdata->quadvarterms[v].adjbilin[i] = invperm[consdata->quadvarterms[v].adjbilin[i]];
1809 
1810  consdata->bilinsorted = TRUE;
1811  assert(consdataCheckBilinTermsSort(consdata));
1812 
1813  /* free temporary memory */
1814  SCIPfreeBufferArray(scip, &invperm);
1815  SCIPfreeBufferArray(scip, &perm);
1816 
1817  return SCIP_OKAY;
1818 }
1819 
1820 /** moves a linear variable from one position to another */
1821 static
1822 void consdataMoveLinearVar(
1823  SCIP_CONSDATA* consdata, /**< constraint data */
1824  int oldpos, /**< position of variable that shall be moved */
1825  int newpos /**< new position of variable */
1826  )
1827 {
1828  assert(consdata != NULL);
1829  assert(oldpos >= 0);
1830  assert(oldpos < consdata->nlinvars);
1831  assert(newpos >= 0);
1832  assert(newpos < consdata->linvarssize);
1833 
1834  if( newpos == oldpos )
1835  return;
1836 
1837  consdata->linvars [newpos] = consdata->linvars [oldpos];
1838  consdata->lincoefs[newpos] = consdata->lincoefs[oldpos];
1839 
1840  if( consdata->lineventdata != NULL )
1841  {
1842  assert(newpos >= consdata->nlinvars || consdata->lineventdata[newpos] == NULL);
1843 
1844  consdata->lineventdata[newpos] = consdata->lineventdata[oldpos];
1845  consdata->lineventdata[newpos]->varidx = newpos;
1846 
1847  consdata->lineventdata[oldpos] = NULL;
1848  }
1849 
1850  consdata->linvarssorted = FALSE;
1851 }
1852 
1853 /** moves a quadratic variable from one position to another */
1854 static
1855 void consdataMoveQuadVarTerm(
1856  SCIP_CONSDATA* consdata, /**< constraint data */
1857  int oldpos, /**< position of variable that shall be moved */
1858  int newpos /**< new position of variable */
1859  )
1860 {
1861  assert(consdata != NULL);
1862  assert(oldpos >= 0);
1863  assert(oldpos < consdata->nquadvars);
1864  assert(newpos >= 0);
1865  assert(newpos < consdata->quadvarssize);
1866 
1867  if( newpos == oldpos )
1868  return;
1869 
1870  assert(newpos >= consdata->nquadvars || consdata->quadvarterms[newpos].eventdata == NULL);
1871 
1872  consdata->quadvarterms[newpos] = consdata->quadvarterms[oldpos];
1873 
1874  if( consdata->quadvarterms[newpos].eventdata != NULL )
1875  {
1876  consdata->quadvarterms[newpos].eventdata->varidx = -newpos-1;
1877  consdata->quadvarterms[oldpos].eventdata = NULL;
1878  }
1879 
1880  consdata->quadvarssorted = FALSE;
1881 }
1882 
1883 /** adds linear coefficient in quadratic constraint */
1884 static
1885 SCIP_RETCODE addLinearCoef(
1886  SCIP* scip, /**< SCIP data structure */
1887  SCIP_CONS* cons, /**< quadratic constraint */
1888  SCIP_VAR* var, /**< variable of constraint entry */
1889  SCIP_Real coef /**< coefficient of constraint entry */
1890  )
1891 {
1892  SCIP_CONSDATA* consdata;
1893  SCIP_Bool transformed;
1894 
1895  assert(scip != NULL);
1896  assert(cons != NULL);
1897  assert(var != NULL);
1898 
1899  /* ignore coefficient if it is nearly zero */
1900  if( SCIPisZero(scip, coef) )
1901  return SCIP_OKAY;
1902 
1903  consdata = SCIPconsGetData(cons);
1904  assert(consdata != NULL);
1905 
1906  /* are we in the transformed problem? */
1907  transformed = SCIPconsIsTransformed(cons);
1908 
1909  /* always use transformed variables in transformed constraints */
1910  if( transformed )
1911  {
1912  SCIP_CALL( SCIPgetTransformedVar(scip, var, &var) );
1913  }
1914  assert(var != NULL);
1915  assert(transformed == SCIPvarIsTransformed(var));
1916 
1917  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars+1) );
1918  consdata->linvars [consdata->nlinvars] = var;
1919  consdata->lincoefs[consdata->nlinvars] = coef;
1920 
1921  ++consdata->nlinvars;
1922 
1923  /* catch variable events */
1924  if( SCIPconsIsEnabled(cons) )
1925  {
1926  SCIP_CONSHDLR* conshdlr;
1927  SCIP_CONSHDLRDATA* conshdlrdata;
1928 
1929  /* get event handler */
1930  conshdlr = SCIPconsGetHdlr(cons);
1931  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1932  assert(conshdlrdata != NULL);
1933  assert(conshdlrdata->eventhdlr != NULL);
1934 
1935  assert(consdata->lineventdata != NULL);
1936  consdata->lineventdata[consdata->nlinvars-1] = NULL;
1937 
1938  /* catch bound change events of variable */
1939  SCIP_CALL( catchLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, consdata->nlinvars-1) );
1940  }
1941 
1942  /* invalidate activity information */
1943  consdata->activity = SCIP_INVALID;
1944  consdata->minlinactivity = SCIP_INVALID;
1945  consdata->maxlinactivity = SCIP_INVALID;
1946  consdata->minlinactivityinf = -1;
1947  consdata->maxlinactivityinf = -1;
1948 
1949  /* invalidate nonlinear row */
1950  if( consdata->nlrow != NULL )
1951  {
1952  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1953  }
1954 
1955  /* install rounding locks for new variable */
1956  SCIP_CALL( lockLinearVariable(scip, cons, var, coef) );
1957 
1958  /* capture new variable */
1959  SCIP_CALL( SCIPcaptureVar(scip, var) );
1960 
1961  consdata->ispropagated = FALSE;
1962  consdata->ispresolved = FALSE;
1963  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
1964  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
1965  if( consdata->nlinvars == 1 )
1966  consdata->linvarssorted = TRUE;
1967  else
1968  consdata->linvarssorted = consdata->linvarssorted && (SCIPvarCompare(consdata->linvars[consdata->nlinvars-2], consdata->linvars[consdata->nlinvars-1]) == -1);
1969  /* always set too FALSE since the new linear variable should be checked if already existing as quad var term */
1970  consdata->linvarsmerged = FALSE;
1971 
1972  return SCIP_OKAY;
1973 }
1974 
1975 /** deletes linear coefficient at given position from quadratic constraint data */
1976 static
1977 SCIP_RETCODE delLinearCoefPos(
1978  SCIP* scip, /**< SCIP data structure */
1979  SCIP_CONS* cons, /**< quadratic constraint */
1980  int pos /**< position of coefficient to delete */
1981  )
1982 {
1983  SCIP_CONSDATA* consdata;
1984  SCIP_VAR* var;
1985  SCIP_Real coef;
1986 
1987  assert(scip != NULL);
1988  assert(cons != NULL);
1989 
1990  consdata = SCIPconsGetData(cons);
1991  assert(consdata != NULL);
1992  assert(0 <= pos && pos < consdata->nlinvars);
1993 
1994  var = consdata->linvars[pos];
1995  coef = consdata->lincoefs[pos];
1996  assert(var != NULL);
1997 
1998  /* remove rounding locks for deleted variable */
1999  SCIP_CALL( unlockLinearVariable(scip, cons, var, coef) );
2000 
2001  /* if we catch variable events, drop the events on the variable */
2002  if( consdata->lineventdata != NULL )
2003  {
2004  SCIP_CONSHDLR* conshdlr;
2005  SCIP_CONSHDLRDATA* conshdlrdata;
2006 
2007  /* get event handler */
2008  conshdlr = SCIPconsGetHdlr(cons);
2009  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2010  assert(conshdlrdata != NULL);
2011  assert(conshdlrdata->eventhdlr != NULL);
2012 
2013  /* drop bound change events of variable */
2014  SCIP_CALL( dropLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2015  }
2016 
2017  /* release variable */
2018  SCIP_CALL( SCIPreleaseVar(scip, &consdata->linvars[pos]) );
2019 
2020  /* move the last variable to the free slot */
2021  consdataMoveLinearVar(consdata, consdata->nlinvars-1, pos);
2022 
2023  --consdata->nlinvars;
2024 
2025  /* invalidate activity */
2026  consdata->activity = SCIP_INVALID;
2027  consdata->minlinactivity = SCIP_INVALID;
2028  consdata->maxlinactivity = SCIP_INVALID;
2029  consdata->minlinactivityinf = -1;
2030  consdata->maxlinactivityinf = -1;
2031 
2032  /* invalidate nonlinear row */
2033  if( consdata->nlrow != NULL )
2034  {
2035  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2036  }
2037 
2038  consdata->ispropagated = FALSE;
2039  consdata->ispresolved = FALSE;
2040 
2041  return SCIP_OKAY;
2042 }
2043 
2044 /** changes linear coefficient value at given position of quadratic constraint */
2045 static
2046 SCIP_RETCODE chgLinearCoefPos(
2047  SCIP* scip, /**< SCIP data structure */
2048  SCIP_CONS* cons, /**< quadratic constraint */
2049  int pos, /**< position of linear coefficient to change */
2050  SCIP_Real newcoef /**< new value of linear coefficient */
2051  )
2052 {
2053  SCIP_CONSHDLR* conshdlr;
2054  SCIP_CONSHDLRDATA* conshdlrdata;
2055  SCIP_CONSDATA* consdata;
2056  SCIP_VAR* var;
2057  SCIP_Real coef;
2058 
2059  assert(scip != NULL);
2060  assert(cons != NULL);
2061  assert(!SCIPconsIsLockedType(cons, SCIP_LOCKTYPE_CONFLICT));
2062  assert(!SCIPisZero(scip, newcoef));
2063 
2064  conshdlrdata = NULL;
2065 
2066  consdata = SCIPconsGetData(cons);
2067  assert(consdata != NULL);
2068  assert(0 <= pos);
2069  assert(pos < consdata->nlinvars);
2070  assert(!SCIPisZero(scip, newcoef));
2071 
2072  var = consdata->linvars[pos];
2073  coef = consdata->lincoefs[pos];
2074  assert(var != NULL);
2075  assert(SCIPconsIsTransformed(cons) == SCIPvarIsTransformed(var));
2076 
2077  /* invalidate activity */
2078  consdata->activity = SCIP_INVALID;
2079  consdata->minlinactivity = SCIP_INVALID;
2080  consdata->maxlinactivity = SCIP_INVALID;
2081  consdata->minlinactivityinf = -1;
2082  consdata->maxlinactivityinf = -1;
2083 
2084  /* invalidate nonlinear row */
2085  if( consdata->nlrow != NULL )
2086  {
2087  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2088  }
2089 
2090  /* if necessary, remove the rounding locks and event catching of the variable */
2091  if( newcoef * coef < 0.0 )
2092  {
2093  if( SCIPconsIsLocked(cons) )
2094  {
2095  assert(SCIPconsIsTransformed(cons));
2096 
2097  /* remove rounding locks for variable with old coefficient */
2098  SCIP_CALL( unlockLinearVariable(scip, cons, var, coef) );
2099  }
2100 
2101  if( consdata->lineventdata[pos] != NULL )
2102  {
2103  /* get event handler */
2104  conshdlr = SCIPconsGetHdlr(cons);
2105  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2106  assert(conshdlrdata != NULL);
2107  assert(conshdlrdata->eventhdlr != NULL);
2108 
2109  /* drop bound change events of variable */
2110  SCIP_CALL( dropLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2111  }
2112  }
2113 
2114  /* change the coefficient */
2115  consdata->lincoefs[pos] = newcoef;
2116 
2117  /* if necessary, install the rounding locks and event catching of the variable again */
2118  if( newcoef * coef < 0.0 )
2119  {
2120  if( SCIPconsIsLocked(cons) )
2121  {
2122  /* install rounding locks for variable with new coefficient */
2123  SCIP_CALL( lockLinearVariable(scip, cons, var, newcoef) );
2124  }
2125 
2126  if( conshdlrdata != NULL )
2127  {
2128  assert(SCIPconsIsEnabled(cons));
2129 
2130  /* catch bound change events of variable */
2131  SCIP_CALL( catchLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2132  }
2133  }
2134 
2135  consdata->ispropagated = FALSE;
2136  consdata->ispresolved = FALSE;
2137 
2138  return SCIP_OKAY;
2139 }
2140 
2141 /** adds quadratic variable term to quadratic constraint */
2142 static
2143 SCIP_RETCODE addQuadVarTerm(
2144  SCIP* scip, /**< SCIP data structure */
2145  SCIP_CONS* cons, /**< quadratic constraint */
2146  SCIP_VAR* var, /**< variable to add */
2147  SCIP_Real lincoef, /**< linear coefficient of variable */
2148  SCIP_Real sqrcoef /**< square coefficient of variable */
2149  )
2150 {
2151  SCIP_CONSDATA* consdata;
2152  SCIP_Bool transformed;
2153  SCIP_QUADVARTERM* quadvarterm;
2154 
2155  assert(scip != NULL);
2156  assert(cons != NULL);
2157  assert(var != NULL);
2158 
2159  consdata = SCIPconsGetData(cons);
2160  assert(consdata != NULL);
2161 
2162  /* are we in the transformed problem? */
2163  transformed = SCIPconsIsTransformed(cons);
2164 
2165  /* always use transformed variables in transformed constraints */
2166  if( transformed )
2167  {
2168  SCIP_CALL( SCIPgetTransformedVar(scip, var, &var) );
2169  }
2170  assert(var != NULL);
2171  assert(transformed == SCIPvarIsTransformed(var));
2172 
2173  SCIP_CALL( consdataEnsureQuadVarTermsSize(scip, consdata, consdata->nquadvars+1) );
2174 
2175  quadvarterm = &consdata->quadvarterms[consdata->nquadvars];
2176  quadvarterm->var = var;
2177  quadvarterm->lincoef = lincoef;
2178  quadvarterm->sqrcoef = sqrcoef;
2179  quadvarterm->adjbilinsize = 0;
2180  quadvarterm->nadjbilin = 0;
2181  quadvarterm->adjbilin = NULL;
2182  quadvarterm->eventdata = NULL;
2183 
2184  ++consdata->nquadvars;
2185 
2186  /* capture variable */
2187  SCIP_CALL( SCIPcaptureVar(scip, var) );
2188 
2189  /* catch variable events, if we do so */
2190  if( SCIPconsIsEnabled(cons) )
2191  {
2192  SCIP_CONSHDLR* conshdlr;
2193  SCIP_CONSHDLRDATA* conshdlrdata;
2194 
2195  /* get event handler */
2196  conshdlr = SCIPconsGetHdlr(cons);
2197  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2198  assert(conshdlrdata != NULL);
2199  assert(conshdlrdata->eventhdlr != NULL);
2200 
2201  /* catch bound change events of variable */
2202  SCIP_CALL( catchQuadVarEvents(scip, conshdlrdata->eventhdlr, cons, consdata->nquadvars-1) );
2203  }
2204 
2205  /* invalidate activity information */
2206  consdata->activity = SCIP_INVALID;
2207  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2208 
2209  /* invalidate nonlinear row */
2210  if( consdata->nlrow != NULL )
2211  {
2212  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2213  }
2214 
2215  /* install rounding locks for new variable */
2216  SCIP_CALL( lockQuadraticVariable(scip, cons, var) );
2217 
2218  consdata->ispropagated = FALSE;
2219  consdata->ispresolved = FALSE;
2220  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
2221  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
2222  if( consdata->nquadvars == 1 )
2223  consdata->quadvarssorted = TRUE;
2224  else
2225  consdata->quadvarssorted = consdata->quadvarssorted &&
2226  (SCIPvarCompare(consdata->quadvarterms[consdata->nquadvars-2].var, consdata->quadvarterms[consdata->nquadvars-1].var) == -1);
2227  /* also set to FALSE if nquadvars == 1, since the new variable should be checked for linearity and other stuff in mergeAndClean ... */
2228  consdata->quadvarsmerged = FALSE;
2229 
2230  consdata->iscurvchecked = FALSE;
2231 
2232  return SCIP_OKAY;
2233 }
2234 
2235 /** deletes quadratic variable term at given position from quadratic constraint data */
2236 static
2237 SCIP_RETCODE delQuadVarTermPos(
2238  SCIP* scip, /**< SCIP data structure */
2239  SCIP_CONS* cons, /**< quadratic constraint */
2240  int pos /**< position of term to delete */
2241  )
2242 {
2243  SCIP_CONSDATA* consdata;
2244  SCIP_VAR* var;
2245 
2246  assert(scip != NULL);
2247  assert(cons != NULL);
2248 
2249  consdata = SCIPconsGetData(cons);
2250  assert(consdata != NULL);
2251  assert(0 <= pos && pos < consdata->nquadvars);
2252 
2253  var = consdata->quadvarterms[pos].var;
2254  assert(var != NULL);
2255  assert(consdata->quadvarterms[pos].nadjbilin == 0);
2256 
2257  /* remove rounding locks for deleted variable */
2258  SCIP_CALL( unlockQuadraticVariable(scip, cons, var) );
2259 
2260  /* if we catch variable events, drop the events on the variable */
2261  if( consdata->quadvarterms[pos].eventdata != NULL )
2262  {
2263  SCIP_CONSHDLR* conshdlr;
2264  SCIP_CONSHDLRDATA* conshdlrdata;
2265 
2266  /* get event handler */
2267  conshdlr = SCIPconsGetHdlr(cons);
2268  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2269  assert(conshdlrdata != NULL);
2270  assert(conshdlrdata->eventhdlr != NULL);
2271 
2272  /* drop bound change events of variable */
2273  SCIP_CALL( dropQuadVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2274  }
2275 
2276  /* release variable */
2277  SCIP_CALL( SCIPreleaseVar(scip, &consdata->quadvarterms[pos].var) );
2278 
2279  /* free adjacency array */
2280  SCIPfreeBlockMemoryArrayNull(scip, &consdata->quadvarterms[pos].adjbilin, consdata->quadvarterms[pos].adjbilinsize);
2281 
2282  /* move the last variable term to the free slot */
2283  consdataMoveQuadVarTerm(consdata, consdata->nquadvars-1, pos);
2284 
2285  --consdata->nquadvars;
2286 
2287  /* invalidate activity */
2288  consdata->activity = SCIP_INVALID;
2289  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2290 
2291  /* invalidate nonlinear row */
2292  if( consdata->nlrow != NULL )
2293  {
2294  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2295  }
2296 
2297  consdata->ispropagated = FALSE;
2298  consdata->ispresolved = FALSE;
2299  consdata->iscurvchecked = FALSE;
2300 
2301  return SCIP_OKAY;
2302 }
2303 
2304 /** replace variable in quadratic variable term at given position of quadratic constraint data
2305  *
2306  * Allows to replace x by coef*y+offset, thereby maintaining linear and square coefficients and bilinear terms.
2307  */
2308 static
2309 SCIP_RETCODE replaceQuadVarTermPos(
2310  SCIP* scip, /**< SCIP data structure */
2311  SCIP_CONS* cons, /**< quadratic constraint */
2312  int pos, /**< position of term to replace */
2313  SCIP_VAR* var, /**< new variable */
2314  SCIP_Real coef, /**< linear coefficient of new variable */
2315  SCIP_Real offset /**< offset of new variable */
2316  )
2317 {
2318  SCIP_CONSDATA* consdata;
2319  SCIP_QUADVARTERM* quadvarterm;
2320  SCIP_EVENTHDLR* eventhdlr;
2321  SCIP_BILINTERM* bilinterm;
2322  SCIP_Real constant;
2323 
2324  int i;
2325  SCIP_VAR* var2;
2326 
2327  consdata = SCIPconsGetData(cons);
2328  assert(consdata != NULL);
2329  assert(pos >= 0);
2330  assert(pos < consdata->nquadvars);
2331 
2332  quadvarterm = &consdata->quadvarterms[pos];
2333 
2334  /* remove rounding locks for old variable */
2335  SCIP_CALL( unlockQuadraticVariable(scip, cons, quadvarterm->var) );
2336 
2337  /* if we catch variable events, drop the events on the old variable */
2338  if( quadvarterm->eventdata != NULL )
2339  {
2340  SCIP_CONSHDLR* conshdlr;
2341  SCIP_CONSHDLRDATA* conshdlrdata;
2342 
2343  /* get event handler */
2344  conshdlr = SCIPconsGetHdlr(cons);
2345  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2346  assert(conshdlrdata != NULL);
2347  assert(conshdlrdata->eventhdlr != NULL);
2348 
2349  eventhdlr = conshdlrdata->eventhdlr;
2350 
2351  /* drop bound change events of variable */
2352  SCIP_CALL( dropQuadVarEvents(scip, eventhdlr, cons, pos) );
2353  }
2354  else
2355  {
2356  eventhdlr = NULL;
2357  }
2358 
2359  /* compute constant and put into lhs/rhs */
2360  constant = quadvarterm->lincoef * offset + quadvarterm->sqrcoef * offset * offset;
2361  if( constant != 0.0 )
2362  {
2363  /* maintain constant part */
2364  if( !SCIPisInfinity(scip, -consdata->lhs) )
2365  consdata->lhs -= constant;
2366  if( !SCIPisInfinity(scip, consdata->rhs) )
2367  consdata->rhs -= constant;
2368  }
2369 
2370  /* update linear and square coefficient */
2371  quadvarterm->lincoef *= coef;
2372  quadvarterm->lincoef += 2.0 * quadvarterm->sqrcoef * coef * offset;
2373  quadvarterm->sqrcoef *= coef * coef;
2374 
2375  /* update bilinear terms */
2376  for( i = 0; i < quadvarterm->nadjbilin; ++i )
2377  {
2378  bilinterm = &consdata->bilinterms[quadvarterm->adjbilin[i]];
2379 
2380  if( bilinterm->var1 == quadvarterm->var )
2381  {
2382  bilinterm->var1 = var;
2383  var2 = bilinterm->var2;
2384  }
2385  else
2386  {
2387  assert(bilinterm->var2 == quadvarterm->var);
2388  bilinterm->var2 = var;
2389  var2 = bilinterm->var1;
2390  }
2391 
2392  if( var == var2 )
2393  {
2394  /* looks like we actually have a square term here */
2395  quadvarterm->lincoef += bilinterm->coef * offset;
2396  quadvarterm->sqrcoef += bilinterm->coef * coef;
2397  /* deleting bilinear terms is expensive, since it requires updating adjacency information
2398  * thus, for now we just set the coefficient to 0.0 and clear in later when the bilinear terms are merged */
2399  bilinterm->coef = 0.0;
2400  continue;
2401  }
2402 
2403  /* swap var1 and var2 if they are in wrong order */
2404  if( SCIPvarCompare(bilinterm->var1, bilinterm->var2) > 0 )
2405  {
2406  SCIP_VAR* tmp;
2407  tmp = bilinterm->var1;
2408  bilinterm->var1 = bilinterm->var2;
2409  bilinterm->var2 = tmp;
2410  }
2411  assert(SCIPvarCompare(bilinterm->var1, bilinterm->var2) == -1);
2412 
2413  if( offset != 0.0 )
2414  {
2415  /* need to find var2 and add offset*bilinterm->coef to linear coefficient */
2416  int var2pos;
2417 
2418  var2pos = 0;
2419  while( consdata->quadvarterms[var2pos].var != var2 )
2420  {
2421  ++var2pos;
2422  assert(var2pos < consdata->nquadvars);
2423  }
2424 
2425  consdata->quadvarterms[var2pos].lincoef += bilinterm->coef * offset;
2426  }
2427 
2428  bilinterm->coef *= coef;
2429  }
2430 
2431  /* release old variable */
2432  SCIP_CALL( SCIPreleaseVar(scip, &quadvarterm->var) );
2433 
2434  /* set new variable */
2435  quadvarterm->var = var;
2436 
2437  /* capture new variable */
2438  SCIP_CALL( SCIPcaptureVar(scip, quadvarterm->var) );
2439 
2440  /* catch variable events, if we do so */
2441  if( eventhdlr != NULL )
2442  {
2443  assert(SCIPconsIsEnabled(cons));
2444 
2445  /* catch bound change events of variable */
2446  SCIP_CALL( catchQuadVarEvents(scip, eventhdlr, cons, pos) );
2447  }
2448 
2449  /* invalidate activity information */
2450  consdata->activity = SCIP_INVALID;
2451  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2452 
2453  /* invalidate nonlinear row */
2454  if( consdata->nlrow != NULL )
2455  {
2456  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2457  }
2458 
2459  /* install rounding locks for new variable */
2460  SCIP_CALL( lockQuadraticVariable(scip, cons, var) );
2461 
2462  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
2463  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
2464  consdata->quadvarssorted = (consdata->nquadvars == 1);
2465  consdata->quadvarsmerged = FALSE;
2466  consdata->bilinsorted &= (quadvarterm->nadjbilin == 0); /*lint !e514*/
2467  consdata->bilinmerged &= (quadvarterm->nadjbilin == 0); /*lint !e514*/
2468 
2469  consdata->ispropagated = FALSE;
2470  consdata->ispresolved = FALSE;
2471  consdata->iscurvchecked = FALSE;
2472 
2473  return SCIP_OKAY;
2474 }
2475 
2476 /** adds a bilinear term to quadratic constraint */
2477 static
2478 SCIP_RETCODE addBilinearTerm(
2479  SCIP* scip, /**< SCIP data structure */
2480  SCIP_CONS* cons, /**< quadratic constraint */
2481  int var1pos, /**< position of first variable in quadratic variables array */
2482  int var2pos, /**< position of second variable in quadratic variables array */
2483  SCIP_Real coef /**< coefficient of bilinear term */
2484  )
2485 {
2486  SCIP_CONSDATA* consdata;
2487  SCIP_BILINTERM* bilinterm;
2488 
2489  assert(scip != NULL);
2490  assert(cons != NULL);
2491 
2492  if( var1pos == var2pos )
2493  {
2494  SCIPerrorMessage("tried to add bilinear term where both variables are the same\n");
2495  return SCIP_INVALIDDATA;
2496  }
2497 
2498  consdata = SCIPconsGetData(cons);
2499  assert(consdata != NULL);
2500 
2501  /* check if the bilinear terms are sorted */
2502  assert(consdataCheckBilinTermsSort(consdata));
2503 
2504  assert(var1pos >= 0);
2505  assert(var1pos < consdata->nquadvars);
2506  assert(var2pos >= 0);
2507  assert(var2pos < consdata->nquadvars);
2508 
2509  SCIP_CALL( consdataEnsureBilinSize(scip, consdata, consdata->nbilinterms + 1) );
2510 
2511  bilinterm = &consdata->bilinterms[consdata->nbilinterms];
2512  if( SCIPvarCompare(consdata->quadvarterms[var1pos].var, consdata->quadvarterms[var2pos].var) < 0 )
2513  {
2514  bilinterm->var1 = consdata->quadvarterms[var1pos].var;
2515  bilinterm->var2 = consdata->quadvarterms[var2pos].var;
2516  }
2517  else
2518  {
2519  bilinterm->var1 = consdata->quadvarterms[var2pos].var;
2520  bilinterm->var2 = consdata->quadvarterms[var1pos].var;
2521  }
2522  bilinterm->coef = coef;
2523 
2524  if( bilinterm->var1 == bilinterm->var2 )
2525  {
2526  SCIPerrorMessage("tried to add bilinear term where both variables are the same, but appear at different positions in quadvarterms array\n");
2527  return SCIP_INVALIDDATA;
2528  }
2529  assert(SCIPvarCompare(bilinterm->var1, bilinterm->var2) == -1);
2530 
2531  SCIP_CALL( consdataEnsureAdjBilinSize(scip, &consdata->quadvarterms[var1pos], consdata->quadvarterms[var1pos].nadjbilin + 1) );
2532  SCIP_CALL( consdataEnsureAdjBilinSize(scip, &consdata->quadvarterms[var2pos], consdata->quadvarterms[var2pos].nadjbilin + 1) );
2533 
2534  consdata->quadvarterms[var1pos].adjbilin[consdata->quadvarterms[var1pos].nadjbilin] = consdata->nbilinterms;
2535  consdata->quadvarterms[var2pos].adjbilin[consdata->quadvarterms[var2pos].nadjbilin] = consdata->nbilinterms;
2536  ++consdata->quadvarterms[var1pos].nadjbilin;
2537  ++consdata->quadvarterms[var2pos].nadjbilin;
2538 
2539  ++consdata->nbilinterms;
2540 
2541  /* invalidate activity information */
2542  consdata->activity = SCIP_INVALID;
2543  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2544 
2545  /* invalidate nonlinear row */
2546  if( consdata->nlrow != NULL )
2547  {
2548  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2549  }
2550 
2551  consdata->ispropagated = FALSE;
2552  consdata->ispresolved = FALSE;
2553  if( consdata->nbilinterms == 1 )
2554  {
2555  consdata->bilinsorted = TRUE;
2556 
2557  /* we have to take care of the bilinear term in mergeAndCleanBilinearTerms() if the coefficient is zero */
2558  consdata->bilinmerged = !SCIPisZero(scip, consdata->bilinterms[0].coef);
2559  }
2560  else
2561  {
2562  consdata->bilinsorted = consdata->bilinsorted
2563  && (bilinTermComp(consdata, consdata->nbilinterms-2, consdata->nbilinterms-1) <= 0);
2564  consdata->bilinmerged = FALSE;
2565  }
2566 
2567  consdata->iscurvchecked = FALSE;
2568 
2569  /* check if the bilinear terms are sorted */
2570  assert(consdataCheckBilinTermsSort(consdata));
2571 
2572  return SCIP_OKAY;
2573 }
2574 
2575 /** removes a set of bilinear terms and updates adjacency information in quad var terms
2576  *
2577  * Note: this function sorts the given array termposs.
2578  */
2579 static
2580 SCIP_RETCODE removeBilinearTermsPos(
2581  SCIP* scip, /**< SCIP data structure */
2582  SCIP_CONS* cons, /**< quadratic constraint */
2583  int nterms, /**< number of terms to delete */
2584  int* termposs /**< indices of terms to delete */
2585  )
2586 {
2587  SCIP_CONSDATA* consdata;
2588  int* newpos;
2589  int i;
2590  int j;
2591  int offset;
2592 
2593  assert(scip != NULL);
2594  assert(cons != NULL);
2595  assert(nterms == 0 || termposs != NULL);
2596 
2597  if( nterms == 0 )
2598  return SCIP_OKAY;
2599 
2600  consdata = SCIPconsGetData(cons);
2601  assert(consdata != NULL);
2602 
2603  SCIPsortInt(termposs, nterms);
2604 
2605  SCIP_CALL( SCIPallocBufferArray(scip, &newpos, consdata->nbilinterms) );
2606 
2607  i = 0;
2608  offset = 0;
2609  for( j = 0; j < consdata->nbilinterms; ++j )
2610  {
2611  /* if j'th term is deleted, increase offset and continue */
2612  if( i < nterms && j == termposs[i] )
2613  {
2614  ++offset;
2615  ++i;
2616  newpos[j] = -1;
2617  continue;
2618  }
2619 
2620  /* otherwise, move it forward and remember new position */
2621  if( offset > 0 )
2622  consdata->bilinterms[j-offset] = consdata->bilinterms[j];
2623  newpos[j] = j - offset;
2624  }
2625  assert(offset == nterms);
2626 
2627  /* update adjacency and activity information in quad var terms */
2628  for( i = 0; i < consdata->nquadvars; ++i )
2629  {
2630  offset = 0;
2631  for( j = 0; j < consdata->quadvarterms[i].nadjbilin; ++j )
2632  {
2633  assert(consdata->quadvarterms[i].adjbilin[j] < consdata->nbilinterms);
2634  if( newpos[consdata->quadvarterms[i].adjbilin[j]] == -1 )
2635  {
2636  /* corresponding bilinear term was deleted, thus increase offset */
2637  ++offset;
2638  }
2639  else
2640  {
2641  /* update index of j'th bilinear term and store at position j-offset */
2642  consdata->quadvarterms[i].adjbilin[j-offset] = newpos[consdata->quadvarterms[i].adjbilin[j]];
2643  }
2644  }
2645  consdata->quadvarterms[i].nadjbilin -= offset;
2646  /* some bilinear term was removed, so invalidate activity bounds */
2647  }
2648 
2649  consdata->nbilinterms -= nterms;
2650 
2651  SCIPfreeBufferArray(scip, &newpos);
2652 
2653  /* some quad vars may be linear now */
2654  consdata->quadvarsmerged = FALSE;
2655 
2656  consdata->ispropagated = FALSE;
2657  consdata->ispresolved = FALSE;
2658  consdata->iscurvchecked = FALSE;
2659 
2660  /* invalidate activity */
2661  consdata->activity = SCIP_INVALID;
2662  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2663 
2664  /* invalidate nonlinear row */
2665  if( consdata->nlrow != NULL )
2666  {
2667  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2668  }
2669 
2670  return SCIP_OKAY;
2671 }
2672 
2673 /** merges quad var terms that correspond to the same variable and does additional cleanup
2674  *
2675  * If a quadratic variable terms is actually linear, makes a linear term out of it
2676  * also replaces squares of binary variables by the binary variables, i.e., adds sqrcoef to lincoef.
2677  */
2678 static
2679 SCIP_RETCODE mergeAndCleanQuadVarTerms(
2680  SCIP* scip, /**< SCIP data structure */
2681  SCIP_CONS* cons /**< quadratic constraint */
2682  )
2683 {
2684  SCIP_QUADVARTERM* quadvarterm;
2685  SCIP_CONSDATA* consdata;
2686  int i;
2687  int j;
2688 
2689  assert(scip != NULL);
2690  assert(cons != NULL);
2691 
2692  consdata = SCIPconsGetData(cons);
2693 
2694  if( consdata->quadvarsmerged )
2695  return SCIP_OKAY;
2696 
2697  if( consdata->nquadvars == 0 )
2698  {
2699  consdata->quadvarsmerged = TRUE;
2700  return SCIP_OKAY;
2701  }
2702 
2703  i = 0;
2704  while( i < consdata->nquadvars )
2705  {
2706  /* make sure quad var terms are sorted (do this in every round, since we may move variables around) */
2707  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
2708 
2709  quadvarterm = &consdata->quadvarterms[i];
2710 
2711  for( j = i+1; j < consdata->nquadvars && consdata->quadvarterms[j].var == quadvarterm->var; ++j )
2712  {
2713  /* add quad var term j to current term i */
2714  quadvarterm->lincoef += consdata->quadvarterms[j].lincoef;
2715  quadvarterm->sqrcoef += consdata->quadvarterms[j].sqrcoef;
2716  if( consdata->quadvarterms[j].nadjbilin > 0 )
2717  {
2718  /* move adjacency information from j to i */
2719  SCIP_CALL( consdataEnsureAdjBilinSize(scip, quadvarterm, quadvarterm->nadjbilin + consdata->quadvarterms[j].nadjbilin) );
2720  BMScopyMemoryArray(&quadvarterm->adjbilin[quadvarterm->nadjbilin], consdata->quadvarterms[j].adjbilin, consdata->quadvarterms[j].nadjbilin); /*lint !e866*/
2721  quadvarterm->nadjbilin += consdata->quadvarterms[j].nadjbilin;
2722  consdata->quadvarterms[j].nadjbilin = 0;
2723  }
2724  consdata->quadvarterms[j].lincoef = 0.0;
2725  consdata->quadvarterms[j].sqrcoef = 0.0;
2726  /* mark that activity information in quadvarterm is not up to date anymore */
2727  }
2728 
2729  /* remove quad var terms i+1..j-1 backwards */
2730  for( j = j-1; j > i; --j )
2731  {
2732  SCIP_CALL( delQuadVarTermPos(scip, cons, j) );
2733  }
2734 
2735  /* for binary variables, x^2 = x
2736  * however, we may destroy convexity of a quadratic term that involves also bilinear terms
2737  * thus, we do this step only if the variable does not appear in any bilinear term */
2738  if( quadvarterm->sqrcoef != 0.0 && SCIPvarIsBinary(quadvarterm->var) && quadvarterm->nadjbilin == 0 )
2739  {
2740  SCIPdebugMsg(scip, "replace square of binary variable by itself: <%s>^2 --> <%s>\n", SCIPvarGetName(quadvarterm->var), SCIPvarGetName(quadvarterm->var));
2741  quadvarterm->lincoef += quadvarterm->sqrcoef;
2742  quadvarterm->sqrcoef = 0.0;
2743 
2744  /* invalidate nonlinear row */
2745  if( consdata->nlrow != NULL )
2746  {
2747  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2748  }
2749  }
2750 
2751  /* if its 0.0 or linear, get rid of it */
2752  if( SCIPisZero(scip, quadvarterm->sqrcoef) && quadvarterm->nadjbilin == 0 )
2753  {
2754  if( !SCIPisZero(scip, quadvarterm->lincoef) )
2755  {
2756  /* seem to be a linear term now, thus add as linear term */
2757  SCIP_CALL( addLinearCoef(scip, cons, quadvarterm->var, quadvarterm->lincoef) );
2758  }
2759  /* remove term at pos i */
2760  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
2761  }
2762  else
2763  {
2764  ++i;
2765  }
2766  }
2767 
2768  consdata->quadvarsmerged = TRUE;
2769  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2770 
2771  return SCIP_OKAY;
2772 }
2773 
2774 /** merges entries with same linear variable into one entry and cleans up entries with coefficient 0.0 */
2775 static
2776 SCIP_RETCODE mergeAndCleanLinearVars(
2777  SCIP* scip, /**< SCIP data structure */
2778  SCIP_CONS* cons /**< quadratic constraint */
2779  )
2780 {
2781  SCIP_CONSDATA* consdata;
2782  SCIP_Real newcoef;
2783  int i;
2784  int j;
2785  int qvarpos;
2786 
2787  assert(scip != NULL);
2788  assert(cons != NULL);
2789 
2790  consdata = SCIPconsGetData(cons);
2791 
2792  if( consdata->linvarsmerged )
2793  return SCIP_OKAY;
2794 
2795  if( consdata->nlinvars == 0 )
2796  {
2797  consdata->linvarsmerged = TRUE;
2798  return SCIP_OKAY;
2799  }
2800 
2801  i = 0;
2802  while( i < consdata->nlinvars )
2803  {
2804  /* make sure linear variables are sorted (do this in every round, since we may move variables around) */
2805  consdataSortLinearVars(consdata);
2806 
2807  /* sum up coefficients that correspond to variable i */
2808  newcoef = consdata->lincoefs[i];
2809  for( j = i+1; j < consdata->nlinvars && consdata->linvars[i] == consdata->linvars[j]; ++j )
2810  newcoef += consdata->lincoefs[j];
2811  /* delete the additional variables in backward order */
2812  for( j = j-1; j > i; --j )
2813  {
2814  SCIP_CALL( delLinearCoefPos(scip, cons, j) );
2815  }
2816 
2817  /* check if there is already a quadratic variable term with this variable */
2818  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, consdata->linvars[i], &qvarpos) );
2819  if( qvarpos >= 0)
2820  {
2821  /* add newcoef to linear coefficient of quadratic variable and mark linear variable as to delete */
2822  assert(qvarpos < consdata->nquadvars);
2823  assert(consdata->quadvarterms[qvarpos].var == consdata->linvars[i]);
2824  consdata->quadvarterms[qvarpos].lincoef += newcoef;
2825  newcoef = 0.0;
2826  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2827  }
2828 
2829  /* delete also entry at position i, if it became zero (or was zero before) */
2830  if( SCIPisZero(scip, newcoef) )
2831  {
2832  SCIP_CALL( delLinearCoefPos(scip, cons, i) );
2833  }
2834  else
2835  {
2836  SCIP_CALL( chgLinearCoefPos(scip, cons, i, newcoef) );
2837  ++i;
2838  }
2839  }
2840 
2841  consdata->linvarsmerged = TRUE;
2842 
2843  return SCIP_OKAY;
2844 }
2845 
2846 /** merges bilinear terms with same variables into a single term, removes bilinear terms with coefficient 0.0 */
2847 static
2848 SCIP_RETCODE mergeAndCleanBilinearTerms(
2849  SCIP* scip, /**< SCIP data structure */
2850  SCIP_CONS* cons /**< quadratic constraint */
2851  )
2852 {
2853  SCIP_CONSDATA* consdata;
2854  SCIP_BILINTERM* bilinterm;
2855  int i;
2856  int j;
2857  int* todelete;
2858  int ntodelete;
2859 
2860  assert(scip != NULL);
2861  assert(cons != NULL);
2862 
2863  consdata = SCIPconsGetData(cons);
2864 
2865  /* check if the bilinear terms are sorted */
2866  assert(consdataCheckBilinTermsSort(consdata));
2867 
2868  if( consdata->bilinmerged )
2869  return SCIP_OKAY;
2870 
2871  if( consdata->nbilinterms == 0 )
2872  {
2873  consdata->bilinmerged = TRUE;
2874  return SCIP_OKAY;
2875  }
2876 
2877  /* alloc memory for array of terms that need to be deleted finally */
2878  ntodelete = 0;
2879  SCIP_CALL( SCIPallocBufferArray(scip, &todelete, consdata->nbilinterms) );
2880 
2881  /* make sure bilinear terms are sorted */
2882  SCIP_CALL( consdataSortBilinTerms(scip, consdata) );
2883 
2884  i = 0;
2885  while( i < consdata->nbilinterms )
2886  {
2887  bilinterm = &consdata->bilinterms[i];
2888 
2889  /* sum up coefficients that correspond to same variables as term i */
2890  for( j = i+1; j < consdata->nbilinterms && bilinterm->var1 == consdata->bilinterms[j].var1 && bilinterm->var2 == consdata->bilinterms[j].var2; ++j )
2891  {
2892  bilinterm->coef += consdata->bilinterms[j].coef;
2893  todelete[ntodelete++] = j;
2894  }
2895 
2896  /* delete also entry at position i, if it became zero (or was zero before) */
2897  if( SCIPisZero(scip, bilinterm->coef) )
2898  {
2899  todelete[ntodelete++] = i;
2900  }
2901 
2902  /* continue with term after the current series */
2903  i = j;
2904  }
2905 
2906  /* delete bilinear terms */
2907  SCIP_CALL( removeBilinearTermsPos(scip, cons, ntodelete, todelete) );
2908 
2909  SCIPfreeBufferArray(scip, &todelete);
2910 
2911  consdata->bilinmerged = TRUE;
2912 
2913  /* check if the bilinear terms are sorted */
2914  assert(consdataCheckBilinTermsSort(consdata));
2915 
2916  return SCIP_OKAY;
2917 }
2918 
2919 /** removes fixes (or aggregated) variables from a quadratic constraint */
2920 static
2921 SCIP_RETCODE removeFixedVariables(
2922  SCIP* scip, /**< SCIP data structure */
2923  SCIP_CONS* cons /**< quadratic constraint */
2924  )
2925 {
2926  SCIP_CONSDATA* consdata;
2927  SCIP_BILINTERM* bilinterm;
2928  SCIP_Real bilincoef;
2929  SCIP_Real coef;
2930  SCIP_Real offset;
2931  SCIP_VAR* var;
2932  SCIP_VAR* var2;
2933  int var2pos;
2934  int i;
2935  int j;
2936  int k;
2937 
2938  SCIP_Bool have_change;
2939 
2940  assert(scip != NULL);
2941  assert(cons != NULL);
2942 
2943  consdata = SCIPconsGetData(cons);
2944 
2945  have_change = FALSE;
2946  i = 0;
2947  while( i < consdata->nlinvars )
2948  {
2949  var = consdata->linvars[i];
2950 
2951  if( SCIPvarIsActive(var) && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
2952  {
2953  ++i;
2954  continue;
2955  }
2956 
2957  have_change = TRUE;
2958 
2959  coef = consdata->lincoefs[i];
2960  offset = 0.0;
2961 
2962  if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
2963  {
2964  offset = coef * (SCIPvarGetLbGlobal(var) + SCIPvarGetUbGlobal(var)) / 2.0;
2965  coef = 0.0;
2966  }
2967  else
2968  {
2969  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &offset) );
2970  }
2971 
2972  SCIPdebugMsg(scip, " linear term %g*<%s> is replaced by %g * <%s> + %g\n", consdata->lincoefs[i], SCIPvarGetName(consdata->linvars[i]),
2973  coef, SCIPvarGetName(var), offset);
2974 
2975  /* delete previous variable (this will move another variable to position i) */
2976  SCIP_CALL( delLinearCoefPos(scip, cons, i) );
2977 
2978  /* put constant part into bounds */
2979  if( offset != 0.0 )
2980  {
2981  if( !SCIPisInfinity(scip, -consdata->lhs) )
2982  consdata->lhs -= offset;
2983  if( !SCIPisInfinity(scip, consdata->rhs) )
2984  consdata->rhs -= offset;
2985  }
2986 
2987  /* nothing left to do if variable had been fixed */
2988  if( coef == 0.0 )
2989  continue;
2990 
2991  /* if GetProbvar gave a linear variable, just add it
2992  * if it's a multilinear variable, add it's disaggregated variables */
2993  if( SCIPvarIsActive(var) )
2994  {
2995  SCIP_CALL( addLinearCoef(scip, cons, var, coef) );
2996  }
2997  else
2998  {
2999  int naggrs;
3000  SCIP_VAR** aggrvars;
3001  SCIP_Real* aggrscalars;
3002  SCIP_Real aggrconstant;
3003 
3004  assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR);
3005 
3006  naggrs = SCIPvarGetMultaggrNVars(var);
3007  aggrvars = SCIPvarGetMultaggrVars(var);
3008  aggrscalars = SCIPvarGetMultaggrScalars(var);
3009  aggrconstant = SCIPvarGetMultaggrConstant(var);
3010 
3011  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars + naggrs) );
3012 
3013  for( j = 0; j < naggrs; ++j )
3014  {
3015  SCIP_CALL( addLinearCoef(scip, cons, aggrvars[j], coef * aggrscalars[j]) );
3016  }
3017 
3018  if( aggrconstant != 0.0 )
3019  {
3020  if( !SCIPisInfinity(scip, -consdata->lhs) )
3021  consdata->lhs -= coef * aggrconstant;
3022  if( !SCIPisInfinity(scip, consdata->rhs) )
3023  consdata->rhs -= coef * aggrconstant;
3024  }
3025  }
3026  }
3027 
3028  i = 0;
3029  while( i < consdata->nquadvars )
3030  {
3031  var = consdata->quadvarterms[i].var;
3032 
3033  if( SCIPvarIsActive(var) && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
3034  {
3035  ++i;
3036  continue;
3037  }
3038 
3039  have_change = TRUE;
3040 
3041  coef = 1.0;
3042  offset = 0.0;
3043 
3044  if( !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
3045  {
3046  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &offset) );
3047  }
3048  else
3049  {
3050  coef = 0.0;
3051  offset = (SCIPvarGetLbGlobal(var) + SCIPvarGetUbGlobal(var)) / 2.0;
3052  }
3053 
3054  SCIPdebugMsg(scip, " quadratic variable <%s> with status %d is replaced by %g * <%s> + %g\n", SCIPvarGetName(consdata->quadvarterms[i].var),
3055  SCIPvarGetStatus(consdata->quadvarterms[i].var), coef, SCIPvarGetName(var), offset);
3056 
3057  /* handle fixed variable */
3058  if( coef == 0.0 )
3059  {
3060  /* if not fixed to 0.0, add to linear coefs of vars in bilinear terms, and deal with linear and square term as constant */
3061  if( offset != 0.0 )
3062  {
3063  for( j = 0; j < consdata->quadvarterms[i].nadjbilin; ++j )
3064  {
3065  bilinterm = &consdata->bilinterms[consdata->quadvarterms[i].adjbilin[j]];
3066 
3067  var2 = bilinterm->var1 == consdata->quadvarterms[i].var ? bilinterm->var2 : bilinterm->var1;
3068  assert(var2 != consdata->quadvarterms[i].var);
3069 
3070  var2pos = 0;
3071  while( consdata->quadvarterms[var2pos].var != var2 )
3072  {
3073  ++var2pos;
3074  assert(var2pos < consdata->nquadvars);
3075  }
3076  consdata->quadvarterms[var2pos].lincoef += bilinterm->coef * offset;
3077  }
3078 
3079  offset = consdata->quadvarterms[i].lincoef * offset + consdata->quadvarterms[i].sqrcoef * offset * offset;
3080  if( !SCIPisInfinity(scip, -consdata->lhs) )
3081  consdata->lhs -= offset;
3082  if( !SCIPisInfinity(scip, consdata->rhs) )
3083  consdata->rhs -= offset;
3084  }
3085 
3086  /* remove bilinear terms */
3087  SCIP_CALL( removeBilinearTermsPos(scip, cons, consdata->quadvarterms[i].nadjbilin, consdata->quadvarterms[i].adjbilin) );
3088 
3089  /* delete quad. var term i */
3090  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
3091 
3092  continue;
3093  }
3094 
3095  assert(var != NULL);
3096 
3097  /* if GetProbvar gave an active variable, replace the quad var term so that it uses the new variable */
3098  if( SCIPvarIsActive(var) )
3099  {
3100  /* replace x by coef*y+offset */
3101  SCIP_CALL( replaceQuadVarTermPos(scip, cons, i, var, coef, offset) );
3102 
3103  continue;
3104  }
3105  else
3106  {
3107  /* if GetProbVar gave a multi-aggregated variable, add new quad var terms and new bilinear terms
3108  * x is replaced by coef * (sum_i a_ix_i + b) + offset
3109  * lcoef * x + scoef * x^2 + bcoef * x * y ->
3110  * (b*coef + offset) * (lcoef + (b*coef + offset) * scoef)
3111  * + sum_i a_i*coef * (lcoef + 2 (b*coef + offset) * scoef) x_i
3112  * + sum_i (a_i*coef)^2 * scoef * x_i^2
3113  * + 2 sum_{i,j, i<j} (a_i a_j coef^2 scoef) x_i x_j
3114  * + bcoef * (b*coef + offset + coef * sum_i a_ix_i) y
3115  */
3116  int naggrs;
3117  SCIP_VAR** aggrvars; /* x_i */
3118  SCIP_Real* aggrscalars; /* a_i */
3119  SCIP_Real aggrconstant; /* b */
3120  int nquadtermsold;
3121 
3122  SCIP_Real lcoef;
3123  SCIP_Real scoef;
3124 
3125  assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR);
3126 
3127  naggrs = SCIPvarGetMultaggrNVars(var);
3128  aggrvars = SCIPvarGetMultaggrVars(var);
3129  aggrscalars = SCIPvarGetMultaggrScalars(var);
3130  aggrconstant = SCIPvarGetMultaggrConstant(var);
3131 
3132  lcoef = consdata->quadvarterms[i].lincoef;
3133  scoef = consdata->quadvarterms[i].sqrcoef;
3134 
3135  nquadtermsold = consdata->nquadvars;
3136 
3137  SCIP_CALL( consdataEnsureQuadVarTermsSize(scip, consdata, consdata->nquadvars + naggrs) );
3138 
3139  /* take care of constant part */
3140  if( aggrconstant != 0.0 || offset != 0.0 )
3141  {
3142  SCIP_Real constant;
3143  constant = (aggrconstant * coef + offset) * (lcoef + (aggrconstant * coef + offset) * scoef);
3144  if( !SCIPisInfinity(scip, -consdata->lhs) )
3145  consdata->lhs -= constant;
3146  if( !SCIPisInfinity(scip, consdata->rhs) )
3147  consdata->rhs -= constant;
3148  }
3149 
3150  /* add x_i's with linear and square coefficients */
3151  for( j = 0; j < naggrs; ++j )
3152  {
3153  SCIP_CALL( addQuadVarTerm(scip, cons, aggrvars[j],
3154  coef * aggrscalars[j] * (lcoef + 2.0 * scoef * (coef * aggrconstant + offset)),
3155  coef * coef * aggrscalars[j] * aggrscalars[j] * scoef) );
3156  }
3157 
3158  /* ensure space for bilinear terms */
3159  SCIP_CALL( consdataEnsureBilinSize(scip, consdata, consdata->nquadvars + (scoef != 0.0 ? (naggrs * (naggrs-1))/2 : 0) + consdata->quadvarterms[j].nadjbilin * naggrs) );
3160 
3161  /* add x_j*x_k's */
3162  if( scoef != 0.0 )
3163  {
3164  for( j = 0; j < naggrs; ++j )
3165  for( k = 0; k < j; ++k )
3166  {
3167  assert(aggrvars[j] != aggrvars[k]);
3168  SCIP_CALL( addBilinearTerm(scip, cons, nquadtermsold + j, nquadtermsold + k,
3169  2.0 * aggrscalars[j] * aggrscalars[k] * coef * coef * scoef) );
3170  }
3171  }
3172 
3173  /* add x_i*y's */
3174  for( k = 0; k < consdata->quadvarterms[i].nadjbilin; ++k )
3175  {
3176  bilinterm = &consdata->bilinterms[consdata->quadvarterms[i].adjbilin[k]];
3177  bilincoef = bilinterm->coef; /* copy coef, as bilinterm pointer may become invalid by realloc in addBilinearTerm() below */
3178  var2 = (bilinterm->var1 == consdata->quadvarterms[i].var) ? bilinterm->var2 : bilinterm->var1;
3179  assert(var2 != consdata->quadvarterms[i].var);
3180 
3181  /* this is not efficient, but we cannot sort the quadratic terms here, since we currently iterate over them */
3182  var2pos = 0;
3183  while( consdata->quadvarterms[var2pos].var != var2 )
3184  {
3185  ++var2pos;
3186  assert(var2pos < consdata->nquadvars);
3187  }
3188 
3189  for( j = 0; j < naggrs; ++j )
3190  {
3191  if( aggrvars[j] == var2 )
3192  { /* x_i == y, so we have a square term here */
3193  consdata->quadvarterms[var2pos].sqrcoef += bilincoef * coef * aggrscalars[j];
3194  }
3195  else
3196  { /* x_i != y, so we need to add a bilinear term here */
3197  SCIP_CALL( addBilinearTerm(scip, cons, nquadtermsold + j, var2pos, bilincoef * coef * aggrscalars[j]) );
3198  }
3199  }
3200 
3201  consdata->quadvarterms[var2pos].lincoef += bilincoef * (aggrconstant * coef + offset);
3202  }
3203 
3204  /* remove bilinear terms */
3205  SCIP_CALL( removeBilinearTermsPos(scip, cons, consdata->quadvarterms[i].nadjbilin, consdata->quadvarterms[i].adjbilin) );
3206 
3207  /* delete quad. var term i */
3208  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
3209  }
3210  }
3211 
3212  consdata->isremovedfixings = TRUE;
3213 
3214  SCIPdebugMsg(scip, "removed fixations from <%s>\n -> ", SCIPconsGetName(cons));
3215  SCIPdebugPrintCons(scip, cons, NULL);
3216 
3217 #ifndef NDEBUG
3218  for( i = 0; i < consdata->nlinvars; ++i )
3219  assert(SCIPvarIsActive(consdata->linvars[i]));
3220 
3221  for( i = 0; i < consdata->nquadvars; ++i )
3222  assert(SCIPvarIsActive(consdata->quadvarterms[i].var));
3223 #endif
3224 
3225  if( !have_change )
3226  return SCIP_OKAY;
3227 
3228  /* some quadratic variable may have been replaced by an already existing linear variable
3229  * in this case, we want the linear variable to be removed, which happens in mergeAndCleanLinearVars
3230  */
3231  consdata->linvarsmerged = FALSE;
3232 
3233  SCIP_CALL( mergeAndCleanBilinearTerms(scip, cons) );
3234  SCIP_CALL( mergeAndCleanQuadVarTerms(scip, cons) );
3235  SCIP_CALL( mergeAndCleanLinearVars(scip, cons) );
3236 
3237 #ifndef NDEBUG
3238  for( i = 0; i < consdata->nbilinterms; ++i )
3239  {
3240  assert(consdata->bilinterms[i].var1 != consdata->bilinterms[i].var2);
3241  assert(consdata->bilinterms[i].coef != 0.0);
3242  assert(SCIPvarCompare(consdata->bilinterms[i].var1, consdata->bilinterms[i].var2) < 0);
3243  }
3244 #endif
3245 
3246  return SCIP_OKAY;
3247 }
3248 
3249 /** create a nonlinear row representation of the constraint and stores them in consdata */
3250 static
3251 SCIP_RETCODE createNlRow(
3252  SCIP* scip, /**< SCIP data structure */
3253  SCIP_CONS* cons /**< quadratic constraint */
3254  )
3255 {
3256  SCIP_CONSDATA* consdata;
3257  int nquadvars; /* number of variables in quadratic terms */
3258  SCIP_VAR** quadvars; /* variables in quadratic terms */
3259  int nquadelems; /* number of quadratic elements (square and bilinear terms) */
3260  SCIP_QUADELEM* quadelems; /* quadratic elements (square and bilinear terms) */
3261  int nquadlinterms; /* number of linear terms using variables that are in quadratic terms */
3262  SCIP_VAR** quadlinvars; /* variables of linear terms using variables that are in quadratic terms */
3263  SCIP_Real* quadlincoefs; /* coefficients of linear terms using variables that are in quadratic terms */
3264  int i;
3265  int idx1;
3266  int idx2;
3267  int lincnt;
3268  int elcnt;
3269  SCIP_VAR* lastvar;
3270  int lastvaridx;
3271  SCIP_EXPRCURV curvature;
3272 
3273  assert(scip != NULL);
3274  assert(cons != NULL);
3275 
3276  consdata = SCIPconsGetData(cons);
3277  assert(consdata != NULL);
3278 
3279  if( consdata->nlrow != NULL )
3280  {
3281  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
3282  }
3283 
3284  nquadvars = consdata->nquadvars;
3285  nquadelems = consdata->nbilinterms;
3286  nquadlinterms = 0;
3287  for( i = 0; i < nquadvars; ++i )
3288  {
3289  if( consdata->quadvarterms[i].sqrcoef != 0.0 )
3290  ++nquadelems;
3291  if( !SCIPisZero(scip, consdata->quadvarterms[i].lincoef) )
3292  ++nquadlinterms;
3293  }
3294 
3295  SCIP_CALL( SCIPallocBufferArray(scip, &quadvars, nquadvars) );
3296  SCIP_CALL( SCIPallocBufferArray(scip, &quadelems, nquadelems) );
3297  SCIP_CALL( SCIPallocBufferArray(scip, &quadlinvars, nquadlinterms) );
3298  SCIP_CALL( SCIPallocBufferArray(scip, &quadlincoefs, nquadlinterms) );
3299 
3300  lincnt = 0;
3301  elcnt = 0;
3302  for( i = 0; i < nquadvars; ++i )
3303  {
3304  quadvars[i] = consdata->quadvarterms[i].var;
3305 
3306  if( consdata->quadvarterms[i].sqrcoef != 0.0 )
3307  {
3308  assert(elcnt < nquadelems);
3309  quadelems[elcnt].idx1 = i;
3310  quadelems[elcnt].idx2 = i;
3311  quadelems[elcnt].coef = consdata->quadvarterms[i].sqrcoef;
3312  ++elcnt;
3313  }
3314 
3315  if( !SCIPisZero(scip, consdata->quadvarterms[i].lincoef) )
3316  {
3317  assert(lincnt < nquadlinterms);
3318  quadlinvars [lincnt] = consdata->quadvarterms[i].var;
3319  quadlincoefs[lincnt] = consdata->quadvarterms[i].lincoef;
3320  ++lincnt;
3321  }
3322  }
3323  assert(lincnt == nquadlinterms);
3324 
3325  /* bilinear terms are sorted first by first variable, then by second variable
3326  * thus, it makes sense to remember the index of the previous first variable for the case a series of bilinear terms with the same first var appears */
3327  lastvar = NULL;
3328  lastvaridx = -1;
3329  for( i = 0; i < consdata->nbilinterms; ++i )
3330  {
3331  if( lastvar == consdata->bilinterms[i].var1 )
3332  {
3333  assert(lastvaridx >= 0);
3334  assert(consdata->quadvarterms[lastvaridx].var == consdata->bilinterms[i].var1);
3335  }
3336  else
3337  {
3338  lastvar = consdata->bilinterms[i].var1;
3339  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, lastvar, &lastvaridx) );
3340  }
3341  idx1 = lastvaridx;
3342 
3343  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, consdata->bilinterms[i].var2, &idx2) );
3344 
3345  assert(elcnt < nquadelems);
3346  quadelems[elcnt].idx1 = MIN(idx1, idx2);
3347  quadelems[elcnt].idx2 = MAX(idx1, idx2);
3348  quadelems[elcnt].coef = consdata->bilinterms[i].coef;
3349  ++elcnt;
3350  }
3351  assert(elcnt == nquadelems);
3352 
3353  /* set curvature for the nonlinear row */
3354  if( consdata->isconcave && consdata->isconvex )
3355  {
3356  assert(consdata->nbilinterms == 0 && consdata->nquadvars == 0);
3357  curvature = SCIP_EXPRCURV_LINEAR;
3358  }
3359  else if( consdata->isconcave )
3360  curvature = SCIP_EXPRCURV_CONCAVE;
3361  else if( consdata->isconvex )
3362  curvature = SCIP_EXPRCURV_CONVEX;
3363  else
3364  curvature = SCIP_EXPRCURV_UNKNOWN;
3365 
3366  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
3367  consdata->nlinvars, consdata->linvars, consdata->lincoefs,
3368  nquadvars, quadvars, nquadelems, quadelems,
3369  NULL, consdata->lhs, consdata->rhs,
3370  curvature) );
3371 
3372  SCIP_CALL( SCIPaddLinearCoefsToNlRow(scip, consdata->nlrow, nquadlinterms, quadlinvars, quadlincoefs) );
3373 
3374  SCIPfreeBufferArray(scip, &quadlincoefs);
3375  SCIPfreeBufferArray(scip, &quadlinvars);
3376  SCIPfreeBufferArray(scip, &quadelems);
3377  SCIPfreeBufferArray(scip, &quadvars);
3378 
3379  return SCIP_OKAY;
3380 }
3381 
3382 /** solve constraint as presolving */
3383 static
3384 SCIP_RETCODE presolveSolve(
3385  SCIP* scip, /**< SCIP data structure */
3386  SCIP_CONS* cons, /**< constraint */
3387  SCIP_RESULT* result, /**< to store result of solve: cutoff, success, or do-not-find */
3388  SCIP_Bool* redundant, /**< to store whether constraint is redundant now (should be deleted) */
3389  int* naggrvars /**< counter on number of variable aggregations */
3390  )
3391 {
3392  SCIP_CONSDATA* consdata;
3393 
3394  assert(scip != NULL);
3395  assert(cons != NULL);
3396  assert(result != NULL);
3397  assert(redundant != NULL);
3398 
3399  *result = SCIP_DIDNOTFIND;
3400  *redundant = FALSE;
3401 
3402  consdata = SCIPconsGetData(cons);
3403  assert(consdata != NULL);
3404 
3405  /* if constraint is an equality with two variables, at least one of them binary,
3406  * and linear after fixing the binary, then we can aggregate the variables */
3407  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) && consdata->nlinvars == 0 && consdata->nquadvars == 2 &&
3408  ((SCIPvarIsBinary(consdata->quadvarterms[0].var) && consdata->quadvarterms[1].sqrcoef == 0.0) ||
3409  (SCIPvarIsBinary(consdata->quadvarterms[1].var) && consdata->quadvarterms[0].sqrcoef == 0.0)) )
3410  {
3411  SCIP_Bool infeasible;
3412  SCIP_Bool aggregated;
3413  SCIP_Real a;
3414  SCIP_Real b;
3415  SCIP_Real c;
3416  SCIP_VAR* x;
3417  SCIP_VAR* y;
3418  int binvaridx;
3419 
3420  /* constraint is a*(x+x^2) + b*y + c*x*y = rhs, with x binary variable
3421  * x = 0 -> b*y == rhs
3422  * x = 1 -> (b+c)*y == rhs - a
3423  *
3424  * if b != 0 and b+c != 0, then y = (rhs-a)/(b+c) * x + rhs/b * (1-x) = ((rhs-a)/(b+c) - rhs/b) * x + rhs/b
3425  */
3426 
3427  binvaridx = (SCIPvarIsBinary(consdata->quadvarterms[0].var) && consdata->quadvarterms[1].sqrcoef == 0.0) ? 0 : 1;
3428 
3429  x = consdata->quadvarterms[binvaridx].var;
3430  a = consdata->quadvarterms[binvaridx].sqrcoef + consdata->quadvarterms[binvaridx].lincoef;
3431 
3432  y = consdata->quadvarterms[1-binvaridx].var;
3433  b = consdata->quadvarterms[1-binvaridx].lincoef;
3434 
3435  assert(consdata->nbilinterms <= 1); /* should actually be 1, since constraint is otherwise linear */
3436  c = (consdata->nbilinterms == 1) ? consdata->bilinterms[0].coef : 0.0;
3437 
3438  if( !SCIPisZero(scip, b) && !SCIPisZero(scip, b+c) )
3439  {
3440  SCIPdebugMsg(scip, "<%s> = 0 -> %g*<%s> = %g and <%s> = 1 -> %g*<%s> = %g\n", SCIPvarGetName(x), b, SCIPvarGetName(y), consdata->rhs,
3441  SCIPvarGetName(x), b+c, SCIPvarGetName(y), consdata->rhs - a);
3442  SCIPdebugMsg(scip, "=> attempt aggregation <%s> = %g*<%s> + %g\n", SCIPvarGetName(y), (consdata->rhs-a)/(b+c) - consdata->rhs/b,
3443  SCIPvarGetName(x), consdata->rhs/b);
3444 
3445  SCIP_CALL( SCIPaggregateVars(scip, x, y, (consdata->rhs-a)/(b+c) - consdata->rhs/b, -1.0, -consdata->rhs/b, &infeasible, redundant, &aggregated) );
3446  if( infeasible )
3447  *result = SCIP_CUTOFF;
3448  else if( *redundant || aggregated )
3449  {
3450  /* aggregated (or were already aggregated), so constraint is now redundant */
3451  *result = SCIP_SUCCESS;
3452  *redundant = TRUE;
3453 
3454  if( aggregated )
3455  ++*naggrvars;
3456  }
3457  }
3458 
3459  /* @todo if b is 0 or b+c is 0, or lhs != rhs, then could replace by varbound constraint */
3460  }
3461 
3462  return SCIP_OKAY;
3463 }
3464 
3465 
3466 /** reformulates products of binary variables as AND constraint
3467  *
3468  * For a product x*y, with x and y binary variables, the product is replaced by a new auxiliary variable z and the constraint z = {x and y} is added.
3469  */
3470 static
3471 SCIP_RETCODE presolveTryAddAND(
3472  SCIP* scip, /**< SCIP data structure */
3473  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3474  SCIP_CONS* cons, /**< constraint */
3475  int* naddconss /**< buffer where to add the number of AND constraints added */
3476  )
3477 {
3478  SCIP_CONSHDLRDATA* conshdlrdata;
3479  SCIP_CONSDATA* consdata;
3480  char name[SCIP_MAXSTRLEN];
3481  SCIP_VAR* vars[2];
3482  SCIP_VAR* auxvar;
3483  SCIP_CONS* andcons;
3484  int i;
3485  int ntodelete;
3486  int* todelete;
3487 
3488  assert(scip != NULL);
3489  assert(conshdlr != NULL);
3490  assert(cons != NULL);
3491  assert(naddconss != NULL);
3492 
3493  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3494  assert(conshdlrdata != NULL);
3495 
3496  /* if no binary variables, then we will find nothing to reformulate here
3497  * (note that this does not count in integer variables with {0,1} bounds...)
3498  */
3499  if( SCIPgetNBinVars(scip) == 0 )
3500  return SCIP_OKAY;
3501 
3502  /* if user does not like AND very much, then return */
3503  if( conshdlrdata->empathy4and < 2 )
3504  return SCIP_OKAY;
3505 
3506  consdata = SCIPconsGetData(cons);
3507  assert(consdata != NULL);
3508 
3509  if( consdata->nbilinterms == 0 )
3510  return SCIP_OKAY;
3511 
3512  /* get array to store indices of bilinear terms that shall be deleted */
3513  SCIP_CALL( SCIPallocBufferArray(scip, &todelete, consdata->nbilinterms) );
3514  ntodelete = 0;
3515 
3516  for( i = 0; i < consdata->nbilinterms; ++i )
3517  {
3518  vars[0] = consdata->bilinterms[i].var1;
3519  if( !SCIPvarIsBinary(vars[0]) )
3520  continue;
3521 
3522  vars[1] = consdata->bilinterms[i].var2;
3523  if( !SCIPvarIsBinary(vars[1]) )
3524  continue;
3525 
3526  /* create auxiliary variable */
3527  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_%s", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]), SCIPconsGetName(cons));
3528  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_BINARY,
3529  SCIPvarIsInitial(vars[0]) || SCIPvarIsInitial(vars[1]), SCIPvarIsRemovable(vars[0]) && SCIPvarIsRemovable(vars[1]), NULL, NULL, NULL, NULL, NULL) );
3530  SCIP_CALL( SCIPaddVar(scip, auxvar) );
3531 #ifdef WITH_DEBUG_SOLUTION
3532  if( SCIPdebugIsMainscip(scip) )
3533  {
3534  SCIP_Real var0val;
3535  SCIP_Real var1val;
3536  SCIP_CALL( SCIPdebugGetSolVal(scip, vars[0], &var0val) );
3537  SCIP_CALL( SCIPdebugGetSolVal(scip, vars[1], &var1val) );
3538  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, var0val * var1val) );
3539  }
3540 #endif
3541 
3542  /* create AND-constraint auxvar = x and y, need to be enforced as not redundant */
3543  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%sAND%s", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]));
3544  SCIP_CALL( SCIPcreateConsAnd(scip, &andcons, name, auxvar, 2, vars,
3545  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
3546  SCIPconsIsSeparated(cons), TRUE, TRUE,
3547  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
3548  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
3549  SCIP_CALL( SCIPaddCons(scip, andcons) );
3550  SCIPdebugMsg(scip, "added AND constraint: ");
3551  SCIPdebugPrintCons(scip, andcons, NULL);
3552  SCIP_CALL( SCIPreleaseCons(scip, &andcons) );
3553  ++*naddconss;
3554 
3555  /* add bilincoef * auxvar to linear terms */
3556  SCIP_CALL( addLinearCoef(scip, cons, auxvar, consdata->bilinterms[i].coef) );
3557  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
3558 
3559  /* remember that we have to delete this bilinear term */
3560  assert(ntodelete < consdata->nbilinterms);
3561  todelete[ntodelete++] = i;
3562  }
3563 
3564  /* remove bilinear terms that have been replaced */
3565  SCIP_CALL( removeBilinearTermsPos(scip, cons, ntodelete, todelete) );
3566  SCIPfreeBufferArray(scip, &todelete);
3567 
3568  return SCIP_OKAY;
3569 }
3570 
3571 /** gets bounds of variable y if x takes a certain value; checks whether x = xval has implications on y */
3572 static
3573 SCIP_RETCODE getImpliedBounds(
3574  SCIP* scip, /**< SCIP data structure */
3575  SCIP_VAR* x, /**< variable which implications to check */
3576  SCIP_Bool xval, /**< value of x to check for (TRUE for 1, FALSE for 0) */
3577  SCIP_VAR* y, /**< variable to check if bounds can be reduced */
3578  SCIP_INTERVAL* resultant /**< buffer to store bounds on y */
3579  )
3580 {
3581  SCIP_VAR** implvars;
3582  SCIP_BOUNDTYPE* impltypes;
3583  SCIP_Real* implbounds;
3584  int nimpls;
3585  int pos;
3586 
3587  assert(scip != NULL);
3588  assert(x != NULL);
3589  assert(y != NULL);
3590  assert(resultant != NULL);
3591 
3592  SCIPintervalSetBounds(resultant, MIN(SCIPvarGetLbGlobal(y), SCIPvarGetUbGlobal(y)), MAX(SCIPvarGetLbGlobal(y), SCIPvarGetUbGlobal(y))); /*lint !e666 */
3593 
3594  if( !SCIPvarIsBinary(x) || !SCIPvarIsActive(x) )
3595  return SCIP_OKAY;
3596 
3597  /* check in cliques for binary to binary implications */
3598  if( SCIPvarIsBinary(y) )
3599  {
3600  resultant->inf = MAX(resultant->inf, MIN(resultant->sup, 0.0));
3601  resultant->sup = MIN(resultant->sup, MAX(resultant->inf, 1.0));
3602 
3603  if( SCIPhaveVarsCommonClique(scip, x, xval, y, TRUE, FALSE) )
3604  {
3605  resultant->sup = MIN(resultant->sup, MAX(resultant->inf, 0.0));
3606  }
3607  else if( SCIPhaveVarsCommonClique(scip, x, xval, y, FALSE, FALSE) )
3608  {
3609  resultant->inf = MAX(resultant->inf, MIN(resultant->sup, 1.0));
3610  }
3611 
3612  return SCIP_OKAY;
3613  }
3614 
3615  /* analyze implications for x = xval */
3616  nimpls = SCIPvarGetNImpls(x, xval);
3617  if( nimpls == 0 )
3618  return SCIP_OKAY;
3619 
3620  implvars = SCIPvarGetImplVars (x, xval);
3621  impltypes = SCIPvarGetImplTypes (x, xval);
3622  implbounds = SCIPvarGetImplBounds(x, xval);
3623 
3624  assert(implvars != NULL);
3625  assert(impltypes != NULL);
3626  assert(implbounds != NULL);
3627 
3628  /* find implications */
3629  if( !SCIPsortedvecFindPtr((void**)implvars, SCIPvarComp, (void*)y, nimpls, &pos) )
3630  return SCIP_OKAY;
3631 
3632  /* if there are several implications on y, go to the first one */
3633  while( pos > 0 && implvars[pos-1] == y )
3634  --pos;
3635 
3636  /* update implied lower and upper bounds on y
3637  * but make sure that resultant will not be empty, due to tolerances
3638  */
3639  while( pos < nimpls && implvars[pos] == y )
3640  {
3641  if( impltypes[pos] == SCIP_BOUNDTYPE_LOWER )
3642  resultant->inf = MAX(resultant->inf, MIN(resultant->sup, implbounds[pos]));
3643  else
3644  resultant->sup = MIN(resultant->sup, MAX(resultant->inf, implbounds[pos]));
3645  ++pos;
3646  }
3647 
3648  assert(resultant->sup >= resultant->inf);
3649 
3650  return SCIP_OKAY;
3651 }
3652 
3653 /** Reformulates products of binary times bounded continuous variables as system of linear inequalities (plus auxiliary variable).
3654  *
3655  * For a product x*y, with y a binary variable and x a continous variable with finite bounds,
3656  * an auxiliary variable z and the inequalities \f$ x^L y \leq z \leq x^U y \f$ and \f$ x - (1-y) x^U \leq z \leq x - (1-y) x^L \f$ are added.
3657  *
3658  * If x is a linear term consisting of more than one variable, it is split up in groups of linear terms of length at most maxnrvar.
3659  * For each product of linear term of length at most maxnrvar with y, an auxiliary z and linear inequalities are added.
3660  *
3661  * If y is a binary variable, the AND constraint \f$ z = x \wedge y \f$ may be added instead of linear constraints.
3662  */
3663 static
3664 SCIP_RETCODE presolveTryAddLinearReform(
3665  SCIP* scip, /**< SCIP data structure */
3666  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3667  SCIP_CONS* cons, /**< constraint */
3668  int* naddconss /**< buffer where to add the number of auxiliary constraints added */
3669  )
3670 { /*lint --e{666} */
3671  SCIP_CONSHDLRDATA* conshdlrdata;
3672  SCIP_CONSDATA* consdata;
3673  SCIP_VAR** xvars;
3674  SCIP_Real* xcoef;
3675  SCIP_INTERVAL xbndszero;
3676  SCIP_INTERVAL xbndsone;
3677  SCIP_INTERVAL act0;
3678  SCIP_INTERVAL act1;
3679  int nxvars;
3680  SCIP_VAR* y;
3681  SCIP_VAR* bvar;
3682  char name[SCIP_MAXSTRLEN];
3683  int nbilinterms;
3684  SCIP_VAR* auxvar;
3685  SCIP_CONS* auxcons;
3686  int i;
3687  int j;
3688  int k;
3689  int bilinidx;
3690  SCIP_Real bilincoef;
3691  SCIP_Real mincoef;
3692  SCIP_Real maxcoef;
3693  int* todelete;
3694  int ntodelete;
3695  int maxnrvar;
3696  SCIP_Bool integral;
3697  SCIP_Longint gcd;
3698  SCIP_Bool auxvarinitial;
3699  SCIP_Bool auxvarremovable;
3700 
3701  assert(scip != NULL);
3702  assert(conshdlr != NULL);
3703  assert(cons != NULL);
3704  assert(naddconss != NULL);
3705 
3706  /* if no binary variables, then we will find nothing to reformulate here
3707  * (note that this does not count in integer variables with {0,1} bounds...)
3708  */
3709  if( SCIPgetNBinVars(scip) == 0 )
3710  return SCIP_OKAY;
3711 
3712  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3713  assert(conshdlrdata != NULL);
3714 
3715  maxnrvar = conshdlrdata->replacebinaryprodlength;
3716  if( maxnrvar == 0 )
3717  return SCIP_OKAY;
3718 
3719  consdata = SCIPconsGetData(cons);
3720  assert(consdata != NULL);
3721 
3722  xvars = NULL;
3723  xcoef = NULL;
3724  todelete = NULL;
3725  gcd = 0;
3726 
3727  for( i = 0; i < consdata->nquadvars; ++i )
3728  {
3729  y = consdata->quadvarterms[i].var;
3730  if( !SCIPvarIsBinary(y) )
3731  continue;
3732 
3733  nbilinterms = consdata->quadvarterms[i].nadjbilin;
3734  if( nbilinterms == 0 )
3735  continue;
3736 
3737  SCIP_CALL( SCIPreallocBufferArray(scip, &xvars, MIN(maxnrvar, nbilinterms)+2) ); /* add 2 for later use when creating linear constraints */
3738  SCIP_CALL( SCIPreallocBufferArray(scip, &xcoef, MIN(maxnrvar, nbilinterms)+2) );
3739 
3740  /* alloc array to store indices of bilinear terms that shall be deleted */
3741  SCIP_CALL( SCIPreallocBufferArray(scip, &todelete, nbilinterms) );
3742  ntodelete = 0;
3743 
3744  auxvarinitial = SCIPvarIsInitial(y);
3745  auxvarremovable = SCIPvarIsRemovable(y);
3746 
3747  /* setup a list of bounded variables x_i with coefficients a_i that are multiplied with binary y: y*(sum_i a_i*x_i)
3748  * and compute range of sum_i a_i*x_i for the cases y = 0 and y = 1
3749  * we may need several rounds if maxnrvar < nbilinterms
3750  */
3751  j = 0;
3752  do
3753  {
3754  nxvars = 0;
3755  SCIPintervalSet(&xbndszero, 0.0);
3756  SCIPintervalSet(&xbndsone, 0.0);
3757 
3758  mincoef = SCIPinfinity(scip);
3759  maxcoef = 0.0;
3760  integral = TRUE;
3761 
3762  /* collect at most maxnrvar variables for x term */
3763  for( ; j < nbilinterms && nxvars < maxnrvar; ++j )
3764  {
3765  bilinidx = consdata->quadvarterms[i].adjbilin[j];
3766  assert(bilinidx >= 0);
3767  assert(bilinidx < consdata->nbilinterms);
3768 
3769  bvar = consdata->bilinterms[bilinidx].var1;
3770  if( bvar == y )
3771  bvar = consdata->bilinterms[bilinidx].var2;
3772  assert(bvar != y);
3773 
3774  /* skip products with unbounded variables */
3775  if( SCIPisInfinity(scip, -SCIPvarGetLbGlobal(bvar)) || SCIPisInfinity(scip, SCIPvarGetUbGlobal(bvar)) )
3776  {
3777  SCIPdebugMsg(scip, "skip reform of <%s><%s> due to unbounded second variable [%g,%g]\n",
3778  SCIPvarGetName(y), SCIPvarGetName(bvar), SCIPvarGetLbGlobal(bvar), SCIPvarGetUbGlobal(bvar));
3779  continue;
3780  }
3781 
3782  /* skip products with non-binary variables if binreformbinaryonly is set */
3783  if( conshdlrdata->binreformbinaryonly && !SCIPvarIsBinary(bvar) )
3784  {
3785  SCIPdebugMsg(scip, "skip reform of <%s><%s> because second variable is not binary\n",
3786  SCIPvarGetName(y), SCIPvarGetName(bvar));
3787  continue;
3788  }
3789 
3790  bilincoef = consdata->bilinterms[bilinidx].coef;
3791  assert(bilincoef != 0.0);
3792 
3793  /* get activity of bilincoef * x if y = 0 */
3794  SCIP_CALL( getImpliedBounds(scip, y, FALSE, bvar, &act0) );
3795  SCIPintervalMulScalar(SCIPinfinity(scip), &act0, act0, bilincoef);
3796 
3797  /* get activity of bilincoef * x if y = 1 */
3798  SCIP_CALL( getImpliedBounds(scip, y, TRUE, bvar, &act1) );
3799  SCIPintervalMulScalar(SCIPinfinity(scip), &act1, act1, bilincoef);
3800 
3801  /* skip products that give rise to very large coefficients (big big-M's) */
3802  if( SCIPfeastol(scip) * REALABS(act0.inf) >= conshdlrdata->binreformmaxcoef || SCIPfeastol(scip) * REALABS(act0.sup) >= conshdlrdata->binreformmaxcoef )
3803  {
3804  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity [%g,%g] for <%s> = 0.0\n",
3805  bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar), SCIPintervalGetInf(act0), SCIPintervalGetSup(act0), SCIPvarGetName(y));
3806  continue;
3807  }
3808  if( SCIPfeastol(scip) * REALABS(act1.inf) >= conshdlrdata->binreformmaxcoef || SCIPfeastol(scip) * REALABS(act1.sup) >= conshdlrdata->binreformmaxcoef )
3809  {
3810  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity [%g,%g] for <%s> = 1.0\n",
3811  bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar), SCIPintervalGetInf(act1), SCIPintervalGetSup(act1), SCIPvarGetName(y));
3812  continue;
3813  }
3814  if( !SCIPisZero(scip, MIN(REALABS(act0.inf), REALABS(act0.sup))) &&
3815  SCIPfeastol(scip) * MAX(REALABS(act0.inf), REALABS(act0.sup)) / MIN(REALABS(act0.inf), REALABS(act0.sup)) >= conshdlrdata->binreformmaxcoef )
3816  {
3817  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity ratio %g for <%s> = 0.0\n", bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar),
3818  MAX(REALABS(act0.inf), REALABS(act0.sup)) / MIN(REALABS(act0.inf), REALABS(act0.sup)), SCIPvarGetName(y));
3819  continue;
3820  }
3821  if( !SCIPisZero(scip, MIN(REALABS(act1.inf), REALABS(act1.sup))) &&
3822  SCIPfeastol(scip) * MAX(REALABS(act1.inf), REALABS(act1.sup)) / MIN(REALABS(act1.inf), REALABS(act1.sup)) >= conshdlrdata->binreformmaxcoef )
3823  {
3824  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity ratio %g for <%s> = 0.0\n", bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar),
3825  MAX(REALABS(act1.inf), REALABS(act1.sup)) / MIN(REALABS(act1.inf), REALABS(act1.sup)), SCIPvarGetName(y));
3826  continue;
3827  }
3828 
3829  /* add bvar to x term */
3830  xvars[nxvars] = bvar;
3831  xcoef[nxvars] = bilincoef;
3832  ++nxvars;
3833 
3834  /* update bounds on x term */
3835  SCIPintervalAdd(SCIPinfinity(scip), &xbndszero, xbndszero, act0);
3836  SCIPintervalAdd(SCIPinfinity(scip), &xbndsone, xbndsone, act1);
3837 
3838  if( REALABS(bilincoef) < mincoef )
3839  mincoef = ABS(bilincoef);
3840  if( REALABS(bilincoef) > maxcoef )
3841  maxcoef = ABS(bilincoef);
3842 
3843  /* update whether all coefficients will be integral and if so, compute their gcd */
3844  integral &= (SCIPvarGetType(bvar) < SCIP_VARTYPE_CONTINUOUS) && SCIPisIntegral(scip, bilincoef); /*lint !e514 */
3845  if( integral )
3846  {
3847  if( nxvars == 1 )
3848  gcd = (SCIP_Longint)SCIPround(scip, REALABS(bilincoef));
3849  else
3850  gcd = SCIPcalcGreComDiv(gcd, (SCIP_Longint)SCIPround(scip, REALABS(bilincoef)));
3851  }
3852 
3853  /* if bvar is initial, then also the auxiliary variable should be initial
3854  * if bvar is not removable, then also the auxiliary variable should not be removable
3855  */
3856  auxvarinitial |= SCIPvarIsInitial(bvar);
3857  auxvarremovable &= SCIPvarIsRemovable(bvar);
3858 
3859  /* remember that we have to remove this bilinear term later */
3860  assert(ntodelete < nbilinterms);
3861  todelete[ntodelete++] = bilinidx;
3862  }
3863 
3864  if( nxvars == 0 ) /* all (remaining) x_j seem to be unbounded */
3865  break;
3866 
3867  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(xbndszero)));
3868  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(xbndszero)));
3869  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(xbndsone)));
3870  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(xbndsone)));
3871 
3872 #ifdef SCIP_DEBUG
3873  if( SCIPintervalGetInf(xbndszero) != SCIPintervalGetInf(xbndsone) || /*lint !e777*/
3874  +SCIPintervalGetSup(xbndszero) != SCIPintervalGetSup(xbndsone) ) /*lint !e777*/
3875  {
3876  SCIPdebugMsg(scip, "got different bounds for y = 0: [%g, %g] and y = 1: [%g, %g]\n", xbndszero.inf, xbndszero.sup, xbndsone.inf, xbndsone.sup);
3877  }
3878 #endif
3879 
3880  if( nxvars == 1 && conshdlrdata->empathy4and >= 1 && SCIPvarIsBinary(xvars[0]) )
3881  {
3882  /* product of two binary variables, replace by auxvar and AND constraint */
3883  /* add auxiliary variable z */
3884  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
3885  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_IMPLINT,
3886  auxvarinitial, auxvarremovable, NULL, NULL, NULL, NULL, NULL) );
3887  SCIP_CALL( SCIPaddVar(scip, auxvar) );
3888 
3889 #ifdef WITH_DEBUG_SOLUTION
3890  if( SCIPdebugIsMainscip(scip) )
3891  {
3892  SCIP_Real var0val;
3893  SCIP_Real var1val;
3894  SCIP_CALL( SCIPdebugGetSolVal(scip, xvars[0], &var0val) );
3895  SCIP_CALL( SCIPdebugGetSolVal(scip, y, &var1val) );
3896  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, var0val * var1val) );
3897  }
3898 #endif
3899 
3900  /* add constraint z = x and y; need to be enforced, as it is not redundant w.r.t. existing constraints */
3901  xvars[1] = y;
3902  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%sAND%s_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
3903  SCIP_CALL( SCIPcreateConsAnd(scip, &auxcons, name, auxvar, 2, xvars,
3904  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
3905  SCIPconsIsSeparated(cons), TRUE, TRUE,
3906  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
3907  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
3908  SCIP_CALL( SCIPaddCons(scip, auxcons) );
3909  SCIPdebugMsg(scip, "added AND constraint: ");
3910  SCIPdebugPrintCons(scip, auxcons, NULL);
3911  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
3912  ++*naddconss;
3913 
3914  /* add linear term coef*auxvar */
3915  SCIP_CALL( addLinearCoef(scip, cons, auxvar, xcoef[0]) );
3916 
3917  /* forget about auxvar */
3918  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
3919  }
3920  else
3921  {
3922  /* product of binary variable with more than one binary or with continuous variables or with binary and user
3923  * did not like AND -> replace by auxvar and linear constraints */
3924  SCIP_Real scale;
3925 
3926  /* scale auxiliary constraint by some nice value,
3927  * if all coefficients are integral, take a value that preserves integrality (-> gcd), so we can make the auxiliary variable impl. integer
3928  */
3929  if( integral )
3930  {
3931  scale = (SCIP_Real)gcd;
3932  assert(scale >= 1.0);
3933  }
3934  else if( nxvars == 1 )
3935  {
3936  /* scaling by the only coefficient gives auxiliary variable = x * y, which thus will be implicit integral provided y is not continuous */
3937  assert(mincoef == maxcoef); /*lint !e777 */
3938  scale = mincoef;
3939  integral = SCIPvarGetType(xvars[0]) < SCIP_VARTYPE_CONTINUOUS;
3940  }
3941  else
3942  {
3943  scale = 1.0;
3944  if( maxcoef < 0.5 )
3945  scale = maxcoef;
3946  if( mincoef > 2.0 )
3947  scale = mincoef;
3948  if( scale != 1.0 )
3949  scale = SCIPselectSimpleValue(scale / 2.0, 1.5 * scale, MAXDNOM);
3950  }
3951  assert(scale > 0.0);
3952  assert(!SCIPisInfinity(scip, scale));
3953 
3954  /* if x-term is always negative for y = 1, negate scale so we get a positive auxiliary variable; maybe this is better sometimes? */
3955  if( !SCIPisPositive(scip, SCIPintervalGetSup(xbndsone)) )
3956  scale = -scale;
3957 
3958  SCIPdebugMsg(scip, "binary reformulation using scale %g, nxvars = %d, integral = %u\n", scale, nxvars, integral);
3959  if( scale != 1.0 )
3960  {
3961  SCIPintervalDivScalar(SCIPinfinity(scip), &xbndszero, xbndszero, scale);
3962  SCIPintervalDivScalar(SCIPinfinity(scip), &xbndsone, xbndsone, scale);
3963  for( k = 0; k < nxvars; ++k )
3964  xcoef[k] /= scale;
3965  }
3966 
3967  /* add auxiliary variable z */
3968  if( nxvars == 1 )
3969  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
3970  else
3971  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_more_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
3972  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, MIN(0., SCIPintervalGetInf(xbndsone)), MAX(0., SCIPintervalGetSup(xbndsone)),
3973  0.0, integral ? SCIP_VARTYPE_IMPLINT : SCIP_VARTYPE_CONTINUOUS,
3974  auxvarinitial, auxvarremovable, NULL, NULL, NULL, NULL, NULL) );
3975  SCIP_CALL( SCIPaddVar(scip, auxvar) );
3976 
3977  /* compute value of auxvar in debug solution */
3978 #ifdef WITH_DEBUG_SOLUTION
3979  if( SCIPdebugIsMainscip(scip) )
3980  {
3981  SCIP_Real debugval;
3982  SCIP_Real varval;
3983 
3984  SCIP_CALL( SCIPdebugGetSolVal(scip, y, &varval) );
3985  if( SCIPisZero(scip, varval) )
3986  {
3987  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, 0.0) );
3988  }
3989  else
3990  {
3991  assert(SCIPisEQ(scip, varval, 1.0));
3992 
3993  debugval = 0.0;
3994  for( k = 0; k < nxvars; ++k )
3995  {
3996  SCIP_CALL( SCIPdebugGetSolVal(scip, xvars[k], &varval) );
3997  debugval += xcoef[k] * varval;
3998  }
3999  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
4000  }
4001  }
4002 #endif
4003 
4004  /* add auxiliary constraints
4005  * it seems to be advantageous to make the varbound constraints initial and the linear constraints not initial
4006  * maybe because it is more likely that a binary variable takes value 0 instead of 1, and thus the varbound constraints
4007  * are more often active, compared to the linear constraints added below
4008  * also, the varbound constraints are more sparse than the linear cons
4009  */
4010  if( SCIPisNegative(scip, SCIPintervalGetInf(xbndsone)) )
4011  {
4012  /* add 0 <= z - xbndsone.inf * y constraint (as varbound constraint), need to be enforced as not redundant */
4013  if( nxvars == 1 )
4014  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_1", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4015  else
4016  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_1", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4017  SCIP_CALL( SCIPcreateConsVarbound(scip, &auxcons, name, auxvar, y, -SCIPintervalGetInf(xbndsone), 0.0, SCIPinfinity(scip),
4018  SCIPconsIsInitial(cons) /*&& conshdlrdata->binreforminitial*/,
4019  SCIPconsIsSeparated(cons), TRUE, TRUE,
4020  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4021  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4022  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4023  SCIPdebugMsg(scip, "added varbound constraint: ");
4024  SCIPdebugPrintCons(scip, auxcons, NULL);
4025  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4026  ++*naddconss;
4027  }
4028  if( SCIPisPositive(scip, SCIPintervalGetSup(xbndsone)) )
4029  {
4030  /* add z - xbndsone.sup * y <= 0 constraint (as varbound constraint), need to be enforced as not redundant */
4031  if( nxvars == 1 )
4032  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_2", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4033  else
4034  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_2", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4035  SCIP_CALL( SCIPcreateConsVarbound(scip, &auxcons, name, auxvar, y, -SCIPintervalGetSup(xbndsone), -SCIPinfinity(scip), 0.0,
4036  SCIPconsIsInitial(cons) /*&& conshdlrdata->binreforminitial*/,
4037  SCIPconsIsSeparated(cons), TRUE, TRUE,
4038  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4039  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4040  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4041  SCIPdebugMsg(scip, "added varbound constraint: ");
4042  SCIPdebugPrintCons(scip, auxcons, NULL);
4043  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4044  ++*naddconss;
4045  }
4046 
4047  /* add xbndszero.inf <= sum_i a_i*x_i + xbndszero.inf * y - z constraint, need to be enforced as not redundant */
4048  xvars[nxvars] = y;
4049  xvars[nxvars+1] = auxvar;
4050  xcoef[nxvars] = SCIPintervalGetInf(xbndszero);
4051  xcoef[nxvars+1] = -1;
4052 
4053  if( nxvars == 1 )
4054  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_3", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4055  else
4056  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_3", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4057  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, name, nxvars+2, xvars, xcoef, SCIPintervalGetInf(xbndszero), SCIPinfinity(scip),
4058  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
4059  SCIPconsIsSeparated(cons), TRUE, TRUE,
4060  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4061  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4062  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4063  SCIPdebugMsg(scip, "added linear constraint: ");
4064  SCIPdebugPrintCons(scip, auxcons, NULL);
4065  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4066  ++*naddconss;
4067 
4068  /* add sum_i a_i*x_i + xbndszero.sup * y - z <= xbndszero.sup constraint, need to be enforced as not redundant */
4069  xcoef[nxvars] = SCIPintervalGetSup(xbndszero);
4070 
4071  if( nxvars == 1 )
4072  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_4", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4073  else
4074  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_4", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4075  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, name, nxvars+2, xvars, xcoef, -SCIPinfinity(scip), SCIPintervalGetSup(xbndszero),
4076  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
4077  SCIPconsIsSeparated(cons), TRUE, TRUE,
4078  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4079  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4080  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4081  SCIPdebugMsg(scip, "added linear constraint: ");
4082  SCIPdebugPrintCons(scip, auxcons, NULL);
4083  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4084  ++*naddconss;
4085 
4086  /* add linear term scale*auxvar to this constraint */
4087  SCIP_CALL( addLinearCoef(scip, cons, auxvar, scale) );
4088 
4089  /* forget about auxvar */
4090  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4091  }
4092  }
4093  while( j < nbilinterms );
4094 
4095  /* remove bilinear terms that have been replaced */
4096  SCIP_CALL( removeBilinearTermsPos(scip, cons, ntodelete, todelete) );
4097  }
4098  SCIPdebugMsg(scip, "resulting quadratic constraint: ");
4099  SCIPdebugPrintCons(scip, cons, NULL);
4100 
4101  SCIPfreeBufferArrayNull(scip, &todelete);
4102  SCIPfreeBufferArrayNull(scip, &xcoef);
4103  SCIPfreeBufferArrayNull(scip, &xvars);
4104 
4105  return SCIP_OKAY;
4106 }
4107 
4108 /** tries to automatically convert a quadratic constraint (or a part of it) into a more specific and more specialized constraint */
4109 static
4110 SCIP_RETCODE presolveUpgrade(
4111  SCIP* scip, /**< SCIP data structure */
4112  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
4113  SCIP_CONS* cons, /**< source constraint to try to convert */
4114  SCIP_Bool* upgraded, /**< buffer to store whether constraint was upgraded */
4115  int* nupgdconss, /**< buffer to increase if constraint was upgraded */
4116  int* naddconss, /**< buffer to increase with number of additional constraints created during upgrade */
4117  SCIP_PRESOLTIMING presoltiming /**< current presolving timing */
4118  )
4119 {
4120  SCIP_CONSHDLRDATA* conshdlrdata;
4121  SCIP_CONSDATA* consdata;
4122  SCIP_VAR* var;
4123  SCIP_Real lincoef;
4124  SCIP_Real quadcoef;
4125  SCIP_Real lb;
4126  SCIP_Real ub;
4127  int nbinlin;
4128  int nbinquad;
4129  int nintlin;
4130  int nintquad;
4131  int nimpllin;
4132  int nimplquad;
4133  int ncontlin;
4134  int ncontquad;
4135  SCIP_Bool integral;
4136  int i;
4137  int j;
4138  SCIP_CONS** upgdconss;
4139  int upgdconsssize;
4140  int nupgdconss_;
4141 
4142  assert(scip != NULL);
4143  assert(conshdlr != NULL);
4144  assert(cons != NULL);
4145  assert(!SCIPconsIsModifiable(cons));
4146  assert(upgraded != NULL);
4147  assert(nupgdconss != NULL);
4148  assert(naddconss != NULL);
4149 
4150  *upgraded = FALSE;
4151 
4152  nupgdconss_ = 0;
4153 
4154  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4155  assert(conshdlrdata != NULL);
4156 
4157  /* if there are no upgrade methods, we can also stop */
4158  if( conshdlrdata->nquadconsupgrades == 0 )
4159  return SCIP_OKAY;
4160 
4161  upgdconsssize = 2;
4162  SCIP_CALL( SCIPallocBufferArray(scip, &upgdconss, upgdconsssize) );
4163 
4164  consdata = SCIPconsGetData(cons);
4165  assert(consdata != NULL);
4166 
4167  /* calculate some statistics on quadratic constraint */
4168  nbinlin = 0;
4169  nbinquad = 0;
4170  nintlin = 0;
4171  nintquad = 0;
4172  nimpllin = 0;
4173  nimplquad = 0;
4174  ncontlin = 0;
4175  ncontquad = 0;
4176  integral = TRUE;
4177  for( i = 0; i < consdata->nlinvars; ++i )
4178  {
4179  var = consdata->linvars[i];
4180  lincoef = consdata->lincoefs[i];
4181  lb = SCIPvarGetLbLocal(var);
4182  ub = SCIPvarGetUbLocal(var);
4183  assert(!SCIPisZero(scip, lincoef));
4184 
4185  switch( SCIPvarGetType(var) )
4186  {
4187  case SCIP_VARTYPE_BINARY:
4188  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4189  integral = integral && SCIPisIntegral(scip, lincoef);
4190  nbinlin++;
4191  break;
4192  case SCIP_VARTYPE_INTEGER:
4193  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4194  integral = integral && SCIPisIntegral(scip, lincoef);
4195  nintlin++;
4196  break;
4197  case SCIP_VARTYPE_IMPLINT:
4198  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4199  integral = integral && SCIPisIntegral(scip, lincoef);
4200  nimpllin++;
4201  break;
4202  case SCIP_VARTYPE_CONTINUOUS:
4203  integral = integral && SCIPisRelEQ(scip, lb, ub) && SCIPisIntegral(scip, lincoef * lb);
4204  ncontlin++;
4205  break;
4206  default:
4207  SCIPerrorMessage("unknown variable type\n");
4208  return SCIP_INVALIDDATA;
4209  }
4210  }
4211 
4212  for( i = 0; i < consdata->nquadvars; ++i )
4213  {
4214  var = consdata->quadvarterms[i].var;
4215  lincoef = consdata->quadvarterms[i].lincoef;
4216  quadcoef = consdata->quadvarterms[i].sqrcoef;
4217  lb = SCIPvarGetLbLocal(var);
4218  ub = SCIPvarGetUbLocal(var);
4219 
4220  switch( SCIPvarGetType(var) )
4221  {
4222  case SCIP_VARTYPE_BINARY:
4223  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4224  integral = integral && SCIPisIntegral(scip, lincoef) && SCIPisIntegral(scip, quadcoef);
4225  nbinquad++;
4226  break;
4227  case SCIP_VARTYPE_INTEGER:
4228  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4229  integral = integral && SCIPisIntegral(scip, lincoef) && SCIPisIntegral(scip, quadcoef);
4230  nintquad++;
4231  break;
4232  case SCIP_VARTYPE_IMPLINT:
4233  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4234  integral = integral && SCIPisIntegral(scip, lincoef) && SCIPisIntegral(scip, quadcoef);
4235  nimplquad++;
4236  break;
4237  case SCIP_VARTYPE_CONTINUOUS:
4238  integral = integral && SCIPisRelEQ(scip, lb, ub) && SCIPisIntegral(scip, lincoef * lb + quadcoef * lb * lb);
4239  ncontquad++;
4240  break;
4241  default:
4242  SCIPerrorMessage("unknown variable type\n");
4243  return SCIP_INVALIDDATA;
4244  }
4245  }
4246 
4247  if( integral )
4248  {
4249  for( i = 0; i < consdata->nbilinterms && integral; ++i )
4250  {
4251  if( SCIPvarGetType(consdata->bilinterms[i].var1) < SCIP_VARTYPE_CONTINUOUS && SCIPvarGetType(consdata->bilinterms[i].var2) < SCIP_VARTYPE_CONTINUOUS )
4252  integral = integral && SCIPisIntegral(scip, consdata->bilinterms[i].coef);
4253  else
4254  integral = FALSE;
4255  }
4256  }
4257 
4258  /* call the upgrading methods */
4259 
4260  SCIPdebugMsg(scip, "upgrading quadratic constraint <%s> (%d upgrade methods):\n",
4261  SCIPconsGetName(cons), conshdlrdata->nquadconsupgrades);
4262  SCIPdebugMsg(scip, " binlin=%d binquad=%d intlin=%d intquad=%d impllin=%d implquad=%d contlin=%d contquad=%d integral=%u\n",
4263  nbinlin, nbinquad, nintlin, nintquad, nimpllin, nimplquad, ncontlin, ncontquad, integral);
4264  SCIPdebugPrintCons(scip, cons, NULL);
4265 
4266  /* try all upgrading methods in priority order in case the upgrading step is enable */
4267  for( i = 0; i < conshdlrdata->nquadconsupgrades; ++i )
4268  {
4269  if( !conshdlrdata->quadconsupgrades[i]->active )
4270  continue;
4271 
4272  SCIP_CALL( conshdlrdata->quadconsupgrades[i]->quadconsupgd(scip, cons,
4273  nbinlin, nbinquad, nintlin, nintquad, nimpllin, nimplquad, ncontlin, ncontquad, integral,
4274  &nupgdconss_, upgdconss, upgdconsssize, presoltiming) );
4275 
4276  while( nupgdconss_ < 0 )
4277  {
4278  /* upgrade function requires more memory: resize upgdconss and call again */
4279  assert(-nupgdconss_ > upgdconsssize);
4280  upgdconsssize = -nupgdconss_;
4281  SCIP_CALL( SCIPreallocBufferArray(scip, &upgdconss, -nupgdconss_) );
4282 
4283  SCIP_CALL( conshdlrdata->quadconsupgrades[i]->quadconsupgd(scip, cons,
4284  nbinlin, nbinquad, nintlin, nintquad, nimpllin, nimplquad, ncontlin, ncontquad, integral,
4285  &nupgdconss_, upgdconss, upgdconsssize, presoltiming) );
4286 
4287  assert(nupgdconss_ != 0);
4288  }
4289 
4290  if( nupgdconss_ > 0 )
4291  {
4292  /* got upgrade */
4293  SCIPdebugPrintCons(scip, cons, NULL);
4294  SCIPdebugMsg(scip, " -> upgraded to %d constraints:\n", nupgdconss_);
4295 
4296  /* add the upgraded constraints to the problem and forget them */
4297  for( j = 0; j < nupgdconss_; ++j )
4298  {
4299  SCIPdebugMsgPrint(scip, "\t");
4300  SCIPdebugPrintCons(scip, upgdconss[j], NULL);
4301 
4302  SCIP_CALL( SCIPaddCons(scip, upgdconss[j]) ); /*lint !e613*/
4303  SCIP_CALL( SCIPreleaseCons(scip, &upgdconss[j]) ); /*lint !e613*/
4304  }
4305 
4306  /* count the first upgrade constraint as constraint upgrade and the remaining ones as added constraints */
4307  *nupgdconss += 1;
4308  *naddconss += nupgdconss_ - 1;
4309  *upgraded = TRUE;
4310 
4311  /* delete upgraded constraint */
4312  SCIPdebugMsg(scip, "delete constraint <%s> after upgrade\n", SCIPconsGetName(cons));
4313  SCIP_CALL( SCIPdelCons(scip, cons) );
4314 
4315  break;
4316  }
4317  }
4318 
4319  SCIPfreeBufferArray(scip, &upgdconss);
4320 
4321  return SCIP_OKAY;
4322 }
4323 
4324 /** helper function for presolveDisaggregate */
4325 static
4326 SCIP_RETCODE presolveDisaggregateMarkComponent(
4327  SCIP* scip, /**< SCIP data structure */
4328  SCIP_CONSDATA* consdata, /**< constraint data */
4329  int quadvaridx, /**< index of quadratic variable to mark */
4330  SCIP_HASHMAP* var2component, /**< variables to components mapping */
4331  int componentnr, /**< the component number to mark to */
4332  int* componentsize /**< buffer to store size of component (incremented by 1) */
4333  )
4334 {
4335  SCIP_QUADVARTERM* quadvarterm;
4336  SCIP_VAR* othervar;
4337  int othervaridx;
4338  int i;
4339 
4340  assert(consdata != NULL);
4341  assert(quadvaridx >= 0);
4342  assert(quadvaridx < consdata->nquadvars);
4343  assert(var2component != NULL);
4344  assert(componentnr >= 0);
4345 
4346  quadvarterm = &consdata->quadvarterms[quadvaridx];
4347 
4348  if( SCIPhashmapExists(var2component, quadvarterm->var) )
4349  {
4350  /* if we saw the variable before, then it should have the same component number */
4351  assert(SCIPhashmapGetImageInt(var2component, quadvarterm->var) == componentnr);
4352  return SCIP_OKAY;
4353  }
4354 
4355  /* assign component number to variable */
4356  SCIP_CALL( SCIPhashmapInsertInt(var2component, quadvarterm->var, componentnr) );
4357  ++*componentsize;
4358 
4359  /* assign same component number to all variables this variable is multiplied with */
4360  for( i = 0; i < quadvarterm->nadjbilin; ++i )
4361  {
4362  othervar = consdata->bilinterms[quadvarterm->adjbilin[i]].var1 == quadvarterm->var ?
4363  consdata->bilinterms[quadvarterm->adjbilin[i]].var2 : consdata->bilinterms[quadvarterm->adjbilin[i]].var1;
4364  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, othervar, &othervaridx) );
4365  assert(othervaridx >= 0);
4366  SCIP_CALL( presolveDisaggregateMarkComponent(scip, consdata, othervaridx, var2component, componentnr, componentsize) );
4367  }
4368 
4369  return SCIP_OKAY;
4370 }
4371 
4372 /** merges components in variables connectivity graph */
4373 static
4374 SCIP_RETCODE presolveDisaggregateMergeComponents(
4375  SCIP* scip, /**< SCIP data structure */
4376  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
4377  SCIP_HASHMAP* var2component, /**< variables to component mapping */
4378  int nvars, /**< number of variables */
4379  int* ncomponents, /**< number of components */
4380  int* componentssize /**< size of components */
4381  )
4382 {
4383  SCIP_CONSHDLRDATA* conshdlrdata;
4384  SCIP_HASHMAPENTRY* entry;
4385  int maxncomponents;
4386  int* oldcompidx;
4387  int* newcompidx;
4388  int i;
4389  int oldcomponent;
4390  int newcomponent;
4391 
4392  assert(scip != NULL);
4393  assert(conshdlr != NULL);
4394  assert(var2component != NULL);
4395  assert(ncomponents != NULL);
4396  assert(componentssize != NULL);
4397 
4398  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4399  assert(conshdlrdata != NULL);
4400 
4401  maxncomponents = conshdlrdata->maxdisaggrsize;
4402  assert(maxncomponents > 0);
4403 
4404  /* if already not too many components, then nothing to do */
4405  if( *ncomponents <= maxncomponents )
4406  return SCIP_OKAY;
4407 
4408  /*
4409  printf("component sizes before:");
4410  for( i = 0; i < *ncomponents; ++i )
4411  printf(" %d", componentssize[i]);
4412  printf("\n");
4413  */
4414 
4415  SCIP_CALL( SCIPallocBufferArray(scip, &oldcompidx, *ncomponents) );
4416  SCIP_CALL( SCIPallocBufferArray(scip, &newcompidx, *ncomponents) );
4417 
4418  for( i = 0; i < *ncomponents; ++i )
4419  oldcompidx[i] = i;
4420 
4421  switch( conshdlrdata->disaggrmergemethod )
4422  {
4423  case 's' :
4424  /* sort components by size, increasing order */
4425  SCIPsortIntInt(componentssize, oldcompidx, *ncomponents);
4426  break;
4427  case 'b' :
4428  case 'm' :
4429  /* sort components by size, decreasing order */
4430  SCIPsortDownIntInt(componentssize, oldcompidx, *ncomponents);
4431  break;
4432  default :
4433  SCIPerrorMessage("invalid value for constraints/quadratic/disaggrmergemethod parameter");
4434  return SCIP_PARAMETERWRONGVAL;
4435  }
4436 
4437  SCIPdebugMsg(scip, "%-30s: % 4d components of size % 4d to % 4d, median: % 4d\n", SCIPgetProbName(scip), *ncomponents, componentssize[0], componentssize[*ncomponents-1], componentssize[*ncomponents/2]);
4438 
4439  if( conshdlrdata->disaggrmergemethod == 'm' )
4440  {
4441  SCIP_Real targetsize;
4442  int count = 0;
4443 
4444  /* a minimal component size we should reach to have all components roughly the same size */
4445  targetsize = nvars / maxncomponents; /*lint !e653*/
4446  for( i = 0; i < *ncomponents; ++i )
4447  {
4448  newcompidx[oldcompidx[i]] = i;
4449  count += componentssize[i];
4450 
4451  /* fill with small components until we reach targetsize
4452  * Since targetsize might be fractional, we also add another component if
4453  * the number of variables remaining (=nvars-count) is larger than
4454  * what we expect to put into the remaining components (=targetsize * (maxncomponents - i-1)).
4455  * Thus, from time to time, a component is made larger than the targetsize to avoid
4456  * having to add much into the last component.
4457  */
4458  while( i < *ncomponents-1 && (componentssize[i] + componentssize[*ncomponents-1] <= targetsize ||
4459  nvars - count > targetsize * (maxncomponents - i)) )
4460  {
4461  /* map last (=smallest) component to component i */
4462  newcompidx[oldcompidx[*ncomponents-1]] = i;
4463 
4464  /* increase size of component i accordingly */
4465  componentssize[i] += componentssize[*ncomponents-1];
4466  count += componentssize[*ncomponents-1];
4467 
4468  /* forget about last component */
4469  --*ncomponents;
4470  }
4471  }
4472  assert(count == nvars);
4473  }
4474  else
4475  {
4476  /* get inverse permutation */
4477  for( i = 0; i < *ncomponents; ++i )
4478  newcompidx[oldcompidx[i]] = i;
4479  }
4480 
4481  /* assign new component numbers to variables, cutting off at maxncomponents */
4482  for( i = 0; i < SCIPhashmapGetNEntries(var2component); ++i )
4483  {
4484  entry = SCIPhashmapGetEntry(var2component, i);
4485  if( entry == NULL )
4486  continue;
4487 
4488  oldcomponent = (int)(size_t)SCIPhashmapEntryGetImage(entry);
4489 
4490  newcomponent = newcompidx[oldcomponent];
4491  if( newcomponent >= maxncomponents )
4492  {
4493  newcomponent = maxncomponents-1;
4494  ++componentssize[maxncomponents-1];
4495  }
4496 
4497  SCIPhashmapEntrySetImage(entry, (void*)(size_t)newcomponent); /*lint !e571*/
4498  }
4499  if( *ncomponents > maxncomponents )
4500  *ncomponents = maxncomponents;
4501 
4502  /*
4503  printf("component sizes after :");
4504  for( i = 0; i < *ncomponents; ++i )
4505  printf(" %d", componentssize[i]);
4506  printf("\n");
4507  */
4508 
4509  SCIPfreeBufferArray(scip, &newcompidx);
4510  SCIPfreeBufferArray(scip, &oldcompidx);
4511 
4512  return SCIP_OKAY;
4513 }
4514 
4515 /** compute the next highest power of 2 for a 32-bit argument
4516  *
4517  * Source: https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
4518  *
4519  * @note Returns 0 for v=0.
4520  */
4521 static
4522 unsigned int nextPowerOf2(
4523  unsigned int v /**< input */
4524  )
4525 {
4526  v--;
4527  v |= v >> 1;
4528  v |= v >> 2;
4529  v |= v >> 4;
4530  v |= v >> 8;
4531  v |= v >> 16;
4532  v++;
4533 
4534  return v;
4535 }
4536 
4537 
4538 /** for quadratic constraints that consists of a sum of quadratic terms, disaggregates the sum into a set of constraints by introducing auxiliary variables
4539  *
4540  * Assume the quadratic constraint can be written in the form
4541  * lhs <= b'x + sum_{k=1..p} q_k(x_k) <= rhs
4542  * where x_k denotes a subset of the variables in x and these subsets are pairwise disjunct
4543  * and q_k(.) is a quadratic form.
4544  * p is selected as large as possible, but to be <= conshdlrdata->maxdisaggrsize.
4545  *
4546  * Without additional scaling, the constraint is disaggregated into
4547  * lhs <= b'x + sum_k c_k z_k <= rhs
4548  * c_k z_k ~ q_k(x)
4549  * where "~" is either "<=", "==", or ">=", depending on whether lhs or rhs are infinite.
4550  * Further, c_k is chosen to be the maximal absolute value of the coefficients of the quadratic terms in q_k(x).
4551  * This is done to ensure that z_k takes values with a similar magnitute as the variables in x_k (better for separation).
4552  *
4553  * However, a solution of this disaggregated system can violate the original constraint by (p+1)*epsilon
4554  * (assuming unscaled violations are used, which is the default).
4555  * Therefore, all constraints are scaled by p+1:
4556  * (p+1)*lhs <= (p+1)*b'x + (p+1) * sum_k c_k z_k <= (p+1) * rhs
4557  * (p+1)*c_k z_k ~ (p+1)*q_k(x)
4558  */
4559 static
4560 SCIP_RETCODE presolveDisaggregate(
4561  SCIP* scip, /**< SCIP data structure */
4562  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
4563  SCIP_CONS* cons, /**< source constraint to try to convert */
4564  int* naddconss /**< pointer to counter of added constraints */
4565  )
4566 {
4567  SCIP_CONSDATA* consdata;
4568  SCIP_HASHMAP* var2component;
4569  int* componentssize;
4570  int ncomponents;
4571  int i;
4572  int comp;
4573  SCIP_CONS** auxconss;
4574  SCIP_VAR** auxvars;
4575  SCIP_Real* auxcoefs;
4576 #ifdef WITH_DEBUG_SOLUTION
4577  SCIP_Real* auxsolvals; /* value of auxiliary variable in debug solution */
4578 #endif
4579  SCIP_Real scale;
4580  char name[SCIP_MAXSTRLEN];
4581 
4582  assert(scip != NULL);
4583  assert(conshdlr != NULL);
4584  assert(cons != NULL);
4585  assert(naddconss != NULL);
4586 
4587  consdata = SCIPconsGetData(cons);
4588  assert(consdata != NULL);
4589 
4590  /* skip if constraint has been already disaggregated */
4591  if( consdata->isdisaggregated )
4592  return SCIP_OKAY;
4593 
4594  consdata->isdisaggregated = TRUE;
4595 
4596  /* make sure there are no quadratic variables without coefficients */
4597  SCIP_CALL( mergeAndCleanBilinearTerms(scip, cons) );
4598  SCIP_CALL( mergeAndCleanQuadVarTerms(scip, cons) );
4599 
4600  if( consdata->nquadvars <= 1 )
4601  return SCIP_OKAY;
4602 
4603  /* sort quadratic variable terms here, so we can later search in it without reordering the array */
4604  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
4605 
4606  /* check how many quadratic terms with non-overlapping variables we have
4607  * in other words, the number of components in the sparsity graph of the quadratic term matrix
4608  */
4609  ncomponents = 0;
4610  SCIP_CALL( SCIPhashmapCreate(&var2component, SCIPblkmem(scip), consdata->nquadvars) );
4611  SCIP_CALL( SCIPallocBufferArray(scip, &componentssize, consdata->nquadvars) );
4612  for( i = 0; i < consdata->nquadvars; ++i )
4613  {
4614  /* if variable was marked already, skip it */
4615  if( SCIPhashmapExists(var2component, (void*)consdata->quadvarterms[i].var) )
4616  continue;
4617 
4618  /* start a new component with variable i */
4619  componentssize[ncomponents] = 0;
4620  SCIP_CALL( presolveDisaggregateMarkComponent(scip, consdata, i, var2component, ncomponents, componentssize + ncomponents) );
4621  ++ncomponents;
4622  }
4623 
4624  assert(ncomponents >= 1);
4625 
4626  /* if there is only one component, we cannot disaggregate
4627  * @todo we could still split the constraint into several while keeping the number of variables sharing several constraints as small as possible
4628  */
4629  if( ncomponents == 1 )
4630  {
4631  SCIPhashmapFree(&var2component);
4632  SCIPfreeBufferArray(scip, &componentssize);
4633  return SCIP_OKAY;
4634  }
4635 
4636  /* merge some components, if necessary */
4637  SCIP_CALL( presolveDisaggregateMergeComponents(scip, conshdlr, var2component, consdata->nquadvars, &ncomponents, componentssize) );
4638 
4639  SCIPfreeBufferArray(scip, &componentssize);
4640 
4641  /* scale all new constraints (ncomponents+1 many) by ncomponents+1 (or its next power of 2), so violations sum up to at most epsilon */
4642  scale = nextPowerOf2((unsigned int)ncomponents + 1);
4643 
4644  SCIP_CALL( SCIPallocBufferArray(scip, &auxconss, ncomponents) );
4645  SCIP_CALL( SCIPallocBufferArray(scip, &auxvars, ncomponents) );
4646  SCIP_CALL( SCIPallocBufferArray(scip, &auxcoefs, ncomponents) );
4647 #ifdef WITH_DEBUG_SOLUTION
4648  SCIP_CALL( SCIPallocClearBufferArray(scip, &auxsolvals, ncomponents) );
4649 #endif
4650 
4651  /* create auxiliary variables and empty constraints for each component */
4652  for( comp = 0; comp < ncomponents; ++comp )
4653  {
4654  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_comp%d", SCIPconsGetName(cons), comp);
4655 
4656  SCIP_CALL( SCIPcreateVar(scip, &auxvars[comp], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
4657  SCIP_VARTYPE_CONTINUOUS, SCIPconsIsInitial(cons), FALSE, NULL, NULL, NULL, NULL, NULL) );
4658 
4659  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &auxconss[comp], name, 0, NULL, NULL, 0, NULL, 0, NULL,
4660  (SCIPisInfinity(scip, -consdata->lhs) ? -SCIPinfinity(scip) : 0.0),
4661  (SCIPisInfinity(scip, consdata->rhs) ? SCIPinfinity(scip) : 0.0),
4662  SCIPconsIsInitial(cons), SCIPconsIsSeparated(cons), TRUE,
4663  TRUE, SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4664  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons)) );
4665 
4666  auxcoefs[comp] = SCIPinfinity(scip);
4667  }
4668 
4669  /* add quadratic variables to each component constraint
4670  * delete adjacency information */
4671  for( i = 0; i < consdata->nquadvars; ++i )
4672  {
4673  assert(SCIPhashmapExists(var2component, consdata->quadvarterms[i].var));
4674 
4675  comp = SCIPhashmapGetImageInt(var2component, consdata->quadvarterms[i].var);
4676  assert(comp >= 0);
4677  assert(comp < ncomponents);
4678 
4679  /* add variable term to corresponding constraint */
4680  SCIP_CALL( SCIPaddQuadVarQuadratic(scip, auxconss[comp], consdata->quadvarterms[i].var, scale * consdata->quadvarterms[i].lincoef, scale * consdata->quadvarterms[i].sqrcoef) );
4681 
4682  /* reduce coefficient of aux variable */
4683  if( !SCIPisZero(scip, consdata->quadvarterms[i].lincoef) && ABS(consdata->quadvarterms[i].lincoef) < auxcoefs[comp] )
4684  auxcoefs[comp] = REALABS(consdata->quadvarterms[i].lincoef);
4685  if( !SCIPisZero(scip, consdata->quadvarterms[i].sqrcoef) && ABS(consdata->quadvarterms[i].sqrcoef) < auxcoefs[comp] )
4686  auxcoefs[comp] = REALABS(consdata->quadvarterms[i].sqrcoef);
4687 
4688  SCIPfreeBlockMemoryArray(scip, &consdata->quadvarterms[i].adjbilin, consdata->quadvarterms[i].adjbilinsize);
4689  consdata->quadvarterms[i].nadjbilin = 0;
4690  consdata->quadvarterms[i].adjbilinsize = 0;
4691 
4692 #ifdef WITH_DEBUG_SOLUTION
4693  if( SCIPdebugIsMainscip(scip) )
4694  {
4695  SCIP_Real debugvarval;
4696 
4697  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata->quadvarterms[i].var, &debugvarval) );
4698  auxsolvals[comp] += consdata->quadvarterms[i].lincoef * debugvarval + consdata->quadvarterms[i].sqrcoef * debugvarval * debugvarval;
4699  }
4700 #endif
4701  }
4702 
4703  /* add bilinear terms to each component constraint */
4704  for( i = 0; i < consdata->nbilinterms; ++i )
4705  {
4706  assert(SCIPhashmapExists(var2component, consdata->bilinterms[i].var1));
4707  assert(SCIPhashmapExists(var2component, consdata->bilinterms[i].var2));
4708 
4709  comp = SCIPhashmapGetImageInt(var2component, consdata->bilinterms[i].var1);
4710  assert(comp == SCIPhashmapGetImageInt(var2component, consdata->bilinterms[i].var2));
4711  assert(!SCIPisZero(scip, consdata->bilinterms[i].coef));
4712 
4713  SCIP_CALL( SCIPaddBilinTermQuadratic(scip, auxconss[comp],
4714  consdata->bilinterms[i].var1, consdata->bilinterms[i].var2, scale * consdata->bilinterms[i].coef) );
4715 
4716  if( ABS(consdata->bilinterms[i].coef) < auxcoefs[comp] )
4717  auxcoefs[comp] = ABS(consdata->bilinterms[i].coef);
4718 
4719 #ifdef WITH_DEBUG_SOLUTION
4720  if( SCIPdebugIsMainscip(scip) )
4721  {
4722  SCIP_Real debugvarval1;
4723  SCIP_Real debugvarval2;
4724 
4725  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata->bilinterms[i].var1, &debugvarval1) );
4726  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata->bilinterms[i].var2, &debugvarval2) );
4727  auxsolvals[comp] += consdata->bilinterms[i].coef * debugvarval1 * debugvarval2;
4728  }
4729 #endif
4730  }
4731 
4732  /* forget about bilinear terms in cons */
4733  SCIPfreeBlockMemoryArray(scip, &consdata->bilinterms, consdata->bilintermssize);
4734  consdata->nbilinterms = 0;
4735  consdata->bilintermssize = 0;
4736 
4737  /* remove quadratic variable terms from cons */
4738  for( i = consdata->nquadvars - 1; i >= 0; --i )
4739  {
4740  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
4741  }
4742  assert(consdata->nquadvars == 0);
4743 
4744  /* scale remaining linear variables and sides by scale */
4745  for( i = 0; i < consdata->nlinvars; ++i )
4746  {
4747  SCIP_CALL( chgLinearCoefPos(scip, cons, i, scale * consdata->lincoefs[i]) );
4748  }
4749  if( !SCIPisInfinity(scip, -consdata->lhs) )
4750  {
4751  consdata->lhs *= scale;
4752  assert(!SCIPisInfinity(scip, -consdata->lhs) );
4753  }
4754  if( !SCIPisInfinity(scip, consdata->rhs) )
4755  {
4756  consdata->rhs *= scale;
4757  assert(!SCIPisInfinity(scip, consdata->rhs) );
4758  }
4759 
4760  /* add auxiliary variables to auxiliary constraints
4761  * add aux vars and constraints to SCIP
4762  * add aux vars to this constraint
4763  * set value of aux vars in debug solution, if any
4764  */
4765  SCIPdebugMsg(scip, "add %d constraints for disaggregation of quadratic constraint <%s>\n", ncomponents, SCIPconsGetName(cons));
4766  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars + ncomponents) );
4767  for( comp = 0; comp < ncomponents; ++comp )
4768  {
4769  SCIP_CONSDATA* auxconsdata;
4770 
4771  SCIP_CALL( SCIPaddLinearVarQuadratic(scip, auxconss[comp], auxvars[comp], -scale * auxcoefs[comp]) );
4772 
4773  SCIP_CALL( SCIPaddVar(scip, auxvars[comp]) );
4774 
4775  SCIP_CALL( SCIPaddCons(scip, auxconss[comp]) );
4776  SCIPdebugPrintCons(scip, auxconss[comp], NULL);
4777 
4778  SCIP_CALL( addLinearCoef(scip, cons, auxvars[comp], scale * auxcoefs[comp]) );
4779 
4780  /* mark that the constraint should not further be disaggregated */
4781  auxconsdata = SCIPconsGetData(auxconss[comp]);
4782  assert(auxconsdata != NULL);
4783  auxconsdata->isdisaggregated = TRUE;
4784 
4785 #ifdef WITH_DEBUG_SOLUTION
4786  if( SCIPdebugIsMainscip(scip) )
4787  {
4788  /* auxvar should take value from auxsolvals in debug solution, but we also scaled auxvar by auxcoefs[comp] */
4789  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvars[comp], auxsolvals[comp] / auxcoefs[comp]) );
4790  }
4791 #endif
4792 
4793  SCIP_CALL( SCIPreleaseCons(scip, &auxconss[comp]) );
4794  SCIP_CALL( SCIPreleaseVar(scip, &auxvars[comp]) );
4795  }
4796  *naddconss += ncomponents;
4797 
4798  SCIPdebugPrintCons(scip, cons, NULL);
4799 
4800  SCIPfreeBufferArray(scip, &auxconss);
4801  SCIPfreeBufferArray(scip, &auxvars);
4802  SCIPfreeBufferArray(scip, &auxcoefs);
4803 #ifdef WITH_DEBUG_SOLUTION
4804  SCIPfreeBufferArray(scip, &auxsolvals);
4805 #endif
4806  SCIPhashmapFree(&var2component);
4807 
4808  return SCIP_OKAY;
4809 }
4810 
4811 #ifdef CHECKIMPLINBILINEAR
4812 /** checks if there are bilinear terms x*y with a binary variable x and an implication x = {0,1} -> y = 0
4813  *
4814  * In this case, the bilinear term can be removed (x=0 case) or replaced by y (x=1 case).
4815  */
4816 static
4817 SCIP_RETCODE presolveApplyImplications(
4818  SCIP* scip, /**< SCIP data structure */
4819  SCIP_CONS* cons, /**< quadratic constraint */
4820  int* nbilinremoved /**< buffer to store number of removed bilinear terms */
4821  )
4822 {
4823  SCIP_CONSDATA* consdata;
4824  SCIP_VAR* x;
4825  SCIP_VAR* y;
4826  SCIP_INTERVAL implbnds;
4827  int i;
4828  int j;
4829  int k;
4830 
4831  assert(scip != NULL);
4832  assert(cons != NULL);
4833  assert(nbilinremoved != NULL);
4834 
4835  *nbilinremoved = 0;
4836 
4837  consdata = SCIPconsGetData(cons);
4838  assert(consdata != NULL);
4839 
4840  SCIPdebugMsg(scip, "apply implications in <%s>\n", SCIPconsGetName(cons));
4841 
4842  /* sort quadvarterms in case we need to search */
4843  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
4844 
4845  for( i = 0; i < consdata->nquadvars; ++i )
4846  {
4847  x = consdata->quadvarterms[i].var;
4848  assert(x != NULL);
4849 
4850  if( consdata->quadvarterms[i].nadjbilin == 0 )
4851  continue;
4852 
4853  if( !SCIPvarIsBinary(x) )
4854  continue;
4855 
4856  if( !SCIPvarIsActive(x) )
4857  continue;
4858 
4859  if( SCIPvarGetNImpls(x, TRUE) == 0 && SCIPvarGetNImpls(x, FALSE) == 0 )
4860  continue;
4861 
4862  for( j = 0; j < consdata->quadvarterms[i].nadjbilin; ++j )
4863  {
4864  k = consdata->quadvarterms[i].adjbilin[j];
4865  assert(k >= 0);
4866  assert(k < consdata->nbilinterms);
4867 
4868  if( consdata->bilinterms[k].coef == 0.0 )
4869  continue;
4870 
4871  y = consdata->bilinterms[k].var1 == x ? consdata->bilinterms[k].var2 : consdata->bilinterms[k].var1;
4872  assert(x != y);
4873 
4874  SCIP_CALL( getImpliedBounds(scip, x, TRUE, y, &implbnds) );
4875  if( SCIPisZero(scip, implbnds.inf) && SCIPisZero(scip, implbnds.sup) )
4876  {
4877  /* if x = 1 implies y = 0, then we can remove the bilinear term x*y, since it is always 0
4878  * we only set the coefficient to 0.0 here and mark the bilinterms as not merged */
4879  SCIPdebugMsg(scip, "remove bilinear term %g<%s><%s> from <%s> due to implication\n", consdata->bilinterms[k].coef, SCIPvarGetName(x), SCIPvarGetName(y), SCIPconsGetName(cons));
4880  consdata->bilinterms[k].coef = 0.0;
4881  consdata->bilinmerged = FALSE;
4882  ++*nbilinremoved;
4883  continue;
4884  }
4885 
4886  SCIP_CALL( getImpliedBounds(scip, x, FALSE, y, &implbnds) );
4887  if( SCIPisZero(scip, implbnds.inf) && SCIPisZero(scip, implbnds.sup) )
4888  {
4889  /* if x = 0 implies y = 0, then we can replace the bilinear term x*y by y
4890  * we only move the coefficient to the linear coef of y here and mark the bilinterms as not merged */
4891  SCIPdebugMsg(scip, "replace bilinear term %g<%s><%s> by %g<%s> in <%s> due to implication\n", consdata->bilinterms[k].coef, SCIPvarGetName(x), SCIPvarGetName(y), consdata->bilinterms[k].coef, SCIPvarGetName(y), SCIPconsGetName(cons));
4892  assert(consdata->quadvarssorted);
4893  SCIP_CALL( SCIPaddQuadVarLinearCoefQuadratic(scip, cons, y, consdata->bilinterms[k].coef) );
4894  consdata->bilinterms[k].coef = 0.0;
4895  consdata->bilinmerged = FALSE;
4896  ++*nbilinremoved;
4897  }
4898  }
4899  }
4900 
4901  if( *nbilinremoved > 0 )
4902  {
4903  SCIP_CALL( mergeAndCleanBilinearTerms(scip, cons) );
4904 
4905  /* invalidate nonlinear row */
4906  if( consdata->nlrow != NULL )
4907  {
4908  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
4909  }
4910 
4911  consdata->ispropagated = FALSE;
4912  consdata->ispresolved = FALSE;
4913  consdata->iscurvchecked = FALSE;
4914  }
4915 
4916  consdata->isimpladded = FALSE;
4917 
4918  return SCIP_OKAY;
4919 }
4920 #endif
4921 
4922 /** checks a quadratic constraint for convexity and/or concavity without checking multivariate functions */
4923 static
4924 void checkCurvatureEasy(
4925  SCIP* scip, /**< SCIP data structure */
4926  SCIP_CONS* cons, /**< quadratic constraint */
4927  SCIP_Bool* determined, /**< pointer to store whether the curvature could be determined */
4928  SCIP_Bool checkmultivariate /**< whether curvature will be checked later on for multivariate functions */
4929  )
4930 {
4931  SCIP_CONSDATA* consdata;
4932  int nquadvars;
4933 
4934  assert(scip != NULL);
4935  assert(cons != NULL);
4936  assert(determined != NULL);
4937 
4938  consdata = SCIPconsGetData(cons);
4939  assert(consdata != NULL);
4940 
4941  nquadvars = consdata->nquadvars;
4942  *determined = TRUE;
4943 
4944  if( consdata->iscurvchecked )
4945  return;
4946 
4947  SCIPdebugMsg(scip, "Checking curvature of constraint <%s> without multivariate functions\n", SCIPconsGetName(cons));
4948 
4949  consdata->maxnonconvexity = 0.0;
4950  if( nquadvars == 1 )
4951  {
4952  assert(consdata->nbilinterms == 0);
4953  consdata->isconvex = !SCIPisNegative(scip, consdata->quadvarterms[0].sqrcoef);
4954  consdata->isconcave = !SCIPisPositive(scip, consdata->quadvarterms[0].sqrcoef);
4955  consdata->iscurvchecked = TRUE;
4956 
4957  if( !SCIPisInfinity(scip, -consdata->lhs) && consdata->quadvarterms[0].sqrcoef > 0.0 )
4958  consdata->maxnonconvexity = consdata->quadvarterms[0].sqrcoef;
4959  if( !SCIPisInfinity(scip, consdata->rhs) && consdata->quadvarterms[0].sqrcoef < 0.0 )
4960  consdata->maxnonconvexity = -consdata->quadvarterms[0].sqrcoef;
4961  }
4962  else if( nquadvars == 0 )
4963  {
4964  consdata->isconvex = TRUE;
4965  consdata->isconcave = TRUE;
4966  consdata->iscurvchecked = TRUE;
4967  }
4968  else if( consdata->nbilinterms == 0 )
4969  {
4970  int v;
4971 
4972  consdata->isconvex = TRUE;
4973  consdata->isconcave = TRUE;
4974 
4975  for( v = nquadvars - 1; v >= 0; --v )
4976  {
4977  consdata->isconvex = consdata->isconvex && !SCIPisNegative(scip, consdata->quadvarterms[v].sqrcoef);
4978  consdata->isconcave = consdata->isconcave && !SCIPisPositive(scip, consdata->quadvarterms[v].sqrcoef);
4979 
4980  if( !SCIPisInfinity(scip, -consdata->lhs) && consdata->quadvarterms[v].sqrcoef > consdata->maxnonconvexity )
4981  consdata->maxnonconvexity = consdata->quadvarterms[0].sqrcoef;
4982  if( !SCIPisInfinity(scip, consdata->rhs) && -consdata->quadvarterms[v].sqrcoef > consdata->maxnonconvexity )
4983  consdata->maxnonconvexity = -consdata->quadvarterms[0].sqrcoef;
4984  }
4985 
4986  consdata->iscurvchecked = TRUE;
4987  }
4988  else if( !checkmultivariate )
4989  {
4990  consdata->isconvex = FALSE;
4991  consdata->isconcave = FALSE;
4992  consdata->iscurvchecked = TRUE;
4993  consdata->maxnonconvexity = SCIPinfinity(scip);
4994  }
4995  else
4996  *determined = FALSE;
4997 }
4998 
4999 /** checks a quadratic constraint for convexity and/or concavity */
5000 static
5001 SCIP_RETCODE checkCurvature(
5002  SCIP* scip, /**< SCIP data structure */
5003  SCIP_CONS* cons, /**< quadratic constraint */
5004  SCIP_Bool checkmultivariate /**< whether curvature should also be checked for multivariate functions */
5005  )
5006 {
5007  SCIP_CONSDATA* consdata;
5008  double* matrix;
5009  SCIP_HASHMAP* var2index;
5010  int i;
5011  int n;
5012  int nn;
5013  int row;
5014  int col;
5015  double* alleigval;
5016  SCIP_Bool determined;
5017 
5018  assert(scip != NULL);
5019  assert(cons != NULL);
5020 
5021  consdata = SCIPconsGetData(cons);
5022  assert(consdata != NULL);
5023 
5024  n = consdata->nquadvars;
5025 
5026  if( consdata->iscurvchecked )
5027  return SCIP_OKAY;
5028 
5029  /* easy checks for curvature detection */
5030  checkCurvatureEasy(scip, cons, &determined, checkmultivariate);
5031 
5032  /* if curvature was already detected stop */
5033  if( determined )
5034  {
5035  return SCIP_OKAY;
5036  }
5037 
5038  SCIPdebugMsg(scip, "Checking curvature of constraint <%s> with multivariate functions\n", SCIPconsGetName(cons));
5039 
5040  if( n == 2 )
5041  {
5042  SCIP_Real tracehalf;
5043  SCIP_Real discriminantroot;
5044 
5045  /* compute eigenvalues by hand */
5046  assert(consdata->nbilinterms == 1);
5047  consdata->isconvex =
5048  consdata->quadvarterms[0].sqrcoef >= 0 &&
5049  consdata->quadvarterms[1].sqrcoef >= 0 &&
5050  4 * consdata->quadvarterms[0].sqrcoef * consdata->quadvarterms[1].sqrcoef >= consdata->bilinterms[0].coef * consdata->bilinterms[0].coef;
5051  consdata->isconcave =
5052  consdata->quadvarterms[0].sqrcoef <= 0 &&
5053  consdata->quadvarterms[1].sqrcoef <= 0 &&
5054  4 * consdata->quadvarterms[0].sqrcoef * consdata->quadvarterms[1].sqrcoef >= consdata->bilinterms[0].coef * consdata->bilinterms[0].coef;
5055 
5056  /* store largest eigenvalue causing nonconvexity according to sides */
5057  tracehalf = (consdata->quadvarterms[0].sqrcoef + consdata->quadvarterms[1].sqrcoef) / 2.0;
5058  discriminantroot = consdata->quadvarterms[0].sqrcoef * consdata->quadvarterms[1].sqrcoef - SQR(consdata->bilinterms[0].coef / 2.0);
5059  discriminantroot = SQR(tracehalf) - discriminantroot;
5060  assert(!SCIPisNegative(scip, discriminantroot));
5061  discriminantroot = SQRT(MAX(0.0, discriminantroot));
5062 
5063  consdata->maxnonconvexity = 0.0;
5064  if( !SCIPisInfinity(scip, -consdata->lhs) )
5065  consdata->maxnonconvexity = MAX(consdata->maxnonconvexity, tracehalf + discriminantroot);
5066  if( !SCIPisInfinity(scip, consdata->rhs) )
5067  consdata->maxnonconvexity = MAX(consdata->maxnonconvexity, discriminantroot - tracehalf);
5068 
5069  consdata->iscurvchecked = TRUE;
5070  return SCIP_OKAY;
5071  }
5072 
5073  /* do not check curvature if n is too large */
5074  nn = n * n;
5075  if( nn < 0 || (unsigned) (int) nn > UINT_MAX / sizeof(SCIP_Real) )
5076  {
5077  SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "cons_quadratic - n is too large to check the curvature\n");
5078  consdata->isconvex = FALSE;
5079  consdata->isconcave = FALSE;
5080  consdata->iscurvchecked = TRUE;
5081  consdata->maxnonconvexity = SCIPinfinity(scip);
5082  return SCIP_OKAY;
5083  }
5084 
5085  /* lower triangular of quadratic term matrix */
5086  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, nn) );
5087  BMSclearMemoryArray(matrix, nn);
5088 
5089  consdata->isconvex = TRUE;
5090  consdata->isconcave = TRUE;
5091  consdata->maxnonconvexity = 0.0;
5092 
5093  SCIP_CALL( SCIPhashmapCreate(&var2index, SCIPblkmem(scip), n) );
5094  for( i = 0; i < n; ++i )
5095  {
5096  if( consdata->quadvarterms[i].nadjbilin > 0 )
5097  {
5098  SCIP_CALL( SCIPhashmapInsertInt(var2index, consdata->quadvarterms[i].var, i) );
5099  matrix[i*n + i] = consdata->quadvarterms[i].sqrcoef;
5100  }
5101  else
5102  {
5103  /* if pure square term, then update maximal nonconvex eigenvalue, as it will not be considered in lapack call below */
5104  if( !SCIPisInfinity(scip, -consdata->lhs) && consdata->quadvarterms[i].sqrcoef > consdata->maxnonconvexity )
5105  consdata->maxnonconvexity = consdata->quadvarterms[i].sqrcoef;
5106  if( !SCIPisInfinity(scip, consdata->rhs) && -consdata->quadvarterms[i].sqrcoef > consdata->maxnonconvexity )
5107  consdata->maxnonconvexity = -consdata->quadvarterms[i].sqrcoef;
5108  }
5109  /* nonzero elements on diagonal tell a lot about convexity/concavity */
5110  if( SCIPisNegative(scip, consdata->quadvarterms[i].sqrcoef) )
5111  consdata->isconvex = FALSE;
5112  if( SCIPisPositive(scip, consdata->quadvarterms[i].sqrcoef) )
5113  consdata->isconcave = FALSE;
5114  }
5115 
5116  /* skip lapack call, if we know already that we are indefinite
5117  * NOTE: this will leave out updating consdata->maxnonconvexity, so that it only provides a lower bound in this case
5118  */
5119  if( !consdata->isconvex && !consdata->isconcave )
5120  {
5121  SCIPfreeBufferArray(scip, &matrix);
5122  SCIPhashmapFree(&var2index);
5123  consdata->iscurvchecked = TRUE;
5124  /* make sure that maxnonconvexity is strictly different from zero if nonconvex
5125  * TODO one could think about doing some eigenvalue estimation here (Gershgorin)
5126  */
5127  consdata->maxnonconvexity = MAX(1000.0, consdata->maxnonconvexity);
5128  return SCIP_OKAY;
5129  }
5130 
5131  if( SCIPisIpoptAvailableIpopt() )
5132  {
5133  for( i = 0; i < consdata->nbilinterms; ++i )
5134  {
5135  assert(SCIPhashmapExists(var2index, consdata->bilinterms[i].var1));
5136  assert(SCIPhashmapExists(var2index, consdata->bilinterms[i].var2));
5137  row = SCIPhashmapGetImageInt(var2index, consdata->bilinterms[i].var1);
5138  col = SCIPhashmapGetImageInt(var2index, consdata->bilinterms[i].var2);
5139  if( row < col )
5140  matrix[row * n + col] = consdata->bilinterms[i].coef/2;
5141  else
5142  matrix[col * n + row] = consdata->bilinterms[i].coef/2;
5143  }
5144 
5145  SCIP_CALL( SCIPallocBufferArray(scip, &alleigval, n) );
5146  /* @todo Can we compute only min and max eigen value?
5147  * @todo Can we estimate the numerical error?
5148  * @todo Trying a cholesky factorization may be much faster.
5149  */
5150  if( LapackDsyev(FALSE, n, matrix, alleigval) != SCIP_OKAY )
5151  {
5152  SCIPwarningMessage(scip, "Failed to compute eigenvalues of quadratic coefficient matrix of constraint %s. Assuming matrix is indefinite.\n", SCIPconsGetName(cons));
5153  consdata->isconvex = FALSE;
5154  consdata->isconcave = FALSE;
5155  }
5156  else
5157  {
5158  /* deconvexification reformulates a stricly convex quadratic function in binaries such that it becomes not-strictly convex
5159  * by adding the -lambda*(x^2-x) terms for lambda the smallest eigenvalue of the matrix
5160  * the result is still a convex form "but less so" (ref. papers by Guignard et.al.), but with hopefully tighter value for the continuous relaxation
5161  */
5162 #ifdef DECONVEXIFY
5163  SCIP_Bool allbinary;
5164  printf("cons <%s>[%g,%g] spectrum = [%g,%g]\n", SCIPconsGetName(cons), consdata->lhs, consdata->rhs, alleigval[0], alleigval[n-1]);
5165 #endif
5166  consdata->isconvex &= !SCIPisNegative(scip, alleigval[0]); /*lint !e514*/
5167  consdata->isconcave &= !SCIPisPositive(scip, alleigval[n-1]); /*lint !e514*/
5168  consdata->iscurvchecked = TRUE;
5169 #ifdef DECONVEXIFY
5170  for( i = 0; i < consdata->nquadvars; ++i )
5171  if( !SCIPvarIsBinary(consdata->quadvarterms[i].var) )
5172  break;
5173  allbinary = i == consdata->nquadvars;
5174 
5175  if( !SCIPisInfinity(scip, consdata->rhs) && alleigval[0] > 0.1 && allbinary )
5176  {
5177  printf("deconvexify cons <%s> by shifting hessian by %g\n", SCIPconsGetName(cons), alleigval[0]);
5178  for( i = 0; i < consdata->nquadvars; ++i )
5179  {
5180  consdata->quadvarterms[i].sqrcoef -= alleigval[0];
5181  consdata->quadvarterms[i].lincoef += alleigval[0];
5182  }
5183  }
5184 
5185  if( !SCIPisInfinity(scip, consdata->lhs) && alleigval[n-1] < -0.1 && allbinary )
5186  {
5187  printf("deconcavify cons <%s> by shifting hessian by %g\n", SCIPconsGetName(cons), alleigval[n-1]);
5188  for( i = 0; i < consdata->nquadvars; ++i )
5189  {
5190  consdata->quadvarterms[i].sqrcoef -= alleigval[n-1];
5191  consdata->quadvarterms[i].lincoef += alleigval[n-1];
5192  }
5193  }
5194 #endif
5195  }
5196 
5197  /* update largest eigenvalue causing nonconvexity according to sides */
5198  if( !SCIPisInfinity(scip, -consdata->lhs) )
5199  consdata->maxnonconvexity = MAX(consdata->maxnonconvexity, alleigval[n-1]);
5200  if( !SCIPisInfinity(scip, consdata->rhs) )
5201  consdata->maxnonconvexity = MAX(consdata->maxnonconvexity, -alleigval[0]);
5202 
5203  SCIPfreeBufferArray(scip, &alleigval);
5204  }
5205  else
5206  {
5207  consdata->isconvex = FALSE;
5208  consdata->isconcave = FALSE;
5209  consdata->iscurvchecked = TRUE; /* set to TRUE since it does not help to repeat this procedure again and again (that will not bring Ipopt in) */
5210  consdata->maxnonconvexity = SCIPinfinity(scip);
5211  }
5212 
5213  SCIPhashmapFree(&var2index);
5214  SCIPfreeBufferArray(scip, &matrix);
5215 
5216  return SCIP_OKAY;
5217 }
5218 
5219 /** check whether indefinite constraint function is factorable and store corresponding coefficients */
5220 static
5221 SCIP_RETCODE checkFactorable(
5222  SCIP* scip, /**< SCIP data structure */
5223  SCIP_CONS* cons /**< constraint */
5224  )
5225 {
5226  SCIP_BILINTERM* bilinterm;
5227  SCIP_CONSDATA* consdata;
5228  SCIP_Real* a;
5229  SCIP_Real* eigvals;
5230  SCIP_Real sigma1;
5231  SCIP_Real sigma2;
5232  SCIP_Bool success;
5233  int n;
5234  int i;
5235  int idx1;
5236  int idx2;
5237  int posidx;
5238  int negidx;
5239 
5240  assert(scip != NULL);
5241  assert(cons != NULL);
5242 
5243  consdata = SCIPconsGetData(cons);
5244  assert(consdata != NULL);
5245  assert(consdata->factorleft == NULL);
5246  assert(consdata->factorright == NULL);
5247 
5248  /* we don't need this if there are no bilinear terms */
5249  if( consdata->nbilinterms == 0 )
5250  return SCIP_OKAY;
5251 
5252  /* write constraint as lhs <= linear + x'^T A x' <= rhs where x' = (x,1) and
5253  * A = ( Q b/2 )
5254  * ( b^T/2 0 )
5255  * compute an eigenvalue factorization of A and check if there are one positive and one negative eigenvalue
5256  * if so, then let sigma1^2 and -sigma2^2 be these eigenvalues and v1 and v2 be the first two rows of the inverse eigenvector matrix
5257  * thus, x'^T A x' = sigma1^2 (v1^T x')^2 - sigma2^2 (v2^T x')^2
5258  * = (sigma1 (v1^T x') - sigma2 (v2^T x')) * (sigma1 (v1^T x') + sigma2 (v2^T x'))
5259  * we then store sigma1 v1^T - sigma2 v2^T as left factor coef, and sigma1 v1^T + sigma2 v2^T as right factor coef
5260  */
5261 
5262  /* if we already know that there are only positive or only negative eigenvalues, then don't try */
5263  if( consdata->iscurvchecked && (consdata->isconvex || consdata->isconcave) )
5264  return SCIP_OKAY;
5265 
5266  n = consdata->nquadvars + 1;
5267 
5268  /* @todo handle case n=3 explicitly */
5269 
5270  /* skip too large matrices */
5271  if( n > 50 )
5272  return SCIP_OKAY;
5273 
5274  /* need routine to compute eigenvalues/eigenvectors */
5275  if( !SCIPisIpoptAvailableIpopt() )
5276  return SCIP_OKAY;
5277 
5278  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
5279 
5280  SCIP_CALL( SCIPallocBufferArray(scip, &a, n*n) );
5281  BMSclearMemoryArray(a, n*n);
5282 
5283  /* set lower triangular entries of A corresponding to bilinear terms */
5284  for( i = 0; i < consdata->nbilinterms; ++i )
5285  {
5286  bilinterm = &consdata->bilinterms[i];
5287 
5288  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var1, &idx1) );
5289  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var2, &idx2) );
5290  assert(idx1 >= 0);
5291  assert(idx2 >= 0);
5292  assert(idx1 != idx2);
5293 
5294  a[MIN(idx1,idx2) * n + MAX(idx1,idx2)] = bilinterm->coef / 2.0;
5295  }
5296 
5297  /* set lower triangular entries of A corresponding to square and linear terms */
5298  for( i = 0; i < consdata->nquadvars; ++i )
5299  {
5300  a[i*n + i] = consdata->quadvarterms[i].sqrcoef;
5301  a[i*n + n-1] = consdata->quadvarterms[i].lincoef / 2.0;
5302  }
5303 
5304  SCIP_CALL( SCIPallocBufferArray(scip, &eigvals, n) );
5305  if( LapackDsyev(TRUE, n, a, eigvals) != SCIP_OKAY )
5306  {
5307  SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors of augmented quadratic form matrix for constraint <%s>.\n", SCIPconsGetName(cons));
5308  goto CLEANUP;
5309  }
5310 
5311  /* check if there is exactly one positive and one negative eigenvalue */
5312  posidx = -1;
5313  negidx = -1;
5314  for( i = 0; i < n; ++i )
5315  {
5316  if( SCIPisPositive(scip, eigvals[i]) )
5317  {
5318  if( posidx == -1 )
5319  posidx = i;
5320  else
5321  break;
5322  }
5323  else if( SCIPisNegative(scip, eigvals[i]) )
5324  {
5325  if( negidx == -1 )
5326  negidx = i;
5327  else
5328  break;
5329  }
5330  }
5331  if( i < n || posidx == -1 || negidx == -1 )
5332  {
5333  SCIPdebugMsg(scip, "Augmented quadratic form of constraint <%s> is not factorable.\n", SCIPconsGetName(cons));
5334  goto CLEANUP;
5335  }
5336  assert(SCIPisPositive(scip, eigvals[posidx]));
5337  assert(SCIPisNegative(scip, eigvals[negidx]));
5338 
5339  /* compute factorleft and factorright */
5340  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->factorleft, consdata->nquadvars + 1) );
5341  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->factorright, consdata->nquadvars + 1) );
5342 
5343  /* eigenvectors are stored in a, inverse eigenvector matrix is transposed of a
5344  * it seems that v1 and v2 are at &a[posidx*n] and &a[negidx*n]
5345  */
5346  sigma1 = sqrt( eigvals[posidx]);
5347  sigma2 = sqrt(-eigvals[negidx]);
5348  for( i = 0; i < n; ++i )
5349  {
5350  consdata->factorleft[i] = sigma1 * a[posidx * n + i] - sigma2 * a[negidx * n + i];
5351  consdata->factorright[i] = sigma1 * a[posidx * n + i] + sigma2 * a[negidx * n + i];
5352  /* set almost-zero elements to zero */
5353  if( SCIPisZero(scip, consdata->factorleft[i]) )
5354  consdata->factorleft[i] = 0.0;
5355  if( SCIPisZero(scip, consdata->factorright[i]) )
5356  consdata->factorright[i] = 0.0;
5357  }
5358 
5359 #ifdef SCIP_DEBUG
5360  SCIPdebugMsg(scip, "constraint <%s> has factorable quadratic form: (%g", SCIPconsGetName(cons), consdata->factorleft[n-1]);
5361  for( i = 0; i < consdata->nquadvars; ++i )
5362  {
5363  if( consdata->factorleft[i] != 0.0 )
5364  SCIPdebugMsgPrint(scip, " %+g<%s>", consdata->factorleft[i], SCIPvarGetName(consdata->quadvarterms[i].var));
5365  }
5366  SCIPdebugMsgPrint(scip, ") * (%g", consdata->factorright[n-1]);
5367  for( i = 0; i < consdata->nquadvars; ++i )
5368  {
5369  if( consdata->factorright[i] != 0.0 )
5370  SCIPdebugMsgPrint(scip, " %+g<%s>", consdata->factorright[i], SCIPvarGetName(consdata->quadvarterms[i].var));
5371  }
5372  SCIPdebugMsgPrint(scip, ")\n");
5373 #endif
5374 
5375  /* check whether factorleft * factorright^T is matrix of augmented quadratic form
5376  * we check here only the nonzero entries from the quadratic form
5377  */
5378  success = TRUE;
5379 
5380  /* check bilinear terms */
5381  for( i = 0; i < consdata->nbilinterms; ++i )
5382  {
5383  bilinterm = &consdata->bilinterms[i];
5384 
5385  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var1, &idx1) );
5386  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var2, &idx2) );
5387 
5388  if( !SCIPisRelEQ(scip, consdata->factorleft[idx1] * consdata->factorright[idx2] + consdata->factorleft[idx2] * consdata->factorright[idx1], bilinterm->coef) )
5389  {
5390  success = FALSE;
5391  break;
5392  }
5393  }
5394 
5395  /* set lower triangular entries of A corresponding to square and linear terms */
5396  for( i = 0; i < consdata->nquadvars; ++i )
5397  {
5398  if( !SCIPisRelEQ(scip, consdata->factorleft[i] * consdata->factorright[i], consdata->quadvarterms[i].sqrcoef) )
5399  {
5400  success = FALSE;
5401  break;
5402  }
5403 
5404  if( !SCIPisRelEQ(scip, consdata->factorleft[n-1] * consdata->factorright[i] + consdata->factorleft[i] * consdata->factorright[n-1], consdata->quadvarterms[i].lincoef) )
5405  {
5406  success = FALSE;
5407  break;
5408  }
5409  }
5410 
5411  if( !success )
5412  {
5413  SCIPdebugMsg(scip, "Factorization not accurate enough. Dropping it.\n");
5414  SCIPfreeBlockMemoryArray(scip, &consdata->factorleft, consdata->nquadvars + 1);
5415  SCIPfreeBlockMemoryArray(scip, &consdata->factorright, consdata->nquadvars + 1);
5416  }
5417 
5418  CLEANUP:
5419  SCIPfreeBufferArray(scip, &eigvals);
5420  SCIPfreeBufferArray(scip, &a);
5421 
5422  return SCIP_OKAY;
5423 }
5424 
5425 /** computes activity and violation of a constraint
5426  *
5427  * If solution violates bounds by more than feastol, the violation is still computed, but *solviolbounds is set to TRUE
5428  */
5429 static
5430 SCIP_RETCODE computeViolation(
5431  SCIP* scip, /**< SCIP data structure */
5432  SCIP_CONS* cons, /**< constraint */
5433  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
5434  SCIP_Bool* solviolbounds /**< buffer to store whether quadratic variables in solution are outside their bounds by more than feastol */
5435  )
5436 { /*lint --e{666}*/
5437  SCIP_CONSDATA* consdata;
5438  SCIP_Real varval;
5439  SCIP_Real varval2;
5440  SCIP_Real absviol;
5441  SCIP_Real relviol;
5442  SCIP_VAR* var;
5443  SCIP_VAR* var2;
5444  int i;
5445  int j;
5446 
5447  assert(scip != NULL);
5448  assert(cons != NULL);
5449  assert(solviolbounds != NULL);
5450 
5451  consdata = SCIPconsGetData(cons);
5452  assert(consdata != NULL);
5453 
5454  *solviolbounds = FALSE;
5455  consdata->activity = 0.0;
5456  consdata->lhsviol = 0.0;
5457  consdata->rhsviol = 0.0;
5458 
5459  for( i = 0; i < consdata->nlinvars; ++i )
5460  {
5461  SCIP_Real activity;
5462 
5463  var = consdata->linvars[i];
5464  varval = SCIPgetSolVal(scip, sol, var);
5465  activity = consdata->lincoefs[i] * varval;
5466 
5467  /* the contribution of a variable with |varval| = +inf is +inf when activity > 0.0, -inf when activity < 0.0, and
5468  * 0.0 otherwise
5469  */
5470  if( SCIPisInfinity(scip, REALABS(varval)) )
5471  {
5472  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
5473  {
5474  consdata->activity = SCIPinfinity(scip);
5475  consdata->rhsviol = SCIPinfinity(scip);
5476  return SCIP_OKAY;
5477  }
5478 
5479  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
5480  {
5481  consdata->activity = -SCIPinfinity(scip);
5482  consdata->lhsviol = SCIPinfinity(scip);
5483  return SCIP_OKAY;
5484  }
5485  }
5486 
5487  consdata->activity += activity;
5488  }
5489 
5490  for( j = 0; j < consdata->nquadvars; ++j )
5491  {
5492  SCIP_Real activity;
5493 
5494  var = consdata->quadvarterms[j].var;
5495  varval = SCIPgetSolVal(scip, sol, var);
5496  activity = (consdata->quadvarterms[j].lincoef + consdata->quadvarterms[j].sqrcoef * varval) * varval;
5497 
5498  /* the contribution of a variable with |varval| = +inf is +inf when activity > 0.0, -inf when activity < 0.0, and
5499  * 0.0 otherwise
5500  */
5501  if( SCIPisInfinity(scip, REALABS(varval)) )
5502  {
5503  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
5504  {
5505  consdata->activity = SCIPinfinity(scip);
5506  consdata->rhsviol = SCIPinfinity(scip);
5507  return SCIP_OKAY;
5508  }
5509 
5510  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
5511  {
5512  consdata->activity = -SCIPinfinity(scip);
5513  consdata->lhsviol = SCIPinfinity(scip);
5514  return SCIP_OKAY;
5515  }
5516  }
5517 
5518  /* project onto local box, in case the LP solution is slightly outside the bounds (which is not our job to enforce) */
5519  if( sol == NULL )
5520  {
5521  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might violated the variable bounds */
5522  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
5523  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
5524  *solviolbounds = TRUE;
5525  else
5526  {
5527  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
5528  activity = (consdata->quadvarterms[j].lincoef + consdata->quadvarterms[j].sqrcoef * varval) * varval;
5529  }
5530  }
5531 
5532  consdata->activity += activity;
5533  }
5534 
5535  for( j = 0; j < consdata->nbilinterms; ++j )
5536  {
5537  SCIP_Real activity;
5538 
5539  var = consdata->bilinterms[j].var1;
5540  var2 = consdata->bilinterms[j].var2;
5541  varval = SCIPgetSolVal(scip, sol, var);
5542  varval2 = SCIPgetSolVal(scip, sol, var2);
5543 
5544  /* project onto local box, in case the LP solution is slightly outside the bounds (which is not our job to enforce) */
5545  if( sol == NULL )
5546  {
5547  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might violated the variable bounds */
5548  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
5549  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
5550  *solviolbounds = TRUE;
5551  else
5552  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
5553 
5554  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might violated the variable bounds */
5555  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var2)) && !SCIPisFeasGE(scip, varval2, SCIPvarGetLbLocal(var2))) ||
5556  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var2)) && !SCIPisFeasLE(scip, varval2, SCIPvarGetUbLocal(var2))) )
5557  *solviolbounds = TRUE;
5558  else
5559  varval2 = MAX(SCIPvarGetLbLocal(var2), MIN(SCIPvarGetUbLocal(var2), varval2));
5560  }
5561 
5562  activity = consdata->bilinterms[j].coef * varval * varval2;
5563 
5564  /* consider var*var2 as a new variable and handle it as it would appear linearly */
5565  if( SCIPisInfinity(scip, REALABS(varval*varval2)) )
5566  {
5567  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
5568  {
5569  consdata->activity = SCIPinfinity(scip);
5570  consdata->rhsviol = SCIPinfinity(scip);
5571  return SCIP_OKAY;
5572  }
5573 
5574  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
5575  {
5576  consdata->activity = -SCIPinfinity(scip);
5577  consdata->lhsviol = SCIPinfinity(scip);
5578  return SCIP_OKAY;
5579  }
5580  }
5581 
5582  consdata->activity += activity;
5583  }
5584 
5585  absviol = 0.0;
5586  relviol = 0.0;
5587  /* compute absolute violation left hand side */
5588  if( consdata->activity < consdata->lhs && !SCIPisInfinity(scip, -consdata->lhs) )
5589  {
5590  consdata->lhsviol = consdata->lhs - consdata->activity;
5591  absviol = consdata->lhsviol;
5592  relviol = SCIPrelDiff(consdata->lhs, consdata->activity);
5593  }
5594  else
5595  consdata->lhsviol = 0.0;
5596 
5597  /* compute absolute violation right hand side */
5598  if( consdata->activity > consdata->rhs && !SCIPisInfinity(scip, consdata->rhs) )
5599  {
5600  consdata->rhsviol = consdata->activity - consdata->rhs;
5601  absviol = consdata->rhsviol;
5602  relviol = SCIPrelDiff(consdata->activity, consdata->rhs);
5603  }
5604  else
5605  consdata->rhsviol = 0.0;
5606 
5607  /* update absolute and relative violation of the solution */
5608  if( sol != NULL )
5609  SCIPupdateSolConsViolation(scip, sol, absviol, relviol);
5610 
5611  return SCIP_OKAY;
5612 }
5613 
5614 /** computes violation of a set of constraints */
5615 static
5616 SCIP_RETCODE computeViolations(
5617  SCIP* scip, /**< SCIP data structure */
5618  SCIP_CONS** conss, /**< constraints */
5619  int nconss, /**< number of constraints */
5620  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
5621  SCIP_Bool* solviolbounds, /**< buffer to store whether quadratic variables in solution are outside their bounds by more than feastol in some constraint */
5622  SCIP_CONS** maxviolcon /**< buffer to store constraint with largest violation, or NULL if solution is feasible */
5623  )
5624 {
5625  SCIP_CONSDATA* consdata;
5626  SCIP_Real viol;
5627  SCIP_Real maxviol;
5628  SCIP_Bool solviolbounds1;
5629  int c;
5630 
5631  assert(scip != NULL);
5632  assert(conss != NULL || nconss == 0);
5633  assert(solviolbounds != NULL);
5634  assert(maxviolcon != NULL);
5635 
5636  *solviolbounds = FALSE;
5637  *maxviolcon = NULL;
5638 
5639  maxviol = 0.0;
5640 
5641  for( c = 0; c < nconss; ++c )
5642  {
5643  assert(conss != NULL);
5644  assert(conss[c] != NULL);
5645 
5646  SCIP_CALL( computeViolation(scip, conss[c], sol, &solviolbounds1) );
5647  *solviolbounds |= solviolbounds1;
5648 
5649  consdata = SCIPconsGetData(conss[c]);
5650  assert(consdata != NULL);
5651 
5652  viol = MAX(consdata->lhsviol, consdata->rhsviol);
5653  if( viol > maxviol && SCIPisGT(scip, viol, SCIPfeastol(scip)) )
5654  {
5655  maxviol = viol;
5656  *maxviolcon = conss[c];
5657  }
5658  }
5659 
5660  return SCIP_OKAY;
5661 }
5662 
5663 
5664 /** index comparison method for bilinear terms */
5665 static
5666 SCIP_DECL_SORTINDCOMP(bilinTermComp2)
5667 { /*lint --e{715}*/
5668  SCIP_BILINTERM* bilinterms = (SCIP_BILINTERM*)dataptr;
5669  int var1cmp;
5670 
5671  assert(bilinterms != NULL);
5672 
5673  var1cmp = SCIPvarCompare(bilinterms[ind1].var1, bilinterms[ind2].var1);
5674  if( var1cmp != 0 )
5675  return var1cmp;
5676 
5677  return SCIPvarCompare(bilinterms[ind1].var2, bilinterms[ind2].var2);
5678 }
5679 
5680 /** volume comparison method for bilinear terms; prioritizes bilinear products with a larger volume */
5681 static
5682 SCIP_DECL_SORTINDCOMP(bilinTermCompVolume)
5683 { /*lint --e{715}*/
5684  SCIP_BILINTERM* bilinterms = (SCIP_BILINTERM*)dataptr;
5685  SCIP_Real vol1;
5686  SCIP_Real vol2;
5687 
5688  assert(bilinterms != NULL);
5689 
5690  vol1 = (SCIPvarGetUbLocal(bilinterms[ind1].var1) - SCIPvarGetLbLocal(bilinterms[ind1].var1))
5691  * (SCIPvarGetUbLocal(bilinterms[ind1].var2) - SCIPvarGetLbLocal(bilinterms[ind1].var2));
5692  vol2 = (SCIPvarGetUbLocal(bilinterms[ind2].var1) - SCIPvarGetLbLocal(bilinterms[ind2].var1))
5693  * (SCIPvarGetUbLocal(bilinterms[ind2].var2) - SCIPvarGetLbLocal(bilinterms[ind2].var2));
5694 
5695  if( vol1 > vol2 )
5696  return -1;
5697  else if( vol1 < vol2 )
5698  return 1;
5699  return bilinTermComp2(dataptr, ind1, ind2);
5700 }
5701 
5702 /** helper function to sort all bilinear terms in the constraint handler data */
5703 static
5704 SCIP_RETCODE sortAllBilinTerms(
5705  SCIP* scip, /**< SCIP data structure */
5706  SCIP_BILINTERM* bilinterms, /**< array containing all bilinear terms */
5707  int nbilinterms, /**< total number of bilinear terms */
5708  SCIP_CONS** bilinconss, /**< array for mapping each term to its constraint */
5709  int* bilinposs /**< array for mapping each term to its position in the corresponding
5710  * bilinconss constraint */
5711  )
5712 {
5713  int* perm;
5714  int i;
5715  int nexti;
5716  int v;
5717  SCIP_BILINTERM bilinterm;
5718  SCIP_CONS* bilincons;
5719  int bilinpos;
5720 
5721  assert(scip != NULL);
5722  assert(bilinterms != NULL);
5723  assert(nbilinterms > 0);
5724  assert(bilinconss != NULL);
5725  assert(bilinposs != NULL);
5726 
5727  /* get temporary memory to store the sorted permutation and the inverse permutation */
5728  SCIP_CALL( SCIPallocBufferArray(scip, &perm, nbilinterms) );
5729 
5730  /* call quicksort */
5731  SCIPsort(perm, bilinTermCompVolume, (void*)bilinterms, nbilinterms);
5732 
5733  /* permute the bilinear terms according to the resulting permutation */
5734  for( v = 0; v < nbilinterms; ++v )
5735  {
5736  if( perm[v] != v )
5737  {
5738  bilinterm = bilinterms[v];
5739  bilincons = bilinconss[v];
5740  bilinpos = bilinposs[v];
5741 
5742  i = v;
5743  do
5744  {
5745  assert(0 <= perm[i] && perm[i] < nbilinterms);
5746  assert(perm[i] != i);
5747 
5748  bilinterms[i] = bilinterms[perm[i]];
5749  bilinconss[i] = bilinconss[perm[i]];
5750  bilinposs[i] = bilinposs[perm[i]];
5751 
5752  nexti = perm[i];
5753  perm[i] = i;
5754  i = nexti;
5755  }
5756  while( perm[i] != v );
5757  bilinterms[i] = bilinterm;
5758  bilinconss[i] = bilincons;
5759  bilinposs[i] = bilinpos;
5760  perm[i] = i;
5761  }
5762  }
5763 
5764  /* free temporary memory */
5765  SCIPfreeBufferArray(scip, &perm);
5766 
5767  return SCIP_OKAY;
5768 }
5769 
5770 /** stores all bilinear terms in the quadratic constraint handler data; in addition, for each bilinear term we store
5771  * the number of nonconvex constraints that require to over- or underestimate this term, which only depends on the
5772  * lhs, rhs, and the bilinear coefficient
5773  */
5774 static
5775 SCIP_RETCODE storeAllBilinearTerms(
5776  SCIP* scip, /**< SCIP data structure */
5777  SCIP_CONSHDLRDATA* conshdlrdata, /**< constraint handler data */
5778  SCIP_CONS** conss, /**< constraints to process */
5779  int nconss /**< number of constraints */
5780  )
5781 {
5782  SCIP_BILINTERM* bilinterms;
5783  SCIP_CONS** bilincons;
5784  int* bilinpos;
5785  int nbilinterms;
5786  int pos;
5787  int c;
5788  int i;
5789 
5790  assert(scip != NULL);
5791  assert(conshdlrdata != NULL);
5792  assert(conss != NULL);
5793 
5794  /* check for all cases for which we don't want to spend time for collecting all bilinear terms */
5795  if( nconss == 0 || conshdlrdata->storedbilinearterms || SCIPgetSubscipDepth(scip) != 0 || SCIPgetDepth(scip) >= 1
5796  || SCIPinProbing(scip) || SCIPinDive(scip) )
5797  return SCIP_OKAY;
5798 
5799  assert(conshdlrdata->bilinestimators == NULL);
5800  assert(conshdlrdata->nbilinterms == 0);
5801 
5802  conshdlrdata->storedbilinearterms = TRUE;
5803  nbilinterms = 0;
5804 
5805  /* count the number of bilinear terms (including duplicates) */
5806  for( c = 0; c < nconss; ++c )
5807  {
5808  SCIP_CONSDATA* consdata = SCIPconsGetData(conss[c]);
5809  assert(consdata != NULL);
5810  nbilinterms += consdata->nbilinterms;
5811  }
5812 
5813  /* no bilinear terms available -> stop */
5814  if( nbilinterms == 0 )
5815  return SCIP_OKAY;
5816 
5817  /* allocate temporary memory for sorting all bilinear terms (including duplicates) */
5818  SCIP_CALL( SCIPallocBufferArray(scip, &bilinterms, nbilinterms) );
5819  SCIP_CALL( SCIPallocBufferArray(scip, &bilincons, nbilinterms) );
5820  SCIP_CALL( SCIPallocBufferArray(scip, &bilinpos, nbilinterms) );
5821 
5822  /* copy all bilinear terms; note that we need separate entries for x*y and y*x */
5823  pos = 0;
5824  for( c = 0; c < nconss; ++c )
5825  {
5826  SCIP_CONSDATA* consdata = SCIPconsGetData(conss[c]);
5827 
5828  /* allocate memory to store the later computed indices of each bilinear term in the bilinterms array of the
5829  * constraint handler data
5830  */
5831  if( consdata->nbilinterms > 0 )
5832  {
5833  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->bilintermsidx, consdata->nbilinterms) );
5834  }
5835 
5836  for( i = 0; i < consdata->nbilinterms; ++i )
5837  {
5838  assert(consdata->bilinterms != NULL);
5839  assert(consdata->bilinterms[i].var1 != consdata->bilinterms[i].var2);
5840 
5841  /* add xy */
5842  bilinterms[pos] = consdata->bilinterms[i];
5843  bilincons[pos] = conss[c];
5844  bilinpos[pos] = i;
5845  ++pos;
5846 
5847  /* invalidate bilinear term index */
5848  assert(consdata->bilintermsidx != NULL);
5849  consdata->bilintermsidx[i] = -1;
5850  }
5851  }
5852  assert(pos == nbilinterms);
5853 
5854  /* sorts all bilinear terms (including duplicates) */
5855  SCIP_CALL( sortAllBilinTerms(scip, bilinterms, nbilinterms, bilincons, bilinpos) );
5856 
5857  /* count the number of bilinear terms without duplicates */
5858  conshdlrdata->nbilinterms = nbilinterms;
5859  for( i = 0; i