source: trunk/PresolveEmpty.cpp @ 173

Last change on this file since 173 was 173, checked in by forrest, 16 years ago

Improvements to presolve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.5 KB
RevLine 
[50]1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <stdio.h>
5#include <math.h>
6
7#include "PresolveMatrix.hpp"
8
9#include "PresolveEmpty.hpp"
10#include "ClpMessage.hpp"
11
12
13const PresolveAction *drop_empty_cols_action::presolve(PresolveMatrix *prob,
14                                                        int *ecols,
15                                                        int necols,
16                                                        const PresolveAction *next)
17{
18  int ncols             = prob->ncols_;
19  CoinBigIndex *mcstrt          = prob->mcstrt_;
20  int *hincol           = prob->hincol_;
21  //  int *hcol         = prob->hcol_;
22
23  // We know we are not going to need row copy again
24  //presolvehlink *clink        = prob->clink_;
25
26  double *clo           = prob->clo_;
27  double *cup           = prob->cup_;
28  double *dcost         = prob->cost_;
29
[173]30  const double ztoldj   = prob->ztoldj_;
[50]31  //int *mrstrt         = prob->mrstrt_;
32  //int *hinrow         = prob->hinrow_;
33  //int *hrow           = prob->hrow_;
34
35  char * integerType     = prob->integerType_;
36  int * originalColumn  = prob->originalColumn_;
37
38  const double maxmin   = prob->maxmin_;
39
40  double * sol = prob->sol_;
41  unsigned char * colstat = prob->colstat_;
42
43  action *actions       = new action[necols];
44  int * colmapping = new int [ncols];
45
46  memset(colmapping,0,ncols*sizeof(int));
47  int i;
48  for (i=necols-1; i>=0; i--) {
49    int jcol = ecols[i];
50    colmapping[jcol]=-1;
51    action &e = actions[i];
52
53    e.jcol      = jcol;
54    e.clo       = clo[jcol];
55    e.cup       = cup[jcol];
56    e.cost      = dcost[jcol];
57
58    // there are no more constraints on this variable,
59    // so we had better be able to compute the answer now
[173]60    if (fabs(dcost[jcol])<ztoldj)
61      dcost[jcol]=0.0;
[50]62    if (dcost[jcol] * maxmin == 0.0) {
63      // hopefully, we can make this non-basic
64      // what does OSL currently do in this case???
65      e.sol = (-PRESOLVE_INF < clo[jcol]
66               ? clo[jcol]
67               : cup[jcol] < PRESOLVE_INF
68               ? cup[jcol]
69               : 0.0);
70    } else if (dcost[jcol] * maxmin > 0.0) {
71      if (-PRESOLVE_INF < clo[jcol])
72        e.sol = clo[jcol];
73      else {
74          prob->messageHandler()->message(CLP_PRESOLVE_COLUMNBOUNDB,
75                                             prob->messages())
76                                               <<jcol
77                                               <<CoinMessageEol;
78        prob->status_ |= 2;
79        break;
80      }
81    } else {
82      if (cup[jcol] < PRESOLVE_INF)
83        e.sol = cup[jcol];
84      else {
85          prob->messageHandler()->message(CLP_PRESOLVE_COLUMNBOUNDA,
86                                             prob->messages())
87                                               <<jcol
88                                               <<CoinMessageEol;
89        prob->status_ |= 2;
90        break;
91      }
92    }
93       
94#if     DEBUG_PRESOLVE
95    if (e.sol * dcost[jcol]) {
96      //printf("\a>>>NON-ZERO COST FOR EMPTY COL %d:  %g\n", jcol, dcost[jcol]);
97    }
98#endif
99    prob->change_bias(e.sol * dcost[jcol]);
100
101
102  }
103  int ncols2=0;
104
105  // now move remaining ones down
106  for (i=0;i<ncols;i++) {
107    if (!colmapping[i]) {
108      mcstrt[ncols2] = mcstrt[i];
109      hincol[ncols2] = hincol[i];
110   
111      clo[ncols2]   = clo[i];
112      cup[ncols2]   = cup[i];
113
114      dcost[ncols2] = dcost[i];
115      if (sol) {
116        sol[ncols2] = sol[i];
117        colstat[ncols2] = colstat[i];
118      }
119
120      integerType[ncols2] = integerType[i];
121      originalColumn[ncols2] = originalColumn[i];
122      ncols2++;
123    }
124  }
125 
126  delete [] colmapping;
127  prob->ncols_ = ncols2;
128
129  return (new drop_empty_cols_action(necols, actions, next));
130}
131
132
133const PresolveAction *drop_empty_cols_action::presolve(PresolveMatrix *prob,
134                                                          const PresolveAction *next)
135{
136  const int *hincol     = prob->hincol_;
137  //  const double *clo = prob->clo_;
138  //  const double *cup = prob->cup_;
139  int ncols             = prob->ncols_;
140  int i;
141  int nempty            = 0;
142  int * empty = new int [ncols];
143
144  // count empty cols
145  for (i=0; i<ncols; i++)
146    if (hincol[i] == 0) {
147#if     DEBUG_PRESOLVE
148      if (nempty==0)
149        printf("UNUSED COLS:  ");
150      printf("%d ", i);
151#endif
152      empty[nempty++] = i;
153    }
154
155  if (nempty) {
156#if     DEBUG_PRESOLVE
157      printf("\ndropped %d cols\n", nempty);
158#endif
159    next = drop_empty_cols_action::presolve(prob, empty, nempty, next);
160  }
161  delete [] empty;
162  return (next);
163}
164
165
166void drop_empty_cols_action::postsolve(PostsolveMatrix *prob) const
167{
168  const int nactions    = nactions_;
169  const action *const actions = actions_;
170
171  int ncols             = prob->ncols_;
172  //  CoinBigIndex nelems               = prob->nelems_;
173
174  CoinBigIndex *mcstrt  = prob->mcstrt_;
175  int *hincol   = prob->hincol_;
176  //  int *hrow = prob->hrow_;
177
178  double *clo   = prob->clo_;
179  double *cup   = prob->cup_;
180
181  double *sol   = prob->sol_;
182  double *cost  = prob->cost_;
183  double *rcosts        = prob->rcosts_;
184  unsigned char *colstat        = prob->colstat_;
185  const double maxmin = prob->maxmin_;
186
187  int ncols2 = ncols+nactions;
188  int * colmapping = new int [ncols2];
189
190  memset(colmapping,0,ncols2*sizeof(int));
191  char *cdone   = prob->cdone_;
[56]192  int action_i;
193  for (action_i = 0; action_i < nactions; action_i++) {
[50]194    const action *e = &actions[action_i];
195    int jcol = e->jcol;
196    colmapping[jcol]=-1;
197  }
198
199  int i;
200
201  // now move remaining ones up
202  for (i=ncols2-1;i>=0;i--) {
203    if (!colmapping[i]) {
204      ncols--;
205      mcstrt[i] = mcstrt[ncols];
206      hincol[i] = hincol[ncols];
207   
208      clo[i]   = clo[ncols];
209      cup[i]   = cup[ncols];
210
211      cost[i] = cost[ncols];
212
213      sol[i] = sol[ncols];
214     
215      rcosts[i] = rcosts[ncols];
216     
217      if (colstat) 
218        colstat[i] = colstat[ncols];
219      cdone[i]   = cdone[ncols];
220    }
221  }
222  assert (!ncols);
223 
224  delete [] colmapping;
225
[56]226  for (action_i = 0; action_i < nactions; action_i++) {
[50]227    const action *e = &actions[action_i];
228    int jcol = e->jcol;
229   
230    // now recreate jcol
231    clo[jcol] = e->clo;
232    cup[jcol] = e->cup;
233    sol[jcol] = e->sol;
234    cost[jcol] = e->cost;
235
236    rcosts[jcol] = maxmin*cost[jcol];
237
238    hincol[jcol] = 0;
239
240    cdone[jcol] = DROP_COL;
241    if (colstat) 
242      prob->setColumnStatusUsingValue(jcol);
243  }
244
245  prob->ncols_ += nactions;
246}
247
248
249
250const PresolveAction *drop_empty_rows_action::presolve(PresolveMatrix *prob,
251                                       const PresolveAction *next)
252{
253  int ncols     = prob->ncols_;
254  CoinBigIndex *mcstrt  = prob->mcstrt_;
255  int *hincol   = prob->hincol_;
256  int *hrow     = prob->hrow_;
257
258  int nrows     = prob->nrows_;
259  // This is done after row copy needed
260  //int *mrstrt = prob->mrstrt_;
261  int *hinrow   = prob->hinrow_;
262  //int *hcol   = prob->hcol_;
263 
264  double *rlo   = prob->rlo_;
265  double *rup   = prob->rup_;
266
267  unsigned char *rowstat        = prob->rowstat_;
268  double *acts  = prob->acts_;
269
270  //presolvehlink *rlink = prob->rlink_;
271 
272
273  int i;
274  int nactions = 0;
275  for (i=0; i<nrows; i++)
276    if (hinrow[i] == 0)
277      nactions++;
278
279  if (nactions == 0)
280    return next;
281  else {
282    action *actions     = new action[nactions];
283    int * rowmapping = new int [nrows];
284
285    nactions = 0;
286    int nrows2=0;
287    for (i=0; i<nrows; i++) {
288      if (hinrow[i] == 0) {
289        action &e = actions[nactions];
290        nactions++;
291
292#if     DEBUG_PRESOLVE
293        if (nactions==1)
294          printf("unused rows:  ");
295
296        printf("%d ", i);
297#endif
298        if (rlo[i] > 0.0 || rup[i] < 0.0) {
[94]299          if (rlo[i]<=prob->feasibilityTolerance_ &&
300              rup[i]>=-prob->feasibilityTolerance_) {
[104]301            rlo[i]=0.0;
302            rup[i]=0.0;
[50]303          } else {
304            prob->status_|= 1;
305          prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
306                                             prob->messages())
307                                               <<i
308                                               <<rlo[i]
309                                               <<rup[i]
310                                               <<CoinMessageEol;
311            break;
312          }
313        }
314        e.row   = i;
315        e.rlo   = rlo[i];
316        e.rup   = rup[i];
317        rowmapping[i]=-1;
318
319      } else {
320        // move down - we want to preserve order
321        rlo[nrows2]=rlo[i];
322        rup[nrows2]=rup[i];
323        if (acts) {
324          acts[nrows2]=acts[i];
325          rowstat[nrows2]=rowstat[i];
326        }
327        rowmapping[i]=nrows2++;
328      }
329    }
330
331    // remap matrix
332    for (i=0;i<ncols;i++) {
333      int j;
334      for (j=mcstrt[i];j<mcstrt[i]+hincol[i];j++) 
335        hrow[j] = rowmapping[hrow[j]];
336    }
337    delete [] rowmapping;
338
339    prob->nrows_ = nrows2;
340
341#if     DEBUG_PRESOLVE
342    if (nactions)
343      printf("\ndropped %d rows\n", nactions);
344#endif
345    return (new drop_empty_rows_action(nactions, actions, next));
346  }
347}
348
349void drop_empty_rows_action::postsolve(PostsolveMatrix *prob) const
350{
351  const int nactions    = nactions_;
352  const action *const actions = actions_;
353
354  int ncols     = prob->ncols_;
355  CoinBigIndex *mcstrt  = prob->mcstrt_;
356  int *hincol   = prob->hincol_;
357
358  int *hrow     = prob->hrow_;
359
360  double *rlo   = prob->rlo_;
361  double *rup   = prob->rup_;
362  unsigned char *rowstat        = prob->rowstat_;
363  double *rowduals = prob->rowduals_;
364  double *acts  = prob->acts_;
365  char *rdone   = prob->rdone_;
366
367  int nrows0    = prob->nrows0_;
368  int nrows     = prob->nrows_;
369
370  int * rowmapping = new int [nrows0];
371
372  int i, action_i;
373  for (i=0; i<nrows0; i++)
374    rowmapping[i] = 0;
375 
376  for (action_i = 0; action_i<nactions; action_i++) {
377    const action *e = &actions[action_i];
378    int hole = e->row;
379    rowmapping[hole]=-1;
380  }
381
382  // move data
383  for (i=nrows0-1; i>=0; i--) {
384    if (!rowmapping[i]) {
385      // not a hole
386      nrows--;
387      rlo[i]=rlo[nrows];
388      rup[i]=rup[nrows];
389      acts[i]=acts[nrows];
390      rowduals[i]=rowduals[nrows];
391      if (rowstat)
392        rowstat[i] = rowstat[nrows];
393    }
394  }
395  assert (!nrows);
396  // set up mapping for matrix
397  for (i=0;i<nrows0;i++) {
398    if (!rowmapping[i])
399      rowmapping[nrows++]=i;
400  }
401
402  for (int j=0; j<ncols; j++) {
403    CoinBigIndex start = mcstrt[j];
404    CoinBigIndex end   = start + hincol[j];
405   
406    for (CoinBigIndex k=start; k<end; ++k) {
407      hrow[k] = rowmapping[hrow[k]];
408    }
409  }
410
411
412  delete [] rowmapping;
413
414  for (action_i = 0; action_i < nactions; action_i++) {
415    const action *e = &actions[action_i];
416    int irow = e->row;
417
418    // Now recreate irow
419    rlo[irow] = e->rlo;
420    rup[irow] = e->rup;
421
422    if (rowstat)
423      prob->setRowStatus(irow,PrePostsolveMatrix::basic);
424    rowduals[irow] = 0.0;       // ???
425    acts[irow] = 0.0;
426
427    // hinrow is not maintained during postsolve
428    //hinrow[irow] = 0;
429
430    rdone[irow] = DROP_ROW;
431  }
432
433  prob->nrows_ = prob->nrows_+nactions;
434}
435
Note: See TracBrowser for help on using the repository browser.