source: branches/devel-1/PresolveForcing.cpp @ 33

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

Presolve nearly ready

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.9 KB
Line 
1#include <stdio.h>
2#include <math.h>
3#include <strings.h>
4
5#include "PresolveMatrix.hpp"
6#include "PresolveEmpty.hpp"    // for DROP_COL/DROP_ROW
7#include "PresolveFixed.hpp"
8#include "PresolveSubst.hpp"
9#include "PresolveUseless.hpp"
10#include "PresolveForcing.hpp"
11#include "ClpMessage.hpp"
12
13#if 0
14inline double max(double x, double y)
15{
16  return (x < y) ? y : x;
17}
18
19inline double min(double x, double y)
20{
21  return (x < y) ? x : y;
22}
23#endif
24static void implied_bounds1(PresolveMatrix * prob,const double *rowels,
25                                const int *mrstrt,
26                                const int *hrow,
27                                const int *hinrow,
28                                const double *clo, const double *cup,
29                                const int *hcol,
30                                int ncols,
31                                const double *rlo, const double *rup,
32                                const char *integerType,
33                                int nrows,
34                                double *ilbound, double *iubound);
35
36
37/*static*/ void implied_bounds(const double *els,
38                           const double *clo, const double *cup,
39                           const int *hcol,
40                           int krs, int kre,
41                           double *maxupp, double *maxdownp,
42                           int jcol,
43                           double rlo, double rup,
44                           double *iclb, double *icub)
45{
46  if (rlo<=-PRESOLVE_INF&&rup>=PRESOLVE_INF) {
47    *iclb = -PRESOLVE_INF;
48    *icub =  PRESOLVE_INF;
49    return;
50  }
51  bool posinf = false;
52  bool neginf = false;
53  double maxup = 0.0;
54  double maxdown = 0.0;
55
56  int jcolk = -1;
57
58  // compute sum of all bounds except for jcol
59  int kk;
60  for (kk=krs; kk<kre; kk++) {
61    if (hcol[kk] == jcol)
62      jcolk = kk;
63
64    // swap jcol with hcol[kre-1];
65    // that is, consider jcol last
66    // this assumes that jcol occurs in this row
67    int k = (hcol[kk] == jcol
68             ? kre-1
69             : kk == kre-1
70             ? jcolk
71             : kk);
72
73    PRESOLVEASSERT(k != -1);    // i.e. jcol had better be in the row
74
75    int col = hcol[k];
76    double coeff = els[k];
77    double lb = clo[col];
78    double ub = cup[col];
79
80    // quick!  compute the implied col bounds before maxup/maxdown are changed
81    if (kk == kre-1) {
82      PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
83
84      double ilb = (rlo - maxup) / coeff;
85      bool finite_ilb = (-PRESOLVE_INF < rlo && !posinf && PRESOLVEFINITE(maxup));
86
87      double iub = (rup - maxdown) / coeff;
88      bool finite_iub = ( rup < PRESOLVE_INF && !neginf && PRESOLVEFINITE(maxdown));
89
90      if (coeff > 0.0) {
91        *iclb = (finite_ilb ? ilb : -PRESOLVE_INF);
92        *icub = (finite_iub ? iub :  PRESOLVE_INF);
93      } else {
94        *iclb = (finite_iub ? iub : -PRESOLVE_INF);
95        *icub = (finite_ilb ? ilb :  PRESOLVE_INF);
96      }
97    }
98
99    if (coeff > 0.0) {
100      if (PRESOLVE_INF <= ub) {
101        posinf = true;
102        if (neginf)
103          break;        // pointless
104      } else
105        maxup += ub * coeff;
106
107      if (lb <= -PRESOLVE_INF) {
108        neginf = true;
109        if (posinf)
110          break;        // pointless
111      } else
112        maxdown += lb * coeff;
113    } else {
114      if (PRESOLVE_INF <= ub) {
115        neginf = true;
116        if (posinf)
117          break;        // pointless
118      } else
119        maxdown += ub * coeff;
120
121      if (lb <= -PRESOLVE_INF) {
122        posinf = true;
123        if (neginf)
124          break;        // pointless
125      } else
126        maxup += lb * coeff;
127    }
128  }
129
130  // If we broke from the loop, then the bounds are infinite.
131  // However, since we put the column whose implied bounds we want
132  // to know at the end, and it doesn't matter if its own bounds
133  // are infinite, don't worry about breaking at the last iteration.
134  if (kk<kre-1) {
135    *iclb = -PRESOLVE_INF;
136    *icub =  PRESOLVE_INF;
137  } else
138    PRESOLVEASSERT(jcolk != -1);
139
140  // store row bounds
141  *maxupp   = (posinf) ?  PRESOLVE_INF : maxup;
142  *maxdownp = (neginf) ? -PRESOLVE_INF : maxdown;
143}
144
145static void implied_row_bounds(const double *els,
146                               const double *clo, const double *cup,
147                               const int *hcol,
148                               int krs, int kre,
149                               double *maxupp, double *maxdownp)
150{
151  int jcol;
152  double iclb, icub;
153
154  implied_bounds(els, clo, cup, hcol, krs, kre, maxupp, maxdownp,
155                 hcol[krs], 0.0, 0.0, &iclb, &icub);
156}
157
158const char *forcing_constraint_action::name() const
159{
160  return ("forcing_constraint_action");
161}
162
163static bool some_col_was_fixed(const int *hcol, int krs, int kre,
164                               const double *clo, 
165                               const double *cup)
166{
167  int k;
168  for (k=krs; k<kre; k++) {
169    int jcol = hcol[k];
170    if (clo[jcol] == cup[jcol])
171      break;
172  }
173  return (k<kre);
174}
175
176
177//
178// It may be the case that the variable bounds are such that no matter what
179// feasible value they take, the constraint cannot be violated;
180// in this case, we obviously don't need to take it into account, and
181// we just drop it as a USELESS constraint.
182//
183// On the other hand, it may be that the only way to satisfy a constraint
184// is to jam all the constraint variables to one of their bounds;
185// in this case, these variables are essentially fixed, which
186// we do with a FORCING constraint.
187// For now, this just tightens the bounds; subsequently the fixed variables
188// will be removed, then the row will be dropped.
189//
190// Since both operations use similar information (the implied row bounds),
191// we combine them into one presolve routine.
192//
193// I claim that these checks could be performed in parallel,
194// that is, the tests could be carried out for all rows in parallel,
195// and then the rows deleted and columns tightened afterward.
196// Obviously, this is true for useless rows.
197// The potential problem is forcing constraints, but I think
198// that is ok.
199// By doing it in parallel rather than sequentially,
200// we may miss transformations due to variables that were fixed
201// by forcing constraints, though.
202//
203// Note that both of these operations will cause problems
204// if the variables in question really need to exceed their bounds in order
205// to make the problem feasible.
206const PresolveAction *forcing_constraint_action::presolve(PresolveMatrix *prob,
207                                          const PresolveAction *next)
208{
209  double *clo   = prob->clo_;
210  double *cup   = prob->cup_;
211
212  const double *rowels  = prob->rowels_;
213  const int *hcol       = prob->hcol_;
214  const int *mrstrt     = prob->mrstrt_;
215  const int *hinrow     = prob->hinrow_;
216  const int nrows       = prob->nrows_;
217
218  const double *rlo     = prob->rlo_;
219  const double *rup     = prob->rup_;
220
221  const char *integerType = prob->integerType_;
222
223  const double tol      = ZTOLDP;
224  const double inftol   = prob->feasibilityTolerance_;
225  const int ncols       = prob->ncols_;
226
227  int *fixed_cols       = new int[ncols];
228  int nfixed_cols       = 0;
229
230  action *actions       = new action [nrows];
231  int nactions = 0;
232
233  int *useless_rows     = new int[nrows];
234  int nuseless_rows     = 0;
235
236  int numberLook = prob->numberRowsToDo_;
237  int iLook;
238  int * look = prob->rowsToDo_;
239
240  for (iLook=0;iLook<numberLook;iLook++) {
241    int irow = look[iLook];
242    if (hinrow[irow] > 0) {
243      int krs = mrstrt[irow];
244      int kre = krs + hinrow[irow];
245
246      double maxup, maxdown;
247      implied_row_bounds(rowels, clo, cup, hcol, krs, kre,
248                         &maxup, &maxdown);
249
250      if (maxup < PRESOLVE_INF && maxup + inftol < rlo[irow]) {
251        /* there is an upper bound and it can't be reached */
252        prob->status_|= 1;
253        prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
254                                             prob->originalModel_->messages())
255                                               <<irow
256                                               <<rlo[irow]
257                                               <<rup[irow]
258                                               <<CoinMessageEol;
259        break;
260      } else if (-PRESOLVE_INF < maxdown && rup[irow] < maxdown - inftol) {
261        /* there is a lower bound and it can't be reached */
262        prob->status_|= 1;
263        prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
264                                             prob->originalModel_->messages())
265                                               <<irow
266                                               <<rlo[irow]
267                                               <<rup[irow]
268                                               <<CoinMessageEol;
269        break;
270      }
271      // ADD TOLERANCE TO THESE TESTS
272      else if ((rlo[irow] <= -PRESOLVE_INF ||
273                (-PRESOLVE_INF < maxdown && rlo[irow] <= maxdown)) &&
274               (rup[irow] >= PRESOLVE_INF ||
275                (maxup < PRESOLVE_INF && rup[irow] >= maxup))) {
276
277        // I'm not sure that these transforms don't intefere with each other
278        // We can get it next time
279        if (some_col_was_fixed(hcol, krs, kre, clo, cup)) {
280          // make sure on next time
281          prob->addRow(irow);
282          continue;
283        }
284
285        // this constraint must always be satisfied - drop it
286        useless_rows[nuseless_rows++] = irow;
287
288      } else if ((maxup < PRESOLVE_INF && fabs(rlo[irow] - maxup) < tol) ||
289                 (-PRESOLVE_INF < maxdown && fabs(rup[irow] - maxdown) < tol)) {
290       
291        // the lower bound can just be reached, or
292        // the upper bound can just be reached;
293        // called a "forcing constraint" in the paper (p. 226)
294        const int lbound_tight = (maxup < PRESOLVE_INF &&
295                                  fabs(rlo[irow] - maxup) < tol);
296
297        // I'm not sure that these transforms don't intefere with each other
298        // We can get it next time
299        if (some_col_was_fixed(hcol, krs, kre, clo, cup)) {
300          // make sure on next time
301          prob->addRow(irow);
302          continue;
303        }
304
305        // out of space - this probably never happens
306        if (nfixed_cols + (kre-krs) >= ncols)
307          break;
308
309        double *bounds = new double[hinrow[irow]];
310        int *rowcols = new int[hinrow[irow]];
311        int lk = krs;   // load fix-to-down in front
312        int uk = kre;   // load fix-to-up in back
313
314        for (int k=krs; k<kre; k++) {
315          int jcol = hcol[k];
316          prob->addCol(jcol);
317          double coeff = rowels[k];
318
319          PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
320
321          // one of the two contributed to maxup - set the other to that
322          if (lbound_tight == (coeff > 0.0)) {
323            --uk;
324            bounds[uk-krs] = clo[jcol];
325            rowcols[uk-krs] = jcol;
326             
327            clo[jcol] = cup[jcol];
328          } else {
329            bounds[lk-krs] = cup[jcol];
330            rowcols[lk-krs] = jcol;
331            ++lk;
332
333            cup[jcol] = clo[jcol];
334          }
335
336          fixed_cols[nfixed_cols++] = jcol;
337        }
338        PRESOLVEASSERT(uk == lk);
339
340        action *f = &actions[nactions];
341        nactions++;
342
343        f->row = irow;
344        f->nlo = lk-krs;
345        f->nup = kre-uk;
346        f->rowcols = rowcols;
347        f->bounds = bounds;
348      }
349    }
350  }
351
352
353  if (nactions) {
354#if     PRESOLVE_SUMMARY
355    printf("NFORCED:  %d\n", nactions);
356#endif
357    next = new forcing_constraint_action(nactions, 
358                                         copyOfArray(actions,nactions), next);
359  }
360  delete [] actions;
361  if (nuseless_rows) {
362    next = useless_constraint_action::presolve(prob,
363                                               useless_rows, nuseless_rows,
364                                               next);
365  }
366  delete[]useless_rows;
367
368  if (nfixed_cols) {
369    next = remove_fixed_action::presolve(prob, fixed_cols, nfixed_cols, next);
370  }
371  delete[]fixed_cols;
372
373  return (next);
374}
375
376void forcing_constraint_action::postsolve(PostsolveMatrix *prob) const
377{
378  const action *const actions = actions_;
379  const int nactions = nactions_;
380
381  const double *colels  = prob->colels_;
382  const int *hrow               = prob->hrow_;
383  const int *mcstrt             = prob->mcstrt_;
384  const int *hincol             = prob->hincol_;
385  const int *link               = prob->link_;
386
387  int free_list = prob->free_list_;
388
389  double *clo   = prob->clo_;
390  double *cup   = prob->cup_;
391
392  const double *sol     = prob->sol_;
393  double *rcosts        = prob->rcosts_;
394
395  //double *acts        = prob->acts_;
396  double *rowduals = prob->rowduals_;
397
398  const double ztoldj   = prob->ztoldj_;
399  const double ztolzb   = prob->ztolzb_;
400
401  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
402
403    const int irow      = f->row;
404    const int nlo       = f->nlo;
405    const int nup       = f->nup;
406    const int ninrow    = nlo + nup;
407    const int *rowcols  = f->rowcols;
408    const double *bounds= f->bounds;
409
410    for (int k=0; k<nlo; k++) {
411      int jcol = rowcols[k];
412      cup[jcol] = bounds[k];
413      PRESOLVEASSERT(prob->getColumnStatus(jcol)!=PrePostsolveMatrix::basic);
414    }
415
416    for (int k=nlo; k<ninrow; k++) {
417      int jcol = rowcols[k];
418      clo[jcol] = bounds[k];
419      PRESOLVEASSERT(prob->getColumnStatus(jcol)!=PrePostsolveMatrix::basic);
420    }
421
422    PRESOLVEASSERT(prob->getRowStatus(irow)==PrePostsolveMatrix::basic);
423    PRESOLVEASSERT(rowduals[irow] == 0.0);
424
425    // this is a lazy implementation.
426    // we tightened the col bounds, then let them be eliminated
427    // by repeated uses of FIX_VARIABLE and a final DROP_ROW.
428    // Therefore, by this point the row has been marked basic,
429    // the rowdual of this row is 0.0,
430    // and the reduced costs for the cols may or may not be ok
431    // for the relaxed column bounds.
432    //
433    // find the one most out of whack and fix it.
434    int whacked = -1;
435    double whack = 0.0;
436    for (int k=0; k<ninrow; k++) {
437      int jcol = rowcols[k];
438      int kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link);
439
440      // choose rowdual to cancel out reduced cost
441      double whack0 = rcosts[jcol] / colels[kk];
442
443      if (((rcosts[jcol] > ztoldj  && !(fabs(sol[jcol] - clo[jcol]) <= ztolzb)) ||
444           (rcosts[jcol] < -ztoldj && !(fabs(sol[jcol] - cup[jcol]) <= ztolzb))) &&
445          fabs(whack0) > fabs(whack)) {
446        whacked = jcol;
447        whack = whack0;
448      }
449    }
450
451    if (whacked != -1) {
452      prob->setColumnStatus(whacked,PrePostsolveMatrix::basic);
453      prob->setRowStatus(irow,PrePostsolveMatrix::atLowerBound);
454      rowduals[irow] = whack;
455
456      for (int k=0; k<ninrow; k++) {
457        int jcol = rowcols[k];
458        int kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link);
459             
460        rcosts[jcol] -= (rowduals[irow] * colels[kk]);
461      }
462    }
463  }
464}
465
466
467
468
469// Determine the maximum and minimum values the constraint sums
470// may take, given the bounds on the variables.
471// If there are infinite terms, record where the first one is,
472// and whether there is more than one.
473// It is possible to compute implied bounds for the (one) variable
474// with no bound.
475static void implied_bounds1(PresolveMatrix * prob, const double *rowels,
476                                const int *mrstrt,
477                                const int *hrow,
478                                const int *hinrow,
479                                const double *clo, const double *cup,
480                                const int *hcol,
481                                int ncols,
482                                const double *rlo, const double *rup,
483                                const char *integerType,
484                                int nrows,
485                                double *ilbound, double *iubound)
486{
487  const double tol = ZTOLDP;
488
489  for (int irow=0; irow<nrows; irow++) {
490    int krs = mrstrt[irow];
491    int kre = krs + hinrow[irow];
492
493    double irlo = rlo[irow];
494    double irup = rup[irow];
495
496    // These are used to set column bounds below.
497    // If there are no (positive) infinite terms,
498    // the loop will range from krs to kre;
499    // if there is just one, it will range over that one variable;
500    // otherwise, it will be empty.
501    int ub_inf_index = -1;
502    int lb_inf_index = -1;
503
504    double maxup = 0.0;
505    double maxdown = 0.0;
506
507    for (int k=krs; k<kre; k++) {
508      int jcol = hcol[k];
509      double coeff = rowels[k];
510      double lb = clo[jcol];
511      double ub = cup[jcol];
512
513      // HAVE TO DEAL WITH BOUNDS OF INTEGER VARIABLES
514      if (coeff > 0.0) {
515        if (PRESOLVE_INF <= ub) {
516          if (ub_inf_index == -1) {
517            ub_inf_index = k;
518          } else {
519            ub_inf_index = -2;
520            if (lb_inf_index == -2)
521              break;    // pointless
522          }
523        } else
524          maxup += ub * coeff;
525
526        if (lb <= -PRESOLVE_INF) {
527          if (lb_inf_index == -1) {
528            lb_inf_index = k;
529          } else {
530            lb_inf_index = -2;
531            if (ub_inf_index == -2)
532              break;    // pointless
533          }
534        } else
535          maxdown += lb * coeff;
536      }
537      else {
538        if (PRESOLVE_INF <= ub) {
539          if (lb_inf_index == -1) {
540            lb_inf_index = k;
541          } else {
542            lb_inf_index = -2;
543            if (ub_inf_index == -2)
544              break;    // pointless
545          }
546        } else
547          maxdown += ub * coeff;
548
549        if (lb <= -PRESOLVE_INF) {
550          if (ub_inf_index == -1) {
551            ub_inf_index = k;
552          } else {
553            ub_inf_index = -2;
554            if (lb_inf_index == -2)
555              break;    // pointless
556          }
557        } else
558          maxup += lb * coeff;
559      }
560    }
561
562    // ub_inf says whether the sum of the "other" ub terms is infinite
563    // in the loop below.
564    // In the case where we only saw one infinite term, the loop
565    // will only cover that case, in which case the other terms
566    // are *not* infinite.
567    // With two or more terms, it is infinite.
568    // If we only saw one infinite term, then
569    if (ub_inf_index == -2)
570      maxup = PRESOLVE_INF;
571
572    if (lb_inf_index == -2)
573      maxdown = -PRESOLVE_INF;
574
575    const bool maxup_finite = PRESOLVEFINITE(maxup);
576    const bool maxdown_finite = PRESOLVEFINITE(maxdown);
577
578    if (ub_inf_index == -1 && maxup_finite && maxup + tol < rlo[irow]) {
579      /* infeasible */
580        prob->status_|= 1;
581        prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
582                                             prob->originalModel_->messages())
583                                               <<irow
584                                               <<rlo[irow]
585                                               <<rup[irow]
586                                               <<CoinMessageEol;
587        break;
588    } else if (lb_inf_index == -1 && maxdown_finite && rup[irow] < maxdown - tol) {
589      /* infeasible */
590        prob->status_|= 1;
591        prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
592                                             prob->originalModel_->messages())
593                                               <<irow
594                                               <<rlo[irow]
595                                               <<rup[irow]
596                                               <<CoinMessageEol;
597        break;
598    }
599
600    for (int k = krs; k<kre; ++k) {
601      int jcol = hcol[k];
602      double coeff = rowels[k];
603
604      // SHOULD GET RID OF THIS
605      if (fabs(coeff) > ZTOLDP &&
606          !integerType[jcol]) {
607        double maxup1 = (ub_inf_index == -1 || ub_inf_index == k
608                         ? maxup
609                         : PRESOLVE_INF);
610        bool maxup_finite1 = (ub_inf_index == -1 || ub_inf_index == k
611                              ? maxup_finite
612                              : false);
613        double maxdown1 = (lb_inf_index == -1 || lb_inf_index == k
614                         ? maxdown
615                         : PRESOLVE_INF);
616        bool maxdown_finite1 = (ub_inf_index == -1 || ub_inf_index == k
617                              ? maxdown_finite
618                              : false);
619
620        double ilb = (irlo - maxup1) / coeff;
621        bool finite_ilb = (-PRESOLVE_INF < irlo && maxup_finite1);
622
623        double iub = (irup - maxdown1) / coeff;
624        bool finite_iub = ( irup < PRESOLVE_INF && maxdown_finite1);
625
626        double ilb1 = (coeff > 0.0
627                       ? (finite_ilb ? ilb : -PRESOLVE_INF)
628                       : (finite_iub ? iub : -PRESOLVE_INF));
629
630        if (ilbound[jcol] < ilb1) {
631          ilbound[jcol] = ilb1;
632          if (jcol == 278001)
633            printf("JCOL LB %g\n", ilb1);
634        }
635      }
636    }
637
638    for (int k = krs; k<kre; ++k) {
639      int jcol = hcol[k];
640      double coeff = rowels[k];
641
642      // SHOULD GET RID OF THIS
643      if (fabs(coeff) > ZTOLDP &&
644          !integerType[jcol]) {
645        double maxup1 = (ub_inf_index == -1 || ub_inf_index == k
646                         ? maxup
647                         : PRESOLVE_INF);
648        bool maxup_finite1 = (ub_inf_index == -1 || ub_inf_index == k
649                              ? maxup_finite
650                              : false);
651        double maxdown1 = (lb_inf_index == -1 || lb_inf_index == k
652                         ? maxdown
653                         : PRESOLVE_INF);
654        bool maxdown_finite1 = (ub_inf_index == -1 || ub_inf_index == k
655                              ? maxdown_finite
656                              : false);
657
658
659        double ilb = (irlo - maxup1) / coeff;
660        bool finite_ilb = (-PRESOLVE_INF < irlo && maxup_finite1);
661
662        double iub = (irup - maxdown1) / coeff;
663        bool finite_iub = ( irup < PRESOLVE_INF && maxdown_finite1);
664
665        double iub1 = (coeff > 0.0
666                       ? (finite_iub ? iub :  PRESOLVE_INF)
667                       : (finite_ilb ? ilb :  PRESOLVE_INF));
668
669        if (iub1 < iubound[jcol]) {
670          iubound[jcol] = iub1;
671          if (jcol == 278001)
672            printf("JCOL UB %g\n", iub1);
673        }
674      }
675    }
676  }
677}
678
679#if 0
680postsolve for implied_bound
681        {
682          double lo0    = pa->clo;
683          double up0    = pa->cup;
684          int irow      = pa->irow;
685          int jcol      = pa->icol;
686          int *rowcols  = pa->rowcols;
687          int ninrow    = pa->ninrow;
688
689          clo[jcol] = lo0;
690          cup[jcol] = up0;
691
692          if ((colstat[jcol] & PRESOLVE_XBASIC) == 0 &&
693              fabs(lo0 - sol[jcol]) > ztolzb &&
694              fabs(up0 - sol[jcol]) > ztolzb) {
695
696            // this non-basic variable is now away from its bound
697            // it is ok just to force it to be basic
698            // informally:  if this variable is at its implied bound,
699            // then the other variables must be at their bounds,
700            // which means the bounds will stop them even if the aren't basic.
701            if (rowstat[irow] & PRESOLVE_XBASIC)
702              rowstat[irow] = 0;
703            else {
704              int k;
705              for (k=0; k<ninrow; k++) {
706                int col = rowcols[k];
707                if (cdone[col] &&
708                    (colstat[col] & PRESOLVE_XBASIC) &&
709                    ((fabs(clo[col] - sol[col]) <= ztolzb && rcosts[col] >= -ztoldj) ||
710                     (fabs(cup[col] - sol[col]) <= ztolzb && rcosts[col] <= ztoldj)))
711                  break;
712              }
713              if (k<ninrow) {
714                int col = rowcols[k];
715                // steal this basic variable
716#if     DEBUG_PRESOLVE
717                printf("PIVOTING ON COL:  %d %d -> %d\n", irow, col, jcol);
718#endif
719                colstat[col] = 0;
720
721                // since all vars were at their bounds, the slack must be 0
722                PRESOLVEASSERT(fabs(acts[irow]) < ZTOLDP);
723                rowstat[irow] = PRESOLVE_XBASIC;
724              }
725              else {
726                // should never happen?
727                abort();
728              }
729
730              // get rid of any remaining basic structurals, since their rcosts
731              // are going to become non-zero in a second.
732              abort();
733              ///////////////////zero_pivot();
734            }
735
736            double rdual_adjust;
737            {
738              int kk = presolve_find_row(irow, mcstrt[jcol], mcstrt[jcol] + hincol[jcol], hrow);
739              // adjust rowdual to cancel out reduced cost
740              // should probably search for col with largest factor
741              rdual_adjust = (rcosts[jcol] / colels[kk]);
742              rowduals[irow] += rdual_adjust;
743              colstat[jcol] = PRESOLVE_XBASIC;
744            }
745
746            for (int k=0; k<ninrow; k++) {
747              int jcol = rowcols[k];
748              int kk = presolve_find_row(irow, mcstrt[jcol], mcstrt[jcol] + hincol[jcol], hrow);
749             
750              rcosts[jcol] -= (rdual_adjust * colels[kk]);
751            }
752
753            {
754              int k;
755              int badbasic = -1;
756
757              // we may have just screwed up the rcost of another basic variable
758              for (k=0; k<ninrow; k++) {
759                int col = rowcols[k];
760                if (col != jcol &&
761                    cdone[col] &&
762                    (colstat[col] & PRESOLVE_XBASIC) &&
763                    !(fabs(rcosts[col]) < ztoldj))
764                  if (badbasic == -1)
765                    badbasic = k;
766                  else
767                    abort();    // two!!  what to do???
768              }
769
770              if (badbasic != -1) {
771                int col = rowcols[badbasic];
772
773                if (fabs(acts[irow]) < ZTOLDP) {
774#if     DEBUG_PRESOLVE
775                  printf("PIVOTING COL TO SLACK!:  %d %d\n", irow, col);
776#endif
777                  colstat[col] = 0;
778                  rowstat[irow] = PRESOLVE_XBASIC;
779                }
780                else
781                  abort();
782              }
783            }
784          }
785        }
786#endif
787
788
789
Note: See TracBrowser for help on using the repository browser.