// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

/*
  Debug compile symbols for CoinPresolve.

  PRESOLVE_CONSISTENCY, PRESOLVE_DEBUG, and PRESOLVE_SUMMARY control
  consistency checking and debugging in the continuous presolve. See the
  comments in CoinPresolvePsdebug.hpp. DO NOT just define the symbols here in
  this file. Unless these symbols are consistent across all presolve code,
  you'll get something between garbage and a core dump.
*/

#include <stdio.h>

#include <cassert>
#include <iostream>

#include "CoinHelperFunctions.hpp"
#include "CoinFinite.hpp"

#include "CoinPackedMatrix.hpp"
#include "CoinWarmStartBasis.hpp"
#include "OsiSolverInterface.hpp"

#include "OsiPresolve.hpp"
#include "CoinPresolveMatrix.hpp"

#if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0 || PRESOLVE_SUMMARY > 0
#include "CoinPresolvePsdebug.hpp"
#include "CoinPresolveMonitor.hpp"
#endif
#include "CoinPresolveEmpty.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolveSingleton.hpp"
#include "CoinPresolveDoubleton.hpp"
#include "CoinPresolveTripleton.hpp"
#include "CoinPresolveZeros.hpp"
#include "CoinPresolveSubst.hpp"
#include "CoinPresolveForcing.hpp"
#include "CoinPresolveDual.hpp"
#include "CoinPresolveTighten.hpp"
#include "CoinPresolveUseless.hpp"
#include "CoinPresolveDupcol.hpp"
#include "CoinPresolveImpliedFree.hpp"
#include "CoinPresolveIsolated.hpp"
#include "CoinMessage.hpp"


OsiPresolve::OsiPresolve() :
  originalModel_(NULL),
  presolvedModel_(NULL),
  nonLinearValue_(0.0),
  originalColumn_(NULL),
  originalRow_(NULL),
  paction_(0),
  ncols_(0),
  nrows_(0),
  nelems_(0),
  presolveActions_(0),
  numberPasses_(5)
{
}

OsiPresolve::~OsiPresolve()
{
  gutsOfDestroy();
}
// Gets rid of presolve actions (e.g.when infeasible)
void 
OsiPresolve::gutsOfDestroy()
{
 const CoinPresolveAction *paction = paction_;
  while (paction) {
    const CoinPresolveAction *next = paction->next;
    delete paction;
    paction = next;
  }
  delete [] originalColumn_;
  delete [] originalRow_;
  paction_=NULL;
  originalColumn_=NULL;
  originalRow_=NULL;
}

/* This version of presolve returns a pointer to a new presolved 
   model.  NULL if infeasible

   doStatus controls activities required to transform an existing
   solution to match the presolved problem. I'd (lh) argue that this should
   default to false, but to maintain previous behaviour it defaults to true.
   Really, this is only useful if you've already optimised before applying
   presolve and also want to work with the solution after presolve.  I think
   that this is the less common case. The more common situation is to apply
   presolve before optimising.
*/
OsiSolverInterface * 
OsiPresolve::presolvedModel(OsiSolverInterface & si,
			    double feasibilityTolerance,
			    bool keepIntegers,
			    int numberPasses,
                            const char * prohibited,
			    bool doStatus,
			    const char * rowProhibited)
{
  ncols_ = si.getNumCols();
  nrows_ = si.getNumRows();
  nelems_ = si.getNumElements();
  numberPasses_ = numberPasses;

  double maxmin = si.getObjSense();
  originalModel_ = &si;
  delete [] originalColumn_;
  originalColumn_ = new int[ncols_];
  delete [] originalRow_;
  originalRow_ = new int[nrows_];
  int i;
  for (i=0;i<ncols_;i++) 
    originalColumn_[i]=i;
  for (i=0;i<nrows_;i++) 
    originalRow_[i]=i;

  // result is 0 - okay, 1 infeasible, -1 go round again
  int result = -1;
  
  // User may have deleted - its their responsibility
  presolvedModel_=NULL;
  // Messages
  CoinMessages msgs = CoinMessage(si.messages().language());
  // Only go round 100 times even if integer preprocessing
  int totalPasses=100;
  while (result == -1) {

    // make new copy
    delete presolvedModel_;
    presolvedModel_ = si.clone();
    totalPasses--;

    // drop integer information if wanted
    if (!keepIntegers) {
      int i;
      for (i=0;i<ncols_;i++)
	presolvedModel_->setContinuous(i);
    }

    
    CoinPresolveMatrix prob(ncols_,
			    maxmin,
			    presolvedModel_,
			    nrows_, nelems_,doStatus,nonLinearValue_,prohibited,
			    rowProhibited);
    // make sure row solution correct
    if (doStatus) {
      double *colels	= prob.colels_;
      int *hrow		= prob.hrow_;
      CoinBigIndex *mcstrt		= prob.mcstrt_;
      int *hincol		= prob.hincol_;
      int ncols		= prob.ncols_;
      
      
      double * csol = prob.sol_;
      double * acts = prob.acts_;
      int nrows = prob.nrows_;

      int colx;

      memset(acts,0,nrows*sizeof(double));
      
      for (colx = 0; colx < ncols; ++colx) {
	double solutionValue = csol[colx];
	for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
	  int row = hrow[i];
	  double coeff = colels[i];
	  acts[row] += solutionValue*coeff;
	}
      }
    }

    // move across feasibility tolerance
    prob.feasibilityTolerance_ = feasibilityTolerance;

/*
  Do presolve. Allow for the possibility that presolve might be ineffective
  (i.e., we're feasible but no postsolve actions are queued.
*/
    paction_ = presolve(&prob) ;
    result = 0 ; 
    // Get rid of useful arrays
    prob.deleteStuff();
/*
  This we don't need to do unless presolve actually reduced the system.
*/
    if (prob.status_==0&&paction_) {
      // Looks feasible but double check to see if anything slipped through
      int n		= prob.ncols_;
      double * lo = prob.clo_;
      double * up = prob.cup_;
      int i;
      
      for (i=0;i<n;i++) {
	if (up[i]<lo[i]) {
	  if (up[i]<lo[i]-1.0e-8) {
	    // infeasible
	    prob.status_=1;
	  } else {
	    up[i]=lo[i];
	  }
	}
      }
      
      n = prob.nrows_;
      lo = prob.rlo_;
      up = prob.rup_;

      for (i=0;i<n;i++) {
	if (up[i]<lo[i]) {
	  if (up[i]<lo[i]-1.0e-8) {
	    // infeasible
	    prob.status_=1;
	  } else {
	    up[i]=lo[i];
	  }
	}
      }
    }
/*
  If we're feasible, load the presolved system into the solver. Presumably we
  could skip model update and copying of status and solution if presolve took
  no action.
*/
    if (prob.status_ == 0) {
    
      prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);

# if PRESOLVE_CONSISTENCY > 0
      if (doStatus)
      { int basicCnt = 0 ;
	int basicColumns = 0;
	int i ;
	CoinPresolveMatrix::Status status ;
	for (i = 0 ; i < prob.ncols_ ; i++)
	{ status = prob.getColumnStatus(i);
	  if (status == CoinPrePostsolveMatrix::basic) basicColumns++ ; }
	basicCnt = basicColumns;
	for (i = 0 ; i < prob.nrows_ ; i++)
	{ status = prob.getRowStatus(i);
	  if (status == CoinPrePostsolveMatrix::basic) basicCnt++ ; }

# if PRESOLVE_DEBUG > 0
	presolve_check_nbasic(&prob) ;
# endif
	if (basicCnt>prob.nrows_) {
	  // Take out slacks
	  double * acts = prob.acts_;
	  double * rlo = prob.rlo_;
	  double * rup = prob.rup_;
	  double infinity = si.getInfinity();
	  for (i = 0 ; i < prob.nrows_ ; i++) {
	    status = prob.getRowStatus(i);
	    if (status == CoinPrePostsolveMatrix::basic) {
	      basicCnt-- ;
	      double down = acts[i]-rlo[i];
	      double up = rup[i]-acts[i];
	      if (CoinMin(up,down)<infinity) {
		if (down<=up)
		  prob.setRowStatus(i,CoinPrePostsolveMatrix::atLowerBound);
		else
		  prob.setRowStatus(i,CoinPrePostsolveMatrix::atUpperBound);
	      } else {
		prob.setRowStatus(i,CoinPrePostsolveMatrix::isFree);
	      }
	    }
	    if (basicCnt==prob.nrows_)
	      break;
	  }
	}
      }
#endif

/*
  Install the status and primal solution, if we've been carrying them along.

  The code that copies status is efficient but brittle. The current definitions
  for CoinWarmStartBasis::Status and CoinPrePostsolveMatrix::Status are in
  one-to-one correspondence. This code will fail if that ever changes.
*/
      if (doStatus) {
	presolvedModel_->setColSolution(prob.sol_);
	CoinWarmStartBasis *basis = 
	  dynamic_cast<CoinWarmStartBasis *>(presolvedModel_->getEmptyWarmStart());
	basis->resize(prob.nrows_,prob.ncols_);
	int i;
	for (i=0;i<prob.ncols_;i++) {
	  CoinWarmStartBasis::Status status = 
	    static_cast<CoinWarmStartBasis::Status> (prob.getColumnStatus(i));
	  basis->setStructStatus(i,status);
	}
	for (i=0;i<prob.nrows_;i++) {
	  CoinWarmStartBasis::Status status = 
	    static_cast<CoinWarmStartBasis::Status> (prob.getRowStatus(i));
	  basis->setArtifStatus(i,status);
	}
	presolvedModel_->setWarmStart(basis);
	delete basis ;
	delete [] prob.sol_;
	delete [] prob.acts_;
	delete [] prob.colstat_;
	prob.sol_=NULL;
	prob.acts_=NULL;
	prob.colstat_=NULL;
      }
/*
  Copy original column and row information from the CoinPresolveMatrix object
  so it'll be available for postsolve.
*/
      int ncolsNow = presolvedModel_->getNumCols();
      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
      delete [] prob.originalColumn_;
      prob.originalColumn_=NULL;
      int nrowsNow = presolvedModel_->getNumRows();
      memcpy(originalRow_,prob.originalRow_,nrowsNow*sizeof(int));
      delete [] prob.originalRow_;
      prob.originalRow_=NULL;

      // now clean up integer variables.  This can modify original
      {
	int numberChanges=0;
	const double * lower0 = originalModel_->getColLower();
	const double * upper0 = originalModel_->getColUpper();
	const double * lower = presolvedModel_->getColLower();
	const double * upper = presolvedModel_->getColUpper();
	for (i=0;i<ncolsNow;i++) {
	  if (!presolvedModel_->isInteger(i))
	    continue;
	  int iOriginal = originalColumn_[i];
	  double lowerValue0 = lower0[iOriginal];
	  double upperValue0 = upper0[iOriginal];
	  double lowerValue = ceil(lower[i]-1.0e-5);
	  double upperValue = floor(upper[i]+1.0e-5);
	  presolvedModel_->setColBounds(i,lowerValue,upperValue);
	  if (lowerValue>upperValue) {
	    numberChanges++;
	    CoinMessageHandler *hdlr = presolvedModel_->messageHandler() ;
	    hdlr->message(COIN_PRESOLVE_COLINFEAS,msgs)
	        << iOriginal << lowerValue << upperValue << CoinMessageEol ;
	    result=1;
	  } else {
	    if (lowerValue>lowerValue0+1.0e-8) {
	      originalModel_->setColLower(iOriginal,lowerValue);
	      numberChanges++;
	    }
	    if (upperValue<upperValue0-1.0e-8) {
	      originalModel_->setColUpper(iOriginal,upperValue);
	      numberChanges++;
	    }
	  }	  
	}
	if (numberChanges) {
	  CoinMessageHandler *hdlr = presolvedModel_->messageHandler() ;
	  hdlr->message(COIN_PRESOLVE_INTEGERMODS,msgs)
	      << numberChanges << CoinMessageEol;
	  // we can't go round again in integer if dupcols
	  if (!result && totalPasses > 0 &&
	      (prob.presolveOptions_&0x80000000) == 0) {
	    result = -1; // round again
	    const CoinPresolveAction *paction = paction_;
	    while (paction) {
	      const CoinPresolveAction *next = paction->next;
	      delete paction;
	      paction = next;
	    }
	    paction_=NULL;
	  }
	}
      }
    } else if (prob.status_ != 0) { 
      // infeasible or unbounded
      result = 1 ;
    }
  }
  if (!result) {
    int nrowsAfter = presolvedModel_->getNumRows();
    int ncolsAfter = presolvedModel_->getNumCols();
    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
    CoinMessageHandler *hdlr = presolvedModel_->messageHandler() ;
    hdlr->message(COIN_PRESOLVE_STATS,msgs)
        << nrowsAfter << -(nrows_-nrowsAfter)
	<< ncolsAfter << -(ncols_-ncolsAfter)
        << nelsAfter << -(nelems_-nelsAfter) << CoinMessageEol ;
  } else {
    gutsOfDestroy();
    delete presolvedModel_;
    presolvedModel_=NULL;
  }
  return presolvedModel_;
}

// Return pointer to presolved model
OsiSolverInterface * 
OsiPresolve::model() const
{
  return presolvedModel_;
}
// Return pointer to original model
OsiSolverInterface * 
OsiPresolve::originalModel() const
{
  return originalModel_;
}

void 
OsiPresolve::postsolve(bool updateStatus)
{
  // Messages
  CoinMessages msgs = CoinMessage(presolvedModel_->messages().language()) ;
  CoinMessageHandler *hdlr = presolvedModel_->messageHandler() ;
  if (!presolvedModel_->isProvenOptimal()) {
    hdlr->message(COIN_PRESOLVE_NONOPTIMAL,msgs) << CoinMessageEol ;
  }

  // this is the size of the original problem
  const int ncols0  = ncols_ ;
  const int nrows0  = nrows_ ;
  const CoinBigIndex nelems0 = nelems_ ;

  // reality check
  assert(ncols0 == originalModel_->getNumCols()) ;
  assert(nrows0 == originalModel_->getNumRows()) ;

  // this is the reduced problem
  int ncols = presolvedModel_->getNumCols() ;
  int nrows = presolvedModel_->getNumRows() ;

  double *acts = new double [nrows0] ;
  double *sol = new double [ncols0] ;
  CoinZeroN(acts,nrows0) ;
  CoinZeroN(sol,ncols0) ;
  
  unsigned char *rowstat = NULL ;
  unsigned char *colstat = NULL ;
  CoinWarmStartBasis *presolvedBasis  = 
    dynamic_cast<CoinWarmStartBasis*>(presolvedModel_->getWarmStart()) ;
  if (!presolvedBasis) updateStatus = false ;
  if (updateStatus) {
    colstat = new unsigned char[ncols0+nrows0] ;
#   ifdef ZEROFAULT
    memset(colstat,0,((ncols0+nrows0)*sizeof(char))) ;
#   endif
    rowstat = colstat+ncols0 ;
    for (int i = 0 ; i < ncols ; i++) {
      colstat[i] = presolvedBasis->getStructStatus(i) ;
    }
    for (int i = 0 ; i < nrows ; i++) {
      rowstat[i] = presolvedBasis->getArtifStatus(i) ;
    }
  } 
  delete presolvedBasis ;

# if PRESOLVE_CONSISTENCY > 0
  if (updateStatus)
  { int basicCnt = 0 ;
    for (int i = 0 ; i < ncols ; i++)
    { if (colstat[i] == CoinWarmStartBasis::basic) basicCnt++ ; }
    for (int i = 0 ; i < nrows ; i++)
    { if (rowstat[i] == CoinWarmStartBasis::basic) basicCnt++ ; }

    assert (basicCnt == nrows) ;
  }
# endif

/*
  Postsolve back to the original problem.  The CoinPostsolveMatrix object
  assumes ownership of sol, acts, colstat, and rowstat.
*/
  CoinPostsolveMatrix prob(presolvedModel_,ncols0,nrows0,nelems0,
			   presolvedModel_->getObjSense(),
			   sol,acts,colstat,rowstat) ;
  postsolve(prob) ;

# if PRESOLVE_CONSISTENCY > 0
  if (updateStatus)
  { int basicCnt = 0 ;
    for (int i = 0 ; i < ncols0 ; i++)
    { if (prob.getColumnStatus(i) == CoinPrePostsolveMatrix::basic)
      basicCnt++ ; }
    for (int i = 0 ; i < nrows0 ; i++)
    { if (prob.getRowStatus(i) == CoinPrePostsolveMatrix::basic)
      basicCnt++ ; }

    assert (basicCnt == nrows0) ;
  }
# endif

  originalModel_->setColSolution(sol) ;
  if (updateStatus) {
    CoinWarmStartBasis *basis = 
      dynamic_cast<CoinWarmStartBasis *>(presolvedModel_->getEmptyWarmStart()) ;
    basis->setSize(ncols0,nrows0) ;
    for (int i = 0 ; i < ncols0 ; i++) {
      CoinWarmStartBasis::Status status =
          static_cast<CoinWarmStartBasis::Status>(prob.getColumnStatus(i)) ;
      assert(status != CoinWarmStartBasis::atLowerBound || originalModel_->getColLower()[i] > -originalModel_->getInfinity()) ;
      assert(status != CoinWarmStartBasis::atUpperBound || originalModel_->getColUpper()[i] <  originalModel_->getInfinity()) ;
      basis->setStructStatus(i,status);
    }

# if PRESOLVE_DEBUG > 0
  /*
    Do a thorough check of row and column solutions. There should be no
    inconsistencies at this point.
  */
  std::cout
    << "Checking solution before transferring basis." << std::endl ;
  presolve_check_sol(&prob,2,2,2) ;
  int errs = 0 ;
# endif
    for (int i = 0 ; i < nrows0 ; i++) {
      CoinWarmStartBasis::Status status =
          static_cast<CoinWarmStartBasis::Status>(prob.getRowStatus(i)) ;
      basis->setArtifStatus(i,status);
    }
    originalModel_->setWarmStart(basis);
    delete basis ;
  } 

}

// return pointer to original columns
const int * 
OsiPresolve::originalColumns() const
{
  return originalColumn_;
}
// return pointer to original rows
const int * 
OsiPresolve::originalRows() const
{
  return originalRow_;
}
// Set pointer to original model
void 
OsiPresolve::setOriginalModel(OsiSolverInterface * model)
{
  originalModel_=model;
}

#if 0
// A lazy way to restrict which transformations are applied
// during debugging.
static int ATOI(const char *name)
{
 return true;
#if	PRESOLVE_DEBUG > 0 || PRESOLVE_SUMMARY > 0
  if (getenv(name)) {
    int val = atoi(getenv(name));
    printf("%s = %d\n", name, val);
    return (val);
  } else {
    if (strcmp(name,"off"))
      return (true);
    else
      return (false);
  }
#else
  return (true);
#endif
}
#endif

#if PRESOLVE_DEBUG > 0
// Anonymous namespace for debug routines
namespace {

/*
  A control routine for debug checks --- keeps down the clutter in doPresolve.
  Each time it's called, it prints a list of transforms applied since the last
  call, then does checks.
*/

void check_and_tell (const CoinPresolveMatrix *const prob,
		     const CoinPresolveAction *first,
		     const CoinPresolveAction *&mark)

{ const CoinPresolveAction *current ;

  if (first != mark)
  { printf("PRESOLVE: applied") ;
    for (current = first ;
	 current != mark && current != 0 ;
	 current = current->next)
    { printf(" %s",current->name()) ; }
    printf("\n") ;

    presolve_check_sol(prob) ;
    presolve_check_nbasic(prob) ;
    mark = first ; }

  return ; }

/*
  At a guess, this code is intended to allow a known solution to be checked
  against presolve progress. Pulled it into the local debug namespace, but
  really should be integrated with CoinPresolvePsdebug.  At the least, needs
  a method to conveniently set debugSolution.
  -- lh, 110605 --
*/
double *debugSolution = NULL ;
int debugNumberColumns = -1 ;
int counter = 1000000 ;

bool break2 (CoinPresolveMatrix *prob) {
  if (counter > 0)
    printf("break2: counter %d\n",counter) ;
  counter-- ;
  if (debugSolution && prob->ncols_ == debugNumberColumns) {
    for (int i = 0 ; i < prob->ncols_ ; i++) {
      double value = debugSolution[i] ;
      if (value < prob->clo_[i]) {
	printf("%d inf %g %g %g\n",i,prob->clo_[i],value,prob->cup_[i]) ;
      } else if (value > prob->cup_[i]) {
	printf("%d inf %g %g %g\n",i,prob->clo_[i],value,prob->cup_[i]) ;
      }
    }
  }
  if (!counter) {
    printf("skipping next and all\n") ;
  }
  return (counter <= 0) ;
}

} // end anonymous namespace for debug routines
#endif

#if PRESOLVE_DEBUG > 0
# define possibleBreak if (break2(prob)) break
# define possibleSkip  if (!break2(prob)) 
#else
# define possibleBreak
# define possibleSkip
#endif

// This is the presolve loop.
// It is a separate virtual function so that it can be easily
// customized by subclassing CoinPresolve.

const CoinPresolveAction *OsiPresolve::presolve(CoinPresolveMatrix *prob)
{
  paction_ = 0 ;

  prob->status_ = 0 ; // say feasible

# if PRESOLVE_DEBUG > 0
  const CoinPresolveAction *pactiond = 0 ;
  presolve_check_sol(prob,2,1,1) ;

  // CoinPresolveMonitor *monitor = new CoinPresolveMonitor(prob,true,22) ;
  CoinPresolveMonitor *monitor = 0 ;
# endif
/*
  Transfer costs off of singleton variables, and also between integer
  variables when advantageous.

  transferCosts is defined in CoinPresolveFixed.cpp
*/
  if ((presolveActions_&0x04) != 0) {
    transferCosts(prob) ;
#   if PRESOLVE_DEBUG > 0
    if (monitor) monitor->checkAndTell(prob) ;
#   endif
  }
/*
  Fix variables before we get into the main transform loop.
*/
  paction_ = make_fixed(prob,paction_) ;

# if PRESOLVE_DEBUG > 0
  check_and_tell(prob,paction_,pactiond) ;
  if (monitor) monitor->checkAndTell(prob) ;
# endif

  // if integers then switch off dual stuff
  // later just do individually
  bool doDualStuff = true ;
  if ((presolveActions_&0x01) == 0) {
    int ncol = presolvedModel_->getNumCols() ;
    for (int i = 0 ; i < ncol ; i++)
      if (presolvedModel_->isInteger(i))
	doDualStuff = false ;
  }
  

# if PRESOLVE_CONSISTENCY > 0
  presolve_links_ok(prob) ;
# endif

/*
  If we're feasible, set up for the main presolve transform loop.
*/
  if (!prob->status_) {
# if 0
/*
  This block is used during debugging. See ATOI to see how it works. Some
  editing will be required to turn it all on.
*/
    bool slackd = ATOI("SLACKD")!=0;
    //bool forcing = ATOI("FORCING")!=0;
    bool doubleton = ATOI("DOUBLETON")!=0;
    bool forcing = ATOI("off")!=0;
    bool ifree = ATOI("off")!=0;
    bool zerocost = ATOI("off")!=0;
    bool dupcol = ATOI("off")!=0;
    bool duprow = ATOI("off")!=0;
    bool dual = ATOI("off")!=0;
# else
# if 1
    // normal operation --- all transforms enabled
    bool slackSingleton = true;
    bool slackd = true;
    bool doubleton = true;
    bool tripleton = true;
    bool forcing = true;
    bool ifree = true;
    bool zerocost = true;
    bool dupcol = true;
    bool duprow = true;
    bool dual = doDualStuff;
# else
    // compile time selection of transforms.
    bool slackSingleton = true;
    bool slackd = false;
    bool doubleton = true;
    bool tripleton = true;
    bool forcing = true;
    bool ifree = false;
    bool zerocost = false;
    bool dupcol = false;
    bool duprow = false;
    bool dual = false;
# endif
# endif
/*
  Process OsiPresolve options. Set corresponding CoinPresolve options and
  control variables here.
*/
    // Switch off some stuff if would annoy set partitioning etc
    if ((presolveActions_&0x02) != 0) {
      doubleton = false;
      tripleton = false;
      ifree = false;
    }
    // stop x+y+z=1
    if ((presolveActions_&0x08) != 0)
      prob->setPresolveOptions(prob->presolveOptions()|0x04) ;
    // switch on stuff which can't be unrolled easily
    if ((presolveActions_&0x10) != 0)
      prob->setPresolveOptions(prob->presolveOptions()|0x10) ;
    // switch on gub stuff (unimplemented as of 110605 -- lh --)
    if ((presolveActions_&0x20) != 0)
      prob->setPresolveOptions(prob->presolveOptions()|0x20) ;
    // allow duplicate column processing for integer columns
    if ((presolveActions_&0x01) != 0)
      prob->setPresolveOptions(prob->presolveOptions()|0x01) ;
/*
  Set [rows,cols]ToDo to process all rows & cols unless there are
  specific prohibitions.
*/
    prob->initColsToDo() ;
    prob->initRowsToDo() ;
/*
  Try to remove duplicate rows and columns.
*/
    if (dupcol) {
      possibleSkip ;
      paction_ = dupcol_action::presolve(prob,paction_) ;
#     if PRESOLVE_DEBUG > 0
      if (monitor) monitor->checkAndTell(prob) ;
#     endif
    }
    if (duprow) {
      possibleSkip ;
      paction_ = duprow_action::presolve(prob,paction_) ;
#     if PRESOLVE_DEBUG > 0
      if (monitor) monitor->checkAndTell(prob) ;
#     endif
    }
/*
  The main loop starts with a minor loop that does inexpensive presolve
  transforms until convergence. At each iteration of this loop,
  next[Rows,Cols]ToDo is copied over to [rows,cols]ToDo.

  Then there's a block to set [rows,cols]ToDo to examine all rows & cols,
  followed by executions of expensive transforms. Then we come back around for
  another iteration of the main loop. [rows,cols]ToDo is not reset as we come
  back around, so we dive into the inexpensive loop set up to process all.

  lastDropped is a count of total number of rows dropped by presolve. Used as
  an additional criterion to end the main presolve loop.
*/
    int lastDropped = 0 ;
    prob->pass_ = 0 ;
    for (int iLoop = 0 ; iLoop < numberPasses_ ; iLoop++) {

#     if PRESOLVE_SUMMARY > 0
      std::cout << "Starting major pass " << (iLoop+1) << std::endl ;
#     endif

      const CoinPresolveAction *const paction0 = paction_ ;
// #define IMPLIED 3
#ifdef IMPLIED
      int fill_level = 3 ;
# define IMPLIED2 1
# if IMPLIED != 3
#   if IMPLIED > 0 && IMPLIED < 11
      fill_level = IMPLIED ;
      printf("** fill_level == %d !\n",fill_level) ;
#   endif
#   if IMPLIED > 11 && IMPLIED < 21
      fill_level = -(IMPLIED-10) ;
      printf("** fill_level == %d !\n",fill_level);
#   endif
# endif
#else
      // look for substitutions with no fill
      int fill_level = 2 ;
#endif
      int whichPass = 0 ;
/*
  Apply inexpensive transforms until convergence or infeasible/unbounded.
*/
      while (true) {
	whichPass++ ;
	prob->pass_++ ;
	const CoinPresolveAction *const paction1 = paction_ ;

	if (slackd) {
	  bool notFinished = true ;
	  while (notFinished) {
	    possibleBreak ;
	    paction_ =
	       slack_doubleton_action::presolve(prob,paction_,notFinished) ;
	  }
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

	if (dual && whichPass == 1) {
	  possibleBreak;
	  // this can also make E rows so do one bit here
	  paction_ = remove_dual_action::presolve(prob,paction_) ;
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

	if (doubleton) {
	  possibleBreak ;
	  paction_ = doubleton_action::presolve(prob,paction_) ;
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

	if (tripleton) {
	  possibleBreak ;
	  paction_ = tripleton_action::presolve(prob,paction_) ;
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

	if (zerocost) {
	  possibleBreak ;
	  paction_ = do_tighten_action::presolve(prob,paction_) ;
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

	if (forcing) {
	  possibleBreak;
	  paction_ = forcing_constraint_action::presolve(prob, paction_);
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

	if (ifree && (whichPass%5) == 1) {
	  possibleBreak ;
	  paction_ = implied_free_action::presolve(prob,paction_,fill_level) ;
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	}

#	if PRESOLVE_CONSISTENCY > 0
	presolve_links_ok(prob) ;
	presolve_no_zeros(prob) ;
	presolve_consistent(prob) ;
#	endif
/*
  Set up for next pass.
  
  Original comment adds: later do faster if many changes i.e. memset and memcpy
*/
        prob->stepRowsToDo() ;

#	if PRESOLVE_DEBUG > 0
	int rowCheck = -1 ;
	bool rowFound = false ;
	for (int i = 0 ; i < prob->numberRowsToDo_ ; i++) {
	  int index = prob->rowsToDo_[i];
	  if (index == rowCheck) {
	    std::cout
	      << "  row " << index << " on list after pass " << whichPass
	      << std::endl ;
	    rowFound = true ;
	  }
	}
	if (!rowFound && rowCheck >= 0)
	  prob->rowsToDo_[prob->numberRowsToDo_++] = rowCheck ;
#	endif

	prob->stepColsToDo() ;

#	if PRESOLVE_DEBUG > 0
	int colCheck = -1 ;
	bool colFound = false ;
	for (int i = 0 ; i < prob->numberNextColsToDo_ ; i++) {
	  int index = prob->colsToDo_[i] ;
	  if (index == colCheck) {
	    std::cout
	      << "  col " << index << " on list after pass " << whichPass
	      << std::endl ;
	    colFound = true ;
	  }
	}
	if (!colFound && colCheck >= 0)
	  prob->colsToDo_[prob->numberColsToDo_++] = colCheck ;
#	endif
/*
  Break if nothing happened (no postsolve actions queued).

  The check for fill_level > 0 is a hack to allow repeating the loop with some
  modified fill level (playing with negative values).
  
  fill_level = 0 (as set in other places) will clearly be a problem.
  -- lh, 110605 --
*/
	if (paction_ == paction1 && fill_level > 0) break ;
      }
/*
  End of inexpensive transform loop.
  Reset [rows,cols]ToDo to process all rows and columns unless there are
  specfic prohibitions.
*/
      prob->initRowsToDo() ;
      prob->initColsToDo() ;
/*
  Try expensive presolve transforms.

  Original comment adds: this caused world.mps to run into numerical
  			 difficulties
*/
#     if PRESOLVE_SUMMARY > 0
      std::cout << "Starting expensive." << std::endl ;
#     endif
/*
  Try and fix variables at upper or lower bound by calculating bounds on the
  dual variables and propagating them to the reduced costs. Every other
  iteration, see if this has created free variables.
*/
      if (dual) {
	for (int itry = 0 ; itry < 5 ; itry++) {
	  const CoinPresolveAction *const paction2 = paction_ ;
	  possibleBreak ;
	  paction_ = remove_dual_action::presolve(prob,paction_) ;
#	  if PRESOLVE_DEBUG > 0
	  check_and_tell(prob,paction_,pactiond) ;
	  if (monitor) monitor->checkAndTell(prob) ;
#	  endif
	  if (prob->status_) break ;
	  if (ifree) {
#ifdef IMPLIED
# if IMPLIED2 == 0
	    int fill_level = 0 ; // switches off substitution
# elif IMPLIED2 != 99
	    int fill_level = IMPLIED2 ;
# endif
#endif
	    if ((itry&1) == 0) {
	      possibleBreak ;
	      paction_ =
	          implied_free_action::presolve(prob,paction_,fill_level) ;
	    }
#	    if PRESOLVE_DEBUG > 0
	    check_and_tell(prob,paction_,pactiond) ;
	    if (monitor) monitor->checkAndTell(prob) ;
#	    endif
	    if (prob->status_) break ;
	  }
	  if (paction_ == paction2) break ;
	}
      } else if (ifree) {
/*
  Just check for free variables.
*/
#ifdef IMPLIED
# if IMPLIED2 == 0
	int fill_level = 0 ; // switches off substitution
# elif IMPLIED2 != 99
	int fill_level = IMPLIED2 ;
# endif
#endif
	possibleBreak ;
	paction_ = implied_free_action::presolve(prob,paction_,fill_level) ;
#	if PRESOLVE_DEBUG > 0
	check_and_tell(prob,paction_,pactiond) ;
	if (monitor) monitor->checkAndTell(prob) ;
#	endif
	if (prob->status_) break ;
      }
/*
  Check if other transformations have produced duplicate rows or columns.
*/
      if (dupcol) {
	possibleBreak ;
	paction_ = dupcol_action::presolve(prob,paction_) ;
#	if PRESOLVE_DEBUG > 0
	check_and_tell(prob,paction_,pactiond) ;
	if (monitor) monitor->checkAndTell(prob) ;
#	endif
	if (prob->status_) break ;
      }
      if (duprow) {
	possibleBreak ;
	paction_ = duprow_action::presolve(prob,paction_) ;
#	if PRESOLVE_DEBUG > 0
	check_and_tell(prob,paction_,pactiond) ;
	if (monitor) monitor->checkAndTell(prob) ;
#	endif
	if (prob->status_) break ;
      }
      // Will trigger abort due to unimplemented postsolve  -- lh, 110605 --
      if ((presolveActions_&0x20) != 0) {
	possibleBreak ;
	paction_ = gubrow_action::presolve(prob,paction_) ;
      }
/*
  Count the number of empty rows and see if we've made progress in this pass.
*/
      bool stopLoop = false ;
      {
	const int *const hinrow = prob->hinrow_ ;
	int numberDropped = 0 ;
	for (int i = 0 ; i < nrows_ ; i++)
	  if (!hinrow[i]) numberDropped++ ;
#	if PRESOLVE_DEBUG > 0
	std::cout
	  << "  " << (numberDropped-lastDropped)
	  << " rows dropped in pass " << iLoop << "." << std::endl ;
#	endif
	if (numberDropped == lastDropped)
	  stopLoop = true ;
	else
	  lastDropped = numberDropped ;
      }
/*
  Check for singleton variables that can act like a logical, allowing a
  row to be transformed from an equality to an inequality.

  The third parameter allows for costs for the existing logicals. This
  is apparently used by clp; consult the clp presolve before implementing
  it here.  -- lh, 110605 --

  Original comment: Do this here as not very loopy
*/
      if (slackSingleton) {
	possibleBreak ;
	paction_ = slack_singleton_action::presolve(prob,paction_,NULL) ;
#	if PRESOLVE_DEBUG > 0
	check_and_tell(prob,paction_,pactiond) ;
	if (monitor) monitor->checkAndTell(prob) ;
#	endif
      }
#     if PRESOLVE_DEBUG > 0
      presolve_check_sol(prob,1) ;
#     endif

      if (paction_ == paction0 || stopLoop) break ;
	  
    } // End of major pass loop
  }
/*
  Final cleanup: drop zero coefficients from the matrix, then drop empty rows
  and columns.
*/
  if (!prob->status_) {
    paction_ = drop_zero_coefficients(prob,paction_) ;
#   if PRESOLVE_DEBUG > 0
    check_and_tell(prob,paction_,pactiond) ;
    if (monitor) monitor->checkAndTell(prob) ;
#   endif

    paction_ = drop_empty_cols_action::presolve(prob,paction_) ;
#   if PRESOLVE_DEBUG > 0
    check_and_tell(prob,paction_,pactiond) ;
#   endif

    paction_ = drop_empty_rows_action::presolve(prob,paction_) ;
#   if PRESOLVE_DEBUG > 0
    check_and_tell(prob,paction_,pactiond) ;
#   endif
  }
/*
  Not feasible? Say something and clean up.
*/
  CoinMessageHandler *hdlr = prob->messageHandler() ;
  CoinMessages msgs = CoinMessage(prob->messages().language());
  if (prob->status_) {
    if (prob->status_ == 1)
      hdlr->message(COIN_PRESOLVE_INFEAS,msgs)
	<< prob->feasibilityTolerance_ << CoinMessageEol ;
    else if (prob->status_ == 2)
      hdlr->message(COIN_PRESOLVE_UNBOUND,msgs) << CoinMessageEol ;
    else
      hdlr->message(COIN_PRESOLVE_INFEASUNBOUND,msgs) << CoinMessageEol ;
    gutsOfDestroy() ;
  }
  return (paction_) ;
}


/*
  We could have implemented this by having each postsolve routine directly
  call the next one, but this makes it easier to add debugging checks.
*/
void OsiPresolve::postsolve (CoinPostsolveMatrix &prob)
{
  const CoinPresolveAction *paction = paction_;

# if PRESOLVE_DEBUG > 0
  std::cout << "Begin POSTSOLVING." << std::endl ;
  if (prob.colstat_)
  { presolve_check_nbasic(&prob) ;
    presolve_check_sol(&prob,2,2,2) ; }
  presolve_check_duals(&prob) ;
# endif
  
  
  while (paction) {
#   if PRESOLVE_DEBUG > 0
    std::cout << "POSTSOLVING " << paction->name() << std::endl ;
#   endif

    paction->postsolve(&prob);
    
#   if PRESOLVE_DEBUG > 0
    if (prob.colstat_) {
      presolve_check_nbasic(&prob) ;
      presolve_check_sol(&prob,2,2,2) ;
    }
#   endif
    paction = paction->next ;
#   if PRESOLVE_DEBUG > 0
    presolve_check_duals(&prob);
#   endif
  }    
# if PRESOLVE_DEBUG > 0
    std::cout << "End POSTSOLVING" << std::endl ;
# endif
  
# if PRESOLVE_DEBUG > 0
  for (int j = 0 ; j < prob.ncols_ ; j++) {
    if (!prob.cdone_[j]) {
      printf("!cdone[%d]\n", j) ;
      abort() ;
    }
  }
  for (int i = 0 ; i < prob.nrows_ ; i++) {
    if (!prob.rdone_[i]) {
      printf("!rdone[%d]\n", i) ;
      abort() ;
    }
  }
  for (int j = 0 ; j < prob.ncols_ ; j++) {
    if (prob.sol_[j] < -1e10 || prob.sol_[j] > 1e10)
      printf("!!!%d %g\n",j,prob.sol_[j]) ;
  }
# endif

  /*
    Put back duals. Flip sign for maximisation problems.
  */
  double maxmin = originalModel_->getObjSense() ;
  if (maxmin < 0.0) {
    double *pi = prob.rowduals_ ;
    for (int i = 0 ; i < nrows_ ; i++) pi[i] = -pi[i] ;
  }
  originalModel_->setRowPrice(prob.rowduals_) ;
}


static inline double getTolerance(const OsiSolverInterface  *si, OsiDblParam key)
{
  double tol;
  if (!si->getDblParam(key,tol)) {
    CoinPresolveAction::throwCoinError("getDblParam failed",
		    "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix") ;
  }
  return (tol) ;
}


// Assumptions:
// 1. nrows>=m.getNumRows()
// 2. ncols>=m.getNumCols()
//
// In presolve, these values are equal.
// In postsolve, they may be inequal, since the reduced problem
// may be smaller, but we need room for the large problem.
// ncols may be larger than si.getNumCols() in postsolve,
// this at that point si will be the reduced problem,
// but we need to reserve enough space for the original problem.

CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const OsiSolverInterface * si,
					     int ncols_in,
					     int nrows_in,
					     CoinBigIndex nelems_in) :
  ncols_(si->getNumCols()),
  nelems_(si->getNumElements()),
  ncols0_(ncols_in),
  nrows0_(nrows_in),
  bulkRatio_(2.0),

  mcstrt_(new CoinBigIndex[ncols_in+1]),
  hincol_(new int[ncols_in+1]),

  cost_(new double[ncols_in]),
  clo_(new double[ncols_in]),
  cup_(new double[ncols_in]),
  rlo_(new double[nrows_in]),
  rup_(new double[nrows_in]),
  originalColumn_(new int[ncols_in]),
  originalRow_(new int[nrows_in]),

  ztolzb_(getTolerance(si, OsiPrimalTolerance)),
  ztoldj_(getTolerance(si, OsiDualTolerance)),

  maxmin_(si->getObjSense()),

  handler_(0),
  defaultHandler_(false),
  messages_()

{
  bulk0_ = static_cast<CoinBigIndex>(bulkRatio_*nelems_in) ;
  hrow_ = new int [bulk0_] ;
  colels_ = new double[bulk0_] ;

  si->getDblParam(OsiObjOffset,originalOffset_);
  int ncols = si->getNumCols();
  int nrows = si->getNumRows();

  setMessageHandler(si->messageHandler()) ;

  CoinDisjointCopyN(si->getColLower(), ncols, clo_);
  CoinDisjointCopyN(si->getColUpper(), ncols, cup_);
  CoinDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
  CoinDisjointCopyN(si->getRowLower(), nrows,  rlo_);
  CoinDisjointCopyN(si->getRowUpper(), nrows,  rup_);
  int i;
  // initialize and clean up bounds
  double infinity = si->getInfinity();
  if (infinity!=COIN_DBL_MAX) {
    for (i=0;i<ncols;i++) {
      if (clo_[i]==-infinity)
	clo_[i]=-COIN_DBL_MAX;
      if (cup_[i]==infinity)
	cup_[i]=COIN_DBL_MAX;
    }
    for (i=0;i<nrows;i++) {
      if (rlo_[i]==-infinity)
	rlo_[i]=-COIN_DBL_MAX;
      if (rup_[i]==infinity)
	rup_[i]=COIN_DBL_MAX;
    }
  }
  for (i=0;i<ncols_in;i++) 
    originalColumn_[i]=i;
  for (i=0;i<nrows_in;i++) 
    originalRow_[i]=i;
  sol_=NULL;
  rowduals_=NULL;
  acts_=NULL;

  rcosts_=NULL;
  colstat_=NULL;
  rowstat_=NULL;
}

// I am not familiar enough with CoinPackedMatrix to be confident
// that I will implement a row-ordered version of toColumnOrderedGapFree
// properly.
static bool isGapFree(const CoinPackedMatrix& matrix)
{
  const CoinBigIndex * start = matrix.getVectorStarts();
  const int * length = matrix.getVectorLengths();
  int i;
  for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
    if (start[i+1] - start[i] != length[i])
      break;
  }
  return (! (i >= 0));
}

CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
				       double maxmin,
				       // end prepost members
				       OsiSolverInterface * si,
				       // rowrep
				       int nrows_in,
				       CoinBigIndex nelems_in,
				       bool doStatus,
				       double nonLinearValue,
                                       const char * prohibited,
				       const char * rowProhibited)
  : CoinPrePostsolveMatrix(si,ncols0_in,nrows_in,nelems_in),
    clink_(new presolvehlink[ncols0_in+1]),
    rlink_(new presolvehlink[nrows_in+1]),
    dobias_(0.0),

    // temporary init
    mrstrt_(new CoinBigIndex[nrows_in+1]),
    hinrow_(new int[nrows_in+1]),
    integerType_(new unsigned char[ncols0_in]),
    tuning_(false),
    startTime_(0.0),
    feasibilityTolerance_(0.0),
    status_(-1),
    maxSubstLevel_(3),
    colsToDo_(new int [ncols0_in]),
    numberColsToDo_(0),
    nextColsToDo_(new int[ncols0_in]),
    numberNextColsToDo_(0),
    rowsToDo_(new int [nrows_in]),
    numberRowsToDo_(0),
    nextRowsToDo_(new int[nrows_in]),
    numberNextRowsToDo_(0),
    presolveOptions_(0)
{

  rowels_ = new double [bulk0_] ;
  hcol_ = new int [bulk0_] ;

  nrows_ = si->getNumRows() ;
  const CoinBigIndex bufsize = static_cast<CoinBigIndex>(bulkRatio_*nelems_in) ;

  // Set up change bits
  rowChanged_ = new unsigned char[nrows_];
  memset(rowChanged_,0,nrows_);
  colChanged_ = new unsigned char[ncols_];
  memset(colChanged_,0,ncols_);
  const CoinPackedMatrix * m1 = si->getMatrixByCol();

  // The coefficient matrix is a big hunk of stuff.
  // Do the copy here to try to avoid running out of memory.

  const CoinBigIndex * start = m1->getVectorStarts();
  const int * length = m1->getVectorLengths();
  const int * row = m1->getIndices();
  const double * element = m1->getElements();
  int icol,nel=0;
  mcstrt_[0]=0;
  for (icol=0;icol<ncols_;icol++) {
    int j;
    for (j=start[icol];j<start[icol]+length[icol];j++) {
      if (fabs(element[j])>ZTOLDP) {
        hrow_[nel]=row[j];
	colels_[nel++]=element[j];
      }
    }
    hincol_[icol]=nel-mcstrt_[icol];
    mcstrt_[icol+1]=nel;
  }

  // same thing for row rep
  CoinPackedMatrix * m = new CoinPackedMatrix();
  m->reverseOrderedCopyOf(*si->getMatrixByCol());
  // do by hand because of zeros m->removeGaps();
  CoinDisjointCopyN(m->getVectorStarts(),  nrows_,  mrstrt_);
  mrstrt_[nrows_] = nelems_;
  CoinDisjointCopyN(m->getVectorLengths(), nrows_,  hinrow_);
  CoinDisjointCopyN(m->getIndices(),       nelems_, hcol_);
  CoinDisjointCopyN(m->getElements(),      nelems_, rowels_);
  start = m->getVectorStarts();
  length = m->getVectorLengths();
  const int * column = m->getIndices();
  element = m->getElements();
  // out zeros
  int irow;
  nel=0;
  mrstrt_[0]=0;
  for (irow=0;irow<nrows_;irow++) {
    int j;
    for (j=start[irow];j<start[irow]+length[irow];j++) {
      if (fabs(element[j])>ZTOLDP) {
        hcol_[nel]=column[j];
        rowels_[nel++]=element[j];
      }
    }
    hinrow_[irow]=nel-mrstrt_[irow];
    mrstrt_[irow+1]=nel;
  }
  nelems_=nel;

  delete m;
  {
    int i;
    for (i=0;i<ncols_;i++) {
      if (si->isInteger(i))  
	integerType_[i] = 1;
      else
	integerType_[i] = 0;
    }
  }

  // Set up prohibited bits if needed
  if (nonLinearValue) {
    anyProhibited_ = true;
    for (icol=0;icol<ncols_;icol++) {
      int j;
      bool nonLinearColumn = false;
      if (cost_[icol]==nonLinearValue)
	nonLinearColumn=true;
      for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
	if (colels_[j]==nonLinearValue) {
	  nonLinearColumn=true;
	  setRowProhibited(hrow_[j]);
	}
      }
      if (nonLinearColumn)
	setColProhibited(icol);
    }
  } else if (prohibited) {
    anyProhibited_ = true;
    for (icol=0;icol<ncols_;icol++) {
      if (prohibited[icol])
	setColProhibited(icol);
    }
  } else {
    anyProhibited_ = false;
  }
  // Any rows special?
  if (rowProhibited) {
    anyProhibited_ = true;
    for (int irow=0;irow<nrows_;irow++) {
      if (rowProhibited[irow])
	setRowProhibited(irow);
    }
  }
  // Go to minimization
  if (maxmin<0.0) {
    for (int i=0;i<ncols_;i++)
      cost_[i]=-cost_[i];
    maxmin_=1.0;
  }
  if (doStatus) {
    // allow for status and solution
    sol_ = new double[ncols_];
    const double *presol ;
    presol = si->getColSolution() ;
    memcpy(sol_,presol,ncols_*sizeof(double));;
    acts_ = new double [nrows_];
    memcpy(acts_,si->getRowActivity(),nrows_*sizeof(double));
    CoinWarmStartBasis * basis  = 
    dynamic_cast<CoinWarmStartBasis*>(si->getWarmStart());
    colstat_ = new unsigned char [nrows_+ncols_];
    rowstat_ = colstat_+ncols_;
    // If basis is NULL then put in all slack basis
    if (basis&&basis->getNumStructural()==ncols_) {
      int i;
      for (i=0;i<ncols_;i++) {
	colstat_[i] = basis->getStructStatus(i);
      }
      for (i=0;i<nrows_;i++) {
	rowstat_[i] = basis->getArtifStatus(i);
      }
    } else {
      int i;
      // no basis
      for (i=0;i<ncols_;i++) {
	colstat_[i] = 3;
      }
      for (i=0;i<nrows_;i++) {
	rowstat_[i] = 1;
      }
    }
    delete basis;
  } 

# if 0
  for (i=0; i<nrows; ++i)
    printf("NR: %6d\n", hinrow[i]);
  for (int i=0; i<ncols; ++i)
    printf("NC: %6d\n", hincol[i]);
# endif

/*
  For building against CoinUtils 2.6, this #if 1 need to be changed into an
  #if 0
*/
# if 0
  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
# else
  presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
  presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
# endif

  // this allows last col/row to expand up to bufsize-1 (22);
  // this must come after the calls to presolve_prefix
  mcstrt_[ncols_] = bufsize-1;
  mrstrt_[nrows_] = bufsize-1;
  // Allocate useful arrays
  initializeStuff();

# if PRESOLVE_CONSISTENCY > 0
  presolve_consistent(this) ;
# endif
}

// avoid compiler warnings about unused variables
#if PRESOLVE_SUMMARY > 0
void CoinPresolveMatrix::update_model(OsiSolverInterface * si,
				      int nrows0, int ncols0,
				      CoinBigIndex nelems0)
#else
void CoinPresolveMatrix::update_model(OsiSolverInterface * si,
				      int /*nrows0*/, int /*ncols0*/,
				      CoinBigIndex /*nelems0*/)
#endif
{
  int nels=0;
  int i;
  if (si->getObjSense() < 0.0) {
    for (int i=0;i<ncols_;i++)
      cost_[i]=-cost_[i];
    dobias_=-dobias_;
    maxmin_ = -1.0;
  }
  for ( i=0; i<ncols_; i++) 
    nels += hincol_[i];
  CoinPackedMatrix m(true,nrows_,ncols_,nels, colels_, hrow_,mcstrt_,hincol_);
  si->loadProblem(m, clo_, cup_, cost_, rlo_, rup_);

  for ( i=0; i<ncols_; i++) {
    if (integerType_[i])
      si->setInteger(i);
    else
      si->setContinuous(i);
  }
  si->setDblParam(OsiObjOffset,originalOffset_-dobias_);

# if PRESOLVE_SUMMARY > 0
  std::cout
    << "New ncol/nrow/nels: "
    << ncols_ << "(-" << ncols0-ncols_ << ") "
    << nrows_ << "(-" << nrows0-nrows_ << ") "
    << si->getNumElements() << "(-" << nelems0-si->getNumElements() << ") "
    << std::endl ;
# endif
}











////////////////  POSTSOLVE

CoinPostsolveMatrix::CoinPostsolveMatrix(OsiSolverInterface*  si,
				       int ncols0_in,
				       int nrows0_in,
				       CoinBigIndex nelems0,
				   
				       double maxmin,
				       // end prepost members

				       double *sol_in,
				       double *acts_in,

				       unsigned char *colstat_in,
				       unsigned char *rowstat_in) :

  CoinPrePostsolveMatrix(si, ncols0_in, nrows0_in, nelems0),
/*
  Used only to mark processed columns and rows so that debugging routines know
  what to check.
*/
# if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
  cdone_(new char[ncols0_in]),
  rdone_(new char[nrows0_in])
# else
  cdone_(0),
  rdone_(0)
# endif

{
/*
  The CoinPrePostsolveMatrix constructor will set bulk0_ to bulkRatio_*nelems0.
  By default, bulkRatio_ is 2. This is certainly larger than absolutely
  necessary, but good for efficiency (minimises the need to compress the bulk
  store). The main storage arrays for the threaded column-major representation
  (hrow_, colels_, link_) should be allocated to this size.
*/
  free_list_ = 0 ;
  maxlink_ = bulk0_ ;
  link_ = new int[maxlink_] ;

  nrows_ = si->getNumRows() ;
  ncols_ = si->getNumCols() ;

  sol_=sol_in;
  rowduals_=NULL;
  acts_=acts_in;

  rcosts_=NULL;
  colstat_=colstat_in;
  rowstat_=rowstat_in;

  // this is the *reduced* model, which is probably smaller
  int ncols1 = ncols_ ;
  int nrows1 = nrows_ ;

  const CoinPackedMatrix * m = si->getMatrixByCol();
#if 0
  if (! isGapFree(*m)) {
    CoinPresolveAction::throwCoinError("Matrix not gap free",
				      "CoinPostsolveMatrix");
  }
#endif
  const CoinBigIndex nelemsr = m->getNumElements();

  if (isGapFree(*m)) {
    CoinDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
    mcstrt_[ncols_] = nelems0;	// points to end of bulk store
    CoinDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
    CoinDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
    CoinDisjointCopyN(m->getElements(),     nelemsr, colels_);
  }
  else
  {
    CoinPackedMatrix* mm = new CoinPackedMatrix(*m);
    if( mm->hasGaps())
      mm->removeGaps();
    assert(nelemsr == mm->getNumElements());
    CoinDisjointCopyN(mm->getVectorStarts(), ncols1, mcstrt_);
    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
    mcstrt_[ncols_] = nelems0;  // points to end of bulk store
    CoinDisjointCopyN(mm->getVectorLengths(),ncols1,  hincol_);
    CoinDisjointCopyN(mm->getIndices(),      nelemsr, hrow_);
    CoinDisjointCopyN(mm->getElements(),     nelemsr, colels_);
  }

# if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
  memset(cdone_, -1, ncols0_);
  memset(rdone_, -1, nrows0_);
# endif

  rowduals_ = new double[nrows0_];
  CoinDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
  rcosts_ = new double[ncols0_];
  CoinDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);

#if PRESOLVE_DEBUG > 0
  // check accuracy of reduced costs (rcosts_ is recalculated reduced costs)
  si->getMatrixByCol()->transposeTimes(rowduals_,rcosts_) ;
  const double *obj = si->getObjCoefficients() ;
  const double *dj = si->getReducedCost() ;
  {
    int i;
    for (i=0;i<ncols1;i++) {
      double newDj = obj[i]-rcosts_[i];
      rcosts_[i]=newDj;
      assert (fabs(newDj-dj[i])<1.0e-1);
    }
  }
  // check reduced costs are 0 for basic variables
  {
    int i;
    for (i=0;i<ncols1;i++)
      if (columnIsBasic(i))
	assert (fabs(rcosts_[i])<1.0e-5);
    for (i=0;i<nrows1;i++)
      if (rowIsBasic(i))
	assert (fabs(rowduals_[i])<1.0e-5);
  }
#endif
/*
  CoinPresolve may, once, have handled both minimisation and maximisation,
  but hard-wired minimisation has crept in.
*/
  if (maxmin<0.0) {
    for (int i = 0 ; i < nrows1 ; i++)
      rowduals_[i] = -rowduals_[i] ;
    for (int j = 0 ; j < ncols1 ; j++) {
      rcosts_[j] = -rcosts_[j] ;
    }
  }

/*
  CoinPresolve requires both column solution and row activity for correct
  operation.
*/
  CoinDisjointCopyN(si->getColSolution(), ncols1, sol_);
  CoinDisjointCopyN(si->getRowActivity(), nrows1, acts_) ;
  si->setDblParam(OsiObjOffset,originalOffset_);

  for (int j=0; j<ncols1; j++) {
    CoinBigIndex kcs = mcstrt_[j];
    CoinBigIndex kce = kcs + hincol_[j];
    for (CoinBigIndex k=kcs; k<kce; ++k) {
      link_[k] = k+1;
    }
    if (kce>0)
      link_[kce-1] = NO_LINK ;
  }
  if (maxlink_>0) {
    int ml = maxlink_;
    for (CoinBigIndex k=nelemsr; k<ml; ++k)
      link_[k] = k+1;
    link_[ml-1] = NO_LINK;
  }
  free_list_ = nelemsr;

# if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
/*
  These are used to track the action of postsolve transforms during debugging.
*/
  CoinFillN(cdone_,ncols1,PRESENT_IN_REDUCED) ;
  CoinZeroN(cdone_+ncols1,ncols0_in-ncols1) ;
  CoinFillN(rdone_,nrows1,PRESENT_IN_REDUCED) ;
  CoinZeroN(rdone_+nrows1,nrows0_in-nrows1) ;
# endif
}