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

Last change on this file since 1900 was 1900, checked in by tkr, 7 years ago

Turning off spurious message during pre-processing

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