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

#include "CoinPragma.hpp"

#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "ClpHelperFunctions.hpp"
#include "ClpSimplexNonlinear.hpp"
#include "ClpSimplexOther.hpp"
#include "ClpSimplexDual.hpp"
#include "ClpFactorization.hpp"
#include "ClpNonLinearCost.hpp"
#include "ClpLinearObjective.hpp"
#include "ClpConstraint.hpp"
#include "ClpPresolve.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_ && !scalingFlag_) && 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 < SHORT_REGION; 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);
    delete nonLinearCost_;
    nonLinearCost_ = new ClpNonLinearCost(this);
    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
    CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_);
    CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_);
    CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_);
    // 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
            CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_);
            CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_);
            CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_);
            // 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;
                         */
          unPerturb();
          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_];
          CoinMemcpyN(dual_, numberRows_, ray_);
          // 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(CoinMin(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
    CoinMemcpyN(status_, (numberColumns_ + numberRows_), saveStatus_);
    CoinMemcpyN(rowActivityWork_, numberRows_, savedSolution_ + numberColumns_);
    CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
  }
  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 == 2 && sequenceIn_ < 0) {
        // does not look good
        double currentObj;
        double thetaObj;
        double predictedObj;
        objective_->stepLength(this, solution_, solution_, 0.0,
          currentObj, thetaObj, predictedObj);
        if (currentObj == predictedObj) {
#ifdef CLP_INVESTIGATE
          printf("looks bad - no change in obj %g\n", currentObj);
#endif
          if (factorization_->pivots())
            result = 3;
          else
            problemStatus_ = 0;
        }
      }
      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, thetaObj, predictedObj);
        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;
#if CLP_DEBUG > 0
        int iLargest = -1;
#endif
        for (iRow = 0; iRow < numberRows_; iRow++) {
          double value = fabs(rhs[iRow]);
          rhs[iRow] = 0.0;
          if (value > largest) {
            largest = value;
#if CLP_DEBUG > 0
            iLargest = iRow;
#endif
          }
        }
#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 > CoinMax(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);
        CoinMemcpyN(work, numberTotal, dArray);
        CoinMemcpyN(work, numberTotal, dArray2);
        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;
#ifdef CLP_DEBUG
          int kPivot = -1;
#endif
          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) {
#ifdef CLP_DEBUG
                kPivot = iPivot;
#endif
                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) {
#ifdef CLP_DEBUG
                kPivot = iPivot;
#endif
                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 && false) {
            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 {
          }
          CoinMemcpyN(zz, number2 * number2, zzz);
          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
          CoinMemcpyN(saveY2, numberNonBasic, saveY);
          CoinMemcpyN(saveS2, numberNonBasic, saveS);
        }
#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
      objTheta2 = CoinMin(objTheta2, 1.0e29);
#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_ < 1.0e-10))
        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;
        } else if (sequenceIn_ < 0 && nTotalPasses > 10) {
          if (objTheta < 1.0e-10) {
            finished = true;
            //printf("zero move\n");
            break;
          }
        }
      }
#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 (%d) passes\n", nPasses, nTotalPasses);
#endif
    if (nTotalPasses >= 1000 || (nTotalPasses > 10 && sequenceIn_ < 0 && theta_ < 1.0e-10))
      returnCode = 2;
    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
    }
    if (theta_ >= 1.0e30) { // odd
      returnCode=2;
      break;
    }
    // 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) {
	assert (sequenceIn_>=0);
        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 = randomNumberGenerator_.randomDouble();
          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_];
  // 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(CoinMin(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;
}
// May use a cut approach for solving any LP
int ClpSimplexNonlinear::primalDualCuts(char *rowsIn, int startUp, int algorithm)
{
  if (!rowsIn) {
    if (algorithm > 0)
      return ClpSimplex::primal(startUp);
    else
      //return static_cast<ClpSimplexDual *>(static_cast<ClpSimplex *>(this))->dual(startUp);
      return ClpSimplex::dual(startUp);
  } else {
    int numberUsed = 0;
    int rowsThreshold = CoinMax(100, numberRows_ / 2);
    //int rowsTry=CoinMax(50,numberRows_/4);
    // Just add this number of rows each time in small problem
    int smallNumberRows = 2 * numberColumns_;
    smallNumberRows = CoinMin(smallNumberRows, numberRows_ / 20);
    // We will need arrays to choose rows to add
    double *weight = new double[numberRows_];
    int *sort = new int[numberRows_ + numberColumns_];
    int *whichColumns = sort + numberRows_;
    int numberSort = 0;
    for (int i = 0; i < numberRows_; i++) {
      if (rowsIn[i])
        numberUsed++;
    }
    if (numberUsed) {
      // normal
    } else {
      // If non slack use that information
      int numberBinding = 0;
      numberPrimalInfeasibilities_ = 0;
      sumPrimalInfeasibilities_ = 0.0;
      memset(rowActivity_, 0, numberRows_ * sizeof(double));
      times(1.0, columnActivity_, rowActivity_);
      for (int i = 0; i < numberRows_; i++) {
        double lowerDifference = rowActivity_[i] - rowLower_[i];
        double upperDifference = rowActivity_[i] - rowUpper_[i];
        if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
          numberPrimalInfeasibilities_++;
          if (lowerDifference < 0.0)
            sumPrimalInfeasibilities_ -= lowerDifference;
          else
            sumPrimalInfeasibilities_ += upperDifference;
          rowsIn[i] = 1;
        } else if (getRowStatus(i) != basic) {
          numberBinding++;
          rowsIn[i] = 1;
        }
      }
      if (numberBinding < rowsThreshold) {
        // try
      } else {
        // random?? - or give up
        // Set up initial list
        numberSort = 0;
        for (int iRow = 0; iRow < numberRows_; iRow++) {
          weight[iRow] = 1.123e50;
          if (rowLower_[iRow] == rowUpper_[iRow]) {
            sort[numberSort++] = iRow;
            weight[iRow] = 0.0;
          }
        }
        numberSort /= 2;
        // and pad out with random rows
        double ratio = ((double)(smallNumberRows - numberSort)) / ((double)numberRows_);
        for (int iRow = 0; iRow < numberRows_; iRow++) {
          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
            sort[numberSort++] = iRow;
        }
        // sort
        CoinSort_2(weight, weight + numberRows_, sort);
        numberSort = CoinMin(numberRows_, smallNumberRows);
        memset(rowsIn, 0, numberRows_);
        for (int iRow = 0; iRow < numberSort; iRow++)
          rowsIn[sort[iRow]] = 1;
      }
    }
    // see if feasible
    numberPrimalInfeasibilities_ = 0;
    memset(rowActivity_, 0, numberRows_ * sizeof(double));
    times(1.0, columnActivity_, rowActivity_);
    for (int i = 0; i < numberRows_; i++) {
      double lowerDifference = rowActivity_[i] - rowLower_[i];
      double upperDifference = rowActivity_[i] - rowUpper_[i];
      if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
        if (lowerDifference < 0.0)
          sumPrimalInfeasibilities_ -= lowerDifference;
        else
          sumPrimalInfeasibilities_ += upperDifference;
        numberPrimalInfeasibilities_++;
      }
    }
    printf("Initial infeasibilities - %g (%d)\n",
      sumPrimalInfeasibilities_, numberPrimalInfeasibilities_);
    // Just do this number of passes
    int maxPass = 50;
    // And take out slack rows until this pass
    int takeOutPass = 30;
    int iPass;

    const CoinBigIndex *start = this->clpMatrix()->getVectorStarts();
    const int *length = this->clpMatrix()->getVectorLengths();
    const int *row = this->clpMatrix()->getIndices();
    problemStatus_ = 1;
    for (iPass = 0; iPass < maxPass; iPass++) {
      printf("Start of pass %d\n", iPass);
      int numberSmallColumns = 0;
      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
        int n = 0;
        for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
          int iRow = row[j];
          if (rowsIn[iRow])
            n++;
        }
        if (n)
          whichColumns[numberSmallColumns++] = iColumn;
      }
      numberSort = 0;
      for (int i = 0; i < numberRows_; i++) {
        if (rowsIn[i])
          sort[numberSort++] = i;
      }
      // Create small problem
      ClpSimplex small(this, numberSort, sort, numberSmallColumns, whichColumns);
      printf("Small model has %d rows, %d columns and %d elements\n",
        small.numberRows(), small.numberColumns(), small.getNumElements());
      small.setFactorizationFrequency(100 + numberSort / 200);
      // Solve
      small.setLogLevel(CoinMax(0, logLevel() - 1));
      if (iPass > 20) {
        if (sumPrimalInfeasibilities_ > 1.0e-1) {
          small.dual();
        } else {
          small.primal(1);
        }
      } else {
        // presolve
#if 0
	ClpSolve::SolveType method;
	ClpSolve::PresolveType presolveType = ClpSolve::presolveOn;
	ClpSolve solveOptions;
	solveOptions.setPresolveType(presolveType, 5);
	if (sumPrimalInfeasibilities_>1.0e-1) 
	  method = ClpSolve::useDual;
	else
	  method = ClpSolve::usePrimalorSprint;
	solveOptions.setSolveType(method);
	small.initialSolve(solveOptions);
#else
#if 1
        ClpPresolve *pinfo = new ClpPresolve();
        ClpSimplex *small2 = pinfo->presolvedModel(small, 1.0e-5);
        assert(small2);
        if (sumPrimalInfeasibilities_ > 1.0e-1) {
          small2->dual();
        } else {
          small2->primal(1);
        }
        pinfo->postsolve(true);
        delete pinfo;
#else
        char *types = new char[small.numberRows() + small.numberColumns()];
        memset(types, 0, small.numberRows() + small.numberColumns());
        if (sumPrimalInfeasibilities_ > 1.0e-1) {
          small.miniSolve(types, types + small.numberRows(),
            -1, 0);
        } else {
          small.miniSolve(types, types + small.numberRows(),
            1, 1);
        }
        delete[] types;
#endif
        if (small.sumPrimalInfeasibilities() > 1.0)
          small.primal(1);
#endif
      }
      bool dualInfeasible = (small.status() == 2);
      // move solution back
      const double *smallSolution = small.primalColumnSolution();
      for (int j = 0; j < numberSmallColumns; j++) {
        int iColumn = whichColumns[j];
        columnActivity_[iColumn] = smallSolution[j];
        this->setColumnStatus(iColumn, small.getColumnStatus(j));
      }
      for (int iRow = 0; iRow < numberSort; iRow++) {
        int kRow = sort[iRow];
        this->setRowStatus(kRow, small.getRowStatus(iRow));
      }
      // compute full solution
      memset(rowActivity_, 0, numberRows_ * sizeof(double));
      times(1.0, columnActivity_, rowActivity_);
      if (iPass != maxPass - 1) {
        // Mark row as not looked at
        for (int iRow = 0; iRow < numberRows_; iRow++)
          weight[iRow] = 1.123e50;
        // Look at rows already in small problem
        int iSort;
        int numberDropped = 0;
        int numberKept = 0;
        int numberBinding = 0;
        numberPrimalInfeasibilities_ = 0;
        sumPrimalInfeasibilities_ = 0.0;
        bool allFeasible = small.numberIterations() == 0;
        for (iSort = 0; iSort < numberSort; iSort++) {
          int iRow = sort[iSort];
          //printf("%d %g %g\n",iRow,rowActivity_[iRow],small.primalRowSolution()[iSort]);
          if (getRowStatus(iRow) == ClpSimplex::basic) {
            // Basic - we can get rid of if early on
            if (iPass < takeOutPass && !dualInfeasible) {
              double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
                rowLower_[iRow] - rowActivity_[iRow]);
              weight[iRow] = -infeasibility;
              if (infeasibility > primalTolerance_ && !allFeasible) {
                numberPrimalInfeasibilities_++;
                sumPrimalInfeasibilities_ += infeasibility;
              } else {
                weight[iRow] = 1.0;
                numberDropped++;
              }
            } else {
              // keep
              weight[iRow] = -1.0e40;
              numberKept++;
            }
          } else {
            // keep
            weight[iRow] = -1.0e50;
            numberKept++;
            numberBinding++;
          }
        }
        // Now rest
        for (int iRow = 0; iRow < numberRows_; iRow++) {
          sort[iRow] = iRow;
          if (weight[iRow] == 1.123e50) {
            // not looked at yet
            double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
              rowLower_[iRow] - rowActivity_[iRow]);
            weight[iRow] = -infeasibility;
            if (infeasibility > primalTolerance_) {
              numberPrimalInfeasibilities_++;
              sumPrimalInfeasibilities_ += infeasibility;
            }
          }
        }
        // sort
        CoinSort_2(weight, weight + numberRows_, sort);
        numberSort = CoinMin(numberRows_, smallNumberRows + numberKept);
        memset(rowsIn, 0, numberRows_);
        for (int iRow = 0; iRow < numberSort; iRow++)
          rowsIn[sort[iRow]] = 1;
        printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
          numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns);
        printf("%d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
          sumPrimalInfeasibilities_);
        if (!numberPrimalInfeasibilities_) {
          problemStatus_ = 0;
          printf("Exiting as looks optimal\n");
          break;
        }
        numberPrimalInfeasibilities_ = 0;
        sumPrimalInfeasibilities_ = 0.0;
        for (iSort = 0; iSort < numberSort; iSort++) {
          if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) {
            numberPrimalInfeasibilities_++;
            sumPrimalInfeasibilities_ += -weight[iSort];
          }
        }
        printf("in small model %d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
          sumPrimalInfeasibilities_);
      } else {
        // out of passes
        problemStatus_ = -1;
      }
    }
    delete[] weight;
    delete[] sort;
    return 0;
  }
}
// A sequential LP method
int ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance,
  int otherOptions)
{
  // 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();
  // See if user wants to use cuts
  char *rowsIn = NULL;
  if ((otherOptions & 1) != 0 || numberPasses < 0) {
    numberPasses = abs(numberPasses);
    rowsIn = new char[numberRows_];
    memset(rowsIn, 0, numberRows_);
  }
  // get feasible
  if (this->status() < 0 || numberPrimalInfeasibilities())
    primalDualCuts(rowsIn, 1, 1);
  // still infeasible
  if (problemStatus_) {
    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];
  CoinMemcpyN(solution, numberColumns, saveSolutionM);
}
#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];
    }
    CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
  }
  // 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) {
        CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
        CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
      }
      // 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;
          }
          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) {
                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) {
                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;
        }
      }
    }
  }
  CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns,
    objective);
  //printf("offset comp %g orig %g\n",offset,objectiveOffset);
  // could do faster
  trueObjective->stepLength(this, solution, changeRegion, 0.0,
    objValue, predictedObj, thetaObj);
#ifdef CLP_INVESTIGATE
  printf("offset comp %g orig %g - obj from stepLength %g\n", offset, objectiveOffset, objValue);
#endif
  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) {
    CoinMemcpyN(solution, numberColumns, saveSolution);
    CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
    CoinMemcpyN(this->dualRowSolution(), numberRows, savePi);
    CoinMemcpyN(status_, numberRows + numberColumns, saveStatus);
#if MULTIPLE > 2
    double *tempSol = saveSolutionM[0];
    for (jNon = 0; jNon < MULTIPLE - 1; jNon++) {
      saveSolutionM[jNon] = saveSolutionM[jNon + 1];
    }
    saveSolutionM[MULTIPLE - 1] = tempSol;
    CoinMemcpyN(solution, numberColumns, tempSol);

#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);
    primalDualCuts(rowsIn, 1, 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
      CoinMemcpyN(safeSolution, numberColumns, solution);
      CoinMemcpyN(solution, numberColumns, saveSolution);
      memset(rowActivity_, 0, numberRows_ * sizeof(double));
      times(1.0, solution, rowActivity_);
      CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
      CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
      CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
      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
      CoinMemcpyN(solution, numberColumns, safeSolution);
    }
    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
    CoinMemcpyN(saveSolution, numberColumns, solution);
    CoinMemcpyN(saveRowSolution, numberRows, rowActivity_);
    CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
    CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
    iPass--;
    assert(exitPass > 0);
    goodMove = -1;
  }
}
#if MULTIPLE > 2
for (jNon = 0; jNon < MULTIPLE; jNon++) {
  delete[] saveSolutionM[jNon];
}
#endif
// restore solution
CoinMemcpyN(saveSolution, numberColumns, solution);
CoinMemcpyN(saveRowSolution, numberRows, rowActivity_);
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
CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns,
  objective);
primalDualCuts(rowsIn, 1, 1);
delete objective_;
delete[] rowsIn;
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;
    CoinBigIndex *addStarts = new CoinBigIndex[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) {
      CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns_,
        objective);
    } else {
      CoinMemcpyN(objective_->gradient(this, solution, offset, true, 2), numberColumns_,
        objective);
    }
    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;
          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;
          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
        CoinMemcpyN(saveSolution, numberColumns2, solution);
        CoinMemcpyN(saveSolution + numberColumns2,
          numberRows, newModel.primalRowSolution());
        CoinMemcpyN(reinterpret_cast< unsigned char * >(saveStatus),
          numberColumns2 + numberRows, newModel.statusArray());
        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];
          CoinMemcpyN(solution, numberColumns2, bestSolution);
        }
      }
      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;
      }
    }
    CoinMemcpyN(solution, numberColumns2, saveSolution);
    CoinMemcpyN(newModel.primalRowSolution(),
      numberRows, saveSolution + numberColumns2);
    CoinMemcpyN(newModel.statusArray(),
      numberColumns2 + numberRows,
      reinterpret_cast< unsigned char * >(saveStatus));

    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) {
    CoinMemcpyN(bestSolution, numberColumns2, solution);
    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();
  CoinMemcpyN(newModel.primalColumnSolution(), numberColumns_, columnActivity_);
  // should do status region
  CoinZeroN(rowActivity_, numberRows_);
  matrix_->times(1.0, columnActivity_, rowActivity_);
  return 0;
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/