1/* $Id: CoinPresolveDual.cpp 1373 2011-01-03 23:57:44Z lou $ */
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//#define PRESOLVE_TIGHTEN_DUALS 1
16//#define PRESOLVE_DEBUG 1
17
18#ifdef PRESOLVE_DEBUG
19#include "CoinPresolvePsdebug.hpp"
20#endif
21
22// this looks for "dominated columns"
23// following ekkredc
24// rdmin == dwrow2
25// rdmax == dwrow
26// this transformation is applied in: k4.mps, lp22.mps
27
28// make inferences using this constraint on reduced cost:
29// at optimality dj>0 ==> var must be at lb
30// dj<0 ==> var must be at ub
31//
32// This implies:
33// there is no lb ==> dj<=0 at optimality
34// there is no ub ==> dj>=0 at optimality
35
36/*
37 This routine looks to be something of a work in progres. See the comment
38 that begins with `Gack!'. And down in the bound propagation loop, why do we
39 only work with variables with u_j = infty? The corresponding section of code
40 for l_j = -infty is ifdef'd away. And why exclude the code protected by
41 PRESOLVE_TIGHTEN_DUALS? And why are we using ekkinf instead of PRESOLVE_INF?
42
43 There is no postsolve action because once we've identified a variable to fix
44 we can invoke make_fixed_action.
45*/
46const CoinPresolveAction *remove_dual_action::presolve(CoinPresolveMatrix *prob,
47 const CoinPresolveAction *next)
48{
49 double startTime = 0.0;
50 int startEmptyRows=0;
51 int startEmptyColumns = 0;
52 if (prob->tuning_) {
53 startTime = CoinCpuTime();
54 startEmptyRows = prob->countEmptyRows();
55 startEmptyColumns = prob->countEmptyCols();
56 }
57 double *colels = prob->colels_;
58 int *hrow = prob->hrow_;
59 CoinBigIndex *mcstrt = prob->mcstrt_;
60 int *hincol = prob->hincol_;
61 int ncols = prob->ncols_;
62
63 double *clo = prob->clo_;
64 double *cup = prob->cup_;
65 unsigned char *colstat = prob->colstat_ ;
66
67 // Used only in `fix if simple' section below. Remove dec'l to avoid
68 // GCC compile warning.
69 // double *rowels = prob->rowels_;
70 int *hcol = prob->hcol_;
71 CoinBigIndex *mrstrt = prob->mrstrt_;
72 int *hinrow = prob->hinrow_;
73 double *csol = prob->sol_;
74 int nrows = prob->nrows_;
75
76 double *rlo = prob->rlo_;
77 double *rup = prob->rup_;
78
79 double *dcost = prob->cost_;
80 const unsigned char *integerType = prob->integerType_;
81
82 const double maxmin = prob->maxmin_;
83
84 const double ekkinf = 1e28;
85 const double ekkinf2 = 1e20;
86 const double ztoldj = prob->ztoldj_;
87
88 CoinRelFltEq relEq(prob->ztolzb_) ;
89
90 double *rdmin = prob->usefulRowDouble_; //new double[nrows];
91 double *rdmax = reinterpret_cast<double *> (prob->usefulRowInt_); //new double[nrows];
92
93# if PRESOLVE_DEBUG
94 std::cout << "Entering remove_dual_action::presolve, " << nrows << " X " << ncols << "." << std::endl ;
95 presolve_check_sol(prob) ;
96 presolve_check_nbasic(prob) ;
97# endif
98
99 // This combines the initialization of rdmin/rdmax to extreme values
100 // (PRESOLVE_INF/-PRESOLVE_INF) with a version of the next loop specialized
101 // for row slacks.
102 // In this case, it is always the case that dprice==0.0 and coeff==1.0.
103 int i;
104 for ( i = 0; i < nrows; i++) {
105 double sup = -rlo[i]; // slack ub; corresponds to cup[j]
106 double slo = -rup[i]; // slack lb; corresponds to clo[j]
107 bool no_lb = (slo <= -ekkinf);
108 bool no_ub = (sup >= ekkinf);
109
110 // here, dj==-row dual
111 // the row slack has no lower bound, but it does have an upper bound,
112 // then the reduced cost must be <= 0, so the row dual must be >= 0
113 rdmin[i] = (no_lb && !no_ub) ? 0.0 : -PRESOLVE_INF;
114
115 rdmax[i] = (no_ub && !no_lb) ? 0.0 : PRESOLVE_INF;
116 }
117
118 // Look for col singletons and update bounds on dual costs
119 // Take the min of the maxes and max of the mins
120 int j;
121 for ( j = 0; j<ncols; j++) {
122 if (integerType[j])
123 continue; // even if it has infinite bound now ....
124 bool no_ub = (cup[j] >= ekkinf);
125 bool no_lb = (clo[j] <= -ekkinf);
126
127 if (hincol[j] == 1 &&
128
129 // we need singleton cols that have exactly one bound
130 (no_ub ^ no_lb)) {
131 int row = hrow[mcstrt[j]];
132 double coeff = colels[mcstrt[j]];
133
134 PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
135
136 // I don't think the sign of dcost[j] matters
137
138 // this row dual would make this col's reduced cost be 0
139 double dprice = maxmin * dcost[j] / coeff;
140
141 // no_ub == !no_lb
142 // no_ub ==> !(dj<0)
143 // no_lb ==> !(dj>0)
144 // I don't think the case where !no_ub has actually been tested
145 if ((coeff > 0.0) == no_ub) {
146 // coeff>0 ==> making the row dual larger would make dj *negative*
147 // ==> dprice is an upper bound on dj if no *ub*
148 // coeff<0 ==> making the row dual larger would make dj *positive*
149 // ==> dprice is an upper bound on dj if no *lb*
150 if (rdmax[row] > dprice) // reduced cost may be positive
151 rdmax[row] = dprice;
152 } else if ((coeff < 0.0) == no_lb) { // no lower bound
153 if (rdmin[row] < dprice) // reduced cost may be negative
154 rdmin[row] = dprice;
155 }
156 }
157 }
158
159 int *fix_cols = prob->usefulColumnInt_; //new int[ncols];
160
161 //int *fixdown_cols = new int[ncols];
162
163#if PRESOLVE_TIGHTEN_DUALS
164 double *djmin = new double[ncols];
165 double *djmax = new double[ncols];
166#endif
167 int nfixup_cols = 0;
168 int nfixdown_cols = ncols;
169
170 int nPass=100;
171 while (nPass-- > 0) {
172 int tightened = 0;
173 /* Perform duality tests */
174 for (int j = 0; j<ncols; j++) {
175 if (hincol[j] > 0) {
176 CoinBigIndex kcs = mcstrt[j];
177 CoinBigIndex kce = kcs + hincol[j];
178 // Number of infinite rows
179 int nflagu = 0;
180 int nflagl = 0;
181 // Number of ordinary rows
182 int nordu = 0;
183 int nordl = 0;
184 double ddjlo = maxmin * dcost[j];
185 double ddjhi = ddjlo;
186
187 for (CoinBigIndex k = kcs; k < kce; k++) {
188 int i = hrow[k];
189 double coeff = colels[k];
190
191 if (coeff > 0.0) {
192 if (rdmin[i] >= -ekkinf2) {
193 ddjhi -= coeff * rdmin[i];
194 nordu ++;
195 } else {
196 nflagu ++;
197 }
198 if (rdmax[i] <= ekkinf2) {
199 ddjlo -= coeff * rdmax[i];
200 nordl ++;
201 } else {
202 nflagl ++;
203 }
204 } else {
205 if (rdmax[i] <= ekkinf2) {
206 ddjhi -= coeff * rdmax[i];
207 nordu ++;
208 } else {
209 nflagu ++;
210 }
211 if (rdmin[i] >= -ekkinf2) {
212 ddjlo -= coeff * rdmin[i];
213 nordl ++;
214 } else {
215 nflagl ++;
216 }
217 }
218 }
219 // See if we may be able to tighten a dual
220 if (!integerType[j]) {
221 if (cup[j]>ekkinf) {
222 // dj can not be negative
223 if (nflagu==1&&ddjhi<-ztoldj) {
224 // We can make bound finite one way
225 for (CoinBigIndex k = kcs; k < kce; k++) {
226 int i = hrow[k];
227 double coeff = colels[k];
228
229 if (coeff > 0.0&&rdmin[i] < -ekkinf2) {
230 // rdmax[i] has upper bound
231 if (ddjhi<rdmax[i]*coeff-ztoldj) {
232 double newValue = ddjhi/coeff;
233 // re-compute lo
234 if (rdmax[i] > ekkinf2 && newValue <= ekkinf2) {
235 nflagl--;
236 ddjlo -= coeff * newValue;
237 } else if (rdmax[i] <= ekkinf2) {
238 ddjlo -= coeff * (newValue-rdmax[i]);
239 }
240 rdmax[i] = newValue;
241 tightened++;
242#if PRESOLVE_DEBUG
243 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
244#endif
245 }
246 } else if (coeff < 0.0 && rdmax[i] > ekkinf2) {
247 // rdmin[i] has lower bound
248 if (ddjhi<rdmin[i]*coeff-ztoldj) {
249 double newValue = ddjhi/coeff;
250 // re-compute lo
251 if (rdmin[i] < -ekkinf2 && newValue >= -ekkinf2) {
252 nflagl--;
253 ddjlo -= coeff * newValue;
254 } else if (rdmax[i] >= -ekkinf2) {
255 ddjlo -= coeff * (newValue-rdmin[i]);
256 }
257 rdmin[i] = newValue;
258 tightened++;
259#if PRESOLVE_DEBUG
260 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
261#endif
262 ddjlo = 0.0;
263 }
264 }
265 }
266 } else if (nflagl==0&&nordl==1&&ddjlo<-ztoldj) {
267 // We may be able to tighten
268 for (CoinBigIndex k = kcs; k < kce; k++) {
269 int i = hrow[k];
270 double coeff = colels[k];
271
272 if (coeff > 0.0) {
273 rdmax[i] += ddjlo/coeff;
274 ddjlo =0.0;
275 tightened++;
276#if PRESOLVE_DEBUG
277 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
278#endif
279 } else if (coeff < 0.0 ) {
280 rdmin[i] += ddjlo/coeff;
281 ddjlo =0.0;
282 tightened++;
283#if PRESOLVE_DEBUG
284 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
285#endif
286 }
287 }
288 }
289 }
290#if 0
291 if (clo[j]<-ekkinf) {
292 // dj can not be positive
293 if (ddjlo>ztoldj&&nflagl==1) {
294 // We can make bound finite one way
295 for (CoinBigIndex k = kcs; k < kce; k++) {
296 int i = hrow[k];
297 double coeff = colels[k];
298
299 if (coeff < 0.0&&rdmin[i] < -ekkinf2) {
300 // rdmax[i] has upper bound
301 if (ddjlo>rdmax[i]*coeff+ztoldj) {
302 double newValue = ddjlo/coeff;
303 // re-compute hi
304 if (rdmax[i] > ekkinf2 && newValue <= ekkinf2) {
305 nflagu--;
306 ddjhi -= coeff * newValue;
307 } else if (rdmax[i] <= ekkinf2) {
308 ddjhi -= coeff * (newValue-rdmax[i]);
309 }
310 rdmax[i] = newValue;
311 tightened++;
312#if PRESOLVE_DEBUG
313 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
314#endif
315 }
316 } else if (coeff > 0.0 && rdmax[i] > ekkinf2) {
317 // rdmin[i] has lower bound
318 if (ddjlo>rdmin[i]*coeff+ztoldj) {
319 double newValue = ddjlo/coeff;
320 // re-compute lo
321 if (rdmin[i] < -ekkinf2 && newValue >= -ekkinf2) {
322 nflagu--;
323 ddjhi -= coeff * newValue;
324 } else if (rdmax[i] >= -ekkinf2) {
325 ddjhi -= coeff * (newValue-rdmin[i]);
326 }
327 rdmin[i] = newValue;
328 tightened++;
329#if PRESOLVE_DEBUG
330 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
331#endif
332 }
333 }
334 }
335 } else if (nflagu==0&&nordu==1&&ddjhi>ztoldj) {
336 // We may be able to tighten
337 for (CoinBigIndex k = kcs; k < kce; k++) {
338 int i = hrow[k];
339 double coeff = colels[k];
340
341 if (coeff < 0.0) {
342 rdmax[i] += ddjhi/coeff;
343 ddjhi =0.0;
344 tightened++;
345#if PRESOLVE_DEBUG
346 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
347#endif
348 } else if (coeff > 0.0 ) {
349 rdmin[i] += ddjhi/coeff;
350 ddjhi =0.0;
351 tightened++;
352#if PRESOLVE_DEBUG
353 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
354#endif
355 }
356 }
357 }
358 }
359#endif
360 }
361
362#if PRESOLVE_TIGHTEN_DUALS
363 djmin[j] = (nflagl ? -PRESOLVE_INF : ddjlo);
364 djmax[j] = (nflagu ? PRESOLVE_INF : ddjhi);
365#endif
366
367 if (ddjlo > ztoldj && nflagl == 0&&!prob->colProhibited2(j)) {
368 // dj>0 at optimality ==> must be at lower bound
369 if (clo[j] <= -ekkinf) {
370 prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDB,
371 prob->messages())
372 <<j
373 <<CoinMessageEol;
374 prob->status_ |= 2;
375 break;
376 } else {
377 fix_cols[--nfixdown_cols] = j;
378 PRESOLVE_DETAIL_PRINT(printf("pre_duallo %dC E\n",j));
379# if PRESOLVE_DEBUG
380 printf("NDUAL: fixing x<%d>",fix_cols[nfixdown_cols]) ;
381 if (csol) printf(" = %g",csol[j]) ;
382 printf(" at lb = %g.\n",clo[j]) ;
383# endif
384 //if (csol[j]-clo[j]>1.0e-7)
385 //printf("down %d row %d nincol %d\n",j,hrow[mcstrt[j]],hincol[j]);
386 // User may have given us feasible solution - move if simple
387 if (csol) {
388# if 0
389 /*
390 Except it's not simple. The net result is that we end up with an
391 excess of basic variables.
392 */
393 if (csol[j]-clo[j]>1.0e-7&&hincol[j]==1) {
394 double value_j = colels[mcstrt[j]];
395 double distance_j = csol[j]-clo[j];
396 int row=hrow[mcstrt[j]];
397 // See if another column can take value
398 for (CoinBigIndex kk=mrstrt[row];kk<mrstrt[row]+hinrow[row];kk++) {
399 int k = hcol[kk];
400 if (colstat[k] == CoinPrePostsolveMatrix::superBasic)
401 continue ;
402
403 if (hincol[k]==1&&k!=j) {
404 double value_k = rowels[kk];
405 double movement;
406 if (value_k*value_j>0.0) {
407 // k needs to increase
408 double distance_k = cup[k]-csol[k];
409 movement = CoinMin((distance_j*value_j)/value_k,distance_k);
410 } else {
411 // k needs to decrease
412 double distance_k = clo[k]-csol[k];
413 movement = CoinMax((distance_j*value_j)/value_k,distance_k);
414 }
415 if (relEq(movement,0)) continue ;
416
417 csol[k] += movement;
418 if (relEq(csol[k],clo[k]))
419 { colstat[k] = CoinPrePostsolveMatrix::atLowerBound ; }
420 else
421 if (relEq(csol[k],cup[k]))
422 { colstat[k] = CoinPrePostsolveMatrix::atUpperBound ; }
423 else
424 if (colstat[k] != CoinPrePostsolveMatrix::isFree)
425 { colstat[k] = CoinPrePostsolveMatrix::basic ; }
426 printf("NDUAL: x<%d> moved %g to %g; ",
427 k,movement,csol[k]) ;
428 printf("lb = %g, ub = %g, status now %s.\n",
429 clo[k],cup[k],columnStatusString(colstat[k])) ;
430 distance_j -= (movement*value_k)/value_j;
431 csol[j] -= (movement*value_k)/value_j;
432 if (distance_j<1.0e-7)
433 break;
434 }
435 }
436 }
437# endif // repair solution.
438
439 csol[j] = clo[j] ; // but the bottom line is we've changed x<j>
440 colstat[j] = CoinPrePostsolveMatrix::atLowerBound ;
441 }
442 }
443 } else if (ddjhi < -ztoldj && nflagu == 0&&!prob->colProhibited2(j)) {
444 // dj<0 at optimality ==> must be at upper bound
445 if (cup[j] >= ekkinf) {
446 prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDA,
447 prob->messages())
448 <<j
449 <<CoinMessageEol;
450 prob->status_ |= 2;
451 break;
452 } else {
453 PRESOLVE_DETAIL_PRINT(printf("pre_dualup %dC E\n",j));
454 fix_cols[nfixup_cols++] = j;
455# if PRESOLVE_DEBUG
456 printf("NDUAL: fixing x<%d>",fix_cols[nfixup_cols-1]) ;
457 if (csol) printf(" = %g",csol[j]) ;
458 printf(" at ub = %g.\n",cup[j]) ;
459# endif
460 // User may have given us feasible solution - move if simple
461 // See comments for `fix at lb' case above.
462 //if (cup[j]-csol[j]>1.0e-7)
463 //printf("up %d row %d nincol %d\n",j,hrow[mcstrt[j]],hincol[j]);
464 if (csol) {
465# if 0
466 // See comments above.
467 if (cup[j]-csol[j]>1.0e-7&&hincol[j]==1) {
468 double value_j = colels[mcstrt[j]];
469 double distance_j = csol[j]-cup[j];
470 int row=hrow[mcstrt[j]];
471 // See if another column can take value
472 for (CoinBigIndex kk=mrstrt[row];kk<mrstrt[row]+hinrow[row];kk++) {
473 int k = hcol[kk];
474 if (colstat[k] == CoinPrePostsolveMatrix::superBasic)
475 continue ;
476
477 if (hincol[k]==1&&k!=j) {
478 double value_k = rowels[kk];
479 double movement;
480 if (value_k*value_j<0.0) {
481 // k needs to increase
482 double distance_k = cup[k]-csol[k];
483 movement = CoinMin((distance_j*value_j)/value_k,distance_k);
484 } else {
485 // k needs to decrease
486 double distance_k = clo[k]-csol[k];
487 movement = CoinMax((distance_j*value_j)/value_k,distance_k);
488 }
489 if (relEq(movement,0)) continue ;
490
491 csol[k] += movement;
492 if (relEq(csol[k],clo[k]))
493 { colstat[k] = CoinPrePostsolveMatrix::atLowerBound ; }
494 else
495 if (relEq(csol[k],cup[k]))
496 { colstat[k] = CoinPrePostsolveMatrix::atUpperBound ; }
497 else
498 if (colstat[k] != CoinPrePostsolveMatrix::isFree)
499 { colstat[k] = CoinPrePostsolveMatrix::basic ; }
500 printf("NDUAL: x<%d> moved %g to %g; ",
501 k,movement,csol[k]) ;
502 printf("lb = %g, ub = %g, status now %s.\n",
503 clo[k],cup[k],columnStatusString(colstat[k])) ;
504 distance_j -= (movement*value_k)/value_j;
505 csol[j] -= (movement*value_k)/value_j;
506 if (distance_j>-1.0e-7)
507 break;
508 }
509 }
510 }
511# endif
512 csol[j] = cup[j] ; // but the bottom line is we've changed x<j>
513 colstat[j] = CoinPrePostsolveMatrix::atUpperBound ;
514 }
515 }
516 }
517 }
518 }
519
520 // I don't know why I stopped doing this.
521#if PRESOLVE_TIGHTEN_DUALS
522 const double *rowels = prob->rowels_;
523 const int *hcol = prob->hcol_;
524 const CoinBigIndex *mrstrt = prob->mrstrt_;
525 int *hinrow = prob->hinrow_;
526 // tighten row dual bounds, as described on p. 229
527 for (int i = 0; i < nrows; i++) {
528 bool no_ub = (rup[i] >= ekkinf);
529 bool no_lb = (rlo[i] <= -ekkinf);
530
531 if ((no_ub ^ no_lb) == true) {
532 CoinBigIndex krs = mrstrt[i];
533 CoinBigIndex kre = krs + hinrow[i];
534 double rmax = rdmax[i];
535 double rmin = rdmin[i];
536
537 // all row columns are non-empty
538 for (CoinBigIndex k=krs; k<kre; k++) {
539 double coeff = rowels[k];
540 int icol = hcol[k];
541 double djmax0 = djmax[icol];
542 double djmin0 = djmin[icol];
543
544 if (no_ub) {
545 // dj must not be negative
546 if (coeff > ZTOLDP2 && djmax0 <PRESOLVE_INF && cup[icol]>=ekkinf) {
547 double bnd = djmax0 / coeff;
548 if (rmax > bnd) {
549#if PRESOLVE_DEBUG
550 printf("MAX TIGHT[%d,%d]: %g --> %g\n", i,hrow[k], rdmax[i], bnd);
551#endif
552 rdmax[i] = rmax = bnd;
553 tightened ++;
554 }
555 } else if (coeff < -ZTOLDP2 && djmax0 <PRESOLVE_INF && cup[icol] >= ekkinf) {
556 double bnd = djmax0 / coeff ;
557 if (rmin < bnd) {
558#if PRESOLVE_DEBUG
559 printf("MIN TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], rdmin[i], bnd);
560#endif
561 rdmin[i] = rmin = bnd;
562 tightened ++;
563 }
564 }
565 } else { // no_lb
566 // dj must not be positive
567 if (coeff > ZTOLDP2 && djmin0 > -PRESOLVE_INF && clo[icol]<=-ekkinf) {
568 double bnd = djmin0 / coeff ;
569 if (rmin < bnd) {
570#if PRESOLVE_DEBUG
571 printf("MIN1 TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], rdmin[i], bnd);
572#endif
573 rdmin[i] = rmin = bnd;
574 tightened ++;
575 }
576 } else if (coeff < -ZTOLDP2 && djmin0 > -PRESOLVE_INF && clo[icol] <= -ekkinf) {
577 double bnd = djmin0 / coeff ;
578 if (rmax > bnd) {
579#if PRESOLVE_DEBUG
580 printf("MAX TIGHT1[%d,%d]: %g --> %g\n", i,hrow[k], rdmax[i], bnd);
581#endif
582 rdmax[i] = rmax = bnd;
583 tightened ++;
584 }
585 }
586 }
587 }
588 }
589 }
590#endif
591
592 if (tightened<100||nfixdown_cols<ncols||nfixup_cols)
593 break;
594#if PRESOLVE_TIGHTEN_DUALS
595 else
596 printf("DUAL TIGHTENED! %d\n", tightened);
597#endif
598 }
599 assert (nfixup_cols<=nfixdown_cols);
600 if (nfixup_cols) {
601#if PRESOLVE_DEBUG
602 printf("NDUAL: %d up", nfixup_cols);
603 for (i = 0 ; i < nfixup_cols ; i++) printf(" %d",fix_cols[i]) ;
604 printf(".\n") ;
605#endif
606 next = make_fixed_action::presolve(prob, fix_cols, nfixup_cols, false, next);
607 }
608
609 if (nfixdown_cols<ncols) {
610 int * fixdown_cols = fix_cols+nfixdown_cols;
611 nfixdown_cols = ncols-nfixdown_cols;
612#if PRESOLVE_DEBUG
613 printf("NDUAL: %d down", nfixdown_cols);
614 for (i = 0 ; i < nfixdown_cols ; i++) printf(" %d",fixdown_cols[i]) ;
615 printf(".\n") ;
616#endif
617 next = make_fixed_action::presolve(prob, fixdown_cols, nfixdown_cols, true, next);
618 }
619 // If dual says so then we can make equality row
620 // Also if cost is in right direction and only one binding row for variable
621 // We may wish to think about giving preference to rows with 2 or 3 elements
622/*
623 Gack! Ok, I can appreciate the thought here, but I'm seriously skeptical
624 about writing canFix[0] before reading rdmin[0]. After that, we should be out
625 of the interference zone for the typical situation where sizeof(double) is
626 twice sizeof(int).
627*/
628 int * canFix = reinterpret_cast<int *> (rdmin);
629 for ( i = 0; i < nrows; i++) {
630 bool no_lb = (rlo[i] <= -ekkinf);
631 bool no_ub = (rup[i] >= ekkinf);
632 canFix[i]=0;
633 if (no_ub && !no_lb ) {
634 if ( rdmin[i]>0.0)
635 canFix[i]=-1;
636 else
637 canFix[i]=-2;
638 } else if (no_lb && !no_ub ) {
639 if (rdmax[i]<0.0)
640 canFix[i]=1;
641 else
642 canFix[i]=2;
643 }
644 }
645 for (j = 0; j<ncols; j++) {
646 if (hincol[j]<=1)
647 continue;
648 if (integerType[j])
649 continue; // even if it has infinite bound now ....
650 CoinBigIndex kcs = mcstrt[j];
651 CoinBigIndex kce = kcs + hincol[j];
652 int bindingUp=-1;
653 int bindingDown=-1;
654 if (cup[j]<ekkinf)
655 bindingUp=-2;
656 if (clo[j]>-ekkinf)
657 bindingDown=-2;
658 for (CoinBigIndex k = kcs; k < kce; k++) {
659 int i = hrow[k];
660 if (abs(canFix[i])!=2) {
661 bindingUp=-2;
662 bindingDown=-2;
663 break;
664 }
665 double coeff = colels[k];
666 if (coeff>0.0) {
667 if (canFix[i]==2) {
668 // binding up
669 if (bindingUp==-1)
670 bindingUp = i;
671 else
672 bindingUp = -2;
673 } else {
674 // binding down
675 if (bindingDown==-1)
676 bindingDown = i;
677 else
678 bindingDown = -2;
679 }
680 } else {
681 if (canFix[i]==2) {
682 // binding down
683 if (bindingDown==-1)
684 bindingDown = i;
685 else
686 bindingDown = -2;
687 } else {
688 // binding up
689 if (bindingUp==-1)
690 bindingUp = i;
691 else
692 bindingUp = -2;
693 }
694 }
695 }
696 double cost = maxmin * dcost[j];
697 if (bindingUp>-2&&cost<=0.0) {
698 // might as well make equality
699 if (bindingUp>=0) {
700 canFix[bindingUp] /= 2; //So -2 goes to -1 etc
701 //printf("fixing row %d to ub by %d\n",bindingUp,j);
702 } else {
703 //printf("binding up row by %d\n",j);
704 }
705 } else if (bindingDown>-2 &&cost>=0.0) {
706 // might as well make equality
707 if (bindingDown>=0) {
708 canFix[bindingDown] /= 2; //So -2 goes to -1 etc
709 //printf("fixing row %d to lb by %d\n",bindingDown,j);
710 } else {
711 //printf("binding down row by %d\n",j);
712 }
713 }
714 }
715 // can't fix if integer and non-unit coefficient
716 //const double *rowels = prob->rowels_;
717 //const int *hcol = prob->hcol_;
718 //const CoinBigIndex *mrstrt = prob->mrstrt_;
719 //int *hinrow = prob->hinrow_;
720 for ( i = 0; i < nrows; i++) {
721 if (abs(canFix[i])==1) {
722 CoinBigIndex krs = mrstrt[i];
723 CoinBigIndex kre = krs + hinrow[i];
724 for (CoinBigIndex k=krs; k<kre; k++) {
725 int icol = hcol[k];
726 if (cup[icol]>clo[icol]&&integerType[icol]) {
727 canFix[i]=0; // not safe
728#ifdef COIN_DEVELOP
729 printf("no dual something CoinPresolveDual row %d col %d\n",
730 i,icol);
731#endif
732 }
733 }
734 }
735 if (canFix[i]==1) {
736 rlo[i]=rup[i];
737 prob->addRow(i);
738 } else if (canFix[i]==-1) {
739 rup[i]=rlo[i];
740 prob->addRow(i);
741 }
742 }
743
744 //delete[]rdmin;
745 //delete[]rdmax;
746
747 //delete[]fixup_cols;
748 //delete[]fixdown_cols;
749
750#if PRESOLVE_TIGHTEN_DUALS
751 delete[]djmin;
752 delete[]djmax;
753#endif
754
755 if (prob->tuning_) {
756 double thisTime=CoinCpuTime();
757 int droppedRows = prob->countEmptyRows() - startEmptyRows;
758 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
759 printf("CoinPresolveDual(1) - %d rows, %d columns dropped in time %g, total %g\n",
760 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
761 }
762
763# if PRESOLVE_DEBUG
764 presolve_check_sol(prob) ;
765 presolve_check_nbasic(prob) ;
766 std::cout << "Leaving remove_dual_action::presolve." << std::endl ;
767# endif
768
769 return (next);
770}
771