// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

//#define	PRESOLVE_CONSISTENCY	1
//#define	PRESOLVE_DEBUG	1

#include <stdio.h>

#include <assert.h>
#include <iostream>

#include "CoinHelperFunctions.hpp"

#include "CoinPackedMatrix.hpp"
#include "ClpPackedMatrix.hpp"
#include "ClpSimplex.hpp"
#ifndef SLIM_CLP
#include "ClpQuadraticObjective.hpp"
#endif

#include "ClpPresolve.hpp"
#include "CoinPresolveMatrix.hpp"

#include "CoinPresolveEmpty.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolvePsdebug.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"



ClpPresolve::ClpPresolve() :
  originalModel_(NULL),
  presolvedModel_(NULL),
  nonLinearValue_(0.0),
  originalColumn_(NULL),
  originalRow_(NULL),
  rowObjective_(NULL),
  paction_(0),   
  ncols_(0),
  nrows_(0),
  nelems_(0),
  numberPasses_(5),
  substitution_(3),
#ifndef CLP_NO_STD
  saveFile_(""),
#endif
  presolveActions_(0)
{
}

ClpPresolve::~ClpPresolve()
{
  destroyPresolve();
}
// Gets rid of presolve actions (e.g.when infeasible)
void 
ClpPresolve::destroyPresolve()
{
  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;
  delete [] rowObjective_;
  rowObjective_=NULL;
}

/* This version of presolve returns a pointer to a new presolved 
   model.  NULL if infeasible
*/
ClpSimplex * 
ClpPresolve::presolvedModel(ClpSimplex & si,
                            double feasibilityTolerance,
                            bool keepIntegers,
                            int numberPasses,
                            bool dropNames,
                            bool doRowObjective)
{
  // Check matrix
  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
					  1.0e20))
    return NULL;
  else
    return gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,dropNames,
                                doRowObjective);
}
#ifndef CLP_NO_STD
/* This version of presolve updates
   model and saves original data to file.  Returns non-zero if infeasible
*/
int
ClpPresolve::presolvedModelToFile(ClpSimplex &si,std::string fileName,
                                  double feasibilityTolerance,
                                  bool keepIntegers,
                                  int numberPasses,
                                  bool doRowObjective)
{
  // Check matrix
  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
					  1.0e20))
    return 2;
  saveFile_=fileName;
  si.saveModel(saveFile_.c_str());
  ClpSimplex * model = gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,true,
                                            doRowObjective);
  if (model==&si) {
    return 0;
  } else {
    si.restoreModel(saveFile_.c_str());
    remove(saveFile_.c_str());
    return 1;
  }
}
#endif
// Return pointer to presolved model
ClpSimplex * 
ClpPresolve::model() const
{
  return presolvedModel_;
}
// Return pointer to original model
ClpSimplex * 
ClpPresolve::originalModel() const
{
  return originalModel_;
}
void 
ClpPresolve::postsolve(bool updateStatus)
{
  // Return at once if no presolved model
  if (!presolvedModel_)
    return;
  // Messages
  CoinMessages messages = originalModel_->coinMessages();
  if (!presolvedModel_->isProvenOptimal()) {
    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
					     messages)
					       <<CoinMessageEol;
  }

  // this is the size of the original problem
  const int ncols0  = ncols_;
  const int nrows0  = nrows_;
  const CoinBigIndex nelems0 = nelems_;
  
  // this is the reduced problem
  int ncols = presolvedModel_->getNumCols();
  int nrows = presolvedModel_->getNumRows();

  double * acts=NULL;
  double * sol =NULL;
  unsigned char * rowstat=NULL;
  unsigned char * colstat = NULL;
#ifndef CLP_NO_STD
  if (saveFile_=="") {
#endif
    // reality check
    assert(ncols0==originalModel_->getNumCols());
    assert(nrows0==originalModel_->getNumRows());
    acts = originalModel_->primalRowSolution();
    sol  = originalModel_->primalColumnSolution();
    if (updateStatus) {
      // postsolve does not know about fixed
      int i;
      for (i=0;i<nrows+ncols;i++) {
        if (presolvedModel_->getColumnStatus(i)==ClpSimplex::isFixed)
          presolvedModel_->setColumnStatus(i,ClpSimplex::atLowerBound);
      }
      unsigned char *status = originalModel_->statusArray();
      if (!status) {
        originalModel_->createStatus();
        status = originalModel_->statusArray();
      }
      rowstat = status + ncols0;
      colstat = status;
      memcpy(colstat, presolvedModel_->statusArray(), ncols);
      memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
    }
#ifndef CLP_NO_STD
  } else {
    // from file
    acts = new double[nrows0];
    sol  = new double[ncols0];
    CoinZeroN(acts,nrows0);
    CoinZeroN(sol,ncols0);
    if (updateStatus) {
      unsigned char *status = new unsigned char [nrows0+ncols0];
      rowstat = status + ncols0;
      colstat = status;
      memcpy(colstat, presolvedModel_->statusArray(), ncols);
      memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
    }
  }
#endif

  // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
  // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
  // cause duplicate free. In case where saveFile == "", as best I can see
  // arrays are owned by originalModel_. fix is to
  // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
  CoinPostsolveMatrix prob(presolvedModel_,
		       ncols0,
		       nrows0,
		       nelems0,
		       presolvedModel_->getObjSense(),
		       // end prepost
		       
		       sol, acts,
		       colstat, rowstat);
    
  postsolve(prob);

#ifndef CLP_NO_STD
  if (saveFile_!="") {
    // From file
    assert (originalModel_==presolvedModel_);
    originalModel_->restoreModel(saveFile_.c_str());
    remove(saveFile_.c_str());
    memcpy(originalModel_->primalRowSolution(),acts,nrows0*sizeof(double));
    // delete [] acts;
    memcpy(originalModel_->primalColumnSolution(),sol,ncols0*sizeof(double));
    // delete [] sol;
    if (updateStatus) {
      memcpy(originalModel_->statusArray(),colstat,nrows0+ncols0);
      // delete [] colstat;
    }
  } else {
#endif
    prob.sol_ = 0 ;
    prob.acts_ = 0 ;
    prob.colstat_ = 0 ;
#ifndef CLP_NO_STD
  }
#endif
  // put back duals
  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
	 nrows_*sizeof(double));
  double maxmin = originalModel_->getObjSense();
  if (maxmin<0.0) {
    // swap signs
    int i;
    double * pi = originalModel_->dualRowSolution();
    for (i=0;i<nrows_;i++)
      pi[i] = -pi[i];
  }
  // Now check solution
  double offset;
  memcpy(originalModel_->dualColumnSolution(),
	 originalModel_->objectiveAsObject()->gradient(originalModel_,
						       originalModel_->primalColumnSolution(),
						       offset,true),ncols_*sizeof(double));
  originalModel_->transposeTimes(-1.0,
				 originalModel_->dualRowSolution(),
				 originalModel_->dualColumnSolution());
  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
			originalModel_->primalRowSolution());
  originalModel_->checkSolutionInternal();
  // Messages
  presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
					    messages)
					      <<originalModel_->objectiveValue()
					      <<originalModel_->sumDualInfeasibilities()
					      <<originalModel_->numberDualInfeasibilities()
					      <<originalModel_->sumPrimalInfeasibilities()
					      <<originalModel_->numberPrimalInfeasibilities()
					       <<CoinMessageEol;
  
  //originalModel_->objectiveValue_=objectiveValue_;
  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
  if (!presolvedModel_->status()) {
    if (!originalModel_->numberDualInfeasibilities()&&
	!originalModel_->numberPrimalInfeasibilities()) {
      originalModel_->setProblemStatus( 0);
    } else {
      originalModel_->setProblemStatus( -1);
      // Say not optimal after presolve
      originalModel_->setSecondaryStatus(7);
      presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
					    messages)
					      <<CoinMessageEol;
    }
  } else {
    originalModel_->setProblemStatus( presolvedModel_->status());
  }
#ifndef CLP_NO_STD
  if (saveFile_!="") 
    presolvedModel_=NULL;
#endif
}

// return pointer to original columns
const int * 
ClpPresolve::originalColumns() const
{
  return originalColumn_;
}
// return pointer to original rows
const int * 
ClpPresolve::originalRows() const
{
  return originalRow_;
}
// Set pointer to original model
void 
ClpPresolve::setOriginalModel(ClpSimplex * 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 || PRESOLVE_SUMMARY
  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
//#define PRESOLVE_DEBUG 1
#if PRESOLVE_DEBUG
void check_sol(CoinPresolveMatrix *prob,double tol)
{
  double *colels	= prob->colels_;
  int *hrow		= prob->hrow_;
  int *mcstrt		= prob->mcstrt_;
  int *hincol		= prob->hincol_;
  int *hinrow		= prob->hinrow_;
  int ncols		= prob->ncols_;


  double * csol = prob->sol_;
  double * acts = prob->acts_;
  double * clo = prob->clo_;
  double * cup = prob->cup_;
  int nrows = prob->nrows_;
  double * rlo = prob->rlo_;
  double * rup = prob->rup_;

  int colx;

  double * rsol = new double[nrows];
  memset(rsol,0,nrows*sizeof(double));

  for (colx = 0; colx < ncols; ++colx) {
    if (1) {
      CoinBigIndex k = mcstrt[colx];
      int nx = hincol[colx];
      double solutionValue = csol[colx];
      for (int i=0; i<nx; ++i) {
	int row = hrow[k];
	double coeff = colels[k];
	k++;
	rsol[row] += solutionValue*coeff;
      }
      if (csol[colx]<clo[colx]-tol) {
	printf("low CSOL:  %d  - %g %g %g\n",
		   colx, clo[colx], csol[colx], cup[colx]);
      } else if (csol[colx]>cup[colx]+tol) {
	printf("high CSOL:  %d  - %g %g %g\n",
		   colx, clo[colx], csol[colx], cup[colx]);
      } 
    }
  }
  int rowx;
  for (rowx = 0; rowx < nrows; ++rowx) {
    if (hinrow[rowx]) {
      if (fabs(rsol[rowx]-acts[rowx])>tol)
	printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
		   rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
      if (rsol[rowx]<rlo[rowx]-tol) {
	printf("low RSOL:  %d - %g %g %g\n",
		   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
      } else if (rsol[rowx]>rup[rowx]+tol ) {
	printf("high RSOL:  %d - %g %g %g\n",
		   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
      } 
    }
  }
  delete [] rsol;
}
#endif
// This is the presolve loop.
// It is a separate virtual function so that it can be easily
// customized by subclassing CoinPresolve.
const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
{
  // Messages
  CoinMessages messages = CoinMessage(prob->messages().language());
  paction_ = 0;

  prob->status_=0; // say feasible
  paction_ = make_fixed(prob, paction_);
  // if integers then switch off dual stuff
  // later just do individually
  bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
  // but allow in some cases
  if ((presolveActions_&512)!=0)
    doDualStuff=true;
  if (prob->anyProhibited()) 
    doDualStuff=false;
  if (!doDual())
    doDualStuff=false;
#if	PRESOLVE_CONSISTENCY
//  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
    presolve_links_ok(prob,false,true) ;
#endif

  if (!prob->status_) {
    bool slackSingleton = doSingletonColumn();
    slackSingleton = false;
    const bool slackd = doSingleton();
    const bool doubleton = doDoubleton();
    const bool tripleton = doTripleton();
    const bool forcing = doForcing();
    const bool ifree = doImpliedFree();
    const bool zerocost = doTighten();
    const bool dupcol = doDupcol();
    const bool duprow = doDuprow();
    const bool dual = doDualStuff;
    
    // some things are expensive so just do once (normally)

    int i;
    // say look at all
    if (!prob->anyProhibited()) {
      for (i=0;i<nrows_;i++) 
	prob->rowsToDo_[i]=i;
      prob->numberRowsToDo_=nrows_;
      for (i=0;i<ncols_;i++) 
	prob->colsToDo_[i]=i;
      prob->numberColsToDo_=ncols_;
    } else {
      // some stuff must be left alone
      prob->numberRowsToDo_=0;
      for (i=0;i<nrows_;i++) 
	if (!prob->rowProhibited(i))
	    prob->rowsToDo_[prob->numberRowsToDo_++]=i;
      prob->numberColsToDo_=0;
      for (i=0;i<ncols_;i++) 
	if (!prob->colProhibited(i))
	    prob->colsToDo_[prob->numberColsToDo_++]=i;
    }


    int iLoop;
#if	PRESOLVE_DEBUG
    check_sol(prob,1.0e0);
#endif
    if (dupcol) {
      //paction_ = dupcol_action::presolve(prob, paction_);
    }

    if ((presolveActions_&16384)!=0)
      prob->setPresolveOptions(prob->presolveOptions()|16384);
    // Check number rows dropped
    int lastDropped=0;
    prob->pass_=0;
    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
    // See if we want statistics
    if ((presolveActions_&0x80000000)!=0)
      printf("Starting major pass %d after %g seconds\n",iLoop+1,CoinCpuTime()-prob->startTime_);
#ifdef PRESOLVE_SUMMARY
      printf("Starting major pass %d\n",iLoop+1);
#endif
      const CoinPresolveAction * const paction0 = paction_;
      // look for substitutions with no fill
      int fill_level=2;
      //fill_level=10;
      //printf("** fill_level == 10 !\n");
      int whichPass=0;
      while (1) {
	whichPass++;
	prob->pass_++;
	const CoinPresolveAction * const paction1 = paction_;

	if (slackd) {
	  bool notFinished = true;
	  while (notFinished) 
	    paction_ = slack_doubleton_action::presolve(prob, paction_,
							notFinished);
	  if (prob->status_)
	    break;
	}
	if (dual&&whichPass==1) {
	  // this can also make E rows so do one bit here
	  paction_ = remove_dual_action::presolve(prob, paction_);
	  if (prob->status_)
	    break;
	}

	if (doubleton) {
	  paction_ = doubleton_action::presolve(prob, paction_);
	  if (prob->status_)
	    break;
	}
	if (tripleton) {
	  paction_ = tripleton_action::presolve(prob, paction_);
	  if (prob->status_)
	    break;
	}

	if (zerocost) {
	  paction_ = do_tighten_action::presolve(prob, paction_);
	  if (prob->status_)
	    break;
	}

	if (forcing) {
	  paction_ = forcing_constraint_action::presolve(prob, paction_);
	  if (prob->status_)
	    break;
	}

	if (ifree&&(whichPass%5)==1) {
	paction_ = implied_free_action::presolve(prob, paction_,fill_level);
	if (prob->status_)
	  break;
	}

#if	PRESOLVE_DEBUG
	check_sol(prob,1.0e0);
#endif

#if	PRESOLVE_CONSISTENCY
//	presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
//			  prob->nrows_);
	presolve_links_ok(prob,false,true) ;
#endif

//#if	PRESOLVE_DEBUG
//	presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
//			  prob->ncols_);
//#endif
//#if	PRESOLVE_CONSISTENCY
//	prob->consistent();
//#endif
#if	PRESOLVE_CONSISTENCY
	presolve_no_zeros(prob,true,false) ;
	presolve_consistent(prob,true) ;
#endif

	  
	// set up for next pass
	// later do faster if many changes i.e. memset and memcpy
	prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
	int kcheck;
	bool found=false;
	kcheck=-1;
	for (i=0;i<prob->numberNextRowsToDo_;i++) {
	  int index = prob->nextRowsToDo_[i];
	  prob->unsetRowChanged(index);
	  prob->rowsToDo_[i] = index;
	  if (index==kcheck) {
	    printf("row %d on list after pass %d\n",kcheck,
		   whichPass);
	    found=true;
	  }
	}
	if (!found&&kcheck>=0)
	  prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
	prob->numberNextRowsToDo_=0;
	prob->numberColsToDo_ = prob->numberNextColsToDo_;
	kcheck=-1;
	found=false;
	for (i=0;i<prob->numberNextColsToDo_;i++) {
	  int index = prob->nextColsToDo_[i];
	  prob->unsetColChanged(index);
	  prob->colsToDo_[i] = index;
	  if (index==kcheck) {
	    printf("col %d on list after pass %d\n",kcheck,
		   whichPass);
	    found=true;
	  }
	}
	if (!found&&kcheck>=0)
	  prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
	prob->numberNextColsToDo_=0;
	if (paction_ == paction1&&fill_level>0)
	  break;
      }
      // say look at all
      int i;
      if (!prob->anyProhibited()) {
	for (i=0;i<nrows_;i++) 
	  prob->rowsToDo_[i]=i;
	prob->numberRowsToDo_=nrows_;
	for (i=0;i<ncols_;i++) 
	  prob->colsToDo_[i]=i;
	prob->numberColsToDo_=ncols_;
      } else {
	// some stuff must be left alone
	prob->numberRowsToDo_=0;
	for (i=0;i<nrows_;i++) 
	  if (!prob->rowProhibited(i))
	    prob->rowsToDo_[prob->numberRowsToDo_++]=i;
	prob->numberColsToDo_=0;
	for (i=0;i<ncols_;i++) 
	  if (!prob->colProhibited(i))
	    prob->colsToDo_[prob->numberColsToDo_++]=i;
      }
      // now expensive things
      // this caused world.mps to run into numerical difficulties
#ifdef PRESOLVE_SUMMARY
      printf("Starting expensive\n");
#endif

      if (dual) {
	int itry;
	for (itry=0;itry<5;itry++) {
	  paction_ = remove_dual_action::presolve(prob, paction_);
	  if (prob->status_)
	    break;
	  const CoinPresolveAction * const paction2 = paction_;
	  if (ifree) {
	    //int fill_level=0; // switches off substitution
	    paction_ = implied_free_action::presolve(prob, paction_,fill_level);
	    if (prob->status_)
	      break;
	  }
	  if (paction_ == paction2)
	    break;
	}
      } else if (ifree) {
	// just free
	int fill_level=0; // switches off substitution
	paction_ = implied_free_action::presolve(prob, paction_,fill_level);
	if (prob->status_)
	  break;
      }
#if	PRESOLVE_DEBUG
      check_sol(prob,1.0e0);
#endif
      if (dupcol) {
        // maybe allow integer columns to be checked
        if ((presolveActions_&512)!=0)
          prob->setPresolveOptions(prob->presolveOptions()|1);
	paction_ = dupcol_action::presolve(prob, paction_);
	if (prob->status_)
	  break;
      }
#if	PRESOLVE_DEBUG
	check_sol(prob,1.0e0);
#endif
      
      if (duprow) {
	paction_ = duprow_action::presolve(prob, paction_);
	if (prob->status_)
	  break;
      }
#if	PRESOLVE_DEBUG
      check_sol(prob,1.0e0);
#endif
      bool stopLoop=false;
      {
	int * hinrow = prob->hinrow_;
	int numberDropped=0;
	for (i=0;i<nrows_;i++) 
	  if (!hinrow[i])
	    numberDropped++;
	
	prob->messageHandler()->message(COIN_PRESOLVE_PASS,
					messages)
					  <<numberDropped<<iLoop+1
					  <<CoinMessageEol;
	//printf("%d rows dropped after pass %d\n",numberDropped,
	//     iLoop+1);
	if (numberDropped==lastDropped)
	  stopLoop=true;
	else
	  lastDropped = numberDropped;
      }
      // Do this here as not very loopy
      if (slackSingleton) {
        // On most passes do not touch costed slacks
        if (paction_ != paction0&&!stopLoop) {
          paction_ = slack_singleton_action::presolve(prob, paction_,NULL);
        } else {
          // do costed if Clp (at end as ruins rest of presolve)
          paction_ = slack_singleton_action::presolve(prob, paction_,rowObjective_);
          stopLoop=true;
        }
      }
#if	PRESOLVE_DEBUG
      check_sol(prob,1.0e0);
#endif
      if (paction_ == paction0||stopLoop)
	break;
	  
    }
  }
  if (!prob->status_) {
    paction_ = drop_zero_coefficients(prob, paction_);
#if	PRESOLVE_DEBUG
	check_sol(prob,1.0e0);
#endif

    paction_ = drop_empty_cols_action::presolve(prob, paction_);
    paction_ = drop_empty_rows_action::presolve(prob, paction_);
#if	PRESOLVE_DEBUG
	check_sol(prob,1.0e0);
#endif
  }
  
  if (prob->status_) {
    if (prob->status_==1) 
	  prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
					     messages)
					       <<prob->feasibilityTolerance_
					       <<CoinMessageEol;
    else if (prob->status_==2)
	  prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
					     messages)
					       <<CoinMessageEol;
    else
	  prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
					     messages)
					       <<CoinMessageEol;
    // get rid of data
    destroyPresolve();
  }
  return (paction_);
}

void check_djs(CoinPostsolveMatrix *prob);


// We could have implemented this by having each postsolve routine
// directly call the next one, but this may make it easier to add debugging checks.
void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
{
  {
    // Check activities
    double *colels	= prob.colels_;
    int *hrow		= prob.hrow_;
    CoinBigIndex *mcstrt		= prob.mcstrt_;
    int *hincol		= prob.hincol_;
    int *link		= prob.link_;
    int ncols		= prob.ncols_;

    char *cdone	= prob.cdone_;

    double * csol = prob.sol_;
    int nrows = prob.nrows_;

    int colx;

    double * rsol = prob.acts_;
    memset(rsol,0,nrows*sizeof(double));

    for (colx = 0; colx < ncols; ++colx) {
      if (cdone[colx]) {
	CoinBigIndex k = mcstrt[colx];
	int nx = hincol[colx];
	double solutionValue = csol[colx];
	for (int i=0; i<nx; ++i) {
	  int row = hrow[k];
	  double coeff = colels[k];
	  k = link[k];
	  rsol[row] += solutionValue*coeff;
	}
      }
    }
  }
  const CoinPresolveAction *paction = paction_;
  //#define PRESOLVE_DEBUG 1
#if	PRESOLVE_DEBUG
  // Check only works after first one
  int checkit=-1;
#endif
  
  while (paction) {
#if	PRESOLVE_DEBUG
    printf("POSTSOLVING %s\n", paction->name());
#endif

    paction->postsolve(&prob);
    
#if	PRESOLVE_DEBUG
  {
    int nr=0;
    int i;
    for (i=0;i<prob.nrows_;i++) {
      if ((prob.rowstat_[i]&7)==1)
        nr++;
    }
    int nc=0;
    for (i=0;i<prob.ncols_;i++) {
      if ((prob.colstat_[i]&7)==1)
        nc++;
    }
    printf("%d rows (%d basic), %d cols (%d basic)\n",prob.nrows_,nr,prob.ncols_,nc);
  }
    checkit++;
    if (prob.colstat_&&checkit>0) {
      presolve_check_nbasic(&prob) ;
      presolve_check_sol(&prob,2,2,1) ;
    }
#endif
    paction = paction->next;
#if	PRESOLVE_DEBUG
//  check_djs(&prob);
    if (checkit>0)
      presolve_check_reduced_costs(&prob) ;
#endif
  }    
#if	PRESOLVE_DEBUG
  if (prob.colstat_) {
    presolve_check_nbasic(&prob) ;
    presolve_check_sol(&prob,2,2,1) ;
  }
#endif
#undef PRESOLVE_DEBUG
  
#if	0 && PRESOLVE_DEBUG
  for (i=0; i<ncols0; i++) {
    if (!cdone[i]) {
      printf("!cdone[%d]\n", i);
      abort();
    }
  }
  
  for (i=0; i<nrows0; i++) {
    if (!rdone[i]) {
      printf("!rdone[%d]\n", i);
      abort();
    }
  }
  
  
  for (i=0; i<ncols0; i++) {
    if (sol[i] < -1e10 || sol[i] > 1e10)
      printf("!!!%d %g\n", i, sol[i]);
    
  }
  
  
#endif
  
#if	0 && PRESOLVE_DEBUG
  // debug check:  make sure we ended up with same original matrix
  {
    int identical = 1;
    
    for (int i=0; i<ncols0; i++) {
      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
      CoinBigIndex kcs0 = &prob->mcstrt0[i];
      CoinBigIndex kcs = mcstrt[i];
      int n = hincol[i];
      for (int k=0; k<n; k++) {
	CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);

	if (k1 == kcs+n) {
	  printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
	  abort();
	}

	if (colels[k1] != &prob->dels0[kcs0+k])
	  printf("BAD COLEL[%d %d %d]:  %g\n",
		 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);

	if (kcs0+k != k1)
	  identical=0;
      }
    }
    printf("identical? %d\n", identical);
  }
#endif
}








static inline double getTolerance(const ClpSimplex  *si, ClpDblParam 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 ClpSimplex * si,
					     int ncols_in,
					     int nrows_in,
					     CoinBigIndex nelems_in,
                                               double bulkRatio) 
  : ncols_(si->getNumCols()),
    nrows_(si->getNumRows()),
    nelems_(si->getNumElements()),
    ncols0_(ncols_in),
    nrows0_(nrows_in),
    bulkRatio_(bulkRatio),
    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, ClpPrimalTolerance)),
    ztoldj_(getTolerance(si, ClpDualTolerance)),
    maxmin_(si->getObjSense()),
    sol_(NULL),
    rowduals_(NULL),
    acts_(NULL),
    rcosts_(NULL),
    colstat_(NULL),
    rowstat_(NULL),
    handler_(NULL),
    defaultHandler_(false),
    messages_()

{
  bulk0_ = (CoinBigIndex) (bulkRatio_*nelems_in);
  hrow_  = new int   [bulk0_];
  colels_ = new double[bulk0_];
  si->getDblParam(ClpObjOffset,originalOffset_);
  int ncols = si->getNumCols();
  int nrows = si->getNumRows();

  setMessageHandler(si->messageHandler()) ;

  ClpDisjointCopyN(si->getColLower(), ncols, clo_);
  ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
  //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
  double offset;
  ClpDisjointCopyN(si->objectiveAsObject()->gradient(si,si->getColSolution(),offset,true), ncols, cost_);
  ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
  ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
  int i;
  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 = matrix.getSizeVectorLengths() - 1;
  // Quick check
  if (matrix.getNumElements()==start[i]) {
    return true;
  } else {
    for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
      if (start[i+1] - start[i] != length[i])
        break;
    }
    return (! (i >= 0));
  }
}
#if	PRESOLVE_DEBUG
static void matrix_bounds_ok(const double *lo, const double *up, int n)
{
  int i;
  for (i=0; i<n; i++) {
    PRESOLVEASSERT(lo[i] <= up[i]);
    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
  }
}
#endif
CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
				     double maxmin,
				     // end prepost members

				     ClpSimplex * si,

				     // rowrep
				     int nrows_in,
				     CoinBigIndex nelems_in,
                                       bool doStatus,
                                       double nonLinearValue,
                                       double bulkRatio) :

  CoinPrePostsolveMatrix(si,
			ncols0_in, nrows_in, nelems_in,bulkRatio),
  clink_(new presolvehlink[ncols0_in+1]),
  rlink_(new presolvehlink[nrows_in+1]),

  dobias_(0.0),


  // temporary init
  integerType_(new unsigned char[ncols0_in]),
  tuning_(false),
  startTime_(0.0),
  feasibilityTolerance_(0.0),
  status_(-1),
  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)

{
  const int bufsize = bulk0_;

  nrows_ = si->getNumRows() ;

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

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

  const CoinBigIndex * start = m->getVectorStarts();
  const int * length = m->getVectorLengths();
  const int * row = m->getIndices();
  const double * element = m->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++) {
      hrow_[nel]=row[j];
      colels_[nel++]=element[j];
    }
    mcstrt_[icol+1]=nel;
  }
  assert(mcstrt_[ncols_] == nelems_);
  ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);

  // same thing for row rep
  CoinPackedMatrix * mRow = new CoinPackedMatrix();
  mRow->reverseOrderedCopyOf(*m);
  mRow->removeGaps();
  mRow->setExtraGap(0.0);

  // Now get rid of matrix
  si->createEmptyMatrix();

  double * el = mRow->getMutableElements();
  int * ind = mRow->getMutableIndices();
  CoinBigIndex * strt = mRow->getMutableVectorStarts();
  int * len = mRow->getMutableVectorLengths();
  // Do carefully to save memory
  rowels_ = new double[bulk0_];
  ClpDisjointCopyN(el,      nelems_, rowels_);
  mRow->nullElementArray();
  delete [] el;
  hcol_ = new int[bulk0_];
  ClpDisjointCopyN(ind,       nelems_, hcol_);
  mRow->nullIndexArray();
  delete [] ind;
  mrstrt_ = new CoinBigIndex[nrows_in+1];
  ClpDisjointCopyN(strt,  nrows_,  mrstrt_);
  mRow->nullStartArray();
  mrstrt_[nrows_] = nelems_;
  delete [] strt;
  hinrow_ = new int[nrows_in+1];
  ClpDisjointCopyN(len, nrows_,  hinrow_);

  delete mRow;
  if (si->integerInformation()) {
    memcpy(integerType_,si->integerInformation(),ncols_*sizeof(char));
  } else {
    ClpFillN<unsigned char>(integerType_, ncols_, (unsigned char) 0);
  }

#ifndef SLIM_CLP
#ifndef NO_RTTI
    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
#else
    ClpQuadraticObjective * quadraticObj = NULL;
    if (si->objectiveAsObject()->type()==2)
      quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
#endif
#endif
  // 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);
    }
#ifndef SLIM_CLP
  } else if (quadraticObj) {
    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    //const int * columnQuadratic = quadratic->getIndices();
    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    const int * columnQuadraticLength = quadratic->getVectorLengths();
    //double * quadraticElement = quadratic->getMutableElements();
    int numberColumns = quadratic->getNumCols();
    anyProhibited_ = true;
    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      if (columnQuadraticLength[iColumn]) {
	setColProhibited(iColumn);
	//printf("%d prohib\n",iColumn);
      }
    }
#endif
  } else {
    anyProhibited_ = false;
  }

  if (doStatus) {
    // allow for status and solution
    sol_ = new double[ncols_];
    memcpy(sol_,si->primalColumnSolution(),ncols_*sizeof(double));;
    acts_ = new double [nrows_];
    memcpy(acts_,si->primalRowSolution(),nrows_*sizeof(double));
    if (!si->statusArray())
      si->createStatus();
    colstat_ = new unsigned char [nrows_+ncols_];
    memcpy(colstat_,si->statusArray(),
	   (nrows_+ncols_)*sizeof(unsigned char));
    rowstat_ = colstat_+ncols_;
  }

  // the original model's fields are now unneeded - free them
  
  si->resize(0,0);

#if	PRESOLVE_DEBUG
  matrix_bounds_ok(rlo_, rup_, nrows_);
  matrix_bounds_ok(clo_, cup_, ncols_);
#endif

#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

  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);

  // 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;

#if	PRESOLVE_CONSISTENCY
//consistent(false);
  presolve_consistent(this,false) ;
#endif
}

void CoinPresolveMatrix::update_model(ClpSimplex * si,
				     int nrows0,
				     int ncols0,
				     CoinBigIndex nelems0)
{
  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
		 clo_, cup_, cost_, rlo_, rup_);
  //delete [] si->integerInformation();
  int numberIntegers=0;
  for (int i=0; i<ncols_; i++) {
    if (integerType_[i])
      numberIntegers++;
  }
  if (numberIntegers) 
    si->copyInIntegerInformation((const char *) integerType_);
  else
    si->copyInIntegerInformation(NULL);

#if	PRESOLVE_SUMMARY
  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
	 ncols_, ncols0-ncols_,
	 nrows_, nrows0-nrows_,
	 si->getNumElements(), nelems0-si->getNumElements());
#endif
  si->setDblParam(ClpObjOffset,originalOffset_-dobias_);

}











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

CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  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,2.0),

  free_list_(0),
  // link, free_list, maxlink
  maxlink_(bulk0_),
  link_(new int[/*maxlink*/ bulk0_]),
      
  cdone_(new char[ncols0_]),
  rdone_(new char[nrows0_in])

{
  bulk0_ = 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->matrix();

  const CoinBigIndex nelemsr = m->getNumElements();
  if (m->getNumElements()&&!isGapFree(*m)) {
    // Odd - gaps
    CoinPackedMatrix mm(*m);
    mm.removeGaps();
    mm.setExtraGap(0.0);
    
    ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
    mcstrt_[ncols1] = nelems0;	// ??    (should point to end of bulk store   -- lh --)
    ClpDisjointCopyN(mm.getVectorLengths(),ncols1,  hincol_);
    ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
    ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
  } else {
    // No gaps
    
    ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
    mcstrt_[ncols1] = nelems0;	// ??    (should point to end of bulk store   -- lh --)
    ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
    ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
    ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
  }



#if	0 && PRESOLVE_DEBUG
  presolve_check_costs(model, &colcopy);
#endif

  // This determines the size of the data structure that contains
  // the matrix being postsolved.  Links are taken from the free_list
  // to recreate matrix entries that were presolved away,
  // and links are added to the free_list when entries created during
  // presolve are discarded.  There is never a need to gc this list.
  // Naturally, it should contain
  // exactly nelems0 entries "in use" when postsolving is done,
  // but I don't know whether the matrix could temporarily get
  // larger during postsolving.  Substitution into more than two
  // rows could do that, in principle.  I am being very conservative
  // here by reserving much more than the amount of space I probably need.
  // If this guess is wrong, check_free_list may be called.
  //  int bufsize = 2*nelems0;

  memset(cdone_, -1, ncols0_);
  memset(rdone_, -1, nrows0_);

  rowduals_ = new double[nrows0_];
  ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);

  rcosts_ = new double[ncols0_];
  ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
  if (maxmin<0.0) {
    // change so will look as if minimize
    int i;
    for (i=0;i<nrows1;i++)
      rowduals_[i] = - rowduals_[i];
    for (i=0;i<ncols1;i++) {
      rcosts_[i] = - rcosts_[i];
    }
  }

  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);

  ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
  si->setDblParam(ClpObjOffset,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;
    }
    link_[kce-1] = NO_LINK ;
  }
  {
    int ml = maxlink_;
    for (CoinBigIndex k=nelemsr; k<ml; ++k)
      link_[k] = k+1;
    if (ml)
      link_[ml-1] = NO_LINK;
  }
  free_list_ = nelemsr;
# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
/*
  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
}
/* This is main part of Presolve */
ClpSimplex * 
ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
				  double feasibilityTolerance,
				  bool keepIntegers,
				  int numberPasses,
				  bool dropNames,
                                  bool doRowObjective)
{
  ncols_ = originalModel->getNumCols();
  nrows_ = originalModel->getNumRows();
  nelems_ = originalModel->getNumElements();
  numberPasses_ = numberPasses;

  double maxmin = originalModel->getObjSense();
  originalModel_ = originalModel;
  delete [] originalColumn_;
  originalColumn_ = new int[ncols_];
  delete [] originalRow_;
  originalRow_ = new int[nrows_];
  // and fill in case returns early
  int i;
  for (i=0;i<ncols_;i++)
    originalColumn_[i]=i;
  for (i=0;i<nrows_;i++)
    originalRow_[i]=i;
  delete [] rowObjective_;
  if (doRowObjective) {
    rowObjective_ = new double [nrows_];
    memset(rowObjective_,0,nrows_*sizeof(double));
  } else {
    rowObjective_=NULL;
  }

  // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
  int result = -1;
  
  // User may have deleted - its their responsibility
  presolvedModel_=NULL;
  // Messages
  CoinMessages messages = originalModel->coinMessages();
  // Only go round 100 times even if integer preprocessing
  int totalPasses=100;
  while (result==-1) {

#ifndef CLP_NO_STD
    // make new copy
    if (saveFile_=="") {
#endif
      delete presolvedModel_;
#ifndef CLP_NO_STD
      // So won't get names
      int lengthNames = originalModel->lengthNames();
      originalModel->setLengthNames(0);
#endif
      presolvedModel_ = new ClpSimplex(*originalModel);
#ifndef CLP_NO_STD
      originalModel->setLengthNames(lengthNames);
    } else {
      presolvedModel_=originalModel;
    }
    presolvedModel_->dropNames();
#endif

    // drop integer information if wanted
    if (!keepIntegers)
      presolvedModel_->deleteIntegerInformation();
    totalPasses--;

    double ratio=2.0;
    if (substitution_>3)
      ratio=substitution_;
    else if (substitution_==2)
      ratio=1.5;
    CoinPresolveMatrix prob(ncols_,
			maxmin,
			presolvedModel_,
			nrows_, nelems_,true,nonLinearValue_,ratio);
    prob.setMaximumSubstitutionLevel(substitution_);
    if (doRowObjective) 
      memset(rowObjective_,0,nrows_*sizeof(double));
    // See if we want statistics
    if ((presolveActions_&0x80000000)!=0)
      prob.statistics();
    // make sure row solution correct
    {
      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
    paction_ = presolve(&prob);

    result =0; 

    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 (prob.status_==0&&paction_) {
      // feasible
    
      prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
      // copy status and solution
      memcpy(presolvedModel_->primalColumnSolution(),
	     prob.sol_,prob.ncols_*sizeof(double));
      memcpy(presolvedModel_->primalRowSolution(),
	     prob.acts_,prob.nrows_*sizeof(double));
      memcpy(presolvedModel_->statusArray(),
	     prob.colstat_,prob.ncols_*sizeof(unsigned char));
      memcpy(presolvedModel_->statusArray()+prob.ncols_,
	     prob.rowstat_,prob.nrows_*sizeof(unsigned char));
      delete [] prob.sol_;
      delete [] prob.acts_;
      delete [] prob.colstat_;
      prob.sol_=NULL;
      prob.acts_=NULL;
      prob.colstat_=NULL;
      
      int ncolsNow = presolvedModel_->getNumCols();
      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
#ifndef SLIM_CLP
#ifndef NO_RTTI
      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
#else
      ClpQuadraticObjective * quadraticObj = NULL;
      if (originalModel->objectiveAsObject()->type()==2)
	quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
#endif
      if (quadraticObj) {
	// set up for subset
	char * mark = new char [ncols_];
	memset(mark,0,ncols_);
	CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
	//const int * columnQuadratic = quadratic->getIndices();
	//const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
	const int * columnQuadraticLength = quadratic->getVectorLengths();
	//double * quadraticElement = quadratic->getMutableElements();
	int numberColumns = quadratic->getNumCols();
	ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
								   ncolsNow,
								   originalColumn_);
	// and modify linear and check
	double * linear = newObj->linearObjective();
	memcpy(linear,presolvedModel_->objective(),ncolsNow*sizeof(double));
	int iColumn;
	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (columnQuadraticLength[iColumn])
	    mark[iColumn]=1;
	}
	// and new
	quadratic = newObj->quadraticObjective();
	columnQuadraticLength = quadratic->getVectorLengths();
	int numberColumns2 = quadratic->getNumCols();
	for ( iColumn=0;iColumn<numberColumns2;iColumn++) {
	  if (columnQuadraticLength[iColumn])
	    mark[originalColumn_[iColumn]]=0;
	}
	presolvedModel_->setObjective(newObj);
	delete newObj;
	// final check
	for ( iColumn=0;iColumn<numberColumns;iColumn++) 
	  if (mark[iColumn])
	    printf("Quadratic column %d modified - may be okay\n",iColumn);
	delete [] mark;
      }
#endif
      delete [] prob.originalColumn_;
      prob.originalColumn_=NULL;
      int nrowsNow = presolvedModel_->getNumRows();
      memcpy(originalRow_,prob.originalRow_,nrowsNow*sizeof(int));
      delete [] prob.originalRow_;
      prob.originalRow_=NULL;
#ifndef CLP_NO_STD
      if (!dropNames&&originalModel->lengthNames()) {
	// Redo names
	int iRow;
	std::vector<std::string> rowNames;
	rowNames.reserve(nrowsNow);
	for (iRow=0;iRow<nrowsNow;iRow++) {
	  int kRow = originalRow_[iRow];
	  rowNames.push_back(originalModel->rowName(kRow));
	}
      
	int iColumn;
	std::vector<std::string> columnNames;
	columnNames.reserve(ncolsNow);
	for (iColumn=0;iColumn<ncolsNow;iColumn++) {
	  int kColumn = originalColumn_[iColumn];
	  columnNames.push_back(originalModel->columnName(kColumn));
	}
	presolvedModel_->copyNames(rowNames,columnNames);
      } else {
	presolvedModel_->setLengthNames(0);
      }
#endif
      if (rowObjective_) {
	int iRow;
        int k=-1;
        int nObj=0;
	for (iRow=0;iRow<nrowsNow;iRow++) {
	  int kRow = originalRow_[iRow];
          assert (kRow>k);
          k=kRow;
          rowObjective_[iRow]=rowObjective_[kRow];
          if (rowObjective_[iRow])
            nObj++;
	}
        if (nObj) {
          printf("%d costed slacks\n",nObj);
          presolvedModel_->setRowObjective(rowObjective_);
        }
      }
      // now clean up integer variables.  This can modify original
      int i;
      const char * information = presolvedModel_->integerInformation();
      if (information) {
	int numberChanges=0;
	double * lower0 = originalModel_->columnLower();
	double * upper0 = originalModel_->columnUpper();
	double * lower = presolvedModel_->columnLower();
	double * upper = presolvedModel_->columnUpper();
	for (i=0;i<ncolsNow;i++) {
	  if (!information[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);
	  lower[i]=lowerValue;
	  upper[i]=upperValue;
	  if (lowerValue>upperValue) {
	    numberChanges++;
	    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
						       messages)
							 <<iOriginal
							 <<lowerValue
							 <<upperValue
							 <<CoinMessageEol;
	    result=1;
	  } else {
	    if (lowerValue>lowerValue0+1.0e-8) {
	      lower0[iOriginal] = lowerValue;
	      numberChanges++;
	    }
	    if (upperValue<upperValue0-1.0e-8) {
	      upper0[iOriginal] = upperValue;
	      numberChanges++;
	    }
	  }	  
	}
	if (numberChanges) {
	  presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
						     messages)
						       <<numberChanges
						       <<CoinMessageEol;
	  if (!result&&totalPasses>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_) {
      // infeasible or unbounded
      result=1;
      originalModel->setProblemStatus(prob.status_);
    } else {
      // no changes - model needs restoring after Lou's changes
#ifndef CLP_NO_STD
      if (saveFile_=="") {
#endif
	delete presolvedModel_;
	presolvedModel_ = new ClpSimplex(*originalModel);
	// but we need to remove gaps
	ClpPackedMatrix* clpMatrix =
	  dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
	if (clpMatrix) {
	  clpMatrix->getPackedMatrix()->removeGaps();
	}
#ifndef CLP_NO_STD
      } else {
	presolvedModel_=originalModel;
      }
      presolvedModel_->dropNames();
#endif
      
      // drop integer information if wanted
      if (!keepIntegers)
	presolvedModel_->deleteIntegerInformation();
      result=2;
    }
  }
  if (result==0||result==2) {
    int nrowsAfter = presolvedModel_->getNumRows();
    int ncolsAfter = presolvedModel_->getNumCols();
    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
					       messages)
						 <<nrowsAfter<< -(nrows_ - nrowsAfter)
						 << ncolsAfter<< -(ncols_ - ncolsAfter)
						 <<nelsAfter<< -(nelems_ - nelsAfter)
						 <<CoinMessageEol;
  } else {
    destroyPresolve();
    if (presolvedModel_!=originalModel_)
      delete presolvedModel_;
    presolvedModel_=NULL;
  }
  return presolvedModel_;
}