source: trunk/PresolveImpliedFree.cpp @ 57

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

Yet more changes for Visual C++

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 KB
Line 
1#include <stdio.h>
2#include <math.h>
3
4#include "PresolveMatrix.hpp"
5#include "PresolveSubst.hpp"
6#include "PresolveIsolated.hpp"
7#include "PresolveImpliedFree.hpp"
8#include "ClpMessage.hpp"
9
10
11
12// If there is a row with a singleton column such that no matter what
13// the values of the other variables are, the constraint forces the singleton
14// column to have a feasible value, then we can drop the column and row,
15// since we just compute the value of the column from the row in postsolve.
16// This seems vaguely similar to the case of a useless constraint, but it
17// is different.  For example, if the singleton column is already free,
18// then this operation will eliminate it, but the constraint is not useless
19// (assuming the constraint is not trivial), since the variables do not imply an
20// upper or lower bound.
21//
22// If the column is not singleton, we can still do something similar if the
23// constraint is an equality constraint.  In that case, we substitute away
24// the variable in the other constraints it appears in.  This introduces
25// new coefficients, but the total number of coefficients never increases
26// if the column has only two constraints, and may not increase much even
27// if there are more.
28//
29// There is nothing to prevent us from substituting away a variable
30// in an equality from the other constraints it appears in, but since
31// that causes fill-in, it wouldn't make sense unless we could then
32// drop the equality itself.  We can't do that if the bounds on the
33// variable in equation aren't implied by the equality.
34// Another way of thinking of this is that there is nothing special
35// about an equality; just like one can't always drop an inequality constraint
36// with a column singleton, one can't always drop an equality.
37//
38// It is possible for two singleton columns to be in the same row.
39// In that case, the other one will become empty.  If it's bounds and
40// costs aren't just right, this signals an unbounded problem.
41// We don't need to check that specially here.
42//
43// invariant:  loosely packed
44const PresolveAction *implied_free_action::presolve(PresolveMatrix *prob,
45                                                     const PresolveAction *next,
46                                                    int & fill_level)
47{
48  double *colels        = prob->colels_;
49  int *hrow     = prob->hrow_;
50  const CoinBigIndex *mcstrt    = prob->mcstrt_;
51  int *hincol   = prob->hincol_;
52  const int ncols       = prob->ncols_;
53
54  const double *clo     = prob->clo_;
55  const double *cup     = prob->cup_;
56
57  const double *rowels  = prob->rowels_;
58  const int *hcol       = prob->hcol_;
59  const CoinBigIndex *mrstrt    = prob->mrstrt_;
60  int *hinrow   = prob->hinrow_;
61  //  const int nrows   = prob->nrows_;
62
63  /*const*/ double *rlo = prob->rlo_;
64  /*const*/ double *rup = prob->rup_;
65
66  double *cost          = prob->cost_;
67
68  presolvehlink *rlink = prob->rlink_;
69  presolvehlink *clink = prob->clink_;
70
71  const char *integerType = prob->integerType_;
72
73  const double tol = prob->feasibilityTolerance_;
74
75  //  int nbounds = 0;
76
77  action *actions       = new action [ncols];
78  int nactions = 0;
79
80  char *implied_free = new char[ncols];
81  bzero(implied_free, ncols*sizeof(char));
82
83  double *ilbound = new double[ncols];
84  double *iubound = new double[ncols];
85
86  double *tclo = new double[ncols];
87  double *tcup = new double[ncols];
88
89#if     PRESOLVE_TRY_TIGHTEN
90  for (int j=0; j<ncols; j++) {
91    ilbound[j] = -PRESOLVE_INF;
92    iubound[j] =  PRESOLVE_INF;
93    tclo[j] = clo[j];
94    tcup[j] = cup[j];
95  }
96
97  {
98    int ntightened;
99    do {
100      implied_bounds1(rowels, mrstrt, hrow, hinrow, clo, cup, hcol, ncols,
101                      rlo, rup, integerType, nrows,
102                      ilbound, iubound);
103
104      ntightened = 0;
105      for (int j=0; j<ncols; j++) {
106        if (tclo[j] < ilbound[j]) {
107          tclo[j] = ilbound[j];
108          ntightened++;
109        }
110        if (tcup[j] > iubound[j]) {
111          tcup[j] = iubound[j];
112          ntightened++;
113        }
114      }
115      printf("NTIGHT: %d\n", ntightened);
116    } while (ntightened);
117  }
118#endif
119
120  int numberLook = prob->numberColsToDo_;
121  int iLook;
122  int * look = prob->colsToDo_;
123  int * look2 = NULL;
124  // if gone from 2 to 3 look at all
125  if (fill_level<0) {
126    look2 = new int[ncols];
127    look=look2;
128    for (iLook=0;iLook<ncols;iLook++) 
129      look[iLook]=iLook;
130    numberLook=ncols;
131 }
132
133  for (iLook=0;iLook<numberLook;iLook++) {
134    int j=look[iLook];
135    if ((hincol[j] >= 1 && hincol[j] <= 3) &&
136        !integerType[j]) {
137      CoinBigIndex kcs = mcstrt[j];
138      CoinBigIndex kce = kcs + hincol[j];
139
140      for (CoinBigIndex k=kcs; k<kce; ++k) {
141        int row = hrow[k];
142        double coeffj = colels[k];
143
144        // if its row is an equality constraint...
145        if (hinrow[row] > 1 &&  // don't bother with singleton rows
146
147            // either this is a singleton col,
148            // or this particular row is an equality
149            (hincol[j] == 1 || fabs(rlo[row] - rup[row]) < tol) &&
150
151            fabs(coeffj) > ZTOLDP) {
152
153          CoinBigIndex krs = mrstrt[row];
154          CoinBigIndex kre = krs + hinrow[row];
155
156          double maxup, maxdown, ilow, iup;
157          implied_bounds(rowels, clo, cup, hcol,
158                         krs, kre,
159                         &maxup, &maxdown,
160                         j, rlo[row], rup[row], &ilow, &iup);
161
162#if     PRESOLVE_TRY_TIGHTEN
163          if ((clo[j] <= ilbound[j] && iubound[j] <= cup[j]) &&
164              ! (clo[j] <= ilow && iup <= cup[j]))
165            printf("TIGHTER:  %6d %6d\n", row, j);
166#endif
167
168          if (maxup < PRESOLVE_INF && maxup + tol < rlo[row]) {
169            /* there is an upper bound and it can't be reached */
170            prob->status_|= 1;
171            prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
172                                             prob->messages())
173                                               <<row
174                                               <<rlo[row]
175                                               <<rup[row]
176                                               <<CoinMessageEol;
177            break;
178          } else if (-PRESOLVE_INF < maxdown && rup[row] < maxdown - tol) {
179            /* there is a lower bound and it can't be reached */
180            prob->status_|= 1;
181            prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
182                                             prob->messages())
183                                               <<row
184                                               <<rlo[row]
185                                               <<rup[row]
186                                               <<CoinMessageEol;
187            break;
188          } else if (clo[j] <= ilow && iup <= cup[j]) {
189
190            // both column bounds implied by the constraints of the problem
191            implied_free[j] = hincol[j];
192            break;
193          }
194        }
195      }
196    }
197  }
198  // implied_free[j] == hincol[j] && hincol[j] > 0 ==> j is implied free
199
200#if 0
201  // DEBUG
202  static int nfree = 0;
203  static int maxfree = atoi(getenv("MAXFREE"));
204#endif
205
206  int isolated_row = -1;
207
208  // first pick off the easy ones
209  // note that this will only deal with columns that were originally
210  // singleton; it will not deal with doubleton columns that become
211  // singletons as a result of dropping rows.
212  for (iLook=0;iLook<numberLook;iLook++) {
213    int j=look[iLook];
214    if (hincol[j] == 1 && implied_free[j] == 1) {
215      CoinBigIndex kcs = mcstrt[j];
216      int row = hrow[kcs];
217      double coeffj = colels[kcs];
218
219      CoinBigIndex krs = mrstrt[row];
220      CoinBigIndex kre = krs + hinrow[row];
221
222#if 0
223      if (nfree >= maxfree)
224        continue;
225      nfree++;
226#endif
227
228      // isolated rows are weird
229      {
230        int n = 0;
231        for (CoinBigIndex k=krs; k<kre; ++k)
232          n += hincol[hcol[k]];
233        if (n==hinrow[row]) {
234          isolated_row = row;
235          break;
236        }
237      }
238
239      const bool nonzero_cost = (cost[j] != 0.0&&fabs(rup[row]-rlo[row])<=tol);
240
241      double *save_costs = nonzero_cost ? new double[hinrow[row]] : NULL;
242
243      {
244        action *s = &actions[nactions++];
245             
246        s->row = row;
247        s->col = j;
248
249        s->clo = clo[j];
250        s->cup = cup[j];
251        s->rlo = rlo[row];
252        s->rup = rup[row];
253
254        s->ninrow = hinrow[row];
255        s->rowels = presolve_duparray(&rowels[krs], hinrow[row]);
256        s->rowcols = presolve_duparray(&hcol[krs], hinrow[row]);
257        s->costs = save_costs;
258      }
259
260      if (nonzero_cost) {
261        double rhs = rlo[row];
262        double costj = cost[j];
263
264#if     DEBUG_PRESOLVE
265        printf("FREE COSTS:  %g  ", costj);
266#endif
267        for (CoinBigIndex k=krs; k<kre; k++) {
268          int jcol = hcol[k];
269          save_costs[k-krs] = cost[jcol];
270
271          if (jcol != j) {
272            double coeff = rowels[k];
273
274#if     DEBUG_PRESOLVE
275            printf("%g %g   ", cost[jcol], coeff/coeffj);
276#endif
277            /*
278             * Similar to eliminating doubleton:
279             *   cost1 x = cost1 (c - b y) / a = (c cost1)/a - (b cost1)/a
280             *   cost[icoly] += cost[icolx] * (-coeff2 / coeff1);
281             */
282            cost[jcol] += costj * (-coeff / coeffj);
283          }
284        }
285#if     DEBUG_PRESOLVE
286        printf("\n");
287
288        /* similar to doubleton */
289        printf("BIAS??????? %g %g %g %g\n",
290               costj * rhs / coeffj,
291               costj, rhs, coeffj);
292#endif
293        prob->change_bias(costj * rhs / coeffj);
294        // ??
295        cost[j] = 0.0;
296      }
297
298      /* remove the row from the columns in the row */
299      for (CoinBigIndex k=krs; k<kre; k++) {
300        int jcol=hcol[k];
301        prob->addCol(jcol);
302        presolve_delete_from_row(jcol, row, mcstrt, hincol, hrow, colels);
303      }
304      PRESOLVE_REMOVE_LINK(rlink, row);
305      hinrow[row] = 0;
306
307      // just to make things squeeky
308      rlo[row] = 0.0;
309      rup[row] = 0.0;
310
311      PRESOLVE_REMOVE_LINK(clink, j);
312      hincol[j] = 0;
313
314      implied_free[j] = 0;      // probably unnecessary
315    }
316  }
317
318  delete [] look2;
319  if (nactions) {
320#if     PRESOLVE_SUMMARY
321    printf("NIMPLIED FREE:  %d\n", nactions);
322#endif
323    next = new implied_free_action(nactions, copyOfArray(actions,nactions), next);
324  }
325  deleteAction(actions);
326
327  delete[]ilbound;
328  delete[]iubound;
329  delete[]tclo;
330  delete[]tcup;
331
332  if (isolated_row != -1) {
333    const PresolveAction *nextX = isolated_constraint_action::presolve(prob, 
334                                                isolated_row, next);
335    if (nextX)
336      next = nextX; // may fail
337  }
338
339  // try more complex ones
340  if (fill_level)
341    next = subst_constraint_action::presolve(prob, implied_free, next,fill_level);
342  delete[]implied_free;
343
344  return (next);
345}
346
347
348
349const char *implied_free_action::name() const
350{
351  return ("implied_free_action");
352}
353
354void implied_free_action::postsolve(PostsolveMatrix *prob) const
355{
356  const action *const actions = actions_;
357  const int nactions = nactions_;
358
359  double *colels        = prob->colels_;
360  int *hrow             = prob->hrow_;
361  CoinBigIndex *mcstrt          = prob->mcstrt_;
362  int *hincol           = prob->hincol_;
363  int *link             = prob->link_;
364
365  double *clo   = prob->clo_;
366  double *cup   = prob->cup_;
367
368  double *rlo   = prob->rlo_;
369  double *rup   = prob->rup_;
370
371  double *sol   = prob->sol_;
372
373  double *rcosts        = prob->rcosts_;
374  double *dcost         = prob->cost_;
375
376  double *acts  = prob->acts_;
377  double *rowduals = prob->rowduals_;
378
379  //  const double ztoldj       = prob->ztoldj_;
380  //  const double ztolzb       = prob->ztolzb_;
381
382  const double maxmin   = prob->maxmin_;
383
384  char *cdone   = prob->cdone_;
385  char *rdone   = prob->rdone_;
386  CoinBigIndex free_list = prob->free_list_;
387
388  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
389
390    int irow = f->row;
391    int icol = f->col;
392         
393    int ninrow = f->ninrow;
394    const double *rowels = f->rowels;
395    const int *rowcols = f->rowcols;
396    const double *save_costs = f->costs;
397
398    // put back coefficients in the row
399    // this includes recreating the singleton column
400    {
401      for (int k = 0; k<ninrow; k++) {
402        int jcol = rowcols[k];
403        double coeff = rowels[k];
404
405        if (save_costs)
406          dcost[jcol] = save_costs[k];
407
408        {
409          CoinBigIndex kk = free_list;
410          free_list = link[free_list];
411
412          check_free_list(free_list);
413
414          link[kk] = mcstrt[jcol];
415          mcstrt[jcol] = kk;
416          colels[kk] = coeff;
417          hrow[kk] = irow;
418        }
419
420        if (jcol == icol) {
421          // initialize the singleton column
422          hincol[jcol] = 1;
423          clo[icol] = f->clo;
424          cup[icol] = f->cup;
425
426          cdone[icol] = IMPLIED_FREE;
427        } else {
428          hincol[jcol]++;
429        }
430      }
431      rdone[irow] = IMPLIED_FREE;
432
433      rlo[irow] = f->rlo;
434      rup[irow] = f->rup;
435    }
436    delete [] save_costs;
437    // coeff has now been initialized
438
439    // compute solution
440    {
441      double act = 0.0;
442      double coeff = 0.0;
443
444      for (int k = 0; k<ninrow; k++)
445        if (rowcols[k] == icol)
446          coeff = rowels[k];
447        else {
448          int jcol = rowcols[k];
449          PRESOLVE_STMT(CoinBigIndex kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link));
450          act += rowels[k] * sol[jcol];
451        }
452           
453      PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
454      // choose rowdual to make this col basic
455      rowduals[irow] = maxmin*dcost[icol] / coeff;
456      rcosts[icol] = 0.0;
457      if ((rlo[irow] < rup[irow] && rowduals[irow] < 0.0)
458          || rlo[irow]< -1.0e20) {
459        if (rlo[irow]<-1.0e20&&rowduals[irow]>=0.0)
460          printf("IMP %g %g %g\n",rlo[irow],rup[irow],rowduals[irow]);
461        sol[icol] = (rup[irow] - act) / coeff;
462        acts[irow] = rup[irow];
463      } else {
464        sol[icol] = (rlo[irow] - act) / coeff;
465        acts[irow] = rlo[irow];
466      }
467
468
469      prob->setRowStatus(irow,PrePostsolveMatrix::atLowerBound);
470      prob->setColumnStatus(icol,PrePostsolveMatrix::basic);
471    }
472  }
473  prob->free_list_ = free_list;
474}
Note: See TracBrowser for help on using the repository browser.