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