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

Last change on this file since 1879 was 1879, checked in by forrest, 7 years ago

fix a few problems with Aboca

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 83.3 KB
Line 
1/* $Id: ClpPresolve.cpp 1879 2012-09-13 16:31:33Z forrest $ */
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#ifndef PRESOLVE_DETAIL
881     if (prob->tuning_) {
882#endif
883       int numberEmptyRows=0;
884       for ( int i=0;i<prob->nrows_;i++) {
885         if (!prob->hinrow_[i]) {
886           PRESOLVE_DETAIL_PRINT(printf("pre_empty row %d\n",i));
887           //printf("pre_empty row %d\n",i);
888           numberEmptyRows++;
889         }
890       }
891       int numberEmptyCols=0;
892       for ( int i=0;i<prob->ncols_;i++) {
893         if (!prob->hincol_[i]) {
894           PRESOLVE_DETAIL_PRINT(printf("pre_empty col %d\n",i));
895           //printf("pre_empty col %d\n",i);
896           numberEmptyCols++;
897         }
898       }
899       printf("CoinPresolve initial state %d empty rows and %d empty columns\n",
900              numberEmptyRows,numberEmptyCols);
901#ifndef PRESOLVE_DETAIL
902     }
903#endif
904     prob->status_ = 0; // say feasible
905     paction_ = make_fixed(prob, paction_);
906     // if integers then switch off dual stuff
907     // later just do individually
908     bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
909     // but allow in some cases
910     if ((presolveActions_ & 512) != 0)
911          doDualStuff = true;
912     if (prob->anyProhibited())
913          doDualStuff = false;
914     if (!doDual())
915          doDualStuff = false;
916#if     PRESOLVE_CONSISTENCY
917//  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
918     presolve_links_ok(prob, false, true) ;
919#endif
920
921     if (!prob->status_) {
922          bool slackSingleton = doSingletonColumn();
923          slackSingleton = true;
924          const bool slackd = doSingleton();
925          const bool doubleton = doDoubleton();
926          const bool tripleton = doTripleton();
927          //#define NO_FORCING
928#ifndef NO_FORCING
929          const bool forcing = doForcing();
930#endif
931          const bool ifree = doImpliedFree();
932          const bool zerocost = doTighten();
933          const bool dupcol = doDupcol();
934          const bool duprow = doDuprow();
935          const bool dual = doDualStuff;
936
937          // some things are expensive so just do once (normally)
938
939          int i;
940          // say look at all
941          if (!prob->anyProhibited()) {
942               for (i = 0; i < nrows_; i++)
943                    prob->rowsToDo_[i] = i;
944               prob->numberRowsToDo_ = nrows_;
945               for (i = 0; i < ncols_; i++)
946                    prob->colsToDo_[i] = i;
947               prob->numberColsToDo_ = ncols_;
948          } else {
949               // some stuff must be left alone
950               prob->numberRowsToDo_ = 0;
951               for (i = 0; i < nrows_; i++)
952                    if (!prob->rowProhibited(i))
953                         prob->rowsToDo_[prob->numberRowsToDo_++] = i;
954               prob->numberColsToDo_ = 0;
955               for (i = 0; i < ncols_; i++)
956                    if (!prob->colProhibited(i))
957                         prob->colsToDo_[prob->numberColsToDo_++] = i;
958          }
959
960            // transfer costs (may want to do it in OsiPresolve)
961            // need a transfer back at end of postsolve transferCosts(prob);
962
963          int iLoop;
964#if     PRESOLVE_DEBUG
965          check_sol(prob, 1.0e0);
966#endif
967          if (dupcol) {
968               // maybe allow integer columns to be checked
969               if ((presolveActions_ & 512) != 0)
970                    prob->setPresolveOptions(prob->presolveOptions() | 1);
971               possibleSkip;
972               paction_ = dupcol_action::presolve(prob, paction_);
973          }
974#ifdef ABC_INHERIT
975          if (doTwoxTwo()) {
976            possibleSkip;
977            paction_ = twoxtwo_action::presolve(prob, paction_);
978          }
979#endif
980          if (duprow) {
981            possibleSkip;
982            int nTightened=tightenDoubletons2(prob);
983            if (nTightened)
984              printf("%d doubletons tightened\n",nTightened);
985            paction_ = duprow_action::presolve(prob, paction_);
986          }
987          if (doGubrow()) {
988            possibleSkip;
989               paction_ = gubrow_action::presolve(prob, paction_);
990          }
991
992          if ((presolveActions_ & 16384) != 0)
993               prob->setPresolveOptions(prob->presolveOptions() | 16384);
994          // Check number rows dropped
995          int lastDropped = 0;
996          prob->pass_ = 0;
997#ifdef ABC_INHERIT
998          int numberRowsStart=nrows_-prob->countEmptyRows();
999          int numberColumnsStart=ncols_-prob->countEmptyCols();
1000          int numberRowsLeft=numberRowsStart;
1001          int numberColumnsLeft=numberColumnsStart;
1002          bool lastPassWasGood=true;
1003#if ABC_NORMAL_DEBUG
1004          printf("Original rows,columns %d,%d starting first pass with %d,%d\n", 
1005                 nrows_,ncols_,numberRowsLeft,numberColumnsLeft);
1006#endif
1007#endif
1008          if (numberPasses_<=5)
1009              prob->presolveOptions_ |= 0x10000; // say more lightweight
1010          for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
1011               // See if we want statistics
1012               if ((presolveActions_ & 0x80000000) != 0)
1013                    printf("Starting major pass %d after %g seconds\n", iLoop + 1, CoinCpuTime() - prob->startTime_);
1014#ifdef PRESOLVE_SUMMARY
1015               printf("Starting major pass %d\n", iLoop + 1);
1016#endif
1017               const CoinPresolveAction * const paction0 = paction_;
1018               // look for substitutions with no fill
1019               //#define IMPLIED 3
1020#ifdef IMPLIED
1021               int fill_level = 3;
1022#define IMPLIED2 99
1023#if IMPLIED!=3
1024#if IMPLIED>2&&IMPLIED<11
1025               fill_level = IMPLIED;
1026               COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
1027#endif
1028#if IMPLIED>11&&IMPLIED<21
1029               fill_level = -(IMPLIED - 10);
1030               COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
1031#endif
1032#endif
1033#else
1034               int fill_level = 2;
1035#endif
1036               int whichPass = 0;
1037               while (1) {
1038                    whichPass++;
1039                    prob->pass_++;
1040                    const CoinPresolveAction * const paction1 = paction_;
1041
1042                    if (slackd) {
1043                         bool notFinished = true;
1044                         while (notFinished) {
1045                           possibleBreak;
1046                              paction_ = slack_doubleton_action::presolve(prob, paction_,
1047                                         notFinished);
1048                         }
1049                         if (prob->status_)
1050                              break;
1051                    }
1052                    if (dual && whichPass == 1) {
1053                         // this can also make E rows so do one bit here
1054                      possibleBreak;
1055                         paction_ = remove_dual_action::presolve(prob, paction_);
1056                         if (prob->status_)
1057                              break;
1058                    }
1059
1060                    if (doubleton) {
1061                      possibleBreak;
1062                         paction_ = doubleton_action::presolve(prob, paction_);
1063                         if (prob->status_)
1064                              break;
1065                    }
1066                    if (tripleton) {
1067                      possibleBreak;
1068                         paction_ = tripleton_action::presolve(prob, paction_);
1069                         if (prob->status_)
1070                              break;
1071                    }
1072
1073                    if (zerocost) {
1074                      possibleBreak;
1075                         paction_ = do_tighten_action::presolve(prob, paction_);
1076                         if (prob->status_)
1077                              break;
1078                    }
1079#ifndef NO_FORCING
1080                    if (forcing) {
1081                      possibleBreak;
1082                         paction_ = forcing_constraint_action::presolve(prob, paction_);
1083                         if (prob->status_)
1084                              break;
1085                    }
1086#endif
1087
1088                    if (ifree && (whichPass % 5) == 1) {
1089                      possibleBreak;
1090                         paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1091                         if (prob->status_)
1092                              break;
1093                    }
1094
1095#if     PRESOLVE_DEBUG
1096                    check_sol(prob, 1.0e0);
1097#endif
1098
1099#if     PRESOLVE_CONSISTENCY
1100//      presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
1101//                        prob->nrows_);
1102                    presolve_links_ok(prob, false, true) ;
1103#endif
1104
1105//#if   PRESOLVE_DEBUG
1106//      presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
1107//                        prob->ncols_);
1108//#endif
1109//#if   PRESOLVE_CONSISTENCY
1110//      prob->consistent();
1111//#endif
1112#if     PRESOLVE_CONSISTENCY
1113                    presolve_no_zeros(prob, true, false) ;
1114                    presolve_consistent(prob, true) ;
1115#endif
1116
1117                    {
1118                      // set up for next pass
1119                      // later do faster if many changes i.e. memset and memcpy
1120                      const int * count = prob->hinrow_;
1121                      const int * nextToDo = prob->nextRowsToDo_;
1122                      int * toDo = prob->rowsToDo_;
1123                      int nNext = prob->numberNextRowsToDo_;
1124                      int n = 0;
1125                      for (int i = 0; i < nNext; i++) {
1126                        int index = nextToDo[i];
1127                        prob->unsetRowChanged(index);
1128                        if (count[index]) 
1129                          toDo[n++] = index;
1130                      }
1131                      prob->numberRowsToDo_ = n;
1132                      prob->numberNextRowsToDo_ = 0;
1133                      count = prob->hincol_;
1134                      nextToDo = prob->nextColsToDo_;
1135                      toDo = prob->colsToDo_;
1136                      nNext = prob->numberNextColsToDo_;
1137                      n = 0;
1138                      for (int i = 0; i < nNext; i++) {
1139                        int index = nextToDo[i];
1140                        prob->unsetColChanged(index);
1141                        if (count[index]) 
1142                          toDo[n++] = index;
1143                      }
1144                      prob->numberColsToDo_ = n;
1145                      prob->numberNextColsToDo_ = 0;
1146                    }
1147                    if (paction_ == paction1 && fill_level > 0)
1148                         break;
1149               }
1150               // say look at all
1151               int i;
1152               if (!prob->anyProhibited()) {
1153                 const int * count = prob->hinrow_;
1154                 int * toDo = prob->rowsToDo_;
1155                 int n = 0;
1156                 for (int i = 0; i < nrows_; i++) {
1157                   prob->unsetRowChanged(i);
1158                   if (count[i]) 
1159                     toDo[n++] = i;
1160                 }
1161                 prob->numberRowsToDo_ = n;
1162                 prob->numberNextRowsToDo_ = 0;
1163                 count = prob->hincol_;
1164                 toDo = prob->colsToDo_;
1165                 n = 0;
1166                 for (int i = 0; i < ncols_; i++) {
1167                   prob->unsetColChanged(i);
1168                   if (count[i]) 
1169                     toDo[n++] = i;
1170                 }
1171                 prob->numberColsToDo_ = n;
1172                 prob->numberNextColsToDo_ = 0;
1173               } else {
1174                    // some stuff must be left alone
1175                    prob->numberRowsToDo_ = 0;
1176                    for (i = 0; i < nrows_; i++)
1177                         if (!prob->rowProhibited(i))
1178                              prob->rowsToDo_[prob->numberRowsToDo_++] = i;
1179                    prob->numberColsToDo_ = 0;
1180                    for (i = 0; i < ncols_; i++)
1181                         if (!prob->colProhibited(i))
1182                              prob->colsToDo_[prob->numberColsToDo_++] = i;
1183               }
1184               // now expensive things
1185               // this caused world.mps to run into numerical difficulties
1186#ifdef PRESOLVE_SUMMARY
1187               printf("Starting expensive\n");
1188#endif
1189
1190               if (dual) {
1191                    int itry;
1192                    for (itry = 0; itry < 5; itry++) {
1193                      possibleBreak;
1194                         paction_ = remove_dual_action::presolve(prob, paction_);
1195                         if (prob->status_)
1196                              break;
1197                         const CoinPresolveAction * const paction2 = paction_;
1198                         if (ifree) {
1199#ifdef IMPLIED
1200#if IMPLIED2 ==0
1201                              int fill_level = 0; // switches off substitution
1202#elif IMPLIED2!=99
1203                              int fill_level = IMPLIED2;
1204#endif
1205#endif
1206                              if ((itry & 1) == 0) {
1207                                possibleBreak;
1208                                   paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1209                              }
1210                              if (prob->status_)
1211                                   break;
1212                         }
1213                         if (paction_ == paction2)
1214                              break;
1215                    }
1216               } else if (ifree) {
1217                    // just free
1218#ifdef IMPLIED
1219#if IMPLIED2 ==0
1220                    int fill_level = 0; // switches off substitution
1221#elif IMPLIED2!=99
1222                    int fill_level = IMPLIED2;
1223#endif
1224#endif
1225                    possibleBreak;
1226                    paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1227                    if (prob->status_)
1228                         break;
1229               }
1230#if     PRESOLVE_DEBUG
1231               check_sol(prob, 1.0e0);
1232#endif
1233               if (dupcol) {
1234                    // maybe allow integer columns to be checked
1235                    if ((presolveActions_ & 512) != 0)
1236                         prob->setPresolveOptions(prob->presolveOptions() | 1);
1237                    possibleBreak;
1238                    paction_ = dupcol_action::presolve(prob, paction_);
1239                    if (prob->status_)
1240                         break;
1241               }
1242#if     PRESOLVE_DEBUG
1243               check_sol(prob, 1.0e0);
1244#endif
1245
1246               if (duprow) {
1247                 possibleBreak;
1248                    paction_ = duprow_action::presolve(prob, paction_);
1249                    if (prob->status_)
1250                         break;
1251               }
1252#if     PRESOLVE_DEBUG
1253               check_sol(prob, 1.0e0);
1254#endif
1255               bool stopLoop = false;
1256               {
1257                    int * hinrow = prob->hinrow_;
1258                    int numberDropped = 0;
1259                    for (i = 0; i < nrows_; i++)
1260                         if (!hinrow[i])
1261                              numberDropped++;
1262
1263                    prob->messageHandler()->message(COIN_PRESOLVE_PASS,
1264                                                    messages)
1265                              << numberDropped << iLoop + 1
1266                              << CoinMessageEol;
1267                    //printf("%d rows dropped after pass %d\n",numberDropped,
1268                    //     iLoop+1);
1269                    if (numberDropped == lastDropped)
1270                         stopLoop = true;
1271                    else
1272                         lastDropped = numberDropped;
1273               }
1274               // Do this here as not very loopy
1275               if (slackSingleton) {
1276                    // On most passes do not touch costed slacks
1277                    if (paction_ != paction0 && !stopLoop) {
1278                      possibleBreak;
1279                         paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
1280                    } else {
1281                         // do costed if Clp (at end as ruins rest of presolve)
1282                      possibleBreak;
1283                         paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
1284                         stopLoop = true;
1285                    }
1286               }
1287#if     PRESOLVE_DEBUG
1288               check_sol(prob, 1.0e0);
1289#endif
1290               if (paction_ == paction0 || stopLoop)
1291                    break;
1292#ifdef ABC_INHERIT
1293               // see whether to stop anyway
1294               int numberRowsNow=nrows_-prob->countEmptyRows();
1295               int numberColumnsNow=ncols_-prob->countEmptyCols();
1296#if ABC_NORMAL_DEBUG
1297               printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n", 
1298                      nrows_,ncols_,numberRowsLeft,numberColumnsLeft,iLoop+1,numberRowsNow,
1299                      numberColumnsNow);
1300#endif
1301               int rowsDeleted=numberRowsLeft-numberRowsNow;
1302               int columnsDeleted=numberColumnsLeft-numberColumnsNow;
1303               if (iLoop>15) {
1304                 if (rowsDeleted*100<numberRowsStart&&
1305                     columnsDeleted*100<numberColumnsStart)
1306                   break;
1307                 lastPassWasGood=true;
1308               } else if (rowsDeleted*100<numberRowsStart&&rowsDeleted<500&&
1309                          columnsDeleted*100<numberColumnsStart&&columnsDeleted<500) {
1310                 if (!lastPassWasGood)
1311                   break;
1312                 else
1313                   lastPassWasGood=false;
1314               } else {
1315                 lastPassWasGood=true;
1316               }
1317               numberRowsLeft=numberRowsNow;
1318               numberColumnsLeft=numberColumnsNow;
1319#endif
1320          }
1321     }
1322     prob->presolveOptions_ &= ~0x10000;
1323     if (!prob->status_) {
1324          paction_ = drop_zero_coefficients(prob, paction_);
1325#if     PRESOLVE_DEBUG
1326          check_sol(prob, 1.0e0);
1327#endif
1328
1329          paction_ = drop_empty_cols_action::presolve(prob, paction_);
1330          paction_ = drop_empty_rows_action::presolve(prob, paction_);
1331#if     PRESOLVE_DEBUG
1332          check_sol(prob, 1.0e0);
1333#endif
1334     }
1335
1336     if (prob->status_) {
1337          if (prob->status_ == 1)
1338               prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
1339                                               messages)
1340                         << prob->feasibilityTolerance_
1341                         << CoinMessageEol;
1342          else if (prob->status_ == 2)
1343               prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
1344                                               messages)
1345                         << CoinMessageEol;
1346          else
1347               prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
1348                                               messages)
1349                         << CoinMessageEol;
1350          // get rid of data
1351          destroyPresolve();
1352     }
1353     return (paction_);
1354}
1355
1356void check_djs(CoinPostsolveMatrix *prob);
1357
1358
1359// We could have implemented this by having each postsolve routine
1360// directly call the next one, but this may make it easier to add debugging checks.
1361void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
1362{
1363     {
1364          // Check activities
1365          double *colels        = prob.colels_;
1366          int *hrow             = prob.hrow_;
1367          CoinBigIndex *mcstrt          = prob.mcstrt_;
1368          int *hincol           = prob.hincol_;
1369          int *link             = prob.link_;
1370          int ncols             = prob.ncols_;
1371
1372          char *cdone   = prob.cdone_;
1373
1374          double * csol = prob.sol_;
1375          int nrows = prob.nrows_;
1376
1377          int colx;
1378
1379          double * rsol = prob.acts_;
1380          memset(rsol, 0, nrows * sizeof(double));
1381
1382          for (colx = 0; colx < ncols; ++colx) {
1383               if (cdone[colx]) {
1384                    CoinBigIndex k = mcstrt[colx];
1385                    int nx = hincol[colx];
1386                    double solutionValue = csol[colx];
1387                    for (int i = 0; i < nx; ++i) {
1388                         int row = hrow[k];
1389                         double coeff = colels[k];
1390                         k = link[k];
1391                         rsol[row] += solutionValue * coeff;
1392                    }
1393               }
1394          }
1395     }
1396     if (prob.maxmin_<0) {
1397       //for (int i=0;i<presolvedModel_->numberRows();i++)
1398       //prob.rowduals_[i]=-prob.rowduals_[i];
1399       for (int i=0;i<ncols_;i++) {
1400         prob.cost_[i]=-prob.cost_[i];
1401         //prob.rcosts_[i]=-prob.rcosts_[i];
1402       }
1403       prob.maxmin_=1.0;
1404     }
1405     const CoinPresolveAction *paction = paction_;
1406     //#define PRESOLVE_DEBUG 1
1407#if     PRESOLVE_DEBUG
1408     // Check only works after first one
1409     int checkit = -1;
1410#endif
1411
1412     while (paction) {
1413#if PRESOLVE_DEBUG
1414          printf("POSTSOLVING %s\n", paction->name());
1415#endif
1416
1417          paction->postsolve(&prob);
1418
1419#if     PRESOLVE_DEBUG
1420#         if 0
1421          /*
1422            This check fails (on exmip1 (!) in osiUnitTest) because clp
1423            enters postsolve with a solution that seems to have incorrect
1424            status for a logical. You can see similar behaviour with
1425            column status --- incorrect entering postsolve.
1426            -- lh, 111207 --
1427          */
1428          {
1429               int nr = 0;
1430               int i;
1431               for (i = 0; i < prob.nrows_; i++) {
1432                    if ((prob.rowstat_[i] & 7) == 1) {
1433                         nr++;
1434                    } else if ((prob.rowstat_[i] & 7) == 2) {
1435                         // at ub
1436                         assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
1437                    } else if ((prob.rowstat_[i] & 7) == 3) {
1438                         // at lb
1439                         assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
1440                    }
1441               }
1442               int nc = 0;
1443               for (i = 0; i < prob.ncols_; i++) {
1444                    if ((prob.colstat_[i] & 7) == 1)
1445                         nc++;
1446               }
1447               printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
1448          }
1449#         endif   // if 0
1450          checkit++;
1451          if (prob.colstat_ && checkit > 0) {
1452               presolve_check_nbasic(&prob) ;
1453               presolve_check_sol(&prob, 2, 2, 1) ;
1454          }
1455#endif
1456          paction = paction->next;
1457#if     PRESOLVE_DEBUG
1458//  check_djs(&prob);
1459          if (checkit > 0)
1460               presolve_check_reduced_costs(&prob) ;
1461#endif
1462     }
1463#if     PRESOLVE_DEBUG
1464     if (prob.colstat_) {
1465          presolve_check_nbasic(&prob) ;
1466          presolve_check_sol(&prob, 2, 2, 1) ;
1467     }
1468#endif
1469#undef PRESOLVE_DEBUG
1470
1471#if     0 && PRESOLVE_DEBUG
1472     for (i = 0; i < ncols0; i++) {
1473          if (!cdone[i]) {
1474               printf("!cdone[%d]\n", i);
1475               abort();
1476          }
1477     }
1478
1479     for (i = 0; i < nrows0; i++) {
1480          if (!rdone[i]) {
1481               printf("!rdone[%d]\n", i);
1482               abort();
1483          }
1484     }
1485
1486
1487     for (i = 0; i < ncols0; i++) {
1488          if (sol[i] < -1e10 || sol[i] > 1e10)
1489               printf("!!!%d %g\n", i, sol[i]);
1490
1491     }
1492
1493
1494#endif
1495
1496#if     0 && PRESOLVE_DEBUG
1497     // debug check:  make sure we ended up with same original matrix
1498     {
1499          int identical = 1;
1500
1501          for (int i = 0; i < ncols0; i++) {
1502               PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
1503               CoinBigIndex kcs0 = &prob->mcstrt0[i];
1504               CoinBigIndex kcs = mcstrt[i];
1505               int n = hincol[i];
1506               for (int k = 0; k < n; k++) {
1507                    CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);
1508
1509                    if (k1 == kcs + n) {
1510                         printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
1511                         abort();
1512                    }
1513
1514                    if (colels[k1] != &prob->dels0[kcs0+k])
1515                         printf("BAD COLEL[%d %d %d]:  %g\n",
1516                                k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
1517
1518                    if (kcs0 + k != k1)
1519                         identical = 0;
1520               }
1521          }
1522          printf("identical? %d\n", identical);
1523     }
1524#endif
1525}
1526
1527
1528
1529
1530
1531
1532
1533
1534static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
1535{
1536     double tol;
1537     if (! si->getDblParam(key, tol)) {
1538          CoinPresolveAction::throwCoinError("getDblParam failed",
1539                                             "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
1540     }
1541     return (tol);
1542}
1543
1544
1545// Assumptions:
1546// 1. nrows>=m.getNumRows()
1547// 2. ncols>=m.getNumCols()
1548//
1549// In presolve, these values are equal.
1550// In postsolve, they may be inequal, since the reduced problem
1551// may be smaller, but we need room for the large problem.
1552// ncols may be larger than si.getNumCols() in postsolve,
1553// this at that point si will be the reduced problem,
1554// but we need to reserve enough space for the original problem.
1555CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
1556          int ncols_in,
1557          int nrows_in,
1558          CoinBigIndex nelems_in,
1559          double bulkRatio)
1560     : ncols_(si->getNumCols()),
1561       nrows_(si->getNumRows()),
1562       nelems_(si->getNumElements()),
1563       ncols0_(ncols_in),
1564       nrows0_(nrows_in),
1565       bulkRatio_(bulkRatio),
1566       mcstrt_(new CoinBigIndex[ncols_in+1]),
1567       hincol_(new int[ncols_in+1]),
1568       cost_(new double[ncols_in]),
1569       clo_(new double[ncols_in]),
1570       cup_(new double[ncols_in]),
1571       rlo_(new double[nrows_in]),
1572       rup_(new double[nrows_in]),
1573       originalColumn_(new int[ncols_in]),
1574       originalRow_(new int[nrows_in]),
1575       ztolzb_(getTolerance(si, ClpPrimalTolerance)),
1576       ztoldj_(getTolerance(si, ClpDualTolerance)),
1577       maxmin_(si->getObjSense()),
1578       sol_(NULL),
1579       rowduals_(NULL),
1580       acts_(NULL),
1581       rcosts_(NULL),
1582       colstat_(NULL),
1583       rowstat_(NULL),
1584       handler_(NULL),
1585       defaultHandler_(false),
1586       messages_()
1587
1588{
1589     bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ * nelems_in);
1590     hrow_  = new int   [bulk0_];
1591     colels_ = new double[bulk0_];
1592     si->getDblParam(ClpObjOffset, originalOffset_);
1593     int ncols = si->getNumCols();
1594     int nrows = si->getNumRows();
1595
1596     setMessageHandler(si->messageHandler()) ;
1597
1598     ClpDisjointCopyN(si->getColLower(), ncols, clo_);
1599     ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
1600     //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
1601     double offset;
1602     ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cost_);
1603     ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
1604     ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
1605     int i;
1606     for (i = 0; i < ncols_in; i++)
1607          originalColumn_[i] = i;
1608     for (i = 0; i < nrows_in; i++)
1609          originalRow_[i] = i;
1610     sol_ = NULL;
1611     rowduals_ = NULL;
1612     acts_ = NULL;
1613
1614     rcosts_ = NULL;
1615     colstat_ = NULL;
1616     rowstat_ = NULL;
1617}
1618
1619// I am not familiar enough with CoinPackedMatrix to be confident
1620// that I will implement a row-ordered version of toColumnOrderedGapFree
1621// properly.
1622static bool isGapFree(const CoinPackedMatrix& matrix)
1623{
1624     const CoinBigIndex * start = matrix.getVectorStarts();
1625     const int * length = matrix.getVectorLengths();
1626     int i = matrix.getSizeVectorLengths() - 1;
1627     // Quick check
1628     if (matrix.getNumElements() == start[i]) {
1629          return true;
1630     } else {
1631          for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
1632               if (start[i+1] - start[i] != length[i])
1633                    break;
1634          }
1635          return (! (i >= 0));
1636     }
1637}
1638#if     PRESOLVE_DEBUG
1639static void matrix_bounds_ok(const double *lo, const double *up, int n)
1640{
1641     int i;
1642     for (i = 0; i < n; i++) {
1643          PRESOLVEASSERT(lo[i] <= up[i]);
1644          PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
1645          PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
1646     }
1647}
1648#endif
1649CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
1650                                       double /*maxmin*/,
1651                                       // end prepost members
1652
1653                                       ClpSimplex * si,
1654
1655                                       // rowrep
1656                                       int nrows_in,
1657                                       CoinBigIndex nelems_in,
1658                                       bool doStatus,
1659                                       double nonLinearValue,
1660                                       double bulkRatio) :
1661
1662     CoinPrePostsolveMatrix(si,
1663                            ncols0_in, nrows_in, nelems_in, bulkRatio),
1664     clink_(new presolvehlink[ncols0_in+1]),
1665     rlink_(new presolvehlink[nrows_in+1]),
1666
1667     dobias_(0.0),
1668
1669
1670     // temporary init
1671     integerType_(new unsigned char[ncols0_in]),
1672     tuning_(false),
1673     startTime_(0.0),
1674     feasibilityTolerance_(0.0),
1675     status_(-1),
1676     colsToDo_(new int [ncols0_in]),
1677     numberColsToDo_(0),
1678     nextColsToDo_(new int[ncols0_in]),
1679     numberNextColsToDo_(0),
1680     rowsToDo_(new int [nrows_in]),
1681     numberRowsToDo_(0),
1682     nextRowsToDo_(new int[nrows_in]),
1683     numberNextRowsToDo_(0),
1684     presolveOptions_(0)
1685{
1686     const int bufsize = bulk0_;
1687
1688     nrows_ = si->getNumRows() ;
1689
1690     // Set up change bits etc
1691     rowChanged_ = new unsigned char[nrows_];
1692     memset(rowChanged_, 0, nrows_);
1693     colChanged_ = new unsigned char[ncols_];
1694     memset(colChanged_, 0, ncols_);
1695     CoinPackedMatrix * m = si->matrix();
1696
1697     // The coefficient matrix is a big hunk of stuff.
1698     // Do the copy here to try to avoid running out of memory.
1699
1700     const CoinBigIndex * start = m->getVectorStarts();
1701     const int * row = m->getIndices();
1702     const double * element = m->getElements();
1703     int icol, nel = 0;
1704     mcstrt_[0] = 0;
1705     ClpDisjointCopyN(m->getVectorLengths(), ncols_,  hincol_);
1706     if (si->getObjSense() < 0.0) {
1707       for (int i=0;i<ncols_;i++)
1708         cost_[i]=-cost_[i];
1709       maxmin_=1.0;
1710     }
1711     for (icol = 0; icol < ncols_; icol++) {
1712          CoinBigIndex j;
1713          for (j = start[icol]; j < start[icol] + hincol_[icol]; j++) {
1714               hrow_[nel] = row[j];
1715               if (fabs(element[j]) > ZTOLDP)
1716                    colels_[nel++] = element[j];
1717          }
1718          mcstrt_[icol+1] = nel;
1719          hincol_[icol] = nel - mcstrt_[icol];
1720     }
1721
1722     // same thing for row rep
1723     CoinPackedMatrix * mRow = new CoinPackedMatrix();
1724     mRow->setExtraGap(0.0);
1725     mRow->setExtraMajor(0.0);
1726     mRow->reverseOrderedCopyOf(*m);
1727     //mRow->removeGaps();
1728     //mRow->setExtraGap(0.0);
1729
1730     // Now get rid of matrix
1731     si->createEmptyMatrix();
1732
1733     double * el = mRow->getMutableElements();
1734     int * ind = mRow->getMutableIndices();
1735     CoinBigIndex * strt = mRow->getMutableVectorStarts();
1736     int * len = mRow->getMutableVectorLengths();
1737     // Do carefully to save memory
1738     rowels_ = new double[bulk0_];
1739     ClpDisjointCopyN(el,      nelems_, rowels_);
1740     mRow->nullElementArray();
1741     delete [] el;
1742     hcol_ = new int[bulk0_];
1743     ClpDisjointCopyN(ind,       nelems_, hcol_);
1744     mRow->nullIndexArray();
1745     delete [] ind;
1746     mrstrt_ = new CoinBigIndex[nrows_in+1];
1747     ClpDisjointCopyN(strt,  nrows_,  mrstrt_);
1748     mRow->nullStartArray();
1749     mrstrt_[nrows_] = nelems_;
1750     delete [] strt;
1751     hinrow_ = new int[nrows_in+1];
1752     ClpDisjointCopyN(len, nrows_,  hinrow_);
1753     if (nelems_ > nel) {
1754          nelems_ = nel;
1755          // Clean any small elements
1756          int irow;
1757          nel = 0;
1758          CoinBigIndex start = 0;
1759          for (irow = 0; irow < nrows_; irow++) {
1760               CoinBigIndex j;
1761               for (j = start; j < start + hinrow_[irow]; j++) {
1762                    hcol_[nel] = hcol_[j];
1763                    if (fabs(rowels_[j]) > ZTOLDP)
1764                         rowels_[nel++] = rowels_[j];
1765               }
1766               start = mrstrt_[irow+1];
1767               mrstrt_[irow+1] = nel;
1768               hinrow_[irow] = nel - mrstrt_[irow];
1769          }
1770     }
1771
1772     delete mRow;
1773     if (si->integerInformation()) {
1774          CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()), ncols_, integerType_);
1775     } else {
1776          ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
1777     }
1778
1779#ifndef SLIM_CLP
1780#ifndef NO_RTTI
1781     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1782#else
1783     ClpQuadraticObjective * quadraticObj = NULL;
1784     if (si->objectiveAsObject()->type() == 2)
1785          quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1786#endif
1787#endif
1788     // Set up prohibited bits if needed
1789     if (nonLinearValue) {
1790          anyProhibited_ = true;
1791          for (icol = 0; icol < ncols_; icol++) {
1792               int j;
1793               bool nonLinearColumn = false;
1794               if (cost_[icol] == nonLinearValue)
1795                    nonLinearColumn = true;
1796               for (j = mcstrt_[icol]; j < mcstrt_[icol+1]; j++) {
1797                    if (colels_[j] == nonLinearValue) {
1798                         nonLinearColumn = true;
1799                         setRowProhibited(hrow_[j]);
1800                    }
1801               }
1802               if (nonLinearColumn)
1803                    setColProhibited(icol);
1804          }
1805#ifndef SLIM_CLP
1806     } else if (quadraticObj) {
1807          CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1808          //const int * columnQuadratic = quadratic->getIndices();
1809          //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1810          const int * columnQuadraticLength = quadratic->getVectorLengths();
1811          //double * quadraticElement = quadratic->getMutableElements();
1812          int numberColumns = quadratic->getNumCols();
1813          anyProhibited_ = true;
1814          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1815               if (columnQuadraticLength[iColumn]) {
1816                    setColProhibited(iColumn);
1817                    //printf("%d prohib\n",iColumn);
1818               }
1819          }
1820#endif
1821     } else {
1822          anyProhibited_ = false;
1823     }
1824
1825     if (doStatus) {
1826          // allow for status and solution
1827          sol_ = new double[ncols_];
1828          CoinMemcpyN(si->primalColumnSolution(), ncols_, sol_);;
1829          acts_ = new double [nrows_];
1830          CoinMemcpyN(si->primalRowSolution(), nrows_, acts_);
1831          if (!si->statusArray())
1832               si->createStatus();
1833          colstat_ = new unsigned char [nrows_+ncols_];
1834          CoinMemcpyN(si->statusArray(),        (nrows_ + ncols_), colstat_);
1835          rowstat_ = colstat_ + ncols_;
1836     }
1837
1838     // the original model's fields are now unneeded - free them
1839
1840     si->resize(0, 0);
1841
1842#if     PRESOLVE_DEBUG
1843     matrix_bounds_ok(rlo_, rup_, nrows_);
1844     matrix_bounds_ok(clo_, cup_, ncols_);
1845#endif
1846
1847#if 0
1848     for (i = 0; i < nrows; ++i)
1849          printf("NR: %6d\n", hinrow[i]);
1850     for (int i = 0; i < ncols; ++i)
1851          printf("NC: %6d\n", hincol[i]);
1852#endif
1853
1854     presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
1855     presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
1856
1857     // this allows last col/row to expand up to bufsize-1 (22);
1858     // this must come after the calls to presolve_prefix
1859     mcstrt_[ncols_] = bufsize - 1;
1860     mrstrt_[nrows_] = bufsize - 1;
1861     // Allocate useful arrays
1862     initializeStuff();
1863
1864#if     PRESOLVE_CONSISTENCY
1865//consistent(false);
1866     presolve_consistent(this, false) ;
1867#endif
1868}
1869
1870// avoid compiler warnings
1871#if PRESOLVE_SUMMARY > 0
1872void CoinPresolveMatrix::update_model(ClpSimplex * si,
1873                                      int nrows0, int ncols0,
1874                                      CoinBigIndex nelems0)
1875#else
1876void CoinPresolveMatrix::update_model(ClpSimplex * si,
1877                                      int /*nrows0*/,
1878                                      int /*ncols0*/,
1879                                      CoinBigIndex /*nelems0*/)
1880#endif
1881{
1882     if (si->getObjSense() < 0.0) {
1883       for (int i=0;i<ncols_;i++)
1884         cost_[i]=-cost_[i];
1885       dobias_=-dobias_;
1886     }
1887     si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1888                     clo_, cup_, cost_, rlo_, rup_);
1889     //delete [] si->integerInformation();
1890     int numberIntegers = 0;
1891     for (int i = 0; i < ncols_; i++) {
1892          if (integerType_[i])
1893               numberIntegers++;
1894     }
1895     if (numberIntegers)
1896          si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
1897     else
1898          si->copyInIntegerInformation(NULL);
1899
1900#if     PRESOLVE_SUMMARY
1901     printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
1902            ncols_, ncols0 - ncols_,
1903            nrows_, nrows0 - nrows_,
1904            si->getNumElements(), nelems0 - si->getNumElements());
1905#endif
1906     si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
1907     if (si->getObjSense() < 0.0) {
1908       // put back
1909       for (int i=0;i<ncols_;i++)
1910         cost_[i]=-cost_[i];
1911       dobias_=-dobias_;
1912       maxmin_=-1.0;
1913     }
1914
1915}
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927////////////////  POSTSOLVE
1928
1929CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
1930          int ncols0_in,
1931          int nrows0_in,
1932          CoinBigIndex nelems0,
1933
1934          double maxmin,
1935          // end prepost members
1936
1937          double *sol_in,
1938          double *acts_in,
1939
1940          unsigned char *colstat_in,
1941          unsigned char *rowstat_in) :
1942     CoinPrePostsolveMatrix(si,
1943                            ncols0_in, nrows0_in, nelems0, 2.0),
1944
1945     free_list_(0),
1946     // link, free_list, maxlink
1947     maxlink_(bulk0_),
1948     link_(new int[/*maxlink*/ bulk0_]),
1949
1950     cdone_(new char[ncols0_]),
1951     rdone_(new char[nrows0_in])
1952
1953{
1954     bulk0_ = maxlink_ ;
1955     nrows_ = si->getNumRows() ;
1956     ncols_ = si->getNumCols() ;
1957
1958     sol_ = sol_in;
1959     rowduals_ = NULL;
1960     acts_ = acts_in;
1961
1962     rcosts_ = NULL;
1963     colstat_ = colstat_in;
1964     rowstat_ = rowstat_in;
1965
1966     // this is the *reduced* model, which is probably smaller
1967     int ncols1 = ncols_ ;
1968     int nrows1 = nrows_ ;
1969
1970     const CoinPackedMatrix * m = si->matrix();
1971
1972     const CoinBigIndex nelemsr = m->getNumElements();
1973     if (m->getNumElements() && !isGapFree(*m)) {
1974          // Odd - gaps
1975          CoinPackedMatrix mm(*m);
1976          mm.removeGaps();
1977          mm.setExtraGap(0.0);
1978
1979          ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
1980          CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
1981          mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
1982          ClpDisjointCopyN(mm.getVectorLengths(), ncols1,  hincol_);
1983          ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
1984          ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
1985     } else {
1986          // No gaps
1987
1988          ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1989          CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
1990          mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
1991          ClpDisjointCopyN(m->getVectorLengths(), ncols1,  hincol_);
1992          ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
1993          ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
1994     }
1995
1996
1997
1998#if     0 && PRESOLVE_DEBUG
1999     presolve_check_costs(model, &colcopy);
2000#endif
2001
2002     // This determines the size of the data structure that contains
2003     // the matrix being postsolved.  Links are taken from the free_list
2004     // to recreate matrix entries that were presolved away,
2005     // and links are added to the free_list when entries created during
2006     // presolve are discarded.  There is never a need to gc this list.
2007     // Naturally, it should contain
2008     // exactly nelems0 entries "in use" when postsolving is done,
2009     // but I don't know whether the matrix could temporarily get
2010     // larger during postsolving.  Substitution into more than two
2011     // rows could do that, in principle.  I am being very conservative
2012     // here by reserving much more than the amount of space I probably need.
2013     // If this guess is wrong, check_free_list may be called.
2014     //  int bufsize = 2*nelems0;
2015
2016     memset(cdone_, -1, ncols0_);
2017     memset(rdone_, -1, nrows0_);
2018
2019     rowduals_ = new double[nrows0_];
2020     ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
2021
2022     rcosts_ = new double[ncols0_];
2023     ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
2024     if (maxmin < 0.0) {
2025          // change so will look as if minimize
2026          int i;
2027          for (i = 0; i < nrows1; i++)
2028               rowduals_[i] = - rowduals_[i];
2029          for (i = 0; i < ncols1; i++) {
2030               rcosts_[i] = - rcosts_[i];
2031          }
2032     }
2033
2034     //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
2035     //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
2036
2037     ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
2038     si->setDblParam(ClpObjOffset, originalOffset_);
2039
2040     for (int j = 0; j < ncols1; j++) {
2041#ifdef COIN_SLOW_PRESOLVE
2042       if (hincol_[j]) {
2043#endif
2044          CoinBigIndex kcs = mcstrt_[j];
2045          CoinBigIndex kce = kcs + hincol_[j];
2046          for (CoinBigIndex k = kcs; k < kce; ++k) {
2047               link_[k] = k + 1;
2048          }
2049          link_[kce-1] = NO_LINK ;
2050#ifdef COIN_SLOW_PRESOLVE
2051       }
2052#endif
2053     }
2054     {
2055          int ml = maxlink_;
2056          for (CoinBigIndex k = nelemsr; k < ml; ++k)
2057               link_[k] = k + 1;
2058          if (ml)
2059               link_[ml-1] = NO_LINK;
2060     }
2061     free_list_ = nelemsr;
2062# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
2063     /*
2064       These are used to track the action of postsolve transforms during debugging.
2065     */
2066     CoinFillN(cdone_, ncols1, PRESENT_IN_REDUCED) ;
2067     CoinZeroN(cdone_ + ncols1, ncols0_in - ncols1) ;
2068     CoinFillN(rdone_, nrows1, PRESENT_IN_REDUCED) ;
2069     CoinZeroN(rdone_ + nrows1, nrows0_in - nrows1) ;
2070# endif
2071}
2072/* This is main part of Presolve */
2073ClpSimplex *
2074ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
2075                                  double feasibilityTolerance,
2076                                  bool keepIntegers,
2077                                  int numberPasses,
2078                                  bool dropNames,
2079                                  bool doRowObjective,
2080                                  const char * prohibitedRows,
2081                                  const char * prohibitedColumns)
2082{
2083     ncols_ = originalModel->getNumCols();
2084     nrows_ = originalModel->getNumRows();
2085     nelems_ = originalModel->getNumElements();
2086     numberPasses_ = numberPasses;
2087
2088     double maxmin = originalModel->getObjSense();
2089     originalModel_ = originalModel;
2090     delete [] originalColumn_;
2091     originalColumn_ = new int[ncols_];
2092     delete [] originalRow_;
2093     originalRow_ = new int[nrows_];
2094     // and fill in case returns early
2095     int i;
2096     for (i = 0; i < ncols_; i++)
2097          originalColumn_[i] = i;
2098     for (i = 0; i < nrows_; i++)
2099          originalRow_[i] = i;
2100     delete [] rowObjective_;
2101     if (doRowObjective) {
2102          rowObjective_ = new double [nrows_];
2103          memset(rowObjective_, 0, nrows_ * sizeof(double));
2104     } else {
2105          rowObjective_ = NULL;
2106     }
2107
2108     // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
2109     int result = -1;
2110
2111     // User may have deleted - its their responsibility
2112     presolvedModel_ = NULL;
2113     // Messages
2114     CoinMessages messages = originalModel->coinMessages();
2115     // Only go round 100 times even if integer preprocessing
2116     int totalPasses = 100;
2117     while (result == -1) {
2118
2119#ifndef CLP_NO_STD
2120          // make new copy
2121          if (saveFile_ == "") {
2122#endif
2123               delete presolvedModel_;
2124#ifndef CLP_NO_STD
2125               // So won't get names
2126               int lengthNames = originalModel->lengthNames();
2127               originalModel->setLengthNames(0);
2128#endif
2129               presolvedModel_ = new ClpSimplex(*originalModel);
2130#ifndef CLP_NO_STD
2131               originalModel->setLengthNames(lengthNames);
2132               presolvedModel_->dropNames();
2133          } else {
2134               presolvedModel_ = originalModel;
2135               if (dropNames)
2136                 presolvedModel_->dropNames();
2137          }
2138#endif
2139
2140          // drop integer information if wanted
2141          if (!keepIntegers)
2142               presolvedModel_->deleteIntegerInformation();
2143          totalPasses--;
2144
2145          double ratio = 2.0;
2146          if (substitution_ > 3)
2147               ratio = substitution_;
2148          else if (substitution_ == 2)
2149               ratio = 1.5;
2150          CoinPresolveMatrix prob(ncols_,
2151                                  maxmin,
2152                                  presolvedModel_,
2153                                  nrows_, nelems_, true, nonLinearValue_, ratio);
2154          if (prohibitedRows) {
2155            prob.setAnyProhibited();
2156            for (int i=0;i<nrows_;i++) {
2157              if (prohibitedRows[i])
2158                prob.setRowProhibited(i);
2159            }
2160          }
2161          if (prohibitedColumns) {
2162            prob.setAnyProhibited();
2163            for (int i=0;i<ncols_;i++) {
2164              if (prohibitedColumns[i])
2165                prob.setColProhibited(i);
2166            }
2167          }
2168          prob.setMaximumSubstitutionLevel(substitution_);
2169          if (doRowObjective)
2170               memset(rowObjective_, 0, nrows_ * sizeof(double));
2171          // See if we want statistics
2172          if ((presolveActions_ & 0x80000000) != 0)
2173               prob.statistics();
2174          // make sure row solution correct
2175          {
2176               double *colels   = prob.colels_;
2177               int *hrow                = prob.hrow_;
2178               CoinBigIndex *mcstrt             = prob.mcstrt_;
2179               int *hincol              = prob.hincol_;
2180               int ncols                = prob.ncols_;
2181
2182
2183               double * csol = prob.sol_;
2184               double * acts = prob.acts_;
2185               int nrows = prob.nrows_;
2186
2187               int colx;
2188
2189               memset(acts, 0, nrows * sizeof(double));
2190
2191               for (colx = 0; colx < ncols; ++colx) {
2192                    double solutionValue = csol[colx];
2193                    for (int i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
2194                         int row = hrow[i];
2195                         double coeff = colels[i];
2196                         acts[row] += solutionValue * coeff;
2197                    }
2198               }
2199          }
2200
2201          // move across feasibility tolerance
2202          prob.feasibilityTolerance_ = feasibilityTolerance;
2203
2204          // Do presolve
2205          paction_ = presolve(&prob);
2206          // Get rid of useful arrays
2207          prob.deleteStuff();
2208
2209          result = 0;
2210
2211          bool fixInfeasibility = (prob.presolveOptions_&16384)!=0;
2212          bool hasSolution = (prob.presolveOptions_&32768)!=0;
2213          if (prob.status_ == 0 && paction_ && (!hasSolution || !fixInfeasibility)) {
2214               // Looks feasible but double check to see if anything slipped through
2215               int n            = prob.ncols_;
2216               double * lo = prob.clo_;
2217               double * up = prob.cup_;
2218               int i;
2219
2220               for (i = 0; i < n; i++) {
2221                    if (up[i] < lo[i]) {
2222                         if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2223                              // infeasible
2224                              prob.status_ = 1;
2225                         } else {
2226                              up[i] = lo[i];
2227                         }
2228                    }
2229               }
2230
2231               n = prob.nrows_;
2232               lo = prob.rlo_;
2233               up = prob.rup_;
2234
2235               for (i = 0; i < n; i++) {
2236                    if (up[i] < lo[i]) {
2237                         if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2238                              // infeasible
2239                              prob.status_ = 1;
2240                         } else {
2241                              up[i] = lo[i];
2242                         }
2243                    }
2244               }
2245          }
2246          if (prob.status_ == 0 && paction_) {
2247               // feasible
2248
2249               prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
2250               // copy status and solution
2251               CoinMemcpyN(          prob.sol_, prob.ncols_, presolvedModel_->primalColumnSolution());
2252               CoinMemcpyN(          prob.acts_, prob.nrows_, presolvedModel_->primalRowSolution());
2253               CoinMemcpyN(          prob.colstat_, prob.ncols_, presolvedModel_->statusArray());
2254               CoinMemcpyN(          prob.rowstat_, prob.nrows_, presolvedModel_->statusArray() + prob.ncols_);
2255               if (fixInfeasibility && hasSolution) {
2256                 // Looks feasible but double check to see if anything slipped through
2257                 int n          = prob.ncols_;
2258                 double * lo = prob.clo_;
2259                 double * up = prob.cup_;
2260                 double * rsol = prob.acts_;
2261                 //memset(prob.acts_,0,prob.nrows_*sizeof(double));
2262                 presolvedModel_->matrix()->times(prob.sol_,rsol);
2263                 int i;
2264                 
2265                 for (i = 0; i < n; i++) {
2266                   double gap=up[i]-lo[i];
2267                   if (rsol[i]<lo[i]-feasibilityTolerance&&fabs(rsol[i]-lo[i])<1.0e-3) {
2268                     lo[i]=rsol[i];
2269                     if (gap<1.0e5)
2270                       up[i]=lo[i]+gap;
2271                   } else if (rsol[i]>up[i]+feasibilityTolerance&&fabs(rsol[i]-up[i])<1.0e-3) {
2272                     up[i]=rsol[i];
2273                     if (gap<1.0e5)
2274                       lo[i]=up[i]-gap;
2275                   }
2276                   if (up[i] < lo[i]) {
2277                     up[i] = lo[i];
2278                   }
2279                 }
2280               }
2281
2282               int n = prob.nrows_;
2283               double * lo = prob.rlo_;
2284               double * up = prob.rup_;
2285
2286               for (i = 0; i < n; i++) {
2287                    if (up[i] < lo[i]) {
2288                         if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2289                              // infeasible
2290                              prob.status_ = 1;
2291                         } else {
2292                              up[i] = lo[i];
2293                         }
2294                    }
2295               }
2296               delete [] prob.sol_;
2297               delete [] prob.acts_;
2298               delete [] prob.colstat_;
2299               prob.sol_ = NULL;
2300               prob.acts_ = NULL;
2301               prob.colstat_ = NULL;
2302
2303               int ncolsNow = presolvedModel_->getNumCols();
2304               CoinMemcpyN(prob.originalColumn_, ncolsNow, originalColumn_);
2305#ifndef SLIM_CLP
2306#ifndef NO_RTTI
2307               ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
2308#else
2309               ClpQuadraticObjective * quadraticObj = NULL;
2310               if (originalModel->objectiveAsObject()->type() == 2)
2311                    quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
2312#endif
2313               if (quadraticObj) {
2314                    // set up for subset
2315                    char * mark = new char [ncols_];
2316                    memset(mark, 0, ncols_);
2317                    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2318                    //const int * columnQuadratic = quadratic->getIndices();
2319                    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
2320                    const int * columnQuadraticLength = quadratic->getVectorLengths();
2321                    //double * quadraticElement = quadratic->getMutableElements();
2322                    int numberColumns = quadratic->getNumCols();
2323                    ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
2324                              ncolsNow,
2325                              originalColumn_);
2326                    // and modify linear and check
2327                    double * linear = newObj->linearObjective();
2328                    CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
2329                    int iColumn;
2330                    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
2331                         if (columnQuadraticLength[iColumn])
2332                              mark[iColumn] = 1;
2333                    }
2334                    // and new
2335                    quadratic = newObj->quadraticObjective();
2336                    columnQuadraticLength = quadratic->getVectorLengths();
2337                    int numberColumns2 = quadratic->getNumCols();
2338                    for ( iColumn = 0; iColumn < numberColumns2; iColumn++) {
2339                         if (columnQuadraticLength[iColumn])
2340                              mark[originalColumn_[iColumn]] = 0;
2341                    }
2342                    presolvedModel_->setObjective(newObj);
2343                    delete newObj;
2344                    // final check
2345                    for ( iColumn = 0; iColumn < numberColumns; iColumn++)
2346                         if (mark[iColumn])
2347                              printf("Quadratic column %d modified - may be okay\n", iColumn);
2348                    delete [] mark;
2349               }
2350#endif
2351               delete [] prob.originalColumn_;
2352               prob.originalColumn_ = NULL;
2353               int nrowsNow = presolvedModel_->getNumRows();
2354               CoinMemcpyN(prob.originalRow_, nrowsNow, originalRow_);
2355               delete [] prob.originalRow_;
2356               prob.originalRow_ = NULL;
2357#ifndef CLP_NO_STD
2358               if (!dropNames && originalModel->lengthNames()) {
2359                    // Redo names
2360                    int iRow;
2361                    std::vector<std::string> rowNames;
2362                    rowNames.reserve(nrowsNow);
2363                    for (iRow = 0; iRow < nrowsNow; iRow++) {
2364                         int kRow = originalRow_[iRow];
2365                         rowNames.push_back(originalModel->rowName(kRow));
2366                    }
2367
2368                    int iColumn;
2369                    std::vector<std::string> columnNames;
2370                    columnNames.reserve(ncolsNow);
2371                    for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
2372                         int kColumn = originalColumn_[iColumn];
2373                         columnNames.push_back(originalModel->columnName(kColumn));
2374                    }
2375                    presolvedModel_->copyNames(rowNames, columnNames);
2376               } else {
2377                    presolvedModel_->setLengthNames(0);
2378               }
2379#endif
2380               if (rowObjective_) {
2381                    int iRow;
2382#ifndef NDEBUG
2383                    int k = -1;
2384#endif
2385                    int nObj = 0;
2386                    for (iRow = 0; iRow < nrowsNow; iRow++) {
2387                         int kRow = originalRow_[iRow];
2388#ifndef NDEBUG
2389                         assert (kRow > k);
2390                         k = kRow;
2391#endif
2392                         rowObjective_[iRow] = rowObjective_[kRow];
2393                         if (rowObjective_[iRow])
2394                              nObj++;
2395                    }
2396                    if (nObj) {
2397                         printf("%d costed slacks\n", nObj);
2398                         presolvedModel_->setRowObjective(rowObjective_);
2399                    }
2400               }
2401               /* now clean up integer variables.  This can modify original
2402                          Don't do if dupcol added columns together */
2403               int i;
2404               const char * information = presolvedModel_->integerInformation();
2405               if ((prob.presolveOptions_ & 0x80000000) == 0 && information) {
2406                    int numberChanges = 0;
2407                    double * lower0 = originalModel_->columnLower();
2408                    double * upper0 = originalModel_->columnUpper();
2409                    double * lower = presolvedModel_->columnLower();
2410                    double * upper = presolvedModel_->columnUpper();
2411                    for (i = 0; i < ncolsNow; i++) {
2412                         if (!information[i])
2413                              continue;
2414                         int iOriginal = originalColumn_[i];
2415                         double lowerValue0 = lower0[iOriginal];
2416                         double upperValue0 = upper0[iOriginal];
2417                         double lowerValue = ceil(lower[i] - 1.0e-5);
2418                         double upperValue = floor(upper[i] + 1.0e-5);
2419                         lower[i] = lowerValue;
2420                         upper[i] = upperValue;
2421                         if (lowerValue > upperValue) {
2422                              numberChanges++;
2423                              presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
2424                                        messages)
2425                                        << iOriginal
2426                                        << lowerValue
2427                                        << upperValue
2428                                        << CoinMessageEol;
2429                              result = 1;
2430                         } else {
2431                              if (lowerValue > lowerValue0 + 1.0e-8) {
2432                                   lower0[iOriginal] = lowerValue;
2433                                   numberChanges++;
2434                              }
2435                              if (upperValue < upperValue0 - 1.0e-8) {
2436                                   upper0[iOriginal] = upperValue;
2437                                   numberChanges++;
2438                              }
2439                         }
2440                    }
2441                    if (numberChanges) {
2442                         presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
2443                                   messages)
2444                                   << numberChanges
2445                                   << CoinMessageEol;
2446                         if (!result && totalPasses > 0) {
2447                              result = -1; // round again
2448                              const CoinPresolveAction *paction = paction_;
2449                              while (paction) {
2450                                   const CoinPresolveAction *next = paction->next;
2451                                   delete paction;
2452                                   paction = next;
2453                              }
2454                              paction_ = NULL;
2455                         }
2456                    }
2457               }
2458          } else if (prob.status_) {
2459               // infeasible or unbounded
2460               result = 1;
2461               // Put status in nelems_!
2462               nelems_ = - prob.status_;
2463               originalModel->setProblemStatus(prob.status_);
2464          } else {
2465               // no changes - model needs restoring after Lou's changes
2466#ifndef CLP_NO_STD
2467               if (saveFile_ == "") {
2468#endif
2469                    delete presolvedModel_;
2470                    presolvedModel_ = new ClpSimplex(*originalModel);
2471                    // but we need to remove gaps
2472                    ClpPackedMatrix* clpMatrix =
2473                         dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
2474                    if (clpMatrix) {
2475                         clpMatrix->getPackedMatrix()->removeGaps();
2476                    }
2477#ifndef CLP_NO_STD
2478               } else {
2479                    presolvedModel_ = originalModel;
2480               }
2481               presolvedModel_->dropNames();
2482#endif
2483
2484               // drop integer information if wanted
2485               if (!keepIntegers)
2486                    presolvedModel_->deleteIntegerInformation();
2487               result = 2;
2488          }
2489     }
2490     if (result == 0 || result == 2) {
2491          int nrowsAfter = presolvedModel_->getNumRows();
2492          int ncolsAfter = presolvedModel_->getNumCols();
2493          CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
2494          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
2495                    messages)
2496                    << nrowsAfter << -(nrows_ - nrowsAfter)
2497                    << ncolsAfter << -(ncols_ - ncolsAfter)
2498                    << nelsAfter << -(nelems_ - nelsAfter)
2499                    << CoinMessageEol;
2500     } else {
2501          destroyPresolve();
2502          if (presolvedModel_ != originalModel_)
2503               delete presolvedModel_;
2504          presolvedModel_ = NULL;
2505     }
2506     return presolvedModel_;
2507}
2508
2509
Note: See TracBrowser for help on using the repository browser.