1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   presol_dualsparsify.c
26   	 * @brief  cancel nonzeros of the constraint matrix based on the columns
27   	 * @author Dieter Weninger
28   	 * @author Leona Gottwald
29   	 * @author Ambros Gleixner
30   	 * @author Weikun Chen
31   	 *
32   	 * This presolver attempts to cancel non-zero entries of the constraint
33   	 * matrix by adding scaled columns to other columns.
34   	 *
35   	 * In more detail, for two columns A_{j.} and A_{k.}, suppose for a given value s, we have
36   	 *
37   	 *                  | A_{j.} | - | A_{j.} - s*A_{k.} | > eta,
38   	 *
39   	 * where eta is an nonnegative integer. Then we introduce a new variable y := s*x_j + x_k
40   	 * and aggregate the variable x_k = y - s*x_j. After aggregation, the column of the variable
41   	 * x_j is A_{j.} + s*A_{j.} which is sparser than A_{j.}. In the case that x_k is nonimplied
42   	 * free variable, we need to add a new constraint l_k <= y - weight*x_j <= u_k into the problem
43   	 * to keep the bounds constraints of variable x_k.
44   	 *
45   	 * Further information can be found in
46   	 * Chen et al. "Two-row and two-column mixed-integer presolve using hasing-based pairing methods".
47   	 *
48   	 * @todo add infrastructure to SCIP for handling aggregated binary variables
49   	 */
50   	
51   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
52   	
53   	#include "blockmemshell/memory.h"
54   	#include "scip/cons_linear.h"
55   	#include "scip/debug.h"
56   	#include "scip/presol_dualsparsify.h"
57   	#include "scip/pub_cons.h"
58   	#include "scip/pub_matrix.h"
59   	#include "scip/pub_message.h"
60   	#include "scip/pub_misc.h"
61   	#include "scip/pub_misc_sort.h"
62   	#include "scip/pub_presol.h"
63   	#include "scip/pub_var.h"
64   	#include "scip/scip_cons.h"
65   	#include "scip/scip_general.h"
66   	#include "scip/scip_mem.h"
67   	#include "scip/scip_message.h"
68   	#include "scip/scip_nlp.h"
69   	#include "scip/scip_numerics.h"
70   	#include "scip/scip_param.h"
71   	#include "scip/scip_presol.h"
72   	#include "scip/scip_pricer.h"
73   	#include "scip/scip_prob.h"
74   	#include "scip/scip_probing.h"
75   	#include "scip/scip_solvingstats.h"
76   	#include "scip/scip_timing.h"
77   	#include "scip/scip_var.h"
78   	#include <string.h>
79   	
80   	#define PRESOL_NAME            "dualsparsify"
81   	#define PRESOL_DESC            "eliminate non-zero coefficients"
82   	
83   	#define PRESOL_PRIORITY           -240000    /**< priority of the presolver (>= 0: before, < 0: after constraint handlers) */
84   	#define PRESOL_MAXROUNDS               -1    /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */
85   	#define PRESOL_TIMING           SCIP_PRESOLTIMING_EXHAUSTIVE /* timing of the presolver (fast, medium, or exhaustive) */
86   	
87   	#define DEFAULT_ENABLECOPY           TRUE    /**< should dualsparsify presolver be copied to sub-SCIPs? */
88   	#define DEFAULT_PRESERVEINTCOEFS    FALSE    /**< should we forbid cancellations that destroy integer coefficients? */
89   	#define DEFAULT_PRESERVEGOODLOCKS   FALSE    /**< should we preserve good locked properties of variables (at most one lock in one direction)? */
90   	#define DEFAULT_MAX_CONT_FILLIN         1    /**< default value for the maximal fillin for continuous variables */
91   	#define DEFAULT_MAX_BIN_FILLIN          1    /**< default value for the maximal fillin for binary variables */
92   	#define DEFAULT_MAX_INT_FILLIN          1    /**< default value for the maximal fillin for integer variables (including binary) */
93   	#define DEFAULT_MAXCONSIDEREDNONZEROS  70    /**< maximal number of considered nonzeros within one column (-1: no limit) */
94   	#define DEFAULT_MINELIMINATEDNONZEROS 100    /**< minimal eleminated nonzeros within one column if we need to add a constraint to the problem */
95   	#define DEFAULT_MAXRETRIEVEFAC      100.0    /**< limit on the number of useless vs. useful hashtable retrieves as a multiple of the number of constraints */
96   	#define DEFAULT_WAITINGFAC            2.0    /**< number of calls to wait until next execution as a multiple of the number of useless calls */
97   	
98   	#define MAXSCALE                   1000.0    /**< maximal allowed scale for cancelling nonzeros */
99   	
100  	
101  	/*
102  	 * Data structures
103  	 */
104  	
105  	/** presolver data */
106  	struct SCIP_PresolData
107  	{
108  	   int                   nfailures;          /**< number of calls to presolver without success */
109  	   int                   nwaitingcalls;      /**< number of presolver calls until next real execution */
110  	   int                   naggregated;        /**< number of aggregated variables */
111  	   int                   maxcontfillin;      /**< maximal fillin for continuous variables */
112  	   int                   maxintfillin;       /**< maximal fillin for integer variables*/
113  	   int                   maxbinfillin;       /**< maximal fillin for binary variables */
114  	   int                   maxconsiderednonzeros;/**< maximal number of considered nonzeros within one column (-1: no limit) */
115  	   int                   mineliminatednonzeros;/**< minimal eliminated nonzeros within one column if we need to add a constraint to the problem */
116  	   SCIP_Real             maxretrievefac;     /**< limit on the number of useless vs. useful hashtable retrieves as a multiple of the number of constraints */
117  	   SCIP_Real             waitingfac;         /**< number of calls to wait until next execution as a multiple of the number of useless calls */
118  	   SCIP_Bool             enablecopy;         /**< should dualsparsify presolver be copied to sub-SCIPs? */
119  	   SCIP_Bool             preserveintcoefs;   /**< should we forbid cancellations that destroy integer coefficients? */
120  	   SCIP_Bool             preservegoodlocks;  /**< should we preserve good locked properties of variables (at most one lock in one direction)? */
121  	};
122  	
123  	/** structure representing a pair of constraints in a column; used for lookup in a hashtable */
124  	struct ColConsPair
125  	{
126  	   int colindex;                             /**< index of the column */
127  	   int consindex1;                           /**< index of the first constraint */
128  	   int consindex2;                           /**< index of the second constraint */
129  	   SCIP_Real conscoef1;                      /**< coefficient of the first constraint */
130  	   SCIP_Real conscoef2;                      /**< coefficient of the second constriant */
131  	};
132  	
133  	typedef struct ColConsPair COLCONSPAIR;
134  	
135  	/*
136  	 * Local methods
137  	 */
138  	
139  	/** returns TRUE iff both keys are equal */
140  	static
141  	SCIP_DECL_HASHKEYEQ(consPairsEqual)
142  	{  /*lint --e{715}*/
143  	   SCIP* scip;
144  	   COLCONSPAIR* conspair1;
145  	   COLCONSPAIR* conspair2;
146  	   SCIP_Real ratio1;
147  	   SCIP_Real ratio2;
148  	
149  	   scip = (SCIP*) userptr;
150  	   conspair1 = (COLCONSPAIR*) key1;
151  	   conspair2 = (COLCONSPAIR*) key2;
152  	
153  	   if( conspair1->consindex1 != conspair2->consindex1 )
154  	      return FALSE;
155  	
156  	   if( conspair1->consindex2 != conspair2->consindex2 )
157  	      return FALSE;
158  	
159  	   ratio1 = conspair1->conscoef2 / conspair1->conscoef1;
160  	   ratio2 = conspair2->conscoef2 / conspair2->conscoef1;
161  	
162  	   return SCIPisEQ(scip, ratio1, ratio2);
163  	}
164  	
165  	/** returns the hash value of the key */
166  	static
167  	SCIP_DECL_HASHKEYVAL(consPairHashval)
168  	{  /*lint --e{715}*/
169  	   COLCONSPAIR* conspair;
170  	
171  	   conspair = (COLCONSPAIR*) key;
172  	
173  	   return SCIPhashThree(conspair->consindex1, conspair->consindex2, SCIPrealHashCode(conspair->conscoef2 / conspair->conscoef1));
174  	}
175  	
176  	/** calculate maximal activity of one row without one specific column */
177  	static
178  	SCIP_Real getMaxActivitySingleRowWithoutCol(
179  	   SCIP*                 scip,               /**< SCIP main data structure */
180  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
181  	   int                   row,                /**< row index */
182  	   int                   col                 /**< column index */
183  	   )
184  	{
185  	   SCIP_Real* valpnt;
186  	   int* rowpnt;
187  	   int* rowend;
188  	   SCIP_Real maxactivity;
189  	   SCIP_Real val;
190  	   SCIP_Real lb;
191  	   SCIP_Real ub;
192  	   int c;
193  	
194  	   assert(scip != NULL);
195  	   assert(matrix != NULL);
196  	
197  	   maxactivity = 0;
198  	
199  	   rowpnt = SCIPmatrixGetRowIdxPtr(matrix, row);
200  	   rowend = rowpnt + SCIPmatrixGetRowNNonzs(matrix, row);
201  	   valpnt = SCIPmatrixGetRowValPtr(matrix, row);
202  	
203  	   for( ; (rowpnt < rowend); rowpnt++, valpnt++ )
204  	   {
205  	      c = *rowpnt;
206  	      val = *valpnt;
207  	
208  	      if( c == col )
209  	         continue;
210  	
211  	      if( val > 0.0 )
212  	      {
213  	         ub = SCIPmatrixGetColUb(matrix, c);
214  	         assert(!SCIPisInfinity(scip, ub));
215  	
216  	         maxactivity += val * ub;
217  	      }
218  	      else if( val < 0.0 )
219  	      {
220  	         lb = SCIPmatrixGetColLb(matrix, c);
221  	         assert(!SCIPisInfinity(scip, -lb));
222  	
223  	         maxactivity += val * lb;
224  	      }
225  	   }
226  	
227  	   return maxactivity;
228  	}
229  	
230  	/** calculate minimal activity of one row without one specific column */
231  	static
232  	SCIP_Real getMinActivitySingleRowWithoutCol(
233  	   SCIP*                 scip,               /**< SCIP main data structure */
234  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
235  	   int                   row,                /**< row index */
236  	   int                   col                 /**< column index */
237  	   )
238  	{
239  	   SCIP_Real* valpnt;
240  	   int* rowpnt;
241  	   int* rowend;
242  	   SCIP_Real minactivity;
243  	   SCIP_Real val;
244  	   SCIP_Real lb;
245  	   SCIP_Real ub;
246  	   int c;
247  	
248  	   assert(scip != NULL);
249  	   assert(matrix != NULL);
250  	
251  	   minactivity = 0;
252  	
253  	   rowpnt = SCIPmatrixGetRowIdxPtr(matrix, row);
254  	   rowend = rowpnt + SCIPmatrixGetRowNNonzs(matrix, row);
255  	   valpnt = SCIPmatrixGetRowValPtr(matrix, row);
256  	
257  	   for( ; (rowpnt < rowend); rowpnt++, valpnt++ )
258  	   {
259  	      c = *rowpnt;
260  	      val = *valpnt;
261  	
262  	      if( c == col )
263  	         continue;
264  	
265  	      if( val > 0.0 )
266  	      {
267  	         lb = SCIPmatrixGetColLb(matrix, c);
268  	         assert(!SCIPisInfinity(scip, -lb));
269  	
270  	         minactivity += val * lb;
271  	      }
272  	      else if( val < 0.0 )
273  	      {
274  	         ub = SCIPmatrixGetColUb(matrix, c);
275  	         assert(!SCIPisInfinity(scip, ub));
276  	
277  	         minactivity += val * ub;
278  	      }
279  	   }
280  	
281  	   return minactivity;
282  	}
283  	
284  	/** get minimal and maximal residual activity without one specified column */
285  	static
286  	void getMinMaxActivityResiduals(
287  	   SCIP*                 scip,               /**< SCIP main data structure */
288  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
289  	   int                   col,                /**< column index */
290  	   int                   row,                /**< row index */
291  	   SCIP_Real             val,                /**< coefficient of this variable in this row */
292  	   SCIP_Real*            minresactivity,     /**< minimum residual activity of this row */
293  	   SCIP_Real*            maxresactivity,     /**< maximum residual activity of this row */
294  	   SCIP_Bool*            isminsettoinfinity, /**< flag indicating if minresactiviy is set to infinity */
295  	   SCIP_Bool*            ismaxsettoinfinity  /**< flag indicating if maxresactiviy is set to infinity */
296  	   )
297  	{
298  	   SCIP_Real lb;
299  	   SCIP_Real ub;
300  	   int nmaxactneginf;
301  	   int nmaxactposinf;
302  	   int nminactneginf;
303  	   int nminactposinf;
304  	   SCIP_Real maxactivity;
305  	   SCIP_Real minactivity;
306  	
307  	   assert(scip != NULL);
308  	   assert(matrix != NULL);
309  	   assert(minresactivity != NULL);
310  	   assert(maxresactivity != NULL);
311  	   assert(isminsettoinfinity != NULL);
312  	   assert(ismaxsettoinfinity != NULL);
313  	
314  	   lb = SCIPmatrixGetColLb(matrix, col);
315  	   ub = SCIPmatrixGetColUb(matrix, col);
316  	
317  	   *isminsettoinfinity = FALSE;
318  	   *ismaxsettoinfinity = FALSE;
319  	
320  	   nmaxactneginf = SCIPmatrixGetRowNMaxActNegInf(matrix, row);
321  	   nmaxactposinf = SCIPmatrixGetRowNMaxActPosInf(matrix, row);
322  	   nminactneginf = SCIPmatrixGetRowNMinActNegInf(matrix, row);
323  	   nminactposinf = SCIPmatrixGetRowNMinActPosInf(matrix, row);
324  	
325  	   maxactivity = SCIPmatrixGetRowMaxActivity(matrix, row);
326  	   minactivity = SCIPmatrixGetRowMinActivity(matrix, row);
327  	
328  	   if( val >= 0.0 )
329  	   {
330  	      if( SCIPisInfinity(scip, ub) )
331  	      {
332  	         assert(nmaxactposinf >= 1);
333  	         if( nmaxactposinf == 1 && nmaxactneginf == 0 )
334  	            *maxresactivity = getMaxActivitySingleRowWithoutCol(scip, matrix, row, col);
335  	         else
336  	         {
337  	            *maxresactivity = SCIPinfinity(scip);
338  	            *ismaxsettoinfinity = TRUE;
339  	         }
340  	      }
341  	      else
342  	      {
343  	         if( (nmaxactneginf + nmaxactposinf) > 0 )
344  	         {
345  	            *maxresactivity = SCIPinfinity(scip);
346  	            *ismaxsettoinfinity = TRUE;
347  	         }
348  	         else
349  	            *maxresactivity = maxactivity - val * ub;
350  	      }
351  	
352  	      if( SCIPisInfinity(scip, -lb) )
353  	      {
354  	         assert(nminactneginf >= 1);
355  	         if( nminactneginf == 1 && nminactposinf == 0 )
356  	            *minresactivity = getMinActivitySingleRowWithoutCol(scip, matrix, row, col);
357  	         else
358  	         {
359  	            *minresactivity = -SCIPinfinity(scip);
360  	            *isminsettoinfinity = TRUE;
361  	         }
362  	      }
363  	      else
364  	      {
365  	         if( (nminactneginf + nminactposinf) > 0 )
366  	         {
367  	            *minresactivity = -SCIPinfinity(scip);
368  	            *isminsettoinfinity = TRUE;
369  	         }
370  	         else
371  	            *minresactivity = minactivity - val * lb;
372  	      }
373  	   }
374  	   else
375  	   {
376  	      if( SCIPisInfinity(scip, -lb) )
377  	      {
378  	         assert(nmaxactneginf >= 1);
379  	         if( nmaxactneginf == 1 && nmaxactposinf == 0 )
380  	            *maxresactivity = getMaxActivitySingleRowWithoutCol(scip, matrix, row, col);
381  	         else
382  	         {
383  	            *maxresactivity = SCIPinfinity(scip);
384  	            *ismaxsettoinfinity = TRUE;
385  	         }
386  	      }
387  	      else
388  	      {
389  	         if( (nmaxactneginf + nmaxactposinf) > 0 )
390  	         {
391  	            *maxresactivity = SCIPinfinity(scip);
392  	            *ismaxsettoinfinity = TRUE;
393  	         }
394  	         else
395  	            *maxresactivity = maxactivity - val * lb;
396  	      }
397  	
398  	      if( SCIPisInfinity(scip, ub) )
399  	      {
400  	         assert(nminactposinf >= 1);
401  	         if( nminactposinf == 1 && nminactneginf == 0 )
402  	            *minresactivity = getMinActivitySingleRowWithoutCol(scip, matrix, row, col);
403  	         else
404  	         {
405  	            *minresactivity = -SCIPinfinity(scip);
406  	            *isminsettoinfinity = TRUE;
407  	         }
408  	      }
409  	      else
410  	      {
411  	         if( (nminactneginf + nminactposinf) > 0 )
412  	         {
413  	            *minresactivity = -SCIPinfinity(scip);
414  	            *isminsettoinfinity = TRUE;
415  	         }
416  	         else
417  	            *minresactivity = minactivity - val * ub;
418  	      }
419  	   }
420  	}
421  	
422  	/** calculate the upper and lower bound of one variable from one row */
423  	static
424  	void getVarBoundsOfRow(
425  	   SCIP*                 scip,               /**< SCIP main data structure */
426  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
427  	   int                   col,                /**< column index of variable */
428  	   int                   row,                /**< row index */
429  	   SCIP_Real             val,                /**< coefficient of this column in this row */
430  	   SCIP_Real*            rowub,              /**< upper bound of row */
431  	   SCIP_Bool*            ubfound,            /**< flag indicating that an upper bound was calculated */
432  	   SCIP_Real*            rowlb,              /**< lower bound of row */
433  	   SCIP_Bool*            lbfound             /**< flag indicating that a lower bound was caluclated */
434  	   )
435  	{
436  	   SCIP_Bool isminsettoinfinity;
437  	   SCIP_Bool ismaxsettoinfinity;
438  	   SCIP_Real minresactivity;
439  	   SCIP_Real maxresactivity;
440  	   SCIP_Real lhs;
441  	   SCIP_Real rhs;
442  	
443  	   assert(rowub != NULL);
444  	   assert(ubfound != NULL);
445  	   assert(rowlb != NULL);
446  	   assert(lbfound != NULL);
447  	
448  	   *rowub = SCIPinfinity(scip);
449  	   *ubfound = FALSE;
450  	   *rowlb = -SCIPinfinity(scip);
451  	   *lbfound = FALSE;
452  	
453  	   getMinMaxActivityResiduals(scip, matrix, col, row, val,
454  	      &minresactivity, &maxresactivity,
455  	      &isminsettoinfinity, &ismaxsettoinfinity);
456  	
457  	   lhs = SCIPmatrixGetRowLhs(matrix, row);
458  	   rhs = SCIPmatrixGetRowRhs(matrix, row);
459  	
460  	   if( val > 0.0 )
461  	   {
462  	      if( !isminsettoinfinity && !SCIPisInfinity(scip, rhs) )
463  	      {
464  	         *rowub = (rhs - minresactivity) / val;
465  	         *ubfound = TRUE;
466  	      }
467  	
468  	      if( !ismaxsettoinfinity && !SCIPisInfinity(scip, -lhs) )
469  	      {
470  	         *rowlb = (lhs - maxresactivity) / val;
471  	         *lbfound = TRUE;
472  	      }
473  	   }
474  	   else
475  	   {
476  	      if( !ismaxsettoinfinity && !SCIPisInfinity(scip, -lhs) )
477  	      {
478  	         *rowub = (lhs - maxresactivity) / val;
479  	         *ubfound = TRUE;
480  	      }
481  	
482  	      if( !isminsettoinfinity && !SCIPisInfinity(scip, rhs) )
483  	      {
484  	         *rowlb = (rhs - minresactivity) / val;
485  	         *lbfound = TRUE;
486  	      }
487  	   }
488  	}
489  	
490  	/** detect implied variable bounds */
491  	static
492  	void getImpliedBounds(
493  	   SCIP*                 scip,               /**< SCIP main data structure */
494  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
495  	   int                   col,                /**< column index for implied free test */
496  	   SCIP_Bool*            ubimplied,          /**< flag indicating an implied upper bound */
497  	   SCIP_Bool*            lbimplied           /**< flag indicating an implied lower bound */
498  	   )
499  	{
500  	   SCIP_Real* valpnt;
501  	   int* colpnt;
502  	   int* colend;
503  	   SCIP_Real impliedub;
504  	   SCIP_Real impliedlb;
505  	   SCIP_Real ub;
506  	   SCIP_Real lb;
507  	
508  	   assert(scip != NULL);
509  	   assert(matrix != NULL);
510  	   assert(ubimplied != NULL);
511  	   assert(lbimplied != NULL);
512  	
513  	   *ubimplied = FALSE;
514  	   impliedub = SCIPinfinity(scip);
515  	
516  	   *lbimplied = FALSE;
517  	   impliedlb = -SCIPinfinity(scip);
518  	
519  	   ub =  SCIPmatrixGetColUb(matrix, col);
520  	   lb =  SCIPmatrixGetColLb(matrix, col);
521  	
522  	   colpnt = SCIPmatrixGetColIdxPtr(matrix, col);
523  	   colend = colpnt + SCIPmatrixGetColNNonzs(matrix, col);
524  	   valpnt = SCIPmatrixGetColValPtr(matrix, col);
525  	   for( ; (colpnt < colend); colpnt++, valpnt++ )
526  	   {
527  	      SCIP_Real rowub;
528  	      SCIP_Bool ubfound;
529  	      SCIP_Real rowlb;
530  	      SCIP_Bool lbfound;
531  	
532  	      getVarBoundsOfRow(scip, matrix, col, *colpnt, *valpnt, &rowub, &ubfound, &rowlb, &lbfound);
533  	
534  	      if( ubfound && (rowub < impliedub) )
535  	         impliedub = rowub;
536  	
537  	      if( lbfound && (rowlb > impliedlb) )
538  	         impliedlb = rowlb;
539  	   }
540  	
541  	   /* we consider +/-inf bounds as implied bounds */
542  	   if( SCIPisInfinity(scip, ub) ||
543  	      (!SCIPisInfinity(scip, ub) && SCIPisLE(scip, impliedub, ub)) )
544  	      *ubimplied = TRUE;
545  	
546  	   if( SCIPisInfinity(scip, -lb) ||
547  	      (!SCIPisInfinity(scip, -lb) && SCIPisGE(scip, impliedlb, lb)) )
548  	      *lbimplied = TRUE;
549  	}
550  	
551  	/** y = weight1 * var[colidx1] + var[colidx2] */
552  	static
553  	SCIP_RETCODE aggregation(
554  	   SCIP*                 scip,               /**< SCIP datastructure */
555  	   SCIP_MATRIX*          matrix,             /**< matrix datastructure */
556  	   SCIP_PRESOLDATA*      presoldata,         /**< presolver data */
557  	   SCIP_VAR**            vars,               /**< the current variables */
558  	   int                   colidx1,            /**< one of the indexes of column to try nonzero cancellation for */
559  	   int                   colidx2,            /**< one of the indexes of column to try nonzero cancellation for */
560  	   SCIP_Bool             isimpliedfree,      /**< is the aggregated variable implied free? */
561  	   SCIP_Real             weight1             /**< weight variable one in the aggregated expression */
562  	   )
563  	{
564  	   SCIP_VAR* tmpvars[2];
565  	   SCIP_Real coefs[2];
566  	   char newvarname[SCIP_MAXSTRLEN];
567  	   char newconsname[SCIP_MAXSTRLEN];
568  	   SCIP_CONS* newcons;
569  	   SCIP_VAR* aggregatedvar;
570  	   SCIP_VAR* newvar;
571  	   SCIP_VARTYPE newvartype;
572  	   SCIP_Real constant;
573  	   SCIP_Real newlb;
574  	   SCIP_Real newub;
575  	   SCIP_Real lhs;
576  	   SCIP_Real rhs;
577  	   SCIP_Bool infeasible;
578  	   SCIP_Bool aggregated;
579  	#ifndef NDEBUG
580  	   if( isimpliedfree )
581  	   {
582  	      SCIP_Bool lbimplied;
583  	      SCIP_Bool ubimplied;
584  	
585  	      getImpliedBounds(scip, matrix, colidx2, &ubimplied, &lbimplied);
586  	      assert(lbimplied && ubimplied);
587  	   }
588  	#endif
589  	
590  	   assert( !SCIPisZero(scip, weight1) );
591  	
592  	   presoldata->naggregated += 1;
593  	   aggregatedvar = vars[colidx2];
594  	
595  	   /* if the variable is implied free, we make sure that the columns bounds are removed,
596  	    * so that subsequent checks for implied bounds do not interfere with the exploitation
597  	    * of this variables implied bounds
598  	    */
599  	   if( isimpliedfree )
600  	   {
601  	      SCIPdebugMsg(scip, "remove column bounds of column %d\n", colidx2);
602  	      SCIPmatrixRemoveColumnBounds(scip, matrix, colidx2);
603  	   }
604  	
605  	   assert(!SCIPdoNotMultaggrVar(scip, aggregatedvar));
606  	
607  	   (void) SCIPsnprintf(newvarname, SCIP_MAXSTRLEN, "dualsparsifyvar_%d", presoldata->naggregated);
608  	
609  	   constant = 0.0;
610  	
611  	   if( weight1 > 0.0 )
612  	   {
613  	      if( SCIPisInfinity(scip, -SCIPvarGetLbGlobal(vars[colidx1])) ||
614  	         SCIPisInfinity(scip, -SCIPvarGetLbGlobal(vars[colidx2])) )
615  	         newlb = -SCIPinfinity(scip);
616  	      else
617  	         newlb = weight1 * SCIPvarGetLbGlobal(vars[colidx1]) + SCIPvarGetLbGlobal(vars[colidx2]);
618  	
619  	      if( SCIPisInfinity(scip, SCIPvarGetUbGlobal(vars[colidx1])) ||
620  	         SCIPisInfinity(scip, SCIPvarGetUbGlobal(vars[colidx2])) )
621  	         newub = SCIPinfinity(scip);
622  	      else
623  	         newub = weight1 * SCIPvarGetUbGlobal(vars[colidx1]) + SCIPvarGetUbGlobal(vars[colidx2]);
624  	   }
625  	   else
626  	   {
627  	      if( SCIPisInfinity(scip, SCIPvarGetUbGlobal(vars[colidx1])) ||
628  	         SCIPisInfinity(scip, -SCIPvarGetLbGlobal(vars[colidx2])) )
629  	         newlb = -SCIPinfinity(scip);
630  	      else
631  	         newlb = weight1 * SCIPvarGetUbGlobal(vars[colidx1]) + SCIPvarGetLbGlobal(vars[colidx2]);
632  	
633  	      if( SCIPisInfinity(scip, SCIPvarGetLbGlobal(vars[colidx1])) ||
634  	         SCIPisInfinity(scip, SCIPvarGetUbGlobal(vars[colidx2])) )
635  	         newub = SCIPinfinity(scip);
636  	      else
637  	         newub = weight1 * SCIPvarGetLbGlobal(vars[colidx1]) + SCIPvarGetUbGlobal(vars[colidx2]);
638  	   }
639  	
640  	   if( SCIPvarIsIntegral(aggregatedvar) )
641  	      newvartype = (SCIPvarGetType(aggregatedvar) == SCIP_VARTYPE_IMPLINT) ?
642  	         SCIP_VARTYPE_IMPLINT : SCIP_VARTYPE_INTEGER;
643  	   else
644  	      newvartype = SCIP_VARTYPE_CONTINUOUS;
645  	
646  	   lhs = SCIPvarGetLbGlobal(vars[colidx2]);
647  	   rhs = SCIPvarGetUbGlobal(vars[colidx2]);
648  	
649  	   SCIP_CALL( SCIPcreateVar(scip, &newvar, newvarname, newlb, newub, 0.0, newvartype,
650  	         SCIPvarIsInitial(aggregatedvar), SCIPvarIsRemovable(aggregatedvar), NULL, NULL, NULL, NULL, NULL) );
651  	   SCIP_CALL( SCIPaddVar(scip, newvar) );
652  	
653  	   /* set the debug solution value for the new variable */
654  	#ifdef WITH_DEBUG_SOLUTION
655  	   if( SCIPdebugIsMainscip(scip) )
656  	   {
657  	      SCIP_Real val1;
658  	      SCIP_Real val2;
659  	
660  	      SCIP_CALL( SCIPdebugGetSolVal(scip, vars[colidx1], &val1) );
661  	      SCIP_CALL( SCIPdebugGetSolVal(scip, vars[colidx2], &val2) );
662  	      SCIP_CALL( SCIPdebugAddSolVal(scip, newvar, weight1 * val1 + val2) );
663  	
664  	      SCIPdebugMsg(scip, "set debug solution value of %s to %g\n", SCIPvarGetName(newvar), weight1 * val1 + val2);
665  	   }
666  	#endif
667  	
668  	   tmpvars[0] = vars[colidx1];
669  	   tmpvars[1] = newvar;
670  	   coefs[0] = -weight1;
671  	   coefs[1] = 1.0;
672  	
673  	   SCIP_CALL( SCIPmultiaggregateVar(scip, aggregatedvar, 2, tmpvars, coefs, constant, &infeasible, &aggregated) );
674  	
675  	   assert(!infeasible);
676  	   assert(aggregated);
677  	
678  	   vars[colidx2] = newvar;
679  	
680  	   /* create a linear constraint that ensures that var[colidx2].lb <= y - weight1 * var[colidx1] <= var[colidx2].ub;
681  	    * note that it might happen that vars[colidx2] is not implied free even though it has infinite bounds because
682  	    * getImpliedBounds() considers infinite bounds to be implied
683  	    */
684  	   if( !isimpliedfree && (!SCIPisInfinity(scip, rhs) || !SCIPisInfinity(scip, -lhs)) )
685  	   {
686  	      SCIPdebugMsg(scip, "create a linear constraint to ensure %g <= %g %s + %g %s <= %g\n", lhs, coefs[0], SCIPvarGetName(tmpvars[0]),
687  	         coefs[1], SCIPvarGetName(tmpvars[1]), rhs);
688  	      (void) SCIPsnprintf(newconsname, SCIP_MAXSTRLEN, "dualsparsifycons_%d", presoldata->naggregated);
689  	
690  	      SCIP_CALL( SCIPcreateConsLinear(scip, &newcons, newconsname, 2, tmpvars, coefs,
691  	            lhs, rhs, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
692  	      SCIP_CALL( SCIPaddCons(scip, newcons) );
693  	
694  	      SCIPdebugPrintCons(scip, newcons, NULL);
695  	
696  	      SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
697  	   }
698  	
699  	   SCIP_CALL( SCIPreleaseVar(scip, &newvar) );
700  	
701  	   return SCIP_OKAY;
702  	}
703  	
704  	/** try nonzero cancellation for given column */
705  	static
706  	SCIP_RETCODE cancelCol(
707  	   SCIP*                 scip,               /**< SCIP datastructure */
708  	   SCIP_MATRIX*          matrix,             /**< the constraint matrix */
709  	   SCIP_PRESOLDATA*      presoldata,         /**< presolver data */
710  	   SCIP_HASHTABLE*       pairtable,          /**< the hashtable containing COLCONSPAIR's of equations */
711  	   SCIP_Bool*            ishashingcols,      /**< array to indicates whether it is impliedfree or not */
712  	   SCIP_VAR**            vars,               /**< array to store the current variables */
713  	   SCIP_Bool*            isblockedvar,       /**< array to indicates whether it is blocked or not */
714  	   int                   colidx,             /**< index of column to try nonzero cancellation for */
715  	   int                   maxcontfillin,      /**< maximal fill-in allowed for continuous variables */
716  	   int                   maxintfillin,       /**< maximal fill-in allowed for integral variables */
717  	   int                   maxbinfillin,       /**< maximal fill-in allowed for binary variables */
718  	   int                   maxconsiderednonzeros, /**< maximal number of nonzeros to consider for cancellation */
719  	   SCIP_Bool             preserveintcoefs,   /**< only perform nonzero cancellation if integrality of coefficients is preserved? */
720  	   SCIP_Longint*         nuseless,           /**< pointer to update number of useless hashtable retrieves */
721  	   int*                  nchgcoefs,          /**< pointer to update number of changed coefficients */
722  	   int*                  ncanceled,          /**< pointer to update number of canceled nonzeros */
723  	   int*                  nfillin,            /**< pointer to update the produced fill-in */
724  	   SCIP_Bool             isaddedcons         /**< whether a linear constraint required to added to keep the validity */
725  	   )
726  	{
727  	   SCIP_VAR* cancelvar;
728  	   SCIP_Real* cancelcolvals;
729  	   SCIP_Real* colvalptr;
730  	   SCIP_Real* tmpvals;
731  	   SCIP_Real* scores;
732  	   int* cancelcolinds;
733  	   int* colidxptr;
734  	   int* tmpinds;
735  	   SCIP_Real bestcancelrate;
736  	   SCIP_Real bestscale;
737  	   SCIP_Real ncols;
738  	   SCIP_Bool colishashing;
739  	   SCIP_Bool swapped = FALSE;
740  	   int cancelcollen;
741  	   int bestnfillin;
742  	   int nretrieves;
743  	   int maxfillin;
744  	   int bestcand;
745  	   int nchgcoef;
746  	
747  	   ncols = SCIPmatrixGetNColumns(matrix);
748  	   colishashing = ishashingcols[colidx];
749  	   cancelcollen = SCIPmatrixGetColNNonzs(matrix, colidx);
750  	   colidxptr = SCIPmatrixGetColIdxPtr(matrix, colidx);
751  	   colvalptr = SCIPmatrixGetColValPtr(matrix, colidx);
752  	   cancelvar = vars[colidx];
753  	
754  	   if( SCIPvarIsIntegral(cancelvar) )
755  	   {
756  	      if( SCIPvarIsBinary(cancelvar) )
757  	         maxfillin = maxbinfillin;
758  	      else
759  	         maxfillin = maxintfillin;
760  	   }
761  	   else
762  	      maxfillin = maxcontfillin;
763  	
764  	   SCIP_CALL( SCIPduplicateBufferArray(scip, &cancelcolinds, colidxptr, cancelcollen) );
765  	   SCIP_CALL( SCIPduplicateBufferArray(scip, &cancelcolvals, colvalptr, cancelcollen) );
766  	   SCIP_CALL( SCIPallocBufferArray(scip, &tmpinds, cancelcollen) );
767  	   SCIP_CALL( SCIPallocBufferArray(scip, &tmpvals, cancelcollen) );
768  	   SCIP_CALL( SCIPallocBufferArray(scip, &scores, cancelcollen) );
769  	
770  	   nchgcoef = 0;
771  	   nretrieves = 0;
772  	   while( TRUE ) /*lint !e716 */
773  	   {
774  	      COLCONSPAIR colconspair;
775  	      int maxlen;
776  	      int i;
777  	      int j;
778  	
779  	      bestcand = -1;
780  	      bestnfillin = 0;
781  	      bestscale = 1.0;
782  	      bestcancelrate = 0.0;
783  	
784  	      /* sort the rows non-decreasingly by number of nonzeros
785  	       * if the number of nonzeros, we use the colindex as tie-breaker
786  	       */
787  	      for( i = 0; i < cancelcollen; ++i )
788  	      {
789  	         tmpinds[i] = i;
790  	         scores[i] = -SCIPmatrixGetRowNNonzs(matrix, cancelcolinds[i]) - 1.0 * cancelcolinds[i] / (ncols);
791  	      }
792  	      SCIPsortRealInt(scores, tmpinds, cancelcollen);
793  	
794  	      maxlen = cancelcollen;
795  	      if( maxconsiderednonzeros >= 0 )
796  	         maxlen = MIN(cancelcollen, maxconsiderednonzeros);
797  	
798  	      for( i = 0; i < maxlen; ++i )
799  	      {
800  	         for( j = i + 1; j < maxlen; ++j )
801  	         {
802  	            COLCONSPAIR* hashingcolconspair;
803  	            SCIP_VAR* hashingcolvar;
804  	            SCIP_Real* hashingcolvals;
805  	            int* hashingcolinds;
806  	            SCIP_Real hashingcollb;
807  	            SCIP_Real hashingcolub;
808  	            SCIP_Real cancelrate;
809  	            SCIP_Real rowlhs;
810  	            SCIP_Real rowrhs;
811  	            SCIP_Real scale;
812  	            SCIP_Bool hashingcolisbin;
813  	            SCIP_Bool abortpair;
814  	            int hashingcollen;
815  	            int ntotfillin;
816  	            int ncancel;
817  	            int a,b;
818  	            int i1,i2;
819  	
820  	            i1 = tmpinds[i];
821  	            i2 = tmpinds[j];
822  	
823  	            assert(cancelcolinds[i] < cancelcolinds[j]);
824  	
825  	            if( cancelcolinds[i1] < cancelcolinds[i2] )
826  	            {
827  	               colconspair.consindex1 = cancelcolinds[i1];
828  	               colconspair.consindex2 = cancelcolinds[i2];
829  	               colconspair.conscoef1 = cancelcolvals[i1];
830  	               colconspair.conscoef2 = cancelcolvals[i2];
831  	            }
832  	            else
833  	            {
834  	               colconspair.consindex1 = cancelcolinds[i2];
835  	               colconspair.consindex2 = cancelcolinds[i1];
836  	               colconspair.conscoef1 = cancelcolvals[i2];
837  	               colconspair.conscoef2 = cancelcolvals[i1];
838  	            }
839  	
840  	            hashingcolconspair = (COLCONSPAIR*)SCIPhashtableRetrieve(pairtable, (void*) &colconspair);
841  	            nretrieves++;
842  	
843  	            if( hashingcolconspair == NULL ||
844  	               hashingcolconspair->colindex == colidx || isblockedvar[hashingcolconspair->colindex] )
845  	               continue;
846  	
847  	            /* if the column we want to cancel is a hashing column (which we stored for canceling other columns),
848  	             * we will only use the hashing columns for canceling with less nonzeros and if the number of nonzeros
849  	             * is equal we use the colindex as tie-breaker to avoid cyclic nonzero cancellation
850  	             */
851  	            hashingcollen = SCIPmatrixGetColNNonzs(matrix, hashingcolconspair->colindex);
852  	            if( colishashing && (cancelcollen < hashingcollen ||
853  	                  (cancelcollen == hashingcollen && colidx < hashingcolconspair->colindex)) )
854  	               continue;
855  	
856  	            hashingcolvals = SCIPmatrixGetColValPtr(matrix, hashingcolconspair->colindex);
857  	            hashingcolinds = SCIPmatrixGetColIdxPtr(matrix, hashingcolconspair->colindex);
858  	            hashingcolvar = vars[hashingcolconspair->colindex];
859  	            hashingcollb = SCIPvarGetLbGlobal(hashingcolvar);
860  	            hashingcolub = SCIPvarGetUbGlobal(hashingcolvar);
861  	            hashingcolisbin = (SCIPvarGetType(hashingcolvar) == SCIP_VARTYPE_BINARY) ||
862  	               (SCIPvarIsIntegral(hashingcolvar) && SCIPisZero(scip, hashingcollb) &&
863  	                  SCIPisEQ(scip, hashingcolub, 1.0));
864  	            scale = -colconspair.conscoef1 / hashingcolconspair->conscoef1;
865  	
866  	            if( SCIPisZero(scip, scale) )
867  	               continue;
868  	
869  	            if( REALABS(scale) > MAXSCALE )
870  	               continue;
871  	
872  	            /* @todo do more reduction if knspsack constraint handler supports downgrading constraint,
873  	             * i.e., converting into a linear constraint
874  	             */
875  	            if( hashingcolisbin )
876  	               continue;
877  	            else if( SCIPvarIsIntegral(hashingcolvar) )
878  	            {
879  	               if( SCIPvarIsIntegral(cancelvar) )
880  	               {
881  	                  /* skip if the hashing variable is an integer variable and
882  	                   * the canceled variable is an implicit integer variable
883  	                   */
884  	                  if( (SCIPvarGetType(hashingcolvar) != SCIP_VARTYPE_IMPLINT) &&
885  	                     (SCIPvarGetType(cancelvar) == SCIP_VARTYPE_IMPLINT) )
886  	                     continue;
887  	
888  	                  /* skip if the scale is non-integral */
889  	                  if( !SCIPisIntegral(scip, scale) )
890  	                     continue;
891  	
892  	                  /* round scale to be exactly integral */
893  	                  scale = floor(scale + 0.5);
894  	               }
895  	               /* skip if the canceled variable is a continuous variable */
896  	               else
897  	                  continue;
898  	            }
899  	
900  	            a = 0;
901  	            b = 0;
902  	            ncancel = 0;
903  	            ntotfillin = 0;
904  	            abortpair = FALSE;
905  	
906  	            while( a < cancelcollen && b < hashingcollen )
907  	            {
908  	               if( cancelcolinds[a] == hashingcolinds[b] )
909  	               {
910  	                  SCIP_Real newcoef;
911  	
912  	                  newcoef = cancelcolvals[a] + scale * hashingcolvals[b];
913  	
914  	                  /* check if coefficient is canceled */
915  	                  if( SCIPisZero(scip, newcoef) )
916  	                  {
917  	                     ++ncancel;
918  	                  }
919  	                  /* otherwise, check if integral coefficients are preserved if the column is integral */
920  	                  else if( (preserveintcoefs && SCIPvarIsIntegral(cancelvar) &&
921  	                        SCIPisIntegral(scip, cancelcolvals[a]) && !SCIPisIntegral(scip, newcoef)) )
922  	                  {
923  	                     abortpair = TRUE;
924  	                     break;
925  	                  }
926  	                  /* finally, check if locks could be modified in a bad way due to flipped signs */
927  	                  else if( COPYSIGN(1.0, newcoef) != COPYSIGN(1.0, cancelcolvals[a]) ) /*lint !e777*/
928  	                  {
929  	                     /* do not flip signs for non-canceled coefficients if this adds a lock to a variable that
930  	                      * had at most one lock in that direction before, except if the other direction gets unlocked
931  	                      */
932  	                     rowrhs = SCIPmatrixGetRowRhs(matrix, cancelcolinds[a]);
933  	                     rowlhs = SCIPmatrixGetRowLhs(matrix, cancelcolinds[a]);
934  	                     if( (cancelcolvals[a] > 0.0 && ! SCIPisInfinity(scip, rowrhs)) ||
935  	                        (cancelcolvals[a] < 0.0 && ! SCIPisInfinity(scip, -rowlhs)) )
936  	                     {
937  	                        /* if we get into this case the variable had a positive coefficient in a <= constraint or
938  	                         * a negative coefficient in a >= constraint, e.g. an uplock. If this was the only uplock
939  	                         * we do not abort their cancelling, otherwise we abort if we had a single or no downlock
940  	                         * and add one
941  	                         */
942  	                        if( presoldata->preservegoodlocks && (SCIPmatrixGetColNUplocks(matrix, colidx) > 1 &&
943  	                              SCIPmatrixGetColNDownlocks(matrix, colidx) <= 1) )
944  	                        {
945  	                           abortpair = TRUE;
946  	                           break;
947  	                        }
948  	                     }
949  	
950  	                     if( (cancelcolvals[a] < 0.0 && ! SCIPisInfinity(scip, rowrhs)) ||
951  	                        (cancelcolvals[a] > 0.0 && ! SCIPisInfinity(scip, -rowlhs)) )
952  	                     {
953  	                        /* symmetric case where the variable had a downlock */
954  	                        if( presoldata->preservegoodlocks && (SCIPmatrixGetColNDownlocks(matrix, colidx) > 1 &&
955  	                              SCIPmatrixGetColNUplocks(matrix, colidx) <= 1) )
956  	                        {
957  	                           abortpair = TRUE;
958  	                           break;
959  	                        }
960  	                     }
961  	                  }
962  	
963  	                  ++a;
964  	                  ++b;
965  	               }
966  	               else if( cancelcolinds[a] < hashingcolinds[b] )
967  	               {
968  	                  ++a;
969  	               }
970  	               else
971  	               {
972  	                  SCIP_Real newcoef;
973  	
974  	                  newcoef = scale * hashingcolvals[b];
975  	                  rowrhs = SCIPmatrixGetRowRhs(matrix, hashingcolinds[b]);
976  	                  rowlhs = SCIPmatrixGetRowLhs(matrix, hashingcolinds[b]);
977  	
978  	                  if( (newcoef > 0.0 && ! SCIPisInfinity(scip, rowrhs)) ||
979  	                     (newcoef < 0.0 && ! SCIPisInfinity(scip, -rowlhs)) )
980  	                  {
981  	                     if( presoldata->preservegoodlocks && SCIPmatrixGetColNUplocks(matrix, colidx) <= 1 )
982  	                     {
983  	                        abortpair = TRUE;
984  	                        ++b;
985  	                        break;
986  	                     }
987  	                  }
988  	
989  	                  if( (newcoef < 0.0 && ! SCIPisInfinity(scip, rowrhs)) ||
990  	                     (newcoef > 0.0 && ! SCIPisInfinity(scip, -rowlhs)) )
991  	                  {
992  	                     if( presoldata->preservegoodlocks && SCIPmatrixGetColNDownlocks(matrix, colidx) <= 1 )
993  	                     {
994  	                        abortpair = TRUE;
995  	                        ++b;
996  	                        break;
997  	                     }
998  	                  }
999  	
1000 	                  ++b;
1001 	
1002 	                  if( ++ntotfillin > maxfillin )
1003 	                  {
1004 	                     abortpair = TRUE;
1005 	                     break;
1006 	                  }
1007 	               }
1008 	            }
1009 	
1010 	            if( abortpair )
1011 	               continue;
1012 	
1013 	            while( b < hashingcollen )
1014 	            {
1015 	               ++b;
1016 	
1017 	               if( ++ntotfillin > maxfillin )
1018 	                  break;
1019 	            }
1020 	CHECKFILLINAGAIN:
1021 	            if( ntotfillin > maxfillin || ntotfillin >= ncancel )
1022 	               continue;
1023 	
1024 	            cancelrate = (ncancel - ntotfillin) / (SCIP_Real) cancelcollen;
1025 	
1026 	            /* if a linear constraint is needed to keep the validity, we require a large nonzero cancellation */
1027 	            if( isaddedcons && (ncancel - ntotfillin < presoldata->mineliminatednonzeros) )
1028 	               continue;
1029 	
1030 	            if( cancelrate > bestcancelrate )
1031 	            {
1032 	               if( ishashingcols[hashingcolconspair->colindex] )
1033 	               {
1034 	                  SCIP_Bool lbimplied;
1035 	                  SCIP_Bool ubimplied;
1036 	
1037 	                  /* recompute whether a variable is still implied free; after some previous multi-aggregations of
1038 	                   * some variables, it might be that other variables that are contained in the same linear rows of the
1039 	                   * matrix are not implied free anymore (see #2971)
1040 	                   */
1041 	                  getImpliedBounds(scip, matrix, hashingcolconspair->colindex, &ubimplied, &lbimplied);
1042 	
1043 	                  if( !lbimplied || !ubimplied )
1044 	                  {
1045 	                     ishashingcols[hashingcolconspair->colindex] = FALSE;
1046 	                     ntotfillin += 2;
1047 	                     goto CHECKFILLINAGAIN;
1048 	                  }
1049 	               }
1050 	
1051 	               bestnfillin = ntotfillin;
1052 	               bestcand = hashingcolconspair->colindex;
1053 	               bestscale = scale;
1054 	               bestcancelrate = cancelrate;
1055 	
1056 	               /* stop looking if the current candidate does not create any fill-in or alter coefficients */
1057 	               if( cancelrate == 1.0 )
1058 	                  break;
1059 	            }
1060 	
1061 	            /* we accept the best candidate immediately if it does not create any fill-in or alter coefficients */
1062 	            if( bestcand != -1 && bestcancelrate == 1.0 )
1063 	               break;
1064 	         }
1065 	      }
1066 	
1067 	      if( bestcand != -1 )
1068 	      {
1069 	         SCIP_Real* hashingcolvals;
1070 	         int* hashingcolinds;
1071 	         int hashingcollen;
1072 	         int tmpcollen;
1073 	         int a;
1074 	         int b;
1075 	
1076 	         SCIPdebugMsg(scip, "cancelcol %d (%s) candidate column %d (%s) (bestcancelrate = %g, bestscale = %g)\n",
1077 	            colidx, SCIPvarGetName(cancelvar), bestcand, SCIPvarGetName(vars[bestcand]), bestcancelrate, bestscale);
1078 	
1079 	         hashingcolvals = SCIPmatrixGetColValPtr(matrix, bestcand);
1080 	         hashingcolinds = SCIPmatrixGetColIdxPtr(matrix, bestcand);
1081 	         hashingcollen = SCIPmatrixGetColNNonzs(matrix, bestcand);
1082 	
1083 	         a = 0;
1084 	         b = 0;
1085 	         tmpcollen = 0;
1086 	
1087 	         while( a < cancelcollen && b < hashingcollen )
1088 	         {
1089 	            if( cancelcolinds[a] == hashingcolinds[b] )
1090 	            {
1091 	               SCIP_Real val = cancelcolvals[a] + bestscale * hashingcolvals[b];
1092 	
1093 	               if( !SCIPisZero(scip, val) )
1094 	               {
1095 	                  tmpinds[tmpcollen] = cancelcolinds[a];
1096 	                  tmpvals[tmpcollen] = val;
1097 	                  ++tmpcollen;
1098 	               }
1099 	               ++nchgcoef;
1100 	
1101 	               ++a;
1102 	               ++b;
1103 	            }
1104 	            else if( cancelcolinds[a] < hashingcolinds[b] )
1105 	            {
1106 	               tmpinds[tmpcollen] = cancelcolinds[a];
1107 	               tmpvals[tmpcollen] = cancelcolvals[a];
1108 	               ++tmpcollen;
1109 	               ++a;
1110 	            }
1111 	            else
1112 	            {
1113 	               tmpinds[tmpcollen] = hashingcolinds[b];
1114 	               tmpvals[tmpcollen] = hashingcolvals[b] * bestscale;
1115 	               ++nchgcoef;
1116 	               ++tmpcollen;
1117 	               ++b;
1118 	            }
1119 	         }
1120 	
1121 	         while( a < cancelcollen )
1122 	         {
1123 	            tmpinds[tmpcollen] = cancelcolinds[a];
1124 	            tmpvals[tmpcollen] = cancelcolvals[a];
1125 	            ++tmpcollen;
1126 	            ++a;
1127 	         }
1128 	
1129 	         while( b < hashingcollen )
1130 	         {
1131 	            tmpinds[tmpcollen] = hashingcolinds[b];
1132 	            tmpvals[tmpcollen] = hashingcolvals[b] * bestscale;
1133 	            ++nchgcoef;
1134 	            ++tmpcollen;
1135 	            ++b;
1136 	         }
1137 	
1138 	         /* update fill-in counter */
1139 	         *nfillin += bestnfillin;
1140 	
1141 	         /* swap the temporary arrays so that the cancelcolinds and cancelcolvals arrays, contain the new
1142 	          * changed column, and the tmpinds and tmpvals arrays can be overwritten in the next iteration
1143 	          */
1144 	         SCIPswapPointers((void**) &tmpinds, (void**) &cancelcolinds);
1145 	         SCIPswapPointers((void**) &tmpvals, (void**) &cancelcolvals);
1146 	         swapped = ! swapped;
1147 	         cancelcollen = tmpcollen;
1148 	         SCIP_CALL( aggregation(scip, matrix, presoldata, vars, colidx, bestcand, ishashingcols[bestcand], -bestscale) );
1149 	
1150 	         /* the newly created variable is now at the position bestcand and is assumed to have the same coefficients.
1151 	          * this is not the case if the variable is not implied free since then a new constraint was added and the
1152 	          * nonzero fillin would not be counted correctly if we do not block this variable
1153 	          */
1154 	         if( !ishashingcols[bestcand] )
1155 	            isblockedvar[bestcand] = TRUE;
1156 	      }
1157 	      else
1158 	         break;
1159 	   }
1160 	
1161 	   if( nchgcoef != 0 )
1162 	   {
1163 	      /* update counters */
1164 	      *nchgcoefs += nchgcoef;
1165 	      *ncanceled += SCIPmatrixGetColNNonzs(matrix, colidx) - cancelcollen;
1166 	
1167 	      isblockedvar[colidx] = TRUE;
1168 	
1169 	      /* if successful, decrease the useless hashtable retrieves counter; the rationale here is that we want to keep
1170 	       * going if, after many useless calls that almost exceeded the budget, we finally reach a useful section; but we
1171 	       * don't allow a negative build-up for the case that the useful section is all at the beginning and we just want
1172 	       * to quit quickly afterwards
1173 	       */
1174 	      *nuseless -= nretrieves;
1175 	      *nuseless = MAX(*nuseless, 0);
1176 	   }
1177 	   else
1178 	   {
1179 	      /* if not successful, increase useless hashtable retrieves counter */
1180 	      *nuseless += nretrieves;
1181 	   }
1182 	
1183 	   SCIPfreeBufferArray(scip, &scores);
1184 	   if( swapped )
1185 	   {
1186 	      SCIPfreeBufferArray(scip, &cancelcolvals);
1187 	      SCIPfreeBufferArray(scip, &cancelcolinds);
1188 	      SCIPfreeBufferArray(scip, &tmpvals);
1189 	      SCIPfreeBufferArray(scip, &tmpinds);
1190 	   }
1191 	   else
1192 	   {
1193 	      SCIPfreeBufferArray(scip, &tmpvals);
1194 	      SCIPfreeBufferArray(scip, &tmpinds);
1195 	      SCIPfreeBufferArray(scip, &cancelcolvals);
1196 	      SCIPfreeBufferArray(scip, &cancelcolinds);
1197 	   }
1198 	
1199 	   return SCIP_OKAY;
1200 	}
1201 	
1202 	/** updates failure counter after one execution */
1203 	static
1204 	void updateFailureStatistic(
1205 	   SCIP_PRESOLDATA*      presoldata,         /**< presolver data */
1206 	   SCIP_Bool             success             /**< was this execution successful? */
1207 	   )
1208 	{
1209 	   assert(presoldata != NULL);
1210 	
1211 	   if( success )
1212 	   {
1213 	      presoldata->nfailures = 0;
1214 	      presoldata->nwaitingcalls = 0;
1215 	   }
1216 	   else
1217 	   {
1218 	      presoldata->nfailures++;
1219 	      presoldata->nwaitingcalls = (int)(presoldata->waitingfac*(SCIP_Real)presoldata->nfailures);
1220 	   }
1221 	}
1222 	
1223 	/*
1224 	 * Callback methods of presolver
1225 	 */
1226 	
1227 	/** copy method for constraint handler plugins (called when SCIP copies plugins) */
1228 	static
1229 	SCIP_DECL_PRESOLCOPY(presolCopyDualsparsify)
1230 	{
1231 	   SCIP_PRESOLDATA* presoldata;
1232 	
1233 	   assert(scip != NULL);
1234 	   assert(presol != NULL);
1235 	   assert(strcmp(SCIPpresolGetName(presol), PRESOL_NAME) == 0);
1236 	
1237 	   /* call inclusion method of presolver if copying is enabled */
1238 	   presoldata = SCIPpresolGetData(presol);
1239 	   assert(presoldata != NULL);
1240 	   if( presoldata->enablecopy )
1241 	   {
1242 	      SCIP_CALL( SCIPincludePresolDualsparsify(scip) );
1243 	   }
1244 	
1245 	   return SCIP_OKAY;
1246 	}
1247 	
1248 	
1249 	/** execution method of presolver */
1250 	static
1251 	SCIP_DECL_PRESOLEXEC(presolExecDualsparsify)
1252 	{  /*lint --e{715}*/
1253 	   SCIP_MATRIX* matrix;
1254 	   int* perm;
1255 	   int* colidxsorted;
1256 	   int* colsparsity;
1257 	   SCIP_Real* scores;
1258 	   COLCONSPAIR* conspairs;
1259 	   SCIP_HASHTABLE* pairtable;
1260 	   SCIP_PRESOLDATA* presoldata;
1261 	   SCIP_Bool* ishashingcols;
1262 	   SCIP_Bool* isblockedvar;
1263 	   SCIP_VAR** vars;
1264 	   SCIP_Longint maxuseless;
1265 	   SCIP_Longint nuseless;
1266 	   SCIP_Bool initialized;
1267 	   SCIP_Bool complete;
1268 	   SCIP_Bool infeasible;
1269 	   int ncols;
1270 	   int c;
1271 	   int i;
1272 	   int j;
1273 	   int conspairssize;
1274 	   int nconspairs;
1275 	   int numcancel;
1276 	   int nfillin;
1277 	
1278 	   assert(result != NULL);
1279 	
1280 	   *result = SCIP_DIDNOTRUN;
1281 	
1282 	   if( SCIPdoNotAggr(scip) )
1283 	      return SCIP_OKAY;
1284 	
1285 	   /* If restart is performed, some cuts will be tranformed into linear constraints.
1286 	    * However, SCIPmatrixCreate() only collects the original constraints (not the constraints transformed from cuts)
1287 	    * For this reason, we only perform this method in the first run of branch-and-cut.
1288 	    * */
1289 	   if( SCIPgetNRuns(scip) > 1 )
1290 	      return SCIP_OKAY;
1291 	
1292 	   presoldata = SCIPpresolGetData(presol);
1293 	
1294 	   if( presoldata->nwaitingcalls > 0 )
1295 	   {
1296 	      presoldata->nwaitingcalls--;
1297 	      SCIPdebugMsg(scip, "skipping dualsparsify: nfailures=%d, nwaitingcalls=%d\n", presoldata->nfailures,
1298 	         presoldata->nwaitingcalls);
1299 	      return SCIP_OKAY;
1300 	   }
1301 	
1302 	   SCIPdebugMsg(scip, "starting dualsparsify. . .\n");
1303 	   *result = SCIP_DIDNOTFIND;
1304 	
1305 	   matrix = NULL;
1306 	   SCIP_CALL( SCIPmatrixCreate(scip, &matrix, TRUE, &initialized, &complete, &infeasible,
1307 	         naddconss, ndelconss, nchgcoefs, nchgbds, nfixedvars) );
1308 	
1309 	   /* if infeasibility was detected during matrix creation or
1310 	    * matrix creation is incomplete, return here.
1311 	    */
1312 	   if( infeasible || !complete )
1313 	   {
1314 	      if( initialized )
1315 	         SCIPmatrixFree(scip, &matrix);
1316 	
1317 	      if( infeasible )
1318 		     *result = SCIP_CUTOFF;
1319 	
1320 	      return SCIP_OKAY;
1321 	   }
1322 	
1323 	   if( !initialized )
1324 	      return SCIP_OKAY;
1325 	
1326 	   ncols = SCIPmatrixGetNColumns(matrix);
1327 	
1328 	   /* sort columns by row indices */
1329 	   for( i = 0; i < ncols; i++ )
1330 	   {
1331 	      int* colpnt = SCIPmatrixGetColIdxPtr(matrix, i);
1332 	      SCIP_Real* valpnt = SCIPmatrixGetColValPtr(matrix, i);
1333 	      SCIPsortIntReal(colpnt, valpnt, SCIPmatrixGetColNNonzs(matrix, i));
1334 	   }
1335 	
1336 	   SCIP_CALL( SCIPallocBufferArray(scip, &scores, SCIPmatrixGetNRows(matrix)) );
1337 	   SCIP_CALL( SCIPallocBufferArray(scip, &perm, SCIPmatrixGetNRows(matrix)) );
1338 	   SCIP_CALL( SCIPallocBufferArray(scip, &ishashingcols, SCIPmatrixGetNColumns(matrix)) );
1339 	   SCIP_CALL( SCIPallocBufferArray(scip, &vars, SCIPmatrixGetNColumns(matrix)) );
1340 	   SCIP_CALL( SCIPallocBufferArray(scip, &isblockedvar, SCIPmatrixGetNColumns(matrix)) );
1341 	
1342 	   /* loop over all columns and create cons pairs */
1343 	   conspairssize = 0;
1344 	   nconspairs = 0;
1345 	   conspairs = NULL;
1346 	   SCIP_CALL( SCIPhashtableCreate(&pairtable, SCIPblkmem(scip), 1,
1347 	         SCIPhashGetKeyStandard, consPairsEqual, consPairHashval, (void*) scip) );
1348 	
1349 	   /* collect implied free variables and their number of nonzeros */
1350 	   for( c = 0; c < ncols; c++ )
1351 	   {
1352 	      SCIP_Bool lbimplied;
1353 	      SCIP_Bool ubimplied;
1354 	      int nnonz;
1355 	
1356 	      vars[c] = SCIPmatrixGetVar(matrix, c);
1357 	
1358 	      /* if the locks do not match do not consider the column for sparsification */
1359 	      if( SCIPmatrixDownlockConflict(matrix, c) || SCIPmatrixUplockConflict(matrix, c) )
1360 	      {
1361 	         isblockedvar[c] = TRUE;
1362 	         ishashingcols[c] = FALSE;
1363 	         continue;
1364 	      }
1365 	
1366 	      /* skip if the variable is not allowed to be multi-aggregated */
1367 	      if( SCIPdoNotMultaggrVar(scip, vars[c]) )
1368 	      {
1369 	         isblockedvar[c] = TRUE;
1370 	         ishashingcols[c] = FALSE;
1371 	         continue;
1372 	      }
1373 	
1374 	      nnonz = SCIPmatrixGetColNNonzs(matrix, c);
1375 	
1376 	      getImpliedBounds(scip, matrix, c, &ubimplied, &lbimplied);
1377 	
1378 	      ishashingcols[c] = FALSE;
1379 	
1380 	      if( lbimplied && ubimplied )
1381 	         ishashingcols[c] = TRUE;
1382 	
1383 	      isblockedvar[c] = FALSE;
1384 	
1385 	      /* only consider implied free variables
1386 	       * skip singleton variables, because either the constraint is redundant
1387 	       * or the variables can be canceled by variable substitution
1388 	       */
1389 	      if( nnonz >= 2 && (lbimplied && ubimplied) )
1390 	      {
1391 	         SCIP_Real* colvals;
1392 	         int* colinds;
1393 	         int failshift;
1394 	         int npairs;
1395 	
1396 	         colinds = SCIPmatrixGetColIdxPtr(matrix, c);
1397 	         colvals = SCIPmatrixGetColValPtr(matrix, c);
1398 	
1399 	         /* sort the rows non-decreasingly by number of nonzeros
1400 	          * if the number of nonzeros is equal, we use the colindex as tie-breaker
1401 	          */
1402 	         for( i = 0; i < nnonz; ++i )
1403 	         {
1404 	            perm[i] = i;
1405 	            scores[i] = -SCIPmatrixGetRowNNonzs(matrix, colinds[i]) - 1.0  *colinds[i] / ncols;
1406 	         }
1407 	         SCIPsortRealInt(scores, perm, nnonz);
1408 	
1409 	         if( presoldata->maxconsiderednonzeros >= 0 )
1410 	            nnonz = MIN(nnonz, presoldata->maxconsiderednonzeros);
1411 	
1412 	         npairs = (nnonz * (nnonz - 1)) / 2;
1413 	         if( nconspairs + npairs > conspairssize )
1414 	         {
1415 	            int newsize = SCIPcalcMemGrowSize(scip, nconspairs + npairs);
1416 	            SCIP_CALL( SCIPreallocBufferArray(scip, &conspairs, newsize) );
1417 	            conspairssize = newsize;
1418 	         }
1419 	
1420 	         /* if we are called after one or more failures, i.e., executions without finding cancellations, then we
1421 	          * shift the section of nonzeros considered; in the case that the maxconsiderednonzeros limit is hit, this
1422 	          * results in different constraint pairs being tried and avoids trying the same useless cancellations
1423 	          * repeatedly
1424 	          */
1425 	         failshift = presoldata->nfailures*presoldata->maxconsiderednonzeros;
1426 	
1427 	         for( i = 0; i < nnonz; ++i )
1428 	         {
1429 	            for( j = i + 1; j < nnonz; ++j )
1430 	            {
1431 	               int i1;
1432 	               int i2;
1433 	
1434 	               assert(nconspairs < conspairssize);
1435 	               assert(conspairs != NULL);
1436 	
1437 	               i1 = perm[(i + failshift) % nnonz];
1438 	               i2 = perm[(j + failshift) % nnonz];
1439 	               /* coverity[var_deref_op] */
1440 	               conspairs[nconspairs].colindex = c;
1441 	
1442 	               if( colinds[i1] < colinds[i2])
1443 	               {
1444 	                  conspairs[nconspairs].consindex1 = colinds[i1];
1445 	                  conspairs[nconspairs].consindex2 = colinds[i2];
1446 	                  conspairs[nconspairs].conscoef1 = colvals[i1];
1447 	                  conspairs[nconspairs].conscoef2 = colvals[i2];
1448 	               }
1449 	               else
1450 	               {
1451 	                  conspairs[nconspairs].consindex1 = colinds[i2];
1452 	                  conspairs[nconspairs].consindex2 = colinds[i1];
1453 	                  conspairs[nconspairs].conscoef1 = colvals[i2];
1454 	                  conspairs[nconspairs].conscoef2 = colvals[i1];
1455 	               }
1456 	               ++nconspairs;
1457 	            }
1458 	         }
1459 	      }
1460 	   }
1461 	
1462 	   /* insert conspairs into hash table */
1463 	   for( c = 0; c < nconspairs; ++c )
1464 	   {
1465 	      COLCONSPAIR* otherconspair;
1466 	      SCIP_Bool insert;
1467 	
1468 	      assert(conspairs != NULL);
1469 	
1470 	      insert = TRUE;
1471 	
1472 	      /* check if this pair is already contained in the hash table;
1473 	       * The loop is required due to the non-transitivity of the hash functions
1474 	       */
1475 	      while( (otherconspair = (COLCONSPAIR*)SCIPhashtableRetrieve(pairtable, (void*) &conspairs[c])) != NULL )
1476 	      {
1477 	         /* if the previous constraint pair has fewer or the same number of nonzeros in the attached column
1478 	          * we keep that pair and skip this one
1479 	          */
1480 	         if( SCIPmatrixGetColNNonzs(matrix, otherconspair->colindex) <=
1481 	            SCIPmatrixGetColNNonzs(matrix, conspairs[c].colindex) )
1482 	         {
1483 	            insert = FALSE;
1484 	            break;
1485 	         }
1486 	
1487 	         /* this pairs column has fewer nonzeros, so remove the other pair from the hash table and loop */
1488 	         SCIP_CALL( SCIPhashtableRemove(pairtable, (void*) otherconspair) );
1489 	      }
1490 	
1491 	      if( insert )
1492 	      {
1493 	         SCIP_CALL( SCIPhashtableInsert(pairtable, (void*) &conspairs[c]) );
1494 	      }
1495 	   }
1496 	
1497 	   /* sort cols according to decreasing sparsity */
1498 	   SCIP_CALL( SCIPallocBufferArray(scip, &colidxsorted, ncols) );
1499 	   SCIP_CALL( SCIPallocBufferArray(scip, &colsparsity, ncols) );
1500 	   for( c = 0; c < ncols; ++c )
1501 	   {
1502 	      colidxsorted[c] = c;
1503 	      colsparsity[c] = -SCIPmatrixGetColNNonzs(matrix, c);
1504 	   }
1505 	   SCIPsortIntInt(colsparsity, colidxsorted, ncols);
1506 	
1507 	   /* loop over the columns and cancel nonzeros until maximum number of retrieves is reached */
1508 	   maxuseless = (SCIP_Longint)(presoldata->maxretrievefac * (SCIP_Real)ncols);
1509 	   nuseless = 0;
1510 	   numcancel = 0;
1511 	   nfillin = 0;
1512 	   for( c = 0; c < ncols && nuseless <= maxuseless && !SCIPisStopped(scip); c++ )
1513 	   {
1514 	      int colidx;
1515 	
1516 	      colidx = colidxsorted[c];
1517 	
1518 	      if( isblockedvar[colidx] )
1519 	         continue;
1520 	
1521 	      /* since the function parameters for the max fillin are unsigned we do not need to handle the
1522 	       * unlimited (-1) case due to implicit conversion rules */
1523 	      SCIP_CALL( cancelCol(scip, matrix, presoldata, pairtable, ishashingcols, vars, isblockedvar, colidx, \
1524 	            presoldata->maxcontfillin == -1 ? INT_MAX : presoldata->maxcontfillin, \
1525 	            presoldata->maxintfillin == -1 ? INT_MAX : presoldata->maxintfillin, \
1526 	            presoldata->maxbinfillin == -1 ? INT_MAX : presoldata->maxbinfillin, \
1527 	            presoldata->maxconsiderednonzeros, presoldata->preserveintcoefs, \
1528 	            &nuseless, nchgcoefs, &numcancel, &nfillin, FALSE) );
1529 	   }
1530 	
1531 	   if( numcancel > 0 )
1532 	   {
1533 	      SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL,
1534 			      "   (%.1fs) dualsparsify: %d nonzeros canceled\n", SCIPgetSolvingTime(scip), numcancel);
1535 	      *result = SCIP_SUCCESS;
1536 	   }
1537 	   else /* do reductions on variables that contain larger nonzero entries */
1538 	   {
1539 	      SCIPhashtableRemoveAll(pairtable);
1540 	      nconspairs = 0;
1541 	
1542 	      /* collect large nonzero entries variables and their number of nonzeros */
1543 	      for( c = 0; c < ncols; c++ )
1544 	      {
1545 	         int nnonz;
1546 	
1547 	         nnonz = SCIPmatrixGetColNNonzs(matrix, c);
1548 	         vars[c] = SCIPmatrixGetVar(matrix, c);
1549 	
1550 	         /* if the locks do not match do not consider the column for sparsification */
1551 	         if( SCIPmatrixDownlockConflict(matrix, c) || SCIPmatrixUplockConflict(matrix, c) )
1552 	         {
1553 	            isblockedvar[c] = TRUE;
1554 	            ishashingcols[c] = FALSE;
1555 	            continue;
1556 	         }
1557 	
1558 	         isblockedvar[c] = FALSE;
1559 	
1560 	         /* only consider nonimplied free variables, i.e., non-hashing columns in the previous step,
1561 	          * with large nonzero entries
1562 	          * skip singleton variables, because either the constraint is redundant
1563 	          * or the variables can be canceled by variables substitution
1564 	          */
1565 	         if( nnonz >= presoldata->mineliminatednonzeros && !ishashingcols[c] )
1566 	         {
1567 	            int* colinds;
1568 	            SCIP_Real* colvals;
1569 	            int npairs;
1570 	            int failshift;
1571 	
1572 	            ishashingcols[c] = TRUE;
1573 	            colinds = SCIPmatrixGetColIdxPtr(matrix, c);
1574 	            colvals = SCIPmatrixGetColValPtr(matrix, c);
1575 	
1576 	            /* sort the rows non-decreasingly by number of nonzeros
1577 	             * if the number of nonzeros, we use the colindex as tie-breaker
1578 	             */
1579 	            for( i = 0; i < nnonz; ++i )
1580 	            {
1581 	               perm[i] = i;
1582 	               scores[i] = -SCIPmatrixGetRowNNonzs(matrix, colinds[i]) - 1.0 * colinds[i] / ncols;
1583 	            }
1584 	            SCIPsortRealInt(scores, perm, nnonz);
1585 	
1586 	            if( presoldata->maxconsiderednonzeros >= 0 )
1587 	               nnonz = MIN(nnonz, presoldata->maxconsiderednonzeros);
1588 	
1589 	            npairs = (nnonz * (nnonz - 1)) / 2;
1590 	            if( nconspairs + npairs > conspairssize )
1591 	            {
1592 	               int newsize = SCIPcalcMemGrowSize(scip, nconspairs + npairs);
1593 	               SCIP_CALL( SCIPreallocBufferArray(scip, &conspairs, newsize) );
1594 	               conspairssize = newsize;
1595 	            }
1596 	
1597 	            /* if we are called after one or more failures, i.e., executions without finding cancellations, then we
1598 	             * shift the section of nonzeros considered; in the case that the maxconsiderednonzeros limit is hit,
1599 	             * this results in different constraint pairs being tried and avoids trying the same useless
1600 	             * cancellations repeatedly
1601 	             */
1602 	            failshift = presoldata->nfailures*presoldata->maxconsiderednonzeros;
1603 	
1604 	            for( i = 0; i < nnonz; ++i )
1605 	            {
1606 	               for( j = i + 1; j < nnonz; ++j )
1607 	               {
1608 	                  int i1;
1609 	                  int i2;
1610 	
1611 	                  assert(nconspairs < conspairssize);
1612 	                  assert(conspairs != NULL);
1613 	
1614 	                  i1 = perm[(i + failshift) % nnonz];
1615 	                  i2 = perm[(j + failshift) % nnonz];
1616 	                  conspairs[nconspairs].colindex = c;
1617 	
1618 	                  if( colinds[i1] < colinds[i2])
1619 	                  {
1620 	                     conspairs[nconspairs].consindex1 = colinds[i1];
1621 	                     conspairs[nconspairs].consindex2 = colinds[i2];
1622 	                     conspairs[nconspairs].conscoef1 = colvals[i1];
1623 	                     conspairs[nconspairs].conscoef2 = colvals[i2];
1624 	                  }
1625 	                  else
1626 	                  {
1627 	                     conspairs[nconspairs].consindex1 = colinds[i2];
1628 	                     conspairs[nconspairs].consindex2 = colinds[i1];
1629 	                     conspairs[nconspairs].conscoef1 = colvals[i2];
1630 	                     conspairs[nconspairs].conscoef2 = colvals[i1];
1631 	                  }
1632 	                  ++nconspairs;
1633 	               }
1634 	            }
1635 	         }
1636 	         else
1637 	         {
1638 	            ishashingcols[c] = FALSE;
1639 	         }
1640 	      }
1641 	
1642 	      /* insert conspairs into hash table */
1643 	      for( c = 0; c < nconspairs; ++c )
1644 	      {
1645 	         SCIP_Bool insert;
1646 	         COLCONSPAIR* otherconspair;
1647 	
1648 	         assert(conspairs != NULL);
1649 	
1650 	         insert = TRUE;
1651 	
1652 	         /* check if this pair is already contained in the hash table;
1653 	          * The loop is required due to the non-transitivity of the hash functions
1654 	          */
1655 	         while( (otherconspair = (COLCONSPAIR*)SCIPhashtableRetrieve(pairtable, (void*) &conspairs[c])) != NULL )
1656 	         {
1657 	            /* if the previous constraint pair has fewer or the same number of nonzeros in the attached column
1658 	             * we keep that pair and skip this one
1659 	             */
1660 	            if( SCIPmatrixGetColNNonzs(matrix, otherconspair->colindex) <=
1661 	               SCIPmatrixGetColNNonzs(matrix, conspairs[c].colindex) )
1662 	            {
1663 	               insert = FALSE;
1664 	               break;
1665 	            }
1666 	
1667 	            /* this pairs column has fewer nonzeros, so remove the other pair from the hash table and loop */
1668 	            SCIP_CALL( SCIPhashtableRemove(pairtable, (void*) otherconspair) );
1669 	         }
1670 	
1671 	         if( insert )
1672 	         {
1673 	            SCIP_CALL( SCIPhashtableInsert(pairtable, (void*) &conspairs[c]) );
1674 	         }
1675 	      }
1676 	
1677 	      /* sort rows according to decreasingly sparsity */
1678 	      assert(colidxsorted != NULL);
1679 	      assert(colsparsity != NULL);
1680 	      for( c = 0; c < ncols; ++c )
1681 	      {
1682 	         colidxsorted[c] = c;
1683 	         colsparsity[c] = -SCIPmatrixGetColNNonzs(matrix, c);
1684 	      }
1685 	      SCIPsortIntInt(colsparsity, colidxsorted, ncols);
1686 	
1687 	      /* loop over the columns and cancel nonzeros until maximum number of retrieves is reached */
1688 	      maxuseless = (SCIP_Longint)(presoldata->maxretrievefac * (SCIP_Real)ncols);
1689 	      nuseless = 0;
1690 	      for( c = 0; c < ncols && nuseless <= maxuseless; c++ )
1691 	      {
1692 	         int colidx;
1693 	         int nnonz;
1694 	
1695 	         colidx = colidxsorted[c];
1696 	         nnonz = SCIPmatrixGetColNNonzs(matrix, colidx);
1697 	
1698 	         if( isblockedvar[colidx] || nnonz < presoldata->mineliminatednonzeros )
1699 	            continue;
1700 	
1701 	         /* since the function parameters for the max fillin are unsigned we do not need to handle the
1702 	          * unlimited (-1) case due to implicit conversion rules */
1703 	         SCIP_CALL( cancelCol(scip, matrix, presoldata, pairtable, ishashingcols, vars, isblockedvar, colidx, \
1704 	               presoldata->maxcontfillin == -1 ? INT_MAX : presoldata->maxcontfillin, \
1705 	               presoldata->maxintfillin == -1 ? INT_MAX : presoldata->maxintfillin, \
1706 	               presoldata->maxbinfillin == -1 ? INT_MAX : presoldata->maxbinfillin, \
1707 	               presoldata->maxconsiderednonzeros, presoldata->preserveintcoefs, \
1708 	               &nuseless, nchgcoefs, &numcancel, &nfillin, TRUE) );
1709 	      }
1710 	
1711 	      if( numcancel > 0 )
1712 	      {
1713 		 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL,
1714 				 "   (%.1fs) dualsparsify: %d nonzeros canceled\n", SCIPgetSolvingTime(scip), numcancel);
1715 	         *result = SCIP_SUCCESS;
1716 	      }
1717 	   }
1718 	
1719 	   updateFailureStatistic(presoldata, numcancel > 0);
1720 	
1721 	   SCIPfreeBufferArrayNull(scip, &conspairs);
1722 	   SCIPfreeBufferArray(scip, &colsparsity);
1723 	   SCIPfreeBufferArray(scip, &colidxsorted);
1724 	
1725 	   SCIPhashtableFree(&pairtable);
1726 	
1727 	   SCIPfreeBufferArray(scip, &isblockedvar);
1728 	   SCIPfreeBufferArray(scip, &vars);
1729 	   SCIPfreeBufferArray(scip, &ishashingcols);
1730 	   SCIPfreeBufferArray(scip, &perm);
1731 	   SCIPfreeBufferArray(scip, &scores);
1732 	
1733 	   SCIPmatrixFree(scip, &matrix);
1734 	
1735 	   return SCIP_OKAY;
1736 	}
1737 	
1738 	/*
1739 	 * presolver specific interface methods
1740 	 */
1741 	
1742 	/** destructor of presolver to free user data (called when SCIP is exiting) */
1743 	static
1744 	SCIP_DECL_PRESOLFREE(presolFreeDualsparsify)
1745 	{  /*lint --e{715}*/
1746 	   SCIP_PRESOLDATA* presoldata;
1747 	
1748 	   /* free presolver data */
1749 	   presoldata = SCIPpresolGetData(presol);
1750 	   assert(presoldata != NULL);
1751 	
1752 	   SCIPfreeBlockMemory(scip, &presoldata);
1753 	   SCIPpresolSetData(presol, NULL);
1754 	
1755 	   return SCIP_OKAY;
1756 	}
1757 	
1758 	/** initialization method of presolver (called after problem was transformed) */
1759 	static
1760 	SCIP_DECL_PRESOLINIT(presolInitDualsparsify)
1761 	{
1762 	   SCIP_PRESOLDATA* presoldata;
1763 	
1764 	   /* set the counters in the init (and not in the initpre) callback such that they persist across restarts */
1765 	   presoldata = SCIPpresolGetData(presol);
1766 	   presoldata->nfailures = 0;
1767 	   presoldata->nwaitingcalls = 0;
1768 	   presoldata->naggregated = 0;
1769 	
1770 	   return SCIP_OKAY;
1771 	}
1772 	
1773 	/** creates the dualsparsify presolver and includes it in SCIP */
1774 	SCIP_RETCODE SCIPincludePresolDualsparsify(
1775 	   SCIP*                 scip                /**< SCIP data structure */
1776 	   )
1777 	{
1778 	   SCIP_PRESOLDATA* presoldata;
1779 	   SCIP_PRESOL* presol;
1780 	
1781 	   /* create dualsparsify presolver data */
1782 	   SCIP_CALL( SCIPallocBlockMemory(scip, &presoldata) );
1783 	
1784 	   /* include presolver */
1785 	   SCIP_CALL( SCIPincludePresolBasic(scip, &presol, PRESOL_NAME, PRESOL_DESC, PRESOL_PRIORITY, PRESOL_MAXROUNDS,
1786 	         PRESOL_TIMING, presolExecDualsparsify, presoldata) );
1787 	
1788 	   SCIP_CALL( SCIPsetPresolCopy(scip, presol, presolCopyDualsparsify) );
1789 	   SCIP_CALL( SCIPsetPresolFree(scip, presol, presolFreeDualsparsify) );
1790 	   SCIP_CALL( SCIPsetPresolInit(scip, presol, presolInitDualsparsify) );
1791 	
1792 	   SCIP_CALL( SCIPaddBoolParam(scip,
1793 	         "presolving/dualsparsify/enablecopy",
1794 	         "should dualsparsify presolver be copied to sub-SCIPs?",
1795 	         &presoldata->enablecopy, TRUE, DEFAULT_ENABLECOPY, NULL, NULL) );
1796 	
1797 	   SCIP_CALL( SCIPaddBoolParam(scip,
1798 	         "presolving/dualsparsify/preserveintcoefs",
1799 	         "should we forbid cancellations that destroy integer coefficients?",
1800 	         &presoldata->preserveintcoefs, TRUE, DEFAULT_PRESERVEINTCOEFS, NULL, NULL) );
1801 	
1802 	   SCIP_CALL( SCIPaddBoolParam(scip,
1803 	         "presolving/dualsparsify/preservegoodlocks",
1804 	         "should we preserve good locked properties of variables (at most one lock in one direction)?",
1805 	         &presoldata->preservegoodlocks, TRUE, DEFAULT_PRESERVEGOODLOCKS, NULL, NULL) );
1806 	
1807 	   SCIP_CALL( SCIPaddIntParam(scip,
1808 	         "presolving/dualsparsify/maxcontfillin",
1809 	         "maximal fillin for continuous variables (-1: unlimited)",
1810 	         &presoldata->maxcontfillin, FALSE, DEFAULT_MAX_CONT_FILLIN, -1, INT_MAX, NULL, NULL) );
1811 	
1812 	   SCIP_CALL( SCIPaddIntParam(scip,
1813 	         "presolving/dualsparsify/maxbinfillin",
1814 	         "maximal fillin for binary variables (-1: unlimited)",
1815 	         &presoldata->maxbinfillin, FALSE, DEFAULT_MAX_BIN_FILLIN, -1, INT_MAX, NULL, NULL) );
1816 	
1817 	   SCIP_CALL( SCIPaddIntParam(scip,
1818 	         "presolving/dualsparsify/maxintfillin",
1819 	         "maximal fillin for integer variables including binaries (-1: unlimited)",
1820 	         &presoldata->maxintfillin, FALSE, DEFAULT_MAX_INT_FILLIN, -1, INT_MAX, NULL, NULL) );
1821 	
1822 	   SCIP_CALL( SCIPaddIntParam(scip,
1823 	         "presolving/dualsparsify/maxconsiderednonzeros",
1824 	         "maximal number of considered nonzeros within one column (-1: no limit)",
1825 	         &presoldata->maxconsiderednonzeros, TRUE, DEFAULT_MAXCONSIDEREDNONZEROS, -1, INT_MAX, NULL, NULL) );
1826 	
1827 	   SCIP_CALL( SCIPaddIntParam(scip,
1828 	         "presolving/dualsparsify/mineliminatednonzeros",
1829 	         "minimal eliminated nonzeros within one column if we need to add a constraint to the problem",
1830 	         &presoldata->mineliminatednonzeros, FALSE, DEFAULT_MINELIMINATEDNONZEROS, 0, INT_MAX, NULL, NULL) );
1831 	
1832 	   SCIP_CALL( SCIPaddRealParam(scip,
1833 	         "presolving/dualsparsify/maxretrievefac",
1834 	         "limit on the number of useless vs. useful hashtable retrieves as a multiple of the number of constraints",
1835 	         &presoldata->maxretrievefac, TRUE, DEFAULT_MAXRETRIEVEFAC, 0.0, SCIP_REAL_MAX, NULL, NULL) );
1836 	
1837 	   SCIP_CALL( SCIPaddRealParam(scip,
1838 	         "presolving/dualsparsify/waitingfac",
1839 	         "number of calls to wait until next execution as a multiple of the number of useless calls",
1840 	         &presoldata->waitingfac, TRUE, DEFAULT_WAITINGFAC, 0.0, SCIP_REAL_MAX, NULL, NULL) );
1841 	
1842 	   return SCIP_OKAY;
1843 	}
1844