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

#include "CoinPragma.hpp"

#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "ClpHelperFunctions.hpp"
#include "ClpSimplexNonlinear.hpp"
#include "ClpFactorization.hpp"
#include "ClpNonLinearCost.hpp"
#include "ClpLinearObjective.hpp"
#include "ClpConstraint.hpp"
#include "ClpQuadraticObjective.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpPrimalColumnPivot.hpp"
#include "ClpMessage.hpp"
#include "ClpEventHandler.hpp"
#include <cfloat>
#include <cassert>
#include <string>
#include <stdio.h>
#include <iostream>
#ifndef NDEBUG
#define CLP_DEBUG 1
#endif
// primal 
int ClpSimplexNonlinear::primal ()
{

  int ifValuesPass=1;
  algorithm_ = +3;

  // save data
  ClpDataSave data = saveData();
  matrix_->refresh(this); // make sure matrix okay

  // Save objective
  ClpObjective * saveObjective=NULL;
  if (objective_->type()>1) {
    // expand to full if quadratic
#ifndef NO_RTTI
    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
#else
    ClpQuadraticObjective * quadraticObj = NULL;
    if (objective_->type()==2)
      quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
#endif
    // for moment only if no scaling
    // May be faster if switched off - but can't see why
    if (!quadraticObj->fullMatrix()&&!rowScale_&&objectiveScale_==1.0) {
      saveObjective = objective_;
      objective_=new ClpQuadraticObjective(*quadraticObj,1);
    }
  }
  double bestObjectiveWhenFlagged=COIN_DBL_MAX;
  int pivotMode=15;
  //pivotMode=20;

  // initialize - maybe values pass and algorithm_ is +1
  // true does something (? perturbs)
  if (!startup(true)) {
    
    // Set average theta
    nonLinearCost_->setAverageTheta(1.0e3);
    int lastCleaned=0; // last time objective or bounds cleaned up
    
    // Say no pivot has occurred (for steepest edge and updates)
    pivotRow_=-2;
    
    // This says whether to restore things etc
    int factorType=0;
    // Start check for cycles
    progress_->startCheck();
    /*
      Status of problem:
      0 - optimal
      1 - infeasible
      2 - unbounded
      -1 - iterating
      -2 - factorization wanted
      -3 - redo checking without factorization
      -4 - looks infeasible
      -5 - looks unbounded
    */
    while (problemStatus_<0) {
      int iRow,iColumn;
      // clear
      for (iRow=0;iRow<4;iRow++) {
	rowArray_[iRow]->clear();
      }    
      
      for (iColumn=0;iColumn<2;iColumn++) {
	columnArray_[iColumn]->clear();
      }    
      
      // give matrix (and model costs and bounds a chance to be
      // refreshed (normally null)
      matrix_->refresh(this);
      // If getting nowhere - why not give it a kick
      // If we have done no iterations - special
      if (lastGoodIteration_==numberIterations_&&factorType)
	factorType=3;
      
      // may factorize, checks if problem finished
      if (objective_->type()>1&&lastFlaggedIteration_>=0&&
	  numberIterations_>lastFlaggedIteration_+507) {
	unflag();
	lastFlaggedIteration_=numberIterations_;
	if (pivotMode>=10) {
	  pivotMode--;
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32) 
	    printf("pivot mode now %d\n",pivotMode);
#endif
	  if (pivotMode==9) 
	    pivotMode=0;	// switch off fast attempt
	}
      }
      statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,
			      bestObjectiveWhenFlagged);
      
      // Say good factorization
      factorType=1;
      
      // Say no pivot has occurred (for steepest edge and updates)
      pivotRow_=-2;

      // exit if victory declared
      if (problemStatus_>=0)
	break;
      
      // test for maximum iterations
      if (hitMaximumIterations()||(ifValuesPass==2&&firstFree_<0)) {
	problemStatus_=3;
	break;
      }

      if (firstFree_<0) {
	if (ifValuesPass) {
	  // end of values pass
	  ifValuesPass=0;
	  int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
	  if (status>=0) {
	    problemStatus_=5;
	    secondaryStatus_=ClpEventHandler::endOfValuesPass;
	    break;
	  }
	}
      }
      // Check event
      {
	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
	if (status>=0) {
	  problemStatus_=5;
	  secondaryStatus_=ClpEventHandler::endOfFactorization;
	  break;
	}
      }
      // Iterate
      whileIterating(pivotMode);
    }
  }
  // if infeasible get real values
  if (problemStatus_==1) {
    infeasibilityCost_=0.0;
    createRim(1+4);
    nonLinearCost_->checkInfeasibilities(0.0);
    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
    // and get good feasible duals
    computeDuals(NULL);
  }
  // correct objective value
  if (numberColumns_)
    objectiveValue_ = nonLinearCost_->feasibleCost()+objective_->nonlinearOffset();
  objectiveValue_ /= (objectiveScale_*rhsScale_);
  // clean up
  unflag();
  finish();
  restoreData(data);
  // restore objective if full
  if (saveObjective) {
    delete objective_;
    objective_=saveObjective;
  }
  return problemStatus_;
}
/*  Refactorizes if necessary 
    Checks if finished.  Updates status.
    lastCleaned refers to iteration at which some objective/feasibility
    cleaning too place.
    
    type - 0 initial so set up save arrays etc
    - 1 normal -if good update save
    - 2 restoring from saved 
*/
void 
ClpSimplexNonlinear::statusOfProblemInPrimal(int & lastCleaned, int type,
						 ClpSimplexProgress * progress,
						 bool doFactorization,
						 double & bestObjectiveWhenFlagged)
{
  int dummy; // for use in generalExpanded
  if (type==2) {
    // trouble - restore solution
    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
	   numberRows_*sizeof(double));
    memcpy(columnActivityWork_,savedSolution_ ,
	   numberColumns_*sizeof(double));
    // restore extra stuff
    matrix_->generalExpanded(this,6,dummy);
    forceFactorization_=1; // a bit drastic but ..
    pivotRow_=-1; // say no weights update
    changeMade_++; // say change made
  }
  int saveThreshold = factorization_->sparseThreshold();
  int tentativeStatus = problemStatus_;
  int numberThrownOut=1; // to loop round on bad factorization in values pass
  while (numberThrownOut) {
    if (problemStatus_>-3||problemStatus_==-4) {
      // factorize
      // later on we will need to recover from singularities
      // also we could skip if first time
      // do weights
      // This may save pivotRow_ for use 
      if (doFactorization)
	primalColumnPivot_->saveWeights(this,1);
      
      if (type&&doFactorization) {
	// is factorization okay?
	int factorStatus = internalFactorize(1);
	if (factorStatus) {
	  if (type!=1||largestPrimalError_>1.0e3
	      ||largestDualError_>1.0e3) {
	    // was ||largestDualError_>1.0e3||objective_->type()>1) {
	    // switch off dense
	    int saveDense = factorization_->denseThreshold();
	    factorization_->setDenseThreshold(0);
	    // make sure will do safe factorization
	    pivotVariable_[0]=-1;
	    internalFactorize(2);
	    factorization_->setDenseThreshold(saveDense);
	    // Go to safe 
	    factorization_->pivotTolerance(0.99);
	    // restore extra stuff
	    matrix_->generalExpanded(this,6,dummy);
	  } else {
	    // no - restore previous basis
	    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
	    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
		   numberRows_*sizeof(double));
	    memcpy(columnActivityWork_,savedSolution_ ,
		   numberColumns_*sizeof(double));
	    // restore extra stuff
	    matrix_->generalExpanded(this,6,dummy);
	    matrix_->generalExpanded(this,5,dummy);
	    forceFactorization_=1; // a bit drastic but ..
	    type = 2;
	    // Go to safe 
	    factorization_->pivotTolerance(0.99);
	    if (internalFactorize(1)!=0)
	       largestPrimalError_=1.0e4; // force other type
	  }
	  // flag incoming
	  if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
	    setFlagged(sequenceIn_);
	    saveStatus_[sequenceIn_]=status_[sequenceIn_];
	  }
	  changeMade_++; // say change made
	}
      }
      if (problemStatus_!=-4)
	problemStatus_=-3;
    }
    // at this stage status is -3 or -5 if looks unbounded
    // get primal and dual solutions
    // put back original costs and then check
    createRim(4);
    // May need to do more if column generation
    dummy=4;
    matrix_->generalExpanded(this,9,dummy);
    numberThrownOut=gutsOfSolution(NULL,NULL,(firstFree_>=0));
    if (numberThrownOut) {
      problemStatus_=tentativeStatus;
      doFactorization=true;
    }
  }
  // Double check reduced costs if no action
  if (progress->lastIterationNumber(0)==numberIterations_) {
    if (primalColumnPivot_->looksOptimal()) {
      numberDualInfeasibilities_ = 0;
      sumDualInfeasibilities_ = 0.0;
    }
  }
  // Check if looping
  int loop;
  if (type!=2) 
    loop = progress->looping();
  else
    loop=-1;
  if (loop>=0) {
    if (!problemStatus_) {
      // declaring victory
      numberPrimalInfeasibilities_ = 0;
      sumPrimalInfeasibilities_ = 0.0;
    } else {
      problemStatus_ = loop; //exit if in loop 
      problemStatus_ = 10; // instead - try other algorithm
    }
    problemStatus_ = 10; // instead - try other algorithm
    return ;
  } else if (loop<-1) {
    // Is it time for drastic measures
    if (nonLinearCost_->numberInfeasibilities()&&progress->badTimes()>5&&
	progress->oddState()<10&&progress->oddState()>=0) {
      progress->newOddState();
      nonLinearCost_->zapCosts();
    }
    // something may have changed
    gutsOfSolution(NULL,NULL,true);
  }
  // If progress then reset costs
  if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
    createRim(4,false); // costs back
    delete nonLinearCost_;
    nonLinearCost_ = new ClpNonLinearCost(this);
    progress->endOddState();
    gutsOfSolution(NULL,NULL,true);
  }
  // Flag to say whether to go to dual to clean up
  bool goToDual=false;
  // really for free variables in
  //if((progressFlag_&2)!=0)
  //problemStatus_=-1;;
  progressFlag_ = 0; //reset progress flag

  handler_->message(CLP_SIMPLEX_STATUS,messages_)
    <<numberIterations_<<nonLinearCost_->feasibleReportCost();
  handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
    <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
  handler_->printing(sumDualInfeasibilities_>0.0)
    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
  handler_->printing(numberDualInfeasibilitiesWithoutFree_
		     <numberDualInfeasibilities_)
		       <<numberDualInfeasibilitiesWithoutFree_;
  handler_->message()<<CoinMessageEol;
  if (!primalFeasible()) {
    nonLinearCost_->checkInfeasibilities(primalTolerance_);
    gutsOfSolution(NULL,NULL,true);
    nonLinearCost_->checkInfeasibilities(primalTolerance_);
  }
  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
  if (trueInfeasibility>1.0) {
    // If infeasibility going up may change weights
    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
    if(progress->lastInfeasibility()<testValue) {
      if (infeasibilityCost_<1.0e14) {
	infeasibilityCost_ *= 1.5;
	if (handler_->logLevel()==63)
	  printf("increasing weight to %g\n",infeasibilityCost_);
	gutsOfSolution(NULL,NULL,true);
      }
    }
  }
  // we may wish to say it is optimal even if infeasible
  bool alwaysOptimal = (specialOptions_&1)!=0;
  // give code benefit of doubt
  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
    // say optimal (with these bounds etc)
    numberDualInfeasibilities_ = 0;
    sumDualInfeasibilities_ = 0.0;
    numberPrimalInfeasibilities_ = 0;
    sumPrimalInfeasibilities_ = 0.0;
  }
  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
  if (dualFeasible()||problemStatus_==-4) {
    // see if extra helps
    if (nonLinearCost_->numberInfeasibilities()&&
	 (nonLinearCost_->sumInfeasibilities()>1.0e-3||sumOfRelaxedPrimalInfeasibilities_)
	&&!alwaysOptimal) {
      //may need infeasiblity cost changed
      // we can see if we can construct a ray
      // make up a new objective
      double saveWeight = infeasibilityCost_;
      // save nonlinear cost as we are going to switch off costs
      ClpNonLinearCost * nonLinear = nonLinearCost_;
      // do twice to make sure Primal solution has settled
      // put non-basics to bounds in case tolerance moved
      // put back original costs
      createRim(4);
      nonLinearCost_->checkInfeasibilities(0.0);
      gutsOfSolution(NULL,NULL,true);

      infeasibilityCost_=1.0e100;
      // put back original costs
      createRim(4);
      nonLinearCost_->checkInfeasibilities(primalTolerance_);
      // may have fixed infeasibilities - double check
      if (nonLinearCost_->numberInfeasibilities()==0) {
	// carry on
	problemStatus_ = -1;
	infeasibilityCost_=saveWeight;
	nonLinearCost_->checkInfeasibilities(primalTolerance_);
      } else {
	nonLinearCost_=NULL;
	// scale
	int i;
	for (i=0;i<numberRows_+numberColumns_;i++) 
	  cost_[i] *= 1.0e-95;
	gutsOfSolution(NULL,NULL,false);
	nonLinearCost_=nonLinear;
	infeasibilityCost_=saveWeight;
	if ((infeasibilityCost_>=1.0e18||
	     numberDualInfeasibilities_==0)&&perturbation_==101) {
	  goToDual=unPerturb(); // stop any further perturbation
	  if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
	    goToDual=false;
	  nonLinearCost_->checkInfeasibilities(primalTolerance_);
	  numberDualInfeasibilities_=1; // carry on
	  problemStatus_=-1;
	}
	if (infeasibilityCost_>=1.0e20||
	    numberDualInfeasibilities_==0) {
	  // we are infeasible - use as ray
	  delete [] ray_;
	  ray_ = new double [numberRows_];
	  memcpy(ray_,dual_,numberRows_*sizeof(double));
	  // and get feasible duals
	  infeasibilityCost_=0.0;
	  createRim(4);
	  nonLinearCost_->checkInfeasibilities(primalTolerance_);
	  gutsOfSolution(NULL,NULL,true);
	  // so will exit
	  infeasibilityCost_=1.0e30;
	  // reset infeasibilities
	  sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
	  numberPrimalInfeasibilities_=
	    nonLinearCost_->numberInfeasibilities();
	}
	if (infeasibilityCost_<1.0e20) {
	  infeasibilityCost_ *= 5.0;
	  changeMade_++; // say change made
	  unflag();
	  handler_->message(CLP_PRIMAL_WEIGHT,messages_)
	    <<infeasibilityCost_
	    <<CoinMessageEol;
	  // put back original costs and then check
	  createRim(4);
	  nonLinearCost_->checkInfeasibilities(0.0);
	  gutsOfSolution(NULL,NULL,true);
	  problemStatus_=-1; //continue
	  goToDual=false;
	} else {
	  // say infeasible
	  problemStatus_ = 1;
	}
      }
    } else {
      // may be optimal
      if (perturbation_==101) {
	goToDual=unPerturb(); // stop any further perturbation
	lastCleaned=-1; // carry on
      }
      bool unflagged = unflag()!=0;
      if ( lastCleaned!=numberIterations_||unflagged) {
	handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
	  <<primalTolerance_
	  <<CoinMessageEol;
	double current = nonLinearCost_->feasibleReportCost();
	if (numberTimesOptimal_<4) {
	  if (bestObjectiveWhenFlagged<=current) {
	    numberTimesOptimal_++;
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("%d times optimal, current %g best %g\n",numberTimesOptimal_,
		     current,bestObjectiveWhenFlagged);
#endif
	  } else {
	    bestObjectiveWhenFlagged=current; 
	  }
	  changeMade_++; // say change made
	  if (numberTimesOptimal_==1) {
	    // better to have small tolerance even if slower
	    factorization_->zeroTolerance(1.0e-15);
	  }
	  lastCleaned=numberIterations_;
	  if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
	    handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
	      <<CoinMessageEol;
	  double oldTolerance = primalTolerance_;
	  primalTolerance_=dblParam_[ClpPrimalTolerance];
	  // put back original costs and then check
	  createRim(4);
	  nonLinearCost_->checkInfeasibilities(oldTolerance);
	  gutsOfSolution(NULL,NULL,true);
	  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
	      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
	    // say optimal (with these bounds etc)
	    numberDualInfeasibilities_ = 0;
	    sumDualInfeasibilities_ = 0.0;
	    numberPrimalInfeasibilities_ = 0;
	    sumPrimalInfeasibilities_ = 0.0;
	  }
	  if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
	    problemStatus_=0;
	  else
	    problemStatus_ = -1;
	} else {
	  problemStatus_=0; // optimal
	  if (lastCleaned<numberIterations_) {
	    handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
	      <<CoinMessageEol;
	  }
	}
      } else {
	problemStatus_=0; // optimal
      }
    }
  } else {
    // see if looks unbounded
    if (problemStatus_==-5) {
      if (nonLinearCost_->numberInfeasibilities()) {
	if (infeasibilityCost_>1.0e18&&perturbation_==101) {
	  // back off weight
	  infeasibilityCost_ = 1.0e13;
	  unPerturb(); // stop any further perturbation
	}
	//we need infeasiblity cost changed
	if (infeasibilityCost_<1.0e20) {
	  infeasibilityCost_ *= 5.0;
	  changeMade_++; // say change made
	  unflag();
	  handler_->message(CLP_PRIMAL_WEIGHT,messages_)
	    <<infeasibilityCost_
	    <<CoinMessageEol;
	  // put back original costs and then check
	  createRim(4);
	  gutsOfSolution(NULL,NULL,true);
	  problemStatus_=-1; //continue
	} else {
	  // say unbounded
	  problemStatus_ = 2;
	}
      } else {
	// say unbounded
	problemStatus_ = 2;
      } 
    } else {
      if(type==3&&problemStatus_!=-5)
	unflag(); // odd
      // carry on
      problemStatus_ = -1;
    }
  }
  // save extra stuff
  matrix_->generalExpanded(this,5,dummy);
  if (type==0||type==1) {
    if (type!=1||!saveStatus_) {
      // create save arrays
      delete [] saveStatus_;
      delete [] savedSolution_;
      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
      savedSolution_ = new double [numberRows_+numberColumns_];
    }
    // save arrays
    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
	   numberRows_*sizeof(double));
    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
  }
  if (doFactorization) {
    // restore weights (if saved) - also recompute infeasibility list
    if (tentativeStatus>-3) 
      primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
    else
      primalColumnPivot_->saveWeights(this,3);
    if (saveThreshold) {
      // use default at present
      factorization_->sparseThreshold(0);
      factorization_->goSparse();
    }
  }
  if (problemStatus_<0&&!changeMade_) {
    problemStatus_=4; // unknown
  }
  lastGoodIteration_ = numberIterations_;
  //if (goToDual) 
  //problemStatus_=10; // try dual
  // Allow matrices to be sorted etc
  int fake=-999; // signal sort
  matrix_->correctSequence(this,fake,fake);
}
/*
  Reasons to come out:
  -1 iterations etc
  -2 inaccuracy 
  -3 slight inaccuracy (and done iterations)
  -4 end of values pass and done iterations
  +0 looks optimal (might be infeasible - but we will investigate)
  +2 looks unbounded
  +3 max iterations 
*/
int
ClpSimplexNonlinear::whileIterating(int & pivotMode)
{
  // Say if values pass
  //int ifValuesPass=(firstFree_>=0) ? 1 : 0;
  int ifValuesPass=1;
  int returnCode=-1;
  // status stays at -1 while iterating, >=0 finished, -2 to invert
  // status -3 to go to top without an invert
  int numberInterior=0;
  int nextUnflag=10;
  int nextUnflagIteration=numberIterations_+10;
  // get two arrays
  double * array1 = new double[2*(numberRows_+numberColumns_)];
  double solutionError=-1.0;
  while (problemStatus_==-1) {
    int result;
    rowArray_[1]->clear();
    //#define CLP_DEBUG
#if CLP_DEBUG > 1
    rowArray_[0]->checkClear();
    rowArray_[1]->checkClear();
    rowArray_[2]->checkClear();
    rowArray_[3]->checkClear();
    columnArray_[0]->checkClear();
#endif
    if (numberInterior>=5) {
      // this can go when we have better minimization
      if (pivotMode<10)
	pivotMode=1;
      unflag();
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32) 
	    printf("interior unflag\n");
#endif
      numberInterior=0;
      nextUnflag=10;
      nextUnflagIteration=numberIterations_+10;
    } else {
      if (numberInterior>nextUnflag&&
	  numberIterations_>nextUnflagIteration) {
	nextUnflagIteration=numberIterations_+10;
	nextUnflag += 10;
	unflag();
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("unflagging as interior\n");
#endif
      }
    } 
    pivotRow_=-1;
    result = pivotColumn(rowArray_[3],rowArray_[0],
			 columnArray_[0],rowArray_[1],pivotMode,solutionError,
			 array1);
    if (result) {
      if (result==3) 
	break; // null vector not accurate
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) {
	double currentObj;
	double thetaObj;
	double predictedObj;
	objective_->stepLength(this,solution_,solution_,0.0,
			       currentObj,predictedObj,thetaObj);
	printf("obj %g after interior move\n",currentObj);
      }
#endif
      // just move and try again
      if (pivotMode<10) {
	pivotMode=result-1;
	numberInterior++;
      }
      continue;
    } else {
      if (pivotMode<10) {
	if (theta_>0.001)
	  pivotMode=0;
	else if (pivotMode==2)
	  pivotMode=1;
      }
      numberInterior=0;
      nextUnflag=10;
      nextUnflagIteration=numberIterations_+10;
    }
    sequenceOut_=-1;
    rowArray_[1]->clear();
    if (sequenceIn_>=0) {
      // we found a pivot column
      assert (!flagged(sequenceIn_));
#ifdef CLP_DEBUG
      if ((handler_->logLevel()&32)) {
	char x = isColumn(sequenceIn_) ? 'C' :'R';
	std::cout<<"pivot column "<<
	  x<<sequenceWithin(sequenceIn_)<<std::endl;
      }
#endif
      // do second half of iteration
      if (pivotRow_<0&&theta_<1.0e-8) {
	assert (sequenceIn_>=0);
	returnCode = pivotResult(ifValuesPass);
      } else {
	// specialized code
	returnCode = pivotNonlinearResult();
	//printf("odd pivrow %d\n",sequenceOut_);
	if (sequenceOut_>=0&&theta_<1.0e-5) {
	  if (getStatus(sequenceOut_)!=isFixed) {
	    if (getStatus(sequenceOut_)==atUpperBound)
	      solution_[sequenceOut_]=upper_[sequenceOut_];
	    else if (getStatus(sequenceOut_)==atLowerBound)
	      solution_[sequenceOut_]=lower_[sequenceOut_];
	    setFlagged(sequenceOut_);
	  }
	}
      }
      if (returnCode<-1&&returnCode>-5) {
	problemStatus_=-2; // 
      } else if (returnCode==-5) {
	// something flagged - continue;
      } else if (returnCode==2) {
	problemStatus_=-5; // looks unbounded
      } else if (returnCode==4) {
	problemStatus_=-2; // looks unbounded but has iterated
      } else if (returnCode!=-1) {
	assert(returnCode==3);
	problemStatus_=3;
      }
    } else {
      // no pivot column
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32)
	printf("** no column pivot\n");
#endif
      if (pivotMode<10) {
	// looks optimal
	primalColumnPivot_->setLooksOptimal(true);
      } else {
	pivotMode--;
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("pivot mode now %d\n",pivotMode);
#endif
	if (pivotMode==9) 
	  pivotMode=0;	// switch off fast attempt
	unflag();
      }
      if (nonLinearCost_->numberInfeasibilities())
	problemStatus_=-4; // might be infeasible 
      returnCode=0;
      break;
    }
  }
  delete [] array1;
  return returnCode;
}
// Creates direction vector
void
ClpSimplexNonlinear::directionVector (CoinIndexedVector * vectorArray,
				      CoinIndexedVector * spare1, CoinIndexedVector * spare2,
				      int pivotMode2,
				      double & normFlagged,double & normUnflagged,
				      int & numberNonBasic)
{
#if CLP_DEBUG > 1
  vectorArray->checkClear();
  spare1->checkClear();
  spare2->checkClear();
#endif
  double *array = vectorArray->denseVector();
  int * index = vectorArray->getIndices();
  int number=0;
  sequenceIn_=-1;
  normFlagged=0.0;
  normUnflagged=1.0;
  double dualTolerance2 = CoinMin(1.0e-8,1.0e-2*dualTolerance_);
  double dualTolerance3 = CoinMin(1.0e-2,1.0e3*dualTolerance_);
  if (!numberNonBasic) {
    //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
    //printf("infeasible\n");
    if (!pivotMode2||pivotMode2>=10) {
      normUnflagged=0.0;
      double bestDj=0.0;
      double bestSuper=0.0;
      double sumSuper=0.0;
      sequenceIn_=-1;
      int nSuper=0;
      for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
	array[iSequence]=0.0;
	if (flagged(iSequence)) {
	  // accumulate  norm
	  switch(getStatus(iSequence)) {
	    
	  case basic:
	  case ClpSimplex::isFixed:
	    break;
	  case atUpperBound:
	    if (dj_[iSequence]>dualTolerance3) 
	      normFlagged += dj_[iSequence]*dj_[iSequence];
	    break;
	  case atLowerBound:
	    if (dj_[iSequence]<-dualTolerance3) 
	      normFlagged += dj_[iSequence]*dj_[iSequence];
	    break;
	  case isFree:
	  case superBasic:
	    if (fabs(dj_[iSequence])>dualTolerance3) 
	      normFlagged += dj_[iSequence]*dj_[iSequence];
	    break;
	  }
	  continue;
	}
	switch(getStatus(iSequence)) {
	  
	case basic:
	case ClpSimplex::isFixed:
	  break;
	case atUpperBound:
	  if (dj_[iSequence]>dualTolerance_) {
	    if (dj_[iSequence]>dualTolerance3) 
	      normUnflagged += dj_[iSequence]*dj_[iSequence];
	    if (pivotMode2<10) {
	      array[iSequence]=-dj_[iSequence];
	      index[number++]=iSequence;
	    } else {
	      if (dj_[iSequence]>bestDj) {
		bestDj=dj_[iSequence];
		sequenceIn_=iSequence;
	      }
	    }
	  }
	  break;
	case atLowerBound:
	  if (dj_[iSequence]<-dualTolerance_) {
	    if (dj_[iSequence]<-dualTolerance3) 
	      normUnflagged += dj_[iSequence]*dj_[iSequence];
	    if (pivotMode2<10) {
	      array[iSequence]=-dj_[iSequence];
	      index[number++]=iSequence;
	    } else {
	      if (-dj_[iSequence]>bestDj) {
		bestDj=-dj_[iSequence];
		sequenceIn_=iSequence;
	      }
	    }
	  }
	  break;
	case isFree:
	case superBasic:
	  //#define ALLSUPER
#define NOSUPER
#ifndef ALLSUPER
	  if (fabs(dj_[iSequence])>dualTolerance_) {
	    if (fabs(dj_[iSequence])>dualTolerance3) 
	      normUnflagged += dj_[iSequence]*dj_[iSequence];
	    nSuper++;
	    bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
	    sumSuper += fabs(dj_[iSequence]);
	  }
	  if (fabs(dj_[iSequence])>dualTolerance2) {
	    array[iSequence]=-dj_[iSequence];
	    index[number++]=iSequence;
	  }
#else
	  array[iSequence]=-dj_[iSequence];
	  index[number++]=iSequence;
	  if (pivotMode2>=10)
	    bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
#endif
	  break;
	}
      }
#ifdef NOSUPER
      // redo
      bestSuper=sumSuper;
      if(sequenceIn_>=0&&bestDj>bestSuper) {
	int j;
	// get rid of superbasics
	for (j=0;j<number;j++) {
	  int iSequence = index[j];
	  array[iSequence]=0.0;
	}
	number=0;
	array[sequenceIn_]=-dj_[sequenceIn_];
	index[number++]=sequenceIn_;
      } else {
	sequenceIn_=-1;
      }
#else
      if (pivotMode2>=10||!nSuper) {
	bool takeBest=true;
	if (pivotMode2==100&&nSuper>1)
	  takeBest=false;
	if(sequenceIn_>=0&&takeBest) {
	  if (fabs(dj_[sequenceIn_])>bestSuper) {
	    array[sequenceIn_]=-dj_[sequenceIn_];
	    index[number++]=sequenceIn_;
	  } else {
	    sequenceIn_=-1;
	  }
	} else {
	  sequenceIn_=-1;
	}
      }
#endif
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) {
	if (sequenceIn_>=0)
	  printf("%d superBasic, chosen %d - dj %g\n",nSuper,sequenceIn_,
		 dj_[sequenceIn_]);
	else
	  printf("%d superBasic - none chosen\n",nSuper);
      }
#endif
    } else {
      double bestDj=0.0;
      double saveDj=0.0;
      if (sequenceOut_>=0) {
	saveDj=dj_[sequenceOut_];
	dj_[sequenceOut_]=0.0;
	switch(getStatus(sequenceOut_)) {
	  
	case basic:
	  sequenceOut_=-1;
	case ClpSimplex::isFixed:
	  break;
	case atUpperBound:
	  if (dj_[sequenceOut_]>dualTolerance_) {
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("after pivot out %d values %g %g %g, dj %g\n",
		     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
		     upper_[sequenceOut_],dj_[sequenceOut_]);
#endif
	  }
	  break;
	case atLowerBound:
	  if (dj_[sequenceOut_]<-dualTolerance_) {
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("after pivot out %d values %g %g %g, dj %g\n",
		     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
		     upper_[sequenceOut_],dj_[sequenceOut_]);
#endif
	  }
	  break;
	case isFree:
	case superBasic:
	  if (dj_[sequenceOut_]>dualTolerance_) {
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("after pivot out %d values %g %g %g, dj %g\n",
		     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
		     upper_[sequenceOut_],dj_[sequenceOut_]);
#endif
	  } else if (dj_[sequenceOut_]<-dualTolerance_) {
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("after pivot out %d values %g %g %g, dj %g\n",
		     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
		     upper_[sequenceOut_],dj_[sequenceOut_]);
#endif
	  }
	  break;
	}
      }
      // Go for dj
      pivotMode2=3;
      for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
	array[iSequence]=0.0;
	if (flagged(iSequence))
	  continue;
	switch(getStatus(iSequence)) {
	  
	case basic:
	case ClpSimplex::isFixed:
	  break;
	case atUpperBound:
	  if (dj_[iSequence]>dualTolerance_) {
	    double distance = CoinMin(1.0e-2,solution_[iSequence]-lower_[iSequence]);
	    double merit = distance*dj_[iSequence];
	    if (pivotMode2==1)
	      merit *= 1.0e-20; // discourage
	    if (pivotMode2==3)
	      merit = fabs(dj_[iSequence]);
	    if (merit>bestDj) {
	      sequenceIn_=iSequence;
	      bestDj=merit;
	    }
	  }
	  break;
	case atLowerBound:
	  if (dj_[iSequence]<-dualTolerance_) {
	    double distance = CoinMin(1.0e-2,upper_[iSequence]-solution_[iSequence]);
	    double merit = -distance*dj_[iSequence];
	    if (pivotMode2==1)
	      merit *= 1.0e-20; // discourage
	    if (pivotMode2==3)
	      merit = fabs(dj_[iSequence]);
	    if (merit>bestDj) {
	      sequenceIn_=iSequence;
	      bestDj=merit;
	    }
	  }
	  break;
	case isFree:
	case superBasic:
	  if (dj_[iSequence]>dualTolerance_) {
	    double distance = CoinMin(1.0e-2,solution_[iSequence]-lower_[iSequence]);
	    distance = CoinMin(solution_[iSequence]-lower_[iSequence],
			   upper_[iSequence]-solution_[iSequence]);
	    double merit = distance*dj_[iSequence];
	    if (pivotMode2==1)
	      merit=distance;
	    if (pivotMode2==3)
	      merit = fabs(dj_[iSequence]);
	    if (merit>bestDj) {
	      sequenceIn_=iSequence;
	      bestDj=merit;
	    }
	  } else if (dj_[iSequence]<-dualTolerance_) {
	    double distance = CoinMin(1.0e-2,upper_[iSequence]-solution_[iSequence]);
	    distance = CoinMin(solution_[iSequence]-lower_[iSequence],
			   upper_[iSequence]-solution_[iSequence]);
	    double merit = -distance*dj_[iSequence];
	    if (pivotMode2==1)
	      merit=distance;
	    if (pivotMode2==3)
	      merit = fabs(dj_[iSequence]);
	    if (merit>bestDj) {
	      sequenceIn_=iSequence;
	      bestDj=merit;
	    }
	  }
	  break;
	}
      }
      if (sequenceOut_>=0) {
	dj_[sequenceOut_]=saveDj;
	sequenceOut_=-1;
      }
      if (sequenceIn_>=0) {
	array[sequenceIn_]=-dj_[sequenceIn_];
	index[number++]=sequenceIn_;
      }
    }
    numberNonBasic = number;
  } else {
    // compute norms
    normUnflagged=0.0;
    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
      if (flagged(iSequence)) {
	// accumulate  norm
	switch(getStatus(iSequence)) {
	  
	case basic:
	case ClpSimplex::isFixed:
	  break;
	case atUpperBound:
	  if (dj_[iSequence]>dualTolerance_) 
	    normFlagged += dj_[iSequence]*dj_[iSequence];
	  break;
	case atLowerBound:
	  if (dj_[iSequence]<-dualTolerance_) 
	    normFlagged += dj_[iSequence]*dj_[iSequence];
	  break;
	case isFree:
	case superBasic:
	  if (fabs(dj_[iSequence])>dualTolerance_) 
	    normFlagged += dj_[iSequence]*dj_[iSequence];
	  break;
	}
      }
    }
    // re-use list
    number=0;
    int j;
    for (j=0;j<numberNonBasic;j++) {
      int iSequence = index[j];
      if (flagged(iSequence))
	continue;
      switch(getStatus(iSequence)) {
	  
      case basic:
      case ClpSimplex::isFixed:
	continue; //abort();
	break;
      case atUpperBound:
	if (dj_[iSequence]>dualTolerance_) {
	  number++;
	  normUnflagged += dj_[iSequence]*dj_[iSequence];
	}
	break;
      case atLowerBound:
	if (dj_[iSequence]<-dualTolerance_) {
	  number++;
	  normUnflagged += dj_[iSequence]*dj_[iSequence];
	}
	break;
      case isFree:
      case superBasic:
	if (fabs(dj_[iSequence])>dualTolerance_) {
	  number++;
	  normUnflagged += dj_[iSequence]*dj_[iSequence];
	}
	break;
      }
      array[iSequence]=-dj_[iSequence];
    }
    // switch to large
    normUnflagged=1.0;
    if (!number) {
      for (j=0;j<numberNonBasic;j++) {
	int iSequence = index[j];
	array[iSequence]=0.0;
      }
      numberNonBasic=0;
    }
    number=numberNonBasic;
  }
  if (number) {
    int j;
    // Now do basic ones - how do I compensate for small basic infeasibilities?
    int iRow;
    for (iRow=0;iRow<numberRows_;iRow++) {
      int iPivot=pivotVariable_[iRow];
      double value=0.0;
      if (solution_[iPivot]>upper_[iPivot]) {
	value = upper_[iPivot]-solution_[iPivot];
      } else if (solution_[iPivot]<lower_[iPivot]) {
	value = lower_[iPivot]-solution_[iPivot];
      }
      //if (value)
      //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
      //     upper_[iPivot]);
      //value=0.0;
      value *= -1.0e0;
      if (value) {
	array[iPivot]=value;
	index[number++]=iPivot;
      }
    }
    double * array2=spare1->denseVector();
    int * index2 = spare1->getIndices();
    int number2=0;
    times(-1.0,array,array2);
    array = array+numberColumns_;
    for (iRow=0;iRow<numberRows_;iRow++) {
      double value = array2[iRow] + array[iRow];
      if (value) {
	array2[iRow]=value;
	index2[number2++]=iRow;
      } else {
	array2[iRow]=0.0;
      }
    }
    array -= numberColumns_;
    spare1->setNumElements(number2);
    // Ftran
    factorization_->updateColumn(spare2,spare1);
    number2=spare1->getNumElements();
    for (j=0;j<number2;j++) {
      int iSequence = index2[j];
      double value = array2[iSequence];
      array2[iSequence]=0.0;
      if (value) {
	int iPivot=pivotVariable_[iSequence];
	double oldValue = array[iPivot];
	if (!oldValue) {
	  array[iPivot] = value;
	  index[number++]=iPivot;
	} else {
	  // something already there
	  array[iPivot] = value+oldValue;
	}
      }
    }
    spare1->setNumElements(0);
  }
  vectorArray->setNumElements(number);
}
#define MINTYPE 1
#if MINTYPE==2
static double 
innerProductIndexed(const double * region1, int size, const double * region2,const int * which)
{
  int i;
  double value=0.0;
  for (i=0;i<size;i++) {
    int j=which[i];
    value += region1[j]*region2[j];
  }
  return value;
}
#endif
/* 
   Row array and column array have direction
   Returns 0 - can do normal iteration (basis change)
   1 - no basis change
*/
int 
ClpSimplexNonlinear::pivotColumn(CoinIndexedVector * longArray,
				 CoinIndexedVector * rowArray,
				 CoinIndexedVector * columnArray,
				 CoinIndexedVector * spare,
				 int & pivotMode,
				 double & solutionError,
				 double * dArray)
{
  // say not optimal
  primalColumnPivot_->setLooksOptimal(false);
  double acceptablePivot=1.0e-10;
  int lastSequenceIn=-1;
  if (pivotMode&&pivotMode<10) {
    acceptablePivot=1.0e-6;
    if (factorization_->pivots())
      acceptablePivot=1.0e-5; // if we have iterated be more strict
  }
  double acceptableBasic=1.0e-7;
  
  int number=longArray->getNumElements();
  int numberTotal = numberRows_+numberColumns_;
  int bestSequence=-1;
  int bestBasicSequence=-1;
  double eps= 1.0e-1;
  eps=1.0e-6;
  double basicTheta = 1.0e30;
  double objTheta=0.0;
  bool finished = false;
  sequenceIn_=-1;
  int nPasses=0;
  int nTotalPasses=0;
  int nBigPasses=0;
  double djNorm0=0.0;
  double djNorm=0.0;
  double normFlagged=0.0;
  double normUnflagged=0.0;
  int localPivotMode=pivotMode;
  bool allFinished=false;
  bool justOne=false;
  int returnCode=1;
  double currentObj;
  double predictedObj;
  double thetaObj;
  objective_->stepLength(this,solution_,solution_,0.0,
					       currentObj,predictedObj,thetaObj);
  double saveObj=currentObj;
#if MINTYPE ==2
  // try Shanno's method
  //would be memory leak
  //double * saveY=new double[numberTotal];
  //double * saveS=new double[numberTotal];
  //double * saveY2=new double[numberTotal];
  //double * saveS2=new double[numberTotal];
  double saveY[100];
  double saveS[100];
  double saveY2[100];
  double saveS2[100];
  double zz[10000];
#endif
  double * dArray2 = dArray+numberTotal;
  // big big loop
  while (!allFinished) {
    double * work=longArray->denseVector();
    int * which=longArray->getIndices();
    allFinished=true;
    // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
    //#define SMALLTHETA1 1.0e-25
    //#define SMALLTHETA2 1.0e-25
#define SMALLTHETA1 1.0e-10
#define SMALLTHETA2 1.0e-10
#define CONJUGATE 2
#if CONJUGATE == 0
    int conjugate = 0;
#elif CONJUGATE == 1
    int conjugate = (pivotMode<10) ? MINTYPE : 0;
#elif CONJUGATE == 2
    int conjugate = MINTYPE;
#else
    int conjugate = MINTYPE;
#endif
    
    //if (pivotMode==1)
    //localPivotMode=11;;
#if CLP_DEBUG > 1
    longArray->checkClear();
#endif
    int numberNonBasic=0;
    int startLocalMode=-1;
    while (!finished) {
      double simpleObjective=COIN_DBL_MAX;
      returnCode=1;
      int iSequence;
      objective_->reducedGradient(this,dj_,false);
      // get direction vector in longArray
      longArray->clear();
      // take out comment to try slightly different idea
      if (nPasses+nTotalPasses>3000||nBigPasses>100) {
	if (factorization_->pivots())
	  returnCode=3;
	break;
      }
      if (!nPasses) {
	numberNonBasic=0;
	nBigPasses++;
      }
      // just do superbasic if in cleanup mode
      int local=localPivotMode;
      if (!local&&pivotMode>=10&&nBigPasses<10) {
	local=100;
      } else if (local==-1||nBigPasses>=10) {
	local=0;
	localPivotMode=0;
      } 
      if (justOne) {
	local=2;
	//local=100;
	justOne=false;
      }
      if (!nPasses)
	startLocalMode=local;
      directionVector(longArray,spare,rowArray,local,
		      normFlagged,normUnflagged,numberNonBasic);
      {
	// check null vector
	double * rhs = spare->denseVector();
#if CLP_DEBUG > 1
	spare->checkClear();
#endif
	int iRow;
	multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,rhs,0.0);
	matrix_->times(1.0,solution_,rhs,rowScale_,columnScale_);
	double largest=0.0;
	int iLargest=-1;
	for (iRow=0;iRow<numberRows_;iRow++) {
	  double value=fabs(rhs[iRow]);
	  rhs[iRow]=0.0;
	  if (value>largest) {
	    largest=value;
	    iLargest=iRow;
	  }
	}
#if CLP_DEBUG > 0
	if ((handler_->logLevel()&32)&&largest>1.0e-8) 
	  printf("largest rhs error %g on row %d\n",largest,iLargest);
#endif
	if (solutionError<0.0) {
	  solutionError=largest;
	} else if (largest>max(1.0e-8,1.0e2*solutionError)&&
		   factorization_->pivots()) {
	  longArray->clear();
	  pivotRow_ = -1;
	  theta_=0.0;
	  return 3;
	}
      }
      if (sequenceIn_>=0)
	lastSequenceIn=sequenceIn_;
#if MINTYPE!=2
      double djNormSave = djNorm;
#endif
      djNorm=0.0;
      int iIndex;
      for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
	int iSequence = which[iIndex];
	double alpha = work[iSequence];
	djNorm += alpha*alpha;
      }
      // go to conjugate gradient if necessary
      if (numberNonBasic&&localPivotMode>=10&&(nPasses>4||sequenceIn_<0)) {
	localPivotMode=0;
	nTotalPasses+= nPasses;
	nPasses=0;
      }
#if CONJUGATE == 2
      conjugate = (localPivotMode<10) ? MINTYPE : 0;
#endif
      //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
      //     nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
      if (!nPasses) {
	//printf("numberNon %d\n",numberNonBasic);
#if MINTYPE==2
	assert (numberNonBasic<100);
	memset(zz,0,numberNonBasic*numberNonBasic*sizeof(double));
	int put=0;
	for (int iVariable=0;iVariable<numberNonBasic;iVariable++) {
	  zz[put]=1.0;
	  put+= numberNonBasic+1;
	}
#endif
	djNorm0 = CoinMax(djNorm, 1.0e-20);
	memcpy(dArray,work,numberTotal*sizeof(double));
	memcpy(dArray2,work,numberTotal*sizeof(double));
	if (sequenceIn_>=0&&numberNonBasic==1) {
	  // see if simple move
	  double objTheta2 =objective_->stepLength(this,solution_,work,1.0e30,
					       currentObj,predictedObj,thetaObj);
	  rowArray->clear();
	  
	  // update the incoming column
	  unpackPacked(rowArray);
	  factorization_->updateColumnFT(spare,rowArray);
	  theta_=0.0;
	  double * work2=rowArray->denseVector();
	  int number=rowArray->getNumElements();
	  int * which2=rowArray->getIndices();
	  int iIndex;
	  bool easyMove=false;
	  double way;
	  if (dj_[sequenceIn_]>0.0)
	    way=-1.0;
	  else
	    way=1.0;
	  double largest=COIN_DBL_MAX;
	  int kPivot=-1;
	  for (iIndex=0;iIndex<number;iIndex++) {
	    int iRow = which2[iIndex];
	    double alpha = way*work2[iIndex];
	    int iPivot = pivotVariable_[iRow];
	    if (alpha<-1.0e-5) {
	      double distance = upper_[iPivot]-solution_[iPivot];
	      if (distance<-largest*alpha) {
		kPivot=iPivot;
		largest=CoinMax(0.0,-distance/alpha);
	      }
	      if (distance<-1.0e-12*alpha) {
		easyMove=true;
		break;
	      }
	    } else if (alpha>1.0e-5) {
	      double distance = solution_[iPivot]-lower_[iPivot];
	      if (distance<largest*alpha) {
		kPivot=iPivot;
		largest=CoinMax(0.0,distance/alpha);
	      }
	      if (distance<1.0e-12*alpha) {
		easyMove=true;
		break;
	      }
	    }
	  }
          // But largest has to match up with change
          assert (work[sequenceIn_]);
          largest /= fabs(work[sequenceIn_]);
	  if (largest<objTheta2) {
	    easyMove=true;
	  } else if (!easyMove) {
	    double objDrop = currentObj-predictedObj;
	    double th = objective_->stepLength(this,solution_,work,largest,
				   currentObj,predictedObj,simpleObjective);
	    simpleObjective = CoinMax(simpleObjective,predictedObj);
	    double easyDrop = currentObj-simpleObjective;
	    if (easyDrop>1.0e-8&&easyDrop>0.5*objDrop) {
	      easyMove=true;
#ifdef CLP_DEBUG
	      if (handler_->logLevel()&32) 
		printf("easy - obj drop %g, easy drop %g\n",objDrop,easyDrop);
#endif
	      if (easyDrop>objDrop) {
		// debug
		printf("****** th %g simple %g\n",th,simpleObjective);
		objective_->stepLength(this,solution_,work,1.0e30,
				   currentObj,predictedObj,simpleObjective);
		objective_->stepLength(this,solution_,work,largest,
				   currentObj,predictedObj,simpleObjective);
	      }
	    }
	  }
	  rowArray->clear();
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32)
	    printf("largest %g piv %d\n",largest,kPivot);
#endif
	  if (easyMove) {
	    valueIn_=solution_[sequenceIn_];
	    dualIn_=dj_[sequenceIn_];
	    lowerIn_=lower_[sequenceIn_];
	    upperIn_=upper_[sequenceIn_];
	    if (dualIn_>0.0)
	      directionIn_ = -1;
	    else 
	      directionIn_ = 1;
	    longArray->clear();
	    pivotRow_ = -1;
	    theta_=0.0;
	    return 0;
	  }
	}
      } else {
#if MINTYPE==1
	if (conjugate) {
	  double djNorm2 = djNorm;
	  if (numberNonBasic&&0) {
	    int iIndex;
	    djNorm2 = 0.0;
	    for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
	      int iSequence = which[iIndex];
	      double alpha = work[iSequence];
	      //djNorm2 += alpha*alpha;
	      double alpha2 = work[iSequence]-dArray2[iSequence];
	      djNorm2 += alpha*alpha2;
	    }
	    //printf("a %.18g b %.18g\n",djNorm,djNorm2);
	  }
	  djNorm=djNorm2;
	  double beta = djNorm2/djNormSave;
	  // reset beta every so often
	  //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1)
	  //beta=0.0;
	  //printf("current norm %g, old %g - beta %g\n",
	  //     djNorm,djNormSave,beta);
	  for (iSequence=0;iSequence<numberTotal;iSequence++) {
	    dArray[iSequence] = work[iSequence] + beta *dArray[iSequence];
	    dArray2[iSequence] = work[iSequence];
	  }
	} else {
	  for (iSequence=0;iSequence<numberTotal;iSequence++)
	    dArray[iSequence] = work[iSequence];
	}
#else
	int number2=numberNonBasic;
	if (number2) {
	  // pack down into dArray
	  int jLast=-1;
	  for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
	    int j=which[iSequence];
	    assert(j>jLast);
	    jLast=j;
	    double value=work[j];
	    dArray[iSequence]=-value;
	  }
	  // see whether to restart
	  // check signs - gradient
	  double g1=innerProduct(dArray,number2,dArray);
	  double g2=innerProduct(dArray,number2,saveY2);
	  // Get differences
	  for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
	    saveY2[iSequence] = dArray[iSequence]-saveY2[iSequence];
	    saveS2[iSequence] = solution_[iSequence]-saveS2[iSequence];
	  }
	  double g3=innerProduct(saveS2,number2,saveY2);
	  printf("inner %g\n",g3);
	  //assert(g3>0);
	  double zzz[10000];
	  int iVariable;
	  g2=1.0e50;// temp
	  if (fabs(g2)>=0.2*fabs(g1)) {
	    // restart
	    double delta = innerProduct(saveY2,number2,saveS2)/
	      innerProduct(saveY2,number2,saveY2);
	    delta=1.0;//temp
	    memset(zz,0,number2*sizeof(double));
	    int put=0;
	    for (iVariable=0;iVariable<number2;iVariable++) {
	      zz[put]=delta;
	      put+= number2+1;
	    }
	  } else {
	  }
	  memcpy(zzz,zz,number2*number2*sizeof(double));
	  double ww[100];
	  // get sk -Hkyk
	  for (iVariable=0;iVariable<number2;iVariable++) {
	    double value=0.0;
	    for (int jVariable=0;jVariable<number2;jVariable++) {
	      value += saveY2[jVariable]*zzz[iVariable+number2*jVariable];
	    }
	    ww[iVariable] = saveS2[iVariable]-value;
	  }
	  double ys = innerProduct(saveY2,number2,saveS2);
	  double multiplier1 = 1.0/ys;
	  double multiplier2= innerProduct(saveY2,number2,ww)/(ys*ys);
#if 1
	  // and second way
	  // Hy
	  double h[100];
	  for (iVariable=0;iVariable<number2;iVariable++) {
	    double value=0.0;
	    for (int jVariable=0;jVariable<number2;jVariable++) {
	      value += saveY2[jVariable]*zzz[iVariable+number2*jVariable];
	    }
	    h[iVariable] = value;
	  }
	  double hh[10000];
	  double yhy1 = innerProduct(h,number2,saveY2)*multiplier1+1.0;
	  yhy1 *= multiplier1;
	  for (iVariable=0;iVariable<number2;iVariable++) {
	    for (int jVariable=0;jVariable<number2;jVariable++) {
	      int put = iVariable+number2*jVariable;
	      double value = zzz[put];
	      value += yhy1*saveS2[iVariable]*saveS2[jVariable];
	      hh[put]=value;
	    }
	  }
	  for (iVariable=0;iVariable<number2;iVariable++) {
	    for (int jVariable=0;jVariable<number2;jVariable++) {
	      int put = iVariable+number2*jVariable;
	      double value = hh[put];
	      value -= multiplier1*(saveS2[iVariable]*h[jVariable]);
	      value -= multiplier1*(saveS2[jVariable]*h[iVariable]);
	      hh[put]=value;
	    }
	  }
#endif
	  // now update H
	  for (iVariable=0;iVariable<number2;iVariable++) {
	    for (int jVariable=0;jVariable<number2;jVariable++) {
	      int put = iVariable+number2*jVariable;
	      double value = zzz[put];
	      value += multiplier1*(ww[iVariable]*saveS2[jVariable]
				    +ww[jVariable]*saveS2[iVariable]);
	      value -= multiplier2*saveS2[iVariable]*saveS2[jVariable];
	      zzz[put]=value;
	    }
	  }
	  //memcpy(zzz,hh,size*sizeof(double));
	  // do search direction
	  memset(dArray,0,numberTotal*sizeof(double));
	  for (iVariable=0;iVariable<numberNonBasic;iVariable++) {
	    double value=0.0;
	    for (int jVariable=0;jVariable<number2;jVariable++) {
	      int k=which[jVariable];
	      value += work[k]*zzz[iVariable+number2*jVariable];
	    }
	    int i=which[iVariable];
	    dArray[i] = value;
	  }
	  // Now fill out dArray
	  {
	    int j;
	    // Now do basic ones 
	    int iRow;
	    CoinIndexedVector * spare1=spare;
	    CoinIndexedVector * spare2=rowArray;
#if CLP_DEBUG > 1
	    spare1->checkClear();
	    spare2->checkClear();
#endif
	    double * array2=spare1->denseVector();
	    int * index2 = spare1->getIndices();
	    int number2=0;
	    times(-1.0,dArray,array2);
	    dArray = dArray+numberColumns_;
	    for (iRow=0;iRow<numberRows_;iRow++) {
	      double value = array2[iRow] + dArray[iRow];
	      if (value) {
		array2[iRow]=value;
		index2[number2++]=iRow;
	      } else {
		array2[iRow]=0.0;
	      }
	    }
	    dArray -= numberColumns_;
	    spare1->setNumElements(number2);
	    // Ftran
	    factorization_->updateColumn(spare2,spare1);
	    number2=spare1->getNumElements();
	    for (j=0;j<number2;j++) {
	      int iSequence = index2[j];
	      double value = array2[iSequence];
	      array2[iSequence]=0.0;
	      if (value) {
		int iPivot=pivotVariable_[iSequence];
		double oldValue = dArray[iPivot];
		dArray[iPivot] = value+oldValue;
	      }
	    }
	    spare1->setNumElements(0);
	  }
	  //assert (innerProductIndexed(dArray,number2,work,which)>0);
	  //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
	  printf ("innerP1 %g innerP2 %g\n",innerProduct(dArray,numberTotal,work),
		  innerProductIndexed(dArray,numberNonBasic,work,which));
	  assert (innerProduct(dArray,numberTotal,work)>0);
#if 1
	  {
	    // check null vector
	    int iRow;
	    double qq[10];
	    memset(qq,0,numberRows_*sizeof(double));
	    multiplyAdd(dArray+numberColumns_,numberRows_,-1.0,qq,0.0);
	    matrix_->times(1.0,dArray,qq,rowScale_,columnScale_);
	    double largest=0.0;
	    int iLargest=-1;
	    for (iRow=0;iRow<numberRows_;iRow++) {
	      double value=fabs(qq[iRow]);
	      if (value>largest) {
		largest=value;
		iLargest=iRow;
	      }
	    }
	    printf("largest null error %g on row %d\n",largest,iLargest);
	    for (iSequence=0;iSequence<numberTotal;iSequence++) {
	      if (getStatus(iSequence)==basic)
		assert (fabs(dj_[iSequence])<1.0e-3);
	    }
	  }
#endif
	  memcpy(saveY,saveY2,numberNonBasic*sizeof(double));
	  memcpy(saveS,saveS2,numberNonBasic*sizeof(double));
	}
#endif
      }
#if MINTYPE==2
      for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
	int j=which[iSequence];
	saveY2[iSequence]=-work[j];
	saveS2[iSequence]=solution_[j];
      }
#endif
      if (djNorm<eps*djNorm0||(nPasses>100&&djNorm<CoinMin(1.0e-1*djNorm0, 1.0e-12))) {
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("dj norm reduced from %g to %g\n",djNorm0,djNorm);
#endif
	if (pivotMode<10||!numberNonBasic) {
	  finished=true;
	} else {
	  localPivotMode=pivotMode;
	  nTotalPasses+= nPasses;
	  nPasses=0;
	  continue;
	}
      }
      //if (nPasses==100||nPasses==50) 
      //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
      //	 normFlagged);
      if (nPasses>100&&djNorm<1.0e-2*normFlagged&&!startLocalMode) {
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n",djNorm0,djNorm,
		 normFlagged);
#endif
	unflag();
	localPivotMode=0;
	nTotalPasses+= nPasses;
	nPasses=0;
	continue;
      }
      if (djNorm>0.99*djNorm0&&nPasses>1500) {
	finished=true;
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("dj norm NOT reduced from %g to %g\n",djNorm0,djNorm);
#endif
	djNorm=1.2345e-20;
      }
      number=longArray->getNumElements();
      if (!numberNonBasic) {
	// looks optimal
	// check if any dj look plausible
	int nSuper=0;
	int nFlagged=0;
	int nNormal=0;
	for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
	  if (flagged(iSequence)) {
	    switch(getStatus(iSequence)) {
	      
	    case basic:
	    case ClpSimplex::isFixed:
	      break;
	    case atUpperBound:
	      if (dj_[iSequence]>dualTolerance_) 
		nFlagged++;
	      break;
	    case atLowerBound:
	      if (dj_[iSequence]<-dualTolerance_) 
		nFlagged++;
	      break;
	    case isFree:
	    case superBasic:
	      if (fabs(dj_[iSequence])>dualTolerance_) 
		nFlagged++;
	      break;
	    }
	    continue;
	  }
	  switch(getStatus(iSequence)) {
	    
	  case basic:
	  case ClpSimplex::isFixed:
	    break;
	  case atUpperBound:
	    if (dj_[iSequence]>dualTolerance_) 
	      nNormal++;
	    break;
	  case atLowerBound:
	    if (dj_[iSequence]<-dualTolerance_) 
	      nNormal++;
	    break;
	  case isFree:
	  case superBasic:
	    if (fabs(dj_[iSequence])>dualTolerance_) 
	      nSuper++;
	    break;
	  }
	}
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("%d super, %d normal, %d flagged\n",
		 nSuper,nNormal,nFlagged);
#endif
	
	int nFlagged2=1;
	if (lastSequenceIn<0&&!nNormal&&!nSuper) {
	  nFlagged2=unflag();
	  if (pivotMode>=10) {
	    pivotMode--;
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("pivot mode now %d\n",pivotMode);
#endif
	    if (pivotMode==9) 
	      pivotMode=0;	// switch off fast attempt
	  }
	} else {
	  lastSequenceIn=-1;
	}
	if (!nFlagged2||(normFlagged+normUnflagged<1.0e-8)) {
	  primalColumnPivot_->setLooksOptimal(true);
	  return 0;
	} else {
	  localPivotMode=-1;
	  nTotalPasses+= nPasses;
	  nPasses=0;
	  finished=false;
	  continue;
	}
      }
      bestSequence=-1;
      bestBasicSequence=-1;
      
      // temp
      nPasses++;
      if (nPasses>2000)
	finished=true;
      double theta = 1.0e30;
      basicTheta = 1.0e30;
      theta_=1.0e30;
      double basicTolerance = 1.0e-4*primalTolerance_;
      for (iSequence=0;iSequence<numberTotal;iSequence++) {
	//if (flagged(iSequence)
	//  continue;
	double alpha = dArray[iSequence];
	Status thisStatus = getStatus(iSequence);
	double oldValue = solution_[iSequence];
	if (thisStatus!=basic) {
	  if (fabs(alpha)>=acceptablePivot) {
	    if (alpha<0.0) {
	      // variable going towards lower bound
	      double bound = lower_[iSequence];
	      oldValue -= bound;
	      if (oldValue+theta*alpha<0.0) {
		bestSequence=iSequence;
		theta = CoinMax(0.0,oldValue/(-alpha));
	      }
	    } else {
	      // variable going towards upper bound
	      double bound = upper_[iSequence];
	      oldValue = bound-oldValue;
	      if (oldValue-theta*alpha<0.0) {
		bestSequence=iSequence;
		theta = CoinMax(0.0,oldValue/alpha);
	      }
	    }
	  }
	} else {
	  if (fabs(alpha)>=acceptableBasic) {
	    if (alpha<0.0) {
	      // variable going towards lower bound
	      double bound = lower_[iSequence];
	      oldValue -= bound;
	      if (oldValue+basicTheta*alpha<-basicTolerance) {
		bestBasicSequence=iSequence;
		basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/(-alpha));
	      }
	    } else {
	      // variable going towards upper bound
	      double bound = upper_[iSequence];
	      oldValue = bound-oldValue;
	      if (oldValue-basicTheta*alpha<-basicTolerance) {
		bestBasicSequence=iSequence;
		basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/alpha);
	      }
	    }
	  }
	}
      }
      theta_ = CoinMin(theta,basicTheta);
      // Now find minimum of function
      double objTheta2 =objective_->stepLength(this,solution_,dArray,CoinMin(theta,basicTheta),
					       currentObj,predictedObj,thetaObj);
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) 
	printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
#endif
#if MINTYPE==1
      if (conjugate) {
	double offset;
	const double * gradient = objective_->gradient(this,
						       dArray, offset,
						       true,0);
	double product=0.0;
	for (iSequence = 0;iSequence<numberColumns_;iSequence++) {
	  double alpha = dArray[iSequence];
	  double value = alpha*gradient[iSequence];
	  product += value;
	}
	//#define INCLUDESLACK
#ifdef INCLUDESLACK
	for (;iSequence<numberColumns_+numberRows_;iSequence++) {
	  double alpha = dArray[iSequence];
	  double value = alpha*cost_[iSequence];
	  product += value;
	}
#endif
	if (product>0.0)
	  objTheta = djNorm/product;
	else
	  objTheta=COIN_DBL_MAX;
#ifndef NDEBUG
	if (product<-1.0e-8&&handler_->logLevel()>1) 
          printf("bad product %g\n",product);
#endif
        product = CoinMax(product,0.0);
      } else {
	objTheta =objTheta2;
      }
#else
      objTheta =objTheta2;
#endif
      // if very small difference then take pivot (depends on djNorm?)
      // distinguish between basic and non-basic
      bool chooseObjTheta=objTheta<theta_;
      if (chooseObjTheta) {
	if (thetaObj<currentObj-1.0e-12&&objTheta+1.0e-10>theta_)
	  chooseObjTheta=false;
	//if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
	//chooseObjTheta=false;
      }
      //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
      if (chooseObjTheta) {
	theta_ = objTheta;
      } else {
	objTheta = CoinMax(objTheta,1.00000001*theta_+1.0e-12);
	//if (theta+1.0e-13>basicTheta) {
	//theta = CoinMax(theta,1.00000001*basicTheta);
	//theta_ = basicTheta;
	//}
      }
      // always out if one variable in and zero move
      if (theta_==basicTheta||(sequenceIn_>=0&&!theta_))
	finished=true; // come out
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) {
	printf("%d pass,",nPasses);
	if (sequenceIn_>=0)
	  printf (" Sin %d,",sequenceIn_);
	if (basicTheta==theta_)
	  printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
	else
	  printf(" basicTheta %g",basicTheta);
	if (theta==theta_)
	  printf(" X(%d) non-basicTheta %g",bestSequence,theta);
	else
	  printf(" non-basicTheta %g",theta);
	printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
	printf(" djNorm %g\n",djNorm);
      }
#endif
      if (handler_->logLevel()>3&&objTheta!=theta_) {
	printf("%d pass obj %g,",nPasses,currentObj);
	if (sequenceIn_>=0)
	  printf (" Sin %d,",sequenceIn_);
	if (basicTheta==theta_)
	  printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
	else
	  printf(" basicTheta %g",basicTheta);
	if (theta==theta_)
	  printf(" X(%d) non-basicTheta %g",bestSequence,theta);
	else
	  printf(" non-basicTheta %g",theta);
	printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
	printf(" djNorm %g\n",djNorm);
      }
      if (objTheta!=theta_) {
	//printf("hit boundary after %d passes\n",nPasses);
	nTotalPasses+= nPasses;
	nPasses=0; // start again
      }
      if (localPivotMode<10||lastSequenceIn==bestSequence) {
	if (theta_==theta&&theta_<basicTheta&&theta_<1.0e-5)
	  setFlagged(bestSequence); // out of active set
      }
      // Update solution
      for (iSequence=0;iSequence<numberTotal;iSequence++) {
	//for (iIndex=0;iIndex<number;iIndex++) {
	
	//int iSequence = which[iIndex];
	double alpha = dArray[iSequence];
	if (alpha) {
	  double value = solution_[iSequence]+theta_*alpha;
	  solution_[iSequence]=value;
	  switch(getStatus(iSequence)) {
	    
	  case basic:
	  case isFixed:
	  case isFree:
	  case atUpperBound:
	  case atLowerBound:
	    nonLinearCost_->setOne(iSequence,value);
	    break;
	  case superBasic:
	    // To get correct action
	    setStatus(iSequence,isFixed);
	    nonLinearCost_->setOne(iSequence,value);
	    //assert (getStatus(iSequence)!=isFixed);
	    break;
	  }
	}
      }
      if (objTheta<SMALLTHETA2&&objTheta==theta_) {
	if (sequenceIn_>=0&&basicTheta<1.0e-9) {
	  // try just one
	  localPivotMode=0;
	  sequenceIn_=-1;
	  nTotalPasses+= nPasses;
	  nPasses=0;
	  //finished=true;
	  //objTheta=0.0;
	  //theta_=0.0;
	}
	//if (objTheta<1.0e-12) 
	//finished=true;
	//printf("zero move\n");
      }
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) {
	objective_->stepLength(this,solution_,work,0.0,
			       currentObj,predictedObj,thetaObj);
	printf("current obj %g after update - simple was %g\n",currentObj,simpleObjective);
      }
#endif
      if (sequenceIn_>=0&&!finished&&objTheta>1.0e-4) {
	// we made some progress - back to normal
	if (localPivotMode<10) {
	  localPivotMode=0;
	  sequenceIn_=-1;
	  nTotalPasses+= nPasses;
	  nPasses=0;
	}
#ifdef CLP_DEBUG
	if (handler_->logLevel()&32) 
	  printf("some progress?\n");
#endif
      }
#if CLP_DEBUG > 1
      longArray->checkClean();
#endif
    }
#ifdef CLP_DEBUG
    if (handler_->logLevel()&32) 
      printf("out of loop after %d passes\n",nPasses);
#endif
    bool ordinaryDj=false;
    //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta) 
    //printf("could try ordinary iteration %g\n",theta_);
    if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-15) {
      //printf("trying ordinary iteration\n");
      ordinaryDj=true;
    }
    if (!basicTheta&&!ordinaryDj) {
      //returnCode=2;
      //objTheta=-1.0; // so we fall through
    }
    assert (theta_<1.0e30); // for now
    // See if we need to pivot
    if (theta_==basicTheta||ordinaryDj) {
      if (!ordinaryDj) {
	// find an incoming column which will force pivot
	int iRow;
	pivotRow_=-1;
	for (iRow=0;iRow<numberRows_;iRow++) {
	  if (pivotVariable_[iRow]==bestBasicSequence) {
	    pivotRow_=iRow;
	    break;
	  }
	}
	assert (pivotRow_>=0);
	// Get good size for pivot
	double acceptablePivot=1.0e-7;
	if (factorization_->pivots()>10)
	  acceptablePivot=1.0e-5; // if we have iterated be more strict
	else if (factorization_->pivots()>5)
	  acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
	// should be dArray but seems better this way!
	double direction = work[bestBasicSequence]>0.0 ? -1.0 : 1.0;
	// create as packed
	rowArray->createPacked(1,&pivotRow_,&direction);
	factorization_->updateColumnTranspose(spare,rowArray);
	// put row of tableau in rowArray and columnArray
	matrix_->transposeTimes(this,-1.0,
				rowArray,spare,columnArray);
	// choose one futhest away from bound which has reasonable pivot
	// If increasing we want negative alpha
	
	double * work2;
	int iSection;
	
	sequenceIn_=-1;
	double bestValue=-1.0;
	double bestDirection=0.0;
	// First pass we take correct direction and large pivots
	// then correct direction
	// then any
	double check[]={1.0e-8,-1.0e-12,-1.0e30};
	double mult[]={100.0,1.0,1.0};
	for (int iPass=0;iPass<3;iPass++) {
	  //if (!bestValue&&iPass==2)
	  //bestValue=-1.0;
	  double acceptable=acceptablePivot*mult[iPass];
	  double checkValue=check[iPass];
	  for (iSection=0;iSection<2;iSection++) {
	    
	    int addSequence;
	    
	    if (!iSection) {
	      work2 = rowArray->denseVector();
	      number = rowArray->getNumElements();
	      which = rowArray->getIndices();
	      addSequence = numberColumns_;
	    } else {
	      work2 = columnArray->denseVector();
	      number = columnArray->getNumElements();
	      which = columnArray->getIndices();
	      addSequence = 0;
	    }
	    int i;
	    
	    for (i=0;i<number;i++) {
	      int iSequence = which[i]+addSequence;
	      if (flagged(iSequence))
		continue;
	      //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
	      //		  upper_[iSequence]-solution_[iSequence]);
	      double alpha=work2[i];
	      // should be dArray but seems better this way!
	      double change=work[iSequence];
	      Status thisStatus = getStatus(iSequence);
	      double direction=0;;
	      switch(thisStatus) {
		
	      case basic:
	      case ClpSimplex::isFixed:
		break;
	      case isFree:
	      case superBasic:
		if (alpha<-acceptable&&change>checkValue)
		  direction=1.0;
		else if (alpha>acceptable&&change<-checkValue)
		  direction=-1.0;
		break;
	      case atUpperBound:
		if (alpha>acceptable&&change<-checkValue)
		  direction=-1.0;
		else if (iPass==2&&alpha<-acceptable&&change<-checkValue)
		  direction=1.0;
		break;
	      case atLowerBound:
		if (alpha<-acceptable&&change>checkValue)
		  direction=1.0;
		else if (iPass==2&&alpha>acceptable&&change>checkValue)
		  direction=-1.0;
		break;
	      }
	      if (direction) {
		if (sequenceIn_!=lastSequenceIn||localPivotMode<10) { 
		  if (CoinMin(solution_[iSequence]-lower_[iSequence],
			  upper_[iSequence]-solution_[iSequence])>bestValue) {
		    bestValue=CoinMin(solution_[iSequence]-lower_[iSequence],
				  upper_[iSequence]-solution_[iSequence]);
		    sequenceIn_=iSequence;
		    bestDirection=direction;
		  }
		} else {
		  // choose
		  bestValue=COIN_DBL_MAX;
		  sequenceIn_=iSequence;
		  bestDirection=direction;
		}
	      }
	    }
	  }
	  if (sequenceIn_>=0&&bestValue>0.0)
	    break;
	}
	if (sequenceIn_>=0) {
	  valueIn_=solution_[sequenceIn_];
	  dualIn_=dj_[sequenceIn_];
	  if (bestDirection<0.0) {
	    // we want positive dj
	    if (dualIn_<=0.0) {
	      //printf("bad dj - xx %g\n",dualIn_);
	      dualIn_=1.0;
	    }
	  } else {
	    // we want negative dj
	    if (dualIn_>=0.0) {
	      //printf("bad dj - xx %g\n",dualIn_);
	      dualIn_=-1.0;
	    }
	  }
	  lowerIn_=lower_[sequenceIn_];
	  upperIn_=upper_[sequenceIn_];
	  if (dualIn_>0.0)
	    directionIn_ = -1;
	  else 
	    directionIn_ = 1;
	} else {
	  //ordinaryDj=true;
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32) {
	    printf("no easy pivot - norm %g mode %d\n",djNorm,localPivotMode);
	    if (rowArray->getNumElements()+columnArray->getNumElements()<12) {
	      for (iSection=0;iSection<2;iSection++) {
		
		int addSequence;
		
		if (!iSection) {
		  work2 = rowArray->denseVector();
		  number = rowArray->getNumElements();
		  which = rowArray->getIndices();
		  addSequence = numberColumns_;
		} else {
		  work2 = columnArray->denseVector();
		  number = columnArray->getNumElements();
		  which = columnArray->getIndices();
		  addSequence = 0;
		}
		int i;
		char section[]={'R','C'};
		for (i=0;i<number;i++) {
		  int iSequence = which[i]+addSequence;
		  if (flagged(iSequence)) {
		    printf("%c%d flagged\n",section[iSection],which[i]);
		    continue;
		  } else {
		    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
			   section[iSection],which[i],status_[iSequence],
			   lower_[iSequence],solution_[iSequence],upper_[iSequence],
			   work2[i],work[iSequence]);
		  }
		}
	      }
	    }
	  }
#endif
	  assert (sequenceIn_<0);
	  justOne=true;
	  allFinished=false; // Round again
	  finished=false;
	  nTotalPasses+= nPasses;
	  nPasses=0;
	  if (djNorm<0.9*djNorm0&&djNorm<1.0e-3&&!localPivotMode) {
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("no pivot - mode %d norms %g %g - unflagging\n",
		     localPivotMode,djNorm0,djNorm);
#endif
	    unflag(); //unflagging
	    returnCode=1;
	  } else {
	    returnCode=2; // do single incoming
	    //returnCode=1;
	  }
	}
      }
      rowArray->clear();
      columnArray->clear();
      longArray->clear();
      if (ordinaryDj) {
	valueIn_=solution_[sequenceIn_];
	dualIn_=dj_[sequenceIn_];
	lowerIn_=lower_[sequenceIn_];
	upperIn_=upper_[sequenceIn_];
	if (dualIn_>0.0)
	  directionIn_ = -1;
	else 
	  directionIn_ = 1;
      }
      if (returnCode==1)
	returnCode=0;
    } else {
      // round again
      longArray->clear();
      if (djNorm<1.0e-3&&!localPivotMode) {
	if (djNorm==1.2345e-20&&djNorm0>1.0e-4) {
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32) 
	    printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n",djNorm0,nPasses,
		   localPivotMode,returnCode);
#endif
	  //if (!localPivotMode)
	  //returnCode=2; // force singleton
	} else {
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32) 
	    printf("unflagging as djNorm %g %g, %d passes\n",djNorm,djNorm0,nPasses);
#endif
	  if (pivotMode>=10) {
	    pivotMode--;
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("pivot mode now %d\n",pivotMode);
#endif
	    if (pivotMode==9) 
	      pivotMode=0;	// switch off fast attempt
	  }
	  bool unflagged = unflag()!=0;
	  if (!unflagged&&djNorm<1.0e-10) {
	    // ? declare victory
	    sequenceIn_=-1;
	    returnCode=0;
	  }
	}
      }
    }
  }
  if (djNorm0<1.0e-12*normFlagged) {
#ifdef CLP_DEBUG
    if (handler_->logLevel()&32) 
      printf("unflagging as djNorm %g %g and flagged norm %g\n",djNorm,djNorm0,normFlagged);
#endif
    unflag();
  }
  if (saveObj-currentObj<1.0e-5&&nTotalPasses>2000) {
    normUnflagged=0.0;
    double dualTolerance3 = CoinMin(1.0e-2,1.0e3*dualTolerance_);
    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
      switch(getStatus(iSequence)) {
	
      case basic:
      case ClpSimplex::isFixed:
	break;
      case atUpperBound:
	if (dj_[iSequence]>dualTolerance3) 
	  normFlagged += dj_[iSequence]*dj_[iSequence];
	break;
      case atLowerBound:
	if (dj_[iSequence]<-dualTolerance3) 
	  normFlagged += dj_[iSequence]*dj_[iSequence];
	break;
      case isFree:
      case superBasic:
	if (fabs(dj_[iSequence])>dualTolerance3) 
	  normFlagged += dj_[iSequence]*dj_[iSequence];
	break;
      }
    }
    if (handler_->logLevel()>2)
      printf("possible optimal  %d %d %g %g\n",
	     nBigPasses,nTotalPasses,saveObj-currentObj,normFlagged);
    if (normFlagged<1.0e-5) {
      unflag();
      primalColumnPivot_->setLooksOptimal(true);
      returnCode=0;
    }
  }
  return returnCode;
}
/* Do last half of an iteration.
   Return codes
   Reasons to come out normal mode
   -1 normal
   -2 factorize now - good iteration
   -3 slight inaccuracy - refactorize - iteration done
   -4 inaccuracy - refactorize - no iteration
   -5 something flagged - go round again
   +2 looks unbounded
   +3 max iterations (iteration done)
   
*/
int 
ClpSimplexNonlinear::pivotNonlinearResult()
{

  int returnCode=-1;

  rowArray_[1]->clear();
  
  // we found a pivot column
  // update the incoming column
  unpackPacked(rowArray_[1]);
  factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
  theta_=0.0;
  double * work=rowArray_[1]->denseVector();
  int number=rowArray_[1]->getNumElements();
  int * which=rowArray_[1]->getIndices();
  bool keepValue=false;
  double saveValue=0.0;
  if (pivotRow_>=0) {
    sequenceOut_=pivotVariable_[pivotRow_];;
    valueOut_ = solution(sequenceOut_);
    keepValue=true;
    saveValue=valueOut_;
    lowerOut_=lower_[sequenceOut_];
    upperOut_=upper_[sequenceOut_];
    int iIndex;
    for (iIndex=0;iIndex<number;iIndex++) {
      
      int iRow = which[iIndex];
      if (iRow==pivotRow_) {
	alpha_ = work[iIndex];
	break;
      }
    }
  } else {
    int iIndex;
    double smallest=COIN_DBL_MAX;
    for (iIndex=0;iIndex<number;iIndex++) {
      int iRow = which[iIndex];
      double alpha = work[iIndex];
      if (fabs(alpha)>1.0e-6) {
	int iPivot = pivotVariable_[iRow];
	double distance = CoinMin(upper_[iPivot]-solution_[iPivot],
			      solution_[iPivot]-lower_[iPivot]);
	if (distance<smallest) {
	  pivotRow_=iRow;
	  alpha_=alpha;
	  smallest=distance;
	}
      }
    }
    if (smallest>primalTolerance_) {
      smallest=COIN_DBL_MAX;
      for (iIndex=0;iIndex<number;iIndex++) {
	int iRow = which[iIndex];
	double alpha = work[iIndex];
	if (fabs(alpha)>1.0e-6) {
	  double distance = CoinDrand48();
	  if (distance<smallest) {
	    pivotRow_=iRow;
	    alpha_=alpha;
	    smallest=distance;
	  }
	}
      }
    }
    assert (pivotRow_>=0);
    sequenceOut_=pivotVariable_[pivotRow_];;
    valueOut_ = solution(sequenceOut_);
    lowerOut_=lower_[sequenceOut_];
    upperOut_=upper_[sequenceOut_];
  }
  double newValue = valueOut_ - theta_*alpha_;
  bool isSuperBasic=false;
  if (valueOut_>=upperOut_-primalTolerance_) {
    directionOut_=-1;      // to upper bound
    upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
    upperOut_ = newValue;
  } else if (valueOut_<=lowerOut_+primalTolerance_) {
    directionOut_=1;      // to lower bound
    lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
  } else {
    lowerOut_=valueOut_;
    upperOut_=valueOut_;
    isSuperBasic=true;
    //printf("XX superbasic out\n");
  }
  dualOut_ = dj_[sequenceOut_];
  double checkValue=1.0e-2;
  if (largestDualError_>1.0e-5)
    checkValue=1.0e-1;
  // if stable replace in basis
  
  int updateStatus = factorization_->replaceColumn(this,
						   rowArray_[2],
						   rowArray_[1],
						   pivotRow_,
						   alpha_);
  
  // if no pivots, bad update but reasonable alpha - take and invert
  if (updateStatus==2&&
      lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
    updateStatus=4;
  if (updateStatus==1||updateStatus==4) {
    // slight error
    if (factorization_->pivots()>5||updateStatus==4) {
      returnCode=-3;
    }
  } else if (updateStatus==2) {
    // major error
    // better to have small tolerance even if slower
    factorization_->zeroTolerance(1.0e-15);
    int maxFactor = factorization_->maximumPivots();
    if (maxFactor>10) {
      if (forceFactorization_<0)
	forceFactorization_= maxFactor;
      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
    } 
    // later we may need to unwind more e.g. fake bounds
    if(lastGoodIteration_ != numberIterations_) {
      clearAll();
      pivotRow_=-1;
      returnCode=-4;
    } else {
      // need to reject something
      char x = isColumn(sequenceIn_) ? 'C' :'R';
      handler_->message(CLP_SIMPLEX_FLAG,messages_)
	<<x<<sequenceWithin(sequenceIn_)
	<<CoinMessageEol;
      setFlagged(sequenceIn_);
      progress_->clearBadTimes();
      lastBadIteration_ = numberIterations_; // say be more cautious
      clearAll();
      pivotRow_=-1;
      sequenceOut_=-1;
      returnCode = -5;
      
    }
    return returnCode;
  } else if (updateStatus==3) {
    // out of memory
    // increase space if not many iterations
    if (factorization_->pivots()<
	0.5*factorization_->maximumPivots()&&
	factorization_->pivots()<200)
      factorization_->areaFactor(
				 factorization_->areaFactor() * 1.1);
    returnCode =-2; // factorize now
  } else if (updateStatus==5) {
    problemStatus_=-2; // factorize now
  }
  
  // update primal solution
  
  double objectiveChange=0.0;
  // after this rowArray_[1] is not empty - used to update djs
  // If pivot row >= numberRows then may be gub
  updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,1);
  
  double oldValue = valueIn_;
  if (directionIn_==-1) {
    // as if from upper bound
    if (sequenceIn_!=sequenceOut_) {
      // variable becoming basic
      valueIn_ -= fabs(theta_);
    } else {
      valueIn_=lowerIn_;
    }
  } else {
    // as if from lower bound
    if (sequenceIn_!=sequenceOut_) {
      // variable becoming basic
      valueIn_ += fabs(theta_);
    } else {
      valueIn_=upperIn_;
    }
  }
  objectiveChange += dualIn_*(valueIn_-oldValue);
  // outgoing
  if (sequenceIn_!=sequenceOut_) {
    if (directionOut_>0) {
      valueOut_ = lowerOut_;
    } else {
      valueOut_ = upperOut_;
    }
    if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
      valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
    else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
      valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
    // may not be exactly at bound and bounds may have changed
    // Make sure outgoing looks feasible
    if (!isSuperBasic)
      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
    solution_[sequenceOut_]=valueOut_;
  }
  // change cost and bounds on incoming if primal
  nonLinearCost_->setOne(sequenceIn_,valueIn_); 
  int whatNext=housekeeping(objectiveChange);
  if (keepValue)
    solution_[sequenceOut_]=saveValue;
  if (isSuperBasic)
    setStatus(sequenceOut_,superBasic);
  //#define CLP_DEBUG
#if CLP_DEBUG > 1
  {
    int ninf= matrix_->checkFeasible(this);
    if (ninf)
      printf("infeas %d\n",ninf);
  }
#endif
  if (whatNext==1) {
    returnCode =-2; // refactorize
  } else if (whatNext==2) {
    // maximum iterations or equivalent
    returnCode=3;
  } else if(numberIterations_ == lastGoodIteration_
	    + 2 * factorization_->maximumPivots()) {
    // done a lot of flips - be safe
    returnCode =-2; // refactorize
  }
  // Check event
  {
    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
    if (status>=0) {
      problemStatus_=5;
      secondaryStatus_=ClpEventHandler::endOfIteration;
      returnCode=4;
    }
  }
  return returnCode;
}
// A sequential LP method
int 
ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance)
{
  // Are we minimizing or maximizing
  double whichWay=optimizationDirection();
  if (whichWay<0.0)
    whichWay=-1.0;
  else if (whichWay>0.0)
    whichWay=1.0;


  int numberColumns = this->numberColumns();
  int numberRows = this->numberRows();
  double * columnLower = this->columnLower();
  double * columnUpper = this->columnUpper();
  double * solution = this->primalColumnSolution();
  
  if (objective_->type()<2) {
    // no nonlinear part
    return ClpSimplex::primal(0);
  }
  // Get list of non linear columns
  char * markNonlinear = new char[numberColumns];
  memset(markNonlinear,0,numberColumns);
  int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
  
  if (!numberNonLinearColumns) {
    delete [] markNonlinear;
    // no nonlinear part
    return ClpSimplex::primal(0);
  }
  int iColumn;
  int * listNonLinearColumn = new int[numberNonLinearColumns];
  numberNonLinearColumns=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if(markNonlinear[iColumn])
      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
  }
  // Replace objective
  ClpObjective * trueObjective = objective_;
  objective_=new ClpLinearObjective(NULL,numberColumns);
  double * objective = this->objective();
  
  // get feasible
  if (this->status()<0||numberPrimalInfeasibilities())
    ClpSimplex::primal(1);
  // still infeasible
  if (numberPrimalInfeasibilities()) {
    delete [] listNonLinearColumn;
    return 0;
  }
  algorithm_=1;
  int jNon;
  int * last[3];
  
  double * trust = new double[numberNonLinearColumns];
  double * trueLower = new double[numberNonLinearColumns];
  double * trueUpper = new double[numberNonLinearColumns];
  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    iColumn=listNonLinearColumn[jNon];
    trust[jNon]=0.5;
    trueLower[jNon]=columnLower[iColumn];
    trueUpper[jNon]=columnUpper[iColumn];
    if (solution[iColumn]<trueLower[jNon])
      solution[iColumn]=trueLower[jNon];
    else if (solution[iColumn]>trueUpper[jNon])
      solution[iColumn]=trueUpper[jNon];
  }
  int saveLogLevel=logLevel();
  int iPass;
  double lastObjective=1.0e31;
  double * saveSolution = new double [numberColumns];
  double * saveRowSolution = new double [numberRows];
  memset(saveRowSolution,0,numberRows*sizeof(double));
  double * savePi = new double [numberRows];
  double * safeSolution = new double [numberColumns];
  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
#define MULTIPLE 0
#if MULTIPLE>2
  // Duplication but doesn't really matter
  double * saveSolutionM[MULTIPLE};
  for (jNon=0;jNon<MULTIPLE;jNon++) {
    saveSolutionM[jNon]=new double[numberColumns];
    memcpy(saveSolutionM,solution,numberColumns*sizeof(double));
  }
#endif
  double targetDrop=1.0e31;
  double objectiveOffset;
  getDblParam(ClpObjOffset,objectiveOffset);
  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
  for (iPass=0;iPass<3;iPass++) {
    last[iPass]=new int[numberNonLinearColumns];
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
      last[iPass][jNon]=0;
  }
  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
  int goodMove=-2;
  char * statusCheck = new char[numberColumns];
  double * changeRegion = new double [numberColumns];
  double offset=0.0;
  double objValue=0.0;
  int exitPass = 2*numberPasses+10;
  for (iPass=0;iPass<numberPasses;iPass++) {
    exitPass--;
    // redo objective
    offset=0.0;
    objValue=-objectiveOffset;
    // make sure x updated
    trueObjective->newXValues();
    double theta=-1.0;
    double maxTheta=COIN_DBL_MAX;
    //maxTheta=1.0;
    if (iPass) {
      int jNon=0;
      for (iColumn=0;iColumn<numberColumns;iColumn++) { 
	changeRegion[iColumn]=solution[iColumn]-saveSolution[iColumn];
	double alpha = changeRegion[iColumn];
	double oldValue = saveSolution[iColumn];
	if (markNonlinear[iColumn]==0) {
	  // linear
	  if (alpha<-1.0e-15) {
	    // variable going towards lower bound
	    double bound = columnLower[iColumn];
	    oldValue -= bound;
	    if (oldValue+maxTheta*alpha<0.0) {
	      maxTheta = CoinMax(0.0,oldValue/(-alpha));
	    }
	  } else if (alpha>1.0e-15) {
	    // variable going towards upper bound
	    double bound = columnUpper[iColumn];
	    oldValue = bound-oldValue;
	    if (oldValue-maxTheta*alpha<0.0) {
	      maxTheta = CoinMax(0.0,oldValue/alpha);
	    }
	  }
	} else {
	  // nonlinear
	  if (alpha<-1.0e-15) {
	    // variable going towards lower bound
	    double bound = trueLower[jNon];
	    oldValue -= bound;
	    if (oldValue+maxTheta*alpha<0.0) {
	      maxTheta = CoinMax(0.0,oldValue/(-alpha));
	    }
	  } else if (alpha>1.0e-15) {
	    // variable going towards upper bound
	    double bound = trueUpper[jNon];
	    oldValue = bound-oldValue;
	    if (oldValue-maxTheta*alpha<0.0) {
	      maxTheta = CoinMax(0.0,oldValue/alpha);
	    }
	  }
	  jNon++;
	}
      }
      // make sure both accurate
      memset(rowActivity_,0,numberRows_*sizeof(double));
      times(1.0,solution,rowActivity_);
      memset(saveRowSolution,0,numberRows_*sizeof(double));
      times(1.0,saveSolution,saveRowSolution);
      for (int iRow=0;iRow<numberRows;iRow++) { 
	double alpha =rowActivity_[iRow]-saveRowSolution[iRow];
	double oldValue = saveRowSolution[iRow];
	if (alpha<-1.0e-15) {
	  // variable going towards lower bound
	  double bound = rowLower_[iRow];
	  oldValue -= bound;
	  if (oldValue+maxTheta*alpha<0.0) {
	    maxTheta = CoinMax(0.0,oldValue/(-alpha));
	  }
	} else if (alpha>1.0e-15) {
	  // variable going towards upper bound
	  double bound = rowUpper_[iRow];
	  oldValue = bound-oldValue;
	  if (oldValue-maxTheta*alpha<0.0) {
	    maxTheta = CoinMax(0.0,oldValue/alpha);
	  }
	}
      }
    } else {
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	changeRegion[iColumn]=0.0;
	saveSolution[iColumn]=solution[iColumn];
      }
      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
    }
    // get current value anyway
    double predictedObj,thetaObj;
    double maxTheta2 = 2.0; // to work out a b c
    double theta2 = trueObjective->stepLength(this,saveSolution,changeRegion,maxTheta2,
					     objValue,predictedObj,thetaObj);
    int lastMoveStatus = goodMove;
    if (goodMove>=0) {
      theta = CoinMin(theta2,maxTheta);
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) 
	printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
	       theta,objValue,thetaObj,predictedObj);
#endif
      if (theta>0.0&&theta<=1.0) {
	// update solution
	double lambda = 1.0-theta;
	for (iColumn=0;iColumn<numberColumns;iColumn++) 
	  solution[iColumn] = lambda * saveSolution[iColumn] 
	    + theta * solution[iColumn];
	memset(rowActivity_,0,numberRows_*sizeof(double));
	times(1.0,solution,rowActivity_);
	if (lambda>0.999) {
	  memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
	  memcpy(status_,saveStatus,numberRows+numberColumns);
	}
	// Do local minimization
#define LOCAL
#ifdef LOCAL
	bool absolutelyOptimal=false;
	int saveScaling = scalingFlag();
	scaling(0);
	int savePerturbation=perturbation_;
	perturbation_=100;
	if (saveLogLevel==1)
	  setLogLevel(0);
	int status=startup(1);
	if (!status) {
	  int numberTotal = numberRows_+numberColumns_;
	  // resize arrays
	  for (int i=0;i<4;i++) {
	    rowArray_[i]->reserve(CoinMax(numberRows_+numberColumns_,rowArray_[i]->capacity()));
	  }
	  CoinIndexedVector * longArray = rowArray_[3];
	  CoinIndexedVector * rowArray = rowArray_[0];
	  //CoinIndexedVector * columnArray = columnArray_[0];
	  CoinIndexedVector * spare = rowArray_[1];
	  double * work=longArray->denseVector();
	  int * which=longArray->getIndices();
	  int nPass=100;
	  //bool conjugate=false;
	  // Put back true bounds
	  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	    int iColumn=listNonLinearColumn[jNon];
	    double value;
	    value=trueLower[jNon];
	    trueLower[jNon]=lower_[iColumn];
	    lower_[iColumn]=value;
	    value=trueUpper[jNon];
	    trueUpper[jNon]=upper_[iColumn];
	    upper_[iColumn]=value;
	    switch(getStatus(iColumn)) {
	      
	    case basic:
	    case isFree:
	    case superBasic:
	      break;
	    case isFixed:
	    case atUpperBound:
	    case atLowerBound:
	      if (solution_[iColumn]>lower_[iColumn]+primalTolerance_) {
		if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
		  setStatus(iColumn,superBasic);
		} else {
		  setStatus(iColumn,atUpperBound);
		}
	      } else {
		if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
		  setStatus(iColumn,atLowerBound);
		} else {
		  setStatus(iColumn,isFixed);
		}
	      }
	      break;
	    }
	  }
	  for (int jPass=0;jPass<nPass;jPass++) {
	    trueObjective->reducedGradient(this,dj_,true);
	    // get direction vector in longArray
	    longArray->clear();
	    double normFlagged=0.0;
	    double normUnflagged=0.0;
	    int numberNonBasic=0;
	    directionVector(longArray,spare,rowArray,0,
			    normFlagged,normUnflagged,numberNonBasic);
	    if (normFlagged+normUnflagged<1.0e-8) {
	      absolutelyOptimal=true;
	      break;  //looks optimal
	    }
	    double djNorm=0.0;
	    int iIndex;
	    for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
	      int iSequence = which[iIndex];
	      double alpha = work[iSequence];
	      djNorm += alpha*alpha;
	    }
	    //if (!jPass)
	    //printf("dj norm %g - %d \n",djNorm,numberNonBasic);
	    //int number=longArray->getNumElements();
	    if (!numberNonBasic) {
	      // looks optimal
	      absolutelyOptimal=true;
	      break;
	    }
	    int bestSequence=-1;
	    double theta = 1.0e30;
	    int iSequence;
	    for (iSequence=0;iSequence<numberTotal;iSequence++) {
	      double alpha = work[iSequence];
	      double oldValue = solution_[iSequence];
	      if (alpha<-1.0e-15) {
		// variable going towards lower bound
		double bound = lower_[iSequence];
		oldValue -= bound;
		if (oldValue+theta*alpha<0.0) {
		  bestSequence=iSequence;
		  theta = CoinMax(0.0,oldValue/(-alpha));
		}
	      } else if (alpha>1.0e-15) {
		// variable going towards upper bound
		double bound = upper_[iSequence];
		oldValue = bound-oldValue;
		if (oldValue-theta*alpha<0.0) {
		  bestSequence=iSequence;
		  theta = CoinMax(0.0,oldValue/alpha);
		}
	      }
	    }
	    // Now find minimum of function
	    double currentObj;
	    double predictedObj;
	    double thetaObj;
	    // need to use true objective
	    double * saveCost = cost_;
	    cost_=NULL;
	    double objTheta =trueObjective->stepLength(this,solution_,work,theta,
						       currentObj,predictedObj,thetaObj);
	    cost_=saveCost;
#ifdef CLP_DEBUG
	    if (handler_->logLevel()&32) 
	      printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
#endif
	    //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
	    //printf("objTheta %g theta %g\n",objTheta,theta);
	    if (theta>objTheta) { 
	      theta=objTheta;
	      thetaObj=predictedObj;
	    }
	    // update one used outside
	    objValue=currentObj;
	    if (theta>1.0e-9&&
		(currentObj-thetaObj<-CoinMax(1.0e-8,1.0e-15*fabs(currentObj))||jPass<5)) {
	      // Update solution
	      for (iSequence=0;iSequence<numberTotal;iSequence++) {
		double alpha = work[iSequence];
		if (alpha) {
		  double value = solution_[iSequence]+theta*alpha;
		  solution_[iSequence]=value;
		  switch(getStatus(iSequence)) {
		    
		  case basic:
		  case isFixed:
		  case isFree:
		    break;
		  case atUpperBound:
		  case atLowerBound:
		  case superBasic:
		    if (value>lower_[iSequence]+primalTolerance_) {
		      if(value<upper_[iSequence]-primalTolerance_) {
			setStatus(iSequence,superBasic);
		      } else {
			setStatus(iSequence,atUpperBound);
		      }
		    } else {
		      if(value<upper_[iSequence]-primalTolerance_) {
			setStatus(iSequence,atLowerBound);
		      } else {
			setStatus(iSequence,isFixed);
		      }
		    }
		    break;
		  }
		}
	      }
	    } else {
	      break;
	    }
	  }
	  // Put back fake bounds
	  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	    int iColumn=listNonLinearColumn[jNon];
	    double value;
	    value=trueLower[jNon];
	    trueLower[jNon]=lower_[iColumn];
	    lower_[iColumn]=value;
	    value=trueUpper[jNon];
	    trueUpper[jNon]=upper_[iColumn];
	    upper_[iColumn]=value;
	  }
	}
	problemStatus_=0;
	finish();
	scaling(saveScaling);
	perturbation_=savePerturbation;
	setLogLevel(saveLogLevel);
#endif
	// redo rowActivity
	memset(rowActivity_,0,numberRows_*sizeof(double));
	times(1.0,solution,rowActivity_);
	if (theta>0.99999&&theta2<1.9&&!absolutelyOptimal) {
	  // If big changes then tighten
	  /*  thetaObj is objvalue + a*2*2 +b*2
	      predictedObj is objvalue + a*theta2*theta2 +b*theta2
	  */
	  double rhs1 = thetaObj-objValue;
	  double rhs2 = predictedObj-objValue;
	  double subtractB = theta2*0.5;
	  double a = (rhs2-subtractB*rhs1)/(theta2*theta2-4.0*subtractB);
	  double b = 0.5*(rhs1 - 4.0*a);
	  if (fabs(a+b)>1.0e-2) {
	    // tighten all
	    goodMove=-1;
	  }
	}
      }
    }
    memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
	   numberColumns*sizeof(double));
    //printf("offset comp %g orig %g\n",offset,objectiveOffset);
    // could do faster
    trueObjective->stepLength(this,solution,changeRegion,0.0,
					     objValue,predictedObj,thetaObj);
    printf("offset comp %g orig %g - obj from stepLength %g\n",offset,objectiveOffset,objValue);
    setDblParam(ClpObjOffset,objectiveOffset+offset);
    int * temp=last[2];
    last[2]=last[1];
    last[1]=last[0];
    last[0]=temp;
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
      iColumn=listNonLinearColumn[jNon];
      double change = solution[iColumn]-saveSolution[iColumn];
      if (change<-1.0e-5) {
	if (fabs(change+trust[jNon])<1.0e-5) 
	  temp[jNon]=-1;
	else
	  temp[jNon]=-2;
      } else if(change>1.0e-5) {
	if (fabs(change-trust[jNon])<1.0e-5) 
	  temp[jNon]=1;
	else
	  temp[jNon]=2;
      } else {
	temp[jNon]=0;
      }
    } 
    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
    double maxDelta=0.0;
    if (goodMove>=0) {
      if (objValue-lastObjective<=1.0e-15*fabs(lastObjective)) 
	goodMove=1;
      else
	goodMove=0;
    } else {
      maxDelta=1.0e10;
    }
    double maxGap=0.0;
    int numberSmaller=0;
    int numberSmaller2=0;
    int numberLarger=0;
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
      iColumn=listNonLinearColumn[jNon];
      maxDelta = CoinMax(maxDelta,
		     fabs(solution[iColumn]-saveSolution[iColumn]));
      if (goodMove>0) {
	if (last[0][jNon]*last[1][jNon]<0) {
	  // halve
	  trust[jNon] *= 0.5;
	  numberSmaller2++;
	} else {
	  if (last[0][jNon]==last[1][jNon]&&
	      last[0][jNon]==last[2][jNon])
	    trust[jNon] = CoinMin(1.5*trust[jNon],1.0e6); 
	  numberLarger++;
	}
      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
	trust[jNon] *= 0.2;
	numberSmaller++;
      }
      maxGap = CoinMax(maxGap,trust[jNon]);
    }
#ifdef CLP_DEBUG
	  if (handler_->logLevel()&32) 
	    std::cout<<"largest gap is "<<maxGap<<" "
		     <<numberSmaller+numberSmaller2<<" reduced ("
		     <<numberSmaller<<" badMove ), "
		     <<numberLarger<<" increased"<<std::endl;
#endif
    if (iPass>10000) {
      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
	trust[jNon] *=0.0001;
    }
    if(lastMoveStatus==-1&&goodMove==-1)
      goodMove=1; // to force solve
    if (goodMove>0) {
      double drop = lastObjective-objValue;
      handler_->message(CLP_SLP_ITER,messages_)
	<<iPass<<objValue-objectiveOffset
	<<drop<<maxDelta
	<<CoinMessageEol;
      if (iPass>20&&drop<1.0e-12*fabs(objValue))
	drop=0.999e-4; // so will exit
      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&theta<0.99999) {
	if (handler_->logLevel()>1) 
	  std::cout<<"Exiting as maxDelta < tolerance and small drop"<<std::endl;
	break;
      }
    } else if (!numberSmaller&&iPass>1) {
      if (handler_->logLevel()>1) 
	  std::cout<<"Exiting as all gaps small"<<std::endl;
	break;
    }
    if (!iPass)
      goodMove=1;
    targetDrop=0.0;
    double * r = this->dualColumnSolution();
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
      iColumn=listNonLinearColumn[jNon];
      columnLower[iColumn]=CoinMax(solution[iColumn]
			       -trust[jNon],
			       trueLower[jNon]);
      columnUpper[iColumn]=CoinMin(solution[iColumn]
			       +trust[jNon],
			       trueUpper[jNon]);
    }
    if (iPass) {
      // get reduced costs
      this->matrix()->transposeTimes(savePi,
				     this->dualColumnSolution());
      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	iColumn=listNonLinearColumn[jNon];
	double dj = objective[iColumn]-r[iColumn];
	r[iColumn]=dj;
	if (dj<-dualTolerance_) 
	  targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
	else if (dj>dualTolerance_)
	  targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
      }
    } else {
      memset(r,0,numberColumns*sizeof(double));
    }
#if 0
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
      iColumn=listNonLinearColumn[jNon];
      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
	columnLower[iColumn]=CoinMax(solution[iColumn],
				 trueLower[jNon]);
	columnUpper[iColumn]=CoinMin(solution[iColumn]
				 +trust[jNon],
				 trueUpper[jNon]);
      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
	columnLower[iColumn]=CoinMax(solution[iColumn]
				 -trust[jNon],
				 trueLower[jNon]);
	columnUpper[iColumn]=CoinMin(solution[iColumn],
				 trueUpper[jNon]);
      } else {
	columnLower[iColumn]=CoinMax(solution[iColumn]
				 -trust[jNon],
				 trueLower[jNon]);
	columnUpper[iColumn]=CoinMin(solution[iColumn]
				 +trust[jNon],
				 trueUpper[jNon]);
      }
    }
#endif
    if (goodMove>0) {
      memcpy(saveSolution,solution,numberColumns*sizeof(double));
      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
      memcpy(saveStatus,status_,numberRows+numberColumns);
#if MULTIPLE>2
      double * tempSol=saveSolutionM[0];
      for (jNon=0;jNon<MULTIPLE-1;jNon++) {
	saveSolutionM[jNon]=saveSolutionM[jNon+1];
      }
      saveSolutionM[MULTIPLE-1]=tempSol;
      memcpy(tempSol,solution,numberColumns*sizeof(double));
      
#endif
      
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) 
	std::cout<<"Pass - "<<iPass
		 <<", target drop is "<<targetDrop
		 <<std::endl;
#endif
      lastObjective = objValue;
      if (targetDrop<CoinMax(1.0e-8,CoinMin(1.0e-6,1.0e-6*fabs(objValue)))&&goodMove&&iPass>3) {
	if (handler_->logLevel()>1) 
	printf("Exiting on target drop %g\n",targetDrop);
	break;
      }
#ifdef CLP_DEBUG
      {
	double * r = this->dualColumnSolution();
	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	  iColumn=listNonLinearColumn[jNon];
	  if (handler_->logLevel()&32) 
	    printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
		   jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
		   r[iColumn],statusCheck[iColumn],columnLower[iColumn],
		   columnUpper[iColumn]);
	}
      }
#endif
      //setLogLevel(63);
      this->scaling(false);
      if (saveLogLevel==1)
	setLogLevel(0);
      ClpSimplex::primal(1);
      algorithm_=1;
      setLogLevel(saveLogLevel);
#ifdef CLP_DEBUG
      if (this->status()) {
	writeMps("xx.mps");
      }
#endif
      if (this->status()==1) {
	// not feasible ! - backtrack and exit
	// use safe solution
	memcpy(solution,safeSolution,numberColumns*sizeof(double));
	memcpy(saveSolution,solution,numberColumns*sizeof(double));
	memset(rowActivity_,0,numberRows_*sizeof(double));
	times(1.0,solution,rowActivity_);
	memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
	memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
	memcpy(status_,saveStatus,numberRows+numberColumns);
	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	  iColumn=listNonLinearColumn[jNon];
	  columnLower[iColumn]=CoinMax(solution[iColumn]
				       -trust[jNon],
				       trueLower[jNon]);
	  columnUpper[iColumn]=CoinMin(solution[iColumn]
				       +trust[jNon],
				       trueUpper[jNon]);
	}
	break;
      } else {
	// save in case problems
	memcpy(safeSolution,solution,numberColumns*sizeof(double));
      }
      if (problemStatus_==3)
	break; // probably user interrupt
      goodMove=1;
    } else {
      // bad pass - restore solution
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32) 
	printf("Backtracking\n");
#endif
      memcpy(solution,saveSolution,numberColumns*sizeof(double));
      memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
      memcpy(status_,saveStatus,numberRows+numberColumns);
      iPass--;
      assert (exitPass>0);
      goodMove=-1;
    }
  }
#if MULTIPLE>2
  for (jNon=0;jNon<MULTIPLE;jNon++) {
    delete [] saveSolutionM[jNon];
  }
#endif
  // restore solution
  memcpy(solution,saveSolution,numberColumns*sizeof(double));
  memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    iColumn=listNonLinearColumn[jNon];
    columnLower[iColumn]=CoinMax(solution[iColumn],
			     trueLower[jNon]);
    columnUpper[iColumn]=CoinMin(solution[iColumn],
			     trueUpper[jNon]);
  }
  delete [] markNonlinear;
  delete [] statusCheck;
  delete [] savePi;
  delete [] saveStatus;
  // redo objective
  memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
	 numberColumns*sizeof(double));
  ClpSimplex::primal(1);
  delete objective_;
  objective_=trueObjective;
  // redo values
  setDblParam(ClpObjOffset,objectiveOffset);
  objectiveValue_ += whichWay*offset;
  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    iColumn=listNonLinearColumn[jNon];
    columnLower[iColumn]= trueLower[jNon];
    columnUpper[iColumn]= trueUpper[jNon];
  }
  delete [] saveSolution;
  delete [] safeSolution;
  delete [] saveRowSolution;
  for (iPass=0;iPass<3;iPass++) 
    delete [] last[iPass];
  delete [] trust;
  delete [] trueUpper;
  delete [] trueLower;
  delete [] listNonLinearColumn;
  delete [] changeRegion;
  // temp
  //setLogLevel(63);
  return 0;
}
/* Primal algorithm for nonlinear constraints
   Using a semi-trust region approach as for pooling problem
   This is in because I have it lying around
   
*/
int 
ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints,
			       int numberPasses, double deltaTolerance)
{
  if (!numberConstraints) {
    // no nonlinear constraints - may be nonlinear objective
    return primalSLP(numberPasses,deltaTolerance);
  }
  // Are we minimizing or maximizing
  double whichWay=optimizationDirection();
  if (whichWay<0.0)
    whichWay=-1.0;
  else if (whichWay>0.0)
    whichWay=1.0;
  // check all matrix for odd rows is empty
  int iConstraint;
  int numberBad=0;
  CoinPackedMatrix * columnCopy = matrix();
  // Get a row copy in standard format
  CoinPackedMatrix copy;
  copy.reverseOrderedCopyOf(*columnCopy);
  // get matrix data pointers
  //const int * column = copy.getIndices();
  //const CoinBigIndex * rowStart = copy.getVectorStarts();
  const int * rowLength = copy.getVectorLengths(); 
  //const double * elementByRow = copy.getElements();
  int numberArtificials=0;
  // We could use nonlinearcost to do segments - maybe later
#define SEGMENTS 3 
  // Penalties may be adjusted by duals
  // Both these should be modified depending on problem
  // Possibly start with big bounds
  //double penalties[]={1.0e-3,1.0e7,1.0e9};
  double penalties[]={1.0e7,1.0e8,1.0e9};
  double bounds[] = {1.0e-2,1.0e2,COIN_DBL_MAX};
  // see how many extra we need
  CoinBigIndex numberExtra=0;
  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
    ClpConstraint * constraint = constraints[iConstraint];
    int iRow = constraint->rowNumber();
    assert (iRow>=0);
    int n = constraint->numberCoefficients() -rowLength[iRow];
    numberExtra += n;
    if (iRow>=numberRows_)
      numberBad++;
    else if (rowLength[iRow]&&n)
      numberBad++;
    if (rowLower_[iRow]>-1.0e20)
      numberArtificials++;
    if (rowUpper_[iRow]<1.0e20)
      numberArtificials++;
  }
  if (numberBad)
    return numberBad;
  ClpObjective * trueObjective = NULL;
  if (objective_->type()>=2) {
    // Replace objective
    trueObjective = objective_;
    objective_=new ClpLinearObjective(NULL,numberColumns_);
  }
  ClpSimplex newModel(*this);
  // we can put back true objective
  if (trueObjective) {
    // Replace objective
    delete objective_;
    objective_=trueObjective;
  }
  int numberColumns2 = numberColumns_;
  int numberSmallGap = numberArtificials;
  if (numberArtificials) {
    numberArtificials *= SEGMENTS;
    numberColumns2 += numberArtificials;
    int * addStarts = new int [numberArtificials+1];
    int * addRow = new int[numberArtificials];
    double * addElement = new double[numberArtificials];
    double * addUpper = new double[numberArtificials];
    addStarts[0]=0;
    double * addCost = new double [numberArtificials];
    numberArtificials=0;
    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
      ClpConstraint * constraint = constraints[iConstraint];
      int iRow = constraint->rowNumber();
      if (rowLower_[iRow]>-1.0e20) {
	for (int k=0;k<SEGMENTS;k++) {
	  addRow[numberArtificials]=iRow;
	  addElement[numberArtificials]=1.0;
	  addCost[numberArtificials]=penalties[k];
	  addUpper[numberArtificials]=bounds[k];
	  numberArtificials++;
	  addStarts[numberArtificials]=numberArtificials;
	}
      }
      if (rowUpper_[iRow]<1.0e20) {
	for (int k=0;k<SEGMENTS;k++) {
	  addRow[numberArtificials]=iRow;
	  addElement[numberArtificials]=-1.0;
	  addCost[numberArtificials]=penalties[k];
	  addUpper[numberArtificials]=bounds[k];
	  numberArtificials++;
	  addStarts[numberArtificials]=numberArtificials;
	}
      }
    }
    newModel.addColumns(numberArtificials,NULL,addUpper,addCost,
		       addStarts,addRow,addElement);
    delete [] addStarts;
    delete [] addRow;
    delete [] addElement;
    delete [] addUpper;
    delete [] addCost;
    //    newModel.primal(1);
  }
  // find nonlinear columns
  int * listNonLinearColumn = new int [numberColumns_+numberSmallGap];
  char * mark = new char[numberColumns_];
  memset(mark,0,numberColumns_);
  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
    ClpConstraint * constraint = constraints[iConstraint];
    constraint->markNonlinear(mark);
  }
  if (trueObjective)
    trueObjective->markNonlinear(mark);
  int iColumn;
  int numberNonLinearColumns=0;
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    if (mark[iColumn])
      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
  }
  // and small gap artificials
  for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) {
    listNonLinearColumn[numberNonLinearColumns++]=iColumn;
  }
  // build row copy of matrix with space for nonzeros
  // Get column copy
  columnCopy = newModel.matrix();
  copy.reverseOrderedCopyOf(*columnCopy);
  // get matrix data pointers
  const int * column = copy.getIndices();
  const CoinBigIndex * rowStart = copy.getVectorStarts();
  rowLength = copy.getVectorLengths(); 
  const double * elementByRow = copy.getElements();
  numberExtra +=copy.getNumElements();
  CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1];
  int * newColumn = new int[numberExtra];
  double * newElement = new double[numberExtra];
  newStarts[0]=0;
  int * backRow = new int [numberRows_];
  int iRow;
  for (iRow=0;iRow<numberRows_;iRow++) 
    backRow[iRow]=-1;
  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
    ClpConstraint * constraint = constraints[iConstraint];
    iRow = constraint->rowNumber();
    backRow[iRow]=iConstraint;
  }
  numberExtra=0;
  for (iRow=0;iRow<numberRows_;iRow++) {
    if (backRow[iRow]<0) {
      // copy normal 
      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
	   j++) {
	newColumn[numberExtra]=column[j];
	newElement[numberExtra++] = elementByRow[j];
      }
    } else {
      ClpConstraint * constraint = constraints[backRow[iRow]];
      assert(iRow == constraint->rowNumber());
      int numberArtificials=0;
      if (rowLower_[iRow]>-1.0e20)
	numberArtificials += SEGMENTS;
      if (rowUpper_[iRow]<1.0e20)
	numberArtificials += SEGMENTS;
      if (numberArtificials==rowLength[iRow]) {
	// all possible
	memset(mark,0,numberColumns_);
	constraint->markNonzero(mark);
	for (int k=0;k<numberColumns_;k++) {
	  if (mark[k]) {
	    newColumn[numberExtra]=k;
	    newElement[numberExtra++] = 1.0;
	  }
	}
	// copy artificials
	for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
	     j++) {
	  newColumn[numberExtra]=column[j];
	  newElement[numberExtra++] = elementByRow[j];
	}
      } else {
	// there already
	// copy
	for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
	     j++) {
	  newColumn[numberExtra]=column[j];
	  assert (elementByRow[j]);
	  newElement[numberExtra++] = elementByRow[j];
	}
      }
    }
    newStarts[iRow+1]=numberExtra;
  }
  delete [] backRow;
  CoinPackedMatrix saveMatrix(false,numberColumns2,numberRows_,
			      numberExtra,newElement,newColumn,newStarts,NULL,0.0,0.0);
  delete [] newStarts;
  delete [] newColumn;
  delete [] newElement;
  delete [] mark;
  // get feasible
  if (whichWay<0.0) {
    newModel.setOptimizationDirection(1.0);
    double * objective = newModel.objective();
    for (int iColumn=0;iColumn<numberColumns_;iColumn++) 
      objective[iColumn] = - objective[iColumn];
  }
  newModel.primal(1);
  // still infeasible
  if (newModel.problemStatus()==1) {
    delete [] listNonLinearColumn;
    return 0;
  } else if (newModel.problemStatus()==2) {
    // unbounded - add bounds
    double * columnLower = newModel.columnLower();
    double * columnUpper = newModel.columnUpper();
    for (int i=0;i<numberColumns_;i++) {
      columnLower[i]= CoinMax(-1.0e8,columnLower[i]);
      columnUpper[i]= CoinMin(1.0e8,columnUpper[i]);
    }
    newModel.primal(1);
  }
  int numberRows = newModel.numberRows();
  double * columnLower = newModel.columnLower();
  double * columnUpper = newModel.columnUpper();
  double * objective = newModel.objective();
  double * rowLower = newModel.rowLower();
  double * rowUpper = newModel.rowUpper(); 
  double * solution = newModel.primalColumnSolution();
  int jNon;
  int * last[3];
  
  double * trust = new double[numberNonLinearColumns];
  double * trueLower = new double[numberNonLinearColumns];
  double * trueUpper = new double[numberNonLinearColumns];
  double objectiveOffset;
  double objectiveOffset2;
  getDblParam(ClpObjOffset,objectiveOffset);
  objectiveOffset *= whichWay;
  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    iColumn=listNonLinearColumn[jNon];
    double upper = columnUpper[iColumn];
    double lower = columnLower[iColumn];
    if (solution[iColumn]<lower)
      solution[iColumn]=lower;
    else if (solution[iColumn]>upper)
      solution[iColumn]=upper;
#if 0
    double large = CoinMax(1000.0,10.0*fabs(solution[iColumn]));
    if (upper>1.0e10)
      upper = solution[iColumn]+large;
    if (lower<-1.0e10)
      lower = solution[iColumn]-large;
#else
    upper = solution[iColumn]+0.5;
    lower = solution[iColumn]-0.5;
#endif
    //columnUpper[iColumn]=upper;
    trust[jNon]=0.05*(1.0+upper-lower);
    trueLower[jNon]=columnLower[iColumn];
    //trueUpper[jNon]=upper;
    trueUpper[jNon]=columnUpper[iColumn];
  }
  bool tryFix=false;
  int iPass;
  double lastObjective=1.0e31;
  double lastGoodObjective=1.0e31;
  double * bestSolution = NULL;
  double * saveSolution = new double [numberColumns2+numberRows];
  char * saveStatus = new char [numberColumns2+numberRows];
  double targetDrop=1.0e31;
  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
  for (iPass=0;iPass<3;iPass++) {
    last[iPass]=new int[numberNonLinearColumns];
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
      last[iPass][jNon]=0;
  }
  int numberZeroPasses=0;
  bool zeroTargetDrop=false;
  double * gradient = new double [numberColumns_];
  bool goneFeasible=false;
  // keep sum of artificials
#define KEEP_SUM 5
  double sumArt[KEEP_SUM];
  for (jNon=0;jNon<KEEP_SUM;jNon++)
    sumArt[jNon]=COIN_DBL_MAX;
#define SMALL_FIX 0.0
  for (iPass=0;iPass<numberPasses;iPass++) {
    objectiveOffset2 = objectiveOffset;
    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
      iColumn=listNonLinearColumn[jNon];
      if (solution[iColumn]<trueLower[jNon])
	solution[iColumn]=trueLower[jNon];
      else if (solution[iColumn]>trueUpper[jNon])
	solution[iColumn]=trueUpper[jNon];
      columnLower[iColumn]=CoinMax(solution[iColumn]
				   -trust[jNon],
				   trueLower[jNon]);
      if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
	columnLower[iColumn]=SMALL_FIX;
      columnUpper[iColumn]=CoinMin(solution[iColumn]
				   +trust[jNon],
				   trueUpper[jNon]);
      if (!trueLower[jNon]) 
	columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
				       columnLower[iColumn]+SMALL_FIX);
      if (!trueLower[jNon]&&tryFix&&
	  columnLower[iColumn]==SMALL_FIX&&
	  columnUpper[iColumn]<3.0*SMALL_FIX) {
	columnLower[iColumn]=0.0;
	solution[iColumn]=0.0;
	columnUpper[iColumn]=0.0;
	printf("fixing %d\n",iColumn);
      }
    }
    // redo matrix
    double offset;
    CoinPackedMatrix newMatrix(saveMatrix);
    // get matrix data pointers
    column = newMatrix.getIndices();
    rowStart = newMatrix.getVectorStarts();
    rowLength = newMatrix.getVectorLengths(); 
    // make sure x updated
    if (numberConstraints)
      constraints[0]->newXValues();
    else
      trueObjective->newXValues();
    double * changeableElement = newMatrix.getMutableElements();
    if (trueObjective) {
      memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
	     numberColumns_*sizeof(double));
    } else {
      memcpy(objective,objective_->gradient(this,solution,offset,true,2),
	     numberColumns_*sizeof(double));
    }
    if (whichWay<0.0) {
      for (int iColumn=0;iColumn<numberColumns_;iColumn++) 
	objective[iColumn] = - objective[iColumn];
    }
    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
      ClpConstraint * constraint = constraints[iConstraint];
      int iRow = constraint->rowNumber();
      double functionValue;
#ifndef NDEBUG
      int numberErrors = 
#endif
      constraint->gradient(&newModel,solution,gradient,functionValue,offset);
      assert (!numberErrors);
      // double dualValue = newModel.dualRowSolution()[iRow];
      int numberCoefficients = constraint->numberCoefficients();
      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+numberCoefficients;j++) {
	int iColumn = column[j];
	changeableElement[j] = gradient[iColumn];
	//objective[iColumn] -= dualValue*gradient[iColumn];
	gradient[iColumn]=0.0;
      }
      for (int k=0;k<numberColumns_;k++)
	assert (!gradient[k]);
      if (rowLower_[iRow]>-1.0e20)
	rowLower[iRow] = rowLower_[iRow] - offset;
      if (rowUpper_[iRow]<1.0e20)
	rowUpper[iRow] = rowUpper_[iRow] - offset;
    }
    // Replace matrix
    // Get a column copy in standard format
    CoinPackedMatrix * columnCopy = new CoinPackedMatrix();;
    columnCopy->reverseOrderedCopyOf(newMatrix);
    newModel.replaceMatrix(columnCopy,true);
    // solve
    newModel.primal(1);
    if (newModel.status()==1) { 
      // Infeasible!
      newModel.allSlackBasis();
      newModel.primal();
      newModel.writeMps("infeas.mps");
      assert(!newModel.status());
    }
    double sumInfeas=0.0;
    int numberInfeas=0;
    for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
      if (solution[iColumn]>1.0e-8) {
	numberInfeas++;
	sumInfeas += solution[iColumn];
      }
    }
    printf("%d artificial infeasibilities - summing to %g\n",
	   numberInfeas,sumInfeas);
    for (jNon=0;jNon<KEEP_SUM-1;jNon++)
      sumArt[jNon]=sumArt[jNon+1];
    sumArt[KEEP_SUM-1]=sumInfeas;
    if (sumInfeas>0.01&&sumInfeas*1.1>sumArt[0]&&penalties[1]<1.0e7) {
      // not doing very well - increase - be more sophisticated later
      lastObjective=COIN_DBL_MAX;
      for (jNon=0;jNon<KEEP_SUM;jNon++)
	sumArt[jNon]=COIN_DBL_MAX;
      //for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) {
      //objective[iColumn+1] *= 1.5;
      //}
      penalties[1] *= 1.5;
      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
	if (trust[jNon]>0.1)
	  trust[jNon] *= 2.0;
	else
	  trust[jNon] = 0.1;
    }
    if (sumInfeas<0.001&&!goneFeasible) {
      goneFeasible=true;
      penalties[0]=1.0e-3;
      penalties[1]=1.0e6;
      printf("Got feasible\n");
    }
    double infValue=0.0;
    double objValue=0.0;
    // make sure x updated
    if (numberConstraints)
      constraints[0]->newXValues();
    else
      trueObjective->newXValues();
    if (trueObjective) {
      objValue=trueObjective->objectiveValue(this,solution);
      printf("objective offset %g\n",offset);
      objectiveOffset2 =objectiveOffset+offset; // ? sign
      newModel.setObjectiveOffset(objectiveOffset2);
    } else {
      objValue=objective_->objectiveValue(this,solution);
    }
    objValue *= whichWay;
    double infPenalty=0.0;
    // This penalty is for target drop
    double infPenalty2=0.0;
    //const int * row = columnCopy->getIndices();
    //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
    //const int * columnLength = columnCopy->getVectorLengths(); 
    //const double * element = columnCopy->getElements();
    double * cost = newModel.objective();
    column = newMatrix.getIndices();
    rowStart = newMatrix.getVectorStarts();
    rowLength = newMatrix.getVectorLengths(); 
    elementByRow = newMatrix.getElements();
    int jColumn = numberColumns_;
    double objectiveAdjustment=0.0;
    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
      ClpConstraint * constraint = constraints[iConstraint];
      int iRow = constraint->rowNumber();
      double functionValue=constraint->functionValue(this,solution);
      double dualValue = newModel.dualRowSolution()[iRow];
      if (numberConstraints<-50)
	printf("For row %d current value is %g (row activity %g) , dual is %g\n",iRow,functionValue,
	       newModel.primalRowSolution()[iRow],
	       dualValue);
      double movement = newModel.primalRowSolution()[iRow]+constraint->offset();
      movement = fabs((movement-functionValue)*dualValue);
      infPenalty2 += movement;
      double sumOfActivities=0.0;
      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
	int iColumn = column[j];
	sumOfActivities += fabs(solution[iColumn]*elementByRow[j]); 
      }
      if (rowLower_[iRow]>-1.0e20) {
	if (functionValue<rowLower_[iRow]-1.0e-5) {
	  double infeasibility = rowLower_[iRow]-functionValue;
	  double thisPenalty=0.0;
	  infValue += infeasibility;
	  double boundMultiplier = 1.0;
	  if (sumOfActivities<0.001)
	    boundMultiplier = 0.1;
	  else if (sumOfActivities>100.0)
	    boundMultiplier=10.0;
	  int k;
	  assert (dualValue>=-1.0e-5);
	  dualValue = CoinMax(dualValue,0.0);
	  for ( k=0;k<SEGMENTS;k++) {
	    if (infeasibility<=0)
	      break;
	    double thisPart = CoinMin(infeasibility,bounds[k]);
	    thisPenalty += thisPart*cost[jColumn+k];
	    infeasibility -= thisPart;
	  }
	  infeasibility = functionValue-rowUpper_[iRow];
	  double newPenalty=0.0;
	  for ( k=0;k<SEGMENTS;k++) {
	    double thisPart = CoinMin(infeasibility,bounds[k]);
	    cost[jColumn+k] = CoinMax(penalties[k],dualValue+1.0e-3);
	    newPenalty += thisPart*cost[jColumn+k];
	    infeasibility -= thisPart;
	  }
	  infPenalty += thisPenalty;
	  objectiveAdjustment += CoinMax(0.0,newPenalty-thisPenalty);
	}
	jColumn += SEGMENTS;
      }
      if (rowUpper_[iRow]<1.0e20) {
	if (functionValue>rowUpper_[iRow]+1.0e-5) {
	  double infeasibility = functionValue-rowUpper_[iRow];
	  double thisPenalty=0.0;
	  infValue += infeasibility;
	  double boundMultiplier = 1.0;
	  if (sumOfActivities<0.001)
	    boundMultiplier = 0.1;
	  else if (sumOfActivities>100.0)
	    boundMultiplier=10.0;
	  int k;
	  dualValue = -dualValue;
	  assert (dualValue>=-1.0e-5);
	  dualValue = CoinMax(dualValue,0.0);
	  for ( k=0;k<SEGMENTS;k++) {
	    if (infeasibility<=0)
	      break;
	    double thisPart = CoinMin(infeasibility,bounds[k]);
	    thisPenalty += thisPart*cost[jColumn+k];
	    infeasibility -= thisPart;
	  }
	  infeasibility = functionValue-rowUpper_[iRow];
	  double newPenalty=0.0;
	  for ( k=0;k<SEGMENTS;k++) {
	    double thisPart = CoinMin(infeasibility,bounds[k]);
	    cost[jColumn+k] = CoinMax(penalties[k],dualValue+1.0e-3);
	    newPenalty += thisPart*cost[jColumn+k];
	    infeasibility -= thisPart;
	  }
	  infPenalty += thisPenalty;
	  objectiveAdjustment += CoinMax(0.0,newPenalty-thisPenalty);
	}
	jColumn += SEGMENTS;
      }
    }
    // adjust last objective value
    lastObjective += objectiveAdjustment;
    if (infValue)
      printf("Sum infeasibilities %g - penalty %g ",infValue,infPenalty);
    if (objectiveOffset2)
      printf("offset2 %g ",objectiveOffset2);
    objValue -= objectiveOffset2;
    printf("True objective %g or maybe %g (with penalty %g) -pen2 %g %g\n",objValue,
	   objValue+objectiveOffset2,objValue+objectiveOffset2+infPenalty,infPenalty2,penalties[1]);
    double useObjValue = objValue+objectiveOffset2+infPenalty;
    objValue += infPenalty+infPenalty2;
    objValue = useObjValue;
    if (iPass) {
      double drop = lastObjective-objValue;
      std::cout<<"True drop was "<<drop<<std::endl;
      if (drop<-0.05*fabs(objValue)-1.0e-4) {
	// pretty bad - go back and halve
	memcpy(solution,saveSolution,numberColumns2*sizeof(double));
	memcpy(newModel.primalRowSolution(),saveSolution+numberColumns2,
	       numberRows*sizeof(double));
	memcpy(newModel.statusArray(),saveStatus,
	       numberColumns2+numberRows);
	for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
	  if (trust[jNon]>0.1)
	    trust[jNon] *= 0.5;
	  else
	    trust[jNon] *= 0.9;
	
	printf("** Halving trust\n");
	objValue=lastObjective;
	continue;
      } else if ((iPass%25)==-1) {
	for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
	  trust[jNon] *= 2.0;
	printf("** Doubling trust\n");
      }
      int * temp=last[2];
      last[2]=last[1];
      last[1]=last[0];
      last[0]=temp;
      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	iColumn=listNonLinearColumn[jNon];
	double change = solution[iColumn]-saveSolution[iColumn];
	if (change<-1.0e-5) {
	  if (fabs(change+trust[jNon])<1.0e-5) 
	    temp[jNon]=-1;
	  else
	    temp[jNon]=-2;
	} else if(change>1.0e-5) {
	  if (fabs(change-trust[jNon])<1.0e-5) 
	    temp[jNon]=1;
	  else
	    temp[jNon]=2;
	} else {
	  temp[jNon]=0;
	}
      } 
      double maxDelta=0.0;
      double smallestTrust=1.0e31;
      double smallestNonLinearGap=1.0e31;
      double smallestGap=1.0e31;
      bool increasing=false;
      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
	double gap = columnUpper[iColumn]-columnLower[iColumn];
	assert (gap>=0.0);
	if (gap)
	  smallestGap = CoinMin(smallestGap,gap);
      }
      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	iColumn=listNonLinearColumn[jNon];
	double gap = columnUpper[iColumn]-columnLower[iColumn];
	assert (gap>=0.0);
	if (gap) {
	  smallestNonLinearGap = CoinMin(smallestNonLinearGap,gap);
	  if (gap<1.0e-7&&iPass==1) {
	    printf("Small gap %d %d %g %g %g\n",
		   jNon,iColumn,columnLower[iColumn],columnUpper[iColumn],
		   gap);
	    //trueUpper[jNon]=trueLower[jNon];
	    //columnUpper[iColumn]=columnLower[iColumn];
	  }
	}
	maxDelta = CoinMax(maxDelta,
		       fabs(solution[iColumn]-saveSolution[iColumn]));
	if (last[0][jNon]*last[1][jNon]<0) {
	  // halve
	  if (trust[jNon]>1.0)
	    trust[jNon] *= 0.5;
	  else
	    trust[jNon] *= 0.7;
	} else {
	  // ? only increase if +=1 ?
	  if (last[0][jNon]==last[1][jNon]&&
	      last[0][jNon]==last[2][jNon]&&
	      last[0][jNon]) {
	    trust[jNon] *= 1.8; 
	    increasing=true;
	  }
	}
	smallestTrust = CoinMin(smallestTrust,trust[jNon]);
      }
      std::cout<<"largest delta is "<<maxDelta
	       <<", smallest trust is "<<smallestTrust
	       <<", smallest gap is "<<smallestGap
	       <<", smallest nonlinear gap is "<<smallestNonLinearGap
	       <<std::endl;
      if (iPass>200) {
	//double useObjValue = objValue+objectiveOffset2+infPenalty;
	if (useObjValue+1.0e-4>	lastGoodObjective&&iPass>250) {
	  std::cout<<"Exiting as objective not changing much"<<std::endl;
	  break;
	} else if (useObjValue<lastGoodObjective) {
	  lastGoodObjective = useObjValue;
	  if (!bestSolution)
	    bestSolution = new double [numberColumns2];
	  memcpy(bestSolution,solution,numberColumns2*sizeof(double));
	}
      }
      if (maxDelta<deltaTolerance&&!increasing&&iPass>100) {
	numberZeroPasses++;
	if (numberZeroPasses==3) {
	  if (tryFix) {
	    std::cout<<"Exiting"<<std::endl;
	    break;
	  } else {
	    tryFix=true;
	    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
	      trust[jNon] = CoinMax(trust[jNon],1.0e-1);
	    numberZeroPasses=0;
	  }
	}
      } else {
	numberZeroPasses=0;
      }
    }
    memcpy(saveSolution,solution,numberColumns2*sizeof(double));
    memcpy(saveSolution+numberColumns2,newModel.primalRowSolution(),
	   numberRows*sizeof(double));
    memcpy(saveStatus,newModel.statusArray(),
	   numberColumns2+numberRows);
    
    targetDrop=infPenalty+infPenalty2;
    if (iPass) {
      // get reduced costs
      const double * pi = newModel.dualRowSolution();
      newModel.matrix()->transposeTimes(pi,
					newModel.dualColumnSolution());
      double * r = newModel.dualColumnSolution();
      for (iColumn=0;iColumn<numberColumns_;iColumn++) 
	r[iColumn] = objective[iColumn]-r[iColumn];
      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
	iColumn=listNonLinearColumn[jNon];
	double dj = r[iColumn];
	if (dj<-1.0e-6) {
	  double drop = -dj*(columnUpper[iColumn]-solution[iColumn]);
	  //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
	  //double drop2 = -dj*(upper-solution[iColumn]);
#if 0
	  if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2&&iPass>100))
	    printf("Big drop %d %g %g %g %g T %g\n",
		   iColumn,columnLower[iColumn],solution[iColumn],
		   columnUpper[iColumn],dj,trueUpper[jNon]);
#endif
	  targetDrop += drop;
	  if (dj<-1.0e-1&&trust[jNon]<1.0e-3
	      &&trueUpper[jNon]-solution[iColumn]>1.0e-2) {
	    trust[jNon] *= 1.5;
	    //printf("Increasing trust on %d to %g\n",
	    //     iColumn,trust[jNon]);
	  }
	} else if (dj>1.0e-6) {
	  double drop = -dj*(columnLower[iColumn]-solution[iColumn]);
	  //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
	  //double drop2 = -dj*(lower-solution[iColumn]);
#if 0
	  if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2))
	    printf("Big drop %d %g %g %g %g T %g\n",
		   iColumn,columnLower[iColumn],solution[iColumn],
		   columnUpper[iColumn],dj,trueLower[jNon]);
#endif
	  targetDrop += drop;
	  if (dj>1.0e-1&&trust[jNon]<1.0e-3
	      &&solution[iColumn]-trueLower[jNon]>1.0e-2) {
	    trust[jNon] *= 1.5;
	    printf("Increasing trust on %d to %g\n",
		   iColumn,trust[jNon]);
	  }
	}
      }
    }
    std::cout<<"Pass - "<<iPass
	     <<", target drop is "<<targetDrop
	     <<std::endl;
    if (iPass>1&&targetDrop<1.0e-5&&zeroTargetDrop)
      break;
    if (iPass>1&&targetDrop<1.0e-5)
      zeroTargetDrop = true;
    else
      zeroTargetDrop = false;
    //if (iPass==5)
    //newModel.setLogLevel(63);
    lastObjective = objValue;
    // take out when ClpPackedMatrix changed
    //newModel.scaling(false);
#if 0
    CoinMpsIO writer;
    writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
		      newModel.getColLower(), newModel.getColUpper(),
		      newModel.getObjCoefficients(),
		      (const char*) 0 /*integrality*/,
		      newModel.getRowLower(), newModel.getRowUpper(),
		      NULL,NULL);
    writer.writeMps("xx.mps");
#endif
  }
  delete [] saveSolution;
  delete [] saveStatus;
  for (iPass=0;iPass<3;iPass++) 
    delete [] last[iPass];
  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    iColumn=listNonLinearColumn[jNon];
    columnLower[iColumn]=trueLower[jNon];
    columnUpper[iColumn]=trueUpper[jNon];
  }
  delete [] trust;
  delete [] trueUpper;
  delete [] trueLower;
  objectiveValue_=newModel.objectiveValue();
  if (bestSolution) {
    memcpy(solution,bestSolution,numberColumns2*sizeof(double));
    delete [] bestSolution;
    printf("restoring objective of %g\n",lastGoodObjective);
    objectiveValue_=lastGoodObjective;
  }
  // Simplest way to get true row activity ?
  double * rowActivity = newModel.primalRowSolution();
  for (iRow=0;iRow<numberRows;iRow++) {
    double difference;
    if (fabs(rowLower_[iRow])<fabs(rowUpper_[iRow]))
      difference = rowLower_[iRow]-rowLower[iRow];
    else
      difference = rowUpper_[iRow]-rowUpper[iRow];
    rowLower[iRow]=rowLower_[iRow];
    rowUpper[iRow]=rowUpper_[iRow];
    if (difference) {
      if (numberRows<50)
	printf("For row %d activity changes from %g to %g\n",
	       iRow,rowActivity[iRow],rowActivity[iRow]+difference);
      rowActivity[iRow]+= difference;
    }
  }
  delete [] listNonLinearColumn;
  delete [] gradient;
  printf("solution still in newModel - do objective etc!\n");
  numberIterations_ =newModel.numberIterations();
  problemStatus_ =newModel.problemStatus();
  secondaryStatus_ =newModel.secondaryStatus();
  memcpy(columnActivity_,newModel.primalColumnSolution(),numberColumns_*sizeof(double));
  // should do status region
  CoinZeroN(rowActivity_,numberRows_);
  matrix_->times(1.0,columnActivity_,rowActivity_);
  return 0;
}