1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file presol_sparsify.c 26 * @ingroup DEFPLUGINS_PRESOL 27 * @brief cancel non-zeros of the constraint matrix 28 * @author Dieter Weninger 29 * @author Leona Gottwald 30 * @author Ambros Gleixner 31 * 32 * This presolver attempts to cancel non-zero entries of the constraint 33 * matrix by adding scaled equalities to other constraints. 34 */ 35 36 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 37 38 #include "blockmemshell/memory.h" 39 #include "scip/cons_linear.h" 40 #include "scip/presol_sparsify.h" 41 #include "scip/pub_cons.h" 42 #include "scip/pub_matrix.h" 43 #include "scip/pub_message.h" 44 #include "scip/pub_misc.h" 45 #include "scip/pub_misc_sort.h" 46 #include "scip/pub_presol.h" 47 #include "scip/pub_var.h" 48 #include "scip/scip_cons.h" 49 #include "scip/scip_general.h" 50 #include "scip/scip_mem.h" 51 #include "scip/scip_message.h" 52 #include "scip/scip_nlp.h" 53 #include "scip/scip_numerics.h" 54 #include "scip/scip_param.h" 55 #include "scip/scip_presol.h" 56 #include "scip/scip_pricer.h" 57 #include "scip/scip_prob.h" 58 #include "scip/scip_probing.h" 59 #include "scip/scip_timing.h" 60 #include <string.h> 61 62 #define PRESOL_NAME "sparsify" 63 #define PRESOL_DESC "eliminate non-zero coefficients" 64 65 #define PRESOL_PRIORITY -24000 /**< priority of the presolver (>= 0: before, < 0: after constraint handlers) */ 66 #define PRESOL_MAXROUNDS -1 /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */ 67 #define PRESOL_TIMING SCIP_PRESOLTIMING_EXHAUSTIVE /* timing of the presolver (fast, medium, or exhaustive) */ 68 69 #define DEFAULT_ENABLECOPY TRUE /**< should sparsify presolver be copied to sub-SCIPs? */ 70 #define DEFAULT_CANCELLINEAR TRUE /**< should we cancel nonzeros in constraints of the linear constraint handler? */ 71 #define DEFAULT_PRESERVEINTCOEFS TRUE /**< should we forbid cancellations that destroy integer coefficients? */ 72 #define DEFAULT_MAX_CONT_FILLIN 0 /**< default value for the maximal fillin for continuous variables */ 73 #define DEFAULT_MAX_BIN_FILLIN 0 /**< default value for the maximal fillin for binary variables */ 74 #define DEFAULT_MAX_INT_FILLIN 0 /**< default value for the maximal fillin for integer variables (including binary) */ 75 #define DEFAULT_MAXNONZEROS -1 /**< maximal support of one equality to be used for cancelling (-1: no limit) */ 76 #define DEFAULT_MAXCONSIDEREDNONZEROS 70 /**< maximal number of considered non-zeros within one row (-1: no limit) */ 77 #define DEFAULT_ROWSORT 'd' /**< order in which to process inequalities ('n'o sorting, 'i'ncreasing nonzeros, 'd'ecreasing nonzeros) */ 78 #define DEFAULT_MAXRETRIEVEFAC 100.0 /**< limit on the number of useless vs. useful hashtable retrieves as a multiple of the number of constraints */ 79 #define DEFAULT_WAITINGFAC 2.0 /**< number of calls to wait until next execution as a multiple of the number of useless calls */ 80 81 #define MAXSCALE 1000.0 /**< maximal allowed scale for cancelling non-zeros */ 82 83 84 /* 85 * Data structures 86 */ 87 88 /** presolver data */ 89 struct SCIP_PresolData 90 { 91 int ncancels; /**< total number of canceled nonzeros (net value, i.e., removed minus added nonzeros) */ 92 int nfillin; /**< total number of added nonzeros */ 93 int nfailures; /**< number of calls to presolver without success */ 94 int nwaitingcalls; /**< number of presolver calls until next real execution */ 95 int maxcontfillin; /**< maximal fillin for continuous variables */ 96 int maxintfillin; /**< maximal fillin for integer variables*/ 97 int maxbinfillin; /**< maximal fillin for binary variables */ 98 int maxnonzeros; /**< maximal support of one equality to be used for cancelling (-1: no limit) */ 99 int maxconsiderednonzeros;/**< maximal number of considered non-zeros within one row (-1: no limit) */ 100 SCIP_Real maxretrievefac; /**< limit on the number of useless vs. useful hashtable retrieves as a multiple of the number of constraints */ 101 SCIP_Real waitingfac; /**< number of calls to wait until next execution as a multiple of the number of useless calls */ 102 char rowsort; /**< order in which to process inequalities ('n'o sorting, 'i'ncreasing nonzeros, 'd'ecreasing nonzeros) */ 103 SCIP_Bool enablecopy; /**< should sparsify presolver be copied to sub-SCIPs? */ 104 SCIP_Bool cancellinear; /**< should we cancel nonzeros in constraints of the linear constraint handler? */ 105 SCIP_Bool preserveintcoefs; /**< should we forbid cancellations that destroy integer coefficients? */ 106 }; 107 108 /** structure representing a pair of variables in a row; used for lookup in a hashtable */ 109 struct RowVarPair 110 { 111 int rowindex; 112 int varindex1; 113 int varindex2; 114 SCIP_Real varcoef1; 115 SCIP_Real varcoef2; 116 }; 117 118 typedef struct RowVarPair ROWVARPAIR; 119 120 /* 121 * Local methods 122 */ 123 124 /** returns TRUE iff both keys are equal */ 125 static 126 SCIP_DECL_HASHKEYEQ(varPairsEqual) 127 { /*lint --e{715}*/ 128 SCIP* scip; 129 ROWVARPAIR* varpair1; 130 ROWVARPAIR* varpair2; 131 SCIP_Real ratio1; 132 SCIP_Real ratio2; 133 134 scip = (SCIP*) userptr; 135 varpair1 = (ROWVARPAIR*) key1; 136 varpair2 = (ROWVARPAIR*) key2; 137 138 if( varpair1->varindex1 != varpair2->varindex1 ) 139 return FALSE; 140 141 if( varpair1->varindex2 != varpair2->varindex2 ) 142 return FALSE; 143 144 ratio1 = varpair1->varcoef2 / varpair1->varcoef1; 145 ratio2 = varpair2->varcoef2 / varpair2->varcoef1; 146 147 return SCIPisEQ(scip, ratio1, ratio2); 148 } 149 150 /** returns the hash value of the key */ 151 static 152 SCIP_DECL_HASHKEYVAL(varPairHashval) 153 { /*lint --e{715}*/ 154 ROWVARPAIR* varpair; 155 156 varpair = (ROWVARPAIR*) key; 157 158 return SCIPhashThree(varpair->varindex1, varpair->varindex2, 159 SCIPrealHashCode(varpair->varcoef2 / varpair->varcoef1)); 160 } 161 162 /** try non-zero cancellation for given row */ 163 static 164 SCIP_RETCODE cancelRow( 165 SCIP* scip, /**< SCIP datastructure */ 166 SCIP_MATRIX* matrix, /**< the constraint matrix */ 167 SCIP_HASHTABLE* pairtable, /**< the hashtable containing ROWVARPAIR's of equations */ 168 int rowidx, /**< index of row to try non-zero cancellation for */ 169 int maxcontfillin, /**< maximal fill-in allowed for continuous variables */ 170 int maxintfillin, /**< maximal fill-in allowed for integral variables */ 171 int maxbinfillin, /**< maximal fill-in allowed for binary variables */ 172 int maxconsiderednonzeros, /**< maximal number of non-zeros to consider for cancellation */ 173 SCIP_Bool preserveintcoefs, /**< only perform non-zero cancellation if integrality of coefficients is preserved? */ 174 SCIP_Longint* nuseless, /**< pointer to update number of useless hashtable retrieves */ 175 int* nchgcoefs, /**< pointer to update number of changed coefficients */ 176 int* ncanceled, /**< pointer to update number of canceled nonzeros */ 177 int* nfillin /**< pointer to update the produced fill-in */ 178 ) 179 { 180 int* cancelrowinds; 181 SCIP_Real* cancelrowvals; 182 SCIP_Real cancellhs; 183 SCIP_Real cancelrhs; 184 SCIP_Real bestcancelrate; 185 int* tmpinds; 186 int* locks; 187 SCIP_Real* tmpvals; 188 int cancelrowlen; 189 int* rowidxptr; 190 SCIP_Real* rowvalptr; 191 int nchgcoef; 192 int nretrieves; 193 int bestnfillin; 194 SCIP_Real mincancelrate; 195 SCIP_Bool rowiseq; 196 SCIP_Bool swapped = FALSE; 197 SCIP_CONS* cancelcons; 198 199 rowiseq = SCIPisEQ(scip, SCIPmatrixGetRowLhs(matrix, rowidx), SCIPmatrixGetRowRhs(matrix, rowidx)); 200 201 cancelrowlen = SCIPmatrixGetRowNNonzs(matrix, rowidx); 202 rowidxptr = SCIPmatrixGetRowIdxPtr(matrix, rowidx); 203 rowvalptr = SCIPmatrixGetRowValPtr(matrix, rowidx); 204 205 cancelcons = SCIPmatrixGetCons(matrix, rowidx); 206 207 mincancelrate = 0.0; 208 209 /* for set packing and logicor constraints, only accept equalities where all modified coefficients are cancelled */ 210 if( SCIPconsGetHdlr(cancelcons) == SCIPfindConshdlr(scip, "setppc") || 211 SCIPconsGetHdlr(cancelcons) == SCIPfindConshdlr(scip, "logicor") ) 212 mincancelrate = 1.0; 213 214 SCIP_CALL( SCIPduplicateBufferArray(scip, &cancelrowinds, rowidxptr, cancelrowlen) ); 215 SCIP_CALL( SCIPduplicateBufferArray(scip, &cancelrowvals, rowvalptr, cancelrowlen) ); 216 SCIP_CALL( SCIPallocBufferArray(scip, &tmpinds, cancelrowlen) ); 217 SCIP_CALL( SCIPallocBufferArray(scip, &tmpvals, cancelrowlen) ); 218 SCIP_CALL( SCIPallocBufferArray(scip, &locks, cancelrowlen) ); 219 220 cancellhs = SCIPmatrixGetRowLhs(matrix, rowidx); 221 cancelrhs = SCIPmatrixGetRowRhs(matrix, rowidx); 222 223 nchgcoef = 0; 224 nretrieves = 0; 225 while( TRUE ) /*lint !e716 */ 226 { 227 SCIP_Real bestscale; 228 int bestcand; 229 int i; 230 int j; 231 ROWVARPAIR rowvarpair; 232 int maxlen; 233 234 bestscale = 1.0; 235 bestcand = -1; 236 bestnfillin = 0; 237 bestcancelrate = 0.0; 238 239 for( i = 0; i < cancelrowlen; ++i ) 240 { 241 tmpinds[i] = i; 242 locks[i] = SCIPmatrixGetColNDownlocks(matrix, cancelrowinds[i]) + SCIPmatrixGetColNUplocks(matrix, cancelrowinds[i]); 243 } 244 245 SCIPsortIntInt(locks, tmpinds, cancelrowlen); 246 247 maxlen = cancelrowlen; 248 if( maxconsiderednonzeros >= 0 ) 249 maxlen = MIN(cancelrowlen, maxconsiderednonzeros); 250 251 for( i = 0; i < maxlen; ++i ) 252 { 253 for( j = i + 1; j < maxlen; ++j ) 254 { 255 int a,b; 256 int ncancel; 257 int ncontfillin; 258 int nintfillin; 259 int nbinfillin; 260 int ntotfillin; 261 int eqrowlen; 262 ROWVARPAIR* eqrowvarpair; 263 SCIP_Real* eqrowvals; 264 int* eqrowinds; 265 SCIP_Real scale; 266 SCIP_Real cancelrate; 267 int i1,i2; 268 SCIP_Bool abortpair; 269 270 i1 = tmpinds[i]; 271 i2 = tmpinds[j]; 272 273 assert(cancelrowinds[i] < cancelrowinds[j]); 274 275 if( cancelrowinds[i1] < cancelrowinds[i2] ) 276 { 277 rowvarpair.varindex1 = cancelrowinds[i1]; 278 rowvarpair.varindex2 = cancelrowinds[i2]; 279 rowvarpair.varcoef1 = cancelrowvals[i1]; 280 rowvarpair.varcoef2 = cancelrowvals[i2]; 281 } 282 else 283 { 284 rowvarpair.varindex1 = cancelrowinds[i2]; 285 rowvarpair.varindex2 = cancelrowinds[i1]; 286 rowvarpair.varcoef1 = cancelrowvals[i2]; 287 rowvarpair.varcoef2 = cancelrowvals[i1]; 288 } 289 290 eqrowvarpair = (ROWVARPAIR*)SCIPhashtableRetrieve(pairtable, (void*) &rowvarpair); 291 nretrieves++; 292 293 if( eqrowvarpair == NULL || eqrowvarpair->rowindex == rowidx ) 294 continue; 295 296 /* if the row we want to cancel is an equality, we will only use equalities 297 * for canceling with less non-zeros and if the number of non-zeros is equal we use the 298 * rowindex as tie-breaker to avoid cyclic non-zero cancellation 299 */ 300 eqrowlen = SCIPmatrixGetRowNNonzs(matrix, eqrowvarpair->rowindex); 301 if( rowiseq && (cancelrowlen < eqrowlen || (cancelrowlen == eqrowlen && rowidx < eqrowvarpair->rowindex)) ) 302 continue; 303 304 eqrowvals = SCIPmatrixGetRowValPtr(matrix, eqrowvarpair->rowindex); 305 eqrowinds = SCIPmatrixGetRowIdxPtr(matrix, eqrowvarpair->rowindex); 306 307 scale = -rowvarpair.varcoef1 / eqrowvarpair->varcoef1; 308 309 if( REALABS(scale) > MAXSCALE ) 310 continue; 311 312 a = 0; 313 b = 0; 314 ncancel = 0; 315 316 ncontfillin = 0; 317 nintfillin = 0; 318 nbinfillin = 0; 319 abortpair = FALSE; 320 while( a < cancelrowlen && b < eqrowlen ) 321 { 322 if( cancelrowinds[a] == eqrowinds[b] ) 323 { 324 SCIP_Real newcoef; 325 326 newcoef = cancelrowvals[a] + scale * eqrowvals[b]; 327 328 /* check if coefficient is cancelled */ 329 if( SCIPisZero(scip, newcoef) ) 330 { 331 ++ncancel; 332 } 333 /* otherwise, check if integral coefficients are preserved if the column is integral */ 334 else if( (preserveintcoefs && SCIPvarIsIntegral(SCIPmatrixGetVar(matrix, cancelrowinds[a])) && 335 SCIPisIntegral(scip, cancelrowvals[a]) && !SCIPisIntegral(scip, newcoef)) ) 336 { 337 abortpair = TRUE; 338 break; 339 } 340 /* finally, check if locks could be modified in a bad way due to flipped signs */ 341 else if( (SCIPisInfinity(scip, cancelrhs) || SCIPisInfinity(scip, -cancellhs)) && 342 COPYSIGN(1.0, newcoef) != COPYSIGN(1.0, cancelrowvals[a]) ) /*lint !e777*/ 343 { 344 /* do not flip signs for non-canceled coefficients if this adds a lock to a variable that had at most one lock 345 * in that direction before, except if the other direction gets unlocked 346 */ 347 if( (cancelrowvals[a] > 0.0 && ! SCIPisInfinity(scip, cancelrhs)) || 348 (cancelrowvals[a] < 0.0 && ! SCIPisInfinity(scip, -cancellhs)) ) 349 { 350 /* if we get into this case the variable had a positive coefficient in a <= constraint or a negative 351 * coefficient in a >= constraint, e.g. an uplock. If this was the only uplock we do not abort their 352 * cancelling, otherwise we abort if we had a single or no downlock and add one 353 */ 354 if( SCIPmatrixGetColNUplocks(matrix, cancelrowinds[a]) > 1 && 355 SCIPmatrixGetColNDownlocks(matrix, cancelrowinds[a]) <= 1 ) 356 { 357 abortpair = TRUE; 358 break; 359 } 360 } 361 362 if( (cancelrowvals[a] < 0.0 && ! SCIPisInfinity(scip, cancelrhs)) || 363 (cancelrowvals[a] > 0.0 && ! SCIPisInfinity(scip, -cancellhs)) ) 364 { 365 /* symmetric case where the variable had a downlock */ 366 if( SCIPmatrixGetColNDownlocks(matrix, cancelrowinds[a]) > 1 && 367 SCIPmatrixGetColNUplocks(matrix, cancelrowinds[a]) <= 1 ) 368 { 369 abortpair = TRUE; 370 break; 371 } 372 } 373 } 374 375 ++a; 376 ++b; 377 } 378 else if( cancelrowinds[a] < eqrowinds[b] ) 379 { 380 ++a; 381 } 382 else 383 { 384 SCIP_Real newcoef; 385 SCIP_VAR* var; 386 387 var = SCIPmatrixGetVar(matrix, eqrowinds[b]); 388 newcoef = scale * eqrowvals[b]; 389 390 if( (newcoef > 0.0 && ! SCIPisInfinity(scip, cancelrhs)) || 391 (newcoef < 0.0 && ! SCIPisInfinity(scip, -cancellhs)) ) 392 { 393 if( SCIPmatrixGetColNUplocks(matrix, eqrowinds[b]) <= 1 ) 394 { 395 abortpair = TRUE; 396 ++b; 397 break; 398 } 399 } 400 401 if( (newcoef < 0.0 && ! SCIPisInfinity(scip, cancelrhs)) || 402 (newcoef > 0.0 && ! SCIPisInfinity(scip, -cancellhs)) ) 403 { 404 if( SCIPmatrixGetColNDownlocks(matrix, eqrowinds[b]) <= 1 ) 405 { 406 abortpair = TRUE; 407 ++b; 408 break; 409 } 410 } 411 412 ++b; 413 414 if( SCIPvarIsIntegral(var) ) 415 { 416 if( SCIPvarIsBinary(var) && ++nbinfillin > maxbinfillin ) 417 { 418 abortpair = TRUE; 419 break; 420 } 421 422 if( ++nintfillin > maxintfillin ) 423 { 424 abortpair = TRUE; 425 break; 426 } 427 } 428 else 429 { 430 if( ++ncontfillin > maxcontfillin ) 431 { 432 abortpair = TRUE; 433 break; 434 } 435 } 436 } 437 } 438 439 if( abortpair ) 440 continue; 441 442 while( b < eqrowlen ) 443 { 444 SCIP_VAR* var = SCIPmatrixGetVar(matrix, eqrowinds[b]); 445 ++b; 446 if( SCIPvarIsIntegral(var) ) 447 { 448 if( SCIPvarIsBinary(var) && ++nbinfillin > maxbinfillin ) 449 break; 450 if( ++nintfillin > maxintfillin ) 451 break; 452 } 453 else 454 { 455 if( ++ncontfillin > maxcontfillin ) 456 break; 457 } 458 } 459 460 if( ncontfillin > maxcontfillin || nbinfillin > maxbinfillin || nintfillin > maxintfillin ) 461 continue; 462 463 ntotfillin = nintfillin + ncontfillin; 464 465 if( ntotfillin >= ncancel ) 466 continue; 467 468 cancelrate = (ncancel - ntotfillin) / (SCIP_Real) eqrowlen; 469 470 if( cancelrate < mincancelrate ) 471 continue; 472 473 if( cancelrate > bestcancelrate ) 474 { 475 bestnfillin = ntotfillin; 476 bestcand = eqrowvarpair->rowindex; 477 bestscale = scale; 478 bestcancelrate = cancelrate; 479 480 /* stop looking if the current candidate does not create any fill-in or alter coefficients */ 481 if( cancelrate == 1.0 ) 482 break; 483 } 484 485 /* we accept the best candidate immediately if it does not create any fill-in or alter coefficients */ 486 if( bestcand != -1 && bestcancelrate == 1.0 ) 487 break; 488 } 489 } 490 491 if( bestcand != -1 ) 492 { 493 int a; 494 int b; 495 SCIP_Real* eqrowvals; 496 int* eqrowinds; 497 int eqrowlen; 498 int tmprowlen; 499 SCIP_Real eqrhs; 500 501 eqrowvals = SCIPmatrixGetRowValPtr(matrix, bestcand); 502 eqrowinds = SCIPmatrixGetRowIdxPtr(matrix, bestcand); 503 eqrowlen = SCIPmatrixGetRowNNonzs(matrix, bestcand); 504 eqrhs = SCIPmatrixGetRowRhs(matrix, bestcand); 505 506 a = 0; 507 b = 0; 508 tmprowlen = 0; 509 510 if( !SCIPisZero(scip, eqrhs) ) 511 { 512 if( !SCIPisInfinity(scip, -cancellhs) ) 513 cancellhs += bestscale * eqrhs; 514 if( !SCIPisInfinity(scip, cancelrhs) ) 515 cancelrhs += bestscale * eqrhs; 516 } 517 518 while( a < cancelrowlen && b < eqrowlen ) 519 { 520 if( cancelrowinds[a] == eqrowinds[b] ) 521 { 522 SCIP_Real val = cancelrowvals[a] + bestscale * eqrowvals[b]; 523 524 if( !SCIPisZero(scip, val) ) 525 { 526 tmpinds[tmprowlen] = cancelrowinds[a]; 527 tmpvals[tmprowlen] = val; 528 ++tmprowlen; 529 } 530 ++nchgcoef; 531 532 ++a; 533 ++b; 534 } 535 else if( cancelrowinds[a] < eqrowinds[b] ) 536 { 537 tmpinds[tmprowlen] = cancelrowinds[a]; 538 tmpvals[tmprowlen] = cancelrowvals[a]; 539 ++tmprowlen; 540 ++a; 541 } 542 else 543 { 544 tmpinds[tmprowlen] = eqrowinds[b]; 545 tmpvals[tmprowlen] = eqrowvals[b] * bestscale; 546 ++nchgcoef; 547 ++tmprowlen; 548 ++b; 549 } 550 } 551 552 while( a < cancelrowlen ) 553 { 554 tmpinds[tmprowlen] = cancelrowinds[a]; 555 tmpvals[tmprowlen] = cancelrowvals[a]; 556 ++tmprowlen; 557 ++a; 558 } 559 560 while( b < eqrowlen ) 561 { 562 tmpinds[tmprowlen] = eqrowinds[b]; 563 tmpvals[tmprowlen] = eqrowvals[b] * bestscale; 564 ++nchgcoef; 565 ++tmprowlen; 566 ++b; 567 } 568 569 /* update fill-in counter */ 570 *nfillin += bestnfillin; 571 572 /* swap the temporary arrays so that the cancelrowinds and cancelrowvals arrays, contain the new 573 * changed row, and the tmpinds and tmpvals arrays can be overwritten in the next iteration 574 */ 575 SCIPswapPointers((void**) &tmpinds, (void**) &cancelrowinds); 576 SCIPswapPointers((void**) &tmpvals, (void**) &cancelrowvals); 577 cancelrowlen = tmprowlen; 578 swapped = !swapped; 579 } 580 else 581 break; 582 } 583 584 if( nchgcoef != 0 ) 585 { 586 SCIP_CONS* cons; 587 SCIP_VAR** consvars; 588 589 int i; 590 591 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, cancelrowlen) ); 592 593 for( i = 0; i < cancelrowlen; ++i ) 594 consvars[i] = SCIPmatrixGetVar(matrix, cancelrowinds[i]); 595 596 /* create sparsified constraint and add it to scip */ 597 SCIP_CALL( SCIPcreateConsLinear(scip, &cons, SCIPmatrixGetRowName(matrix, rowidx), cancelrowlen, consvars, cancelrowvals, 598 cancellhs, cancelrhs, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) ); 599 SCIP_CALL( SCIPdelCons(scip, SCIPmatrixGetCons(matrix, rowidx)) ); 600 SCIP_CALL( SCIPaddCons(scip, cons) ); 601 602 #ifdef SCIP_MORE_DEBUG 603 SCIPdebugMsg(scip, "########\n"); 604 SCIPdebugMsg(scip, "old:\n"); 605 SCIPmatrixPrintRow(scip, matrix, rowidx); 606 SCIPdebugMsg(scip, "new:\n"); 607 SCIPdebugPrintCons(scip, cons, NULL); 608 SCIPdebugMsg(scip, "########\n"); 609 #endif 610 611 SCIP_CALL( SCIPreleaseCons(scip, &cons) ); 612 613 /* update counters */ 614 *nchgcoefs += nchgcoef; 615 *ncanceled += SCIPmatrixGetRowNNonzs(matrix, rowidx) - cancelrowlen; 616 617 /* if successful, decrease the useless hashtable retrieves counter; the rationale here is that we want to keep 618 * going if, after many useless calls that almost exceeded the budget, we finally reach a useful section; but we 619 * don't allow a negative build-up for the case that the useful section is all at the beginning and we just want 620 * to quit quickly afterwards 621 */ 622 *nuseless -= nretrieves; 623 *nuseless = MAX(*nuseless, 0); 624 625 SCIPfreeBufferArray(scip, &consvars); 626 } 627 else 628 { 629 /* if not successful, increase useless hashtable retrieves counter */ 630 *nuseless += nretrieves; 631 } 632 633 SCIPfreeBufferArray(scip, &locks); 634 if( !swapped ) 635 { 636 SCIPfreeBufferArray(scip, &tmpvals); 637 SCIPfreeBufferArray(scip, &tmpinds); 638 SCIPfreeBufferArray(scip, &cancelrowvals); 639 SCIPfreeBufferArray(scip, &cancelrowinds); 640 } 641 else 642 { 643 SCIPfreeBufferArray(scip, &cancelrowvals); 644 SCIPfreeBufferArray(scip, &cancelrowinds); 645 SCIPfreeBufferArray(scip, &tmpvals); 646 SCIPfreeBufferArray(scip, &tmpinds); 647 } 648 649 return SCIP_OKAY; 650 } 651 652 /** updates failure counter after one execution */ 653 static 654 void updateFailureStatistic( 655 SCIP_PRESOLDATA* presoldata, /**< presolver data */ 656 SCIP_Bool success /**< was this execution successful? */ 657 ) 658 { 659 assert(presoldata != NULL); 660 661 if( success ) 662 { 663 presoldata->nfailures = 0; 664 presoldata->nwaitingcalls = 0; 665 } 666 else 667 { 668 presoldata->nfailures++; 669 presoldata->nwaitingcalls = (int)(presoldata->waitingfac*(SCIP_Real)presoldata->nfailures); 670 } 671 } 672 673 674 /* 675 * Callback methods of presolver 676 */ 677 678 /** copy method for constraint handler plugins (called when SCIP copies plugins) */ 679 static 680 SCIP_DECL_PRESOLCOPY(presolCopySparsify) 681 { 682 SCIP_PRESOLDATA* presoldata; 683 684 assert(scip != NULL); 685 assert(presol != NULL); 686 assert(strcmp(SCIPpresolGetName(presol), PRESOL_NAME) == 0); 687 688 /* call inclusion method of presolver if copying is enabled */ 689 presoldata = SCIPpresolGetData(presol); 690 assert(presoldata != NULL); 691 if( presoldata->enablecopy ) 692 { 693 SCIP_CALL( SCIPincludePresolSparsify(scip) ); 694 } 695 696 return SCIP_OKAY; 697 } 698 699 /** execution method of presolver */ 700 static 701 SCIP_DECL_PRESOLEXEC(presolExecSparsify) 702 { /*lint --e{715}*/ 703 SCIP_MATRIX* matrix; 704 SCIP_Bool initialized; 705 SCIP_Bool complete; 706 SCIP_Bool infeasible; 707 int nrows; 708 int r; 709 int i; 710 int j; 711 int numcancel; 712 int oldnchgcoefs; 713 int nfillin; 714 int* locks; 715 int* perm; 716 int* rowidxsorted; 717 int* rowsparsity; 718 SCIP_HASHTABLE* pairtable; 719 ROWVARPAIR* varpairs; 720 int nvarpairs; 721 int varpairssize; 722 SCIP_PRESOLDATA* presoldata; 723 SCIP_Longint maxuseless; 724 SCIP_Longint nuseless; 725 SCIP_CONSHDLR* linearhdlr; 726 727 assert(result != NULL); 728 729 *result = SCIP_DIDNOTRUN; 730 731 if( (SCIPgetStage(scip) != SCIP_STAGE_PRESOLVING) || SCIPinProbing(scip) || SCIPisNLPEnabled(scip) ) 732 return SCIP_OKAY; 733 734 if( SCIPisStopped(scip) || SCIPgetNActivePricers(scip) > 0 ) 735 return SCIP_OKAY; 736 737 presoldata = SCIPpresolGetData(presol); 738 739 if( presoldata->nwaitingcalls > 0 ) 740 { 741 presoldata->nwaitingcalls--; 742 SCIPdebugMsg(scip, "skipping sparsify: nfailures=%d, nwaitingcalls=%d\n", presoldata->nfailures, 743 presoldata->nwaitingcalls); 744 return SCIP_OKAY; 745 } 746 747 /* if we want to cancel only from specialized constraints according to the parameter, then we can skip execution if 748 * only linear constraints are present 749 */ 750 linearhdlr = SCIPfindConshdlr(scip, "linear"); 751 if( !presoldata->cancellinear && linearhdlr != NULL && SCIPconshdlrGetNConss(linearhdlr) >= SCIPgetNConss(scip) ) 752 { 753 SCIPdebugMsg(scip, "skipping sparsify: only linear constraints found\n"); 754 return SCIP_OKAY; 755 } 756 757 SCIPdebugMsg(scip, "starting sparsify. . .\n"); 758 *result = SCIP_DIDNOTFIND; 759 760 matrix = NULL; 761 SCIP_CALL( SCIPmatrixCreate(scip, &matrix, TRUE, &initialized, &complete, &infeasible, 762 naddconss, ndelconss, nchgcoefs, nchgbds, nfixedvars) ); 763 764 /* if infeasibility was detected during matrix creation, return here */ 765 if( infeasible ) 766 { 767 if( initialized ) 768 SCIPmatrixFree(scip, &matrix); 769 770 *result = SCIP_CUTOFF; 771 return SCIP_OKAY; 772 } 773 774 /* we only work on pure MIPs currently */ 775 if( initialized && complete ) 776 { 777 nrows = SCIPmatrixGetNRows(matrix); 778 779 /* sort rows by column indices */ 780 for( i = 0; i < nrows; i++ ) 781 { 782 int* rowpnt = SCIPmatrixGetRowIdxPtr(matrix, i); 783 SCIP_Real* valpnt = SCIPmatrixGetRowValPtr(matrix, i); 784 SCIPsortIntReal(rowpnt, valpnt, SCIPmatrixGetRowNNonzs(matrix, i)); 785 } 786 787 SCIP_CALL( SCIPallocBufferArray(scip, &locks, SCIPmatrixGetNColumns(matrix)) ); 788 SCIP_CALL( SCIPallocBufferArray(scip, &perm, SCIPmatrixGetNColumns(matrix)) ); 789 790 /* loop over all rows and create var pairs */ 791 numcancel = 0; 792 nfillin = 0; 793 varpairssize = 0; 794 nvarpairs = 0; 795 varpairs = NULL; 796 SCIP_CALL( SCIPhashtableCreate(&pairtable, SCIPblkmem(scip), 1, SCIPhashGetKeyStandard, varPairsEqual, varPairHashval, (void*) scip) ); 797 798 /* collect equalities and their number of non-zeros */ 799 for( r = 0; r < nrows; r++ ) 800 { 801 int nnonz; 802 803 nnonz = SCIPmatrixGetRowNNonzs(matrix, r); 804 805 /* consider equalities with support at most maxnonzeros; skip singleton equalities, because these are faster 806 * processed by trivial presolving 807 */ 808 if( nnonz >= 2 && (presoldata->maxnonzeros < 0 || nnonz <= presoldata->maxnonzeros) 809 && SCIPisEQ(scip, SCIPmatrixGetRowRhs(matrix, r), SCIPmatrixGetRowLhs(matrix, r)) ) 810 { 811 int* rowinds; 812 SCIP_Real* rowvals; 813 int npairs; 814 int failshift; 815 816 rowinds = SCIPmatrixGetRowIdxPtr(matrix, r); 817 rowvals = SCIPmatrixGetRowValPtr(matrix, r); 818 819 for( i = 0; i < nnonz; ++i ) 820 { 821 perm[i] = i; 822 locks[i] = SCIPmatrixGetColNDownlocks(matrix, rowinds[i]) + SCIPmatrixGetColNUplocks(matrix, rowinds[i]); 823 } 824 825 SCIPsortIntInt(locks, perm, nnonz); 826 827 if( presoldata->maxconsiderednonzeros >= 0 ) 828 nnonz = MIN(nnonz, presoldata->maxconsiderednonzeros); 829 830 npairs = (nnonz * (nnonz - 1)) / 2; 831 if( nvarpairs + npairs > varpairssize ) 832 { 833 int newsize = SCIPcalcMemGrowSize(scip, nvarpairs + npairs); 834 SCIP_CALL( SCIPreallocBufferArray(scip, &varpairs, newsize) ); 835 varpairssize = newsize; 836 } 837 838 /* if we are called after one or more failures, i.e., executions without finding cancellations, then we 839 * shift the section of nonzeros considered; in the case that the maxconsiderednonzeros limit is hit, this 840 * results in different variable pairs being tried and avoids trying the same useless cancellations 841 * repeatedly 842 */ 843 failshift = presoldata->nfailures*presoldata->maxconsiderednonzeros; 844 845 for( i = 0; i < nnonz; ++i ) 846 { 847 for( j = i + 1; j < nnonz; ++j ) 848 { 849 int i1; 850 int i2; 851 852 assert(nvarpairs < varpairssize); 853 assert(varpairs != NULL); 854 855 i1 = perm[(i + failshift) % nnonz]; 856 i2 = perm[(j + failshift) % nnonz]; 857 varpairs[nvarpairs].rowindex = r; 858 859 if( rowinds[i1] < rowinds[i2]) 860 { 861 varpairs[nvarpairs].varindex1 = rowinds[i1]; 862 varpairs[nvarpairs].varindex2 = rowinds[i2]; 863 varpairs[nvarpairs].varcoef1 = rowvals[i1]; 864 varpairs[nvarpairs].varcoef2 = rowvals[i2]; 865 } 866 else 867 { 868 varpairs[nvarpairs].varindex1 = rowinds[i2]; 869 varpairs[nvarpairs].varindex2 = rowinds[i1]; 870 varpairs[nvarpairs].varcoef1 = rowvals[i2]; 871 varpairs[nvarpairs].varcoef2 = rowvals[i1]; 872 } 873 ++nvarpairs; 874 } 875 } 876 } 877 } 878 879 /* insert varpairs into hash table */ 880 for( r = 0; r < nvarpairs; ++r ) 881 { 882 SCIP_Bool insert; 883 ROWVARPAIR* othervarpair; 884 885 assert(varpairs != NULL); 886 887 insert = TRUE; 888 889 /* check if this pair is already contained in the hash table; 890 * The loop is required due to the non-transitivity of the hash functions 891 */ 892 while( (othervarpair = (ROWVARPAIR*)SCIPhashtableRetrieve(pairtable, (void*) &varpairs[r])) != NULL ) 893 { 894 /* if the previous variable pair has fewer or the same number of non-zeros in the attached row 895 * we keep that pair and skip this one 896 */ 897 if( SCIPmatrixGetRowNNonzs(matrix, othervarpair->rowindex) <= SCIPmatrixGetRowNNonzs(matrix, varpairs[r].rowindex) ) 898 { 899 insert = FALSE; 900 break; 901 } 902 903 /* this pairs row has fewer non-zeros, so remove the other pair from the hash table and loop */ 904 SCIP_CALL( SCIPhashtableRemove(pairtable, (void*) othervarpair) ); 905 } 906 907 if( insert ) 908 { 909 /* prevent the insertion of too many variable pairs into the hashtable. 910 * a safety margin of factor 4 is built into the 8=2*4. 911 */ 912 if( ((SCIP_Longint)SCIPhashtableGetNEntries(pairtable) * (SCIP_Longint)(8 * sizeof(void*))) > (SCIP_Longint)INT_MAX ) 913 { 914 break; 915 } 916 917 SCIP_CALL( SCIPhashtableInsert(pairtable, (void*) &varpairs[r]) ); 918 } 919 } 920 921 /* sort rows according to parameter value */ 922 if( presoldata->rowsort == 'i' || presoldata->rowsort == 'd' ) 923 { 924 SCIP_CALL( SCIPallocBufferArray(scip, &rowidxsorted, nrows) ); 925 SCIP_CALL( SCIPallocBufferArray(scip, &rowsparsity, nrows) ); 926 for( r = 0; r < nrows; ++r ) 927 rowidxsorted[r] = r; 928 if( presoldata->rowsort == 'i' ) 929 { 930 for( r = 0; r < nrows; ++r ) 931 rowsparsity[r] = SCIPmatrixGetRowNNonzs(matrix, r); 932 } 933 else if( presoldata->rowsort == 'd' ) 934 { 935 for( r = 0; r < nrows; ++r ) 936 rowsparsity[r] = -SCIPmatrixGetRowNNonzs(matrix, r); 937 } 938 SCIPsortIntInt(rowsparsity, rowidxsorted, nrows); 939 } 940 else 941 { 942 assert(presoldata->rowsort == 'n'); 943 rowidxsorted = NULL; 944 rowsparsity = NULL; 945 } 946 947 /* loop over the rows and cancel non-zeros until maximum number of retrieves is reached */ 948 maxuseless = (SCIP_Longint)(presoldata->maxretrievefac * (SCIP_Real)nrows); 949 nuseless = 0; 950 oldnchgcoefs = *nchgcoefs; 951 for( r = 0; r < nrows && nuseless <= maxuseless && !SCIPisStopped(scip); r++ ) 952 { 953 int rowidx; 954 955 rowidx = rowidxsorted != NULL ? rowidxsorted[r] : r; 956 957 /* check whether we want to cancel only from specialized constraints; one reasoning behind this may be that 958 * cancelling fractional coefficients requires more numerical care than is currently implemented in method 959 * cancelRow() 960 */ 961 assert(SCIPmatrixGetCons(matrix, rowidx) != NULL); 962 if( !presoldata->cancellinear && SCIPconsGetHdlr(SCIPmatrixGetCons(matrix, rowidx)) == linearhdlr ) 963 continue; 964 965 /* since the function parameters for the max fillin are unsigned we do not need to handle the 966 * unlimited (-1) case due to implicit conversion rules */ 967 SCIP_CALL( cancelRow(scip, matrix, pairtable, rowidx, \ 968 presoldata->maxcontfillin == -1 ? INT_MAX : presoldata->maxcontfillin, \ 969 presoldata->maxintfillin == -1 ? INT_MAX : presoldata->maxintfillin, \ 970 presoldata->maxbinfillin == -1 ? INT_MAX : presoldata->maxbinfillin, \ 971 presoldata->maxconsiderednonzeros, presoldata->preserveintcoefs, \ 972 &nuseless, nchgcoefs, &numcancel, &nfillin) ); 973 } 974 975 SCIPfreeBufferArrayNull(scip, &rowsparsity); 976 SCIPfreeBufferArrayNull(scip, &rowidxsorted); 977 978 SCIPhashtableFree(&pairtable); 979 SCIPfreeBufferArrayNull(scip, &varpairs); 980 981 SCIPfreeBufferArray(scip, &perm); 982 SCIPfreeBufferArray(scip, &locks); 983 984 /* update result */ 985 presoldata->ncancels += numcancel; 986 presoldata->nfillin += nfillin; 987 988 if( numcancel > 0 ) 989 { 990 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 991 " (%.1fs) sparsify %s: %d/%d (%.1f%%) nonzeros canceled" 992 " - in total %d canceled nonzeros, %d changed coefficients, %d added nonzeros\n", 993 SCIPgetSolvingTime(scip), (nuseless > maxuseless ? "aborted" : "finished"), numcancel, 994 SCIPmatrixGetNNonzs(matrix), 100.0*(SCIP_Real)numcancel/(SCIP_Real)SCIPmatrixGetNNonzs(matrix), 995 presoldata->ncancels, SCIPpresolGetNChgCoefs(presol) + *nchgcoefs - oldnchgcoefs, presoldata->nfillin); 996 *result = SCIP_SUCCESS; 997 } 998 999 updateFailureStatistic(presoldata, numcancel > 0); 1000 1001 SCIPdebugMsg(scip, "sparsify failure statistic: nfailures=%d, nwaitingcalls=%d\n", presoldata->nfailures, 1002 presoldata->nwaitingcalls); 1003 } 1004 /* if matrix construction fails once, we do not ever want to be called again */ 1005 else 1006 { 1007 updateFailureStatistic(presoldata, FALSE); 1008 presoldata->nwaitingcalls = INT_MAX; 1009 } 1010 1011 SCIPmatrixFree(scip, &matrix); 1012 1013 return SCIP_OKAY; 1014 } 1015 1016 /* 1017 * presolver specific interface methods 1018 */ 1019 1020 /** destructor of presolver to free user data (called when SCIP is exiting) */ 1021 static 1022 SCIP_DECL_PRESOLFREE(presolFreeSparsify) 1023 { /*lint --e{715}*/ 1024 SCIP_PRESOLDATA* presoldata; 1025 1026 /* free presolver data */ 1027 presoldata = SCIPpresolGetData(presol); 1028 assert(presoldata != NULL); 1029 1030 SCIPfreeBlockMemory(scip, &presoldata); 1031 SCIPpresolSetData(presol, NULL); 1032 1033 return SCIP_OKAY; 1034 } 1035 1036 /** initialization method of presolver (called after problem was transformed) */ 1037 static 1038 SCIP_DECL_PRESOLINIT(presolInitSparsify) 1039 { 1040 SCIP_PRESOLDATA* presoldata; 1041 1042 /* set the counters in the init (and not in the initpre) callback such that they persist across restarts */ 1043 presoldata = SCIPpresolGetData(presol); 1044 presoldata->ncancels = 0; 1045 presoldata->nfillin = 0; 1046 presoldata->nfailures = 0; 1047 presoldata->nwaitingcalls = 0; 1048 1049 return SCIP_OKAY; 1050 } 1051 1052 /** creates the sparsify presolver and includes it in SCIP */ 1053 SCIP_RETCODE SCIPincludePresolSparsify( 1054 SCIP* scip /**< SCIP data structure */ 1055 ) 1056 { 1057 SCIP_PRESOLDATA* presoldata; 1058 SCIP_PRESOL* presol; 1059 1060 /* create sparsify presolver data */ 1061 SCIP_CALL( SCIPallocBlockMemory(scip, &presoldata) ); 1062 1063 /* include presolver */ 1064 SCIP_CALL( SCIPincludePresolBasic(scip, &presol, PRESOL_NAME, PRESOL_DESC, PRESOL_PRIORITY, PRESOL_MAXROUNDS, 1065 PRESOL_TIMING, presolExecSparsify, presoldata) ); 1066 1067 SCIP_CALL( SCIPsetPresolCopy(scip, presol, presolCopySparsify) ); 1068 SCIP_CALL( SCIPsetPresolFree(scip, presol, presolFreeSparsify) ); 1069 SCIP_CALL( SCIPsetPresolInit(scip, presol, presolInitSparsify) ); 1070 1071 SCIP_CALL( SCIPaddBoolParam(scip, 1072 "presolving/sparsify/enablecopy", 1073 "should sparsify presolver be copied to sub-SCIPs?", 1074 &presoldata->enablecopy, TRUE, DEFAULT_ENABLECOPY, NULL, NULL) ); 1075 1076 SCIP_CALL( SCIPaddBoolParam(scip, 1077 "presolving/sparsify/cancellinear", 1078 "should we cancel nonzeros in constraints of the linear constraint handler?", 1079 &presoldata->cancellinear, TRUE, DEFAULT_CANCELLINEAR, NULL, NULL) ); 1080 1081 SCIP_CALL( SCIPaddBoolParam(scip, 1082 "presolving/sparsify/preserveintcoefs", 1083 "should we forbid cancellations that destroy integer coefficients?", 1084 &presoldata->preserveintcoefs, TRUE, DEFAULT_PRESERVEINTCOEFS, NULL, NULL) ); 1085 1086 SCIP_CALL( SCIPaddIntParam(scip, 1087 "presolving/sparsify/maxcontfillin", 1088 "maximal fillin for continuous variables (-1: unlimited)", 1089 &presoldata->maxcontfillin, FALSE, DEFAULT_MAX_CONT_FILLIN, -1, INT_MAX, NULL, NULL) ); 1090 1091 SCIP_CALL( SCIPaddIntParam(scip, 1092 "presolving/sparsify/maxbinfillin", 1093 "maximal fillin for binary variables (-1: unlimited)", 1094 &presoldata->maxbinfillin, FALSE, DEFAULT_MAX_BIN_FILLIN, -1, INT_MAX, NULL, NULL) ); 1095 1096 SCIP_CALL( SCIPaddIntParam(scip, 1097 "presolving/sparsify/maxintfillin", 1098 "maximal fillin for integer variables including binaries (-1: unlimited)", 1099 &presoldata->maxintfillin, FALSE, DEFAULT_MAX_INT_FILLIN, -1, INT_MAX, NULL, NULL) ); 1100 1101 SCIP_CALL( SCIPaddIntParam(scip, 1102 "presolving/sparsify/maxnonzeros", 1103 "maximal support of one equality to be used for cancelling (-1: no limit)", 1104 &presoldata->maxnonzeros, TRUE, DEFAULT_MAXNONZEROS, -1, INT_MAX, NULL, NULL) ); 1105 1106 SCIP_CALL( SCIPaddIntParam(scip, 1107 "presolving/sparsify/maxconsiderednonzeros", 1108 "maximal number of considered non-zeros within one row (-1: no limit)", 1109 &presoldata->maxconsiderednonzeros, TRUE, DEFAULT_MAXCONSIDEREDNONZEROS, -1, INT_MAX, NULL, NULL) ); 1110 1111 SCIP_CALL( SCIPaddCharParam(scip, 1112 "presolving/sparsify/rowsort", 1113 "order in which to process inequalities ('n'o sorting, 'i'ncreasing nonzeros, 'd'ecreasing nonzeros)", 1114 &presoldata->rowsort, TRUE, DEFAULT_ROWSORT, "nid", NULL, NULL) ); 1115 1116 SCIP_CALL( SCIPaddRealParam(scip, 1117 "presolving/sparsify/maxretrievefac", 1118 "limit on the number of useless vs. useful hashtable retrieves as a multiple of the number of constraints", 1119 &presoldata->maxretrievefac, TRUE, DEFAULT_MAXRETRIEVEFAC, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 1120 1121 SCIP_CALL( SCIPaddRealParam(scip, 1122 "presolving/sparsify/waitingfac", 1123 "number of calls to wait until next execution as a multiple of the number of useless calls", 1124 &presoldata->waitingfac, TRUE, DEFAULT_WAITINGFAC, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 1125 1126 return SCIP_OKAY; 1127 } 1128