source: branches/devel-1/PresolveFixed.cpp @ 41

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

Testing presolve on nasty burlington problem

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.5 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 "PresolveMatrix.hpp"
9#include "PresolveFixed.hpp"
10
11const char *remove_fixed_action::name() const
12{
13  return ("remove_fixed_action");
14}
15
16/*
17 * invariant:  both reps are loosely packed.
18 * coefficients of both reps remain consistent.
19 *
20 * Note that this concerns variables whose column bounds determine that
21 * they are slack; this does NOT concern singleton row constraints that
22 * determine that the relevant variable is slack.
23 *
24 * Invariant:  col and row rep are consistent
25 */
26const remove_fixed_action *remove_fixed_action::presolve(PresolveMatrix *prob,
27                                                     int *fcols,
28                                                     int nfcols,
29                                                     const PresolveAction *next)
30{
31  double *colels        = prob->colels_;
32  int *hrow             = prob->hrow_;
33  int *mcstrt           = prob->mcstrt_;
34  int *hincol           = prob->hincol_;
35
36  double *rowels        = prob->rowels_;
37  int *hcol             = prob->hcol_;
38  int *mrstrt           = prob->mrstrt_;
39  int *hinrow           = prob->hinrow_;
40  int nrows             = prob->nrows_;
41
42  double *clo   = prob->clo_;
43  double *cup   = prob->cup_;
44  double *rlo   = prob->rlo_;
45  double *rup   = prob->rup_;
46  double *acts  = prob->acts_;
47
48  double *dcost = prob->cost_;
49
50  presolvehlink *clink = prob->clink_;
51
52  action *actions       = new  action[nfcols];
53
54  for (int ckc=0; ckc<nfcols; ckc++) {
55    int j = fcols[ckc];
56
57    PRESOLVEASSERT(/*hincol[j] > 0 &&*/ cup[j] == clo[j]);
58
59    double sol = clo[j];
60    int kcs = mcstrt[j];
61    int kce = kcs + hincol[j];
62    int k;
63
64    {
65      action &f = actions[ckc];
66     
67      f.col = j;
68      f.sol = sol;
69
70      f.nincol = hincol[j];
71      f.colels = presolve_duparray(&colels[kcs], hincol[j]); // inefficient
72      f.colrows = presolve_duparray(&hrow[kcs], hincol[j]);  // inefficient
73    }
74    // the bias is updated when the empty column is removed
75    //prob->change_bias(sol * dcost[j]);
76
77    for (k=kcs; k<kce; k++) {
78      int row = hrow[k];
79      double coeff = colels[k];
80
81      rlo[row] -= sol * coeff;
82      rup[row] -= sol * coeff;
83      acts[row] -= sol * coeff;
84
85      // remove this column from all rows it occurs in in the row rep
86      presolve_delete_from_row(row, j, mrstrt, hinrow, hcol, rowels);
87
88      // mark
89      prob->addRow(row);
90      int krs = mrstrt[row];
91      int kre = krs + hinrow[row];
92      for (int k=krs; k<kre; k++) {
93        int jcol = hcol[k];
94        prob->addCol(jcol);
95      }
96      // remove the column entirely from the col rep
97      PRESOLVE_REMOVE_LINK(clink, j);
98      hincol[j] = 0;
99    }
100  }
101#if     PRESOLVE_SUMMARY
102  printf("NFIXED:  %d\n", nfcols);
103#endif
104
105#if 0
106  remove_fixed_action * nextAction =  new
107    remove_fixed_action(nfcols, actions, next);
108  delete [] actions;
109  return nextAction;
110#else
111  return (new remove_fixed_action(nfcols, actions, next));
112#endif
113}
114
115remove_fixed_action::remove_fixed_action(int nactions,
116                                         const action *actions,
117                                         const PresolveAction *next) :
118  PresolveAction(next),
119  nactions_(nactions),
120  actions_(actions)
121{
122}
123remove_fixed_action::~remove_fixed_action()
124{
125  for (int i=nactions_-1; i>=0; i--) {
126    delete[]actions_[i].colrows;
127    delete[]actions_[i].colels;
128  }
129  delete[]actions_;
130}
131
132/*
133 * Say we determined that cup - clo <= ztolzb, so we fixed sol at clo.
134 * This involved subtracting clo*coeff from ub/lb for each row the
135 * variable occurred in.
136 * Now when we put the variable back in, by construction the variable
137 * is within tolerance, the non-slacks are unchanged, and the
138 * distances of the affected slacks from their bounds should remain
139 * unchanged (ignoring roundoff errors).
140 * It may be that by adding the term back in, the affected constraints
141 * now aren't as accurate due to round-off errors; this could happen
142 * if only one summand and the slack in the original formulation were large
143 * (and naturally had opposite signs), and the new term in the constraint
144 * is about the size of the old slack, so the new slack becomes about 0.
145 * It may be that there is catastrophic cancellation in the summation,
146 * so it might not compute to 0.
147 */
148void remove_fixed_action::postsolve(PostsolveMatrix *prob) const
149{
150  const action *const actions = actions_;
151  const int nactions    = nactions_;
152
153  double *colels        = prob->colels_;
154  int *hrow             = prob->hrow_;
155  int *mcstrt           = prob->mcstrt_;
156  int *hincol           = prob->hincol_;
157  int *link             = prob->link_;
158  int ncols             = prob->ncols_;
159  int free_list         = prob->free_list_;
160
161  double *clo   = prob->clo_;
162  double *cup   = prob->cup_;
163  double *rlo   = prob->rlo_;
164  double *rup   = prob->rup_;
165
166  double *sol   = prob->sol_;
167  double *dcost = prob->cost_;
168  double *rcosts        = prob->rcosts_;
169
170  double *acts  = prob->acts_;
171  double *rowduals = prob->rowduals_;
172
173  unsigned char *colstat        = prob->colstat_;
174  unsigned char *rowstat        = prob->rowstat_;
175
176  const double maxmin   = prob->maxmin_;
177
178  char *cdone   = prob->cdone_;
179  char *rdone   = prob->rdone_;
180
181  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
182    int icol = f->col;
183    const double thesol = f->sol;
184
185    cdone[icol] = FIXED_VARIABLE;
186
187    sol[icol] = thesol;
188    clo[icol] = thesol;
189    cup[icol] = thesol;
190
191    {
192      int cs = -11111;
193      int nc = f->nincol;
194      double dj = maxmin * dcost[icol];
195
196      for (int i=0; i<nc; ++i) {
197        int row = f->colrows[i];
198        double coeff = f->colels[i];
199
200        // pop free_list
201        int k = free_list;
202        free_list = link[free_list];
203       
204        check_free_list(free_list);
205
206        // restore
207        hrow[k] = row;
208        colels[k] = coeff;
209        link[k] = cs;
210        cs = k;
211
212        PRESOLVEASSERT(rdone[row]);
213
214        if (-PRESOLVE_INF < rlo[row])
215          rlo[row] += coeff * thesol;
216        if (rup[row] < PRESOLVE_INF)
217          rup[row] += coeff * thesol;
218        acts[row] += coeff * thesol;
219
220        dj -= rowduals[row] * coeff;
221      }
222      mcstrt[icol] = cs;
223
224      rcosts[icol] = dj;
225    }
226    hincol[icol] = f->nincol;
227
228    /* the bounds in the reduced problem were tightened.
229     * that means that this variable may not have been basic
230     * because it didn't have to be,
231     * but now it may have to.
232     * no - the bounds aren't changed by this operation.
233     */
234    if (colstat)
235      prob->setColumnStatus(icol,PrePostsolveMatrix::atUpperBound);
236  }
237
238  prob->free_list_ = free_list;
239}
240
241
242const PresolveAction *remove_fixed(PresolveMatrix *prob,
243                                    const PresolveAction *next)
244{
245  int ncols     = prob->ncols_;
246  int *fcols    = new int[ncols];
247  int nfcols    = 0;
248
249  int *hincol           = prob->hincol_;
250
251  double *clo   = prob->clo_;
252  double *cup   = prob->cup_;
253
254  for (int i=0; i<ncols; i++)
255    if (hincol[i] > 0 && clo[i] == cup[i])
256      fcols[nfcols++] = i;
257
258  next = remove_fixed_action::presolve(prob, fcols, nfcols, next);
259  delete[]fcols;
260  return (next);
261}
262
263
264
265const char *make_fixed_action::name() const
266{
267  return ("make_fixed_action");
268}
269
270const PresolveAction *make_fixed_action::presolve(PresolveMatrix *prob,
271                                                     int *fcols,
272                                                     int nfcols,
273                                                   bool fix_to_lower,
274                                                   const PresolveAction *next)
275{
276  double *clo   = prob->clo_;
277  double *cup   = prob->cup_;
278
279  action *actions       = new  action[nfcols];
280
281  for (int ckc=0; ckc<nfcols; ckc++) {
282    int j = fcols[ckc];
283
284    action &f = actions[ckc];
285     
286    if (fix_to_lower) {
287      f.bound = cup[j];
288      cup[j] = clo[j];
289    } else {
290      f.bound = clo[j];
291      clo[j] = cup[j];
292    }
293  }
294
295  // this is unusual in that the make_fixed_action transform
296  // contains within it a remove_fixed_action transform
297  // bad idea?
298  return (new make_fixed_action(nfcols, actions, fix_to_lower,
299                                remove_fixed_action::presolve(prob,
300                                                              fcols, nfcols,
301                                                              0),
302                                next));
303}
304
305void make_fixed_action::postsolve(PostsolveMatrix *prob) const
306{
307  const action *const actions = actions_;
308  const int nactions    = nactions_;
309  const bool fix_to_lower       = fix_to_lower_;
310
311  double *clo   = prob->clo_;
312  double *cup   = prob->cup_;
313
314  faction_->postsolve(prob);
315
316  for (int cnt = faction_->nactions_-1; cnt>=0; cnt--) {
317    const action *f = &actions[cnt];
318    const remove_fixed_action::action *f1 = &faction_->actions_[cnt];
319
320    int icol = f1->col;
321    if (fix_to_lower) {
322      cup[icol] = f->bound;
323    } else {
324      clo[icol] = f->bound;
325    }
326  }
327}
328
329
330const PresolveAction *make_fixed(PresolveMatrix *prob,
331                                  const PresolveAction *next)
332{
333  int ncols     = prob->ncols_;
334  int *fcols    = new int[ncols];
335  int nfcols    = 0;
336
337  int *hincol           = prob->hincol_;
338
339  double *clo   = prob->clo_;
340  double *cup   = prob->cup_;
341
342  for (int i=0; i<ncols; i++)
343    if (hincol[i] > 0 && fabs(cup[i] - clo[i]) < ZTOLDP) 
344      fcols[nfcols++] = i;
345
346  next = make_fixed_action::presolve(prob, fcols, nfcols,
347                                     true, // arbitrary
348                                     next);
349  delete[]fcols;
350  return (next);
351}
352
353
Note: See TracBrowser for help on using the repository browser.