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