source: trunk/PresolveDual.cpp @ 55

Last change on this file since 55 was 55, checked in by forrest, 17 years ago

Making it a bit more likely Visual Studio will work

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1#include <stdio.h>
2#include <math.h>
3
4#include "PresolveMatrix.hpp"
5#include "PresolveFixed.hpp"
6#include "PresolveDual.hpp"
7#include "ClpMessage.hpp"
8
9// this looks for "dominated columns"
10// following ekkredc
11// rdmin == dwrow2
12// rdmax == dwrow
13// this transformation is applied in: k4.mps, lp22.mps
14
15// make inferences using this constraint on reduced cost:
16//      at optimality   dj>0 ==> var must be at lb
17//                      dj<0 ==> var must be at ub
18//
19// This implies:
20//      there is no lb ==> dj<=0 at optimality
21//      there is no ub ==> dj>=0 at optimality
22const PresolveAction *remove_dual_action::presolve(PresolveMatrix *prob,
23                                   const PresolveAction *next)
24{
25  double *colels        = prob->colels_;
26  int *hrow             = prob->hrow_;
27  CoinBigIndex *mcstrt          = prob->mcstrt_;
28  int *hincol           = prob->hincol_;
29  int ncols             = prob->ncols_;
30
31  double *clo   = prob->clo_;
32  double *cup   = prob->cup_;
33
34  //  double *rowels    = prob->rowels_;
35  //  int *hcol         = prob->hcol_;
36  //  CoinBigIndex *mrstrt              = prob->mrstrt_;
37  //  int *hinrow               = prob->hinrow_;
38  int nrows             = prob->nrows_;
39
40  double *rlo   = prob->rlo_;
41  double *rup   = prob->rup_;
42
43  double *dcost = prob->cost_;
44
45  const double maxmin   = prob->maxmin_;
46
47  const double ekkinf = 1e28;
48  const double ekkinf2 = 1e20;
49
50  double *rdmin = new double[nrows];
51  double *rdmax = new double[nrows];
52
53  // This combines the initialization of rdmin/rdmax to extreme values
54  // (PRESOLVE_INF/-PRESOLVE_INF) with a version of the next loop specialized
55  // for row slacks.
56  // In this case, it is always the case that dprice==0.0 and coeff==1.0.
57  for (int i = 0; i < nrows; i++) {
58    double sup = -rlo[i];       // slack ub; corresponds to cup[j]
59    double slo = -rup[i];       // slack lb; corresponds to clo[j]
60    bool no_lb = (slo <= -ekkinf);
61    bool no_ub = (sup >= ekkinf);
62
63    // here, dj==-row dual
64    // the row slack has no lower bound, but it does have an upper bound,
65    // then the reduced cost must be <= 0, so the row dual must be >= 0
66    rdmin[i] = (no_lb && !no_ub) ? 0.0 : -PRESOLVE_INF;
67
68    rdmax[i] = (no_ub && !no_lb) ? 0.0 :  PRESOLVE_INF;
69  }
70
71  // Look for col singletons and update bounds on dual costs
72  // Take the min of the maxes and max of the mins
73  for (int j = 0; j<ncols; j++) {
74    bool no_ub = (cup[j] >= ekkinf);
75    bool no_lb = (clo[j] <= -ekkinf);
76   
77    if (hincol[j] == 1 &&
78
79        // we need singleton cols that have exactly one bound
80        (no_ub ^ no_lb)) {
81      int row = hrow[mcstrt[j]];
82      double coeff = colels[mcstrt[j]];
83
84      PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
85
86      // I don't think the sign of dcost[j] matters
87
88      // this row dual would make this col's reduced cost be 0
89      double dprice = maxmin * dcost[j] / coeff;
90
91      // no_ub == !no_lb
92      // no_ub ==> !(dj<0)
93      // no_lb ==> !(dj>0)
94      // I don't think the case where !no_ub has actually been tested
95      if ((coeff > 0.0) == no_ub) {
96        // coeff>0 ==> making the row dual larger would make dj *negative*
97        //              ==> dprice is an upper bound on dj if no *ub*
98        // coeff<0 ==> making the row dual larger would make dj *positive*
99        //              ==> dprice is an upper bound on dj if no *lb*
100        if (rdmax[row] > dprice)        // reduced cost may be positive
101          rdmax[row] = dprice;
102      } else {                  // no lower bound
103        if (rdmin[row] < dprice)        // reduced cost may be negative
104          rdmin[row] = dprice;
105      }
106    }
107  }
108
109  int *fixup_cols       = new int[ncols];
110
111  int *fixdown_cols     = new int[ncols];
112
113#if     PRESOLVE_TIGHTEN_DUALS
114  double *djmin = new double[ncols];
115  double *djmax = new double[ncols];
116#endif
117  int nfixup_cols       = 0;
118  int nfixdown_cols     = 0;
119
120  while (1) {
121    /* Perform duality tests */
122    for (int j = 0; j<ncols; j++)
123      if (hincol[j] > 0) {
124
125        //if (cup[j] - clo[j] <= ztolzb && cup[j] - clo[j] > ZTOLDP) {
126        //  cup[j] = clo[j];
127        //}
128        CoinBigIndex kcs = mcstrt[j];
129        CoinBigIndex kce = kcs + hincol[j];
130        int iflagu = 0;
131        int iflagl = 0;
132        double ddjlo = maxmin * dcost[j];
133        double ddjhi = ddjlo;
134
135        for (CoinBigIndex k = kcs; k < kce; k++) {
136          int i = hrow[k];
137          double coeff = colels[k];
138
139          if (coeff > 0.0) {
140            ddjhi -= coeff * rdmin[i];
141            if (rdmin[i] >= ekkinf2) {
142              iflagu = 1;
143              ddjhi = 0.0;
144            }
145            ddjlo -= coeff * rdmax[i];
146            if (rdmax[i] <= -ekkinf2) {
147              iflagl = 1;
148              ddjlo = 0.0;
149            }
150          } else {
151            ddjhi -= coeff * rdmax[i];
152            if (rdmax[i] <= -ekkinf2) {
153              iflagu = 1;
154              ddjhi = 0.0;
155            }
156            ddjlo -= coeff * rdmin[i];
157            if (rdmin[i] >= ekkinf2) {
158              iflagl = 1;
159              ddjlo = 0.0;
160            }
161          }
162        }
163
164#if     PRESOLVE_TIGHTEN_DUALS
165        djmin[j] = (iflagl ?  PRESOLVE_INF : ddjlo);
166        djmax[j] = (iflagu ? -PRESOLVE_INF : ddjhi);
167#endif
168
169        if (ddjlo > ZTOLDP && iflagl == 0) {
170          // dj>0 at optimality ==> must be at lower bound
171          if (clo[j] <= -ekkinf) {
172            prob->messageHandler()->message(CLP_PRESOLVE_COLUMNBOUNDB,
173                                             prob->messages())
174                                               <<j
175                                               <<CoinMessageEol;
176            prob->status_ |= 2;
177            break;
178          } else {
179            fixdown_cols[nfixdown_cols++] = j;
180          }
181        } else if (ddjhi < -ZTOLDP && iflagu == 0) {
182          // dj<0 at optimality ==> must be at upper bound
183          if (cup[j] >= ekkinf) {
184            prob->messageHandler()->message(CLP_PRESOLVE_COLUMNBOUNDA,
185                                             prob->messages())
186                                               <<j
187                                               <<CoinMessageEol;
188            prob->status_ |= 2;
189            break;
190          } else {
191            fixup_cols[nfixup_cols++] = j;
192          }
193        }
194      }
195
196    int tightened = 0;
197
198    // I don't know why I stopped doing this.
199#if     PRESOLVE_TIGHTEN_DUALS
200    // tighten row dual bounds, as described on p. 229
201    for (int i = 0; i < nrows; i++) {
202      bool no_ub = (rup[i] >= ekkinf);
203      bool no_lb = (rlo[i] <= -ekkinf);
204     
205      if ((no_ub ^ no_lb) == true) {
206        CoinBigIndex krs = mrstrt[i];
207        CoinBigIndex kre = krs + hinrow[i];
208        double rmax  = rdmax[i];
209        double rmin  = rdmin[i];
210
211        // it may be that we could just use rmin/rmax, but I'm not sure,
212        // since I don't update djmax0/djmin0
213        double rmax0 = rmax;
214        double rmin0 = rmin;
215
216        // all row columns are non-empty
217        for (CoinBigIndex k=krs; k<kre; k++) {
218          double coeff = rowels[k];
219          double djmax0 = djmax[hcol[k]];
220          double djmin0 = djmin[hcol[k]];
221
222          if (no_ub) {
223            if (coeff > ZTOLDP && djmax0 > -PRESOLVE_INF && rmin0 < PRESOLVE_INF) {
224              double bnd = djmax0 / coeff + rmin0;
225              if (rmax < bnd) {
226#if     DEBUG_PRESOLVE
227                printf("MAX TIGHT[%d,%d]:  %g --> %g\n", i,hrow[k], rdmax[i], bnd);
228#endif
229                rdmax[i] = rmax = bnd;
230                tightened = 1;
231              }
232            } else if (coeff < -ZTOLDP && djmax0 > -PRESOLVE_INF && rmax0 > -PRESOLVE_INF) {
233              double bnd = djmax0 / coeff + rmax0;
234              if (rmin > bnd) {
235#if     DEBUG_PRESOLVE
236                printf("MIN TIGHT[%d,%d]:  %g --> %g\n", i, hrow[k], rdmin[i], bnd);
237#endif
238                rdmin[i] = rmin = bnd;
239                tightened = 1;
240              }         
241            }
242          } else {      // no_lb
243            if (coeff > ZTOLDP && djmin0 < PRESOLVE_INF && rmax0 > -PRESOLVE_INF) {
244              double bnd = djmin0 / coeff + rmax0;
245              if (rmin > bnd) {
246#if     DEBUG_PRESOLVE
247                printf("MIN1 TIGHT[%d,%d]:  %g --> %g\n", i, hrow[k], rdmin[i], bnd);
248#endif
249                rdmin[i] = rmin = bnd;
250                tightened = 1;
251              }
252            } else if (coeff < -ZTOLDP && djmin0 < PRESOLVE_INF && rmin0 < PRESOLVE_INF) {
253              double bnd = djmin0 / coeff + rmin0;
254              if (rmax < bnd) {
255#if     DEBUG_PRESOLVE
256                printf("MAX TIGHT1[%d,%d]:  %g --> %g\n", i,hrow[k], rdmax[i], bnd);
257#endif
258                rdmax[i] = rmax = bnd;
259                tightened = 1;
260              }
261            }
262          }
263        }
264      }
265    }
266#endif
267
268    if (!tightened)
269      break;
270#if     PRESOLVE_TIGHTEN_DUALS
271    else
272      printf("DUAL TIGHTENED!  %d\n", nactions);
273#endif
274  }
275
276  if (nfixup_cols) {
277#if     DEBUG_PRESOLVE
278    printf("NDUAL:  %d\n", nfixup_cols);
279#endif
280    next = make_fixed_action::presolve(prob, fixup_cols, nfixup_cols, false, next);
281  }
282
283  if (nfixdown_cols) {
284#if     DEBUG_PRESOLVE
285    printf("NDUAL:  %d\n", nfixdown_cols);
286#endif
287    next = make_fixed_action::presolve(prob, fixdown_cols, nfixdown_cols, true, next);
288  }
289
290  delete[]rdmin;
291  delete[]rdmax;
292
293  delete[]fixup_cols;
294  delete[]fixdown_cols;
295
296#if     PRESOLVE_TIGHTEN_DUALS
297  delete[]djmin;
298  delete[]djmax;
299#endif
300
301  return (next);
302}
Note: See TracBrowser for help on using the repository browser.