source: trunk/PresolveEmpty.cpp @ 104

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

Changes for free variables and unbounded

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