/* $Id: ClpPredictorCorrector.cpp 1959 2013-06-14 15:43:10Z stefan $ */
// 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;
}