source: trunk/PresolveForcing.cpp @ 95

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

Was trying something for dual - didn't really work

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