1 /* $Id: CoinPresolveDual.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include <stdio.h>
7 #include <math.h>
8
9 #include "CoinPresolveMatrix.hpp"
10 #include "CoinPresolveFixed.hpp"
11 #include "CoinPresolveDual.hpp"
12 #include "CoinMessage.hpp"
13 #include "CoinHelperFunctions.hpp"
14 #include "CoinFloatEqual.hpp"
15
16 /*
17 Define PRESOLVE_DEBUG and PRESOLVE_CONSISTENCY as compile flags! If not
18 uniformly on across all uses of presolve code you'll get something between
19 garbage and a core dump. See comments in CoinPresolvePsdebug.hpp
20 */
21 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
22 #include "CoinPresolvePsdebug.hpp"
23 #endif
24
25 /*
26 Set to > 0 to enable bound propagation on dual variables? This has been
27 disabled for a good while. See notes in code.
28
29 #define PRESOLVE_TIGHTEN_DUALS 1
30 */
31 /*
32 Guards incorrect code that attempts to adjust a column solution. See
33 comments with the code. Could possibly be fixed, which is why it hasn't
34 been chopped out.
35
36 #define REMOVE_DUAL_ACTION_REPAIR_SOLN 1
37 */
38
39 /*
40 In this transform we're looking to prove bounds on the duals y<i> and
41 reduced costs cbar<j>. We can use this in two ways.
42
43 First, if we can prove a bound on the reduced cost cbar<j> with strict
44 inequality,
45 * cbar<j> > 0 ==> x<j> NBLB at optimality
46 * cbar<j> < 0 ==> x<j> NBUB at optimality
47 If the required bound is not finite, the problem is unbounded. Andersen &
48 Andersen call this a dominated column.
49
50 Second, suppose we can show that cbar<i> = -y<i> is strictly nonzero
51 for the logical s<i> associated with some inequality i. (There must be
52 exactly one finite row bound, so that the logical has exactly one finite
53 bound which is 0). If cbar<i> demands the logical be nonbasic at bound,
54 it's zero, and we can convert the inequality to an equality.
55
56 Not based on duals and reduced costs, but in the same vein, if we can
57 identify an architectural x<j> that'll accomplish the same goal (bring
58 one constraint tight when moved in a favourable direction), we can make
59 that constraint an equality. The conditions are different, however:
60 * Moving x<j> in the favourable direction tightens exactly one
61 constraint.
62 * The bound (l<j> or u<j>) in the favourable direction is infinite.
63 * If x<j> is entangled in any other constraints, moving x<j> in the
64 favourable direction loosens those constraints.
65 A bit of thought and linear algebra is required to convince yourself that in
66 the circumstances, cbar<j> = c<j>, hence we need only look at c<j> to
67 determine the favourable direction.
68
69
70 To get from bounds on the duals to bounds on the reduced costs, note that
71 cbar<j> = c<j> - (SUM{P}a<ij>y<i> + SUM{M}a<ij>y<i>)
72 for P = {a<ij> > 0} and M = {a<ij> < 0}. Then
73 cbarmin<j> = c<j> - (SUM{P}a<ij>ymax<i> + SUM{M}a<ij>ymin<i>)
74 cbarmax<j> = c<j> - (SUM{P}a<ij>ymin<i> + SUM{M}a<ij>ymax<i>)
75
76 As A&A note, the reverse implication also holds:
77 * if l<j> = -infty, cbar<j> <= 0 at optimality
78 * if u<j> = infty, cbar<j> >= 0 at optimality
79
80 We can use this to run bound propagation on the duals in an attempt to
81 tighten bounds and force more reduced costs to strict inequality. Suppose
82 u<j> = infty. Then cbar<j> >= 0. It must be possible to achieve this, so
83 cbarmax<j> >= 0:
84 0 <= c<j> - (SUM{P}a<ij>ymin<i> + SUM{M}a<ij>ymax<i>)
85 Solve for y<t>, a<tj> > 0:
86 y<t> <= 1/a<tj>[c<j> - (SUM{P\t}a<ij>ymin<i> + SUM{M}a<ij>ymax<i>)]
87 If a<tj> < 0,
88 y<t> >= 1/a<tj>[c<j> - (SUM{P}a<ij>ymin<i> + SUM{M\t}a<ij>ymax<i>)]
89 For l<j> = -infty, cbar<j> <= 0, hence cbarmin <= 0:
90 0 >= c<j> - (SUM{P}a<ij>ymax<i> + SUM{M}a<ij>ymin<i>)
91 Solve for y<t>, a<tj> > 0:
92 y<t> >= 1/a<tj>[c<j> - (SUM{P\t}a<ij>ymax<i> + SUM{M}a<ij>ymin<i>)]
93 If a<tj> < 0,
94 y<t> <= 1/a<tj>[c<j> - (SUM{P}a<ij>ymax<i> + SUM{M\t}a<ij>ymin<i>)]
95
96 We can get initial bounds on ymin<i> and ymax<i> from column singletons
97 x<j>, where
98 cbar<j> = c<j> - a<ij>y<i>
99 If u<j> = infty, then at optimality
100 0 <= cbar<j> = c<j> - a<ij>y<i>
101 For a<ij> > 0 we have
102 y<i> <= c<j>/a<ij>
103 and for a<ij> < 0 we have
104 y<i> >= c<j>/a<ij>
105 We can do a similar calculation for l<j> = -infty, so that
106 0 >= cbar<j> = c<j> - a<ij>y<i>
107 For a<ij> > 0 we have
108 y<i> >= c<j>/a<ij>
109 and for a<ij> < 0 we have
110 y<i> <= c<j>/a<ij>
111 A logical is a column singleton with c<j> = 0.0 and |a<ij>| = 1.0. One or
112 both bounds can be finite, depending on constraint representation.
113
114 This code is hardwired for minimisation.
115
116 Extensive rework of the second part of the routine Fall 2011 to add a
117 postsolve transform to fix bug #67, incorrect status for logicals.
118 -- lh, 111208 --
119 */
120 /*
121 Original comments from 040916:
122
123 This routine looks to be something of a work in progress. Down in the
124 bound propagation loop, why do we only work with variables with u_j =
125 infty? The corresponding section of code for l_j = -infty is ifdef'd away.
126 l<j> = -infty is uncommon; perhaps it's not worth the effort?
127
128 Why exclude the code protected by PRESOLVE_TIGHTEN_DUALS? Why are we
129 using ekkinf instead of PRESOLVE_INF?
130 */
131 const CoinPresolveAction
132 *
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)133 remove_dual_action::presolve(CoinPresolveMatrix *prob,
134 const CoinPresolveAction *next)
135 {
136 #if PRESOLVE_DEBUG > 0
137 std::cout
138 << "Entering remove_dual_action::presolve, " << prob->nrows_
139 << "x" << prob->ncols_ << "." << std::endl;
140 #endif
141 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
142 presolve_check_sol(prob, 2, 1, 1);
143 presolve_check_nbasic(prob);
144 #endif
145 #if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
146 int startEmptyRows = 0;
147 int startEmptyColumns = 0;
148 startEmptyRows = prob->countEmptyRows();
149 startEmptyColumns = prob->countEmptyCols();
150 #if COIN_PRESOLVE_TUNING > 0
151 double startTime = 0.0;
152 if (prob->tuning_)
153 startTime = CoinCpuTime();
154 #endif
155 #endif
156
157 // column-major representation
158 const int ncols = prob->ncols_;
159 const CoinBigIndex *const mcstrt = prob->mcstrt_;
160 const int *const hincol = prob->hincol_;
161 const int *const hrow = prob->hrow_;
162 const double *const colels = prob->colels_;
163 const double *const cost = prob->cost_;
164
165 // column type, bounds, solution, and status
166 const unsigned char *const integerType = prob->integerType_;
167 const double *const clo = prob->clo_;
168 const double *const cup = prob->cup_;
169 double *const csol = prob->sol_;
170 unsigned char *&colstat = prob->colstat_;
171
172 // row-major representation
173 const int nrows = prob->nrows_;
174 const CoinBigIndex *const mrstrt = prob->mrstrt_;
175 const int *const hinrow = prob->hinrow_;
176 const int *const hcol = prob->hcol_;
177 #if 1 //REMOVE_DUAL_ACTION_REPAIR_SOLN > 0
178 const double *const rowels = prob->rowels_;
179 #endif
180
181 // row bounds
182 double *const rlo = prob->rlo_;
183 double *const rup = prob->rup_;
184
185 // tolerances
186 const double ekkinf2 = PRESOLVE_SMALL_INF;
187 const double ekkinf = ekkinf2 * 1.0e8;
188 const double ztolcbarj = prob->ztoldj_;
189 const CoinRelFltEq relEq(prob->ztolzb_);
190
191 /*
192 Grab one of the preallocated scratch arrays to hold min and max values of
193 row duals.
194 */
195 double *ymin = prob->usefulRowDouble_;
196 double *ymax = ymin + nrows;
197 /* use arrays to say which bound active
198 1 - lower
199 2 - upper
200 4 - fixed
201 8 - don't touch
202 */
203 char *active = reinterpret_cast< char * >(prob->usefulColumnInt_);
204 memset(active, 0, ncols);
205 int nOneBound = 0;
206 int numberLook = prob->numberColsToDo_;
207 int *look = prob->colsToDo_;
208 for (int iLook = 0; iLook < numberLook; iLook++) {
209 int j = look[iLook];
210 char type;
211 if (cup[j] >= ekkinf) {
212 if (clo[j] <= -ekkinf)
213 type = 0;
214 else
215 type = 1;
216 } else {
217 if (clo[j] <= -ekkinf)
218 type = 2;
219 else if (clo[j] < cup[j])
220 type = 3;
221 else
222 type = 4;
223 }
224 if (type == 1 || type == 2)
225 nOneBound++;
226 active[j] = type;
227 }
228 int nFreed = 0;
229 #define USE_ACTIVE 1
230 //#define PRESOLVE_DEBUG 2
231 for (int i = 0; i < nrows; i++) {
232 const CoinBigIndex krs = mrstrt[i];
233 const CoinBigIndex kre = krs + hinrow[i];
234 int positive = 0;
235 int negative = 0;
236 CoinBigIndex onlyPositive = -1;
237 CoinBigIndex onlyNegative = -1;
238 for (CoinBigIndex k = krs; k < kre; k++) {
239 const double coeff = rowels[k];
240 const int icol = hcol[k];
241 char type = active[icol];
242 //const double lb = clo[icol] ;
243 //const double ub = cup[icol] ;
244 if (type == 0 || (type & 8) != 0) {
245 // free or already used
246 positive = 2;
247 negative = 2;
248 break;
249 } else if (type == 1) {
250 // lower bound active
251 if (coeff > 0.0) {
252 positive++;
253 } else {
254 negative++;
255 }
256 } else if (type == 2) {
257 // upper bound active
258 if (coeff < 0.0) {
259 positive++;
260 } else {
261 negative++;
262 }
263 } else if (type == 3) {
264 // both bounds active
265 if (coeff > 0.0) {
266 onlyPositive = k;
267 positive++;
268 } else {
269 onlyNegative = k;
270 negative++;
271 }
272 } else {
273 // fixed
274 }
275 }
276 bool doneSomething = false;
277 if (onlyPositive >= 0 && positive == 1) {
278 const double coeff = rowels[onlyPositive];
279 const int icol = hcol[onlyPositive];
280 const double lb = clo[icol];
281 const double ub = cup[icol];
282 // Get possible without icol
283 double maxUp = 0.0;
284 double maxDown = 0.0;
285 #if USE_ACTIVE > 2
286 printf("%d pos coeff %g bounds %g,%g rhs %g,%g - %d negative",
287 icol, coeff, lb, ub, rlo[i], rup[i], negative);
288 printf(", row is %g <= ", rlo[i]);
289 #endif
290 for (CoinBigIndex k = krs; k < kre; k++) {
291 const int icol2 = hcol[k];
292 #if USE_ACTIVE > 2
293 printf("(%d,%g - bounds %g,%g) ", icol2, rowels[k],
294 clo[icol2], cup[icol2]);
295 #endif
296 if (icol2 == icol)
297 continue;
298 const double coeff2 = rowels[k];
299 const double lb2 = clo[icol2];
300 const double ub2 = cup[icol2];
301 if (coeff2 > 0.0) {
302 if (ub2 > 1.0e30)
303 maxUp = COIN_DBL_MAX;
304 else
305 maxUp += coeff2 * ub2;
306 if (lb2 < -1.0e30)
307 maxDown = -COIN_DBL_MAX;
308 else
309 maxDown += coeff2 * lb2;
310 } else {
311 if (lb2 < -1.0e30)
312 maxUp = COIN_DBL_MAX;
313 else
314 maxUp += coeff2 * lb2;
315 if (ub2 > 1.0e30)
316 maxDown = -COIN_DBL_MAX;
317 else
318 maxDown += coeff2 * ub2;
319 }
320 }
321 double impliedLower = (rlo[i] - maxUp) / coeff;
322 double impliedUpper = (rup[i] - maxDown) / coeff;
323 #if USE_ACTIVE > 2
324 printf("<= %g - implied l,u %g,%g\n", rup[i], impliedLower, impliedUpper);
325 #endif
326 if (impliedLower > lb - 1.0e-7) {
327 #if USE_ACTIVE == 2
328 printf("%d pos coeff %g bounds %g,%g maxDown %g rhs %g,%g - %d negative\n",
329 icol, coeff, lb, ub, maxDown, rlo[i], rup[i], negative);
330 #endif
331 #if USE_ACTIVE > 1
332 printf("can take off lb as implied lower %g > %g\n", impliedLower, lb);
333 #endif
334 doneSomething = true;
335 active[icol] = 2 + 8;
336 nFreed++;
337 #define TRY_UPPER 3
338 #if TRY_UPPER > 1
339 } else if (impliedUpper < ub + 1.0e-7) {
340 #if 0
341 printf("%d pos coeff %g bounds %g,%g maxDown %g, maxUp %g rhs %g,%g\n",
342 icol,coeff,lb,ub,maxDown,maxUp,rlo[i],rup[i]);
343 printf("row is %g <= ",rlo[i]);
344 for (CoinBigIndex k = krs ; k < kre ; k++) {
345 const int icol2 = hcol[k] ;
346 printf("(%d,%g - bounds %g,%g) ",icol2,rowels[k],
347 clo[icol2],cup[icol2]);
348 }
349 printf("<= %g - implied l,u %g,%g\n",rup[i],impliedLower,impliedUpper);
350 printf("can take off ub as implied upper %g < %g\n",impliedUpper,ub);
351 #endif
352 #if TRY_UPPER > 2
353 doneSomething = true;
354 active[icol] = 1 + 8;
355 nFreed++;
356 #endif
357 #endif
358 } else {
359 #if 0
360 printf("cant take off lb as implied lower %g - ",impliedLower);
361 printf("row is %g <= ",rlo[i]);
362 for (CoinBigIndex k = krs ; k < kre ; k++) {
363 const double coeff = rowels[k] ;
364 const int icol = hcol[k] ;
365 printf("(%d,%g) ",icol,coeff);
366 }
367 printf("<= %g\n",rup[i]);
368 #endif
369 }
370 }
371 if (onlyNegative >= 0 && negative == 1 && !doneSomething) {
372 const double coeff = rowels[onlyNegative];
373 const int icol = hcol[onlyNegative];
374 const double lb = clo[icol];
375 const double ub = cup[icol];
376 // Get possible without icol
377 double maxUp = 0.0;
378 double maxDown = 0.0;
379 #if USE_ACTIVE > 2
380 printf("%d neg coeff %g bounds %g,%g rhs %g,%g - %d positive",
381 icol, coeff, lb, ub, rlo[i], rup[i], positive);
382 printf(", row is %g <= ", rlo[i]);
383 #endif
384 for (CoinBigIndex k = krs; k < kre; k++) {
385 const int icol2 = hcol[k];
386 #if USE_ACTIVE > 2
387 printf("(%d,%g - bounds %g,%g) ", icol2, rowels[k],
388 clo[icol2], cup[icol2]);
389 #endif
390 if (icol2 == icol)
391 continue;
392 const double coeff2 = rowels[k];
393 const double lb2 = clo[icol2];
394 const double ub2 = cup[icol2];
395 if (coeff2 > 0.0) {
396 if (ub2 > 1.0e30)
397 maxUp = COIN_DBL_MAX;
398 else
399 maxUp += coeff2 * ub2;
400 if (lb2 < -1.0e30)
401 maxDown = -COIN_DBL_MAX;
402 else
403 maxDown += coeff2 * lb2;
404 } else {
405 if (lb2 < -1.0e30)
406 maxUp = COIN_DBL_MAX;
407 else
408 maxUp += coeff2 * lb2;
409 if (ub2 > 1.0e30)
410 maxDown = -COIN_DBL_MAX;
411 else
412 maxDown += coeff2 * ub2;
413 }
414 }
415 double impliedLower = (rup[i] - maxDown) / coeff;
416 double impliedUpper = (rlo[i] - maxUp) / coeff;
417 #if USE_ACTIVE > 2
418 printf("<= %g - implied l,u %g,%g\n", rup[i], impliedLower, impliedUpper);
419 #endif
420 int nOne = 0;
421 for (CoinBigIndex k = krs; k < kre; k++) {
422 const double coeff = rowels[k];
423 const int icol = hcol[k];
424 if (hincol[icol] == 1 && coeff > 0.0)
425 nOne++;
426 }
427 if (nOne >= hinrow[i] - 1) {
428 for (CoinBigIndex k = krs; k < kre; k++) {
429 const int icol = hcol[k];
430 const double coeff = rowels[k];
431 if (hincol[icol] == 1 && coeff > 0.0) {
432 #if USE_ACTIVE > 2
433 printf("(%d has one entry) ", icol);
434 #endif
435 if (active[icol] == 3 || clo[icol] < -1.0e30)
436 nOne = 0; // no good
437 }
438 }
439 #if USE_ACTIVE > 2
440 printf("\n");
441 #endif
442 }
443 if (impliedLower > lb - 1.0e-7 /*&& nOne*/) {
444 #if USE_ACTIVE > 1
445 printf("second type can take off lb of %g on column %d as implied lower %g - effective upper %g, coeff %g on row %d\n", lb, icol, impliedLower, maxUp, coeff, i);
446 #endif
447 active[icol] = 2 + 8;
448 nFreed++;
449 #if 0 //TRY_UPPER
450 } else if (impliedLower>lb-1.0e-7 ) {
451 printf("second type can't take off lb of %g on column %d as implied lower %g - effective upper %g, coeff %g on row %d\n",lb,icol,impliedLower,maxUp,coeff,i);
452 #endif
453 #if TRY_UPPER > 1
454 } else if (impliedUpper < ub + 1.0e-7) {
455 #if 0
456 printf("%d neg coeff %g bounds %g,%g maxDown %g, maxUp %g rhs %g,%g\n",
457 icol,coeff,lb,ub,maxDown,maxUp,rlo[i],rup[i]);
458 printf("row is %g <= ",rlo[i]);
459 for (CoinBigIndex k = krs ; k < kre ; k++) {
460 const int icol2 = hcol[k] ;
461 printf("(%d,%g - bounds %g,%g) ",icol2,rowels[k],
462 clo[icol2],cup[icol2]);
463 }
464 printf("<= %g - implied l,u %g,%g\n",rup[i],impliedLower,impliedUpper);
465 printf("can take off ub as implied upper %g < %g\n",impliedUpper,ub);
466 #endif
467 #if TRY_UPPER > 2
468 doneSomething = true;
469 active[icol] = 1 + 8;
470 nFreed++;
471 #endif
472 #endif
473 }
474 } else if (!positive && false) {
475 printf("all negative %g <= ", rlo[i]);
476 for (CoinBigIndex k = krs; k < kre; k++) {
477 const double coeff = rowels[k];
478 const int icol = hcol[k];
479 printf("(%d,%g) ", icol, coeff);
480 }
481 printf("<= %g\n", rup[i]);
482 } else if (!negative && false) {
483 printf("all positive %g <= ", rlo[i]);
484 for (CoinBigIndex k = krs; k < kre; k++) {
485 const double coeff = rowels[k];
486 const int icol = hcol[k];
487 printf("(%d,%g) ", icol, coeff);
488 }
489 printf("<= %g\n", rup[i]);
490 }
491 }
492 //printf("%d had only one bound - %d added\n",nOneBound,nFreed);
493 /*
494 Initialise row dual min/max. The defaults are +/- infty, but if we know
495 that the logical has no upper (lower) bound, it can only be nonbasic at
496 the other bound. Given minimisation, we can conclude:
497 * <= constraint ==> [0,infty] ==> NBLB ==> cbar<i> > 0 ==> y<i> < 0
498 * >= constraint ==> [-infty,0] ==> NBUB ==> cbar<i> < 0 ==> y<i> > 0
499 Range constraints are not helpful here, the dual can be either sign. There's
500 no calculation because we assume the objective coefficient c<i> = 0 and the
501 coefficient a<ii> = 1.0, hence cbar<i> = -y<i>.
502 */
503 for (int i = 0; i < nrows; i++) {
504 const bool no_lb = (rup[i] >= ekkinf);
505 const bool no_ub = (rlo[i] <= -ekkinf);
506
507 ymin[i] = ((no_lb && !no_ub) ? 0.0 : (-PRESOLVE_INF));
508 ymax[i] = ((no_ub && !no_lb) ? 0.0 : PRESOLVE_INF);
509 }
510 /*
511 We can do a similar calculation with singleton columns where the variable
512 has only one finite bound, but we have to work a bit harder, as the cost
513 coefficient c<j> is not necessarily 0.0 and a<ij> is not necessarily 1.0.
514
515 cbar<j> = c<j> - y<i>a<ij>, hence y<i> = c<j>/a<ij> at cbar<j> = 0. The
516 question is whether this is an upper or lower bound on y<i>.
517 * x<j> NBLB ==> cbar<j> >= 0
518 a<ij> > 0 ==> increasing y<i> decreases cbar<j> ==> upper bound
519 a<ij> < 0 ==> increasing y<i> increases cbar<j> ==> lower bound
520 * x<j> NBUB ==> cbar<j> <= 0
521 a<ij> > 0 ==> increasing y<i> decreases cbar<j> ==> lower bound
522 a<ij> < 0 ==> increasing y<i> increases cbar<j> ==> upper bound
523 The condition below (simple test for equality to choose the bound) looks
524 a bit odd, but a bit of boolean algebra should convince you it's correct.
525 We have a bound; the only question is whether it's upper or lower.
526
527 Skip integer variables; it's far too likely that we'll tighten infinite
528 bounds elsewhere in presolve.
529
530 NOTE: If bound propagation is applied to continuous variables, the same
531 hazard will apply. -- lh, 110611 --
532 */
533 for (int j = 0; j < ncols; j++) {
534 if (integerType[j])
535 continue;
536 if (hincol[j] != 1)
537 continue;
538 #ifndef USE_ACTIVE
539 const bool no_ub = (cup[j] >= ekkinf);
540 const bool no_lb = (clo[j] <= -ekkinf);
541 #else
542 const bool no_ub = (active[j] & 6) == 0;
543 const bool no_lb = (active[j] & 5) == 0;
544 #endif
545 if (no_ub != no_lb) {
546 const int &i = hrow[mcstrt[j]];
547 double aij = colels[mcstrt[j]];
548 PRESOLVEASSERT(fabs(aij) > ZTOLDP);
549 const double yzero = cost[j] / aij;
550 if ((aij > 0.0) == no_ub) {
551 if (ymax[i] > yzero)
552 ymax[i] = yzero;
553 } else {
554 if (ymin[i] < yzero)
555 ymin[i] = yzero;
556 }
557 }
558 }
559 int nfixup_cols = 0;
560 int nfixdown_cols = ncols;
561 // Grab another work array, sized to ncols_
562 int *fix_cols = prob->usefulColumnInt_;
563 #if PRESOLVE_TIGHTEN_DUALS > 0
564 double *cbarmin = new double[ncols];
565 double *cbarmax = new double[ncols];
566 #endif
567 /*
568 Now we have (admittedly weak) bounds on the dual for each row. We can use
569 these to calculate upper and lower bounds on cbar<j>. Open loops to
570 take multiple passes over all columns.
571 */
572 int nPass = 0;
573 while (nPass++ < 100) {
574 int tightened = 0;
575 for (int j = 0; j < ncols; j++) {
576 if (hincol[j] <= 0)
577 continue;
578 /*
579 Calculate min cbar<j> and max cbar<j> for the column by calculating the
580 contribution to c<j>-ya<j> using ymin<i> and ymax<i> from above.
581 */
582 const CoinBigIndex &kcs = mcstrt[j];
583 const CoinBigIndex kce = kcs + hincol[j];
584 // Number of infinite rows
585 int nflagu = 0;
586 int nflagl = 0;
587 // Number of ordinary rows
588 int nordu = 0;
589 int nordl = 0;
590 double cbarjmin = cost[j];
591 double cbarjmax = cbarjmin;
592 for (CoinBigIndex k = kcs; k < kce; k++) {
593 const int &i = hrow[k];
594 const double &aij = colels[k];
595 const double mindelta = aij * ymin[i];
596 const double maxdelta = aij * ymax[i];
597
598 if (aij > 0.0) {
599 if (ymin[i] >= -ekkinf2) {
600 cbarjmax -= mindelta;
601 nordu++;
602 } else {
603 nflagu++;
604 }
605 if (ymax[i] <= ekkinf2) {
606 cbarjmin -= maxdelta;
607 nordl++;
608 } else {
609 nflagl++;
610 }
611 } else {
612 if (ymax[i] <= ekkinf2) {
613 cbarjmax -= maxdelta;
614 nordu++;
615 } else {
616 nflagu++;
617 }
618 if (ymin[i] >= -ekkinf2) {
619 cbarjmin -= mindelta;
620 nordl++;
621 } else {
622 nflagl++;
623 }
624 }
625 }
626 /*
627 See if we can tighten bounds on a dual y<t>. See the comments at the head of
628 the file for the linear algebra.
629
630 At this point, I don't understand the restrictions on propagation. Neither
631 are necessary. The net effect is to severely restrict the circumstances
632 where we'll propagate. Requiring nflagu == 1 excludes the case where all
633 duals have finite bounds (unlikely?). And why cbarjmax < -ztolcbarj?
634
635 In any event, we're looking for the row that's contributing the infinity. If
636 a<tj> > 0, the contribution is due to ymin<t> and we tighten ymax<t>;
637 a<tj> < 0, ymax<t> and ymin<t>, respectively.
638
639 The requirement cbarjmax < (ymax[i]*aij-ztolcbarj) ensures we don't propagate
640 tiny changes. If we make a change, it will affect cbarjmin. Make the
641 adjustment immediately.
642
643 Continuous variables only.
644 */
645 if (!integerType[j]) {
646 #ifndef USE_ACTIVE
647 const bool no_cub = (cup[j] >= ekkinf);
648 //const bool no_clb = (clo[j] <= -ekkinf) ;
649 #else
650 const bool no_cub = (active[j] & 6) == 0;
651 //const bool no_clb = (active[j]&5)==0 ;
652 #endif
653 if (no_cub) {
654 if (nflagu == 1 && cbarjmax < -ztolcbarj) {
655 for (CoinBigIndex k = kcs; k < kce; k++) {
656 const int i = hrow[k];
657 const double aij = colels[k];
658 if (aij > 0.0 && ymin[i] < -ekkinf2) {
659 if (cbarjmax < (ymax[i] * aij - ztolcbarj)) {
660 const double newValue = cbarjmax / aij;
661 if (ymax[i] > ekkinf2 && newValue <= ekkinf2) {
662 nflagl--;
663 cbarjmin -= aij * newValue;
664 } else if (ymax[i] <= ekkinf2) {
665 cbarjmin -= aij * (newValue - ymax[i]);
666 }
667 #if PRESOLVE_DEBUG > 1
668 std::cout
669 << "NDUAL(infu/inf): u(" << j << ") = " << cup[j]
670 << ": max y(" << i << ") was " << ymax[i]
671 << " now " << newValue << "." << std::endl;
672 #endif
673 ymax[i] = newValue;
674 tightened++;
675 }
676 } else if (aij < 0.0 && ymax[i] > ekkinf2) {
677 if (cbarjmax < (ymin[i] * aij - ztolcbarj)) {
678 const double newValue = cbarjmax / aij;
679 if (ymin[i] < -ekkinf2 && newValue >= -ekkinf2) {
680 nflagl--;
681 cbarjmin -= aij * newValue;
682 } else if (ymin[i] >= -ekkinf2) {
683 cbarjmin -= aij * (newValue - ymin[i]);
684 }
685 #if PRESOLVE_DEBUG > 1
686 std::cout
687 << "NDUAL(infu/inf): u(" << j << ") = " << cup[j]
688 << ": min y(" << i << ") was " << ymin[i]
689 << " now " << newValue << "." << std::endl;
690 #endif
691 ymin[i] = newValue;
692 tightened++;
693 // Huh? asymmetric
694 // cbarjmin = 0.0 ;
695 }
696 }
697 }
698 } else if (nflagl == 0 && nordl == 1 && cbarjmin < -ztolcbarj) {
699 /*
700 This is a column singleton. Why are we doing this? It's not like changes
701 to other y will affect this.
702 */
703 /* In fact a bug as ymax/min changed - but as I don't
704 wish to totally understand coding will let it abort
705 if not a singleton */
706 if (kce != kcs + 1)
707 abort();
708 #if 0
709 for (CoinBigIndex k = kcs; k < kce; k++) {
710 const int i = hrow[k] ;
711 const double aij = colels[k] ;
712 if (aij > 0.0) {
713 #if PRESOLVE_DEBUG > 1
714 std::cout
715 << "NDUAL(infu/sing): u(" << j << ") = " << cup[j]
716 << ": max y(" << i << ") was " << ymax[i]
717 << " now " << ymax[i]+cbarjmin/aij << "." << std::endl ;
718 #endif
719 ymax[i] += cbarjmin/aij ;
720 cbarjmin = 0.0 ;
721 tightened++ ;
722 } else if (aij < 0.0 ) {
723 #if PRESOLVE_DEBUG > 1
724 std::cout
725 << "NDUAL(infu/sing): u(" << j << ") = " << cup[j]
726 << ": min y(" << i << ") was " << ymin[i]
727 << " now " << ymin[i]+cbarjmin/aij << "." << std::endl ;
728 #endif
729 ymin[i] += cbarjmin/aij ;
730 cbarjmin = 0.0 ;
731 tightened++ ;
732 }
733 }
734 #endif
735 }
736 } // end u<j> = infty
737 #if PROCESS_INFINITE_LB
738 /*
739 Unclear why this section is commented out, except for the possibility that
740 bounds of -infty < x < something are rare and it likely suffered fromm the
741 same errors. Consider the likelihood that this whole block needs to be
742 edited to sway min/max, & similar.
743 */
744 if (no_clb) {
745 // cbarj can not be positive
746 if (cbarjmin > ztolcbarj && nflagl == 1) {
747 // We can make bound finite one way
748 for (CoinBigIndex k = kcs; k < kce; k++) {
749 const int i = hrow[k];
750 const double coeff = colels[k];
751
752 if (coeff < 0.0 && ymin[i] < -ekkinf2) {
753 // ymax[i] has upper bound
754 if (cbarjmin > ymax[i] * coeff + ztolcbarj) {
755 const double newValue = cbarjmin / coeff;
756 // re-compute hi
757 if (ymax[i] > ekkinf2 && newValue <= ekkinf2) {
758 nflagu--;
759 cbarjmax -= coeff * newValue;
760 } else if (ymax[i] <= ekkinf2) {
761 cbarjmax -= coeff * (newValue - ymax[i]);
762 }
763 ymax[i] = newValue;
764 tightened++;
765 #if PRESOLVE_DEBUG > 1
766 printf("Col %d, row %d max pi now %g\n", j, i, ymax[i]);
767 #endif
768 }
769 } else if (coeff > 0.0 && ymax[i] > ekkinf2) {
770 // ymin[i] has lower bound
771 if (cbarjmin > ymin[i] * coeff + ztolcbarj) {
772 const double newValue = cbarjmin / coeff;
773 // re-compute lo
774 if (ymin[i] < -ekkinf2 && newValue >= -ekkinf2) {
775 nflagu--;
776 cbarjmax -= coeff * newValue;
777 } else if (ymin[i] >= -ekkinf2) {
778 cbarjmax -= coeff * (newValue - ymin[i]);
779 }
780 ymin[i] = newValue;
781 tightened++;
782 #if PRESOLVE_DEBUG > 1
783 printf("Col %d, row %d min pi now %g\n", j, i, ymin[i]);
784 #endif
785 }
786 }
787 }
788 } else if (nflagu == 0 && nordu == 1 && cbarjmax > ztolcbarj) {
789 // We may be able to tighten
790 for (CoinBigIndex k = kcs; k < kce; k++) {
791 const int i = hrow[k];
792 const double coeff = colels[k];
793
794 if (coeff < 0.0) {
795 ymax[i] += cbarjmax / coeff;
796 cbarjmax = 0.0;
797 tightened++;
798 #if PRESOLVE_DEBUG > 1
799 printf("Col %d, row %d max pi now %g\n", j, i, ymax[i]);
800 #endif
801 } else if (coeff > 0.0) {
802 ymin[i] += cbarjmax / coeff;
803 cbarjmax = 0.0;
804 tightened++;
805 #if PRESOLVE_DEBUG > 1
806 printf("Col %d, row %d min pi now %g\n", j, i, ymin[i]);
807 #endif
808 }
809 }
810 }
811 } // end l<j> < -infty
812 #endif // PROCESS_INFINITE_LB
813 }
814 /*
815 That's the end of propagation of bounds for dual variables.
816 */
817 #if PRESOLVE_TIGHTEN_DUALS > 0
818 cbarmin[j] = (nflagl ? (-PRESOLVE_INF) : cbarjmin);
819 cbarmax[j] = (nflagu ? PRESOLVE_INF : cbarjmax);
820 #endif
821 /*
822 If cbarmin<j> > 0 (strict inequality) then x<j> NBLB at optimality. If l<j>
823 is -infinity, notify the user and set the status to unbounded.
824 */
825 if (cbarjmin > ztolcbarj && nflagl == 0 && !prob->colProhibited2(j)) {
826 if (clo[j] <= -ekkinf) {
827 CoinMessageHandler *msghdlr = prob->messageHandler();
828 msghdlr->message(COIN_PRESOLVE_COLUMNBOUNDB, prob->messages())
829 << j << CoinMessageEol;
830 prob->status_ |= 2;
831 break;
832 } else {
833 fix_cols[--nfixdown_cols] = j;
834 #if PRESOLVE_DEBUG > 1
835 std::cout << "NDUAL(fix l): fix x(" << j << ")";
836 if (csol)
837 std::cout << " = " << csol[j];
838 std::cout
839 << " at l(" << j << ") = " << clo[j] << "; cbar("
840 << j << ") > " << cbarjmin << "." << std::endl;
841 #endif
842 if (csol) {
843 #if REMOVE_DUAL_ACTION_REPAIR_SOLN > 0
844 /*
845 Original comment: User may have given us feasible solution - move if simple
846
847 Except it's not simple. The net result is that we end up with an excess of
848 basic variables. Mark x<j> NBLB and let the client recalculate the solution
849 after establishing the basis.
850 */
851 if (csol[j] - clo[j] > 1.0e-7 && hincol[j] == 1) {
852 double value_j = colels[mcstrt[j]];
853 double distance_j = csol[j] - clo[j];
854 int row = hrow[mcstrt[j]];
855 // See if another column can take value
856 for (CoinBigIndex kk = mrstrt[row]; kk < mrstrt[row] + hinrow[row]; kk++) {
857 const int k = hcol[kk];
858 if (colstat[k] == CoinPrePostsolveMatrix::superBasic)
859 continue;
860
861 if (hincol[k] == 1 && k != j) {
862 const double value_k = rowels[kk];
863 double movement;
864 if (value_k * value_j > 0.0) {
865 // k needs to increase
866 double distance_k = cup[k] - csol[k];
867 movement = CoinMin((distance_j * value_j) / value_k, distance_k);
868 } else {
869 // k needs to decrease
870 double distance_k = clo[k] - csol[k];
871 movement = CoinMax((distance_j * value_j) / value_k, distance_k);
872 }
873 if (relEq(movement, 0))
874 continue;
875
876 csol[k] += movement;
877 if (relEq(csol[k], clo[k])) {
878 colstat[k] = CoinPrePostsolveMatrix::atLowerBound;
879 } else if (relEq(csol[k], cup[k])) {
880 colstat[k] = CoinPrePostsolveMatrix::atUpperBound;
881 } else if (colstat[k] != CoinPrePostsolveMatrix::isFree) {
882 colstat[k] = CoinPrePostsolveMatrix::basic;
883 }
884 printf("NDUAL: x<%d> moved %g to %g; ",
885 k, movement, csol[k]);
886 printf("lb = %g, ub = %g, status now %s.\n",
887 clo[k], cup[k], columnStatusString(k));
888 distance_j -= (movement * value_k) / value_j;
889 csol[j] -= (movement * value_k) / value_j;
890 if (distance_j < 1.0e-7)
891 break;
892 }
893 }
894 }
895 #endif // repair solution.
896
897 csol[j] = clo[j];
898 colstat[j] = CoinPrePostsolveMatrix::atLowerBound;
899 }
900 }
901 } else if (cbarjmax < -ztolcbarj && nflagu == 0 && !prob->colProhibited2(j)) {
902 /*
903 If cbarmax<j> < 0 (strict inequality) then x<j> NBUB at optimality. If u<j>
904 is infinity, notify the user and set the status to unbounded.
905 */
906 if (cup[j] >= ekkinf) {
907 CoinMessageHandler *msghdlr = prob->messageHandler();
908 msghdlr->message(COIN_PRESOLVE_COLUMNBOUNDA, prob->messages())
909 << j << CoinMessageEol;
910 prob->status_ |= 2;
911 break;
912 } else {
913 fix_cols[nfixup_cols++] = j;
914 #if PRESOLVE_DEBUG > 1
915 std::cout << "NDUAL(fix u): fix x(" << j << ")";
916 if (csol)
917 std::cout << " = " << csol[j];
918 std::cout
919 << " at u(" << j << ") = " << cup[j] << "; cbar("
920 << j << ") < " << cbarjmax << "." << std::endl;
921 #endif
922 if (csol) {
923 #if 0
924 // See comments above for 'fix at lb'.
925 if (cup[j]-csol[j] > 1.0e-7 && hincol[j] == 1) {
926 double value_j = colels[mcstrt[j]] ;
927 double distance_j = csol[j]-cup[j] ;
928 int row = hrow[mcstrt[j]] ;
929 // See if another column can take value
930 for (CoinBigIndex kk = mrstrt[row] ; kk < mrstrt[row]+hinrow[row] ; kk++) {
931 const int k = hcol[kk] ;
932 if (colstat[k] == CoinPrePostsolveMatrix::superBasic)
933 continue ;
934
935 if (hincol[k] == 1 && k != j) {
936 const double value_k = rowels[kk] ;
937 double movement ;
938 if (value_k*value_j<0.0) {
939 // k needs to increase
940 double distance_k = cup[k]-csol[k] ;
941 movement = CoinMin((distance_j*value_j)/value_k,distance_k) ;
942 } else {
943 // k needs to decrease
944 double distance_k = clo[k]-csol[k] ;
945 movement = CoinMax((distance_j*value_j)/value_k,distance_k) ;
946 }
947 if (relEq(movement,0)) continue ;
948
949 csol[k] += movement ;
950 if (relEq(csol[k],clo[k]))
951 { colstat[k] = CoinPrePostsolveMatrix::atLowerBound ; }
952 else
953 if (relEq(csol[k],cup[k]))
954 { colstat[k] = CoinPrePostsolveMatrix::atUpperBound ; }
955 else
956 if (colstat[k] != CoinPrePostsolveMatrix::isFree)
957 { colstat[k] = CoinPrePostsolveMatrix::basic ; }
958 printf("NDUAL: x<%d> moved %g to %g; ",
959 k,movement,csol[k]) ;
960 printf("lb = %g, ub = %g, status now %s.\n",
961 clo[k],cup[k],columnStatusString(k)) ;
962 distance_j -= (movement*value_k)/value_j ;
963 csol[j] -= (movement*value_k)/value_j ;
964 if (distance_j>-1.0e-7)
965 break ;
966 }
967 }
968 }
969 #endif
970 csol[j] = cup[j];
971 colstat[j] = CoinPrePostsolveMatrix::atUpperBound;
972 }
973 }
974 } // end cbar<j> < 0
975 }
976 /*
977 That's the end of this walk through the columns.
978 */
979 // I don't know why I stopped doing this.
980 #if PRESOLVE_TIGHTEN_DUALS > 0
981 const double *rowels = prob->rowels_;
982 const int *hcol = prob->hcol_;
983 const CoinBigIndex *mrstrt = prob->mrstrt_;
984 int *hinrow = prob->hinrow_;
985 // tighten row dual bounds, as described on p. 229
986 for (int i = 0; i < nrows; i++) {
987 const bool no_ub = (rup[i] >= ekkinf);
988 const bool no_lb = (rlo[i] <= -ekkinf);
989
990 if ((no_ub ^ no_lb) == true) {
991 const CoinBigIndex krs = mrstrt[i];
992 const CoinBigIndex kre = krs + hinrow[i];
993 const double rmax = ymax[i];
994 const double rmin = ymin[i];
995
996 // all row columns are non-empty
997 for (CoinBigIndex k = krs; k < kre; k++) {
998 const double coeff = rowels[k];
999 const int icol = hcol[k];
1000 const double cbarmax0 = cbarmax[icol];
1001 const double cbarmin0 = cbarmin[icol];
1002
1003 if (no_ub) {
1004 // cbarj must not be negative
1005 if (coeff > ZTOLDP2 && cbarjmax0 < PRESOLVE_INF && cup[icol] >= ekkinf) {
1006 const double bnd = cbarjmax0 / coeff;
1007 if (rmax > bnd) {
1008 #if PRESOLVE_DEBUG > 1
1009 printf("MAX TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], ymax[i], bnd);
1010 #endif
1011 ymax[i] = rmax = bnd;
1012 tightened++;
1013 ;
1014 }
1015 } else if (coeff < -ZTOLDP2 && cbarjmax0 < PRESOLVE_INF && cup[icol] >= ekkinf) {
1016 const double bnd = cbarjmax0 / coeff;
1017 if (rmin < bnd) {
1018 #if PRESOLVE_DEBUG > 1
1019 printf("MIN TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], ymin[i], bnd);
1020 #endif
1021 ymin[i] = rmin = bnd;
1022 tightened++;
1023 ;
1024 }
1025 }
1026 } else { // no_lb
1027 // cbarj must not be positive
1028 if (coeff > ZTOLDP2 && cbarmin0 > -PRESOLVE_INF && clo[icol] <= -ekkinf) {
1029 const double bnd = cbarmin0 / coeff;
1030 if (rmin < bnd) {
1031 #if PRESOLVE_DEBUG > 1
1032 printf("MIN1 TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], ymin[i], bnd);
1033 #endif
1034 ymin[i] = rmin = bnd;
1035 tightened++;
1036 ;
1037 }
1038 } else if (coeff < -ZTOLDP2 && cbarmin0 > -PRESOLVE_INF && clo[icol] <= -ekkinf) {
1039 const double bnd = cbarmin0 / coeff;
1040 if (rmax > bnd) {
1041 #if PRESOLVE_DEBUG > 1
1042 printf("MAX TIGHT1[%d,%d]: %g --> %g\n", i, hrow[k], ymax[i], bnd);
1043 #endif
1044 ymax[i] = rmax = bnd;
1045 tightened++;
1046 ;
1047 }
1048 }
1049 }
1050 }
1051 }
1052 }
1053 #endif // PRESOLVE_TIGHTEN_DUALS
1054 /*
1055 Is it productive to continue with another pass? Essentially, we need lots of
1056 tightening but no fixing. If we fixed any variables, break and process them.
1057 */
1058 #if PRESOLVE_DEBUG > 1
1059 std::cout
1060 << "NDUAL: pass " << nPass << ": fixed " << (ncols - nfixdown_cols)
1061 << " down, " << nfixup_cols << " up, tightened " << tightened
1062 << " duals." << std::endl;
1063 #endif
1064 if (tightened < 100 || nfixdown_cols < ncols || nfixup_cols)
1065 break;
1066 }
1067 assert(nfixup_cols <= nfixdown_cols);
1068 /*
1069 Process columns fixed at upper bound.
1070 */
1071 if (nfixup_cols) {
1072 #if PRESOLVE_DEBUG > 1
1073 std::cout << "NDUAL(upper):";
1074 for (int k = 0; k < nfixup_cols; k++)
1075 std::cout << " " << fix_cols[k];
1076 std::cout << "." << std::endl;
1077 #endif
1078 next = make_fixed_action::presolve(prob, fix_cols, nfixup_cols, false, next);
1079 }
1080 /*
1081 Process columns fixed at lower bound.
1082 */
1083 if (nfixdown_cols < ncols) {
1084 int *fixdown_cols = fix_cols + nfixdown_cols;
1085 nfixdown_cols = ncols - nfixdown_cols;
1086 #if PRESOLVE_DEBUG > 1
1087 std::cout << "NDUAL(lower):";
1088 for (int k = 0; k < nfixdown_cols; k++)
1089 std::cout << " " << fixdown_cols[k];
1090 std::cout << "." << std::endl;
1091 #endif
1092 next = make_fixed_action::presolve(prob, fixdown_cols,
1093 nfixdown_cols, true, next);
1094 }
1095 /*
1096 Now look for variables that, when moved in the favourable direction
1097 according to reduced cost, will naturally tighten an inequality to an
1098 equality. We can convert that inequality to an equality. See the comments
1099 at the head of the routine.
1100
1101 Start with logicals. Open a loop and look for suitable rows. At the
1102 end of this loop, rows marked with +/-1 will be forced to equality by
1103 their logical. Rows marked with +/-2 are inequalities that (perhaps)
1104 can be forced to equality using architecturals. Rows marked with 0 are
1105 not suitable (range or nonbinding).
1106
1107 Note that usefulRowInt_ is 3*nrows_; we'll use the second partition below.
1108 */
1109 int *canFix = prob->usefulRowInt_;
1110 for (int i = 0; i < nrows; i++) {
1111 const bool no_rlb = (rlo[i] <= -ekkinf);
1112 const bool no_rub = (rup[i] >= ekkinf);
1113 canFix[i] = 0;
1114 if (no_rub && !no_rlb) {
1115 if (ymin[i] > 0.0)
1116 canFix[i] = -1;
1117 else
1118 canFix[i] = -2;
1119 } else if (no_rlb && !no_rub) {
1120 if (ymax[i] < 0.0)
1121 canFix[i] = 1;
1122 else
1123 canFix[i] = 2;
1124 }
1125 #if PRESOLVE_DEBUG > 1
1126 if (abs(canFix[i]) == 1) {
1127 std::cout
1128 << "NDUAL(eq): candidate row (" << i << ") (" << rlo[i]
1129 << "," << rup[i] << "), logical must be "
1130 << ((canFix[i] == -1) ? "NBUB" : "NBLB") << "." << std::endl;
1131 }
1132 #endif
1133 }
1134 /*
1135 Can we do a similar trick with architectural variables? Here, we're looking
1136 for x<j> such that
1137 (1) Exactly one of l<j> or u<j> is infinite.
1138 (2) Moving x<j> towards the infinite bound can tighten exactly one
1139 constraint i to equality. If x<j> is entangled with other constraints,
1140 moving x<j> towards the infinite bound will loosen those constraints.
1141 (3) Moving x<j> towards the infinite bound is a good idea according to the
1142 cost c<j> (note we don't have to consider reduced cost here).
1143 If we can find a suitable x<j>, constraint i can become an equality.
1144
1145 This is yet another instance of bound propagation, but we're looking for
1146 a very specific pattern: A variable that can be increased arbitrarily
1147 in all rows it's entangled with, except for one, which bounds it. And
1148 we're going to push the variable so as to make that row an equality.
1149 But note what we're *not* doing: No actual comparison of finite bound
1150 values to the amount necessary to force an equality. So no worries about
1151 accuracy, the bane of bound propagation.
1152
1153 Open a loop to scan columns. bindingUp and bindingDown indicate the result
1154 of the analysis; -1 says `possible', -2 is ruled out. Test first for
1155 condition (1). Column singletons are presumably handled elsewhere. Integer
1156 variables need not apply. If both bounds are finite, no need to look
1157 further.
1158 */
1159 for (int j = 0; j < ncols; j++) {
1160 if (hincol[j] <= 1)
1161 continue;
1162 if (integerType[j])
1163 continue;
1164 int bindingUp = -1;
1165 int bindingDown = -1;
1166 if (cup[j] < ekkinf)
1167 bindingUp = -2;
1168 if (clo[j] > -ekkinf)
1169 bindingDown = -2;
1170 if (bindingUp == -2 && bindingDown == -2)
1171 continue;
1172 /*
1173 Open a loop to walk the column and check for condition (2).
1174
1175 The test for |canFix[i]| != 2 is a non-interference check. We don't want to
1176 mess with constraints where we've already decided to use the logical to
1177 force equality. Nor do we want to deal with range or nonbinding constraints.
1178 */
1179 const CoinBigIndex &kcs = mcstrt[j];
1180 const CoinBigIndex kce = kcs + hincol[j];
1181 for (CoinBigIndex k = kcs; k < kce; k++) {
1182 const int &i = hrow[k];
1183 if (abs(canFix[i]) != 2) {
1184 bindingUp = -2;
1185 bindingDown = -2;
1186 break;
1187 }
1188 double aij = colels[k];
1189 /*
1190 For a<ij> > 0 in a <= constraint (canFix = 2), the up direction is
1191 binding. For a >= constraint, it'll be the down direction. If the relevant
1192 binding code is still -1, set it to the index of the row. Similarly for
1193 a<ij> < 0.
1194
1195 If this is the second or subsequent binding constraint in that direction,
1196 set binding[Up,Down] to -2 (we don't want to get into the business of
1197 calculating which constraint is actually binding).
1198 */
1199 if (aij > 0.0) {
1200 if (canFix[i] == 2) {
1201 if (bindingUp == -1)
1202 bindingUp = i;
1203 else
1204 bindingUp = -2;
1205 } else {
1206 if (bindingDown == -1)
1207 bindingDown = i;
1208 else
1209 bindingDown = -2;
1210 }
1211 } else {
1212 if (canFix[i] == 2) {
1213 if (bindingDown == -1)
1214 bindingDown = i;
1215 else
1216 bindingDown = -2;
1217 } else {
1218 if (bindingUp == -1)
1219 bindingUp = i;
1220 else
1221 bindingUp = -2;
1222 }
1223 }
1224 }
1225 if (bindingUp == -2 && bindingDown == -2)
1226 continue;
1227 /*
1228 If bindingUp > -2, then either no constraint provided a bound (-1) or
1229 there's a single constraint (0 <= i < m) that bounds x<j>. If we have
1230 just one binding constraint, check that the reduced cost is favourable
1231 (c<j> <= 0 for x<j> NBUB at optimum for minimisation). If so, declare
1232 that we will force the row to equality (canFix[i] = +/-1). Note that we
1233 don't adjust the primal solution value for x<j>.
1234
1235 If no constraint provided a bound, we might be headed for unboundedness,
1236 but leave that for some other code to determine.
1237 */
1238 double cj = cost[j];
1239 if (bindingUp > -2 && cj <= 0.0) {
1240 if (bindingUp >= 0) {
1241 canFix[bindingUp] /= 2;
1242 #if PRESOLVE_DEBUG > 1
1243 std::cout
1244 << "NDUAL(eq): candidate row (" << bindingUp << ") ("
1245 << rlo[bindingUp] << "," << rup[bindingUp] << "), "
1246 << " increasing x(" << j << "), cbar = " << cj << "."
1247 << std::endl;
1248 } else {
1249 std::cout
1250 << "NDUAL(eq): no binding upper bound for x(" << j
1251 << "), cbar = " << cj << "." << std::endl;
1252 #endif
1253 }
1254 } else if (bindingDown > -2 && cj >= 0.0) {
1255 if (bindingDown >= 0) {
1256 canFix[bindingDown] /= 2;
1257 #if PRESOLVE_DEBUG > 1
1258 std::cout
1259 << "NDUAL(eq): candidate row (" << bindingDown << ") ("
1260 << rlo[bindingDown] << "," << rup[bindingDown] << "), "
1261 << " decreasing x(" << j << "), cbar = " << cj << "."
1262 << std::endl;
1263 } else {
1264 std::cout
1265 << "NDUAL(eq): no binding lower bound for x(" << j
1266 << "), cbar = " << cj << "." << std::endl;
1267 #endif
1268 }
1269 }
1270 }
1271 /*
1272 We have candidate rows. We've avoided scanning full rows until now,
1273 but there's one remaining hazard: if the row contains unfixed integer
1274 variables then we don't want to just pin the row to a fixed rhs; that
1275 might prevent us from achieving integrality. Scan canFix, count and
1276 record suitable rows (use the second partition of usefulRowInt_).
1277 */
1278 // get min max stuff - we can re-use ymin/ymax
1279 int *infCount = prob->usefulRowInt_ + 2 * nrows;
1280 {
1281 // copied from CglProbing
1282 int i, j;
1283 CoinBigIndex k, krs, kre;
1284 int iflagu, iflagl;
1285 double dmaxup, dmaxdown;
1286
1287 for (i = 0; i < nrows; ++i) {
1288 if (rlo[i] == rup[i]) {
1289 infCount[i] = 10 | (10 << 16);
1290 } else if (rlo[i] > -1.0e20 || rup[i] < 1.0e20) {
1291 iflagu = 0;
1292 iflagl = 0;
1293 dmaxup = 0.0;
1294 dmaxdown = 0.0;
1295 krs = mrstrt[i];
1296 kre = mrstrt[i] + hinrow[i];
1297
1298 /* ------------------------------------------------------------*/
1299 /* Compute L(i) and U(i) */
1300 /* ------------------------------------------------------------*/
1301 for (k = krs; k < kre; ++k) {
1302 double value = rowels[k];
1303 j = hcol[k];
1304 if (value > 0.0) {
1305 if (cup[j] < 1.0e12)
1306 dmaxup += cup[j] * value;
1307 else
1308 ++iflagu;
1309 if (clo[j] > -1.0e12)
1310 dmaxdown += clo[j] * value;
1311 else
1312 ++iflagl;
1313 } else if (value < 0.0) {
1314 if (cup[j] < 1.0e12)
1315 dmaxdown += cup[j] * value;
1316 else
1317 ++iflagl;
1318 if (clo[j] > -1.0e12)
1319 dmaxup += clo[j] * value;
1320 else
1321 ++iflagu;
1322 }
1323 }
1324 iflagl = CoinMin(iflagl, 2);
1325 iflagu = CoinMin(iflagu, 2);
1326 infCount[i] = iflagl | (iflagu << 16);
1327 ymax[i] = dmaxup;
1328 ymin[i] = dmaxdown;
1329 } else {
1330 infCount[i] = 10 | (10 << 16);
1331 }
1332 }
1333 }
1334 for (int j = 0; j < ncols; j++) {
1335 if (hincol[j] == 1 && cost[j] > -1.0e10) {
1336 // can we make equality row
1337 double coeff = colels[mcstrt[j]];
1338 int irow = hrow[mcstrt[j]];
1339 int iflagl = infCount[irow] & 63;
1340 int iflagu = infCount[irow] >> 16;
1341 if (iflagu > 1 && iflagl > 1)
1342 continue;
1343 double dmaxdown = ymin[irow];
1344 double dmaxup = ymax[irow];
1345 #if USE_ACTIVE > 1
1346 printf("singleton col %d has %g on row %d - rhs %g,%g infs %d,%d maxes %g,%g\n",
1347 j, coeff, irow, rlo[irow], rup[irow],
1348 iflagl, iflagu,
1349 dmaxdown, dmaxup);
1350 #endif
1351 if (coeff > 0.0) {
1352 if (clo[j] > -1.0e12) {
1353 dmaxdown -= coeff * clo[j];
1354 if (iflagl)
1355 dmaxdown = -1.0e60;
1356 } else if (iflagl > 1) {
1357 dmaxdown = -1.0e60;
1358 }
1359 if (cup[j] < 1.0e12) {
1360 dmaxup -= coeff * cup[j];
1361 if (iflagu)
1362 dmaxup = 1.0e60;
1363 } else if (iflagu > 1) {
1364 dmaxup = 1.0e60;
1365 }
1366 } else {
1367 }
1368 double rhs;
1369 if (cost[j] > 0.0) {
1370 // we want as small as possible
1371 if (coeff > 0) {
1372 if (rlo[irow] > -1.0e12)
1373 rhs = rlo[irow] / coeff;
1374 else
1375 rhs = -COIN_DBL_MAX;
1376 } else {
1377 if (rup[irow] < 1.0e12)
1378 rhs = -rup[irow] / coeff;
1379 else
1380 rhs = COIN_DBL_MAX;
1381 }
1382 } else {
1383 // we want as large as possible
1384 if (coeff > 0) {
1385 if (rup[irow] < 1.0e12) {
1386 rhs = rup[irow];
1387 if (cup[j] * coeff + dmaxdown > rhs - 1.0e-7 && clo[j] * coeff + dmaxup < rhs + 1.0e-7) {
1388 #if 0 //USE_ACTIVE<2
1389 printf("singleton col %d has %g on row %d - rhs %g,%g infs %d,%d maxes %g,%g\n",
1390 j,coeff,irow,rlo[irow],rup[irow],
1391 iflagl,iflagu,
1392 dmaxdown,dmaxup);
1393 #endif
1394 //printf("making rup equality\n");
1395 canFix[irow] = 1;
1396 infCount[irow] = 10 | (10 << 16);
1397 }
1398 rhs /= coeff;
1399 } else {
1400 if (cup[j] * coeff + dmaxdown > 1.0e30) {
1401 //printf("unbounded?\n");
1402 }
1403 }
1404 } else {
1405 if (rlo[irow] > -1.0e12) {
1406 rhs = rlo[irow];
1407 if (cup[j] * coeff + dmaxdown < rhs + 1.0e-7
1408 && clo[j] * coeff + dmaxup > rhs - 1.0e-7) {
1409 #if USE_ACTIVE > 1
1410 printf("singleton col %d has %g on row %d - rhs %g,%g infs %d,%d maxes %g,%g\n",
1411 j, coeff, irow, rlo[irow], rup[irow],
1412 iflagl, iflagu,
1413 dmaxdown, dmaxup);
1414 printf("making rlo equality\n");
1415 #endif
1416 canFix[irow] = -11;
1417 infCount[irow] = 10 | (10 << 16);
1418 }
1419 rhs /= coeff;
1420 } else {
1421 if (cup[j] * coeff + dmaxdown > 1.0e30) {
1422 //printf("unbounded?\n");
1423 }
1424 }
1425 }
1426 }
1427 #if USE_ACTIVE > 1
1428 printf("adjusted maxes %g,%g - rhs %g\n", dmaxdown, dmaxup, rhs);
1429 #endif
1430 }
1431 }
1432 #if 1 // PRESOLVE_DEBUG > 0
1433 int makeEqCandCnt = 0;
1434 for (int i = 0; i < nrows; i++) {
1435 if (abs(canFix[i]) == 1)
1436 makeEqCandCnt++;
1437 }
1438 #endif
1439 int makeEqCnt = nrows;
1440 for (int i = 0; i < nrows; i++) {
1441 if (abs(canFix[i]) == 1 && hinrow[i] > 1) {
1442 const CoinBigIndex krs = mrstrt[i];
1443 const CoinBigIndex kre = krs + hinrow[i];
1444 int nBinary = 0;
1445 int nOther = 0;
1446 int nSameCoeff = 0;
1447 double rhs = canFix[i] == 1 ? rup[i] : rlo[i];
1448 double firstCoeff = 0.0;
1449 int j = -1;
1450 for (CoinBigIndex k = krs; k < kre; k++) {
1451 j = hcol[k];
1452 double coeff = rowels[k];
1453 if (cup[j] > clo[j]) {
1454 if (!firstCoeff) {
1455 firstCoeff = coeff;
1456 nSameCoeff = 1;
1457 } else {
1458 if (coeff == firstCoeff)
1459 nSameCoeff++;
1460 else if (coeff != -firstCoeff)
1461 nOther = 1;
1462 }
1463 if (prob->colProhibited2(j)) {
1464 nBinary = 1;
1465 nOther = 1;
1466 break;
1467 } else if (integerType[j]) {
1468 if (clo[j] == 0.0 && cup[j] == 1.0) {
1469 nBinary++;
1470 } else {
1471 nBinary = 1;
1472 nOther = 1;
1473 break;
1474 }
1475 } else {
1476 nOther++;
1477 }
1478 } else {
1479 rhs -= coeff * clo[j];
1480 }
1481 }
1482 bool canFixThis = true;
1483 if (nBinary) {
1484 if (nOther) {
1485 canFixThis = false;
1486 #if PRESOLVE_DEBUG > 1
1487 std::cout
1488 << "NDUAL(eq): cannot convert row " << i << " to equality; "
1489 << "unfixed integer variable x(" << j << ")." << std::endl;
1490 #endif
1491 } else {
1492 // may be able to if all easy integer
1493 if (!rhs) {
1494 // for now just two of opposite sign
1495 if (nSameCoeff == 1 && nBinary == 2) {
1496 // ok
1497 #if PRESOLVE_DEBUG > 1
1498 printf("int eq row is %g <= ", rlo[i]);
1499 for (CoinBigIndex k = krs; k < kre; k++) {
1500 const double coeff = rowels[k];
1501 const int icol = hcol[k];
1502 printf("(%d,%g - bounds %g,%g) ", icol, coeff, clo[icol], cup[icol]);
1503 }
1504 printf("<= %g\n", rup[i]);
1505 #endif
1506 } else {
1507 canFixThis = false;
1508 #if PRESOLVE_DEBUG > 1
1509 std::cout
1510 << "NDUAL(eq2): cannot convert row " << i << " to equality; "
1511 << "unfixed integer variable x(" << j << ")." << std::endl;
1512 #endif
1513 }
1514 } else if (rhs == 1.0 && canFix[i] == 1 && nSameCoeff == nBinary && firstCoeff == 1.0) {
1515 // ok
1516 #if PRESOLVE_DEBUG > 1
1517 printf("int2 eq row is %g <= ", rlo[i]);
1518 for (CoinBigIndex k = krs; k < kre; k++) {
1519 const double coeff = rowels[k];
1520 const int icol = hcol[k];
1521 printf("(%d,%g - bounds %g,%g) ", icol, coeff, clo[icol], cup[icol]);
1522 }
1523 printf("<= %g\n", rup[i]);
1524 #endif
1525 } else {
1526 canFixThis = false;
1527 #if PRESOLVE_DEBUG > 1
1528 std::cout
1529 << "NDUAL(eq3): cannot convert row " << i << " to equality; "
1530 << "unfixed integer variable x(" << j << ")." << std::endl;
1531 #endif
1532 }
1533 }
1534 }
1535 if (canFixThis)
1536 canFix[makeEqCnt++] = i;
1537 }
1538 }
1539 makeEqCnt -= nrows;
1540 #if PRESOLVE_DEBUG > 0
1541 if ((makeEqCandCnt - makeEqCnt) > 0) {
1542 std::cout
1543 << "NDUAL(eq): rejected " << (makeEqCandCnt - makeEqCnt)
1544 << " rows due to unfixed integer variables." << std::endl;
1545 }
1546 #endif
1547 /*
1548 If we've identified inequalities to convert, do the conversion, record
1549 the information needed to restore bounds in postsolve, and finally create
1550 the postsolve object.
1551 */
1552 if (makeEqCnt > 0) {
1553 action *bndRecords = new action[makeEqCnt];
1554 for (int k = 0; k < makeEqCnt; k++) {
1555 const int i = canFix[k + nrows];
1556 #if PRESOLVE_DEBUG > 1
1557 std::cout << "NDUAL(eq): forcing row " << i << " to equality;";
1558 if (canFix[i] == -1)
1559 std::cout << " dropping b = " << rup[i] << " to " << rlo[i];
1560 else
1561 std::cout << " raising blow = " << rlo[i] << " to " << rup[i];
1562 std::cout << "." << std::endl;
1563 #endif
1564 action &bndRec = bndRecords[(k)];
1565 bndRec.rlo_ = rlo[i];
1566 bndRec.rup_ = rup[i];
1567 bndRec.ndx_ = i;
1568 if (canFix[i] == 1) {
1569 rlo[i] = rup[i];
1570 prob->addRow(i);
1571 } else if (canFix[i] == -1) {
1572 rup[i] = rlo[i];
1573 prob->addRow(i);
1574 }
1575 }
1576 next = new remove_dual_action(makeEqCnt, bndRecords, next);
1577 }
1578
1579 #if PRESOLVE_TIGHTEN_DUALS > 0
1580 delete[] cbarmin;
1581 delete[] cbarmax;
1582 #endif
1583 #undef PRESOLVE_DEBUG
1584 #if COIN_PRESOLVE_TUNING > 0
1585 double thisTime = 0.0;
1586 if (prob->tuning_)
1587 thisTime = CoinCpuTime();
1588 #endif
1589 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
1590 presolve_check_sol(prob, 2, 1, 1);
1591 presolve_check_nbasic(prob);
1592 #endif
1593 #if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
1594 int droppedRows = prob->countEmptyRows() - startEmptyRows;
1595 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
1596 std::cout
1597 << "Leaving remove_dual_action::presolve, dropped " << droppedRows
1598 << " rows, " << droppedColumns << " columns, forced "
1599 << makeEqCnt << " equalities";
1600 #if COIN_PRESOLVE_TUNING > 0
1601 if (prob->tuning_)
1602 std::cout << " in " << (thisTime - prob->startTime_) << "s";
1603 #endif
1604 std::cout << "." << std::endl;
1605 #endif
1606
1607 return (next);
1608 }
1609
1610 /*
1611 Postsolve: replace the original row bounds.
1612
1613 The catch here is that each constraint was an equality in the presolved
1614 problem, with a logical s<i> that had l<i> = u<i> = 0. We're about to
1615 convert the equality back to an inequality. One row bound will go to
1616 infinity, as will one of the bounds of the logical. We may need to patch the
1617 basis. The logical for a <= constraint cannot be NBUB, and the logical for a
1618 >= constraint cannot be NBLB.
1619 */
postsolve(CoinPostsolveMatrix * prob) const1620 void remove_dual_action::postsolve(CoinPostsolveMatrix *prob) const
1621 {
1622 const action *const &bndRecords = actions_;
1623 const int &numRecs = nactions_;
1624
1625 double *&rlo = prob->rlo_;
1626 double *&rup = prob->rup_;
1627 unsigned char *&rowstat = prob->rowstat_;
1628
1629 #if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
1630 #if PRESOLVE_DEBUG > 0
1631 std::cout
1632 << "Entering remove_dual_action::postsolve, " << numRecs
1633 << " bounds to restore." << std::endl;
1634 #endif
1635 presolve_check_threads(prob);
1636 presolve_check_sol(prob, 2, 2, 2);
1637 presolve_check_nbasic(prob);
1638 #endif
1639
1640 /*
1641 For each record, restore the row bounds. If we have status arrays, check
1642 the status of the logical and adjust if necessary.
1643
1644 In spite of the fact that the status array is an unsigned char array,
1645 we still need to use getRowStatus to make sure we're only looking at the
1646 bottom three bits. Why is this an issue? Because the status array isn't
1647 necessarily cleared to zeros, and setRowStatus carefully changes only
1648 the bottom three bits!
1649 */
1650 for (int k = 0; k < numRecs; k++) {
1651 const action &bndRec = bndRecords[k];
1652 const int &i = bndRec.ndx_;
1653 const double &rloi = bndRec.rlo_;
1654 const double &rupi = bndRec.rup_;
1655
1656 #if PRESOLVE_DEBUG > 1
1657 std::cout << "NDUAL(eq): row(" << i << ")";
1658 if (rlo[i] != rloi)
1659 std::cout << " LB " << rlo[i] << " -> " << rloi;
1660 if (rup[i] != rupi)
1661 std::cout << " UB " << rup[i] << " -> " << rupi;
1662 #endif
1663
1664 rlo[i] = rloi;
1665 rup[i] = rupi;
1666 if (rowstat) {
1667 unsigned char stati = prob->getRowStatus(i);
1668 if (stati == CoinPresolveMatrix::atUpperBound) {
1669 if (rloi <= -PRESOLVE_INF) {
1670 rowstat[i] = CoinPresolveMatrix::atLowerBound;
1671 #if PRESOLVE_DEBUG > 1
1672 std::cout
1673 << ", status forced to "
1674 << statusName(static_cast< CoinPresolveMatrix::Status >(rowstat[i]));
1675 #endif
1676 }
1677 } else if (stati == CoinPresolveMatrix::atLowerBound) {
1678 if (rupi >= PRESOLVE_INF) {
1679 rowstat[i] = CoinPresolveMatrix::atUpperBound;
1680 #if PRESOLVE_DEBUG > 1
1681 std::cout
1682 << ", status forced to "
1683 << statusName(static_cast< CoinPresolveMatrix::Status >(rowstat[i]));
1684 #endif
1685 }
1686 }
1687 #if PRESOLVE_DEBUG > 2
1688 else if (stati == CoinPresolveMatrix::basic) {
1689 std::cout << ", status is basic.";
1690 } else if (stati == CoinPresolveMatrix::isFree) {
1691 std::cout << ", status is free?!";
1692 } else {
1693 unsigned int tmp = static_cast< unsigned int >(stati);
1694 std::cout << ", status is invalid (" << tmp << ")!";
1695 }
1696 #endif
1697 }
1698 #if PRESOLVE_DEBUG > 1
1699 std::cout << "." << std::endl;
1700 #endif
1701 }
1702
1703 #if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
1704 presolve_check_threads(prob);
1705 presolve_check_sol(prob, 2, 2, 2);
1706 presolve_check_nbasic(prob);
1707 #if PRESOLVE_DEBUG > 0
1708 std::cout << "Leaving remove_dual_action::postsolve." << std::endl;
1709 #endif
1710 #endif
1711
1712 return;
1713 }
1714
1715 /*
1716 Destructor
1717 */
~remove_dual_action()1718 remove_dual_action::~remove_dual_action()
1719 {
1720 deleteAction(actions_, action *);
1721 }
1722
1723 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1724 */
1725