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   	/**@file    presol_dualinfer.c
25   	 * @ingroup DEFPLUGINS_PRESOL
26   	 * @brief   dual inference presolver
27   	 * @author  Dieter Weninger
28   	 * @author  Patrick Gemander
29   	 *
30   	 * This presolver does bound strengthening on continuous variables (columns) for getting bounds on dual variables y.
31   	 * The bounds of the dual variables are then used to fix primal variables or change the side of constraints.
32   	 * For ranged rows one needs to decide which side (rhs or lhs) determines the equality.
33   	 *
34   	 * We distinguish two cases concerning complementary slackness:
35   	 * i)  reduced cost fixing:       c_j - sup_y(y^T A_{.j}) > 0 => x_j = l_j
36   	 *                                c_j - inf_y(y^T A_{.j}) < 0 => x_j = u_j
37   	 * ii) positive dual lower bound: y_i > 0 =>  A_{i.}x = b_i
38   	 *
39   	 * Further information on this presolving approach are given in
40   	 * Achterberg et al. "Presolve reductions in mixed integer programming"
41   	 * and for a two-column extension in
42   	 * Chen et al. "Two-row and two-column mixed-integer presolve using hasing-based pairing methods".
43   	 */
44   	
45   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
46   	
47   	#include <stdio.h>
48   	#include <assert.h>
49   	#include <string.h>
50   	
51   	#include "scip/scipdefplugins.h"
52   	#include "scip/pub_matrix.h"
53   	#include "blockmemshell/memory.h"
54   	#include "scip/cons_linear.h"
55   	#include "scip/presol_dualinfer.h"
56   	#include "scip/pub_cons.h"
57   	#include "scip/pub_matrix.h"
58   	#include "scip/pub_message.h"
59   	#include "scip/pub_presol.h"
60   	#include "scip/pub_var.h"
61   	#include "scip/scip_general.h"
62   	#include "scip/scip_mem.h"
63   	#include "scip/scip_message.h"
64   	#include "scip/scip_numerics.h"
65   	#include "scip/scip_presol.h"
66   	#include "scip/scip_prob.h"
67   	#include "scip/scip_probing.h"
68   	#include "scip/scip_var.h"
69   	
70   	#define PRESOL_NAME             "dualinfer"
71   	#define PRESOL_DESC             "exploit dual information for fixings and side changes"
72   	#define PRESOL_PRIORITY                (-3000) /**< priority of the presolver (>= 0: before, < 0: after constraint handlers) */
73   	#define PRESOL_MAXROUNDS                0    /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */
74   	#define PRESOL_TIMING                   SCIP_PRESOLTIMING_EXHAUSTIVE /* timing of the presolver (fast, medium, or exhaustive) */
75   	
76   	#define DEFAULT_TWOCOLUMN_COMBINE       TRUE /**< should two column convex combination be used per default */
77   	#define DEFAULT_MAXLOOPS_DUALBNDSTR     12   /**< default maximal number of loops for dual bound strengthening */
78   	#define DEFAULT_MAXCONSIDEREDNONZEROS   100  /**< default maximal number of considered non-zeros within one row */
79   	#define DEFAULT_MAXRETRIEVEFAILS        1000 /**< default maximal number of consecutive useless hashtable retrieves */
80   	#define DEFAULT_MAXCOMBINEFAILS         1000 /**< default maximal number of consecutive useless row combines */
81   	#define DEFAULT_MAXHASHFAC              10   /**< default maximal number of hashlist entries as multiple of number of rows in the problem */
82   	#define DEFAULT_MAXPAIRFAC              1    /**< default maximal number of processed row pairs as multiple of the number of rows in the problem */
83   	#define DEFAULT_MAXROWSUPPORT           3    /**< default maximal number of non-zeros in one row for turning an inequality into an equality */
84   	
85   	
86   	/*
87   	 * Data structures
88   	 */
89   	
90   	/** control parameters */
91   	struct SCIP_PresolData
92   	{
93   	   SCIP_Bool usetwocolcombine;               /**< use convex combination of two columns */
94   	   int maxdualbndloops;                      /**< default number of dual bound strengthening loops */
95   	   int maxpairfac;                           /**< maximal number of processed row pairs as multiple of the number of rows in the problem (-1: no limit) */
96   	   int maxhashfac;                           /**< maximal number of hashlist entries as multiple of number of rows in the problem (-1: no limit) */
97   	   int maxretrievefails;                     /**< maximal number of consecutive useless hashtable retrieves */
98   	   int maxcombinefails;                      /**< maximal number of consecutive useless row combines */
99   	   int maxconsiderednonzeros;                /**< maximal number of considered non-zeros within one row (-1: no limit) */
100  	   int maxrowsupport;                        /**< maximal number of non-zeros in one row for turning an inequality into an equality */
101  	};
102  	
103  	/** type of variable fixing direction */
104  	enum Fixingdirection
105  	{
106  	   FIXATLB = -1,                             /** fix variable at its lower bound */
107  	   NOFIX   =  0,                             /** no fixing */
108  	   FIXATUB =  1                              /** fix variable at its upper bound */
109  	};
110  	typedef enum Fixingdirection FIXINGDIRECTION;
111  	
112  	/** type of constraint side change */
113  	enum SideChange
114  	{
115  	   RHSTOLHS = -1,                            /** set rhs to value of lhs */
116  	   NOCHANGE =  0,                            /** no side change */
117  	   LHSTORHS =  1                             /** set lhs to value of rhs */
118  	};
119  	typedef enum SideChange SIDECHANGE;
120  	
121  	/** Signum for convex-combined variable coefficients \f$(\lambda * A_{ri} + (1 - \lambda) * A_{si})\f$
122  	 *  UP  - Coefficient changes from negative to positive for increasing lambda
123  	 *  DN  - Coefficient changes from positive to negative for increasing lambda
124  	 *  POS - Coefficient is positive for all lambda in (0,1)
125  	 *  NEG - Coefficient is negative for all lambda in (0,1)
126  	 */
127  	enum signum {UP, DN, POS, NEG};
128  	
129  	/** structure representing a pair of column indices; used for lookup in a hashtable */
130  	struct ColPair
131  	{
132  	   int col1idx;                              /**< first row index */
133  	   int col2idx;                              /**< second row index */
134  	};
135  	typedef struct ColPair COLPAIR;
136  	
137  	/*
138  	 * Local methods
139  	 */
140  	
141  	/** encode contents of a colpair as void* pointer */
142  	static
143  	void*
144  	encodeColPair(
145  	   COLPAIR*              colpair             /**< pointer to colpair */
146  	   )
147  	{
148  	   uint64_t a;
149  	   uint64_t b;
150  	
151  	   assert(colpair->col1idx >= 0);
152  	   assert(colpair->col2idx >= 0);
153  	
154  	   a = (uint64_t)(long)colpair->col1idx;
155  	   b = (uint64_t)(long)colpair->col2idx;
156  	   return (void*)((a << 32) | b);
157  	}
158  	
159  	/** compute single positive int hashvalue for two ints */
160  	static
161  	int
162  	hashIndexPair(
163  	   int                   idx1,               /**< first integer index */
164  	   int                   idx2                /**< second integer index */
165  	   )
166  	{
167  	   uint32_t hash = SCIPhashTwo(idx1, idx2);
168  	   return (int)(hash>>1);
169  	}
170  	
171  	/** add hash/rowidx pair to hashlist/rowidxlist **/
172  	static
173  	SCIP_RETCODE addEntry(
174  	   SCIP*                 scip,               /**< SCIP datastructure */
175  	   int*                  pos,                /**< position of last entry added */
176  	   int*                  listsize,           /**< size of hashlist and rowidxlist */
177  	   int**                 hashlist,           /**< block memory array containing hashes */
178  	   int**                 colidxlist,         /**< block memory array containing column indices */
179  	   int                   hash,               /**< hash to be inserted */
180  	   int                   colidx              /**< column index to be inserted */
181  	   )
182  	{
183  	   if( (*pos) >= (*listsize) )
184  	   {
185  	      int newsize  = SCIPcalcMemGrowSize(scip, (*pos) + 1);
186  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, hashlist, (*listsize), newsize) );
187  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, colidxlist, (*listsize), newsize) );
188  	      (*listsize) = newsize;
189  	   }
190  	
191  	   (*hashlist)[(*pos)] = hash;
192  	   (*colidxlist)[(*pos)] = colidx;
193  	   (*pos)++;
194  	
195  	   return SCIP_OKAY;
196  	}
197  	
198  	/** Within a sorted list, get next block with same value
199  	 *  E.g. for [h1, h1, h1, h2, h2, h2, h2, h3,...] and end = 0
200  	 *  returns start = 0, end = 3
201  	 *  and on a second call with end = 3 on the same list
202  	 *  returns start = 3, end = 7.
203  	 */
204  	static
205  	void findNextBlock(
206  	   const int*            list,               /**< list of integers */
207  	   int                   len,                /**< length of list */
208  	   int*                  start,              /**< variable to contain start index of found block */
209  	   int*                  end                 /**< variable to contain end index of found block */
210  	   )
211  	{
212  	   int i;
213  	   (*start) = (*end);
214  	   i = (*end) + 1;
215  	   while( i < len && list[i] == list[i - 1] )
216  	      i++;
217  	
218  	   (*end) = i;
219  	}
220  	
221  	/**
222  	 * The algorithm described in Belotti P. "Bound reduction using pairs of linear inequalities"
223  	 * tries to derive upper and lower bounds for all variables via convex combinations of linear inequalities
224  	 * We apply Belotti's algorithm to pairs of columns of continuous variables.
225  	 */
226  	static
227  	SCIP_RETCODE combineCols(
228  	   SCIP*                 scip,               /**< SCIP datastructure */
229  	   int*                  row1idxptr,         /**< indices specifying bound positions in lbs and ubs for first row */
230  	   int*                  row2idxptr,         /**< indices specifying bound positions in lbs und ubs for second row */
231  	   SCIP_Real*            row1valptr,         /**< first row coefficients */
232  	   SCIP_Real*            row2valptr,         /**< second row coefficients */
233  	   SCIP_Real             b1,                 /**< rhs of first row */
234  	   SCIP_Real             b2,                 /**< rhs of second row*/
235  	   int                   row1len,            /**< length of first row (e.g. row1idxptr and row1valptr)*/
236  	   int                   row2len,            /**< length of second row (e.g. row2idxptr and row2valptr)*/
237  	   int                   ncols,              /**< length of bound arrays lbs and ubs */
238  	   SCIP_Bool             swaprow1,           /**< should the sense of the first row be swapped to <= ? */
239  	   SCIP_Bool             swaprow2,           /**< should the sense of the second row be swapped to <= ? */
240  	   SCIP_Real*            lbs,                /**< lower bound array */
241  	   SCIP_Real*            ubs,                /**< upper bound array */
242  	   SCIP_Bool*            success             /**< we return (success || found better bounds") */
243  	   )
244  	{
245  	   int i;
246  	   int j;
247  	   int nvars;
248  	   int* varinds;
249  	   int nbreakpoints;
250  	   SCIP_Real* breakpoints;
251  	   int idx;
252  	   int idx1;
253  	   int idx2;
254  	   SCIP_Real* row1coefs;
255  	   SCIP_Real* row2coefs;
256  	   enum signum* signs;
257  	   int ninfs;
258  	   int l1infs;
259  	   SCIP_Real  l1;
260  	   SCIP_Real  l2;
261  	   SCIP_Real* newlbs;
262  	   SCIP_Real* newubs;
263  	   SCIP_Real coef;
264  	   int sign;
265  	   int shift;
266  	
267  	   SCIP_CALL( SCIPallocBufferArray(scip, &row1coefs, ncols) );
268  	   SCIP_CALL( SCIPallocBufferArray(scip, &row2coefs, ncols) );
269  	
270  	   SCIPsortIntReal(row1idxptr, row1valptr, row1len);
271  	   SCIPsortIntReal(row2idxptr, row2valptr, row2len);
272  	
273  	   /* swap rows if necessary */
274  	   if( swaprow1 )
275  	   {
276  	      for( i = 0; i < row1len; i++ )
277  	         row1coefs[row1idxptr[i]] = -row1valptr[i];
278  	      b1 = -b1;
279  	   }
280  	   else
281  	   {
282  	      for( i = 0; i < row1len; i++ )
283  	         row1coefs[row1idxptr[i]] = row1valptr[i];
284  	   }
285  	
286  	   if( swaprow2 )
287  	   {
288  	      for( i = 0; i < row2len; i++ )
289  	         row2coefs[row2idxptr[i]] = -row2valptr[i];
290  	      b2 = -b2;
291  	   }
292  	   else
293  	   {
294  	      for( i = 0; i < row2len; i++ )
295  	         row2coefs[row2idxptr[i]] = row2valptr[i];
296  	   }
297  	
298  	   SCIP_CALL( SCIPallocBufferArray(scip, &varinds, ncols) );
299  	   SCIP_CALL( SCIPallocBufferArray(scip, &signs, ncols) );
300  	   SCIP_CALL( SCIPallocBufferArray(scip, &breakpoints, ncols) );
301  	
302  	   /* calculate cancellation breakpoints and sign behaviour */
303  	   i = 0;
304  	   j = 0;
305  	   nvars = 0;
306  	   nbreakpoints = 0;
307  	   while( i < row1len && j < row2len )
308  	   {
309  	      assert(i + 1 == row1len || row1idxptr[i] < row1idxptr[i + 1]);
310  	      assert(j + 1 == row2len || row2idxptr[j] < row2idxptr[j + 1]);
311  	
312  	      idx1 = row1idxptr[i];
313  	      idx2 = row2idxptr[j];
314  	
315  	      /* We use 2.0 as default value for "no cancellation". For cancellations, this will be replaced by values in (0,1).
316  	       * A value larger than 1.0 is used because we sort the array and want to put non-cancellations to the end. */
317  	      breakpoints[nvars] = 2.0;
318  	
319  	      if( idx1 == idx2 )
320  	      {
321  	         if( (SCIPisNegative(scip, row1coefs[idx1]) && SCIPisPositive(scip, row2coefs[idx2])) ||
322  	            (SCIPisPositive(scip, row1coefs[idx1]) && SCIPisNegative(scip, row2coefs[idx2])) )
323  	         {
324  	            if( SCIPisNegative(scip, row2coefs[idx2]) )
325  	               signs[idx1] = UP;
326  	            else
327  	               signs[idx1] = DN;
328  	
329  	            breakpoints[nvars] = row2coefs[idx2] / (row2coefs[idx2] - row1coefs[idx1]);
330  	            nbreakpoints++;
331  	         }
332  	         else if( SCIPisPositive(scip, row1coefs[idx1]) )
333  	            signs[idx1] = POS;
334  	         else
335  	            signs[idx1] = NEG;
336  	
337  	         varinds[nvars] = idx1;
338  	         i++;
339  	         j++;
340  	      }
341  	      else if( idx1 < idx2 )
342  	      {
343  	         if( SCIPisPositive(scip, row1coefs[idx1]) )
344  	            signs[idx1] = POS;
345  	         else
346  	            signs[idx1] = NEG;
347  	
348  	         /* We will access this entry later on, so we explicitly write a zero here */
349  	         row2coefs[idx1] = 0.0;
350  	
351  	         varinds[nvars] = idx1;
352  	         i++;
353  	      }
354  	      else
355  	      {
356  	         assert(idx1 > idx2);
357  	         if( SCIPisPositive(scip, row2coefs[idx2]) )
358  	            signs[idx2] = POS;
359  	         else
360  	            signs[idx2] = NEG;
361  	
362  	         /* We will access this entry later on, so we explicitly write a zero here */
363  	         row1coefs[idx2] = 0.0;
364  	
365  	         varinds[nvars] = idx2;
366  	         j++;
367  	      }
368  	      nvars++;
369  	   }
370  	
371  	   while( i < row1len )
372  	   {
373  	      idx1 = row1idxptr[i];
374  	
375  	      if( SCIPisPositive(scip, row1coefs[idx1]) )
376  	         signs[idx1] = POS;
377  	      else
378  	         signs[idx1] = NEG;
379  	
380  	      /* We will access this entry later on, so we explicitly write a zero here */
381  	      row2coefs[idx1] = 0.0;
382  	
383  	      varinds[nvars] = idx1;
384  	      breakpoints[nvars] = 2.0;
385  	      nvars++;
386  	      i++;
387  	   }
388  	
389  	   while( j < row2len )
390  	   {
391  	      idx2 = row2idxptr[j];
392  	
393  	      if( SCIPisPositive(scip, row2coefs[idx2]) )
394  	         signs[idx2] = POS;
395  	      else
396  	         signs[idx2] = NEG;
397  	
398  	      /* We will access this entry later on, so we explicitly write a zero here */
399  	      row1coefs[idx2] = 0.0;
400  	
401  	      varinds[nvars] = idx2;
402  	      breakpoints[nvars] = 2.0;
403  	      nvars++;
404  	      j++;
405  	   }
406  	
407  	   SCIPsortRealInt(breakpoints, varinds, nvars);
408  	
409  	   /* The obvious preconditions for bound tightenings are met, so we try to calculate new bounds. */
410  	   if( nbreakpoints >= 1 )
411  	   {
412  	      SCIP_CALL( SCIPallocBufferArray(scip, &newlbs, nvars) );
413  	      SCIP_CALL( SCIPallocBufferArray(scip, &newubs, nvars) );
414  	
415  	      for( i = 0; i < nvars; i++)
416  	      {
417  	         idx = varinds[i];
418  	         newlbs[i] = lbs[idx];
419  	         newubs[i] = ubs[idx];
420  	      }
421  	
422  	      /* calculate activity contributions of each row */
423  	      l1 = b1;
424  	      l2 = b2;
425  	      l1infs = 0;
426  	      ninfs = 0;
427  	      for( i = 0; i < nvars; i++ )
428  	      {
429  	         idx = varinds[i];
430  	         if( !SCIPisZero(scip, row2coefs[idx]) )
431  	         {
432  	            if( SCIPisNegative(scip, row2coefs[idx]) )
433  	            {
434  	               if( !SCIPisInfinity(scip, -lbs[idx]) )
435  	               {
436  	                  l1 -= row1coefs[idx] * lbs[idx];
437  	                  l2 -= row2coefs[idx] * lbs[idx];
438  	               }
439  	               else
440  	                  ninfs++;
441  	            }
442  	            else
443  	            {
444  	               /* coefficient of second row is positive */
445  	               if( !SCIPisInfinity(scip, ubs[idx]) )
446  	               {
447  	                  l1 -= row1coefs[idx] * ubs[idx];
448  	                  l2 -= row2coefs[idx] * ubs[idx];
449  	               }
450  	               else
451  	                  ninfs++;
452  	            }
453  	         }
454  	         else
455  	         {
456  	            /* since row2coefs[idx] is zero, we have to choose the bound using row1coefs[idx] */
457  	            assert(!SCIPisZero(scip, row1coefs[idx]) && SCIPisZero(scip, row2coefs[idx]));
458  	            if( SCIPisNegative(scip, row1coefs[idx]) )
459  	            {
460  	               if( !SCIPisInfinity(scip, -lbs[idx]) )
461  	                  l1 -= row1coefs[idx] * lbs[idx];
462  	               else
463  	                  l1infs++;
464  	            }
465  	            else
466  	            {
467  	               /* coefficient of first row is positive */
468  	               if( !SCIPisInfinity(scip, ubs[idx]) )
469  	                  l1 -= row1coefs[idx] * ubs[idx];
470  	               else
471  	                  l1infs++;
472  	            }
473  	         }
474  	      }
475  	
476  	      /* Calculate bounds for lambda = 0 */
477  	#ifdef SCIP_MORE_DEBUG
478  	      SCIPdebugMsg(scip, "lambda = 0, l1 = %g, l2 = %g, ninfs = %d\n", i, breakpoints[i], l1, l2, ninfs);
479  	#endif
480  	
481  	      if( ninfs <= 1 )
482  	      {
483  	#ifdef SCIP_MORE_DEBUG
484  	         SCIP_Real oldlb;
485  	         SCIP_Real oldub;
486  	#endif
487  	         for( i = 0; i < nvars; i++ )
488  	         {
489  	#ifdef SCIP_MORE_DEBUG
490  	            oldlb = newlbs[i];
491  	            oldub = newubs[i];
492  	#endif
493  	            idx = varinds[i];
494  	            if( SCIPisPositive(scip, row2coefs[idx]) )
495  	            {
496  	               if( ninfs == 0 )
497  	                  newlbs[i] = MAX(newlbs[i], (l2 + row2coefs[idx] * ubs[idx]) / row2coefs[idx]);
498  	               else if( SCIPisInfinity(scip, ubs[idx]) )
499  	                  newlbs[i] = MAX(newlbs[i], l2 / row2coefs[idx]);
500  	            }
501  	            else if ( SCIPisNegative(scip, row2coefs[idx]) )
502  	            {
503  	               if( ninfs == 0 )
504  	                  newubs[i] = MIN(newubs[i], (l2 + row2coefs[idx] * lbs[idx]) / row2coefs[idx]);
505  	               else if( SCIPisInfinity(scip, -lbs[idx]) )
506  	                  newubs[i] = MIN(newubs[i], l2 / row2coefs[idx]);
507  	            }
508  	#ifdef SCIP_MORE_DEBUG
509  	            if( !SCIPisEQ(scip, oldlb, newlbs[i]) || !SCIPisEQ(scip, oldub, newubs[i]) )
510  	               SCIPdebugMsg(scip, "%g <= %g <= var_%d <= %g <= %g\n", oldlb, newlbs[i], i, newubs[i], oldub);
511  	#endif
512  	         }
513  	      }
514  	
515  	      ninfs += l1infs;
516  	
517  	      i = 0;
518  	      while( i < nbreakpoints )
519  	      {
520  	         int nnewinfs;
521  	         SCIP_Real l1update;
522  	         SCIP_Real l2update;
523  	         SCIP_Bool updated;
524  	
525  	         /* determine number of infinities and compute update for l1 and l2 */
526  	         shift = 0;
527  	         nnewinfs = 0;
528  	         l1update = 0.0;
529  	         l2update = 0.0;
530  	         updated = FALSE;
531  	         j = i;
532  	         while( !updated )
533  	         {
534  	            idx = varinds[j];
535  	            assert(signs[idx] == UP || signs[idx] == DN);
536  	            if( signs[idx] == UP )
537  	               sign = 1;
538  	            else
539  	               sign = -1;
540  	
541  	            if( !SCIPisInfinity(scip, -lbs[idx]) )
542  	            {
543  	               l1update += sign * row1coefs[idx] * lbs[idx];
544  	               l2update += sign * row2coefs[idx] * lbs[idx];
545  	            }
546  	            else
547  	            {
548  	               if( signs[idx] == UP  )
549  	                  ninfs--;
550  	               else
551  	                  nnewinfs++;
552  	            }
553  	
554  	            if( !SCIPisInfinity(scip, ubs[idx]) )
555  	            {
556  	               l1update -= sign * row1coefs[idx] * ubs[idx];
557  	               l2update -= sign * row2coefs[idx] * ubs[idx];
558  	            }
559  	            else
560  	            {
561  	               if( signs[idx] == UP  )
562  	                  nnewinfs++;
563  	               else
564  	                  ninfs--;
565  	            }
566  	
567  	            if( signs[idx] == UP )
568  	               signs[idx] = POS;
569  	            else
570  	               signs[idx] = NEG;
571  	
572  	            if( j + 1 >= nbreakpoints || !SCIPisEQ(scip, breakpoints[j], breakpoints[j + 1]) )
573  	               updated = TRUE;
574  	
575  	            shift++;
576  	            j++;
577  	         }
578  	
579  	#ifdef SCIP_MORE_DEBUG
580  	         SCIPdebugMsg(scip, "lambda_%d = %g, l1 = %g, l2 = %g, ninfs = %d\n", i, breakpoints[i], l1, l2, ninfs);
581  	#endif
582  	
583  	         assert(ninfs >= 0);
584  	
585  	         /* if more than one infinity destroys our bounds we cannot tighten anything */
586  	         if( ninfs <= 1 )
587  	         {
588  	            /* check for bounds to be tightened */
589  	            for( j = 0; j < nvars; j++ )
590  	            {
591  	#ifdef SCIP_MORE_DEBUG
592  	               SCIP_Real oldlb;
593  	               SCIP_Real oldub;
594  	#endif
595  	
596  	               /* catch the special case where the entire remaining constraint is cancelled */
597  	               if( j >= nvars )
598  	                  break;
599  	
600  	#ifdef SCIP_MORE_DEBUG
601  	               oldlb = newlbs[j];
602  	               oldub = newubs[j];
603  	#endif
604  	
605  	               idx = varinds[j];
606  	               coef = breakpoints[i] * row1coefs[idx] + (1 - breakpoints[i]) * row2coefs[idx];
607  	               assert(!SCIPisEQ(scip, breakpoints[i], 2.0));
608  	
609  	               /* skip if the coefficient is too close to zero as it becomes numerically unstable */
610  	               if( SCIPisZero(scip, coef) )
611  	                  continue;
612  	
613  	               if( signs[idx] == POS || signs[idx] == DN )
614  	               {
615  	                  if( ninfs == 0 )
616  	                     newlbs[j] = MAX(newlbs[j], (breakpoints[i] * l1 + (1 - breakpoints[i]) * l2 + coef * ubs[idx]) / coef);
617  	                  else if( SCIPisInfinity(scip, ubs[idx]) )
618  	                     newlbs[j] = MAX(newlbs[j], (breakpoints[i] * l1 + (1 - breakpoints[i]) * l2) / coef);
619  	               }
620  	               else if ( signs[idx] == NEG || signs[idx] == UP )
621  	               {
622  	                  if( ninfs == 0 )
623  	                     newubs[j] = MIN(newubs[j], (breakpoints[i] * l1 + (1 - breakpoints[i]) * l2 + coef * lbs[idx]) / coef);
624  	                  else if( SCIPisInfinity(scip, -lbs[idx]) )
625  	                     newubs[j] = MIN(newubs[j], (breakpoints[i] * l1 + (1 - breakpoints[i]) * l2) / coef);
626  	               }
627  	#ifdef SCIP_MORE_DEBUG
628  	               if( !SCIPisEQ(scip, oldlb, newlbs[j]) || !SCIPisEQ(scip, oldub, newubs[j]) )
629  	                  SCIPdebugMsg(scip, "%g <= %g <= var_%d <= %g <= %g\n", oldlb, newlbs[j], j, newubs[j], oldub);
630  	#endif
631  	            }
632  	         }
633  	
634  	         i += shift;
635  	         ninfs += nnewinfs;
636  	         l1 += l1update;
637  	         l2 += l2update;
638  	      }
639  	
640  	      /* check infinities in first row */
641  	      ninfs = 0;
642  	      for( i = 0; i < nvars; i++ )
643  	      {
644  	         idx = varinds[i];
645  	         if( (SCIPisPositive(scip, row1coefs[idx]) && SCIPisInfinity(scip, ubs[idx]))
646  	            || (SCIPisNegative(scip, row1coefs[idx]) && SCIPisInfinity(scip, -lbs[idx])) )
647  	            ninfs++;
648  	      }
649  	
650  	      /* calculate bounds for lambda = 1 */
651  	#ifdef SCIP_MORE_DEBUG
652  	      SCIPdebugMsg(scip, "lambda = 1, l1 = %g, l2 = %g, ninfs = %d\n", i, breakpoints[i], l1, l2, ninfs);
653  	#endif
654  	      if( ninfs <= 1 )
655  	      {
656  	#ifdef SCIP_MORE_DEBUG
657  	         SCIP_Real oldlb;
658  	         SCIP_Real oldub;
659  	#endif
660  	         for( i = 0; i < nvars; i++ )
661  	         {
662  	#ifdef SCIP_MORE_DEBUG
663  	            oldlb = newlbs[i];
664  	            oldub = newubs[i];
665  	#endif
666  	            idx = varinds[i];
667  	            if( SCIPisPositive(scip, row1coefs[idx]) )
668  	            {
669  	               if( ninfs == 0 )
670  	                  newlbs[i] = MAX(newlbs[i], (l1 + row1coefs[idx] * ubs[idx]) / row1coefs[idx]);
671  	               else if( SCIPisInfinity(scip, ubs[idx]) )
672  	                  newlbs[i] = MAX(newlbs[i], l1 / row1coefs[idx]);
673  	            }
674  	            else if ( SCIPisNegative(scip, row1coefs[idx]) )
675  	            {
676  	               if( ninfs == 0 )
677  	                  newubs[i] = MIN(newubs[i], (l1 + row1coefs[idx] * lbs[idx]) / row1coefs[idx]);
678  	               else if( SCIPisInfinity(scip, -lbs[idx]) )
679  	                  newubs[i] = MIN(newubs[i], l1 / row1coefs[idx]);
680  	            }
681  	#ifdef SCIP_MORE_DEBUG
682  	            if( !SCIPisEQ(scip, oldlb, newlbs[i]) || !SCIPisEQ(scip, oldub, newubs[i]) )
683  	               SCIPdebugMsg(scip, "%g <= %g <= var_%i <= %g <= %g\n", oldlb, newlbs[i], i, newubs[i], oldub);
684  	#endif
685  	         }
686  	      }
687  	
688  	      /* update bound arrays and determine success */
689  	      for( i = 0; i < nvars; i++ )
690  	      {
691  	         idx = varinds[i];
692  	
693  	         assert(SCIPisLE(scip, lbs[idx], newlbs[i]));
694  	         assert(SCIPisGE(scip, ubs[idx], newubs[i]));
695  	
696  	         if( SCIPisGT(scip, newlbs[i], lbs[idx]) || SCIPisLT(scip, newubs[i], ubs[idx]) )
697  	         {
698  	            (*success) = TRUE;
699  	
700  	            lbs[idx] = newlbs[i];
701  	            ubs[idx] = newubs[i];
702  	         }
703  	      }
704  	      SCIPfreeBufferArray(scip, &newubs);
705  	      SCIPfreeBufferArray(scip, &newlbs);
706  	   }
707  	
708  	   SCIPfreeBufferArray(scip, &breakpoints);
709  	   SCIPfreeBufferArray(scip, &signs);
710  	   SCIPfreeBufferArray(scip, &varinds);
711  	   SCIPfreeBufferArray(scip, &row2coefs);
712  	   SCIPfreeBufferArray(scip, &row1coefs);
713  	
714  	   return SCIP_OKAY;
715  	}
716  	
717  	/** get minimal and maximal residual activities without one specific column */
718  	static
719  	void getMinMaxActivityResiduals(
720  	   SCIP*                 scip,               /**< SCIP main data structure */
721  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
722  	   int                   withoutcol,         /**< exclude this column index */
723  	   int                   row,                /**< row index */
724  	   SCIP_Real*            lbs,                /**< lower bounds */
725  	   SCIP_Real*            ubs,                /**< upper bounds */
726  	   SCIP_Real*            minresactivity,     /**< minimum residual activity of this row */
727  	   SCIP_Real*            maxresactivity,     /**< maximum residual activity of this row */
728  	   SCIP_Bool*            isminsettoinfinity, /**< flag indicating if minresactiviy is set to infinity */
729  	   SCIP_Bool*            ismaxsettoinfinity  /**< flag indicating if maxresactiviy is set to infinity */
730  	   )
731  	{
732  	   SCIP_Real coef;
733  	   int* rowpnt;
734  	   int* rowend;
735  	   SCIP_Real* valpnt;
736  	   int nmaxactneginf;
737  	   int nmaxactposinf;
738  	   int nminactneginf;
739  	   int nminactposinf;
740  	   SCIP_Real maxresact;
741  	   SCIP_Real minresact;
742  	   int col;
743  	
744  	   assert(scip != NULL);
745  	   assert(matrix != NULL);
746  	   assert(minresactivity != NULL);
747  	   assert(maxresactivity != NULL);
748  	   assert(isminsettoinfinity != NULL);
749  	   assert(ismaxsettoinfinity != NULL);
750  	
751  	   *isminsettoinfinity = FALSE;
752  	   *ismaxsettoinfinity = FALSE;
753  	
754  	   nmaxactneginf = 0;
755  	   nmaxactposinf = 0;
756  	   nminactneginf = 0;
757  	   nminactposinf = 0;
758  	   maxresact = 0;
759  	   minresact = 0;
760  	
761  	   rowpnt = SCIPmatrixGetRowIdxPtr(matrix, row);
762  	   rowend = rowpnt + SCIPmatrixGetRowNNonzs(matrix, row);
763  	   valpnt = SCIPmatrixGetRowValPtr(matrix, row);
764  	
765  	   for( ; rowpnt < rowend; rowpnt++, valpnt++ )
766  	   {
767  	      col = *rowpnt;
768  	
769  	      if( col == withoutcol )
770  	         continue;
771  	
772  	      coef = *valpnt;
773  	
774  	      /* positive coefficient */
775  	      if( coef > 0.0 )
776  	      {
777  	         if( SCIPisInfinity(scip, ubs[col]) )
778  	            nmaxactposinf++;
779  	         else
780  	            maxresact += coef * ubs[col];
781  	
782  	         if( SCIPisInfinity(scip, -lbs[col]) )
783  	            nminactneginf++;
784  	         else
785  	            minresact += coef * lbs[col];
786  	      }
787  	      else /* negative coefficient */
788  	      {
789  	         if( SCIPisInfinity(scip, -lbs[col]) )
790  	            nmaxactneginf++;
791  	         else
792  	            maxresact += coef * lbs[col];
793  	
794  	         if( SCIPisInfinity(scip, ubs[col]) )
795  	            nminactposinf++;
796  	         else
797  	            minresact += coef * ubs[col];
798  	      }
799  	   }
800  	
801  	   if( (nmaxactneginf + nmaxactposinf) > 0 )
802  	      *ismaxsettoinfinity = TRUE;
803  	   else
804  	      *maxresactivity = maxresact;
805  	
806  	   if( (nminactneginf + nminactposinf) > 0 )
807  	      *isminsettoinfinity = TRUE;
808  	   else
809  	      *minresactivity = minresact;
810  	}
811  	
812  	/** calculate the upper and lower bound of one variable from one row */
813  	static
814  	void getVarBoundsOfRow(
815  	   SCIP*                 scip,               /**< SCIP main data structure */
816  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
817  	   int                   col,                /**< column index of variable */
818  	   int                   row,                /**< row index */
819  	   SCIP_Real             val,                /**< coefficient of this column in this row */
820  	   SCIP_Real*            lbs,                /**< lower bounds */
821  	   SCIP_Real*            ubs,                /**< upper bounds */
822  	   SCIP_Real*            rowub,              /**< upper bound of row */
823  	   SCIP_Bool*            ubfound,            /**< flag indicating that an upper bound was calculated */
824  	   SCIP_Real*            rowlb,              /**< lower bound of row */
825  	   SCIP_Bool*            lbfound             /**< flag indicating that a lower bound was caluclated */
826  	   )
827  	{
828  	   SCIP_Bool isminsettoinfinity;
829  	   SCIP_Bool ismaxsettoinfinity;
830  	   SCIP_Real minresactivity;
831  	   SCIP_Real maxresactivity;
832  	   SCIP_Real lhs;
833  	   SCIP_Real rhs;
834  	
835  	   assert(rowub != NULL);
836  	   assert(ubfound != NULL);
837  	   assert(rowlb != NULL);
838  	   assert(lbfound != NULL);
839  	
840  	   *rowub = SCIPinfinity(scip);
841  	   *ubfound = FALSE;
842  	   *rowlb = -SCIPinfinity(scip);
843  	   *lbfound = FALSE;
844  	
845  	   getMinMaxActivityResiduals(scip, matrix, col, row, lbs, ubs,
846  	      &minresactivity, &maxresactivity,
847  	      &isminsettoinfinity, &ismaxsettoinfinity);
848  	
849  	   lhs = SCIPmatrixGetRowLhs(matrix, row);
850  	   rhs = SCIPmatrixGetRowRhs(matrix, row);
851  	
852  	   if( val > 0.0 )
853  	   {
854  	      if( !isminsettoinfinity && !SCIPisInfinity(scip, rhs) )
855  	      {
856  	         *rowub = (rhs - minresactivity) / val; // maybe one wants some kind of numerical guard of check that values is not too small for all these
857  	         *ubfound = TRUE;
858  	      }
859  	
860  	      if( !ismaxsettoinfinity && !SCIPisInfinity(scip, -lhs) )
861  	      {
862  	         *rowlb = (lhs - maxresactivity) / val;
863  	         *lbfound = TRUE;
864  	      }
865  	   }
866  	   else
867  	   {
868  	      if( !ismaxsettoinfinity && !SCIPisInfinity(scip, -lhs) )
869  	      {
870  	         *rowub = (lhs - maxresactivity) / val;
871  	         *ubfound = TRUE;
872  	      }
873  	
874  	      if( !isminsettoinfinity && !SCIPisInfinity(scip, rhs) )
875  	      {
876  	         *rowlb = (rhs - minresactivity) / val;
877  	         *lbfound = TRUE;
878  	      }
879  	   }
880  	}
881  	
882  	
883  	/** detect implied variable bounds */
884  	static
885  	void getImpliedBounds(
886  	   SCIP*                 scip,               /**< SCIP main data structure */
887  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
888  	   int                   col,                /**< column index for implied free test */
889  	   SCIP_Real*            lbs,                /**< lower bounds */
890  	   SCIP_Real*            ubs,                /**< upper bounds */
891  	   SCIP_Bool*            ubimplied,          /**< flag indicating an implied upper bound */
892  	   SCIP_Bool*            lbimplied           /**< flag indicating an implied lower bound */
893  	   )
894  	{
895  	   SCIP_Real impliedub;
896  	   SCIP_Real impliedlb;
897  	   int* colpnt;
898  	   int* colend;
899  	   SCIP_Real* valpnt;
900  	
901  	   assert(scip != NULL);
902  	   assert(matrix != NULL);
903  	   assert(lbs != NULL);
904  	   assert(ubs != NULL);
905  	   assert(ubimplied != NULL);
906  	   assert(lbimplied != NULL);
907  	
908  	   *ubimplied = FALSE;
909  	   impliedub = SCIPinfinity(scip);
910  	
911  	   *lbimplied = FALSE;
912  	   impliedlb = -SCIPinfinity(scip);
913  	
914  	   colpnt = SCIPmatrixGetColIdxPtr(matrix, col);
915  	   colend = colpnt + SCIPmatrixGetColNNonzs(matrix, col);
916  	   valpnt = SCIPmatrixGetColValPtr(matrix, col);
917  	   for( ; (colpnt < colend); colpnt++, valpnt++ )
918  	   {
919  	      SCIP_Real rowub;
920  	      SCIP_Bool ubfound;
921  	      SCIP_Real rowlb;
922  	      SCIP_Bool lbfound;
923  	
924  	      getVarBoundsOfRow(scip, matrix, col, *colpnt, *valpnt, lbs, ubs,
925  	         &rowub, &ubfound, &rowlb, &lbfound);
926  	
927  	      if( ubfound && (rowub < impliedub) )
928  	         impliedub = rowub;
929  	
930  	      if( lbfound && (rowlb > impliedlb) )
931  	         impliedlb = rowlb;
932  	   }
933  	
934  	   /* we consider +/-inf bounds as implied bounds */
935  	   if( SCIPisInfinity(scip, ubs[col]) ||
936  	      (!SCIPisInfinity(scip, ubs[col]) && SCIPisLE(scip, impliedub, ubs[col])) )
937  	      *ubimplied = TRUE;
938  	
939  	   if( SCIPisInfinity(scip, -lbs[col]) ||
940  	      (!SCIPisInfinity(scip, -lbs[col]) && SCIPisGE(scip, impliedlb, lbs[col])) )
941  	      *lbimplied = TRUE;
942  	}
943  	
944  	
945  	/** calculate minimal column activity from one variable without one row */
946  	static
947  	SCIP_Real getMinColActWithoutRow(
948  	   SCIP*                 scip,               /**< SCIP main data structure */
949  	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
950  	   int                   col,                /**< column index */
951  	   int                   withoutrow,         /**< exclude this row index */
952  	   SCIP_Real*            lbdual,             /**< lower bounds of dual variables */
953  	   SCIP_Real*            ubdual              /**< upper bounds of dual variables */
954  	   )
955  	{
956  	   SCIP_Real* valpnt;
957  	   int* colpnt;
958  	   int* colend;
959  	   SCIP_Real val;
960  	   SCIP_Real mincolactivity;
961  	   int row;
962  	
963  	   assert(scip != NULL);
964  	   assert(matrix != NULL);
965  	   assert(lbdual != NULL);
966  	   assert(ubdual != NULL);
967  	
968  	   mincolactivity = 0;
969  	
970  	   colpnt = SCIPmatrixGetColIdxPtr(matrix, col);
971  	   colend = colpnt + SCIPmatrixGetColNNonzs(matrix, col);
972  	   valpnt = SCIPmatrixGetColValPtr(matrix, col);
973  	
974  	   for( ; colpnt < colend; colpnt++, valpnt++ )
975  	   {
976  	      row = *colpnt;
977  	      val = *valpnt;
978  	
979  	      if( row == withoutrow )
980  	         continue;
981  	
982  	      if( val > 0.0 )
983  	      {
984  	         assert(!SCIPisInfinity(scip, -lbdual[row]));
985  	         mincolactivity += val * lbdual[row];
986  	      }
987  	      else if( val < 0.0 )
988  	      {
989  	         assert(!SCIPisInfinity(scip, ubdual[row]));
990  	         mincolactivity += val * ubdual[row];
991  	      }
992  	   }
993  	
994  	   return mincolactivity;
995  	}
996  	
997  	
998  	/** In the primal the residual activity of a constraint w.r.t. a variable is the activity of the constraint without the variable.
999  	 *  This function does the same but in the dual.
1000 	 *  It computes the residual activity of column 'col' w.r.t. variable 'row'
1001 	 */
1002 	static
1003 	void calcMinColActResidual(
1004 	   SCIP*                 scip,               /**< SCIP main data structure */
1005 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1006 	   int                   col,                /**< column index */
1007 	   int                   row,                /**< row index */
1008 	   SCIP_Real             val,                /**< matrix coefficient */
1009 	   SCIP_Real*            lbdual,             /**< lower bounds of the dual variables */
1010 	   SCIP_Real*            ubdual,             /**< upper bounds of the dual variables */
1011 	   const SCIP_Real*      mincolact,          /**< minimal column activities */
1012 	   const int*            mincolactinf,       /**< number of infinite contributions to minimal column activity */
1013 	   SCIP_Real*            mincolresact        /**< minimal residual column activity */
1014 	   )
1015 	{
1016 	   assert(scip != NULL);
1017 	   assert(matrix != NULL);
1018 	   assert(lbdual != NULL);
1019 	   assert(ubdual != NULL);
1020 	   assert(mincolact != NULL);
1021 	   assert(mincolactinf != NULL);
1022 	   assert(mincolresact != NULL);
1023 	
1024 	   *mincolresact = -SCIPinfinity(scip);
1025 	
1026 	   if( val > 0.0 )
1027 	   {
1028 	      if( SCIPisInfinity(scip, -lbdual[row]) )
1029 	      {
1030 	         assert(mincolactinf[col] >= 1);
1031 	         if( mincolactinf[col] == 1 )
1032 	            *mincolresact = getMinColActWithoutRow(scip, matrix, col, row, lbdual, ubdual);
1033 	         else
1034 	            *mincolresact = -SCIPinfinity(scip);
1035 	      }
1036 	      else
1037 	      {
1038 	         if( mincolactinf[col] > 0 )
1039 	            *mincolresact = -SCIPinfinity(scip);
1040 	         else
1041 	            *mincolresact = mincolact[col] - val * lbdual[row];
1042 	      }
1043 	   }
1044 	   else if( val < 0.0 )
1045 	   {
1046 	      if( SCIPisInfinity(scip, ubdual[row]) )
1047 	      {
1048 	         assert(mincolactinf[col] >= 1);
1049 	         if( mincolactinf[col] == 1 )
1050 	            *mincolresact = getMinColActWithoutRow(scip, matrix, col, row, lbdual, ubdual);
1051 	         else
1052 	            *mincolresact = -SCIPinfinity(scip);
1053 	      }
1054 	      else
1055 	      {
1056 	         if( mincolactinf[col] > 0 )
1057 	            *mincolresact = -SCIPinfinity(scip);
1058 	         else
1059 	            *mincolresact = mincolact[col] - val * ubdual[row];
1060 	      }
1061 	   }
1062 	}
1063 	
1064 	/** calculate minimal column activity of one column */
1065 	static
1066 	void calcMinColActivity(
1067 	   SCIP*                 scip,               /**< SCIP main data structure */
1068 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1069 	   int                   col,                /**< column for activity calculations */
1070 	   SCIP_Real*            lbdual,             /**< lower bounds of dual variables */
1071 	   SCIP_Real*            ubdual,             /**< upper bounds of dual variables */
1072 	   SCIP_Real*            mincolact,          /**< minimal column activities */
1073 	   int*                  mincolactinf        /**< number of -inf contributions to minimal column activity */
1074 	   )
1075 	{
1076 	   SCIP_Real* valpnt;
1077 	   int* colpnt;
1078 	   int* colend;
1079 	   SCIP_Real val;
1080 	   int row;
1081 	
1082 	   assert(scip != NULL);
1083 	   assert(matrix != NULL);
1084 	   assert(lbdual != NULL);
1085 	   assert(ubdual != NULL);
1086 	   assert(mincolact != NULL);
1087 	   assert(mincolactinf != NULL);
1088 	
1089 	   /* init activities */
1090 	   mincolact[col] = 0.0;
1091 	   mincolactinf[col] = 0;
1092 	
1093 	   colpnt = SCIPmatrixGetColIdxPtr(matrix, col);
1094 	   colend = colpnt + SCIPmatrixGetColNNonzs(matrix, col);
1095 	   valpnt = SCIPmatrixGetColValPtr(matrix, col);
1096 	
1097 	   /* calculate column activities */
1098 	   for( ; colpnt < colend; colpnt++, valpnt++ )
1099 	   {
1100 	      row = *colpnt;
1101 	      val = *valpnt;
1102 	
1103 	      if( val > 0.0 )
1104 	      {
1105 	         if(SCIPisInfinity(scip, -lbdual[row]))
1106 	            mincolactinf[col]++;
1107 	         else
1108 	            mincolact[col] += val * lbdual[row];
1109 	      }
1110 	      else if( val < 0.0 )
1111 	      {
1112 	         if(SCIPisInfinity(scip, ubdual[row]))
1113 	            mincolactinf[col]++;
1114 	         else
1115 	            mincolact[col] += val * ubdual[row];
1116 	      }
1117 	   }
1118 	
1119 	   /* update column activities if infinity counters are greater 0 */
1120 	   if( mincolactinf[col] > 0 )
1121 	      mincolact[col] = -SCIPinfinity(scip);
1122 	}
1123 	
1124 	/** calculate maximal column activity of one column */
1125 	static
1126 	void calcMaxColActivity(
1127 	   SCIP*                 scip,               /**< SCIP main data structure */
1128 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1129 	   int                   col,                /**< column for activity calculations */
1130 	   SCIP_Real*            lbdual,             /**< lower bounds of dual variables */
1131 	   SCIP_Real*            ubdual,             /**< upper bounds of dual variables */
1132 	   SCIP_Real*            maxcolact,          /**< minimal column activities */
1133 	   int*                  maxcolactinf        /**< number of -inf contributions to minimal column activity */
1134 	   )
1135 	{
1136 	   SCIP_Real* valpnt;
1137 	   int* colpnt;
1138 	   int* colend;
1139 	   SCIP_Real val;
1140 	   int row;
1141 	
1142 	   assert(scip != NULL);
1143 	   assert(matrix != NULL);
1144 	   assert(lbdual != NULL);
1145 	   assert(ubdual != NULL);
1146 	   assert(maxcolact != NULL);
1147 	   assert(maxcolactinf != NULL);
1148 	
1149 	   /* init activities */
1150 	   maxcolact[col] = 0.0;
1151 	   maxcolactinf[col] = 0;
1152 	
1153 	   colpnt = SCIPmatrixGetColIdxPtr(matrix, col);
1154 	   colend = colpnt + SCIPmatrixGetColNNonzs(matrix, col);
1155 	   valpnt = SCIPmatrixGetColValPtr(matrix, col);
1156 	
1157 	   /* calculate column activities */
1158 	   for( ; colpnt < colend; colpnt++, valpnt++ )
1159 	   {
1160 	      row = *colpnt;
1161 	      val = *valpnt;
1162 	
1163 	      if( val > 0.0 )
1164 	      {
1165 	         if(SCIPisInfinity(scip, ubdual[row]))
1166 	            maxcolactinf[col]++;
1167 	         else
1168 	            maxcolact[col] += val * ubdual[row];
1169 	      }
1170 	      else if( val < 0.0 )
1171 	      {
1172 	         if(SCIPisInfinity(scip, -lbdual[row]))
1173 	            maxcolactinf[col]++;
1174 	         else
1175 	            maxcolact[col] += val * lbdual[row];
1176 	      }
1177 	   }
1178 	
1179 	   /* update column activities if infinity counters are greater 0 */
1180 	   if( maxcolactinf[col] > 0 )
1181 	      maxcolact[col] = SCIPinfinity(scip);
1182 	}
1183 	
1184 	
1185 	/** update minimal/maximal column activity infinity counters */
1186 	static
1187 	void infinityCountUpdate(
1188 	   SCIP*                 scip,               /**< SCIP main data structure */
1189 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1190 	   int                   row,                /**< row index */
1191 	   SCIP_Real*            lbdual,             /**< lower bounds of dual variables */
1192 	   SCIP_Real*            ubdual,             /**< upper bounds of dual variables */
1193 	   const SCIP_Bool*      isubimplied,        /**< flags indicating of the upper bound is implied */
1194 	   SCIP_Real*            mincolact,          /**< minimal column activities */
1195 	   int*                  mincolactinf,       /**< number of infinity contributions to minimal column activity */
1196 	   SCIP_Bool             ubinfchange,        /**< flag indicating if the upper bound has changed from infinity to a finite value */
1197 	   SCIP_Bool             lbinfchange         /**< flag indicating if the lower bound has changed from -infinity to a finite value */
1198 	   )
1199 	{
1200 	   SCIP_Real* valpnt;
1201 	   int* rowpnt;
1202 	   int* rowend;
1203 	   SCIP_Real val;
1204 	   int col;
1205 	
1206 	   rowpnt = SCIPmatrixGetRowIdxPtr(matrix, row);
1207 	   rowend = rowpnt + SCIPmatrixGetRowNNonzs(matrix, row);
1208 	   valpnt = SCIPmatrixGetRowValPtr(matrix, row);
1209 	
1210 	   /* look at all column entries present within row and update the
1211 	    * corresponding infinity counters. if one counter gets to zero,
1212 	    * then calculate this column activity new.
1213 	    */
1214 	
1215 	   for(; (rowpnt < rowend); rowpnt++, valpnt++ )
1216 	   {
1217 	      col = *rowpnt;
1218 	      val = *valpnt;
1219 	
1220 	      if( isubimplied[col] )
1221 	      {
1222 	         if( val < 0 )
1223 	         {
1224 	            if( ubinfchange )
1225 	            {
1226 	               assert(mincolactinf[col] > 0);
1227 	               mincolactinf[col]--;
1228 	            }
1229 	         }
1230 	         else if( val > 0 )
1231 	         {
1232 	            if( lbinfchange )
1233 	            {
1234 	               assert(mincolactinf[col] > 0);
1235 	               mincolactinf[col]--;
1236 	            }
1237 	         }
1238 	
1239 	         if( mincolactinf[col] == 0 )
1240 	            calcMinColActivity(scip, matrix, col, lbdual, ubdual, mincolact, mincolactinf);
1241 	      }
1242 	   }
1243 	}
1244 	
1245 	#ifdef SCIP_DEBUG
1246 	/** use LP calculations for determining the best dual variable bounds from a specific row index */
1247 	static
1248 	SCIP_RETCODE determineBestBounds(
1249 	   SCIP*                 scip,               /**< SCIP main data structure */
1250 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1251 	   int                   row,                /**< row index for dual bound calculations */
1252 	   SCIP_Bool             solveLP,            /**< flag indicating to solve subscip LP */
1253 	   SCIP_Real*            lowerbnddual,       /**< lower bound of dual variable */
1254 	   SCIP_Real*            upperbnddual        /**< upper bound of dual variable */
1255 	   )
1256 	{
1257 	   int i;
1258 	   int nrows;
1259 	   int ncols;
1260 	   int numberconvars;
1261 	   SCIP_VAR* var;
1262 	   SCIP_VAR** variables;
1263 	   SCIP_VAR** tmpvars;
1264 	   SCIP_Real* tmpcoef;
1265 	   SCIP_CONS** constraints;
1266 	   int numDualVars;
1267 	   SCIP* subscip;
1268 	   SCIP_RETCODE retcode;
1269 	   char name[SCIP_MAXSTRLEN+3];
1270 	   int fillcnt;
1271 	   int* colpnt;
1272 	   int* colend;
1273 	   SCIP_Real* valpnt;
1274 	   int* colmap;
1275 	
1276 	   *lowerbnddual = -SCIPinfinity(scip);
1277 	   *upperbnddual = SCIPinfinity(scip);
1278 	
1279 	   nrows = SCIPmatrixGetNRows(matrix);
1280 	   assert(0 <= row && row < nrows);
1281 	   ncols = SCIPmatrixGetNColumns(matrix);
1282 	
1283 	   SCIP_CALL( SCIPcreate(&subscip) );
1284 	   SCIP_CALL( SCIPcreateProbBasic(subscip, "subscip") );
1285 	   SCIP_CALL( SCIPincludeDefaultPlugins(subscip) );
1286 	
1287 	   /* avoid recursive calls */
1288 	   SCIP_CALL( SCIPsetIntParam(subscip, "presolving/dualinfer/maxrounds", 0) );
1289 	   SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) );
1290 	   SCIP_CALL( SCIPsetBoolParam(subscip, "misc/catchctrlc", TRUE) );
1291 	
1292 	   SCIP_CALL( SCIPallocBufferArray(scip, &colmap, ncols) );
1293 	   numberconvars = 0;
1294 	   for(i = 0; i < ncols; i++)
1295 	   {
1296 	      var = SCIPmatrixGetVar(matrix, i);
1297 	      if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT )
1298 	      {
1299 	         colmap[i] = numberconvars; /* start numbering with 0 */
1300 	         numberconvars++;
1301 	      }
1302 	      else
1303 	         colmap[i] = -1;
1304 	   }
1305 	   numDualVars = nrows + 2 * numberconvars;
1306 	
1307 	   /* create dual variables */
1308 	   SCIP_CALL( SCIPallocBufferArray(scip, &variables, numDualVars) );
1309 	   for( i = 0; i < nrows; i++ )
1310 	   {
1311 	      variables[i] = NULL;
1312 	      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "y%d", i);
1313 	      if( !SCIPmatrixIsRowRhsInfinity(matrix, i ) )
1314 	      {
1315 	         /* dual variable for equation or ranged row */
1316 	         SCIP_CALL( SCIPcreateVarBasic(subscip, &variables[i], name,
1317 	               -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
1318 	      }
1319 	      else
1320 	      {
1321 	         /* dual variable for >= inequality */
1322 	         SCIP_CALL( SCIPcreateVarBasic(subscip, &variables[i], name,
1323 	               0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
1324 	      }
1325 	      SCIP_CALL( SCIPaddVar(subscip, variables[i]) );
1326 	      assert( variables[i] != NULL );
1327 	   }
1328 	
1329 	   /* in addition, we introduce dual variables for the bounds,
1330 	      because we treat each continuous variable as a free variable */
1331 	   fillcnt = nrows;
1332 	   for( i = 0; i < numberconvars; i++ )
1333 	   {
1334 	      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "ylb%d", fillcnt);
1335 	      SCIP_CALL( SCIPcreateVarBasic(subscip, &variables[fillcnt], name,
1336 	            0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
1337 	      SCIP_CALL( SCIPaddVar(subscip, variables[fillcnt]) );
1338 	      assert( variables[fillcnt] != NULL );
1339 	      fillcnt++;
1340 	
1341 	      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "yub%d", fillcnt);
1342 	      SCIP_CALL( SCIPcreateVarBasic(subscip, &variables[fillcnt], name,
1343 	            0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
1344 	      SCIP_CALL( SCIPaddVar(subscip, variables[fillcnt]) );
1345 	      assert( variables[fillcnt] != NULL );
1346 	      fillcnt++;
1347 	   }
1348 	   assert(numDualVars == fillcnt);
1349 	
1350 	   SCIP_CALL( SCIPallocBufferArray(scip, &tmpvars, numDualVars) );
1351 	   SCIP_CALL( SCIPallocBufferArray(scip, &tmpcoef, numDualVars) );
1352 	
1353 	   SCIP_CALL( SCIPallocBufferArray(scip, &constraints, numberconvars) );
1354 	   for( i = 0; i <numberconvars; i++)
1355 	      constraints[i] = NULL;
1356 	
1357 	   for(i = 0; i < ncols; i++)
1358 	   {
1359 	      var = SCIPmatrixGetVar(matrix, i);
1360 	      if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT )
1361 	      {
1362 	         SCIP_Real objval = SCIPvarGetObj(var);
1363 	         int cidx = colmap[i];
1364 	         assert(0 <= cidx && cidx < numberconvars);
1365 	
1366 	         colpnt = SCIPmatrixGetColIdxPtr(matrix, i);
1367 	         colend = colpnt + SCIPmatrixGetColNNonzs(matrix, i);
1368 	         valpnt = SCIPmatrixGetColValPtr(matrix, i);
1369 	         fillcnt = 0;
1370 	         for( ; colpnt < colend; colpnt++, valpnt++ )
1371 	         {
1372 	            assert(0 <= *colpnt && *colpnt < nrows);
1373 	            assert(variables[*colpnt] != NULL);
1374 	            tmpvars[fillcnt] = variables[*colpnt];
1375 	            tmpcoef[fillcnt] = *valpnt;
1376 	            fillcnt++;
1377 	         }
1378 	
1379 	         /* consider dual variable for a lower bound */
1380 	         if(SCIPisGT(scip, SCIPvarGetLbGlobal(var), -SCIPinfinity(scip)))
1381 	         {
1382 	            assert(variables[nrows + 2 * cidx] != NULL);
1383 	            tmpvars[fillcnt] = variables[nrows + 2 * cidx];
1384 	            tmpcoef[fillcnt] = 1.0;
1385 	            fillcnt++;
1386 	         }
1387 	
1388 	         /* consider dual variable for an upper bound */
1389 	         if(SCIPisLT(scip, SCIPvarGetUbGlobal(var), SCIPinfinity(scip)))
1390 	         {
1391 	            assert(variables[nrows + 2 * cidx + 1] != NULL);
1392 	            tmpvars[fillcnt] = variables[nrows + 2 * cidx + 1];
1393 	            tmpcoef[fillcnt] = -1.0;
1394 	            fillcnt++;
1395 	         }
1396 	
1397 	         /* because we treat the continuous columns as free variable,
1398 	            we need here an equality */
1399 	         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "c%d", cidx);
1400 	         SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &constraints[cidx], name,
1401 	               fillcnt, tmpvars, tmpcoef, objval, objval) );
1402 	         SCIP_CALL( SCIPaddCons(subscip, constraints[cidx]) );
1403 	      }
1404 	   }
1405 	
1406 	   /* determine lower dual bound via a minimization problem */
1407 	   SCIP_CALL( SCIPsetObjsense(subscip,SCIP_OBJSENSE_MINIMIZE) );
1408 	   SCIP_CALL( SCIPchgVarObj(subscip, variables[row], 1.0) );
1409 	   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "dbg_min_%s.lp", SCIPvarGetName(variables[row]));
1410 	   SCIP_CALL( SCIPwriteOrigProblem(subscip, name, "lp", FALSE) );
1411 	   if( solveLP )
1412 	   {
1413 	      retcode = SCIPsolve(subscip);
1414 	      if( retcode != SCIP_OKAY )
1415 	         SCIPwarningMessage(scip, "Error subscip: <%d>\n", retcode);
1416 	      else
1417 	      {
1418 	         if( SCIPgetStatus(subscip) == SCIP_STATUS_OPTIMAL )
1419 	         {
1420 	            SCIP_SOL* sol;
1421 	            SCIP_Bool feasible;
1422 	            sol = SCIPgetBestSol(subscip);
1423 	            SCIP_CALL( SCIPcheckSolOrig(subscip, sol, &feasible, TRUE, TRUE) );
1424 	
1425 	            if(feasible)
1426 	               *lowerbnddual = SCIPgetSolOrigObj(subscip, sol);
1427 	         }
1428 	      }
1429 	      SCIP_CALL( SCIPfreeTransform(subscip) );
1430 	   }
1431 	
1432 	   /* determine upper dual bound via a maximization problem */
1433 	   SCIP_CALL( SCIPsetObjsense(subscip,SCIP_OBJSENSE_MAXIMIZE) );
1434 	   SCIP_CALL( SCIPchgVarObj(subscip, variables[row], 1.0) );
1435 	   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "dbg_max_%s.lp", SCIPvarGetName(variables[row]));
1436 	   SCIP_CALL( SCIPwriteOrigProblem(subscip, name, "lp", FALSE) );
1437 	   if( solveLP )
1438 	   {
1439 	      retcode = SCIPsolve(subscip);
1440 	      if( retcode != SCIP_OKAY )
1441 	         SCIPwarningMessage(scip, "Error subscip: <%d>\n", retcode);
1442 	      else
1443 	      {
1444 	         if( SCIPgetStatus(subscip) == SCIP_STATUS_OPTIMAL )
1445 	         {
1446 	            SCIP_SOL* sol;
1447 	            SCIP_Bool feasible;
1448 	            sol = SCIPgetBestSol(subscip);
1449 	            SCIP_CALL( SCIPcheckSolOrig(subscip, sol, &feasible, TRUE, TRUE) );
1450 	
1451 	            if(feasible)
1452 	               *upperbnddual = SCIPgetSolOrigObj(subscip, sol);
1453 	         }
1454 	      }
1455 	      SCIP_CALL( SCIPfreeTransform(subscip) );
1456 	   }
1457 	
1458 	   /* release variables and constraints */
1459 	   for( i = 0; i < numDualVars; i++ )
1460 	   {
1461 	      if(variables[i] != NULL)
1462 	         SCIP_CALL( SCIPreleaseVar(subscip, &variables[i]) );
1463 	   }
1464 	   for( i = 0; i < numberconvars; i++ )
1465 	   {
1466 	      if(constraints[i] != NULL)
1467 	         SCIP_CALL( SCIPreleaseCons(subscip, &constraints[i]) );
1468 	   }
1469 	
1470 	   SCIPfreeBufferArray(scip, &constraints);
1471 	   SCIPfreeBufferArray(scip, &tmpcoef);
1472 	   SCIPfreeBufferArray(scip, &tmpvars);
1473 	   SCIPfreeBufferArray(scip, &variables);
1474 	   SCIP_CALL( SCIPfree(&subscip) );
1475 	   SCIPfreeBufferArray(scip, &colmap);
1476 	
1477 	   return SCIP_OKAY;
1478 	}
1479 	#endif
1480 	
1481 	/** update bounds of the dual variables */
1482 	static
1483 	void updateDualBounds(
1484 	   SCIP*                 scip,               /**< SCIP main data structure */
1485 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1486 	   SCIP_Real             objval,             /**< objective function value */
1487 	   SCIP_Real             val,                /**< matrix coefficient */
1488 	   int                   row,                /**< row index */
1489 	   SCIP_Real             mincolresact,       /**< minimal column residual activity */
1490 	   SCIP_Real*            lbdual,             /**< dual lower bounds */
1491 	   SCIP_Real*            ubdual,             /**< dual upper bounds */
1492 	   int*                  boundchanges,       /**< counter for the number of bound changes */
1493 	   SCIP_Bool*            ubinfchange,        /**< flag indicating an upper bound change from infinite to finite */
1494 	   SCIP_Bool*            lbinfchange         /**< flag indicating a lower bound change from infinite to finite */
1495 	   )
1496 	{
1497 	   SCIP_Real newlbdual;
1498 	   SCIP_Real newubdual;
1499 	
1500 	   assert(scip != NULL);
1501 	   assert(matrix != NULL);
1502 	   assert(lbdual != NULL);
1503 	   assert(ubdual != NULL);
1504 	   assert(boundchanges != NULL);
1505 	   assert(ubinfchange != NULL);
1506 	   assert(lbinfchange != NULL);
1507 	
1508 	   *ubinfchange = FALSE;
1509 	   *lbinfchange = FALSE;
1510 	
1511 	   if( !SCIPisInfinity(scip, -mincolresact) )
1512 	   {
1513 	      if( val > 0 )
1514 	      {
1515 	         newubdual = (objval - mincolresact) / val;
1516 	
1517 	         if( newubdual < ubdual[row] )
1518 	         {
1519 	            /* accept the new upper bound only if the numerics are reliable */
1520 	            if( SCIPisLE(scip,lbdual[row],newubdual) )
1521 	            {
1522 	               if( SCIPisInfinity(scip, ubdual[row]) )
1523 	                  *ubinfchange = TRUE;
1524 	
1525 	               ubdual[row] = newubdual;
1526 	               (*boundchanges)++;
1527 	            }
1528 	         }
1529 	      }
1530 	      else if( val < 0 )
1531 	      {
1532 	         newlbdual = (objval - mincolresact) / val;
1533 	
1534 	         if( newlbdual > lbdual[row] )
1535 	         {
1536 	            /* accept the new lower bound only if the numerics are reliable */
1537 	            if( SCIPisLE(scip,newlbdual,ubdual[row]) )
1538 	            {
1539 	               if( SCIPisInfinity(scip, -lbdual[row]) )
1540 	                  *lbinfchange = TRUE;
1541 	
1542 	               lbdual[row] = newlbdual;
1543 	               (*boundchanges)++;
1544 	            }
1545 	         }
1546 	      }
1547 	   }
1548 	}
1549 	
1550 	/** dual bound strengthening */
1551 	static
1552 	SCIP_RETCODE dualBoundStrengthening(
1553 	   SCIP*                 scip,               /**< SCIP main data structure */
1554 	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
1555 	   SCIP_PRESOLDATA*      presoldata,         /**< presolver data structure */
1556 	   FIXINGDIRECTION*      varstofix,          /**< array holding information for later upper/lower bound fixing */
1557 	   int*                  npossiblefixings,   /**< number of possible fixings */
1558 	   SIDECHANGE*           sidestochange,      /**< array holding if this is an implied equality */
1559 	   int*                  npossiblesidechanges/**< number of possible equality changes */
1560 	   )
1561 	{
1562 	   SCIP_Real* lbdual;
1563 	   SCIP_Real* ubdual;
1564 	   SCIP_Real* mincolact;
1565 	   int* mincolactinf;
1566 	   SCIP_Real* maxcolact;
1567 	   int* maxcolactinf;
1568 	   int* colpnt;
1569 	   int* colend;
1570 	   SCIP_Real* valpnt;
1571 	   int boundchanges;
1572 	   int loops;
1573 	   int i;
1574 	   int j;
1575 	   int k;
1576 	   int nrows;
1577 	   int ncols;
1578 	   SCIP_Bool* isubimplied;
1579 	   SCIP_Bool* islbimplied;
1580 	   SCIP_Real* tmplbs;
1581 	   SCIP_Real* tmpubs;
1582 	   SCIP_VAR* var;
1583 	   int* implubvars;
1584 	   int nimplubvars;
1585 	
1586 	   SCIP_Longint maxhashes;
1587 	   int maxlen;
1588 	   int pospp;
1589 	   int listsizepp;
1590 	   int posmm;
1591 	   int listsizemm;
1592 	   int pospm;
1593 	   int listsizepm;
1594 	   int posmp;
1595 	   int listsizemp;
1596 	
1597 	   int* hashlistpp;
1598 	   int* hashlistmm;
1599 	   int* hashlistpm;
1600 	   int* hashlistmp;
1601 	
1602 	   int* colidxlistpp;
1603 	   int* colidxlistmm;
1604 	   int* colidxlistpm;
1605 	   int* colidxlistmp;
1606 	
1607 	   int block1start;
1608 	   int block1end;
1609 	   int block2start;
1610 	   int block2end;
1611 	
1612 	   SCIP_HASHSET* pairhashset;
1613 	   SCIP_Real* colvalptr;
1614 	   int* colidxptr;
1615 	
1616 	   assert(scip != NULL);
1617 	   assert(matrix != NULL);
1618 	   assert(varstofix != NULL);
1619 	   assert(npossiblefixings != NULL);
1620 	   assert(sidestochange != NULL);
1621 	   assert(npossiblesidechanges != NULL);
1622 	
1623 	   nrows = SCIPmatrixGetNRows(matrix);
1624 	   ncols = SCIPmatrixGetNColumns(matrix);
1625 	
1626 	   SCIP_CALL( SCIPallocBufferArray(scip, &tmplbs, ncols) );
1627 	   SCIP_CALL( SCIPallocBufferArray(scip, &tmpubs, ncols) );
1628 	   for( i = 0; i < ncols; i++ )
1629 	   {
1630 	      var = SCIPmatrixGetVar(matrix, i);
1631 	      tmplbs[i] = SCIPvarGetLbLocal(var);
1632 	      tmpubs[i] = SCIPvarGetUbLocal(var);
1633 	   }
1634 	
1635 	   /* verify which bounds of continuous variables are implied */
1636 	   SCIP_CALL( SCIPallocBufferArray(scip, &isubimplied, ncols) );
1637 	   SCIP_CALL( SCIPallocBufferArray(scip, &islbimplied, ncols) );
1638 	   SCIP_CALL( SCIPallocBufferArray(scip, &implubvars, ncols) );
1639 	   nimplubvars = 0;
1640 	   for( i = 0; i < ncols; i++ )
1641 	   {
1642 	      var = SCIPmatrixGetVar(matrix, i);
1643 	
1644 	      if( SCIPmatrixUplockConflict(matrix, i) || SCIPmatrixDownlockConflict(matrix, i) ||
1645 	         (SCIPvarGetType(var) != SCIP_VARTYPE_CONTINUOUS && SCIPvarGetType(var) != SCIP_VARTYPE_IMPLINT) )
1646 	      {
1647 	         /* we don't care about integral variables or variables that have conflicting locks */
1648 	         isubimplied[i] = FALSE;
1649 	         islbimplied[i] = FALSE;
1650 	      }
1651 	      else
1652 	      {
1653 	         getImpliedBounds(scip, matrix, i, tmplbs, tmpubs, &(isubimplied[i]), &(islbimplied[i]));
1654 	
1655 	         /* if a continuous variable has a not implied upper bound we can
1656 	          * not use this variable (column) for propagating dual bounds.
1657 	          * not implied lowers bound can usually be treated.
1658 	          */
1659 	
1660 	         /* collect continuous variables with implied upper bound */
1661 	         if( isubimplied[i] )
1662 	         {
1663 	            implubvars[nimplubvars] = i;
1664 	            nimplubvars++;
1665 	         }
1666 	
1667 	         /* reset implied bounds for further detections of other implied bounds */
1668 	         if( isubimplied[i] )
1669 	            tmpubs[i] = SCIPinfinity(scip);
1670 	
1671 	         if( islbimplied[i] )
1672 	            tmplbs[i] = -SCIPinfinity(scip);
1673 	      }
1674 	   }
1675 	
1676 	   /* initialize bounds of the dual variables */
1677 	   SCIP_CALL( SCIPallocBufferArray(scip, &lbdual, nrows) );
1678 	   SCIP_CALL( SCIPallocBufferArray(scip, &ubdual, nrows) );
1679 	   for( i = 0; i < nrows; i++ )
1680 	   {
1681 	      if( !SCIPmatrixIsRowRhsInfinity(matrix, i) )
1682 	      {
1683 	         /* dual free variable for equation or ranged row */
1684 	         lbdual[i] = -SCIPinfinity(scip);
1685 	         ubdual[i] = SCIPinfinity(scip);
1686 	      }
1687 	      else
1688 	      {
1689 	         /* dual variable for >= inequality */
1690 	         lbdual[i] = 0.0;
1691 	         ubdual[i] = SCIPinfinity(scip);
1692 	      }
1693 	   }
1694 	
1695 	   /* run convex combination on pairs of continuous variables (columns) using Belotti's algorithm */
1696 	   if( nimplubvars >= 2 && presoldata->usetwocolcombine )
1697 	   {
1698 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &hashlistpp, ncols) );
1699 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &hashlistmm, ncols) );
1700 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &hashlistpm, ncols) );
1701 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &hashlistmp, ncols) );
1702 	
1703 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &colidxlistpp, ncols) );
1704 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &colidxlistmm, ncols) );
1705 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &colidxlistpm, ncols) );
1706 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &colidxlistmp, ncols) );
1707 	
1708 	      pospp = 0;
1709 	      posmm = 0;
1710 	      pospm = 0;
1711 	      posmp = 0;
1712 	      listsizepp = ncols;
1713 	      listsizemm = ncols;
1714 	      listsizepm = ncols;
1715 	      listsizemp = ncols;
1716 	      maxhashes = presoldata->maxhashfac == -1 ? SCIP_LONGINT_MAX : (((SCIP_Longint)ncols) * presoldata->maxhashfac);
1717 	
1718 	      for( i = 0; i < nimplubvars; i++)
1719 	      {
1720 	         if( ((SCIP_Longint)pospp) + posmm + pospm + posmp > maxhashes )
1721 	            break;
1722 	
1723 	         colvalptr = SCIPmatrixGetColValPtr(matrix, implubvars[i]);
1724 	         colidxptr = SCIPmatrixGetColIdxPtr(matrix, implubvars[i]);
1725 	         maxlen = MIN(presoldata->maxconsiderednonzeros, SCIPmatrixGetColNNonzs(matrix, implubvars[i])); /*lint !e666*/
1726 	         for( j = 0; j < maxlen; j++)
1727 	         {
1728 	            for( k = j + 1; k < maxlen; k++)
1729 	            {
1730 	               if( SCIPisPositive(scip, colvalptr[j]) )
1731 	               {
1732 	                  if(SCIPisPositive(scip, colvalptr[k]) )
1733 	                  {
1734 	                     SCIP_CALL( addEntry(scip, &pospp, &listsizepp, &hashlistpp, &colidxlistpp,
1735 	                        hashIndexPair(colidxptr[j], colidxptr[k]), implubvars[i]) );
1736 	                  }
1737 	                  else
1738 	                  {
1739 	                     SCIP_CALL( addEntry(scip, &pospm, &listsizepm, &hashlistpm, &colidxlistpm,
1740 	                        hashIndexPair(colidxptr[j], colidxptr[k]), implubvars[i]) );
1741 	                  }
1742 	               }
1743 	               else
1744 	               {
1745 	                  if(SCIPisPositive(scip, colvalptr[k]) )
1746 	                  {
1747 	                     SCIP_CALL( addEntry(scip, &posmp, &listsizemp, &hashlistmp, &colidxlistmp,
1748 	                        hashIndexPair(colidxptr[j], colidxptr[k]), implubvars[i]) );
1749 	                  }
1750 	                  else
1751 	                  {
1752 	                     SCIP_CALL( addEntry(scip, &posmm, &listsizemm, &hashlistmm, &colidxlistmm,
1753 	                        hashIndexPair(colidxptr[j], colidxptr[k]), implubvars[i]) );
1754 	                  }
1755 	               }
1756 	            }
1757 	         }
1758 	      }
1759 	#ifdef SCIP_MORE_DEBUG
1760 	      SCIPdebugMsg(scip, "hashlist sizes: pp %d, mm %d, pm %d, mp %d \n", pospp, posmm, pospm, posmp);
1761 	#endif
1762 	      SCIPsortIntInt(hashlistpp, colidxlistpp, pospp);
1763 	      SCIPsortIntInt(hashlistmm, colidxlistmm, posmm);
1764 	      SCIPsortIntInt(hashlistpm, colidxlistpm, pospm);
1765 	      SCIPsortIntInt(hashlistmp, colidxlistmp, posmp);
1766 	
1767 	      SCIP_CALL( SCIPhashsetCreate(&pairhashset, SCIPblkmem(scip), 1) );
1768 	
1769 	      /* Process pp and mm lists */
1770 	      if( pospp > 0 && posmm > 0 )
1771 	      {
1772 	         SCIP_Longint ncombines;
1773 	         SCIP_Longint maxcombines;
1774 	         SCIP_Bool finished;
1775 	         SCIP_Bool success;
1776 	         int combinefails;
1777 	         int retrievefails;
1778 	         COLPAIR colpair;
1779 	
1780 	         finished = FALSE;
1781 	         block1start = 0;
1782 	         block1end = 0;
1783 	         block2start = 0;
1784 	         block2end = 0;
1785 	
1786 	         maxcombines = presoldata->maxpairfac == -1 ? SCIP_LONGINT_MAX : (((SCIP_Longint)ncols) * presoldata->maxpairfac);
1787 	
1788 	         ncombines = 0;
1789 	         combinefails = 0;
1790 	         retrievefails = 0;
1791 	         findNextBlock(hashlistpp, pospp, &block1start, &block1end);
1792 	         findNextBlock(hashlistmm, posmm, &block2start, &block2end);
1793 	#ifdef SCIP_MORE_DEBUG
1794 	         SCIPdebugMsg(scip, "processing pp and mm\n");
1795 	#endif
1796 	
1797 	         // same as in the rworowbnd presolver - both while loops to basically the same with one using pp and mm and the other pm and mp
1798 	         // I would write an additional function and remove the code duplication
1799 	         while( !finished )
1800 	         {
1801 	            if( hashlistpp[block1start] == hashlistmm[block2start] )
1802 	            {
1803 	               for( i = block1start; i < block1end; i++ )
1804 	               {
1805 	                  for( j = block2start; j < block2end; j++ )
1806 	                  {
1807 	                     if( colidxlistpp[i] != colidxlistmm[j] )
1808 	                     {
1809 	                        colpair.col1idx = MIN(colidxlistpp[i], colidxlistmm[j]);
1810 	                        colpair.col2idx = MAX(colidxlistpp[i], colidxlistmm[j]);
1811 	
1812 	                        if( !SCIPhashsetExists(pairhashset, encodeColPair(&colpair)) )
1813 	                        {
1814 	                           int* colpnt1       = SCIPmatrixGetColIdxPtr(matrix, colpair.col1idx);
1815 	                           SCIP_Real* valpnt1 = SCIPmatrixGetColValPtr(matrix, colpair.col1idx);
1816 	                           int* colpnt2       = SCIPmatrixGetColIdxPtr(matrix, colpair.col2idx);
1817 	                           SCIP_Real* valpnt2 = SCIPmatrixGetColValPtr(matrix, colpair.col2idx);
1818 	                           SCIP_Real obj1     = SCIPvarGetObj(SCIPmatrixGetVar(matrix, colpair.col1idx));
1819 	                           SCIP_Real obj2     = SCIPvarGetObj(SCIPmatrixGetVar(matrix, colpair.col2idx));
1820 	                           int collen1        = SCIPmatrixGetColNNonzs(matrix, colpair.col1idx);
1821 	                           int collen2        = SCIPmatrixGetColNNonzs(matrix, colpair.col2idx);
1822 	
1823 	                           success = FALSE;
1824 	
1825 	                           SCIP_CALL( combineCols(scip, colpnt1, colpnt2, valpnt1, valpnt2, obj1, obj2, collen1,
1826 	                                    collen2, nrows, TRUE, TRUE, lbdual, ubdual, &success) );
1827 	
1828 	                           if( success )
1829 	                              combinefails = 0;
1830 	                           else
1831 	                              combinefails++;
1832 	
1833 	                           SCIP_CALL( SCIPhashsetInsert(pairhashset, SCIPblkmem(scip), encodeColPair(&colpair)) );
1834 	                           ncombines++;
1835 	                           if( ncombines >= maxcombines || combinefails >= presoldata->maxcombinefails )
1836 	                              finished = TRUE;
1837 	#ifdef SCIP_MORE_DEBUG
1838 	                           SCIPdebugMsg(scip, "pm/mp: %d retrievefails before reset, %d combines\n", retrievefails, ncombines);
1839 	#endif
1840 	                           retrievefails = 0;
1841 	                        }
1842 	                        else if( retrievefails < presoldata->maxretrievefails )
1843 	                           retrievefails++;
1844 	                        else
1845 	                           finished = TRUE;
1846 	                     }
1847 	                     if( finished )
1848 	                        break;
1849 	                  }
1850 	                  if( finished )
1851 	                     break;
1852 	               }
1853 	
1854 	               if( block1end < pospp && block2end < posmm )
1855 	               {
1856 	                  findNextBlock(hashlistpp, pospp, &block1start, &block1end);
1857 	                  findNextBlock(hashlistmm, posmm, &block2start, &block2end);
1858 	               }
1859 	               else
1860 	                  finished = TRUE;
1861 	            }
1862 	            else if( hashlistpp[block1start] < hashlistmm[block2start] && block1end < pospp )
1863 	               findNextBlock(hashlistpp, pospp, &block1start, &block1end);
1864 	            else if( hashlistpp[block1start] > hashlistmm[block2start] && block2end < posmm )
1865 	               findNextBlock(hashlistmm, posmm, &block2start, &block2end);
1866 	            else
1867 	               finished = TRUE;
1868 	         }
1869 	      }
1870 	
1871 	      /* Process pm and mp lists */
1872 	      if( pospm > 0 && posmp > 0 )
1873 	      {
1874 	         SCIP_Longint maxcombines;
1875 	         SCIP_Longint ncombines;
1876 	         SCIP_Bool finished;
1877 	         SCIP_Bool success;
1878 	         int combinefails;
1879 	         int retrievefails;
1880 	         COLPAIR colpair;
1881 	
1882 	         finished = FALSE;
1883 	         block1start = 0;
1884 	         block1end = 0;
1885 	         block2start = 0;
1886 	         block2end = 0;
1887 	
1888 	         maxcombines = presoldata->maxpairfac == -1 ? SCIP_LONGINT_MAX : (((SCIP_Longint)ncols) * presoldata->maxpairfac);
1889 	
1890 	         ncombines = 0;
1891 	         combinefails = 0;
1892 	         retrievefails = 0;
1893 	         findNextBlock(hashlistpm, pospm, &block1start, &block1end);
1894 	         findNextBlock(hashlistmp, posmp, &block2start, &block2end);
1895 	#ifdef SCIP_MORE_DEBUG
1896 	         SCIPdebugMsg(scip, "processing pm and mp\n");
1897 	#endif
1898 	
1899 	         while( !finished )
1900 	         {
1901 	            if( hashlistpm[block1start] == hashlistmp[block2start] )
1902 	            {
1903 	               for( i = block1start; i < block1end; i++ )
1904 	               {
1905 	                  for( j = block2start; j < block2end; j++ )
1906 	                  {
1907 	                     if( colidxlistpm[i] != colidxlistmp[j] )
1908 	                     {
1909 	                        colpair.col1idx = MIN(colidxlistpm[i], colidxlistmp[j]);
1910 	                        colpair.col2idx = MAX(colidxlistpm[i], colidxlistmp[j]);
1911 	
1912 	                        if( !SCIPhashsetExists(pairhashset, encodeColPair(&colpair)) )
1913 	                        {
1914 	                           int* colpnt1       = SCIPmatrixGetColIdxPtr(matrix, colpair.col1idx);
1915 	                           SCIP_Real* valpnt1 = SCIPmatrixGetColValPtr(matrix, colpair.col1idx);
1916 	                           int* colpnt2       = SCIPmatrixGetColIdxPtr(matrix, colpair.col2idx);
1917 	                           SCIP_Real* valpnt2 = SCIPmatrixGetColValPtr(matrix, colpair.col2idx);
1918 	                           SCIP_Real obj1     = SCIPvarGetObj(SCIPmatrixGetVar(matrix, colpair.col1idx));
1919 	                           SCIP_Real obj2     = SCIPvarGetObj(SCIPmatrixGetVar(matrix, colpair.col2idx));
1920 	                           int collen1        = SCIPmatrixGetColNNonzs(matrix, colpair.col1idx);
1921 	                           int collen2        = SCIPmatrixGetColNNonzs(matrix, colpair.col2idx);
1922 	
1923 	                           success = FALSE;
1924 	
1925 	                           SCIP_CALL( combineCols(scip, colpnt1, colpnt2, valpnt1, valpnt2, obj1, obj2, collen1,
1926 	                                    collen2, nrows, TRUE, TRUE, lbdual, ubdual, &success) );
1927 	
1928 	                           if( success )
1929 	                              combinefails = 0;
1930 	                           else
1931 	                              combinefails++;
1932 	
1933 	                           SCIP_CALL( SCIPhashsetInsert(pairhashset, SCIPblkmem(scip), encodeColPair(&colpair)) );
1934 	                           ncombines++;
1935 	                           if( ncombines >= maxcombines || combinefails >= presoldata->maxcombinefails )
1936 	                              finished = TRUE;
1937 	
1938 	                           retrievefails = 0;
1939 	                        }
1940 	                        else if( retrievefails < presoldata->maxretrievefails )
1941 	                           retrievefails++;
1942 	                        else
1943 	                           finished = TRUE;
1944 	                     }
1945 	                     if( finished )
1946 	                        break;
1947 	                  }
1948 	                  if( finished )
1949 	                     break;
1950 	               }
1951 	
1952 	               if( block1end < pospm && block2end < posmp )
1953 	               {
1954 	                  findNextBlock(hashlistpm, pospm, &block1start, &block1end);
1955 	                  findNextBlock(hashlistmp, posmp, &block2start, &block2end);
1956 	               }
1957 	               else
1958 	                  finished = TRUE;
1959 	            }
1960 	            else if( hashlistpm[block1start] < hashlistmp[block2start] && block1end < pospm )
1961 	               findNextBlock(hashlistpm, pospm, &block1start, &block1end);
1962 	            else if( hashlistpm[block1start] > hashlistmp[block2start] && block2end < posmp )
1963 	               findNextBlock(hashlistmp, posmp, &block2start, &block2end);
1964 	            else
1965 	               finished = TRUE;
1966 	         }
1967 	      }
1968 	
1969 	      SCIPhashsetFree(&pairhashset, SCIPblkmem(scip));
1970 	      SCIPfreeBlockMemoryArray(scip, &colidxlistmp, listsizemp);
1971 	      SCIPfreeBlockMemoryArray(scip, &colidxlistpm, listsizepm);
1972 	      SCIPfreeBlockMemoryArray(scip, &colidxlistmm, listsizemm);
1973 	      SCIPfreeBlockMemoryArray(scip, &colidxlistpp, listsizepp);
1974 	      SCIPfreeBlockMemoryArray(scip, &hashlistmp, listsizemp);
1975 	      SCIPfreeBlockMemoryArray(scip, &hashlistpm, listsizepm);
1976 	      SCIPfreeBlockMemoryArray(scip, &hashlistmm, listsizemm);
1977 	      SCIPfreeBlockMemoryArray(scip, &hashlistpp, listsizepp);
1978 	
1979 	#ifdef SCIP_MORE_DEBUG
1980 	      SCIPdebugMsg(scip, "CombCols:\n");
1981 	      for( i = 0; i < nrows; i++ )
1982 	      {
1983 	         assert(SCIPisLE(scip,lbdual[i],ubdual[i]));
1984 	         SCIPdebugMsg(scip, "y%d=[%g,%g]\n",i,lbdual[i],ubdual[i]);
1985 	      }
1986 	      SCIPdebugMsg(scip,"\n");
1987 	#endif
1988 	   }
1989 	
1990 	   SCIP_CALL( SCIPallocBufferArray(scip, &mincolact, ncols) );
1991 	   SCIP_CALL( SCIPallocBufferArray(scip, &mincolactinf, ncols) );
1992 	
1993 	   /* apply dual bound strengthening */
1994 	   loops = 0;
1995 	   boundchanges = 1;
1996 	   while( 0 < boundchanges && loops < presoldata->maxdualbndloops )
1997 	   {
1998 	      loops++;
1999 	      boundchanges = 0;
2000 	
2001 	      for( i = 0; i < nimplubvars; i++ )
2002 	      {
2003 	         assert(SCIPvarGetType(SCIPmatrixGetVar(matrix, implubvars[i])) == SCIP_VARTYPE_CONTINUOUS ||
2004 	            SCIPvarGetType(SCIPmatrixGetVar(matrix, implubvars[i])) == SCIP_VARTYPE_IMPLINT);
2005 	         calcMinColActivity(scip, matrix, implubvars[i], lbdual, ubdual, mincolact, mincolactinf);
2006 	      }
2007 	
2008 	      for( i = 0; i < nimplubvars; i++ )
2009 	      {
2010 	         SCIP_Real objval;
2011 	         SCIP_Bool ubinfchange;
2012 	         SCIP_Bool lbinfchange;
2013 	         int col;
2014 	
2015 	         col = implubvars[i];
2016 	         var = SCIPmatrixGetVar(matrix, col);
2017 	
2018 	         objval = SCIPvarGetObj(var);
2019 	         colpnt = SCIPmatrixGetColIdxPtr(matrix, col);
2020 	         colend = colpnt + SCIPmatrixGetColNNonzs(matrix, col);
2021 	         valpnt = SCIPmatrixGetColValPtr(matrix, col);
2022 	
2023 	         for( ; colpnt < colend; colpnt++, valpnt++ )
2024 	         {
2025 	            int row;
2026 	            SCIP_Real val;
2027 	            SCIP_Real mincolresact;
2028 	
2029 	            row = *colpnt;
2030 	            val = *valpnt;
2031 	
2032 	            calcMinColActResidual(scip, matrix, col, row, val, lbdual, ubdual,
2033 	               mincolact, mincolactinf, &mincolresact);
2034 	
2035 	            updateDualBounds(scip, matrix, objval, val, row, mincolresact,
2036 	               lbdual, ubdual, &boundchanges, &ubinfchange, &lbinfchange);
2037 	
2038 	            if( ubinfchange || lbinfchange )
2039 	               infinityCountUpdate(scip, matrix, row, lbdual, ubdual, isubimplied,
2040 	                  mincolact, mincolactinf, ubinfchange, lbinfchange);
2041 	         }
2042 	      }
2043 	   }
2044 	
2045 	#ifdef SCIP_MORE_DEBUG
2046 	   SCIPdebugMsg(scip, "BndStr:\n");
2047 	   for( i = 0; i < nrows; i++ )
2048 	   {
2049 	      assert(SCIPisLE(scip,lbdual[i],ubdual[i]));
2050 	      SCIPdebugMsg(scip, "y%d=[%g,%g]\n",i,lbdual[i],ubdual[i]);
2051 	   }
2052 	   SCIPdebugMsg(scip,"\n");
2053 	#endif
2054 	
2055 	   SCIP_CALL( SCIPallocBufferArray(scip, &maxcolact, ncols) );
2056 	   SCIP_CALL( SCIPallocBufferArray(scip, &maxcolactinf, ncols) );
2057 	
2058 	   /* calculate final minimal and maximal column activities */
2059 	   for( i = 0; i < ncols; i++ )
2060 	   {
2061 	      calcMinColActivity(scip, matrix, i, lbdual, ubdual, mincolact, mincolactinf);
2062 	      calcMaxColActivity(scip, matrix, i, lbdual, ubdual, maxcolact, maxcolactinf);
2063 	   }
2064 	
2065 	   for( i = 0; i < ncols; i++ )
2066 	   {
2067 	      SCIP_Real objval;
2068 	
2069 	      var = SCIPmatrixGetVar(matrix, i);
2070 	
2071 	      /* do not fix variables if the locks do not match */
2072 	      if( SCIPmatrixUplockConflict(matrix, i) || SCIPmatrixDownlockConflict(matrix, i) )
2073 	         continue;
2074 	
2075 	      objval = SCIPvarGetObj(var);
2076 	
2077 	      /* c_j - sup(y^T A_{.j}) > 0 => fix x_j to its lower bound */
2078 	      if( SCIPisGT(scip, objval, maxcolact[i]) && varstofix[i] == NOFIX )
2079 	      {
2080 	         if( SCIPisGT(scip, SCIPvarGetLbGlobal(var), -SCIPinfinity(scip)) )
2081 	         {
2082 	            varstofix[i] = FIXATLB;
2083 	            (*npossiblefixings)++;
2084 	         }
2085 	      }
2086 	
2087 	      /* c_j - inf(y^T A_{.j}) < 0 => fix x_j to its upper bound */
2088 	      if( SCIPisLT(scip, objval, mincolact[i]) && varstofix[i] == NOFIX )
2089 	      {
2090 	         if( SCIPisLT(scip, SCIPvarGetUbGlobal(var), SCIPinfinity(scip)) )
2091 	         {
2092 	            varstofix[i] = FIXATUB;
2093 	            (*npossiblefixings)++;
2094 	         }
2095 	      }
2096 	   }
2097 	
2098 	   for( i = 0; i < nrows; i++ )
2099 	   {
2100 	      /* implied equality: y_i > 0 =>  A_{i.}x - b_i = 0 */
2101 	      if( SCIPmatrixIsRowRhsInfinity(matrix, i) )
2102 	      {
2103 	         if( SCIPisGT(scip, lbdual[i], 0.0) && (sidestochange[i] == NOCHANGE) )
2104 	         {
2105 	            /* change >= inequality to equality */
2106 	            sidestochange[i] = RHSTOLHS;
2107 	            (*npossiblesidechanges)++;
2108 	         }
2109 	      }
2110 	      else
2111 	      {
2112 	         if( !SCIPmatrixIsRowRhsInfinity(matrix, i) &&
2113 	            !SCIPisEQ(scip,SCIPmatrixGetRowLhs(matrix, i),SCIPmatrixGetRowRhs(matrix, i)) )
2114 	         {
2115 	            /* for ranged rows we have to decide which side (lhs or rhs) determines the equality */
2116 	            if( SCIPisGT(scip, lbdual[i], 0.0) && sidestochange[i]==NOCHANGE )
2117 	            {
2118 	               sidestochange[i] = RHSTOLHS;
2119 	               (*npossiblesidechanges)++;
2120 	            }
2121 	
2122 	            if( SCIPisLT(scip, ubdual[i], 0.0) && sidestochange[i]==NOCHANGE)
2123 	            {
2124 	               sidestochange[i] = LHSTORHS;
2125 	               (*npossiblesidechanges)++;
2126 	            }
2127 	         }
2128 	      }
2129 	   }
2130 	
2131 	   SCIPfreeBufferArray(scip, &maxcolactinf);
2132 	   SCIPfreeBufferArray(scip, &maxcolact);
2133 	   SCIPfreeBufferArray(scip, &mincolactinf);
2134 	   SCIPfreeBufferArray(scip, &mincolact);
2135 	
2136 	   SCIPfreeBufferArray(scip, &ubdual);
2137 	   SCIPfreeBufferArray(scip, &lbdual);
2138 	   SCIPfreeBufferArray(scip, &implubvars);
2139 	   SCIPfreeBufferArray(scip, &islbimplied);
2140 	   SCIPfreeBufferArray(scip, &isubimplied);
2141 	   SCIPfreeBufferArray(scip, &tmpubs);
2142 	   SCIPfreeBufferArray(scip, &tmplbs);
2143 	
2144 	   return SCIP_OKAY;
2145 	}
2146 	
2147 	/*
2148 	 * Callback methods of presolver
2149 	 */
2150 	
2151 	/** copy method for constraint handler plugins (called when SCIP copies plugins) */
2152 	static
2153 	SCIP_DECL_PRESOLCOPY(presolCopyDualinfer)
2154 	{  /*lint --e{715}*/
2155 	   assert(scip != NULL);
2156 	   assert(presol != NULL);
2157 	   assert(strcmp(SCIPpresolGetName(presol), PRESOL_NAME) == 0);
2158 	
2159 	   /* call inclusion method of presolver */
2160 	   SCIP_CALL( SCIPincludePresolDualinfer(scip) );
2161 	
2162 	   return SCIP_OKAY;
2163 	}
2164 	
2165 	/** destructor of presolver to free user data (called when SCIP is exiting) */
2166 	static
2167 	SCIP_DECL_PRESOLFREE(presolFreeDualinfer)
2168 	{  /*lint --e{715}*/
2169 	   SCIP_PRESOLDATA* presoldata;
2170 	
2171 	   /* free presolver data */
2172 	   presoldata = SCIPpresolGetData(presol);
2173 	   assert(presoldata != NULL);
2174 	
2175 	   SCIPfreeBlockMemory(scip, &presoldata);
2176 	   SCIPpresolSetData(presol, NULL);
2177 	
2178 	   return SCIP_OKAY;
2179 	}
2180 	
2181 	/** execution method of presolver */
2182 	static
2183 	SCIP_DECL_PRESOLEXEC(presolExecDualinfer)
2184 	{  /*lint --e{715}*/
2185 	   SCIP_MATRIX* matrix;
2186 	   SCIP_Bool initialized;
2187 	   SCIP_Bool complete;
2188 	   SCIP_Bool infeasible;
2189 	   SCIP_PRESOLDATA* presoldata;
2190 	   FIXINGDIRECTION* varstofix;
2191 	   int npossiblefixings;
2192 	   int nconvarsfixed;
2193 	   int nintvarsfixed;
2194 	   int nbinvarsfixed;
2195 	   SIDECHANGE* sidestochange;
2196 	   int npossiblesidechanges;
2197 	   int nsideschanged;
2198 	   int i;
2199 	   int nrows;
2200 	   int ncols;
2201 	   SCIP_VAR* var;
2202 	
2203 	   assert(result != NULL);
2204 	   *result = SCIP_DIDNOTRUN;
2205 	
2206 	   if( SCIPgetNContVars(scip) + SCIPgetNImplVars(scip) == 0 )
2207 	      return SCIP_OKAY;
2208 	
2209 	   /* the reductions made in this presolver apply to all optimal solutions because of complementary slackness */
2210 	   if( !SCIPallowWeakDualReds(scip) )
2211 	      return SCIP_OKAY;
2212 	
2213 	   *result = SCIP_DIDNOTFIND;
2214 	
2215 	   if( SCIPisInfinity(scip, -SCIPgetPseudoObjval(scip)) )
2216 	   {
2217 	      SCIPdebugMsg(scip, "DualInfer not executed because condition of existing dual solution is not fulfilled.\n");
2218 	      return SCIP_OKAY;
2219 	   }
2220 	
2221 	   presoldata = SCIPpresolGetData(presol);
2222 	   assert(presoldata != NULL);
2223 	
2224 	   matrix = NULL;
2225 	   SCIP_CALL( SCIPmatrixCreate(scip, &matrix, TRUE, &initialized, &complete, &infeasible,
2226 	      naddconss, ndelconss, nchgcoefs, nchgbds, nfixedvars) );
2227 	
2228 	   /* if infeasibility was detected during matrix creation, return here */
2229 	   if( infeasible )
2230 	   {
2231 	      if( initialized )
2232 	         SCIPmatrixFree(scip, &matrix);
2233 	
2234 	      *result = SCIP_CUTOFF;
2235 	      return SCIP_OKAY;
2236 	   }
2237 	
2238 	   if( !initialized )
2239 	      return SCIP_OKAY;
2240 	
2241 	   npossiblefixings = 0;
2242 	   nconvarsfixed = 0;
2243 	   nintvarsfixed = 0;
2244 	   nbinvarsfixed = 0;
2245 	   npossiblesidechanges = 0;
2246 	   nsideschanged = 0;
2247 	
2248 	   nrows = SCIPmatrixGetNRows(matrix);
2249 	   ncols = SCIPmatrixGetNColumns(matrix);
2250 	
2251 	   SCIP_CALL( SCIPallocBufferArray(scip, &varstofix, ncols) );
2252 	   SCIP_CALL( SCIPallocBufferArray(scip, &sidestochange, nrows) );
2253 	
2254 	   BMSclearMemoryArray(varstofix, ncols);
2255 	   BMSclearMemoryArray(sidestochange, nrows);
2256 	
2257 	   SCIP_CALL( dualBoundStrengthening(scip, matrix, presoldata,
2258 	         varstofix, &npossiblefixings, sidestochange, &npossiblesidechanges) );
2259 	
2260 	   if( npossiblefixings > 0 )
2261 	   {
2262 	      for( i = ncols - 1; i >= 0; --i )
2263 	      {
2264 	         SCIP_Bool fixed;
2265 	
2266 	         var = SCIPmatrixGetVar(matrix, i);
2267 	
2268 	         /* there should be no fixings for variables with inconsistent locks */
2269 	         assert(varstofix[i] == NOFIX || (!SCIPmatrixUplockConflict(matrix, i) && !SCIPmatrixDownlockConflict(matrix, i)));
2270 	
2271 	         fixed = FALSE;
2272 	
2273 	         if( varstofix[i] == FIXATLB )
2274 	         {
2275 	            SCIP_Real lb;
2276 	            lb = SCIPvarGetLbLocal(var);
2277 	
2278 	            /* fix at lower bound */
2279 	            SCIP_CALL( SCIPfixVar(scip, var, lb, &infeasible, &fixed) );
2280 	            if( infeasible )
2281 	            {
2282 	               SCIPdebugMsg(scip, " -> infeasible fixing\n");
2283 	               *result = SCIP_CUTOFF;
2284 	               break;
2285 	            }
2286 	            assert(fixed);
2287 	            (*nfixedvars)++;
2288 	            *result = SCIP_SUCCESS;
2289 	         }
2290 	         else if( varstofix[i] == FIXATUB )
2291 	         {
2292 	            SCIP_Real ub;
2293 	            ub = SCIPvarGetUbLocal(var);
2294 	
2295 	            /* fix at upper bound */
2296 	            SCIP_CALL( SCIPfixVar(scip, var, ub, &infeasible, &fixed) );
2297 	            if( infeasible )
2298 	            {
2299 	               SCIPdebugMsg(scip, " -> infeasible fixing\n");
2300 	               *result = SCIP_CUTOFF;
2301 	               break;
2302 	            }
2303 	            assert(fixed);
2304 	            (*nfixedvars)++;
2305 	            *result = SCIP_SUCCESS;
2306 	         }
2307 	
2308 	         /* keep a small statistic which types of variables are fixed */
2309 	         if( fixed )
2310 	         {
2311 	            if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT )
2312 	               nconvarsfixed++;
2313 	            else if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY )
2314 	               nbinvarsfixed++;
2315 	            else
2316 	               nintvarsfixed++;
2317 	         }
2318 	      }
2319 	   }
2320 	
2321 	   if( npossiblesidechanges > 0 )
2322 	   {
2323 	      for( i = 0; i < nrows; i++ )
2324 	      {
2325 	         SCIP_CONS* cons;
2326 	         SCIP_CONSHDLR* conshdlr;
2327 	         const char* conshdlrname;
2328 	
2329 	         if( sidestochange[i] == NOCHANGE )
2330 	            continue;
2331 	
2332 	         if( presoldata->maxrowsupport < SCIPmatrixGetRowNNonzs(matrix, i) )
2333 	            continue;
2334 	
2335 	         cons = SCIPmatrixGetCons(matrix,i);
2336 	         conshdlr = SCIPconsGetHdlr(cons);
2337 	         conshdlrname = SCIPconshdlrGetName(conshdlr);
2338 	
2339 	         if( strcmp(conshdlrname, "linear") == 0 )
2340 	         {
2341 	            SCIP_Real lhs;
2342 	            SCIP_Real rhs;
2343 	            SCIP_Real matrixlhs;
2344 	            SCIP_Real matrixrhs;
2345 	
2346 	            lhs = SCIPgetLhsLinear(scip, cons);
2347 	            rhs = SCIPgetRhsLinear(scip, cons);
2348 	            matrixlhs = SCIPmatrixGetRowLhs(matrix, i);
2349 	            matrixrhs = SCIPmatrixGetRowRhs(matrix, i);
2350 	
2351 	            assert(!SCIPisEQ(scip, matrixlhs, matrixrhs));
2352 	
2353 	            /* when creating the matrix, constraints are multiplied if necessary by (-1)
2354 	               * to ensure that the following representation is obtained:
2355 	               * infty >= a x >= b
2356 	               * or
2357 	               * c >= ax >= b (ranged rows)
2358 	               */
2359 	
2360 	            /* for ranged constraints we have to distinguish between both sides */
2361 	            if( sidestochange[i] == RHSTOLHS )
2362 	            {
2363 	               if( SCIPisEQ(scip, matrixlhs, lhs) )
2364 	               {
2365 	                  /* change rhs to lhs */
2366 	                  SCIP_CALL( SCIPchgRhsLinear(scip, cons, matrixlhs) );
2367 	               }
2368 	               else
2369 	               {
2370 	                  /* consider multiplication by (-1) in the matrix */
2371 	                  SCIP_CALL( SCIPchgLhsLinear(scip, cons, -matrixlhs) );
2372 	               }
2373 	
2374 	               nsideschanged++;
2375 	               (*nchgsides)++;
2376 	            }
2377 	            else if( sidestochange[i] == LHSTORHS )
2378 	            {
2379 	               if( SCIPisEQ(scip, matrixrhs, rhs) )
2380 	               {
2381 	                  /* change lhs to rhs */
2382 	                  SCIP_CALL( SCIPchgLhsLinear(scip, cons, matrixrhs) );
2383 	               }
2384 	               else
2385 	               {
2386 	                  /* consider multiplication by (-1) in the matrix */
2387 	                  SCIP_CALL( SCIPchgRhsLinear(scip, cons, -matrixrhs) );
2388 	               }
2389 	
2390 	               nsideschanged++;
2391 	               (*nchgsides)++;
2392 	            }
2393 	         }
2394 	      }
2395 	   }
2396 	
2397 	   SCIPfreeBufferArray(scip, &sidestochange);
2398 	   SCIPfreeBufferArray(scip, &varstofix);
2399 	
2400 	   if( (nconvarsfixed + nintvarsfixed + nbinvarsfixed) > 0 || npossiblesidechanges > 0 )
2401 	   {
2402 	      SCIPdebugMsg(scip, "### fixed vars [cont: %d, int: %d, bin: %d], changed sides [%d]\n",
2403 	         nconvarsfixed, nintvarsfixed, nbinvarsfixed, nsideschanged);
2404 	   }
2405 	
2406 	   SCIPmatrixFree(scip, &matrix);
2407 	
2408 	   return SCIP_OKAY;
2409 	}
2410 	
2411 	
2412 	/*
2413 	 * presolver specific interface methods
2414 	 */
2415 	
2416 	/** creates the dual inference presolver and includes it in SCIP */
2417 	SCIP_RETCODE SCIPincludePresolDualinfer(
2418 	   SCIP*                 scip                /**< SCIP data structure */
2419 	   )
2420 	{
2421 	   SCIP_PRESOL* presol;
2422 	   SCIP_PRESOLDATA* presoldata;
2423 	
2424 	   /* create presolver data */
2425 	   SCIP_CALL( SCIPallocBlockMemory(scip, &presoldata) );
2426 	
2427 	   /* include presolver */
2428 	   SCIP_CALL( SCIPincludePresolBasic(scip, &presol, PRESOL_NAME, PRESOL_DESC, PRESOL_PRIORITY, PRESOL_MAXROUNDS,
2429 	         PRESOL_TIMING, presolExecDualinfer, presoldata) );
2430 	   SCIP_CALL( SCIPsetPresolCopy(scip, presol, presolCopyDualinfer) );
2431 	   SCIP_CALL( SCIPsetPresolFree(scip, presol, presolFreeDualinfer) );
2432 	
2433 	   SCIP_CALL( SCIPaddBoolParam(scip,
2434 	         "presolving/dualinfer/twocolcombine",
2435 	         "use convex combination of columns for determining dual bounds",
2436 	         &presoldata->usetwocolcombine, FALSE, DEFAULT_TWOCOLUMN_COMBINE, NULL, NULL) );
2437 	
2438 	   SCIP_CALL( SCIPaddIntParam(scip,
2439 	         "presolving/dualinfer/maxdualbndloops",
2440 	         "maximal number of dual bound strengthening loops",
2441 	         &presoldata->maxdualbndloops, FALSE, DEFAULT_MAXLOOPS_DUALBNDSTR, -1, INT_MAX, NULL, NULL) );
2442 	
2443 	   SCIP_CALL( SCIPaddIntParam(scip,
2444 	         "presolving/dualinfer/maxconsiderednonzeros",
2445 	         "maximal number of considered non-zeros within one column (-1: no limit)",
2446 	         &presoldata->maxconsiderednonzeros, TRUE, DEFAULT_MAXCONSIDEREDNONZEROS, -1, INT_MAX, NULL, NULL) );
2447 	
2448 	   SCIP_CALL( SCIPaddIntParam(scip,
2449 	         "presolving/dualinfer/maxretrievefails",
2450 	         "maximal number of consecutive useless hashtable retrieves",
2451 	         &presoldata->maxretrievefails, TRUE, DEFAULT_MAXRETRIEVEFAILS, -1, INT_MAX, NULL, NULL) );
2452 	
2453 	   SCIP_CALL( SCIPaddIntParam(scip,
2454 	         "presolving/dualinfer/maxcombinefails",
2455 	         "maximal number of consecutive useless column combines",
2456 	         &presoldata->maxcombinefails, TRUE, DEFAULT_MAXCOMBINEFAILS, -1, INT_MAX, NULL, NULL) );
2457 	
2458 	   SCIP_CALL( SCIPaddIntParam(scip,
2459 	         "presolving/dualinfer/maxhashfac",
2460 	         "Maximum number of hashlist entries as multiple of number of columns in the problem (-1: no limit)",
2461 	         &presoldata->maxhashfac, TRUE, DEFAULT_MAXHASHFAC, -1, INT_MAX, NULL, NULL) );
2462 	
2463 	   SCIP_CALL( SCIPaddIntParam(scip,
2464 	         "presolving/dualinfer/maxpairfac",
2465 	         "Maximum number of processed column pairs as multiple of the number of columns in the problem (-1: no limit)",
2466 	         &presoldata->maxpairfac, TRUE, DEFAULT_MAXPAIRFAC, -1, INT_MAX, NULL, NULL) );
2467 	
2468 	    SCIP_CALL( SCIPaddIntParam(scip,
2469 	         "presolving/dualinfer/maxrowsupport",
2470 	         "Maximum number of row's non-zeros for changing inequality to equality",
2471 	         &presoldata->maxrowsupport, FALSE, DEFAULT_MAXROWSUPPORT, 2, INT_MAX, NULL, NULL) );
2472 	
2473 	   return SCIP_OKAY;
2474 	}
2475