source: branches/devel-1/PresolveSingleton.cpp @ 29

Last change on this file since 29 was 29, checked in by forrest, 18 years ago

Presolve (no changes to Makefile)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <stdio.h>
5#include <math.h>
6#include <strings.h>
7
8#include "CoinHelperFunctions.hpp"
9#include "PresolveMatrix.hpp"
10
11#include "PresolveEmpty.hpp"    // for DROP_COL/DROP_ROW
12#include "PresolveFixed.hpp"
13#include "PresolveSingleton.hpp"
14#include "ClpMessage.hpp"
15
16/*
17 * Transfers singleton row bound information to the corresponding column bounds.
18 * What I refer to as a row singleton would be called a doubleton
19 * in the paper, since my terminology doesn't refer to the slacks.
20 * In terms of the paper, we transfer the bounds of the slack onto
21 * the variable (vii) and then "substitute" the slack out of the problem
22 * (which is a noop).
23 */
24const PresolveAction *
25slack_doubleton_action::presolve(PresolveMatrix *prob,
26                                 const PresolveAction *next,
27                                 bool & notFinished)
28{
29  double *colels        = prob->colels_;
30  int *hrow             = prob->hrow_;
31  int *mcstrt           = prob->mcstrt_;
32  int *hincol           = prob->hincol_;
33  int ncols             = prob->ncols_;
34
35  double *clo           = prob->clo_;
36  double *cup           = prob->cup_;
37
38  double *rowels        = prob->rowels_;
39  const int *hcol       = prob->hcol_;
40  const int *mrstrt     = prob->mrstrt_;
41  int *hinrow           = prob->hinrow_;
42  int nrows             = prob->nrows_;
43
44  double *rlo           = prob->rlo_;
45  double *rup           = prob->rup_;
46
47  // If rowstat exists then all do
48  unsigned char *rowstat        = prob->rowstat_;
49  double *acts  = prob->acts_;
50  double * sol = prob->sol_;
51  unsigned char * colstat = prob->colstat_;
52
53  const char *integerType = prob->integerType_;
54
55  const double ztolzb   = prob->ztolzb_;
56
57  int numberLook = prob->numberRowsToDo_;
58  int iLook;
59  int * look = prob->rowsToDo_;
60
61  action actions[MAX_SLACK_DOUBLETONS];
62  int nactions = 0;
63  notFinished = false;
64
65  int fixed_cols[ncols];
66  int nfixed_cols       = 0;
67
68  // this loop can apparently create new singleton rows;
69  // I haven't bothered to detect this.
70  // wasfor (int irow=0; irow<nrows; irow++)
71  for (iLook=0;iLook<numberLook;iLook++) {
72    int irow = look[iLook];
73    if (hinrow[irow] == 1) {
74      int jcol = hcol[mrstrt[irow]];
75      double coeff = rowels[mrstrt[irow]];
76      double lo = rlo[irow];
77      double up = rup[irow];
78      double acoeff = fabs(coeff);
79      const bool singleton_col = (hincol[jcol] == 1);
80
81      if (acoeff < ZTOLDP)
82        continue;
83
84      // don't bother with fixed cols
85      if (fabs(cup[jcol] - clo[jcol]) < ztolzb)
86        continue;
87
88      {
89        // put column on stack of things to do next time
90        prob->addCol(jcol);
91        action *s = &actions[nactions];
92        nactions++;
93
94        s->col = jcol;
95        s->clo = clo[jcol];
96        s->cup = cup[jcol];
97
98        s->row = irow;
99        s->rlo = rlo[irow];
100        s->rup = rup[irow];
101
102        s->coeff = coeff;
103      }
104
105      if (coeff < 0.0) {
106        swap(lo, up);
107        lo = -lo;
108        up = -up;
109      }
110     
111      if (lo <= -PRESOLVE_INF)
112        lo = -PRESOLVE_INF;
113      else {
114        lo /= acoeff;
115        if (lo <= -PRESOLVE_INF)
116          lo = -PRESOLVE_INF;
117      }
118
119      if (up > PRESOLVE_INF)
120        up = PRESOLVE_INF;
121      else {
122        up /= acoeff;
123        if (up > PRESOLVE_INF)
124          up = PRESOLVE_INF;
125      }
126     
127      if (clo[jcol] < lo)
128        clo[jcol] = lo;
129       
130      if (cup[jcol] > up) 
131        cup[jcol] = up;
132
133      if (fabs(cup[jcol] - clo[jcol]) < ZTOLDP) {
134        fixed_cols[nfixed_cols++] = jcol;
135      }
136
137      if (lo > up) {
138        if (lo <= up + prob->feasibilityTolerance_) {
139          // If close to integer then go there
140          double nearest = floor(lo+0.5);
141          if (fabs(nearest-lo)<2.0*prob->feasibilityTolerance_) {
142            lo = nearest;
143            up = nearest;
144          } else {
145            lo = up;
146          }
147          clo[jcol] = lo;
148          cup[jcol] = up;
149        } else {
150          prob->status_ |= 1;
151          prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_COLINFEAS,
152                                             prob->originalModel_->messages())
153                                               <<jcol
154                                               <<lo
155                                               <<up
156                                               <<CoinMessageEol;
157          break;
158        }
159      }
160
161#if     DEBUG_PRESOLVE
162      printf("SINGLETON R-%d C-%d\n", irow, jcol);
163#endif
164
165      // eliminate the row entirely from the row rep
166      // rlinks don't seem to be used PRESOLVE_REMOVE_LINK(rlink, irow);
167      hinrow[irow] = 0;
168
169      // just to make things squeeky
170      rlo[irow] = 0.0;
171      rup[irow] = 0.0;
172
173      if (rowstat) {
174        // update solution and basis
175        int basisChoice=0;
176        int numberBasic=0;
177        if (prob->columnIsBasic(jcol))
178          numberBasic++;
179        if (prob->rowIsBasic(irow))
180          numberBasic++;
181        if (sol[jcol]<=clo[jcol]+ztolzb) {
182          sol[jcol] =clo[jcol];
183          prob->setColumnStatus(jcol,PrePostsolveMatrix::atLowerBound);
184        } else if (sol[jcol]>=cup[jcol]-ztolzb) {
185          sol[jcol] =cup[jcol];
186          prob->setColumnStatus(jcol,PrePostsolveMatrix::atUpperBound);
187        } else {
188          basisChoice=1;
189        }
190        if (numberBasic>1||basisChoice)
191          prob->setColumnStatus(jcol,PrePostsolveMatrix::basic);
192      }
193
194      // remove the row from this col in the col rep
195      presolve_delete_from_row(jcol, irow, mcstrt, hincol, hrow, colels);
196
197      if (nactions >= MAX_SLACK_DOUBLETONS) {
198        notFinished=true;
199        break;
200      }
201    }
202  }
203
204  if (nactions) {
205#if     PRESOLVE_SUMMARY
206    printf("SINGLETON ROWS:  %d\n", nactions);
207#endif
208    action *save_actions = new action[nactions];
209    CoinDisjointCopyN(actions, nactions, save_actions);
210    next = new slack_doubleton_action(nactions, save_actions, next);
211
212    if (nfixed_cols)
213      next = make_fixed_action::presolve(prob, fixed_cols, nfixed_cols,
214                                         true, // arbitrary
215                                         next);
216  }
217  return (next);
218}
219
220void slack_doubleton_action::postsolve(PostsolveMatrix *prob) const
221{
222  const action *const actions = actions_;
223  const int nactions = nactions_;
224
225  double *colels        = prob->colels_;
226  int *hrow             = prob->hrow_;
227  int *mcstrt           = prob->mcstrt_;
228  int *hincol           = prob->hincol_;
229  int *link             = prob->link_;
230  int ncols             = prob->ncols_;
231
232  double *clo           = prob->clo_;
233  double *cup           = prob->cup_;
234
235  double *rlo           = prob->rlo_;
236  double *rup           = prob->rup_;
237
238  double *sol           = prob->sol_;
239  double *rcosts        = prob->rcosts_;
240
241  double *acts          = prob->acts_;
242  double *rowduals      = prob->rowduals_;
243
244  unsigned char *colstat                = prob->colstat_;
245  unsigned char *rowstat                = prob->rowstat_;
246
247  char *cdone           = prob->cdone_;
248  char *rdone           = prob->rdone_;
249  int free_list         = prob->free_list_;
250
251  const double ztolzb   = prob->ztolzb_;
252
253  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
254    int irow = f->row;
255    double lo0 = f->clo;
256    double up0 = f->cup;
257    double coeff = f->coeff;
258    int jcol = f->col;
259
260    /* the column was in the reduced problem */
261    PRESOLVEASSERT(cdone[jcol] && rdone[irow]==DROP_ROW);
262
263    rlo[irow] = f->rlo;
264    rup[irow] = f->rup;
265
266    clo[jcol] = lo0;
267    cup[jcol] = up0;
268
269    acts[irow] = coeff * sol[jcol];
270
271    // copy col to end to make room for new element
272    {
273      int k = free_list;
274      free_list = link[free_list];
275
276      check_free_list(free_list);
277
278      hrow[k] = irow;
279      colels[k] = coeff;
280      link[k] = mcstrt[jcol];
281      mcstrt[jcol] = k;
282    }
283    hincol[jcol]++;     // right?
284
285    /*
286     * Have to compute rowstat and rowduals, since we are adding the row.
287     * that satisfy complentarity slackness.
288     *
289     * We may also have to modify colstat and rcosts if bounds
290     * have been relaxed.
291     */
292    if (!colstat) {
293      // ????
294      rowduals[irow] = 0.0;
295    } else {
296      if (prob->columnIsBasic(jcol)) {
297        /* variable is already basic, so slack in this row must be basic */
298        prob->setRowStatus(irow,PrePostsolveMatrix::basic);
299       
300        rowduals[irow] = 0.0;
301        /* hence no reduced costs change */
302        /* since the column was basic, it doesn't matter if it is now
303           away from its bounds. */
304        /* the slack is basic and its reduced cost is 0 */
305      } else if ((fabs(sol[jcol] - lo0) <= ztolzb &&
306                  rcosts[jcol] >= 0) ||
307                 
308                 (fabs(sol[jcol] - up0) <= ztolzb &&
309                  rcosts[jcol] <= 0)) {
310        /* up against its bound but the rcost is still ok */
311        prob->setRowStatus(irow,PrePostsolveMatrix::basic);
312       
313        rowduals[irow] = 0.0;
314        /* hence no reduced costs change */
315      } else if (! (fabs(sol[jcol] - lo0) <= ztolzb) &&
316                 ! (fabs(sol[jcol] - up0) <= ztolzb)) {
317        /* variable not marked basic,
318         * so was at a bound in the reduced problem,
319         * but its bounds were tighter in the reduced problem,
320         * so now it has to be made basic.
321         */
322        prob->setColumnStatus(jcol,PrePostsolveMatrix::basic);
323        prob->setRowStatusUsingValue(irow);
324
325        /* choose a value for rowduals[irow] that forces rcosts[jcol] to 0.0 */
326        /* new reduced cost = 0 = old reduced cost - new dual * coeff */
327        rowduals[irow] = rcosts[jcol] / coeff;
328        rcosts[jcol] = 0.0;
329
330        /* this value is provably of the right sign for the slack */
331        /* SHOULD SHOW THIS */
332      } else {
333        /* the solution is at a bound, but the rcost is wrong */
334
335        prob->setColumnStatus(jcol,PrePostsolveMatrix::basic);
336        prob->setRowStatusUsingValue(irow);
337
338        /* choose a value for rowduals[irow] that forcesd rcosts[jcol] to 0.0 */
339        rowduals[irow] = rcosts[jcol] / coeff;
340        rcosts[jcol] = 0.0;
341
342        /* this value is provably of the right sign for the slack */
343        /*rcosts[irow] = -rowduals[irow];*/
344      }
345    }
346
347    rdone[irow] = SLACK_DOUBLETON;
348  }
349  prob->free_list_ = free_list;
350}
351
Note: See TracBrowser for help on using the repository browser.