// 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).

#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cfloat>
#include <cassert>
#include <iostream>
//#define CGL_DEBUG 1
//#ifdef NDEBUG
//#undef NDEBUG
//#endif
#include "CoinPragma.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedVector.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinIndexedVector.hpp"
#include "OsiRowCutDebugger.hpp"
#define COIN_HAS_CLP_GOMORY
#ifdef COIN_HAS_CLP_GOMORY
#include "OsiClpSolverInterface.hpp"
#endif
#include "CoinFactorization.hpp"
#undef CLP_OSL
#if 1
#define CLP_OSL 1
#if CLP_OSL!=1&&CLP_OSL!=3
#undef CLP_OSL
#else
#include "CoinOslFactorization.hpp"
#endif
#endif
#include "CoinWarmStartBasis.hpp"
#include "CglGomory.hpp"
#include "CoinFinite.hpp"
#ifdef CGL_DEBUG_GOMORY
int gomory_try=CGL_DEBUG_GOMORY;
#endif
//-------------------------------------------------------------------
// Generate Gomory cuts
//------------------------------------------------------------------- 
void CglGomory::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
			     const CglTreeInfo info)
{
#ifdef CGL_DEBUG_GOMORY
  gomory_try++;
#endif
  // Get basic problem information
  int numberColumns=si.getNumCols(); 
  
  // get integer variables and basis
  char * intVar = new char[numberColumns];
  int i;
  CoinWarmStart * warmstart = si.getWarmStart();
  CoinWarmStartBasis* warm =
    dynamic_cast<CoinWarmStartBasis*>(warmstart);
  const double * colUpper = si.getColUpper();
  const double * colLower = si.getColLower();
  //#define CLP_INVESTIGATE2
#ifndef CLP_INVESTIGATE2
  if ((info.options&16)!=0)
#endif
    printf("%d %d %d\n",info.inTree,info.options,info.pass);
  for (i=0;i<numberColumns;i++) {
    if (si.isInteger(i)) {
      if (colUpper[i]>colLower[i]+0.5) {
	if (fabs(colUpper[i]-1.0)<1.0e-12&&fabs(colLower[i])<1.0e-12) {
	  intVar[i]=1; //0-1
	} else  if (colLower[i]>=0.0) {
	  intVar[i] = 2; // other
	} else {
	  // negative bounds - I am not sure works
	  intVar[i] = 0;
	}
      } else {
	intVar[i] = 0;
      }
    } else {
      intVar[i]=0;
    }
  }
  const OsiSolverInterface * useSolver=&si;
#ifdef COIN_HAS_CLP_GOMORY
  double * objective = NULL;
  OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
  int numberOriginalRows = -1;
  if (clpSolver) {
    useSolver = originalSolver_;
    assert (gomoryType_);
    // check simplex is plausible
    if (!clpSolver->getNumRows()||numberColumns!=clpSolver->getNumCols()) {
      delete originalSolver_;
      originalSolver_=si.clone();
      clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
      assert (clpSolver);
      useSolver = originalSolver_;
    }
    ClpSimplex * simplex = clpSolver->getModelPtr();
    numberOriginalRows = simplex->numberRows();
    int numberRows = si.getNumRows();
    assert (numberOriginalRows<=numberRows);
    // only do if different (unless type 2x)
    int gomoryType = gomoryType_%10;
    int whenToDo = gomoryType_/10;
    if (whenToDo==2 ||(numberRows>numberOriginalRows && whenToDo==1
		       && (info.options&512)==0) ||
	((info.options&1024)!=0 && (info.options&512)==0
	 && numberTimesStalled_<3)) {
      // bounds
      memcpy(simplex->columnLower(),colLower,numberColumns*sizeof(double));
      memcpy(simplex->columnUpper(),colUpper,numberColumns*sizeof(double));
      double * obj = simplex->objective();
      objective = CoinCopyOfArray(obj,numberColumns);
      const double * pi = si.getRowPrice();
      const CoinPackedMatrix * rowCopy = si.getMatrixByRow();
      const int * column = rowCopy->getIndices();
      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
      const int * rowLength = rowCopy->getVectorLengths(); 
      const double * rowElements = rowCopy->getElements();
      const double * rowLower = si.getRowLower();
      const double * rowUpper = si.getRowUpper();
      int numberCopy;
      int numberAdd;
      double * rowLower2 = NULL;
      double * rowUpper2 = NULL;
      int * column2 = NULL;
      CoinBigIndex * rowStart2 = NULL;
      double * rowElements2 = NULL;
      char * copy = new char [numberRows-numberOriginalRows];
      memset(copy,0,numberRows-numberOriginalRows);
      if (gomoryType==2) {
	numberCopy=0;
	numberAdd=0;
	for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
	  bool simple = true;
	  for (int k=rowStart[iRow];
	       k<rowStart[iRow]+rowLength[iRow];k++) {
	    double value = rowElements[k];
	    if (value!=floor(value+0.5)) {
	      simple=false;
	      break;
	    }
	  }
	  if (simple) {
	    numberCopy++;
	    numberAdd+=rowLength[iRow];
	    copy[iRow-numberOriginalRows]=1;
	  }
	}
	if (numberCopy) {
	  //printf("Using %d rows out of %d\n",numberCopy,numberRows-numberOriginalRows);
	  rowLower2 = new double [numberCopy];
	  rowUpper2 = new double [numberCopy];
	  rowStart2 = new CoinBigIndex [numberCopy+1];
	  rowStart2[0]=0;
	  column2 = new int [numberAdd];
	  rowElements2 = new double [numberAdd];
	}
      }
      numberCopy=0;
      numberAdd=0;
      const double * rowSolution = si.getRowActivity();
      double offset=0.0;
      for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
	if (!copy[iRow-numberOriginalRows]) {
	  double value = pi[iRow];
	  offset += rowSolution[iRow]*value;
	  for (int k=rowStart[iRow];
	       k<rowStart[iRow]+rowLength[iRow];k++) {
	    int iColumn=column[k];
	    obj[iColumn] -= value*rowElements[k];
	  }
	} else {
	  rowLower2[numberCopy]=rowLower[iRow];
	  rowUpper2[numberCopy]=rowUpper[iRow];
	  for (int k=rowStart[iRow];
	       k<rowStart[iRow]+rowLength[iRow];k++) {
	    column2[numberAdd]=column[k];
	    rowElements2[numberAdd++]=rowElements[k];
	  }
	  numberCopy++;
	  rowStart2[numberCopy]=numberAdd;
	}
      }
#if 0
      CoinThreadRandom randomNumberGenerator;
      const double * solution = si.getColSolution();
      for (int i=0;i<numberColumns;i++) {
	if (intVar[i]==1) {
	  double randomNumber = randomNumberGenerator.randomDouble();
	  //randomNumber = 0.001*floor(randomNumber*1000.0);
	  if (solution[i]>0.5)
	    obj[i] -= randomNumber*0.001*fabs(obj[i]);
	  else
	    obj[i] += randomNumber*0.001*fabs(obj[i]);
	}
      }
#endif
      if (numberCopy) {
	clpSolver->addRows(numberCopy,
			   rowStart2,column2,rowElements2,
			   rowLower2,rowUpper2);
	delete [] rowLower2 ;
	delete [] rowUpper2 ;
	delete [] column2 ;
	delete [] rowStart2 ;
	delete [] rowElements2 ;
      }
      delete [] copy;
      memcpy(simplex->primalColumnSolution(),si.getColSolution(),
	     numberColumns*sizeof(double));
      warm->resize(numberOriginalRows,numberColumns);
      clpSolver->setBasis(*warm);
      delete warm;
      simplex->setDualObjectiveLimit(COIN_DBL_MAX);
      simplex->setLogLevel(0);
      simplex->primal(1);
      // check basis
      int numberTotal=simplex->numberRows()+simplex->numberColumns();
      int superbasic=0;
      for (int i=0;i<numberTotal;i++) {
	if (simplex->getStatus(i)==ClpSimplex::superBasic)
	  superbasic++;
      }
      if (superbasic) {
	//printf("%d superbasic!\n",superbasic);
	simplex->dual();
	superbasic=0;
	for (int i=0;i<numberTotal;i++) {
	  if (simplex->getStatus(i)==ClpSimplex::superBasic)
	    superbasic++;
	}
	assert (!superbasic);
      }
      //printf("Trying - %d its status %d objs %g %g - with offset %g\n",
      //     simplex->numberIterations(),simplex->status(),
      //     simplex->objectiveValue(),si.getObjValue(),simplex->objectiveValue()+offset);
      //simplex->setLogLevel(0);
      warm=simplex->getBasis();
      warmstart=warm;
      if (simplex->status()) {
	//printf("BAD status %d\n",simplex->status());
	//clpSolver->writeMps("clp");
	//si.writeMps("si");
	delete [] objective;
	objective=NULL;
	useSolver=&si;
      }
    } else {
      // don't do
      delete warmstart;
      warmstart=NULL;
      if ((info.options&1024)==0)
	numberTimesStalled_=0;
    }
  }
#endif
#ifdef CGL_DEBUG
  const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
  if (debugger&&!debugger->onOptimalPath(si))
    debugger = NULL;
#else
  const OsiRowCutDebugger * debugger = NULL;
#endif
  int numberRowCutsBefore = cs.sizeRowCuts();

  if (warmstart)
    generateCuts(debugger, cs, *useSolver->getMatrixByCol(), 
		 *useSolver->getMatrixByRow(),
		 useSolver->getColSolution(),
		 useSolver->getColLower(), useSolver->getColUpper(), 
		 useSolver->getRowLower(), useSolver->getRowUpper(),
		 intVar,warm,info);
#ifdef COIN_HAS_CLP_GOMORY
  if (objective) {
    ClpSimplex * simplex = clpSolver->getModelPtr();
    memcpy(simplex->objective(),objective,numberColumns*sizeof(double));
    delete [] objective;
    // take out locally useless cuts
    const double * solution = si.getColSolution();
    double primalTolerance = 1.0e-7;
    int numberRowCutsAfter = cs.sizeRowCuts();
    for (int k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) {
      const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
      double sum = 0.0;
      int n = thisCut->row().getNumElements();
      const int * column = thisCut->row().getIndices();
      const double * element = thisCut->row().getElements();
      assert (n);
      for (int i = 0; i < n; i++) {
	double value = element[i];
	sum += value * solution[column[i]];
      }
      if (sum > thisCut->ub() + primalTolerance) {
	sum = sum - thisCut->ub();
      } else if (sum < thisCut->lb() - primalTolerance) {
	sum = thisCut->lb() - sum;
      } else {
	sum = 0.0;
      }
      if (!sum) {
	// take out
	cs.eraseRowCut(k);
      }
    }
#ifdef CLP_INVESTIGATE2
    printf("OR %p pass %d inTree %c - %d cuts (but %d deleted)\n",
       originalSolver_,info.pass,info.inTree?'Y':'N',
       numberRowCutsAfter-numberRowCutsBefore,
       numberRowCutsAfter-cs.sizeRowCuts());
#endif
  }
#endif

  delete warmstart;
  delete [] intVar;
  if ((!info.inTree&&((info.options&4)==4||((info.options&8)&&!info.pass)))
      ||(info.options&16)!=0) {
    int limit = maximumLengthOfCutInTree();
    int numberRowCutsAfter = cs.sizeRowCuts();
    for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++) {
      int length = cs.rowCutPtr(i)->row().getNumElements();
      if (length<=limit)
	cs.rowCutPtr(i)->setGloballyValid();
    }
  }
  if ((gomoryType_%10)==2) {
    // back to original
    assert(clpSolver);
    int numberRows = clpSolver->getNumRows();
    if (numberRows>numberOriginalRows) {
      int numberDelete = numberRows-numberOriginalRows;
      int * delRow = new int [numberDelete];
      for (int i=0;i<numberDelete;i++)
	delRow[i]=i+numberOriginalRows;
      clpSolver->deleteRows(numberDelete,delRow);
      delete [] delRow;
    }
  }
}

// Returns value - floor but allowing for small errors
inline double above_integer(double value) {
  double value2=floor(value);
  double value3=floor(value+0.5);
  if (fabs(value3-value)<1.0e-7*(fabs(value3)+1.0))
    return 0.0;
  return value-value2;
}
//-------------------------------------------------------------------
// Returns the greatest common denominator of two 
// positive integers, a and b, found using Euclid's algorithm 
//-------------------------------------------------------------------
static int gcd(int a, int b) 
{
  int remainder = -1;
#if CGL_DEBUG>1
  printf("gcd of %d and %d\n",a,b);
  int nLoop=0;
#endif
  // make sure a<=b (will always remain so)
  if(a > b) {
    // Swap a and b
    int temp = a;
    a = b;
    b = temp;
  }
  // if zero then gcd is nonzero (zero may occur in rhs of packed)
  if (!a) {
    if (b) {
      return b;
    } else {
      printf("**** gcd given two zeros!!\n");
      abort();
    }
  }
  while (remainder) {

#if CGL_DEBUG>1
    nLoop++;
    if (nLoop>50) {
      abort();
      return -1;
    }
#endif
    remainder = b % a;
    b = a;
    a = remainder;
  }
#if CGL_DEBUG>1
  printf("=> %d\n",b);
#endif
  return b;
}

//-------------------------------------------------------------------
// Returns the nearest rational with denominator < maxDenominator
//-------------------------------------------------------------------
typedef struct {
  int numerator;
  int denominator;
} Rational;
inline Rational nearestRational(double value, int maxDenominator) 
{
  Rational tryThis;
  Rational tryA;
  Rational tryB;
  double integerPart;

#if CGL_DEBUG>1
  printf("Rational of %g is ",value);
#endif
  int nLoop=0;

  tryA.numerator=0;
  tryA.denominator=1;
  tryB.numerator=1;
  tryB.denominator=0;

  if (fabs(value)<1.0e-10)
    return tryA;
  integerPart = floor(value);
  value -= integerPart;
  tryThis.numerator = tryB.numerator* static_cast<int> (integerPart) + tryA.numerator;
  tryThis.denominator = tryB.denominator* static_cast<int> (integerPart) + tryA.denominator;
  tryA = tryB;
  tryB = tryThis;

  while (value>1.0e-10 && tryB.denominator <=maxDenominator) {
    nLoop++;
    if (nLoop>50) {
      Rational bad;
      bad.numerator=-1;
      bad.denominator=-1;
#if CGL_DEBUG>1
      printf(" *** bad rational\n");
#endif
      return bad;
    }
    value = 1.0/value;
    integerPart = floor(value+1.0e-10);
    value -= integerPart;
    tryThis.numerator = tryB.numerator* static_cast<int> (integerPart) + tryA.numerator;
    tryThis.denominator = tryB.denominator* static_cast<int>(integerPart) + tryA.denominator;
    tryA = tryB;
    tryB = tryThis;
  }
  if (tryB.denominator <= maxDenominator) {
#if CGL_DEBUG>1
    printf("%d/%d\n",tryB.numerator,tryB.denominator);
#endif
    return tryB;
  } else {
#if CGL_DEBUG>1
    printf("%d/%d\n",tryA.numerator,tryA.denominator);
#endif
    return tryA;
  }
}
// Does actual work - returns number of cuts
int
CglGomory::generateCuts( 
#ifdef CGL_DEBUG
			const OsiRowCutDebugger * debugger, 
#else
			const OsiRowCutDebugger * , 
#endif
                         OsiCuts & cs,
                         const CoinPackedMatrix & columnCopy,
                         const CoinPackedMatrix & rowCopy,
                         const double * colsol,
                         const double * colLower, const double * colUpper,
                         const double * rowLower, const double * rowUpper,
			 const char * intVar,
                         const CoinWarmStartBasis* warm,
                         const CglTreeInfo info)
{
  int infoOptions=info.options;
  bool globalCuts = (infoOptions&16)!=0;
  double testFixed = (!globalCuts) ? 1.0e-8 : -1.0;
  // get what to look at
  double away = info.inTree ? away_ : CoinMin(away_,awayAtRoot_);
  int numberRows=columnCopy.getNumRows();
  int numberColumns=columnCopy.getNumCols(); 
  int numberElements=columnCopy.getNumElements();
  // Allow bigger length on initial matrix (if special setting)
  //if (limit==512&&!info.inTree&&!info.pass)
  //limit=1024;
  // Start of code to create a factorization from warm start (A) ====
  // check factorization is okay
  CoinFactorization factorization;
#ifdef CLP_OSL
  CoinOslFactorization * factorization2=NULL;
  if (alternateFactorization_) {
    factorization2 = new CoinOslFactorization();
  }
#endif
  // We can either set increasing rows so ...IsBasic gives pivot row
  // or we can just increment iBasic one by one
  // for now let ...iBasic give pivot row
  int status=-100;
  // probably could use pivotVariables from OsiSimplexModel
  int * rowIsBasic = new int[numberRows];
  int * columnIsBasic = new int[numberColumns];
  int i;
  int numberBasic=0;
  for (i=0;i<numberRows;i++) {
    if (warm->getArtifStatus(i) == CoinWarmStartBasis::basic) {
      rowIsBasic[i]=1;
      numberBasic++;
    } else {
      rowIsBasic[i]=-1;
    }
  }
  for (i=0;i<numberColumns;i++) {
    if (warm->getStructStatus(i) == CoinWarmStartBasis::basic) {
      columnIsBasic[i]=1;
      numberBasic++;
    } else {
      columnIsBasic[i]=-1;
    }
  }
  //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
  while (status<-98) {
#ifdef CLP_OSL
    if (!alternateFactorization_) {
#endif
      status=factorization.factorize(columnCopy,
				     rowIsBasic, columnIsBasic);
      if (status==-99) factorization.areaFactor(factorization.areaFactor() * 2.0);
#ifdef CLP_OSL
    } else {
      double areaFactor=1.0;
      status=factorization2->factorize(columnCopy,
				      rowIsBasic, columnIsBasic,areaFactor);
      if (status==-99) areaFactor *= 2.0;
    }
#endif
  } 
  if (status) {
#ifdef COIN_DEVELOP
    std::cout<<"Bad factorization of basis - status "<<status<<std::endl;
#endif
    delete [] rowIsBasic;
    delete [] columnIsBasic;
    return -1;
  }
  // End of creation of factorization (A) ====
  
#ifdef CLP_OSL
  double relaxation = !alternateFactorization_ ? factorization.conditionNumber() :
    factorization2->conditionNumber();
#else
  double relaxation = factorization.conditionNumber();
#endif
#ifdef COIN_DEVELOP_z
  if (relaxation>1.0e49)
    printf("condition %g\n",relaxation);
#endif
  relaxation *= conditionNumberMultiplier_;
  double bounds[2]={-COIN_DBL_MAX,0.0};
  int iColumn,iRow;

  const int * column = rowCopy.getIndices();
  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
  const int * rowLength = rowCopy.getVectorLengths(); 
  const double * rowElements = rowCopy.getElements();
  const int * row = columnCopy.getIndices();
  const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
  const int * columnLength = columnCopy.getVectorLengths(); 
  const double * columnElements = columnCopy.getElements();

  // we need to do book-keeping for variables at ub
  double tolerance = 1.0e-7;
  bool * swap= new bool [numberColumns];
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (columnIsBasic[iColumn]<0&&
	colUpper[iColumn]-colsol[iColumn]<=tolerance) {
      swap[iColumn]=true;
    } else {
      swap[iColumn]=false;
    }
  }

  // get row activities (could use solver but lets do here )
  double * rowActivity = new double [numberRows];
  CoinFillN(rowActivity,numberRows,0.0);
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    double value = colsol[iColumn];
    int k;
    for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];k++) {
      iRow = row[k];
      rowActivity[iRow] += columnElements[k]*value;
    }
  }
  /* we need to mark rows:
     0) equality
     1) slack at lb (activity at ub)
     2) slack at ub (activity at lb)
     and 4 bit is set if slack must be integer

  */
  int * rowType = new int[numberRows];
  for (iRow=0;iRow<numberRows;iRow++) {
    if (rowIsBasic[iRow]<0&&rowUpper[iRow]>rowLower[iRow]+1.0e-7) {
      int type=0;
      double rhs=0.0;
      if (rowActivity[iRow]>=rowUpper[iRow]-1.0e-7) {
	type=1;
	rhs=rowUpper[iRow];
      } else if (rowActivity[iRow]<=rowLower[iRow]+1.0e-7) {
	type=2;
	rhs=rowLower[iRow];
      } else {
	// probably large rhs
	if (rowActivity[iRow]-rowLower[iRow]<
	    rowUpper[iRow]-rowActivity[iRow])
	  rowType[iRow]=2;
	else
	  rowType[iRow]=1;
#ifdef CGL_DEBUG
	assert (CoinMin(rowUpper[iRow]-rowActivity[iRow],
		    rowActivity[iRow]-rowUpper[iRow])<1.0e-5);
	//abort();
        continue;
#else
	continue;
#endif
      }
      if (above_integer(rhs)<1.0e-10) {
	// could be integer slack
	bool allInteger=true;
	int k;
	for (k=rowStart[iRow];
	     k<rowStart[iRow]+rowLength[iRow];k++) {
	  int iColumn=column[k];
	  if (!intVar[iColumn]||above_integer(rowElements[k])>1.0e-10) {
	    // not integer slacks
	    allInteger=false;
	    break;
	  }
	}
	if (allInteger) {
	  type |= 4;
	}
      }
      rowType[iRow]=type;
    } else {
      // row is equality or basic
      rowType[iRow]=0;
    }
  }

  // Start of code to create work arrays for factorization (B) ====
  // two vectors for updating (one is work)
  CoinIndexedVector work;
  CoinIndexedVector array;
  // make sure large enough
  work.reserve(numberRows);
  array.reserve(numberRows);
  int * arrayRows = array.getIndices();
  double * arrayElements = array.denseVector();
  // End of code to create work arrays (B) ====

  int numberAdded=0;
  // we also need somewhere to accumulate cut
  CoinIndexedVector cutVector;
  cutVector.reserve(numberColumns+1);
  int * cutIndex = cutVector.getIndices();
  double * cutElement = cutVector.denseVector(); 
  // and for packed form (as not necessarily in order)
  // also space for sort
  bool doSorted = (infoOptions&256)!=0;
  int lengthArray = static_cast<int>(numberColumns+1+((numberColumns+1)*sizeof(int))/sizeof(double));
  if (doSorted)
    lengthArray+=numberColumns;
  double * packed = new double[lengthArray]; 
  double * sort = packed+numberColumns+1;
  int * which = reinterpret_cast<int *>(doSorted ? (sort+numberColumns): (sort));
  double tolerance1=1.0e-6;
  double tolerance2=0.9;
  double tolerance3=1.0e-4;
  double tolerance6=1.0e-6;
  double tolerance9=1.0e-4;
#define MORE_GOMORY_CUTS 1
#ifdef CLP_INVESTIGATE2
  int saveLimit = info.inTree ? 50 : 1000;
#else
#if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
  int saveLimit;
#endif  
#endif  
  // get limit on length of cut
  int limit = 0;
  if (!limit_)
    dynamicLimitInTree_ = CoinMax(50,numberColumns/40);
  if (!info.inTree) {
    limit = limitAtRoot_;
    if (!info.pass) {
      tolerance1=1.0;
      tolerance2=1.0e-2;
      tolerance3=1.0e-6;
      tolerance6=1.0e-7;
      tolerance9=1.0e-5;
      if (!limit)
	limit=numberColumns;
    } else {
      if((infoOptions&32)==0/*&&numberTimesStalled_<3*/) {
	if (!limit) {
	  if(numberElements>8*numberColumns)
	    limit=numberColumns;
	  else
	    limit = CoinMax(1000,numberRows/4);
	}
      } else {
	limit=numberColumns;
	numberTimesStalled_++;
      } 
    }
  } else {
    limit = limit_;
    if (!limit) {
      if (!info.pass) 
	limit = dynamicLimitInTree_;
      else
	limit=50;
    }
  }
  // If big - allow for rows
  if (limit>=numberColumns)
    limit += numberRows;
#ifdef CLP_INVESTIGATE2
  if (limit>saveLimit&&!info.inTree&&(infoOptions&512)==0) 
    printf("Gomory limit changed from %d to %d, inTree %c, pass %d, r %d,c %d,e %d\n",
	   saveLimit,limit,info.inTree ? 'Y' : 'N',info.pass,
	   numberRows,numberColumns,numberElements);
#endif
  int nCandidates=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    // This returns pivot row for columns or -1 if not basic (C) ====
    int iBasic=columnIsBasic[iColumn];
    if (iBasic>=0&&intVar[iColumn]) {
      double reducedValue=above_integer(colsol[iColumn]);
      //printf("col %d bas %d val %.18g\n",iColumn,iBasic,colsol[iColumn]);
      if(intVar[iColumn]&&reducedValue<1.0-away&&reducedValue>away) {
	if (doSorted)
	  sort[nCandidates]=fabs(0.5-reducedValue);
	which[nCandidates++]=iColumn;
      }
    }
  }
  int nTotalEls=COIN_INT_MAX;
  if (doSorted) {
    CoinSort_2(sort,sort+nCandidates,which);
    int nElsNow = columnCopy.getNumElements();
    int nAdd;
    int nAdd2;
    int nReasonable;
    int depth=info.level;
    if (depth<2) {
      nAdd=10000;
      if (info.pass>0)
	nAdd = CoinMin(nAdd,nElsNow+2*numberRows);
      nAdd2 = 5*numberColumns;
      nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
      if (!depth&&!info.pass) {
	// allow more
	nAdd += nElsNow/2;
	nAdd2 += nElsNow/2;
	nReasonable += nElsNow/2;
	limit=numberRows+numberColumns;
      }
    } else {
      nAdd = 200;
      nAdd2 = 2*numberColumns;
      nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
    }
    nTotalEls=nReasonable;
  }
#ifdef MORE_GOMORY_CUTS
  int saveTotalEls=nTotalEls;
#endif
#if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
  saveLimit=limit;
  if (doSorted)
    limit=numberRows+numberColumns;
  OsiCuts longCuts;
#endif
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  OsiCuts secondaryCuts;
#endif
  for (int kColumn=0;kColumn<nCandidates;kColumn++) {
    if (nTotalEls<=0)
      break;  // Got enough
    iColumn=which[kColumn];
    double reducedValue=above_integer(colsol[iColumn]);;
    // This returns pivot row for columns or -1 if not basic (C) ====
    int iBasic=columnIsBasic[iColumn];
    double ratio=reducedValue/(1.0-reducedValue);
    if (iBasic>=0) {
  // Debug code below computes tableau column of basic ====
      int j;
#ifdef CGL_DEBUG
      {
	// put column into array
	array.setVector(columnLength[iColumn],row+columnStart[iColumn],
			columnElements+columnStart[iColumn]);
	// get column in tableau
#ifdef CLP_OSL
	if (!alternateFactorization_)
#endif
	  factorization.updateColumn ( &work, &array );
#ifdef CLP_OSL
	else
	  factorization2->updateColumn ( &work, &array );
#endif
	int nn=0;
	int numberInArray=array.getNumElements();
	for (j=0;j<numberInArray;j++) {
	  int indexValue=arrayRows[j];
	  double value=arrayElements[indexValue];
	  if (fabs(value)>1.0e-5) {
	    assert (fabs(value-1.0)<1.0e-7);
	    assert (indexValue==iBasic);
	    nn++;
	  }
	}
	assert (nn==1);
	array.clear();
	work.checkClear();
      }
#endif
      array.clear();
      assert(intVar[iColumn]&&reducedValue<1.0-away&&reducedValue>away);
      {
#ifdef CGL_DEBUG
	cutVector.checkClear();
#endif
	// get row of tableau
	double one =1.0;
	array.setVector(1,&iBasic,&one);
	int numberNonInteger=0;
	//Code below computes tableau row ====
	// get pi
#ifdef CLP_OSL
	if (!alternateFactorization_)
#endif
	  factorization.updateColumnTranspose ( &work, &array );
#ifdef CLP_OSL
	else
	  factorization2->updateColumnTranspose ( &work, &array );
#endif
	int numberInArray=array.getNumElements();
#ifdef CGL_DEBUG
	// check pivot on iColumn
	{
	  double value=0.0;
	  int k;
	  // add in row of tableau
	  for (k=columnStart[iColumn];
	       k<columnStart[iColumn]+columnLength[iColumn];k++) {
	    iRow = row[k];
	    value += columnElements[k]*arrayElements[iRow];
	  }
	  // should be 1
	  assert (fabs(value-1.0) < 1.0e-7);
	}
#endif
	double largestFactor=0.0;
	for (j=0;j<numberInArray;j++) {
	  int indexValue=arrayRows[j];
	  double value=arrayElements[indexValue];
	  largestFactor = CoinMax(largestFactor,fabs(value));
	}
	//reducedValue=colsol[iColumn];
	// coding from pg 130 of Wolsey 
	// adjustment to rhs
	double rhs=0.0;
	int number=0;
#ifdef CGL_DEBUG_GOMORY
	    if (!gomory_try)
	      printf("start for basic column %d\n",iColumn);
#endif
	// columns
	for (j=0;j<numberColumns;j++) {
	  if (columnIsBasic[j]<0&&colUpper[j]>colLower[j]+testFixed) {
	    double value=0.0;
	    int k;
	    // add in row of tableau
	    for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
	      iRow = row[k];
	      double value2 = columnElements[k]*arrayElements[iRow];
	      largestFactor = CoinMax(largestFactor,fabs(value2));
	      value += value2;
	    }
#ifdef CGL_DEBUG_GOMORY
	    if (!gomory_try&&value)
	      printf("col %d alpha %g colsol %g swap %c bounds %g %g\n",
		     j,value,colsol[j],swap[j] ? 'Y' : 'N',
		     colLower[j],colUpper[j]);
#endif
	    // value is entry in tableau row end (C) ====
	    if (fabs(value)<1.0e-16) {
	      // small value
	      continue;
	    } else {
	      // left in to stop over compilation?
	      //if (iColumn==-52) printf("for basic %d, column %d has alpha %g, colsol %g\n",
	      //		      iColumn,j,value,colsol[j]);
#if CGL_DEBUG>1
	      if (iColumn==52) printf("for basic %d, column %d has alpha %g, colsol %g\n",
				      iColumn,j,value,colsol[j]);
#endif
	      // deal with bounds
	      if (swap[j]) {
		//reducedValue -= value*colUpper[j];
		// negate
		value = - value;
	      } else {
		//reducedValue -= value*colLower[j];
	      }
#if CGL_DEBUG>1
	      if (iColumn==52) printf("%d value %g reduced %g int %d rhs %g swap %d\n",
				      j,value,reducedValue,intVar[j],rhs,swap[j]);
#endif
	      double coefficient;
	      if (intVar[j]) {
		// integer
		coefficient = above_integer(value);
		if (coefficient > reducedValue) {
		  coefficient = ratio * (1.0-coefficient);
		} 
	      } else {
		// continuous
		numberNonInteger++;
		if (value > 0.0) {
		  coefficient = value;
		} else {
		  //??? sign wrong in book
		  coefficient = -ratio*value;
		}
	      }
	      if (swap[j]) {
		// negate
		coefficient = - coefficient;
		rhs += colUpper[j]*coefficient;
	      } else {
		rhs += colLower[j]*coefficient;
	      }
	      if (fabs(coefficient)>= COIN_INDEXED_TINY_ELEMENT) {
		cutElement[j] = coefficient;
		cutIndex[number++]=j;
		// If too many - break from loop
		if (number>limit) 
		  break;
	      }
	    }
	  } else {
	    // basic
	    continue;
	  }
	}
	cutVector.setNumElements(number);
	// If too many - just clear vector and skip
	if (number>limit) {
	  cutVector.clear();
	  continue;
	}
	//check will be cut
	//reducedValue=above_integer(reducedValue);
	rhs += reducedValue;
	double violation = reducedValue;
#ifdef CGL_DEBUG
	std::cout<<"cut has violation of "<<violation
		 <<" value "<<colsol[iColumn]<<std::endl;
#endif
	// now do slacks part
	for (j=0;j<numberInArray;j++) {
	  iRow=arrayRows[j];
	  double value = arrayElements[iRow];
	  int type=rowType[iRow];
	  if (type&&fabs(value)>=1.0e-16) {
	    if ((type&1)==0) {
	      // negate to get correct coefficient
	      value = - value;
	    }
	    double coefficient;
	    if ((type&4)!=0) {
	      // integer
	      coefficient = above_integer(value);
	      if (coefficient > reducedValue) {
		coefficient = ratio * (1.0-coefficient);
	      } 
	    } else {
	      numberNonInteger++;
	      // continuous
	      if (value > 0.0) {
		coefficient = value;
	      } else {
		coefficient = -ratio*value;
	      }
	    }
	    if ((type&1)!=0) {
	      // slack at ub - treat as +1.0
	      rhs -= coefficient*rowUpper[iRow];
	    } else {
	      // negate yet again ?
	      coefficient = - coefficient;
	      rhs -= coefficient*rowLower[iRow];
	    }
	    int k;
	    for (k=rowStart[iRow];
		 k<rowStart[iRow]+rowLength[iRow];k++) {
	      int jColumn=column[k];
	      double value=rowElements[k];
	      cutVector.quickAdd(jColumn,-coefficient*value);
	    }
	  }
	}
	//check again and pack down
	// also change signs
	// also zero cutElement
	double sum=0.0;
	rhs = - rhs;
	int n = cutVector.getNumElements();
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
	double violation2=violation;
#endif
	number=0;
	for (j=0;j<n;j++) {
	  int jColumn =cutIndex[j];
	  double value=-cutElement[jColumn];
	  cutElement[jColumn]=0.0;
	  if (fabs(value)>1.0e-8) {
	    sum+=value*colsol[jColumn];
	    packed[number]=value;
	    cutIndex[number++]=jColumn;
          } else {
#define LARGE_BOUND 1.0e20
            // small - adjust rhs if rhs reasonable
            if (value>0.0&&colLower[jColumn]>-LARGE_BOUND) {
              rhs -= value*colLower[jColumn];
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
	      // weaken violation
	      violation2 -= fabs(value*(colsol[jColumn]-colLower[jColumn]));
#endif
            } else if (value<0.0&&colUpper[jColumn]<LARGE_BOUND) {
              rhs -= value*colUpper[jColumn];
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
	      // weaken violation
	      violation2 -= fabs(value*(colsol[jColumn]-colUpper[jColumn]));
#endif
            } else if (fabs(value)>1.0e-13) {
              // take anyway
              sum+=value*colsol[jColumn];
              packed[number]=value;
              cutIndex[number++]=jColumn;
            } 
          }
	}
	// Final test on number
	//if (number>limit)
	//continue;
	// say zeroed out
	cutVector.setNumElements(0);
	bool accurate2=false;
	double difference=fabs((sum-rhs)-violation);
	double useTolerance;
	if (tolerance1>0.99) {
	  // use absolute
	  useTolerance = tolerance;
	} else {
	  double rhs2=CoinMax(fabs(rhs),10.0);
	  useTolerance=rhs2*0.1*tolerance1;
	}
	bool accurate = (difference<useTolerance);
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
	double difference2=fabs((sum-rhs)-violation2);
#if MORE_GOMORY_CUTS==1
	if (difference2<useTolerance&&doSorted) 
#else
	if (difference2<useTolerance&&doSorted&&number<saveLimit) 
#endif
	  accurate2=true;
#endif
	if (sum >rhs+tolerance2*away&&
	    (accurate||accurate2)) {
	  //#ifdef CGL_DEBUG
#ifdef CGL_DEBUG
#if CGL_DEBUG<=1
	  if (number<=-10) {
#endif
	    for (j=0;j<number;j++) {
	      std::cout<<" ["<<cutIndex[j]<<","<<packed[j]<<"]";
	    }
	    std::cout<<" <= "<<rhs<<std::endl;
#if CGL_DEBUG<=1
	  }
#endif
#endif
	  if (!numberNonInteger&&number) {
#ifdef CGL_DEBUG
	    assert (sizeof(Rational)==sizeof(double));
#endif
	    Rational * cleaned = reinterpret_cast<Rational *> (cutElement);
	    int * xInt = reinterpret_cast<int *> (cutElement);
	    // cut should have an integer slack so try and simplify
	    // add in rhs and put in cutElements (remember to zero later)
	    cutIndex[number]=numberColumns+1;
	    packed[number]=rhs;
	    int numberNonSmall=0;
	    int lcm = 1;
	    
	    for (j=0;j<number+1;j++) {
	      double value=above_integer(fabs(packed[j]));
	      if (fabs(value)<tolerance3) {
		// too small
		continue;
	      } else {
		numberNonSmall++;
	      }
	      
	      cleaned[j]=nearestRational(value,100000);
	      if (cleaned[j].denominator<0) {
		// bad rational
		lcm=-1;
		break;
	      }
	      int thisGcd = gcd(lcm,cleaned[j].denominator);
	      // may need to check for overflow
	      lcm /= thisGcd;
	      lcm *= cleaned[j].denominator;
	    }
	    if (lcm>0&&numberNonSmall) {
	      double multiplier = lcm;
	      int nOverflow = 0; 
	      for (j=0; j<number+1;j++) {
		double value = fabs(packed[j]);
		double dxInt = value*multiplier;
		xInt[j]= static_cast<int> (dxInt+0.5); 
#if CGL_DEBUG>1
		printf("%g => %g   \n",value,dxInt);
#endif
		if (dxInt>1.0e9||fabs(dxInt-xInt[j])> 1.0e-8) {
		  nOverflow++;
		  break;
		}
	      }
	      
	      if (nOverflow){
#ifdef CGL_DEBUG
		printf("Gomory Scaling: Warning: Overflow detected \n");
#endif
		numberNonInteger=-1;
	      } else {
		
		// find greatest common divisor of the elements
		j=0;
		while (!xInt[j])
		  j++; // skip zeros
		int thisGcd = gcd(xInt[j],xInt[j+1]);
		j++;
		for (;j<number+1;j++) {
		  if (xInt[j])
		    thisGcd = gcd(thisGcd,xInt[j]);
		}
#if 0
		// Check nothing too illegal - FIX this
		for (j=0;j<number+1;j++) {
		  double old = lcm*packed[j];
		  int newOne;
		  if (old>0.0)
		    newOne=xInt[j]/thisGcd;
		  else
		    newOne=-xInt[j]/thisGcd;
		  if (fabs(((double) newOne)-old)>
		      1.0e-10*(fabs(newOne)+1.0)) {
		    // say no good - first see if happens
		    printf("Fix this test 456 - just skip\n");
		    abort();
		  }
		} 
#endif		  
#if CGL_DEBUG>1
		printf("The gcd of xInt is %i\n",thisGcd);    
#endif
		
		// construct new cut by dividing through by gcd and 
		double minMultiplier=1.0e100;
		double maxMultiplier=0.0;
		for (j=0; j<number+1;j++) {
		  double old=packed[j];
		  if (old>0.0) {
		    packed[j]=xInt[j]/thisGcd;
		  } else {
		    packed[j]=-xInt[j]/thisGcd;
		  }
#if CGL_DEBUG>1
		  printf("%g => %g   \n",old,packed[j]);
#endif
		  if (packed[j]) {
		    if (fabs(packed[j])>maxMultiplier*fabs(old))
		      maxMultiplier = packed[j]/old;
		    if (fabs(packed[j])<minMultiplier*fabs(old))
		      minMultiplier = packed[j]/old;
		  }
		}
		rhs = packed[number];
#ifdef CGL_DEBUG
		printf("min, max multipliers - %g, %g\n",
		       minMultiplier,maxMultiplier);
#endif
		assert(maxMultiplier/minMultiplier>0.9999&&maxMultiplier/minMultiplier<1.0001);
	      }
	    }
	    // erase cutElement
	    CoinFillN(cutElement,number+1,0.0);
	  } else {
	    // relax rhs a tiny bit
	    rhs += 1.0e-8;
	    // relax if lots of elements for mixed gomory
	    if (number>=20) {
	      rhs  += 1.0e-7*(static_cast<double> (number/20));
	    }
	  }
	  // Take off tiny elements
	  // for first pass reject
#define TINY_ELEMENT 1.0e-12
	  {
	    int i,number2=number;
	    number=0;
	    double largest=0.0;
	    double smallest=1.0e30;
	    for (i=0;i<number2;i++) {
	      double value=fabs(packed[i]);
	      if (value<TINY_ELEMENT) {
		int iColumn = cutIndex[i];
		if (colUpper[iColumn]-colLower[iColumn]<LARGE_BOUND) {
		  // weaken cut
		  if (packed[i]>0.0) 
		    rhs -= value*colLower[iColumn];
		  else
		    rhs += value*colUpper[iColumn];
		} else {
		  // throw away
		  number=limit+1;
		  numberNonInteger=1;
		  break;
		}
	      } else {
		int iColumn = cutIndex[i];
		if (colUpper[iColumn]!=colLower[iColumn]||globalCuts) {
		  largest=CoinMax(largest,value);
		  smallest=CoinMin(smallest,value);
		  cutIndex[number]=cutIndex[i];
		  packed[number++]=packed[i];
		} else {
		  // fixed so subtract out
		  rhs -= packed[i]*colLower[iColumn];
		}
	      }
	    }
	    if (largest>1.0e10*smallest) {
	      number=limit+1; //reject
	      numberNonInteger=1;
	    } else if (largest>1.0e9*smallest) {
#ifdef CLP_INVESTIGATE2
	      printf("WOuld reject %g %g ratio %g\n",smallest,largest,
		     smallest/largest);
#endif
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
	      accurate=false;
#endif
	    }
	  }
	  if (number<limit||!numberNonInteger) {
	    bounds[1]=rhs;
	    if (number>50&&numberNonInteger)
	      bounds[1] = rhs+tolerance6+1.0e-8*fabs(rhs); // weaken
	    double test = CoinMin(largestFactor*largestFactorMultiplier_,
				  relaxation);
	    if (number>5&&numberNonInteger&&test>1.0e-20) {
#ifdef CLP_INVESTIGATE2
	      printf("relaxing rhs by %g - largestFactor was %g, rel %g\n",
	         CoinMin(test*fabs(rhs),tolerance9),largestFactor,relaxation);
#endif
	      //bounds[1] = CoinMax(bounds[1],
	      //		  rhs+CoinMin(test*fabs(rhs),tolerance9)); // weaken
	      bounds[1] = bounds[1]+CoinMin(test*fabs(rhs),tolerance9); // weaken
	    }
#ifdef MORE_GOMORY_CUTS
	    if (accurate) {
#else
	    {
#endif
	      OsiRowCut rc;
	      rc.setRow(number,cutIndex,packed,false);
	      rc.setLb(bounds[0]);
	      rc.setUb(bounds[1]);   
#if MORE_GOMORY_CUTS<2
	      nTotalEls -= number;
	      cs.insert(rc);
#else
	      if(number<saveLimit) {
		nTotalEls -= number;
		cs.insert(rc);
	      } else {
		longCuts.insert(rc);
	      }
#endif
	      //printf("nTot %d kCol %d iCol %d ibasic %d\n",
	      //     nTotalEls,kColumn,iColumn,iBasic);
	      numberAdded++;
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
	    } else if (accurate2) {
	      OsiRowCut rc;
	      rc.setRow(number,cutIndex,packed,false);
	      rc.setLb(bounds[0]);
	      rc.setUb(bounds[1]);   
	      secondaryCuts.insert(rc);
#endif
	    }
	  } else {
#ifdef CGL_DEBUG
	    std::cout<<"cut has "<<number<<" entries - skipped"
		     <<std::endl;
	    if (!number)
	      std::cout<<"******* Empty cut - infeasible"<<std::endl;
#endif
	  }
	} else {
	  // why dropped?
#ifdef CGL_DEBUG
	  std::cout<<"first violation "<<violation<<" now "
		   <<sum-rhs<<" why?, rhs= "<<rhs<<std::endl;
	  
	  for (j=0;j<number;j++) {
	    int jColumn =cutIndex[j];
	    double value=packed[j];
	    std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
		     <<") ";
	  }
	  std::cout<<std::endl;
	  //abort();
#endif
	}
      }
    } else {
      // not basic
#if CGL_DEBUG>1
      {
	// put column into array
	array.setVector(columnLength[iColumn],row+columnStart[iColumn],
			columnElements+columnStart[iColumn]);
	// get column in tableau
#ifdef CLP_OSL
	if (!alternateFactorization_)
#endif
	  factorization.updateColumn ( &work, &array );
#ifdef CLP_OSL
	else
	  factorization2->updateColumn ( &work, &array );
#endif
	int numberInArray=array.getNumElements();
	printf("non-basic %d\n",iColumn);
	for (int j=0;j<numberInArray;j++) {
	  int indexValue=arrayRows[j];
	  double value=arrayElements[indexValue];
	  if (fabs(value)>1.0e-6) {
	    printf("%d %g\n",indexValue,value);
	  }
	}
      }
#endif
    }
  }
#ifdef CLP_OSL
  delete factorization2;
#endif

  delete [] rowActivity;
  delete [] swap;
  delete [] rowType;
  delete [] packed;
  delete [] rowIsBasic;
  delete [] columnIsBasic;
#ifdef MORE_GOMORY_CUTS
#if MORE_GOMORY_CUTS==1
  int numberInaccurate = secondaryCuts.sizeRowCuts();
#ifdef CLP_INVESTIGATE2
  int numberOrdinary = numberAdded-numberInaccurate;
  if (!info.inTree&&(infoOptions&512)==0) 
    printf("Gomory has %d ordinary and %d less accurate cuts(%d els)\n",
	   numberOrdinary,numberInaccurate,saveTotalEls-nTotalEls);
#endif
#elif MORE_GOMORY_CUTS==2
  int numberLong = longCuts.sizeRowCuts();
#ifdef CLP_INVESTIGATE2
  int numberOrdinary = numberAdded-numberLong;
  if (!info.inTree&&(infoOptions&512)==0) 
    printf("Gomory has %d ordinary and %d long cuts(%d els)\n",
	   numberOrdinary,numberLong,saveTotalEls-nTotalEls);
#endif
#elif MORE_GOMORY_CUTS==3
  int numberLong = longCuts.sizeRowCuts();
  int numberInaccurate = secondaryCuts.sizeRowCuts();
#ifdef CLP_INVESTIGATE2
  int numberOrdinary = numberAdded-numberLong-numberInaccurate;
  if (!info.inTree&&(infoOptions&512)==0) 
    printf("Gomory has %d ordinary, %d long and %d less accurate cuts(%d els)\n",
	   numberOrdinary,numberLong,numberInaccurate,saveTotalEls-nTotalEls);
#endif
#endif
  if (doSorted&&limit<numberColumns) {
    // Just half
    nTotalEls -= saveTotalEls/2;
#if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
    while (nTotalEls>0) {
      for (int i=0;i<numberLong;i++) {
	nTotalEls -= longCuts.rowCutPtr(i)->row().getNumElements();
	cs.insert(longCuts.rowCut(i));
	numberAdded ++;
	if (nTotalEls<=0)
	  break;
      }
      break;
    }
#endif
#if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
    while (nTotalEls>0) {
      for (int i=0;i<numberInaccurate;i++) {
	nTotalEls -= secondaryCuts.rowCutPtr(i)->row().getNumElements();
	cs.insert(secondaryCuts.rowCut(i));
	numberAdded ++;
	if (nTotalEls<=0)
	  break;
      }
      break;
    }
#endif
  }
#else
#ifdef CLP_INVESTIGATE2
  if (!info.inTree&&(infoOptions&512)==0) 
    printf("Gomory added %d cuts(%d els)\n",numberAdded,saveTotalEls-nTotalEls);
#endif
#endif
  return numberAdded;
}
// Limit stuff
void CglGomory::setLimit(int limit)
{
  if (limit>=0)
    limit_=limit; 
}
int CglGomory::getLimit() const
{
  return limit_;
}
// Limit stuff at root
void CglGomory::setLimitAtRoot(int limit)
{
  if (limit>=0)
    limitAtRoot_=limit;
}
int CglGomory::getLimitAtRoot() const
{
  return limitAtRoot_;
}
// Return maximum length of cut in tree
int 
CglGomory::maximumLengthOfCutInTree() const
{
  if (limit_)
    return limit_;
  else
    return dynamicLimitInTree_;
}

// Away stuff
void CglGomory::setAway(double value)
{
  if (value>0.0&&value<=0.5)
    away_=value;
}
double CglGomory::getAway() const
{
  return away_;
}

// Away stuff at root
void CglGomory::setAwayAtRoot(double value)
{
  if (value>0.0&&value<=0.5)
    awayAtRoot_=value;
}
double CglGomory::getAwayAtRoot() const
{
  return awayAtRoot_;
}

// ConditionNumberMultiplier stuff
void CglGomory::setConditionNumberMultiplier(double value)
{
  if (value>=0.0)
    conditionNumberMultiplier_=value;
}
double CglGomory::getConditionNumberMultiplier() const
{
  return conditionNumberMultiplier_;
}

// LargestFactorMultiplier stuff
void CglGomory::setLargestFactorMultiplier(double value)
{
  if (value>=0.0)
    largestFactorMultiplier_=value;
}
double CglGomory::getLargestFactorMultiplier() const
{
  return largestFactorMultiplier_;
}

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
CglGomory::CglGomory ()
:
CglCutGenerator(),
away_(0.05),
awayAtRoot_(0.05),
conditionNumberMultiplier_(1.0e-18),
largestFactorMultiplier_(1.0e-13),
originalSolver_(NULL),
limit_(50),
limitAtRoot_(0),
dynamicLimitInTree_(-1),
alternateFactorization_(0),
gomoryType_(0)
{

}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
CglGomory::CglGomory (const CglGomory & source) :
  CglCutGenerator(source),
  away_(source.away_),
  awayAtRoot_(source.awayAtRoot_),
  conditionNumberMultiplier_(source.conditionNumberMultiplier_),
  largestFactorMultiplier_(source.largestFactorMultiplier_),
  originalSolver_(NULL),
  limit_(source.limit_),
  limitAtRoot_(source.limitAtRoot_),
  dynamicLimitInTree_(source.dynamicLimitInTree_),
  alternateFactorization_(source.alternateFactorization_),
  gomoryType_(source.gomoryType_)
{ 
  if (source.originalSolver_)
    originalSolver_ = source.originalSolver_->clone();
}

//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
CglCutGenerator *
CglGomory::clone() const
{
  return new CglGomory(*this);
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
CglGomory::~CglGomory ()
{
  delete originalSolver_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
CglGomory &
CglGomory::operator=(const CglGomory& rhs)
{
  if (this != &rhs) {
    CglCutGenerator::operator=(rhs);
    away_=rhs.away_;
    awayAtRoot_=rhs.awayAtRoot_;
    conditionNumberMultiplier_ = rhs.conditionNumberMultiplier_;
    largestFactorMultiplier_ = rhs.largestFactorMultiplier_;
    limit_=rhs.limit_;
    limitAtRoot_=rhs.limitAtRoot_;
    dynamicLimitInTree_ = rhs.dynamicLimitInTree_;
    alternateFactorization_=rhs.alternateFactorization_; 
    gomoryType_ = rhs.gomoryType_;
    delete originalSolver_;
    if (rhs.originalSolver_)
      originalSolver_ = rhs.originalSolver_->clone();
    else
      originalSolver_=NULL;
  }
  return *this;
}
// This can be used to refresh any information
void 
CglGomory::refreshSolver(OsiSolverInterface * solver)
{
  int numberColumns=solver->getNumCols(); 
  const double * colUpper = solver->getColUpper();
  const double * colLower = solver->getColLower();
  canDoGlobalCuts_ = true;
  if (originalSolver_) {
    delete originalSolver_;
    originalSolver_ = solver->clone();
  }
  for (int i=0;i<numberColumns;i++) {
    if (solver->isInteger(i)) {
      if (colUpper[i]>colLower[i]+1.0) {
	canDoGlobalCuts_ = false;
	break;
      }
    }
  }
}
// Pass in a copy of original solver (clone it)
void 
CglGomory::passInOriginalSolver(OsiSolverInterface * solver)
{
  delete originalSolver_;
  if (solver) {
    if (!gomoryType_)
      gomoryType_=1;
    originalSolver_ = solver->clone();
  } else {
    gomoryType_=0;
    originalSolver_=NULL;
  }
}
// Does actual work - returns number of cuts
int
CglGomory::generateCuts( const OsiRowCutDebugger * debugger, 
                         OsiCuts & cs,
                         const CoinPackedMatrix & columnCopy,
                         const double * colsol,
                         const double * colLower, const double * colUpper,
                         const double * rowLower, const double * rowUpper,
			 const char * intVar,
                         const CoinWarmStartBasis* warm,
                         const CglTreeInfo info)
{
  CoinPackedMatrix rowCopy;
  rowCopy.reverseOrderedCopyOf(columnCopy);
  return generateCuts( debugger, cs, columnCopy, rowCopy,
		       colsol, colLower, colUpper,
		       rowLower, rowUpper, intVar, warm, info);
}
// Create C++ lines to get to current state
std::string
CglGomory::generateCpp( FILE * fp) 
{
  CglGomory other;
  fprintf(fp,"0#include \"CglGomory.hpp\"\n");
  fprintf(fp,"3  CglGomory gomory;\n");
  if (limit_!=other.limit_)
    fprintf(fp,"3  gomory.setLimit(%d);\n",limit_);
  else
    fprintf(fp,"4  gomory.setLimit(%d);\n",limit_);
  if (limitAtRoot_!=other.limitAtRoot_)
    fprintf(fp,"3  gomory.setLimitAtRoot(%d);\n",limitAtRoot_);
  else
    fprintf(fp,"4  gomory.setLimitAtRoot(%d);\n",limitAtRoot_);
  if (away_!=other.away_)
    fprintf(fp,"3  gomory.setAway(%g);\n",away_);
  else
    fprintf(fp,"4  gomory.setAway(%g);\n",away_);
  if (awayAtRoot_!=other.awayAtRoot_)
    fprintf(fp,"3  gomory.setAwayAtRoot(%g);\n",awayAtRoot_);
  else
    fprintf(fp,"4  gomory.setAwayAtRoot(%g);\n",awayAtRoot_);
  if (getAggressiveness()!=other.getAggressiveness())
    fprintf(fp,"3  gomory.setAggressiveness(%d);\n",getAggressiveness());
  else
    fprintf(fp,"4  gomory.setAggressiveness(%d);\n",getAggressiveness());
  return "gomory";
}