source: branches/devel-1/PresolveDual.cpp @ 2351

Last change on this file since 2351 was 49, checked in by ladanyi, 17 years ago

everything compiles

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