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   scip_dcmp.c
26   	 * @ingroup OTHER_CFILES
27   	 * @brief  methods for working with decompositions
28   	 * @author Gregor Hendel
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include <assert.h>
34   	#include "scip/struct_dcmp.h"
35   	#include "scip/debug.h"
36   	#include "scip/dcmp.h"
37   	#include "scip/mem.h"
38   	#include "scip/pub_misc.h"
39   	#include "scip/pub_var.h"
40   	#include "scip/scip_cons.h"
41   	#include "scip/scip_prob.h"
42   	#include "scip/scip_var.h"
43   	#include "scip/scip_mem.h"
44   	#include "scip/struct_scip.h"
45   	#include "scip/pub_cons.h"
46   	#include "scip/pub_dcmp.h"
47   	#include "scip/scip_dcmp.h"
48   	#include "scip/scip_general.h"
49   	#include "scip/scip_param.h"
50   	#include "scip/scip_var.h"
51   	#include "scip/scip_datastructures.h"
52   	#include "scip/scip_message.h"
53   	
54   	
55   	#define LABEL_UNASSIGNED INT_MIN /* label constraints or variables as unassigned. Only for internal use */
56   	
57   	/** count occurrences of label in array, starting from pos */
58   	static
59   	int countLabelFromPos(
60   	   int*                  labels,             /**< array of labels */
61   	   int                   pos,                /**< position to start counting from */
62   	   int                   nlabels             /**< the number of labels */
63   	   )
64   	{
65   	   int endpos = pos;
66   	   int currlabel;
67   	
68   	   assert(labels != NULL);
69   	   assert(pos < nlabels);
70   	
71   	   currlabel = labels[pos];
72   	
73   	   do
74   	   {
75   	      endpos++;
76   	   }
77   	   while( endpos < nlabels && labels[endpos] == currlabel );
78   	
79   	   return endpos - pos;
80   	}
81   	
82   	/** raises an error if the condition is not TRUE */
83   	static
84   	SCIP_RETCODE ensureCondition(
85   	   SCIP_Bool             condition           /**< some condition that must hold */
86   	   )
87   	{
88   	   return condition ? SCIP_OKAY : SCIP_ERROR;
89   	}
90   	
91   	/** get variable buffer storage size for the buffer to be large enough to hold all variables */
92   	static
93   	int getVarbufSize(
94   	   SCIP*                 scip                /**< SCIP data structure */
95   	   )
96   	{
97   	   int norigvars;
98   	   int ntransvars;
99   	
100  	   norigvars = SCIPgetNOrigVars(scip);
101  	   ntransvars = SCIPgetNVars(scip);
102  	
103  	   return 2 * MAX(norigvars, ntransvars);
104  	}
105  	
106  	/** get variables and constraints of the original or transformed problem, to which the decomposition corresponds */
107  	static
108  	void getDecompVarsConssData(
109  	   SCIP*                 scip,               /**< SCIP data structure */
110  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
111  	   SCIP_VAR***           vars,               /**< pointer to store original/transformed variables array, or NULL */
112  	   SCIP_CONS***          conss,              /**< pointer to store original/transformed constraints array, or NULL */
113  	   int*                  nvars,              /**< pointer to store original/transformed variables array's length, or NULL */
114  	   int*                  nconss              /**< pointer to store original/transformed constraints array's length, or NULL */
115  	   )
116  	{
117  	   SCIP_Bool original;
118  	   assert(scip != NULL);
119  	   assert(decomp != NULL);
120  	
121  	   original = SCIPdecompIsOriginal(decomp);
122  	
123  	   if( vars )
124  	      *vars = original ? SCIPgetOrigVars(scip) : SCIPgetVars(scip);
125  	
126  	   if( nvars )
127  	         *nvars = original ? SCIPgetNOrigVars(scip) : SCIPgetNVars(scip);
128  	
129  	   if( conss )
130  	      *conss = original ? SCIPgetOrigConss(scip) : SCIPgetConss(scip);
131  	
132  	   if( nconss )
133  	      *nconss = original ? SCIPgetNOrigConss(scip) : SCIPgetNConss(scip);
134  	}
135  	
136  	/** query the constraints active variables and their labels */
137  	static
138  	SCIP_RETCODE decompGetConsVarsAndLabels(
139  	   SCIP*                 scip,               /**< SCIP data structure */
140  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
141  	   SCIP_CONS*            cons,               /**< the constraint */
142  	   SCIP_VAR**            varbuf,             /**< variable buffer array */
143  	   int*                  labelbuf,           /**< buffer to store labels, or NULL if not needed */
144  	   int                   bufsize,            /**< size of buffer arrays */
145  	   int*                  nvars,              /**< pointer to store number of variables */
146  	   int*                  requiredsize,       /**< pointer to store required size */
147  	   SCIP_Bool*            success             /**< pointer to store whether variables and labels were successfully inserted */
148  	   )
149  	{
150  	   SCIP_Bool success2;
151  	
152  	   assert(scip != NULL);
153  	   assert(decomp != NULL);
154  	   assert(cons != NULL);
155  	   assert(varbuf != NULL);
156  	   assert(nvars != NULL);
157  	   assert(requiredsize != NULL);
158  	   assert(success != NULL);
159  	
160  	   *success = FALSE;
161  	   *requiredsize = 0;
162  	   *nvars = 0;
163  	   SCIP_CALL( SCIPgetConsNVars(scip, cons, nvars, &success2) );
164  	
165  	   /* the constraint does not have the corresponding callback */
166  	   if( ! success2 )
167  	   {
168  	      return SCIP_OKAY;
169  	   }
170  	
171  	   if( bufsize < *nvars )
172  	   {
173  	      *requiredsize = *nvars;
174  	
175  	      return SCIP_OKAY;
176  	   }
177  	
178  	   SCIP_CALL( SCIPgetConsVars(scip, cons, varbuf, bufsize, &success2) );
179  	   /* the constraint does not have the corresponding callback */
180  	   if( ! success2 )
181  	   {
182  	      return SCIP_OKAY;
183  	   }
184  	
185  	   if( ! SCIPdecompIsOriginal(decomp) )
186  	   {
187  	      SCIP_CALL( SCIPgetActiveVars(scip, varbuf, nvars, bufsize, requiredsize) );
188  	
189  	      if( *requiredsize > bufsize )
190  	         return SCIP_OKAY;
191  	   }
192  	   else
193  	   {
194  	      int v;
195  	      for( v = 0; v < *nvars; ++v )
196  	      {
197  	         assert(SCIPvarIsActive(varbuf[v]) || SCIPvarIsNegated(varbuf[v]));
198  	
199  	         /* some constraint handlers such as indicator may already return inactive variables */
200  	         if( SCIPvarIsNegated(varbuf[v]) )
201  	            varbuf[v] = SCIPvarGetNegatedVar(varbuf[v]);
202  	      }
203  	   }
204  	
205  	   /* get variables labels, if requested */
206  	   if( labelbuf != NULL )
207  	   {
208  	      SCIPdecompGetVarsLabels(decomp, varbuf, labelbuf, *nvars);
209  	   }
210  	
211  	   *success = TRUE;
212  	
213  	   return SCIP_OKAY;
214  	}
215  	
216  	/** creates a decomposition */
217  	SCIP_RETCODE SCIPcreateDecomp(
218  	   SCIP*                 scip,               /**< SCIP data structure */
219  	   SCIP_DECOMP**         decomp,             /**< pointer to store the decomposition data structure */
220  	   int                   nblocks,            /**< the number of blocks (without the linking block) */
221  	   SCIP_Bool             original,           /**< is this a decomposition in the original (TRUE) or transformed space? */
222  	   SCIP_Bool             benderslabels       /**< should the variables be labeled for the application of Benders' decomposition */
223  	   )
224  	{
225  	   assert(scip != NULL);
226  	
227  	   SCIP_CALL( SCIPdecompCreate(decomp, SCIPblkmem(scip), nblocks, original, benderslabels) );
228  	
229  	   return SCIP_OKAY;
230  	}
231  	
232  	/** frees a decomposition */
233  	void SCIPfreeDecomp(
234  	   SCIP*                 scip,               /**< SCIP data structure */
235  	   SCIP_DECOMP**         decomp              /**< pointer to free the decomposition data structure */
236  	   )
237  	{
238  	   assert(scip != NULL);
239  	
240  	   SCIPdecompFree(decomp, SCIPblkmem(scip));
241  	}
242  	
243  	/** adds decomposition to SCIP */
244  	SCIP_RETCODE SCIPaddDecomp(
245  	   SCIP*                 scip,               /**< SCIP data structure */
246  	   SCIP_DECOMP*          decomp              /**< decomposition to add */
247  	   )
248  	{
249  	   assert(scip != NULL);
250  	   assert(decomp != NULL);
251  	
252  	   SCIP_CALL( SCIPcheckStage(scip, "SCIPaddDecomp", FALSE, SCIPdecompIsOriginal(decomp), SCIPdecompIsOriginal(decomp),
253  	      SCIPdecompIsOriginal(decomp), SCIPdecompIsOriginal(decomp), TRUE, TRUE, TRUE, !SCIPdecompIsOriginal(decomp),
254  	      !SCIPdecompIsOriginal(decomp), !SCIPdecompIsOriginal(decomp), FALSE, FALSE, FALSE) );
255  	
256  	   SCIP_CALL( SCIPdecompstoreAdd(scip->decompstore, decomp) );
257  	
258  	   return SCIP_OKAY;
259  	}
260  	
261  	/** gets available user decompositions for either the original or transformed problem */
262  	void SCIPgetDecomps(
263  	   SCIP*                 scip,               /**< SCIP data structure */
264  	   SCIP_DECOMP***        decomps,            /**< pointer to store decompositions array */
265  	   int*                  ndecomps,           /**< pointer to store number of decompositions */
266  	   SCIP_Bool             original            /**< should the decompositions for the original problem be returned? */
267  	   )
268  	{
269  	   assert(scip != NULL);
270  	
271  	   SCIP_CALL_ABORT( SCIPcheckStage(scip, "SCIPgetDecomps", FALSE, original, original, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE) );
272  	
273  	   if( decomps != NULL )
274  	      *decomps = original ? SCIPdecompstoreGetOrigDecomps(scip->decompstore) : SCIPdecompstoreGetDecomps(scip->decompstore);
275  	
276  	   if( ndecomps != NULL )
277  	      *ndecomps = original ? SCIPdecompstoreGetNOrigDecomps(scip->decompstore) : SCIPdecompstoreGetNDecomps(scip->decompstore);
278  	}
279  	
280  	/** returns TRUE if the constraint \p cons contains only linking variables in decomposition \p decomp */
281  	SCIP_RETCODE SCIPhasConsOnlyLinkVars(
282  	   SCIP*                 scip,               /**< SCIP data structure */
283  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
284  	   SCIP_CONS*            cons,               /**< the constraint */
285  	   SCIP_Bool*            hasonlylinkvars     /**< will be set to TRUE if this constraint has only linking variables */
286  	   )
287  	{
288  	   SCIP_VAR** consvars;
289  	   int nvars;
290  	   int i;
291  	   SCIP_Bool success;
292  	
293  	   assert(scip != NULL);
294  	   assert(cons != NULL);
295  	   assert(decomp != NULL);
296  	   assert(hasonlylinkvars != NULL);
297  	
298  	   SCIP_CALL( SCIPgetConsNVars(scip, cons, &nvars, &success) );
299  	   SCIP_CALL( ensureCondition(success) );
300  	
301  	   SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nvars) );
302  	
303  	   SCIP_CALL( SCIPgetConsVars(scip, cons, consvars, nvars, &success) );
304  	   SCIP_CALL( ensureCondition(success) );
305  	
306  	   if( ! SCIPdecompIsOriginal(decomp) )
307  	   {
308  	      int requiredsize;
309  	      SCIP_CALL( SCIPgetActiveVars(scip, consvars, &nvars, nvars, &requiredsize) );
310  	      assert(requiredsize <= nvars);
311  	   }
312  	
313  	   *hasonlylinkvars = TRUE;
314  	   /* check if variables are all linking variables */
315  	   for( i = 0; i < nvars && *hasonlylinkvars; ++i )
316  	   {
317  	      int label;
318  	
319  	      SCIPdecompGetVarsLabels(decomp, &consvars[i], &label, 1);
320  	
321  	      *hasonlylinkvars = (label == SCIP_DECOMP_LINKVAR);
322  	   }
323  	
324  	   SCIPfreeBufferArray(scip, &consvars);
325  	
326  	   return SCIP_OKAY;
327  	}
328  	
329  	
330  	/** computes constraint labels from variable labels
331  	 *
332  	 *  Existing labels for the constraints are simply overridden
333  	 *
334  	 *  The computed labels depend on the flag SCIPdecompUseBendersLabels() of the decomposition. If the flag is set
335  	 *  to FALSE, the labeling assigns
336  	 *
337  	 *  - label i, if only variables labeled i are present in the constraint (and optionally linking variables)
338  	 *  - SCIP_DECOMP_LINKCONS, if there are either only variables labeled with SCIP_DECOMP_LINKVAR present, or
339  	 *    if there are variables with more than one block label.
340  	 *
341  	 *  If the flag is set to TRUE, the assignment is the same, unless variables from 2 named blocks occur in the same
342  	 *  constraint, which is an invalid labeling for the Benders case.
343  	 */
344  	SCIP_RETCODE SCIPcomputeDecompConsLabels(
345  	   SCIP*                 scip,               /**< SCIP data structure */
346  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
347  	   SCIP_CONS**           conss,              /**< array of constraints */
348  	   int                   nconss              /**< number of constraints */
349  	   )
350  	{
351  	   SCIP_VAR** varbuffer;
352  	   int c;
353  	   int varbufsize;
354  	   int* varlabels;
355  	   int* conslabels;
356  	   SCIP_Bool benderserror;
357  	   SCIP_Bool benderslabels;
358  	
359  	   assert(decomp != NULL);
360  	
361  	   if( nconss == 0 )
362  	      return SCIP_OKAY;
363  	
364  	   varbufsize = getVarbufSize(scip);
365  	
366  	   SCIP_CALL( SCIPallocBufferArray(scip, &varbuffer, varbufsize) );
367  	   SCIP_CALL( SCIPallocBufferArray(scip, &varlabels, varbufsize) );
368  	   SCIP_CALL( SCIPallocBufferArray(scip, &conslabels, nconss) );
369  	
370  	   benderslabels = SCIPdecompUseBendersLabels(decomp);
371  	   benderserror = FALSE;
372  	
373  	   /* assign label to each individual constraint */
374  	   for( c = 0; c < nconss && ! benderserror; ++c )
375  	   {
376  	      int nconsvars;
377  	      int v;
378  	      int nlinkingvars = 0;
379  	      int conslabel;
380  	      int requiredsize;
381  	
382  	      SCIP_Bool success;
383  	
384  	      SCIP_CALL( decompGetConsVarsAndLabels(scip, decomp, conss[c], varbuffer, varlabels,
385  	            varbufsize, &nconsvars, &requiredsize, &success) );
386  	      SCIP_CALL( ensureCondition(success) );
387  	
388  	      /* loop over variable labels to compute the constraint label */
389  	      conslabel = LABEL_UNASSIGNED;
390  	      for( v = 0; v < nconsvars; ++v )
391  	      {
392  	         int varlabel = varlabels[v];
393  	
394  	         /* count the number of linking variables, and keep track if there are two variables with different labels */
395  	         if( varlabel == SCIP_DECOMP_LINKVAR )
396  	            ++nlinkingvars;
397  	         else if( conslabel == LABEL_UNASSIGNED )
398  	            conslabel = varlabel;
399  	         else if( conslabel != varlabel )
400  	         {
401  	            /* there must not be two variables from different named blocks in a single constraint, since the presence
402  	             * of named block variables forbids this constraint from the master (linking) block
403  	             */
404  	            if( benderslabels )
405  	               benderserror = TRUE;
406  	
407  	            conslabel = SCIP_DECOMP_LINKCONS;
408  	            break;
409  	         }
410  	      }
411  	
412  	      assert(nlinkingvars == nconsvars || conslabel != LABEL_UNASSIGNED);
413  	      assert(v == nconsvars || conslabel == SCIP_DECOMP_LINKCONS);
414  	      SCIP_UNUSED(nlinkingvars);
415  	
416  	      /* if there are only linking variables, the constraint is unassigned */
417  	      if( conslabel == LABEL_UNASSIGNED )
418  	         conslabel = SCIP_DECOMP_LINKCONS;
419  	      conslabels[c] = conslabel;
420  	   }
421  	
422  	   SCIP_CALL( SCIPdecompSetConsLabels(decomp, conss, conslabels, nconss) );
423  	
424  	   SCIPfreeBufferArray(scip, &conslabels);
425  	   SCIPfreeBufferArray(scip, &varlabels);
426  	   SCIPfreeBufferArray(scip, &varbuffer);
427  	
428  	   /* throw an error and inform the user if the variable block decomposition does not allow a benders constraint labeling */
429  	   if( benderserror )
430  	   {
431  	      SCIPerrorMessage("Error in constraint label computation; variables from multiple named blocks in a single constraint\n");
432  	
433  	      return SCIP_INVALIDDATA;
434  	   }
435  	
436  	   return SCIP_OKAY;
437  	}
438  	
439  	/** creates a decomposition of the variables from a labeling of the constraints
440  	 *
441  	 *  NOTE: by default, the variable labeling is based on a Dantzig-Wolfe decomposition. This means that constraints in named
442  	 *  blocks have have precedence over linking constraints. If a variable exists in constraints from
443  	 *  two or more named blocks, then this variable is marked as a linking variable.
444  	 *  If a variable occurs in exactly one named block i>=0, it is assigned label i.
445  	 *  Variables which are only in linking constraints are unlabeled. However, SCIPdecompGetVarsLabels() will
446  	 *  label them as linking variables.
447  	 *
448  	 *  If the variables should be labeled for the application of Benders' decomposition, the decomposition must be
449  	 *  flagged explicitly via SCIPdecompSetUseBendersLabels().
450  	 *  With this setting, the presence in linking constraints takes precedence over the presence in named blocks.
451  	 *  Now, a variable is considered linking if it is present in at least one linking constraint and an arbitrary
452  	 *  number of constraints from named blocks.
453  	 */
454  	SCIP_RETCODE SCIPcomputeDecompVarsLabels(
455  	   SCIP*                 scip,               /**< SCIP data structure */
456  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
457  	   SCIP_CONS**           conss,              /**< array of constraints */
458  	   int                   nconss              /**< number of constraints */
459  	   )
460  	{
461  	   int c;
462  	   int* conslabels;
463  	   SCIP_VAR** varbuffer;
464  	   int varbufsize;
465  	   SCIP_Bool benderslabels;
466  	
467  	   assert(scip != NULL);
468  	   assert(decomp != NULL);
469  	   assert(conss != NULL);
470  	
471  	   /* make the buffer array large enough such that we do not have to reallocate buffer */
472  	   varbufsize = getVarbufSize(scip);
473  	
474  	   /* allocate buffer to store constraint variables and labels */
475  	   SCIP_CALL( SCIPallocBufferArray(scip, &conslabels, nconss) );
476  	   SCIP_CALL( SCIPallocBufferArray(scip, &varbuffer, varbufsize) );
477  	
478  	   /* query constraint labels */
479  	   SCIPdecompGetConsLabels(decomp, conss, conslabels, nconss);
480  	
481  	   benderslabels = SCIPdecompUseBendersLabels(decomp);
482  	   /* iterate over constraints and query the corresponding constraint labels */
483  	   for( c = 0; c < nconss; ++c )
484  	   {
485  	      int conslabel;
486  	      int v;
487  	      int nconsvars;
488  	      SCIP_Bool success;
489  	      int requiredsize;
490  	      int newvarlabel;
491  	
492  	      conslabel = conslabels[c];
493  	
494  	      if( conslabel == SCIP_DECOMP_LINKCONS )
495  	      {
496  	         /* skip linking constraints unless Benders labeling is used */
497  	         if( ! benderslabels )
498  	            continue;
499  	         else
500  	            newvarlabel = SCIP_DECOMP_LINKVAR;
501  	      }
502  	      else
503  	         newvarlabel = conslabel;
504  	
505  	      SCIP_CALL( decompGetConsVarsAndLabels(scip, decomp, conss[c], varbuffer, NULL,
506  	            varbufsize, &nconsvars, &requiredsize, &success) );
507  	      SCIP_CALL( ensureCondition(success) );
508  	
509  	      /* each variable in this constraint gets the constraint label unless it already has a different label -> make it a linking variable */
510  	      for( v = 0; v < nconsvars; ++v )
511  	      {
512  	         SCIP_VAR* var = varbuffer[v];
513  	
514  	         assert(SCIPvarIsActive(var) || (SCIPdecompIsOriginal(decomp) && SCIPvarIsNegated(var)));
515  	
516  	         /* some constraint handlers such as indicator may already return inactive variables */
517  	         if( SCIPvarIsNegated(var) )
518  	            var = SCIPvarGetNegatedVar(var);
519  	
520  	         if( SCIPhashmapExists(decomp->var2block, (void *)var) )
521  	         {
522  	            int varlabel = SCIPhashmapGetImageInt(decomp->var2block, (void *)var);
523  	
524  	            /* store the label linking variable explicitly to distinguish it from the default */
525  	            if( varlabel != SCIP_DECOMP_LINKVAR && varlabel != newvarlabel )
526  	               SCIP_CALL( SCIPhashmapSetImageInt(decomp->var2block, (void *)var, SCIP_DECOMP_LINKVAR) );
527  	         }
528  	         else
529  	         {
530  	            SCIP_CALL( SCIPhashmapInsertInt(decomp->var2block, (void *)var, newvarlabel) );
531  	         }
532  	      }
533  	   }
534  	
535  	   SCIPfreeBufferArray(scip, &varbuffer);
536  	   SCIPfreeBufferArray(scip, &conslabels);
537  	
538  	   return SCIP_OKAY;
539  	}
540  	
541  	/** assigns linking constraints to blocks
542  	 *
543  	 * Each linking constraint is assigned to the most frequent block among its variables.
544  	 * Variables of other blocks are relabeled as linking variables.
545  	 * Constraints that have only linking variables are skipped.
546  	 *
547  	 * @note: In contrast to SCIPcomputeDecompConsLabels(), this method potentially relabels variables.
548  	 */
549  	SCIP_RETCODE SCIPassignDecompLinkConss(
550  	   SCIP*                 scip,               /**< SCIP data structure */
551  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
552  	   SCIP_CONS**           conss,              /**< array of linking constraints that should be reassigned */
553  	   int                   nconss,             /**< number of constraints */
554  	   int*                  nskipconss          /**< pointer to store the number of constraints that were skipped, or NULL */
555  	   )
556  	{
557  	   SCIP_VAR** vars;
558  	   SCIP_VAR** allvars;
559  	   int* varslabels;
560  	   int requiredsize;
561  	   int nvars;
562  	   int nconsvars;
563  	   int varbufsize;
564  	   int c;
565  	   int nskipconsslocal;
566  	   int defaultlabel;
567  	
568  	   assert(scip != NULL);
569  	   assert(decomp != NULL);
570  	
571  	   nvars = SCIPgetNVars(scip);
572  	   varbufsize = getVarbufSize(scip);
573  	
574  	   SCIP_CALL( SCIPallocBufferArray(scip, &varslabels, varbufsize) );
575  	   SCIP_CALL( SCIPallocBufferArray(scip, &vars, varbufsize) );
576  	
577  	   /* get one label as default label */
578  	   allvars = SCIPdecompIsOriginal(decomp) ? SCIPgetOrigVars(scip) : SCIPgetVars(scip);
579  	   SCIPdecompGetVarsLabels(decomp, allvars, varslabels, nvars);
580  	   for( c = 0; c < nvars; c++ )
581  	   {
582  	      if( varslabels[c] != SCIP_DECOMP_LINKVAR )
583  	      {
584  	         defaultlabel = varslabels[c];
585  	         break;
586  	      }
587  	   }
588  	
589  	   nskipconsslocal = 0;
590  	   for( c = 0; c < nconss; c++ )
591  	   {
592  	      SCIP_Bool success;
593  	
594  	      SCIP_CALL( decompGetConsVarsAndLabels(scip, decomp, conss[c], vars, varslabels, varbufsize,
595  	               &nconsvars, &requiredsize, &success) );
596  	      SCIP_CALL( ensureCondition(success) );
597  	
598  	      SCIPsortIntPtr(varslabels, (void **)vars, nconsvars);
599  	      /* constraint contains no (active) variables */
600  	      if( nconsvars == 0 )
601  	      {
602  	         SCIP_CALL( SCIPdecompSetConsLabels(decomp, &conss[c], &defaultlabel, 1) );
603  	      }
604  	      /* constraint contains only linking variables */
605  	      else if( varslabels[nconsvars - 1] == SCIP_DECOMP_LINKVAR )
606  	      {
607  	         nskipconsslocal++;
608  	
609  	         continue;
610  	      }
611  	      else
612  	      {
613  	         int startposs[2];
614  	         int endposs[2];
615  	         int nlinkvars;
616  	         int block;
617  	         int maxnblockvars;
618  	         int curr;
619  	         int v;
620  	         int p;
621  	
622  	         /* count linking variables */
623  	         if( varslabels[0] == SCIP_DECOMP_LINKVAR )
624  	         {
625  	            nlinkvars = countLabelFromPos(varslabels, 0, nconsvars);
626  	         }
627  	         else
628  	         {
629  	            nlinkvars = 0;
630  	         }
631  	
632  	         assert(nlinkvars < nconsvars);
633  	
634  	         curr = nlinkvars;
635  	         /* find the most frequent block label among the nonlinking variables */
636  	         maxnblockvars = 0;
637  	         block = -1;
638  	         do
639  	         {
640  	            int nblockvars = countLabelFromPos(varslabels, curr, nconsvars);
641  	            if( nblockvars > maxnblockvars )
642  	            {
643  	               maxnblockvars = nblockvars;
644  	               block = curr;
645  	            }
646  	            curr += nblockvars;
647  	         }
648  	         while( curr < nconsvars );
649  	
650  	         /* reassign all variables from other blocks as linking variables */
651  	         startposs[0] = nlinkvars;
652  	         endposs[0] = block;
653  	         startposs[1] = block + maxnblockvars;
654  	         endposs[1] = nconsvars;
655  	
656  	         /* loop over all variables before (p==0) and after (p==1) the most frequent block label */
657  	         for( p = 0; p < 2; ++p )
658  	         {
659  	            /* relabel */
660  	            for( v = startposs[p]; v < endposs[p]; ++v )
661  	               varslabels[v] = SCIP_DECOMP_LINKVAR;
662  	
663  	            /* set labels in the decomposition */
664  	            SCIP_CALL( SCIPdecompSetVarsLabels(decomp, &vars[startposs[p]], &varslabels[startposs[p]], endposs[p] - startposs[p]) );
665  	         }
666  	
667  	         SCIP_CALL( SCIPdecompSetConsLabels(decomp, &conss[c], &varslabels[block], 1) );
668  	      }
669  	   }
670  	
671  	   if( nskipconss != NULL )
672  	      *nskipconss = nskipconsslocal;
673  	
674  	   SCIPfreeBufferArray(scip, &vars);
675  	   SCIPfreeBufferArray(scip, &varslabels);
676  	
677  	   return SCIP_OKAY;
678  	}
679  	
680  	/** return position of a label in decomp array */
681  	static
682  	int findLabelIdx(
683  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
684  	   int                   label               /**< the label */
685  	   )
686  	{
687  	   int pos;
688  	
689  	   (void)SCIPsortedvecFindInt(decomp->labels, label, decomp->nblocks + 1, &pos);
690  	
691  	   return pos;
692  	}
693  	
694  	/** compute decomposition modularity (comparison of within block edges against a random decomposition) */
695  	static
696  	SCIP_RETCODE computeModularity(
697  	   SCIP*                 scip,               /**< SCIP data structure */
698  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
699  	   SCIP_Real*            modularity          /**< pointer to store modularity value */
700  	   )
701  	{
702  	   SCIP_CONS** conss;
703  	   SCIP_VAR** varbuf;
704  	   int* varslabels;
705  	   int* conslabels;
706  	   int* totaldegrees; /* the total degree for every block */
707  	   int* withinedges; /* the number of edges within each block */
708  	   int nnonzeroes = 0;
709  	   int varbufsize;
710  	   int nconss;
711  	   int c;
712  	   int b;
713  	
714  	   /* allocate buffer arrays to hold constraint and variable labels, and store within-edges and total community degrees */
715  	   getDecompVarsConssData(scip, decomp, NULL, &conss, NULL, &nconss);
716  	   varbufsize = getVarbufSize(scip);
717  	
718  	   SCIP_CALL( SCIPallocBufferArray(scip, &conslabels, nconss) );
719  	   SCIP_CALL( SCIPallocBufferArray(scip, &varslabels, varbufsize) );
720  	   SCIP_CALL( SCIPallocBufferArray(scip, &varbuf, varbufsize) );
721  	
722  	   SCIP_CALL( SCIPallocClearBufferArray(scip, &totaldegrees, decomp->nblocks + 1) );
723  	   SCIP_CALL( SCIPallocClearBufferArray(scip, &withinedges, decomp->nblocks + 1) );
724  	
725  	   SCIPdecompGetConsLabels(decomp, conss, conslabels, nconss);
726  	
727  	   /*
728  	    * loop over all nonzeros, consider the labels of the incident nodes (cons and variable)
729  	    * and increase the corresponding counters
730  	    */
731  	   for( c = 0; c < nconss; ++c )
732  	   {
733  	      int nconsvars;
734  	      int conslabel = conslabels[c];
735  	      int blockpos;
736  	      int varblockstart;
737  	      int requiredsize;
738  	      SCIP_Bool success;
739  	
740  	      /* linking constraints do not contribute to the modularity */
741  	      if( conslabel == SCIP_DECOMP_LINKCONS )
742  	         continue;
743  	
744  	      /* find the position of the constraint label. Constraints of the border always belong to the first block at index 0 */
745  	      blockpos = findLabelIdx(decomp, conslabel);
746  	
747  	      SCIP_CALL( decompGetConsVarsAndLabels(scip, decomp, conss[c], varbuf, varslabels,
748  	               varbufsize, &nconsvars, &requiredsize, &success) );
749  	      SCIP_CALL( ensureCondition(success) );
750  	
751  	      SCIPsortInt(varslabels, nconsvars);
752  	
753  	      /* count occurrences of labels (blocks) in the sorted labels array */
754  	      varblockstart = 0;
755  	      while( varblockstart < nconsvars )
756  	      {
757  	         int varblockpos;
758  	         int nblockvars = countLabelFromPos(varslabels, varblockstart, nconsvars);
759  	
760  	         varblockpos = findLabelIdx(decomp, varslabels[varblockstart]);
761  	
762  	         /* don't consider linking variables for modularity statistics */
763  	         if( varslabels[varblockstart] != SCIP_DECOMP_LINKVAR )
764  	         {
765  	            /* increase the number of within edges for variable and constraints from the same block */
766  	            if( varblockpos == blockpos )
767  	               withinedges[varblockpos] += nblockvars;
768  	
769  	            /* increase the total degrees and nonzero (edge) counts; it is intended that the total degrees sum up
770  	             * to twice the number of edges
771  	             */
772  	            totaldegrees[blockpos] += nblockvars;
773  	            totaldegrees[varblockpos] += nblockvars;
774  	            nnonzeroes += nblockvars;
775  	         }
776  	
777  	         varblockstart += nblockvars;
778  	      }
779  	   }
780  	
781  	/* ensure that total degrees sum up to twice the number of edges */
782  	#ifndef NDEBUG
783  	   {
784  	      int totaldegreesum = 0;
785  	      for( b = 1; b < decomp->nblocks + 1; ++b )
786  	         totaldegreesum += totaldegrees[b];
787  	
788  	      assert(totaldegreesum == 2 * nnonzeroes);
789  	   }
790  	#endif
791  	
792  	   /* compute modularity */
793  	   *modularity = 0.0;
794  	   nnonzeroes = MAX(nnonzeroes, 1);
795  	   for( b = 1; b < decomp->nblocks + 1; ++b )
796  	   {
797  	      SCIP_Real expectedval;
798  	      expectedval = totaldegrees[b] / (2.0 * nnonzeroes);
799  	      expectedval = SQR(expectedval);
800  	      *modularity += (withinedges[b] / (SCIP_Real)nnonzeroes) - expectedval;
801  	   }
802  	
803  	   SCIPfreeBufferArray(scip, &withinedges);
804  	   SCIPfreeBufferArray(scip, &totaldegrees);
805  	   SCIPfreeBufferArray(scip, &varbuf);
806  	   SCIPfreeBufferArray(scip, &varslabels);
807  	   SCIPfreeBufferArray(scip, &conslabels);
808  	
809  	   return SCIP_OKAY;
810  	}
811  	
812  	/** compute the area score */
813  	static
814  	void computeAreaScore(
815  	   SCIP*                 scip,               /**< SCIP data structure */
816  	   SCIP_DECOMP*          decomp              /**< decomposition data structure */
817  	   )
818  	{
819  	   SCIP_Real areascore = 1.0;
820  	   int nvars;
821  	   int nconss;
822  	
823  	   getDecompVarsConssData(scip, decomp, NULL, NULL, &nvars, &nconss);
824  	
825  	   if( nvars > 0 && nconss > 0 )
826  	   {
827  	      int nlinkconss = decomp->consssize[0];
828  	      int nlinkvars = decomp->varssize[0];
829  	      SCIP_Real factor = 1.0 / ((SCIP_Real)nvars * nconss);
830  	
831  	      int i;
832  	
833  	      /* compute diagonal contributions to the area score */
834  	      for( i = 1; i < decomp->nblocks + 1; ++i )
835  	      {
836  	         areascore -= (factor * decomp->consssize[i]) * decomp->varssize[i];
837  	      }
838  	
839  	      areascore -= ((SCIP_Real)nlinkconss * nvars + (SCIP_Real)nconss * nlinkvars - (SCIP_Real)nlinkconss * nlinkvars) * factor;
840  	   }
841  	
842  	   decomp->areascore = areascore;
843  	}
844  	
845  	/** build the block decomposition graph */
846  	static
847  	SCIP_RETCODE buildBlockGraph(
848  	   SCIP*                 scip,               /**< SCIP data structure */
849  	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
850  	   int                   maxgraphedge        /**< maximum number of edges in block graph computation (-1: no limit, 0: disable block graph computation) */
851  	   )
852  	{
853  	   SCIP_VAR** vars;
854  	   SCIP_CONS** conss;
855  	   SCIP_VAR** consvars;
856  	   SCIP_DIGRAPH* blocklinkingvargraph;
857  	   SCIP_DIGRAPH* blockgraph = NULL;
858  	   int* varlabels;
859  	   int* conslabels;
860  	   SCIP_CONS** consscopy; /* working copy of the constraints */
861  	   int* linkvaridx;
862  	   int* succnodesblk;
863  	   int* succnodesvar;
864  	   SCIP_Bool success;
865  	   int varbufsize;
866  	   int nvars;
867  	   int nconss;
868  	   int nblocks;
869  	   int nlinkingvars = 0;
870  	   int nconsvars;
871  	   int conspos;
872  	   int tempmin;
873  	   int tempmax;
874  	   int nblockgraphedges;
875  	   int blocknodeidx;
876  	   int i;
877  	   int j;
878  	   int v;
879  	   int n;
880  	
881  	   if( maxgraphedge == -1 )
882  	      maxgraphedge = INT_MAX;
883  	
884  	   assert(scip != NULL);
885  	   assert(decomp != NULL);
886  	   assert(decomp->statscomplete == FALSE);
887  	
888  	   /* capture the trivial case that no linking variables are present */
889  	   if( decomp->varssize[0] == 0 || decomp->nblocks == 0 )
890  	   {
891  	      decomp->mindegree = 0;
892  	      decomp->maxdegree = 0;
893  	      decomp->nedges = 0;
894  	      decomp->ncomponents = SCIPdecompGetNBlocks(decomp);
895  	      decomp->narticulations = 0;
896  	      decomp->statscomplete = TRUE;
897  	
898  	      return SCIP_OKAY;
899  	   }
900  	
901  	   getDecompVarsConssData(scip, decomp, &vars, &conss, &nvars, &nconss);
902  	   varbufsize = getVarbufSize(scip);
903  	   nblocks = SCIPdecompGetNBlocks(decomp);
904  	
905  	   SCIP_CALL( SCIPallocBufferArray(scip, &conslabels, nconss) );
906  	   SCIP_CALL( SCIPallocBufferArray(scip, &varlabels, varbufsize) );
907  	   SCIP_CALL( SCIPallocBufferArray(scip, &linkvaridx, varbufsize) );
908  	   SCIP_CALL( SCIPallocBufferArray(scip, &consvars, varbufsize) );
909  	
910  	   /* store variable and constraint labels in buffer arrays */
911  	   SCIPdecompGetConsLabels(decomp, conss, conslabels, nconss);
912  	   SCIPdecompGetVarsLabels(decomp, vars, varlabels, nvars);
913  	
914  	   /* create a mapping of all linking variables to 0,..,nlinkingvars -1 and store it in array linkvaridx */
915  	   for( v = 0; v < nvars; ++v )
916  	   {
917  	      if( varlabels[v] == SCIP_DECOMP_LINKVAR )
918  	      {
919  	         linkvaridx[v] = nlinkingvars;
920  	         assert(SCIPvarGetProbindex(vars[v]) == v);
921  	         ++nlinkingvars;
922  	      }
923  	      else
924  	         linkvaridx[v] = -1;
925  	   }
926  	
927  	   /* create a bipartite graph composed of block and linking var nodes */
928  	   SCIP_CALL( SCIPcreateDigraph(scip, &blocklinkingvargraph, nblocks + nlinkingvars) );/* nblocks does not include the linking constraints block */
929  	
930  	   /* initialize position to start after the linking constraints, which we skip anyway */
931  	   SCIP_CALL( SCIPduplicateBufferArray(scip, &consscopy, conss, nconss) );
932  	   SCIPsortIntPtr(conslabels, (void**)consscopy, nconss);
933  	   if( conslabels[0] == SCIP_DECOMP_LINKCONS )
934  	      conspos = countLabelFromPos(conslabels, 0, nconss);
935  	   else
936  	      /* no linking constraints present */
937  	      conspos = 0;
938  	
939  	   blocknodeidx = -1;
940  	   /* loop over each block */
941  	   while( conspos < nconss )
942  	   {
943  	      SCIP_Bool* adjacent;
944  	      int* adjacentidxs;
945  	      int nblockconss = countLabelFromPos(conslabels, conspos, nconss);
946  	      int nblocklinkingvars = 0;
947  	      int c;
948  	
949  	      ++blocknodeidx;
950  	      /* allocate buffer storage to store all linking variable indices adjacent to this block */
951  	      SCIP_CALL( SCIPallocCleanBufferArray(scip, &adjacent, nlinkingvars) );
952  	      SCIP_CALL( SCIPallocBufferArray(scip, &adjacentidxs, nlinkingvars) );
953  	
954  	      /* loop over the constraints of this block; stop if the block vertex has maximum degree */
955  	      for( c = conspos; c < conspos + nblockconss && nblocklinkingvars < nlinkingvars; ++c )
956  	      {
957  	         int requiredsize;
958  	         SCIP_CONS* cons = consscopy[c];
959  	         assert(conslabels[c] != SCIP_DECOMP_LINKCONS);
960  	
961  	         SCIP_CALL( decompGetConsVarsAndLabels(scip, decomp, cons, consvars, varlabels,
962  	               varbufsize, &nconsvars, &requiredsize, &success) );
963  	         SCIP_CALL( ensureCondition(success) );
964  	
965  	         /* search for linking variables that are not connected so far; stop as soon as block vertex has max degree */
966  	         for( j = 0; j < nconsvars && nblocklinkingvars < nlinkingvars; ++j )
967  	         {
968  	            int linkingvarnodeidx;
969  	            /* consider only linking variables */
970  	            if( varlabels[j] != SCIP_DECOMP_LINKVAR )
971  	               continue;
972  	
973  	            linkingvarnodeidx = linkvaridx[SCIPvarGetProbindex(consvars[j])];
974  	            assert(linkingvarnodeidx >= 0);
975  	
976  	            if( !adjacent[linkingvarnodeidx] )
977  	            {
978  	               adjacent[linkingvarnodeidx] = TRUE;
979  	               adjacentidxs[nblocklinkingvars++] = linkingvarnodeidx;
980  	            }
981  	         }
982  	      }
983  	
984  	      /* connect block and linking variables in the digraph */
985  	      assert(blocknodeidx == findLabelIdx(decomp, conslabels[conspos]) - 1);
986  	      for( i = 0; i < nblocklinkingvars; ++i )
987  	      {
988  	         SCIP_CALL( SCIPdigraphAddArc(blocklinkingvargraph, blocknodeidx, nblocks + adjacentidxs[i], NULL) );
989  	         SCIP_CALL( SCIPdigraphAddArc(blocklinkingvargraph, nblocks + adjacentidxs[i], blocknodeidx, NULL) );
990  	      }
991  	
992  	      /* clean up the adjacent array before freeing */
993  	      for( i = 0; i < nblocklinkingvars; ++i )
994  	         adjacent[adjacentidxs[i]] = FALSE;
995  	
996  	      /* check that adjacent has been entirely cleaned up */
997  	#ifndef NDEBUG
998  	      for( i = 0; i < nlinkingvars; ++i )
999  	         assert(adjacent[i] == FALSE);
1000 	#endif
1001 	
1002 	      SCIPfreeBufferArray(scip, &adjacentidxs);
1003 	      SCIPfreeCleanBufferArray(scip, &adjacent);
1004 	
1005 	      conspos += nblockconss;
1006 	   }
1007 	
1008 	   SCIPfreeBufferArray(scip, &consscopy);
1009 	
1010 	   assert(SCIPdigraphGetNNodes(blocklinkingvargraph) > 0);
1011 	   /* check first if any of the linking variables is connected with all blocks -> block graph is complete and connected */
1012 	   for( n = nblocks; n < SCIPdigraphGetNNodes(blocklinkingvargraph); ++n )
1013 	   {
1014 	      int nsuccvar;
1015 	      nsuccvar = (int) SCIPdigraphGetNSuccessors(blocklinkingvargraph, n);
1016 	
1017 	      if( nsuccvar == nblocks )
1018 	      {
1019 	         decomp->nedges = nblocks * (nblocks - 1) / 2;
1020 	         decomp->mindegree = decomp->maxdegree = nblocks - 1;
1021 	         decomp->narticulations = 0;
1022 	         decomp->ncomponents = 1;
1023 	         decomp->statscomplete = TRUE;
1024 	
1025 	         goto TERMINATE;
1026 	      }
1027 	   }
1028 	
1029 	   /* from the information of the above bipartite graph, build the block-decomposition graph: nodes -> blocks and double-direction arcs -> linking variables */
1030 	   SCIP_CALL( SCIPcreateDigraph(scip, &blockgraph, nblocks) );
1031 	
1032 	   /* we count the number of block graph edges manually, because SCIPdigraphGetNArcs() iterates over all nodes */
1033 	   nblockgraphedges = 0;
1034 	   for( n = 0; n < nblocks - 1 && nblockgraphedges < maxgraphedge; ++n )
1035 	   {
1036 	      SCIP_Bool* adjacent; /* an array to mark the adjacent blocks to the current block */
1037 	      int* adjacentidxs;
1038 	      int nsuccblk;
1039 	      int nadjacentblks = 0;
1040 	      SCIP_CALL( SCIPallocCleanBufferArray(scip, &adjacent, nblocks) );
1041 	      SCIP_CALL( SCIPallocBufferArray(scip, &adjacentidxs, nblocks) );
1042 	
1043 	      assert(n < SCIPdigraphGetNNodes(blocklinkingvargraph));
1044 	
1045 	      /* loop over the connected linking variables to the current block and their connected blocks to update the adjacency array */
1046 	      nsuccblk = SCIPdigraphGetNSuccessors(blocklinkingvargraph, n);
1047 	      succnodesblk = SCIPdigraphGetSuccessors(blocklinkingvargraph, n);
1048 	      for( i = 0; i < nsuccblk && nadjacentblks < nblocks - (n + 1); ++i )
1049 	      {
1050 	         int startpos;
1051 	         int nsuccvar;
1052 	
1053 	         assert(succnodesblk[i] < SCIPdigraphGetNNodes(blocklinkingvargraph));
1054 	
1055 	         nsuccvar = SCIPdigraphGetNSuccessors(blocklinkingvargraph, succnodesblk[i]);
1056 	         succnodesvar = SCIPdigraphGetSuccessors(blocklinkingvargraph, succnodesblk[i]);
1057 	
1058 	         /* previously visited blocks can be skipped in this step */
1059 	         if( ! SCIPsortedvecFindInt(succnodesvar, n, nsuccvar, &startpos) )
1060 	            SCIPABORT();
1061 	         for( j = startpos + 1; j < nsuccvar; ++j )
1062 	         {
1063 	            assert( succnodesvar[j] > n );
1064 	            if( !adjacent[succnodesvar[j]] )
1065 	            {
1066 	               adjacent[succnodesvar[j]] = TRUE;
1067 	               adjacentidxs[nadjacentblks] = succnodesvar[j];
1068 	               ++nadjacentblks;
1069 	            }
1070 	         }
1071 	      }
1072 	
1073 	      /* double-direction arcs are added in this step between the current block and its adjacent block nodes */
1074 	      for( i = 0; i < nadjacentblks && nblockgraphedges < maxgraphedge; ++i )
1075 	      {
1076 	          SCIP_CALL( SCIPdigraphAddArc(blockgraph, n, adjacentidxs[i], NULL) );
1077 	          SCIP_CALL( SCIPdigraphAddArc(blockgraph, adjacentidxs[i], n, NULL) );
1078 	
1079 	          ++nblockgraphedges;
1080 	      }
1081 	
1082 	      /* clean up the adjacent array and free it */
1083 	      for( i = 0; i < nadjacentblks; ++i )
1084 	         adjacent[adjacentidxs[i]] = FALSE;
1085 	
1086 	      SCIPfreeBufferArray(scip, &adjacentidxs);
1087 	      SCIPfreeCleanBufferArray(scip, &adjacent);
1088 	   }
1089 	
1090 	   assert(SCIPdigraphGetNNodes(blockgraph) > 0);
1091 	
1092 	   /* Get the number of edges in the block-decomposition graph.*/
1093 	   decomp->nedges = nblockgraphedges;
1094 	   decomp->statscomplete = nblockgraphedges < maxgraphedge;
1095 	
1096 	   /* Get the minimum and maximum degree of the block-decomposition graph */
1097 	   tempmin = (int) SCIPdigraphGetNSuccessors(blockgraph, 0);
1098 	   tempmax = (int) SCIPdigraphGetNSuccessors(blockgraph, 0);
1099 	   for( n = 1; n < SCIPdigraphGetNNodes(blockgraph); ++n )
1100 	   {
1101 	      int nsuccblk = SCIPdigraphGetNSuccessors(blockgraph, n);
1102 	
1103 	      if( nsuccblk < tempmin )
1104 	         tempmin = nsuccblk;
1105 	      else if( nsuccblk > tempmax )
1106 	         tempmax = nsuccblk;
1107 	   }
1108 	
1109 	   decomp->mindegree = tempmin;
1110 	   decomp->maxdegree = tempmax;
1111 	
1112 	   /* Get the number of connected components in the block-decomposition graph.*/
1113 	   SCIP_CALL( SCIPdigraphComputeUndirectedComponents(blockgraph, -1, NULL, NULL) );
1114 	   decomp->ncomponents = SCIPdigraphGetNComponents(blockgraph);
1115 	
1116 	   /* Get the number of articulation points in the block-decomposition graph using DFS.*/
1117 	   SCIP_CALL( SCIPdigraphGetArticulationPoints(blockgraph, NULL, &decomp->narticulations) );
1118 	
1119 	TERMINATE:
1120 	   SCIPfreeBufferArray(scip, &consvars);
1121 	   SCIPfreeBufferArray(scip, &linkvaridx);
1122 	   SCIPfreeBufferArray(scip, &varlabels);
1123 	   SCIPfreeBufferArray(scip, &conslabels);
1124 	
1125 	   /* blockgraph has probably not been allocated */
1126 	   if( blockgraph != NULL)
1127 	      SCIPdigraphFree(&blockgraph);
1128 	
1129 	   SCIPdigraphFree(&blocklinkingvargraph);
1130 	
1131 	   return SCIP_OKAY;
1132 	}
1133 	
1134 	/** computes decomposition statistics and store them in the decomposition object */
1135 	SCIP_RETCODE SCIPcomputeDecompStats(
1136 	   SCIP*                 scip,               /**< SCIP data structure */
1137 	   SCIP_DECOMP*          decomp,             /**< decomposition data structure */
1138 	   SCIP_Bool             uselimits           /**< respect user limits on potentially expensive graph statistics? */
1139 	   )
1140 	{
1141 	   SCIP_VAR** vars;
1142 	   SCIP_CONS** conss;
1143 	   SCIP_CONS** conssarray;
1144 	   SCIP_VAR** varsarray;
1145 	   int* varslabels;
1146 	   int* conslabels;
1147 	   int nvars;
1148 	   int nconss;
1149 	   int varblockstart;
1150 	   int consblockstart;
1151 	   int currlabelidx;
1152 	   int varidx;
1153 	   int considx;
1154 	   int i;
1155 	   int maxgraphedge;
1156 	   SCIP_Bool disablemeasures;
1157 	
1158 	   assert(scip != NULL);
1159 	   assert(decomp != NULL);
1160 	
1161 	   getDecompVarsConssData(scip, decomp, &vars, &conss, &nvars, &nconss);
1162 	
1163 	   /* return if problem is empty
1164 	    *
1165 	    * TODO ensure that statistics reflect this correctly
1166 	    */
1167 	   if( nvars == 0 || nconss == 0 )
1168 	   {
1169 	      decomp->nblocks = 0;
1170 	      decomp->varssize[0] = nvars;
1171 	      decomp->consssize[0] = nconss;
1172 	      decomp->labels[0] = SCIP_DECOMP_LINKVAR;
1173 	      return SCIP_OKAY;
1174 	   }
1175 	
1176 	   decomp->statscomplete = FALSE;
1177 	
1178 	   /* store variable and constraint labels in buffer arrays */
1179 	   SCIP_CALL( SCIPduplicateBufferArray(scip, &conssarray, conss, nconss) );
1180 	   SCIP_CALL( SCIPallocBufferArray(scip, &conslabels, nconss) );
1181 	   SCIP_CALL( SCIPduplicateBufferArray(scip, &varsarray, vars, nvars) );
1182 	   SCIP_CALL( SCIPallocBufferArray(scip, &varslabels, nvars) );
1183 	
1184 	   SCIPdecompGetVarsLabels(decomp, varsarray, varslabels, nvars);
1185 	   SCIPdecompGetConsLabels(decomp, conssarray, conslabels, nconss);
1186 	
1187 	   /* sort both buffer arrays for quick counting */
1188 	   SCIPsortIntPtr(varslabels, (void**)varsarray, nvars);
1189 	   SCIPsortIntPtr(conslabels, (void**)conssarray, nconss);
1190 	
1191 	   /* the first label is always LINKVAR, even if Benders' variable labels are used. We can ignore the variables
1192 	    * labelled as LINKCONS since this label is only required when computing the variable labels for Benders'
1193 	    * decomposition.
1194 	    */
1195 	   decomp->labels[0] = SCIP_DECOMP_LINKVAR;
1196 	
1197 	   /* treating the linking variables first */
1198 	   if( varslabels[0] == SCIP_DECOMP_LINKVAR )
1199 	      decomp->varssize[0] = countLabelFromPos(varslabels, 0, nvars);
1200 	   else
1201 	      decomp->varssize[0] = 0;
1202 	
1203 	   /* count border constraints and store their number */
1204 	   if( conslabels[0] == SCIP_DECOMP_LINKCONS )
1205 	      decomp->consssize[0] = countLabelFromPos(conslabels, 0, nconss);
1206 	   else
1207 	      decomp->consssize[0] = 0;
1208 	
1209 	   /* merge labels (except for border at position 0) since neither variable nor constraint labels by themselves need to be complete */
1210 	   currlabelidx = 1;
1211 	   varidx = decomp->varssize[0];
1212 	   considx = decomp->consssize[0];
1213 	
1214 	   while( varidx < nvars || considx < nconss )
1215 	   {
1216 	      int varlabel;
1217 	      int conslabel;
1218 	
1219 	      varlabel = varidx < nvars ? varslabels[varidx] : INT_MAX;
1220 	      conslabel = considx < nconss ? conslabels[considx] : INT_MAX;
1221 	
1222 	      assert(currlabelidx < decomp->memsize);
1223 	      /* store the smaller of the two current labels */
1224 	      decomp->labels[currlabelidx] = MIN(varlabel, conslabel);
1225 	
1226 	      /* a strictly larger variable label means that there are no variables for the current label */
1227 	      if( varlabel <= conslabel )
1228 	         decomp->varssize[currlabelidx] = countLabelFromPos(varslabels, varidx, nvars);
1229 	      else
1230 	         decomp->varssize[currlabelidx] = 0;
1231 	
1232 	      /* the same for constraint labels */
1233 	      if( conslabel <= varlabel )
1234 	         decomp->consssize[currlabelidx] = countLabelFromPos(conslabels, considx, nconss);
1235 	      else
1236 	         decomp->consssize[currlabelidx] = 0;
1237 	
1238 	      /* increase indices appropriately */
1239 	      varidx += decomp->varssize[currlabelidx];
1240 	      considx += decomp->consssize[currlabelidx];
1241 	
1242 	      currlabelidx++;
1243 	   }
1244 	
1245 	   SCIPdebugMsg(scip, "Counted %d different labels (should be %d)\n", currlabelidx, decomp->nblocks + 1);
1246 	
1247 	   /* strip the remaining, unused blocks */
1248 	   if( currlabelidx < decomp->nblocks + 1 )
1249 	      decomp->nblocks = currlabelidx - 1;
1250 	
1251 	   /* delete empty blocks from statistics, relabel the corresponding constraints/variables as linking */
1252 	   varblockstart = decomp->varssize[0];
1253 	   consblockstart = decomp->consssize[0];
1254 	
1255 	   for( i = 1; i < decomp->nblocks + 1; ++i )
1256 	   {
1257 	      assert(MAX(decomp->varssize[i], decomp->consssize[i]) > 0);
1258 	      /* relabel constraint blocks as linking, if there are no corresponding variables */
1259 	      if( decomp->varssize[i] == 0 )
1260 	      {
1261 	         int nblockconss = decomp->consssize[i];
1262 	         int c;
1263 	         /* relabel these constraints as linking */
1264 	         for( c = consblockstart; c < consblockstart + nblockconss; ++c )
1265 	            conslabels[c] = SCIP_DECOMP_LINKCONS;
1266 	
1267 	         SCIP_CALL( SCIPdecompSetConsLabels(decomp, &conssarray[consblockstart], &conslabels[consblockstart], nblockconss) );
1268 	
1269 	         /* increase number of linking constraints */
1270 	         decomp->consssize[0] += nblockconss;
1271 	      }
1272 	
1273 	      /* same for constraints */
1274 	      if( decomp->consssize[i] == 0 )
1275 	      {
1276 	         int nblockvars = decomp->varssize[i];
1277 	         int v;
1278 	
1279 	         /* relabel the variables as linking variables */
1280 	         for( v = varblockstart; v < varblockstart + nblockvars; ++v )
1281 	            varslabels[v] = SCIP_DECOMP_LINKVAR;
1282 	
1283 	         SCIP_CALL( SCIPdecompSetVarsLabels(decomp, &varsarray[varblockstart], &varslabels[varblockstart], nblockvars) );
1284 	
1285 	         /* increase number of linking variables */
1286 	         decomp->varssize[0] += nblockvars;
1287 	      }
1288 	
1289 	      varblockstart += decomp->varssize[i];
1290 	      consblockstart += decomp->consssize[i];
1291 	   }
1292 	
1293 	   currlabelidx = 1;
1294 	
1295 	   /* delete empty blocks; they are no longer present */
1296 	   for( i = 1; i < decomp->nblocks + 1; ++i )
1297 	   {
1298 	      /* keep only nonempty blocks */
1299 	      if( decomp->varssize[i] > 0 && decomp->consssize[i] > 0 )
1300 	      {
1301 	         decomp->labels[currlabelidx] = decomp->labels[i];
1302 	         decomp->varssize[currlabelidx] = decomp->varssize[i];
1303 	         decomp->consssize[currlabelidx] = decomp->consssize[i];
1304 	
1305 	         currlabelidx++;
1306 	      }
1307 	   }
1308 	
1309 	   decomp->nblocks = currlabelidx - 1;
1310 	
1311 	   decomp->idxsmallestblock = decomp->idxlargestblock = -1;
1312 	   /* now that indices are fixed, store indices with largest and smallest number of constraints */
1313 	   for( i = 1; i < decomp->nblocks + 1; ++i )
1314 	   {
1315 	      if( decomp->idxsmallestblock == -1 )
1316 	         decomp->idxsmallestblock = decomp->idxlargestblock = i;
1317 	      else if( decomp->consssize[decomp->idxsmallestblock] > decomp->consssize[i] )
1318 	         decomp->idxsmallestblock = i;
1319 	      else if( decomp->consssize[decomp->idxlargestblock] < decomp->consssize[i] )
1320 	         decomp->idxlargestblock = i;
1321 	   }
1322 	
1323 	   /* compute more involved statistics such as the area score, the modularity, and the block graph statistics */
1324 	   SCIP_CALL( SCIPgetBoolParam(scip, "decomposition/disablemeasures", &disablemeasures) );
1325 	   if( !disablemeasures )
1326 	   {
1327 	      SCIP_CALL( computeModularity(scip, decomp, &decomp->modularity) );
1328 	      computeAreaScore(scip, decomp);
1329 	   }
1330 	
1331 	   if( uselimits )
1332 	   {
1333 	      SCIP_CALL( SCIPgetIntParam(scip, "decomposition/maxgraphedge", &maxgraphedge) );
1334 	   }
1335 	   else
1336 	      maxgraphedge = -1;
1337 	
1338 	   /* do not start computation of the block graph if maxgraphedge is set to 0 */
1339 	   if( maxgraphedge != 0 )
1340 	   {
1341 	      SCIP_CALL( buildBlockGraph(scip, decomp, maxgraphedge) );
1342 	   }
1343 	
1344 	   SCIPfreeBufferArray(scip, &varslabels);
1345 	   SCIPfreeBufferArray(scip, &varsarray);
1346 	   SCIPfreeBufferArray(scip, &conslabels);
1347 	   SCIPfreeBufferArray(scip, &conssarray);
1348 	
1349 	   return SCIP_OKAY;
1350 	}
1351