Scippy

SCIP

Solving Constraint Integer Programs

presol_domcol.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-2014 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 email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file presol_domcol.c
17  * @ingroup PRESOLVERS
18  * @brief dominated column presolver
19  * @author Dieter Weninger
20  * @author Gerald Gamrath
21  *
22  * This presolver looks for dominance relations between variable pairs.
23  * From a dominance relation and certain bound/clique-constellations
24  * variable fixings mostly at the lower bound of the dominated variable can be derived.
25  * Additionally it is possible to improve bounds by predictive bound strengthening.
26  *
27  * @todo Also run on general CIPs, if the number of locks of the investigated variables comes only from (upgraded)
28  * linear constraints.
29  *
30  * @todo Instead of choosing variables for comparison out of one row, we should try to use 'hashvalues' for columns that
31  * indicate in which constraint type (<=, >=, or ranged row / ==) they are existing. Then sort the variables (and
32  * corresponding data) after the ranged row/equation hashvalue and only try to derive dominance on variables with
33  * the same hashvalue on ranged row/equation and fitting hashvalues for the other constraint types.
34  *
35  */
36 
37 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
38 
39 #include <stdio.h>
40 #include <assert.h>
41 #include <string.h>
42 
43 /* includes necessary for matrix data structure */
44 #include "scip/cons_knapsack.h"
45 #include "scip/cons_linear.h"
46 #include "scip/cons_logicor.h"
47 #include "scip/cons_setppc.h"
48 #include "scip/cons_varbound.h"
49 
50 #include "presol_domcol.h"
51 
52 #define PRESOL_NAME "domcol"
53 #define PRESOL_DESC "dominated column presolver"
54 #define PRESOL_PRIORITY 20000000 /**< priority of the presolver (>= 0: before, < 0: after constraint handlers) */
55 #define PRESOL_MAXROUNDS -1 /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */
56 #define PRESOL_DELAY TRUE /**< should presolver be delayed, if other presolvers found reductions? */
57 
58 #define DEFAULT_NUMMINPAIRS 1024 /**< minimal number of pair comparisons */
59 #define DEFAULT_NUMMAXPAIRS 1048576 /**< maximal number of pair comparisons */
60 
61 #define DEFAULT_PREDBNDSTR FALSE /**< should predictive bound strengthening be applied? */
62 #define DEFAULT_SINGCOLSTUFF TRUE /**< should singleton columns stuffing be applied? */
63 
64 /*
65  * Data structures
66  */
67 
68 /** control parameters */
69 struct SCIP_PresolData
70 {
71  int numminpairs; /**< minimal number of pair comparisons */
72  int nummaxpairs; /**< maximal number of pair comparisons */
73  int numcurrentpairs; /**< current number of pair comparisons */
74 
75  SCIP_Bool predbndstr; /**< flag indicating if predictive bound strengthening should be applied */
76  SCIP_Bool singcolstuffing; /**< flag indicating if singleton columns stuffing should be applied */
77 #ifdef SCIP_STATISTIC
78  int nfixingsstuffing; /**< number of fixings done by singleton columns stuffing */
79  int nfixingsdomcol; /**< number of fixings only from dominated columns */
80 #endif
81 };
82 
83 /** type of fixing direction */
84 enum Fixingdirection
85 {
86  FIXATLB = -1, /**< fix variable at lower bound */
87  NOFIX = 0, /**< do not fix variable */
88  FIXATUB = 1 /**< fix variable at upper bound */
89 };
90 typedef enum Fixingdirection FIXINGDIRECTION;
91 
92 
93 /********************************************************************/
94 /************** matrix data structure and functions *****************/
95 
96 /** constraint matrix data structure in column and row major format */
97 struct ConstraintMatrix
98 {
99  SCIP_Real* colmatval; /**< coefficients in column major format */
100  int* colmatind; /**< row indexes in column major format */
101  int* colmatbeg; /**< column storage offset */
102  int* colmatcnt; /**< number of row entries per column */
103  int ncols; /**< complete number of columns */
104  SCIP_Real* lb; /**< lower bound per variable */
105  SCIP_Real* ub; /**< upper bound per variable */
106  int* nuplocks; /**< number of up locks per variable */
107  int* ndownlocks; /**< number of down locks per variable */
108 
109  SCIP_VAR** vars; /**< variables pointer */
110 
111  SCIP_Real* rowmatval; /**< coefficients in row major format */
112  int* rowmatind; /**< column indexed in row major format */
113  int* rowmatbeg; /**< row storage offset */
114  int* rowmatcnt; /**< number of column entries per row */
115 #ifdef SCIP_DEBUG
116  const char** rowname; /**< name of row */
117 #endif
118  int nrows; /**< complete number of rows */
119  SCIP_Real* lhs; /**< left hand side per row */
120  SCIP_Real* rhs; /**< right hand side per row */
121  SCIP_Bool* isrhsinfinite; /**< is right hand side infinity */
122  int nnonzs; /**< sparsity counter */
123  SCIP_Real* minactivity; /**< min activity per row */
124  SCIP_Real* maxactivity; /**< max activity per row */
125  int* minactivityneginf; /**< min activity negative infinity counter */
126  int* minactivityposinf; /**< min activity positive infinity counter */
127  int* maxactivityneginf; /**< max activity negative infinity counter */
128  int* maxactivityposinf; /**< max activity positive infinity counter */
129 };
130 typedef struct ConstraintMatrix CONSTRAINTMATRIX;
131 
132 /** transforms given variables, scalars, and constant to the corresponding active variables, scalars, and constant */
133 static
134 SCIP_RETCODE getActiveVariables(
135  SCIP* scip, /**< SCIP data structure */
136  SCIP_VAR*** vars, /**< vars array to get active variables for */
137  SCIP_Real** scalars, /**< scalars a_1, ..., a_n in linear sum a_1*x_1 + ... + a_n*x_n + c */
138  int* nvars, /**< pointer to number of variables and values in vars and vals array */
139  SCIP_Real* constant /**< pointer to constant c in linear sum a_1*x_1 + ... + a_n*x_n + c */
140  )
141 {
142  int requiredsize;
143 
144  assert(scip != NULL);
145  assert(vars != NULL);
146  assert(scalars != NULL);
147  assert(*vars != NULL);
148  assert(*scalars != NULL);
149  assert(nvars != NULL);
150  assert(constant != NULL);
151 
152  SCIP_CALL( SCIPgetProbvarLinearSum(scip, *vars, *scalars, nvars, *nvars, constant, &requiredsize, TRUE) );
153 
154  if( requiredsize > *nvars )
155  {
156  SCIP_CALL( SCIPreallocBufferArray(scip, vars, requiredsize) );
157  SCIP_CALL( SCIPreallocBufferArray(scip, scalars, requiredsize) );
158 
159  /* call function a second time with enough memory */
160  SCIP_CALL( SCIPgetProbvarLinearSum(scip, *vars, *scalars, nvars, requiredsize, constant, &requiredsize, TRUE) );
161  assert(requiredsize <= *nvars);
162  }
163 
164  return SCIP_OKAY;
165 }
166 
167 /** add one row to the constraint matrix */
168 static
169 SCIP_RETCODE addRow(
170  SCIP* scip, /**< SCIP data structure */
171  CONSTRAINTMATRIX* matrix, /**< constraint matrix */
172 #ifdef SCIP_DEBUG
173  const char* name, /**< name of constraint corresponding to row */
174 #endif
175  SCIP_VAR** vars, /**< variables of this row */
176  SCIP_Real* vals, /**< coefficients of this row */
177  int nvars, /**< number of variables of this row */
178  SCIP_Real lhs, /**< left hand side */
179  SCIP_Real rhs /**< right hand side */
180  )
181 {
182  int j;
183  int probindex;
184  int rowidx;
185  SCIP_Real factor;
186  SCIP_Bool rangedorequality;
187 
188  assert(vars != NULL);
189  assert(vals != NULL);
190 
191  rowidx = matrix->nrows;
192  rangedorequality = FALSE;
193 
194  if( SCIPisInfinity(scip, -lhs) )
195  {
196  factor = -1.0;
197  matrix->lhs[rowidx] = -rhs;
198  matrix->rhs[rowidx] = SCIPinfinity(scip);
199  matrix->isrhsinfinite[rowidx] = TRUE;
200  }
201  else
202  {
203  factor = 1.0;
204  matrix->lhs[rowidx] = lhs;
205  matrix->rhs[rowidx] = rhs;
206  matrix->isrhsinfinite[rowidx] = SCIPisInfinity(scip, matrix->rhs[rowidx]);
207 
208  if( !SCIPisInfinity(scip, rhs) )
209  rangedorequality = TRUE;
210  }
211  assert(!SCIPisInfinity(scip, -matrix->lhs[rowidx]));
212 
213 
214 #ifdef SCIP_DEBUG
215  matrix->rowname[rowidx] = name;
216 #endif
217 
218  matrix->rowmatbeg[rowidx] = matrix->nnonzs;
219 
220  /* = or ranged */
221  if( rangedorequality )
222  {
223  assert(factor > 0);
224 
225  for( j = 0; j < nvars; j++ )
226  {
227  matrix->rowmatval[matrix->nnonzs] = factor * vals[j];
228  probindex = SCIPvarGetProbindex(vars[j]);
229  assert(matrix->vars[probindex] == vars[j]);
230 
231  matrix->nuplocks[probindex]++;
232  matrix->ndownlocks[probindex]++;
233 
234  assert(0 <= probindex && probindex < matrix->ncols);
235  matrix->rowmatind[matrix->nnonzs] = probindex;
236  (matrix->nnonzs)++;
237  }
238  }
239  /* >= or <= */
240  else
241  {
242  for( j = 0; j < nvars; j++ )
243  {
244  /* due to the factor, <= constraints will be transfered to >= */
245  matrix->rowmatval[matrix->nnonzs] = factor * vals[j];
246  probindex = SCIPvarGetProbindex(vars[j]);
247  assert(matrix->vars[probindex] == vars[j]);
248 
249  if( matrix->rowmatval[matrix->nnonzs] > 0 )
250  matrix->ndownlocks[probindex]++;
251  else
252  {
253  assert(matrix->rowmatval[matrix->nnonzs] < 0);
254  matrix->nuplocks[probindex]++;
255  }
256 
257  assert(0 <= probindex && probindex < matrix->ncols);
258  matrix->rowmatind[matrix->nnonzs] = probindex;
259  (matrix->nnonzs)++;
260  }
261  }
262 
263  matrix->rowmatcnt[rowidx] = matrix->nnonzs - matrix->rowmatbeg[rowidx];
264 
265  ++(matrix->nrows);
266 
267  return SCIP_OKAY;
268 }
269 
270 /** add one constraint to matrix */
271 static
272 SCIP_RETCODE addConstraint(
273  SCIP* scip, /**< current scip instance */
274  CONSTRAINTMATRIX* matrix, /**< constraint matrix */
275 #ifdef SCIP_DEBUG
276  const char* name, /**< name of constraint corresponding to row */
277 #endif
278  SCIP_VAR** vars, /**< variables of this constraint */
279  SCIP_Real* vals, /**< variable coefficients of this constraint */
280  int nvars, /**< number of variables */
281  SCIP_Real lhs, /**< left hand side */
282  SCIP_Real rhs /**< right hand side */
283  )
284 {
285  SCIP_VAR** activevars;
286  SCIP_Real* activevals;
287  SCIP_Real activeconstant;
288  int nactivevars;
289  int v;
290 
291  assert(scip != NULL);
292  assert(matrix != NULL);
293  assert(vars != NULL || nvars == 0);
294  assert(SCIPisLE(scip, lhs, rhs));
295 
296  /* constraint is redundant */
297  if( SCIPisInfinity(scip, -lhs) && SCIPisInfinity(scip, rhs) )
298  return SCIP_OKAY;
299 
300  /* we do not add empty constraints to the matrix */
301  if( nvars == 0 )
302  return SCIP_OKAY;
303 
304  activevars = NULL;
305  activevals = NULL;
306  nactivevars = nvars;
307  activeconstant = 0.0;
308 
309  /* duplicate variable and value array */
310  SCIP_CALL( SCIPduplicateBufferArray(scip, &activevars, vars, nactivevars ) );
311  if( vals != NULL )
312  {
313  SCIP_CALL( SCIPduplicateBufferArray(scip, &activevals, vals, nactivevars ) );
314  }
315  else
316  {
317  SCIP_CALL( SCIPallocBufferArray(scip, &activevals, nactivevars) );
318 
319  for( v = 0; v < nactivevars; v++ )
320  activevals[v] = 1.0;
321  }
322 
323  /* retransform given variables to active variables */
324  SCIP_CALL( getActiveVariables(scip, &activevars, &activevals, &nactivevars, &activeconstant) );
325 
326  /* adapt left and right hand side */
327  if( !SCIPisInfinity(scip, -lhs) )
328  lhs -= activeconstant;
329  if( !SCIPisInfinity(scip, rhs) )
330  rhs -= activeconstant;
331 
332  /* add single row to matrix */
333  if( nactivevars > 0 )
334  {
335 #ifdef SCIP_DEBUG
336  SCIP_CALL( addRow(scip, matrix, name, activevars, activevals, nactivevars, lhs, rhs) );
337 #else
338  SCIP_CALL( addRow(scip, matrix, activevars, activevals, nactivevars, lhs, rhs) );
339 #endif
340  }
341 
342  /* free buffer arrays */
343  SCIPfreeBufferArray(scip, &activevals);
344  SCIPfreeBufferArray(scip, &activevars);
345 
346  return SCIP_OKAY;
347 }
348 
349 /** transform row major format into column major format */
350 static
351 SCIP_RETCODE setColumnMajorFormat(
352  SCIP* scip, /**< current scip instance */
353  CONSTRAINTMATRIX* matrix /**< constraint matrix */
354  )
355 {
356  int colidx;
357  int i;
358  int* rowpnt;
359  int* rowend;
360  SCIP_Real* valpnt;
361  int* fillidx;
362 
363  assert(scip != NULL);
364  assert(matrix != NULL);
365  assert(matrix->colmatval != NULL);
366  assert(matrix->colmatind != NULL);
367  assert(matrix->colmatbeg != NULL);
368  assert(matrix->colmatcnt != NULL);
369  assert(matrix->rowmatval != NULL);
370  assert(matrix->rowmatind != NULL);
371  assert(matrix->rowmatbeg != NULL);
372  assert(matrix->rowmatcnt != NULL);
373 
374  SCIP_CALL( SCIPallocBufferArray(scip, &fillidx, matrix->ncols) );
375  BMSclearMemoryArray(fillidx, matrix->ncols);
376  BMSclearMemoryArray(matrix->colmatcnt, matrix->ncols);
377 
378  for( i = 0; i < matrix->nrows; i++ )
379  {
380  rowpnt = matrix->rowmatind + matrix->rowmatbeg[i];
381  rowend = rowpnt + matrix->rowmatcnt[i];
382  for( ; rowpnt < rowend; rowpnt++ )
383  {
384  colidx = *rowpnt;
385  (matrix->colmatcnt[colidx])++;
386  }
387  }
388 
389  matrix->colmatbeg[0] = 0;
390  for( i = 0; i < matrix->ncols-1; i++ )
391  {
392  matrix->colmatbeg[i+1] = matrix->colmatbeg[i] + matrix->colmatcnt[i];
393  }
394 
395  for( i = 0; i < matrix->nrows; i++ )
396  {
397  rowpnt = matrix->rowmatind + matrix->rowmatbeg[i];
398  rowend = rowpnt + matrix->rowmatcnt[i];
399  valpnt = matrix->rowmatval + matrix->rowmatbeg[i];
400 
401  for( ; rowpnt < rowend; rowpnt++, valpnt++ )
402  {
403  assert(*rowpnt < matrix->ncols);
404  colidx = *rowpnt;
405  matrix->colmatval[matrix->colmatbeg[colidx] + fillidx[colidx]] = *valpnt;
406  matrix->colmatind[matrix->colmatbeg[colidx] + fillidx[colidx]] = i;
407  fillidx[colidx]++;
408  }
409  }
410 
411  SCIPfreeBufferArray(scip, &fillidx);
412 
413  return SCIP_OKAY;
414 }
415 
416 /** calculate min/max activity per row */
417 static
418 SCIP_RETCODE calcActivityBounds(
419  SCIP* scip, /**< current scip instance */
420  CONSTRAINTMATRIX* matrix /**< constraint matrix */
421  )
422 {
423  SCIP_Real val;
424  int* rowpnt;
425  int* rowend;
426  SCIP_Real* valpnt;
427  int col;
428  int row;
429 
430  assert(scip != NULL);
431  assert(matrix != NULL);
432 
433  for( row = 0; row < matrix->nrows; row++ )
434  {
435  matrix->minactivity[row] = 0;
436  matrix->maxactivity[row] = 0;
437  matrix->minactivityneginf[row] = 0;
438  matrix->minactivityposinf[row] = 0;
439  matrix->maxactivityneginf[row] = 0;
440  matrix->maxactivityposinf[row] = 0;
441 
442  rowpnt = matrix->rowmatind + matrix->rowmatbeg[row];
443  rowend = rowpnt + matrix->rowmatcnt[row];
444  valpnt = matrix->rowmatval + matrix->rowmatbeg[row];
445 
446  for( ; rowpnt < rowend; rowpnt++, valpnt++ )
447  {
448  /* get column index */
449  col = *rowpnt;
450 
451  /* get variable coefficient */
452  val = *valpnt;
453  assert(!SCIPisZero(scip, val));
454 
455  assert(matrix->ncols > col);
456 
457  assert(!SCIPisInfinity(scip, matrix->lb[col]));
458  assert(!SCIPisInfinity(scip, -matrix->ub[col]));
459 
460  /* positive coefficient */
461  if( val > 0.0 )
462  {
463  if( SCIPisInfinity(scip, matrix->ub[col]) )
464  matrix->maxactivityposinf[row]++;
465  else
466  matrix->maxactivity[row] += val * matrix->ub[col];
467 
468  if( SCIPisInfinity(scip, -matrix->lb[col]) )
469  matrix->minactivityneginf[row]++;
470  else
471  matrix->minactivity[row] += val * matrix->lb[col];
472  }
473  /* negative coefficient */
474  else
475  {
476  if( SCIPisInfinity(scip, -matrix->lb[col]) )
477  matrix->maxactivityneginf[row]++;
478  else
479  matrix->maxactivity[row] += val * matrix->lb[col];
480 
481  if( SCIPisInfinity(scip, matrix->ub[col]) )
482  matrix->minactivityposinf[row]++;
483  else
484  matrix->minactivity[row] += val * matrix->ub[col];
485  }
486  }
487  }
488 
489  return SCIP_OKAY;
490 }
491 
492 /** initialize matrix */
493 static
494 SCIP_RETCODE initMatrix(
495  SCIP* scip, /**< current scip instance */
496  CONSTRAINTMATRIX** matrixptr, /**< pointer to constraint matrix object to be initialized */
497  SCIP_Bool* initialized /**< was the initialization successful? */
498  )
499 {
500  CONSTRAINTMATRIX* matrix;
501  SCIP_CONSHDLR** conshdlrs;
502  const char* conshdlrname;
503  SCIP_Bool stopped;
504  SCIP_VAR** vars;
505  SCIP_VAR* var;
506  SCIP_CONS* cons;
507  int nconshdlrs;
508  int nconss;
509  int nnonzstmp;
510  int nvars;
511  int c;
512  int i;
513  int v;
514 
515  nnonzstmp = 0;
516 
517  assert(scip != NULL);
518  assert(matrixptr != NULL);
519  assert(initialized != NULL);
520 
521  *initialized = FALSE;
522 
523  /* return if no variables or constraints are present */
524  if( SCIPgetNVars(scip) == 0 || SCIPgetNConss(scip) == 0 )
525  return SCIP_OKAY;
526 
527  /* loop over all constraint handlers and collect the number of checked constraints */
528  nconshdlrs = SCIPgetNConshdlrs(scip);
529  conshdlrs = SCIPgetConshdlrs(scip);
530  nconss = 0;
531 
532  for( i = 0; i < nconshdlrs; ++i )
533  {
534  int nconshdlrconss;
535 
536  nconshdlrconss = SCIPconshdlrGetNCheckConss(conshdlrs[i]);
537 
538  if( nconshdlrconss > 0 )
539  {
540  conshdlrname = SCIPconshdlrGetName(conshdlrs[i]);
541 
542  if( (strcmp(conshdlrname, "linear") != 0) && (strcmp(conshdlrname, "setppc") != 0)
543  && (strcmp(conshdlrname, "logicor") != 0) && (strcmp(conshdlrname, "knapsack") != 0)
544  && (strcmp(conshdlrname, "varbound") != 0) )
545  {
546  SCIPdebugMessage("unsupported constraint type <%s>: aborting domcol presolver\n", conshdlrname);
547  break;
548  }
549  }
550 
551  nconss += nconshdlrconss;
552  }
553 
554  /* do nothing if we have unsupported constraint types or no checked constraints */
555  if( i < nconshdlrs || nconss == 0 )
556  return SCIP_OKAY;
557 
558  stopped = FALSE;
559 
560  vars = SCIPgetVars(scip);
561  nvars = SCIPgetNVars(scip);
562 
563  /* approximate number of nonzeros by taking for each variable the number of up- and downlocks;
564  * this counts nonzeros in equalities twice, but can be at most two times as high as the exact number
565  */
566  for( i = nvars - 1; i >= 0; --i )
567  {
568  nnonzstmp += SCIPvarGetNLocksDown(vars[i]);
569  nnonzstmp += SCIPvarGetNLocksUp(vars[i]);
570  }
571 
572  /* do nothing if we have no entries */
573  if( nnonzstmp == 0 )
574  return SCIP_OKAY;
575 
576  /* build the matrix structure */
577  SCIP_CALL( SCIPallocBuffer(scip, matrixptr) );
578  matrix = *matrixptr;
579 
580  /* copy vars array and set number of variables */
581  SCIP_CALL( SCIPduplicateBufferArray(scip, &matrix->vars, SCIPgetVars(scip), nvars) );
582  matrix->ncols = nvars;
583 
584  matrix->nrows = 0;
585  matrix->nnonzs = 0;
586 
587  /* allocate memory for columns */
588  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatval, nnonzstmp) );
589  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatind, nnonzstmp) );
590  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatbeg, matrix->ncols) );
591  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatcnt, matrix->ncols) );
592  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->lb, matrix->ncols) );
593  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->ub, matrix->ncols) );
594  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->nuplocks, matrix->ncols) );
595  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->ndownlocks, matrix->ncols) );
596 
597  BMSclearMemoryArray(matrix->nuplocks, matrix->ncols);
598  BMSclearMemoryArray(matrix->ndownlocks, matrix->ncols);
599 
600  for( v = 0; v < matrix->ncols; v++ )
601  {
602  var = matrix->vars[v];
603  assert(var != NULL);
604 
605  matrix->lb[v] = SCIPvarGetLbGlobal(var);
606  matrix->ub[v] = SCIPvarGetUbGlobal(var);
607  }
608 
609  /* allocate memory for rows */
610  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatval, nnonzstmp) );
611  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatind, nnonzstmp) );
612  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatbeg, nconss) );
613  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatcnt, nconss) );
614  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->lhs, nconss) );
615  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rhs, nconss) );
616 
617  /* allocate memory for status of finiteness of row sides */
618  SCIP_CALL( SCIPallocClearMemoryArray(scip, &matrix->isrhsinfinite, nconss) );
619 
620 #ifdef SCIP_DEBUG
621  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowname, nconss) );
622 #endif
623 
624  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->minactivity, nconss) );
625  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->maxactivity, nconss) );
626  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->minactivityneginf, nconss) );
627  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->minactivityposinf, nconss) );
628  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->maxactivityneginf, nconss) );
629  SCIP_CALL( SCIPallocBufferArray(scip, &matrix->maxactivityposinf, nconss) );
630 
631  /* loop a second time over constraints handlers and add constraints to the matrix */
632  for( i = 0; i < nconshdlrs; ++i )
633  {
634  SCIP_CONS** conshdlrconss;
635  int nconshdlrconss;
636 
637  if( SCIPisStopped(scip) )
638  {
639  stopped = TRUE;
640  break;
641  }
642 
643  conshdlrname = SCIPconshdlrGetName(conshdlrs[i]);
644  conshdlrconss = SCIPconshdlrGetCheckConss(conshdlrs[i]);
645  nconshdlrconss = SCIPconshdlrGetNCheckConss(conshdlrs[i]);
646 
647  if( strcmp(conshdlrname, "linear") == 0 )
648  {
649  for( c = 0; c < nconshdlrconss && (c % 1000 != 0 || !SCIPisStopped(scip)); ++c )
650  {
651  cons = conshdlrconss[c];
652  assert(SCIPconsIsTransformed(cons));
653 
654 #ifdef SCIP_DEBUG
655  SCIP_CALL( addConstraint(scip, matrix, SCIPconsGetName(cons), SCIPgetVarsLinear(scip, cons),
656  SCIPgetValsLinear(scip, cons), SCIPgetNVarsLinear(scip, cons),
657  SCIPgetLhsLinear(scip, cons), SCIPgetRhsLinear(scip, cons)) );
658 #else
659  SCIP_CALL( addConstraint(scip, matrix, SCIPgetVarsLinear(scip, cons),
660  SCIPgetValsLinear(scip, cons), SCIPgetNVarsLinear(scip, cons),
661  SCIPgetLhsLinear(scip, cons), SCIPgetRhsLinear(scip, cons)) );
662 #endif
663  }
664  }
665  else if( strcmp(conshdlrname, "setppc") == 0 )
666  {
667  for( c = 0; c < nconshdlrconss && (c % 1000 != 0 || !SCIPisStopped(scip)); ++c )
668  {
669  SCIP_Real lhs;
670  SCIP_Real rhs;
671 
672  cons = conshdlrconss[c];
673  assert(SCIPconsIsTransformed(cons));
674 
675  switch( SCIPgetTypeSetppc(scip, cons) )
676  {
677  case SCIP_SETPPCTYPE_PARTITIONING :
678  lhs = 1.0;
679  rhs = 1.0;
680  break;
681  case SCIP_SETPPCTYPE_PACKING :
682  lhs = -SCIPinfinity(scip);
683  rhs = 1.0;
684  break;
685  case SCIP_SETPPCTYPE_COVERING :
686  lhs = 1.0;
687  rhs = SCIPinfinity(scip);
688  break;
689  default:
690  return SCIP_ERROR;
691  }
692 
693 #ifdef SCIP_DEBUG
694  SCIP_CALL( addConstraint(scip, matrix, SCIPconsGetName(cons), SCIPgetVarsSetppc(scip, cons), NULL,
695  SCIPgetNVarsSetppc(scip, cons), lhs, rhs) );
696 #else
697  SCIP_CALL( addConstraint(scip, matrix, SCIPgetVarsSetppc(scip, cons), NULL,
698  SCIPgetNVarsSetppc(scip, cons), lhs, rhs) );
699 #endif
700  }
701  }
702  else if( strcmp(conshdlrname, "logicor") == 0 )
703  {
704  for( c = 0; c < nconshdlrconss && (c % 1000 != 0 || !SCIPisStopped(scip)); ++c )
705  {
706  cons = conshdlrconss[c];
707  assert(SCIPconsIsTransformed(cons));
708 
709 #ifdef SCIP_DEBUG
710  SCIP_CALL( addConstraint(scip, matrix, SCIPconsGetName(cons), SCIPgetVarsLogicor(scip, cons),
711  NULL, SCIPgetNVarsLogicor(scip, cons), 1.0, SCIPinfinity(scip)) );
712 #else
713  SCIP_CALL( addConstraint(scip, matrix, SCIPgetVarsLogicor(scip, cons),
714  NULL, SCIPgetNVarsLogicor(scip, cons), 1.0, SCIPinfinity(scip)) );
715 #endif
716  }
717  }
718  else if( strcmp(conshdlrname, "knapsack") == 0 )
719  {
720  if( nconshdlrconss > 0 )
721  {
722  SCIP_Real* consvals;
723  int valssize;
724 
725  valssize = 100;
726  SCIP_CALL( SCIPallocBufferArray(scip, &consvals, valssize) );
727 
728  for( c = 0; c < nconshdlrconss && (c % 1000 != 0 || !SCIPisStopped(scip)); ++c )
729  {
730  SCIP_Longint* weights;
731 
732  cons = conshdlrconss[c];
733  assert(SCIPconsIsTransformed(cons));
734 
735  weights = SCIPgetWeightsKnapsack(scip, cons);
736  nvars = SCIPgetNVarsKnapsack(scip, cons);
737 
738  if( nvars > valssize )
739  {
740  valssize = (int) (1.5 * nvars);
741  SCIP_CALL( SCIPreallocBufferArray(scip, &consvals, valssize) );
742  }
743 
744  for( v = 0; v < nvars; v++ )
745  consvals[v] = (SCIP_Real)weights[v];
746 
747 #ifdef SCIP_DEBUG
748  SCIP_CALL( addConstraint(scip, matrix, SCIPconsGetName(cons), SCIPgetVarsKnapsack(scip, cons), consvals,
749  SCIPgetNVarsKnapsack(scip, cons), -SCIPinfinity(scip),
750  (SCIP_Real)SCIPgetCapacityKnapsack(scip, cons)) );
751 #else
752  SCIP_CALL( addConstraint(scip, matrix, SCIPgetVarsKnapsack(scip, cons), consvals,
753  SCIPgetNVarsKnapsack(scip, cons), -SCIPinfinity(scip),
754  (SCIP_Real)SCIPgetCapacityKnapsack(scip, cons)) );
755 #endif
756  }
757 
758  SCIPfreeBufferArray(scip, &consvals);
759  }
760  }
761  else if( strcmp(conshdlrname, "varbound") == 0 )
762  {
763  if( nconshdlrconss > 0 )
764  {
765  SCIP_VAR** consvars;
766  SCIP_Real* consvals;
767 
768  SCIP_CALL( SCIPallocBufferArray(scip, &consvars, 2) );
769  SCIP_CALL( SCIPallocBufferArray(scip, &consvals, 2) );
770  consvals[0] = 1.0;
771 
772  for( c = 0; c < nconshdlrconss && (c % 1000 != 0 || !SCIPisStopped(scip)); ++c )
773  {
774  cons = conshdlrconss[c];
775  assert(SCIPconsIsTransformed(cons));
776 
777  consvars[0] = SCIPgetVarVarbound(scip, cons);
778  consvars[1] = SCIPgetVbdvarVarbound(scip, cons);
779 
780  consvals[1] = SCIPgetVbdcoefVarbound(scip, cons);
781 
782 #ifdef SCIP_DEBUG
783  SCIP_CALL( addConstraint(scip, matrix, SCIPconsGetName(cons), consvars, consvals, 2, SCIPgetLhsVarbound(scip, cons),
784  SCIPgetRhsVarbound(scip, cons)) );
785 #else
786  SCIP_CALL( addConstraint(scip, matrix, consvars, consvals, 2, SCIPgetLhsVarbound(scip, cons),
787  SCIPgetRhsVarbound(scip, cons)) );
788 #endif
789  }
790 
791  SCIPfreeBufferArray(scip, &consvals);
792  SCIPfreeBufferArray(scip, &consvars);
793  }
794  }
795 #ifndef NDEBUG
796  else
797  {
798  assert(nconshdlrconss == 0);
799  }
800 #endif
801  }
802  assert(matrix->nrows <= nconss);
803  assert(matrix->nnonzs <= nnonzstmp);
804 
805  if( !stopped )
806  {
807  /* calculate row activity bounds */
808  SCIP_CALL( calcActivityBounds(scip, matrix) );
809 
810  /* transform row major format into column major format */
811  SCIP_CALL( setColumnMajorFormat(scip, matrix) );
812 
813  *initialized = TRUE;
814  }
815 
816  return SCIP_OKAY;
817 }
818 
819 
820 /** frees the constraint matrix */
821 static
822 void freeMatrix(
823  SCIP* scip, /**< current SCIP instance */
824  CONSTRAINTMATRIX** matrix /**< constraint matrix object */
825  )
826 {
827  assert(scip != NULL);
828  assert(matrix != NULL);
829 
830  if( (*matrix) != NULL )
831  {
832  assert((*matrix)->colmatval != NULL);
833  assert((*matrix)->colmatind != NULL);
834  assert((*matrix)->colmatbeg != NULL);
835  assert((*matrix)->colmatcnt != NULL);
836  assert((*matrix)->lb != NULL);
837  assert((*matrix)->ub != NULL);
838  assert((*matrix)->nuplocks != NULL);
839  assert((*matrix)->ndownlocks != NULL);
840 
841  assert((*matrix)->rowmatval != NULL);
842  assert((*matrix)->rowmatind != NULL);
843  assert((*matrix)->rowmatbeg != NULL);
844  assert((*matrix)->rowmatcnt != NULL);
845  assert((*matrix)->lhs != NULL);
846  assert((*matrix)->rhs != NULL);
847 #ifdef SCIP_DEBUG
848  assert((*matrix)->rowname != NULL);
849 #endif
850 
851  SCIPfreeBufferArray(scip, &((*matrix)->maxactivityposinf));
852  SCIPfreeBufferArray(scip, &((*matrix)->maxactivityneginf));
853  SCIPfreeBufferArray(scip, &((*matrix)->minactivityposinf));
854  SCIPfreeBufferArray(scip, &((*matrix)->minactivityneginf));
855  SCIPfreeBufferArray(scip, &((*matrix)->maxactivity));
856  SCIPfreeBufferArray(scip, &((*matrix)->minactivity));
857 
858 #ifdef SCIP_DEBUG
859  SCIPfreeBufferArray(scip, &((*matrix)->rowname));
860 #endif
861 
862  SCIPfreeMemoryArray(scip, &((*matrix)->isrhsinfinite));
863 
864  SCIPfreeBufferArray(scip, &((*matrix)->rhs));
865  SCIPfreeBufferArray(scip, &((*matrix)->lhs));
866  SCIPfreeBufferArray(scip, &((*matrix)->rowmatcnt));
867  SCIPfreeBufferArray(scip, &((*matrix)->rowmatbeg));
868  SCIPfreeBufferArray(scip, &((*matrix)->rowmatind));
869  SCIPfreeBufferArray(scip, &((*matrix)->rowmatval));
870 
871  SCIPfreeBufferArray(scip, &((*matrix)->ndownlocks));
872  SCIPfreeBufferArray(scip, &((*matrix)->nuplocks));
873  SCIPfreeBufferArray(scip, &((*matrix)->ub));
874  SCIPfreeBufferArray(scip, &((*matrix)->lb));
875  SCIPfreeBufferArray(scip, &((*matrix)->colmatcnt));
876  SCIPfreeBufferArray(scip, &((*matrix)->colmatbeg));
877  SCIPfreeBufferArray(scip, &((*matrix)->colmatind));
878  SCIPfreeBufferArray(scip, &((*matrix)->colmatval));
879 
880  (*matrix)->nrows = 0;
881  (*matrix)->ncols = 0;
882  (*matrix)->nnonzs = 0;
883 
884  SCIPfreeBufferArrayNull(scip, &((*matrix)->vars));
885 
886  SCIPfreeBuffer(scip, matrix);
887  }
888 }
889 
890 /************** end matrix data structure and functions *************/
891 /********************************************************************/
892 
893 
894 /*
895  * Local methods
896  */
897 #ifdef SCIP_DEBUG
898 /** print a row from the constraint matrix */
899 static
900 void printRow(
901  SCIP* scip, /**< SCIP main data structure */
902  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
903  int row /**< row index for printing */
904  )
905 {
906  int* rowpnt;
907  int* rowend;
908  int c;
909  SCIP_Real val;
910  SCIP_Real* valpnt;
911  char relation;
912 
913  relation='-';
914  if( !matrix->isrhsinfinite &&
915  SCIPisEQ(scip, matrix->lhs[row], matrix->rhs[row]))
916  {
917  relation='e';
918  }
919  else if( !matrix->isrhsinfinite &&
920  !SCIPisEQ(scip, matrix->lhs[row], matrix->rhs[row]))
921  {
922  relation='r';
923  }
924  else
925  {
926  relation='g';
927  }
928 
929  rowpnt = matrix->rowmatind + matrix->rowmatbeg[row];
930  rowend = rowpnt + matrix->rowmatcnt[row];
931  valpnt = matrix->rowmatval + matrix->rowmatbeg[row];
932 
933  SCIPdebugPrintf("\n(L:%g) [%c] %g <=", (matrix->minactivityposinf[row] + matrix->minactivityneginf[row] > 0) ? -SCIPinfinity(scip) : matrix->minactivity[row], relation, matrix->lhs[row]);
934  for(; (rowpnt < rowend); rowpnt++, valpnt++)
935  {
936  c = *rowpnt;
937  val = *valpnt;
938  SCIPdebugPrintf(" %g{%s[idx:%d][bnd:%g,%g]}",
939  val, SCIPvarGetName(matrix->vars[c]), c,
940  SCIPvarGetLbGlobal(matrix->vars[c]), SCIPvarGetUbGlobal(matrix->vars[c]));
941  }
942  SCIPdebugPrintf(" <= %g (U:%g)", (matrix->maxactivityposinf[row] + matrix->maxactivityneginf[row] > 0) ? SCIPinfinity(scip) : matrix->rhs[row], matrix->maxactivity[row]);
943 }
944 
945 /** print all rows from the constraint matrix containing a variable */
946 static
947 SCIP_RETCODE printRowsOfCol(
948  SCIP* scip, /**< SCIP main data structure */
949  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
950  int col /**< column index for printing */
951  )
952 {
953  int numrows;
954  int* rows;
955  int i;
956  int* colpnt;
957  int* colend;
958 
959  numrows = matrix->colmatcnt[col];
960 
961  SCIP_CALL( SCIPallocBufferArray(scip, &rows, numrows) );
962 
963  colpnt = matrix->colmatind + matrix->colmatbeg[col];
964  colend = colpnt + matrix->colmatcnt[col];
965  for( i = 0; (colpnt < colend); colpnt++, i++ )
966  {
967  rows[i] = *colpnt;
968  }
969 
970  SCIPdebugPrintf("\n-------");
971  SCIPdebugPrintf("\ncol %d number rows: %d",col,numrows);
972  for( i = 0; i < numrows; i++ )
973  {
974  printRow(scip, matrix, rows[i]);
975  }
976  SCIPdebugPrintf("\n-------");
977 
978  SCIPfreeBufferArray(scip, &rows);
979 
980  return SCIP_OKAY;
981 }
982 
983 /** print information about a dominance relation */
984 static
985 SCIP_RETCODE printDomRelInfo(
986  SCIP* scip, /**< SCIP main data structure */
987  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
988  SCIP_VAR* dominatingvar, /**< dominating variable */
989  int dominatingidx, /**< index of dominating variable */
990  SCIP_VAR* dominatedvar, /**< dominated variable */
991  int dominatedidx, /**< index of dominated variable */
992  SCIP_Real dominatingub, /**< predicted upper bound of dominating variable */
993  SCIP_Real dominatingwclb /**< worst case lower bound of dominating variable */
994  )
995 {
996  char type;
997 
998  assert(SCIPvarGetType(dominatingvar)==SCIPvarGetType(dominatedvar));
999 
1000  switch(SCIPvarGetType(dominatingvar))
1001  {
1002  case SCIP_VARTYPE_CONTINUOUS:
1003  type='C';
1004  break;
1005  case SCIP_VARTYPE_BINARY:
1006  type='B';
1007  break;
1008  case SCIP_VARTYPE_INTEGER:
1009  case SCIP_VARTYPE_IMPLINT:
1010  type='I';
1011  break;
1012  default:
1013  SCIPerrorMessage("unknown variable type\n");
1014  SCIPABORT();
1015  return SCIP_INVALIDDATA; /*lint !e527*/
1016  }
1017 
1018  SCIPdebugPrintf("\n\n### [%c], obj:%g->%g,\t%s[idx:%d](nrows:%d)->%s[idx:%d](nrows:%d)\twclb=%g, ub'=%g, ub=%g",
1019  type, SCIPvarGetObj(dominatingvar), SCIPvarGetObj(dominatedvar),
1020  SCIPvarGetName(dominatingvar), dominatingidx, matrix->colmatcnt[dominatingidx],
1021  SCIPvarGetName(dominatedvar), dominatedidx, matrix->colmatcnt[dominatedidx],
1022  dominatingwclb, dominatingub, SCIPvarGetUbGlobal(dominatingvar));
1023 
1024  SCIP_CALL( printRowsOfCol(scip, matrix, dominatingidx) );
1025 
1026  return SCIP_OKAY;
1027 }
1028 #endif
1029 
1030 
1031 /** get minimum/maximum residual activity for the specified variable and with another variable set to its upper bound */
1032 static
1033 void getActivityResidualsUpperBound(
1034  SCIP* scip,
1035  CONSTRAINTMATRIX* matrix,
1036  int row,
1037  int col,
1038  SCIP_Real coef,
1039  int upperboundcol,
1040  SCIP_Real upperboundcoef,
1041  SCIP_Real* minresactivity,
1042  SCIP_Real* maxresactivity,
1043  SCIP_Bool* success
1044  )
1045 {
1046  SCIP_VAR* var;
1047  SCIP_VAR* upperboundvar;
1048  SCIP_Real lb;
1049  SCIP_Real ub;
1050  SCIP_Real lbupperboundvar;
1051  SCIP_Real ubupperboundvar;
1052  SCIP_Real tmpmaxact;
1053  SCIP_Real tmpminact;
1054  int maxactinf;
1055  int minactinf;
1056 
1057  assert(scip != NULL);
1058  assert(matrix != NULL);
1059  assert(row < matrix->nrows);
1060  assert(minresactivity != NULL);
1061  assert(maxresactivity != NULL);
1062 
1063  var = matrix->vars[col];
1064  upperboundvar = matrix->vars[upperboundcol];
1065  assert(var != NULL);
1066  assert(upperboundvar != NULL);
1067 
1068  lb = SCIPvarGetLbGlobal(var);
1069  ub = SCIPvarGetUbGlobal(var);
1070 
1071  maxactinf = matrix->maxactivityposinf[row] + matrix->maxactivityneginf[row];
1072  minactinf = matrix->minactivityposinf[row] + matrix->minactivityneginf[row];
1073 
1074  assert(!SCIPisInfinity(scip, lb));
1075  assert(!SCIPisInfinity(scip, -ub));
1076 
1077  lbupperboundvar = SCIPvarGetLbGlobal(upperboundvar);
1078  ubupperboundvar = SCIPvarGetUbGlobal(upperboundvar);
1079 
1080  assert(!SCIPisInfinity(scip, lbupperboundvar));
1081  assert(!SCIPisInfinity(scip, -ubupperboundvar));
1082 
1083  tmpminact = matrix->minactivity[row];
1084  tmpmaxact = matrix->maxactivity[row];
1085 
1086  /* first, adjust min and max activity such that upperboundvar is always set to its upper bound */
1087 
1088  /* abort if the upper bound is infty */
1089  if( SCIPisInfinity(scip, ubupperboundvar) )
1090  {
1091  *success = FALSE;
1092  return;
1093  }
1094  else
1095  {
1096  /* coef > 0 --> the lower bound is used for the minactivity --> update the minactivity */
1097  if( upperboundcoef > 0.0 )
1098  {
1099  if( SCIPisInfinity(scip, -lbupperboundvar) )
1100  {
1101  assert(minactinf > 0);
1102  minactinf--;
1103  tmpminact += (upperboundcoef * ubupperboundvar);
1104  }
1105  else
1106  {
1107  tmpminact = tmpminact - (upperboundcoef * lbupperboundvar) + (upperboundcoef * ubupperboundvar);
1108  }
1109  }
1110  /* coef < 0 --> the lower bound is used for the maxactivity --> update the maxactivity */
1111  else
1112  {
1113  if( SCIPisInfinity(scip, -lbupperboundvar) )
1114  {
1115  assert(maxactinf > 0);
1116  maxactinf--;
1117  tmpmaxact += (upperboundcoef * ubupperboundvar);
1118  }
1119  else
1120  {
1121  tmpmaxact = tmpmaxact - (upperboundcoef * lbupperboundvar) + (upperboundcoef * ubupperboundvar);
1122  }
1123  }
1124  }
1125 
1126  *success = TRUE;
1127 
1128  /* the coefficient is positive, so the upper bound contributed to the maxactivity and the lower bound to the minactivity */
1129  if( coef >= 0.0 )
1130  {
1131  if( SCIPisInfinity(scip, ub) )
1132  {
1133  assert(maxactinf >= 1);
1134  if( maxactinf == 1 )
1135  {
1136  *maxresactivity = tmpmaxact;
1137  }
1138  else
1139  *maxresactivity = SCIPinfinity(scip);
1140  }
1141  else
1142  {
1143  if( maxactinf > 0 )
1144  *maxresactivity = SCIPinfinity(scip);
1145  else
1146  *maxresactivity = tmpmaxact - coef * ub;
1147  }
1148 
1149  if( SCIPisInfinity(scip, -lb) )
1150  {
1151  assert(minactinf >= 1);
1152  if( minactinf == 1 )
1153  {
1154  *minresactivity = tmpminact;
1155  }
1156  else
1157  *minresactivity = -SCIPinfinity(scip);
1158  }
1159  else
1160  {
1161  if( minactinf > 0 )
1162  *minresactivity = -SCIPinfinity(scip);
1163  else
1164  *minresactivity = tmpminact - coef * lb;
1165  }
1166  }
1167  /* the coefficient is negative, so the lower bound contributed to the maxactivity and the upper bound to the minactivity */
1168  else
1169  {
1170  if( SCIPisInfinity(scip, -lb) )
1171  {
1172  assert(maxactinf >= 1);
1173  if( maxactinf == 1 )
1174  {
1175  *maxresactivity = tmpmaxact;
1176  }
1177  else
1178  *maxresactivity = SCIPinfinity(scip);
1179  }
1180  else
1181  {
1182  if( maxactinf > 0 )
1183  *maxresactivity = SCIPinfinity(scip);
1184  else
1185  *maxresactivity = tmpmaxact - coef * lb;
1186  }
1187 
1188  if( SCIPisInfinity(scip, ub) )
1189  {
1190  assert(minactinf >= 1);
1191  if( minactinf == 1 )
1192  {
1193  *minresactivity = tmpminact;
1194  }
1195  else
1196  *minresactivity = -SCIPinfinity(scip);
1197  }
1198  else
1199  {
1200  if( minactinf > 0 )
1201  *minresactivity = -SCIPinfinity(scip);
1202  else
1203  *minresactivity = tmpminact - coef * ub;
1204  }
1205  }
1206 }
1207 
1208 /** get minimum/maximum residual activity for the specified variable and with another variable set to its lower bound */
1209 static
1210 void getActivityResidualsLowerBound(
1211  SCIP* scip, /**< SCIP main data structure */
1212  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
1213  int row, /**< row index */
1214  int col, /**< column index */
1215  SCIP_Real coef, /**< coefficient of the column in this row */
1216  int lowerboundcol, /**< column index of variable to set to its lower bound */
1217  SCIP_Real lowerboundcoef, /**< coefficient in this row of the column to be set to its lower bound */
1218  SCIP_Real* minresactivity, /**< minimum residual activity of this row */
1219  SCIP_Real* maxresactivity, /**< maximum residual activity of this row */
1220  SCIP_Bool* success /**< pointer to store whether the computation was successful */
1221  )
1222 {
1223  SCIP_VAR* var;
1224  SCIP_VAR* lowerboundvar;
1225  SCIP_Real lb;
1226  SCIP_Real ub;
1227  SCIP_Real lblowerboundvar;
1228  SCIP_Real ublowerboundvar;
1229  SCIP_Real tmpmaxact;
1230  SCIP_Real tmpminact;
1231  int maxactinf;
1232  int minactinf;
1233 
1234  assert(scip != NULL);
1235  assert(matrix != NULL);
1236  assert(row < matrix->nrows);
1237  assert(minresactivity != NULL);
1238  assert(maxresactivity != NULL);
1239 
1240  var = matrix->vars[col];
1241  lowerboundvar = matrix->vars[lowerboundcol];
1242  assert(var != NULL);
1243  assert(lowerboundvar != NULL);
1244 
1245  lb = SCIPvarGetLbGlobal(var);
1246  ub = SCIPvarGetUbGlobal(var);
1247 
1248  maxactinf = matrix->maxactivityposinf[row] + matrix->maxactivityneginf[row];
1249  minactinf = matrix->minactivityposinf[row] + matrix->minactivityneginf[row];
1250 
1251  assert(!SCIPisInfinity(scip, lb));
1252  assert(!SCIPisInfinity(scip, -ub));
1253 
1254  lblowerboundvar = SCIPvarGetLbGlobal(lowerboundvar);
1255  ublowerboundvar = SCIPvarGetUbGlobal(lowerboundvar);
1256 
1257  assert(!SCIPisInfinity(scip, lblowerboundvar));
1258  assert(!SCIPisInfinity(scip, -ublowerboundvar));
1259 
1260  tmpminact = matrix->minactivity[row];
1261  tmpmaxact = matrix->maxactivity[row];
1262 
1263  /* first, adjust min and max activity such that lowerboundvar is always set to its lower bound */
1264 
1265  /* abort if the lower bound is -infty */
1266  if( SCIPisInfinity(scip, -lblowerboundvar) )
1267  {
1268  *success = FALSE;
1269  return;
1270  }
1271  else
1272  {
1273  /* coef > 0 --> the upper bound is used for the maxactivity --> update the maxactivity */
1274  if( lowerboundcoef > 0.0 )
1275  {
1276  if( SCIPisInfinity(scip, ublowerboundvar) )
1277  {
1278  assert(maxactinf > 0);
1279  maxactinf--;
1280  tmpmaxact += (lowerboundcoef * lblowerboundvar);
1281  }
1282  else
1283  {
1284  tmpmaxact = tmpmaxact - (lowerboundcoef * ublowerboundvar) + (lowerboundcoef * lblowerboundvar);
1285  }
1286  }
1287  /* coef < 0 --> the upper bound is used for the minactivity --> update the minactivity */
1288  else
1289  {
1290  if( SCIPisInfinity(scip, ublowerboundvar) )
1291  {
1292  assert(minactinf > 0);
1293  minactinf--;
1294  tmpminact += (lowerboundcoef * lblowerboundvar);
1295  }
1296  else
1297  {
1298  tmpminact = tmpminact - (lowerboundcoef * ublowerboundvar) + (lowerboundcoef * lblowerboundvar);
1299  }
1300  }
1301  }
1302 
1303  *success = TRUE;
1304 
1305  /* the coefficient is positive, so the upper bound contributed to the maxactivity and the lower bound to the minactivity */
1306  if( coef >= 0.0 )
1307  {
1308  if( SCIPisInfinity(scip, ub) )
1309  {
1310  assert(maxactinf >= 1);
1311  if( maxactinf == 1 )
1312  {
1313  *maxresactivity = tmpmaxact;
1314  }
1315  else
1316  *maxresactivity = SCIPinfinity(scip);
1317  }
1318  else
1319  {
1320  if( maxactinf > 0 )
1321  *maxresactivity = SCIPinfinity(scip);
1322  else
1323  *maxresactivity = tmpmaxact - coef * ub;
1324  }
1325 
1326  if( SCIPisInfinity(scip, -lb) )
1327  {
1328  assert(minactinf >= 1);
1329  if( minactinf == 1 )
1330  {
1331  *minresactivity = tmpminact;
1332  }
1333  else
1334  *minresactivity = -SCIPinfinity(scip);
1335  }
1336  else
1337  {
1338  if( minactinf > 0 )
1339  *minresactivity = -SCIPinfinity(scip);
1340  else
1341  *minresactivity = tmpminact - coef * lb;
1342  }
1343  }
1344  /* the coefficient is negative, so the lower bound contributed to the maxactivity and the upper bound to the minactivity */
1345  else
1346  {
1347  if( SCIPisInfinity(scip, -lb) )
1348  {
1349  assert(maxactinf >= 1);
1350  if( maxactinf == 1 )
1351  {
1352  *maxresactivity = tmpmaxact;
1353  }
1354  else
1355  *maxresactivity = SCIPinfinity(scip);
1356  }
1357  else
1358  {
1359  if( maxactinf > 0 )
1360  *maxresactivity = SCIPinfinity(scip);
1361  else
1362  *maxresactivity = tmpmaxact - coef * lb;
1363  }
1364 
1365  if( SCIPisInfinity(scip, ub) )
1366  {
1367  assert(minactinf >= 1);
1368  if( minactinf == 1 )
1369  {
1370  *minresactivity = tmpminact;
1371  }
1372  else
1373  *minresactivity = -SCIPinfinity(scip);
1374  }
1375  else
1376  {
1377  if( minactinf > 0 )
1378  *minresactivity = -SCIPinfinity(scip);
1379  else
1380  *minresactivity = tmpminact - coef * ub;
1381  }
1382  }
1383 }
1384 
1385 /** Calculate bounds of the dominated variable by rowbound analysis.
1386  * We use a special kind of predictive rowbound analysis by first setting the dominating variable to its upper bound.
1387  */
1388 static
1389 SCIP_RETCODE calcVarBoundsDominated(
1390  SCIP* scip, /**< SCIP main data structure */
1391  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
1392  int row, /**< current row index */
1393  int coldominating, /**< column index of dominating variable */
1394  SCIP_Real valdominating, /**< row coefficient of dominating variable */
1395  int coldominated, /**< column index of dominated variable */
1396  SCIP_Real valdominated, /**< row coefficient of dominated variable */
1397  SCIP_Bool* ubcalculated, /**< was a upper bound calculated? */
1398  SCIP_Real* calculatedub, /**< predicted upper bound */
1399  SCIP_Bool* wclbcalculated, /**< was a lower worst case bound calculated? */
1400  SCIP_Real* calculatedwclb, /**< predicted worst case lower bound */
1401  SCIP_Bool* lbcalculated, /**< was a lower bound calculated? */
1402  SCIP_Real* calculatedlb, /**< predicted lower bound */
1403  SCIP_Bool* wcubcalculated, /**< was a worst case upper bound calculated? */
1404  SCIP_Real* calculatedwcub /**< calculated worst case upper bound */
1405  )
1406 {
1407  SCIP_Real minresactivity;
1408  SCIP_Real maxresactivity;
1409  SCIP_Real lhs;
1410  SCIP_Real rhs;
1411  SCIP_Bool success;
1412 
1413  assert(scip != NULL);
1414  assert(matrix != NULL);
1415  assert(0 <= row && row < matrix->nrows);
1416  assert(0 <= coldominating && coldominating < matrix->ncols);
1417  assert(0 <= coldominated && coldominated < matrix->ncols);
1418 
1419  assert(ubcalculated != NULL);
1420  assert(calculatedub != NULL);
1421  assert(wclbcalculated != NULL);
1422  assert(calculatedwclb != NULL);
1423  assert(lbcalculated != NULL);
1424  assert(calculatedlb != NULL);
1425  assert(wcubcalculated != NULL);
1426  assert(calculatedwcub != NULL);
1427 
1428  assert(!SCIPisZero(scip, valdominated));
1429  assert(matrix->vars[coldominated] != NULL);
1430 
1431  *ubcalculated = FALSE;
1432  *wclbcalculated = FALSE;
1433  *lbcalculated = FALSE;
1434  *wcubcalculated = FALSE;
1435 
1436  /* no rowbound analysis for multiaggregated variables, which should not exist, because the matrix only consists of
1437  * active variables
1438  */
1439  assert(SCIPvarGetStatus(matrix->vars[coldominating]) != SCIP_VARSTATUS_MULTAGGR);
1440  assert(SCIPvarGetStatus(matrix->vars[coldominated]) != SCIP_VARSTATUS_MULTAGGR);
1441 
1442  lhs = matrix->lhs[row];
1443  rhs = matrix->rhs[row];
1444  assert(!SCIPisInfinity(scip, lhs));
1445  assert(!SCIPisInfinity(scip, -rhs));
1446 
1447  /* get minimum/maximum activity of this row without the dominated variable */
1448  getActivityResidualsUpperBound(scip, matrix, row, coldominated, valdominated, coldominating, valdominating,
1449  &minresactivity, &maxresactivity, &success);
1450 
1451  if( !success )
1452  return SCIP_OKAY;
1453 
1454  assert(!SCIPisInfinity(scip, minresactivity));
1455  assert(!SCIPisInfinity(scip, -maxresactivity));
1456 
1457  *calculatedub = SCIPinfinity(scip);
1458  *calculatedwclb = -SCIPinfinity(scip);
1459  *calculatedlb = -SCIPinfinity(scip);
1460  *calculatedwcub = SCIPinfinity(scip);
1461 
1462  /* predictive rowbound analysis */
1463 
1464  if( valdominated > 0.0 )
1465  {
1466  /* lower bound calculation */
1467  if( !SCIPisInfinity(scip, maxresactivity) )
1468  {
1469  *calculatedlb = (lhs - maxresactivity)/valdominated;
1470  *lbcalculated = TRUE;
1471  }
1472 
1473  /* worst case calculation of lower bound */
1474  if( !SCIPisInfinity(scip, -minresactivity) )
1475  {
1476  *calculatedwclb = (lhs - minresactivity)/valdominated;
1477  *wclbcalculated = TRUE;
1478  }
1479  else
1480  {
1481  /* worst case lower bound is infinity */
1482  *calculatedwclb = SCIPinfinity(scip);
1483  *wclbcalculated = TRUE;
1484  }
1485 
1486  /* we can only calculate upper bounds, of the right hand side is finite */
1487  if( !matrix->isrhsinfinite[row] )
1488  {
1489  /* upper bound calculation */
1490  if( !SCIPisInfinity(scip, -minresactivity) )
1491  {
1492  *calculatedub = (rhs - minresactivity)/valdominated;
1493  *ubcalculated = TRUE;
1494  }
1495 
1496  /* worst case calculation of upper bound */
1497  if( !SCIPisInfinity(scip, maxresactivity) )
1498  {
1499  *calculatedwcub = (rhs - maxresactivity)/valdominated;
1500  *wcubcalculated = TRUE;
1501  }
1502  else
1503  {
1504  /* worst case upper bound is -infinity */
1505  *calculatedwcub = -SCIPinfinity(scip);
1506  *wcubcalculated = TRUE;
1507  }
1508  }
1509  }
1510  else
1511  {
1512  /* upper bound calculation */
1513  if( !SCIPisInfinity(scip, maxresactivity) )
1514  {
1515  *calculatedub = (lhs - maxresactivity)/valdominated;
1516  *ubcalculated = TRUE;
1517  }
1518 
1519  /* worst case calculation of upper bound */
1520  if( !SCIPisInfinity(scip, -minresactivity) )
1521  {
1522  *calculatedwcub = (lhs - minresactivity)/valdominated;
1523  *wcubcalculated = TRUE;
1524  }
1525  else
1526  {
1527  /* worst case upper bound is -infinity */
1528  *calculatedwcub = -SCIPinfinity(scip);
1529  *wcubcalculated = TRUE;
1530  }
1531 
1532  /* we can only calculate lower bounds, of the right hand side is finite */
1533  if( !matrix->isrhsinfinite[row] )
1534  {
1535  /* lower bound calculation */
1536  if( !SCIPisInfinity(scip, -minresactivity) )
1537  {
1538  *calculatedlb = (rhs - minresactivity)/valdominated;
1539  *lbcalculated = TRUE;
1540  }
1541 
1542  /* worst case calculation of lower bound */
1543  if( !SCIPisInfinity(scip, maxresactivity) )
1544  {
1545  *calculatedwclb = (rhs - maxresactivity)/valdominated;
1546  *wclbcalculated = TRUE;
1547  }
1548  else
1549  {
1550  /* worst case lower bound is infinity */
1551  *calculatedwclb = SCIPinfinity(scip);
1552  *wclbcalculated = TRUE;
1553  }
1554  }
1555  }
1556 
1557  return SCIP_OKAY;
1558 }
1559 
1560 /** Calculate bounds of the dominating variable by rowbound analysis.
1561  * We use a special kind of predictive rowbound analysis by first setting the dominated variable to its lower bound.
1562  */
1563 static
1564 SCIP_RETCODE calcVarBoundsDominating(
1565  SCIP* scip, /**< SCIP main data structure */
1566  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
1567  int row, /**< current row index */
1568  int coldominating, /**< column index of dominating variable */
1569  SCIP_Real valdominating, /**< row coefficient of dominating variable */
1570  int coldominated, /**< column index of dominated variable */
1571  SCIP_Real valdominated, /**< row coefficient of dominated variable */
1572  SCIP_Bool* ubcalculated, /**< was a upper bound calculated? */
1573  SCIP_Real* calculatedub, /**< predicted upper bound */
1574  SCIP_Bool* wclbcalculated, /**< was a lower worst case bound calculated? */
1575  SCIP_Real* calculatedwclb, /**< predicted worst case lower bound */
1576  SCIP_Bool* lbcalculated, /**< was a lower bound calculated? */
1577  SCIP_Real* calculatedlb, /**< predicted lower bound */
1578  SCIP_Bool* wcubcalculated, /**< was a worst case upper bound calculated? */
1579  SCIP_Real* calculatedwcub /**< calculated worst case upper bound */
1580  )
1581 {
1582  SCIP_Real minresactivity;
1583  SCIP_Real maxresactivity;
1584  SCIP_Real lhs;
1585  SCIP_Real rhs;
1586  SCIP_Bool success;
1587 
1588  assert(scip != NULL);
1589  assert(matrix != NULL);
1590  assert(0 <= row && row < matrix->nrows);
1591  assert(0 <= coldominating && coldominating < matrix->ncols);
1592  assert(0 <= coldominated && coldominated < matrix->ncols);
1593 
1594  assert(ubcalculated != NULL);
1595  assert(calculatedub != NULL);
1596  assert(wclbcalculated != NULL);
1597  assert(calculatedwclb != NULL);
1598  assert(lbcalculated != NULL);
1599  assert(calculatedlb != NULL);
1600  assert(wcubcalculated != NULL);
1601  assert(calculatedwcub != NULL);
1602 
1603  assert(!SCIPisZero(scip, valdominating));
1604  assert(matrix->vars[coldominating] != NULL);
1605 
1606  *ubcalculated = FALSE;
1607  *wclbcalculated = FALSE;
1608  *lbcalculated = FALSE;
1609  *wcubcalculated = FALSE;
1610 
1611  /* no rowbound analysis for multiaggregated variables, which should not exist, because the matrix only consists of
1612  * active variables
1613  */
1614  assert(SCIPvarGetStatus(matrix->vars[coldominating]) != SCIP_VARSTATUS_MULTAGGR);
1615  assert(SCIPvarGetStatus(matrix->vars[coldominated]) != SCIP_VARSTATUS_MULTAGGR);
1616 
1617  lhs = matrix->lhs[row];
1618  rhs = matrix->rhs[row];
1619  assert(!SCIPisInfinity(scip, lhs));
1620  assert(!SCIPisInfinity(scip, -rhs));
1621 
1622  /* get minimum/maximum activity of this row without the dominating variable */
1623  getActivityResidualsLowerBound(scip, matrix, row, coldominating, valdominating, coldominated, valdominated,
1624  &minresactivity, &maxresactivity, &success);
1625 
1626  if( !success )
1627  return SCIP_OKAY;
1628 
1629  assert(!SCIPisInfinity(scip, minresactivity));
1630  assert(!SCIPisInfinity(scip, -maxresactivity));
1631 
1632  *calculatedub = SCIPinfinity(scip);
1633  *calculatedwclb = -SCIPinfinity(scip);
1634  *calculatedlb = -SCIPinfinity(scip);
1635  *calculatedwcub = SCIPinfinity(scip);
1636 
1637  /* predictive rowbound analysis */
1638 
1639  if( valdominating > 0.0 )
1640  {
1641  /* lower bound calculation */
1642  if( !SCIPisInfinity(scip, maxresactivity) )
1643  {
1644  *calculatedlb = (lhs - maxresactivity)/valdominating;
1645  *lbcalculated = TRUE;
1646  }
1647 
1648  /* worst case calculation of lower bound */
1649  if( !SCIPisInfinity(scip, -minresactivity) )
1650  {
1651  *calculatedwclb = (lhs - minresactivity)/valdominating;
1652  *wclbcalculated = TRUE;
1653  }
1654  else
1655  {
1656  /* worst case lower bound is infinity */
1657  *calculatedwclb = SCIPinfinity(scip);
1658  *wclbcalculated = TRUE;
1659  }
1660 
1661  /* we can only calculate upper bounds, of the right hand side is finite */
1662  if( !matrix->isrhsinfinite[row] )
1663  {
1664  /* upper bound calculation */
1665  if( !SCIPisInfinity(scip, -minresactivity) )
1666  {
1667  *calculatedub = (rhs - minresactivity)/valdominating;
1668  *ubcalculated = TRUE;
1669  }
1670 
1671  /* worst case calculation of upper bound */
1672  if( !SCIPisInfinity(scip, maxresactivity) )
1673  {
1674  *calculatedwcub = (rhs - maxresactivity)/valdominating;
1675  *wcubcalculated = TRUE;
1676  }
1677  else
1678  {
1679  /* worst case upper bound is -infinity */
1680  *calculatedwcub = -SCIPinfinity(scip);
1681  *wcubcalculated = TRUE;
1682  }
1683  }
1684  }
1685  else
1686  {
1687  /* upper bound calculation */
1688  if( !SCIPisInfinity(scip, maxresactivity) )
1689  {
1690  *calculatedub = (lhs - maxresactivity)/valdominating;
1691  *ubcalculated = TRUE;
1692  }
1693 
1694  /* worst case calculation of upper bound */
1695  if( !SCIPisInfinity(scip, -minresactivity) )
1696  {
1697  *calculatedwcub = (lhs - minresactivity)/valdominating;
1698  *wcubcalculated = TRUE;
1699  }
1700  else
1701  {
1702  /* worst case upper bound is -infinity */
1703  *calculatedwcub = -SCIPinfinity(scip);
1704  *wcubcalculated = TRUE;
1705  }
1706 
1707  /* we can only calculate lower bounds, of the right hand side is finite */
1708  if( !matrix->isrhsinfinite[row] )
1709  {
1710  /* lower bound calculation */
1711  if( !SCIPisInfinity(scip, -minresactivity) )
1712  {
1713  *calculatedlb = (rhs - minresactivity)/valdominating;
1714  *lbcalculated = TRUE;
1715  }
1716 
1717  /* worst case calculation of lower bound */
1718  if( !SCIPisInfinity(scip, maxresactivity) )
1719  {
1720  *calculatedwclb = (rhs - maxresactivity)/valdominating;
1721  *wclbcalculated = TRUE;
1722  }
1723  else
1724  {
1725  /* worst case lower bound is infinity */
1726  *calculatedwclb = SCIPinfinity(scip);
1727  *wclbcalculated = TRUE;
1728  }
1729  }
1730  }
1731 
1732  return SCIP_OKAY;
1733 }
1734 
1735 /** try to find new variable bounds and update them when they are better then the old bounds */
1736 static
1737 SCIP_RETCODE updateBounds(
1738  SCIP* scip, /**< SCIP main data structure */
1739  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
1740  int row, /**< row index */
1741  int col1, /**< dominating variable index */
1742  SCIP_Real val1, /**< dominating variable coefficient */
1743  int col2, /**< dominated variable index */
1744  SCIP_Real val2, /**< dominated variable coefficient */
1745  SCIP_Bool predictdominating, /**< flag indicating if bounds of dominating or dominated variable are predicted */
1746  SCIP_Real* upperbound, /**< predicted upper bound */
1747  SCIP_Real* wclowerbound, /**< predicted worst case lower bound */
1748  SCIP_Real* lowerbound, /**< predicted lower bound */
1749  SCIP_Real* wcupperbound /**< predicted worst case upper bound */
1750  )
1751 {
1752  SCIP_Bool ubcalculated;
1753  SCIP_Bool wclbcalculated;
1754  SCIP_Bool lbcalculated;
1755  SCIP_Bool wcubcalculated;
1756  SCIP_Real newub;
1757  SCIP_Real newwclb;
1758  SCIP_Real newlb;
1759  SCIP_Real newwcub;
1760 
1761  assert(scip != NULL);
1762  assert(matrix != NULL);
1763  assert(row < matrix->nrows);
1764  assert(col1 < matrix->ncols);
1765  assert(col2 < matrix->ncols);
1766  assert(upperbound != NULL);
1767  assert(wclowerbound != NULL);
1768  assert(lowerbound != NULL);
1769  assert(wcupperbound != NULL);
1770 
1771  newub = SCIPinfinity(scip);
1772  newlb = -SCIPinfinity(scip);
1773  newwcub = newub;
1774  newwclb = newlb;
1775 
1776  if( predictdominating )
1777  {
1778  /* do predictive rowbound analysis for the dominating variable */
1779  SCIP_CALL( calcVarBoundsDominating(scip, matrix, row, col1, val1, col2, val2,
1780  &ubcalculated, &newub, &wclbcalculated, &newwclb,
1781  &lbcalculated, &newlb, &wcubcalculated, &newwcub) );
1782  }
1783  else
1784  {
1785  /* do predictive rowbound analysis for the dominated variable */
1786  SCIP_CALL( calcVarBoundsDominated(scip, matrix, row, col1, val1, col2, val2,
1787  &ubcalculated, &newub, &wclbcalculated, &newwclb,
1788  &lbcalculated, &newlb, &wcubcalculated, &newwcub) );
1789  }
1790 
1791  /* update bounds in case if they are better */
1792  if( ubcalculated )
1793  {
1794  if( newub < *upperbound )
1795  *upperbound = newub;
1796  }
1797  if( wclbcalculated )
1798  {
1799  if( newwclb > *wclowerbound )
1800  *wclowerbound = newwclb;
1801  }
1802  if( lbcalculated )
1803  {
1804  if( newlb > *lowerbound )
1805  *lowerbound = newlb;
1806  }
1807  if( wcubcalculated )
1808  {
1809  if( newwcub < *wcupperbound )
1810  *wcupperbound = newwcub;
1811  }
1812 
1813  return SCIP_OKAY;
1814 }
1815 
1816 /** detect parallel columns by using the algorithm of Bixby and Wagner
1817  * see paper: "A note on Detecting Simple Redundancies in Linear Systems", June 1986
1818  */
1819 static
1820 SCIP_RETCODE detectParallelCols(
1821  SCIP* scip, /**< SCIP main data structure */
1822  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
1823  int* pclass, /**< parallel column classes */
1824  SCIP_Bool* varineq /**< indicating if variable is within an equation */
1825  )
1826 {
1827  SCIP_Real* valpnt;
1828  SCIP_Real* values;
1829  SCIP_Real* scale;
1830  int* classsizes;
1831  int* pcset;
1832  int* rowpnt;
1833  int* rowend;
1834  int* colindices;
1835  int* pcs;
1836  SCIP_Real startval;
1837  SCIP_Real aij;
1838  int startpc;
1839  int startk;
1840  int startt;
1841  int pcsetfill;
1842  int colidx;
1843  int k;
1844  int t;
1845  int m;
1846  int i;
1847  int r;
1848  int newpclass;
1849  int pc;
1850 
1851  assert(scip != NULL);
1852  assert(matrix != NULL);
1853  assert(pclass != NULL);
1854  assert(varineq != NULL);
1855 
1856  SCIP_CALL( SCIPallocBufferArray(scip, &classsizes, matrix->ncols) );
1857  SCIP_CALL( SCIPallocBufferArray(scip, &scale, matrix->ncols) );
1858  SCIP_CALL( SCIPallocBufferArray(scip, &pcset, matrix->ncols) );
1859  SCIP_CALL( SCIPallocBufferArray(scip, &values, matrix->ncols) );
1860  SCIP_CALL( SCIPallocBufferArray(scip, &colindices, matrix->ncols) );
1861  SCIP_CALL( SCIPallocBufferArray(scip, &pcs, matrix->ncols) );
1862 
1863  /* init */
1864  BMSclearMemoryArray(scale, matrix->ncols);
1865  BMSclearMemoryArray(pclass, matrix->ncols);
1866  BMSclearMemoryArray(classsizes, matrix->ncols);
1867  classsizes[0] = matrix->ncols;
1868  pcsetfill = 0;
1869  for( t = 1; t < matrix->ncols; ++t )
1870  pcset[pcsetfill++] = t;
1871 
1872  /* loop over all rows */
1873  for( r = 0; r < matrix->nrows; ++r )
1874  {
1875  /* we consider only equations or ranged rows */
1876  if( !matrix->isrhsinfinite[r] )
1877  {
1878  rowpnt = matrix->rowmatind + matrix->rowmatbeg[r];
1879  rowend = rowpnt + matrix->rowmatcnt[r];
1880  valpnt = matrix->rowmatval + matrix->rowmatbeg[r];
1881 
1882  i = 0;
1883  for( ; (rowpnt < rowend); rowpnt++, valpnt++ )
1884  {
1885  aij = *valpnt;
1886  colidx = *rowpnt;
1887 
1888 #ifdef SCIP_DEBUG
1889  if( SCIPisZero(scip, aij) )
1890  {
1891  SCIPdebugMessage("Matrix coefficient is very small !\n");
1892  }
1893 #endif
1894  /* remember variable was part of an equation or ranged row */
1895  varineq[colidx] = TRUE;
1896 
1897  if( scale[colidx] == 0.0 )
1898  scale[colidx] = aij;
1899  assert(scale[colidx] != 0.0);
1900 
1901  colindices[i] = colidx;
1902  values[i] = aij / scale[colidx];
1903  pc = pclass[colidx];
1904  assert(pc < matrix->ncols);
1905 
1906  /* update class sizes and pclass set */
1907  assert(classsizes[pc] > 0);
1908  classsizes[pc]--;
1909  if( classsizes[pc] == 0 )
1910  {
1911  assert(pcsetfill < matrix->ncols);
1912  pcset[pcsetfill++] = pc;
1913  }
1914  pcs[i] = pc;
1915 
1916  i++;
1917  }
1918 
1919  /* sort on the pclass values */
1920  if( i > 1 )
1921  {
1922  SCIPsortIntIntReal(pcs, colindices, values, i);
1923  }
1924 
1925  k = 0;
1926  while( TRUE ) /*lint !e716*/
1927  {
1928  assert(k < i);
1929  startpc = pcs[k];
1930  startk = k;
1931 
1932  /* find pclass-sets */
1933  while( k < i && pcs[k] == startpc )
1934  k++;
1935 
1936  /* sort on the A values which have equal pclass values */
1937  if( k - startk > 1 )
1938  SCIPsortRealInt(&(values[startk]), &(colindices[startk]), k - startk);
1939 
1940  t = 0;
1941  while( TRUE ) /*lint !e716*/
1942  {
1943  assert(startk + t < i);
1944  startval = values[startk + t];
1945  startt = t;
1946 
1947  /* find A-sets */
1948  while( t < k - startk && SCIPisEQ(scip, startval, values[startk + t]) )
1949  t++;
1950 
1951  /* get new pclass */
1952  newpclass = pcset[0];
1953  assert(pcsetfill > 0);
1954  pcset[0] = pcset[--pcsetfill];
1955 
1956  /* renumbering */
1957  for( m = startk + startt; m < startk + t; m++ )
1958  {
1959  assert(m < i);
1960  assert(colindices[m] < matrix->ncols);
1961  assert(newpclass < matrix->ncols);
1962 
1963  pclass[colindices[m]] = newpclass;
1964  classsizes[newpclass]++;
1965  }
1966 
1967  if( t == k - startk )
1968  break;
1969  }
1970 
1971  if( k == matrix->rowmatcnt[r] )
1972  break;
1973  }
1974  }
1975  }
1976 
1977  SCIPfreeBufferArray(scip, &pcs);
1978  SCIPfreeBufferArray(scip, &colindices);
1979  SCIPfreeBufferArray(scip, &values);
1980  SCIPfreeBufferArray(scip, &pcset);
1981  SCIPfreeBufferArray(scip, &scale);
1982  SCIPfreeBufferArray(scip, &classsizes);
1983 
1984  return SCIP_OKAY;
1985 }
1986 
1987 
1988 /** try to improve variable bounds by predictive bound strengthening */
1989 static
1990 SCIP_RETCODE predBndStr(
1991  SCIP* scip, /**< SCIP main data structure */
1992  SCIP_VAR* dominatingvar, /**< dominating variable */
1993  int dominatingidx, /**< column index of the dominating variable */
1994  SCIP_Real dominatingub, /**< predicted upper bound of the dominating variable */
1995  SCIP_Real dominatinglb, /**< predicted lower bound of the dominating variable */
1996  SCIP_Real dominatingwcub, /**< predicted worst case upper bound of the dominating variable */
1997  SCIP_VAR* dominatedvar, /**< dominated variable */
1998  int dominatedidx, /**< column index of the dominated variable */
1999  SCIP_Real dominatedub, /**< predicted upper bound of the dominated variable */
2000  SCIP_Real dominatedwclb, /**< predicted worst case lower bound of the dominated variable */
2001  SCIP_Real dominatedlb, /**< predicted lower bound of the dominated variable */
2002  FIXINGDIRECTION* varstofix, /**< array holding fixing information */
2003  int* nchgbds /**< count number of bound changes */
2004  )
2005 {
2006  /* we compare only variables from the same type */
2007  if( !(SCIPvarGetType(dominatingvar) == SCIPvarGetType(dominatedvar) ||
2008  SCIPvarIsBinary(dominatingvar) == SCIPvarIsBinary(dominatedvar) ||
2009  (SCIPvarGetType(dominatingvar) == SCIP_VARTYPE_INTEGER && SCIPvarGetType(dominatedvar) == SCIP_VARTYPE_IMPLINT) ||
2010  (SCIPvarGetType(dominatedvar) == SCIP_VARTYPE_INTEGER && SCIPvarGetType(dominatingvar) == SCIP_VARTYPE_IMPLINT)) )
2011  {
2012  return SCIP_OKAY;
2013  }
2014 
2015  if( varstofix[dominatingidx] == NOFIX )
2016  {
2017  /* assume x dominates y (x->y). we get this bound from a positive coefficient
2018  * of the dominating variable. because x->y the dominated variable y has
2019  * a positive coefficient too. thus y contributes to the minactivity with its
2020  * lower bound. but this case is considered within predictive bound analysis.
2021  * thus the dominating upper bound is a common upper bound.
2022  */
2023  if( !SCIPisInfinity(scip, dominatingub) )
2024  {
2025  SCIP_Real newub;
2026  SCIP_Real oldub;
2027  SCIP_Real lb;
2028 
2029  newub = dominatingub;
2030  oldub = SCIPvarGetUbGlobal(dominatingvar);
2031  lb = SCIPvarGetLbGlobal(dominatingvar);
2032 
2033  /* if( SCIPvarGetType(dominatingvar) != SCIP_VARTYPE_CONTINUOUS )
2034  newub = SCIPceil(scip, newub); */
2035 
2036  if( SCIPisLE(scip, lb, newub) && SCIPisLT(scip, newub, oldub) )
2037  {
2038  SCIPdebugMessage("[ub]\tupper bound for dominating variable <%s> changed: [%.17f,%.17f] -> [%.17f,%.17f]\n",
2039  SCIPvarGetName(dominatingvar), lb, oldub, lb, newub);
2040  SCIP_CALL( SCIPchgVarUb(scip, dominatingvar, newub) );
2041  (*nchgbds)++;
2042  }
2043  }
2044 
2045  /* assume x dominates y (x->y). we get this lower bound of the dominating variable
2046  * from a negative coefficient within a <= relation. if y has a positive coefficient
2047  * we get a common lower bound of x from predictive bound analysis. if y has a
2048  * negative coefficient we get an improved lower bound of x because the minactivity
2049  * is greater. for discrete variables we have to round down the lower bound.
2050  */
2051  if( !SCIPisInfinity(scip, -dominatinglb) )
2052  {
2053  SCIP_Real newlb;
2054  SCIP_Real oldlb;
2055  SCIP_Real ub;
2056 
2057  newlb = dominatinglb;
2058  oldlb = SCIPvarGetLbGlobal(dominatingvar);
2059  ub = SCIPvarGetUbGlobal(dominatingvar);
2060 
2061  if( SCIPvarGetType(dominatingvar) != SCIP_VARTYPE_CONTINUOUS )
2062  newlb = SCIPfloor(scip, newlb);
2063 
2064  if( SCIPisLT(scip, oldlb, newlb) && SCIPisLE(scip, newlb, ub) )
2065  {
2066  SCIPdebugMessage("[lb]\tlower bound of dominating variable <%s> changed: [%.17f,%.17f] -> [%.17f,%.17f]\n",
2067  SCIPvarGetName(dominatingvar), oldlb, ub, newlb, ub);
2068  SCIP_CALL( SCIPchgVarLb(scip, dominatingvar, newlb) );
2069  (*nchgbds)++;
2070  }
2071  }
2072 
2073  /* assume x dominates y (x->y). we get this bound from a positive coefficient
2074  * of x within a <= relation. from x->y it follows, that y has a positive
2075  * coefficient in this row too. the worst case upper bound of x is an estimation
2076  * of how great x can be in every case. if the objective coefficient of x is
2077  * negative we get thus a lower bound of x. for discrete variables we have
2078  * to round down the lower bound before.
2079  */
2080  if( !SCIPisInfinity(scip, dominatingwcub) && SCIPisNegative(scip, SCIPvarGetObj(dominatingvar)))
2081  {
2082  SCIP_Real newlb;
2083  SCIP_Real oldlb;
2084  SCIP_Real ub;
2085 
2086  newlb = dominatingwcub;
2087  oldlb = SCIPvarGetLbGlobal(dominatingvar);
2088  ub = SCIPvarGetUbGlobal(dominatingvar);
2089 
2090  if( SCIPvarGetType(dominatingvar) != SCIP_VARTYPE_CONTINUOUS )
2091  newlb = SCIPfloor(scip, newlb);
2092 
2093  if( SCIPisLT(scip, oldlb, newlb) && SCIPisLE(scip, newlb, ub) )
2094  {
2095  SCIPdebugMessage("[wcub]\tlower bound of dominating variable <%s> changed: [%.17f,%.17f] -> [%.17f,%.17f]\n",
2096  SCIPvarGetName(dominatingvar), oldlb, ub, newlb, ub);
2097  SCIP_CALL( SCIPchgVarLb(scip, dominatingvar, newlb) );
2098  (*nchgbds)++;
2099  }
2100  }
2101  }
2102 
2103  if( varstofix[dominatedidx] == NOFIX )
2104  {
2105  /* assume x dominates y (x->y). we get this bound for a positive coefficient of y
2106  * within a <= relation. if x has a negative coefficient we get a common upper
2107  * bound of y. if x has a positive coefficient we get an improved upper bound
2108  * of y because the minactivity is greater.
2109  */
2110  if( !SCIPisInfinity(scip, dominatedub) )
2111  {
2112  SCIP_Real newub;
2113  SCIP_Real oldub;
2114  SCIP_Real lb;
2115 
2116  newub = dominatedub;
2117  oldub = SCIPvarGetUbGlobal(dominatedvar);
2118  lb = SCIPvarGetLbGlobal(dominatedvar);
2119 
2120  if( SCIPisLE(scip, lb, newub) && SCIPisLT(scip, newub, oldub) )
2121  {
2122  SCIPdebugMessage("[ub]\tupper bound of dominated variable <%s> changed: [%.17f,%.17f] -> [%.17f,%.17f]\n",
2123  SCIPvarGetName(dominatedvar), lb, oldub, lb, newub);
2124  SCIP_CALL( SCIPchgVarUb(scip, dominatedvar, newub) );
2125  (*nchgbds)++;
2126  }
2127  }
2128 
2129  /* assume x dominates y (x->y). we get this bound only from a negative
2130  * coefficient of y within a <= relation. because of x->y then x has a negative
2131  * coefficient too. the worst case lower bound is an estimation what property
2132  * the dominated variable must have if the dominating variable is at its upper bound.
2133  * to get an upper bound of the dominated variable we need to consider a positive
2134  * objective coefficient. for discrete variables we have to round up the upper bound.
2135  */
2136  if( !SCIPisInfinity(scip, -dominatedwclb) && SCIPisPositive(scip, SCIPvarGetObj(dominatedvar)) )
2137  {
2138  SCIP_Real newub;
2139  SCIP_Real oldub;
2140  SCIP_Real lb;
2141 
2142  newub = dominatedwclb;
2143  oldub = SCIPvarGetUbGlobal(dominatedvar);
2144  lb = SCIPvarGetLbGlobal(dominatedvar);
2145 
2146  if( SCIPvarGetType(dominatedvar) != SCIP_VARTYPE_CONTINUOUS )
2147  newub = SCIPceil(scip, newub);
2148 
2149  if( SCIPisLE(scip, lb, newub) && SCIPisLT(scip, newub, oldub) )
2150  {
2151  SCIPdebugMessage("[wclb]\tupper bound of dominated variable <%s> changed: [%.17f,%.17f] -> [%.17f,%.17f]\n",
2152  SCIPvarGetName(dominatedvar), lb, oldub, lb, newub);
2153  SCIP_CALL( SCIPchgVarUb(scip, dominatedvar, newub) );
2154  (*nchgbds)++;
2155  }
2156  }
2157 
2158  /* assume x dominates y (x->y). we get a lower bound of y from a negative
2159  * coefficient within a <= relation. but if x->y then x has a negative
2160  * coefficient too and contributes with its upper bound to the minactivity.
2161  * thus in all we have a common lower bound calculation and no rounding is necessary.
2162  */
2163  if( !SCIPisInfinity(scip, -dominatedlb) )
2164  {
2165  SCIP_Real newlb;
2166  SCIP_Real oldlb;
2167  SCIP_Real ub;
2168 
2169  newlb = dominatedlb;
2170  oldlb = SCIPvarGetLbGlobal(dominatedvar);
2171  ub = SCIPvarGetUbGlobal(dominatedvar);
2172 
2173  if( SCIPisLT(scip, oldlb, newlb) && SCIPisLE(scip, newlb, ub) )
2174  {
2175  SCIPdebugMessage("[lb]\tlower bound of dominated variable <%s> changed: [%.17f,%.17f] -> [%.17f,%.17f]\n",
2176  SCIPvarGetName(dominatedvar), oldlb, ub, newlb, ub);
2177  SCIP_CALL( SCIPchgVarLb(scip, dominatedvar, newlb) );
2178  (*nchgbds)++;
2179  }
2180  }
2181  }
2182 
2183  return SCIP_OKAY;
2184 }
2185 
2186 /** try to find variable fixings */
2187 static
2188 SCIP_RETCODE findFixings(
2189  SCIP* scip, /**< SCIP main data structure */
2190  CONSTRAINTMATRIX* matrix, /**< constraint matrix structure */
2191  SCIP_VAR* dominatingvar, /**< dominating variable */
2192  int dominatingidx, /**< column index of the dominating variable */
2193  SCIP_Real dominatingub, /**< predicted upper bound of the dominating variable */
2194  SCIP_Real dominatingwclb, /**< predicted worst case lower bound of the dominating variable */
2195  SCIP_Real dominatinglb, /**< predicted lower bound of the dominating variable */
2196  SCIP_Real dominatingwcub, /**< predicted worst case upper bound of the dominating variable */
2197  SCIP_VAR* dominatedvar, /**< dominated variable */
2198  int dominatedidx, /**< column index of the dominated variable */
2199  FIXINGDIRECTION* varstofix, /**< array holding fixing information */
2200  SCIP_Bool onlybinvars, /**< flag indicating only binary variables are present */
2201  SCIP_Bool onlyoneone, /**< when onlybinvars is TRUE, flag indicates if both binary variables are in clique */
2202  int* nfixings /**< counter for possible fixings */
2203  )
2204 {
2205  /* we compare only variables from the same type */
2206  if( !(SCIPvarGetType(dominatingvar) == SCIPvarGetType(dominatedvar) ||
2207  SCIPvarIsBinary(dominatingvar) == SCIPvarIsBinary(dominatedvar) ||
2208  (SCIPvarGetType(dominatingvar) == SCIP_VARTYPE_INTEGER && SCIPvarGetType(dominatedvar) == SCIP_VARTYPE_IMPLINT) ||
2209  (SCIPvarGetType(dominatedvar) == SCIP_VARTYPE_INTEGER && SCIPvarGetType(dominatingvar) == SCIP_VARTYPE_IMPLINT)) )
2210  {
2211  return SCIP_OKAY;
2212  }
2213 
2214  if( varstofix[dominatedidx] == NOFIX && matrix->colmatcnt[dominatingidx] == 1
2215  && matrix->colmatcnt[dominatedidx] == 1 )
2216  {
2217  /* We have a x->y dominance relation and only one equality constraint
2218  * where the dominating variable x with an infinity upper bound and the
2219  * dominated variable y are present. Then the dominated variable y
2220  * can be fixed at its lower bound.
2221  */
2222  int row;
2223  row = *(matrix->colmatind + matrix->colmatbeg[dominatedidx]);
2224 
2225  if( SCIPisEQ(scip, matrix->lhs[row], matrix->rhs[row]) &&
2226  SCIPisInfinity(scip, SCIPvarGetUbGlobal(dominatingvar)) )
2227  {
2228  varstofix[dominatedidx] = FIXATLB;
2229  (*nfixings)++;
2230 
2231  return SCIP_OKAY;
2232  }
2233  }
2234 
2235  if( varstofix[dominatedidx] == NOFIX && !SCIPisNegative(scip, SCIPvarGetObj(dominatedvar)) )
2236  {
2237  if( !SCIPisInfinity(scip, -dominatingwclb) &&
2238  SCIPisLE(scip, dominatingwclb, SCIPvarGetUbGlobal(dominatingvar)) )
2239  {
2240  /* we have a x->y dominance relation with a positive obj coefficient
2241  * of the dominated variable y. we need to secure feasibility
2242  * by testing if the predicted lower worst case bound is less equal the
2243  * current upper bound. it is possible, that the lower worst case bound
2244  * is infinity and the upper bound of the dominating variable x is
2245  * infinity too.
2246  */
2247  varstofix[dominatedidx] = FIXATLB;
2248  (*nfixings)++;
2249  }
2250  }
2251 
2252  if( varstofix[dominatedidx] == NOFIX && !SCIPisInfinity(scip, dominatingub) &&
2253  SCIPisLE(scip, dominatingub, SCIPvarGetUbGlobal(dominatingvar)) )
2254  {
2255  /* we have a x->y dominance relation with an arbitrary obj coefficient
2256  * of the dominating variable x. in all cases we have to look
2257  * if the predicted upper bound of the dominating variable is great enough.
2258  * by testing, that the predicted upper bound is not infinity we avoid problems
2259  * with x->y e.g.
2260  * min -x -y
2261  * s.t. -x -y <= -1
2262  * 0<=x<=1, 0<=y<=1
2263  * where y is not at their lower bound.
2264  */
2265  varstofix[dominatedidx] = FIXATLB;
2266  (*nfixings)++;
2267  }
2268 
2269  if( varstofix[dominatingidx] == NOFIX && !SCIPisPositive(scip, SCIPvarGetObj(dominatingvar)) )
2270  {
2271  /* we have a x->y dominance relation with a negative obj coefficient
2272  * of the dominating variable x. if the worst case upper bound is
2273  * greater equal than upper bound, we fix x at the upper bound
2274  */
2275  if( !SCIPisInfinity(scip, dominatingwcub) &&
2276  SCIPisGE(scip, dominatingwcub, SCIPvarGetUbGlobal(dominatingvar)) )
2277  {
2278  varstofix[dominatingidx] = FIXATUB;
2279  (*nfixings)++;
2280  }
2281  }
2282 
2283  if( varstofix[dominatingidx] == NOFIX && !SCIPisInfinity(scip, -dominatinglb) &&
2284  SCIPisGE(scip, dominatinglb, SCIPvarGetUbGlobal(dominatingvar)) )
2285  {
2286  /* we have a x->y dominance relation with an arbitrary obj coefficient
2287  * of the dominating variable x. if the predicted lower bound is greater
2288  * equal than upper bound, we fix x at the upper bound.
2289  */
2290  varstofix[dominatingidx] = FIXATUB;
2291  (*nfixings)++;
2292  }
2293 
2294  if( onlybinvars )
2295  {
2296  if( varstofix[dominatedidx] == NOFIX && (onlyoneone || SCIPvarsHaveCommonClique(dominatingvar, TRUE, dominatedvar, TRUE, TRUE)) )
2297  {
2298  /* We have a (1->1)-clique with dominance relation (x->y) (x dominates y).
2299  * From this dominance relation, we know (1->0) is possible and not worse than (0->1)
2300  * concerning the objective function. It follows that only (1->0) or (0->0) are possible,
2301  * but in both cases y has the value 0 => y=0.
2302  */
2303  varstofix[dominatedidx] = FIXATLB;
2304  (*nfixings)++;
2305  }
2306 
2307  if( varstofix[dominatingidx] == NOFIX && SCIPvarsHaveCommonClique(dominatingvar, FALSE, dominatedvar, FALSE, TRUE) )
2308  {
2309  /* We have a (0->0)-clique with dominance relation x->y (x dominates y).
2310  * From this dominance relation, we know (1->0) is possible and not worse than (0->1)
2311  * concerning the objective function. It follows that only (1->0) or (1->1) are possible,
2312  * but in both cases x has the value 1 => x=1
2313  */
2314  varstofix[dominatingidx] = FIXATUB;
2315  (*nfixings)++;
2316  }
2317  }
2318  else
2319  assert(!onlyoneone);
2320 
2321  return SCIP_OKAY;
2322 }
2323 
2324 
2325 /** find dominance relation between variable pairs */
2326 static
2327 SCIP_RETCODE findDominancePairs(
2328  SCIP* scip, /**< SCIP main data structure */
2329  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
2330  SCIP_PRESOLDATA* presoldata, /**< presolver data */
2331  int* searchcols, /**< indexes of variables for pair comparisons */
2332  int searchsize, /**< number of variables for pair comparisons */
2333  SCIP_Bool onlybinvars, /**< flag indicating searchcols has only binary variable indexes */
2334  FIXINGDIRECTION* varstofix, /**< array holding information for later upper/lower bound fixing */
2335  int* nfixings, /**< found number of possible fixings */
2336  SCIP_Longint* ndomrelations, /**< found number of dominance relations */
2337  int* nchgbds /**< number of changed bounds */
2338  )
2339 {
2340  SCIP_Real* vals1;
2341  SCIP_Real* vals2;
2342  SCIP_Real tmpupperbounddominatingcol1;
2343  SCIP_Real tmpupperbounddominatingcol2;
2344  SCIP_Real tmpwclowerbounddominatingcol1;
2345  SCIP_Real tmpwclowerbounddominatingcol2;
2346  SCIP_Real tmplowerbounddominatingcol1;
2347  SCIP_Real tmplowerbounddominatingcol2;
2348  SCIP_Real tmpwcupperbounddominatingcol1;
2349  SCIP_Real tmpwcupperbounddominatingcol2;
2350  int* rows1;
2351  int* rows2;
2352  int nrows1;
2353  int nrows2;
2354  SCIP_Real tmpupperbounddominatedcol1;
2355  SCIP_Real tmpupperbounddominatedcol2;
2356  SCIP_Real tmpwclowerbounddominatedcol1;
2357  SCIP_Real tmpwclowerbounddominatedcol2;
2358  SCIP_Real tmplowerbounddominatedcol1;
2359  SCIP_Real tmplowerbounddominatedcol2;
2360  SCIP_Real tmpwcupperbounddominatedcol1;
2361  SCIP_Real tmpwcupperbounddominatedcol2;
2362  SCIP_Real obj1;
2363  SCIP_Bool col1domcol2;
2364  SCIP_Bool col2domcol1;
2365  SCIP_Bool onlyoneone;
2366  int cnt1;
2367  int cnt2;
2368  int col1;
2369  int col2;
2370  int r1;
2371  int r2;
2372  int paircnt;
2373  int oldnfixings;
2374 
2375  assert(scip != NULL);
2376  assert(matrix != NULL);
2377  assert(presoldata != NULL);
2378  assert(searchcols != NULL);
2379  assert(varstofix != NULL);
2380  assert(nfixings != NULL);
2381  assert(ndomrelations != NULL);
2382  assert(nchgbds != NULL);
2383 
2384  paircnt = 0;
2385  oldnfixings = *nfixings;
2386 
2387  /* pair comparisons */
2388  for( cnt1 = 0; cnt1 < searchsize; cnt1++ )
2389  {
2390  /* get index of the first variable */
2391  col1 = searchcols[cnt1];
2392 
2393  if( varstofix[col1] == FIXATLB )
2394  continue;
2395 
2396  obj1 = SCIPvarGetObj(matrix->vars[col1]);
2397 
2398  for( cnt2 = cnt1+1; cnt2 < searchsize; cnt2++ )
2399  {
2400  /* get index of the second variable */
2401  col2 = searchcols[cnt2];
2402 
2403  onlyoneone = FALSE;
2404 
2405  /* we always have minimize as obj sense */
2406 
2407  /* column 1 dominating column 2 */
2408  col1domcol2 = (obj1 <= SCIPvarGetObj(matrix->vars[col2]));
2409 
2410  /* column 2 dominating column 1 */
2411  col2domcol1 = (SCIPvarGetObj(matrix->vars[col2]) <= obj1);
2412 
2413  /* search only if nothing was found yet */
2414  col1domcol2 = col1domcol2 && (varstofix[col2] == NOFIX);
2415  col2domcol1 = col2domcol1 && (varstofix[col1] == NOFIX);
2416 
2417  /* we only search for a dominance relation if the lower bounds are not negative */
2418  if( !onlybinvars )
2419  {
2420  if( SCIPisLT(scip, SCIPvarGetLbGlobal(matrix->vars[col1]), 0.0) ||
2421  SCIPisLT(scip, SCIPvarGetLbGlobal(matrix->vars[col2]), 0.0) )
2422  {
2423  col1domcol2 = FALSE;
2424  col2domcol1 = FALSE;
2425  }
2426  }
2427 
2428  /* pair comparison control */
2429  if( paircnt == presoldata->numcurrentpairs )
2430  {
2431  assert(*nfixings >= oldnfixings);
2432  if( *nfixings == oldnfixings )
2433  {
2434  /* not enough fixings found, decrement number of comparisons */
2435  presoldata->numcurrentpairs >>= 1; /*lint !e702*/
2436  if( presoldata->numcurrentpairs < presoldata->numminpairs )
2437  presoldata->numcurrentpairs = presoldata->numminpairs;
2438 
2439  /* stop searching in this row */
2440  return SCIP_OKAY;
2441  }
2442  oldnfixings = *nfixings;
2443  paircnt = 0;
2444 
2445  /* increment number of comparisons */
2446  presoldata->numcurrentpairs <<= 1; /*lint !e701*/
2447  if( presoldata->numcurrentpairs > presoldata->nummaxpairs )
2448  presoldata->numcurrentpairs = presoldata->nummaxpairs;
2449  }
2450  paircnt++;
2451 
2452  if( !col1domcol2 && !col2domcol1 )
2453  continue;
2454 
2455  /* get the data for both columns */
2456  vals1 = matrix->colmatval + matrix->colmatbeg[col1];
2457  rows1 = matrix->colmatind + matrix->colmatbeg[col1];
2458  nrows1 = matrix->colmatcnt[col1];
2459  r1 = 0;
2460  vals2 = matrix->colmatval + matrix->colmatbeg[col2];
2461  rows2 = matrix->colmatind + matrix->colmatbeg[col2];
2462  nrows2 = matrix->colmatcnt[col2];
2463  r2 = 0;
2464 
2465  /* do we have a obj constant? */
2466  if( nrows1 == 0 || nrows2 == 0 )
2467  {
2468  col1domcol2 = FALSE;
2469  col2domcol1 = FALSE;
2470  continue;
2471  }
2472 
2473  /* initialize temporary bounds of dominating variable */
2474  tmpupperbounddominatingcol1 = SCIPinfinity(scip);
2475  tmpupperbounddominatingcol2 = tmpupperbounddominatingcol1;
2476  tmpwclowerbounddominatingcol1 = -SCIPinfinity(scip);
2477  tmpwclowerbounddominatingcol2 = tmpwclowerbounddominatingcol1;
2478  tmplowerbounddominatingcol1 = -SCIPinfinity(scip);
2479  tmplowerbounddominatingcol2 = tmplowerbounddominatingcol1;
2480  tmpwcupperbounddominatingcol1 = SCIPinfinity(scip);
2481  tmpwcupperbounddominatingcol2 = tmpwcupperbounddominatingcol1;
2482 
2483  /* initialize temporary bounds of dominated variable */
2484  tmpupperbounddominatedcol1 = SCIPinfinity(scip);
2485  tmpupperbounddominatedcol2 = tmpupperbounddominatedcol1;
2486  tmpwclowerbounddominatedcol1 = -SCIPinfinity(scip);
2487  tmpwclowerbounddominatedcol2 = tmpwclowerbounddominatedcol1;
2488  tmplowerbounddominatedcol1 = -SCIPinfinity(scip);
2489  tmplowerbounddominatedcol2 = tmplowerbounddominatedcol1;
2490  tmpwcupperbounddominatedcol1 = SCIPinfinity(scip);
2491  tmpwcupperbounddominatedcol2 = tmpwcupperbounddominatedcol1;
2492 
2493  /* compare rows of this column pair */
2494  while( (col1domcol2 || col2domcol1) && (r1 < nrows1 || r2 < nrows2))
2495  {
2496  assert((r1 >= nrows1-1) || (rows1[r1] < rows1[r1+1]));
2497  assert((r2 >= nrows2-1) || (rows2[r2] < rows2[r2+1]));
2498 
2499  /* there is a nonredundant row containing column 1 but not column 2 */
2500  if( r1 < nrows1 && (r2 == nrows2 || rows1[r1] < rows2[r2]) )
2501  {
2502  /* dominance depends on the type of relation */
2503  if( !matrix->isrhsinfinite[rows1[r1]] )
2504  {
2505  /* no dominance relation for equations or ranged rows */
2506  col2domcol1 = FALSE;
2507  col1domcol2 = FALSE;
2508  }
2509  else
2510  {
2511  /* >= relation, larger coefficients dominate smaller ones */
2512  if( vals1[r1] > 0.0 )
2513  col2domcol1 = FALSE;
2514  else
2515  col1domcol2 = FALSE;
2516  }
2517 
2518  r1++;
2519  }
2520  /* there is a nonredundant row containing column 2, but not column 1 */
2521  else if( r2 < nrows2 && (r1 == nrows1 || rows1[r1] > rows2[r2]) )
2522  {
2523  /* dominance depends on the type of relation */
2524  if( !matrix->isrhsinfinite[rows2[r2]] )
2525  {
2526  /* no dominance relation for equations or ranged rows */
2527  col2domcol1 = FALSE;
2528  col1domcol2 = FALSE;
2529  }
2530  else
2531  {
2532  /* >= relation, larger coefficients dominate smaller ones */
2533  if( vals2[r2] < 0.0 )
2534  col2domcol1 = FALSE;
2535  else
2536  col1domcol2 = FALSE;
2537  }
2538 
2539  r2++;
2540  }
2541  /* if both columns appear in a common row, compare the coefficients */
2542  else
2543  {
2544  assert(r1 < nrows1 && r2 < nrows2);
2545  assert(rows1[r1] == rows2[r2]);
2546 
2547  /* dominance depends on the type of inequality */
2548  if( !matrix->isrhsinfinite[rows1[r1]] )
2549  {
2550  /* no dominance relation if coefficients differ for equations or ranged rows */
2551  if( !SCIPisEQ(scip, vals1[r1], vals2[r2]) )
2552  {
2553  col2domcol1 = FALSE;
2554  col1domcol2 = FALSE;
2555  }
2556 
2557  if( onlybinvars )
2558  {
2559  if( !onlyoneone && (matrix->minactivityposinf[rows1[r1]] + matrix->minactivityneginf[rows1[r1]] == 0)
2560  && SCIPisFeasGE(scip, matrix->minactivity[rows1[r1]] + MIN(vals1[r1], vals2[r2]), matrix->rhs[rows1[r1]]) )
2561  {
2562  onlyoneone = TRUE;
2563  }
2564  }
2565  }
2566  else
2567  {
2568  /* >= relation, larger coefficients dominate smaller ones */
2569  if( vals1[r1] > vals2[r2] )
2570  col2domcol1 = FALSE;
2571  else if( vals1[r1] < vals2[r2] )
2572  col1domcol2 = FALSE;
2573  }
2574 
2575  /* we claim the same sign for the coefficients to achieve monotonically
2576  * decreasing predictive bound functions
2577  */
2578  if( !onlyoneone && ((vals1[r1] < 0 && vals2[r2] < 0) || (vals1[r1] > 0 && vals2[r2] > 0)) )
2579  {
2580  if( col1domcol2 )
2581  {
2582  /* update bounds of dominating variable for column 1 */
2583  SCIP_CALL( updateBounds(scip, matrix, rows1[r1],
2584  col1, vals1[r1], col2, vals2[r2], TRUE,
2585  &tmpupperbounddominatingcol1, &tmpwclowerbounddominatingcol1,
2586  &tmplowerbounddominatingcol1, &tmpwcupperbounddominatingcol1) );
2587 
2588  /* update bounds of dominated variable for column 1 */
2589  SCIP_CALL( updateBounds(scip, matrix, rows1[r1],
2590  col1, vals1[r1], col2, vals2[r2], FALSE,
2591  &tmpupperbounddominatedcol1, &tmpwclowerbounddominatedcol1,
2592  &tmplowerbounddominatedcol1, &tmpwcupperbounddominatedcol1) );
2593  }
2594 
2595  if( col2domcol1 )
2596  {
2597  /* update bounds of dominating variable for column 2 */
2598  SCIP_CALL( updateBounds(scip, matrix, rows2[r2],
2599  col2, vals2[r2], col1, vals1[r1], TRUE,
2600  &tmpupperbounddominatingcol2, &tmpwclowerbounddominatingcol2,
2601  &tmplowerbounddominatingcol2, &tmpwcupperbounddominatingcol2) );
2602 
2603  /* update bounds of dominated variable for column 2 */
2604  SCIP_CALL( updateBounds(scip, matrix, rows2[r2],
2605  col2, vals2[r2], col1, vals1[r1], FALSE,
2606  &tmpupperbounddominatedcol2, &tmpwclowerbounddominatedcol2,
2607  &tmplowerbounddominatedcol2, &tmpwcupperbounddominatedcol2) );
2608  }
2609  }
2610 
2611  r1++;
2612  r2++;
2613  }
2614  }
2615 
2616  /* a column is only dominated, if there are no more rows in which it is contained */
2617  col1domcol2 = col1domcol2 && r2 == nrows2;
2618  col2domcol1 = col2domcol1 && r1 == nrows1;
2619 
2620  if( !col1domcol2 && !col2domcol1 )
2621  continue;
2622 
2623  /* no dominance relation for left equations or ranged rows */
2624  while( r1 < nrows1 )
2625  {
2626  if( !matrix->isrhsinfinite[rows1[r1]] )
2627  {
2628  col2domcol1 = FALSE;
2629  col1domcol2 = FALSE;
2630  break;
2631  }
2632  r1++;
2633  }
2634  if( !col1domcol2 && !col2domcol1 )
2635  continue;
2636  while( r2 < nrows2 )
2637  {
2638  if( !matrix->isrhsinfinite[rows2[r2]] )
2639  {
2640  col2domcol1 = FALSE;
2641  col1domcol2 = FALSE;
2642  break;
2643  }
2644  r2++;
2645  }
2646 
2647  if( col1domcol2 || col2domcol1 )
2648  (*ndomrelations)++;
2649 
2650  if( col1domcol2 && col2domcol1 )
2651  {
2652  /* prefer the variable as dominating variable with the greater upper bound */
2653  if( SCIPisGE(scip, SCIPvarGetUbGlobal(matrix->vars[col1]), SCIPvarGetUbGlobal(matrix->vars[col2])) )
2654  col2domcol1 = FALSE;
2655  else
2656  col1domcol2 = FALSE;
2657  }
2658 
2659  /* use dominance relation and clique/bound-information
2660  * to find variable fixings. additionally try to strengthen
2661  * variable bounds by predictive bound strengthening.
2662  */
2663  if( col1domcol2 )
2664  {
2665  SCIP_CALL( findFixings(scip, matrix, matrix->vars[col1], col1,
2666  tmpupperbounddominatingcol1, tmpwclowerbounddominatingcol1,
2667  tmplowerbounddominatingcol1, tmpwcupperbounddominatingcol1,
2668  matrix->vars[col2], col2,
2669  varstofix, onlybinvars, onlyoneone, nfixings) );
2670 
2671  if( presoldata->predbndstr )
2672  {
2673  SCIP_CALL( predBndStr(scip, matrix->vars[col1], col1,
2674  tmpupperbounddominatingcol1,
2675  tmplowerbounddominatingcol1, tmpwcupperbounddominatingcol1,
2676  matrix->vars[col2], col2,
2677  tmpupperbounddominatedcol1, tmpwclowerbounddominatedcol1,
2678  tmplowerbounddominatedcol1,
2679  varstofix, nchgbds) );
2680  }
2681  }
2682  else if( col2domcol1 )
2683  {
2684  SCIP_CALL( findFixings(scip, matrix, matrix->vars[col2], col2,
2685  tmpupperbounddominatingcol2, tmpwclowerbounddominatingcol2,
2686  tmplowerbounddominatingcol2, tmpwcupperbounddominatingcol2,
2687  matrix->vars[col1], col1,
2688  varstofix, onlybinvars, onlyoneone, nfixings) );
2689 
2690  if( presoldata->predbndstr )
2691  {
2692  SCIP_CALL( predBndStr(scip, matrix->vars[col2], col2,
2693  tmpupperbounddominatingcol2,
2694  tmplowerbounddominatingcol2, tmpwcupperbounddominatingcol2,
2695  matrix->vars[col1], col1,
2696  tmpupperbounddominatedcol2, tmpwclowerbounddominatedcol2,
2697  tmplowerbounddominatedcol2,
2698  varstofix, nchgbds) );
2699  }
2700  }
2701  if( varstofix[col1] == FIXATLB )
2702  break;
2703  }
2704  }
2705 
2706  return SCIP_OKAY;
2707 }
2708 
2709 /** try to fix singleton column continuous variables */
2710 static
2711 SCIP_RETCODE singletonColumnStuffing(
2712  SCIP* scip, /**< SCIP main data structure */
2713  CONSTRAINTMATRIX* matrix, /**< matrix containing the constraints */
2714  SCIP_Bool* varsprocessed, /**< array indicating that this variable has been processed */
2715  FIXINGDIRECTION* varstofix, /**< array holding fixing information */
2716  int* nfixings /**< number of possible fixings */
2717  )
2718 {
2719  SCIP_VAR* var;
2720  SCIP_Real* valpnt;
2721  SCIP_Real* colratios;
2722  SCIP_Real* colcoeffs;
2723  SCIP_Bool* rowprocessed;
2724  int* rowpnt;
2725  int* rowend;
2726  int* colindices;
2727  int* colnozerolb;
2728  SCIP_Real constant1;
2729  SCIP_Real constant2;
2730  SCIP_Real val;
2731  SCIP_Real value;
2732  SCIP_Real boundoffset;
2733  SCIP_Bool tryfixing;
2734  int idx;
2735  int col;
2736  int row;
2737  int fillcnt;
2738  int colidx;
2739  int k;
2740 
2741  assert(scip != NULL);
2742  assert(matrix != NULL);
2743  assert(varsprocessed != NULL);
2744  assert(varstofix != NULL);
2745  assert(nfixings != NULL);
2746 
2747  SCIP_CALL( SCIPallocBufferArray(scip, &colindices, matrix->ncols) );
2748  SCIP_CALL( SCIPallocBufferArray(scip, &colratios, matrix->ncols) );
2749  SCIP_CALL( SCIPallocBufferArray(scip, &colcoeffs, matrix->ncols) );
2750 
2751  SCIP_CALL( SCIPallocBufferArray(scip, &colnozerolb, matrix->ncols) );
2752  BMSclearMemoryArray(colnozerolb, matrix->ncols);
2753 
2754  SCIP_CALL( SCIPallocBufferArray(scip, &rowprocessed, matrix->nrows) );
2755  BMSclearMemoryArray(rowprocessed, matrix->nrows);
2756 
2757  for( col = 0; col < matrix->ncols; col++ )
2758  {
2759  /* we look only at rows with minimal one continuous singleton column */
2760  if( matrix->colmatcnt[col] == 1 && SCIPvarGetType(matrix->vars[col]) == SCIP_VARTYPE_CONTINUOUS )
2761  {
2762  row = *(matrix->colmatind + matrix->colmatbeg[col]);
2763  if( rowprocessed[row] )
2764  continue;
2765 
2766  rowprocessed[row] = TRUE;
2767 
2768  if( matrix->isrhsinfinite[row] )
2769  {
2770  /* case for >= relation: coeff<0 and obj<0 */
2771  fillcnt = 0;
2772  tryfixing = TRUE;
2773  constant1 = 0.0;
2774  constant2 = 0.0;
2775 
2776  rowpnt = matrix->rowmatind + matrix->rowmatbeg[row];
2777  rowend = rowpnt + matrix->rowmatcnt[row];
2778  valpnt = matrix->rowmatval + matrix->rowmatbeg[row];
2779 
2780  for( ; (rowpnt < rowend); rowpnt++, valpnt++ )
2781  {
2782  val = *valpnt;
2783  colidx = *rowpnt;
2784  var = matrix->vars[colidx];
2785  assert(var != NULL);
2786  assert(val != 0.0);
2787 
2788  if( SCIPisGE(scip, SCIPvarGetLbGlobal(var), 0.0) )
2789  {
2790  if( SCIPisNegative(scip, val) )
2791  {
2792  /* do we have a continuous singleton column */
2793  if( matrix->colmatcnt[colidx] == 1 && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
2794  {
2795  if( SCIPisPositive(scip, SCIPvarGetLbGlobal(var)) )
2796  {
2797  constant1 += val * SCIPvarGetLbGlobal(var);
2798  constant2 += val * SCIPvarGetLbGlobal(var);
2799  colnozerolb[fillcnt] = 1;
2800  }
2801 
2802  if( SCIPisNegative(scip, SCIPvarGetObj(matrix->vars[colidx])) )
2803  {
2804  colratios[fillcnt] = SCIPvarGetObj(var) / val;
2805  colindices[fillcnt] = colidx;
2806  colcoeffs[fillcnt] = val;
2807  fillcnt++;
2808  }
2809  }
2810  else
2811  {
2812  /* discrete variables or variables in more than one row are present */
2813  if( SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) )
2814  {
2815  tryfixing = FALSE;
2816  break;
2817  }
2818  constant1 += val * SCIPvarGetUbGlobal(var);
2819  constant2 += val * SCIPvarGetLbGlobal(var);
2820  }
2821  }
2822  else if( SCIPisPositive(scip, val) )
2823  {
2824  constant1 += val * SCIPvarGetLbGlobal(var);
2825  constant2 += val * SCIPvarGetUbGlobal(var);
2826  }
2827  }
2828  else
2829  {
2830  tryfixing = FALSE;
2831  break;
2832  }
2833  }
2834 
2835  if( tryfixing )
2836  {
2837  SCIPsortRealRealIntInt(colratios, colcoeffs, colindices, colnozerolb, fillcnt);
2838 
2839  /* try to fix continuous singleton columns by their ratio */
2840  for( k = fillcnt - 1; k >= 0; k-- )
2841  {
2842  boundoffset = 0.0;
2843  idx = colindices[k];
2844  value = colcoeffs[k] * SCIPvarGetUbGlobal(matrix->vars[idx]);
2845 
2846  if( colnozerolb[k] )
2847  {
2848  boundoffset = colcoeffs[k] * SCIPvarGetLbGlobal(matrix->vars[idx]);
2849  }
2850 
2851  assert(SCIPisNegative(scip, SCIPvarGetObj(matrix->vars[idx])));
2852 
2853  if( matrix->colmatcnt[idx] == 1 && SCIPvarGetType(matrix->vars[idx]) == SCIP_VARTYPE_CONTINUOUS )
2854  {
2855  if( SCIPisGE(scip, value, matrix->lhs[row] - constant1 + boundoffset) )
2856  {
2857  varstofix[idx] = FIXATUB;
2858  varsprocessed[idx] = TRUE;
2859  (*nfixings)++;
2860  }
2861  else if( SCIPisGE(scip, matrix->lhs[row], constant2) )
2862  {
2863  varstofix[idx] = FIXATLB;
2864  varsprocessed[idx] = TRUE;
2865  (*nfixings)++;
2866  }
2867 
2868  constant1 += value;
2869  if( colnozerolb[k] )
2870  constant1 -= boundoffset;
2871 
2872  constant2 += value;
2873  if( colnozerolb[k] )
2874  constant2 -= boundoffset;
2875  }
2876  }
2877  }
2878  }
2879 
2880 
2881  if( matrix->isrhsinfinite[row] )
2882  {
2883  /* case for >= relation: coeff>0 and obj>0 */
2884  fillcnt = 0;
2885  tryfixing = TRUE;
2886  constant1 = 0.0;
2887  constant2 = 0.0;
2888 
2889  rowpnt = matrix->rowmatind + matrix->rowmatbeg[row];
2890  rowend = rowpnt + matrix->rowmatcnt[row];
2891  valpnt = matrix->rowmatval + matrix->rowmatbeg[row];
2892 
2893  for( ; (rowpnt < rowend); rowpnt++, valpnt++ )
2894  {
2895  val = *valpnt;
2896  colidx = *rowpnt;
2897  var = matrix->vars[colidx];
2898  assert(var != NULL);
2899  assert(val != 0.0);
2900 
2901  if( SCIPisGE(scip, SCIPvarGetLbGlobal(var), 0.0) )
2902  {
2903  if( SCIPisPositive(scip, val) )
2904  {
2905  /* do we have a continuous singleton column */
2906  if( matrix->colmatcnt[colidx] == 1 && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS &&
2907  SCIPisPositive(scip, SCIPvarGetObj(var)) )
2908  {
2909  if( SCIPisPositive(scip, SCIPvarGetLbGlobal(var)) )
2910  {
2911  constant1 += val * SCIPvarGetLbGlobal(var);
2912  constant2 += val * SCIPvarGetLbGlobal(var);
2913  colnozerolb[fillcnt] = 1;
2914  }
2915 
2916  colratios[fillcnt] = SCIPvarGetObj(var) / val;
2917  colindices[fillcnt] = colidx;
2918  colcoeffs[fillcnt] = val;
2919  fillcnt++;
2920  }
2921  else
2922  {
2923  /* discrete variables or variables which are present within
2924  * more than one row are estimated at their upper bound
2925  */
2926  if( SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) )
2927  {
2928  tryfixing = FALSE;
2929  break;
2930  }
2931  constant1 += val * SCIPvarGetUbGlobal(var);
2932  constant2 += val * SCIPvarGetLbGlobal(var);
2933  }
2934  }
2935  else if( SCIPisNegative(scip, val) )
2936  {
2937  /* consider lower bound for negative coefficients */
2938  constant1 += val * SCIPvarGetLbGlobal(var);
2939  constant2 += val * SCIPvarGetUbGlobal(var);
2940  }
2941  }
2942  else
2943  {
2944  tryfixing = FALSE;
2945  break;
2946  }
2947  }
2948 
2949  if( tryfixing )
2950  {
2951  SCIPsortRealRealIntInt(colratios, colcoeffs, colindices, colnozerolb, fillcnt);
2952 
2953  /* try to fix continuous singleton columns by their ratio */
2954  for( k = 0; k < fillcnt; k++ )
2955  {
2956  boundoffset = 0.0;
2957  idx = colindices[k];
2958  value = colcoeffs[k] * SCIPvarGetUbGlobal(matrix->vars[idx]);
2959 
2960  if( colnozerolb[k] )
2961  {
2962  boundoffset = colcoeffs[k] * SCIPvarGetLbGlobal(matrix->vars[idx]);
2963  }
2964 
2965  assert(SCIPisPositive(scip, SCIPvarGetObj(matrix->vars[idx])));
2966 
2967  if( matrix->colmatcnt[idx] == 1 && SCIPvarGetType(matrix->vars[idx]) == SCIP_VARTYPE_CONTINUOUS )
2968  {
2969  if( SCIPisLE(scip, value, matrix->lhs[row] - constant1 + boundoffset) )
2970  {
2971  varstofix[idx] = FIXATUB;
2972  varsprocessed[idx] = TRUE;
2973  (*nfixings)++;
2974  }
2975  else if( SCIPisLE(scip, matrix->lhs[row], constant2) )
2976  {
2977  varstofix[idx] = FIXATLB;
2978  varsprocessed[idx] = TRUE;
2979  (*nfixings)++;
2980  }
2981 
2982  constant1 += value;
2983  if( colnozerolb[k] )
2984  constant1 -= boundoffset;
2985 
2986  constant2 += value;
2987  if( colnozerolb[k] )
2988  constant2 -= boundoffset;
2989  }
2990  }
2991  }
2992  }
2993  }
2994  }
2995 
2996  SCIPfreeBufferArray(scip, &rowprocessed);
2997  SCIPfreeBufferArray(scip, &colnozerolb);
2998  SCIPfreeBufferArray(scip, &colcoeffs);
2999  SCIPfreeBufferArray(scip, &colratios);
3000  SCIPfreeBufferArray(scip, &colindices);
3001 
3002  return SCIP_OKAY;
3003 }
3004 
3005 /*
3006  * Callback methods of presolver
3007  */
3008 
3009 #ifdef SCIP_STATISTIC
3010 /** presolving deinitialization method of presolver (called after presolving has been finished) */
3011 static
3012 SCIP_DECL_PRESOLEXITPRE(presolExitpreDomcol)
3013 { /*lint --e{715}*/
3014  SCIP_PRESOLDATA* presoldata;
3015 
3016  /* free presolver data */
3017  presoldata = SCIPpresolGetData(presol);
3018  assert(presoldata != NULL);
3019 
3020  printf("### Stuffing: %d fixings ###\n",presoldata->nfixingsstuffing);
3021  printf("### Domcol: %d fixings ###\n", presoldata->nfixingsdomcol);
3022 
3023  presoldata->nfixingsstuffing = 0;
3024  presoldata->nfixingsdomcol = 0;
3025 
3026  return SCIP_OKAY;
3027 }
3028 #endif
3029 
3030 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
3031 static
3032 SCIP_DECL_PRESOLCOPY(presolCopyDomcol)
3033 { /*lint --e{715}*/
3034  assert(scip != NULL);
3035  assert(presol != NULL);
3036  assert(strcmp(SCIPpresolGetName(presol), PRESOL_NAME) == 0);
3037 
3038  /* call inclusion method of presolver */
3039  SCIP_CALL( SCIPincludePresolDomcol(scip) );
3040 
3041  return SCIP_OKAY;
3042 }
3043 
3044 /** destructor of presolver to free user data (called when SCIP is exiting) */
3045 static
3046 SCIP_DECL_PRESOLFREE(presolFreeDomcol)
3047 { /*lint --e{715}*/
3048  SCIP_PRESOLDATA* presoldata;
3049 
3050  /* free presolver data */
3051  presoldata = SCIPpresolGetData(presol);
3052  assert(presoldata != NULL);
3053 
3054  SCIPfreeMemory(scip, &presoldata);
3055  SCIPpresolSetData(presol, NULL);
3056 
3057  return SCIP_OKAY;
3058 }
3059 
3060 /** execution method of presolver */
3061 static
3062 SCIP_DECL_PRESOLEXEC(presolExecDomcol)
3063 { /*lint --e{715}*/
3064  SCIP_PRESOLDATA* presoldata;
3065  CONSTRAINTMATRIX* matrix;
3066  SCIP_Bool initialized;
3067 
3068  assert(result != NULL);
3069  *result = SCIP_DIDNOTRUN;
3070 
3071  /* do no dominated column presolving in case of probing and nonlinear processing
3072  * @todo SCIPisNLPEnabled() always returns FALSE during presolve, since the necessary flag is set after presolve (in exitpre, currently)
3073  */
3074  if( (SCIPgetStage(scip) != SCIP_STAGE_PRESOLVING) || SCIPinProbing(scip) || SCIPisNLPEnabled(scip) )
3075  return SCIP_OKAY;
3076 
3077  *result = SCIP_DIDNOTFIND;
3078 
3079  presoldata = SCIPpresolGetData(presol);
3080  assert(presoldata != NULL);
3081 
3082  /* initialize constraint matrix */
3083  matrix = NULL;
3084  initialized = FALSE;
3085  SCIP_CALL( initMatrix(scip, &matrix, &initialized) );
3086 
3087  if( initialized )
3088  {
3089  int nfixings;
3090  int nfixingsstuffing;
3091  SCIP_Longint ndomrelations;
3092  int v;
3093  int r;
3094  FIXINGDIRECTION* varstofix;
3095  SCIP_Bool* varsprocessed;
3096  int nvars;
3097  int nrows;
3098  int* rowidxsorted;
3099  int* rowsparsity;
3100  int varcount;
3101  int* consearchcols;
3102  int* intsearchcols;
3103  int* binsearchcols;
3104  int nconfill;
3105  int nintfill;
3106  int nbinfill;
3107 #ifdef SCIP_DEBUG
3108  int nconvarsfixed = 0;
3109  int nintvarsfixed = 0;
3110  int nbinvarsfixed = 0;
3111 #endif
3112  int* pclass;
3113  int* colidx;
3114  int pclassstart;
3115  int pc;
3116  SCIP_Bool* varineq;
3117 
3118  assert(SCIPgetNVars(scip) == matrix->ncols);
3119 
3120  nfixings = 0;
3121  nfixingsstuffing = 0;
3122  ndomrelations = 0;
3123  nvars = matrix->ncols;
3124  nrows = matrix->nrows;
3125 
3126  SCIP_CALL( SCIPallocBufferArray(scip, &varstofix, nvars) );
3127  SCIP_CALL( SCIPallocBufferArray(scip, &varsprocessed, nvars) );
3128  BMSclearMemoryArray(varstofix, nvars);
3129  BMSclearMemoryArray(varsprocessed, nvars);
3130 
3131  SCIP_CALL( SCIPallocBufferArray(scip, &consearchcols, nvars) );
3132  SCIP_CALL( SCIPallocBufferArray(scip, &intsearchcols, nvars) );
3133  SCIP_CALL( SCIPallocBufferArray(scip, &binsearchcols, nvars) );
3134 
3135  SCIP_CALL( SCIPallocBufferArray(scip, &rowidxsorted, nrows) );
3136  SCIP_CALL( SCIPallocBufferArray(scip, &rowsparsity, nrows) );
3137  for( r = 0; r < nrows; ++r )
3138  {
3139  rowidxsorted[r] = r;
3140  rowsparsity[r] = matrix->rowmatcnt[r];
3141  }
3142 
3143  SCIP_CALL( SCIPallocBufferArray(scip, &pclass, nvars) );
3144  SCIP_CALL( SCIPallocBufferArray(scip, &colidx, nvars) );
3145  SCIP_CALL( SCIPallocBufferArray(scip, &varineq, nvars) );
3146  for( v = 0; v < nvars; v++ )
3147  {
3148  colidx[v] = v;
3149  varineq[v] = FALSE;
3150  }
3151 
3152  /* init pair comparision control */
3153  presoldata->numcurrentpairs = presoldata->nummaxpairs;
3154 
3155 
3156  /* before doing dominated column presolving we stuff singleton continuous columns.
3157  * this sometimes helps to do a more effective predictive bound analysis.
3158  */
3159  if( presoldata->singcolstuffing )
3160  {
3161  SCIP_CALL( singletonColumnStuffing(scip, matrix, varsprocessed, varstofix, &nfixingsstuffing) );
3162 #ifdef SCIP_STATISTIC
3163  presoldata->nfixingsstuffing += nfixingsstuffing;
3164 #endif
3165  }
3166 
3167 
3168  /* 1.stage: we search for dominance relations only within parallel columns
3169  * concerning equalities and ranged rows
3170  */
3171 
3172  SCIP_CALL( detectParallelCols(scip, matrix, pclass, varineq) );
3173  SCIPsortIntInt(pclass, colidx, nvars);
3174 
3175  varcount = 0;
3176 
3177  pc = 0;
3178  while( pc < nvars )
3179  {
3180  int varidx;
3181 
3182  varidx = 0;
3183  nconfill = 0;
3184  nintfill = 0;
3185  nbinfill = 0;
3186 
3187  pclassstart = pclass[pc];
3188  while( pc < nvars && pclassstart == pclass[pc] )
3189  {
3190  varidx = colidx[pc];
3191 
3192  /* we only regard variables which were not processed yet and are present within equalities or ranged rows */
3193  if( !varsprocessed[varidx] && varineq[varidx] )
3194  {
3195  /* we search only for dominance relations between the same variable type */
3196  if( SCIPvarGetType(matrix->vars[varidx]) == SCIP_VARTYPE_CONTINUOUS )
3197  {
3198  consearchcols[nconfill++] = varidx;
3199  }
3200  else if( SCIPvarIsBinary(matrix->vars[varidx]) )
3201  {
3202  binsearchcols[nbinfill++] = varidx;
3203  }
3204  else
3205  {
3206  assert(SCIPvarGetType(matrix->vars[varidx]) == SCIP_VARTYPE_INTEGER || SCIPvarGetType(matrix->vars[varidx]) == SCIP_VARTYPE_IMPLINT);
3207  intsearchcols[nintfill++] = varidx;
3208  }
3209  }
3210  ++pc;
3211  }
3212 
3213  /* search for dominance relations between continuous variables */
3214  if( nconfill > 1 )
3215  {
3216  SCIP_CALL( findDominancePairs(scip, matrix, presoldata, consearchcols, nconfill, FALSE,
3217  varstofix, &nfixings, &ndomrelations, nchgbds) );
3218 
3219  for( v = 0; v < nconfill; ++v )
3220  varsprocessed[consearchcols[v]] = TRUE;
3221 
3222  varcount += nconfill;
3223  }
3224  else if( nconfill == 1 )
3225  {
3226  if( varineq[varidx] )
3227  varsprocessed[consearchcols[0]] = TRUE;
3228  }
3229 
3230  /* search for dominance relations between integer and impl-integer variables */
3231  if( nintfill > 1 )
3232  {
3233  SCIP_CALL( findDominancePairs(scip, matrix, presoldata, intsearchcols, nintfill, FALSE,
3234  varstofix, &nfixings, &ndomrelations, nchgbds) );
3235 
3236  for( v = 0; v < nintfill; ++v )
3237  varsprocessed[intsearchcols[v]] = TRUE;
3238 
3239  varcount += nintfill;
3240  }
3241  else if( nintfill == 1 )
3242  {
3243  if( varineq[varidx] )
3244  varsprocessed[intsearchcols[0]] = TRUE;
3245  }
3246 
3247  /* search for dominance relations between binary variables */
3248  if( nbinfill > 1 )
3249  {
3250  SCIP_CALL( findDominancePairs(scip, matrix, presoldata, binsearchcols, nbinfill, TRUE,
3251  varstofix, &nfixings, &ndomrelations, nchgbds) );
3252 
3253  for( v = 0; v < nbinfill; ++v )
3254  varsprocessed[binsearchcols[v]] = TRUE;
3255 
3256  varcount += nbinfill;
3257  }
3258  else if( nbinfill == 1 )
3259  {
3260  if( varineq[varidx] )
3261  varsprocessed[binsearchcols[0]] = TRUE;
3262  }
3263 
3264  /* break if no vars are left */
3265  if( varcount >= nvars )
3266  break;
3267  }
3268 
3269 
3270  /* 2.stage: we search for dominance relations of the left columns
3271  * by row-sparsity order
3272  */
3273 
3274  /* sort rows per sparsity monotonically increasing */
3275  SCIPsortIntInt(rowsparsity, rowidxsorted, nrows);
3276 
3277  for( r = 0; r < nrows; ++r )
3278  {
3279  int rowidx;
3280  int* rowpnt;
3281  int* rowend;
3282 
3283  /* break if the time limit was reached; since the check is expensive,
3284  * we only check all 1000 constraints
3285  */
3286  if( (r % 1000 == 0) && SCIPisStopped(scip) )
3287  break;
3288 
3289  rowidx = rowidxsorted[r];
3290  rowpnt = matrix->rowmatind + matrix->rowmatbeg[rowidx];
3291  rowend = rowpnt + matrix->rowmatcnt[rowidx];
3292 
3293  if( matrix->rowmatcnt[rowidx] == 1 )
3294  continue;
3295 
3296  nconfill = 0;
3297  nintfill = 0;
3298  nbinfill = 0;
3299 
3300  for( ; rowpnt < rowend; rowpnt++ )
3301  {
3302  if( !(varsprocessed[*rowpnt]) )
3303  {
3304  int varidx;
3305  varidx = *rowpnt;
3306 
3307  /* we search only for dominance relations between the same variable type */
3308  if( SCIPvarGetType(matrix->vars[varidx]) == SCIP_VARTYPE_CONTINUOUS )
3309  {
3310  consearchcols[nconfill++] = varidx;
3311  }
3312  else if( SCIPvarIsBinary(matrix->vars[varidx]) )
3313  {
3314  binsearchcols[nbinfill++] = varidx;
3315  }
3316  else
3317  {
3318  assert(SCIPvarGetType(matrix->vars[varidx]) == SCIP_VARTYPE_INTEGER || SCIPvarGetType(matrix->vars[varidx]) == SCIP_VARTYPE_IMPLINT);
3319  intsearchcols[nintfill++] = varidx;
3320  }
3321  }
3322  }
3323 
3324  /* search for dominance relations between continuous variables */
3325  if( nconfill > 1 )
3326  {
3327  SCIP_CALL( findDominancePairs(scip, matrix, presoldata, consearchcols, nconfill, FALSE,
3328  varstofix, &nfixings, &ndomrelations, nchgbds) );
3329 
3330  for( v = 0; v < nconfill; ++v )
3331  varsprocessed[consearchcols[v]] = TRUE;
3332 
3333  varcount += nconfill;
3334  }
3335 
3336  /* search for dominance relations between integer and impl-integer variables */
3337  if( nintfill > 1 )
3338  {
3339  SCIP_CALL( findDominancePairs(scip, matrix, presoldata, intsearchcols, nintfill, FALSE,
3340  varstofix, &nfixings, &ndomrelations, nchgbds) );
3341 
3342  for( v = 0; v < nintfill; ++v )
3343  varsprocessed[intsearchcols[v]] = TRUE;
3344 
3345  varcount += nintfill;
3346  }
3347 
3348  /* search for dominance relations between binary variables */
3349  if( nbinfill > 1 )
3350  {
3351  SCIP_CALL( findDominancePairs(scip, matrix, presoldata, binsearchcols, nbinfill, TRUE,
3352  varstofix, &nfixings, &ndomrelations, nchgbds) );
3353 
3354  for( v = 0; v < nbinfill; ++v )
3355  varsprocessed[binsearchcols[v]] = TRUE;
3356 
3357  varcount += nbinfill;
3358  }
3359 
3360  /* break if no vars are left */
3361  if( varcount >= nvars )
3362  break;
3363  }
3364 
3365 
3366 #ifdef SCIP_STATISTIC
3367  presoldata->nfixingsdomcol += nfixings;
3368 #endif
3369 
3370  if( (nfixings + nfixingsstuffing) > 0 )
3371  {
3372  int oldnfixedvars = *nfixedvars;
3373 
3374  /* look for fixable variables */
3375  for( v = matrix->ncols - 1; v >= 0; --v )
3376  {
3377  SCIP_Bool infeasible;
3378  SCIP_Bool fixed;
3379  SCIP_VAR* var;
3380 
3381  if( SCIPvarGetNLocksUp(matrix->vars[v]) != matrix->nuplocks[v] ||
3382  SCIPvarGetNLocksDown(matrix->vars[v]) != matrix->ndownlocks[v] )
3383  {
3384  /* no fixing, maybe countsols constraint handler did additional locks */
3385  continue;
3386  }
3387 
3388  if( varstofix[v] == FIXATLB )
3389  {
3390  SCIP_Real lb;
3391 
3392  var = matrix->vars[v];
3393  lb = SCIPvarGetLbGlobal(var);
3394  assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || SCIPisFeasIntegral(scip, lb));
3395 
3396  SCIPdebugMessage("fixing dominated variable <%s> to its lower bound %g\n", SCIPvarGetName(var), lb);
3397 
3398  /* fix at lower bound */
3399  SCIP_CALL( SCIPfixVar(scip, var, lb, &infeasible, &fixed) );
3400  if( infeasible )
3401  {
3402  SCIPdebugMessage(" -> infeasible fixing\n");
3403  *result = SCIP_CUTOFF;
3404 
3405  break;
3406  }
3407  assert(fixed);
3408  (*nfixedvars)++;
3409 
3410 #ifdef SCIP_DEBUG
3411  if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
3412  {
3413  nconvarsfixed++;
3414  }
3415  else if( SCIPvarIsBinary(var) )
3416  {
3417  nbinvarsfixed++;
3418  }
3419  else
3420  {
3421  assert(SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER || SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT);
3422  nintvarsfixed++;
3423  }
3424 #endif
3425  }
3426  else if( varstofix[v] == FIXATUB )
3427  {
3428  SCIP_Real ub;
3429 
3430  var = matrix->vars[v];
3431  ub = SCIPvarGetUbGlobal(var);
3432  assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || SCIPisFeasIntegral(scip, ub));
3433 
3434  SCIPdebugMessage("fixing dominated variable <%s> to its upper bound %g\n", SCIPvarGetName(var), ub);
3435 
3436  /* fix at upper bound */
3437  SCIP_CALL( SCIPfixVar(scip, var, ub, &infeasible, &fixed) );
3438  if( infeasible )
3439  {
3440  SCIPdebugMessage(" -> infeasible fixing\n");
3441  *result = SCIP_CUTOFF;
3442 
3443  break;
3444  }
3445  assert(fixed);
3446 
3447  (*nfixedvars)++;
3448 
3449 #ifdef SCIP_DEBUG
3450  if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
3451  {
3452  nconvarsfixed++;
3453  }
3454  else if( SCIPvarIsBinary(var) )
3455  {
3456  nbinvarsfixed++;
3457  }
3458  else
3459  {
3460  assert(SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER || SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT);
3461  nintvarsfixed++;
3462  }
3463 #endif
3464  }
3465  }
3466 
3467  if( *result != SCIP_CUTOFF && *nfixedvars > oldnfixedvars )
3468  *result = SCIP_SUCCESS;
3469  }
3470 
3471  SCIPfreeBufferArray(scip, &varineq);
3472  SCIPfreeBufferArray(scip, &colidx);
3473  SCIPfreeBufferArray(scip, &pclass);
3474  SCIPfreeBufferArray(scip, &rowsparsity);
3475  SCIPfreeBufferArray(scip, &rowidxsorted);
3476  SCIPfreeBufferArray(scip, &binsearchcols);
3477  SCIPfreeBufferArray(scip, &intsearchcols);
3478  SCIPfreeBufferArray(scip, &consearchcols);
3479  SCIPfreeBufferArray(scip, &varsprocessed);
3480  SCIPfreeBufferArray(scip, &varstofix);
3481 
3482 #ifdef SCIP_DEBUG
3483  if( (nconvarsfixed + nintvarsfixed + nbinvarsfixed) > 0 )
3484  {
3485  SCIPdebugMessage("### %d vars [%lld dom] ===>>> fixed [cont: %d, int: %d, bin: %d], %scutoff detected\n",
3486  matrix->ncols, ndomrelations, nconvarsfixed, nintvarsfixed, nbinvarsfixed, (*result != SCIP_CUTOFF) ? "no " : "");
3487  }
3488 #endif
3489  }
3490 
3491  freeMatrix(scip, &matrix);
3492 
3493  return SCIP_OKAY;
3494 }
3495 
3496 /*
3497  * presolver specific interface methods
3498  */
3499 
3500 /** creates the domcol presolver and includes it in SCIP */
3501 SCIP_RETCODE SCIPincludePresolDomcol(
3502  SCIP* scip /**< SCIP data structure */
3503  )
3504 {
3505  SCIP_PRESOLDATA* presoldata;
3506  SCIP_PRESOL* presol;
3507 
3508  /* create domcol presolver data */
3509  SCIP_CALL( SCIPallocMemory(scip, &presoldata) );
3510 
3511  /* include presolver */
3512  SCIP_CALL( SCIPincludePresolBasic(scip, &presol, PRESOL_NAME, PRESOL_DESC, PRESOL_PRIORITY, PRESOL_MAXROUNDS,
3513  PRESOL_DELAY, presolExecDomcol, presoldata) );
3514  SCIP_CALL( SCIPsetPresolCopy(scip, presol, presolCopyDomcol) );
3515  SCIP_CALL( SCIPsetPresolFree(scip, presol, presolFreeDomcol) );
3516 
3517 #ifdef SCIP_STATISTIC
3518  presoldata->nfixingsstuffing = 0;
3519  presoldata->nfixingsdomcol = 0;
3520  SCIP_CALL( SCIPsetPresolExitpre(scip, presol, presolExitpreDomcol) );
3521 #endif
3522 
3523  SCIP_CALL( SCIPaddIntParam(scip,
3524  "presolving/domcol/numminpairs",
3525  "minimal number of pair comparisons",
3526  &presoldata->numminpairs, FALSE, DEFAULT_NUMMINPAIRS, 100, DEFAULT_NUMMAXPAIRS, NULL, NULL) );
3527 
3528  SCIP_CALL( SCIPaddIntParam(scip,
3529  "presolving/domcol/nummaxpairs",
3530  "maximal number of pair comparisons",
3531  &presoldata->nummaxpairs, FALSE, DEFAULT_NUMMAXPAIRS, DEFAULT_NUMMINPAIRS, 1000000000, NULL, NULL) );
3532 
3533  SCIP_CALL( SCIPaddBoolParam(scip,
3534  "presolving/domcol/predbndstr",
3535  "should predictive bound strengthening be applied?",
3536  &presoldata->predbndstr, FALSE, DEFAULT_PREDBNDSTR, NULL, NULL) );
3537 
3538  SCIP_CALL( SCIPaddBoolParam(scip,
3539  "presolving/domcol/singcolstuffing",
3540  "should singleton columns stuffing be applied?",
3541  &presoldata->singcolstuffing, FALSE, DEFAULT_SINGCOLSTUFF, NULL, NULL) );
3542 
3543  return SCIP_OKAY;
3544 }
3545