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

/*
   Implements crude primal dual predictor corrector algorithm

 */
//#define SOME_DEBUG

#include "CoinPragma.hpp"
#include <math.h>

#include "CoinHelperFunctions.hpp"
#include "ClpPredictorCorrector.hpp"
#include "ClpEventHandler.hpp"
#include "CoinPackedMatrix.hpp"
#include "ClpMessage.hpp"
#include "ClpCholeskyBase.hpp"
#include "ClpHelperFunctions.hpp"
#include "ClpQuadraticObjective.hpp"
#include <cfloat>
#include <cassert>
#include <string>
#include <cstdio>
#include <iostream>
#if 0
static int yyyyyy = 0;
void ClpPredictorCorrector::saveSolution(std::string fileName)
{
     FILE * fp = fopen(fileName.c_str(), "wb");
     if (fp) {
          int numberRows = numberRows_;
          int numberColumns = numberColumns_;
          fwrite(&numberRows, sizeof(int), 1, fp);
          fwrite(&numberColumns, sizeof(int), 1, fp);
          CoinWorkDouble dsave[20];
          memset(dsave, 0, sizeof(dsave));
          fwrite(dsave, sizeof(CoinWorkDouble), 20, fp);
          int msave[20];
          memset(msave, 0, sizeof(msave));
          msave[0] = numberIterations_;
          fwrite(msave, sizeof(int), 20, fp);
          fwrite(dual_, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(errorRegion_, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(rhsFixRegion_, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(solution_, sizeof(CoinWorkDouble), numberColumns, fp);
          fwrite(solution_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(diagonal_, sizeof(CoinWorkDouble), numberColumns, fp);
          fwrite(diagonal_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(wVec_, sizeof(CoinWorkDouble), numberColumns, fp);
          fwrite(wVec_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(zVec_, sizeof(CoinWorkDouble), numberColumns, fp);
          fwrite(zVec_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(upperSlack_, sizeof(CoinWorkDouble), numberColumns, fp);
          fwrite(upperSlack_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
          fwrite(lowerSlack_, sizeof(CoinWorkDouble), numberColumns, fp);
          fwrite(lowerSlack_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
          fclose(fp);
     } else {
          std::cout << "Unable to open file " << fileName << std::endl;
     }
}
#endif
// Could change on CLP_LONG_CHOLESKY or COIN_LONG_WORK?
static CoinWorkDouble eScale = 1.0e27;
static CoinWorkDouble eBaseCaution = 1.0e-12;
static CoinWorkDouble eBase = 1.0e-12;
static CoinWorkDouble eDiagonal = 1.0e25;
static CoinWorkDouble eDiagonalCaution = 1.0e18;
static CoinWorkDouble eExtra = 1.0e-12;

// main function

int ClpPredictorCorrector::solve()
{
  problemStatus_ = -1;
  algorithm_ = 1;
  //create all regions
  if (!createWorkingData()) {
    problemStatus_ = 4;
    return 2;
  }
#if COIN_LONG_WORK
  // reallocate some regions
  double *dualSave = dual_;
  dual_ = reinterpret_cast< double * >(new CoinWorkDouble[numberRows_]);
  double *reducedCostSave = reducedCost_;
  reducedCost_ = reinterpret_cast< double * >(new CoinWorkDouble[numberColumns_]);
#endif
  //diagonalPerturbation_=1.0e-25;
  ClpMatrixBase *saveMatrix = NULL;
  // If quadratic then make copy so we can actually scale or normalize
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (objective_->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
#endif
  /* If modeSwitch is :
        0 - normal
        1 - bit switch off centering
        2 - bit always do type 2
        4 - accept corrector nearly always
     */
  int modeSwitch = 0;
  //if (quadraticObj)
  //modeSwitch |= 1; // switch off centring for now
  //if (quadraticObj)
  //modeSwitch |=4;
  ClpObjective *saveObjective = NULL;
  if (quadraticObj) {
    // check KKT is on
    if (!cholesky_->kkt()) {
      //No!
      handler_->message(CLP_BARRIER_KKT, messages_)
        << CoinMessageEol;
      return -1;
    }
    saveObjective = objective_;
    // We are going to make matrix full rather half
    objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
  }
  bool allowIncreasingGap = (modeSwitch & 4) != 0;
  // If scaled then really scale matrix
  if (scalingFlag_ > 0 && rowScale_) {
    saveMatrix = matrix_;
    matrix_ = matrix_->scaledColumnCopy(this);
  }
  //initializeFeasible(); - this just set fixed flag
  smallestInfeasibility_ = COIN_DBL_MAX;
  int i;
  for (i = 0; i < LENGTH_HISTORY; i++)
    historyInfeasibility_[i] = COIN_DBL_MAX;

  //bool firstTime=true;
  //firstFactorization(true);
  int returnCode = cholesky_->order(this);
  if (returnCode || cholesky_->symbolic()) {
    COIN_DETAIL_PRINT(printf("Error return from symbolic - probably not enough memory\n"));
    problemStatus_ = 4;
    //delete all temporary regions
    deleteWorkingData();
    if (saveMatrix) {
      // restore normal copy
      delete matrix_;
      matrix_ = saveMatrix;
    }
    // Restore quadratic objective if necessary
    if (saveObjective) {
      delete objective_;
      objective_ = saveObjective;
    }
    return -1;
  }
  mu_ = 1.0e10;
  diagonalScaleFactor_ = 1.0;
  //set iterations
  numberIterations_ = -1;
  int numberTotal = numberRows_ + numberColumns_;
  //initialize solution here
  if (createSolution() < 0) {
    COIN_DETAIL_PRINT(printf("Not enough memory\n"));
    problemStatus_ = 4;
    //delete all temporary regions
    deleteWorkingData();
    if (saveMatrix) {
      // restore normal copy
      delete matrix_;
      matrix_ = saveMatrix;
    }
    return -1;
  }
  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
  // Could try centering steps without any original step i.e. just center
  //firstFactorization(false);
  CoinZeroN(dualArray, numberRows_);
  multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
  matrix_->times(1.0, solution_, errorRegion_);
  maximumRHSError_ = maximumAbsElement(errorRegion_, numberRows_);
  maximumBoundInfeasibility_ = maximumRHSError_;
  //CoinWorkDouble maximumDualError_=COIN_DBL_MAX;
  //initialize
  actualDualStep_ = 0.0;
  actualPrimalStep_ = 0.0;
  gonePrimalFeasible_ = false;
  goneDualFeasible_ = false;
  //bool hadGoodSolution=false;
  diagonalNorm_ = solutionNorm_;
  mu_ = solutionNorm_;
  int numberFixed = updateSolution(-COIN_DBL_MAX);
  int numberFixedTotal = numberFixed;
  //int numberRows_DroppedBefore=0;
  //CoinWorkDouble extra=eExtra;
  //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX;
  //constants for infeas interior point
  const CoinWorkDouble beta2 = 0.99995;
  const CoinWorkDouble tau = 0.00002;
  CoinWorkDouble lastComplementarityGap = COIN_DBL_MAX * 1.0e-20;
  CoinWorkDouble lastStep = 1.0;
  // use to see if to take affine
  CoinWorkDouble checkGap = COIN_DBL_MAX;
  int lastGoodIteration = 0;
  CoinWorkDouble bestObjectiveGap = COIN_DBL_MAX;
  CoinWorkDouble bestObjective = COIN_DBL_MAX;
  int bestKilled = -1;
  int saveIteration = -1;
  int saveIteration2 = -1;
  bool sloppyOptimal = false;
  // this just to be used to exit
  bool sloppyOptimal2 = false;
  CoinWorkDouble *savePi = NULL;
  CoinWorkDouble *savePrimal = NULL;
  CoinWorkDouble *savePi2 = NULL;
  CoinWorkDouble *savePrimal2 = NULL;
  // Extra regions for centering
  CoinWorkDouble *saveX = new CoinWorkDouble[numberTotal];
  CoinWorkDouble *saveY = new CoinWorkDouble[numberRows_];
  CoinWorkDouble *saveZ = new CoinWorkDouble[numberTotal];
  CoinWorkDouble *saveW = new CoinWorkDouble[numberTotal];
  CoinWorkDouble *saveSL = new CoinWorkDouble[numberTotal];
  CoinWorkDouble *saveSU = new CoinWorkDouble[numberTotal];
  // Save smallest mu used in primal dual moves
  CoinWorkDouble objScale = optimizationDirection_ / (rhsScale_ * objectiveScale_);
  while (problemStatus_ < 0) {
    //#define FULL_DEBUG
#ifdef FULL_DEBUG
    {
      int i;
      printf("row    pi          artvec       rhsfx\n");
      for (i = 0; i < numberRows_; i++) {
        printf("%d %g %g %g\n", i, dual_[i], errorRegion_[i], rhsFixRegion_[i]);
      }
      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
      for (i = 0; i < numberColumns_ + numberRows_; i++) {
        printf(" %d %g %g %g %g %g %g\n", i, solution_[i], diagonal_[i], wVec_[i],
          zVec_[i], upperSlack_[i], lowerSlack_[i]);
      }
    }
#endif
    complementarityGap_ = complementarityGap(numberComplementarityPairs_,
      numberComplementarityItems_, 0);
    handler_->message(CLP_BARRIER_ITERATION, messages_)
      << numberIterations_
      << static_cast< double >(primalObjective_ * objScale - dblParam_[ClpObjOffset])
      << static_cast< double >(dualObjective_ * objScale - dblParam_[ClpObjOffset])
      << static_cast< double >(complementarityGap_)
      << numberFixedTotal
      << cholesky_->rank()
      << CoinMessageEol;
    // Check event
    {
      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
      if (status >= 0) {
        problemStatus_ = 5;
        secondaryStatus_ = ClpEventHandler::endOfIteration;
        break;
      }
    }
#if 0
          if (numberIterations_ == -1) {
               saveSolution("xxx.sav");
               if (yyyyyy)
                    exit(99);
          }
#endif
    // move up history
    for (i = 1; i < LENGTH_HISTORY; i++)
      historyInfeasibility_[i - 1] = historyInfeasibility_[i];
    historyInfeasibility_[LENGTH_HISTORY - 1] = complementarityGap_;
    // switch off saved if changes
    //if (saveIteration+10<numberIterations_&&
    //complementarityGap_*2.0<historyInfeasibility_[0])
    //saveIteration=-1;
    lastStep = CoinMin(actualPrimalStep_, actualDualStep_);
    CoinWorkDouble goodGapChange;
    //#define KEEP_GOING_IF_FIXED 5
#ifndef KEEP_GOING_IF_FIXED
#define KEEP_GOING_IF_FIXED 10000
#endif
    if (!sloppyOptimal) {
      goodGapChange = 0.93;
    } else {
      goodGapChange = 0.7;
      if (numberFixed > KEEP_GOING_IF_FIXED)
        goodGapChange = 0.99; // make more likely to carry on
    }
    CoinWorkDouble gapO;
    CoinWorkDouble lastGood = bestObjectiveGap;
    if (gonePrimalFeasible_ && goneDualFeasible_) {
      CoinWorkDouble largestObjective;
      if (CoinAbs(primalObjective_) > CoinAbs(dualObjective_)) {
        largestObjective = CoinAbs(primalObjective_);
      } else {
        largestObjective = CoinAbs(dualObjective_);
      }
      if (largestObjective < 1.0) {
        largestObjective = 1.0;
      }
      gapO = CoinAbs(primalObjective_ - dualObjective_) / largestObjective;
      handler_->message(CLP_BARRIER_OBJECTIVE_GAP, messages_)
        << static_cast< double >(gapO)
        << CoinMessageEol;
      //start saving best
      bool saveIt = false;
      if (gapO < bestObjectiveGap) {
        bestObjectiveGap = gapO;
#ifndef SAVE_ON_OBJ
        saveIt = true;
#endif
      }
      if (primalObjective_ < bestObjective) {
        bestObjective = primalObjective_;
#ifdef SAVE_ON_OBJ
        saveIt = true;
#endif
      }
      if (numberFixedTotal > bestKilled) {
        bestKilled = numberFixedTotal;
#if KEEP_GOING_IF_FIXED < 10
        saveIt = true;
#endif
      }
      if (saveIt) {
#if KEEP_GOING_IF_FIXED < 10
        COIN_DETAIL_PRINT(printf("saving\n"));
#endif
        saveIteration = numberIterations_;
        if (!savePi) {
          savePi = new CoinWorkDouble[numberRows_];
          savePrimal = new CoinWorkDouble[numberTotal];
        }
        CoinMemcpyN(dualArray, numberRows_, savePi);
        CoinMemcpyN(solution_, numberTotal, savePrimal);
      } else if (gapO > 2.0 * bestObjectiveGap) {
        //maybe be more sophisticated e.g. re-initialize having
        //fixed variables and dropped rows
        //std::cout <<" gap increasing "<<std::endl;
      }
      //std::cout <<"could stop"<<std::endl;
      //gapO=0.0;
      if (CoinAbs(primalObjective_ - dualObjective_) < dualTolerance()) {
        gapO = 0.0;
      }
    } else {
      gapO = COIN_DBL_MAX;
      if (saveIteration >= 0) {
        handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
          << CoinMessageEol;
        CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
        // save alternate
        if (numberFixedTotal > bestKilled && maximumBoundInfeasibility_ < 1.0e-6 && scaledRHSError < 1.0e-2) {
          bestKilled = numberFixedTotal;
#if KEEP_GOING_IF_FIXED < 10
          COIN_DETAIL_PRINT(printf("saving alternate\n"));
#endif
          saveIteration2 = numberIterations_;
          if (!savePi2) {
            savePi2 = new CoinWorkDouble[numberRows_];
            savePrimal2 = new CoinWorkDouble[numberTotal];
          }
          CoinMemcpyN(dualArray, numberRows_, savePi2);
          CoinMemcpyN(solution_, numberTotal, savePrimal2);
        }
        if (sloppyOptimal2) {
          // vaguely optimal
          if (maximumBoundInfeasibility_ > 1.0e-2 || scaledRHSError > 1.0e-2 || maximumDualError_ > objectiveNorm_ * 1.0e-2) {
            handler_->message(CLP_BARRIER_EXIT2, messages_)
              << saveIteration
              << CoinMessageEol;
            problemStatus_ = 0; // benefit of doubt
            break;
          }
        } else {
          // not close to optimal but check if getting bad
          CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
          if ((maximumBoundInfeasibility_ > 1.0e-1 || scaledRHSError > 1.0e-1 || maximumDualError_ > objectiveNorm_ * 1.0e-1)
            && (numberIterations_ > 50
                 && complementarityGap_ > 0.9 * historyInfeasibility_[0])) {
            handler_->message(CLP_BARRIER_EXIT2, messages_)
              << saveIteration
              << CoinMessageEol;
            break;
          }
          if (complementarityGap_ > 0.95 * checkGap && bestObjectiveGap < 1.0e-3 && (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
            handler_->message(CLP_BARRIER_EXIT2, messages_)
              << saveIteration
              << CoinMessageEol;
            break;
          }
        }
      }
      if (complementarityGap_ > 0.5 * checkGap && primalObjective_ > bestObjective + 1.0e-9 && (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
        handler_->message(CLP_BARRIER_EXIT2, messages_)
          << saveIteration
          << CoinMessageEol;
        break;
      }
    }
    // See if we should be thinking about exit if diverging
    double relativeMultiplier = 1.0 + fabs(primalObjective_) + fabs(dualObjective_);
    // Quadratic coding is rubbish so be more forgiving?
    if (quadraticObj)
      relativeMultiplier *= 5.0;
    if (gapO < 1.0e-5 + 1.0e-9 * relativeMultiplier
      || complementarityGap_ < 0.1 + 1.0e-9 * relativeMultiplier)
      sloppyOptimal2 = true;
    if ((gapO < 1.0e-6 || (gapO < 1.0e-4 && complementarityGap_ < 0.1)) && !sloppyOptimal) {
      sloppyOptimal = true;
      sloppyOptimal2 = true;
      handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL, messages_)
        << numberIterations_ << static_cast< double >(complementarityGap_)
        << CoinMessageEol;
    }
    int numberBack = quadraticObj ? 10 : 5;
    //tryJustPredictor=true;
    //printf("trying just predictor\n");
    //}
    if (complementarityGap_ >= 1.05 * lastComplementarityGap) {
      handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
        << static_cast< double >(complementarityGap_) << "increasing"
        << CoinMessageEol;
      if (saveIteration >= 0 && sloppyOptimal2) {
        handler_->message(CLP_BARRIER_EXIT2, messages_)
          << saveIteration
          << CoinMessageEol;
        break;
      } else if (numberIterations_ - lastGoodIteration >= numberBack && complementarityGap_ < 1.0e-6) {
        break; // not doing very well - give up
      }
    } else if (complementarityGap_ < goodGapChange * lastComplementarityGap) {
      lastGoodIteration = numberIterations_;
      lastComplementarityGap = complementarityGap_;
    } else if (numberIterations_ - lastGoodIteration >= numberBack && complementarityGap_ < 1.0e-3) {
      handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
        << static_cast< double >(complementarityGap_) << "not decreasing"
        << CoinMessageEol;
      if (gapO > 0.75 * lastGood && numberFixed < KEEP_GOING_IF_FIXED) {
        break;
      }
    } else if (numberIterations_ - lastGoodIteration >= 2 && complementarityGap_ < 1.0e-6) {
      handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
        << static_cast< double >(complementarityGap_) << "not decreasing"
        << CoinMessageEol;
      break;
    }
    if (numberIterations_ > maximumBarrierIterations_ || hitMaximumIterations()) {
      handler_->message(CLP_BARRIER_STOPPING, messages_)
        << CoinMessageEol;
      problemStatus_ = 3;
      onStopped(); // set secondary status
      break;
    }
    if (gapO < targetGap_) {
      problemStatus_ = 0;
      handler_->message(CLP_BARRIER_EXIT, messages_)
        << " "
        << CoinMessageEol;
      break; //finished
    }
    if (complementarityGap_ < 1.0e-12) {
      problemStatus_ = 0;
      handler_->message(CLP_BARRIER_EXIT, messages_)
        << "- small complementarity gap"
        << CoinMessageEol;
      break; //finished
    }
    if (complementarityGap_ < 1.0e-10 && gapO < 1.0e-10) {
      problemStatus_ = 0;
      handler_->message(CLP_BARRIER_EXIT, messages_)
        << "- objective gap and complementarity gap both small"
        << CoinMessageEol;
      break; //finished
    }
    if (gapO < 1.0e-9) {
      CoinWorkDouble value = gapO * complementarityGap_;
      value *= actualPrimalStep_;
      value *= actualDualStep_;
      //std::cout<<value<<std::endl;
      if (value < 1.0e-17 && numberIterations_ > lastGoodIteration) {
        problemStatus_ = 0;
        handler_->message(CLP_BARRIER_EXIT, messages_)
          << "- objective gap and complementarity gap both smallish and small steps"
          << CoinMessageEol;
        break; //finished
      }
    }
    CoinWorkDouble nextGap = COIN_DBL_MAX;
    int nextNumber = 0;
    int nextNumberItems = 0;
    worstDirectionAccuracy_ = 0.0;
    int newDropped = 0;
    //Predictor step
    //prepare for cholesky.  Set up scaled diagonal in deltaX
    //  ** for efficiency may be better if scale factor known before
    CoinWorkDouble norm2 = 0.0;
    CoinWorkDouble maximumValue;
    getNorms(diagonal_, numberTotal, maximumValue, norm2);
    diagonalNorm_ = CoinSqrt(norm2 / numberComplementarityPairs_);
    diagonalScaleFactor_ = 1.0;
    CoinWorkDouble maximumAllowable = eScale;
    //scale so largest is less than allowable ? could do better
    CoinWorkDouble factor = 0.5;
    while (maximumValue > maximumAllowable) {
      diagonalScaleFactor_ *= factor;
      maximumValue *= factor;
    } /* endwhile */
    if (diagonalScaleFactor_ != 1.0) {
      handler_->message(CLP_BARRIER_SCALING, messages_)
        << "diagonal" << static_cast< double >(diagonalScaleFactor_)
        << CoinMessageEol;
      diagonalNorm_ *= diagonalScaleFactor_;
    }
    multiplyAdd(NULL, numberTotal, 0.0, diagonal_,
      diagonalScaleFactor_);
    int *rowsDroppedThisTime = new int[numberRows_];
    newDropped = cholesky_->factorize(diagonal_, rowsDroppedThisTime);
    if (newDropped) {
      if (newDropped == -1) {
        COIN_DETAIL_PRINT(printf("Out of memory\n"));
        problemStatus_ = 4;
        //delete all temporary regions
        deleteWorkingData();
        if (saveMatrix) {
          // restore normal copy
          delete matrix_;
          matrix_ = saveMatrix;
        }
        return -1;
      } else {
#ifndef NDEBUG
        //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
        //assert(!newDropped2);
#endif
        if (newDropped < 0 && 0) {
          //replace dropped
          newDropped = -newDropped;
          //off 1 to allow for reset all
          newDropped--;
          //set all bits false
          cholesky_->resetRowsDropped();
        }
      }
    }
    delete[] rowsDroppedThisTime;
    if (cholesky_->status()) {
      std::cout << "bad cholesky?" << std::endl;
      abort();
    }
    int phase = 0; // predictor, corrector , primal dual
    CoinWorkDouble directionAccuracy = 0.0;
    bool doCorrector = true;
    bool goodMove = true;
    //set up for affine direction
    setupForSolve(phase);
    if ((modeSwitch & 2) == 0) {
      directionAccuracy = findDirectionVector(phase);
      if (directionAccuracy > worstDirectionAccuracy_) {
        worstDirectionAccuracy_ = directionAccuracy;
      }
      if (saveIteration > 0 && directionAccuracy > 1.0) {
        handler_->message(CLP_BARRIER_EXIT2, messages_)
          << saveIteration
          << CoinMessageEol;
        break;
      }
      findStepLength(phase);
      nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
      debugMove(0, actualPrimalStep_, actualDualStep_);
      debugMove(0, 1.0e-2, 1.0e-2);
    }
    CoinWorkDouble affineGap = nextGap;
    int bestPhase = 0;
    CoinWorkDouble bestNextGap = nextGap;
    // ?
    bestNextGap = CoinMax(nextGap, 0.8 * complementarityGap_);
    if (quadraticObj)
      bestNextGap = CoinMax(nextGap, 0.99 * complementarityGap_);
    if (complementarityGap_ > 1.0e-4 * numberComplementarityPairs_) {
      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
      CoinWorkDouble part1 = nextGap / numberComplementarityPairs_;
      part1 = nextGap / numberComplementarityItems_;
      CoinWorkDouble part2 = nextGap / complementarityGap_;
      mu_ = part1 * part2 * part2;
#if 0
               CoinWorkDouble papermu = complementarityGap_ / numberComplementarityPairs_;
               CoinWorkDouble affmu = nextGap / nextNumber;
               CoinWorkDouble sigma = pow(affmu / papermu, 3);
               printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
                      mu_, papermu, affmu, sigma, sigma * papermu);
#endif
      //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
      //					    (CoinWorkDouble) numberComplementarityPairs_));
    } else {
      CoinWorkDouble phi;
      if (numberComplementarityPairs_ <= 5000) {
        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 2.0);
      } else {
        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 1.5);
        if (phi < 500.0 * 500.0) {
          phi = 500.0 * 500.0;
        }
      }
      mu_ = complementarityGap_ / phi;
    }
    //save information
    CoinWorkDouble product = affineProduct();
#if 0
          //can we do corrector step?
          CoinWorkDouble xx = complementarityGap_ * (beta2 - tau) + product;
          if (xx > 0.0) {
               CoinWorkDouble saveMu = mu_;
               CoinWorkDouble mu2 = numberComplementarityPairs_;
               mu2 = xx / mu2;
               if (mu2 > mu_) {
                    //std::cout<<" could increase to "<<mu2<<std::endl;
                    //was mu2=mu2*0.25;
                    mu2 = mu2 * 0.99;
                    if (mu2 < mu_) {
                         mu_ = mu2;
                         //std::cout<<" changing mu to "<<mu_<<std::endl;
                    } else {
                         //std::cout<<std::endl;
                    }
               } else {
                    //std::cout<<" should decrease to "<<mu2<<std::endl;
                    mu_ = 0.5 * mu2;
                    //std::cout<<" changing mu to "<<mu_<<std::endl;
               }
               handler_->message(CLP_BARRIER_MU, messages_)
                         << saveMu << mu_
                         << CoinMessageEol;
          } else {
               //std::cout<<" bad by any standards"<<std::endl;
          }
#endif
    if (complementarityGap_ * (beta2 - tau) + product - mu_ * numberComplementarityPairs_ < 0.0 && 0) {
#ifdef SOME_DEBUG
      printf("failed 1 product %.18g mu %.18g - %.18g < 0.0, nextGap %.18g\n", product, mu_,
        complementarityGap_ * (beta2 - tau) + product - mu_ * numberComplementarityPairs_,
        nextGap);
#endif
      doCorrector = false;
      if (nextGap > 0.9 * complementarityGap_ || 1) {
        goodMove = false;
        bestNextGap = COIN_DBL_MAX;
      }
      //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_;
      //floatNumber = 1.0*numberComplementarityItems_;
      //mu_=nextGap/floatNumber;
      handler_->message(CLP_BARRIER_INFO, messages_)
        << "no corrector step"
        << CoinMessageEol;
    } else {
      phase = 1;
    }
    // If bad gap - try standard primal dual
    if (nextGap > complementarityGap_ * 1.001)
      goodMove = false;
    if ((modeSwitch & 2) != 0)
      goodMove = false;
    if (goodMove && doCorrector) {
      CoinMemcpyN(deltaX_, numberTotal, saveX);
      CoinMemcpyN(deltaY_, numberRows_, saveY);
      CoinMemcpyN(deltaZ_, numberTotal, saveZ);
      CoinMemcpyN(deltaW_, numberTotal, saveW);
      CoinMemcpyN(deltaSL_, numberTotal, saveSL);
      CoinMemcpyN(deltaSU_, numberTotal, saveSU);
#ifdef HALVE
      CoinWorkDouble savePrimalStep = actualPrimalStep_;
      CoinWorkDouble saveDualStep = actualDualStep_;
      CoinWorkDouble saveMu = mu_;
#endif
      //set up for next step
      setupForSolve(phase);
      CoinWorkDouble directionAccuracy2 = findDirectionVector(phase);
      if (directionAccuracy2 > worstDirectionAccuracy_) {
        worstDirectionAccuracy_ = directionAccuracy2;
      }
      CoinWorkDouble testValue = 1.0e2 * directionAccuracy;
      if (1.0e2 * projectionTolerance_ > testValue) {
        testValue = 1.0e2 * projectionTolerance_;
      }
      if (primalTolerance() > testValue) {
        testValue = primalTolerance();
      }
      if (maximumRHSError_ > testValue) {
        testValue = maximumRHSError_;
      }
      if (directionAccuracy2 > testValue && numberIterations_ >= -77) {
        goodMove = false;
#ifdef SOME_DEBUG
        printf("accuracy %g phase 1 failed, test value %g\n",
          directionAccuracy2, testValue);
#endif
      }
      if (goodMove) {
        phase = 1;
        CoinWorkDouble norm = findStepLength(phase);
        nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
        debugMove(1, actualPrimalStep_, actualDualStep_);
        //debugMove(1,1.0e-7,1.0e-7);
        goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
        if (norm < 0)
          goodMove = false;
        if (!goodMove) {
#ifdef SOME_DEBUG
          printf("checkGoodMove failed\n");
#endif
        }
      }
#ifdef HALVE
      int nHalve = 0;
      // relax test
      bestNextGap = CoinMax(bestNextGap, 0.9 * complementarityGap_);
      while (!goodMove) {
        mu_ = saveMu;
        actualPrimalStep_ = savePrimalStep;
        actualDualStep_ = saveDualStep;
        int i;
        //printf("halve %d\n",nHalve);
        nHalve++;
        const CoinWorkDouble lambda = 0.5;
        for (i = 0; i < numberRows_; i++)
          deltaY_[i] = lambda * deltaY_[i] + (1.0 - lambda) * saveY[i];
        for (i = 0; i < numberTotal; i++) {
          deltaX_[i] = lambda * deltaX_[i] + (1.0 - lambda) * saveX[i];
          deltaZ_[i] = lambda * deltaZ_[i] + (1.0 - lambda) * saveZ[i];
          deltaW_[i] = lambda * deltaW_[i] + (1.0 - lambda) * saveW[i];
          deltaSL_[i] = lambda * deltaSL_[i] + (1.0 - lambda) * saveSL[i];
          deltaSU_[i] = lambda * deltaSU_[i] + (1.0 - lambda) * saveSU[i];
        }
        CoinMemcpyN(saveX, numberTotal, deltaX_);
        CoinMemcpyN(saveY, numberRows_, deltaY_);
        CoinMemcpyN(saveZ, numberTotal, deltaZ_);
        CoinMemcpyN(saveW, numberTotal, deltaW_);
        CoinMemcpyN(saveSL, numberTotal, deltaSL_);
        CoinMemcpyN(saveSU, numberTotal, deltaSU_);
        findStepLength(1);
        nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
        goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
        if (nHalve > 10)
          break;
        //assert (goodMove);
      }
      if (nHalve && handler_->logLevel() > 2)
        printf("halved %d times\n", nHalve);
#endif
    }
    //bestPhase=-1;
    //goodMove=false;
    if (!goodMove) {
      // Just primal dual step
      CoinWorkDouble floatNumber;
      floatNumber = 2.0 * numberComplementarityPairs_;
      //floatNumber = numberComplementarityItems_;
      CoinWorkDouble saveMu = mu_; // use one from predictor corrector
      mu_ = complementarityGap_ / floatNumber;
      // If going well try small mu
      mu_ *= CoinSqrt((1.0 - lastStep) / (1.0 + 10.0 * lastStep));
      CoinWorkDouble mu1 = mu_;
      CoinWorkDouble phi;
      if (numberComplementarityPairs_ <= 500) {
        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 2.0);
      } else {
        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 1.5);
        if (phi < 500.0 * 500.0) {
          phi = 500.0 * 500.0;
        }
      }
      mu_ = complementarityGap_ / phi;
      //printf("pd mu %g, alternate %g, smallest\n",mu_,mu1);
      mu_ = CoinSqrt(mu_ * mu1);
      mu_ = mu1;
      if ((numberIterations_ & 1) == 0 || numberIterations_ < 10)
        mu_ = saveMu;
      // Try simpler
      floatNumber = numberComplementarityItems_;
      mu_ = 0.5 * complementarityGap_ / floatNumber;
      //if ((modeSwitch&2)==0) {
      //if ((numberIterations_&1)==0)
      //  mu_ *= 0.5;
      //} else {
      //mu_ *= 0.8;
      //}
      //set up for next step
      setupForSolve(2);
      findDirectionVector(2);
      CoinWorkDouble norm = findStepLength(2);
      // just for debug
      bestNextGap = complementarityGap_ * 1.0005;
      //bestNextGap=COIN_DBL_MAX;
      nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
      debugMove(2, actualPrimalStep_, actualDualStep_);
      //debugMove(2,1.0e-7,1.0e-7);
      checkGoodMove(false, bestNextGap, allowIncreasingGap);
      if ((nextGap > 0.9 * complementarityGap_ && bestPhase == 0 && affineGap < nextGap
            && (numberIterations_ > 80 || (numberIterations_ > 20 && quadraticObj)))
        || norm < 0.0) {
        // Back to affine
        phase = 0;
        setupForSolve(phase);
        directionAccuracy = findDirectionVector(phase);
        findStepLength(phase);
        nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
        bestNextGap = complementarityGap_;
        //checkGoodMove(false, bestNextGap,allowIncreasingGap);
      }
    }
    if (!goodMove)
      mu_ = nextGap / (static_cast< CoinWorkDouble >(nextNumber) * 1.1);
    //if (quadraticObj)
    //goodMove=true;
    //goodMove=false; //TEMP
    // Do centering steps
    int numberTries = 0;
    int numberGoodTries = 0;
#ifdef COIN_DETAIL
    CoinWorkDouble nextCenterGap = 0.0;
    CoinWorkDouble originalDualStep = actualDualStep_;
    CoinWorkDouble originalPrimalStep = actualPrimalStep_;
#endif
    if (actualDualStep_ > 0.9 && actualPrimalStep_ > 0.9)
      goodMove = false; // don't bother
    if ((modeSwitch & 1) != 0)
      goodMove = false;
    while (goodMove && numberTries < 5) {
      goodMove = false;
      numberTries++;
      CoinMemcpyN(deltaX_, numberTotal, saveX);
      CoinMemcpyN(deltaY_, numberRows_, saveY);
      CoinMemcpyN(deltaZ_, numberTotal, saveZ);
      CoinMemcpyN(deltaW_, numberTotal, saveW);
      CoinWorkDouble savePrimalStep = actualPrimalStep_;
      CoinWorkDouble saveDualStep = actualDualStep_;
      CoinWorkDouble saveMu = mu_;
      setupForSolve(3);
      findDirectionVector(3);
      findStepLength(3);
      debugMove(3, actualPrimalStep_, actualDualStep_);
      //debugMove(3,1.0e-7,1.0e-7);
      CoinWorkDouble xGap = complementarityGap(nextNumber, nextNumberItems, 3);
      // If one small then that's the one that counts
      CoinWorkDouble checkDual = saveDualStep;
      CoinWorkDouble checkPrimal = savePrimalStep;
      if (checkDual > 5.0 * checkPrimal) {
        checkDual = 2.0 * checkPrimal;
      } else if (checkPrimal > 5.0 * checkDual) {
        checkPrimal = 2.0 * checkDual;
      }
      if (actualPrimalStep_ < checkPrimal || actualDualStep_ < checkDual || (xGap > nextGap && xGap > 0.9 * complementarityGap_)) {
        //if (actualPrimalStep_<=checkPrimal||
        //actualDualStep_<=checkDual) {
#ifdef SOME_DEBUG
        printf("PP rejected gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
          actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
#endif
        mu_ = saveMu;
        actualPrimalStep_ = savePrimalStep;
        actualDualStep_ = saveDualStep;
        CoinMemcpyN(saveX, numberTotal, deltaX_);
        CoinMemcpyN(saveY, numberRows_, deltaY_);
        CoinMemcpyN(saveZ, numberTotal, deltaZ_);
        CoinMemcpyN(saveW, numberTotal, deltaW_);
      } else {
#ifdef SOME_DEBUG
        printf("PPphase 3 gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
          actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
#endif
        numberGoodTries++;
#ifdef COIN_DETAIL
        nextCenterGap = xGap;
#endif
        // See if big enough change
        if (actualPrimalStep_ < 1.01 * checkPrimal || actualDualStep_ < 1.01 * checkDual) {
          // stop now
        } else {
          // carry on
          goodMove = true;
        }
      }
    }
    if (numberGoodTries && handler_->logLevel() > 1) {
      COIN_DETAIL_PRINT(printf("%d centering steps moved from (gap %.18g, dual %.18g, primal %.18g) to (gap %.18g, dual %.18g, primal %.18g)\n",
        numberGoodTries, static_cast< double >(nextGap), static_cast< double >(originalDualStep),
        static_cast< double >(originalPrimalStep),
        static_cast< double >(nextCenterGap), static_cast< double >(actualDualStep_),
        static_cast< double >(actualPrimalStep_)));
    }
    // save last gap
    checkGap = complementarityGap_;
    numberFixed = updateSolution(nextGap);
    numberFixedTotal += numberFixed;
  } /* endwhile */
  delete[] saveX;
  delete[] saveY;
  delete[] saveZ;
  delete[] saveW;
  delete[] saveSL;
  delete[] saveSU;
  if (savePi) {
    if (numberIterations_ - saveIteration > 20 && numberIterations_ - saveIteration2 < 5) {
#if KEEP_GOING_IF_FIXED < 10
      std::cout << "Restoring2 from iteration " << saveIteration2 << std::endl;
#endif
      CoinMemcpyN(savePi2, numberRows_, dualArray);
      CoinMemcpyN(savePrimal2, numberTotal, solution_);
    } else {
#if KEEP_GOING_IF_FIXED < 10
      std::cout << "Restoring from iteration " << saveIteration << std::endl;
#endif
      CoinMemcpyN(savePi, numberRows_, dualArray);
      CoinMemcpyN(savePrimal, numberTotal, solution_);
    }
    delete[] savePi;
    delete[] savePrimal;
  }
  delete[] savePi2;
  delete[] savePrimal2;
  //recompute slacks
  // Split out solution
  CoinZeroN(rowActivity_, numberRows_);
  CoinMemcpyN(solution_, numberColumns_, columnActivity_);
  matrix_->times(1.0, columnActivity_, rowActivity_);
  //unscale objective
  multiplyAdd(NULL, numberTotal, 0.0, cost_, scaleFactor_);
  multiplyAdd(NULL, numberRows_, 0, dualArray, scaleFactor_);
  checkSolution();
  //CoinMemcpyN(reducedCost_,numberColumns_,dj_);
  // If quadratic use last solution
  // Restore quadratic objective if necessary
  if (saveObjective) {
    delete objective_;
    objective_ = saveObjective;
    objectiveValue_ = 0.5 * (primalObjective_ + dualObjective_);
  }
  handler_->message(CLP_BARRIER_END, messages_)
    << static_cast< double >(sumPrimalInfeasibilities_)
    << static_cast< double >(sumDualInfeasibilities_)
    << static_cast< double >(complementarityGap_)
    << static_cast< double >(objectiveValue())
    << CoinMessageEol;
  //#ifdef SOME_DEBUG
  if (handler_->logLevel() > 1)
    COIN_DETAIL_PRINT(printf("ENDRUN status %d after %d iterations\n", problemStatus_, numberIterations_));
    //#endif
    //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
    //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
    //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
    //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
    //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
#if COIN_LONG_WORK
  // put back dual
  delete[] dual_;
  delete[] reducedCost_;
  dual_ = dualSave;
  reducedCost_ = reducedCostSave;
#endif
  //delete all temporary regions
  deleteWorkingData();
#if KEEP_GOING_IF_FIXED < 10
#if 0 //ndef NDEBUG
     {
          static int kk = 0;
          char name[20];
          sprintf(name, "save.sol.%d", kk);
          kk++;
          printf("saving to file %s\n", name);
          FILE * fp = fopen(name, "wb");
          int numberWritten;
          numberWritten = fwrite(&numberColumns_, sizeof(int), 1, fp);
          assert (numberWritten == 1);
          numberWritten = fwrite(columnActivity_, sizeof(double), numberColumns_, fp);
          assert (numberWritten == numberColumns_);
          fclose(fp);
     }
#endif
#endif
  if (saveMatrix) {
    // restore normal copy
    delete matrix_;
    matrix_ = saveMatrix;
  }
  return problemStatus_;
}
// findStepLength.
//phase  - 0 predictor
//         1 corrector
//         2 primal dual
CoinWorkDouble ClpPredictorCorrector::findStepLength(int phase)
{
  CoinWorkDouble directionNorm = 0.0;
  CoinWorkDouble maximumPrimalStep = COIN_DBL_MAX * 1.0e-20;
  CoinWorkDouble maximumDualStep = COIN_DBL_MAX;
  int numberTotal = numberRows_ + numberColumns_;
  CoinWorkDouble tolerance = 1.0e-12;
#ifdef SOME_DEBUG
  int chosenPrimalSequence = -1;
  int chosenDualSequence = -1;
  bool lowPrimal = false;
  bool lowDual = false;
#endif
  // If done many iterations then allow to hit boundary
  CoinWorkDouble hitTolerance;
  //printf("objective norm %g\n",objectiveNorm_);
  if (numberIterations_ < 80 || !gonePrimalFeasible_)
    hitTolerance = COIN_DBL_MAX;
  else
    hitTolerance = CoinMax(1.0e3, 1.0e-3 * objectiveNorm_);
  int iColumn;
  //printf("dual value %g\n",dual_[0]);
  //printf("     X     dX      lS     dlS     uS     dUs    dj    Z dZ     t   dT\n");
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      CoinWorkDouble directionElement = deltaX_[iColumn];
      if (directionNorm < CoinAbs(directionElement)) {
        directionNorm = CoinAbs(directionElement);
      }
      if (lowerBound(iColumn)) {
        CoinWorkDouble delta = -deltaSL_[iColumn];
        CoinWorkDouble z1 = deltaZ_[iColumn];
        CoinWorkDouble newZ = zVec_[iColumn] + z1;
        if (zVec_[iColumn] > tolerance) {
          if (zVec_[iColumn] < -z1 * maximumDualStep) {
            maximumDualStep = -zVec_[iColumn] / z1;
#ifdef SOME_DEBUG
            chosenDualSequence = iColumn;
            lowDual = true;
#endif
          }
        }
        if (lowerSlack_[iColumn] < maximumPrimalStep * delta) {
          CoinWorkDouble newStep = lowerSlack_[iColumn] / delta;
          if (newStep > 0.2 || newZ < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] < hitTolerance) {
            maximumPrimalStep = newStep;
#ifdef SOME_DEBUG
            chosenPrimalSequence = iColumn;
            lowPrimal = true;
#endif
          } else {
            //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
          }
        }
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble delta = -deltaSU_[iColumn];
        ;
        CoinWorkDouble w1 = deltaW_[iColumn];
        CoinWorkDouble newT = wVec_[iColumn] + w1;
        if (wVec_[iColumn] > tolerance) {
          if (wVec_[iColumn] < -w1 * maximumDualStep) {
            maximumDualStep = -wVec_[iColumn] / w1;
#ifdef SOME_DEBUG
            chosenDualSequence = iColumn;
            lowDual = false;
#endif
          }
        }
        if (upperSlack_[iColumn] < maximumPrimalStep * delta) {
          CoinWorkDouble newStep = upperSlack_[iColumn] / delta;
          if (newStep > 0.2 || newT < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] > -hitTolerance) {
            maximumPrimalStep = newStep;
#ifdef SOME_DEBUG
            chosenPrimalSequence = iColumn;
            lowPrimal = false;
#endif
          } else {
            //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
          }
        }
      }
    }
  }
#ifdef SOME_DEBUG
  printf("new step - phase %d, norm %.18g, dual step %.18g, primal step %.18g\n",
    phase, directionNorm, maximumDualStep, maximumPrimalStep);
  if (lowDual)
    printf("ld %d %g %g => %g (dj %g,sol %g) ",
      chosenDualSequence, zVec_[chosenDualSequence],
      deltaZ_[chosenDualSequence], zVec_[chosenDualSequence] + maximumDualStep * deltaZ_[chosenDualSequence], dj_[chosenDualSequence],
      solution_[chosenDualSequence]);
  else
    printf("ud %d %g %g => %g (dj %g,sol %g) ",
      chosenDualSequence, wVec_[chosenDualSequence],
      deltaW_[chosenDualSequence], wVec_[chosenDualSequence] + maximumDualStep * deltaW_[chosenDualSequence], dj_[chosenDualSequence],
      solution_[chosenDualSequence]);
  if (lowPrimal)
    printf("lp %d %g %g => %g (dj %g,sol %g)\n",
      chosenPrimalSequence, lowerSlack_[chosenPrimalSequence],
      deltaSL_[chosenPrimalSequence], lowerSlack_[chosenPrimalSequence] + maximumPrimalStep * deltaSL_[chosenPrimalSequence],
      dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
  else
    printf("up %d %g %g => %g (dj %g,sol %g)\n",
      chosenPrimalSequence, upperSlack_[chosenPrimalSequence],
      deltaSU_[chosenPrimalSequence], upperSlack_[chosenPrimalSequence] + maximumPrimalStep * deltaSU_[chosenPrimalSequence],
      dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
#endif
  actualPrimalStep_ = stepLength_ * maximumPrimalStep;
  if (phase >= 0 && actualPrimalStep_ > 1.0) {
    actualPrimalStep_ = 1.0;
  }
  actualDualStep_ = stepLength_ * maximumDualStep;
  if (phase >= 0 && actualDualStep_ > 1.0) {
    actualDualStep_ = 1.0;
  }
  // See if quadratic objective
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (objective_->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
#endif
  if (quadraticObj) {
    // Use smaller unless very small
    CoinWorkDouble smallerStep = CoinMin(actualDualStep_, actualPrimalStep_);
    if (smallerStep > 0.0001) {
      actualDualStep_ = smallerStep;
      actualPrimalStep_ = smallerStep;
    }
  }
#define OFFQ
#ifndef OFFQ
  if (quadraticObj) {
    // Don't bother if phase 0 or 3 or large gap
    //if ((phase==1||phase==2||phase==0)&&maximumDualError_>0.1*complementarityGap_
    //&&smallerStep>0.001) {
    if ((phase == 1 || phase == 2 || phase == 0 || phase == 3)) {
      // minimize complementarity + norm*dual inf ? primal inf
      // at first - just check better - if not
      // Complementarity gap will be a*change*change + b*change +c
      CoinWorkDouble a = 0.0;
      CoinWorkDouble b = 0.0;
      CoinWorkDouble c = 0.0;
      /* SQUARE of dual infeasibility will be:
               square of dj - ......
               */
      CoinWorkDouble aq = 0.0;
      CoinWorkDouble bq = 0.0;
      CoinWorkDouble cq = 0.0;
      CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
      CoinWorkDouble *linearDjChange = new CoinWorkDouble[numberTotal];
      CoinZeroN(linearDjChange, numberColumns_);
      multiplyAdd(deltaY_, numberRows_, 1.0, linearDjChange + numberColumns_, 0.0);
      matrix_->transposeTimes(-1.0, deltaY_, linearDjChange);
      CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
      const int *columnQuadratic = quadratic->getIndices();
      const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
      const int *columnQuadraticLength = quadratic->getVectorLengths();
      CoinWorkDouble *quadraticElement = quadratic->getMutableElements();
      for (iColumn = 0; iColumn < numberTotal; iColumn++) {
        CoinWorkDouble oldPrimal = solution_[iColumn];
        if (!flagged(iColumn)) {
          if (lowerBound(iColumn)) {
            CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
            c += lowerSlack_[iColumn] * zVec_[iColumn];
            b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
            a += deltaZ_[iColumn] * change;
          }
          if (upperBound(iColumn)) {
            CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
            c += upperSlack_[iColumn] * wVec_[iColumn];
            b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
            a += deltaW_[iColumn] * change;
          }
          // new djs are dj_ + change*value
          CoinWorkDouble djChange = linearDjChange[iColumn];
          if (iColumn < numberColumns_) {
            for (CoinBigIndex j = columnQuadraticStart[iColumn];
                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
              int jColumn = columnQuadratic[j];
              CoinWorkDouble changeJ = deltaX_[jColumn];
              CoinWorkDouble elementValue = quadraticElement[j];
              djChange += changeJ * elementValue;
            }
          }
          CoinWorkDouble gammaTerm = gamma2;
          if (primalR_) {
            gammaTerm += primalR_[iColumn];
          }
          djChange += gammaTerm;
          // and dual infeasibility
          CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] + gammaTerm * solution_[iColumn];
          CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
          cq += oldInf * oldInf;
          bq += 2.0 * oldInf * changeInf;
          aq += changeInf * changeInf;
        } else {
          // fixed
          if (lowerBound(iColumn)) {
            c += lowerSlack_[iColumn] * zVec_[iColumn];
          }
          if (upperBound(iColumn)) {
            c += upperSlack_[iColumn] * wVec_[iColumn];
          }
          // new djs are dj_ + change*value
          CoinWorkDouble djChange = linearDjChange[iColumn];
          if (iColumn < numberColumns_) {
            for (CoinBigIndex j = columnQuadraticStart[iColumn];
                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
              int jColumn = columnQuadratic[j];
              CoinWorkDouble changeJ = deltaX_[jColumn];
              CoinWorkDouble elementValue = quadraticElement[j];
              djChange += changeJ * elementValue;
            }
          }
          CoinWorkDouble gammaTerm = gamma2;
          if (primalR_) {
            gammaTerm += primalR_[iColumn];
          }
          djChange += gammaTerm;
          // and dual infeasibility
          CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] + gammaTerm * solution_[iColumn];
          CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
          cq += oldInf * oldInf;
          bq += 2.0 * oldInf * changeInf;
          aq += changeInf * changeInf;
        }
      }
      delete[] linearDjChange;
      // ? We want to minimize complementarityGap + solutionNorm_*square of inf ??
      // maybe use inf and do line search
      // To check see if matches at current step
      CoinWorkDouble step = actualPrimalStep_;
      //Current gap + solutionNorm_ * CoinSqrt (sum square inf)
      CoinWorkDouble multiplier = solutionNorm_;
      multiplier *= 0.01;
      multiplier = 1.0;
      CoinWorkDouble currentInf = multiplier * CoinSqrt(cq);
      CoinWorkDouble nextInf = multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
      CoinWorkDouble allowedIncrease = 1.4;
#ifdef SOME_DEBUG
      printf("lin %g %g %g -> %g\n", a, b, c,
        c + b * step + a * step * step);
      printf("quad %g %g %g -> %g\n", aq, bq, cq,
        cq + bq * step + aq * step * step);
      debugMove(7, step, step);
      printf("current dualInf %g, with step of %g is %g\n",
        currentInf, step, nextInf);
#endif
      if (b > -1.0e-6) {
        if (phase != 0)
          directionNorm = -1.0;
      }
      if ((phase == 1 || phase == 2 || phase == 0 || phase == 3) && nextInf > 0.1 * complementarityGap_ && nextInf > currentInf * allowedIncrease) {
        //cq = CoinMax(cq,10.0);
        // convert to (x+q)*(x+q) = w
        CoinWorkDouble q = bq / (1.0 * aq);
        CoinWorkDouble w = CoinMax(q * q + (cq / aq) * (allowedIncrease - 1.0), 0.0);
        w = CoinSqrt(w);
        CoinWorkDouble stepX = w - q;
        step = stepX;
        nextInf = multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
#ifdef SOME_DEBUG
        printf("with step of %g dualInf is %g\n",
          step, nextInf);
#endif
        actualDualStep_ = CoinMin(step, actualDualStep_);
        actualPrimalStep_ = CoinMin(step, actualPrimalStep_);
      }
    }
  } else {
    // probably pointless as linear
    // minimize complementarity
    // Complementarity gap will be a*change*change + b*change +c
    CoinWorkDouble a = 0.0;
    CoinWorkDouble b = 0.0;
    CoinWorkDouble c = 0.0;
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      CoinWorkDouble oldPrimal = solution_[iColumn];
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
          c += lowerSlack_[iColumn] * zVec_[iColumn];
          b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
          a += deltaZ_[iColumn] * change;
        }
        if (upperBound(iColumn)) {
          CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
          c += upperSlack_[iColumn] * wVec_[iColumn];
          b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
          a += deltaW_[iColumn] * change;
        }
      } else {
        // fixed
        if (lowerBound(iColumn)) {
          c += lowerSlack_[iColumn] * zVec_[iColumn];
        }
        if (upperBound(iColumn)) {
          c += upperSlack_[iColumn] * wVec_[iColumn];
        }
      }
    }
    // ? We want to minimize complementarityGap;
    // maybe use inf and do line search
    // To check see if matches at current step
    CoinWorkDouble step = CoinMin(actualPrimalStep_, actualDualStep_);
    CoinWorkDouble next = c + b * step + a * step * step;
#ifdef SOME_DEBUG
    printf("lin %g %g %g -> %g\n", a, b, c,
      c + b * step + a * step * step);
    debugMove(7, step, step);
#endif
    if (b > -1.0e-6) {
      if (phase == 0) {
#ifdef SOME_DEBUG
        printf("*** odd phase 0 direction\n");
#endif
      } else {
        directionNorm = -1.0;
      }
    }
    // and with ratio
    a = 0.0;
    b = 0.0;
    CoinWorkDouble ratio = actualDualStep_ / actualPrimalStep_;
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      CoinWorkDouble oldPrimal = solution_[iColumn];
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
          b += lowerSlack_[iColumn] * deltaZ_[iColumn] * ratio + zVec_[iColumn] * change;
          a += deltaZ_[iColumn] * change * ratio;
        }
        if (upperBound(iColumn)) {
          CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
          b += upperSlack_[iColumn] * deltaW_[iColumn] * ratio + wVec_[iColumn] * change;
          a += deltaW_[iColumn] * change * ratio;
        }
      }
    }
    // ? We want to minimize complementarityGap;
    // maybe use inf and do line search
    // To check see if matches at current step
    step = actualPrimalStep_;
    CoinWorkDouble next2 = c + b * step + a * step * step;
    if (next2 > next) {
      actualPrimalStep_ = CoinMin(actualPrimalStep_, actualDualStep_);
      actualDualStep_ = actualPrimalStep_;
    }
#ifdef SOME_DEBUG
    printf("linb %g %g %g -> %g\n", a, b, c,
      c + b * step + a * step * step);
    debugMove(7, actualPrimalStep_, actualDualStep_);
#endif
    if (b > -1.0e-6) {
      if (phase == 0) {
#ifdef SOME_DEBUG
        printf("*** odd phase 0 direction\n");
#endif
      } else {
        directionNorm = -1.0;
      }
    }
  }
#else
  //actualPrimalStep_ =0.5*actualDualStep_;
#endif
#ifdef FULL_DEBUG
  if (phase == 3) {
    CoinWorkDouble minBeta = 0.1 * mu_;
    CoinWorkDouble maxBeta = 10.0 * mu_;
    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
            printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
              iColumn, primalValue, dualValue, gapProduct, delta2Z_[iColumn]);
        }
        if (upperBound(iColumn)) {
          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
              iColumn, primalValue, dualValue, gapProduct, delta2W_[iColumn]);
        }
      }
    }
  }
#endif
#ifdef SOME_DEBUG_not
  {
    CoinWorkDouble largestL = 0.0;
    CoinWorkDouble smallestL = COIN_DBL_MAX;
    CoinWorkDouble largestU = 0.0;
    CoinWorkDouble smallestU = COIN_DBL_MAX;
    CoinWorkDouble sumL = 0.0;
    CoinWorkDouble sumU = 0.0;
    int nL = 0;
    int nU = 0;
    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          largestL = CoinMax(largestL, gapProduct);
          smallestL = CoinMin(smallestL, gapProduct);
          nL++;
          sumL += gapProduct;
        }
        if (upperBound(iColumn)) {
          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          largestU = CoinMax(largestU, gapProduct);
          smallestU = CoinMin(smallestU, gapProduct);
          nU++;
          sumU += gapProduct;
        }
      }
    }
    CoinWorkDouble mu = (sumL + sumU) / (static_cast< CoinWorkDouble >(nL + nU));

    CoinWorkDouble minBeta = 0.1 * mu;
    CoinWorkDouble maxBeta = 10.0 * mu;
    int nBL = 0;
    int nAL = 0;
    int nBU = 0;
    int nAU = 0;
    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          if (gapProduct < minBeta)
            nBL++;
          else if (gapProduct > maxBeta)
            nAL++;
          //if (gapProduct<0.1*minBeta)
          //printf("Lsmall one %d dual %g primal %g\n",iColumn,
          //   dualValue,primalValue);
        }
        if (upperBound(iColumn)) {
          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          if (gapProduct < minBeta)
            nBU++;
          else if (gapProduct > maxBeta)
            nAU++;
          //if (gapProduct<0.1*minBeta)
          //printf("Usmall one %d dual %g primal %g\n",iColumn,
          //   dualValue,primalValue);
        }
      }
    }
    printf("phase %d new mu %.18g new gap %.18g\n", phase, mu, sumL + sumU);
    printf("          %d lower, smallest %.18g, %d below - largest %.18g, %d above\n",
      nL, smallestL, nBL, largestL, nAL);
    printf("          %d upper, smallest %.18g, %d below - largest %.18g, %d above\n",
      nU, smallestU, nBU, largestU, nAU);
  }
#endif
  return directionNorm;
}
/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
void ClpPredictorCorrector::solveSystem(CoinWorkDouble *region1, CoinWorkDouble *region2,
  const CoinWorkDouble *region1In, const CoinWorkDouble *region2In,
  const CoinWorkDouble *saveRegion1, const CoinWorkDouble *saveRegion2,
  bool gentleRefine)
{
  int iRow;
  int numberTotal = numberRows_ + numberColumns_;
  if (region2In) {
    // normal
    for (iRow = 0; iRow < numberRows_; iRow++)
      region2[iRow] = region2In[iRow];
  } else {
    // initial solution - (diagonal is 1 or 0)
    CoinZeroN(region2, numberRows_);
  }
  int iColumn;
  if (cholesky_->type() < 20) {
    // not KKT
    for (iColumn = 0; iColumn < numberTotal; iColumn++)
      region1[iColumn] = region1In[iColumn] * diagonal_[iColumn];
    multiplyAdd(region1 + numberColumns_, numberRows_, -1.0, region2, 1.0);
    matrix_->times(1.0, region1, region2);
    CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
    CoinWorkDouble scale = 1.0;
    CoinWorkDouble unscale = 1.0;
    if (maximumRHS > 1.0e-30) {
      if (maximumRHS <= 0.5) {
        CoinWorkDouble factor = 2.0;
        while (maximumRHS <= 0.5) {
          maximumRHS *= factor;
          scale *= factor;
        } /* endwhile */
      } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
        CoinWorkDouble factor = 0.5;
        while (maximumRHS >= 2.0) {
          maximumRHS *= factor;
          scale *= factor;
        } /* endwhile */
      }
      unscale = diagonalScaleFactor_ / scale;
    } else {
      //effectively zero
      scale = 0.0;
      unscale = 0.0;
    }
    multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
    cholesky_->solve(region2);
    multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
    multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns_, 0.0);
    CoinZeroN(region1, numberColumns_);
    matrix_->transposeTimes(1.0, region2, region1);
    for (iColumn = 0; iColumn < numberTotal; iColumn++)
      region1[iColumn] = (region1[iColumn] - region1In[iColumn]) * diagonal_[iColumn];
  } else {
    for (iColumn = 0; iColumn < numberTotal; iColumn++)
      region1[iColumn] = region1In[iColumn];
    cholesky_->solveKKT(region1, region2, diagonal_, diagonalScaleFactor_);
  }
  if (saveRegion2) {
    //refine?
    CoinWorkDouble scaleX = 1.0;
    if (gentleRefine)
      scaleX = 0.8;
    multiplyAdd(saveRegion2, numberRows_, 1.0, region2, scaleX);
    assert(saveRegion1);
    multiplyAdd(saveRegion1, numberTotal, 1.0, region1, scaleX);
  }
}
// findDirectionVector.
CoinWorkDouble ClpPredictorCorrector::findDirectionVector(const int phase)
{
  CoinWorkDouble projectionTolerance = projectionTolerance_;
  //temporary
  //projectionTolerance=1.0e-15;
  CoinWorkDouble errorCheck = 0.9 * maximumRHSError_ / solutionNorm_;
  if (errorCheck > primalTolerance()) {
    if (errorCheck < projectionTolerance) {
      projectionTolerance = errorCheck;
    }
  } else {
    if (primalTolerance() < projectionTolerance) {
      projectionTolerance = primalTolerance();
    }
  }
  CoinWorkDouble *newError = new CoinWorkDouble[numberRows_];
  int numberTotal = numberRows_ + numberColumns_;
  //if flagged then entries zero so can do
  // For KKT separate out
  CoinWorkDouble *region1Save = NULL; //for refinement
  int iColumn;
  if (cholesky_->type() < 20) {
    int iColumn;
    for (iColumn = 0; iColumn < numberTotal; iColumn++)
      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
    multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
    matrix_->times(1.0, deltaX_, deltaY_);
  } else {
    // regions in will be workArray and newError
    // regions out deltaX_ and deltaY_
    multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, newError, 0.0);
    matrix_->times(-1.0, solution_, newError);
    // This is inefficient but just for now get values which will be in deltay
    int iColumn;
    for (iColumn = 0; iColumn < numberTotal; iColumn++)
      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
    multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
    matrix_->times(1.0, deltaX_, deltaY_);
  }
  bool goodSolve = false;
  CoinWorkDouble *regionSave = NULL; //for refinement
  int numberTries = 0;
  CoinWorkDouble relativeError = COIN_DBL_MAX;
  CoinWorkDouble tryError = 1.0e31;
  CoinWorkDouble saveMaximum = 0.0;
  double firstError = 0.0;
  double lastError2 = 0.0;
  while (!goodSolve && numberTries < 30) {
    CoinWorkDouble lastError = relativeError;
    goodSolve = true;
    CoinWorkDouble maximumRHS;
    maximumRHS = CoinMax(maximumAbsElement(deltaY_, numberRows_), 1.0e-12);
    if (!numberTries)
      saveMaximum = maximumRHS;
    if (cholesky_->type() < 20) {
      // no kkt
      CoinWorkDouble scale = 1.0;
      CoinWorkDouble unscale = 1.0;
      if (maximumRHS > 1.0e-30) {
        if (maximumRHS <= 0.5) {
          CoinWorkDouble factor = 2.0;
          while (maximumRHS <= 0.5) {
            maximumRHS *= factor;
            scale *= factor;
          } /* endwhile */
        } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
          CoinWorkDouble factor = 0.5;
          while (maximumRHS >= 2.0) {
            maximumRHS *= factor;
            scale *= factor;
          } /* endwhile */
        }
        unscale = diagonalScaleFactor_ / scale;
      } else {
        //effectively zero
        scale = 0.0;
        unscale = 0.0;
      }
      //printf("--putting scales to 1.0\n");
      //scale=1.0;
      //unscale=1.0;
      multiplyAdd(NULL, numberRows_, 0.0, deltaY_, scale);
      cholesky_->solve(deltaY_);
      multiplyAdd(NULL, numberRows_, 0.0, deltaY_, unscale);
#if 0
               {
                    printf("deltay\n");
                    for (int i = 0; i < numberRows_; i++)
                         printf("%d %.18g\n", i, deltaY_[i]);
               }
               exit(66);
#endif
      if (numberTries) {
        //refine?
        CoinWorkDouble scaleX = 1.0;
        if (lastError > 1.0e-5)
          scaleX = 0.8;
        multiplyAdd(regionSave, numberRows_, 1.0, deltaY_, scaleX);
      }
      //CoinZeroN(newError,numberRows_);
      multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
      CoinZeroN(deltaX_, numberColumns_);
      matrix_->transposeTimes(1.0, deltaY_, deltaX_);
      //if flagged then entries zero so can do
      for (iColumn = 0; iColumn < numberTotal; iColumn++)
        deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
          - workArray_[iColumn];
    } else {
      // KKT
      solveSystem(deltaX_, deltaY_,
        workArray_, newError, region1Save, regionSave, lastError > 1.0e-5);
    }
    multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, newError, 0.0);
    matrix_->times(1.0, deltaX_, newError);
    numberTries++;

    //now add in old Ax - doing extra checking
    CoinWorkDouble maximumRHSError = 0.0;
    CoinWorkDouble maximumRHSChange = 0.0;
    int iRow;
    char *dropped = cholesky_->rowsDropped();
    for (iRow = 0; iRow < numberRows_; iRow++) {
      if (!dropped[iRow]) {
        CoinWorkDouble newValue = newError[iRow];
        CoinWorkDouble oldValue = errorRegion_[iRow];
        //severity of errors depend on signs
        //**later                                                             */
        if (CoinAbs(newValue) > maximumRHSChange) {
          maximumRHSChange = CoinAbs(newValue);
        }
        CoinWorkDouble result = newValue + oldValue;
        if (CoinAbs(result) > maximumRHSError) {
          maximumRHSError = CoinAbs(result);
        }
        newError[iRow] = result;
      } else {
        CoinWorkDouble newValue = newError[iRow];
        CoinWorkDouble oldValue = errorRegion_[iRow];
        if (CoinAbs(newValue) > maximumRHSChange) {
          maximumRHSChange = CoinAbs(newValue);
        }
        CoinWorkDouble result = newValue + oldValue;
        newError[iRow] = result;
        //newError[iRow]=0.0;
        //assert(deltaY_[iRow]==0.0);
        deltaY_[iRow] = 0.0;
      }
    }
    relativeError = maximumRHSError / solutionNorm_;
    relativeError = maximumRHSError / saveMaximum;
    if (relativeError > tryError)
      relativeError = tryError;
    if (numberTries == 1)
      firstError = relativeError;
    if (relativeError < lastError) {
      lastError2 = relativeError;
      maximumRHSChange_ = maximumRHSChange;
      if (relativeError > projectionTolerance && numberTries <= 3) {
        //try and refine
        goodSolve = false;
      }
      //*** extra test here
      if (!goodSolve) {
        if (!regionSave) {
          regionSave = new CoinWorkDouble[numberRows_];
          if (cholesky_->type() >= 20)
            region1Save = new CoinWorkDouble[numberTotal];
        }
        CoinMemcpyN(deltaY_, numberRows_, regionSave);
        if (cholesky_->type() < 20) {
          // not KKT
          multiplyAdd(newError, numberRows_, -1.0, deltaY_, 0.0);
        } else {
          // KKT
          CoinMemcpyN(deltaX_, numberTotal, region1Save);
          // and back to input region
          CoinMemcpyN(deltaY_, numberRows_, newError);
        }
      }
    } else {
      //std::cout <<" worse residual = "<<relativeError;
      //bring back previous
      relativeError = lastError;
      if (regionSave) {
        CoinMemcpyN(regionSave, numberRows_, deltaY_);
        if (cholesky_->type() < 20) {
          // not KKT
          multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
          CoinZeroN(deltaX_, numberColumns_);
          matrix_->transposeTimes(1.0, deltaY_, deltaX_);
          //if flagged then entries zero so can do
          for (iColumn = 0; iColumn < numberTotal; iColumn++)
            deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
              - workArray_[iColumn];
        } else {
          // KKT
          CoinMemcpyN(region1Save, numberTotal, deltaX_);
        }
      } else {
        // disaster
        CoinFillN(deltaX_, numberTotal, static_cast< CoinWorkDouble >(1.0));
        CoinFillN(deltaY_, numberRows_, static_cast< CoinWorkDouble >(1.0));
        COIN_DETAIL_PRINT(printf("bad cholesky\n"));
      }
    }
  } /* endwhile */
  if (firstError > 1.0e-8 || numberTries > 1) {
    handler_->message(CLP_BARRIER_ACCURACY, messages_)
      << phase << numberTries << static_cast< double >(firstError)
      << static_cast< double >(lastError2)
      << CoinMessageEol;
  }
  delete[] regionSave;
  delete[] region1Save;
  delete[] newError;
  // now rest
  CoinWorkDouble extra = eExtra;
  //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
  //CoinZeroN(deltaW_,numberColumns_);
  //matrix_->transposeTimes(-1.0,deltaY_,deltaW_);

  for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    deltaSU_[iColumn] = 0.0;
    deltaSL_[iColumn] = 0.0;
    deltaZ_[iColumn] = 0.0;
    CoinWorkDouble dd = deltaW_[iColumn];
    deltaW_[iColumn] = 0.0;
    if (!flagged(iColumn)) {
      CoinWorkDouble deltaX = deltaX_[iColumn];
      if (lowerBound(iColumn)) {
        CoinWorkDouble zValue = rhsZ_[iColumn];
        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
        CoinWorkDouble slack = lowerSlack_[iColumn] + extra;
        deltaSL_[iColumn] = -rhsL_[iColumn] + deltaX;
        deltaZ_[iColumn] = (gHat - zVec_[iColumn] * deltaX) / slack;
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble wValue = rhsW_[iColumn];
        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
        CoinWorkDouble slack = upperSlack_[iColumn] + extra;
        deltaSU_[iColumn] = rhsU_[iColumn] - deltaX;
        deltaW_[iColumn] = (hHat + wVec_[iColumn] * deltaX) / slack;
      }
      if (0) {
        // different way of calculating
        CoinWorkDouble gamma2 = gamma_ * gamma_;
        CoinWorkDouble dZ = 0.0;
        CoinWorkDouble dW = 0.0;
        CoinWorkDouble zValue = rhsZ_[iColumn];
        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
        CoinWorkDouble slackL = lowerSlack_[iColumn] + extra;
        CoinWorkDouble wValue = rhsW_[iColumn];
        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
        CoinWorkDouble slackU = upperSlack_[iColumn] + extra;
        CoinWorkDouble q = rhsC_[iColumn] + gamma2 * deltaX + dd;
        if (primalR_)
          q += deltaX * primalR_[iColumn];
        dW = (gHat + hHat - slackL * q + (wValue - zValue) * deltaX) / (slackL + slackU);
        dZ = dW + q;
        //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
        //deltaW_[iColumn],dZ,dW);
        if (lowerBound(iColumn)) {
          if (upperBound(iColumn)) {
            //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
            //deltaW_[iColumn],dZ,dW);
            deltaW_[iColumn] = dW;
            deltaZ_[iColumn] = dZ;
          } else {
            // just lower
            //printf("L %d old %g new %g\n",iColumn,deltaZ_[iColumn],
            //dZ);
          }
        } else {
          assert(upperBound(iColumn));
          //printf("U %d old %g new %g\n",iColumn,deltaW_[iColumn],
          //dW);
        }
      }
    }
  }
#if 0
     CoinWorkDouble * check = new CoinWorkDouble[numberTotal];
     // Check out rhsC_
     multiplyAdd(deltaY_, numberRows_, -1.0, check + numberColumns_, 0.0);
     CoinZeroN(check, numberColumns_);
     matrix_->transposeTimes(1.0, deltaY_, check);
     quadraticDjs(check, deltaX_, -1.0);
     for (iColumn = 0; iColumn < numberTotal; iColumn++) {
          check[iColumn] += deltaZ_[iColumn] - deltaW_[iColumn];
          if (CoinAbs(check[iColumn] - rhsC_[iColumn]) > 1.0e-3)
               printf("rhsC %d %g %g\n", iColumn, check[iColumn], rhsC_[iColumn]);
     }
     // Check out rhsZ_
     for (iColumn = 0; iColumn < numberTotal; iColumn++) {
          check[iColumn] += lowerSlack_[iColumn] * deltaZ_[iColumn] +
                            zVec_[iColumn] * deltaSL_[iColumn];
          if (CoinAbs(check[iColumn] - rhsZ_[iColumn]) > 1.0e-3)
               printf("rhsZ %d %g %g\n", iColumn, check[iColumn], rhsZ_[iColumn]);
     }
     // Check out rhsW_
     for (iColumn = 0; iColumn < numberTotal; iColumn++) {
          check[iColumn] += upperSlack_[iColumn] * deltaW_[iColumn] +
                            wVec_[iColumn] * deltaSU_[iColumn];
          if (CoinAbs(check[iColumn] - rhsW_[iColumn]) > 1.0e-3)
               printf("rhsW %d %g %g\n", iColumn, check[iColumn], rhsW_[iColumn]);
     }
     delete [] check;
#endif
  return relativeError;
}
// createSolution.  Creates solution from scratch
int ClpPredictorCorrector::createSolution()
{
  int numberTotal = numberRows_ + numberColumns_;
  int iColumn;
  CoinWorkDouble tolerance = primalTolerance();
  // See if quadratic objective
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (objective_->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
#endif
  if (!quadraticObj) {
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      if (upper_[iColumn] - lower_[iColumn] > tolerance)
        clearFixed(iColumn);
      else
        setFixed(iColumn);
    }
  } else {
    // try leaving fixed
    for (iColumn = 0; iColumn < numberTotal; iColumn++)
      clearFixed(iColumn);
  }

  CoinWorkDouble maximumObjective = 0.0;
  CoinWorkDouble objectiveNorm2 = 0.0;
  getNorms(cost_, numberTotal, maximumObjective, objectiveNorm2);
  if (!maximumObjective) {
    maximumObjective = 1.0; // objective all zero
  }
  objectiveNorm2 = CoinSqrt(objectiveNorm2) / static_cast< CoinWorkDouble >(numberTotal);
  objectiveNorm_ = maximumObjective;
  scaleFactor_ = 1.0;
  if (maximumObjective > 0.0) {
    if (maximumObjective < 1.0) {
      scaleFactor_ = maximumObjective;
    } else if (maximumObjective > 1.0e4) {
      scaleFactor_ = maximumObjective / 1.0e4;
    }
  }
  if (scaleFactor_ != 1.0) {
    objectiveNorm2 *= scaleFactor_;
    multiplyAdd(NULL, numberTotal, 0.0, cost_, 1.0 / scaleFactor_);
    objectiveNorm_ = maximumObjective / scaleFactor_;
  }
  // See if quadratic objective
  if (quadraticObj) {
    // If scaled then really scale matrix
    CoinWorkDouble scaleFactor = scaleFactor_ * optimizationDirection_ * objectiveScale_ * rhsScale_;
    if ((scalingFlag_ > 0 && rowScale_) || scaleFactor != 1.0) {
      CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
      const int *columnQuadratic = quadratic->getIndices();
      const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
      const int *columnQuadraticLength = quadratic->getVectorLengths();
      double *quadraticElement = quadratic->getMutableElements();
      int numberColumns = quadratic->getNumCols();
      CoinWorkDouble scale = 1.0 / scaleFactor;
      if (scalingFlag_ > 0 && rowScale_) {
        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
          CoinWorkDouble scaleI = columnScale_[iColumn] * scale;
          for (CoinBigIndex j = columnQuadraticStart[iColumn];
               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
            int jColumn = columnQuadratic[j];
            CoinWorkDouble scaleJ = columnScale_[jColumn];
            quadraticElement[j] *= scaleI * scaleJ;
            objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
          }
        }
      } else {
        // not scaled
        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
          for (CoinBigIndex j = columnQuadraticStart[iColumn];
               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
            quadraticElement[j] *= scale;
            objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
          }
        }
      }
    }
  }
  baseObjectiveNorm_ = objectiveNorm_;
  //accumulate fixed in dj region (as spare)
  //accumulate primal solution in primal region
  //DZ in lowerDual
  //DW in upperDual
  CoinWorkDouble infiniteCheck = 1.0e40;
  //CoinWorkDouble     fakeCheck=1.0e10;
  //use deltaX region for work region
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    CoinWorkDouble primalValue = solution_[iColumn];
    clearFlagged(iColumn);
    clearFixedOrFree(iColumn);
    clearLowerBound(iColumn);
    clearUpperBound(iColumn);
    clearFakeLower(iColumn);
    clearFakeUpper(iColumn);
    if (!fixed(iColumn)) {
      dj_[iColumn] = 0.0;
      diagonal_[iColumn] = 1.0;
      deltaX_[iColumn] = 1.0;
      CoinWorkDouble lowerValue = lower_[iColumn];
      CoinWorkDouble upperValue = upper_[iColumn];
      if (lowerValue > -infiniteCheck) {
        if (upperValue < infiniteCheck) {
          //upper and lower bounds
          setLowerBound(iColumn);
          setUpperBound(iColumn);
          if (lowerValue >= 0.0) {
            solution_[iColumn] = lowerValue;
          } else if (upperValue <= 0.0) {
            solution_[iColumn] = upperValue;
          } else {
            solution_[iColumn] = 0.0;
          }
        } else {
          //just lower bound
          setLowerBound(iColumn);
          if (lowerValue >= 0.0) {
            solution_[iColumn] = lowerValue;
          } else {
            solution_[iColumn] = 0.0;
          }
        }
      } else {
        if (upperValue < infiniteCheck) {
          //just upper bound
          setUpperBound(iColumn);
          if (upperValue <= 0.0) {
            solution_[iColumn] = upperValue;
          } else {
            solution_[iColumn] = 0.0;
          }
        } else {
          //free
          setFixedOrFree(iColumn);
          solution_[iColumn] = 0.0;
          //std::cout<<" free "<<i<<std::endl;
        }
      }
    } else {
      setFlagged(iColumn);
      setFixedOrFree(iColumn);
      setLowerBound(iColumn);
      setUpperBound(iColumn);
      dj_[iColumn] = primalValue;
      ;
      solution_[iColumn] = lower_[iColumn];
      diagonal_[iColumn] = 0.0;
      deltaX_[iColumn] = 0.0;
    }
  }
  //   modify fixed RHS
  multiplyAdd(dj_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 0.0);
  //   create plausible RHS?
  matrix_->times(-1.0, dj_, rhsFixRegion_);
  multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, errorRegion_, 0.0);
  matrix_->times(-1.0, solution_, errorRegion_);
  rhsNorm_ = maximumAbsElement(errorRegion_, numberRows_);
  if (rhsNorm_ < 1.0) {
    rhsNorm_ = 1.0;
  }
  int *rowsDropped = new int[numberRows_];
  int returnCode = cholesky_->factorize(diagonal_, rowsDropped);
  if (returnCode == -1) {
    COIN_DETAIL_PRINT(printf("Out of memory\n"));
    problemStatus_ = 4;
    return -1;
  }
  if (cholesky_->status()) {
    std::cout << "singular on initial cholesky?" << std::endl;
    cholesky_->resetRowsDropped();
    //cholesky_->factorize(rowDropped_);
    //if (cholesky_->status()) {
    //std::cout << "bad cholesky??? (after retry)" <<std::endl;
    //abort();
    //}
  }
  delete[] rowsDropped;
  if (cholesky_->type() < 20) {
    // not KKT
    cholesky_->solve(errorRegion_);
    //create information for solution
    multiplyAdd(errorRegion_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
    CoinZeroN(deltaX_, numberColumns_);
    matrix_->transposeTimes(1.0, errorRegion_, deltaX_);
  } else {
    // KKT
    // reverse sign on solution
    multiplyAdd(NULL, numberRows_ + numberColumns_, 0.0, solution_, -1.0);
    solveSystem(deltaX_, errorRegion_, solution_, NULL, NULL, NULL, false);
  }
  CoinWorkDouble initialValue = 1.0e2;
  if (rhsNorm_ * 1.0e-2 > initialValue) {
    initialValue = rhsNorm_ * 1.0e-2;
  }
  //initialValue = CoinMax(1.0,rhsNorm_);
  CoinWorkDouble smallestBoundDifference = COIN_DBL_MAX;
  CoinWorkDouble *fakeSolution = deltaX_;
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      if (lower_[iColumn] - fakeSolution[iColumn] > initialValue) {
        initialValue = lower_[iColumn] - fakeSolution[iColumn];
      }
      if (fakeSolution[iColumn] - upper_[iColumn] > initialValue) {
        initialValue = fakeSolution[iColumn] - upper_[iColumn];
      }
      if (upper_[iColumn] - lower_[iColumn] < smallestBoundDifference) {
        smallestBoundDifference = upper_[iColumn] - lower_[iColumn];
      }
    }
  }
  solutionNorm_ = 1.0e-12;
  handler_->message(CLP_BARRIER_SAFE, messages_)
    << static_cast< double >(initialValue) << static_cast< double >(objectiveNorm_)
    << CoinMessageEol;
  CoinWorkDouble extra = 1.0e-10;
  CoinWorkDouble largeGap = 1.0e15;
  //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_;
  CoinWorkDouble safeObjectiveValue = objectiveNorm_ + 1.0;
  CoinWorkDouble safeFree = 1.0e-1 * initialValue;
  //printf("normal safe dual value of %g, primal value of %g\n",
  // safeObjectiveValue,initialValue);
  //safeObjectiveValue=CoinMax(2.0,1.0e-1*safeObjectiveValue);
  //initialValue=CoinMax(100.0,1.0e-1*initialValue);;
  //printf("temp safe dual value of %g, primal value of %g\n",
  // safeObjectiveValue,initialValue);
  CoinWorkDouble zwLarge = 1.0e2 * initialValue;
  //zwLarge=1.0e40;
  if (cholesky_->choleskyCondition() < 0.0 && cholesky_->type() < 20) {
    // looks bad - play safe
    initialValue *= 10.0;
    safeObjectiveValue *= 10.0;
    safeFree *= 10.0;
  }
  CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
  // First do primal side
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      CoinWorkDouble lowerValue = lower_[iColumn];
      CoinWorkDouble upperValue = upper_[iColumn];
      CoinWorkDouble newValue;
      CoinWorkDouble setPrimal = initialValue;
      if (quadraticObj) {
        // perturb primal solution a bit
        //fakeSolution[iColumn]  *= 0.002*CoinDrand48()+0.999;
      }
      if (lowerBound(iColumn)) {
        if (upperBound(iColumn)) {
          //upper and lower bounds
          if (upperValue - lowerValue > 2.0 * setPrimal) {
            CoinWorkDouble fakeValue = fakeSolution[iColumn];
            if (fakeValue < lowerValue + setPrimal) {
              fakeValue = lowerValue + setPrimal;
            }
            if (fakeValue > upperValue - setPrimal) {
              fakeValue = upperValue - setPrimal;
            }
            newValue = fakeValue;
          } else {
            newValue = 0.5 * (upperValue + lowerValue);
          }
        } else {
          //just lower bound
          CoinWorkDouble fakeValue = fakeSolution[iColumn];
          if (fakeValue < lowerValue + setPrimal) {
            fakeValue = lowerValue + setPrimal;
          }
          newValue = fakeValue;
        }
      } else {
        if (upperBound(iColumn)) {
          //just upper bound
          CoinWorkDouble fakeValue = fakeSolution[iColumn];
          if (fakeValue > upperValue - setPrimal) {
            fakeValue = upperValue - setPrimal;
          }
          newValue = fakeValue;
        } else {
          //free
          newValue = fakeSolution[iColumn];
          if (newValue >= 0.0) {
            if (newValue < safeFree) {
              newValue = safeFree;
            }
          } else {
            if (newValue > -safeFree) {
              newValue = -safeFree;
            }
          }
        }
      }
      solution_[iColumn] = newValue;
    } else {
      // fixed
      lowerSlack_[iColumn] = 0.0;
      upperSlack_[iColumn] = 0.0;
      solution_[iColumn] = lower_[iColumn];
      zVec_[iColumn] = 0.0;
      wVec_[iColumn] = 0.0;
      diagonal_[iColumn] = 0.0;
    }
  }
  solutionNorm_ = maximumAbsElement(solution_, numberTotal);
  // Set bounds and do dj including quadratic
  largeGap = CoinMax(1.0e7, 1.02 * solutionNorm_);
  CoinPackedMatrix *quadratic = NULL;
  const int *columnQuadratic = NULL;
  const CoinBigIndex *columnQuadraticStart = NULL;
  const int *columnQuadraticLength = NULL;
  const double *quadraticElement = NULL;
  if (quadraticObj) {
    quadratic = quadraticObj->quadraticObjective();
    columnQuadratic = quadratic->getIndices();
    columnQuadraticStart = quadratic->getVectorStarts();
    columnQuadraticLength = quadratic->getVectorLengths();
    quadraticElement = quadratic->getElements();
  }
  CoinWorkDouble quadraticNorm = 0.0;
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      CoinWorkDouble primalValue = solution_[iColumn];
      CoinWorkDouble lowerValue = lower_[iColumn];
      CoinWorkDouble upperValue = upper_[iColumn];
      // Do dj
      CoinWorkDouble reducedCost = cost_[iColumn];
      if (lowerBound(iColumn)) {
        reducedCost += linearPerturbation_;
      }
      if (upperBound(iColumn)) {
        reducedCost -= linearPerturbation_;
      }
      if (quadraticObj && iColumn < numberColumns_) {
        for (CoinBigIndex j = columnQuadraticStart[iColumn];
             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
          int jColumn = columnQuadratic[j];
          CoinWorkDouble valueJ = solution_[jColumn];
          CoinWorkDouble elementValue = quadraticElement[j];
          reducedCost += valueJ * elementValue;
        }
        quadraticNorm = CoinMax(quadraticNorm, CoinAbs(reducedCost));
      }
      dj_[iColumn] = reducedCost;
      if (primalValue > lowerValue + largeGap && primalValue < upperValue - largeGap) {
        clearFixedOrFree(iColumn);
        setLowerBound(iColumn);
        setUpperBound(iColumn);
        lowerValue = CoinMax(lowerValue, primalValue - largeGap);
        upperValue = CoinMin(upperValue, primalValue + largeGap);
        lower_[iColumn] = lowerValue;
        upper_[iColumn] = upperValue;
      }
    }
  }
  safeObjectiveValue = CoinMax(safeObjectiveValue, quadraticNorm);
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      CoinWorkDouble primalValue = solution_[iColumn];
      CoinWorkDouble lowerValue = lower_[iColumn];
      CoinWorkDouble upperValue = upper_[iColumn];
      CoinWorkDouble reducedCost = dj_[iColumn];
      CoinWorkDouble low = 0.0;
      CoinWorkDouble high = 0.0;
      if (lowerBound(iColumn)) {
        if (upperBound(iColumn)) {
          //upper and lower bounds
          if (upperValue - lowerValue > 2.0 * initialValue) {
            low = primalValue - lowerValue;
            high = upperValue - primalValue;
          } else {
            low = initialValue;
            high = initialValue;
          }
          CoinWorkDouble s = low + extra;
          CoinWorkDouble ratioZ;
          if (s < zwLarge) {
            ratioZ = 1.0;
          } else {
            ratioZ = CoinSqrt(zwLarge / s);
          }
          CoinWorkDouble t = high + extra;
          CoinWorkDouble ratioT;
          if (t < zwLarge) {
            ratioT = 1.0;
          } else {
            ratioT = CoinSqrt(zwLarge / t);
          }
          //modify s and t
          if (s > largeGap) {
            s = largeGap;
          }
          if (t > largeGap) {
            t = largeGap;
          }
          //modify if long long way away from bound
          if (reducedCost >= 0.0) {
            zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
            zVec_[iColumn] = CoinMax(reducedCost, safeObjectiveValue * ratioZ);
            wVec_[iColumn] = safeObjectiveValue * ratioT;
          } else {
            zVec_[iColumn] = safeObjectiveValue * ratioZ;
            wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
            wVec_[iColumn] = CoinMax(-reducedCost, safeObjectiveValue * ratioT);
          }
          CoinWorkDouble gammaTerm = gamma2;
          if (primalR_)
            gammaTerm += primalR_[iColumn];
          diagonal_[iColumn] = (t * s) / (s * wVec_[iColumn] + t * zVec_[iColumn] + gammaTerm * t * s);
        } else {
          //just lower bound
          low = primalValue - lowerValue;
          high = 0.0;
          CoinWorkDouble s = low + extra;
          CoinWorkDouble ratioZ;
          if (s < zwLarge) {
            ratioZ = 1.0;
          } else {
            ratioZ = CoinSqrt(zwLarge / s);
          }
          //modify s
          if (s > largeGap) {
            s = largeGap;
          }
          if (reducedCost >= 0.0) {
            zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
            zVec_[iColumn] = CoinMax(reducedCost, safeObjectiveValue * ratioZ);
            wVec_[iColumn] = 0.0;
          } else {
            zVec_[iColumn] = safeObjectiveValue * ratioZ;
            wVec_[iColumn] = 0.0;
          }
          CoinWorkDouble gammaTerm = gamma2;
          if (primalR_)
            gammaTerm += primalR_[iColumn];
          diagonal_[iColumn] = s / (zVec_[iColumn] + s * gammaTerm);
        }
      } else {
        if (upperBound(iColumn)) {
          //just upper bound
          low = 0.0;
          high = upperValue - primalValue;
          CoinWorkDouble t = high + extra;
          CoinWorkDouble ratioT;
          if (t < zwLarge) {
            ratioT = 1.0;
          } else {
            ratioT = CoinSqrt(zwLarge / t);
          }
          //modify t
          if (t > largeGap) {
            t = largeGap;
          }
          if (reducedCost >= 0.0) {
            zVec_[iColumn] = 0.0;
            wVec_[iColumn] = safeObjectiveValue * ratioT;
          } else {
            zVec_[iColumn] = 0.0;
            wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
            wVec_[iColumn] = CoinMax(-reducedCost, safeObjectiveValue * ratioT);
          }
          CoinWorkDouble gammaTerm = gamma2;
          if (primalR_)
            gammaTerm += primalR_[iColumn];
          diagonal_[iColumn] = t / (wVec_[iColumn] + t * gammaTerm);
        }
      }
      lowerSlack_[iColumn] = low;
      upperSlack_[iColumn] = high;
    }
  }
#if 0
     if (solution_[0] > 0.0) {
          for (int i = 0; i < numberTotal; i++)
               printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
                      diagonal_[i], CoinAbs(dj_[i]),
                      lowerSlack_[i], zVec_[i],
                      upperSlack_[i], wVec_[i]);
     } else {
          for (int i = 0; i < numberTotal; i++)
               printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
                      diagonal_[i], CoinAbs(dj_[i]),
                      upperSlack_[i], wVec_[i],
                      lowerSlack_[i], zVec_[i] );
     }
     exit(66);
#endif
  return 0;
}
// complementarityGap.  Computes gap
//phase 0=as is , 1 = after predictor , 2 after corrector
CoinWorkDouble ClpPredictorCorrector::complementarityGap(int &numberComplementarityPairs,
  int &numberComplementarityItems,
  const int phase)
{
  CoinWorkDouble gap = 0.0;
  //seems to be same coding for phase = 1 or 2
  numberComplementarityPairs = 0;
  numberComplementarityItems = 0;
  int numberTotal = numberRows_ + numberColumns_;
  CoinWorkDouble toleranceGap = 0.0;
  CoinWorkDouble largestGap = 0.0;
  CoinWorkDouble smallestGap = COIN_DBL_MAX;
  //seems to be same coding for phase = 1 or 2
  int numberNegativeGaps = 0;
  CoinWorkDouble sumNegativeGap = 0.0;
  CoinWorkDouble largeGap = 1.0e2 * solutionNorm_;
  if (largeGap < 1.0e10) {
    largeGap = 1.0e10;
  }
  largeGap = 1.0e30;
  CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
  CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance];
  dualTolerance = dualTolerance / scaleFactor_;
  for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!fixedOrFree(iColumn)) {
      numberComplementarityPairs++;
      //can collapse as if no lower bound both zVec and deltaZ 0.0
      if (lowerBound(iColumn)) {
        numberComplementarityItems++;
        CoinWorkDouble dualValue;
        CoinWorkDouble primalValue;
        if (!phase) {
          dualValue = zVec_[iColumn];
          primalValue = lowerSlack_[iColumn];
        } else {
          CoinWorkDouble change;
          change = solution_[iColumn] + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
          dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
          primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
        }
        //reduce primalValue
        if (primalValue > largeGap) {
          primalValue = largeGap;
        }
        CoinWorkDouble gapProduct = dualValue * primalValue;
        if (gapProduct < 0.0) {
          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
          //primalValue<<endl;
          numberNegativeGaps++;
          sumNegativeGap -= gapProduct;
          gapProduct = 0.0;
        }
        gap += gapProduct;
        //printf("l %d prim %g dual %g totalGap %g\n",
        //   iColumn,primalValue,dualValue,gap);
        if (gapProduct > largestGap) {
          largestGap = gapProduct;
        }
        smallestGap = CoinMin(smallestGap, gapProduct);
        if (dualValue > dualTolerance && primalValue > primalTolerance) {
          toleranceGap += dualValue * primalValue;
        }
      }
      if (upperBound(iColumn)) {
        numberComplementarityItems++;
        CoinWorkDouble dualValue;
        CoinWorkDouble primalValue;
        if (!phase) {
          dualValue = wVec_[iColumn];
          primalValue = upperSlack_[iColumn];
        } else {
          CoinWorkDouble change;
          change = upper_[iColumn] - solution_[iColumn] - deltaX_[iColumn] - upperSlack_[iColumn];
          dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
          primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
        }
        //reduce primalValue
        if (primalValue > largeGap) {
          primalValue = largeGap;
        }
        CoinWorkDouble gapProduct = dualValue * primalValue;
        if (gapProduct < 0.0) {
          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
          //primalValue<<endl;
          numberNegativeGaps++;
          sumNegativeGap -= gapProduct;
          gapProduct = 0.0;
        }
        gap += gapProduct;
        //printf("u %d prim %g dual %g totalGap %g\n",
        //   iColumn,primalValue,dualValue,gap);
        if (gapProduct > largestGap) {
          largestGap = gapProduct;
        }
        if (dualValue > dualTolerance && primalValue > primalTolerance) {
          toleranceGap += dualValue * primalValue;
        }
      }
    }
  }
  //if (numberIterations_>4)
  //exit(9);
  if (!phase && numberNegativeGaps) {
    handler_->message(CLP_BARRIER_NEGATIVE_GAPS, messages_)
      << numberNegativeGaps << static_cast< double >(sumNegativeGap)
      << CoinMessageEol;
  }

  //in case all free!
  if (!numberComplementarityPairs) {
    numberComplementarityPairs = 1;
  }
#ifdef SOME_DEBUG
  printf("with d,p steps %g,%g gap %g - smallest %g, largest %g, pairs %d\n",
    actualDualStep_, actualPrimalStep_,
    gap, smallestGap, largestGap, numberComplementarityPairs);
#endif
  return gap;
}
// setupForSolve.
//phase 0=affine , 1 = corrector , 2 = primal-dual
void ClpPredictorCorrector::setupForSolve(const int phase)
{
  CoinWorkDouble extra = eExtra;
  int numberTotal = numberRows_ + numberColumns_;
  int iColumn;
#ifdef SOME_DEBUG
  printf("phase %d in setupForSolve, mu %.18g\n", phase, mu_);
#endif
  CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
  switch (phase) {
  case 0:
    CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
    if (delta_ || dualR_) {
      // add in regularization
      CoinWorkDouble delta2 = delta_ * delta_;
      for (int iRow = 0; iRow < numberRows_; iRow++) {
        rhsB_[iRow] -= delta2 * dualArray[iRow];
        if (dualR_)
          rhsB_[iRow] -= dualR_[iRow] * dualArray[iRow];
      }
    }
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      rhsC_[iColumn] = 0.0;
      rhsU_[iColumn] = 0.0;
      rhsL_[iColumn] = 0.0;
      rhsZ_[iColumn] = 0.0;
      rhsW_[iColumn] = 0.0;
      if (!flagged(iColumn)) {
        rhsC_[iColumn] = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn];
        rhsC_[iColumn] += gamma2 * solution_[iColumn];
        if (primalR_)
          rhsC_[iColumn] += primalR_[iColumn] * solution_[iColumn];
        if (lowerBound(iColumn)) {
          rhsZ_[iColumn] = -zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
          rhsL_[iColumn] = CoinMax(0.0, (lower_[iColumn] + lowerSlack_[iColumn]) - solution_[iColumn]);
        }
        if (upperBound(iColumn)) {
          rhsW_[iColumn] = -wVec_[iColumn] * (upperSlack_[iColumn] + extra);
          rhsU_[iColumn] = CoinMin(0.0, (upper_[iColumn] - upperSlack_[iColumn]) - solution_[iColumn]);
        }
      }
    }
#if 0
          for (int i = 0; i < 3; i++) {
               if (!CoinAbs(rhsZ_[i]))
                    rhsZ_[i] = 0.0;
               if (!CoinAbs(rhsW_[i]))
                    rhsW_[i] = 0.0;
               if (!CoinAbs(rhsU_[i]))
                    rhsU_[i] = 0.0;
               if (!CoinAbs(rhsL_[i]))
                    rhsL_[i] = 0.0;
          }
          if (solution_[0] > 0.0) {
               for (int i = 0; i < 3; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, solution_[i],
                           diagonal_[i], dj_[i],
                           lowerSlack_[i], zVec_[i],
                           upperSlack_[i], wVec_[i]);
               for (int i = 0; i < 3; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, rhsC_[i],
                           rhsZ_[i], rhsL_[i],
                           rhsW_[i], rhsU_[i]);
          } else {
               for (int i = 0; i < 3; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, solution_[i],
                           diagonal_[i], dj_[i],
                           lowerSlack_[i], zVec_[i],
                           upperSlack_[i], wVec_[i]);
               for (int i = 0; i < 3; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, rhsC_[i],
                           rhsZ_[i], rhsL_[i],
                           rhsW_[i], rhsU_[i]);
          }
#endif
    break;
  case 1:
    // could be stored in delta2?
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      rhsZ_[iColumn] = 0.0;
      rhsW_[iColumn] = 0.0;
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra)
            - deltaZ_[iColumn] * deltaX_[iColumn];
          // To bring in line with OSL
          rhsZ_[iColumn] += deltaZ_[iColumn] * rhsL_[iColumn];
        }
        if (upperBound(iColumn)) {
          rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra)
            + deltaW_[iColumn] * deltaX_[iColumn];
          // To bring in line with OSL
          rhsW_[iColumn] -= deltaW_[iColumn] * rhsU_[iColumn];
        }
      }
    }
#if 0
          for (int i = 0; i < numberTotal; i++) {
               if (!CoinAbs(rhsZ_[i]))
                    rhsZ_[i] = 0.0;
               if (!CoinAbs(rhsW_[i]))
                    rhsW_[i] = 0.0;
               if (!CoinAbs(rhsU_[i]))
                    rhsU_[i] = 0.0;
               if (!CoinAbs(rhsL_[i]))
                    rhsL_[i] = 0.0;
          }
          if (solution_[0] > 0.0) {
               for (int i = 0; i < numberTotal; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
                           diagonal_[i], CoinAbs(dj_[i]),
                           lowerSlack_[i], zVec_[i],
                           upperSlack_[i], wVec_[i]);
               for (int i = 0; i < numberTotal; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(rhsC_[i]),
                           rhsZ_[i], rhsL_[i],
                           rhsW_[i], rhsU_[i]);
          } else {
               for (int i = 0; i < numberTotal; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
                           diagonal_[i], CoinAbs(dj_[i]),
                           upperSlack_[i], wVec_[i],
                           lowerSlack_[i], zVec_[i] );
               for (int i = 0; i < numberTotal; i++)
                    printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(rhsC_[i]),
                           rhsW_[i], rhsU_[i],
                           rhsZ_[i], rhsL_[i]);
          }
          exit(66);
#endif
    break;
  case 2:
    CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      rhsZ_[iColumn] = 0.0;
      rhsW_[iColumn] = 0.0;
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
        }
        if (upperBound(iColumn)) {
          rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra);
        }
      }
    }
    break;
  case 3: {
    CoinWorkDouble minBeta = 0.1 * mu_;
    CoinWorkDouble maxBeta = 10.0 * mu_;
    CoinWorkDouble dualStep = CoinMin(1.0, actualDualStep_ + 0.1);
    CoinWorkDouble primalStep = CoinMin(1.0, actualPrimalStep_ + 0.1);
#ifdef SOME_DEBUG
    printf("good complementarity range %g to %g\n", minBeta, maxBeta);
#endif
    //minBeta=0.0;
    //maxBeta=COIN_DBL_MAX;
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
          CoinWorkDouble dualValue = zVec_[iColumn] + dualStep * deltaZ_[iColumn];
          CoinWorkDouble primalValue = lowerSlack_[iColumn] + primalStep * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          if (gapProduct > 0.0 && dualValue < 0.0)
            gapProduct = -gapProduct;
#ifdef FULL_DEBUG
          delta2Z_[iColumn] = gapProduct;
          if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
            printf("lower %d primal %g, dual %g, gap %g\n",
              iColumn, primalValue, dualValue, gapProduct);
#endif
          CoinWorkDouble value = 0.0;
          if (gapProduct < minBeta) {
            value = 2.0 * (minBeta - gapProduct);
            value = (mu_ - gapProduct);
            value = (minBeta - gapProduct);
            assert(value > 0.0);
          } else if (gapProduct > maxBeta) {
            value = CoinMax(maxBeta - gapProduct, -maxBeta);
            assert(value < 0.0);
          }
          rhsZ_[iColumn] += value;
        }
        if (upperBound(iColumn)) {
          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
          CoinWorkDouble dualValue = wVec_[iColumn] + dualStep * deltaW_[iColumn];
          CoinWorkDouble primalValue = upperSlack_[iColumn] + primalStep * change;
          CoinWorkDouble gapProduct = dualValue * primalValue;
          if (gapProduct > 0.0 && dualValue < 0.0)
            gapProduct = -gapProduct;
#ifdef FULL_DEBUG
          delta2W_[iColumn] = gapProduct;
          if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
            printf("upper %d primal %g, dual %g, gap %g\n",
              iColumn, primalValue, dualValue, gapProduct);
#endif
          CoinWorkDouble value = 0.0;
          if (gapProduct < minBeta) {
            value = (minBeta - gapProduct);
            assert(value > 0.0);
          } else if (gapProduct > maxBeta) {
            value = CoinMax(maxBeta - gapProduct, -maxBeta);
            assert(value < 0.0);
          }
          rhsW_[iColumn] += value;
        }
      }
    }
  } break;
  } /* endswitch */
  if (cholesky_->type() < 20) {
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      CoinWorkDouble value = rhsC_[iColumn];
      CoinWorkDouble zValue = rhsZ_[iColumn];
      CoinWorkDouble wValue = rhsW_[iColumn];
#if 0
#if 1
               if (phase == 0) {
                    // more accurate
                    value = dj[iColumn];
                    zValue = 0.0;
                    wValue = 0.0;
               } else if (phase == 2) {
                    // more accurate
                    value = dj[iColumn];
                    zValue = mu_;
                    wValue = mu_;
               }
#endif
               assert (rhsL_[iColumn] >= 0.0);
               assert (rhsU_[iColumn] <= 0.0);
               if (lowerBound(iColumn)) {
                    value += (-zVec_[iColumn] * rhsL_[iColumn] - zValue) /
                             (lowerSlack_[iColumn] + extra);
               }
               if (upperBound(iColumn)) {
                    value += (wValue - wVec_[iColumn] * rhsU_[iColumn]) /
                             (upperSlack_[iColumn] + extra);
               }
#else
      if (lowerBound(iColumn)) {
        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
        value -= gHat / (lowerSlack_[iColumn] + extra);
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
        value += hHat / (upperSlack_[iColumn] + extra);
      }
#endif
      workArray_[iColumn] = diagonal_[iColumn] * value;
    }
#if 0
          if (solution_[0] > 0.0) {
               for (int i = 0; i < numberTotal; i++)
                    printf("%d %.18g\n", i, workArray_[i]);
          } else {
               for (int i = 0; i < numberTotal; i++)
                    printf("%d %.18g\n", i, workArray_[i]);
          }
          exit(66);
#endif
  } else {
    // KKT
    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
      CoinWorkDouble value = rhsC_[iColumn];
      CoinWorkDouble zValue = rhsZ_[iColumn];
      CoinWorkDouble wValue = rhsW_[iColumn];
      if (lowerBound(iColumn)) {
        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
        value -= gHat / (lowerSlack_[iColumn] + extra);
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
        value += hHat / (upperSlack_[iColumn] + extra);
      }
      workArray_[iColumn] = value;
    }
  }
}
//method: sees if looks plausible change in complementarity
bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
  CoinWorkDouble &bestNextGap,
  bool allowIncreasingGap)
{
  const CoinWorkDouble beta3 = 0.99997;
  bool goodMove = false;
  int nextNumber;
  int nextNumberItems;
  int numberTotal = numberRows_ + numberColumns_;
  CoinWorkDouble returnGap = bestNextGap;
  CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (objective_->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
#endif
  if (nextGap > bestNextGap && nextGap > 0.9 * complementarityGap_ && doCorrector
    && !quadraticObj && !allowIncreasingGap) {
#ifdef SOME_DEBUG
    printf("checkGood phase 1 next gap %.18g, phase 0 %.18g, old gap %.18g\n",
      nextGap, bestNextGap, complementarityGap_);
#endif
    return false;
  } else {
    returnGap = nextGap;
  }
  CoinWorkDouble step;
  if (actualDualStep_ > actualPrimalStep_) {
    step = actualDualStep_;
  } else {
    step = actualPrimalStep_;
  }
  CoinWorkDouble testValue = 1.0 - step * (1.0 - beta3);
  //testValue=0.0;
  testValue *= complementarityGap_;
  if (nextGap < testValue) {
    //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
    goodMove = true;
  } else if (doCorrector) {
    //if (actualDualStep_<actualPrimalStep_) {
    //step=actualDualStep_;
    //} else {
    //step=actualPrimalStep_;
    //}
    CoinWorkDouble gap = bestNextGap;
    goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
    if (goodMove)
      returnGap = gap;
  } else {
    goodMove = true;
  }
  if (goodMove)
    goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
  // Say good if small
  //if (quadraticObj) {
  if (CoinMax(actualDualStep_, actualPrimalStep_) < 1.0e-6)
    goodMove = true;
  if (!goodMove) {
    //try smaller of two
    if (actualDualStep_ < actualPrimalStep_) {
      step = actualDualStep_;
    } else {
      step = actualPrimalStep_;
    }
    if (step > 1.0) {
      step = 1.0;
    }
    actualPrimalStep_ = step;
    //if (quadraticObj)
    //actualPrimalStep_ *=0.5;
    actualDualStep_ = step;
    goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
    int pass = 0;
    while (!goodMove) {
      pass++;
      CoinWorkDouble gap = bestNextGap;
      goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
      if (goodMove || pass > 3) {
        returnGap = gap;
        break;
      }
      if (step < 1.0e-4) {
        break;
      }
      step *= 0.5;
      actualPrimalStep_ = step;
      //if (quadraticObj)
      //actualPrimalStep_ *=0.5;
      actualDualStep_ = step;
    } /* endwhile */
    if (doCorrector) {
      //say bad move if both small
      if (numberIterations_ & 1) {
        if (actualPrimalStep_ < 1.0e-2 && actualDualStep_ < 1.0e-2) {
          goodMove = false;
        }
      } else {
        if (actualPrimalStep_ < 1.0e-5 && actualDualStep_ < 1.0e-5) {
          goodMove = false;
        }
        if (actualPrimalStep_ * actualDualStep_ < 1.0e-20) {
          goodMove = false;
        }
      }
    }
  }
  if (goodMove) {
    //compute delta in objectives
    CoinWorkDouble deltaObjectivePrimal = 0.0;
    CoinWorkDouble deltaObjectiveDual = innerProduct(deltaY_, numberRows_,
      rhsFixRegion_);
    CoinWorkDouble error = 0.0;
    CoinWorkDouble *workArray = workArray_;
    CoinZeroN(workArray, numberColumns_);
    CoinMemcpyN(deltaY_, numberRows_, workArray + numberColumns_);
    matrix_->transposeTimes(-1.0, deltaY_, workArray);
    //CoinWorkDouble sumPerturbCost=0.0;
    for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
      if (!flagged(iColumn)) {
        if (lowerBound(iColumn)) {
          //sumPerturbCost+=deltaX_[iColumn];
          deltaObjectiveDual += deltaZ_[iColumn] * lower_[iColumn];
        }
        if (upperBound(iColumn)) {
          //sumPerturbCost-=deltaX_[iColumn];
          deltaObjectiveDual -= deltaW_[iColumn] * upper_[iColumn];
        }
        CoinWorkDouble change = CoinAbs(workArray_[iColumn] - deltaZ_[iColumn] + deltaW_[iColumn]);
        error = CoinMax(change, error);
      }
      deltaObjectivePrimal += cost_[iColumn] * deltaX_[iColumn];
    }
    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
    CoinWorkDouble testValue;
    if (error > 0.0) {
      testValue = 1.0e1 * CoinMax(maximumDualError_, 1.0e-12) / error;
    } else {
      testValue = 1.0e1;
    }
    // If quadratic then primal step may compensate
    if (testValue < actualDualStep_ && !quadraticObj) {
      handler_->message(CLP_BARRIER_REDUCING, messages_)
        << "dual" << static_cast< double >(actualDualStep_)
        << static_cast< double >(testValue)
        << CoinMessageEol;
      actualDualStep_ = testValue;
    }
  }
  if (maximumRHSError_ < 1.0e1 * solutionNorm_ * primalTolerance()
    && maximumRHSChange_ > 1.0e-16 * solutionNorm_) {
    //check change in AX not too much
    //??? could be dropped row going infeasible
    CoinWorkDouble ratio = 1.0e1 * CoinMax(maximumRHSError_, 1.0e-12) / maximumRHSChange_;
    if (ratio < actualPrimalStep_) {
      handler_->message(CLP_BARRIER_REDUCING, messages_)
        << "primal" << static_cast< double >(actualPrimalStep_)
        << static_cast< double >(ratio)
        << CoinMessageEol;
      if (ratio > 1.0e-6) {
        actualPrimalStep_ = ratio;
      } else {
        actualPrimalStep_ = ratio;
        //std::cout <<"sign we should be stopping"<<std::endl;
      }
    }
  }
  if (goodMove)
    bestNextGap = returnGap;
  return goodMove;
}
//:  checks for one step size
bool ClpPredictorCorrector::checkGoodMove2(CoinWorkDouble move,
  CoinWorkDouble &bestNextGap,
  bool allowIncreasingGap)
{
  CoinWorkDouble complementarityMultiplier = 1.0 / numberComplementarityPairs_;
  const CoinWorkDouble gamma = 1.0e-8;
  const CoinWorkDouble gammap = 1.0e-8;
  CoinWorkDouble gammad = 1.0e-8;
  int nextNumber;
  int nextNumberItems;
  CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
  if (nextGap > bestNextGap && !allowIncreasingGap)
    return false;
  CoinWorkDouble lowerBoundGap = gamma * nextGap * complementarityMultiplier;
  bool goodMove = true;
  int iColumn;
  for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    if (!flagged(iColumn)) {
      if (lowerBound(iColumn)) {
        CoinWorkDouble part1 = lowerSlack_[iColumn] + actualPrimalStep_ * deltaSL_[iColumn];
        CoinWorkDouble part2 = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
        if (part1 * part2 < lowerBoundGap) {
          goodMove = false;
          break;
        }
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble part1 = upperSlack_[iColumn] + actualPrimalStep_ * deltaSU_[iColumn];
        CoinWorkDouble part2 = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
        if (part1 * part2 < lowerBoundGap) {
          goodMove = false;
          break;
        }
      }
    }
  }
  CoinWorkDouble *nextDj = NULL;
  CoinWorkDouble maximumDualError = maximumDualError_;
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (objective_->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
#endif
  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
  if (quadraticObj) {
    // change gammad
    gammad = 1.0e-4;
    CoinWorkDouble gamma2 = gamma_ * gamma_;
    nextDj = new CoinWorkDouble[numberColumns_];
    CoinWorkDouble *nextSolution = new CoinWorkDouble[numberColumns_];
    // put next primal into nextSolution
    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
      if (!flagged(iColumn)) {
        nextSolution[iColumn] = solution_[iColumn] + actualPrimalStep_ * deltaX_[iColumn];
      } else {
        nextSolution[iColumn] = solution_[iColumn];
      }
    }
    // do reduced costs
    CoinMemcpyN(cost_, numberColumns_, nextDj);
    matrix_->transposeTimes(-1.0, dualArray, nextDj);
    matrix_->transposeTimes(-actualDualStep_, deltaY_, nextDj);
    quadraticDjs(nextDj, nextSolution, 1.0);
    delete[] nextSolution;
    CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
    const int *columnQuadraticLength = quadratic->getVectorLengths();
    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
      if (!fixedOrFree(iColumn)) {
        CoinWorkDouble newZ = 0.0;
        CoinWorkDouble newW = 0.0;
        if (lowerBound(iColumn)) {
          newZ = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
        }
        if (upperBound(iColumn)) {
          newW = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
        }
        if (columnQuadraticLength[iColumn]) {
          CoinWorkDouble gammaTerm = gamma2;
          if (primalR_)
            gammaTerm += primalR_[iColumn];
          //CoinWorkDouble dualInfeasibility=
          //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
          //+gammaTerm*solution_[iColumn];
          CoinWorkDouble newInfeasibility = nextDj[iColumn] - newZ + newW
            + gammaTerm * (solution_[iColumn] + actualPrimalStep_ * deltaX_[iColumn]);
          maximumDualError = CoinMax(maximumDualError, newInfeasibility);
          //if (CoinAbs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
          //if (dualInfeasibility*newInfeasibility<0.0) {
          //  printf("%d current %g next %g\n",iColumn,dualInfeasibility,
          //       newInfeasibility);
          //  goodMove=false;
          //}
          //}
        }
      }
    }
    delete[] nextDj;
  }
  //      Satisfy g_p(alpha)?
  if (rhsNorm_ > solutionNorm_) {
    solutionNorm_ = rhsNorm_;
  }
  CoinWorkDouble errorCheck = maximumRHSError_ / solutionNorm_;
  if (errorCheck < maximumBoundInfeasibility_) {
    errorCheck = maximumBoundInfeasibility_;
  }
  // scale back move
  move = CoinMin(move, 0.95);
  //scale
  if ((1.0 - move) * errorCheck > primalTolerance()) {
    if (nextGap < gammap * (1.0 - move) * errorCheck) {
      goodMove = false;
    }
  }
  //      Satisfy g_d(alpha)?
  errorCheck = maximumDualError / objectiveNorm_;
  if ((1.0 - move) * errorCheck > dualTolerance()) {
    if (nextGap < gammad * (1.0 - move) * errorCheck) {
      goodMove = false;
    }
  }
  if (goodMove)
    bestNextGap = nextGap;
  return goodMove;
}
// updateSolution.  Updates solution at end of iteration
//returns number fixed
int ClpPredictorCorrector::updateSolution(CoinWorkDouble /*nextGap*/)
{
  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
  int numberTotal = numberRows_ + numberColumns_;
  //update pi
  multiplyAdd(deltaY_, numberRows_, actualDualStep_, dualArray, 1.0);
  CoinZeroN(errorRegion_, numberRows_);
  CoinZeroN(rhsFixRegion_, numberRows_);
  CoinWorkDouble maximumRhsInfeasibility = 0.0;
  CoinWorkDouble maximumBoundInfeasibility = 0.0;
  CoinWorkDouble maximumDualError = 1.0e-12;
  CoinWorkDouble primalObjectiveValue = 0.0;
  CoinWorkDouble dualObjectiveValue = 0.0;
  CoinWorkDouble solutionNorm = 1.0e-12;
  int numberKilled = 0;
  CoinWorkDouble freeMultiplier = 1.0e6;
  CoinWorkDouble trueNorm = diagonalNorm_ / diagonalScaleFactor_;
  if (freeMultiplier < trueNorm) {
    freeMultiplier = trueNorm;
  }
  if (freeMultiplier > 1.0e12) {
    freeMultiplier = 1.0e12;
  }
  freeMultiplier = 0.5 / freeMultiplier;
  CoinWorkDouble condition = CoinAbs(cholesky_->choleskyCondition());
  bool caution;
  if ((condition < 1.0e10 && trueNorm < 1.0e12) || numberIterations_ < 20) {
    caution = false;
  } else {
    caution = true;
  }
  CoinWorkDouble extra = eExtra;
  const CoinWorkDouble largeFactor = 1.0e2;
  CoinWorkDouble largeGap = largeFactor * solutionNorm_;
  if (largeGap < largeFactor) {
    largeGap = largeFactor;
  }
  CoinWorkDouble dualFake = 0.0;
  CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
  dualTolerance = dualTolerance / scaleFactor_;
  if (dualTolerance < 1.0e-12) {
    dualTolerance = 1.0e-12;
  }
  CoinWorkDouble offsetObjective = 0.0;
  CoinWorkDouble killTolerance = primalTolerance();
  //CoinWorkDouble nextMu = nextGap/(static_cast<CoinWorkDouble>(2*numberComplementarityPairs_));
  //printf("using gap of %g\n",nextMu);
  //largest allowable ratio of lowerSlack/zVec (etc)
  CoinWorkDouble epsilonBase;
  CoinWorkDouble diagonalLimit;
  if (!caution) {
    epsilonBase = eBase;
    diagonalLimit = eDiagonal;
  } else {
    epsilonBase = eBaseCaution;
    diagonalLimit = eDiagonalCaution;
  }
  CoinWorkDouble maximumDJInfeasibility = 0.0;
  int numberIncreased = 0;
  int numberDecreased = 0;
  CoinWorkDouble largestDiagonal = 0.0;
  CoinWorkDouble smallestDiagonal = 1.0e50;
  CoinWorkDouble largeGap2 = CoinMax(1.0e7, 1.0e2 * solutionNorm_);
  //largeGap2 = 1.0e9;
  // When to start looking at killing (factor0
  CoinWorkDouble killFactor;
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (objective_->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
#endif
#ifndef CLP_CAUTION
#define KILL_ITERATION 50
#else
#if CLP_CAUTION < 1
#define KILL_ITERATION 50
#else
#define KILL_ITERATION 100
#endif
#endif
  if (!quadraticObj || 1) {
    if (numberIterations_ < KILL_ITERATION) {
      killFactor = 1.0;
    } else if (numberIterations_ < 2 * KILL_ITERATION) {
      killFactor = 5.0;
      stepLength_ = CoinMax(stepLength_, 0.9995);
    } else if (numberIterations_ < 4 * KILL_ITERATION) {
      killFactor = 20.0;
      stepLength_ = CoinMax(stepLength_, 0.99995);
    } else {
      killFactor = 1.0e2;
      stepLength_ = CoinMax(stepLength_, 0.999995);
    }
  } else {
    killFactor = 1.0;
  }
  // put next primal into deltaSL_
  int iColumn;
  int iRow;
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    CoinWorkDouble thisWeight = deltaX_[iColumn];
    CoinWorkDouble newPrimal = solution_[iColumn] + 1.0 * actualPrimalStep_ * thisWeight;
    deltaSL_[iColumn] = newPrimal;
  }
#if 0
     // nice idea but doesn't work
     multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
     matrix_->times(1.0, solution_, errorRegion_);
     multiplyAdd(deltaSL_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 0.0);
     matrix_->times(1.0, deltaSL_, rhsFixRegion_);
     CoinWorkDouble newNorm =  maximumAbsElement(deltaSL_, numberTotal);
     CoinWorkDouble tol = newNorm * primalTolerance();
     bool goneInf = false;
     for (iRow = 0; iRow < numberRows_; iRow++) {
          CoinWorkDouble value = errorRegion_[iRow];
          CoinWorkDouble valueNew = rhsFixRegion_[iRow];
          if (CoinAbs(value) < tol && CoinAbs(valueNew) > tol) {
               printf("row %d old %g new %g\n", iRow, value, valueNew);
               goneInf = true;
          }
     }
     if (goneInf) {
          actualPrimalStep_ *= 0.5;
          for (iColumn = 0; iColumn < numberTotal; iColumn++) {
               CoinWorkDouble thisWeight = deltaX_[iColumn];
               CoinWorkDouble newPrimal = solution_[iColumn] + 1.0 * actualPrimalStep_ * thisWeight;
               deltaSL_[iColumn] = newPrimal;
          }
     }
     CoinZeroN(errorRegion_, numberRows_);
     CoinZeroN(rhsFixRegion_, numberRows_);
#endif
  // do reduced costs
  CoinMemcpyN(dualArray, numberRows_, dj_ + numberColumns_);
  CoinMemcpyN(cost_, numberColumns_, dj_);
  CoinWorkDouble quadraticOffset = quadraticDjs(dj_, deltaSL_, 1.0);
  // Save modified costs for fixed variables
  CoinMemcpyN(dj_, numberColumns_, deltaSU_);
  matrix_->transposeTimes(-1.0, dualArray, dj_);
  CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
  CoinWorkDouble gammaOffset = 0.0;
#if 0
     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     const int * columnLength = matrix_->getVectorLengths();
     const int * row = matrix_->getIndices();
     const double * element = matrix_->getElements();
#endif
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      CoinWorkDouble reducedCost = dj_[iColumn];
      bool thisKilled = false;
      CoinWorkDouble zValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
      CoinWorkDouble wValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
      zVec_[iColumn] = zValue;
      wVec_[iColumn] = wValue;
      CoinWorkDouble thisWeight = deltaX_[iColumn];
      CoinWorkDouble oldPrimal = solution_[iColumn];
      CoinWorkDouble newPrimal = solution_[iColumn] + actualPrimalStep_ * thisWeight;
      CoinWorkDouble dualObjectiveThis = 0.0;
      CoinWorkDouble sUpper = extra;
      CoinWorkDouble sLower = extra;
      CoinWorkDouble kill;
      if (CoinAbs(newPrimal) > 1.0e4) {
        kill = killTolerance * 1.0e-4 * newPrimal;
      } else {
        kill = killTolerance;
      }
      kill *= 1.0e-3; //be conservative
      CoinWorkDouble smallerSlack = COIN_DBL_MAX;
      bool fakeOldBounds = false;
      bool fakeNewBounds = false;
      CoinWorkDouble trueLower;
      CoinWorkDouble trueUpper;
      if (iColumn < numberColumns_) {
        trueLower = columnLower_[iColumn];
        trueUpper = columnUpper_[iColumn];
      } else {
        trueLower = rowLower_[iColumn - numberColumns_];
        trueUpper = rowUpper_[iColumn - numberColumns_];
      }
      if (oldPrimal > trueLower + largeGap2 && oldPrimal < trueUpper - largeGap2)
        fakeOldBounds = true;
      if (newPrimal > trueLower + largeGap2 && newPrimal < trueUpper - largeGap2)
        fakeNewBounds = true;
      if (fakeOldBounds) {
        if (fakeNewBounds) {
          lower_[iColumn] = newPrimal - largeGap2;
          lowerSlack_[iColumn] = largeGap2;
          upper_[iColumn] = newPrimal + largeGap2;
          upperSlack_[iColumn] = largeGap2;
        } else {
          lower_[iColumn] = trueLower;
          setLowerBound(iColumn);
          lowerSlack_[iColumn] = CoinMax(newPrimal - trueLower, 1.0);
          upper_[iColumn] = trueUpper;
          setUpperBound(iColumn);
          upperSlack_[iColumn] = CoinMax(trueUpper - newPrimal, 1.0);
        }
      } else if (fakeNewBounds) {
        lower_[iColumn] = newPrimal - largeGap2;
        lowerSlack_[iColumn] = largeGap2;
        upper_[iColumn] = newPrimal + largeGap2;
        upperSlack_[iColumn] = largeGap2;
        // so we can just have one test
        fakeOldBounds = true;
      }
      CoinWorkDouble lowerBoundInfeasibility = 0.0;
      CoinWorkDouble upperBoundInfeasibility = 0.0;
      //double saveNewPrimal = newPrimal;
      if (lowerBound(iColumn)) {
        CoinWorkDouble oldSlack = lowerSlack_[iColumn];
        CoinWorkDouble newSlack;
        newSlack = lowerSlack_[iColumn] + actualPrimalStep_ * (oldPrimal - oldSlack + thisWeight - lower_[iColumn]);
        if (fakeOldBounds)
          newSlack = lowerSlack_[iColumn];
        CoinWorkDouble epsilon = CoinAbs(newSlack) * epsilonBase;
        epsilon = CoinMin(epsilon, 1.0e-5);
        //epsilon=1.0e-14;
        //make sure reasonable
        if (zValue < epsilon) {
          zValue = epsilon;
        }
        CoinWorkDouble feasibleSlack = newPrimal - lower_[iColumn];
        if (feasibleSlack > 0.0 && newSlack > 0.0) {
          CoinWorkDouble larger;
          if (newSlack > feasibleSlack) {
            larger = newSlack;
          } else {
            larger = feasibleSlack;
          }
          if (CoinAbs(feasibleSlack - newSlack) < 1.0e-6 * larger) {
            newSlack = feasibleSlack;
          }
        }
        if (zVec_[iColumn] > dualTolerance) {
          dualObjectiveThis += lower_[iColumn] * zVec_[iColumn];
        }
        lowerSlack_[iColumn] = newSlack;
        if (newSlack < smallerSlack) {
          smallerSlack = newSlack;
        }
        lowerBoundInfeasibility = CoinAbs(newPrimal - lowerSlack_[iColumn] - lower_[iColumn]);
        if (lowerSlack_[iColumn] <= kill * killFactor && CoinAbs(newPrimal - lower_[iColumn]) <= kill * killFactor) {
          CoinWorkDouble step = CoinMin(actualPrimalStep_ * 1.1, 1.0);
          CoinWorkDouble newPrimal2 = solution_[iColumn] + step * thisWeight;
          if (newPrimal2 < newPrimal && dj_[iColumn] > 1.0e-5 && numberIterations_ > 50 - 40) {
            newPrimal = lower_[iColumn];
            lowerSlack_[iColumn] = 0.0;
            //printf("fixing %d to lower\n",iColumn);
          }
        }
        if (lowerSlack_[iColumn] <= kill && CoinAbs(newPrimal - lower_[iColumn]) <= kill) {
          //may be better to leave at value?
          newPrimal = lower_[iColumn];
          lowerSlack_[iColumn] = 0.0;
          thisKilled = true;
          //cout<<j<<" l "<<reducedCost<<" "<<zVec_[iColumn]<<endl;
        } else {
          sLower += lowerSlack_[iColumn];
        }
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble oldSlack = upperSlack_[iColumn];
        CoinWorkDouble newSlack;
        newSlack = upperSlack_[iColumn] + actualPrimalStep_ * (-oldPrimal - oldSlack - thisWeight + upper_[iColumn]);
        if (fakeOldBounds)
          newSlack = upperSlack_[iColumn];
        CoinWorkDouble epsilon = CoinAbs(newSlack) * epsilonBase;
        epsilon = CoinMin(epsilon, 1.0e-5);
        //make sure reasonable
        //epsilon=1.0e-14;
        if (wValue < epsilon) {
          wValue = epsilon;
        }
        CoinWorkDouble feasibleSlack = upper_[iColumn] - newPrimal;
        if (feasibleSlack > 0.0 && newSlack > 0.0) {
          CoinWorkDouble larger;
          if (newSlack > feasibleSlack) {
            larger = newSlack;
          } else {
            larger = feasibleSlack;
          }
          if (CoinAbs(feasibleSlack - newSlack) < 1.0e-6 * larger) {
            newSlack = feasibleSlack;
          }
        }
        if (wVec_[iColumn] > dualTolerance) {
          dualObjectiveThis -= upper_[iColumn] * wVec_[iColumn];
        }
        upperSlack_[iColumn] = newSlack;
        if (newSlack < smallerSlack) {
          smallerSlack = newSlack;
        }
        upperBoundInfeasibility = CoinAbs(newPrimal + upperSlack_[iColumn] - upper_[iColumn]);
        if (upperSlack_[iColumn] <= kill * killFactor && CoinAbs(newPrimal - upper_[iColumn]) <= kill * killFactor) {
          CoinWorkDouble step = CoinMin(actualPrimalStep_ * 1.1, 1.0);
          CoinWorkDouble newPrimal2 = solution_[iColumn] + step * thisWeight;
          if (newPrimal2 > newPrimal && dj_[iColumn] < -1.0e-5 && numberIterations_ > 50 - 40) {
            newPrimal = upper_[iColumn];
            upperSlack_[iColumn] = 0.0;
            //printf("fixing %d to upper\n",iColumn);
          }
        }
        if (upperSlack_[iColumn] <= kill && CoinAbs(newPrimal - upper_[iColumn]) <= kill) {
          //may be better to leave at value?
          newPrimal = upper_[iColumn];
          upperSlack_[iColumn] = 0.0;
          thisKilled = true;
        } else {
          sUpper += upperSlack_[iColumn];
        }
      }
#if 0
               if (newPrimal != saveNewPrimal && iColumn < numberColumns_) {
                    // adjust slacks
                    double movement = newPrimal - saveNewPrimal;
                    for (CoinBigIndex j = columnStart[iColumn];
                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                         int iRow = row[j];
                         double slackMovement = element[j] * movement;
                         solution_[iRow+numberColumns_] += slackMovement; // sign?
                    }
               }
#endif
      solution_[iColumn] = newPrimal;
      if (CoinAbs(newPrimal) > solutionNorm) {
        solutionNorm = CoinAbs(newPrimal);
      }
      if (!thisKilled) {
        CoinWorkDouble gammaTerm = gamma2;
        if (primalR_) {
          gammaTerm += primalR_[iColumn];
          quadraticOffset += newPrimal * newPrimal * primalR_[iColumn];
        }
        CoinWorkDouble dualInfeasibility = reducedCost - zVec_[iColumn] + wVec_[iColumn] + gammaTerm * newPrimal;
        if (CoinAbs(dualInfeasibility) > dualTolerance) {
#if 0
                         if (dualInfeasibility > 0.0) {
                              // To improve we could reduce t and/or increase z
                              if (lowerBound(iColumn)) {
                                   CoinWorkDouble complementarity = zVec_[iColumn] * lowerSlack_[iColumn];
                                   if (complementarity < nextMu) {
                                        CoinWorkDouble change =
                                             CoinMin(dualInfeasibility,
                                                     (nextMu - complementarity) / lowerSlack_[iColumn]);
                                        dualInfeasibility -= change;
                                        COIN_DETAIL_PRINT(printf("%d lb locomp %g - dual inf from %g to %g\n",
                                               iColumn, complementarity, dualInfeasibility + change,
								 dualInfeasibility));
                                        zVec_[iColumn] += change;
                                        zValue = CoinMax(zVec_[iColumn], 1.0e-12);
                                   }
                              }
                              if (upperBound(iColumn)) {
                                   CoinWorkDouble complementarity = wVec_[iColumn] * upperSlack_[iColumn];
                                   if (complementarity > nextMu) {
                                        CoinWorkDouble change =
                                             CoinMin(dualInfeasibility,
                                                     (complementarity - nextMu) / upperSlack_[iColumn]);
                                        dualInfeasibility -= change;
                                        COIN_DETAIL_PRINT(printf("%d ub hicomp %g - dual inf from %g to %g\n",
                                               iColumn, complementarity, dualInfeasibility + change,
								 dualInfeasibility));
                                        wVec_[iColumn] -= change;
                                        wValue = CoinMax(wVec_[iColumn], 1.0e-12);
                                   }
                              }
                         } else {
                              // To improve we could reduce z and/or increase t
                              if (lowerBound(iColumn)) {
                                   CoinWorkDouble complementarity = zVec_[iColumn] * lowerSlack_[iColumn];
                                   if (complementarity > nextMu) {
                                        CoinWorkDouble change =
                                             CoinMax(dualInfeasibility,
                                                     (nextMu - complementarity) / lowerSlack_[iColumn]);
                                        dualInfeasibility -= change;
                                        COIN_DETAIL_PRINT(printf("%d lb hicomp %g - dual inf from %g to %g\n",
                                               iColumn, complementarity, dualInfeasibility + change,
								 dualInfeasibility));
                                        zVec_[iColumn] += change;
                                        zValue = CoinMax(zVec_[iColumn], 1.0e-12);
                                   }
                              }
                              if (upperBound(iColumn)) {
                                   CoinWorkDouble complementarity = wVec_[iColumn] * upperSlack_[iColumn];
                                   if (complementarity < nextMu) {
                                        CoinWorkDouble change =
                                             CoinMax(dualInfeasibility,
                                                     (complementarity - nextMu) / upperSlack_[iColumn]);
                                        dualInfeasibility -= change;
                                        COIN_DETAIL_PRINT(printf("%d ub locomp %g - dual inf from %g to %g\n",
                                               iColumn, complementarity, dualInfeasibility + change,
								 dualInfeasibility));
                                        wVec_[iColumn] -= change;
                                        wValue = CoinMax(wVec_[iColumn], 1.0e-12);
                                   }
                              }
                         }
#endif
          dualFake += newPrimal * dualInfeasibility;
        }
        if (lowerBoundInfeasibility > maximumBoundInfeasibility) {
          maximumBoundInfeasibility = lowerBoundInfeasibility;
        }
        if (upperBoundInfeasibility > maximumBoundInfeasibility) {
          maximumBoundInfeasibility = upperBoundInfeasibility;
        }
        dualInfeasibility = CoinAbs(dualInfeasibility);
        if (dualInfeasibility > maximumDualError) {
          //printf("bad dual %d %g\n",iColumn,
          // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
          maximumDualError = dualInfeasibility;
        }
        dualObjectiveValue += dualObjectiveThis;
        gammaOffset += newPrimal * newPrimal;
        if (sLower > largeGap) {
          sLower = largeGap;
        }
        if (sUpper > largeGap) {
          sUpper = largeGap;
        }
#if 1
        CoinWorkDouble divisor = sLower * wValue + sUpper * zValue + gammaTerm * sLower * sUpper;
        CoinWorkDouble diagonalValue = (sUpper * sLower) / divisor;
#else
        CoinWorkDouble divisor = sLower * wValue + sUpper * zValue + gammaTerm * sLower * sUpper;
        CoinWorkDouble diagonalValue2 = (sUpper * sLower) / divisor;
        CoinWorkDouble diagonalValue;
        if (!lowerBound(iColumn)) {
          diagonalValue = wValue / sUpper + gammaTerm;
        } else if (!upperBound(iColumn)) {
          diagonalValue = zValue / sLower + gammaTerm;
        } else {
          diagonalValue = zValue / sLower + wValue / sUpper + gammaTerm;
        }
        diagonalValue = 1.0 / diagonalValue;
#endif
        diagonal_[iColumn] = diagonalValue;
        //FUDGE
        if (diagonalValue > diagonalLimit) {
#ifdef COIN_DEVELOP
          std::cout << "large diagonal " << diagonalValue << std::endl;
#endif
          diagonal_[iColumn] = diagonalLimit;
        }
#ifdef COIN_DEVELOP
        if (diagonalValue < 1.0e-10) {
          //std::cout<<"small diagonal "<<diagonalValue<<std::endl;
        }
#endif
        if (diagonalValue > largestDiagonal) {
          largestDiagonal = diagonalValue;
        }
        if (diagonalValue < smallestDiagonal) {
          smallestDiagonal = diagonalValue;
        }
        deltaX_[iColumn] = 0.0;
      } else {
        numberKilled++;
        if (solution_[iColumn] != lower_[iColumn] && solution_[iColumn] != upper_[iColumn]) {
          COIN_DETAIL_PRINT(printf("%d %g %g %g\n", iColumn, static_cast< double >(lower_[iColumn]),
            static_cast< double >(solution_[iColumn]), static_cast< double >(upper_[iColumn])));
        }
        diagonal_[iColumn] = 0.0;
        zVec_[iColumn] = 0.0;
        wVec_[iColumn] = 0.0;
        setFlagged(iColumn);
        setFixedOrFree(iColumn);
        deltaX_[iColumn] = newPrimal;
        offsetObjective += newPrimal * deltaSU_[iColumn];
      }
    } else {
      deltaX_[iColumn] = solution_[iColumn];
      diagonal_[iColumn] = 0.0;
      offsetObjective += solution_[iColumn] * deltaSU_[iColumn];
      if (upper_[iColumn] - lower_[iColumn] > 1.0e-5) {
        if (solution_[iColumn] < lower_[iColumn] + 1.0e-8 && dj_[iColumn] < -1.0e-8) {
          if (-dj_[iColumn] > maximumDJInfeasibility) {
            maximumDJInfeasibility = -dj_[iColumn];
          }
        }
        if (solution_[iColumn] > upper_[iColumn] - 1.0e-8 && dj_[iColumn] > 1.0e-8) {
          if (dj_[iColumn] > maximumDJInfeasibility) {
            maximumDJInfeasibility = dj_[iColumn];
          }
        }
      }
    }
    primalObjectiveValue += solution_[iColumn] * cost_[iColumn];
  }
  handler_->message(CLP_BARRIER_DIAGONAL, messages_)
    << static_cast< double >(largestDiagonal) << static_cast< double >(smallestDiagonal)
    << CoinMessageEol;
#if 0
     // If diagonal wild - kill some
     if (largestDiagonal > 1.0e17 * smallestDiagonal) {
          CoinWorkDouble killValue = largestDiagonal * 1.0e-17;
          for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
               if (CoinAbs(diagonal_[iColumn]) < killValue)
                    diagonal_[iolumn] = 0.0;
          }
     }
#endif
  // update rhs region
  multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 1.0);
  matrix_->times(1.0, deltaX_, rhsFixRegion_);
  primalObjectiveValue += 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
  if (quadraticOffset) {
    //  printf("gamma offset %g %g, quadoffset %g\n",gammaOffset,gamma2*gammaOffset,quadraticOffset);
  }

  dualObjectiveValue += offsetObjective + dualFake;
  dualObjectiveValue -= 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
  if (numberIncreased || numberDecreased) {
    handler_->message(CLP_BARRIER_SLACKS, messages_)
      << numberIncreased << numberDecreased
      << CoinMessageEol;
  }
  if (maximumDJInfeasibility) {
    handler_->message(CLP_BARRIER_DUALINF, messages_)
      << static_cast< double >(maximumDJInfeasibility)
      << CoinMessageEol;
  }
  // Need to rethink (but it is only for printing)
  sumPrimalInfeasibilities_ = maximumRhsInfeasibility;
  sumDualInfeasibilities_ = maximumDualError;
  maximumBoundInfeasibility_ = maximumBoundInfeasibility;
  //compute error and fixed RHS
  multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
  matrix_->times(1.0, solution_, errorRegion_);
  maximumDualError_ = maximumDualError;
  maximumBoundInfeasibility_ = maximumBoundInfeasibility;
  solutionNorm_ = solutionNorm;
  //finish off objective computation
  primalObjective_ = primalObjectiveValue * scaleFactor_;
  CoinWorkDouble dualValue2 = innerProduct(dualArray, numberRows_,
    rhsFixRegion_);
  dualObjectiveValue -= dualValue2;
  dualObjective_ = dualObjectiveValue * scaleFactor_;
  if (numberKilled) {
    handler_->message(CLP_BARRIER_KILLED, messages_)
      << numberKilled
      << CoinMessageEol;
  }
  CoinWorkDouble maximumRHSError1 = 0.0;
  CoinWorkDouble maximumRHSError2 = 0.0;
  CoinWorkDouble primalOffset = 0.0;
  char *dropped = cholesky_->rowsDropped();
  for (iRow = 0; iRow < numberRows_; iRow++) {
    CoinWorkDouble value = errorRegion_[iRow];
    if (!dropped[iRow]) {
      if (CoinAbs(value) > maximumRHSError1) {
        maximumRHSError1 = CoinAbs(value);
      }
    } else {
      if (CoinAbs(value) > maximumRHSError2) {
        maximumRHSError2 = CoinAbs(value);
      }
      primalOffset += value * dualArray[iRow];
    }
  }
  primalObjective_ -= primalOffset * scaleFactor_;
  if (maximumRHSError1 > maximumRHSError2) {
    maximumRHSError_ = maximumRHSError1;
  } else {
    maximumRHSError_ = maximumRHSError1; //note change
    if (maximumRHSError2 > primalTolerance()) {
      handler_->message(CLP_BARRIER_ABS_DROPPED, messages_)
        << static_cast< double >(maximumRHSError2)
        << CoinMessageEol;
    }
  }
  objectiveNorm_ = maximumAbsElement(dualArray, numberRows_);
  if (objectiveNorm_ < 1.0e-12) {
    objectiveNorm_ = 1.0e-12;
  }
  if (objectiveNorm_ < baseObjectiveNorm_) {
    //std::cout<<" base "<<baseObjectiveNorm_<<" "<<objectiveNorm_<<std::endl;
    if (objectiveNorm_ < baseObjectiveNorm_ * 1.0e-4) {
      objectiveNorm_ = baseObjectiveNorm_ * 1.0e-4;
    }
  }
  bool primalFeasible = true;
  if (maximumRHSError_ > primalTolerance() || maximumDualError_ > dualTolerance / scaleFactor_) {
    handler_->message(CLP_BARRIER_ABS_ERROR, messages_)
      << static_cast< double >(maximumRHSError_) << static_cast< double >(maximumDualError_)
      << CoinMessageEol;
  }
  if (rhsNorm_ > solutionNorm_) {
    solutionNorm_ = rhsNorm_;
  }
  CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
  bool dualFeasible = true;
#if KEEP_GOING_IF_FIXED > 5
  if (maximumBoundInfeasibility_ > primalTolerance() || scaledRHSError > primalTolerance())
    primalFeasible = false;
#else
  if (maximumBoundInfeasibility_ > primalTolerance() || scaledRHSError > CoinMax(CoinMin(100.0 * primalTolerance(), 1.0e-5), primalTolerance()))
    primalFeasible = false;
#endif
  // relax dual test if obj big and gap smallish
  CoinWorkDouble gap = CoinAbs(primalObjective_ - dualObjective_);
  CoinWorkDouble sizeObj = CoinMin(CoinAbs(primalObjective_), CoinAbs(dualObjective_)) + 1.0e-50;
  //printf("gap %g sizeObj %g ratio %g comp %g\n",
  //     gap,sizeObj,gap/sizeObj,complementarityGap_);
  if (numberIterations_ > 100 && gap / sizeObj < 1.0e-9 && complementarityGap_ < 1.0e-7 * sizeObj)
    dualTolerance *= 1.0e2;
  if (maximumDualError_ > objectiveNorm_ * dualTolerance)
    dualFeasible = false;
  if (!primalFeasible || !dualFeasible) {
    handler_->message(CLP_BARRIER_FEASIBLE, messages_)
      << static_cast< double >(maximumBoundInfeasibility_) << static_cast< double >(scaledRHSError)
      << static_cast< double >(maximumDualError_ / objectiveNorm_)
      << CoinMessageEol;
  }
  if (!gonePrimalFeasible_) {
    gonePrimalFeasible_ = primalFeasible;
  } else if (!primalFeasible) {
    gonePrimalFeasible_ = primalFeasible;
    if (!numberKilled) {
      handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
        << CoinMessageEol;
    }
  }
  if (!goneDualFeasible_) {
    goneDualFeasible_ = dualFeasible;
  } else if (!dualFeasible) {
    handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
      << CoinMessageEol;
    goneDualFeasible_ = dualFeasible;
  }
  //objectiveValue();
  if (solutionNorm_ > 1.0e40) {
    std::cout << "primal off to infinity" << std::endl;
    abort();
  }
  if (objectiveNorm_ > 1.0e40) {
    std::cout << "dual off to infinity" << std::endl;
    abort();
  }
  handler_->message(CLP_BARRIER_STEP, messages_)
    << static_cast< double >(actualPrimalStep_)
    << static_cast< double >(actualDualStep_)
    << static_cast< double >(mu_)
    << CoinMessageEol;
  numberIterations_++;
  return numberKilled;
}
//  Save info on products of affine deltaSU*deltaW and deltaSL*deltaZ
CoinWorkDouble
ClpPredictorCorrector::affineProduct()
{
  CoinWorkDouble product = 0.0;
  //IF zVec starts as 0 then deltaZ always zero
  //(remember if free then zVec not 0)
  //I think free can be done with careful use of boundSlacks to zero
  //out all we want
  for (int iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    CoinWorkDouble w3 = deltaZ_[iColumn] * deltaX_[iColumn];
    CoinWorkDouble w4 = -deltaW_[iColumn] * deltaX_[iColumn];
    if (lowerBound(iColumn)) {
      w3 += deltaZ_[iColumn] * (solution_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn]);
      product += w3;
    }
    if (upperBound(iColumn)) {
      w4 += deltaW_[iColumn] * (-solution_[iColumn] - upperSlack_[iColumn] + upper_[iColumn]);
      product += w4;
    }
  }
  return product;
}
//See exactly what would happen given current deltas
void ClpPredictorCorrector::debugMove(int /*phase*/,
  CoinWorkDouble primalStep, CoinWorkDouble dualStep)
{
#ifndef SOME_DEBUG
  return;
#endif
  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
  int numberTotal = numberRows_ + numberColumns_;
  CoinWorkDouble *dualNew = ClpCopyOfArray(dualArray, numberRows_);
  CoinWorkDouble *errorRegionNew = new CoinWorkDouble[numberRows_];
  CoinWorkDouble *rhsFixRegionNew = new CoinWorkDouble[numberRows_];
  CoinWorkDouble *primalNew = ClpCopyOfArray(solution_, numberTotal);
  CoinWorkDouble *djNew = new CoinWorkDouble[numberTotal];
  //update pi
  multiplyAdd(deltaY_, numberRows_, dualStep, dualNew, 1.0);
  // do reduced costs
  CoinMemcpyN(dualNew, numberRows_, djNew + numberColumns_);
  CoinMemcpyN(cost_, numberColumns_, djNew);
  matrix_->transposeTimes(-1.0, dualNew, djNew);
  // update x
  int iColumn;
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn))
      primalNew[iColumn] += primalStep * deltaX_[iColumn];
  }
  CoinWorkDouble quadraticOffset = quadraticDjs(djNew, primalNew, 1.0);
  CoinZeroN(errorRegionNew, numberRows_);
  CoinZeroN(rhsFixRegionNew, numberRows_);
  CoinWorkDouble maximumBoundInfeasibility = 0.0;
  CoinWorkDouble maximumDualError = 1.0e-12;
  CoinWorkDouble primalObjectiveValue = 0.0;
  CoinWorkDouble dualObjectiveValue = 0.0;
  CoinWorkDouble solutionNorm = 1.0e-12;
  const CoinWorkDouble largeFactor = 1.0e2;
  CoinWorkDouble largeGap = largeFactor * solutionNorm_;
  if (largeGap < largeFactor) {
    largeGap = largeFactor;
  }
  CoinWorkDouble dualFake = 0.0;
  CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
  dualTolerance = dualTolerance / scaleFactor_;
  if (dualTolerance < 1.0e-12) {
    dualTolerance = 1.0e-12;
  }
  CoinWorkDouble newGap = 0.0;
  CoinWorkDouble offsetObjective = 0.0;
  CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
  CoinWorkDouble gammaOffset = 0.0;
  CoinWorkDouble maximumDjInfeasibility = 0.0;
  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    if (!flagged(iColumn)) {
      CoinWorkDouble reducedCost = djNew[iColumn];
      CoinWorkDouble zValue = zVec_[iColumn] + dualStep * deltaZ_[iColumn];
      CoinWorkDouble wValue = wVec_[iColumn] + dualStep * deltaW_[iColumn];
      CoinWorkDouble thisWeight = deltaX_[iColumn];
      CoinWorkDouble oldPrimal = solution_[iColumn];
      CoinWorkDouble newPrimal = primalNew[iColumn];
      CoinWorkDouble lowerBoundInfeasibility = 0.0;
      CoinWorkDouble upperBoundInfeasibility = 0.0;
      if (lowerBound(iColumn)) {
        CoinWorkDouble oldSlack = lowerSlack_[iColumn];
        CoinWorkDouble newSlack = lowerSlack_[iColumn] + primalStep * (oldPrimal - oldSlack + thisWeight - lower_[iColumn]);
        if (zValue > dualTolerance) {
          dualObjectiveValue += lower_[iColumn] * zVec_[iColumn];
        }
        lowerBoundInfeasibility = CoinAbs(newPrimal - newSlack - lower_[iColumn]);
        newGap += newSlack * zValue;
      }
      if (upperBound(iColumn)) {
        CoinWorkDouble oldSlack = upperSlack_[iColumn];
        CoinWorkDouble newSlack = upperSlack_[iColumn] + primalStep * (-oldPrimal - oldSlack - thisWeight + upper_[iColumn]);
        if (wValue > dualTolerance) {
          dualObjectiveValue -= upper_[iColumn] * wVec_[iColumn];
        }
        upperBoundInfeasibility = CoinAbs(newPrimal + newSlack - upper_[iColumn]);
        newGap += newSlack * wValue;
      }
      if (CoinAbs(newPrimal) > solutionNorm) {
        solutionNorm = CoinAbs(newPrimal);
      }
      CoinWorkDouble gammaTerm = gamma2;
      if (primalR_) {
        gammaTerm += primalR_[iColumn];
        quadraticOffset += newPrimal * newPrimal * primalR_[iColumn];
      }
      CoinWorkDouble dualInfeasibility = reducedCost - zValue + wValue + gammaTerm * newPrimal;
      if (CoinAbs(dualInfeasibility) > dualTolerance) {
        dualFake += newPrimal * dualInfeasibility;
      }
      if (lowerBoundInfeasibility > maximumBoundInfeasibility) {
        maximumBoundInfeasibility = lowerBoundInfeasibility;
      }
      if (upperBoundInfeasibility > maximumBoundInfeasibility) {
        maximumBoundInfeasibility = upperBoundInfeasibility;
      }
      dualInfeasibility = CoinAbs(dualInfeasibility);
      if (dualInfeasibility > maximumDualError) {
        //printf("bad dual %d %g\n",iColumn,
        // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
        maximumDualError = dualInfeasibility;
      }
      gammaOffset += newPrimal * newPrimal;
      djNew[iColumn] = 0.0;
    } else {
      offsetObjective += primalNew[iColumn] * cost_[iColumn];
      if (upper_[iColumn] - lower_[iColumn] > 1.0e-5) {
        if (primalNew[iColumn] < lower_[iColumn] + 1.0e-8 && djNew[iColumn] < -1.0e-8) {
          if (-djNew[iColumn] > maximumDjInfeasibility) {
            maximumDjInfeasibility = -djNew[iColumn];
          }
        }
        if (primalNew[iColumn] > upper_[iColumn] - 1.0e-8 && djNew[iColumn] > 1.0e-8) {
          if (djNew[iColumn] > maximumDjInfeasibility) {
            maximumDjInfeasibility = djNew[iColumn];
          }
        }
      }
      djNew[iColumn] = primalNew[iColumn];
    }
    primalObjectiveValue += solution_[iColumn] * cost_[iColumn];
  }
  // update rhs region
  multiplyAdd(djNew + numberColumns_, numberRows_, -1.0, rhsFixRegionNew, 1.0);
  matrix_->times(1.0, djNew, rhsFixRegionNew);
  primalObjectiveValue += 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
  dualObjectiveValue += offsetObjective + dualFake;
  dualObjectiveValue -= 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
  // Need to rethink (but it is only for printing)
  //compute error and fixed RHS
  multiplyAdd(primalNew + numberColumns_, numberRows_, -1.0, errorRegionNew, 0.0);
  matrix_->times(1.0, primalNew, errorRegionNew);
  //finish off objective computation
  CoinWorkDouble primalObjectiveNew = primalObjectiveValue * scaleFactor_;
  CoinWorkDouble dualValue2 = innerProduct(dualNew, numberRows_,
    rhsFixRegionNew);
  dualObjectiveValue -= dualValue2;
  //CoinWorkDouble dualObjectiveNew=dualObjectiveValue*scaleFactor_;
  CoinWorkDouble maximumRHSError1 = 0.0;
  CoinWorkDouble maximumRHSError2 = 0.0;
  CoinWorkDouble primalOffset = 0.0;
  char *dropped = cholesky_->rowsDropped();
  int iRow;
  for (iRow = 0; iRow < numberRows_; iRow++) {
    CoinWorkDouble value = errorRegionNew[iRow];
    if (!dropped[iRow]) {
      if (CoinAbs(value) > maximumRHSError1) {
        maximumRHSError1 = CoinAbs(value);
      }
    } else {
      if (CoinAbs(value) > maximumRHSError2) {
        maximumRHSError2 = CoinAbs(value);
      }
      primalOffset += value * dualNew[iRow];
    }
  }
  primalObjectiveNew -= primalOffset * scaleFactor_;
  //CoinWorkDouble maximumRHSError;
  if (maximumRHSError1 > maximumRHSError2) {
    ; //maximumRHSError = maximumRHSError1;
  } else {
    //maximumRHSError = maximumRHSError1; //note change
    if (maximumRHSError2 > primalTolerance()) {
      handler_->message(CLP_BARRIER_ABS_DROPPED, messages_)
        << static_cast< double >(maximumRHSError2)
        << CoinMessageEol;
    }
  }
  /*printf("PH %d %g, %g new comp %g, b %g, p %g, d %g\n",phase,
      primalStep,dualStep,newGap,maximumBoundInfeasibility,
      maximumRHSError,maximumDualError);
     if (handler_->logLevel()>1)
       printf("       objs %g %g\n",
        primalObjectiveNew,dualObjectiveNew);
     if (maximumDjInfeasibility) {
       printf(" max dj error on fixed %g\n",
        maximumDjInfeasibility);
        } */
  delete[] dualNew;
  delete[] errorRegionNew;
  delete[] rhsFixRegionNew;
  delete[] primalNew;
  delete[] djNew;
}

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