source: trunk/Clp/src/ClpPresolve.cpp @ 1896

Last change on this file since 1896 was 1896, checked in by lou, 7 years ago

Minimal changes to presolve driver for testRedundant, fill_level.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 83.5 KB
Line 
1/* $Id: ClpPresolve.cpp 1896 2012-11-29 19:34:35Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6//#define       PRESOLVE_CONSISTENCY    1
7//#define       PRESOLVE_DEBUG  1
8
9#include <stdio.h>
10
11#include <cassert>
12#include <iostream>
13
14#include "CoinHelperFunctions.hpp"
15#if CLP_HAS_ABC
16#include "CoinAbcCommon.hpp"
17#endif
18
19#include "CoinPackedMatrix.hpp"
20#include "ClpPackedMatrix.hpp"
21#include "ClpSimplex.hpp"
22#include "ClpSimplexOther.hpp"
23#ifndef SLIM_CLP
24#include "ClpQuadraticObjective.hpp"
25#endif
26
27#include "ClpPresolve.hpp"
28#include "CoinPresolveMatrix.hpp"
29
30#include "CoinPresolveEmpty.hpp"
31#include "CoinPresolveFixed.hpp"
32#include "CoinPresolvePsdebug.hpp"
33#include "CoinPresolveSingleton.hpp"
34#include "CoinPresolveDoubleton.hpp"
35#include "CoinPresolveTripleton.hpp"
36#include "CoinPresolveZeros.hpp"
37#include "CoinPresolveSubst.hpp"
38#include "CoinPresolveForcing.hpp"
39#include "CoinPresolveDual.hpp"
40#include "CoinPresolveTighten.hpp"
41#include "CoinPresolveUseless.hpp"
42#include "CoinPresolveDupcol.hpp"
43#include "CoinPresolveImpliedFree.hpp"
44#include "CoinPresolveIsolated.hpp"
45#include "CoinMessage.hpp"
46
47
48
49ClpPresolve::ClpPresolve() :
50     originalModel_(NULL),
51     presolvedModel_(NULL),
52     nonLinearValue_(0.0),
53     originalColumn_(NULL),
54     originalRow_(NULL),
55     rowObjective_(NULL),
56     paction_(0),
57     ncols_(0),
58     nrows_(0),
59     nelems_(0),
60#ifdef ABC_INHERIT
61     numberPasses_(20),
62#else
63     numberPasses_(5),
64#endif
65     substitution_(3),
66#ifndef CLP_NO_STD
67     saveFile_(""),
68#endif
69     presolveActions_(0)
70{
71}
72
73ClpPresolve::~ClpPresolve()
74{
75     destroyPresolve();
76}
77// Gets rid of presolve actions (e.g.when infeasible)
78void
79ClpPresolve::destroyPresolve()
80{
81     const CoinPresolveAction *paction = paction_;
82     while (paction) {
83          const CoinPresolveAction *next = paction->next;
84          delete paction;
85          paction = next;
86     }
87     delete [] originalColumn_;
88     delete [] originalRow_;
89     paction_ = NULL;
90     originalColumn_ = NULL;
91     originalRow_ = NULL;
92     delete [] rowObjective_;
93     rowObjective_ = NULL;
94}
95
96/* This version of presolve returns a pointer to a new presolved
97   model.  NULL if infeasible
98*/
99ClpSimplex *
100ClpPresolve::presolvedModel(ClpSimplex & si,
101                            double feasibilityTolerance,
102                            bool keepIntegers,
103                            int numberPasses,
104                            bool dropNames,
105                            bool doRowObjective,
106                            const char * prohibitedRows,
107                            const char * prohibitedColumns)
108{
109     // Check matrix
110     int checkType = ((si.specialOptions() & 128) != 0) ? 14 : 15;
111     if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
112                                             1.0e20,checkType))
113          return NULL;
114     else
115          return gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
116                                      doRowObjective,
117                                      prohibitedRows,
118                                      prohibitedColumns);
119}
120#ifndef CLP_NO_STD
121/* This version of presolve updates
122   model and saves original data to file.  Returns non-zero if infeasible
123*/
124int
125ClpPresolve::presolvedModelToFile(ClpSimplex &si, std::string fileName,
126                                  double feasibilityTolerance,
127                                  bool keepIntegers,
128                                  int numberPasses,
129                                  bool dropNames,
130                                  bool doRowObjective)
131{
132     // Check matrix
133     if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
134                                             1.0e20))
135          return 2;
136     saveFile_ = fileName;
137     si.saveModel(saveFile_.c_str());
138     ClpSimplex * model = gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
139                          doRowObjective);
140     if (model == &si) {
141          return 0;
142     } else {
143          si.restoreModel(saveFile_.c_str());
144          remove(saveFile_.c_str());
145          return 1;
146     }
147}
148#endif
149// Return pointer to presolved model
150ClpSimplex *
151ClpPresolve::model() const
152{
153     return presolvedModel_;
154}
155// Return pointer to original model
156ClpSimplex *
157ClpPresolve::originalModel() const
158{
159     return originalModel_;
160}
161// Return presolve status (0,1,2)
162int 
163ClpPresolve::presolveStatus() const
164{
165  if (nelems_>=0) {
166    // feasible (or not done yet)
167    return 0;
168  } else {
169    int presolveStatus = - nelems_;
170    // If both infeasible and unbounded - say infeasible
171    if (presolveStatus>2)
172      presolveStatus = 1;
173    return presolveStatus;
174  }
175}
176void
177ClpPresolve::postsolve(bool updateStatus)
178{
179     // Return at once if no presolved model
180     if (!presolvedModel_)
181          return;
182     // Messages
183     CoinMessages messages = originalModel_->coinMessages();
184     if (!presolvedModel_->isProvenOptimal()) {
185          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
186                    messages)
187                    << CoinMessageEol;
188     }
189
190     // this is the size of the original problem
191     const int ncols0  = ncols_;
192     const int nrows0  = nrows_;
193     const CoinBigIndex nelems0 = nelems_;
194
195     // this is the reduced problem
196     int ncols = presolvedModel_->getNumCols();
197     int nrows = presolvedModel_->getNumRows();
198
199     double * acts = NULL;
200     double * sol = NULL;
201     unsigned char * rowstat = NULL;
202     unsigned char * colstat = NULL;
203#ifndef CLP_NO_STD
204     if (saveFile_ == "") {
205#endif
206          // reality check
207          assert(ncols0 == originalModel_->getNumCols());
208          assert(nrows0 == originalModel_->getNumRows());
209          acts = originalModel_->primalRowSolution();
210          sol  = originalModel_->primalColumnSolution();
211          if (updateStatus) {
212               // postsolve does not know about fixed
213               int i;
214               for (i = 0; i < nrows + ncols; i++) {
215                    if (presolvedModel_->getColumnStatus(i) == ClpSimplex::isFixed)
216                         presolvedModel_->setColumnStatus(i, ClpSimplex::atLowerBound);
217               }
218               unsigned char *status = originalModel_->statusArray();
219               if (!status) {
220                    originalModel_->createStatus();
221                    status = originalModel_->statusArray();
222               }
223               rowstat = status + ncols0;
224               colstat = status;
225               CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
226               CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
227          }
228#ifndef CLP_NO_STD
229     } else {
230          // from file
231          acts = new double[nrows0];
232          sol  = new double[ncols0];
233          CoinZeroN(acts, nrows0);
234          CoinZeroN(sol, ncols0);
235          if (updateStatus) {
236               unsigned char *status = new unsigned char [nrows0+ncols0];
237               rowstat = status + ncols0;
238               colstat = status;
239               CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
240               CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
241          }
242     }
243#endif
244
245     // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
246     // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
247     // cause duplicate free. In case where saveFile == "", as best I can see
248     // arrays are owned by originalModel_. fix is to
249     // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
250     CoinPostsolveMatrix prob(presolvedModel_,
251                              ncols0,
252                              nrows0,
253                              nelems0,
254                              presolvedModel_->getObjSense(),
255                              // end prepost
256
257                              sol, acts,
258                              colstat, rowstat);
259
260     postsolve(prob);
261
262#ifndef CLP_NO_STD
263     if (saveFile_ != "") {
264          // From file
265          assert (originalModel_ == presolvedModel_);
266          originalModel_->restoreModel(saveFile_.c_str());
267          remove(saveFile_.c_str());
268          CoinMemcpyN(acts, nrows0, originalModel_->primalRowSolution());
269          // delete [] acts;
270          CoinMemcpyN(sol, ncols0, originalModel_->primalColumnSolution());
271          // delete [] sol;
272          if (updateStatus) {
273               CoinMemcpyN(colstat, nrows0 + ncols0, originalModel_->statusArray());
274               // delete [] colstat;
275          }
276     } else {
277#endif
278          prob.sol_ = 0 ;
279          prob.acts_ = 0 ;
280          prob.colstat_ = 0 ;
281#ifndef CLP_NO_STD
282     }
283#endif
284     // put back duals
285     CoinMemcpyN(prob.rowduals_,        nrows_, originalModel_->dualRowSolution());
286     double maxmin = originalModel_->getObjSense();
287     if (maxmin < 0.0) {
288          // swap signs
289          int i;
290          double * pi = originalModel_->dualRowSolution();
291          for (i = 0; i < nrows_; i++)
292               pi[i] = -pi[i];
293     }
294     // Now check solution
295     double offset;
296     CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
297                 originalModel_->primalColumnSolution(), offset, true),
298                 ncols_, originalModel_->dualColumnSolution());
299     originalModel_->clpMatrix()->transposeTimes(-1.0,
300                                    originalModel_->dualRowSolution(),
301                                    originalModel_->dualColumnSolution());
302     memset(originalModel_->primalRowSolution(), 0, nrows_ * sizeof(double));
303     originalModel_->clpMatrix()->times(1.0, 
304                                        originalModel_->primalColumnSolution(),
305                                        originalModel_->primalRowSolution());
306     originalModel_->checkSolutionInternal();
307     if (originalModel_->sumDualInfeasibilities() > 1.0e-1) {
308          // See if we can fix easily
309          static_cast<ClpSimplexOther *> (originalModel_)->cleanupAfterPostsolve();
310     }
311     // Messages
312     presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
313               messages)
314               << originalModel_->objectiveValue()
315               << originalModel_->sumDualInfeasibilities()
316               << originalModel_->numberDualInfeasibilities()
317               << originalModel_->sumPrimalInfeasibilities()
318               << originalModel_->numberPrimalInfeasibilities()
319               << CoinMessageEol;
320
321     //originalModel_->objectiveValue_=objectiveValue_;
322     originalModel_->setNumberIterations(presolvedModel_->numberIterations());
323     if (!presolvedModel_->status()) {
324          if (!originalModel_->numberDualInfeasibilities() &&
325                    !originalModel_->numberPrimalInfeasibilities()) {
326               originalModel_->setProblemStatus( 0);
327          } else {
328               originalModel_->setProblemStatus( -1);
329               // Say not optimal after presolve
330               originalModel_->setSecondaryStatus(7);
331               presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
332                         messages)
333                         << CoinMessageEol;
334          }
335     } else {
336          originalModel_->setProblemStatus( presolvedModel_->status());
337          // but not if close to feasible
338          if( originalModel_->sumPrimalInfeasibilities()<1.0e-1) {
339               originalModel_->setProblemStatus( -1);
340               // Say not optimal after presolve
341               originalModel_->setSecondaryStatus(7);
342          }
343     }
344#ifndef CLP_NO_STD
345     if (saveFile_ != "")
346          presolvedModel_ = NULL;
347#endif
348}
349
350// return pointer to original columns
351const int *
352ClpPresolve::originalColumns() const
353{
354     return originalColumn_;
355}
356// return pointer to original rows
357const int *
358ClpPresolve::originalRows() const
359{
360     return originalRow_;
361}
362// Set pointer to original model
363void
364ClpPresolve::setOriginalModel(ClpSimplex * model)
365{
366     originalModel_ = model;
367}
368#if 0
369// A lazy way to restrict which transformations are applied
370// during debugging.
371static int ATOI(const char *name)
372{
373     return true;
374#if     PRESOLVE_DEBUG || PRESOLVE_SUMMARY
375     if (getenv(name)) {
376          int val = atoi(getenv(name));
377          printf("%s = %d\n", name, val);
378          return (val);
379     } else {
380          if (strcmp(name, "off"))
381               return (true);
382          else
383               return (false);
384     }
385#else
386     return (true);
387#endif
388}
389#endif
390//#define PRESOLVE_DEBUG 1
391#if PRESOLVE_DEBUG
392void check_sol(CoinPresolveMatrix *prob, double tol)
393{
394     double *colels     = prob->colels_;
395     int *hrow          = prob->hrow_;
396     int *mcstrt                = prob->mcstrt_;
397     int *hincol                = prob->hincol_;
398     int *hinrow                = prob->hinrow_;
399     int ncols          = prob->ncols_;
400
401
402     double * csol = prob->sol_;
403     double * acts = prob->acts_;
404     double * clo = prob->clo_;
405     double * cup = prob->cup_;
406     int nrows = prob->nrows_;
407     double * rlo = prob->rlo_;
408     double * rup = prob->rup_;
409
410     int colx;
411
412     double * rsol = new double[nrows];
413     memset(rsol, 0, nrows * sizeof(double));
414
415     for (colx = 0; colx < ncols; ++colx) {
416          if (1) {
417               CoinBigIndex k = mcstrt[colx];
418               int nx = hincol[colx];
419               double solutionValue = csol[colx];
420               for (int i = 0; i < nx; ++i) {
421                    int row = hrow[k];
422                    double coeff = colels[k];
423                    k++;
424                    rsol[row] += solutionValue * coeff;
425               }
426               if (csol[colx] < clo[colx] - tol) {
427                    printf("low CSOL:  %d  - %g %g %g\n",
428                           colx, clo[colx], csol[colx], cup[colx]);
429               } else if (csol[colx] > cup[colx] + tol) {
430                    printf("high CSOL:  %d  - %g %g %g\n",
431                           colx, clo[colx], csol[colx], cup[colx]);
432               }
433          }
434     }
435     int rowx;
436     for (rowx = 0; rowx < nrows; ++rowx) {
437          if (hinrow[rowx]) {
438               if (fabs(rsol[rowx] - acts[rowx]) > tol)
439                    printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
440                           rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
441               if (rsol[rowx] < rlo[rowx] - tol) {
442                    printf("low RSOL:  %d - %g %g %g\n",
443                           rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
444               } else if (rsol[rowx] > rup[rowx] + tol ) {
445                    printf("high RSOL:  %d - %g %g %g\n",
446                           rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
447               }
448          }
449     }
450     delete [] rsol;
451}
452#endif
453static int tightenDoubletons2(CoinPresolveMatrix * prob)
454{
455  // column-major representation
456  const int ncols = prob->ncols_ ;
457  const CoinBigIndex *const mcstrt = prob->mcstrt_ ;
458  const int *const hincol = prob->hincol_ ;
459  const int *const hrow = prob->hrow_ ;
460  double * colels = prob->colels_ ;
461  double * cost = prob->cost_ ;
462
463  // column type, bounds, solution, and status
464  const unsigned char *const integerType = prob->integerType_ ;
465  double * clo = prob->clo_ ;
466  double * cup = prob->cup_ ;
467  // row-major representation
468  //const int nrows = prob->nrows_ ;
469  const CoinBigIndex *const mrstrt = prob->mrstrt_ ;
470  const int *const hinrow = prob->hinrow_ ;
471  const int *const hcol = prob->hcol_ ;
472  double * rowels = prob->rowels_ ;
473
474  // row bounds
475  double *const rlo = prob->rlo_ ;
476  double *const rup = prob->rup_ ;
477
478  // tolerances
479  //const double ekkinf2 = PRESOLVE_SMALL_INF ;
480  //const double ekkinf = ekkinf2*1.0e8 ;
481  //const double ztolcbarj = prob->ztoldj_ ;
482  //const CoinRelFltEq relEq(prob->ztolzb_) ;
483  int numberChanged=0;
484  double bound[2];
485  double alpha[2]={0.0,0.0};
486  double offset=0.0;
487
488  for (int icol=0;icol<ncols;icol++) {
489    if (hincol[icol]==2) {
490      CoinBigIndex start=mcstrt[icol];
491      int row0 = hrow[start];
492      if (hinrow[row0]!=2)
493        continue;
494      int row1 = hrow[start+1];
495      if (hinrow[row1]!=2)
496        continue;
497      double element0 = colels[start];
498      double rowUpper0=rup[row0];
499      bool swapSigns0=false;
500      if (rlo[row0]>-1.0e30) {
501        if (rup[row0]>1.0e30) {
502          swapSigns0=true;
503          rowUpper0=-rlo[row0];
504          element0=-element0;
505        } else {
506          // range or equality
507          continue;
508        }
509      } else if (rup[row0]>1.0e30) {
510        // free
511        continue;
512      }
513#if 0
514      // skip here for speed
515      // skip if no cost (should be able to get rid of)
516      if (!cost[icol]) {
517        printf("should be able to get rid of %d with no cost\n",icol);
518        continue;
519      }
520      // skip if negative cost for now
521      if (cost[icol]<0.0) {
522        printf("code for negative cost\n");
523        continue;
524      }
525#endif
526      double element1 = colels[start+1];
527      double rowUpper1=rup[row1];
528      bool swapSigns1=false;
529      if (rlo[row1]>-1.0e30) {
530        if (rup[row1]>1.0e30) {
531          swapSigns1=true;
532          rowUpper1=-rlo[row1];
533          element1=-element1;
534        } else {
535          // range or equality
536          continue;
537        }
538      } else if (rup[row1]>1.0e30) {
539        // free
540        continue;
541      }
542      double lowerX=clo[icol];
543      double upperX=cup[icol];
544      int otherCol=-1;
545      CoinBigIndex startRow=mrstrt[row0];
546      for (CoinBigIndex j=startRow;j<startRow+2;j++) {
547        int jcol=hcol[j];
548        if (jcol!=icol) {
549          alpha[0]=swapSigns0 ? -rowels[j] :rowels[j];
550          otherCol=jcol;
551        }
552      }
553      startRow=mrstrt[row1];
554      bool possible=true;
555      for (CoinBigIndex j=startRow;j<startRow+2;j++) {
556        int jcol=hcol[j];
557        if (jcol!=icol) {
558          if (jcol==otherCol) {
559            alpha[1]=swapSigns1 ? -rowels[j] :rowels[j];
560          } else {
561            possible=false;
562          }
563        }
564      }
565      if (possible) {
566        // skip if no cost (should be able to get rid of)
567        if (!cost[icol]) {
568          PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n",icol));
569          continue;
570        }
571        // skip if negative cost for now
572        if (cost[icol]<0.0) {
573          PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
574          continue;
575        }
576        bound[0]=clo[otherCol];
577        bound[1]=cup[otherCol];
578        double lowestLowest=COIN_DBL_MAX;
579        double highestLowest=-COIN_DBL_MAX;
580        double lowestHighest=COIN_DBL_MAX;
581        double highestHighest=-COIN_DBL_MAX;
582        int binding0=0;
583        int binding1=0;
584        for (int k=0;k<2;k++) {
585          bool infLow0=false;
586          bool infLow1=false;
587          double sum0=0.0;
588          double sum1=0.0;
589          double value=bound[k];
590          if (fabs(value)<1.0e30) {
591            sum0+=alpha[0]*value;
592            sum1+=alpha[1]*value;
593          } else {
594            if (alpha[0]>0.0) {
595              if (value<0.0)
596                infLow0 =true;
597            } else if (alpha[0]<0.0) {
598              if (value>0.0)
599                infLow0 =true;
600            }
601            if (alpha[1]>0.0) {
602              if (value<0.0)
603                infLow1 =true;
604            } else if (alpha[1]<0.0) {
605              if (value>0.0)
606                infLow1 =true;
607            }
608          }
609          /* Got sums
610           */
611          double thisLowest0=-COIN_DBL_MAX;
612          double thisHighest0=COIN_DBL_MAX;
613          if (element0>0.0) {
614            // upper bound unless inf&2 !=0
615            if (!infLow0)
616              thisHighest0 = (rowUpper0-sum0)/element0;
617          } else {
618            // lower bound unless inf&2 !=0
619            if (!infLow0)
620              thisLowest0 = (rowUpper0-sum0)/element0;
621          }
622          double thisLowest1=-COIN_DBL_MAX;
623          double thisHighest1=COIN_DBL_MAX;
624          if (element1>0.0) {
625            // upper bound unless inf&2 !=0
626            if (!infLow1)
627              thisHighest1 = (rowUpper1-sum1)/element1;
628          } else {
629            // lower bound unless inf&2 !=0
630            if (!infLow1)
631              thisLowest1 = (rowUpper1-sum1)/element1;
632          }
633          if (thisLowest0>thisLowest1+1.0e-12) {
634            if (thisLowest0>lowerX+1.0e-12)
635              binding0|= 1<<k;
636          } else if (thisLowest1>thisLowest0+1.0e-12) {
637            if (thisLowest1>lowerX+1.0e-12)
638              binding1|= 1<<k;
639            thisLowest0=thisLowest1;
640          }
641          if (thisHighest0<thisHighest1-1.0e-12) {
642            if (thisHighest0<upperX-1.0e-12)
643              binding0|= 1<<k;
644          } else if (thisHighest1<thisHighest0-1.0e-12) {
645            if (thisHighest1<upperX-1.0e-12)
646              binding1|= 1<<k;
647            thisHighest0=thisHighest1;
648          }
649          lowestLowest=CoinMin(lowestLowest,thisLowest0);
650          highestHighest=CoinMax(highestHighest,thisHighest0);
651          lowestHighest=CoinMin(lowestHighest,thisHighest0);
652          highestLowest=CoinMax(highestLowest,thisLowest0);
653        }
654        // see if any good
655        //#define PRINT_VALUES
656        if (!binding0||!binding1) {
657          PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n",icol));
658        } else {
659#ifdef PRINT_VALUES
660          printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
661                 icol,lowerX,upperX,lowestLowest,lowestHighest,
662                 highestLowest,highestHighest);
663#endif
664          // if integer adjust
665          if (integerType[icol]) {
666            lowestLowest=ceil(lowestLowest-1.0e-5);
667            highestLowest=ceil(highestLowest-1.0e-5);
668            lowestHighest=floor(lowestHighest+1.0e-5);
669            highestHighest=floor(highestHighest+1.0e-5);
670          }
671          // if costed may be able to adjust
672          if (cost[icol]>=0.0) {
673            if (highestLowest<upperX&&highestLowest>=lowerX&&highestHighest<1.0e30) {
674              highestHighest=CoinMin(highestHighest,highestLowest);
675            }
676          }
677          if (cost[icol]<=0.0) {
678            if (lowestHighest>lowerX&&lowestHighest<=upperX&&lowestHighest>-1.0e30) {
679              lowestLowest=CoinMax(lowestLowest,lowestHighest);
680            }
681          }
682#if 1
683          if (lowestLowest>lowerX+1.0e-8) {
684#ifdef PRINT_VALUES
685            printf("Can increase lower bound on %d from %g to %g\n",
686                   icol,lowerX,lowestLowest);
687#endif
688            lowerX=lowestLowest;
689          }
690          if (highestHighest<upperX-1.0e-8) {
691#ifdef PRINT_VALUES
692            printf("Can decrease upper bound on %d from %g to %g\n",
693                   icol,upperX,highestHighest);
694#endif
695            upperX=highestHighest;
696           
697          }
698#endif
699          // see if we can move costs
700          double xValue;
701          double yValue0;
702          double yValue1;
703          double newLower=COIN_DBL_MAX;
704          double newUpper=-COIN_DBL_MAX;
705          double ranges0[2];
706          double ranges1[2];
707          double costEqual;
708          double slope[2];
709          assert (binding0+binding1==3);
710          // get where equal
711          xValue=(rowUpper0*element1-rowUpper1*element0)/(alpha[0]*element1-alpha[1]*element0);
712          yValue0=(rowUpper0-xValue*alpha[0])/element0;
713          yValue1=(rowUpper1-xValue*alpha[1])/element1;
714          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
715          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
716          double xValueEqual=xValue;
717          double yValueEqual=yValue0;
718          costEqual = xValue*cost[otherCol]+yValueEqual*cost[icol];
719          if (binding0==1) {
720            ranges0[0]=bound[0];
721            ranges0[1]=yValue0;
722            ranges1[0]=yValue0;
723            ranges1[1]=bound[1];
724            // take x 1.0 down
725            double x=xValue-1.0;
726            double y=(rowUpper0-x*alpha[0])/element0;
727            double costTotal = x*cost[otherCol]+y*cost[icol];
728            slope[0] = costEqual-costTotal;
729            // take x 1.0 up
730            x=xValue+1.0;
731            y=(rowUpper1-x*alpha[1])/element0;
732            costTotal = x*cost[otherCol]+y*cost[icol];
733            slope[1] = costTotal-costEqual;
734          } else {
735            ranges1[0]=bound[0];
736            ranges1[1]=yValue0;
737            ranges0[0]=yValue0;
738            ranges0[1]=bound[1];
739            // take x 1.0 down
740            double x=xValue-1.0;
741            double y=(rowUpper1-x*alpha[1])/element0;
742            double costTotal = x*cost[otherCol]+y*cost[icol];
743            slope[1] = costEqual-costTotal;
744            // take x 1.0 up
745            x=xValue+1.0;
746            y=(rowUpper0-x*alpha[0])/element0;
747            costTotal = x*cost[otherCol]+y*cost[icol];
748            slope[0] = costTotal-costEqual;
749          }
750#ifdef PRINT_VALUES
751          printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
752                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
753          printf("Cost at equality %g for constraint 0 ranges %g -> %g slope %g for constraint 1 ranges %g -> %g slope %g\n",
754                 costEqual,ranges0[0],ranges0[1],slope[0],ranges1[0],ranges1[1],slope[1]);
755#endif
756          xValue=bound[0];
757          yValue0=(rowUpper0-xValue*alpha[0])/element0;
758          yValue1=(rowUpper1-xValue*alpha[1])/element1;
759#ifdef PRINT_VALUES
760          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
761                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
762#endif
763          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
764          // cost>0 so will be at lower
765          //double yValueAtBound0=newLower;
766          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
767          xValue=bound[1];
768          yValue0=(rowUpper0-xValue*alpha[0])/element0;
769          yValue1=(rowUpper1-xValue*alpha[1])/element1;
770#ifdef PRINT_VALUES
771          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
772                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
773#endif
774          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
775          // cost>0 so will be at lower
776          //double yValueAtBound1=newLower;
777          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
778          lowerX=CoinMax(lowerX,newLower-1.0e-12*fabs(newLower));
779          upperX=CoinMin(upperX,newUpper+1.0e-12*fabs(newUpper));
780          // Now make duplicate row
781          // keep row 0 so need to adjust costs so same
782#ifdef PRINT_VALUES
783          printf("Costs for x %g,%g,%g are %g,%g,%g\n",
784                 xValueEqual-1.0,xValueEqual,xValueEqual+1.0,
785                 costEqual-slope[0],costEqual,costEqual+slope[1]);
786#endif
787          double costOther=cost[otherCol]+slope[1];
788          double costThis=cost[icol]+slope[1]*(element0/alpha[0]);
789          xValue=xValueEqual;
790          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
791          double thisOffset=costEqual-(costOther*xValue+costThis*yValue0);
792          offset += thisOffset;
793#ifdef PRINT_VALUES
794          printf("new cost at equal %g\n",costOther*xValue+costThis*yValue0+thisOffset);
795#endif
796          xValue=xValueEqual-1.0;
797          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
798#ifdef PRINT_VALUES
799          printf("new cost at -1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
800#endif
801          assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual-slope[0]))<1.0e-5);
802          xValue=xValueEqual+1.0;
803          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
804#ifdef PRINT_VALUES
805          printf("new cost at +1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
806#endif
807          assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual+slope[1]))<1.0e-5);
808          numberChanged++;
809          //      continue;
810          cost[otherCol] = costOther;
811          cost[icol] = costThis;
812          clo[icol]=lowerX;
813          cup[icol]=upperX;
814          int startCol[2];
815          int endCol[2];
816          startCol[0]=mcstrt[icol];
817          endCol[0]=startCol[0]+2;
818          startCol[1]=mcstrt[otherCol];
819          endCol[1]=startCol[1]+hincol[otherCol];
820          double values[2]={0.0,0.0};
821          for (int k=0;k<2;k++) {
822            for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
823              if (hrow[i]==row0)
824                values[k]=colels[i];
825            }
826            for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
827              if (hrow[i]==row1)
828                colels[i]=values[k];
829            }
830          }
831          for (CoinBigIndex i=mrstrt[row1];i<mrstrt[row1]+2;i++) {
832            if (hcol[i]==icol)
833              rowels[i]=values[0];
834            else
835              rowels[i]=values[1];
836          }
837        }
838      }
839    }
840  }
841#if ABC_NORMAL_DEBUG>0
842  if (offset)
843    printf("Cost offset %g\n",offset);
844#endif
845  return numberChanged;
846}
847//#define COIN_PRESOLVE_BUG
848#ifdef COIN_PRESOLVE_BUG
849static int counter=1000000;
850static int startEmptyRows=0;
851static int startEmptyColumns=0;
852static bool break2(CoinPresolveMatrix *prob)
853{
854  int droppedRows = prob->countEmptyRows() - startEmptyRows ;
855  int droppedColumns =  prob->countEmptyCols() - startEmptyColumns;
856  startEmptyRows=prob->countEmptyRows();
857  startEmptyColumns=prob->countEmptyCols();
858  printf("Dropped %d rows and %d columns - current empty %d, %d\n",droppedRows,
859         droppedColumns,startEmptyRows,startEmptyColumns);
860  counter--;
861  if (!counter) {
862    printf("skipping next and all\n");
863  }
864  return (counter<=0);
865}
866#define possibleBreak if (break2(prob)) break
867#define possibleSkip  if (!break2(prob))
868#else
869#define possibleBreak
870#define possibleSkip
871#endif
872// This is the presolve loop.
873// It is a separate virtual function so that it can be easily
874// customized by subclassing CoinPresolve.
875const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
876{
877     // Messages
878     CoinMessages messages = CoinMessage(prob->messages().language());
879     paction_ = 0;
880     prob->maxSubstLevel_ = 3 ;
881#ifndef PRESOLVE_DETAIL
882     if (prob->tuning_) {
883#endif
884       int numberEmptyRows=0;
885       for ( int i=0;i<prob->nrows_;i++) {
886         if (!prob->hinrow_[i]) {
887           PRESOLVE_DETAIL_PRINT(printf("pre_empty row %d\n",i));
888           //printf("pre_empty row %d\n",i);
889           numberEmptyRows++;
890         }
891       }
892       int numberEmptyCols=0;
893       for ( int i=0;i<prob->ncols_;i++) {
894         if (!prob->hincol_[i]) {
895           PRESOLVE_DETAIL_PRINT(printf("pre_empty col %d\n",i));
896           //printf("pre_empty col %d\n",i);
897           numberEmptyCols++;
898         }
899       }
900       printf("CoinPresolve initial state %d empty rows and %d empty columns\n",
901              numberEmptyRows,numberEmptyCols);
902#ifndef PRESOLVE_DETAIL
903     }
904#endif
905     prob->status_ = 0; // say feasible
906     paction_ = make_fixed(prob, paction_);
907     paction_ = testRedundant(prob,paction_) ;
908     // if integers then switch off dual stuff
909     // later just do individually
910     bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
911     // but allow in some cases
912     if ((presolveActions_ & 512) != 0)
913          doDualStuff = true;
914     if (prob->anyProhibited())
915          doDualStuff = false;
916     if (!doDual())
917          doDualStuff = false;
918#if     PRESOLVE_CONSISTENCY
919//  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
920     presolve_links_ok(prob, false, true) ;
921#endif
922
923     if (!prob->status_) {
924          bool slackSingleton = doSingletonColumn();
925          slackSingleton = true;
926          const bool slackd = doSingleton();
927          const bool doubleton = doDoubleton();
928          const bool tripleton = doTripleton();
929          //#define NO_FORCING
930#ifndef NO_FORCING
931          const bool forcing = doForcing();
932#endif
933          const bool ifree = doImpliedFree();
934          const bool zerocost = doTighten();
935          const bool dupcol = doDupcol();
936          const bool duprow = doDuprow();
937          const bool dual = doDualStuff;
938
939          // some things are expensive so just do once (normally)
940
941          int i;
942          // say look at all
943          if (!prob->anyProhibited()) {
944               for (i = 0; i < nrows_; i++)
945                    prob->rowsToDo_[i] = i;
946               prob->numberRowsToDo_ = nrows_;
947               for (i = 0; i < ncols_; i++)
948                    prob->colsToDo_[i] = i;
949               prob->numberColsToDo_ = ncols_;
950          } else {
951               // some stuff must be left alone
952               prob->numberRowsToDo_ = 0;
953               for (i = 0; i < nrows_; i++)
954                    if (!prob->rowProhibited(i))
955                         prob->rowsToDo_[prob->numberRowsToDo_++] = i;
956               prob->numberColsToDo_ = 0;
957               for (i = 0; i < ncols_; i++)
958                    if (!prob->colProhibited(i))
959                         prob->colsToDo_[prob->numberColsToDo_++] = i;
960          }
961
962            // transfer costs (may want to do it in OsiPresolve)
963            // need a transfer back at end of postsolve transferCosts(prob);
964
965          int iLoop;
966#if     PRESOLVE_DEBUG
967          check_sol(prob, 1.0e0);
968#endif
969          if (dupcol) {
970               // maybe allow integer columns to be checked
971               if ((presolveActions_ & 512) != 0)
972                    prob->setPresolveOptions(prob->presolveOptions() | 1);
973               possibleSkip;
974               paction_ = dupcol_action::presolve(prob, paction_);
975          }
976#ifdef ABC_INHERIT
977          if (doTwoxTwo()) {
978            possibleSkip;
979            paction_ = twoxtwo_action::presolve(prob, paction_);
980          }
981#endif
982          if (duprow) {
983            possibleSkip;
984            int nTightened=tightenDoubletons2(prob);
985            if (nTightened)
986              printf("%d doubletons tightened\n",nTightened);
987            paction_ = duprow_action::presolve(prob, paction_);
988          }
989          if (doGubrow()) {
990            possibleSkip;
991               paction_ = gubrow_action::presolve(prob, paction_);
992          }
993
994          if ((presolveActions_ & 16384) != 0)
995               prob->setPresolveOptions(prob->presolveOptions() | 16384);
996          // Check number rows dropped
997          int lastDropped = 0;
998          prob->pass_ = 0;
999#ifdef ABC_INHERIT
1000          int numberRowsStart=nrows_-prob->countEmptyRows();
1001          int numberColumnsStart=ncols_-prob->countEmptyCols();
1002          int numberRowsLeft=numberRowsStart;
1003          int numberColumnsLeft=numberColumnsStart;
1004          bool lastPassWasGood=true;
1005#if ABC_NORMAL_DEBUG
1006          printf("Original rows,columns %d,%d starting first pass with %d,%d\n", 
1007                 nrows_,ncols_,numberRowsLeft,numberColumnsLeft);
1008#endif
1009#endif
1010          if (numberPasses_<=5)
1011              prob->presolveOptions_ |= 0x10000; // say more lightweight
1012          for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
1013               // See if we want statistics
1014               if ((presolveActions_ & 0x80000000) != 0)
1015                    printf("Starting major pass %d after %g seconds\n", iLoop + 1, CoinCpuTime() - prob->startTime_);
1016#ifdef PRESOLVE_SUMMARY
1017               printf("Starting major pass %d\n", iLoop + 1);
1018#endif
1019               const CoinPresolveAction * const paction0 = paction_;
1020               // look for substitutions with no fill
1021               //#define IMPLIED 3
1022#ifdef IMPLIED
1023               int fill_level = 3;
1024#define IMPLIED2 99
1025#if IMPLIED!=3
1026#if IMPLIED>2&&IMPLIED<11
1027               fill_level = IMPLIED;
1028               COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
1029#endif
1030#if IMPLIED>11&&IMPLIED<21
1031               fill_level = -(IMPLIED - 10);
1032               COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
1033#endif
1034#endif
1035#else
1036               int fill_level = prob->maxSubstLevel_;
1037#endif
1038               int whichPass = 0;
1039               while (1) {
1040                    whichPass++;
1041                    prob->pass_++;
1042                    const CoinPresolveAction * const paction1 = paction_;
1043
1044                    if (slackd) {
1045                         bool notFinished = true;
1046                         while (notFinished) {
1047                           possibleBreak;
1048                              paction_ = slack_doubleton_action::presolve(prob, paction_,
1049                                         notFinished);
1050                         }
1051                         if (prob->status_)
1052                              break;
1053                    }
1054                    if (dual && whichPass == 1) {
1055                         // this can also make E rows so do one bit here
1056                      possibleBreak;
1057                         paction_ = remove_dual_action::presolve(prob, paction_);
1058                         if (prob->status_)
1059                              break;
1060                    }
1061
1062                    if (doubleton) {
1063                      possibleBreak;
1064                         paction_ = doubleton_action::presolve(prob, paction_);
1065                         if (prob->status_)
1066                              break;
1067                    }
1068                    if (tripleton) {
1069                      possibleBreak;
1070                         paction_ = tripleton_action::presolve(prob, paction_);
1071                         if (prob->status_)
1072                              break;
1073                    }
1074
1075                    if (zerocost) {
1076                      possibleBreak;
1077                         paction_ = do_tighten_action::presolve(prob, paction_);
1078                         if (prob->status_)
1079                              break;
1080                    }
1081#ifndef NO_FORCING
1082                    if (forcing) {
1083                      possibleBreak;
1084                         paction_ = forcing_constraint_action::presolve(prob, paction_);
1085                         if (prob->status_)
1086                              break;
1087                    }
1088#endif
1089
1090                    if (ifree && (whichPass % 5) == 1) {
1091                      possibleBreak;
1092                         paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1093                         if (prob->status_)
1094                              break;
1095                    }
1096
1097#if     PRESOLVE_DEBUG
1098                    check_sol(prob, 1.0e0);
1099#endif
1100
1101#if     PRESOLVE_CONSISTENCY
1102//      presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
1103//                        prob->nrows_);
1104                    presolve_links_ok(prob, false, true) ;
1105#endif
1106
1107//#if   PRESOLVE_DEBUG
1108//      presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
1109//                        prob->ncols_);
1110//#endif
1111//#if   PRESOLVE_CONSISTENCY
1112//      prob->consistent();
1113//#endif
1114#if     PRESOLVE_CONSISTENCY
1115                    presolve_no_zeros(prob, true, false) ;
1116                    presolve_consistent(prob, true) ;
1117#endif
1118
1119                    {
1120                      // set up for next pass
1121                      // later do faster if many changes i.e. memset and memcpy
1122                      const int * count = prob->hinrow_;
1123                      const int * nextToDo = prob->nextRowsToDo_;
1124                      int * toDo = prob->rowsToDo_;
1125                      int nNext = prob->numberNextRowsToDo_;
1126                      int n = 0;
1127                      for (int i = 0; i < nNext; i++) {
1128                        int index = nextToDo[i];
1129                        prob->unsetRowChanged(index);
1130                        if (count[index]) 
1131                          toDo[n++] = index;
1132                      }
1133                      prob->numberRowsToDo_ = n;
1134                      prob->numberNextRowsToDo_ = 0;
1135                      count = prob->hincol_;
1136                      nextToDo = prob->nextColsToDo_;
1137                      toDo = prob->colsToDo_;
1138                      nNext = prob->numberNextColsToDo_;
1139                      n = 0;
1140                      for (int i = 0; i < nNext; i++) {
1141                        int index = nextToDo[i];
1142                        prob->unsetColChanged(index);
1143                        if (count[index]) 
1144                          toDo[n++] = index;
1145                      }
1146                      prob->numberColsToDo_ = n;
1147                      prob->numberNextColsToDo_ = 0;
1148                    }
1149                    if (paction_ == paction1 && fill_level > 0)
1150                         break;
1151               }
1152               // say look at all
1153               int i;
1154               if (!prob->anyProhibited()) {
1155                 const int * count = prob->hinrow_;
1156                 int * toDo = prob->rowsToDo_;
1157                 int n = 0;
1158                 for (int i = 0; i < nrows_; i++) {
1159                   prob->unsetRowChanged(i);
1160                   if (count[i]) 
1161                     toDo[n++] = i;
1162                 }
1163                 prob->numberRowsToDo_ = n;
1164                 prob->numberNextRowsToDo_ = 0;
1165                 count = prob->hincol_;
1166                 toDo = prob->colsToDo_;
1167                 n = 0;
1168                 for (int i = 0; i < ncols_; i++) {
1169                   prob->unsetColChanged(i);
1170                   if (count[i]) 
1171                     toDo[n++] = i;
1172                 }
1173                 prob->numberColsToDo_ = n;
1174                 prob->numberNextColsToDo_ = 0;
1175               } else {
1176                    // some stuff must be left alone
1177                    prob->numberRowsToDo_ = 0;
1178                    for (i = 0; i < nrows_; i++)
1179                         if (!prob->rowProhibited(i))
1180                              prob->rowsToDo_[prob->numberRowsToDo_++] = i;
1181                    prob->numberColsToDo_ = 0;
1182                    for (i = 0; i < ncols_; i++)
1183                         if (!prob->colProhibited(i))
1184                              prob->colsToDo_[prob->numberColsToDo_++] = i;
1185               }
1186               // now expensive things
1187               // this caused world.mps to run into numerical difficulties
1188#ifdef PRESOLVE_SUMMARY
1189               printf("Starting expensive\n");
1190#endif
1191
1192               if (dual) {
1193                    int itry;
1194                    for (itry = 0; itry < 5; itry++) {
1195                      possibleBreak;
1196                         paction_ = remove_dual_action::presolve(prob, paction_);
1197                         if (prob->status_)
1198                              break;
1199                         const CoinPresolveAction * const paction2 = paction_;
1200                         if (ifree) {
1201#ifdef IMPLIED
1202#if IMPLIED2 ==0
1203                              int fill_level = 0; // switches off substitution
1204#elif IMPLIED2!=99
1205                              int fill_level = IMPLIED2;
1206#endif
1207#endif
1208                              if ((itry & 1) == 0) {
1209                                possibleBreak;
1210                                   paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1211                              }
1212                              if (prob->status_)
1213                                   break;
1214                         }
1215                         if (paction_ == paction2)
1216                              break;
1217                    }
1218               } else if (ifree) {
1219                    // just free
1220#ifdef IMPLIED
1221#if IMPLIED2 ==0
1222                    int fill_level = 0; // switches off substitution
1223#elif IMPLIED2!=99
1224                    int fill_level = IMPLIED2;
1225#endif
1226#endif
1227                    possibleBreak;
1228                    paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1229                    if (prob->status_)
1230                         break;
1231               }
1232#if     PRESOLVE_DEBUG
1233               check_sol(prob, 1.0e0);
1234#endif
1235               if (dupcol) {
1236                    // maybe allow integer columns to be checked
1237                    if ((presolveActions_ & 512) != 0)
1238                         prob->setPresolveOptions(prob->presolveOptions() | 1);
1239                    possibleBreak;
1240                    paction_ = dupcol_action::presolve(prob, paction_);
1241                    if (prob->status_)
1242                         break;
1243               }
1244#if     PRESOLVE_DEBUG
1245               check_sol(prob, 1.0e0);
1246#endif
1247
1248               if (duprow) {
1249                 possibleBreak;
1250                    paction_ = duprow_action::presolve(prob, paction_);
1251                    if (prob->status_)
1252                         break;
1253               }
1254               // Marginally slower on netlib if this call is enabled.
1255               // paction_ = testRedundant(prob,paction_) ;
1256#if     PRESOLVE_DEBUG
1257               check_sol(prob, 1.0e0);
1258#endif
1259               bool stopLoop = false;
1260               {
1261                    int * hinrow = prob->hinrow_;
1262                    int numberDropped = 0;
1263                    for (i = 0; i < nrows_; i++)
1264                         if (!hinrow[i])
1265                              numberDropped++;
1266
1267                    prob->messageHandler()->message(COIN_PRESOLVE_PASS,
1268                                                    messages)
1269                              << numberDropped << iLoop + 1
1270                              << CoinMessageEol;
1271                    //printf("%d rows dropped after pass %d\n",numberDropped,
1272                    //     iLoop+1);
1273                    if (numberDropped == lastDropped)
1274                         stopLoop = true;
1275                    else
1276                         lastDropped = numberDropped;
1277               }
1278               // Do this here as not very loopy
1279               if (slackSingleton) {
1280                    // On most passes do not touch costed slacks
1281                    if (paction_ != paction0 && !stopLoop) {
1282                      possibleBreak;
1283                         paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
1284                    } else {
1285                         // do costed if Clp (at end as ruins rest of presolve)
1286                      possibleBreak;
1287                         paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
1288                         stopLoop = true;
1289                    }
1290               }
1291#if     PRESOLVE_DEBUG
1292               check_sol(prob, 1.0e0);
1293#endif
1294               if (paction_ == paction0 || stopLoop)
1295                    break;
1296#ifdef ABC_INHERIT
1297               // see whether to stop anyway
1298               int numberRowsNow=nrows_-prob->countEmptyRows();
1299               int numberColumnsNow=ncols_-prob->countEmptyCols();
1300#if ABC_NORMAL_DEBUG
1301               printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n", 
1302                      nrows_,ncols_,numberRowsLeft,numberColumnsLeft,iLoop+1,numberRowsNow,
1303                      numberColumnsNow);
1304#endif
1305               int rowsDeleted=numberRowsLeft-numberRowsNow;
1306               int columnsDeleted=numberColumnsLeft-numberColumnsNow;
1307               if (iLoop>15) {
1308                 if (rowsDeleted*100<numberRowsStart&&
1309                     columnsDeleted*100<numberColumnsStart)
1310                   break;
1311                 lastPassWasGood=true;
1312               } else if (rowsDeleted*100<numberRowsStart&&rowsDeleted<500&&
1313                          columnsDeleted*100<numberColumnsStart&&columnsDeleted<500) {
1314                 if (!lastPassWasGood)
1315                   break;
1316                 else
1317                   lastPassWasGood=false;
1318               } else {
1319                 lastPassWasGood=true;
1320               }
1321               numberRowsLeft=numberRowsNow;
1322               numberColumnsLeft=numberColumnsNow;
1323#endif
1324          }
1325     }
1326     prob->presolveOptions_ &= ~0x10000;
1327     if (!prob->status_) {
1328          paction_ = drop_zero_coefficients(prob, paction_);
1329#if     PRESOLVE_DEBUG
1330          check_sol(prob, 1.0e0);
1331#endif
1332
1333          paction_ = drop_empty_cols_action::presolve(prob, paction_);
1334          paction_ = drop_empty_rows_action::presolve(prob, paction_);
1335#if     PRESOLVE_DEBUG
1336          check_sol(prob, 1.0e0);
1337#endif
1338     }
1339
1340     if (prob->status_) {
1341          if (prob->status_ == 1)
1342               prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
1343                                               messages)
1344                         << prob->feasibilityTolerance_
1345                         << CoinMessageEol;
1346          else if (prob->status_ == 2)
1347               prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
1348                                               messages)
1349                         << CoinMessageEol;
1350          else
1351               prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
1352                                               messages)
1353                         << CoinMessageEol;
1354          // get rid of data
1355          destroyPresolve();
1356     }
1357     return (paction_);
1358}
1359
1360void check_djs(CoinPostsolveMatrix *prob);
1361
1362
1363// We could have implemented this by having each postsolve routine
1364// directly call the next one, but this may make it easier to add debugging checks.
1365void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
1366{
1367     {
1368          // Check activities
1369          double *colels        = prob.colels_;
1370          int *hrow             = prob.hrow_;
1371          CoinBigIndex *mcstrt          = prob.mcstrt_;
1372          int *hincol           = prob.hincol_;
1373          int *link             = prob.link_;
1374          int ncols             = prob.ncols_;
1375
1376          char *cdone   = prob.cdone_;
1377
1378          double * csol = prob.sol_;
1379          int nrows = prob.nrows_;
1380
1381          int colx;
1382
1383          double * rsol = prob.acts_;
1384          memset(rsol, 0, nrows * sizeof(double));
1385
1386          for (colx = 0; colx < ncols; ++colx) {
1387               if (cdone[colx]) {
1388                    CoinBigIndex k = mcstrt[colx];
1389                    int nx = hincol[colx];
1390                    double solutionValue = csol[colx];
1391                    for (int i = 0; i < nx; ++i) {
1392                         int row = hrow[k];
1393                         double coeff = colels[k];
1394                         k = link[k];
1395                         rsol[row] += solutionValue * coeff;
1396                    }
1397               }
1398          }
1399     }
1400     if (prob.maxmin_<0) {
1401       //for (int i=0;i<presolvedModel_->numberRows();i++)
1402       //prob.rowduals_[i]=-prob.rowduals_[i];
1403       for (int i=0;i<ncols_;i++) {
1404         prob.cost_[i]=-prob.cost_[i];
1405         //prob.rcosts_[i]=-prob.rcosts_[i];
1406       }
1407       prob.maxmin_=1.0;
1408     }
1409     const CoinPresolveAction *paction = paction_;
1410     //#define PRESOLVE_DEBUG 1
1411#if     PRESOLVE_DEBUG
1412     // Check only works after first one
1413     int checkit = -1;
1414#endif
1415
1416     while (paction) {
1417#if PRESOLVE_DEBUG
1418          printf("POSTSOLVING %s\n", paction->name());
1419#endif
1420
1421          paction->postsolve(&prob);
1422
1423#if     PRESOLVE_DEBUG
1424#         if 0
1425          /*
1426            This check fails (on exmip1 (!) in osiUnitTest) because clp
1427            enters postsolve with a solution that seems to have incorrect
1428            status for a logical. You can see similar behaviour with
1429            column status --- incorrect entering postsolve.
1430            -- lh, 111207 --
1431          */
1432          {
1433               int nr = 0;
1434               int i;
1435               for (i = 0; i < prob.nrows_; i++) {
1436                    if ((prob.rowstat_[i] & 7) == 1) {
1437                         nr++;
1438                    } else if ((prob.rowstat_[i] & 7) == 2) {
1439                         // at ub
1440                         assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
1441                    } else if ((prob.rowstat_[i] & 7) == 3) {
1442                         // at lb
1443                         assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
1444                    }
1445               }
1446               int nc = 0;
1447               for (i = 0; i < prob.ncols_; i++) {
1448                    if ((prob.colstat_[i] & 7) == 1)
1449                         nc++;
1450               }
1451               printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
1452          }
1453#         endif   // if 0
1454          checkit++;
1455          if (prob.colstat_ && checkit > 0) {
1456               presolve_check_nbasic(&prob) ;
1457               presolve_check_sol(&prob, 2, 2, 1) ;
1458          }
1459#endif
1460          paction = paction->next;
1461#if     PRESOLVE_DEBUG
1462//  check_djs(&prob);
1463          if (checkit > 0)
1464               presolve_check_reduced_costs(&prob) ;
1465#endif
1466     }
1467#if     PRESOLVE_DEBUG
1468     if (prob.colstat_) {
1469          presolve_check_nbasic(&prob) ;
1470          presolve_check_sol(&prob, 2, 2, 1) ;
1471     }
1472#endif
1473#undef PRESOLVE_DEBUG
1474
1475#if     0 && PRESOLVE_DEBUG
1476     for (i = 0; i < ncols0; i++) {
1477          if (!cdone[i]) {
1478               printf("!cdone[%d]\n", i);
1479               abort();
1480          }
1481     }
1482
1483     for (i = 0; i < nrows0; i++) {
1484          if (!rdone[i]) {
1485               printf("!rdone[%d]\n", i);
1486               abort();
1487          }
1488     }
1489
1490
1491     for (i = 0; i < ncols0; i++) {
1492          if (sol[i] < -1e10 || sol[i] > 1e10)
1493               printf("!!!%d %g\n", i, sol[i]);
1494
1495     }
1496
1497
1498#endif
1499
1500#if     0 && PRESOLVE_DEBUG
1501     // debug check:  make sure we ended up with same original matrix
1502     {
1503          int identical = 1;
1504
1505          for (int i = 0; i < ncols0; i++) {
1506               PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
1507               CoinBigIndex kcs0 = &prob->mcstrt0[i];
1508               CoinBigIndex kcs = mcstrt[i];
1509               int n = hincol[i];
1510               for (int k = 0; k < n; k++) {
1511                    CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);
1512
1513                    if (k1 == kcs + n) {
1514                         printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
1515                         abort();
1516                    }
1517
1518                    if (colels[k1] != &prob->dels0[kcs0+k])
1519                         printf("BAD COLEL[%d %d %d]:  %g\n",
1520                                k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
1521
1522                    if (kcs0 + k != k1)
1523                         identical = 0;
1524               }
1525          }
1526          printf("identical? %d\n", identical);
1527     }
1528#endif
1529}
1530
1531
1532
1533
1534
1535
1536
1537
1538static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
1539{
1540     double tol;
1541     if (! si->getDblParam(key, tol)) {
1542          CoinPresolveAction::throwCoinError("getDblParam failed",
1543                                             "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
1544     }
1545     return (tol);
1546}
1547
1548
1549// Assumptions:
1550// 1. nrows>=m.getNumRows()
1551// 2. ncols>=m.getNumCols()
1552//
1553// In presolve, these values are equal.
1554// In postsolve, they may be inequal, since the reduced problem
1555// may be smaller, but we need room for the large problem.
1556// ncols may be larger than si.getNumCols() in postsolve,
1557// this at that point si will be the reduced problem,
1558// but we need to reserve enough space for the original problem.
1559CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
1560          int ncols_in,
1561          int nrows_in,
1562          CoinBigIndex nelems_in,
1563          double bulkRatio)
1564     : ncols_(si->getNumCols()),
1565       nrows_(si->getNumRows()),
1566       nelems_(si->getNumElements()),
1567       ncols0_(ncols_in),
1568       nrows0_(nrows_in),
1569       bulkRatio_(bulkRatio),
1570       mcstrt_(new CoinBigIndex[ncols_in+1]),
1571       hincol_(new int[ncols_in+1]),
1572       cost_(new double[ncols_in]),
1573       clo_(new double[ncols_in]),
1574       cup_(new double[ncols_in]),
1575       rlo_(new double[nrows_in]),
1576       rup_(new double[nrows_in]),
1577       originalColumn_(new int[ncols_in]),
1578       originalRow_(new int[nrows_in]),
1579       ztolzb_(getTolerance(si, ClpPrimalTolerance)),
1580       ztoldj_(getTolerance(si, ClpDualTolerance)),
1581       maxmin_(si->getObjSense()),
1582       sol_(NULL),
1583       rowduals_(NULL),
1584       acts_(NULL),
1585       rcosts_(NULL),
1586       colstat_(NULL),
1587       rowstat_(NULL),
1588       handler_(NULL),
1589       defaultHandler_(false),
1590       messages_()
1591
1592{
1593     bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ * nelems_in);
1594     hrow_  = new int   [bulk0_];
1595     colels_ = new double[bulk0_];
1596     si->getDblParam(ClpObjOffset, originalOffset_);
1597     int ncols = si->getNumCols();
1598     int nrows = si->getNumRows();
1599
1600     setMessageHandler(si->messageHandler()) ;
1601
1602     ClpDisjointCopyN(si->getColLower(), ncols, clo_);
1603     ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
1604     //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
1605     double offset;
1606     ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cost_);
1607     ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
1608     ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
1609     int i;
1610     for (i = 0; i < ncols_in; i++)
1611          originalColumn_[i] = i;
1612     for (i = 0; i < nrows_in; i++)
1613          originalRow_[i] = i;
1614     sol_ = NULL;
1615     rowduals_ = NULL;
1616     acts_ = NULL;
1617
1618     rcosts_ = NULL;
1619     colstat_ = NULL;
1620     rowstat_ = NULL;
1621}
1622
1623// I am not familiar enough with CoinPackedMatrix to be confident
1624// that I will implement a row-ordered version of toColumnOrderedGapFree
1625// properly.
1626static bool isGapFree(const CoinPackedMatrix& matrix)
1627{
1628     const CoinBigIndex * start = matrix.getVectorStarts();
1629     const int * length = matrix.getVectorLengths();
1630     int i = matrix.getSizeVectorLengths() - 1;
1631     // Quick check
1632     if (matrix.getNumElements() == start[i]) {
1633          return true;
1634     } else {
1635          for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
1636               if (start[i+1] - start[i] != length[i])
1637                    break;
1638          }
1639          return (! (i >= 0));
1640     }
1641}
1642#if     PRESOLVE_DEBUG
1643static void matrix_bounds_ok(const double *lo, const double *up, int n)
1644{
1645     int i;
1646     for (i = 0; i < n; i++) {
1647          PRESOLVEASSERT(lo[i] <= up[i]);
1648          PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
1649          PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
1650     }
1651}
1652#endif
1653CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
1654                                       double /*maxmin*/,
1655                                       // end prepost members
1656
1657                                       ClpSimplex * si,
1658
1659                                       // rowrep
1660                                       int nrows_in,
1661                                       CoinBigIndex nelems_in,
1662                                       bool doStatus,
1663                                       double nonLinearValue,
1664                                       double bulkRatio) :
1665
1666     CoinPrePostsolveMatrix(si,
1667                            ncols0_in, nrows_in, nelems_in, bulkRatio),
1668     clink_(new presolvehlink[ncols0_in+1]),
1669     rlink_(new presolvehlink[nrows_in+1]),
1670
1671     dobias_(0.0),
1672
1673
1674     // temporary init
1675     integerType_(new unsigned char[ncols0_in]),
1676     tuning_(false),
1677     startTime_(0.0),
1678     feasibilityTolerance_(0.0),
1679     status_(-1),
1680     colsToDo_(new int [ncols0_in]),
1681     numberColsToDo_(0),
1682     nextColsToDo_(new int[ncols0_in]),
1683     numberNextColsToDo_(0),
1684     rowsToDo_(new int [nrows_in]),
1685     numberRowsToDo_(0),
1686     nextRowsToDo_(new int[nrows_in]),
1687     numberNextRowsToDo_(0),
1688     presolveOptions_(0)
1689{
1690     const int bufsize = bulk0_;
1691
1692     nrows_ = si->getNumRows() ;
1693
1694     // Set up change bits etc
1695     rowChanged_ = new unsigned char[nrows_];
1696     memset(rowChanged_, 0, nrows_);
1697     colChanged_ = new unsigned char[ncols_];
1698     memset(colChanged_, 0, ncols_);
1699     CoinPackedMatrix * m = si->matrix();
1700
1701     // The coefficient matrix is a big hunk of stuff.
1702     // Do the copy here to try to avoid running out of memory.
1703
1704     const CoinBigIndex * start = m->getVectorStarts();
1705     const int * row = m->getIndices();
1706     const double * element = m->getElements();
1707     int icol, nel = 0;
1708     mcstrt_[0] = 0;
1709     ClpDisjointCopyN(m->getVectorLengths(), ncols_,  hincol_);
1710     if (si->getObjSense() < 0.0) {
1711       for (int i=0;i<ncols_;i++)
1712         cost_[i]=-cost_[i];
1713       maxmin_=1.0;
1714     }
1715     for (icol = 0; icol < ncols_; icol++) {
1716          CoinBigIndex j;
1717          for (j = start[icol]; j < start[icol] + hincol_[icol]; j++) {
1718               hrow_[nel] = row[j];
1719               if (fabs(element[j]) > ZTOLDP)
1720                    colels_[nel++] = element[j];
1721          }
1722          mcstrt_[icol+1] = nel;
1723          hincol_[icol] = nel - mcstrt_[icol];
1724     }
1725
1726     // same thing for row rep
1727     CoinPackedMatrix * mRow = new CoinPackedMatrix();
1728     mRow->setExtraGap(0.0);
1729     mRow->setExtraMajor(0.0);
1730     mRow->reverseOrderedCopyOf(*m);
1731     //mRow->removeGaps();
1732     //mRow->setExtraGap(0.0);
1733
1734     // Now get rid of matrix
1735     si->createEmptyMatrix();
1736
1737     double * el = mRow->getMutableElements();
1738     int * ind = mRow->getMutableIndices();
1739     CoinBigIndex * strt = mRow->getMutableVectorStarts();
1740     int * len = mRow->getMutableVectorLengths();
1741     // Do carefully to save memory
1742     rowels_ = new double[bulk0_];
1743     ClpDisjointCopyN(el,      nelems_, rowels_);
1744     mRow->nullElementArray();
1745     delete [] el;
1746     hcol_ = new int[bulk0_];
1747     ClpDisjointCopyN(ind,       nelems_, hcol_);
1748     mRow->nullIndexArray();
1749     delete [] ind;
1750     mrstrt_ = new CoinBigIndex[nrows_in+1];
1751     ClpDisjointCopyN(strt,  nrows_,  mrstrt_);
1752     mRow->nullStartArray();
1753     mrstrt_[nrows_] = nelems_;
1754     delete [] strt;
1755     hinrow_ = new int[nrows_in+1];
1756     ClpDisjointCopyN(len, nrows_,  hinrow_);
1757     if (nelems_ > nel) {
1758          nelems_ = nel;
1759          // Clean any small elements
1760          int irow;
1761          nel = 0;
1762          CoinBigIndex start = 0;
1763          for (irow = 0; irow < nrows_; irow++) {
1764               CoinBigIndex j;
1765               for (j = start; j < start + hinrow_[irow]; j++) {
1766                    hcol_[nel] = hcol_[j];
1767                    if (fabs(rowels_[j]) > ZTOLDP)
1768                         rowels_[nel++] = rowels_[j];
1769               }
1770               start = mrstrt_[irow+1];
1771               mrstrt_[irow+1] = nel;
1772               hinrow_[irow] = nel - mrstrt_[irow];
1773          }
1774     }
1775
1776     delete mRow;
1777     if (si->integerInformation()) {
1778          CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()), ncols_, integerType_);
1779     } else {
1780          ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
1781     }
1782
1783#ifndef SLIM_CLP
1784#ifndef NO_RTTI
1785     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1786#else
1787     ClpQuadraticObjective * quadraticObj = NULL;
1788     if (si->objectiveAsObject()->type() == 2)
1789          quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1790#endif
1791#endif
1792     // Set up prohibited bits if needed
1793     if (nonLinearValue) {
1794          anyProhibited_ = true;
1795          for (icol = 0; icol < ncols_; icol++) {
1796               int j;
1797               bool nonLinearColumn = false;
1798               if (cost_[icol] == nonLinearValue)
1799                    nonLinearColumn = true;
1800               for (j = mcstrt_[icol]; j < mcstrt_[icol+1]; j++) {
1801                    if (colels_[j] == nonLinearValue) {
1802                         nonLinearColumn = true;
1803                         setRowProhibited(hrow_[j]);
1804                    }
1805               }
1806               if (nonLinearColumn)
1807                    setColProhibited(icol);
1808          }
1809#ifndef SLIM_CLP
1810     } else if (quadraticObj) {
1811          CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1812          //const int * columnQuadratic = quadratic->getIndices();
1813          //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1814          const int * columnQuadraticLength = quadratic->getVectorLengths();
1815          //double * quadraticElement = quadratic->getMutableElements();
1816          int numberColumns = quadratic->getNumCols();
1817          anyProhibited_ = true;
1818          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1819               if (columnQuadraticLength[iColumn]) {
1820                    setColProhibited(iColumn);
1821                    //printf("%d prohib\n",iColumn);
1822               }
1823          }
1824#endif
1825     } else {
1826          anyProhibited_ = false;
1827     }
1828
1829     if (doStatus) {
1830          // allow for status and solution
1831          sol_ = new double[ncols_];
1832          CoinMemcpyN(si->primalColumnSolution(), ncols_, sol_);;
1833          acts_ = new double [nrows_];
1834          CoinMemcpyN(si->primalRowSolution(), nrows_, acts_);
1835          if (!si->statusArray())
1836               si->createStatus();
1837          colstat_ = new unsigned char [nrows_+ncols_];
1838          CoinMemcpyN(si->statusArray(),        (nrows_ + ncols_), colstat_);
1839          rowstat_ = colstat_ + ncols_;
1840     }
1841
1842     // the original model's fields are now unneeded - free them
1843
1844     si->resize(0, 0);
1845
1846#if     PRESOLVE_DEBUG
1847     matrix_bounds_ok(rlo_, rup_, nrows_);
1848     matrix_bounds_ok(clo_, cup_, ncols_);
1849#endif
1850
1851#if 0
1852     for (i = 0; i < nrows; ++i)
1853          printf("NR: %6d\n", hinrow[i]);
1854     for (int i = 0; i < ncols; ++i)
1855          printf("NC: %6d\n", hincol[i]);
1856#endif
1857
1858     presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
1859     presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
1860
1861     // this allows last col/row to expand up to bufsize-1 (22);
1862     // this must come after the calls to presolve_prefix
1863     mcstrt_[ncols_] = bufsize - 1;
1864     mrstrt_[nrows_] = bufsize - 1;
1865     // Allocate useful arrays
1866     initializeStuff();
1867
1868#if     PRESOLVE_CONSISTENCY
1869//consistent(false);
1870     presolve_consistent(this, false) ;
1871#endif
1872}
1873
1874// avoid compiler warnings
1875#if PRESOLVE_SUMMARY > 0
1876void CoinPresolveMatrix::update_model(ClpSimplex * si,
1877                                      int nrows0, int ncols0,
1878                                      CoinBigIndex nelems0)
1879#else
1880void CoinPresolveMatrix::update_model(ClpSimplex * si,
1881                                      int /*nrows0*/,
1882                                      int /*ncols0*/,
1883                                      CoinBigIndex /*nelems0*/)
1884#endif
1885{
1886     if (si->getObjSense() < 0.0) {
1887       for (int i=0;i<ncols_;i++)
1888         cost_[i]=-cost_[i];
1889       dobias_=-dobias_;
1890     }
1891     si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1892                     clo_, cup_, cost_, rlo_, rup_);
1893     //delete [] si->integerInformation();
1894     int numberIntegers = 0;
1895     for (int i = 0; i < ncols_; i++) {
1896          if (integerType_[i])
1897               numberIntegers++;
1898     }
1899     if (numberIntegers)
1900          si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
1901     else
1902          si->copyInIntegerInformation(NULL);
1903
1904#if     PRESOLVE_SUMMARY
1905     printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
1906            ncols_, ncols0 - ncols_,
1907            nrows_, nrows0 - nrows_,
1908            si->getNumElements(), nelems0 - si->getNumElements());
1909#endif
1910     si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
1911     if (si->getObjSense() < 0.0) {
1912       // put back
1913       for (int i=0;i<ncols_;i++)
1914         cost_[i]=-cost_[i];
1915       dobias_=-dobias_;
1916       maxmin_=-1.0;
1917     }
1918
1919}
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931////////////////  POSTSOLVE
1932
1933CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
1934          int ncols0_in,
1935          int nrows0_in,
1936          CoinBigIndex nelems0,
1937
1938          double maxmin,
1939          // end prepost members
1940
1941          double *sol_in,
1942          double *acts_in,
1943
1944          unsigned char *colstat_in,
1945          unsigned char *rowstat_in) :
1946     CoinPrePostsolveMatrix(si,
1947                            ncols0_in, nrows0_in, nelems0, 2.0),
1948
1949     free_list_(0),
1950     // link, free_list, maxlink
1951     maxlink_(bulk0_),
1952     link_(new int[/*maxlink*/ bulk0_]),
1953
1954     cdone_(new char[ncols0_]),
1955     rdone_(new char[nrows0_in])
1956
1957{
1958     bulk0_ = maxlink_ ;
1959     nrows_ = si->getNumRows() ;
1960     ncols_ = si->getNumCols() ;
1961
1962     sol_ = sol_in;
1963     rowduals_ = NULL;
1964     acts_ = acts_in;
1965
1966     rcosts_ = NULL;
1967     colstat_ = colstat_in;
1968     rowstat_ = rowstat_in;
1969
1970     // this is the *reduced* model, which is probably smaller
1971     int ncols1 = ncols_ ;
1972     int nrows1 = nrows_ ;
1973
1974     const CoinPackedMatrix * m = si->matrix();
1975
1976     const CoinBigIndex nelemsr = m->getNumElements();
1977     if (m->getNumElements() && !isGapFree(*m)) {
1978          // Odd - gaps
1979          CoinPackedMatrix mm(*m);
1980          mm.removeGaps();
1981          mm.setExtraGap(0.0);
1982
1983          ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
1984          CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
1985          mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
1986          ClpDisjointCopyN(mm.getVectorLengths(), ncols1,  hincol_);
1987          ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
1988          ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
1989     } else {
1990          // No gaps
1991
1992          ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1993          CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
1994          mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
1995          ClpDisjointCopyN(m->getVectorLengths(), ncols1,  hincol_);
1996          ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
1997          ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
1998     }
1999
2000
2001
2002#if     0 && PRESOLVE_DEBUG
2003     presolve_check_costs(model, &colcopy);
2004#endif
2005
2006     // This determines the size of the data structure that contains
2007     // the matrix being postsolved.  Links are taken from the free_list
2008     // to recreate matrix entries that were presolved away,
2009     // and links are added to the free_list when entries created during
2010     // presolve are discarded.  There is never a need to gc this list.
2011     // Naturally, it should contain
2012     // exactly nelems0 entries "in use" when postsolving is done,
2013     // but I don't know whether the matrix could temporarily get
2014     // larger during postsolving.  Substitution into more than two
2015     // rows could do that, in principle.  I am being very conservative
2016     // here by reserving much more than the amount of space I probably need.
2017     // If this guess is wrong, check_free_list may be called.
2018     //  int bufsize = 2*nelems0;
2019
2020     memset(cdone_, -1, ncols0_);
2021     memset(rdone_, -1, nrows0_);
2022
2023     rowduals_ = new double[nrows0_];
2024     ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
2025
2026     rcosts_ = new double[ncols0_];
2027     ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
2028     if (maxmin < 0.0) {
2029          // change so will look as if minimize
2030          int i;
2031          for (i = 0; i < nrows1; i++)
2032               rowduals_[i] = - rowduals_[i];
2033          for (i = 0; i < ncols1; i++) {
2034               rcosts_[i] = - rcosts_[i];
2035          }
2036     }
2037
2038     //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
2039     //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
2040
2041     ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
2042     si->setDblParam(ClpObjOffset, originalOffset_);
2043
2044     for (int j = 0; j < ncols1; j++) {
2045#ifdef COIN_SLOW_PRESOLVE
2046       if (hincol_[j]) {
2047#endif
2048          CoinBigIndex kcs = mcstrt_[j];
2049          CoinBigIndex kce = kcs + hincol_[j];
2050          for (CoinBigIndex k = kcs; k < kce; ++k) {
2051               link_[k] = k + 1;
2052          }
2053          link_[kce-1] = NO_LINK ;
2054#ifdef COIN_SLOW_PRESOLVE
2055       }
2056#endif
2057     }
2058     {
2059          int ml = maxlink_;
2060          for (CoinBigIndex k = nelemsr; k < ml; ++k)
2061               link_[k] = k + 1;
2062          if (ml)
2063               link_[ml-1] = NO_LINK;
2064     }
2065     free_list_ = nelemsr;
2066# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
2067     /*
2068       These are used to track the action of postsolve transforms during debugging.
2069     */
2070     CoinFillN(cdone_, ncols1, PRESENT_IN_REDUCED) ;
2071     CoinZeroN(cdone_ + ncols1, ncols0_in - ncols1) ;
2072     CoinFillN(rdone_, nrows1, PRESENT_IN_REDUCED) ;
2073     CoinZeroN(rdone_ + nrows1, nrows0_in - nrows1) ;
2074# endif
2075}
2076/* This is main part of Presolve */
2077ClpSimplex *
2078ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
2079                                  double feasibilityTolerance,
2080                                  bool keepIntegers,
2081                                  int numberPasses,
2082                                  bool dropNames,
2083                                  bool doRowObjective,
2084                                  const char * prohibitedRows,
2085                                  const char * prohibitedColumns)
2086{
2087     ncols_ = originalModel->getNumCols();
2088     nrows_ = originalModel->getNumRows();
2089     nelems_ = originalModel->getNumElements();
2090     numberPasses_ = numberPasses;
2091
2092     double maxmin = originalModel->getObjSense();
2093     originalModel_ = originalModel;
2094     delete [] originalColumn_;
2095     originalColumn_ = new int[ncols_];
2096     delete [] originalRow_;
2097     originalRow_ = new int[nrows_];
2098     // and fill in case returns early
2099     int i;
2100     for (i = 0; i < ncols_; i++)
2101          originalColumn_[i] = i;
2102     for (i = 0; i < nrows_; i++)
2103          originalRow_[i] = i;
2104     delete [] rowObjective_;
2105     if (doRowObjective) {
2106          rowObjective_ = new double [nrows_];
2107          memset(rowObjective_, 0, nrows_ * sizeof(double));
2108     } else {
2109          rowObjective_ = NULL;
2110     }
2111
2112     // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
2113     int result = -1;
2114
2115     // User may have deleted - its their responsibility
2116     presolvedModel_ = NULL;
2117     // Messages
2118     CoinMessages messages = originalModel->coinMessages();
2119     // Only go round 100 times even if integer preprocessing
2120     int totalPasses = 100;
2121     while (result == -1) {
2122
2123#ifndef CLP_NO_STD
2124          // make new copy
2125          if (saveFile_ == "") {
2126#endif
2127               delete presolvedModel_;
2128#ifndef CLP_NO_STD
2129               // So won't get names
2130               int lengthNames = originalModel->lengthNames();
2131               originalModel->setLengthNames(0);
2132#endif
2133               presolvedModel_ = new ClpSimplex(*originalModel);
2134#ifndef CLP_NO_STD
2135               originalModel->setLengthNames(lengthNames);
2136               presolvedModel_->dropNames();
2137          } else {
2138               presolvedModel_ = originalModel;
2139               if (dropNames)
2140                 presolvedModel_->dropNames();
2141          }
2142#endif
2143
2144          // drop integer information if wanted
2145          if (!keepIntegers)
2146               presolvedModel_->deleteIntegerInformation();
2147          totalPasses--;
2148
2149          double ratio = 2.0;
2150          if (substitution_ > 3)
2151               ratio = substitution_;
2152          else if (substitution_ == 2)
2153               ratio = 1.5;
2154          CoinPresolveMatrix prob(ncols_,
2155                                  maxmin,
2156                                  presolvedModel_,
2157                                  nrows_, nelems_, true, nonLinearValue_, ratio);
2158          if (prohibitedRows) {
2159            prob.setAnyProhibited();
2160            for (int i=0;i<nrows_;i++) {
2161              if (prohibitedRows[i])
2162                prob.setRowProhibited(i);
2163            }
2164          }
2165          if (prohibitedColumns) {
2166            prob.setAnyProhibited();
2167            for (int i=0;i<ncols_;i++) {
2168              if (prohibitedColumns[i])
2169                prob.setColProhibited(i);
2170            }
2171          }
2172          prob.setMaximumSubstitutionLevel(substitution_);
2173          if (doRowObjective)
2174               memset(rowObjective_, 0, nrows_ * sizeof(double));
2175          // See if we want statistics
2176          if ((presolveActions_ & 0x80000000) != 0)
2177               prob.statistics();
2178          // make sure row solution correct
2179          {
2180               double *colels   = prob.colels_;
2181               int *hrow                = prob.hrow_;
2182               CoinBigIndex *mcstrt             = prob.mcstrt_;
2183               int *hincol              = prob.hincol_;
2184               int ncols                = prob.ncols_;
2185
2186
2187               double * csol = prob.sol_;
2188               double * acts = prob.acts_;
2189               int nrows = prob.nrows_;
2190
2191               int colx;
2192
2193               memset(acts, 0, nrows * sizeof(double));
2194
2195               for (colx = 0; colx < ncols; ++colx) {
2196                    double solutionValue = csol[colx];
2197                    for (int i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
2198                         int row = hrow[i];
2199                         double coeff = colels[i];
2200                         acts[row] += solutionValue * coeff;
2201                    }
2202               }
2203          }
2204
2205          // move across feasibility tolerance
2206          prob.feasibilityTolerance_ = feasibilityTolerance;
2207
2208          // Do presolve
2209          paction_ = presolve(&prob);
2210          // Get rid of useful arrays
2211          prob.deleteStuff();
2212
2213          result = 0;
2214
2215          bool fixInfeasibility = (prob.presolveOptions_&16384)!=0;
2216          bool hasSolution = (prob.presolveOptions_&32768)!=0;
2217          if (prob.status_ == 0 && paction_ && (!hasSolution || !fixInfeasibility)) {
2218               // Looks feasible but double check to see if anything slipped through
2219               int n            = prob.ncols_;
2220               double * lo = prob.clo_;
2221               double * up = prob.cup_;
2222               int i;
2223
2224               for (i = 0; i < n; i++) {
2225                    if (up[i] < lo[i]) {
2226                         if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2227                              // infeasible
2228                              prob.status_ = 1;
2229                         } else {
2230                              up[i] = lo[i];
2231                         }
2232                    }
2233               }
2234
2235               n = prob.nrows_;
2236               lo = prob.rlo_;
2237               up = prob.rup_;
2238
2239               for (i = 0; i < n; i++) {
2240                    if (up[i] < lo[i]) {
2241                         if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2242                              // infeasible
2243                              prob.status_ = 1;
2244                         } else {
2245                              up[i] = lo[i];
2246                         }
2247                    }
2248               }
2249          }
2250          if (prob.status_ == 0 && paction_) {
2251               // feasible
2252
2253               prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
2254               // copy status and solution
2255               CoinMemcpyN(          prob.sol_, prob.ncols_, presolvedModel_->primalColumnSolution());
2256               CoinMemcpyN(          prob.acts_, prob.nrows_, presolvedModel_->primalRowSolution());
2257               CoinMemcpyN(          prob.colstat_, prob.ncols_, presolvedModel_->statusArray());
2258               CoinMemcpyN(          prob.rowstat_, prob.nrows_, presolvedModel_->statusArray() + prob.ncols_);
2259               if (fixInfeasibility && hasSolution) {
2260                 // Looks feasible but double check to see if anything slipped through
2261                 int n          = prob.ncols_;
2262                 double * lo = prob.clo_;
2263                 double * up = prob.cup_;
2264                 double * rsol = prob.acts_;
2265                 //memset(prob.acts_,0,prob.nrows_*sizeof(double));
2266                 presolvedModel_->matrix()->times(prob.sol_,rsol);
2267                 int i;
2268                 
2269                 for (i = 0; i < n; i++) {
2270                   double gap=up[i]-lo[i];
2271                   if (rsol[i]<lo[i]-feasibilityTolerance&&fabs(rsol[i]-lo[i])<1.0e-3) {
2272                     lo[i]=rsol[i];
2273                     if (gap<1.0e5)
2274                       up[i]=lo[i]+gap;
2275                   } else if (rsol[i]>up[i]+feasibilityTolerance&&fabs(rsol[i]-up[i])<1.0e-3) {
2276                     up[i]=rsol[i];
2277                     if (gap<1.0e5)
2278                       lo[i]=up[i]-gap;
2279                   }
2280                   if (up[i] < lo[i]) {
2281                     up[i] = lo[i];
2282                   }
2283                 }
2284               }
2285
2286               int n = prob.nrows_;
2287               double * lo = prob.rlo_;
2288               double * up = prob.rup_;
2289
2290               for (i = 0; i < n; i++) {
2291                    if (up[i] < lo[i]) {
2292                         if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2293                              // infeasible
2294                              prob.status_ = 1;
2295                         } else {
2296                              up[i] = lo[i];
2297                         }
2298                    }
2299               }
2300               delete [] prob.sol_;
2301               delete [] prob.acts_;
2302               delete [] prob.colstat_;
2303               prob.sol_ = NULL;
2304               prob.acts_ = NULL;
2305               prob.colstat_ = NULL;
2306
2307               int ncolsNow = presolvedModel_->getNumCols();
2308               CoinMemcpyN(prob.originalColumn_, ncolsNow, originalColumn_);
2309#ifndef SLIM_CLP
2310#ifndef NO_RTTI
2311               ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
2312#else
2313               ClpQuadraticObjective * quadraticObj = NULL;
2314               if (originalModel->objectiveAsObject()->type() == 2)
2315                    quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
2316#endif
2317               if (quadraticObj) {
2318                    // set up for subset
2319                    char * mark = new char [ncols_];
2320                    memset(mark, 0, ncols_);
2321                    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2322                    //const int * columnQuadratic = quadratic->getIndices();
2323                    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
2324                    const int * columnQuadraticLength = quadratic->getVectorLengths();
2325                    //double * quadraticElement = quadratic->getMutableElements();
2326                    int numberColumns = quadratic->getNumCols();
2327                    ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
2328                              ncolsNow,
2329                              originalColumn_);
2330                    // and modify linear and check
2331                    double * linear = newObj->linearObjective();
2332                    CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
2333                    int iColumn;
2334                    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
2335                         if (columnQuadraticLength[iColumn])
2336                              mark[iColumn] = 1;
2337                    }
2338                    // and new
2339                    quadratic = newObj->quadraticObjective();
2340                    columnQuadraticLength = quadratic->getVectorLengths();
2341                    int numberColumns2 = quadratic->getNumCols();
2342                    for ( iColumn = 0; iColumn < numberColumns2; iColumn++) {
2343                         if (columnQuadraticLength[iColumn])
2344                              mark[originalColumn_[iColumn]] = 0;
2345                    }
2346                    presolvedModel_->setObjective(newObj);
2347                    delete newObj;
2348                    // final check
2349                    for ( iColumn = 0; iColumn < numberColumns; iColumn++)
2350                         if (mark[iColumn])
2351                              printf("Quadratic column %d modified - may be okay\n", iColumn);
2352                    delete [] mark;
2353               }
2354#endif
2355               delete [] prob.originalColumn_;
2356               prob.originalColumn_ = NULL;
2357               int nrowsNow = presolvedModel_->getNumRows();
2358               CoinMemcpyN(prob.originalRow_, nrowsNow, originalRow_);
2359               delete [] prob.originalRow_;
2360               prob.originalRow_ = NULL;
2361#ifndef CLP_NO_STD
2362               if (!dropNames && originalModel->lengthNames()) {
2363                    // Redo names
2364                    int iRow;
2365                    std::vector<std::string> rowNames;
2366                    rowNames.reserve(nrowsNow);
2367                    for (iRow = 0; iRow < nrowsNow; iRow++) {
2368                         int kRow = originalRow_[iRow];
2369                         rowNames.push_back(originalModel->rowName(kRow));
2370                    }
2371
2372                    int iColumn;
2373                    std::vector<std::string> columnNames;
2374                    columnNames.reserve(ncolsNow);
2375                    for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
2376                         int kColumn = originalColumn_[iColumn];
2377                         columnNames.push_back(originalModel->columnName(kColumn));
2378                    }
2379                    presolvedModel_->copyNames(rowNames, columnNames);
2380               } else {
2381                    presolvedModel_->setLengthNames(0);
2382               }
2383#endif
2384               if (rowObjective_) {
2385                    int iRow;
2386#ifndef NDEBUG
2387                    int k = -1;
2388#endif
2389                    int nObj = 0;
2390                    for (iRow = 0; iRow < nrowsNow; iRow++) {
2391                         int kRow = originalRow_[iRow];
2392#ifndef NDEBUG
2393                         assert (kRow > k);
2394                         k = kRow;
2395#endif
2396                         rowObjective_[iRow] = rowObjective_[kRow];
2397                         if (rowObjective_[iRow])
2398                              nObj++;
2399                    }
2400                    if (nObj) {
2401                         printf("%d costed slacks\n", nObj);
2402                         presolvedModel_->setRowObjective(rowObjective_);
2403                    }
2404               }
2405               /* now clean up integer variables.  This can modify original
2406                          Don't do if dupcol added columns together */
2407               int i;
2408               const char * information = presolvedModel_->integerInformation();
2409               if ((prob.presolveOptions_ & 0x80000000) == 0 && information) {
2410                    int numberChanges = 0;
2411                    double * lower0 = originalModel_->columnLower();
2412                    double * upper0 = originalModel_->columnUpper();
2413                    double * lower = presolvedModel_->columnLower();
2414                    double * upper = presolvedModel_->columnUpper();
2415                    for (i = 0; i < ncolsNow; i++) {
2416                         if (!information[i])
2417                              continue;
2418                         int iOriginal = originalColumn_[i];
2419                         double lowerValue0 = lower0[iOriginal];
2420                         double upperValue0 = upper0[iOriginal];
2421                         double lowerValue = ceil(lower[i] - 1.0e-5);
2422                         double upperValue = floor(upper[i] + 1.0e-5);
2423                         lower[i] = lowerValue;
2424                         upper[i] = upperValue;
2425                         if (lowerValue > upperValue) {
2426                              numberChanges++;
2427                              presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
2428                                        messages)
2429                                        << iOriginal
2430                                        << lowerValue
2431                                        << upperValue
2432                                        << CoinMessageEol;
2433                              result = 1;
2434                         } else {
2435                              if (lowerValue > lowerValue0 + 1.0e-8) {
2436                                   lower0[iOriginal] = lowerValue;
2437                                   numberChanges++;
2438                              }
2439                              if (upperValue < upperValue0 - 1.0e-8) {
2440                                   upper0[iOriginal] = upperValue;
2441                                   numberChanges++;
2442                              }
2443                         }
2444                    }
2445                    if (numberChanges) {
2446                         presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
2447                                   messages)
2448                                   << numberChanges
2449                                   << CoinMessageEol;
2450                         if (!result && totalPasses > 0) {
2451                              result = -1; // round again
2452                              const CoinPresolveAction *paction = paction_;
2453                              while (paction) {
2454                                   const CoinPresolveAction *next = paction->next;
2455                                   delete paction;
2456                                   paction = next;
2457                              }
2458                              paction_ = NULL;
2459                         }
2460                    }
2461               }
2462          } else if (prob.status_) {
2463               // infeasible or unbounded
2464               result = 1;
2465               // Put status in nelems_!
2466               nelems_ = - prob.status_;
2467               originalModel->setProblemStatus(prob.status_);
2468          } else {
2469               // no changes - model needs restoring after Lou's changes
2470#ifndef CLP_NO_STD
2471               if (saveFile_ == "") {
2472#endif
2473                    delete presolvedModel_;
2474                    presolvedModel_ = new ClpSimplex(*originalModel);
2475                    // but we need to remove gaps
2476                    ClpPackedMatrix* clpMatrix =
2477                         dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
2478                    if (clpMatrix) {
2479                         clpMatrix->getPackedMatrix()->removeGaps();
2480                    }
2481#ifndef CLP_NO_STD
2482               } else {
2483                    presolvedModel_ = originalModel;
2484               }
2485               presolvedModel_->dropNames();
2486#endif
2487
2488               // drop integer information if wanted
2489               if (!keepIntegers)
2490                    presolvedModel_->deleteIntegerInformation();
2491               result = 2;
2492          }
2493     }
2494     if (result == 0 || result == 2) {
2495          int nrowsAfter = presolvedModel_->getNumRows();
2496          int ncolsAfter = presolvedModel_->getNumCols();
2497          CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
2498          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
2499                    messages)
2500                    << nrowsAfter << -(nrows_ - nrowsAfter)
2501                    << ncolsAfter << -(ncols_ - ncolsAfter)
2502                    << nelsAfter << -(nelems_ - nelsAfter)
2503                    << CoinMessageEol;
2504     } else {
2505          destroyPresolve();
2506          if (presolvedModel_ != originalModel_)
2507               delete presolvedModel_;
2508          presolvedModel_ = NULL;
2509     }
2510     return presolvedModel_;
2511}
2512
2513
Note: See TracBrowser for help on using the repository browser.