1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   heur_dps.c
26   	 * @ingroup DEFPLUGINS_HEUR
27   	 * @brief  dynamic partition search
28   	 * @author Katrin Halbig
29   	 *
30   	 * The dynamic partition search (DPS) is a construction heuristic which additionally needs a
31   	 * user decomposition with linking constraints only.
32   	 *
33   	 * This heuristic splits the problem into several sub-SCIPs according to the given decomposition. Thereby the linking constraints
34   	 * with their right-hand and left-hand sides are also split. DPS searches for a partition of the sides on the blocks
35   	 * so that a feasible solution is obtained.
36   	 * For each block the parts of the original linking constraints are extended by slack variables. Moreover, the objective function
37   	 * is replaced by the sum of these additional variables weighted by penalty parameters lambda. If all blocks have an optimal solution
38   	 * of zero, the algorithm terminates with a feasible solution for the main problem. Otherwise, the partition and the penalty parameters
39   	 * are updated, and the sub-SCIPs are solved again.
40   	 *
41   	 * A detailed description can be found in
42   	 * K. Halbig, A. Göß and D. Weninger (2023). Exploiting user-supplied Decompositions inside Heuristics. https://optimization-online.org/?p=23386
43   	 */
44   	
45   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
46   	
47   	#include "scip/heur_dps.h"
48   	#include "scip/pub_dcmp.h"
49   	#include "scip/pub_heur.h"
50   	#include "scip/pub_misc.h"
51   	#include "scip/scipdefplugins.h"
52   	#include "scip/scip_cons.h"
53   	#include "scip/scip_dcmp.h"
54   	#include "scip/scip_general.h"
55   	#include "scip/scip_heur.h"
56   	#include "scip/scip_mem.h"
57   	#include "scip/scip_message.h"
58   	#include "scip/scip_param.h"
59   	#include "scip/scip_prob.h"
60   	
61   	
62   	#define HEUR_NAME             "dps"
63   	#define HEUR_DESC             "primal heuristic for decomposable MIPs"
64   	#define HEUR_DISPCHAR         SCIP_HEURDISPCHAR_LNS
65   	#define HEUR_PRIORITY         75000
66   	#define HEUR_FREQ             -1
67   	#define HEUR_FREQOFS          0
68   	#define HEUR_MAXDEPTH         -1
69   	#define HEUR_TIMING           SCIP_HEURTIMING_BEFORENODE | SCIP_HEURTIMING_AFTERNODE
70   	#define HEUR_USESSUBSCIP      TRUE  /**< does the heuristic use a secondary SCIP instance? */
71   	
72   	#define DEFAULT_MAXIT         50    /**< maximum number of iterations */
73   	#define DEFAULT_PENALTY       100.0 /**< multiplier for absolute increase of penalty parameters */
74   	
75   	/* event handler properties */
76   	#define EVENTHDLR_NAME        "Dps"
77   	#define EVENTHDLR_DESC        "event handler for " HEUR_NAME " heuristic"
78   	
79   	/*
80   	 * Data structures
81   	 */
82   	
83   	/** primal heuristic data */
84   	struct SCIP_HeurData
85   	{
86   	   SCIP_CONS**           linkingconss;       /**< linking constraints */
87   	   int                   nlinking;           /**< number of linking constraints */
88   	   int                   nblocks;            /**< number of blocks */
89   	   int                   maxit;              /**< maximal number of iterations */
90   	   int                   timing;             /**< should the heuristic run before or after the processing of the node?
91   	                                                  (0: before, 1: after, 2: both)*/
92   	   SCIP_Real             maxlinkscore;       /**< maximal linking score of used decomposition (equivalent to percentage of linking constraints) */
93   	   SCIP_Real             penalty;            /**< multiplier for absolute increase of penalty parameters */
94   	   SCIP_Bool             reoptimize;         /**< should the problem get reoptimized with the original objective function? */
95   	   SCIP_Bool             reuse;              /**< should solutions get reused in subproblems? */
96   	   SCIP_Bool             reoptlimits;        /**< should strict limits for reoptimization be set? */
97   	};
98   	
99   	/** data related to one block */
100  	struct Blockproblem
101  	{
102  	   SCIP*                 blockscip;          /**< SCIP data structure */
103  	   SCIP_VAR**            slackvars;          /**< slack variables */
104  	   SCIP_CONS**           linkingconss;       /**< linking constraints */
105  	   int*                  linkingindices;     /**< indices of linking constraints in original problem */
106  	   int                   nlinking;           /**< number of linking constraints */
107  	   int                   nblockvars;         /**< number of block variables */
108  	   int                   nslackvars;         /**< number of slack variables */
109  	   SCIP_Real*            origobj;            /**< original objective coefficients */
110  	};
111  	typedef struct Blockproblem BLOCKPROBLEM;
112  	
113  	/** data related to one linking constraint */
114  	struct Linking
115  	{
116  	   SCIP_CONS*            linkingcons;        /**< corresponding linking constraint of original problem */
117  	   SCIP_CONS**           blockconss;         /**< linking constraints of the blocks */
118  	   SCIP_VAR**            slacks;             /**< slackvars of the blocks */
119  	   SCIP_Real*            minactivity;        /**< minimal activity of constraint for each block */
120  	   SCIP_Real*            maxactivity;        /**< maximal activity of constraint for each block */
121  	   SCIP_Real*            currentrhs;         /**< current partition of rhs */
122  	   SCIP_Real*            currentlhs;         /**< current partition of lhs */
123  	   int*                  blocknumbers;       /**< number of the blocks */
124  	   int                   nblocks;            /**< number of blocks in which this linking constraint participates; dimension of arrays */
125  	   int                   nslacks;            /**< number of slack variables */
126  	   int                   nslacksperblock;    /**< 2, if ranged constraint; 1, if only rhs or lhs */
127  	   int                   lastviolations;     /**< number of iterations in which the constraint was violated in succession */
128  	   SCIP_Bool             hasrhs;             /**< has linking constraint finite right-hand side? */
129  	   SCIP_Bool             haslhs;             /**< has linking constraint finite left-hand side? */
130  	};
131  	typedef struct Linking LINKING;
132  	
133  	/*
134  	 * Local methods
135  	 */
136  	
137  	/** assigns linking variables to last block
138  	 *
139  	 * The labels are copied to newdecomp and the linking variables are assigned to the last block (i.e., highest block label).
140  	 * Constraint labels and statistics are recomputed.
141  	 */
142  	static
143  	SCIP_RETCODE assignLinking(
144  	   SCIP*                 scip,               /**< SCIP data structure */
145  	   SCIP_DECOMP*          newdecomp,          /**< decomposition with assigned linking variables */
146  	   SCIP_VAR**            vars,               /**< sorted array of variables */
147  	   SCIP_CONS**           conss,              /**< sorted array of constraints */
148  	   int*                  varlabels,          /**< sorted array of variable labels */
149  	   int*                  conslabels,         /**< sorted array of constraint labels */
150  	   int                   nvars,              /**< number of variables */
151  	   int                   nconss,             /**< number of constraints */
152  	   int                   nlinkvars           /**< number of linking variables */
153  	   )
154  	{
155  	   int newlabel;
156  	   int v;
157  	
158  	   assert(scip != NULL);
159  	   assert(newdecomp != NULL);
160  	   assert(vars != NULL);
161  	   assert(conss != NULL);
162  	   assert(varlabels != NULL);
163  	   assert(conslabels != NULL);
164  	
165  	   /* copy the labels */
166  	   SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, vars, varlabels, nvars) );
167  	   SCIP_CALL( SCIPdecompSetConsLabels(newdecomp, conss, conslabels, nconss) );
168  	
169  	   /* assign linking variables */
170  	   newlabel = varlabels[nvars - 1]; /* take always label of last block */
171  	   assert(newlabel >= 0);
172  	   for( v = 0; v < nlinkvars; v++ )
173  	   {
174  	      SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, &vars[v], &newlabel, 1) );
175  	   }
176  	   SCIPdebugMsg(scip, "assigned %d linking variables\n", nlinkvars);
177  	
178  	   /* recompute constraint labels and statistics */
179  	   SCIP_CALL( SCIPcomputeDecompConsLabels(scip, newdecomp, conss, nconss) );
180  	   SCIP_CALL( SCIPcomputeDecompStats(scip, newdecomp, TRUE) );
181  	   nlinkvars = SCIPdecompGetNBorderVars(newdecomp);
182  	
183  	   /* get new labels and sort */
184  	   SCIPdecompGetConsLabels(newdecomp, conss, conslabels, nconss);
185  	   SCIPdecompGetVarsLabels(newdecomp, vars, varlabels, nvars);
186  	   SCIPsortIntPtr(conslabels, (void**)conss, nconss);
187  	   SCIPsortIntPtr(varlabels, (void**)vars, nvars);
188  	
189  	   /* After assigning the linking variables, blocks can have zero constraints.
190  	    * So the remaining variables are labeled as linking in SCIPcomputeDecompStats().
191  	    * We assign this variables to the same label as above.
192  	    */
193  	   if( nlinkvars >= 1 )
194  	   {
195  	      assert(varlabels[0] == SCIP_DECOMP_LINKVAR);
196  	      SCIPdebugMsg(scip, "assign again %d linking variables\n", nlinkvars);
197  	
198  	      for( v = 0; v < nlinkvars; v++ )
199  	      {
200  	         SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, &vars[v], &newlabel, 1) );
201  	      }
202  	      SCIP_CALL( SCIPcomputeDecompConsLabels(scip, newdecomp, conss, nconss) );
203  	      SCIP_CALL( SCIPcomputeDecompStats(scip, newdecomp, TRUE) );
204  	
205  	      SCIPdecompGetConsLabels(newdecomp, conss, conslabels, nconss);
206  	      SCIPdecompGetVarsLabels(newdecomp, vars, varlabels, nvars);
207  	      SCIPsortIntPtr(conslabels, (void**)conss, nconss);
208  	      SCIPsortIntPtr(varlabels, (void**)vars, nvars);
209  	   }
210  	   assert(varlabels[0] != SCIP_DECOMP_LINKVAR);
211  	
212  	   return SCIP_OKAY;
213  	}
214  	
215  	/** creates a sub-SCIP and sets parameters */
216  	static
217  	SCIP_RETCODE createSubscip(
218  	   SCIP*                 scip,               /**< main SCIP data structure */
219  	   SCIP**                subscip             /**< pointer to store created sub-SCIP */
220  	   )
221  	{
222  	   SCIP_Real infvalue;
223  	
224  	   assert(scip != NULL);
225  	   assert(subscip != NULL);
226  	
227  	   /* create a new SCIP instance */
228  	   SCIP_CALL( SCIPcreate(subscip) );
229  	   SCIP_CALL( SCIPincludeDefaultPlugins(*subscip) );
230  	
231  	   /* copy value for infinity */
232  	   SCIP_CALL( SCIPgetRealParam(scip, "numerics/infinity", &infvalue) );
233  	   SCIP_CALL( SCIPsetRealParam(*subscip, "numerics/infinity", infvalue) );
234  	
235  	   SCIP_CALL( SCIPcopyLimits(scip, *subscip) );
236  	
237  	   /* avoid recursive calls */
238  	   SCIP_CALL( SCIPsetSubscipsOff(*subscip, TRUE) );
239  	
240  	   /* disable cutting plane separation */
241  	   SCIP_CALL( SCIPsetSeparating(*subscip, SCIP_PARAMSETTING_OFF, TRUE) );
242  	
243  	   /* disable expensive presolving */
244  	   SCIP_CALL( SCIPsetPresolving(*subscip, SCIP_PARAMSETTING_FAST, TRUE) );
245  	
246  	   /* disable expensive techniques */
247  	   SCIP_CALL( SCIPsetIntParam(*subscip, "misc/usesymmetry", 0) );
248  	
249  	   /* do not abort subproblem on CTRL-C */
250  	   SCIP_CALL( SCIPsetBoolParam(*subscip, "misc/catchctrlc", FALSE) );
251  	
252  	#ifdef SCIP_DEBUG
253  	   /* for debugging, enable full output */
254  	   SCIP_CALL( SCIPsetIntParam(*subscip, "display/verblevel", 5) );
255  	   SCIP_CALL( SCIPsetIntParam(*subscip, "display/freq", 100000000) );
256  	#else
257  	   /* disable statistic timing inside sub SCIP and output to console */
258  	   SCIP_CALL( SCIPsetIntParam(*subscip, "display/verblevel", 0) );
259  	   SCIP_CALL( SCIPsetBoolParam(*subscip, "timing/statistictiming", FALSE) );
260  	#endif
261  	
262  	   /* speed up sub-SCIP by not checking dual LP feasibility */
263  	   SCIP_CALL( SCIPsetBoolParam(*subscip, "lp/checkdualfeas", FALSE) );
264  	
265  	   return SCIP_OKAY;
266  	}
267  	
268  	/** copies the given variables and constraints to the given sub-SCIP */
269  	static
270  	SCIP_RETCODE copyToSubscip(
271  	   SCIP*                 scip,               /**< source SCIP */
272  	   SCIP*                 subscip,            /**< target SCIP */
273  	   const char*           name,               /**< name for copied problem */
274  	   SCIP_VAR**            vars,               /**< array of variables to copy */
275  	   SCIP_CONS**           conss,              /**< array of constraints to copy */
276  	   SCIP_HASHMAP*         varsmap,            /**< hashmap for copied variables */
277  	   SCIP_HASHMAP*         conssmap,           /**< hashmap for copied constraints */
278  	   int                   nvars,              /**< number of variables to copy */
279  	   int                   nconss,             /**< number of constraints to copy */
280  	   SCIP_Bool*            success             /**< was copying successful? */
281  	   )
282  	{
283  	   SCIP_CONS* newcons;
284  	   SCIP_VAR* newvar;
285  	   int i;
286  	
287  	   assert(scip != NULL);
288  	   assert(subscip != NULL);
289  	   assert(vars != NULL);
290  	   assert(conss != NULL);
291  	   assert(varsmap != NULL);
292  	   assert(conssmap != NULL);
293  	   assert(success != NULL);
294  	
295  	   SCIPdebugMsg(scip, "copyToSubscip\n");
296  	
297  	   /* create problem in sub-SCIP */
298  	   SCIP_CALL( SCIPcreateProb(subscip, name, NULL, NULL, NULL, NULL, NULL, NULL, NULL) );
299  	
300  	   /* copy variables */
301  	   for( i = 0; i < nvars; ++i )
302  	   {
303  	      SCIP_CALL( SCIPgetVarCopy(scip, subscip, vars[i], &newvar, varsmap, conssmap, FALSE, success) );
304  	
305  	      /* abort if variable was not successfully copied */
306  	      if( !(*success) )
307  	      {
308  	         SCIPwarningMessage(scip, "Abort heuristic dps since not all variables were successfully copied.\n");
309  	         SCIPABORT();
310  	         return SCIP_OKAY;
311  	      }
312  	   }
313  	   assert(nvars == SCIPgetNOrigVars(subscip));
314  	
315  	   /* copy constraints */
316  	   for( i = 0; i < nconss; ++i )
317  	   {
318  	      assert(conss[i] != NULL);
319  	      assert(!SCIPconsIsModifiable(conss[i]));
320  	      assert(SCIPconsIsActive(conss[i]));
321  	      assert(!SCIPconsIsDeleted(conss[i]));
322  	
323  	      /* copy the constraint */
324  	      SCIP_CALL( SCIPgetConsCopy(scip, subscip, conss[i], &newcons, SCIPconsGetHdlr(conss[i]), varsmap, conssmap, NULL,
325  	            SCIPconsIsInitial(conss[i]), SCIPconsIsSeparated(conss[i]), SCIPconsIsEnforced(conss[i]),
326  	            SCIPconsIsChecked(conss[i]), SCIPconsIsPropagated(conss[i]), FALSE, FALSE,
327  	            SCIPconsIsDynamic(conss[i]), SCIPconsIsRemovable(conss[i]), FALSE, FALSE, success) );
328  	
329  	      /* abort if constraint was not successfully copied */
330  	      if( !(*success) )
331  	         return SCIP_OKAY;
332  	
333  	      SCIP_CALL( SCIPaddCons(subscip, newcons) );
334  	      SCIP_CALL( SCIPreleaseCons(subscip, &newcons) );
335  	   }
336  	
337  	   /* block constraint contains variables which are not part of this block
338  	    *
339  	    * todo: maybe they are part of the block, but it is not recognized, because they are, for example, negated or aggregated.
340  	    */
341  	   if( nvars != SCIPgetNOrigVars(subscip) )
342  	      *success = FALSE;
343  	
344  	   return SCIP_OKAY;
345  	}
346  	
347  	/** creates the subscip for a given block */
348  	static
349  	SCIP_RETCODE createBlockproblem(
350  	   SCIP*                 scip,               /**< SCIP data structure */
351  	   BLOCKPROBLEM*         blockproblem,       /**< blockproblem that should be created */
352  	   LINKING**             linkings,           /**< linkings that will be (partially) initialized */
353  	   SCIP_CONS**           conss,              /**< sorted array of constraints of this block */
354  	   SCIP_VAR**            vars,               /**< sorted array of variables of this block */
355  	   int                   nconss,             /**< number of constraints of this block */
356  	   int                   nvars,              /**< number of variables of this block */
357  	   SCIP_CONS**           linkingconss,       /**< linking constraints in the original problem */
358  	   int                   nlinking,           /**< number of linking constraints in the original problem */
359  	   int                   blocknumber,        /**< number of block that should be created */
360  	   SCIP_Bool*            success             /**< pointer to store whether creation was successful */
361  	   )
362  	{
363  	   char name[SCIP_MAXSTRLEN];
364  	   SCIP_HASHMAP* varsmap;
365  	   SCIP_HASHMAP* conssmap;
366  	   SCIP_VAR** consvars; /* all vars in original linking cons */
367  	   SCIP_Real* consvals;
368  	   int nconsvars;
369  	   SCIP_VAR** blockvars; /* vars of current linking cons of current block */
370  	   SCIP_Real* blockvals;
371  	   int nblockvars;
372  	   SCIP_VAR** subvars; /* all vars of subscip */
373  	   int maxnconsvars; /* current size of arrays */
374  	   int c;
375  	   int v;
376  	
377  	   assert(scip != NULL);
378  	   assert(blockproblem != NULL);
379  	   assert(conss != NULL);
380  	   assert(vars != NULL);
381  	   assert(blockproblem->blockscip != NULL);
382  	
383  	   maxnconsvars = 20; /* start size; increase size if necessary */
384  	
385  	   SCIPdebugMsg(scip, "Create blockproblem %d\n", blocknumber);
386  	
387  	   /* create the variable/constraint mapping hash map */
388  	   SCIP_CALL( SCIPhashmapCreate(&varsmap, SCIPblkmem(scip), nvars) );
389  	   SCIP_CALL( SCIPhashmapCreate(&conssmap, SCIPblkmem(scip), nconss) );
390  	
391  	   /* get name of the original problem and add "comp_nr" */
392  	   (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_comp_%d", SCIPgetProbName(scip), blocknumber);
393  	
394  	   SCIP_CALL( copyToSubscip(scip, blockproblem->blockscip, name, vars, conss, varsmap, conssmap,
395  	                            nvars, nconss, success) );
396  	   if( !(*success) )
397  	   {
398  	      SCIPdebugMsg(scip, "Copy to subscip failed\n");
399  	      SCIPhashmapFree(&conssmap);
400  	      SCIPhashmapFree(&varsmap);
401  	
402  	      return SCIP_OKAY;
403  	   }
404  	
405  	   /* save number of variables that have a corresponding variable in original problem*/
406  	   blockproblem->nblockvars = SCIPgetNVars(blockproblem->blockscip);
407  	
408  	   /* save original objective and set objective to zero */
409  	   subvars = SCIPgetVars(blockproblem->blockscip);
410  	   for( v = 0; v < nvars; v++ )
411  	   {
412  	      blockproblem->origobj[v] = SCIPvarGetObj(subvars[v]);
413  	      SCIP_CALL( SCIPchgVarObj(blockproblem->blockscip, subvars[v], 0.0) );
414  	   }
415  	
416  	   /* allocate memory */
417  	   SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &blockvars, nvars + 2) ); /* two entries for the slack variables */
418  	   SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &blockvals, nvars + 2) );
419  	   SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &consvars, maxnconsvars) );
420  	   SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &consvals, maxnconsvars) );
421  	
422  	   /* find and add parts of linking constraints */
423  	   SCIPdebugMsg(scip, "add parts of linking constraints\n");
424  	   for( c = 0; c < nlinking; c++ )
425  	   {
426  	      const char* conshdlrname;
427  	      char consname[SCIP_MAXSTRLEN];
428  	      SCIP_CONS* newcons;
429  	      SCIP_Real rhs;
430  	      SCIP_Real lhs;
431  	      SCIP_Real minact;
432  	      SCIP_Real maxact;
433  	      SCIP_Bool mininfinite;
434  	      SCIP_Bool maxinfinite;
435  	
436  	      assert(linkingconss[c] != NULL);
437  	
438  	      newcons = NULL;
439  	
440  	#ifdef SCIP_MORE_DEBUG
441  	      SCIPdebugMsg(scip, "consider constraint %s\n", SCIPconsGetName(linkingconss[c]));
442  	      SCIPdebugPrintCons(scip, linkingconss[c], NULL);
443  	#endif
444  	
445  	      nblockvars = 0;
446  	
447  	      /* every constraint with linear representation is allowed */
448  	      conshdlrname = SCIPconshdlrGetName(SCIPconsGetHdlr(linkingconss[c]));
449  	      if( !( (strcmp(conshdlrname, "linear") == 0) || (strcmp(conshdlrname, "setppc") == 0)
450  	            || (strcmp(conshdlrname, "logicor") == 0) || (strcmp(conshdlrname, "knapsack") == 0)
451  	            || (strcmp(conshdlrname, "varbound") == 0) ) )
452  	      {
453  	         SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "Heuristic %s cannot handle linking constraints of type %s\n", HEUR_NAME, conshdlrname);
454  	         /* TODO which other types can we handle/transform in a linear constraint? */
455  	
456  	         *success = FALSE;
457  	         break; /* releases memory and breaks heuristic */
458  	      }
459  	
460  	      SCIP_CALL( SCIPgetConsNVars(scip, linkingconss[c], &nconsvars, success) );
461  	
462  	      /* reallocate memory if we have more variables than maxnconsvars */
463  	      if( nconsvars > maxnconsvars )
464  	      {
465  	         int newsize;
466  	
467  	         /* calculate new size */
468  	         newsize = SCIPcalcMemGrowSize(scip, MAX(2 * maxnconsvars, nconsvars)); /* at least double size */
469  	
470  	         SCIP_CALL( SCIPreallocBufferArray(blockproblem->blockscip, &consvars, newsize) );
471  	         SCIP_CALL( SCIPreallocBufferArray(blockproblem->blockscip, &consvals, newsize) );
472  	         maxnconsvars = newsize;
473  	      }
474  	
475  	      SCIP_CALL( SCIPgetConsVars(scip, linkingconss[c], consvars, nconsvars, success) );
476  	      SCIP_CALL( SCIPgetConsVals(scip, linkingconss[c], consvals, nconsvars, success) );
477  	
478  	      if( !(*success) )
479  	      {
480  	         SCIPdebugMsg(scip, "Create blockproblem failed\n");
481  	         break; /* releases memory and breaks heuristic */
482  	      }
483  	
484  	      /* check if constraint contains variables of this block */
485  	      for( v = 0; v < nconsvars; v++ )
486  	      {
487  	         if( SCIPhashmapExists(varsmap, (void*)consvars[v]) )
488  	         {
489  	            blockvars[nblockvars] = (SCIP_VAR*) SCIPhashmapGetImage(varsmap, (void*)consvars[v]);
490  	            blockvals[nblockvars] = consvals[v];
491  	            ++nblockvars;
492  	         }
493  	         /* handle negated variables*/
494  	         else if( SCIPvarGetStatus(consvars[v]) == SCIP_VARSTATUS_NEGATED)
495  	         {
496  	            if( SCIPhashmapExists(varsmap, (void*)SCIPvarGetNegationVar(consvars[v])) ) /* negation exists in this block */
497  	            {
498  	               /* save negated variable */
499  	               SCIP_VAR* origblockvar = (SCIP_VAR*) SCIPhashmapGetImage(varsmap, (void*)SCIPvarGetNegationVar(consvars[v]));
500  	               SCIP_VAR* negblockvar = NULL;
501  	               SCIP_CALL( SCIPgetNegatedVar(blockproblem->blockscip, origblockvar, &negblockvar) );
502  	               blockvars[nblockvars] = negblockvar;
503  	               blockvals[nblockvars] = consvals[v];
504  	               ++nblockvars;
505  	            }
506  	         }
507  	      }
508  	
509  	      /* continue with next linking constraint if it has no part in current block */
510  	      if( nblockvars == 0 )
511  	         continue;
512  	
513  	      /* get rhs and/or lhs */
514  	      rhs = SCIPconsGetRhs(scip, linkingconss[c], success);
515  	      if( !(*success) )
516  	      {
517  	         SCIPdebugMsg(scip, "Create blockproblem failed\n");
518  	         return SCIP_OKAY;
519  	      }
520  	      lhs = SCIPconsGetLhs(scip, linkingconss[c], success);
521  	      if( !(*success) )
522  	      {
523  	         SCIPdebugMsg(scip, "Create blockproblem failed\n");
524  	         return SCIP_OKAY;
525  	      }
526  	      assert(!SCIPisInfinity(scip, rhs) || !SCIPisInfinity(scip, -lhs)); /* at least one side bounded */
527  	      assert(SCIPisLE(scip, lhs, rhs));
528  	
529  	      if( !SCIPisInfinity(scip, rhs) )
530  	         linkings[c]->hasrhs = TRUE;
531  	      if( !SCIPisInfinity(scip, -lhs) )
532  	         linkings[c]->haslhs = TRUE;
533  	      if( !SCIPisInfinity(scip, rhs) && !SCIPisInfinity(scip, -lhs))
534  	         linkings[c]->nslacksperblock = 2;
535  	      else
536  	         linkings[c]->nslacksperblock = 1;
537  	
538  	      /* add slack variable for rhs */
539  	      if( linkings[c]->hasrhs )
540  	      {
541  	         /* slack variable z_r >= 0 */
542  	         char varname[SCIP_MAXSTRLEN];
543  	         (void)SCIPsnprintf(varname, SCIP_MAXSTRLEN, "z_r_%s", SCIPconsGetName(linkingconss[c]));
544  	         SCIP_CALL( SCIPcreateVarBasic(blockproblem->blockscip, &blockvars[nblockvars], varname,
545  	                                          0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) );
546  	         blockvals[nblockvars] = -1.0;
547  	         SCIP_CALL( SCIPaddVar(blockproblem->blockscip, blockvars[nblockvars]) );
548  	#ifdef SCIP_MORE_DEBUG
549  	         SCIPdebugMsg(scip, "Add variable %s\n", SCIPvarGetName(blockvars[nblockvars]));
550  	#endif
551  	         linkings[c]->slacks[linkings[c]->nslacks] = blockvars[nblockvars];
552  	         blockproblem->slackvars[blockproblem->nslackvars] = blockvars[nblockvars];
553  	         ++blockproblem->nslackvars;
554  	         ++linkings[c]->nslacks;
555  	         ++nblockvars;
556  	      }
557  	
558  	      /* add slack variable for lhs */
559  	      if( linkings[c]->haslhs )
560  	      {
561  	         /* slack variable z_l >= 0 */
562  	         char varname[SCIP_MAXSTRLEN];
563  	         (void)SCIPsnprintf(varname, SCIP_MAXSTRLEN, "z_l_%s", SCIPconsGetName(linkingconss[c]));
564  	         SCIP_CALL( SCIPcreateVarBasic(blockproblem->blockscip, &blockvars[nblockvars], varname,
565  	                                          0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) );
566  	         blockvals[nblockvars] = 1.0;
567  	         SCIP_CALL( SCIPaddVar(blockproblem->blockscip, blockvars[nblockvars]) );
568  	#ifdef SCIP_MORE_DEBUG
569  	         SCIPdebugMsg(scip, "Add variable %s\n", SCIPvarGetName(blockvars[nblockvars]));
570  	#endif
571  	         linkings[c]->slacks[linkings[c]->nslacks] = blockvars[nblockvars];
572  	         blockproblem->slackvars[blockproblem->nslackvars] = blockvars[nblockvars];
573  	         ++blockproblem->nslackvars;
574  	         ++linkings[c]->nslacks;
575  	         ++nblockvars;
576  	      }
577  	
578  	      /* add linking constraint with slack variable */
579  	      (void)SCIPsnprintf(consname, SCIP_MAXSTRLEN, "%s", SCIPconsGetName(linkingconss[c]));
580  	      SCIP_CALL( SCIPcreateConsBasicLinear(blockproblem->blockscip, &newcons, consname, nblockvars, blockvars, blockvals, lhs, rhs) );
581  	      SCIP_CALL( SCIPaddCons(blockproblem->blockscip, newcons) );
582  	#ifdef SCIP_MORE_DEBUG
583  	      SCIPdebugMsg(blockproblem->blockscip, "add constraint %s\n", SCIPconsGetName(newcons));
584  	      SCIPdebugPrintCons(blockproblem->blockscip, newcons, NULL);
585  	#endif
586  	
587  	      blockproblem->linkingconss[blockproblem->nlinking] = newcons;
588  	      linkings[c]->blockconss[linkings[c]->nblocks] = newcons;
589  	      linkings[c]->blocknumbers[linkings[c]->nblocks] = blocknumber;
590  	      blockproblem->linkingindices[blockproblem->nlinking] = c;
591  	
592  	      /* calculate minimal und maximal activity (exclude slack variables) */
593  	      minact = 0;
594  	      maxact = 0;
595  	      mininfinite = FALSE;
596  	      maxinfinite = FALSE;
597  	      for( v = 0; v < nblockvars - linkings[c]->nslacksperblock && (!mininfinite || !maxinfinite); v++ )
598  	      {
599  	         SCIP_Real lb;
600  	         SCIP_Real ub;
601  	         lb = SCIPvarGetLbGlobal(blockvars[v]);
602  	         ub = SCIPvarGetUbGlobal(blockvars[v]);
603  	
604  	         if( blockvals[v] >= 0.0 )
605  	         {
606  	            mininfinite = (mininfinite || SCIPisInfinity(scip, -lb));
607  	            maxinfinite = (maxinfinite || SCIPisInfinity(scip, ub));
608  	            if( !mininfinite )
609  	               minact += blockvals[v] * lb;
610  	            if( !maxinfinite )
611  	               maxact += blockvals[v] * ub;
612  	         }
613  	         else
614  	         {
615  	            mininfinite = (mininfinite || SCIPisInfinity(scip, ub));
616  	            maxinfinite = (maxinfinite || SCIPisInfinity(scip, -lb));
617  	            if( !mininfinite )
618  	               minact += blockvals[v] * ub;
619  	            if( !maxinfinite )
620  	               maxact += blockvals[v] * lb;
621  	         }
622  	      }
623  	
624  	      if( mininfinite )
625  	         linkings[c]->minactivity[linkings[c]->nblocks] = -SCIPinfinity(scip);
626  	      else
627  	         linkings[c]->minactivity[linkings[c]->nblocks] = minact;
628  	      if( maxinfinite )
629  	         linkings[c]->maxactivity[linkings[c]->nblocks] = SCIPinfinity(scip);
630  	      else
631  	         linkings[c]->maxactivity[linkings[c]->nblocks] = maxact;
632  	      assert(SCIPisLE(scip, linkings[c]->minactivity[linkings[c]->nblocks], linkings[c]->maxactivity[linkings[c]->nblocks]));
633  	
634  	      linkings[c]->nblocks++;
635  	      blockproblem->nlinking++;
636  	
637  	      for( v = 1; v <= linkings[c]->nslacksperblock; v++ )
638  	      {
639  	         SCIP_CALL( SCIPreleaseVar(blockproblem->blockscip, &blockvars[nblockvars - v]) );
640  	      }
641  	
642  	      SCIP_CALL( SCIPreleaseCons(blockproblem->blockscip, &newcons) );
643  	   }
644  	   assert(blockproblem->nlinking <= nlinking);
645  	
646  	   /* free memory */
647  	   SCIPfreeBufferArray(blockproblem->blockscip, &consvals);
648  	   SCIPfreeBufferArray(blockproblem->blockscip, &consvars);
649  	   SCIPfreeBufferArray(blockproblem->blockscip, &blockvals);
650  	   SCIPfreeBufferArray(blockproblem->blockscip, &blockvars);
651  	
652  	   SCIPhashmapFree(&conssmap);
653  	   SCIPhashmapFree(&varsmap);
654  	
655  	   return SCIP_OKAY;
656  	}
657  	
658  	/** creates data structures and splits problem into blocks */
659  	static
660  	SCIP_RETCODE createAndSplitProblem(
661  	   SCIP*                 scip,               /**< SCIP data structure */
662  	   SCIP_HEURDATA*        heurdata,           /**< heuristic data */
663  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
664  	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
665  	   LINKING**             linkings,           /**< array of linking data structures */
666  	   SCIP_VAR**            vars,               /**< sorted array of variables */
667  	   SCIP_CONS**           conss,              /**< sorted array of constraints */
668  	   SCIP_Bool*            success             /**< pointer to store whether splitting was successful */
669  	   )
670  	{
671  	   int* nconssblock;
672  	   int* nvarsblock;
673  	   int conssoffset;
674  	   int varsoffset;
675  	   int i;   /* blocknumber */
676  	
677  	   assert(scip != NULL);
678  	   assert(heurdata != NULL);
679  	   assert(vars != NULL);
680  	   assert(conss != NULL);
681  	
682  	   SCIP_CALL( SCIPallocBufferArray(scip, &nvarsblock, heurdata->nblocks + 1) );
683  	   SCIP_CALL( SCIPallocBufferArray(scip, &nconssblock, heurdata->nblocks + 1) );
684  	   SCIP_CALL( SCIPdecompGetVarsSize(decomp, nvarsblock, heurdata->nblocks + 1) );
685  	   SCIP_CALL( SCIPdecompGetConssSize(decomp, nconssblock, heurdata->nblocks + 1) );
686  	   assert(0 == nvarsblock[0]);
687  	
688  	   varsoffset = 0;
689  	   conssoffset = 0;
690  	
691  	   for( i = 0; i < heurdata->nblocks; i++)
692  	   {
693  	      conssoffset += nconssblock[i];
694  	      varsoffset += nvarsblock[i];
695  	
696  	      SCIP_CALL( createBlockproblem(scip, blockproblem[i], linkings, &conss[conssoffset], &vars[varsoffset], nconssblock[i+1], nvarsblock[i+1],
697  	                                    heurdata->linkingconss, heurdata->nlinking, i, success) );
698  	      if( !(*success) )
699  	         break;
700  	   }
701  	
702  	   SCIPfreeBufferArray(scip, &nconssblock);
703  	   SCIPfreeBufferArray(scip, &nvarsblock);
704  	
705  	   return SCIP_OKAY;
706  	}
707  	
708  	/** rounds partition for one linking constraint to integer value if variables and coefficients are integer
709  	 *
710  	 *  changes only currentrhs/currentlhs
711  	 */
712  	static
713  	SCIP_RETCODE roundPartition(
714  	   SCIP*                 scip,               /**< SCIP data structure */
715  	   LINKING*              linking,            /**< one linking data structure */
716  	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
717  	   SCIP_Bool             roundbyrhs          /**< round by right hand side? */
718  	   )
719  	{
720  	   SCIP_Real* fracPart;
721  	   int* sorting;
722  	   int* isinteger;
723  	   SCIP_Real sumbefor; /* includes value at idx */
724  	   SCIP_Real sumafter;
725  	   SCIP_Real diff;
726  	   int nnonintblocks; /* number of non integer blocks */
727  	   int idx;
728  	   int b;
729  	   int i;
730  	   int k;
731  	
732  	   assert(scip != NULL);
733  	   assert(linking != NULL);
734  	   assert(blockproblem != NULL);
735  	
736  	   nnonintblocks = 0;
737  	   idx = 0;
738  	
739  	   SCIP_CALL( SCIPallocBufferArray(scip, &fracPart, linking->nblocks) );
740  	   SCIP_CALL( SCIPallocBufferArray(scip, &sorting, linking->nblocks) );
741  	   SCIP_CALL( SCIPallocBufferArray(scip, &isinteger, linking->nblocks) );
742  	
743  	   /* get integer blocks and fractional parts */
744  	   for( b = 0; b < linking->nblocks; b++ )
745  	   {
746  	      SCIP* subscip;
747  	      SCIP_CONS* blockcons;
748  	      SCIP_VAR** blockvars;
749  	      SCIP_Real* blockvals;
750  	      int nblockvars;
751  	      int length; /* number of block variables without slack variables */
752  	      SCIP_Bool success;
753  	
754  	      subscip = blockproblem[linking->blocknumbers[b]]->blockscip;
755  	      blockcons = linking->blockconss[b];
756  	      sorting[b] = b; /* store current sorting to sort back */
757  	
758  	      SCIP_CALL( SCIPgetConsNVars(subscip, blockcons, &nblockvars, &success) );
759  	      assert(success);
760  	      SCIP_CALL( SCIPallocBufferArray(scip, &blockvars, nblockvars) );
761  	      SCIP_CALL( SCIPallocBufferArray(scip, &blockvals, nblockvars) );
762  	
763  	      SCIP_CALL( SCIPgetConsVars(subscip, blockcons, blockvars, nblockvars, &success) );
764  	      assert(success);
765  	      SCIP_CALL( SCIPgetConsVals(subscip, blockcons, blockvals, nblockvars, &success) );
766  	      assert(success);
767  	
768  	      /* get number of block variables in this constraint without slack variables */
769  	      length = nblockvars - linking->nslacksperblock;
770  	
771  	      /* get blocks with integer value */
772  	      isinteger[b] = 1;
773  	      for( i = 0; i < length; i++ )
774  	      {
775  	         if( SCIPvarGetType(blockvars[i]) == SCIP_VARTYPE_CONTINUOUS ||  !SCIPisIntegral(scip, blockvals[i]) )
776  	         {
777  	            isinteger[b] = 0;
778  	            nnonintblocks++;
779  	            break;
780  	         }
781  	      }
782  	
783  	      /* get fractional part of blockconstraints */
784  	      if( roundbyrhs )
785  	         fracPart[b] = linking->currentrhs[b] - floor(linking->currentrhs[b]); /* do not use SCIPfrac()! */
786  	      else
787  	         fracPart[b] = linking->currentlhs[b] - floor(linking->currentlhs[b]);
788  	
789  	      SCIPfreeBufferArray(scip, &blockvals);
790  	      SCIPfreeBufferArray(scip, &blockvars);
791  	   }
792  	
793  	   /* sort non-integer blocks to the front */
794  	   SCIPsortIntIntReal(isinteger, sorting, fracPart, linking->nblocks);
795  	
796  	   /* sort by fractional part */
797  	   SCIPsortRealInt(fracPart, sorting, nnonintblocks);
798  	   SCIPsortRealInt(&fracPart[nnonintblocks], &sorting[nnonintblocks], linking->nblocks - nnonintblocks);
799  	
800  	   /* detect blocks for rounding down and rounding up:
801  	    * integer blocks with small fractional parts are rounded down
802  	    * integer blocks with big fractional parts are rounded up
803  	    */
804  	
805  	   sumbefor = 0.0;
806  	   sumafter = 0.0;
807  	
808  	   for( i = 0; i < linking->nblocks - nnonintblocks; i++ )
809  	      sumafter += 1 - fracPart[nnonintblocks + i];
810  	
811  	   for( i = 0; i < linking->nblocks - nnonintblocks; i++ )
812  	   {
813  	      sumbefor += fracPart[nnonintblocks + i];
814  	      sumafter -= 1 - fracPart[nnonintblocks + i];
815  	
816  	      if( sumbefor >= sumafter )
817  	      {
818  	         for( k = 0; k <= i; k++ )
819  	            fracPart[nnonintblocks + k] = -fracPart[nnonintblocks + k];
820  	
821  	         for( k = i + 1; k < linking->nblocks - nnonintblocks; k++ )
822  	            fracPart[nnonintblocks + k] = 1 - fracPart[nnonintblocks + k];
823  	
824  	         idx = i;
825  	         break;
826  	      }
827  	   }
828  	   diff = sumbefor - sumafter;
829  	   assert(SCIPisGE(scip, diff, 0.0));
830  	
831  	   /* add difference to last non integer block */
832  	   for( i = nnonintblocks - 1; i >= 0; i-- )
833  	   {
834  	      if( SCIPisGT(scip, diff, 0.0) )
835  	      {
836  	         fracPart[i] = diff;
837  	         diff = 0;
838  	      }
839  	      else
840  	         fracPart[i] = 0;
841  	   }
842  	
843  	   /* add difference to last rounded down block if no non integer block exists */
844  	   if( SCIPisGT(scip, diff, 0.0))
845  	   {
846  	      assert(nnonintblocks == 0);
847  	      fracPart[idx] += diff;
848  	   }
849  	
850  	   /* sort back */
851  	   SCIPsortIntReal(sorting, fracPart, linking->nblocks);
852  	
853  	   /* round partition
854  	    * if we have a ranged constraint, both sides get rounded in the same way
855  	    */
856  	   for( b = 0; b < linking->nblocks; b++ )
857  	   {
858  	      if( linking->hasrhs )
859  	         linking->currentrhs[b] += fracPart[b];
860  	
861  	      if( linking->haslhs )
862  	         linking->currentlhs[b] += fracPart[b];
863  	   }
864  	
865  	   SCIPfreeBufferArray(scip, &isinteger);
866  	   SCIPfreeBufferArray(scip, &sorting);
867  	   SCIPfreeBufferArray(scip, &fracPart);
868  	
869  	   return SCIP_OKAY;
870  	}
871  	
872  	/** calculates initial partition and sets rhs/lhs in blockproblems */
873  	static
874  	SCIP_RETCODE initCurrent(
875  	   SCIP*                 scip,               /**< SCIP data structure of main scip */
876  	   LINKING**             linkings,           /**< array of linking data structures */
877  	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
878  	   SCIP_HEURTIMING       heurtiming,         /**< current point in the node solving process */
879  	   int                   nlinking,           /**< number of linking constraints */
880  	   SCIP_Bool*            success             /**< pointer to store whether initialization was successful */
881  	   )
882  	{
883  	   LINKING* linking;
884  	   SCIP_Real rhs;
885  	   SCIP_Real lhs;
886  	   SCIP_Real residualrhs;
887  	   SCIP_Real residuallhs;
888  	   SCIP_Real goalvalue;
889  	   int c;
890  	   int b;
891  	
892  	   assert(scip != NULL);
893  	   assert(linkings != NULL);
894  	   assert(blockproblem != NULL);
895  	   assert(nlinking > 0);
896  	
897  	   SCIPdebugMsg(scip, "initialize partition\n");
898  	
899  	   /* for each linking constraint the rhs/lhs is split between the blocks */
900  	   for( c = 0; c < nlinking; c++ )
901  	   {
902  	      linking = linkings[c];
903  	      rhs = SCIPconsGetRhs(scip, linking->linkingcons, success);
904  	      assert(*success);
905  	      lhs = SCIPconsGetLhs(scip, linking->linkingcons, success);
906  	      assert(*success);
907  	      residualrhs = rhs;
908  	      residuallhs = lhs;
909  	
910  	      /* LP value for each block plus part of remaining amount if sum is not equal to rhs/lhs */
911  	      if( (heurtiming & SCIP_HEURTIMING_AFTERNODE) && (linking->hasrhs || linking->haslhs) )
912  	      {
913  	         SCIP_Real sumrhs = 0.0;
914  	         SCIP_Real sumlhs = 0.0;
915  	         for( b = 0; b < linking->nblocks; b++ )
916  	         {
917  	            SCIP_VAR** consvars;
918  	            SCIP_Real* consvals;
919  	            SCIP_CONS* cons;
920  	            int nconsvars;
921  	            SCIP_Real lpvalue;
922  	            int i;
923  	
924  	            /* get variables of linking constraint of one block */
925  	            cons = linking->blockconss[b];
926  	            SCIP_CALL( SCIPgetConsNVars(blockproblem[linking->blocknumbers[b]]->blockscip, cons, &nconsvars, success) );
927  	            SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nconsvars) );
928  	            SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nconsvars) );
929  	            SCIP_CALL( SCIPgetConsVars(blockproblem[linking->blocknumbers[b]]->blockscip, cons, consvars, nconsvars, success) );
930  	            SCIP_CALL( SCIPgetConsVals(blockproblem[linking->blocknumbers[b]]->blockscip, cons, consvals, nconsvars, success) );
931  	            assert(success);
932  	
933  	            /* calculate value of partition variable in lp solution */
934  	            lpvalue = 0.0;
935  	            for( i = 0; i < nconsvars - linking->nslacksperblock; i++ )
936  	            {
937  	               SCIP_VAR* origvar;
938  	               SCIP_Real varlpvalue;
939  	
940  	               origvar = SCIPfindVar(scip, SCIPvarGetName(consvars[i]));
941  	               if( origvar == NULL ) /* e.g. variable is negated */
942  	               {
943  	                  *success = FALSE;
944  	                  return SCIP_OKAY;
945  	               }
946  	               varlpvalue = SCIPvarGetLPSol(origvar);
947  	               lpvalue += varlpvalue * consvals[i];
948  	            }
949  	            sumrhs += lpvalue;
950  	            sumlhs += lpvalue;
951  	
952  	            /* set currentrhs/lhs at lp value */
953  	            if( linking->hasrhs )
954  	               linking->currentrhs[b] = lpvalue;
955  	            if( linking->haslhs )
956  	               linking->currentlhs[b] = lpvalue;
957  	
958  	            SCIPfreeBufferArray(scip, &consvars);
959  	            SCIPfreeBufferArray(scip, &consvals);
960  	         }
961  	         assert(SCIPisLE(scip, sumrhs, rhs));
962  	         assert(SCIPisGE(scip, sumlhs, lhs));
963  	
964  	         /* distribute remaining amount */
965  	         if( !SCIPisEQ(scip, rhs, sumrhs) && linking->hasrhs )
966  	         {
967  	            SCIP_Real diff;
968  	            SCIP_Real part;
969  	            SCIP_Real residual;
970  	            diff = rhs - sumrhs;
971  	            residual = 0.0;
972  	
973  	            for( b = 0; b < linking->nblocks; b++ )
974  	            {
975  	               goalvalue = linking->currentrhs[b] + ( diff / linking->nblocks ) + residual;
976  	               part = MIN(goalvalue, linking->maxactivity[b]);
977  	               residual = goalvalue - part;
978  	               linking->currentrhs[b] = part;
979  	            }
980  	            if( !SCIPisZero(scip, residual) )
981  	               linking->currentrhs[0] += residual;
982  	         }
983  	         if( !SCIPisEQ(scip, lhs, sumlhs) && linking->haslhs )
984  	         {
985  	            SCIP_Real diff;
986  	            SCIP_Real part;
987  	            SCIP_Real residual;
988  	            diff = sumlhs - lhs; /* always positive*/
989  	            residual = 0.0;
990  	
991  	            for( b = 0; b < linking->nblocks; b++ )
992  	            {
993  	               goalvalue = linking->currentlhs[b] - ( diff / linking->nblocks ) + residual;
994  	               part = MAX(goalvalue, linking->minactivity[b]);
995  	               residual = goalvalue - part;
996  	               linking->currentlhs[b] = part;
997  	            }
998  	            if( !SCIPisZero(scip, residual) )
999  	               linking->currentlhs[0] += residual;
1000 	         }
1001 	      }
1002 	      /* equal parts for each block with respect to minimal/maximal activity */
1003 	      else if( linking->hasrhs || linking->haslhs )
1004 	      {
1005 	         if( linking->hasrhs )
1006 	         {
1007 	            for( b = 0; b < linking->nblocks; b++ )
1008 	            {
1009 	               goalvalue = residualrhs / (linking->nblocks - b);
1010 	               linking->currentrhs[b] = MIN(MAX(goalvalue, linking->minactivity[b]), linking->maxactivity[b]);
1011 	               residualrhs -= linking->currentrhs[b];
1012 	            }
1013 	            /* add residual partition to first block */
1014 	            linking->currentrhs[0] += residualrhs;
1015 	         }
1016 	         if( linking->haslhs )
1017 	         {
1018 	            for( b = 0; b < linking->nblocks; b++ )
1019 	            {
1020 	               goalvalue = residuallhs / (linking->nblocks - b);
1021 	               linking->currentlhs[b] = MIN(MAX(goalvalue, linking->minactivity[b]), linking->maxactivity[b]);
1022 	               residuallhs -= linking->currentlhs[b];
1023 	            }
1024 	            /* add residual partition to first block */
1025 	            linking->currentlhs[0] += residuallhs;
1026 	         }
1027 	      }
1028 	      else
1029 	      {
1030 	         assert(linking->nblocks == 0 && !SCIPconsIsChecked(linking->linkingcons));
1031 	      }
1032 	
1033 	      SCIP_CALL( roundPartition(scip, linking, blockproblem, linking->hasrhs) );
1034 	
1035 	      /* set sides in blockproblem at initial partition */
1036 	      for( b = 0; b < linking->nblocks; b++ )
1037 	      {
1038 	         if( linking->hasrhs )
1039 	         {
1040 	            SCIP_CALL( SCIPchgRhsLinear(blockproblem[linking->blocknumbers[b]]->blockscip,
1041 	                                       linking->blockconss[b], linking->currentrhs[b]) );
1042 	#ifdef SCIP_MORE_DEBUG
1043 	            SCIPdebugMsg(scip, "change rhs of %s in block %d to %f\n",
1044 	                        SCIPconsGetName(linking->linkingcons), linking->blocknumbers[b], linking->currentrhs[b]);
1045 	#endif
1046 	         }
1047 	         if( linking->haslhs )
1048 	         {
1049 	            SCIP_CALL( SCIPchgLhsLinear(blockproblem[linking->blocknumbers[b]]->blockscip,
1050 	                                       linking->blockconss[b], linking->currentlhs[b]) );
1051 	#ifdef SCIP_MORE_DEBUG
1052 	            SCIPdebugMsg(scip, "change lhs of %s in block %d to %f\n",
1053 	                        SCIPconsGetName(linking->linkingcons), linking->blocknumbers[b], linking->currentlhs[b]);
1054 	#endif
1055 	         }
1056 	      }
1057 	   }
1058 	
1059 	   return SCIP_OKAY;
1060 	}
1061 	
1062 	/** calculates shift */
1063 	static
1064 	SCIP_RETCODE calculateShift(
1065 	   SCIP*                 scip,               /**< SCIP data structure of main scip */
1066 	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
1067 	   LINKING*              linking,            /**< one linking data structure */
1068 	   SCIP_Real**           shift,              /**< pointer to store direction vector */
1069 	   int*                  nviolatedblocksrhs, /**< pointer to store number of blocks which violate rhs */
1070 	   int*                  nviolatedblockslhs, /**< pointer to store number of blocks which violate lhs */
1071 	   SCIP_Bool*            update              /**< should we update the partition? */
1072 	   )
1073 	{
(1) Event var_decl: Declaring variable "nonviolatedblocksrhs" without initializer.
Also see events: [uninit_use]
1074 	   int* nonviolatedblocksrhs;
1075 	   int* nonviolatedblockslhs;
1076 	   SCIP_Real sumviols = 0.0;
1077 	   int v;
1078 	
(2) Event cond_false: Condition "linking->hasrhs", taking false branch.
1079 	   if( linking->hasrhs )
1080 	   {
1081 	      SCIP_CALL( SCIPallocBufferArray(scip, &nonviolatedblocksrhs, linking->nblocks) );
(3) Event if_end: End of if statement.
1082 	   }
(4) Event cond_true: Condition "linking->haslhs", taking true branch.
1083 	   if( linking->haslhs )
1084 	   {
(5) Event cond_false: Condition "(nonviolatedblockslhs = BMSallocBufferMemoryArray_call(SCIPbuffer(scip), (size_t)(ptrdiff_t)linking->nblocks, 4UL /* sizeof (*nonviolatedblockslhs) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/heur_dps.c", 1085)) == NULL", taking false branch.
(6) Event cond_false: Condition "(_restat_ = (((nonviolatedblockslhs = BMSallocBufferMemoryArray_call(SCIPbuffer(scip), (size_t)(ptrdiff_t)linking->nblocks, 4UL /* sizeof (*nonviolatedblockslhs) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/heur_dps.c", 1085)) == NULL) ? SCIP_NOMEMORY : SCIP_OKAY)) != SCIP_OKAY", taking false branch.
(7) Event if_end: End of if statement.
1085 	      SCIP_CALL( SCIPallocBufferArray(scip, &nonviolatedblockslhs, linking->nblocks) );
1086 	   }
1087 	
1088 	   /* get violated blocks and calculate shift */
(8) Event cond_true: Condition "v < linking->nblocks", taking true branch.
(16) Event loop_begin: Jumped back to beginning of loop.
(17) Event cond_false: Condition "v < linking->nblocks", taking false branch.
1089 	   for( v = 0; v < linking->nblocks; v++ )
1090 	   {
1091 	      SCIP* subscip;
1092 	      SCIP_SOL* subsol;
1093 	      SCIP_Real slackval;
1094 	
1095 	      subscip = blockproblem[linking->blocknumbers[v]]->blockscip;
1096 	      subsol = SCIPgetBestSol(subscip);
1097 	
1098 	      /* if we have ranged constraints, the slack variables of the rhs are in front;
1099 	       * get slack variable of block; save violated blocks
1100 	       */
(9) Event cond_false: Condition "linking->hasrhs", taking false branch.
1101 	      if( linking->hasrhs )
1102 	      {
1103 	         slackval = SCIPgetSolVal(subscip, subsol, linking->slacks[v * linking->nslacksperblock]);
1104 	
1105 	         /* block is violated */
1106 	         if( SCIPisPositive(scip, slackval) )
1107 	         {
1108 	            (*nviolatedblocksrhs)++;
1109 	
1110 	            (*shift)[v] += slackval;
1111 	            sumviols += slackval;
1112 	         }
1113 	         else
1114 	         {
1115 	            nonviolatedblocksrhs[v - *nviolatedblocksrhs] = v; /*lint !e644*/
1116 	         }
(10) Event if_end: End of if statement.
1117 	      }
(11) Event cond_true: Condition "linking->haslhs", taking true branch.
1118 	      if( linking->haslhs )
1119 	      {
1120 	         slackval = SCIPgetSolVal(subscip, subsol, linking->slacks[(v * linking->nslacksperblock) + linking->nslacksperblock - 1]);
1121 	
1122 	         /* block is violated */
(12) Event cond_true: Condition "slackval > scip->set->num_epsilon", taking true branch.
1123 	         if( SCIPisPositive(scip, slackval) )
1124 	         {
1125 	            (*nviolatedblockslhs)++;
1126 	
1127 	            (*shift)[v] -= slackval;
1128 	            sumviols -= slackval;
(13) Event if_fallthrough: Falling through to end of if statement.
1129 	         }
1130 	         else
1131 	         {
1132 	            nonviolatedblockslhs[v - *nviolatedblockslhs] = v; /*lint !e644*/
(14) Event if_end: End of if statement.
1133 	         }
1134 	      }
(15) Event loop: Jumping back to the beginning of the loop.
(18) Event loop_end: Reached end of loop.
1135 	   }
1136 	
1137 	   /* linking constraint is in no block violated or always violated -> do not update partition */
(19) Event cond_false: Condition "*nviolatedblocksrhs + *nviolatedblockslhs == 0", taking false branch.
(20) Event cond_false: Condition "linking->nblocks == *nviolatedblocksrhs", taking false branch.
(21) Event cond_false: Condition "linking->nblocks == *nviolatedblockslhs", taking false branch.
1138 	   if( *nviolatedblocksrhs + *nviolatedblockslhs == 0 ||
1139 	         linking->nblocks == *nviolatedblocksrhs || linking->nblocks == *nviolatedblockslhs )
1140 	   {
1141 	      *update = FALSE;
1142 	      /* free memory */
1143 	      if( linking->haslhs )
1144 	         SCIPfreeBufferArray(scip, &nonviolatedblockslhs);
1145 	      if( linking->hasrhs )
1146 	         SCIPfreeBufferArray(scip, &nonviolatedblocksrhs);
1147 	      return SCIP_OKAY;
(22) Event if_end: End of if statement.
1148 	   }
1149 	
1150 	   /* set values of non violated blocks */
(23) Event cond_true: Condition "sumviols > scip->set->num_epsilon", taking true branch.
1151 	   if( SCIPisPositive(scip, sumviols) )
1152 	   {
1153 	      /* rhs of original linking constraint is violated */
1154 	      SCIP_Real residual = sumviols;
1155 	      SCIP_Real part;
1156 	      SCIP_Real shift_tmp;
1157 	
1158 	      assert(linking->hasrhs);
1159 	      assert(*nviolatedblocksrhs != 0);
1160 	
1161 	      /* substract from each non violated block the same amount with respect to minimal/maximal activity,
1162 	       * so that the shift is zero in sum
1163 	       */
(24) Event cond_true: Condition "v < linking->nblocks - *nviolatedblocksrhs", taking true branch.
1164 	      for( v = 0; v < linking->nblocks - *nviolatedblocksrhs; v++ )
1165 	      {
(25) Event uninit_use: Using uninitialized value "nonviolatedblocksrhs".
Also see events: [var_decl]
1166 	         part = linking->currentrhs[nonviolatedblocksrhs[v]] - residual/(linking->nblocks - *nviolatedblocksrhs - v);
1167 	         part = MIN(MAX(part, linking->minactivity[nonviolatedblocksrhs[v]]), linking->maxactivity[nonviolatedblocksrhs[v]]);
1168 	         shift_tmp = part - linking->currentrhs[nonviolatedblocksrhs[v]];
1169 	         residual += shift_tmp;
1170 	         (*shift)[nonviolatedblocksrhs[v]] += shift_tmp;
1171 	      }
1172 	      if( !SCIPisZero(scip, residual) )
1173 	      {
1174 	         /* assign residual to ... */
1175 	         if( linking->nblocks - *nviolatedblocksrhs == 1 )
1176 	            (*shift)[nonviolatedblocksrhs[0] == 0 ? 1 : 0] -= residual; /* first violated block */
1177 	         else
1178 	            (*shift)[nonviolatedblocksrhs[0]] -= residual; /* first nonviolated block */
1179 	      }
1180 	   }
1181 	   if( SCIPisNegative(scip, sumviols) )
1182 	   {
1183 	      /* lhs of original linking constraint is violated */
1184 	      SCIP_Real residual = sumviols;
1185 	      SCIP_Real part;
1186 	      SCIP_Real shift_tmp;
1187 	
1188 	      assert(linking->haslhs);
1189 	      assert(*nviolatedblockslhs != 0);
1190 	
1191 	      /* add to each non violated block the same amount with respect to minimal/maximal activity,
1192 	       * so that the shift is zero in sum
1193 	       */
1194 	      for( v = 0; v < linking->nblocks - *nviolatedblockslhs; v++ )
1195 	      {
1196 	         part = linking->currentlhs[nonviolatedblockslhs[v]] - residual/(linking->nblocks - *nviolatedblockslhs - v);
1197 	         part = MIN(MAX(part, linking->minactivity[nonviolatedblockslhs[v]]), linking->maxactivity[nonviolatedblockslhs[v]]);
1198 	         shift_tmp = part - linking->currentlhs[nonviolatedblockslhs[v]];
1199 	         residual += shift_tmp;
1200 	         (*shift)[nonviolatedblockslhs[v]] += shift_tmp;
1201 	      }
1202 	      if( !SCIPisZero(scip, residual) )
1203 	      {
1204 	         /* assign residual to ... */
1205 	         if( linking->nblocks - *nviolatedblockslhs == 1 )
1206 	            (*shift)[nonviolatedblockslhs[0] == 0 ? 1 : 0] -= residual; /* first violated block */
1207 	         else
1208 	            (*shift)[nonviolatedblockslhs[0]] -= residual; /* first nonviolated block */
1209 	      }
1210 	   }
1211 	
1212 	   *update = TRUE;
1213 	
1214 	   /* free memory */
1215 	   if( linking->haslhs )
1216 	      SCIPfreeBufferArray(scip, &nonviolatedblockslhs);
1217 	   if( linking->hasrhs )
1218 	      SCIPfreeBufferArray(scip, &nonviolatedblocksrhs);
1219 	
1220 	   return SCIP_OKAY;
1221 	}
1222 	
1223 	/** update partition */
1224 	static
1225 	SCIP_RETCODE updatePartition(
1226 	   SCIP*                 scip,               /**< SCIP data structure of main scip */
1227 	   LINKING**             linkings,           /**< linking data structure */
1228 	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
1229 	   int**                 nviolatedblocksrhs, /**< pointer to store number of blocks which violate rhs */
1230 	   int**                 nviolatedblockslhs, /**< pointer to store number of blocks which violate lhs */
1231 	   int                   nlinking,           /**< number of linking constraints */
1232 	   int                   nblocks,            /**< number of blocks */
1233 	   int                   iteration,          /**< number of iteration */
1234 	   SCIP_Bool*            oneupdate           /**< is at least partition of one constraint updated? */
1235 	   )
1236 	{
1237 	   SCIP_Real* shift; /* direction vector */
1238 	   int v;
1239 	   int c;
1240 	   SCIP_Bool update; /* is partition of current constraint updated? */
1241 	
1242 	   assert(scip != NULL);
1243 	   assert(linkings != NULL);
1244 	   assert(blockproblem != NULL);
1245 	   assert(iteration >= 0);
1246 	   assert(!*oneupdate);
1247 	
1248 	   SCIP_CALL( SCIPallocBufferArray(scip, &shift, nblocks) ); /* allocate memory only once */
1249 	
1250 	   for( c = 0; c < nlinking; c++ )
1251 	   {
1252 	      LINKING* linking; /* one linking data structure */
1253 	
1254 	      linking = linkings[c];
1255 	      (*nviolatedblocksrhs)[c] = 0;
1256 	      (*nviolatedblockslhs)[c] = 0;
1257 	      update = TRUE;
1258 	      BMSclearMemoryArray(shift, linking->nblocks);
1259 	      assert(nblocks >= linking->nblocks);
1260 	
1261 	      SCIP_CALL( calculateShift(scip, blockproblem, linking, &shift, &(*nviolatedblocksrhs)[c], &(*nviolatedblockslhs)[c], &update) );
1262 	
1263 	      if( !update )
1264 	         continue;
1265 	
1266 	      *oneupdate = TRUE;
1267 	
1268 	#ifdef SCIP_DEBUG
1269 	      SCIP_Real sum = 0.0;
1270 	      /* check sum of shift; must be zero */
1271 	      for( v = 0; v < linking->nblocks; v++ )
1272 	         sum += shift[v];
1273 	      assert(SCIPisFeasZero(scip, sum));
1274 	#endif
1275 	
1276 	      /* add shift to both sides */
1277 	      for( v = 0; v < linking->nblocks; v++)
1278 	      {
1279 	         if( linking->hasrhs )
1280 	            linking->currentrhs[v] += shift[v];
1281 	
1282 	         if( linking->haslhs )
1283 	            linking->currentlhs[v] += shift[v];
1284 	      }
1285 	
1286 	      SCIP_CALL( roundPartition(scip, linking, blockproblem, ((linking->hasrhs && ((*nviolatedblocksrhs)[c] != 0)) || !linking->haslhs)) );
1287 	
1288 	      /* set sides in blockproblems to new partition */
1289 	      for( v = 0; v < linking->nblocks; v++ )
1290 	      {
1291 	         SCIP* subscip;
1292 	         subscip = blockproblem[linking->blocknumbers[v]]->blockscip;
1293 	
1294 	         if( linking->hasrhs )
1295 	         {
1296 	            SCIP_CALL( SCIPchgRhsLinear(subscip, linking->blockconss[v], linking->currentrhs[v]) );
1297 	         }
1298 	         if( linking->haslhs )
1299 	         {
1300 	            SCIP_CALL( SCIPchgLhsLinear(subscip, linking->blockconss[v], linking->currentlhs[v]) );
1301 	         }
1302 	      }
1303 	   }
1304 	
1305 	   /* free memory */
1306 	   SCIPfreeBufferArray(scip, &shift);
1307 	
1308 	   return SCIP_OKAY;
1309 	}
1310 	
1311 	/** update penalty parameters lambda
1312 	 *
1313 	 * if a linking constraint is violated two times in succession, the corresponding penalty parameter is increased in each block
1314 	 */
1315 	static
1316 	SCIP_RETCODE updateLambda(
1317 	   SCIP*                 scip,               /**< SCIP data structure */
1318 	   SCIP_HEURDATA*        heurdata,           /**< heuristic data */
1319 	   LINKING**             linkings,           /**< array of linking data structures */
1320 	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
1321 	   int*                  nviolatedblocksrhs, /**< number of blocks which violate rhs */
1322 	   int*                  nviolatedblockslhs, /**< number of blocks which violate lhs */
1323 	   int                   nlinking            /**< number of linking constraints */
1324 	   )
1325 	{
1326 	   SCIP_VAR* slackvar;
1327 	   SCIP_Real old_obj;
1328 	   SCIP_Real new_obj;
1329 	   int c;
1330 	   int b;
1331 	
1332 	   assert(scip != NULL);
1333 	   assert(linkings != NULL);
1334 	   assert(blockproblem != NULL);
1335 	
1336 	   for( c = 0; c < nlinking; c++ )
1337 	   {
1338 	      assert(nviolatedblocksrhs[c] >= 0);
1339 	      assert(nviolatedblockslhs[c] >= 0);
1340 	
1341 	      /* skip constraint if it is not violated */
1342 	      if( nviolatedblocksrhs[c] + nviolatedblockslhs[c] == 0 )
1343 	      {
1344 	         linkings[c]->lastviolations = 0; /* reset flag */
1345 	         continue;
1346 	      }
1347 	
1348 	      /* add number of violated blocks multiplied with parameter "penalty" to lambda (initial value is 1) */
1349 	      for( b = 0; b < linkings[c]->nblocks; b++ )
1350 	      {
1351 	         SCIP* subscip;
1352 	         subscip = blockproblem[linkings[c]->blocknumbers[b]]->blockscip;
1353 	         assert(subscip != NULL);
1354 	
1355 	         if( linkings[c]->hasrhs && (nviolatedblocksrhs[c] >= 1) && (linkings[c]->lastviolations >= 1) )
1356 	         {
1357 	            slackvar = linkings[c]->slacks[b * linkings[c]->nslacksperblock];
1358 	            old_obj = SCIPvarGetObj(slackvar);
1359 	            new_obj = old_obj + heurdata->penalty * nviolatedblocksrhs[c];
1360 	
1361 	            SCIP_CALL( SCIPchgVarObj(subscip, slackvar, new_obj) );
1362 	         }
1363 	         if( linkings[c]->haslhs && (nviolatedblockslhs[c] >= 1) && (linkings[c]->lastviolations >= 1) )
1364 	         {
1365 	            slackvar = linkings[c]->slacks[b * linkings[c]->nslacksperblock + linkings[c]->nslacksperblock - 1];
1366 	            old_obj = SCIPvarGetObj(slackvar);
1367 	            new_obj = old_obj + heurdata->penalty * nviolatedblockslhs[c];
1368 	
1369 	            SCIP_CALL( SCIPchgVarObj(subscip, slackvar, new_obj) );
1370 	         }
1371 	      }
1372 	
1373 	      /* update number of violations in the last iterations */
1374 	      linkings[c]->lastviolations += 1;
1375 	   }
1376 	
1377 	   return SCIP_OKAY;
1378 	}
1379 	
1380 	/** computes feasible solution from last stored solution for each block and adds it to the solution storage */
1381 	static
1382 	SCIP_RETCODE reuseSolution(
1383 	   LINKING**             linkings,           /**< array of linking data structures */
1384 	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
1385 	   int                   nblocks             /**< number of blocks */
1386 	   )
1387 	{
1388 	   SCIP_SOL** sols;
1389 	   SCIP_SOL* sol; /* solution of block that will be repaired */
1390 	   SCIP_SOL* newsol;
1391 	   SCIP_VAR** blockvars;
1392 	   SCIP_Real* blockvals;
1393 	   int nsols;
1394 	   int nvars;
1395 	   int b;
1396 	   int c;
1397 	   int i;
1398 	   SCIP_Bool success;
1399 	
1400 	   assert(linkings != NULL);
1401 	   assert(blockproblem != NULL);
1402 	
1403 	   for( b = 0; b < nblocks; b++ )
1404 	   {
1405 	      SCIP* subscip;
1406 	
1407 	      subscip = blockproblem[b]->blockscip;
1408 	      nsols = SCIPgetNSols(subscip);
1409 	
1410 	      /* no solution in solution candidate storage found */
1411 	      if( nsols == 0 )
1412 	         return SCIP_OKAY;
1413 	
1414 	      /* take last solution */
1415 	      sols = SCIPgetSols(subscip);
1416 	      sol = sols[nsols - 1];
1417 	
1418 	      /* copy the solution */
1419 	      nvars = SCIPgetNVars(subscip);
1420 	      blockvars = SCIPgetVars(subscip);
1421 	      SCIP_CALL( SCIPallocBufferArray(subscip, &blockvals, nvars) );
1422 	      SCIP_CALL( SCIPgetSolVals(subscip, sol, nvars, blockvars, blockvals) );
1423 	      SCIP_CALL( SCIPcreateOrigSol(subscip, &newsol, NULL) );
1424 	      SCIP_CALL( SCIPsetSolVals(subscip, newsol, nvars, blockvars, blockvals) );
1425 	
1426 	      /* correct each coupling constraint:
1427 	      * lhs <= orig_var - z_r + z_l <= rhs
1428 	      * adapt slack variables so that constraint is feasible
1429 	      */
1430 	      for( c = 0; c < blockproblem[b]->nlinking; c++ )
1431 	      {
1432 	         LINKING* linking; /* linking data structure of this constraint */
1433 	         SCIP_VAR* rvar; /* slack variable z_r */
1434 	         SCIP_VAR* lvar; /* slack variable z_l */
1435 	         SCIP_Real rval; /* value of slack variable z_r */
1436 	         SCIP_Real lval; /* value of slack variable z_l */
1437 	         SCIP_Real activitycons; /* activity of constraint*/
1438 	         SCIP_Real activity; /* activity of constraint without slack variables */
1439 	         SCIP_Real rhs; /* current right hand side */
1440 	         SCIP_Real lhs; /* current left hand side */
1441 	
1442 	         linking = linkings[blockproblem[b]->linkingindices[c]];
1443 	         rhs = SCIPgetRhsLinear(subscip, blockproblem[b]->linkingconss[c]);
1444 	         lhs = SCIPgetLhsLinear(subscip, blockproblem[b]->linkingconss[c]);
1445 	         assert(SCIPisGE(subscip, rhs, lhs));
1446 	
1447 	         activitycons = SCIPgetActivityLinear(subscip, blockproblem[b]->linkingconss[c], sol);
1448 	
1449 	         /* get slack variables and subtract their value from the activity;
1450 	          * calculate and set values of slack variables
1451 	          */
1452 	         for( i = 0; i < linking->nblocks; i++ )
1453 	         {
1454 	            if( linking->blocknumbers[i] == b )
1455 	            {
1456 	               if( linking->hasrhs && linking->haslhs )
1457 	               {
1458 	                  rvar = linking->slacks[2 * i];
1459 	                  lvar = linking->slacks[2 * i + 1];
1460 	                  rval = SCIPgetSolVal(subscip, sol, rvar);
1461 	                  lval = SCIPgetSolVal(subscip, sol, lvar);
1462 	                  activity = activitycons + rval - lval;
1463 	                  SCIP_CALL( SCIPsetSolVal(subscip, newsol, rvar, MAX(0.0, activity - rhs)) );
1464 	                  SCIP_CALL( SCIPsetSolVal(subscip, newsol, lvar, MAX(0.0, lhs - activity)) );
1465 	               }
1466 	               else if( linking->hasrhs )
1467 	               {
1468 	                  rvar = linking->slacks[i];
1469 	                  rval = SCIPgetSolVal(subscip, sol, rvar);
1470 	                  activity = activitycons + rval;
1471 	                  SCIP_CALL( SCIPsetSolVal(subscip, newsol, rvar, MAX(0.0, activity - rhs)) );
1472 	               }
1473 	               else /* linking->haslhs */
1474 	               {
1475 	                  assert(linking->haslhs);
1476 	                  lvar = linking->slacks[i];
1477 	                  lval = SCIPgetSolVal(subscip, sol, lvar);
1478 	                  activity = activitycons - lval;
1479 	                  SCIP_CALL( SCIPsetSolVal(subscip, newsol, lvar, MAX(0.0, lhs - activity)) );
1480 	               }
1481 	               break;
1482 	            }
1483 	         }
1484 	      }
1485 	
1486 	      SCIPdebugMsg(subscip, "Try adding solution with objective value %.2f\n", SCIPgetSolOrigObj(subscip, newsol));
1487 	      SCIP_CALL( SCIPaddSolFree(subscip, &newsol, &success) );
1488 	
1489 	      if( !success )
1490 	         SCIPdebugMsg(subscip, "Correcting solution failed\n"); /* maybe not better than old solutions */
1491 	      else
1492 	         SCIPdebugMsg(subscip, "Correcting solution successful\n");
1493 	
1494 	      SCIPfreeBufferArray(subscip, &blockvals);
1495 	   }
1496 	
1497 	   return SCIP_OKAY;
1498 	}
1499 	
1500 	/** reoptimizes the heuristic solution with original objective function */
1501 	static
1502 	SCIP_RETCODE reoptimize(
1503 	   SCIP*                 scip,               /**< SCIP data structure */
1504 	   SCIP_HEUR*            heur,               /**< pointer to heuristic */
1505 	   SCIP_SOL*             sol,                /**< heuristic solution */
1506 	   BLOCKPROBLEM**        blockproblem,       /**< array of blockproblem data structures */
1507 	   int                   nblocks,            /**< number of blockproblems */
1508 	   SCIP_Bool             limits,             /**< should strict limits be set? */
1509 	   SCIP_SOL**            newsol,             /**< pointer to store improved solution */
1510 	   SCIP_Bool*            success             /**< pointer to store whether reoptimization was successful */
1511 	   )
1512 	{
1513 	   SCIP_Real time;
1514 	   SCIP_Real timesubscip;
1515 	   SCIP_Bool check;
1516 	   int b;
1517 	   int v;
1518 	
1519 	   assert(scip != NULL);
1520 	   assert(heur != NULL);
1521 	   assert(sol != NULL);
1522 	   assert(blockproblem != NULL);
1523 	
1524 	   *success = FALSE;
1525 	   check = FALSE;
1526 	
1527 	   /* get used time */
1528 	   timesubscip = SCIPgetTotalTime(blockproblem[0]->blockscip);
1529 	
1530 	   /* for each blockproblem:
1531 	    * - change back to original objective function
1532 	    * - fix slack variables to zero
1533 	    * - set limits and solve problem
1534 	    */
1535 	   for( b = 0; b < nblocks; b++ )
1536 	   {
1537 	      SCIP* subscip;
1538 	      SCIP_VAR** blockvars;
1539 	      SCIP_Real infvalue;
1540 	      int nvars;
1541 	      int nsols;
1542 	
1543 	      subscip = blockproblem[b]->blockscip;
1544 	      blockvars = SCIPgetOrigVars(subscip);
1545 	      nvars = SCIPgetNOrigVars(subscip);
1546 	      nsols = 0;
1547 	
1548 	      /* in order to change objective function */
1549 	      SCIP_CALL( SCIPfreeTransform(subscip) );
1550 	
1551 	      /* change back to original objective function */
1552 	      for( v = 0; v < blockproblem[b]->nblockvars; v++ )
1553 	      {
1554 	         SCIP_CALL( SCIPchgVarObj(subscip, blockvars[v], blockproblem[b]->origobj[v]) );
1555 	      }
1556 	
1557 	      /* fix slack variables to zero */
1558 	      for( v = blockproblem[b]->nblockvars; v < nvars; v++ )
1559 	      {
1560 	         SCIP_CALL( SCIPchgVarUb(subscip, blockvars[v], 0.0) );
1561 	         SCIP_CALL( SCIPchgVarLb(subscip, blockvars[v], 0.0) );
1562 	      }
1563 	
1564 	      /* do not abort subproblem on CTRL-C */
1565 	      SCIP_CALL( SCIPsetBoolParam(subscip, "misc/catchctrlc", FALSE) );
1566 	
1567 	      /* forbid recursive call of heuristics and separators solving sub-SCIPs */
1568 	      SCIP_CALL( SCIPsetSubscipsOff(subscip, TRUE) );
1569 	
1570 	#ifdef SCIP_DEBUG
1571 	      /* for debugging, enable full output */
1572 	      SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) );
1573 	      SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 100000000) );
1574 	#else
1575 	      /* disable statistic timing inside sub SCIP and output to console */
1576 	      SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) );
1577 	      SCIP_CALL( SCIPsetBoolParam(subscip, "timing/statistictiming", FALSE) );
1578 	#endif
1579 	
1580 	      /* disable cutting plane separation */
1581 	      SCIP_CALL( SCIPsetSeparating(subscip, SCIP_PARAMSETTING_OFF, TRUE) );
1582 	
1583 	      /* disable expensive presolving */
1584 	      SCIP_CALL( SCIPsetPresolving(subscip, SCIP_PARAMSETTING_FAST, TRUE) );
1585 	
1586 	      /* disable expensive techniques */
1587 	      SCIP_CALL( SCIPsetIntParam(subscip, "misc/usesymmetry", 0) );
1588 	
1589 	      /* speed up sub-SCIP by not checking dual LP feasibility */
1590 	      SCIP_CALL( SCIPsetBoolParam(subscip, "lp/checkdualfeas", FALSE) );
1591 	
1592 	      /* copy value for infinity */
1593 	      SCIP_CALL( SCIPgetRealParam(scip, "numerics/infinity", &infvalue) );
1594 	      SCIP_CALL( SCIPsetRealParam(subscip, "numerics/infinity", infvalue) );
1595 	
1596 	      /* set limits; do not use more time in each subproblem than the heuristic has already used for first solution */
1597 	      SCIP_CALL( SCIPcopyLimits(scip, subscip) );
1598 	      if( limits )
1599 	      {
1600 	         SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", 1LL) );
1601 	         SCIP_CALL( SCIPgetRealParam(subscip, "limits/time", &time) );
1602 	         if( timesubscip <  time - 1.0 )
1603 	         {
1604 	            SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", timesubscip + 1.0) );
1605 	         }
1606 	         SCIP_CALL( SCIPtransformProb(subscip) );
1607 	         nsols = SCIPgetNSols(subscip);
1608 	         /* one improving solution */
1609 	         SCIP_CALL( SCIPsetIntParam(subscip, "limits/bestsol", nsols + 1) );
1610 	      }
1611 	      else
1612 	      {
1613 	         /* only 50% of remaining time */
1614 	         SCIP_CALL( SCIPgetRealParam(subscip, "limits/time", &time) );
1615 	         SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", time * 0.5) );
1616 	         /* set big node limits */
1617 	         SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", 1000LL) );
1618 	         SCIP_CALL( SCIPsetLongintParam(subscip, "limits/stallnodes", 100LL) );
1619 	      }
1620 	
1621 	      /* reoptimize problem */
1622 	      SCIP_CALL_ABORT( SCIPsolve(subscip) );
1623 	
1624 	      if( SCIPgetNSols(subscip) == 0 )
1625 	      {
1626 	         /* we found no solution */
1627 	         return SCIP_OKAY;
1628 	      }
1629 	      else if( SCIPgetStatus(subscip) == SCIP_STATUS_BESTSOLLIMIT ||
1630 	               SCIPgetStatus(subscip) == SCIP_STATUS_OPTIMAL ||
1631 	               SCIPgetNSols(subscip) > nsols )
1632 	      {
1633 	         check = TRUE;
1634 	      }
1635 	   }
1636 	
1637 	   if( !check )
1638 	   {
1639 	      /* we have no better solution */
1640 	      return SCIP_OKAY;
1641 	   }
1642 	
1643 	   /* create sol of main scip */
1644 	   SCIP_CALL( SCIPcreateSol(scip, newsol, heur) );
1645 	
1646 	   /* copy solution to main scip */
1647 	   for( b = 0; b < nblocks; b++ )
1648 	   {
1649 	      SCIP_SOL* blocksol;
1650 	      SCIP_VAR** blockvars;
1651 	      SCIP_Real* blocksolvals;
1652 	      int nblockvars;
1653 	
1654 	      /* get solution of block variables (without slack variables) */
1655 	      blocksol = SCIPgetBestSol(blockproblem[b]->blockscip);
1656 	      blockvars = SCIPgetOrigVars(blockproblem[b]->blockscip);
1657 	      nblockvars = blockproblem[b]->nblockvars;
1658 	
1659 	      SCIP_CALL( SCIPallocBufferArray(scip, &blocksolvals, nblockvars) );
1660 	      SCIP_CALL( SCIPgetSolVals(blockproblem[b]->blockscip, blocksol, nblockvars, blockvars, blocksolvals) );
1661 	
1662 	      for( v = 0; v < nblockvars; v++ )
1663 	      {
1664 	         SCIP_VAR* origvar;
1665 	
1666 	         origvar = SCIPfindVar(scip, SCIPvarGetName(blockvars[v]));
1667 	         SCIP_CALL( SCIPsetSolVal(scip, *newsol, origvar, blocksolvals[v]) );
1668 	      }
1669 	
1670 	      SCIPfreeBufferArray(scip, &blocksolvals);
1671 	   }
1672 	
1673 	   *success = TRUE;
1674 	
1675 	   return SCIP_OKAY;
1676 	}
1677 	
1678 	
1679 	/* ---------------- Callback methods of event handler ---------------- */
1680 	
1681 	/* exec the event handler
1682 	 *
1683 	 * interrupt solution process of sub-SCIP if dual bound is greater than zero and a solution is available
1684 	 */
1685 	static
1686 	SCIP_DECL_EVENTEXEC(eventExecDps)
1687 	{
1688 	   assert(eventhdlr != NULL);
1689 	   assert(eventdata != NULL);
1690 	   assert(strcmp(SCIPeventhdlrGetName(eventhdlr), EVENTHDLR_NAME) == 0);
1691 	   assert(event != NULL);
1692 	   assert(SCIPeventGetType(event) & SCIP_EVENTTYPE_LPSOLVED);
1693 	
1694 	   SCIPdebugMsg(scip, "dual bound: %0.2f\n", SCIPgetDualbound(scip));
1695 	
1696 	   if( SCIPisFeasGT(scip, SCIPgetDualbound(scip), 0.0) && SCIPgetNSols(scip) >= 1 )
1697 	   {
1698 	      SCIPdebugMsg(scip, "DPS: interrupt subscip\n");
1699 	      SCIP_CALL( SCIPinterruptSolve(scip) );
1700 	   }
1701 	
1702 	   return SCIP_OKAY;
1703 	}
1704 	
1705 	
1706 	/*
1707 	 * Callback methods of primal heuristic
1708 	 */
1709 	
1710 	/** copy method for primal heuristic plugins (called when SCIP copies plugins) */
1711 	static
1712 	SCIP_DECL_HEURCOPY(heurCopyDps)
1713 	{  /*lint --e{715}*/
1714 	   assert(scip != NULL);
1715 	   assert(heur != NULL);
1716 	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1717 	
1718 	   /* call inclusion method of primal heuristic */
1719 	   SCIP_CALL( SCIPincludeHeurDps(scip) );
1720 	
1721 	   return SCIP_OKAY;
1722 	}
1723 	
1724 	/** destructor of primal heuristic to free user data (called when SCIP is exiting) */
1725 	static
1726 	SCIP_DECL_HEURFREE(heurFreeDps)
1727 	{  /*lint --e{715}*/
1728 	   SCIP_HEURDATA* heurdata;
1729 	
1730 	   assert(heur != NULL);
1731 	   assert(scip != NULL);
1732 	
1733 	   /* free heuristic data */
1734 	   heurdata = SCIPheurGetData(heur);
1735 	   assert(heurdata != NULL);
1736 	
1737 	   SCIPfreeBlockMemory(scip, &heurdata);
1738 	   SCIPheurSetData(heur, NULL);
1739 	
1740 	   return SCIP_OKAY;
1741 	}
1742 	
1743 	/** execution method of primal heuristic */
1744 	static
1745 	SCIP_DECL_HEUREXEC(heurExecDps)
1746 	{  /*lint --e{715}*/
1747 	   SCIP_DECOMP** alldecomps;
1748 	   SCIP_DECOMP* decomp;
1749 	   SCIP_DECOMP* assigneddecomp;
1750 	   SCIP_VAR** vars;
1751 	   SCIP_CONS** conss;
1752 	   SCIP_VAR** sortedvars;
1753 	   SCIP_CONS** sortedconss;
1754 	   SCIP_HEURDATA* heurdata;
1755 	   BLOCKPROBLEM** blockproblem;
1756 	   LINKING** linkings;
1757 	   int* sortedvarlabels;
1758 	   int* sortedconslabels;
1759 	   SCIP_EVENTHDLR* eventhdlr; /* event handler */
1760 	   SCIP_Real memory; /* in MB */
1761 	   SCIP_Real timelimit;
1762 	   SCIP_Real allslacksval;
1763 	   SCIP_Real blocksolval;
1764 	   SCIP_STATUS status;
1765 	   SCIP_Bool avoidmemout;
1766 	   SCIP_Bool disablemeasures;
1767 	   int maxgraphedge;
1768 	   int ndecomps;
1769 	   int nvars;
1770 	   int nconss;
1771 	   int nblocks;
1772 	   SCIP_Bool success;
1773 	   int b;
1774 	   int c;
1775 	   int k;
1776 	
1777 	   assert( heur != NULL );
1778 	   assert( scip != NULL );
1779 	   assert( result != NULL );
1780 	
1781 	   heurdata = SCIPheurGetData(heur);
1782 	   assert(heurdata != NULL);
1783 	
1784 	   assigneddecomp = NULL;
1785 	   blockproblem = NULL;
1786 	   linkings = NULL;
1787 	   eventhdlr = NULL;
1788 	
1789 	   *result = SCIP_DIDNOTRUN;
1790 	
1791 	   if( !((heurtiming & SCIP_HEURTIMING_BEFORENODE) && (heurdata->timing == 0 || heurdata->timing == 2))
1792 	         && !((heurtiming & SCIP_HEURTIMING_AFTERNODE) && (heurdata->timing == 1 || heurdata->timing == 2)) )
1793 	   {
1794 	      return SCIP_OKAY;
1795 	   }
1796 	
1797 	   /* call heuristic only once */
1798 	   if( SCIPheurGetNCalls(heur) > 0 )
1799 	      return SCIP_OKAY;
1800 	
1801 	   /* -------------------------------------------------------------------- */
1802 	   SCIPdebugMsg(scip, "initialize dps heuristic\n");
1803 	
1804 	   /* take the first transformed decomposition */
1805 	   SCIPgetDecomps(scip, &alldecomps, &ndecomps, FALSE);
1806 	   if( ndecomps == 0 )
1807 	      return SCIP_OKAY;
1808 	
1809 	   decomp = alldecomps[0];
1810 	   assert(decomp != NULL);
1811 	   SCIPdebugMsg(scip, "First transformed decomposition is selected\n");
1812 	
1813 	   nblocks = SCIPdecompGetNBlocks(decomp);
1814 	   nconss = SCIPgetNConss(scip);
1815 	   nvars = SCIPgetNVars(scip);
1816 	
1817 	   /* if problem has no constraints, no variables or less than two blocks, return */
1818 	   if( nconss == 0 || nvars == 0 || nblocks <= 1 )
1819 	   {
1820 	      SCIPdebugMsg(scip, "problem has no constraints, no variables or less than two blocks\n");
1821 	      return SCIP_OKAY;
1822 	   }
1823 	
1824 	   /* estimate required memory for all blocks and terminate if not enough memory is available */
1825 	   SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &memory) );
1826 	   SCIP_CALL( SCIPgetBoolParam(scip, "misc/avoidmemout", &avoidmemout) );
1827 	   if( avoidmemout && (((SCIPgetMemUsed(scip) + SCIPgetMemExternEstim(scip))/1048576.0) * (nblocks/4.0 + 2) >= memory) )
1828 	   {
1829 	      SCIPdebugMsg(scip, "The estimated memory usage for %d blocks is too large.\n", nblocks);
1830 	      return SCIP_OKAY;
1831 	   }
1832 	
1833 	   /* we do not need the block decomposition graph and expensive measures of the decomposition statistics */
1834 	   SCIP_CALL( SCIPgetIntParam(scip, "decomposition/maxgraphedge", &maxgraphedge) );
1835 	   if( !SCIPisParamFixed(scip, "decomposition/maxgraphedge") )
1836 	   {
1837 	      SCIP_CALL( SCIPsetIntParam(scip, "decomposition/maxgraphedge", 0) );
1838 	   }
1839 	   SCIP_CALL( SCIPgetBoolParam(scip, "decomposition/disablemeasures", &disablemeasures) );
1840 	   if( !SCIPisParamFixed(scip, "decomposition/disablemeasures") )
1841 	   {
1842 	      SCIP_CALL( SCIPsetBoolParam(scip, "decomposition/disablemeasures", TRUE) );
1843 	   }
1844 	
1845 	   vars = SCIPgetVars(scip);
1846 	   conss = SCIPgetConss(scip);
1847 	   SCIP_CALL( SCIPduplicateBufferArray(scip, &sortedvars, vars, nvars) );
1848 	   SCIP_CALL( SCIPduplicateBufferArray(scip, &sortedconss, conss, nconss) );
1849 	   SCIP_CALL( SCIPallocBufferArray(scip, &sortedvarlabels, nvars) );
1850 	   SCIP_CALL( SCIPallocBufferArray(scip, &sortedconslabels, nconss) );
1851 	
1852 	   /* get labels and sort in increasing order */
1853 	   SCIPdecompGetVarsLabels(decomp, sortedvars, sortedvarlabels, nvars);
1854 	   SCIPdecompGetConsLabels(decomp, sortedconss, sortedconslabels, nconss);
1855 	   SCIPsortIntPtr(sortedvarlabels, (void**)sortedvars, nvars);
1856 	   SCIPsortIntPtr(sortedconslabels, (void**)sortedconss, nconss);
1857 	
1858 	   if( sortedvarlabels[0] == SCIP_DECOMP_LINKVAR )
1859 	   {
1860 	      /* create new decomposition; don't change the decompositions in the decompstore */
1861 	      SCIP_CALL( SCIPdecompCreate(&assigneddecomp, SCIPblkmem(scip), nblocks, SCIPdecompIsOriginal(decomp), SCIPdecompUseBendersLabels(decomp)) );
1862 	
1863 	      SCIP_CALL( assignLinking(scip, assigneddecomp, sortedvars, sortedconss, sortedvarlabels, sortedconslabels, nvars, nconss, SCIPdecompGetNBorderVars(decomp)) );
1864 	      assert(SCIPdecompGetNBlocks(decomp) >= SCIPdecompGetNBlocks(assigneddecomp));
1865 	      decomp = assigneddecomp;
1866 	
1867 	      /* number of blocks can get smaller */
1868 	      nblocks = SCIPdecompGetNBlocks(decomp);
1869 	   }
1870 	   else
1871 	   {
1872 	      /* The decomposition statistics were computed during transformation of the decomposition store.
1873 	       * Since propagators can have changed the number of constraints/variables,
1874 	       * the statistics are no longer up-to-date and have to be recomputed.
1875 	       */
1876 	      SCIP_CALL( SCIPcomputeDecompStats(scip, decomp, TRUE) );
1877 	      nblocks = SCIPdecompGetNBlocks(decomp);
1878 	   }
1879 	
1880 	   /* reset parameters */
1881 	   SCIP_CALL( SCIPsetIntParam(scip, "decomposition/maxgraphedge", maxgraphedge) );
1882 	   SCIP_CALL( SCIPsetBoolParam(scip, "decomposition/disablemeasures", disablemeasures) );
1883 	
1884 	#ifdef SCIP_DEBUG
1885 	      char buffer[SCIP_MAXSTRLEN];
1886 	      SCIPdebugMsg(scip, "DPS used decomposition:\n%s\n", SCIPdecompPrintStats(decomp, buffer));
1887 	#endif
1888 	
1889 	   /* check properties of used decomposition */
1890 	   if( heurdata->maxlinkscore != 1.0 )
1891 	   {
1892 	      SCIP_Real linkscore;
1893 	      int nlinkconss;
1894 	
1895 	      nlinkconss = SCIPdecompGetNBorderConss(decomp);
1896 	
1897 	      linkscore = (SCIP_Real)nlinkconss / (SCIP_Real)nconss;
1898 	
1899 	      if( linkscore > heurdata->maxlinkscore )
1900 	      {
1901 	         SCIPdebugMsg(scip, "decomposition has not the required properties\n");
1902 	         goto TERMINATE;
1903 	      }
1904 	   }
1905 	
1906 	   if( sortedvarlabels[0] == SCIP_DECOMP_LINKVAR ||
1907 	         sortedconslabels[0] != SCIP_DECOMP_LINKCONS ||
1908 	         nblocks <= 1 )
1909 	   {
1910 	      SCIPdebugMsg(scip, "Problem has linking variables or no linking constraints or less than two blocks\n");
1911 	      goto TERMINATE;
1912 	   }
1913 	
1914 	   /* initialize heurdata */
1915 	   heurdata->linkingconss = sortedconss;
1916 	   heurdata->nlinking = SCIPdecompGetNBorderConss(decomp);
1917 	   heurdata->nblocks = nblocks;
1918 	
1919 	   /* allocate memory for blockproblems and initialize partially */
1920 	   SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem, nblocks) );
1921 	   for( b = 0; b < nblocks; b++ )
1922 	   {
1923 	      SCIP_CALL( SCIPallocBlockMemory(scip, &(blockproblem[b])) ); /*lint !e866*/
1924 	      SCIP_CALL( createSubscip(scip, &blockproblem[b]->blockscip) );
1925 	
1926 	      SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->linkingconss, heurdata->nlinking) );
1927 	      SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->linkingindices, heurdata->nlinking) );
1928 	      SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->slackvars, heurdata->nlinking * 2) ); /* maximum two slacks per linking constraint */
1929 	      SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->origobj, nvars) );
1930 	      blockproblem[b]->nblockvars = 0;
1931 	      blockproblem[b]->nlinking = 0;
1932 	      blockproblem[b]->nslackvars = 0;
1933 	   }
1934 	
1935 	   /* allocate memory for linkings and initialize partially */
1936 	   SCIP_CALL( SCIPallocBufferArray(scip, &linkings, heurdata->nlinking) );
1937 	   for( c = 0; c < heurdata->nlinking; c++ )
1938 	   {
1939 	      SCIP_CALL( SCIPallocBlockMemory(scip, &(linkings[c])) ); /*lint !e866*/
1940 	      SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->blockconss, heurdata->nblocks) );
1941 	      SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->slacks, heurdata->nblocks*2) ); /* maximum two slacks per block */
1942 	      SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->blocknumbers, heurdata->nblocks) );
1943 	      SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->minactivity, heurdata->nblocks) );
1944 	      SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->maxactivity, heurdata->nblocks) );
1945 	
1946 	      linkings[c]->linkingcons = heurdata->linkingconss[c];
1947 	      linkings[c]->currentrhs = NULL;
1948 	      linkings[c]->currentlhs = NULL;
1949 	      linkings[c]->nblocks = 0;
1950 	      linkings[c]->nslacks = 0;
1951 	      linkings[c]->nslacksperblock = 0;
1952 	      linkings[c]->lastviolations = 0;
1953 	      linkings[c]->hasrhs = FALSE;
1954 	      linkings[c]->haslhs = FALSE;
1955 	   }
1956 	
1957 	   SCIP_CALL( createAndSplitProblem(scip, heurdata, decomp, blockproblem, linkings, sortedvars, sortedconss, &success) );
1958 	   if( !success )
1959 	   {
1960 	      SCIPdebugMsg(scip, "Create and split problem failed\n");
1961 	      goto TERMINATE;
1962 	   }
1963 	
1964 	   /* allocate memory for current partition*/
1965 	   for( c = 0; c < heurdata->nlinking; c++ )
1966 	   {
1967 	      if( linkings[c]->hasrhs )
1968 	      {
1969 	         SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->currentrhs, linkings[c]->nblocks ) );
1970 	      }
1971 	
1972 	      if( linkings[c]->haslhs )
1973 	      {
1974 	         SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->currentlhs, linkings[c]->nblocks ) );
1975 	      }
1976 	   }
1977 	
1978 	   /* initialize partition */
1979 	   SCIP_CALL( initCurrent(scip, linkings, blockproblem, heurtiming, heurdata->nlinking, &success) );
1980 	   if( !success )
1981 	   {
1982 	      SCIPdebugMsg(scip, "Initialization of partition failed\n");
1983 	      goto TERMINATE;
1984 	   }
1985 	
1986 	   /* ------------------------------------------------------------------------ */
1987 	   SCIPdebugMsg(scip, "Start heuristik DPS\n");
1988 	   *result = SCIP_DIDNOTFIND;
1989 	
1990 	   for( k = 0; k < heurdata->maxit; k++ )
1991 	   {
1992 	      /* do not exceed the timelimit */
1993 	      SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1994 	      if( (timelimit - SCIPgetSolvingTime(scip)) <= 0 )
1995 	         goto TERMINATE;
1996 	
1997 	      /* solve the subproblems */
1998 	      allslacksval = 0.0;
1999 	      for( b = 0; b < heurdata->nblocks; b++ )
2000 	      {
2001 	         SCIP* subscip;
2002 	         subscip = blockproblem[b]->blockscip;
2003 	
2004 	         /* update time and memory limit of subproblem */
2005 	         SCIP_CALL( SCIPcopyLimits(scip, subscip) );
2006 	
2007 	         /* create event handler for LP events */
2008 	         if( k == 0 )
2009 	         {
2010 	            SCIP_CALL( SCIPincludeEventhdlrBasic(subscip, &eventhdlr, EVENTHDLR_NAME, EVENTHDLR_DESC, eventExecDps, NULL) );
2011 	            if( eventhdlr == NULL )
2012 	            {
2013 	               SCIPerrorMessage("event handler for " HEUR_NAME " heuristic not found.\n");
2014 	               return SCIP_PLUGINNOTFOUND;
2015 	            }
2016 	         }
2017 	
2018 	         /* catch LP events of sub-SCIP */
2019 	         SCIP_CALL( SCIPtransformProb(subscip) );
2020 	         SCIP_CALL( SCIPcatchEvent(subscip, SCIP_EVENTTYPE_LPSOLVED, eventhdlr, (SCIP_EVENTDATA*) heurdata, NULL) );
2021 	
2022 	         SCIPdebugMsg(scip, "Solve blockproblem %d\n", b);
2023 	         SCIP_CALL_ABORT( SCIPsolve(subscip) );
2024 	
2025 	         /* drop LP events of sub-SCIP */
2026 	         SCIP_CALL( SCIPdropEvent(subscip, SCIP_EVENTTYPE_LPSOLVED, eventhdlr, (SCIP_EVENTDATA*) heurdata, -1) );
2027 	
2028 	         /* get status and objective value if available */
2029 	         status = SCIPgetStatus(subscip);
2030 	         if( status == SCIP_STATUS_INFEASIBLE )
2031 	         {
2032 	            SCIPdebugMsg(scip, "Subproblem is infeasible\n");
2033 	            goto TERMINATE;
2034 	         }
2035 	         else if( status == SCIP_STATUS_UNBOUNDED )
2036 	         {
2037 	            SCIPdebugMsg(scip, "Subproblem is unbounded\n");
2038 	            goto TERMINATE;
2039 	         }
2040 	         else if( SCIPgetNSols(subscip) >= 1 )
2041 	         {
2042 	            blocksolval = SCIPgetPrimalbound(subscip);
2043 	
2044 	            if( status == SCIP_STATUS_TIMELIMIT && !SCIPisZero(scip, blocksolval) )
2045 	            {
2046 	               SCIPdebugMsg(scip, "Subproblem reached timelimit without optimal solution\n");
2047 	               goto TERMINATE;
2048 	            }
2049 	            SCIPdebugMsg(scip, "Solution value: %f\n", blocksolval);
2050 	            allslacksval += blocksolval;
2051 	         }
2052 	         else
2053 	         {
2054 	            SCIPdebugMsg(scip, "No subproblem solution available\n");
2055 	            goto TERMINATE;
2056 	         }
2057 	      }
2058 	
2059 	      /* all slack variables are zero -> we found a feasible solution */
2060 	      if( SCIPisZero(scip, allslacksval) )
2061 	      {
2062 	         SCIP_SOL* newsol;
2063 	
2064 	         SCIPdebugMsg(scip, "Feasible solution found after %i iterations\n", k);
2065 	
2066 	         /* create new solution */
2067 	         SCIP_CALL( SCIPcreateSol(scip, &newsol, heur) );
2068 	         for( b = 0; b < heurdata->nblocks; b++ )
2069 	         {
2070 	            SCIP_SOL* blocksol;
2071 	            SCIP_VAR** blockvars;
2072 	            SCIP_Real* blocksolvals;
2073 	            int nblockvars;
2074 	
2075 	            /* get solution of block variables (without slack variables) */
2076 	            blocksol = SCIPgetBestSol(blockproblem[b]->blockscip);
2077 	            blockvars = SCIPgetOrigVars(blockproblem[b]->blockscip);
2078 	            nblockvars = blockproblem[b]->nblockvars;
2079 	
2080 	            SCIP_CALL( SCIPallocBufferArray(scip, &blocksolvals, nblockvars) );
2081 	            SCIP_CALL( SCIPgetSolVals(blockproblem[b]->blockscip, blocksol, nblockvars, blockvars, blocksolvals) );
2082 	
2083 	            for( c = 0; c < nblockvars; c++ )
2084 	            {
2085 	               SCIP_VAR* origvar;
2086 	
2087 	               origvar = SCIPfindVar(scip, SCIPvarGetName(blockvars[c]));
2088 	               SCIP_CALL( SCIPsetSolVal(scip, newsol, origvar, blocksolvals[c]) );
2089 	            }
2090 	
2091 	            SCIPfreeBufferArray(scip, &blocksolvals);
2092 	         }
2093 	
2094 	         /* if reoptimization is activated, fix partition and reoptimize with original objective function */
2095 	         if( heurdata->reoptimize )
2096 	         {
2097 	            SCIP_SOL* improvedsol = NULL;
2098 	            SCIP_CALL( reoptimize(scip, heur, newsol, blockproblem, heurdata->nblocks, heurdata->reoptlimits, &improvedsol, &success) );
2099 	            assert(improvedsol != NULL || success == FALSE);
2100 	
2101 	            if( success )
2102 	            {
2103 	               SCIP_CALL( SCIPtrySolFree(scip, &improvedsol, TRUE, FALSE, TRUE, TRUE, TRUE, &success) );
2104 	               if( success )
2105 	               {
2106 	                  SCIPdebugMsg(scip, "Reoptimizing solution successful\n");
2107 	                  *result = SCIP_FOUNDSOL;
2108 	               }
2109 	            }
2110 	         }
2111 	
2112 	         /* if reoptimization is turned off or reoptimization found no solution, try initial solution */
2113 	         if( *result != SCIP_FOUNDSOL )
2114 	         {
2115 	            SCIPdebugMsg(scip, "Solution has value: %.2f\n", SCIPgetSolOrigObj(scip, newsol));
2116 	            SCIP_CALL( SCIPtrySolFree(scip, &newsol, TRUE, FALSE, TRUE, TRUE, TRUE, &success) );
2117 	            if( success )
2118 	            {
2119 	               SCIPdebugMsg(scip, "Solution copy successful\n");
2120 	               *result = SCIP_FOUNDSOL;
2121 	            }
2122 	         }
2123 	         else
2124 	         {
2125 	            SCIP_CALL( SCIPfreeSol(scip, &newsol) );
2126 	         }
2127 	
2128 	         goto TERMINATE;
2129 	      }
2130 	      /* current partition is not feasible:
2131 	       * - update partition
2132 	       * - update penalty parameters lambda
2133 	       * - reuse last solution
2134 	       */
2135 	      else
2136 	      {
2137 	         int* nviolatedblocksrhs; /* number of blocks which violate rhs for each linking constraint */
2138 	         int* nviolatedblockslhs; /* number of blocks which violate lhs for each linking constraint */
2139 	         SCIP_Bool update = FALSE;
2140 	
2141 	         SCIP_CALL( SCIPallocBufferArray(scip, &nviolatedblocksrhs, heurdata->nlinking) );
2142 	         SCIP_CALL( SCIPallocBufferArray(scip, &nviolatedblockslhs, heurdata->nlinking) );
2143 	         BMSclearMemoryArray(nviolatedblocksrhs, heurdata->nlinking);
2144 	         BMSclearMemoryArray(nviolatedblockslhs, heurdata->nlinking);
2145 	
2146 	         SCIPdebugMsg(scip, "update partition\n");
2147 	         SCIP_CALL( updatePartition(scip, linkings, blockproblem, &nviolatedblocksrhs, &nviolatedblockslhs,
2148 	                                    heurdata->nlinking, nblocks, k, &update) );
2149 	
2150 	         /* terminate if partition is not updated */
2151 	         if( !update )
2152 	         {
2153 	            SCIPfreeBufferArray(scip, &nviolatedblockslhs);
2154 	            SCIPfreeBufferArray(scip, &nviolatedblocksrhs);
2155 	            goto TERMINATE;
2156 	         }
2157 	
2158 	         /* in order to change objective function */
2159 	         for( b = 0; b < heurdata->nblocks; b++ )
2160 	         {
2161 	            SCIP_CALL( SCIPfreeTransform(blockproblem[b]->blockscip) );
2162 	         }
2163 	
2164 	         SCIPdebugMsg(scip, "update lambda\n");
2165 	         SCIP_CALL( updateLambda(scip, heurdata, linkings, blockproblem, nviolatedblocksrhs, nviolatedblockslhs, heurdata->nlinking) );
2166 	
2167 	         if( heurdata->reuse )
2168 	         {
2169 	            /* reuse old solution in each block if available */
2170 	            SCIP_CALL( reuseSolution(linkings, blockproblem, nblocks) );
2171 	         }
2172 	
2173 	         SCIPfreeBufferArray(scip, &nviolatedblockslhs);
2174 	         SCIPfreeBufferArray(scip, &nviolatedblocksrhs);
2175 	      }
2176 	   }
2177 	   SCIPdebugMsg(scip, "maximum number of iterations reached\n");
2178 	
2179 	   /* ------------------------------------------------------------------------ */
2180 	   /* free memory */
2181 	TERMINATE:
2182 	   if( linkings != NULL )
2183 	   {
2184 	      for( c = heurdata->nlinking - 1; c >= 0; c-- )
2185 	      {
2186 	         if( linkings[c]->currentlhs != NULL )
2187 	            SCIPfreeBufferArray(scip, &(linkings[c])->currentlhs);
2188 	
2189 	         if( linkings[c]->currentrhs != NULL )
2190 	            SCIPfreeBufferArray(scip, &(linkings[c])->currentrhs);
2191 	      }
2192 	
2193 	      for( c = heurdata->nlinking - 1; c >= 0; c-- )
2194 	      {
2195 	         linkings[c]->linkingcons = NULL;
2196 	         SCIPfreeBufferArray(scip, &(linkings[c])->maxactivity);
2197 	         SCIPfreeBufferArray(scip, &(linkings[c])->minactivity);
2198 	         SCIPfreeBufferArray(scip, &(linkings[c])->blocknumbers);
2199 	         SCIPfreeBufferArray(scip, &(linkings[c])->slacks);
2200 	         SCIPfreeBufferArray(scip, &(linkings[c])->blockconss);
2201 	         SCIPfreeBlockMemory(scip, &(linkings[c])); /*lint !e866*/
2202 	      }
2203 	      SCIPfreeBufferArray(scip, &linkings);
2204 	   }
2205 	
2206 	   if( blockproblem != NULL )
2207 	   {
2208 	      for( b = nblocks - 1; b >= 0; b-- )
2209 	      {
2210 	         SCIPfreeBufferArray(scip, &(blockproblem[b])->origobj);
2211 	         SCIPfreeBufferArray(scip, &(blockproblem[b])->slackvars);
2212 	         SCIPfreeBufferArray(scip, &(blockproblem[b])->linkingindices);
2213 	         SCIPfreeBufferArray(scip, &(blockproblem[b])->linkingconss);
2214 	         SCIP_CALL( SCIPfree(&blockproblem[b]->blockscip) );
2215 	         SCIPfreeBlockMemory(scip, &(blockproblem[b])); /*lint !e866*/
2216 	      }
2217 	      SCIPfreeBufferArray(scip, &blockproblem);
2218 	   }
2219 	
2220 	   if( assigneddecomp != NULL )
2221 	      SCIPdecompFree(&assigneddecomp, SCIPblkmem(scip));
2222 	
2223 	   if( sortedconslabels != NULL )
2224 	      SCIPfreeBufferArray(scip, &sortedconslabels);
2225 	
2226 	   if( sortedvarlabels != NULL )
2227 	      SCIPfreeBufferArray(scip, &sortedvarlabels);
2228 	
2229 	   if( sortedconss != NULL )
2230 	      SCIPfreeBufferArray(scip, &sortedconss);
2231 	
2232 	   if( sortedvars != NULL )
2233 	      SCIPfreeBufferArray(scip, &sortedvars);
2234 	
2235 	   SCIPdebugMsg(scip, "Leave DPS heuristic\n");
2236 	
2237 	   return SCIP_OKAY;
2238 	}
2239 	
2240 	
2241 	/*
2242 	 * primal heuristic specific interface methods
2243 	 */
2244 	
2245 	/** creates the dps primal heuristic and includes it in SCIP */
2246 	SCIP_RETCODE SCIPincludeHeurDps(
2247 	   SCIP*                 scip                /**< SCIP data structure */
2248 	   )
2249 	{
2250 	   SCIP_HEURDATA* heurdata;
2251 	   SCIP_HEUR* heur;
2252 	
2253 	   /* create dps primal heuristic data */
2254 	   SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
2255 	
2256 	   heur = NULL;
2257 	
2258 	   /* include primal heuristic */
2259 	   SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
2260 	         HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
2261 	         HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecDps, heurdata) );
2262 	
2263 	   assert(heur != NULL);
2264 	
2265 	   /* set non fundamental callbacks via setter functions */
2266 	   SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyDps) );
2267 	   SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeDps) );
2268 	
2269 	   /* add dps primal heuristic parameters */
2270 	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxiterations",
2271 	   "maximal number of iterations", &heurdata->maxit, FALSE, DEFAULT_MAXIT, 1, INT_MAX, NULL, NULL) );
2272 	
2273 	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxlinkscore",
2274 	   "maximal linking score of used decomposition (equivalent to percentage of linking constraints)",
2275 	   &heurdata->maxlinkscore, FALSE, 1.0, 0.0, 1.0, NULL, NULL) );
2276 	
2277 	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/penalty",
2278 	   "multiplier for absolute increase of penalty parameters (0: no increase)",
2279 	   &heurdata->penalty, FALSE, DEFAULT_PENALTY, 0.0, SCIP_REAL_MAX, NULL, NULL) );
2280 	
2281 	   SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reoptimize",
2282 	   "should the problem get reoptimized with the original objective function?", &heurdata->reoptimize, FALSE, FALSE, NULL, NULL) );
2283 	
2284 	   SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reuse",
2285 	   "should solutions get reused in subproblems?", &heurdata->reuse, FALSE, FALSE, NULL, NULL) );
2286 	
2287 	   SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reoptlimits",
2288 	   "should strict limits for reoptimization be set?", &heurdata->reoptlimits, FALSE, TRUE, NULL, NULL) );
2289 	
2290 	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/timing",
2291 	      "should the heuristic run before or after the processing of the node? (0: before, 1: after, 2: both)",
2292 	      &heurdata->timing, FALSE, 0, 0, 2, NULL, NULL) );
2293 	
2294 	   return SCIP_OKAY;
2295 	}
2296