Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (3 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpPredictorCorrector.cpp

    r1951 r2385  
    7474// main function
    7575
    76 int ClpPredictorCorrector::solve ( )
     76int ClpPredictorCorrector::solve()
    7777{
    78      problemStatus_ = -1;
    79      algorithm_ = 1;
    80      //create all regions
    81      if (!createWorkingData()) {
    82           problemStatus_ = 4;
    83           return 2;
    84      }
     78  problemStatus_ = -1;
     79  algorithm_ = 1;
     80  //create all regions
     81  if (!createWorkingData()) {
     82    problemStatus_ = 4;
     83    return 2;
     84  }
    8585#if COIN_LONG_WORK
    86      // reallocate some regions
    87      double * dualSave = dual_;
    88      dual_ = reinterpret_cast<double *>(new CoinWorkDouble[numberRows_]);
    89      double * reducedCostSave = reducedCost_;
    90      reducedCost_ = reinterpret_cast<double *>(new CoinWorkDouble[numberColumns_]);
    91 #endif
    92      //diagonalPerturbation_=1.0e-25;
    93      ClpMatrixBase * saveMatrix = NULL;
    94      // If quadratic then make copy so we can actually scale or normalize
     86  // reallocate some regions
     87  double *dualSave = dual_;
     88  dual_ = reinterpret_cast< double * >(new CoinWorkDouble[numberRows_]);
     89  double *reducedCostSave = reducedCost_;
     90  reducedCost_ = reinterpret_cast< double * >(new CoinWorkDouble[numberColumns_]);
     91#endif
     92  //diagonalPerturbation_=1.0e-25;
     93  ClpMatrixBase *saveMatrix = NULL;
     94  // If quadratic then make copy so we can actually scale or normalize
    9595#ifndef NO_RTTI
    96      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     96  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
    9797#else
    98      ClpQuadraticObjective * quadraticObj = NULL;
    99      if (objective_->type() == 2)
    100           quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    101 #endif
    102      /* If modeSwitch is :
     98  ClpQuadraticObjective *quadraticObj = NULL;
     99  if (objective_->type() == 2)
     100    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     101#endif
     102  /* If modeSwitch is :
    103103        0 - normal
    104104        1 - bit switch off centering
     
    106106        4 - accept corrector nearly always
    107107     */
    108      int modeSwitch = 0;
    109      //if (quadraticObj)
    110      //modeSwitch |= 1; // switch off centring for now
    111      //if (quadraticObj)
    112      //modeSwitch |=4;
    113      ClpObjective * saveObjective = NULL;
    114      if (quadraticObj) {
    115           // check KKT is on
    116           if (!cholesky_->kkt()) {
    117                //No!
    118                handler_->message(CLP_BARRIER_KKT, messages_)
    119                          << CoinMessageEol;
    120                return -1;
    121           }
    122           saveObjective = objective_;
    123           // We are going to make matrix full rather half
    124           objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
    125      }
    126      bool allowIncreasingGap = (modeSwitch & 4) != 0;
    127      // If scaled then really scale matrix
    128      if (scalingFlag_ > 0 && rowScale_) {
    129           saveMatrix = matrix_;
    130           matrix_ = matrix_->scaledColumnCopy(this);
    131      }
    132      //initializeFeasible(); - this just set fixed flag
    133      smallestInfeasibility_ = COIN_DBL_MAX;
    134      int i;
    135      for (i = 0; i < LENGTH_HISTORY; i++)
    136           historyInfeasibility_[i] = COIN_DBL_MAX;
     108  int modeSwitch = 0;
     109  //if (quadraticObj)
     110  //modeSwitch |= 1; // switch off centring for now
     111  //if (quadraticObj)
     112  //modeSwitch |=4;
     113  ClpObjective *saveObjective = NULL;
     114  if (quadraticObj) {
     115    // check KKT is on
     116    if (!cholesky_->kkt()) {
     117      //No!
     118      handler_->message(CLP_BARRIER_KKT, messages_)
     119        << CoinMessageEol;
     120      return -1;
     121    }
     122    saveObjective = objective_;
     123    // We are going to make matrix full rather half
     124    objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
     125  }
     126  bool allowIncreasingGap = (modeSwitch & 4) != 0;
     127  // If scaled then really scale matrix
     128  if (scalingFlag_ > 0 && rowScale_) {
     129    saveMatrix = matrix_;
     130    matrix_ = matrix_->scaledColumnCopy(this);
     131  }
     132  //initializeFeasible(); - this just set fixed flag
     133  smallestInfeasibility_ = COIN_DBL_MAX;
     134  int i;
     135  for (i = 0; i < LENGTH_HISTORY; i++)
     136    historyInfeasibility_[i] = COIN_DBL_MAX;
    137137
    138      //bool firstTime=true;
    139      //firstFactorization(true);
    140      int returnCode = cholesky_->order(this);
    141      if (returnCode || cholesky_->symbolic()) {
    142        COIN_DETAIL_PRINT(printf("Error return from symbolic - probably not enough memory\n"));
    143           problemStatus_ = 4;
    144           //delete all temporary regions
    145           deleteWorkingData();
    146           if (saveMatrix) {
    147                // restore normal copy
    148                delete matrix_;
    149                matrix_ = saveMatrix;
    150           }
    151           // Restore quadratic objective if necessary
    152           if (saveObjective) {
    153                delete objective_;
    154                objective_ = saveObjective;
    155           }
    156           return -1;
    157      }
    158      mu_ = 1.0e10;
    159      diagonalScaleFactor_ = 1.0;
    160      //set iterations
    161      numberIterations_ = -1;
    162      int numberTotal = numberRows_ + numberColumns_;
    163      //initialize solution here
    164      if(createSolution() < 0) {
    165        COIN_DETAIL_PRINT(printf("Not enough memory\n"));
    166           problemStatus_ = 4;
    167           //delete all temporary regions
    168           deleteWorkingData();
    169           if (saveMatrix) {
    170                // restore normal copy
    171                delete matrix_;
    172                matrix_ = saveMatrix;
    173           }
    174           return -1;
    175      }
    176      CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    177      // Could try centering steps without any original step i.e. just center
    178      //firstFactorization(false);
    179      CoinZeroN(dualArray, numberRows_);
    180      multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
    181      matrix_->times(1.0, solution_, errorRegion_);
    182      maximumRHSError_ = maximumAbsElement(errorRegion_, numberRows_);
    183      maximumBoundInfeasibility_ = maximumRHSError_;
    184      //CoinWorkDouble maximumDualError_=COIN_DBL_MAX;
    185      //initialize
    186      actualDualStep_ = 0.0;
    187      actualPrimalStep_ = 0.0;
    188      gonePrimalFeasible_ = false;
    189      goneDualFeasible_ = false;
    190      //bool hadGoodSolution=false;
    191      diagonalNorm_ = solutionNorm_;
    192      mu_ = solutionNorm_;
    193      int numberFixed = updateSolution(-COIN_DBL_MAX);
    194      int numberFixedTotal = numberFixed;
    195      //int numberRows_DroppedBefore=0;
    196      //CoinWorkDouble extra=eExtra;
    197      //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX;
    198      //constants for infeas interior point
    199      const CoinWorkDouble beta2 = 0.99995;
    200      const CoinWorkDouble tau   = 0.00002;
    201      CoinWorkDouble lastComplementarityGap = COIN_DBL_MAX * 1.0e-20;
    202      CoinWorkDouble lastStep = 1.0;
    203      // use to see if to take affine
    204      CoinWorkDouble checkGap = COIN_DBL_MAX;
    205      int lastGoodIteration = 0;
    206      CoinWorkDouble bestObjectiveGap = COIN_DBL_MAX;
    207      CoinWorkDouble bestObjective = COIN_DBL_MAX;
    208      int bestKilled = -1;
    209      int saveIteration = -1;
    210      int saveIteration2 = -1;
    211      bool sloppyOptimal = false;
    212      // this just to be used to exit
    213      bool sloppyOptimal2 = false;
    214      CoinWorkDouble * savePi = NULL;
    215      CoinWorkDouble * savePrimal = NULL;
    216      CoinWorkDouble * savePi2 = NULL;
    217      CoinWorkDouble * savePrimal2 = NULL;
    218      // Extra regions for centering
    219      CoinWorkDouble * saveX = new CoinWorkDouble[numberTotal];
    220      CoinWorkDouble * saveY = new CoinWorkDouble[numberRows_];
    221      CoinWorkDouble * saveZ = new CoinWorkDouble[numberTotal];
    222      CoinWorkDouble * saveW = new CoinWorkDouble[numberTotal];
    223      CoinWorkDouble * saveSL = new CoinWorkDouble[numberTotal];
    224      CoinWorkDouble * saveSU = new CoinWorkDouble[numberTotal];
    225      // Save smallest mu used in primal dual moves
    226      CoinWorkDouble objScale = optimizationDirection_ /
    227                                (rhsScale_ * objectiveScale_);
    228      while (problemStatus_ < 0) {
    229           //#define FULL_DEBUG
     138  //bool firstTime=true;
     139  //firstFactorization(true);
     140  int returnCode = cholesky_->order(this);
     141  if (returnCode || cholesky_->symbolic()) {
     142    COIN_DETAIL_PRINT(printf("Error return from symbolic - probably not enough memory\n"));
     143    problemStatus_ = 4;
     144    //delete all temporary regions
     145    deleteWorkingData();
     146    if (saveMatrix) {
     147      // restore normal copy
     148      delete matrix_;
     149      matrix_ = saveMatrix;
     150    }
     151    // Restore quadratic objective if necessary
     152    if (saveObjective) {
     153      delete objective_;
     154      objective_ = saveObjective;
     155    }
     156    return -1;
     157  }
     158  mu_ = 1.0e10;
     159  diagonalScaleFactor_ = 1.0;
     160  //set iterations
     161  numberIterations_ = -1;
     162  int numberTotal = numberRows_ + numberColumns_;
     163  //initialize solution here
     164  if (createSolution() < 0) {
     165    COIN_DETAIL_PRINT(printf("Not enough memory\n"));
     166    problemStatus_ = 4;
     167    //delete all temporary regions
     168    deleteWorkingData();
     169    if (saveMatrix) {
     170      // restore normal copy
     171      delete matrix_;
     172      matrix_ = saveMatrix;
     173    }
     174    return -1;
     175  }
     176  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
     177  // Could try centering steps without any original step i.e. just center
     178  //firstFactorization(false);
     179  CoinZeroN(dualArray, numberRows_);
     180  multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
     181  matrix_->times(1.0, solution_, errorRegion_);
     182  maximumRHSError_ = maximumAbsElement(errorRegion_, numberRows_);
     183  maximumBoundInfeasibility_ = maximumRHSError_;
     184  //CoinWorkDouble maximumDualError_=COIN_DBL_MAX;
     185  //initialize
     186  actualDualStep_ = 0.0;
     187  actualPrimalStep_ = 0.0;
     188  gonePrimalFeasible_ = false;
     189  goneDualFeasible_ = false;
     190  //bool hadGoodSolution=false;
     191  diagonalNorm_ = solutionNorm_;
     192  mu_ = solutionNorm_;
     193  int numberFixed = updateSolution(-COIN_DBL_MAX);
     194  int numberFixedTotal = numberFixed;
     195  //int numberRows_DroppedBefore=0;
     196  //CoinWorkDouble extra=eExtra;
     197  //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX;
     198  //constants for infeas interior point
     199  const CoinWorkDouble beta2 = 0.99995;
     200  const CoinWorkDouble tau = 0.00002;
     201  CoinWorkDouble lastComplementarityGap = COIN_DBL_MAX * 1.0e-20;
     202  CoinWorkDouble lastStep = 1.0;
     203  // use to see if to take affine
     204  CoinWorkDouble checkGap = COIN_DBL_MAX;
     205  int lastGoodIteration = 0;
     206  CoinWorkDouble bestObjectiveGap = COIN_DBL_MAX;
     207  CoinWorkDouble bestObjective = COIN_DBL_MAX;
     208  int bestKilled = -1;
     209  int saveIteration = -1;
     210  int saveIteration2 = -1;
     211  bool sloppyOptimal = false;
     212  // this just to be used to exit
     213  bool sloppyOptimal2 = false;
     214  CoinWorkDouble *savePi = NULL;
     215  CoinWorkDouble *savePrimal = NULL;
     216  CoinWorkDouble *savePi2 = NULL;
     217  CoinWorkDouble *savePrimal2 = NULL;
     218  // Extra regions for centering
     219  CoinWorkDouble *saveX = new CoinWorkDouble[numberTotal];
     220  CoinWorkDouble *saveY = new CoinWorkDouble[numberRows_];
     221  CoinWorkDouble *saveZ = new CoinWorkDouble[numberTotal];
     222  CoinWorkDouble *saveW = new CoinWorkDouble[numberTotal];
     223  CoinWorkDouble *saveSL = new CoinWorkDouble[numberTotal];
     224  CoinWorkDouble *saveSU = new CoinWorkDouble[numberTotal];
     225  // Save smallest mu used in primal dual moves
     226  CoinWorkDouble objScale = optimizationDirection_ / (rhsScale_ * objectiveScale_);
     227  while (problemStatus_ < 0) {
     228    //#define FULL_DEBUG
    230229#ifdef FULL_DEBUG
    231           {
    232                int i;
    233                printf("row    pi          artvec       rhsfx\n");
    234                for (i = 0; i < numberRows_; i++) {
    235                     printf("%d %g %g %g\n", i, dual_[i], errorRegion_[i], rhsFixRegion_[i]);
    236                }
    237                printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
    238                for (i = 0; i < numberColumns_ + numberRows_; i++) {
    239                     printf(" %d %g %g %g %g %g %g\n", i, solution_[i], diagonal_[i], wVec_[i],
    240                            zVec_[i], upperSlack_[i], lowerSlack_[i]);
    241                }
    242           }
    243 #endif
    244           complementarityGap_ = complementarityGap(numberComplementarityPairs_,
    245                                 numberComplementarityItems_, 0);
    246           handler_->message(CLP_BARRIER_ITERATION, messages_)
    247                     << numberIterations_
    248                     << static_cast<double>(primalObjective_ * objScale - dblParam_[ClpObjOffset])
    249                     << static_cast<double>(dualObjective_ * objScale - dblParam_[ClpObjOffset])
    250                     << static_cast<double>(complementarityGap_)
    251                     << numberFixedTotal
    252                     << cholesky_->rank()
    253                     << CoinMessageEol;
    254           // Check event
    255           {
    256             int status = eventHandler_->event(ClpEventHandler::endOfIteration);
    257             if (status >= 0) {
    258               problemStatus_ = 5;
    259               secondaryStatus_ = ClpEventHandler::endOfIteration;
    260               break;
    261             }
    262           }
     230    {
     231      int i;
     232      printf("row    pi          artvec       rhsfx\n");
     233      for (i = 0; i < numberRows_; i++) {
     234        printf("%d %g %g %g\n", i, dual_[i], errorRegion_[i], rhsFixRegion_[i]);
     235      }
     236      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
     237      for (i = 0; i < numberColumns_ + numberRows_; i++) {
     238        printf(" %d %g %g %g %g %g %g\n", i, solution_[i], diagonal_[i], wVec_[i],
     239          zVec_[i], upperSlack_[i], lowerSlack_[i]);
     240      }
     241    }
     242#endif
     243    complementarityGap_ = complementarityGap(numberComplementarityPairs_,
     244      numberComplementarityItems_, 0);
     245    handler_->message(CLP_BARRIER_ITERATION, messages_)
     246      << numberIterations_
     247      << static_cast< double >(primalObjective_ * objScale - dblParam_[ClpObjOffset])
     248      << static_cast< double >(dualObjective_ * objScale - dblParam_[ClpObjOffset])
     249      << static_cast< double >(complementarityGap_)
     250      << numberFixedTotal
     251      << cholesky_->rank()
     252      << CoinMessageEol;
     253    // Check event
     254    {
     255      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
     256      if (status >= 0) {
     257        problemStatus_ = 5;
     258        secondaryStatus_ = ClpEventHandler::endOfIteration;
     259        break;
     260      }
     261    }
    263262#if 0
    264263          if (numberIterations_ == -1) {
     
    268267          }
    269268#endif
    270           // move up history
    271           for (i = 1; i < LENGTH_HISTORY; i++)
    272                historyInfeasibility_[i-1] = historyInfeasibility_[i];
    273           historyInfeasibility_[LENGTH_HISTORY-1] = complementarityGap_;
    274           // switch off saved if changes
    275           //if (saveIteration+10<numberIterations_&&
    276           //complementarityGap_*2.0<historyInfeasibility_[0])
    277           //saveIteration=-1;
    278           lastStep = CoinMin(actualPrimalStep_, actualDualStep_);
    279           CoinWorkDouble goodGapChange;
    280           //#define KEEP_GOING_IF_FIXED 5
     269    // move up history
     270    for (i = 1; i < LENGTH_HISTORY; i++)
     271      historyInfeasibility_[i - 1] = historyInfeasibility_[i];
     272    historyInfeasibility_[LENGTH_HISTORY - 1] = complementarityGap_;
     273    // switch off saved if changes
     274    //if (saveIteration+10<numberIterations_&&
     275    //complementarityGap_*2.0<historyInfeasibility_[0])
     276    //saveIteration=-1;
     277    lastStep = CoinMin(actualPrimalStep_, actualDualStep_);
     278    CoinWorkDouble goodGapChange;
     279    //#define KEEP_GOING_IF_FIXED 5
    281280#ifndef KEEP_GOING_IF_FIXED
    282281#define KEEP_GOING_IF_FIXED 10000
    283282#endif
    284           if (!sloppyOptimal) {
    285                goodGapChange = 0.93;
    286           } else {
    287                goodGapChange = 0.7;
    288                if (numberFixed > KEEP_GOING_IF_FIXED)
    289                     goodGapChange = 0.99; // make more likely to carry on
    290           }
    291           CoinWorkDouble gapO;
    292           CoinWorkDouble lastGood = bestObjectiveGap;
    293           if (gonePrimalFeasible_ && goneDualFeasible_) {
    294                CoinWorkDouble largestObjective;
    295                if (CoinAbs(primalObjective_) > CoinAbs(dualObjective_)) {
    296                     largestObjective = CoinAbs(primalObjective_);
    297                } else {
    298                     largestObjective = CoinAbs(dualObjective_);
    299                }
    300                if (largestObjective < 1.0) {
    301                     largestObjective = 1.0;
    302                }
    303                gapO = CoinAbs(primalObjective_ - dualObjective_) / largestObjective;
    304                handler_->message(CLP_BARRIER_OBJECTIVE_GAP, messages_)
    305                          << static_cast<double>(gapO)
    306                          << CoinMessageEol;
    307                //start saving best
    308                bool saveIt = false;
    309                if (gapO < bestObjectiveGap) {
    310                     bestObjectiveGap = gapO;
     283    if (!sloppyOptimal) {
     284      goodGapChange = 0.93;
     285    } else {
     286      goodGapChange = 0.7;
     287      if (numberFixed > KEEP_GOING_IF_FIXED)
     288        goodGapChange = 0.99; // make more likely to carry on
     289    }
     290    CoinWorkDouble gapO;
     291    CoinWorkDouble lastGood = bestObjectiveGap;
     292    if (gonePrimalFeasible_ && goneDualFeasible_) {
     293      CoinWorkDouble largestObjective;
     294      if (CoinAbs(primalObjective_) > CoinAbs(dualObjective_)) {
     295        largestObjective = CoinAbs(primalObjective_);
     296      } else {
     297        largestObjective = CoinAbs(dualObjective_);
     298      }
     299      if (largestObjective < 1.0) {
     300        largestObjective = 1.0;
     301      }
     302      gapO = CoinAbs(primalObjective_ - dualObjective_) / largestObjective;
     303      handler_->message(CLP_BARRIER_OBJECTIVE_GAP, messages_)
     304        << static_cast< double >(gapO)
     305        << CoinMessageEol;
     306      //start saving best
     307      bool saveIt = false;
     308      if (gapO < bestObjectiveGap) {
     309        bestObjectiveGap = gapO;
    311310#ifndef SAVE_ON_OBJ
    312                     saveIt = true;
    313 #endif
    314                }
    315                if (primalObjective_ < bestObjective) {
    316                     bestObjective = primalObjective_;
     311        saveIt = true;
     312#endif
     313      }
     314      if (primalObjective_ < bestObjective) {
     315        bestObjective = primalObjective_;
    317316#ifdef SAVE_ON_OBJ
    318                     saveIt = true;
    319 #endif
    320                }
    321                if (numberFixedTotal > bestKilled) {
    322                     bestKilled = numberFixedTotal;
    323 #if KEEP_GOING_IF_FIXED<10
    324                     saveIt = true;
    325 #endif
    326                }
    327                if (saveIt) {
    328 #if KEEP_GOING_IF_FIXED<10
    329                  COIN_DETAIL_PRINT(printf("saving\n"));
    330 #endif
    331                     saveIteration = numberIterations_;
    332                     if (!savePi) {
    333                          savePi = new CoinWorkDouble[numberRows_];
    334                          savePrimal = new CoinWorkDouble [numberTotal];
    335                     }
    336                     CoinMemcpyN(dualArray, numberRows_, savePi);
    337                     CoinMemcpyN(solution_, numberTotal, savePrimal);
    338                } else if(gapO > 2.0 * bestObjectiveGap) {
    339                     //maybe be more sophisticated e.g. re-initialize having
    340                     //fixed variables and dropped rows
    341                     //std::cout <<" gap increasing "<<std::endl;
    342                }
    343                //std::cout <<"could stop"<<std::endl;
    344                //gapO=0.0;
    345                if (CoinAbs(primalObjective_ - dualObjective_) < dualTolerance()) {
    346                     gapO = 0.0;
    347                }
    348           } else {
    349                gapO = COIN_DBL_MAX;
    350                if (saveIteration >= 0) {
    351                     handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
    352                               << CoinMessageEol;
    353                     CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
    354                     // save alternate
    355                     if (numberFixedTotal > bestKilled &&
    356                               maximumBoundInfeasibility_ < 1.0e-6 &&
    357                               scaledRHSError < 1.0e-2) {
    358                          bestKilled = numberFixedTotal;
    359 #if KEEP_GOING_IF_FIXED<10
    360                          COIN_DETAIL_PRINT(printf("saving alternate\n"));
    361 #endif
    362                          saveIteration2 = numberIterations_;
    363                          if (!savePi2) {
    364                               savePi2 = new CoinWorkDouble[numberRows_];
    365                               savePrimal2 = new CoinWorkDouble [numberTotal];
    366                          }
    367                          CoinMemcpyN(dualArray, numberRows_, savePi2);
    368                          CoinMemcpyN(solution_, numberTotal, savePrimal2);
    369                     }
    370                     if (sloppyOptimal2) {
    371                          // vaguely optimal
    372                          if (maximumBoundInfeasibility_ > 1.0e-2 ||
    373                                    scaledRHSError > 1.0e-2 ||
    374                                    maximumDualError_ > objectiveNorm_ * 1.0e-2) {
    375                               handler_->message(CLP_BARRIER_EXIT2, messages_)
    376                                         << saveIteration
    377                                         << CoinMessageEol;
    378                               problemStatus_ = 0; // benefit of doubt
    379                               break;
    380                          }
    381                     } else {
    382                          // not close to optimal but check if getting bad
    383                          CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
    384                          if ((maximumBoundInfeasibility_ > 1.0e-1 ||
    385                                    scaledRHSError > 1.0e-1 ||
    386                                    maximumDualError_ > objectiveNorm_ * 1.0e-1)
    387                                    && (numberIterations_ > 50
    388                                        && complementarityGap_ > 0.9 * historyInfeasibility_[0])) {
    389                               handler_->message(CLP_BARRIER_EXIT2, messages_)
    390                                         << saveIteration
    391                                         << CoinMessageEol;
    392                               break;
    393                          }
    394                          if (complementarityGap_ > 0.95 * checkGap && bestObjectiveGap < 1.0e-3 &&
    395                                    (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
    396                               handler_->message(CLP_BARRIER_EXIT2, messages_)
    397                                         << saveIteration
    398                                         << CoinMessageEol;
    399                               break;
    400                          }
    401                     }
    402                }
    403                if (complementarityGap_ > 0.5 * checkGap && primalObjective_ >
    404                          bestObjective + 1.0e-9 &&
    405                          (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
    406                     handler_->message(CLP_BARRIER_EXIT2, messages_)
    407                               << saveIteration
    408                               << CoinMessageEol;
    409                     break;
    410                }
    411           }
    412           // See if we should be thinking about exit if diverging
    413           double relativeMultiplier = 1.0 + fabs(primalObjective_) + fabs(dualObjective_);
    414           // Quadratic coding is rubbish so be more forgiving?
    415           if (quadraticObj)
    416             relativeMultiplier *= 5.0;
    417           if (gapO < 1.0e-5 + 1.0e-9 * relativeMultiplier
    418               || complementarityGap_ < 0.1 + 1.0e-9 * relativeMultiplier)
    419               sloppyOptimal2 = true;
    420           if ((gapO < 1.0e-6 || (gapO < 1.0e-4 && complementarityGap_ < 0.1)) && !sloppyOptimal) {
    421                sloppyOptimal = true;
    422                sloppyOptimal2 = true;
    423                handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL, messages_)
    424                          << numberIterations_ << static_cast<double>(complementarityGap_)
    425                          << CoinMessageEol;
    426           }
    427           int numberBack = quadraticObj ? 10 : 5;
    428           //tryJustPredictor=true;
    429           //printf("trying just predictor\n");
    430           //}
    431           if (complementarityGap_ >= 1.05 * lastComplementarityGap) {
    432                handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
    433                          << static_cast<double>(complementarityGap_) << "increasing"
    434                          << CoinMessageEol;
    435                if (saveIteration >= 0 && sloppyOptimal2) {
    436                     handler_->message(CLP_BARRIER_EXIT2, messages_)
    437                               << saveIteration
    438                               << CoinMessageEol;
    439                     break;
    440                } else if (numberIterations_ - lastGoodIteration >= numberBack &&
    441                           complementarityGap_ < 1.0e-6) {
    442                     break; // not doing very well - give up
    443                }
    444           } else if (complementarityGap_ < goodGapChange * lastComplementarityGap) {
    445                lastGoodIteration = numberIterations_;
    446                lastComplementarityGap = complementarityGap_;
    447           } else if (numberIterations_ - lastGoodIteration >= numberBack &&
    448                      complementarityGap_ < 1.0e-3) {
    449                handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
    450                          << static_cast<double>(complementarityGap_) << "not decreasing"
    451                          << CoinMessageEol;
    452                if (gapO > 0.75 * lastGood && numberFixed < KEEP_GOING_IF_FIXED) {
    453                     break;
    454                }
    455           } else if (numberIterations_ - lastGoodIteration >= 2 &&
    456                      complementarityGap_ < 1.0e-6) {
    457                handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
    458                          << static_cast<double>(complementarityGap_) << "not decreasing"
    459                          << CoinMessageEol;
    460                break;
    461           }
    462           if (numberIterations_ > maximumBarrierIterations_ || hitMaximumIterations()) {
    463                handler_->message(CLP_BARRIER_STOPPING, messages_)
    464                          << CoinMessageEol;
    465                problemStatus_ = 3;
    466                onStopped(); // set secondary status
    467                break;
    468           }
    469           if (gapO < targetGap_) {
    470                problemStatus_ = 0;
    471                handler_->message(CLP_BARRIER_EXIT, messages_)
    472                          << " "
    473                          << CoinMessageEol;
    474                break;//finished
    475           }
    476           if (complementarityGap_ < 1.0e-12) {
    477                problemStatus_ = 0;
    478                handler_->message(CLP_BARRIER_EXIT, messages_)
    479                          << "- small complementarity gap"
    480                          << CoinMessageEol;
    481                break;//finished
    482           }
    483           if (complementarityGap_ < 1.0e-10 && gapO < 1.0e-10) {
    484                problemStatus_ = 0;
    485                handler_->message(CLP_BARRIER_EXIT, messages_)
    486                          << "- objective gap and complementarity gap both small"
    487                          << CoinMessageEol;
    488                break;//finished
    489           }
    490           if (gapO < 1.0e-9) {
    491                CoinWorkDouble value = gapO * complementarityGap_;
    492                value *= actualPrimalStep_;
    493                value *= actualDualStep_;
    494                //std::cout<<value<<std::endl;
    495                if (value < 1.0e-17 && numberIterations_ > lastGoodIteration) {
    496                     problemStatus_ = 0;
    497                     handler_->message(CLP_BARRIER_EXIT, messages_)
    498                               << "- objective gap and complementarity gap both smallish and small steps"
    499                               << CoinMessageEol;
    500                     break;//finished
    501                }
    502           }
    503           CoinWorkDouble nextGap = COIN_DBL_MAX;
    504           int nextNumber = 0;
    505           int nextNumberItems = 0;
    506           worstDirectionAccuracy_ = 0.0;
    507           int newDropped = 0;
    508           //Predictor step
    509           //prepare for cholesky.  Set up scaled diagonal in deltaX
    510           //  ** for efficiency may be better if scale factor known before
    511           CoinWorkDouble norm2 = 0.0;
    512           CoinWorkDouble maximumValue;
    513           getNorms(diagonal_, numberTotal, maximumValue, norm2);
    514           diagonalNorm_ = CoinSqrt(norm2 / numberComplementarityPairs_);
    515           diagonalScaleFactor_ = 1.0;
    516           CoinWorkDouble maximumAllowable = eScale;
    517           //scale so largest is less than allowable ? could do better
    518           CoinWorkDouble factor = 0.5;
    519           while (maximumValue > maximumAllowable) {
    520                diagonalScaleFactor_ *= factor;
    521                maximumValue *= factor;
    522           } /* endwhile */
    523           if (diagonalScaleFactor_ != 1.0) {
    524                handler_->message(CLP_BARRIER_SCALING, messages_)
    525                          << "diagonal" << static_cast<double>(diagonalScaleFactor_)
    526                          << CoinMessageEol;
    527                diagonalNorm_ *= diagonalScaleFactor_;
    528           }
    529           multiplyAdd(NULL, numberTotal, 0.0, diagonal_,
    530                       diagonalScaleFactor_);
    531           int * rowsDroppedThisTime = new int [numberRows_];
    532           newDropped = cholesky_->factorize(diagonal_, rowsDroppedThisTime);
    533           if (newDropped) {
    534                if (newDropped == -1) {
    535                  COIN_DETAIL_PRINT(printf("Out of memory\n"));
    536                     problemStatus_ = 4;
    537                     //delete all temporary regions
    538                     deleteWorkingData();
    539                     if (saveMatrix) {
    540                          // restore normal copy
    541                          delete matrix_;
    542                          matrix_ = saveMatrix;
    543                     }
    544                     return -1;
    545                } else {
     317        saveIt = true;
     318#endif
     319      }
     320      if (numberFixedTotal > bestKilled) {
     321        bestKilled = numberFixedTotal;
     322#if KEEP_GOING_IF_FIXED < 10
     323        saveIt = true;
     324#endif
     325      }
     326      if (saveIt) {
     327#if KEEP_GOING_IF_FIXED < 10
     328        COIN_DETAIL_PRINT(printf("saving\n"));
     329#endif
     330        saveIteration = numberIterations_;
     331        if (!savePi) {
     332          savePi = new CoinWorkDouble[numberRows_];
     333          savePrimal = new CoinWorkDouble[numberTotal];
     334        }
     335        CoinMemcpyN(dualArray, numberRows_, savePi);
     336        CoinMemcpyN(solution_, numberTotal, savePrimal);
     337      } else if (gapO > 2.0 * bestObjectiveGap) {
     338        //maybe be more sophisticated e.g. re-initialize having
     339        //fixed variables and dropped rows
     340        //std::cout <<" gap increasing "<<std::endl;
     341      }
     342      //std::cout <<"could stop"<<std::endl;
     343      //gapO=0.0;
     344      if (CoinAbs(primalObjective_ - dualObjective_) < dualTolerance()) {
     345        gapO = 0.0;
     346      }
     347    } else {
     348      gapO = COIN_DBL_MAX;
     349      if (saveIteration >= 0) {
     350        handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
     351          << CoinMessageEol;
     352        CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
     353        // save alternate
     354        if (numberFixedTotal > bestKilled && maximumBoundInfeasibility_ < 1.0e-6 && scaledRHSError < 1.0e-2) {
     355          bestKilled = numberFixedTotal;
     356#if KEEP_GOING_IF_FIXED < 10
     357          COIN_DETAIL_PRINT(printf("saving alternate\n"));
     358#endif
     359          saveIteration2 = numberIterations_;
     360          if (!savePi2) {
     361            savePi2 = new CoinWorkDouble[numberRows_];
     362            savePrimal2 = new CoinWorkDouble[numberTotal];
     363          }
     364          CoinMemcpyN(dualArray, numberRows_, savePi2);
     365          CoinMemcpyN(solution_, numberTotal, savePrimal2);
     366        }
     367        if (sloppyOptimal2) {
     368          // vaguely optimal
     369          if (maximumBoundInfeasibility_ > 1.0e-2 || scaledRHSError > 1.0e-2 || maximumDualError_ > objectiveNorm_ * 1.0e-2) {
     370            handler_->message(CLP_BARRIER_EXIT2, messages_)
     371              << saveIteration
     372              << CoinMessageEol;
     373            problemStatus_ = 0; // benefit of doubt
     374            break;
     375          }
     376        } else {
     377          // not close to optimal but check if getting bad
     378          CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
     379          if ((maximumBoundInfeasibility_ > 1.0e-1 || scaledRHSError > 1.0e-1 || maximumDualError_ > objectiveNorm_ * 1.0e-1)
     380            && (numberIterations_ > 50
     381                 && complementarityGap_ > 0.9 * historyInfeasibility_[0])) {
     382            handler_->message(CLP_BARRIER_EXIT2, messages_)
     383              << saveIteration
     384              << CoinMessageEol;
     385            break;
     386          }
     387          if (complementarityGap_ > 0.95 * checkGap && bestObjectiveGap < 1.0e-3 && (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
     388            handler_->message(CLP_BARRIER_EXIT2, messages_)
     389              << saveIteration
     390              << CoinMessageEol;
     391            break;
     392          }
     393        }
     394      }
     395      if (complementarityGap_ > 0.5 * checkGap && primalObjective_ > bestObjective + 1.0e-9 && (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
     396        handler_->message(CLP_BARRIER_EXIT2, messages_)
     397          << saveIteration
     398          << CoinMessageEol;
     399        break;
     400      }
     401    }
     402    // See if we should be thinking about exit if diverging
     403    double relativeMultiplier = 1.0 + fabs(primalObjective_) + fabs(dualObjective_);
     404    // Quadratic coding is rubbish so be more forgiving?
     405    if (quadraticObj)
     406      relativeMultiplier *= 5.0;
     407    if (gapO < 1.0e-5 + 1.0e-9 * relativeMultiplier
     408      || complementarityGap_ < 0.1 + 1.0e-9 * relativeMultiplier)
     409      sloppyOptimal2 = true;
     410    if ((gapO < 1.0e-6 || (gapO < 1.0e-4 && complementarityGap_ < 0.1)) && !sloppyOptimal) {
     411      sloppyOptimal = true;
     412      sloppyOptimal2 = true;
     413      handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL, messages_)
     414        << numberIterations_ << static_cast< double >(complementarityGap_)
     415        << CoinMessageEol;
     416    }
     417    int numberBack = quadraticObj ? 10 : 5;
     418    //tryJustPredictor=true;
     419    //printf("trying just predictor\n");
     420    //}
     421    if (complementarityGap_ >= 1.05 * lastComplementarityGap) {
     422      handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
     423        << static_cast< double >(complementarityGap_) << "increasing"
     424        << CoinMessageEol;
     425      if (saveIteration >= 0 && sloppyOptimal2) {
     426        handler_->message(CLP_BARRIER_EXIT2, messages_)
     427          << saveIteration
     428          << CoinMessageEol;
     429        break;
     430      } else if (numberIterations_ - lastGoodIteration >= numberBack && complementarityGap_ < 1.0e-6) {
     431        break; // not doing very well - give up
     432      }
     433    } else if (complementarityGap_ < goodGapChange * lastComplementarityGap) {
     434      lastGoodIteration = numberIterations_;
     435      lastComplementarityGap = complementarityGap_;
     436    } else if (numberIterations_ - lastGoodIteration >= numberBack && complementarityGap_ < 1.0e-3) {
     437      handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
     438        << static_cast< double >(complementarityGap_) << "not decreasing"
     439        << CoinMessageEol;
     440      if (gapO > 0.75 * lastGood && numberFixed < KEEP_GOING_IF_FIXED) {
     441        break;
     442      }
     443    } else if (numberIterations_ - lastGoodIteration >= 2 && complementarityGap_ < 1.0e-6) {
     444      handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
     445        << static_cast< double >(complementarityGap_) << "not decreasing"
     446        << CoinMessageEol;
     447      break;
     448    }
     449    if (numberIterations_ > maximumBarrierIterations_ || hitMaximumIterations()) {
     450      handler_->message(CLP_BARRIER_STOPPING, messages_)
     451        << CoinMessageEol;
     452      problemStatus_ = 3;
     453      onStopped(); // set secondary status
     454      break;
     455    }
     456    if (gapO < targetGap_) {
     457      problemStatus_ = 0;
     458      handler_->message(CLP_BARRIER_EXIT, messages_)
     459        << " "
     460        << CoinMessageEol;
     461      break; //finished
     462    }
     463    if (complementarityGap_ < 1.0e-12) {
     464      problemStatus_ = 0;
     465      handler_->message(CLP_BARRIER_EXIT, messages_)
     466        << "- small complementarity gap"
     467        << CoinMessageEol;
     468      break; //finished
     469    }
     470    if (complementarityGap_ < 1.0e-10 && gapO < 1.0e-10) {
     471      problemStatus_ = 0;
     472      handler_->message(CLP_BARRIER_EXIT, messages_)
     473        << "- objective gap and complementarity gap both small"
     474        << CoinMessageEol;
     475      break; //finished
     476    }
     477    if (gapO < 1.0e-9) {
     478      CoinWorkDouble value = gapO * complementarityGap_;
     479      value *= actualPrimalStep_;
     480      value *= actualDualStep_;
     481      //std::cout<<value<<std::endl;
     482      if (value < 1.0e-17 && numberIterations_ > lastGoodIteration) {
     483        problemStatus_ = 0;
     484        handler_->message(CLP_BARRIER_EXIT, messages_)
     485          << "- objective gap and complementarity gap both smallish and small steps"
     486          << CoinMessageEol;
     487        break; //finished
     488      }
     489    }
     490    CoinWorkDouble nextGap = COIN_DBL_MAX;
     491    int nextNumber = 0;
     492    int nextNumberItems = 0;
     493    worstDirectionAccuracy_ = 0.0;
     494    int newDropped = 0;
     495    //Predictor step
     496    //prepare for cholesky.  Set up scaled diagonal in deltaX
     497    //  ** for efficiency may be better if scale factor known before
     498    CoinWorkDouble norm2 = 0.0;
     499    CoinWorkDouble maximumValue;
     500    getNorms(diagonal_, numberTotal, maximumValue, norm2);
     501    diagonalNorm_ = CoinSqrt(norm2 / numberComplementarityPairs_);
     502    diagonalScaleFactor_ = 1.0;
     503    CoinWorkDouble maximumAllowable = eScale;
     504    //scale so largest is less than allowable ? could do better
     505    CoinWorkDouble factor = 0.5;
     506    while (maximumValue > maximumAllowable) {
     507      diagonalScaleFactor_ *= factor;
     508      maximumValue *= factor;
     509    } /* endwhile */
     510    if (diagonalScaleFactor_ != 1.0) {
     511      handler_->message(CLP_BARRIER_SCALING, messages_)
     512        << "diagonal" << static_cast< double >(diagonalScaleFactor_)
     513        << CoinMessageEol;
     514      diagonalNorm_ *= diagonalScaleFactor_;
     515    }
     516    multiplyAdd(NULL, numberTotal, 0.0, diagonal_,
     517      diagonalScaleFactor_);
     518    int *rowsDroppedThisTime = new int[numberRows_];
     519    newDropped = cholesky_->factorize(diagonal_, rowsDroppedThisTime);
     520    if (newDropped) {
     521      if (newDropped == -1) {
     522        COIN_DETAIL_PRINT(printf("Out of memory\n"));
     523        problemStatus_ = 4;
     524        //delete all temporary regions
     525        deleteWorkingData();
     526        if (saveMatrix) {
     527          // restore normal copy
     528          delete matrix_;
     529          matrix_ = saveMatrix;
     530        }
     531        return -1;
     532      } else {
    546533#ifndef NDEBUG
    547                     //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
    548                     //assert(!newDropped2);
    549 #endif
    550                     if (newDropped < 0 && 0) {
    551                          //replace dropped
    552                          newDropped = -newDropped;
    553                          //off 1 to allow for reset all
    554                          newDropped--;
    555                          //set all bits false
    556                          cholesky_->resetRowsDropped();
    557                     }
    558                }
    559           }
    560           delete [] rowsDroppedThisTime;
    561           if (cholesky_->status()) {
    562                std::cout << "bad cholesky?" << std::endl;
    563                abort();
    564           }
    565           int phase = 0; // predictor, corrector , primal dual
    566           CoinWorkDouble directionAccuracy = 0.0;
    567           bool doCorrector = true;
    568           bool goodMove = true;
    569           //set up for affine direction
    570           setupForSolve(phase);
    571           if ((modeSwitch & 2) == 0) {
    572                directionAccuracy = findDirectionVector(phase);
    573                if (directionAccuracy > worstDirectionAccuracy_) {
    574                     worstDirectionAccuracy_ = directionAccuracy;
    575                }
    576                if (saveIteration > 0 && directionAccuracy > 1.0) {
    577                     handler_->message(CLP_BARRIER_EXIT2, messages_)
    578                               << saveIteration
    579                               << CoinMessageEol;
    580                     break;
    581                }
    582                findStepLength(phase);
    583                nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
    584                debugMove(0, actualPrimalStep_, actualDualStep_);
    585                debugMove(0, 1.0e-2, 1.0e-2);
    586           }
    587           CoinWorkDouble affineGap = nextGap;
    588           int bestPhase = 0;
    589           CoinWorkDouble bestNextGap = nextGap;
    590           // ?
    591           bestNextGap = CoinMax(nextGap, 0.8 * complementarityGap_);
    592           if (quadraticObj)
    593                bestNextGap = CoinMax(nextGap, 0.99 * complementarityGap_);
    594           if (complementarityGap_ > 1.0e-4 * numberComplementarityPairs_) {
    595                //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
    596                CoinWorkDouble part1 = nextGap / numberComplementarityPairs_;
    597                part1 = nextGap / numberComplementarityItems_;
    598                CoinWorkDouble part2 = nextGap / complementarityGap_;
    599                mu_ = part1 * part2 * part2;
     534        //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
     535        //assert(!newDropped2);
     536#endif
     537        if (newDropped < 0 && 0) {
     538          //replace dropped
     539          newDropped = -newDropped;
     540          //off 1 to allow for reset all
     541          newDropped--;
     542          //set all bits false
     543          cholesky_->resetRowsDropped();
     544        }
     545      }
     546    }
     547    delete[] rowsDroppedThisTime;
     548    if (cholesky_->status()) {
     549      std::cout << "bad cholesky?" << std::endl;
     550      abort();
     551    }
     552    int phase = 0; // predictor, corrector , primal dual
     553    CoinWorkDouble directionAccuracy = 0.0;
     554    bool doCorrector = true;
     555    bool goodMove = true;
     556    //set up for affine direction
     557    setupForSolve(phase);
     558    if ((modeSwitch & 2) == 0) {
     559      directionAccuracy = findDirectionVector(phase);
     560      if (directionAccuracy > worstDirectionAccuracy_) {
     561        worstDirectionAccuracy_ = directionAccuracy;
     562      }
     563      if (saveIteration > 0 && directionAccuracy > 1.0) {
     564        handler_->message(CLP_BARRIER_EXIT2, messages_)
     565          << saveIteration
     566          << CoinMessageEol;
     567        break;
     568      }
     569      findStepLength(phase);
     570      nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
     571      debugMove(0, actualPrimalStep_, actualDualStep_);
     572      debugMove(0, 1.0e-2, 1.0e-2);
     573    }
     574    CoinWorkDouble affineGap = nextGap;
     575    int bestPhase = 0;
     576    CoinWorkDouble bestNextGap = nextGap;
     577    // ?
     578    bestNextGap = CoinMax(nextGap, 0.8 * complementarityGap_);
     579    if (quadraticObj)
     580      bestNextGap = CoinMax(nextGap, 0.99 * complementarityGap_);
     581    if (complementarityGap_ > 1.0e-4 * numberComplementarityPairs_) {
     582      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
     583      CoinWorkDouble part1 = nextGap / numberComplementarityPairs_;
     584      part1 = nextGap / numberComplementarityItems_;
     585      CoinWorkDouble part2 = nextGap / complementarityGap_;
     586      mu_ = part1 * part2 * part2;
    600587#if 0
    601588               CoinWorkDouble papermu = complementarityGap_ / numberComplementarityPairs_;
     
    605592                      mu_, papermu, affmu, sigma, sigma * papermu);
    606593#endif
    607                //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
    608                //                                           (CoinWorkDouble) numberComplementarityPairs_));
    609           } else {
    610                CoinWorkDouble phi;
    611                if (numberComplementarityPairs_ <= 5000) {
    612                     phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 2.0);
    613                } else {
    614                     phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 1.5);
    615                     if (phi < 500.0 * 500.0) {
    616                          phi = 500.0 * 500.0;
    617                     }
    618                }
    619                mu_ = complementarityGap_ / phi;
    620           }
    621           //save information
    622           CoinWorkDouble product = affineProduct();
     594      //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
     595      //                                            (CoinWorkDouble) numberComplementarityPairs_));
     596    } else {
     597      CoinWorkDouble phi;
     598      if (numberComplementarityPairs_ <= 5000) {
     599        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 2.0);
     600      } else {
     601        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 1.5);
     602        if (phi < 500.0 * 500.0) {
     603          phi = 500.0 * 500.0;
     604        }
     605      }
     606      mu_ = complementarityGap_ / phi;
     607    }
     608    //save information
     609    CoinWorkDouble product = affineProduct();
    623610#if 0
    624611          //can we do corrector step?
     
    650637          }
    651638#endif
    652           if (complementarityGap_*(beta2 - tau) + product - mu_ * numberComplementarityPairs_ < 0.0 && 0) {
     639    if (complementarityGap_ * (beta2 - tau) + product - mu_ * numberComplementarityPairs_ < 0.0 && 0) {
    653640#ifdef SOME_DEBUG
    654                printf("failed 1 product %.18g mu %.18g - %.18g < 0.0, nextGap %.18g\n", product, mu_,
    655                       complementarityGap_*(beta2 - tau) + product - mu_ * numberComplementarityPairs_,
    656                       nextGap);
    657 #endif
    658                doCorrector = false;
    659                if (nextGap > 0.9 * complementarityGap_ || 1) {
    660                     goodMove = false;
    661                     bestNextGap = COIN_DBL_MAX;
    662                }
    663                //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_;
    664                //floatNumber = 1.0*numberComplementarityItems_;
    665                //mu_=nextGap/floatNumber;
    666                handler_->message(CLP_BARRIER_INFO, messages_)
    667                          << "no corrector step"
    668                          << CoinMessageEol;
    669           } else {
    670                phase = 1;
    671           }
    672           // If bad gap - try standard primal dual
    673           if (nextGap > complementarityGap_ * 1.001)
    674                goodMove = false;
    675           if ((modeSwitch & 2) != 0)
    676                goodMove = false;
    677           if (goodMove && doCorrector) {
    678                CoinMemcpyN(deltaX_, numberTotal, saveX);
    679                CoinMemcpyN(deltaY_, numberRows_, saveY);
    680                CoinMemcpyN(deltaZ_, numberTotal, saveZ);
    681                CoinMemcpyN(deltaW_, numberTotal, saveW);
    682                CoinMemcpyN(deltaSL_, numberTotal, saveSL);
    683                CoinMemcpyN(deltaSU_, numberTotal, saveSU);
     641      printf("failed 1 product %.18g mu %.18g - %.18g < 0.0, nextGap %.18g\n", product, mu_,
     642        complementarityGap_ * (beta2 - tau) + product - mu_ * numberComplementarityPairs_,
     643        nextGap);
     644#endif
     645      doCorrector = false;
     646      if (nextGap > 0.9 * complementarityGap_ || 1) {
     647        goodMove = false;
     648        bestNextGap = COIN_DBL_MAX;
     649      }
     650      //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_;
     651      //floatNumber = 1.0*numberComplementarityItems_;
     652      //mu_=nextGap/floatNumber;
     653      handler_->message(CLP_BARRIER_INFO, messages_)
     654        << "no corrector step"
     655        << CoinMessageEol;
     656    } else {
     657      phase = 1;
     658    }
     659    // If bad gap - try standard primal dual
     660    if (nextGap > complementarityGap_ * 1.001)
     661      goodMove = false;
     662    if ((modeSwitch & 2) != 0)
     663      goodMove = false;
     664    if (goodMove && doCorrector) {
     665      CoinMemcpyN(deltaX_, numberTotal, saveX);
     666      CoinMemcpyN(deltaY_, numberRows_, saveY);
     667      CoinMemcpyN(deltaZ_, numberTotal, saveZ);
     668      CoinMemcpyN(deltaW_, numberTotal, saveW);
     669      CoinMemcpyN(deltaSL_, numberTotal, saveSL);
     670      CoinMemcpyN(deltaSU_, numberTotal, saveSU);
    684671#ifdef HALVE
    685                CoinWorkDouble savePrimalStep = actualPrimalStep_;
    686                CoinWorkDouble saveDualStep = actualDualStep_;
    687                CoinWorkDouble saveMu = mu_;
    688 #endif
    689                //set up for next step
    690                setupForSolve(phase);
    691                CoinWorkDouble directionAccuracy2 = findDirectionVector(phase);
    692                if (directionAccuracy2 > worstDirectionAccuracy_) {
    693                     worstDirectionAccuracy_ = directionAccuracy2;
    694                }
    695                CoinWorkDouble testValue = 1.0e2 * directionAccuracy;
    696                if (1.0e2 * projectionTolerance_ > testValue) {
    697                     testValue = 1.0e2 * projectionTolerance_;
    698                }
    699                if (primalTolerance() > testValue) {
    700                     testValue = primalTolerance();
    701                }
    702                if (maximumRHSError_ > testValue) {
    703                     testValue = maximumRHSError_;
    704                }
    705                if (directionAccuracy2 > testValue && numberIterations_ >= -77) {
    706                     goodMove = false;
     672      CoinWorkDouble savePrimalStep = actualPrimalStep_;
     673      CoinWorkDouble saveDualStep = actualDualStep_;
     674      CoinWorkDouble saveMu = mu_;
     675#endif
     676      //set up for next step
     677      setupForSolve(phase);
     678      CoinWorkDouble directionAccuracy2 = findDirectionVector(phase);
     679      if (directionAccuracy2 > worstDirectionAccuracy_) {
     680        worstDirectionAccuracy_ = directionAccuracy2;
     681      }
     682      CoinWorkDouble testValue = 1.0e2 * directionAccuracy;
     683      if (1.0e2 * projectionTolerance_ > testValue) {
     684        testValue = 1.0e2 * projectionTolerance_;
     685      }
     686      if (primalTolerance() > testValue) {
     687        testValue = primalTolerance();
     688      }
     689      if (maximumRHSError_ > testValue) {
     690        testValue = maximumRHSError_;
     691      }
     692      if (directionAccuracy2 > testValue && numberIterations_ >= -77) {
     693        goodMove = false;
    707694#ifdef SOME_DEBUG
    708                     printf("accuracy %g phase 1 failed, test value %g\n",
    709                            directionAccuracy2, testValue);
    710 #endif
    711                }
    712                if (goodMove) {
    713                     phase = 1;
    714                     CoinWorkDouble norm = findStepLength(phase);
    715                     nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
    716                     debugMove(1, actualPrimalStep_, actualDualStep_);
    717                     //debugMove(1,1.0e-7,1.0e-7);
    718                     goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
    719                     if (norm < 0)
    720                          goodMove = false;
    721                     if (!goodMove) {
     695        printf("accuracy %g phase 1 failed, test value %g\n",
     696          directionAccuracy2, testValue);
     697#endif
     698      }
     699      if (goodMove) {
     700        phase = 1;
     701        CoinWorkDouble norm = findStepLength(phase);
     702        nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
     703        debugMove(1, actualPrimalStep_, actualDualStep_);
     704        //debugMove(1,1.0e-7,1.0e-7);
     705        goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
     706        if (norm < 0)
     707          goodMove = false;
     708        if (!goodMove) {
    722709#ifdef SOME_DEBUG
    723                          printf("checkGoodMove failed\n");
    724 #endif
    725                     }
    726                }
     710          printf("checkGoodMove failed\n");
     711#endif
     712        }
     713      }
    727714#ifdef HALVE
    728                int nHalve = 0;
    729                // relax test
    730                bestNextGap = CoinMax(bestNextGap, 0.9 * complementarityGap_);
    731                while (!goodMove) {
    732                     mu_ = saveMu;
    733                     actualPrimalStep_ = savePrimalStep;
    734                     actualDualStep_ = saveDualStep;
    735                     int i;
    736                     //printf("halve %d\n",nHalve);
    737                     nHalve++;
    738                     const CoinWorkDouble lambda = 0.5;
    739                     for (i = 0; i < numberRows_; i++)
    740                          deltaY_[i] = lambda * deltaY_[i] + (1.0 - lambda) * saveY[i];
    741                     for (i = 0; i < numberTotal; i++) {
    742                          deltaX_[i] = lambda * deltaX_[i] + (1.0 - lambda) * saveX[i];
    743                          deltaZ_[i] = lambda * deltaZ_[i] + (1.0 - lambda) * saveZ[i];
    744                          deltaW_[i] = lambda * deltaW_[i] + (1.0 - lambda) * saveW[i];
    745                          deltaSL_[i] = lambda * deltaSL_[i] + (1.0 - lambda) * saveSL[i];
    746                          deltaSU_[i] = lambda * deltaSU_[i] + (1.0 - lambda) * saveSU[i];
    747                     }
    748                     CoinMemcpyN(saveX, numberTotal, deltaX_);
    749                     CoinMemcpyN(saveY, numberRows_, deltaY_);
    750                     CoinMemcpyN(saveZ, numberTotal, deltaZ_);
    751                     CoinMemcpyN(saveW, numberTotal, deltaW_);
    752                     CoinMemcpyN(saveSL, numberTotal, deltaSL_);
    753                     CoinMemcpyN(saveSU, numberTotal, deltaSU_);
    754                     findStepLength(1);
    755                     nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
    756                     goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
    757                     if (nHalve > 10)
    758                          break;
    759                     //assert (goodMove);
    760                }
    761                if (nHalve && handler_->logLevel() > 2)
    762                     printf("halved %d times\n", nHalve);
    763 #endif
    764           }
    765           //bestPhase=-1;
    766           //goodMove=false;
    767           if (!goodMove) {
    768                // Just primal dual step
    769                CoinWorkDouble floatNumber;
    770                floatNumber = 2.0 * numberComplementarityPairs_;
    771                //floatNumber = numberComplementarityItems_;
    772                CoinWorkDouble saveMu = mu_; // use one from predictor corrector
    773                mu_ = complementarityGap_ / floatNumber;
    774                // If going well try small mu
    775                mu_ *= CoinSqrt((1.0 - lastStep) / (1.0 + 10.0 * lastStep));
    776                CoinWorkDouble mu1 = mu_;
    777                CoinWorkDouble phi;
    778                if (numberComplementarityPairs_ <= 500) {
    779                     phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 2.0);
    780                } else {
    781                     phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 1.5);
    782                     if (phi < 500.0 * 500.0) {
    783                          phi = 500.0 * 500.0;
    784                     }
    785                }
    786                mu_ = complementarityGap_ / phi;
    787                //printf("pd mu %g, alternate %g, smallest\n",mu_,mu1);
    788                mu_ = CoinSqrt(mu_ * mu1);
    789                mu_ = mu1;
    790                if ((numberIterations_ & 1) == 0 || numberIterations_ < 10)
    791                     mu_ = saveMu;
    792                // Try simpler
    793                floatNumber = numberComplementarityItems_;
    794                mu_ = 0.5 * complementarityGap_ / floatNumber;
    795                //if ((modeSwitch&2)==0) {
    796                //if ((numberIterations_&1)==0)
    797                //  mu_ *= 0.5;
    798                //} else {
    799                //mu_ *= 0.8;
    800                //}
    801                //set up for next step
    802                setupForSolve(2);
    803                findDirectionVector(2);
    804                CoinWorkDouble norm = findStepLength(2);
    805                // just for debug
    806                bestNextGap = complementarityGap_ * 1.0005;
    807                //bestNextGap=COIN_DBL_MAX;
    808                nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
    809                debugMove(2, actualPrimalStep_, actualDualStep_);
    810                //debugMove(2,1.0e-7,1.0e-7);
    811                checkGoodMove(false, bestNextGap, allowIncreasingGap);
    812                if ((nextGap > 0.9 * complementarityGap_ && bestPhase == 0 && affineGap < nextGap
    813                          && (numberIterations_ > 80 || (numberIterations_ > 20 && quadraticObj))) || norm < 0.0) {
    814                     // Back to affine
    815                     phase = 0;
    816                     setupForSolve(phase);
    817                     directionAccuracy = findDirectionVector(phase);
    818                     findStepLength(phase);
    819                     nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
    820                     bestNextGap = complementarityGap_;
    821                     //checkGoodMove(false, bestNextGap,allowIncreasingGap);
    822                }
    823           }
    824           if (!goodMove)
    825                mu_ = nextGap / (static_cast<CoinWorkDouble> (nextNumber) * 1.1);
    826           //if (quadraticObj)
    827           //goodMove=true;
    828           //goodMove=false; //TEMP
    829           // Do centering steps
    830           int numberTries = 0;
    831           int numberGoodTries = 0;
     715      int nHalve = 0;
     716      // relax test
     717      bestNextGap = CoinMax(bestNextGap, 0.9 * complementarityGap_);
     718      while (!goodMove) {
     719        mu_ = saveMu;
     720        actualPrimalStep_ = savePrimalStep;
     721        actualDualStep_ = saveDualStep;
     722        int i;
     723        //printf("halve %d\n",nHalve);
     724        nHalve++;
     725        const CoinWorkDouble lambda = 0.5;
     726        for (i = 0; i < numberRows_; i++)
     727          deltaY_[i] = lambda * deltaY_[i] + (1.0 - lambda) * saveY[i];
     728        for (i = 0; i < numberTotal; i++) {
     729          deltaX_[i] = lambda * deltaX_[i] + (1.0 - lambda) * saveX[i];
     730          deltaZ_[i] = lambda * deltaZ_[i] + (1.0 - lambda) * saveZ[i];
     731          deltaW_[i] = lambda * deltaW_[i] + (1.0 - lambda) * saveW[i];
     732          deltaSL_[i] = lambda * deltaSL_[i] + (1.0 - lambda) * saveSL[i];
     733          deltaSU_[i] = lambda * deltaSU_[i] + (1.0 - lambda) * saveSU[i];
     734        }
     735        CoinMemcpyN(saveX, numberTotal, deltaX_);
     736        CoinMemcpyN(saveY, numberRows_, deltaY_);
     737        CoinMemcpyN(saveZ, numberTotal, deltaZ_);
     738        CoinMemcpyN(saveW, numberTotal, deltaW_);
     739        CoinMemcpyN(saveSL, numberTotal, deltaSL_);
     740        CoinMemcpyN(saveSU, numberTotal, deltaSU_);
     741        findStepLength(1);
     742        nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
     743        goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
     744        if (nHalve > 10)
     745          break;
     746        //assert (goodMove);
     747      }
     748      if (nHalve && handler_->logLevel() > 2)
     749        printf("halved %d times\n", nHalve);
     750#endif
     751    }
     752    //bestPhase=-1;
     753    //goodMove=false;
     754    if (!goodMove) {
     755      // Just primal dual step
     756      CoinWorkDouble floatNumber;
     757      floatNumber = 2.0 * numberComplementarityPairs_;
     758      //floatNumber = numberComplementarityItems_;
     759      CoinWorkDouble saveMu = mu_; // use one from predictor corrector
     760      mu_ = complementarityGap_ / floatNumber;
     761      // If going well try small mu
     762      mu_ *= CoinSqrt((1.0 - lastStep) / (1.0 + 10.0 * lastStep));
     763      CoinWorkDouble mu1 = mu_;
     764      CoinWorkDouble phi;
     765      if (numberComplementarityPairs_ <= 500) {
     766        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 2.0);
     767      } else {
     768        phi = pow(static_cast< CoinWorkDouble >(numberComplementarityPairs_), 1.5);
     769        if (phi < 500.0 * 500.0) {
     770          phi = 500.0 * 500.0;
     771        }
     772      }
     773      mu_ = complementarityGap_ / phi;
     774      //printf("pd mu %g, alternate %g, smallest\n",mu_,mu1);
     775      mu_ = CoinSqrt(mu_ * mu1);
     776      mu_ = mu1;
     777      if ((numberIterations_ & 1) == 0 || numberIterations_ < 10)
     778        mu_ = saveMu;
     779      // Try simpler
     780      floatNumber = numberComplementarityItems_;
     781      mu_ = 0.5 * complementarityGap_ / floatNumber;
     782      //if ((modeSwitch&2)==0) {
     783      //if ((numberIterations_&1)==0)
     784      //  mu_ *= 0.5;
     785      //} else {
     786      //mu_ *= 0.8;
     787      //}
     788      //set up for next step
     789      setupForSolve(2);
     790      findDirectionVector(2);
     791      CoinWorkDouble norm = findStepLength(2);
     792      // just for debug
     793      bestNextGap = complementarityGap_ * 1.0005;
     794      //bestNextGap=COIN_DBL_MAX;
     795      nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
     796      debugMove(2, actualPrimalStep_, actualDualStep_);
     797      //debugMove(2,1.0e-7,1.0e-7);
     798      checkGoodMove(false, bestNextGap, allowIncreasingGap);
     799      if ((nextGap > 0.9 * complementarityGap_ && bestPhase == 0 && affineGap < nextGap
     800            && (numberIterations_ > 80 || (numberIterations_ > 20 && quadraticObj)))
     801        || norm < 0.0) {
     802        // Back to affine
     803        phase = 0;
     804        setupForSolve(phase);
     805        directionAccuracy = findDirectionVector(phase);
     806        findStepLength(phase);
     807        nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
     808        bestNextGap = complementarityGap_;
     809        //checkGoodMove(false, bestNextGap,allowIncreasingGap);
     810      }
     811    }
     812    if (!goodMove)
     813      mu_ = nextGap / (static_cast< CoinWorkDouble >(nextNumber) * 1.1);
     814    //if (quadraticObj)
     815    //goodMove=true;
     816    //goodMove=false; //TEMP
     817    // Do centering steps
     818    int numberTries = 0;
     819    int numberGoodTries = 0;
    832820#ifdef COIN_DETAIL
    833           CoinWorkDouble nextCenterGap = 0.0;
    834           CoinWorkDouble originalDualStep = actualDualStep_;
    835           CoinWorkDouble originalPrimalStep = actualPrimalStep_;
    836 #endif
    837           if (actualDualStep_ > 0.9 && actualPrimalStep_ > 0.9)
    838                goodMove = false; // don't bother
    839           if ((modeSwitch & 1) != 0)
    840                goodMove = false;
    841           while (goodMove && numberTries < 5) {
    842                goodMove = false;
    843                numberTries++;
    844                CoinMemcpyN(deltaX_, numberTotal, saveX);
    845                CoinMemcpyN(deltaY_, numberRows_, saveY);
    846                CoinMemcpyN(deltaZ_, numberTotal, saveZ);
    847                CoinMemcpyN(deltaW_, numberTotal, saveW);
    848                CoinWorkDouble savePrimalStep = actualPrimalStep_;
    849                CoinWorkDouble saveDualStep = actualDualStep_;
    850                CoinWorkDouble saveMu = mu_;
    851                setupForSolve(3);
    852                findDirectionVector(3);
    853                findStepLength(3);
    854                debugMove(3, actualPrimalStep_, actualDualStep_);
    855                //debugMove(3,1.0e-7,1.0e-7);
    856                CoinWorkDouble xGap = complementarityGap(nextNumber, nextNumberItems, 3);
    857                // If one small then that's the one that counts
    858                CoinWorkDouble checkDual = saveDualStep;
    859                CoinWorkDouble checkPrimal = savePrimalStep;
    860                if (checkDual > 5.0 * checkPrimal) {
    861                     checkDual = 2.0 * checkPrimal;
    862                } else if (checkPrimal > 5.0 * checkDual) {
    863                     checkPrimal = 2.0 * checkDual;
    864                }
    865                if (actualPrimalStep_ < checkPrimal ||
    866                          actualDualStep_ < checkDual ||
    867                          (xGap > nextGap && xGap > 0.9 * complementarityGap_)) {
    868                     //if (actualPrimalStep_<=checkPrimal||
    869                     //actualDualStep_<=checkDual) {
     821    CoinWorkDouble nextCenterGap = 0.0;
     822    CoinWorkDouble originalDualStep = actualDualStep_;
     823    CoinWorkDouble originalPrimalStep = actualPrimalStep_;
     824#endif
     825    if (actualDualStep_ > 0.9 && actualPrimalStep_ > 0.9)
     826      goodMove = false; // don't bother
     827    if ((modeSwitch & 1) != 0)
     828      goodMove = false;
     829    while (goodMove && numberTries < 5) {
     830      goodMove = false;
     831      numberTries++;
     832      CoinMemcpyN(deltaX_, numberTotal, saveX);
     833      CoinMemcpyN(deltaY_, numberRows_, saveY);
     834      CoinMemcpyN(deltaZ_, numberTotal, saveZ);
     835      CoinMemcpyN(deltaW_, numberTotal, saveW);
     836      CoinWorkDouble savePrimalStep = actualPrimalStep_;
     837      CoinWorkDouble saveDualStep = actualDualStep_;
     838      CoinWorkDouble saveMu = mu_;
     839      setupForSolve(3);
     840      findDirectionVector(3);
     841      findStepLength(3);
     842      debugMove(3, actualPrimalStep_, actualDualStep_);
     843      //debugMove(3,1.0e-7,1.0e-7);
     844      CoinWorkDouble xGap = complementarityGap(nextNumber, nextNumberItems, 3);
     845      // If one small then that's the one that counts
     846      CoinWorkDouble checkDual = saveDualStep;
     847      CoinWorkDouble checkPrimal = savePrimalStep;
     848      if (checkDual > 5.0 * checkPrimal) {
     849        checkDual = 2.0 * checkPrimal;
     850      } else if (checkPrimal > 5.0 * checkDual) {
     851        checkPrimal = 2.0 * checkDual;
     852      }
     853      if (actualPrimalStep_ < checkPrimal || actualDualStep_ < checkDual || (xGap > nextGap && xGap > 0.9 * complementarityGap_)) {
     854        //if (actualPrimalStep_<=checkPrimal||
     855        //actualDualStep_<=checkDual) {
    870856#ifdef SOME_DEBUG
    871                     printf("PP rejected gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
    872                            actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
    873 #endif
    874                     mu_ = saveMu;
    875                     actualPrimalStep_ = savePrimalStep;
    876                     actualDualStep_ = saveDualStep;
    877                     CoinMemcpyN(saveX, numberTotal, deltaX_);
    878                     CoinMemcpyN(saveY, numberRows_, deltaY_);
    879                     CoinMemcpyN(saveZ, numberTotal, deltaZ_);
    880                     CoinMemcpyN(saveW, numberTotal, deltaW_);
    881                } else {
     857        printf("PP rejected gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
     858          actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
     859#endif
     860        mu_ = saveMu;
     861        actualPrimalStep_ = savePrimalStep;
     862        actualDualStep_ = saveDualStep;
     863        CoinMemcpyN(saveX, numberTotal, deltaX_);
     864        CoinMemcpyN(saveY, numberRows_, deltaY_);
     865        CoinMemcpyN(saveZ, numberTotal, deltaZ_);
     866        CoinMemcpyN(saveW, numberTotal, deltaW_);
     867      } else {
    882868#ifdef SOME_DEBUG
    883                     printf("PPphase 3 gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
    884                            actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
    885 #endif
    886                     numberGoodTries++;
     869        printf("PPphase 3 gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
     870          actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
     871#endif
     872        numberGoodTries++;
    887873#ifdef COIN_DETAIL
    888                     nextCenterGap = xGap;
    889 #endif
    890                     // See if big enough change
    891                     if (actualPrimalStep_ < 1.01 * checkPrimal ||
    892                               actualDualStep_ < 1.01 * checkDual) {
    893                          // stop now
    894                     } else {
    895                          // carry on
    896                          goodMove = true;
    897                     }
    898                }
    899           }
    900           if (numberGoodTries && handler_->logLevel() > 1) {
    901                COIN_DETAIL_PRINT(printf("%d centering steps moved from (gap %.18g, dual %.18g, primal %.18g) to (gap %.18g, dual %.18g, primal %.18g)\n",
    902                       numberGoodTries, static_cast<double>(nextGap), static_cast<double>(originalDualStep),
    903                       static_cast<double>(originalPrimalStep),
    904                       static_cast<double>(nextCenterGap), static_cast<double>(actualDualStep_),
    905                                         static_cast<double>(actualPrimalStep_)));
    906           }
    907           // save last gap
    908           checkGap = complementarityGap_;
    909           numberFixed = updateSolution(nextGap);
    910           numberFixedTotal += numberFixed;
    911      } /* endwhile */
    912      delete [] saveX;
    913      delete [] saveY;
    914      delete [] saveZ;
    915      delete [] saveW;
    916      delete [] saveSL;
    917      delete [] saveSU;
    918      if (savePi) {
    919           if (numberIterations_ - saveIteration > 20 &&
    920                     numberIterations_ - saveIteration2 < 5) {
    921 #if KEEP_GOING_IF_FIXED<10
    922                std::cout << "Restoring2 from iteration " << saveIteration2 << std::endl;
    923 #endif
    924                CoinMemcpyN(savePi2, numberRows_, dualArray);
    925                CoinMemcpyN(savePrimal2, numberTotal, solution_);
    926           } else {
    927 #if KEEP_GOING_IF_FIXED<10
    928                std::cout << "Restoring from iteration " << saveIteration << std::endl;
    929 #endif
    930                CoinMemcpyN(savePi, numberRows_, dualArray);
    931                CoinMemcpyN(savePrimal, numberTotal, solution_);
    932           }
    933           delete [] savePi;
    934           delete [] savePrimal;
    935      }
    936      delete [] savePi2;
    937      delete [] savePrimal2;
    938      //recompute slacks
    939      // Split out solution
    940      CoinZeroN(rowActivity_, numberRows_);
    941      CoinMemcpyN(solution_, numberColumns_, columnActivity_);
    942      matrix_->times(1.0, columnActivity_, rowActivity_);
    943      //unscale objective
    944      multiplyAdd(NULL, numberTotal, 0.0, cost_, scaleFactor_);
    945      multiplyAdd(NULL, numberRows_, 0, dualArray, scaleFactor_);
    946      checkSolution();
    947      //CoinMemcpyN(reducedCost_,numberColumns_,dj_);
    948      // If quadratic use last solution
    949      // Restore quadratic objective if necessary
    950      if (saveObjective) {
    951           delete objective_;
    952           objective_ = saveObjective;
    953           objectiveValue_ = 0.5 * (primalObjective_ + dualObjective_);
    954      }
    955      handler_->message(CLP_BARRIER_END, messages_)
    956                << static_cast<double>(sumPrimalInfeasibilities_)
    957                << static_cast<double>(sumDualInfeasibilities_)
    958                << static_cast<double>(complementarityGap_)
    959                << static_cast<double>(objectiveValue())
    960                << CoinMessageEol;
    961      //#ifdef SOME_DEBUG
    962      if (handler_->logLevel() > 1)
    963        COIN_DETAIL_PRINT(printf("ENDRUN status %d after %d iterations\n", problemStatus_, numberIterations_));
    964      //#endif
    965      //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
    966      //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
    967      //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
    968      //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
    969      //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
     874        nextCenterGap = xGap;
     875#endif
     876        // See if big enough change
     877        if (actualPrimalStep_ < 1.01 * checkPrimal || actualDualStep_ < 1.01 * checkDual) {
     878          // stop now
     879        } else {
     880          // carry on
     881          goodMove = true;
     882        }
     883      }
     884    }
     885    if (numberGoodTries && handler_->logLevel() > 1) {
     886      COIN_DETAIL_PRINT(printf("%d centering steps moved from (gap %.18g, dual %.18g, primal %.18g) to (gap %.18g, dual %.18g, primal %.18g)\n",
     887        numberGoodTries, static_cast< double >(nextGap), static_cast< double >(originalDualStep),
     888        static_cast< double >(originalPrimalStep),
     889        static_cast< double >(nextCenterGap), static_cast< double >(actualDualStep_),
     890        static_cast< double >(actualPrimalStep_)));
     891    }
     892    // save last gap
     893    checkGap = complementarityGap_;
     894    numberFixed = updateSolution(nextGap);
     895    numberFixedTotal += numberFixed;
     896  } /* endwhile */
     897  delete[] saveX;
     898  delete[] saveY;
     899  delete[] saveZ;
     900  delete[] saveW;
     901  delete[] saveSL;
     902  delete[] saveSU;
     903  if (savePi) {
     904    if (numberIterations_ - saveIteration > 20 && numberIterations_ - saveIteration2 < 5) {
     905#if KEEP_GOING_IF_FIXED < 10
     906      std::cout << "Restoring2 from iteration " << saveIteration2 << std::endl;
     907#endif
     908      CoinMemcpyN(savePi2, numberRows_, dualArray);
     909      CoinMemcpyN(savePrimal2, numberTotal, solution_);
     910    } else {
     911#if KEEP_GOING_IF_FIXED < 10
     912      std::cout << "Restoring from iteration " << saveIteration << std::endl;
     913#endif
     914      CoinMemcpyN(savePi, numberRows_, dualArray);
     915      CoinMemcpyN(savePrimal, numberTotal, solution_);
     916    }
     917    delete[] savePi;
     918    delete[] savePrimal;
     919  }
     920  delete[] savePi2;
     921  delete[] savePrimal2;
     922  //recompute slacks
     923  // Split out solution
     924  CoinZeroN(rowActivity_, numberRows_);
     925  CoinMemcpyN(solution_, numberColumns_, columnActivity_);
     926  matrix_->times(1.0, columnActivity_, rowActivity_);
     927  //unscale objective
     928  multiplyAdd(NULL, numberTotal, 0.0, cost_, scaleFactor_);
     929  multiplyAdd(NULL, numberRows_, 0, dualArray, scaleFactor_);
     930  checkSolution();
     931  //CoinMemcpyN(reducedCost_,numberColumns_,dj_);
     932  // If quadratic use last solution
     933  // Restore quadratic objective if necessary
     934  if (saveObjective) {
     935    delete objective_;
     936    objective_ = saveObjective;
     937    objectiveValue_ = 0.5 * (primalObjective_ + dualObjective_);
     938  }
     939  handler_->message(CLP_BARRIER_END, messages_)
     940    << static_cast< double >(sumPrimalInfeasibilities_)
     941    << static_cast< double >(sumDualInfeasibilities_)
     942    << static_cast< double >(complementarityGap_)
     943    << static_cast< double >(objectiveValue())
     944    << CoinMessageEol;
     945  //#ifdef SOME_DEBUG
     946  if (handler_->logLevel() > 1)
     947    COIN_DETAIL_PRINT(printf("ENDRUN status %d after %d iterations\n", problemStatus_, numberIterations_));
     948    //#endif
     949    //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
     950    //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
     951    //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
     952    //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
     953    //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
    970954#if COIN_LONG_WORK
    971      // put back dual
    972      delete [] dual_;
    973      delete [] reducedCost_;
    974      dual_ = dualSave;
    975      reducedCost_ = reducedCostSave;
    976 #endif
    977      //delete all temporary regions
    978      deleteWorkingData();
    979 #if KEEP_GOING_IF_FIXED<10
     955  // put back dual
     956  delete[] dual_;
     957  delete[] reducedCost_;
     958  dual_ = dualSave;
     959  reducedCost_ = reducedCostSave;
     960#endif
     961  //delete all temporary regions
     962  deleteWorkingData();
     963#if KEEP_GOING_IF_FIXED < 10
    980964#if 0 //ndef NDEBUG
    981965     {
     
    995979#endif
    996980#endif
    997      if (saveMatrix) {
    998           // restore normal copy
    999           delete matrix_;
    1000           matrix_ = saveMatrix;
    1001      }
    1002      return problemStatus_;
     981  if (saveMatrix) {
     982    // restore normal copy
     983    delete matrix_;
     984    matrix_ = saveMatrix;
     985  }
     986  return problemStatus_;
    1003987}
    1004988// findStepLength.
     
    1006990//         1 corrector
    1007991//         2 primal dual
    1008 CoinWorkDouble ClpPredictorCorrector::findStepLength( int phase)
     992CoinWorkDouble ClpPredictorCorrector::findStepLength(int phase)
    1009993{
    1010      CoinWorkDouble directionNorm = 0.0;
    1011      CoinWorkDouble maximumPrimalStep = COIN_DBL_MAX * 1.0e-20;
    1012      CoinWorkDouble maximumDualStep = COIN_DBL_MAX;
    1013      int numberTotal = numberRows_ + numberColumns_;
    1014      CoinWorkDouble tolerance = 1.0e-12;
     994  CoinWorkDouble directionNorm = 0.0;
     995  CoinWorkDouble maximumPrimalStep = COIN_DBL_MAX * 1.0e-20;
     996  CoinWorkDouble maximumDualStep = COIN_DBL_MAX;
     997  int numberTotal = numberRows_ + numberColumns_;
     998  CoinWorkDouble tolerance = 1.0e-12;
    1015999#ifdef SOME_DEBUG
    1016      int chosenPrimalSequence = -1;
    1017      int chosenDualSequence = -1;
    1018      bool lowPrimal = false;
    1019      bool lowDual = false;
    1020 #endif
    1021      // If done many iterations then allow to hit boundary
    1022      CoinWorkDouble hitTolerance;
    1023      //printf("objective norm %g\n",objectiveNorm_);
    1024      if (numberIterations_ < 80 || !gonePrimalFeasible_)
    1025           hitTolerance = COIN_DBL_MAX;
    1026      else
    1027           hitTolerance = CoinMax(1.0e3, 1.0e-3 * objectiveNorm_);
    1028      int iColumn;
    1029      //printf("dual value %g\n",dual_[0]);
    1030      //printf("     X     dX      lS     dlS     uS     dUs    dj    Z dZ     t   dT\n");
    1031      for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    1032           if (!flagged(iColumn)) {
    1033                CoinWorkDouble directionElement = deltaX_[iColumn];
    1034                if (directionNorm < CoinAbs(directionElement)) {
    1035                     directionNorm = CoinAbs(directionElement);
    1036                }
    1037                if (lowerBound(iColumn)) {
    1038                     CoinWorkDouble delta = - deltaSL_[iColumn];
    1039                     CoinWorkDouble z1 = deltaZ_[iColumn];
    1040                     CoinWorkDouble newZ = zVec_[iColumn] + z1;
    1041                     if (zVec_[iColumn] > tolerance) {
    1042                          if (zVec_[iColumn] < -z1 * maximumDualStep) {
    1043                               maximumDualStep = -zVec_[iColumn] / z1;
     1000  int chosenPrimalSequence = -1;
     1001  int chosenDualSequence = -1;
     1002  bool lowPrimal = false;
     1003  bool lowDual = false;
     1004#endif
     1005  // If done many iterations then allow to hit boundary
     1006  CoinWorkDouble hitTolerance;
     1007  //printf("objective norm %g\n",objectiveNorm_);
     1008  if (numberIterations_ < 80 || !gonePrimalFeasible_)
     1009    hitTolerance = COIN_DBL_MAX;
     1010  else
     1011    hitTolerance = CoinMax(1.0e3, 1.0e-3 * objectiveNorm_);
     1012  int iColumn;
     1013  //printf("dual value %g\n",dual_[0]);
     1014  //printf("     X     dX      lS     dlS     uS     dUs    dj    Z dZ     t   dT\n");
     1015  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     1016    if (!flagged(iColumn)) {
     1017      CoinWorkDouble directionElement = deltaX_[iColumn];
     1018      if (directionNorm < CoinAbs(directionElement)) {
     1019        directionNorm = CoinAbs(directionElement);
     1020      }
     1021      if (lowerBound(iColumn)) {
     1022        CoinWorkDouble delta = -deltaSL_[iColumn];
     1023        CoinWorkDouble z1 = deltaZ_[iColumn];
     1024        CoinWorkDouble newZ = zVec_[iColumn] + z1;
     1025        if (zVec_[iColumn] > tolerance) {
     1026          if (zVec_[iColumn] < -z1 * maximumDualStep) {
     1027            maximumDualStep = -zVec_[iColumn] / z1;
    10441028#ifdef SOME_DEBUG
    1045                               chosenDualSequence = iColumn;
    1046                               lowDual = true;
    1047 #endif
    1048                          }
    1049                     }
    1050                     if (lowerSlack_[iColumn] < maximumPrimalStep * delta) {
    1051                          CoinWorkDouble newStep = lowerSlack_[iColumn] / delta;
    1052                          if (newStep > 0.2 || newZ < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] < hitTolerance) {
    1053                               maximumPrimalStep = newStep;
     1029            chosenDualSequence = iColumn;
     1030            lowDual = true;
     1031#endif
     1032          }
     1033        }
     1034        if (lowerSlack_[iColumn] < maximumPrimalStep * delta) {
     1035          CoinWorkDouble newStep = lowerSlack_[iColumn] / delta;
     1036          if (newStep > 0.2 || newZ < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] < hitTolerance) {
     1037            maximumPrimalStep = newStep;
    10541038#ifdef SOME_DEBUG
    1055                               chosenPrimalSequence = iColumn;
    1056                               lowPrimal = true;
    1057 #endif
    1058                          } else {
    1059                               //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
    1060                          }
    1061                     }
    1062                }
    1063                if (upperBound(iColumn)) {
    1064                     CoinWorkDouble delta = - deltaSU_[iColumn];;
    1065                     CoinWorkDouble w1 = deltaW_[iColumn];
    1066                     CoinWorkDouble newT = wVec_[iColumn] + w1;
    1067                     if (wVec_[iColumn] > tolerance) {
    1068                          if (wVec_[iColumn] < -w1 * maximumDualStep) {
    1069                               maximumDualStep = -wVec_[iColumn] / w1;
     1039            chosenPrimalSequence = iColumn;
     1040            lowPrimal = true;
     1041#endif
     1042          } else {
     1043            //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
     1044          }
     1045        }
     1046      }
     1047      if (upperBound(iColumn)) {
     1048        CoinWorkDouble delta = -deltaSU_[iColumn];
     1049        ;
     1050        CoinWorkDouble w1 = deltaW_[iColumn];
     1051        CoinWorkDouble newT = wVec_[iColumn] + w1;
     1052        if (wVec_[iColumn] > tolerance) {
     1053          if (wVec_[iColumn] < -w1 * maximumDualStep) {
     1054            maximumDualStep = -wVec_[iColumn] / w1;
    10701055#ifdef SOME_DEBUG
    1071                               chosenDualSequence = iColumn;
    1072                               lowDual = false;
    1073 #endif
    1074                          }
    1075                     }
    1076                     if (upperSlack_[iColumn] < maximumPrimalStep * delta) {
    1077                          CoinWorkDouble newStep = upperSlack_[iColumn] / delta;
    1078                          if (newStep > 0.2 || newT < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] > -hitTolerance) {
    1079                               maximumPrimalStep = newStep;
     1056            chosenDualSequence = iColumn;
     1057            lowDual = false;
     1058#endif
     1059          }
     1060        }
     1061        if (upperSlack_[iColumn] < maximumPrimalStep * delta) {
     1062          CoinWorkDouble newStep = upperSlack_[iColumn] / delta;
     1063          if (newStep > 0.2 || newT < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] > -hitTolerance) {
     1064            maximumPrimalStep = newStep;
    10801065#ifdef SOME_DEBUG
    1081                               chosenPrimalSequence = iColumn;
    1082                               lowPrimal = false;
    1083 #endif
    1084                          } else {
    1085                               //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
    1086                          }
    1087                     }
    1088                }
    1089           }
    1090      }
     1066            chosenPrimalSequence = iColumn;
     1067            lowPrimal = false;
     1068#endif
     1069          } else {
     1070            //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
     1071          }
     1072        }
     1073      }
     1074    }
     1075  }
    10911076#ifdef SOME_DEBUG
    1092      printf("new step - phase %d, norm %.18g, dual step %.18g, primal step %.18g\n",
    1093             phase, directionNorm, maximumDualStep, maximumPrimalStep);
    1094      if (lowDual)
    1095           printf("ld %d %g %g => %g (dj %g,sol %g) ",
    1096                  chosenDualSequence, zVec_[chosenDualSequence],
    1097                  deltaZ_[chosenDualSequence], zVec_[chosenDualSequence] +
    1098                  maximumDualStep * deltaZ_[chosenDualSequence], dj_[chosenDualSequence],
    1099                  solution_[chosenDualSequence]);
    1100      else
    1101           printf("ud %d %g %g => %g (dj %g,sol %g) ",
    1102                  chosenDualSequence, wVec_[chosenDualSequence],
    1103                  deltaW_[chosenDualSequence], wVec_[chosenDualSequence] +
    1104                  maximumDualStep * deltaW_[chosenDualSequence], dj_[chosenDualSequence],
    1105                  solution_[chosenDualSequence]);
    1106      if (lowPrimal)
    1107           printf("lp %d %g %g => %g (dj %g,sol %g)\n",
    1108                  chosenPrimalSequence, lowerSlack_[chosenPrimalSequence],
    1109                  deltaSL_[chosenPrimalSequence], lowerSlack_[chosenPrimalSequence] +
    1110                  maximumPrimalStep * deltaSL_[chosenPrimalSequence],
    1111                  dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
    1112      else
    1113           printf("up %d %g %g => %g (dj %g,sol %g)\n",
    1114                  chosenPrimalSequence, upperSlack_[chosenPrimalSequence],
    1115                  deltaSU_[chosenPrimalSequence], upperSlack_[chosenPrimalSequence] +
    1116                  maximumPrimalStep * deltaSU_[chosenPrimalSequence],
    1117                  dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
    1118 #endif
    1119      actualPrimalStep_ = stepLength_ * maximumPrimalStep;
    1120      if (phase >= 0 && actualPrimalStep_ > 1.0) {
    1121           actualPrimalStep_ = 1.0;
    1122      }
    1123      actualDualStep_ = stepLength_ * maximumDualStep;
    1124      if (phase >= 0 && actualDualStep_ > 1.0) {
    1125           actualDualStep_ = 1.0;
    1126      }
    1127      // See if quadratic objective
     1077  printf("new step - phase %d, norm %.18g, dual step %.18g, primal step %.18g\n",
     1078    phase, directionNorm, maximumDualStep, maximumPrimalStep);
     1079  if (lowDual)
     1080    printf("ld %d %g %g => %g (dj %g,sol %g) ",
     1081      chosenDualSequence, zVec_[chosenDualSequence],
     1082      deltaZ_[chosenDualSequence], zVec_[chosenDualSequence] + maximumDualStep * deltaZ_[chosenDualSequence], dj_[chosenDualSequence],
     1083      solution_[chosenDualSequence]);
     1084  else
     1085    printf("ud %d %g %g => %g (dj %g,sol %g) ",
     1086      chosenDualSequence, wVec_[chosenDualSequence],
     1087      deltaW_[chosenDualSequence], wVec_[chosenDualSequence] + maximumDualStep * deltaW_[chosenDualSequence], dj_[chosenDualSequence],
     1088      solution_[chosenDualSequence]);
     1089  if (lowPrimal)
     1090    printf("lp %d %g %g => %g (dj %g,sol %g)\n",
     1091      chosenPrimalSequence, lowerSlack_[chosenPrimalSequence],
     1092      deltaSL_[chosenPrimalSequence], lowerSlack_[chosenPrimalSequence] + maximumPrimalStep * deltaSL_[chosenPrimalSequence],
     1093      dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
     1094  else
     1095    printf("up %d %g %g => %g (dj %g,sol %g)\n",
     1096      chosenPrimalSequence, upperSlack_[chosenPrimalSequence],
     1097      deltaSU_[chosenPrimalSequence], upperSlack_[chosenPrimalSequence] + maximumPrimalStep * deltaSU_[chosenPrimalSequence],
     1098      dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
     1099#endif
     1100  actualPrimalStep_ = stepLength_ * maximumPrimalStep;
     1101  if (phase >= 0 && actualPrimalStep_ > 1.0) {
     1102    actualPrimalStep_ = 1.0;
     1103  }
     1104  actualDualStep_ = stepLength_ * maximumDualStep;
     1105  if (phase >= 0 && actualDualStep_ > 1.0) {
     1106    actualDualStep_ = 1.0;
     1107  }
     1108  // See if quadratic objective
    11281109#ifndef NO_RTTI
    1129      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1110  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
    11301111#else
    1131      ClpQuadraticObjective * quadraticObj = NULL;
    1132      if (objective_->type() == 2)
    1133           quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    1134 #endif
    1135      if (quadraticObj) {
    1136           // Use smaller unless very small
    1137           CoinWorkDouble smallerStep = CoinMin(actualDualStep_, actualPrimalStep_);
    1138           if (smallerStep > 0.0001) {
    1139                actualDualStep_ = smallerStep;
    1140                actualPrimalStep_ = smallerStep;
    1141           }
    1142      }
     1112  ClpQuadraticObjective *quadraticObj = NULL;
     1113  if (objective_->type() == 2)
     1114    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     1115#endif
     1116  if (quadraticObj) {
     1117    // Use smaller unless very small
     1118    CoinWorkDouble smallerStep = CoinMin(actualDualStep_, actualPrimalStep_);
     1119    if (smallerStep > 0.0001) {
     1120      actualDualStep_ = smallerStep;
     1121      actualPrimalStep_ = smallerStep;
     1122    }
     1123  }
    11431124#define OFFQ
    11441125#ifndef OFFQ
    1145      if (quadraticObj) {
    1146           // Don't bother if phase 0 or 3 or large gap
    1147           //if ((phase==1||phase==2||phase==0)&&maximumDualError_>0.1*complementarityGap_
    1148           //&&smallerStep>0.001) {
    1149           if ((phase == 1 || phase == 2 || phase == 0 || phase == 3)) {
    1150                // minimize complementarity + norm*dual inf ? primal inf
    1151                // at first - just check better - if not
    1152                // Complementarity gap will be a*change*change + b*change +c
    1153                CoinWorkDouble a = 0.0;
    1154                CoinWorkDouble b = 0.0;
    1155                CoinWorkDouble c = 0.0;
    1156                /* SQUARE of dual infeasibility will be:
     1126  if (quadraticObj) {
     1127    // Don't bother if phase 0 or 3 or large gap
     1128    //if ((phase==1||phase==2||phase==0)&&maximumDualError_>0.1*complementarityGap_
     1129    //&&smallerStep>0.001) {
     1130    if ((phase == 1 || phase == 2 || phase == 0 || phase == 3)) {
     1131      // minimize complementarity + norm*dual inf ? primal inf
     1132      // at first - just check better - if not
     1133      // Complementarity gap will be a*change*change + b*change +c
     1134      CoinWorkDouble a = 0.0;
     1135      CoinWorkDouble b = 0.0;
     1136      CoinWorkDouble c = 0.0;
     1137      /* SQUARE of dual infeasibility will be:
    11571138               square of dj - ......
    11581139               */
    1159                CoinWorkDouble aq = 0.0;
    1160                CoinWorkDouble bq = 0.0;
    1161                CoinWorkDouble cq = 0.0;
    1162                CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
    1163                CoinWorkDouble * linearDjChange = new CoinWorkDouble[numberTotal];
    1164                CoinZeroN(linearDjChange, numberColumns_);
    1165                multiplyAdd(deltaY_, numberRows_, 1.0, linearDjChange + numberColumns_, 0.0);
    1166                matrix_->transposeTimes(-1.0, deltaY_, linearDjChange);
    1167                CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    1168                const int * columnQuadratic = quadratic->getIndices();
    1169                const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    1170                const int * columnQuadraticLength = quadratic->getVectorLengths();
    1171                CoinWorkDouble * quadraticElement = quadratic->getMutableElements();
    1172                for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    1173                     CoinWorkDouble oldPrimal = solution_[iColumn];
    1174                     if (!flagged(iColumn)) {
    1175                          if (lowerBound(iColumn)) {
    1176                               CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
    1177                               c += lowerSlack_[iColumn] * zVec_[iColumn];
    1178                               b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
    1179                               a += deltaZ_[iColumn] * change;
    1180                          }
    1181                          if (upperBound(iColumn)) {
    1182                               CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
    1183                               c += upperSlack_[iColumn] * wVec_[iColumn];
    1184                               b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
    1185                               a += deltaW_[iColumn] * change;
    1186                          }
    1187                          // new djs are dj_ + change*value
    1188                          CoinWorkDouble djChange = linearDjChange[iColumn];
    1189                          if (iColumn < numberColumns_) {
    1190                               for (CoinBigIndex j = columnQuadraticStart[iColumn];
    1191                                         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    1192                                    int jColumn = columnQuadratic[j];
    1193                                    CoinWorkDouble changeJ = deltaX_[jColumn];
    1194                                    CoinWorkDouble elementValue = quadraticElement[j];
    1195                                    djChange += changeJ * elementValue;
    1196                               }
    1197                          }
    1198                          CoinWorkDouble gammaTerm = gamma2;
    1199                          if (primalR_) {
    1200                               gammaTerm += primalR_[iColumn];
    1201                          }
    1202                          djChange += gammaTerm;
    1203                          // and dual infeasibility
    1204                          CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] +
    1205                                                  gammaTerm * solution_[iColumn];
    1206                          CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
    1207                          cq += oldInf * oldInf;
    1208                          bq += 2.0 * oldInf * changeInf;
    1209                          aq += changeInf * changeInf;
    1210                     } else {
    1211                          // fixed
    1212                          if (lowerBound(iColumn)) {
    1213                               c += lowerSlack_[iColumn] * zVec_[iColumn];
    1214                          }
    1215                          if (upperBound(iColumn)) {
    1216                               c += upperSlack_[iColumn] * wVec_[iColumn];
    1217                          }
    1218                          // new djs are dj_ + change*value
    1219                          CoinWorkDouble djChange = linearDjChange[iColumn];
    1220                          if (iColumn < numberColumns_) {
    1221                               for (CoinBigIndex j = columnQuadraticStart[iColumn];
    1222                                         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    1223                                    int jColumn = columnQuadratic[j];
    1224                                    CoinWorkDouble changeJ = deltaX_[jColumn];
    1225                                    CoinWorkDouble elementValue = quadraticElement[j];
    1226                                    djChange += changeJ * elementValue;
    1227                               }
    1228                          }
    1229                          CoinWorkDouble gammaTerm = gamma2;
    1230                          if (primalR_) {
    1231                               gammaTerm += primalR_[iColumn];
    1232                          }
    1233                          djChange += gammaTerm;
    1234                          // and dual infeasibility
    1235                          CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] +
    1236                                                  gammaTerm * solution_[iColumn];
    1237                          CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
    1238                          cq += oldInf * oldInf;
    1239                          bq += 2.0 * oldInf * changeInf;
    1240                          aq += changeInf * changeInf;
    1241                     }
    1242                }
    1243                delete [] linearDjChange;
    1244                // ? We want to minimize complementarityGap + solutionNorm_*square of inf ??
    1245                // maybe use inf and do line search
    1246                // To check see if matches at current step
    1247                CoinWorkDouble step = actualPrimalStep_;
    1248                //Current gap + solutionNorm_ * CoinSqrt (sum square inf)
    1249                CoinWorkDouble multiplier = solutionNorm_;
    1250                multiplier *= 0.01;
    1251                multiplier = 1.0;
    1252                CoinWorkDouble currentInf =  multiplier * CoinSqrt(cq);
    1253                CoinWorkDouble nextInf =         multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
    1254                CoinWorkDouble allowedIncrease = 1.4;
     1140      CoinWorkDouble aq = 0.0;
     1141      CoinWorkDouble bq = 0.0;
     1142      CoinWorkDouble cq = 0.0;
     1143      CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
     1144      CoinWorkDouble *linearDjChange = new CoinWorkDouble[numberTotal];
     1145      CoinZeroN(linearDjChange, numberColumns_);
     1146      multiplyAdd(deltaY_, numberRows_, 1.0, linearDjChange + numberColumns_, 0.0);
     1147      matrix_->transposeTimes(-1.0, deltaY_, linearDjChange);
     1148      CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
     1149      const int *columnQuadratic = quadratic->getIndices();
     1150      const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     1151      const int *columnQuadraticLength = quadratic->getVectorLengths();
     1152      CoinWorkDouble *quadraticElement = quadratic->getMutableElements();
     1153      for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     1154        CoinWorkDouble oldPrimal = solution_[iColumn];
     1155        if (!flagged(iColumn)) {
     1156          if (lowerBound(iColumn)) {
     1157            CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
     1158            c += lowerSlack_[iColumn] * zVec_[iColumn];
     1159            b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
     1160            a += deltaZ_[iColumn] * change;
     1161          }
     1162          if (upperBound(iColumn)) {
     1163            CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
     1164            c += upperSlack_[iColumn] * wVec_[iColumn];
     1165            b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
     1166            a += deltaW_[iColumn] * change;
     1167          }
     1168          // new djs are dj_ + change*value
     1169          CoinWorkDouble djChange = linearDjChange[iColumn];
     1170          if (iColumn < numberColumns_) {
     1171            for (CoinBigIndex j = columnQuadraticStart[iColumn];
     1172                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     1173              int jColumn = columnQuadratic[j];
     1174              CoinWorkDouble changeJ = deltaX_[jColumn];
     1175              CoinWorkDouble elementValue = quadraticElement[j];
     1176              djChange += changeJ * elementValue;
     1177            }
     1178          }
     1179          CoinWorkDouble gammaTerm = gamma2;
     1180          if (primalR_) {
     1181            gammaTerm += primalR_[iColumn];
     1182          }
     1183          djChange += gammaTerm;
     1184          // and dual infeasibility
     1185          CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] + gammaTerm * solution_[iColumn];
     1186          CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
     1187          cq += oldInf * oldInf;
     1188          bq += 2.0 * oldInf * changeInf;
     1189          aq += changeInf * changeInf;
     1190        } else {
     1191          // fixed
     1192          if (lowerBound(iColumn)) {
     1193            c += lowerSlack_[iColumn] * zVec_[iColumn];
     1194          }
     1195          if (upperBound(iColumn)) {
     1196            c += upperSlack_[iColumn] * wVec_[iColumn];
     1197          }
     1198          // new djs are dj_ + change*value
     1199          CoinWorkDouble djChange = linearDjChange[iColumn];
     1200          if (iColumn < numberColumns_) {
     1201            for (CoinBigIndex j = columnQuadraticStart[iColumn];
     1202                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     1203              int jColumn = columnQuadratic[j];
     1204              CoinWorkDouble changeJ = deltaX_[jColumn];
     1205              CoinWorkDouble elementValue = quadraticElement[j];
     1206              djChange += changeJ * elementValue;
     1207            }
     1208          }
     1209          CoinWorkDouble gammaTerm = gamma2;
     1210          if (primalR_) {
     1211            gammaTerm += primalR_[iColumn];
     1212          }
     1213          djChange += gammaTerm;
     1214          // and dual infeasibility
     1215          CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] + gammaTerm * solution_[iColumn];
     1216          CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
     1217          cq += oldInf * oldInf;
     1218          bq += 2.0 * oldInf * changeInf;
     1219          aq += changeInf * changeInf;
     1220        }
     1221      }
     1222      delete[] linearDjChange;
     1223      // ? We want to minimize complementarityGap + solutionNorm_*square of inf ??
     1224      // maybe use inf and do line search
     1225      // To check see if matches at current step
     1226      CoinWorkDouble step = actualPrimalStep_;
     1227      //Current gap + solutionNorm_ * CoinSqrt (sum square inf)
     1228      CoinWorkDouble multiplier = solutionNorm_;
     1229      multiplier *= 0.01;
     1230      multiplier = 1.0;
     1231      CoinWorkDouble currentInf = multiplier * CoinSqrt(cq);
     1232      CoinWorkDouble nextInf = multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
     1233      CoinWorkDouble allowedIncrease = 1.4;
    12551234#ifdef SOME_DEBUG
    1256                printf("lin %g %g %g -> %g\n", a, b, c,
    1257                       c + b * step + a * step * step);
    1258                printf("quad %g %g %g -> %g\n", aq, bq, cq,
    1259                       cq + bq * step + aq * step * step);
    1260                debugMove(7, step, step);
    1261                printf ("current dualInf %g, with step of %g is %g\n",
    1262                        currentInf, step, nextInf);
    1263 #endif
    1264                if (b > -1.0e-6) {
    1265                     if (phase != 0)
    1266                          directionNorm = -1.0;
    1267                }
    1268                if ((phase == 1 || phase == 2 || phase == 0 || phase == 3) && nextInf > 0.1 * complementarityGap_ &&
    1269                          nextInf > currentInf * allowedIncrease) {
    1270                     //cq = CoinMax(cq,10.0);
    1271                     // convert to (x+q)*(x+q) = w
    1272                     CoinWorkDouble q = bq / (1.0 * aq);
    1273                     CoinWorkDouble w = CoinMax(q * q + (cq / aq) * (allowedIncrease - 1.0), 0.0);
    1274                     w = CoinSqrt(w);
    1275                     CoinWorkDouble stepX = w - q;
    1276                     step = stepX;
    1277                     nextInf =
    1278                          multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
     1235      printf("lin %g %g %g -> %g\n", a, b, c,
     1236        c + b * step + a * step * step);
     1237      printf("quad %g %g %g -> %g\n", aq, bq, cq,
     1238        cq + bq * step + aq * step * step);
     1239      debugMove(7, step, step);
     1240      printf("current dualInf %g, with step of %g is %g\n",
     1241        currentInf, step, nextInf);
     1242#endif
     1243      if (b > -1.0e-6) {
     1244        if (phase != 0)
     1245          directionNorm = -1.0;
     1246      }
     1247      if ((phase == 1 || phase == 2 || phase == 0 || phase == 3) && nextInf > 0.1 * complementarityGap_ && nextInf > currentInf * allowedIncrease) {
     1248        //cq = CoinMax(cq,10.0);
     1249        // convert to (x+q)*(x+q) = w
     1250        CoinWorkDouble q = bq / (1.0 * aq);
     1251        CoinWorkDouble w = CoinMax(q * q + (cq / aq) * (allowedIncrease - 1.0), 0.0);
     1252        w = CoinSqrt(w);
     1253        CoinWorkDouble stepX = w - q;
     1254        step = stepX;
     1255        nextInf = multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
    12791256#ifdef SOME_DEBUG
    1280                     printf ("with step of %g dualInf is %g\n",
    1281                             step, nextInf);
    1282 #endif
    1283                     actualDualStep_ = CoinMin(step, actualDualStep_);
    1284                     actualPrimalStep_ = CoinMin(step, actualPrimalStep_);
    1285                }
    1286           }
    1287      } else {
    1288           // probably pointless as linear
    1289           // minimize complementarity
    1290           // Complementarity gap will be a*change*change + b*change +c
    1291           CoinWorkDouble a = 0.0;
    1292           CoinWorkDouble b = 0.0;
    1293           CoinWorkDouble c = 0.0;
    1294           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    1295                CoinWorkDouble oldPrimal = solution_[iColumn];
    1296                if (!flagged(iColumn)) {
    1297                     if (lowerBound(iColumn)) {
    1298                          CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
    1299                          c += lowerSlack_[iColumn] * zVec_[iColumn];
    1300                          b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
    1301                          a += deltaZ_[iColumn] * change;
    1302                     }
    1303                     if (upperBound(iColumn)) {
    1304                          CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
    1305                          c += upperSlack_[iColumn] * wVec_[iColumn];
    1306                          b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
    1307                          a += deltaW_[iColumn] * change;
    1308                     }
    1309                } else {
    1310                     // fixed
    1311                     if (lowerBound(iColumn)) {
    1312                          c += lowerSlack_[iColumn] * zVec_[iColumn];
    1313                     }
    1314                     if (upperBound(iColumn)) {
    1315                          c += upperSlack_[iColumn] * wVec_[iColumn];
    1316                     }
    1317                }
    1318           }
    1319           // ? We want to minimize complementarityGap;
    1320           // maybe use inf and do line search
    1321           // To check see if matches at current step
    1322           CoinWorkDouble step = CoinMin(actualPrimalStep_, actualDualStep_);
    1323           CoinWorkDouble next = c + b * step + a * step * step;
     1257        printf("with step of %g dualInf is %g\n",
     1258          step, nextInf);
     1259#endif
     1260        actualDualStep_ = CoinMin(step, actualDualStep_);
     1261        actualPrimalStep_ = CoinMin(step, actualPrimalStep_);
     1262      }
     1263    }
     1264  } else {
     1265    // probably pointless as linear
     1266    // minimize complementarity
     1267    // Complementarity gap will be a*change*change + b*change +c
     1268    CoinWorkDouble a = 0.0;
     1269    CoinWorkDouble b = 0.0;
     1270    CoinWorkDouble c = 0.0;
     1271    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     1272      CoinWorkDouble oldPrimal = solution_[iColumn];
     1273      if (!flagged(iColumn)) {
     1274        if (lowerBound(iColumn)) {
     1275          CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
     1276          c += lowerSlack_[iColumn] * zVec_[iColumn];
     1277          b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
     1278          a += deltaZ_[iColumn] * change;
     1279        }
     1280        if (upperBound(iColumn)) {
     1281          CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
     1282          c += upperSlack_[iColumn] * wVec_[iColumn];
     1283          b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
     1284          a += deltaW_[iColumn] * change;
     1285        }
     1286      } else {
     1287        // fixed
     1288        if (lowerBound(iColumn)) {
     1289          c += lowerSlack_[iColumn] * zVec_[iColumn];
     1290        }
     1291        if (upperBound(iColumn)) {
     1292          c += upperSlack_[iColumn] * wVec_[iColumn];
     1293        }
     1294      }
     1295    }
     1296    // ? We want to minimize complementarityGap;
     1297    // maybe use inf and do line search
     1298    // To check see if matches at current step
     1299    CoinWorkDouble step = CoinMin(actualPrimalStep_, actualDualStep_);
     1300    CoinWorkDouble next = c + b * step + a * step * step;
    13241301#ifdef SOME_DEBUG
    1325           printf("lin %g %g %g -> %g\n", a, b, c,
    1326                  c + b * step + a * step * step);
    1327           debugMove(7, step, step);
    1328 #endif
    1329           if (b > -1.0e-6) {
    1330                if (phase == 0) {
     1302    printf("lin %g %g %g -> %g\n", a, b, c,
     1303      c + b * step + a * step * step);
     1304    debugMove(7, step, step);
     1305#endif
     1306    if (b > -1.0e-6) {
     1307      if (phase == 0) {
    13311308#ifdef SOME_DEBUG
    1332                     printf("*** odd phase 0 direction\n");
    1333 #endif
    1334                } else {
    1335                     directionNorm = -1.0;
    1336                }
    1337           }
    1338           // and with ratio
    1339           a = 0.0;
    1340           b = 0.0;
    1341           CoinWorkDouble ratio = actualDualStep_ / actualPrimalStep_;
    1342           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    1343                CoinWorkDouble oldPrimal = solution_[iColumn];
    1344                if (!flagged(iColumn)) {
    1345                     if (lowerBound(iColumn)) {
    1346                          CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
    1347                          b += lowerSlack_[iColumn] * deltaZ_[iColumn] * ratio + zVec_[iColumn] * change;
    1348                          a += deltaZ_[iColumn] * change * ratio;
    1349                     }
    1350                     if (upperBound(iColumn)) {
    1351                          CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
    1352                          b += upperSlack_[iColumn] * deltaW_[iColumn] * ratio + wVec_[iColumn] * change;
    1353                          a += deltaW_[iColumn] * change * ratio;
    1354                     }
    1355                }
    1356           }
    1357           // ? We want to minimize complementarityGap;
    1358           // maybe use inf and do line search
    1359           // To check see if matches at current step
    1360           step = actualPrimalStep_;
    1361           CoinWorkDouble next2 = c + b * step + a * step * step;
    1362           if (next2 > next) {
    1363                actualPrimalStep_ = CoinMin(actualPrimalStep_, actualDualStep_);
    1364                actualDualStep_ = actualPrimalStep_;
    1365           }
     1309        printf("*** odd phase 0 direction\n");
     1310#endif
     1311      } else {
     1312        directionNorm = -1.0;
     1313      }
     1314    }
     1315    // and with ratio
     1316    a = 0.0;
     1317    b = 0.0;
     1318    CoinWorkDouble ratio = actualDualStep_ / actualPrimalStep_;
     1319    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     1320      CoinWorkDouble oldPrimal = solution_[iColumn];
     1321      if (!flagged(iColumn)) {
     1322        if (lowerBound(iColumn)) {
     1323          CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
     1324          b += lowerSlack_[iColumn] * deltaZ_[iColumn] * ratio + zVec_[iColumn] * change;
     1325          a += deltaZ_[iColumn] * change * ratio;
     1326        }
     1327        if (upperBound(iColumn)) {
     1328          CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
     1329          b += upperSlack_[iColumn] * deltaW_[iColumn] * ratio + wVec_[iColumn] * change;
     1330          a += deltaW_[iColumn] * change * ratio;
     1331        }
     1332      }
     1333    }
     1334    // ? We want to minimize complementarityGap;
     1335    // maybe use inf and do line search
     1336    // To check see if matches at current step
     1337    step = actualPrimalStep_;
     1338    CoinWorkDouble next2 = c + b * step + a * step * step;
     1339    if (next2 > next) {
     1340      actualPrimalStep_ = CoinMin(actualPrimalStep_, actualDualStep_);
     1341      actualDualStep_ = actualPrimalStep_;
     1342    }
    13661343#ifdef SOME_DEBUG
    1367           printf("linb %g %g %g -> %g\n", a, b, c,
    1368                  c + b * step + a * step * step);
    1369           debugMove(7, actualPrimalStep_, actualDualStep_);
    1370 #endif
    1371           if (b > -1.0e-6) {
    1372                if (phase == 0) {
     1344    printf("linb %g %g %g -> %g\n", a, b, c,
     1345      c + b * step + a * step * step);
     1346    debugMove(7, actualPrimalStep_, actualDualStep_);
     1347#endif
     1348    if (b > -1.0e-6) {
     1349      if (phase == 0) {
    13731350#ifdef SOME_DEBUG
    1374                     printf("*** odd phase 0 direction\n");
    1375 #endif
    1376                } else {
    1377                     directionNorm = -1.0;
    1378                }
    1379           }
    1380      }
     1351        printf("*** odd phase 0 direction\n");
     1352#endif
     1353      } else {
     1354        directionNorm = -1.0;
     1355      }
     1356    }
     1357  }
    13811358#else
    1382      //actualPrimalStep_ =0.5*actualDualStep_;
     1359  //actualPrimalStep_ =0.5*actualDualStep_;
    13831360#endif
    13841361#ifdef FULL_DEBUG
    1385      if (phase == 3) {
    1386           CoinWorkDouble minBeta = 0.1 * mu_;
    1387           CoinWorkDouble maxBeta = 10.0 * mu_;
    1388           for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    1389                if (!flagged(iColumn)) {
    1390                     if (lowerBound(iColumn)) {
    1391                          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
    1392                          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
    1393                          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
    1394                          CoinWorkDouble gapProduct = dualValue * primalValue;
    1395                          if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
    1396                               printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
    1397                                      iColumn, primalValue, dualValue, gapProduct, delta2Z_[iColumn]);
    1398                     }
    1399                     if (upperBound(iColumn)) {
    1400                          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
    1401                          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
    1402                          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
    1403                          CoinWorkDouble gapProduct = dualValue * primalValue;
    1404                          if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
    1405                               printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
    1406                                      iColumn, primalValue, dualValue, gapProduct, delta2W_[iColumn]);
    1407                     }
    1408                }
    1409           }
    1410      }
     1362  if (phase == 3) {
     1363    CoinWorkDouble minBeta = 0.1 * mu_;
     1364    CoinWorkDouble maxBeta = 10.0 * mu_;
     1365    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
     1366      if (!flagged(iColumn)) {
     1367        if (lowerBound(iColumn)) {
     1368          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     1369          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
     1370          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
     1371          CoinWorkDouble gapProduct = dualValue * primalValue;
     1372          if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
     1373            printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
     1374              iColumn, primalValue, dualValue, gapProduct, delta2Z_[iColumn]);
     1375        }
     1376        if (upperBound(iColumn)) {
     1377          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
     1378          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
     1379          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
     1380          CoinWorkDouble gapProduct = dualValue * primalValue;
     1381          if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
     1382            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
     1383              iColumn, primalValue, dualValue, gapProduct, delta2W_[iColumn]);
     1384        }
     1385      }
     1386    }
     1387  }
    14111388#endif
    14121389#ifdef SOME_DEBUG_not
    1413      {
    1414           CoinWorkDouble largestL = 0.0;
    1415           CoinWorkDouble smallestL = COIN_DBL_MAX;
    1416           CoinWorkDouble largestU = 0.0;
    1417           CoinWorkDouble smallestU = COIN_DBL_MAX;
    1418           CoinWorkDouble sumL = 0.0;
    1419           CoinWorkDouble sumU = 0.0;
    1420           int nL = 0;
    1421           int nU = 0;
    1422           for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    1423                if (!flagged(iColumn)) {
    1424                     if (lowerBound(iColumn)) {
    1425                          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
    1426                          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
    1427                          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
    1428                          CoinWorkDouble gapProduct = dualValue * primalValue;
    1429                          largestL = CoinMax(largestL, gapProduct);
    1430                          smallestL = CoinMin(smallestL, gapProduct);
    1431                          nL++;
    1432                          sumL += gapProduct;
    1433                     }
    1434                     if (upperBound(iColumn)) {
    1435                          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
    1436                          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
    1437                          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
    1438                          CoinWorkDouble gapProduct = dualValue * primalValue;
    1439                          largestU = CoinMax(largestU, gapProduct);
    1440                          smallestU = CoinMin(smallestU, gapProduct);
    1441                          nU++;
    1442                          sumU += gapProduct;
    1443                     }
    1444                }
    1445           }
    1446           CoinWorkDouble mu = (sumL + sumU) / (static_cast<CoinWorkDouble> (nL + nU));
     1390  {
     1391    CoinWorkDouble largestL = 0.0;
     1392    CoinWorkDouble smallestL = COIN_DBL_MAX;
     1393    CoinWorkDouble largestU = 0.0;
     1394    CoinWorkDouble smallestU = COIN_DBL_MAX;
     1395    CoinWorkDouble sumL = 0.0;
     1396    CoinWorkDouble sumU = 0.0;
     1397    int nL = 0;
     1398    int nU = 0;
     1399    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
     1400      if (!flagged(iColumn)) {
     1401        if (lowerBound(iColumn)) {
     1402          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     1403          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
     1404          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
     1405          CoinWorkDouble gapProduct = dualValue * primalValue;
     1406          largestL = CoinMax(largestL, gapProduct);
     1407          smallestL = CoinMin(smallestL, gapProduct);
     1408          nL++;
     1409          sumL += gapProduct;
     1410        }
     1411        if (upperBound(iColumn)) {
     1412          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
     1413          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
     1414          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
     1415          CoinWorkDouble gapProduct = dualValue * primalValue;
     1416          largestU = CoinMax(largestU, gapProduct);
     1417          smallestU = CoinMin(smallestU, gapProduct);
     1418          nU++;
     1419          sumU += gapProduct;
     1420        }
     1421      }
     1422    }
     1423    CoinWorkDouble mu = (sumL + sumU) / (static_cast< CoinWorkDouble >(nL + nU));
    14471424
    1448           CoinWorkDouble minBeta = 0.1 * mu;
    1449           CoinWorkDouble maxBeta = 10.0 * mu;
    1450           int nBL = 0;
    1451           int nAL = 0;
    1452           int nBU = 0;
    1453           int nAU = 0;
    1454           for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    1455                if (!flagged(iColumn)) {
    1456                     if (lowerBound(iColumn)) {
    1457                          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
    1458                          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
    1459                          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
    1460                          CoinWorkDouble gapProduct = dualValue * primalValue;
    1461                          if (gapProduct < minBeta)
    1462                               nBL++;
    1463                          else if (gapProduct > maxBeta)
    1464                               nAL++;
    1465                          //if (gapProduct<0.1*minBeta)
    1466                          //printf("Lsmall one %d dual %g primal %g\n",iColumn,
    1467                          //   dualValue,primalValue);
    1468                     }
    1469                     if (upperBound(iColumn)) {
    1470                          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
    1471                          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
    1472                          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
    1473                          CoinWorkDouble gapProduct = dualValue * primalValue;
    1474                          if (gapProduct < minBeta)
    1475                               nBU++;
    1476                          else if (gapProduct > maxBeta)
    1477                               nAU++;
    1478                          //if (gapProduct<0.1*minBeta)
    1479                          //printf("Usmall one %d dual %g primal %g\n",iColumn,
    1480                          //   dualValue,primalValue);
    1481                     }
    1482                }
    1483           }
    1484           printf("phase %d new mu %.18g new gap %.18g\n", phase, mu, sumL + sumU);
    1485           printf("          %d lower, smallest %.18g, %d below - largest %.18g, %d above\n",
    1486                  nL, smallestL, nBL, largestL, nAL);
    1487           printf("          %d upper, smallest %.18g, %d below - largest %.18g, %d above\n",
    1488                  nU, smallestU, nBU, largestU, nAU);
    1489      }
    1490 #endif
    1491      return directionNorm;
     1425    CoinWorkDouble minBeta = 0.1 * mu;
     1426    CoinWorkDouble maxBeta = 10.0 * mu;
     1427    int nBL = 0;
     1428    int nAL = 0;
     1429    int nBU = 0;
     1430    int nAU = 0;
     1431    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
     1432      if (!flagged(iColumn)) {
     1433        if (lowerBound(iColumn)) {
     1434          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     1435          CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
     1436          CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
     1437          CoinWorkDouble gapProduct = dualValue * primalValue;
     1438          if (gapProduct < minBeta)
     1439            nBL++;
     1440          else if (gapProduct > maxBeta)
     1441            nAL++;
     1442          //if (gapProduct<0.1*minBeta)
     1443          //printf("Lsmall one %d dual %g primal %g\n",iColumn,
     1444          //   dualValue,primalValue);
     1445        }
     1446        if (upperBound(iColumn)) {
     1447          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
     1448          CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
     1449          CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
     1450          CoinWorkDouble gapProduct = dualValue * primalValue;
     1451          if (gapProduct < minBeta)
     1452            nBU++;
     1453          else if (gapProduct > maxBeta)
     1454            nAU++;
     1455          //if (gapProduct<0.1*minBeta)
     1456          //printf("Usmall one %d dual %g primal %g\n",iColumn,
     1457          //   dualValue,primalValue);
     1458        }
     1459      }
     1460    }
     1461    printf("phase %d new mu %.18g new gap %.18g\n", phase, mu, sumL + sumU);
     1462    printf("          %d lower, smallest %.18g, %d below - largest %.18g, %d above\n",
     1463      nL, smallestL, nBL, largestL, nAL);
     1464    printf("          %d upper, smallest %.18g, %d below - largest %.18g, %d above\n",
     1465      nU, smallestU, nBU, largestU, nAU);
     1466  }
     1467#endif
     1468  return directionNorm;
    14921469}
    14931470/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
    1494 void
    1495 ClpPredictorCorrector::solveSystem(CoinWorkDouble * region1, CoinWorkDouble * region2,
    1496                                    const CoinWorkDouble * region1In, const CoinWorkDouble * region2In,
    1497                                    const CoinWorkDouble * saveRegion1, const CoinWorkDouble * saveRegion2,
    1498                                    bool gentleRefine)
     1471void ClpPredictorCorrector::solveSystem(CoinWorkDouble *region1, CoinWorkDouble *region2,
     1472  const CoinWorkDouble *region1In, const CoinWorkDouble *region2In,
     1473  const CoinWorkDouble *saveRegion1, const CoinWorkDouble *saveRegion2,
     1474  bool gentleRefine)
    14991475{
    1500      int iRow;
    1501      int numberTotal = numberRows_ + numberColumns_;
    1502      if (region2In) {
    1503           // normal
    1504           for (iRow = 0; iRow < numberRows_; iRow++)
    1505                region2[iRow] = region2In[iRow];
    1506      } else {
    1507           // initial solution - (diagonal is 1 or 0)
    1508           CoinZeroN(region2, numberRows_);
    1509      }
    1510      int iColumn;
    1511      if (cholesky_->type() < 20) {
    1512           // not KKT
    1513           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1514                region1[iColumn] = region1In[iColumn] * diagonal_[iColumn];
    1515           multiplyAdd(region1 + numberColumns_, numberRows_, -1.0, region2, 1.0);
    1516           matrix_->times(1.0, region1, region2);
    1517           CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
    1518           CoinWorkDouble scale = 1.0;
    1519           CoinWorkDouble unscale = 1.0;
    1520           if (maximumRHS > 1.0e-30) {
    1521                if (maximumRHS <= 0.5) {
    1522                     CoinWorkDouble factor = 2.0;
    1523                     while (maximumRHS <= 0.5) {
    1524                          maximumRHS *= factor;
    1525                          scale *= factor;
    1526                     } /* endwhile */
    1527                } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
    1528                     CoinWorkDouble factor = 0.5;
    1529                     while (maximumRHS >= 2.0) {
    1530                          maximumRHS *= factor;
    1531                          scale *= factor;
    1532                     } /* endwhile */
    1533                }
    1534                unscale = diagonalScaleFactor_ / scale;
    1535           } else {
    1536                //effectively zero
    1537                scale = 0.0;
    1538                unscale = 0.0;
    1539           }
    1540           multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
    1541           cholesky_->solve(region2);
    1542           multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
    1543           multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns_, 0.0);
    1544           CoinZeroN(region1, numberColumns_);
    1545           matrix_->transposeTimes(1.0, region2, region1);
    1546           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1547                region1[iColumn] = (region1[iColumn] - region1In[iColumn]) * diagonal_[iColumn];
    1548      } else {
    1549           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1550                region1[iColumn] = region1In[iColumn];
    1551           cholesky_->solveKKT(region1, region2, diagonal_, diagonalScaleFactor_);
    1552      }
    1553      if (saveRegion2) {
    1554           //refine?
    1555           CoinWorkDouble scaleX = 1.0;
    1556           if (gentleRefine)
    1557                scaleX = 0.8;
    1558           multiplyAdd(saveRegion2, numberRows_, 1.0, region2, scaleX);
    1559           assert (saveRegion1);
    1560           multiplyAdd(saveRegion1, numberTotal, 1.0, region1, scaleX);
    1561      }
     1476  int iRow;
     1477  int numberTotal = numberRows_ + numberColumns_;
     1478  if (region2In) {
     1479    // normal
     1480    for (iRow = 0; iRow < numberRows_; iRow++)
     1481      region2[iRow] = region2In[iRow];
     1482  } else {
     1483    // initial solution - (diagonal is 1 or 0)
     1484    CoinZeroN(region2, numberRows_);
     1485  }
     1486  int iColumn;
     1487  if (cholesky_->type() < 20) {
     1488    // not KKT
     1489    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1490      region1[iColumn] = region1In[iColumn] * diagonal_[iColumn];
     1491    multiplyAdd(region1 + numberColumns_, numberRows_, -1.0, region2, 1.0);
     1492    matrix_->times(1.0, region1, region2);
     1493    CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
     1494    CoinWorkDouble scale = 1.0;
     1495    CoinWorkDouble unscale = 1.0;
     1496    if (maximumRHS > 1.0e-30) {
     1497      if (maximumRHS <= 0.5) {
     1498        CoinWorkDouble factor = 2.0;
     1499        while (maximumRHS <= 0.5) {
     1500          maximumRHS *= factor;
     1501          scale *= factor;
     1502        } /* endwhile */
     1503      } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
     1504        CoinWorkDouble factor = 0.5;
     1505        while (maximumRHS >= 2.0) {
     1506          maximumRHS *= factor;
     1507          scale *= factor;
     1508        } /* endwhile */
     1509      }
     1510      unscale = diagonalScaleFactor_ / scale;
     1511    } else {
     1512      //effectively zero
     1513      scale = 0.0;
     1514      unscale = 0.0;
     1515    }
     1516    multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
     1517    cholesky_->solve(region2);
     1518    multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
     1519    multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns_, 0.0);
     1520    CoinZeroN(region1, numberColumns_);
     1521    matrix_->transposeTimes(1.0, region2, region1);
     1522    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1523      region1[iColumn] = (region1[iColumn] - region1In[iColumn]) * diagonal_[iColumn];
     1524  } else {
     1525    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1526      region1[iColumn] = region1In[iColumn];
     1527    cholesky_->solveKKT(region1, region2, diagonal_, diagonalScaleFactor_);
     1528  }
     1529  if (saveRegion2) {
     1530    //refine?
     1531    CoinWorkDouble scaleX = 1.0;
     1532    if (gentleRefine)
     1533      scaleX = 0.8;
     1534    multiplyAdd(saveRegion2, numberRows_, 1.0, region2, scaleX);
     1535    assert(saveRegion1);
     1536    multiplyAdd(saveRegion1, numberTotal, 1.0, region1, scaleX);
     1537  }
    15621538}
    15631539// findDirectionVector.
    15641540CoinWorkDouble ClpPredictorCorrector::findDirectionVector(const int phase)
    15651541{
    1566      CoinWorkDouble projectionTolerance = projectionTolerance_;
    1567      //temporary
    1568      //projectionTolerance=1.0e-15;
    1569      CoinWorkDouble errorCheck = 0.9 * maximumRHSError_ / solutionNorm_;
    1570      if (errorCheck > primalTolerance()) {
    1571           if (errorCheck < projectionTolerance) {
    1572                projectionTolerance = errorCheck;
    1573           }
    1574      } else {
    1575           if (primalTolerance() < projectionTolerance) {
    1576                projectionTolerance = primalTolerance();
    1577           }
    1578      }
    1579      CoinWorkDouble * newError = new CoinWorkDouble [numberRows_];
    1580      int numberTotal = numberRows_ + numberColumns_;
    1581      //if flagged then entries zero so can do
    1582      // For KKT separate out
    1583      CoinWorkDouble * region1Save = NULL; //for refinement
    1584      int iColumn;
    1585      if (cholesky_->type() < 20) {
    1586           int iColumn;
    1587           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1588                deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
    1589           multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
    1590           matrix_->times(1.0, deltaX_, deltaY_);
    1591      } else {
    1592           // regions in will be workArray and newError
    1593           // regions out deltaX_ and deltaY_
    1594           multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, newError, 0.0);
    1595           matrix_->times(-1.0, solution_, newError);
    1596           // This is inefficient but just for now get values which will be in deltay
    1597           int iColumn;
    1598           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1599                deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
    1600           multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
    1601           matrix_->times(1.0, deltaX_, deltaY_);
    1602      }
    1603      bool goodSolve = false;
    1604      CoinWorkDouble * regionSave = NULL; //for refinement
    1605      int numberTries = 0;
    1606      CoinWorkDouble relativeError = COIN_DBL_MAX;
    1607      CoinWorkDouble tryError = 1.0e31;
    1608      CoinWorkDouble saveMaximum = 0.0;
    1609      double firstError = 0.0;
    1610      double lastError2 = 0.0;
    1611      while (!goodSolve && numberTries < 30) {
    1612           CoinWorkDouble lastError = relativeError;
    1613           goodSolve = true;
    1614           CoinWorkDouble maximumRHS;
    1615           maximumRHS = CoinMax(maximumAbsElement(deltaY_, numberRows_), 1.0e-12);
    1616           if (!numberTries)
    1617                saveMaximum = maximumRHS;
    1618           if (cholesky_->type() < 20) {
    1619                // no kkt
    1620                CoinWorkDouble scale = 1.0;
    1621                CoinWorkDouble unscale = 1.0;
    1622                if (maximumRHS > 1.0e-30) {
    1623                     if (maximumRHS <= 0.5) {
    1624                          CoinWorkDouble factor = 2.0;
    1625                          while (maximumRHS <= 0.5) {
    1626                               maximumRHS *= factor;
    1627                               scale *= factor;
    1628                          } /* endwhile */
    1629                     } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
    1630                          CoinWorkDouble factor = 0.5;
    1631                          while (maximumRHS >= 2.0) {
    1632                               maximumRHS *= factor;
    1633                               scale *= factor;
    1634                          } /* endwhile */
    1635                     }
    1636                     unscale = diagonalScaleFactor_ / scale;
    1637                } else {
    1638                     //effectively zero
    1639                     scale = 0.0;
    1640                     unscale = 0.0;
    1641                }
    1642                //printf("--putting scales to 1.0\n");
    1643                //scale=1.0;
    1644                //unscale=1.0;
    1645                multiplyAdd(NULL, numberRows_, 0.0, deltaY_, scale);
    1646                cholesky_->solve(deltaY_);
    1647                multiplyAdd(NULL, numberRows_, 0.0, deltaY_, unscale);
     1542  CoinWorkDouble projectionTolerance = projectionTolerance_;
     1543  //temporary
     1544  //projectionTolerance=1.0e-15;
     1545  CoinWorkDouble errorCheck = 0.9 * maximumRHSError_ / solutionNorm_;
     1546  if (errorCheck > primalTolerance()) {
     1547    if (errorCheck < projectionTolerance) {
     1548      projectionTolerance = errorCheck;
     1549    }
     1550  } else {
     1551    if (primalTolerance() < projectionTolerance) {
     1552      projectionTolerance = primalTolerance();
     1553    }
     1554  }
     1555  CoinWorkDouble *newError = new CoinWorkDouble[numberRows_];
     1556  int numberTotal = numberRows_ + numberColumns_;
     1557  //if flagged then entries zero so can do
     1558  // For KKT separate out
     1559  CoinWorkDouble *region1Save = NULL; //for refinement
     1560  int iColumn;
     1561  if (cholesky_->type() < 20) {
     1562    int iColumn;
     1563    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1564      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
     1565    multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
     1566    matrix_->times(1.0, deltaX_, deltaY_);
     1567  } else {
     1568    // regions in will be workArray and newError
     1569    // regions out deltaX_ and deltaY_
     1570    multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, newError, 0.0);
     1571    matrix_->times(-1.0, solution_, newError);
     1572    // This is inefficient but just for now get values which will be in deltay
     1573    int iColumn;
     1574    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1575      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
     1576    multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
     1577    matrix_->times(1.0, deltaX_, deltaY_);
     1578  }
     1579  bool goodSolve = false;
     1580  CoinWorkDouble *regionSave = NULL; //for refinement
     1581  int numberTries = 0;
     1582  CoinWorkDouble relativeError = COIN_DBL_MAX;
     1583  CoinWorkDouble tryError = 1.0e31;
     1584  CoinWorkDouble saveMaximum = 0.0;
     1585  double firstError = 0.0;
     1586  double lastError2 = 0.0;
     1587  while (!goodSolve && numberTries < 30) {
     1588    CoinWorkDouble lastError = relativeError;
     1589    goodSolve = true;
     1590    CoinWorkDouble maximumRHS;
     1591    maximumRHS = CoinMax(maximumAbsElement(deltaY_, numberRows_), 1.0e-12);
     1592    if (!numberTries)
     1593      saveMaximum = maximumRHS;
     1594    if (cholesky_->type() < 20) {
     1595      // no kkt
     1596      CoinWorkDouble scale = 1.0;
     1597      CoinWorkDouble unscale = 1.0;
     1598      if (maximumRHS > 1.0e-30) {
     1599        if (maximumRHS <= 0.5) {
     1600          CoinWorkDouble factor = 2.0;
     1601          while (maximumRHS <= 0.5) {
     1602            maximumRHS *= factor;
     1603            scale *= factor;
     1604          } /* endwhile */
     1605        } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
     1606          CoinWorkDouble factor = 0.5;
     1607          while (maximumRHS >= 2.0) {
     1608            maximumRHS *= factor;
     1609            scale *= factor;
     1610          } /* endwhile */
     1611        }
     1612        unscale = diagonalScaleFactor_ / scale;
     1613      } else {
     1614        //effectively zero
     1615        scale = 0.0;
     1616        unscale = 0.0;
     1617      }
     1618      //printf("--putting scales to 1.0\n");
     1619      //scale=1.0;
     1620      //unscale=1.0;
     1621      multiplyAdd(NULL, numberRows_, 0.0, deltaY_, scale);
     1622      cholesky_->solve(deltaY_);
     1623      multiplyAdd(NULL, numberRows_, 0.0, deltaY_, unscale);
    16481624#if 0
    16491625               {
     
    16541630               exit(66);
    16551631#endif
    1656                if (numberTries) {
    1657                     //refine?
    1658                     CoinWorkDouble scaleX = 1.0;
    1659                     if (lastError > 1.0e-5)
    1660                          scaleX = 0.8;
    1661                     multiplyAdd(regionSave, numberRows_, 1.0, deltaY_, scaleX);
    1662                }
    1663                //CoinZeroN(newError,numberRows_);
    1664                multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
    1665                CoinZeroN(deltaX_, numberColumns_);
    1666                matrix_->transposeTimes(1.0, deltaY_, deltaX_);
    1667                //if flagged then entries zero so can do
    1668                for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1669                     deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
    1670                                        - workArray_[iColumn];
     1632      if (numberTries) {
     1633        //refine?
     1634        CoinWorkDouble scaleX = 1.0;
     1635        if (lastError > 1.0e-5)
     1636          scaleX = 0.8;
     1637        multiplyAdd(regionSave, numberRows_, 1.0, deltaY_, scaleX);
     1638      }
     1639      //CoinZeroN(newError,numberRows_);
     1640      multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
     1641      CoinZeroN(deltaX_, numberColumns_);
     1642      matrix_->transposeTimes(1.0, deltaY_, deltaX_);
     1643      //if flagged then entries zero so can do
     1644      for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1645        deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
     1646          - workArray_[iColumn];
     1647    } else {
     1648      // KKT
     1649      solveSystem(deltaX_, deltaY_,
     1650        workArray_, newError, region1Save, regionSave, lastError > 1.0e-5);
     1651    }
     1652    multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, newError, 0.0);
     1653    matrix_->times(1.0, deltaX_, newError);
     1654    numberTries++;
     1655
     1656    //now add in old Ax - doing extra checking
     1657    CoinWorkDouble maximumRHSError = 0.0;
     1658    CoinWorkDouble maximumRHSChange = 0.0;
     1659    int iRow;
     1660    char *dropped = cholesky_->rowsDropped();
     1661    for (iRow = 0; iRow < numberRows_; iRow++) {
     1662      if (!dropped[iRow]) {
     1663        CoinWorkDouble newValue = newError[iRow];
     1664        CoinWorkDouble oldValue = errorRegion_[iRow];
     1665        //severity of errors depend on signs
     1666        //**later                                                             */
     1667        if (CoinAbs(newValue) > maximumRHSChange) {
     1668          maximumRHSChange = CoinAbs(newValue);
     1669        }
     1670        CoinWorkDouble result = newValue + oldValue;
     1671        if (CoinAbs(result) > maximumRHSError) {
     1672          maximumRHSError = CoinAbs(result);
     1673        }
     1674        newError[iRow] = result;
     1675      } else {
     1676        CoinWorkDouble newValue = newError[iRow];
     1677        CoinWorkDouble oldValue = errorRegion_[iRow];
     1678        if (CoinAbs(newValue) > maximumRHSChange) {
     1679          maximumRHSChange = CoinAbs(newValue);
     1680        }
     1681        CoinWorkDouble result = newValue + oldValue;
     1682        newError[iRow] = result;
     1683        //newError[iRow]=0.0;
     1684        //assert(deltaY_[iRow]==0.0);
     1685        deltaY_[iRow] = 0.0;
     1686      }
     1687    }
     1688    relativeError = maximumRHSError / solutionNorm_;
     1689    relativeError = maximumRHSError / saveMaximum;
     1690    if (relativeError > tryError)
     1691      relativeError = tryError;
     1692    if (numberTries == 1)
     1693      firstError = relativeError;
     1694    if (relativeError < lastError) {
     1695      lastError2 = relativeError;
     1696      maximumRHSChange_ = maximumRHSChange;
     1697      if (relativeError > projectionTolerance && numberTries <= 3) {
     1698        //try and refine
     1699        goodSolve = false;
     1700      }
     1701      //*** extra test here
     1702      if (!goodSolve) {
     1703        if (!regionSave) {
     1704          regionSave = new CoinWorkDouble[numberRows_];
     1705          if (cholesky_->type() >= 20)
     1706            region1Save = new CoinWorkDouble[numberTotal];
     1707        }
     1708        CoinMemcpyN(deltaY_, numberRows_, regionSave);
     1709        if (cholesky_->type() < 20) {
     1710          // not KKT
     1711          multiplyAdd(newError, numberRows_, -1.0, deltaY_, 0.0);
     1712        } else {
     1713          // KKT
     1714          CoinMemcpyN(deltaX_, numberTotal, region1Save);
     1715          // and back to input region
     1716          CoinMemcpyN(deltaY_, numberRows_, newError);
     1717        }
     1718      }
     1719    } else {
     1720      //std::cout <<" worse residual = "<<relativeError;
     1721      //bring back previous
     1722      relativeError = lastError;
     1723      if (regionSave) {
     1724        CoinMemcpyN(regionSave, numberRows_, deltaY_);
     1725        if (cholesky_->type() < 20) {
     1726          // not KKT
     1727          multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
     1728          CoinZeroN(deltaX_, numberColumns_);
     1729          matrix_->transposeTimes(1.0, deltaY_, deltaX_);
     1730          //if flagged then entries zero so can do
     1731          for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1732            deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
     1733              - workArray_[iColumn];
     1734        } else {
     1735          // KKT
     1736          CoinMemcpyN(region1Save, numberTotal, deltaX_);
     1737        }
     1738      } else {
     1739        // disaster
     1740        CoinFillN(deltaX_, numberTotal, static_cast< CoinWorkDouble >(1.0));
     1741        CoinFillN(deltaY_, numberRows_, static_cast< CoinWorkDouble >(1.0));
     1742        COIN_DETAIL_PRINT(printf("bad cholesky\n"));
     1743      }
     1744    }
     1745  } /* endwhile */
     1746  if (firstError > 1.0e-8 || numberTries > 1) {
     1747    handler_->message(CLP_BARRIER_ACCURACY, messages_)
     1748      << phase << numberTries << static_cast< double >(firstError)
     1749      << static_cast< double >(lastError2)
     1750      << CoinMessageEol;
     1751  }
     1752  delete[] regionSave;
     1753  delete[] region1Save;
     1754  delete[] newError;
     1755  // now rest
     1756  CoinWorkDouble extra = eExtra;
     1757  //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
     1758  //CoinZeroN(deltaW_,numberColumns_);
     1759  //matrix_->transposeTimes(-1.0,deltaY_,deltaW_);
     1760
     1761  for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
     1762    deltaSU_[iColumn] = 0.0;
     1763    deltaSL_[iColumn] = 0.0;
     1764    deltaZ_[iColumn] = 0.0;
     1765    CoinWorkDouble dd = deltaW_[iColumn];
     1766    deltaW_[iColumn] = 0.0;
     1767    if (!flagged(iColumn)) {
     1768      CoinWorkDouble deltaX = deltaX_[iColumn];
     1769      if (lowerBound(iColumn)) {
     1770        CoinWorkDouble zValue = rhsZ_[iColumn];
     1771        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
     1772        CoinWorkDouble slack = lowerSlack_[iColumn] + extra;
     1773        deltaSL_[iColumn] = -rhsL_[iColumn] + deltaX;
     1774        deltaZ_[iColumn] = (gHat - zVec_[iColumn] * deltaX) / slack;
     1775      }
     1776      if (upperBound(iColumn)) {
     1777        CoinWorkDouble wValue = rhsW_[iColumn];
     1778        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
     1779        CoinWorkDouble slack = upperSlack_[iColumn] + extra;
     1780        deltaSU_[iColumn] = rhsU_[iColumn] - deltaX;
     1781        deltaW_[iColumn] = (hHat + wVec_[iColumn] * deltaX) / slack;
     1782      }
     1783      if (0) {
     1784        // different way of calculating
     1785        CoinWorkDouble gamma2 = gamma_ * gamma_;
     1786        CoinWorkDouble dZ = 0.0;
     1787        CoinWorkDouble dW = 0.0;
     1788        CoinWorkDouble zValue = rhsZ_[iColumn];
     1789        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
     1790        CoinWorkDouble slackL = lowerSlack_[iColumn] + extra;
     1791        CoinWorkDouble wValue = rhsW_[iColumn];
     1792        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
     1793        CoinWorkDouble slackU = upperSlack_[iColumn] + extra;
     1794        CoinWorkDouble q = rhsC_[iColumn] + gamma2 * deltaX + dd;
     1795        if (primalR_)
     1796          q += deltaX * primalR_[iColumn];
     1797        dW = (gHat + hHat - slackL * q + (wValue - zValue) * deltaX) / (slackL + slackU);
     1798        dZ = dW + q;
     1799        //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
     1800        //deltaW_[iColumn],dZ,dW);
     1801        if (lowerBound(iColumn)) {
     1802          if (upperBound(iColumn)) {
     1803            //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
     1804            //deltaW_[iColumn],dZ,dW);
     1805            deltaW_[iColumn] = dW;
     1806            deltaZ_[iColumn] = dZ;
    16711807          } else {
    1672                // KKT
    1673                solveSystem(deltaX_, deltaY_,
    1674                            workArray_, newError, region1Save, regionSave, lastError > 1.0e-5);
    1675           }
    1676           multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, newError, 0.0);
    1677           matrix_->times(1.0, deltaX_, newError);
    1678           numberTries++;
    1679 
    1680           //now add in old Ax - doing extra checking
    1681           CoinWorkDouble maximumRHSError = 0.0;
    1682           CoinWorkDouble maximumRHSChange = 0.0;
    1683           int iRow;
    1684           char * dropped = cholesky_->rowsDropped();
    1685           for (iRow = 0; iRow < numberRows_; iRow++) {
    1686                if (!dropped[iRow]) {
    1687                     CoinWorkDouble newValue = newError[iRow];
    1688                     CoinWorkDouble oldValue = errorRegion_[iRow];
    1689                     //severity of errors depend on signs
    1690                     //**later                                                             */
    1691                     if (CoinAbs(newValue) > maximumRHSChange) {
    1692                          maximumRHSChange = CoinAbs(newValue);
    1693                     }
    1694                     CoinWorkDouble result = newValue + oldValue;
    1695                     if (CoinAbs(result) > maximumRHSError) {
    1696                          maximumRHSError = CoinAbs(result);
    1697                     }
    1698                     newError[iRow] = result;
    1699                } else {
    1700                     CoinWorkDouble newValue = newError[iRow];
    1701                     CoinWorkDouble oldValue = errorRegion_[iRow];
    1702                     if (CoinAbs(newValue) > maximumRHSChange) {
    1703                          maximumRHSChange = CoinAbs(newValue);
    1704                     }
    1705                     CoinWorkDouble result = newValue + oldValue;
    1706                     newError[iRow] = result;
    1707                     //newError[iRow]=0.0;
    1708                     //assert(deltaY_[iRow]==0.0);
    1709                     deltaY_[iRow] = 0.0;
    1710                }
    1711           }
    1712           relativeError = maximumRHSError / solutionNorm_;
    1713           relativeError = maximumRHSError / saveMaximum;
    1714           if (relativeError > tryError)
    1715                relativeError = tryError;
    1716           if (numberTries == 1)
    1717                firstError = relativeError;
    1718           if (relativeError < lastError) {
    1719                lastError2 = relativeError;
    1720                maximumRHSChange_ = maximumRHSChange;
    1721                if (relativeError > projectionTolerance && numberTries <= 3) {
    1722                     //try and refine
    1723                     goodSolve = false;
    1724                }
    1725                //*** extra test here
    1726                if (!goodSolve) {
    1727                     if (!regionSave) {
    1728                          regionSave = new CoinWorkDouble [numberRows_];
    1729                          if (cholesky_->type() >= 20)
    1730                               region1Save = new CoinWorkDouble [numberTotal];
    1731                     }
    1732                     CoinMemcpyN(deltaY_, numberRows_, regionSave);
    1733                     if (cholesky_->type() < 20) {
    1734                          // not KKT
    1735                          multiplyAdd(newError, numberRows_, -1.0, deltaY_, 0.0);
    1736                     } else {
    1737                          // KKT
    1738                          CoinMemcpyN(deltaX_, numberTotal, region1Save);
    1739                          // and back to input region
    1740                          CoinMemcpyN(deltaY_, numberRows_, newError);
    1741                     }
    1742                }
    1743           } else {
    1744                //std::cout <<" worse residual = "<<relativeError;
    1745                //bring back previous
    1746                relativeError = lastError;
    1747                if (regionSave) {
    1748                     CoinMemcpyN(regionSave, numberRows_, deltaY_);
    1749                     if (cholesky_->type() < 20) {
    1750                          // not KKT
    1751                          multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
    1752                          CoinZeroN(deltaX_, numberColumns_);
    1753                          matrix_->transposeTimes(1.0, deltaY_, deltaX_);
    1754                          //if flagged then entries zero so can do
    1755                          for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1756                               deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
    1757                                                  - workArray_[iColumn];
    1758                     } else {
    1759                          // KKT
    1760                          CoinMemcpyN(region1Save, numberTotal, deltaX_);
    1761                     }
    1762                } else {
    1763                     // disaster
    1764                     CoinFillN(deltaX_, numberTotal, static_cast<CoinWorkDouble>(1.0));
    1765                     CoinFillN(deltaY_, numberRows_, static_cast<CoinWorkDouble>(1.0));
    1766                     COIN_DETAIL_PRINT(printf("bad cholesky\n"));
    1767                }
    1768           }
    1769      } /* endwhile */
    1770      if (firstError > 1.0e-8 || numberTries > 1) {
    1771           handler_->message(CLP_BARRIER_ACCURACY, messages_)
    1772                     << phase << numberTries << static_cast<double>(firstError)
    1773                     << static_cast<double>(lastError2)
    1774                     << CoinMessageEol;
    1775      }
    1776      delete [] regionSave;
    1777      delete [] region1Save;
    1778      delete [] newError;
    1779      // now rest
    1780      CoinWorkDouble extra = eExtra;
    1781      //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
    1782      //CoinZeroN(deltaW_,numberColumns_);
    1783      //matrix_->transposeTimes(-1.0,deltaY_,deltaW_);
    1784 
    1785      for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    1786           deltaSU_[iColumn] = 0.0;
    1787           deltaSL_[iColumn] = 0.0;
    1788           deltaZ_[iColumn] = 0.0;
    1789           CoinWorkDouble dd = deltaW_[iColumn];
    1790           deltaW_[iColumn] = 0.0;
    1791           if (!flagged(iColumn)) {
    1792                CoinWorkDouble deltaX = deltaX_[iColumn];
    1793                if (lowerBound(iColumn)) {
    1794                     CoinWorkDouble zValue = rhsZ_[iColumn];
    1795                     CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
    1796                     CoinWorkDouble slack = lowerSlack_[iColumn] + extra;
    1797                     deltaSL_[iColumn] = -rhsL_[iColumn] + deltaX;
    1798                     deltaZ_[iColumn] = (gHat - zVec_[iColumn] * deltaX) / slack;
    1799                }
    1800                if (upperBound(iColumn)) {
    1801                     CoinWorkDouble wValue = rhsW_[iColumn];
    1802                     CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
    1803                     CoinWorkDouble slack = upperSlack_[iColumn] + extra;
    1804                     deltaSU_[iColumn] = rhsU_[iColumn] - deltaX;
    1805                     deltaW_[iColumn] = (hHat + wVec_[iColumn] * deltaX) / slack;
    1806                }
    1807                if (0) {
    1808                     // different way of calculating
    1809                     CoinWorkDouble gamma2 = gamma_ * gamma_;
    1810                     CoinWorkDouble dZ = 0.0;
    1811                     CoinWorkDouble dW = 0.0;
    1812                     CoinWorkDouble zValue = rhsZ_[iColumn];
    1813                     CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
    1814                     CoinWorkDouble slackL = lowerSlack_[iColumn] + extra;
    1815                     CoinWorkDouble wValue = rhsW_[iColumn];
    1816                     CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
    1817                     CoinWorkDouble slackU = upperSlack_[iColumn] + extra;
    1818                     CoinWorkDouble q = rhsC_[iColumn] + gamma2 * deltaX + dd;
    1819                     if (primalR_)
    1820                          q += deltaX * primalR_[iColumn];
    1821                     dW = (gHat + hHat - slackL * q + (wValue - zValue) * deltaX) / (slackL + slackU);
    1822                     dZ = dW + q;
    1823                     //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
    1824                     //deltaW_[iColumn],dZ,dW);
    1825                     if (lowerBound(iColumn)) {
    1826                          if (upperBound(iColumn)) {
    1827                               //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
    1828                               //deltaW_[iColumn],dZ,dW);
    1829                               deltaW_[iColumn] = dW;
    1830                               deltaZ_[iColumn] = dZ;
    1831                          } else {
    1832                               // just lower
    1833                               //printf("L %d old %g new %g\n",iColumn,deltaZ_[iColumn],
    1834                               //dZ);
    1835                          }
    1836                     } else {
    1837                          assert (upperBound(iColumn));
    1838                          //printf("U %d old %g new %g\n",iColumn,deltaW_[iColumn],
    1839                          //dW);
    1840                     }
    1841                }
    1842           }
    1843      }
     1808            // just lower
     1809            //printf("L %d old %g new %g\n",iColumn,deltaZ_[iColumn],
     1810            //dZ);
     1811          }
     1812        } else {
     1813          assert(upperBound(iColumn));
     1814          //printf("U %d old %g new %g\n",iColumn,deltaW_[iColumn],
     1815          //dW);
     1816        }
     1817      }
     1818    }
     1819  }
    18441820#if 0
    18451821     CoinWorkDouble * check = new CoinWorkDouble[numberTotal];
     
    18701846     delete [] check;
    18711847#endif
    1872      return relativeError;
     1848  return relativeError;
    18731849}
    18741850// createSolution.  Creates solution from scratch
    18751851int ClpPredictorCorrector::createSolution()
    18761852{
    1877      int numberTotal = numberRows_ + numberColumns_;
    1878      int iColumn;
    1879      CoinWorkDouble tolerance = primalTolerance();
    1880      // See if quadratic objective
     1853  int numberTotal = numberRows_ + numberColumns_;
     1854  int iColumn;
     1855  CoinWorkDouble tolerance = primalTolerance();
     1856  // See if quadratic objective
    18811857#ifndef NO_RTTI
    1882      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1858  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
    18831859#else
    1884      ClpQuadraticObjective * quadraticObj = NULL;
    1885      if (objective_->type() == 2)
    1886           quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    1887 #endif
    1888      if (!quadraticObj) {
    1889           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    1890                if (upper_[iColumn] - lower_[iColumn] > tolerance)
    1891                     clearFixed(iColumn);
    1892                else
    1893                     setFixed(iColumn);
    1894           }
    1895      } else {
    1896           // try leaving fixed
    1897           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    1898                clearFixed(iColumn);
    1899      }
     1860  ClpQuadraticObjective *quadraticObj = NULL;
     1861  if (objective_->type() == 2)
     1862    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     1863#endif
     1864  if (!quadraticObj) {
     1865    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     1866      if (upper_[iColumn] - lower_[iColumn] > tolerance)
     1867        clearFixed(iColumn);
     1868      else
     1869        setFixed(iColumn);
     1870    }
     1871  } else {
     1872    // try leaving fixed
     1873    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     1874      clearFixed(iColumn);
     1875  }
    19001876
    1901      CoinWorkDouble maximumObjective = 0.0;
    1902      CoinWorkDouble objectiveNorm2 = 0.0;
    1903      getNorms(cost_, numberTotal, maximumObjective, objectiveNorm2);
    1904      if (!maximumObjective) {
    1905           maximumObjective = 1.0; // objective all zero
    1906      }
    1907      objectiveNorm2 = CoinSqrt(objectiveNorm2) / static_cast<CoinWorkDouble> (numberTotal);
    1908      objectiveNorm_ = maximumObjective;
    1909      scaleFactor_ = 1.0;
    1910      if (maximumObjective > 0.0) {
    1911           if (maximumObjective < 1.0) {
    1912                scaleFactor_ = maximumObjective;
    1913           } else if (maximumObjective > 1.0e4) {
    1914                scaleFactor_ = maximumObjective / 1.0e4;
    1915           }
    1916      }
    1917      if (scaleFactor_ != 1.0) {
    1918           objectiveNorm2 *= scaleFactor_;
    1919           multiplyAdd(NULL, numberTotal, 0.0, cost_, 1.0 / scaleFactor_);
    1920           objectiveNorm_ = maximumObjective / scaleFactor_;
    1921      }
    1922      // See if quadratic objective
    1923      if (quadraticObj) {
    1924           // If scaled then really scale matrix
    1925           CoinWorkDouble scaleFactor =
    1926                scaleFactor_ * optimizationDirection_ * objectiveScale_ *
    1927                rhsScale_;
    1928           if ((scalingFlag_ > 0 && rowScale_) || scaleFactor != 1.0) {
    1929                CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    1930                const int * columnQuadratic = quadratic->getIndices();
    1931                const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    1932                const int * columnQuadraticLength = quadratic->getVectorLengths();
    1933                double * quadraticElement = quadratic->getMutableElements();
    1934                int numberColumns = quadratic->getNumCols();
    1935                CoinWorkDouble scale = 1.0 / scaleFactor;
    1936                if (scalingFlag_ > 0 && rowScale_) {
    1937                     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1938                          CoinWorkDouble scaleI = columnScale_[iColumn] * scale;
    1939                          for (CoinBigIndex j = columnQuadraticStart[iColumn];
    1940                                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    1941                               int jColumn = columnQuadratic[j];
    1942                               CoinWorkDouble scaleJ = columnScale_[jColumn];
    1943                               quadraticElement[j] *= scaleI * scaleJ;
    1944                               objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
    1945                          }
    1946                     }
    1947                } else {
    1948                     // not scaled
    1949                     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1950                          for (CoinBigIndex j = columnQuadraticStart[iColumn];
    1951                                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    1952                               quadraticElement[j] *= scale;
    1953                               objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
    1954                          }
    1955                     }
    1956                }
    1957           }
    1958      }
    1959      baseObjectiveNorm_ = objectiveNorm_;
    1960      //accumulate fixed in dj region (as spare)
    1961      //accumulate primal solution in primal region
    1962      //DZ in lowerDual
    1963      //DW in upperDual
    1964      CoinWorkDouble infiniteCheck = 1.0e40;
    1965      //CoinWorkDouble     fakeCheck=1.0e10;
    1966      //use deltaX region for work region
    1967      for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    1968           CoinWorkDouble primalValue = solution_[iColumn];
    1969           clearFlagged(iColumn);
    1970           clearFixedOrFree(iColumn);
    1971           clearLowerBound(iColumn);
    1972           clearUpperBound(iColumn);
    1973           clearFakeLower(iColumn);
    1974           clearFakeUpper(iColumn);
    1975           if (!fixed(iColumn)) {
    1976                dj_[iColumn] = 0.0;
    1977                diagonal_[iColumn] = 1.0;
    1978                deltaX_[iColumn] = 1.0;
    1979                CoinWorkDouble lowerValue = lower_[iColumn];
    1980                CoinWorkDouble upperValue = upper_[iColumn];
    1981                if (lowerValue > -infiniteCheck) {
    1982                     if (upperValue < infiniteCheck) {
    1983                          //upper and lower bounds
    1984                          setLowerBound(iColumn);
    1985                          setUpperBound(iColumn);
    1986                          if (lowerValue >= 0.0) {
    1987                               solution_[iColumn] = lowerValue;
    1988                          } else if (upperValue <= 0.0) {
    1989                               solution_[iColumn] = upperValue;
    1990                          } else {
    1991                               solution_[iColumn] = 0.0;
    1992                          }
    1993                     } else {
    1994                          //just lower bound
    1995                          setLowerBound(iColumn);
    1996                          if (lowerValue >= 0.0) {
    1997                               solution_[iColumn] = lowerValue;
    1998                          } else {
    1999                               solution_[iColumn] = 0.0;
    2000                          }
    2001                     }
    2002                } else {
    2003                     if (upperValue < infiniteCheck) {
    2004                          //just upper bound
    2005                          setUpperBound(iColumn);
    2006                          if (upperValue <= 0.0) {
    2007                               solution_[iColumn] = upperValue;
    2008                          } else {
    2009                               solution_[iColumn] = 0.0;
    2010                          }
    2011                     } else {
    2012                          //free
    2013                          setFixedOrFree(iColumn);
    2014                          solution_[iColumn] = 0.0;
    2015                          //std::cout<<" free "<<i<<std::endl;
    2016                     }
    2017                }
     1877  CoinWorkDouble maximumObjective = 0.0;
     1878  CoinWorkDouble objectiveNorm2 = 0.0;
     1879  getNorms(cost_, numberTotal, maximumObjective, objectiveNorm2);
     1880  if (!maximumObjective) {
     1881    maximumObjective = 1.0; // objective all zero
     1882  }
     1883  objectiveNorm2 = CoinSqrt(objectiveNorm2) / static_cast< CoinWorkDouble >(numberTotal);
     1884  objectiveNorm_ = maximumObjective;
     1885  scaleFactor_ = 1.0;
     1886  if (maximumObjective > 0.0) {
     1887    if (maximumObjective < 1.0) {
     1888      scaleFactor_ = maximumObjective;
     1889    } else if (maximumObjective > 1.0e4) {
     1890      scaleFactor_ = maximumObjective / 1.0e4;
     1891    }
     1892  }
     1893  if (scaleFactor_ != 1.0) {
     1894    objectiveNorm2 *= scaleFactor_;
     1895    multiplyAdd(NULL, numberTotal, 0.0, cost_, 1.0 / scaleFactor_);
     1896    objectiveNorm_ = maximumObjective / scaleFactor_;
     1897  }
     1898  // See if quadratic objective
     1899  if (quadraticObj) {
     1900    // If scaled then really scale matrix
     1901    CoinWorkDouble scaleFactor = scaleFactor_ * optimizationDirection_ * objectiveScale_ * rhsScale_;
     1902    if ((scalingFlag_ > 0 && rowScale_) || scaleFactor != 1.0) {
     1903      CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
     1904      const int *columnQuadratic = quadratic->getIndices();
     1905      const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     1906      const int *columnQuadraticLength = quadratic->getVectorLengths();
     1907      double *quadraticElement = quadratic->getMutableElements();
     1908      int numberColumns = quadratic->getNumCols();
     1909      CoinWorkDouble scale = 1.0 / scaleFactor;
     1910      if (scalingFlag_ > 0 && rowScale_) {
     1911        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1912          CoinWorkDouble scaleI = columnScale_[iColumn] * scale;
     1913          for (CoinBigIndex j = columnQuadraticStart[iColumn];
     1914               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     1915            int jColumn = columnQuadratic[j];
     1916            CoinWorkDouble scaleJ = columnScale_[jColumn];
     1917            quadraticElement[j] *= scaleI * scaleJ;
     1918            objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
     1919          }
     1920        }
     1921      } else {
     1922        // not scaled
     1923        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1924          for (CoinBigIndex j = columnQuadraticStart[iColumn];
     1925               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     1926            quadraticElement[j] *= scale;
     1927            objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
     1928          }
     1929        }
     1930      }
     1931    }
     1932  }
     1933  baseObjectiveNorm_ = objectiveNorm_;
     1934  //accumulate fixed in dj region (as spare)
     1935  //accumulate primal solution in primal region
     1936  //DZ in lowerDual
     1937  //DW in upperDual
     1938  CoinWorkDouble infiniteCheck = 1.0e40;
     1939  //CoinWorkDouble     fakeCheck=1.0e10;
     1940  //use deltaX region for work region
     1941  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     1942    CoinWorkDouble primalValue = solution_[iColumn];
     1943    clearFlagged(iColumn);
     1944    clearFixedOrFree(iColumn);
     1945    clearLowerBound(iColumn);
     1946    clearUpperBound(iColumn);
     1947    clearFakeLower(iColumn);
     1948    clearFakeUpper(iColumn);
     1949    if (!fixed(iColumn)) {
     1950      dj_[iColumn] = 0.0;
     1951      diagonal_[iColumn] = 1.0;
     1952      deltaX_[iColumn] = 1.0;
     1953      CoinWorkDouble lowerValue = lower_[iColumn];
     1954      CoinWorkDouble upperValue = upper_[iColumn];
     1955      if (lowerValue > -infiniteCheck) {
     1956        if (upperValue < infiniteCheck) {
     1957          //upper and lower bounds
     1958          setLowerBound(iColumn);
     1959          setUpperBound(iColumn);
     1960          if (lowerValue >= 0.0) {
     1961            solution_[iColumn] = lowerValue;
     1962          } else if (upperValue <= 0.0) {
     1963            solution_[iColumn] = upperValue;
    20181964          } else {
    2019                setFlagged(iColumn);
    2020                setFixedOrFree(iColumn);
    2021                setLowerBound(iColumn);
    2022                setUpperBound(iColumn);
    2023                dj_[iColumn] = primalValue;;
    2024                solution_[iColumn] = lower_[iColumn];
    2025                diagonal_[iColumn] = 0.0;
    2026                deltaX_[iColumn] = 0.0;
    2027           }
    2028      }
    2029      //   modify fixed RHS
    2030      multiplyAdd(dj_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 0.0);
    2031      //   create plausible RHS?
    2032      matrix_->times(-1.0, dj_, rhsFixRegion_);
    2033      multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, errorRegion_, 0.0);
    2034      matrix_->times(-1.0, solution_, errorRegion_);
    2035      rhsNorm_ = maximumAbsElement(errorRegion_, numberRows_);
    2036      if (rhsNorm_ < 1.0) {
    2037           rhsNorm_ = 1.0;
    2038      }
    2039      int * rowsDropped = new int [numberRows_];
    2040      int returnCode = cholesky_->factorize(diagonal_, rowsDropped);
    2041      if (returnCode == -1) {
    2042        COIN_DETAIL_PRINT(printf("Out of memory\n"));
    2043           problemStatus_ = 4;
    2044           return -1;
    2045      }
    2046      if (cholesky_->status()) {
    2047           std::cout << "singular on initial cholesky?" << std::endl;
    2048           cholesky_->resetRowsDropped();
    2049           //cholesky_->factorize(rowDropped_);
    2050           //if (cholesky_->status()) {
    2051           //std::cout << "bad cholesky??? (after retry)" <<std::endl;
    2052           //abort();
    2053           //}
    2054      }
    2055      delete [] rowsDropped;
    2056      if (cholesky_->type() < 20) {
    2057           // not KKT
    2058           cholesky_->solve(errorRegion_);
    2059           //create information for solution
    2060           multiplyAdd(errorRegion_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
    2061           CoinZeroN(deltaX_, numberColumns_);
    2062           matrix_->transposeTimes(1.0, errorRegion_, deltaX_);
    2063      } else {
    2064           // KKT
    2065           // reverse sign on solution
    2066           multiplyAdd(NULL, numberRows_ + numberColumns_, 0.0, solution_, -1.0);
    2067           solveSystem(deltaX_, errorRegion_, solution_, NULL, NULL, NULL, false);
    2068      }
    2069      CoinWorkDouble initialValue = 1.0e2;
    2070      if (rhsNorm_ * 1.0e-2 > initialValue) {
    2071           initialValue = rhsNorm_ * 1.0e-2;
    2072      }
    2073      //initialValue = CoinMax(1.0,rhsNorm_);
    2074      CoinWorkDouble smallestBoundDifference = COIN_DBL_MAX;
    2075      CoinWorkDouble * fakeSolution = deltaX_;
    2076      for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
    2077           if (!flagged(iColumn)) {
    2078                if (lower_[iColumn] - fakeSolution[iColumn] > initialValue) {
    2079                     initialValue = lower_[iColumn] - fakeSolution[iColumn];
    2080                }
    2081                if (fakeSolution[iColumn] - upper_[iColumn] > initialValue) {
    2082                     initialValue = fakeSolution[iColumn] - upper_[iColumn];
    2083                }
    2084                if (upper_[iColumn] - lower_[iColumn] < smallestBoundDifference) {
    2085                     smallestBoundDifference = upper_[iColumn] - lower_[iColumn];
    2086                }
    2087           }
    2088      }
    2089      solutionNorm_ = 1.0e-12;
    2090      handler_->message(CLP_BARRIER_SAFE, messages_)
    2091                << static_cast<double>(initialValue) << static_cast<double>(objectiveNorm_)
    2092                << CoinMessageEol;
    2093      CoinWorkDouble extra = 1.0e-10;
    2094      CoinWorkDouble largeGap = 1.0e15;
    2095      //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_;
    2096      CoinWorkDouble safeObjectiveValue = objectiveNorm_ + 1.0;
    2097      CoinWorkDouble safeFree = 1.0e-1 * initialValue;
    2098      //printf("normal safe dual value of %g, primal value of %g\n",
    2099      // safeObjectiveValue,initialValue);
    2100      //safeObjectiveValue=CoinMax(2.0,1.0e-1*safeObjectiveValue);
    2101      //initialValue=CoinMax(100.0,1.0e-1*initialValue);;
    2102      //printf("temp safe dual value of %g, primal value of %g\n",
    2103      // safeObjectiveValue,initialValue);
    2104      CoinWorkDouble zwLarge = 1.0e2 * initialValue;
    2105      //zwLarge=1.0e40;
    2106      if (cholesky_->choleskyCondition() < 0.0 && cholesky_->type() < 20) {
    2107           // looks bad - play safe
    2108           initialValue *= 10.0;
    2109           safeObjectiveValue *= 10.0;
    2110           safeFree *= 10.0;
    2111      }
    2112      CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
    2113      // First do primal side
    2114      for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
    2115           if (!flagged(iColumn)) {
    2116                CoinWorkDouble lowerValue = lower_[iColumn];
    2117                CoinWorkDouble upperValue = upper_[iColumn];
    2118                CoinWorkDouble newValue;
    2119                CoinWorkDouble setPrimal = initialValue;
    2120                if (quadraticObj) {
    2121                     // perturb primal solution a bit
    2122                     //fakeSolution[iColumn]  *= 0.002*CoinDrand48()+0.999;
    2123                }
    2124                if (lowerBound(iColumn)) {
    2125                     if (upperBound(iColumn)) {
    2126                          //upper and lower bounds
    2127                          if (upperValue - lowerValue > 2.0 * setPrimal) {
    2128                               CoinWorkDouble fakeValue = fakeSolution[iColumn];
    2129                               if (fakeValue < lowerValue + setPrimal) {
    2130                                    fakeValue = lowerValue + setPrimal;
    2131                               }
    2132                               if (fakeValue > upperValue - setPrimal) {
    2133                                    fakeValue = upperValue - setPrimal;
    2134                               }
    2135                               newValue = fakeValue;
    2136                          } else {
    2137                               newValue = 0.5 * (upperValue + lowerValue);
    2138                          }
    2139                     } else {
    2140                          //just lower bound
    2141                          CoinWorkDouble fakeValue = fakeSolution[iColumn];
    2142                          if (fakeValue < lowerValue + setPrimal) {
    2143                               fakeValue = lowerValue + setPrimal;
    2144                          }
    2145                          newValue = fakeValue;
    2146                     }
    2147                } else {
    2148                     if (upperBound(iColumn)) {
    2149                          //just upper bound
    2150                          CoinWorkDouble fakeValue = fakeSolution[iColumn];
    2151                          if (fakeValue > upperValue - setPrimal) {
    2152                               fakeValue = upperValue - setPrimal;
    2153                          }
    2154                          newValue = fakeValue;
    2155                     } else {
    2156                          //free
    2157                          newValue = fakeSolution[iColumn];
    2158                          if (newValue >= 0.0) {
    2159                               if (newValue < safeFree) {
    2160                                    newValue = safeFree;
    2161                               }
    2162                          } else {
    2163                               if (newValue > -safeFree) {
    2164                                    newValue = -safeFree;
    2165                               }
    2166                          }
    2167                     }
    2168                }
    2169                solution_[iColumn] = newValue;
     1965            solution_[iColumn] = 0.0;
     1966          }
     1967        } else {
     1968          //just lower bound
     1969          setLowerBound(iColumn);
     1970          if (lowerValue >= 0.0) {
     1971            solution_[iColumn] = lowerValue;
    21701972          } else {
    2171                // fixed
    2172                lowerSlack_[iColumn] = 0.0;
    2173                upperSlack_[iColumn] = 0.0;
    2174                solution_[iColumn] = lower_[iColumn];
    2175                zVec_[iColumn] = 0.0;
    2176                wVec_[iColumn] = 0.0;
    2177                diagonal_[iColumn] = 0.0;
    2178           }
    2179      }
    2180      solutionNorm_ =  maximumAbsElement(solution_, numberTotal);
    2181      // Set bounds and do dj including quadratic
    2182      largeGap = CoinMax(1.0e7, 1.02 * solutionNorm_);
    2183      CoinPackedMatrix * quadratic = NULL;
    2184      const int * columnQuadratic = NULL;
    2185      const CoinBigIndex * columnQuadraticStart = NULL;
    2186      const int * columnQuadraticLength = NULL;
    2187      const double * quadraticElement = NULL;
    2188      if (quadraticObj) {
    2189           quadratic = quadraticObj->quadraticObjective();
    2190           columnQuadratic = quadratic->getIndices();
    2191           columnQuadraticStart = quadratic->getVectorStarts();
    2192           columnQuadraticLength = quadratic->getVectorLengths();
    2193           quadraticElement = quadratic->getElements();
    2194      }
    2195      CoinWorkDouble quadraticNorm = 0.0;
    2196      for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
    2197           if (!flagged(iColumn)) {
    2198                CoinWorkDouble primalValue = solution_[iColumn];
    2199                CoinWorkDouble lowerValue = lower_[iColumn];
    2200                CoinWorkDouble upperValue = upper_[iColumn];
    2201                // Do dj
    2202                CoinWorkDouble reducedCost = cost_[iColumn];
    2203                if (lowerBound(iColumn)) {
    2204                     reducedCost += linearPerturbation_;
    2205                }
    2206                if (upperBound(iColumn)) {
    2207                     reducedCost -= linearPerturbation_;
    2208                }
    2209                if (quadraticObj && iColumn < numberColumns_) {
    2210                     for (CoinBigIndex j = columnQuadraticStart[iColumn];
    2211                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    2212                          int jColumn = columnQuadratic[j];
    2213                          CoinWorkDouble valueJ = solution_[jColumn];
    2214                          CoinWorkDouble elementValue = quadraticElement[j];
    2215                          reducedCost += valueJ * elementValue;
    2216                     }
    2217                     quadraticNorm = CoinMax(quadraticNorm, CoinAbs(reducedCost));
    2218                }
    2219                dj_[iColumn] = reducedCost;
    2220                if (primalValue > lowerValue + largeGap && primalValue < upperValue - largeGap) {
    2221                     clearFixedOrFree(iColumn);
    2222                     setLowerBound(iColumn);
    2223                     setUpperBound(iColumn);
    2224                     lowerValue = CoinMax(lowerValue, primalValue - largeGap);
    2225                     upperValue = CoinMin(upperValue, primalValue + largeGap);
    2226                     lower_[iColumn] = lowerValue;
    2227                     upper_[iColumn] = upperValue;
    2228                }
    2229           }
    2230      }
    2231      safeObjectiveValue = CoinMax(safeObjectiveValue, quadraticNorm);
    2232      for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
    2233           if (!flagged(iColumn)) {
    2234                CoinWorkDouble primalValue = solution_[iColumn];
    2235                CoinWorkDouble lowerValue = lower_[iColumn];
    2236                CoinWorkDouble upperValue = upper_[iColumn];
    2237                CoinWorkDouble reducedCost = dj_[iColumn];
    2238                CoinWorkDouble low = 0.0;
    2239                CoinWorkDouble high = 0.0;
    2240                if (lowerBound(iColumn)) {
    2241                     if (upperBound(iColumn)) {
    2242                          //upper and lower bounds
    2243                          if (upperValue - lowerValue > 2.0 * initialValue) {
    2244                               low = primalValue - lowerValue;
    2245                               high = upperValue - primalValue;
    2246                          } else {
    2247                               low = initialValue;
    2248                               high = initialValue;
    2249                          }
    2250                          CoinWorkDouble s = low + extra;
    2251                          CoinWorkDouble ratioZ;
    2252                          if (s < zwLarge) {
    2253                               ratioZ = 1.0;
    2254                          } else {
    2255                               ratioZ = CoinSqrt(zwLarge / s);
    2256                          }
    2257                          CoinWorkDouble t = high + extra;
    2258                          CoinWorkDouble ratioT;
    2259                          if (t < zwLarge) {
    2260                               ratioT = 1.0;
    2261                          } else {
    2262                               ratioT = CoinSqrt(zwLarge / t);
    2263                          }
    2264                          //modify s and t
    2265                          if (s > largeGap) {
    2266                               s = largeGap;
    2267                          }
    2268                          if (t > largeGap) {
    2269                               t = largeGap;
    2270                          }
    2271                          //modify if long long way away from bound
    2272                          if (reducedCost >= 0.0) {
    2273                               zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
    2274                               zVec_[iColumn] = CoinMax(reducedCost, safeObjectiveValue * ratioZ);
    2275                               wVec_[iColumn] = safeObjectiveValue * ratioT;
    2276                          } else {
    2277                               zVec_[iColumn] = safeObjectiveValue * ratioZ;
    2278                               wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
    2279                               wVec_[iColumn] = CoinMax(-reducedCost , safeObjectiveValue * ratioT);
    2280                          }
    2281                          CoinWorkDouble gammaTerm = gamma2;
    2282                          if (primalR_)
    2283                               gammaTerm += primalR_[iColumn];
    2284                          diagonal_[iColumn] = (t * s) /
    2285                                               (s * wVec_[iColumn] + t * zVec_[iColumn] + gammaTerm * t * s);
    2286                     } else {
    2287                          //just lower bound
    2288                          low = primalValue - lowerValue;
    2289                          high = 0.0;
    2290                          CoinWorkDouble s = low + extra;
    2291                          CoinWorkDouble ratioZ;
    2292                          if (s < zwLarge) {
    2293                               ratioZ = 1.0;
    2294                          } else {
    2295                               ratioZ = CoinSqrt(zwLarge / s);
    2296                          }
    2297                          //modify s
    2298                          if (s > largeGap) {
    2299                               s = largeGap;
    2300                          }
    2301                          if (reducedCost >= 0.0) {
    2302                               zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
    2303                               zVec_[iColumn] = CoinMax(reducedCost , safeObjectiveValue * ratioZ);
    2304                               wVec_[iColumn] = 0.0;
    2305                          } else {
    2306                               zVec_[iColumn] = safeObjectiveValue * ratioZ;
    2307                               wVec_[iColumn] = 0.0;
    2308                          }
    2309                          CoinWorkDouble gammaTerm = gamma2;
    2310                          if (primalR_)
    2311                               gammaTerm += primalR_[iColumn];
    2312                          diagonal_[iColumn] = s / (zVec_[iColumn] + s * gammaTerm);
    2313                     }
    2314                } else {
    2315                     if (upperBound(iColumn)) {
    2316                          //just upper bound
    2317                          low = 0.0;
    2318                          high = upperValue - primalValue;
    2319                          CoinWorkDouble t = high + extra;
    2320                          CoinWorkDouble ratioT;
    2321                          if (t < zwLarge) {
    2322                               ratioT = 1.0;
    2323                          } else {
    2324                               ratioT = CoinSqrt(zwLarge / t);
    2325                          }
    2326                          //modify t
    2327                          if (t > largeGap) {
    2328                               t = largeGap;
    2329                          }
    2330                          if (reducedCost >= 0.0) {
    2331                               zVec_[iColumn] = 0.0;
    2332                               wVec_[iColumn] = safeObjectiveValue * ratioT;
    2333                          } else {
    2334                               zVec_[iColumn] = 0.0;
    2335                               wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
    2336                               wVec_[iColumn] = CoinMax(-reducedCost , safeObjectiveValue * ratioT);
    2337                          }
    2338                          CoinWorkDouble gammaTerm = gamma2;
    2339                          if (primalR_)
    2340                               gammaTerm += primalR_[iColumn];
    2341                          diagonal_[iColumn] =  t / (wVec_[iColumn] + t * gammaTerm);
    2342                     }
    2343                }
    2344                lowerSlack_[iColumn] = low;
    2345                upperSlack_[iColumn] = high;
    2346           }
    2347      }
     1973            solution_[iColumn] = 0.0;
     1974          }
     1975        }
     1976      } else {
     1977        if (upperValue < infiniteCheck) {
     1978          //just upper bound
     1979          setUpperBound(iColumn);
     1980          if (upperValue <= 0.0) {
     1981            solution_[iColumn] = upperValue;
     1982          } else {
     1983            solution_[iColumn] = 0.0;
     1984          }
     1985        } else {
     1986          //free
     1987          setFixedOrFree(iColumn);
     1988          solution_[iColumn] = 0.0;
     1989          //std::cout<<" free "<<i<<std::endl;
     1990        }
     1991      }
     1992    } else {
     1993      setFlagged(iColumn);
     1994      setFixedOrFree(iColumn);
     1995      setLowerBound(iColumn);
     1996      setUpperBound(iColumn);
     1997      dj_[iColumn] = primalValue;
     1998      ;
     1999      solution_[iColumn] = lower_[iColumn];
     2000      diagonal_[iColumn] = 0.0;
     2001      deltaX_[iColumn] = 0.0;
     2002    }
     2003  }
     2004  //   modify fixed RHS
     2005  multiplyAdd(dj_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 0.0);
     2006  //   create plausible RHS?
     2007  matrix_->times(-1.0, dj_, rhsFixRegion_);
     2008  multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, errorRegion_, 0.0);
     2009  matrix_->times(-1.0, solution_, errorRegion_);
     2010  rhsNorm_ = maximumAbsElement(errorRegion_, numberRows_);
     2011  if (rhsNorm_ < 1.0) {
     2012    rhsNorm_ = 1.0;
     2013  }
     2014  int *rowsDropped = new int[numberRows_];
     2015  int returnCode = cholesky_->factorize(diagonal_, rowsDropped);
     2016  if (returnCode == -1) {
     2017    COIN_DETAIL_PRINT(printf("Out of memory\n"));
     2018    problemStatus_ = 4;
     2019    return -1;
     2020  }
     2021  if (cholesky_->status()) {
     2022    std::cout << "singular on initial cholesky?" << std::endl;
     2023    cholesky_->resetRowsDropped();
     2024    //cholesky_->factorize(rowDropped_);
     2025    //if (cholesky_->status()) {
     2026    //std::cout << "bad cholesky??? (after retry)" <<std::endl;
     2027    //abort();
     2028    //}
     2029  }
     2030  delete[] rowsDropped;
     2031  if (cholesky_->type() < 20) {
     2032    // not KKT
     2033    cholesky_->solve(errorRegion_);
     2034    //create information for solution
     2035    multiplyAdd(errorRegion_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
     2036    CoinZeroN(deltaX_, numberColumns_);
     2037    matrix_->transposeTimes(1.0, errorRegion_, deltaX_);
     2038  } else {
     2039    // KKT
     2040    // reverse sign on solution
     2041    multiplyAdd(NULL, numberRows_ + numberColumns_, 0.0, solution_, -1.0);
     2042    solveSystem(deltaX_, errorRegion_, solution_, NULL, NULL, NULL, false);
     2043  }
     2044  CoinWorkDouble initialValue = 1.0e2;
     2045  if (rhsNorm_ * 1.0e-2 > initialValue) {
     2046    initialValue = rhsNorm_ * 1.0e-2;
     2047  }
     2048  //initialValue = CoinMax(1.0,rhsNorm_);
     2049  CoinWorkDouble smallestBoundDifference = COIN_DBL_MAX;
     2050  CoinWorkDouble *fakeSolution = deltaX_;
     2051  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2052    if (!flagged(iColumn)) {
     2053      if (lower_[iColumn] - fakeSolution[iColumn] > initialValue) {
     2054        initialValue = lower_[iColumn] - fakeSolution[iColumn];
     2055      }
     2056      if (fakeSolution[iColumn] - upper_[iColumn] > initialValue) {
     2057        initialValue = fakeSolution[iColumn] - upper_[iColumn];
     2058      }
     2059      if (upper_[iColumn] - lower_[iColumn] < smallestBoundDifference) {
     2060        smallestBoundDifference = upper_[iColumn] - lower_[iColumn];
     2061      }
     2062    }
     2063  }
     2064  solutionNorm_ = 1.0e-12;
     2065  handler_->message(CLP_BARRIER_SAFE, messages_)
     2066    << static_cast< double >(initialValue) << static_cast< double >(objectiveNorm_)
     2067    << CoinMessageEol;
     2068  CoinWorkDouble extra = 1.0e-10;
     2069  CoinWorkDouble largeGap = 1.0e15;
     2070  //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_;
     2071  CoinWorkDouble safeObjectiveValue = objectiveNorm_ + 1.0;
     2072  CoinWorkDouble safeFree = 1.0e-1 * initialValue;
     2073  //printf("normal safe dual value of %g, primal value of %g\n",
     2074  // safeObjectiveValue,initialValue);
     2075  //safeObjectiveValue=CoinMax(2.0,1.0e-1*safeObjectiveValue);
     2076  //initialValue=CoinMax(100.0,1.0e-1*initialValue);;
     2077  //printf("temp safe dual value of %g, primal value of %g\n",
     2078  // safeObjectiveValue,initialValue);
     2079  CoinWorkDouble zwLarge = 1.0e2 * initialValue;
     2080  //zwLarge=1.0e40;
     2081  if (cholesky_->choleskyCondition() < 0.0 && cholesky_->type() < 20) {
     2082    // looks bad - play safe
     2083    initialValue *= 10.0;
     2084    safeObjectiveValue *= 10.0;
     2085    safeFree *= 10.0;
     2086  }
     2087  CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
     2088  // First do primal side
     2089  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2090    if (!flagged(iColumn)) {
     2091      CoinWorkDouble lowerValue = lower_[iColumn];
     2092      CoinWorkDouble upperValue = upper_[iColumn];
     2093      CoinWorkDouble newValue;
     2094      CoinWorkDouble setPrimal = initialValue;
     2095      if (quadraticObj) {
     2096        // perturb primal solution a bit
     2097        //fakeSolution[iColumn]  *= 0.002*CoinDrand48()+0.999;
     2098      }
     2099      if (lowerBound(iColumn)) {
     2100        if (upperBound(iColumn)) {
     2101          //upper and lower bounds
     2102          if (upperValue - lowerValue > 2.0 * setPrimal) {
     2103            CoinWorkDouble fakeValue = fakeSolution[iColumn];
     2104            if (fakeValue < lowerValue + setPrimal) {
     2105              fakeValue = lowerValue + setPrimal;
     2106            }
     2107            if (fakeValue > upperValue - setPrimal) {
     2108              fakeValue = upperValue - setPrimal;
     2109            }
     2110            newValue = fakeValue;
     2111          } else {
     2112            newValue = 0.5 * (upperValue + lowerValue);
     2113          }
     2114        } else {
     2115          //just lower bound
     2116          CoinWorkDouble fakeValue = fakeSolution[iColumn];
     2117          if (fakeValue < lowerValue + setPrimal) {
     2118            fakeValue = lowerValue + setPrimal;
     2119          }
     2120          newValue = fakeValue;
     2121        }
     2122      } else {
     2123        if (upperBound(iColumn)) {
     2124          //just upper bound
     2125          CoinWorkDouble fakeValue = fakeSolution[iColumn];
     2126          if (fakeValue > upperValue - setPrimal) {
     2127            fakeValue = upperValue - setPrimal;
     2128          }
     2129          newValue = fakeValue;
     2130        } else {
     2131          //free
     2132          newValue = fakeSolution[iColumn];
     2133          if (newValue >= 0.0) {
     2134            if (newValue < safeFree) {
     2135              newValue = safeFree;
     2136            }
     2137          } else {
     2138            if (newValue > -safeFree) {
     2139              newValue = -safeFree;
     2140            }
     2141          }
     2142        }
     2143      }
     2144      solution_[iColumn] = newValue;
     2145    } else {
     2146      // fixed
     2147      lowerSlack_[iColumn] = 0.0;
     2148      upperSlack_[iColumn] = 0.0;
     2149      solution_[iColumn] = lower_[iColumn];
     2150      zVec_[iColumn] = 0.0;
     2151      wVec_[iColumn] = 0.0;
     2152      diagonal_[iColumn] = 0.0;
     2153    }
     2154  }
     2155  solutionNorm_ = maximumAbsElement(solution_, numberTotal);
     2156  // Set bounds and do dj including quadratic
     2157  largeGap = CoinMax(1.0e7, 1.02 * solutionNorm_);
     2158  CoinPackedMatrix *quadratic = NULL;
     2159  const int *columnQuadratic = NULL;
     2160  const CoinBigIndex *columnQuadraticStart = NULL;
     2161  const int *columnQuadraticLength = NULL;
     2162  const double *quadraticElement = NULL;
     2163  if (quadraticObj) {
     2164    quadratic = quadraticObj->quadraticObjective();
     2165    columnQuadratic = quadratic->getIndices();
     2166    columnQuadraticStart = quadratic->getVectorStarts();
     2167    columnQuadraticLength = quadratic->getVectorLengths();
     2168    quadraticElement = quadratic->getElements();
     2169  }
     2170  CoinWorkDouble quadraticNorm = 0.0;
     2171  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2172    if (!flagged(iColumn)) {
     2173      CoinWorkDouble primalValue = solution_[iColumn];
     2174      CoinWorkDouble lowerValue = lower_[iColumn];
     2175      CoinWorkDouble upperValue = upper_[iColumn];
     2176      // Do dj
     2177      CoinWorkDouble reducedCost = cost_[iColumn];
     2178      if (lowerBound(iColumn)) {
     2179        reducedCost += linearPerturbation_;
     2180      }
     2181      if (upperBound(iColumn)) {
     2182        reducedCost -= linearPerturbation_;
     2183      }
     2184      if (quadraticObj && iColumn < numberColumns_) {
     2185        for (CoinBigIndex j = columnQuadraticStart[iColumn];
     2186             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     2187          int jColumn = columnQuadratic[j];
     2188          CoinWorkDouble valueJ = solution_[jColumn];
     2189          CoinWorkDouble elementValue = quadraticElement[j];
     2190          reducedCost += valueJ * elementValue;
     2191        }
     2192        quadraticNorm = CoinMax(quadraticNorm, CoinAbs(reducedCost));
     2193      }
     2194      dj_[iColumn] = reducedCost;
     2195      if (primalValue > lowerValue + largeGap && primalValue < upperValue - largeGap) {
     2196        clearFixedOrFree(iColumn);
     2197        setLowerBound(iColumn);
     2198        setUpperBound(iColumn);
     2199        lowerValue = CoinMax(lowerValue, primalValue - largeGap);
     2200        upperValue = CoinMin(upperValue, primalValue + largeGap);
     2201        lower_[iColumn] = lowerValue;
     2202        upper_[iColumn] = upperValue;
     2203      }
     2204    }
     2205  }
     2206  safeObjectiveValue = CoinMax(safeObjectiveValue, quadraticNorm);
     2207  for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2208    if (!flagged(iColumn)) {
     2209      CoinWorkDouble primalValue = solution_[iColumn];
     2210      CoinWorkDouble lowerValue = lower_[iColumn];
     2211      CoinWorkDouble upperValue = upper_[iColumn];
     2212      CoinWorkDouble reducedCost = dj_[iColumn];
     2213      CoinWorkDouble low = 0.0;
     2214      CoinWorkDouble high = 0.0;
     2215      if (lowerBound(iColumn)) {
     2216        if (upperBound(iColumn)) {
     2217          //upper and lower bounds
     2218          if (upperValue - lowerValue > 2.0 * initialValue) {
     2219            low = primalValue - lowerValue;
     2220            high = upperValue - primalValue;
     2221          } else {
     2222            low = initialValue;
     2223            high = initialValue;
     2224          }
     2225          CoinWorkDouble s = low + extra;
     2226          CoinWorkDouble ratioZ;
     2227          if (s < zwLarge) {
     2228            ratioZ = 1.0;
     2229          } else {
     2230            ratioZ = CoinSqrt(zwLarge / s);
     2231          }
     2232          CoinWorkDouble t = high + extra;
     2233          CoinWorkDouble ratioT;
     2234          if (t < zwLarge) {
     2235            ratioT = 1.0;
     2236          } else {
     2237            ratioT = CoinSqrt(zwLarge / t);
     2238          }
     2239          //modify s and t
     2240          if (s > largeGap) {
     2241            s = largeGap;
     2242          }
     2243          if (t > largeGap) {
     2244            t = largeGap;
     2245          }
     2246          //modify if long long way away from bound
     2247          if (reducedCost >= 0.0) {
     2248            zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
     2249            zVec_[iColumn] = CoinMax(reducedCost, safeObjectiveValue * ratioZ);
     2250            wVec_[iColumn] = safeObjectiveValue * ratioT;
     2251          } else {
     2252            zVec_[iColumn] = safeObjectiveValue * ratioZ;
     2253            wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
     2254            wVec_[iColumn] = CoinMax(-reducedCost, safeObjectiveValue * ratioT);
     2255          }
     2256          CoinWorkDouble gammaTerm = gamma2;
     2257          if (primalR_)
     2258            gammaTerm += primalR_[iColumn];
     2259          diagonal_[iColumn] = (t * s) / (s * wVec_[iColumn] + t * zVec_[iColumn] + gammaTerm * t * s);
     2260        } else {
     2261          //just lower bound
     2262          low = primalValue - lowerValue;
     2263          high = 0.0;
     2264          CoinWorkDouble s = low + extra;
     2265          CoinWorkDouble ratioZ;
     2266          if (s < zwLarge) {
     2267            ratioZ = 1.0;
     2268          } else {
     2269            ratioZ = CoinSqrt(zwLarge / s);
     2270          }
     2271          //modify s
     2272          if (s > largeGap) {
     2273            s = largeGap;
     2274          }
     2275          if (reducedCost >= 0.0) {
     2276            zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
     2277            zVec_[iColumn] = CoinMax(reducedCost, safeObjectiveValue * ratioZ);
     2278            wVec_[iColumn] = 0.0;
     2279          } else {
     2280            zVec_[iColumn] = safeObjectiveValue * ratioZ;
     2281            wVec_[iColumn] = 0.0;
     2282          }
     2283          CoinWorkDouble gammaTerm = gamma2;
     2284          if (primalR_)
     2285            gammaTerm += primalR_[iColumn];
     2286          diagonal_[iColumn] = s / (zVec_[iColumn] + s * gammaTerm);
     2287        }
     2288      } else {
     2289        if (upperBound(iColumn)) {
     2290          //just upper bound
     2291          low = 0.0;
     2292          high = upperValue - primalValue;
     2293          CoinWorkDouble t = high + extra;
     2294          CoinWorkDouble ratioT;
     2295          if (t < zwLarge) {
     2296            ratioT = 1.0;
     2297          } else {
     2298            ratioT = CoinSqrt(zwLarge / t);
     2299          }
     2300          //modify t
     2301          if (t > largeGap) {
     2302            t = largeGap;
     2303          }
     2304          if (reducedCost >= 0.0) {
     2305            zVec_[iColumn] = 0.0;
     2306            wVec_[iColumn] = safeObjectiveValue * ratioT;
     2307          } else {
     2308            zVec_[iColumn] = 0.0;
     2309            wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
     2310            wVec_[iColumn] = CoinMax(-reducedCost, safeObjectiveValue * ratioT);
     2311          }
     2312          CoinWorkDouble gammaTerm = gamma2;
     2313          if (primalR_)
     2314            gammaTerm += primalR_[iColumn];
     2315          diagonal_[iColumn] = t / (wVec_[iColumn] + t * gammaTerm);
     2316        }
     2317      }
     2318      lowerSlack_[iColumn] = low;
     2319      upperSlack_[iColumn] = high;
     2320    }
     2321  }
    23482322#if 0
    23492323     if (solution_[0] > 0.0) {
     
    23622336     exit(66);
    23632337#endif
    2364      return 0;
     2338  return 0;
    23652339}
    23662340// complementarityGap.  Computes gap
    23672341//phase 0=as is , 1 = after predictor , 2 after corrector
    2368 CoinWorkDouble ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
    2369           int & numberComplementarityItems,
    2370           const int phase)
     2342CoinWorkDouble ClpPredictorCorrector::complementarityGap(int &numberComplementarityPairs,
     2343  int &numberComplementarityItems,
     2344  const int phase)
    23712345{
    2372      CoinWorkDouble gap = 0.0;
    2373      //seems to be same coding for phase = 1 or 2
    2374      numberComplementarityPairs = 0;
    2375      numberComplementarityItems = 0;
    2376      int numberTotal = numberRows_ + numberColumns_;
    2377      CoinWorkDouble toleranceGap = 0.0;
    2378      CoinWorkDouble largestGap = 0.0;
    2379      CoinWorkDouble smallestGap = COIN_DBL_MAX;
    2380      //seems to be same coding for phase = 1 or 2
    2381      int numberNegativeGaps = 0;
    2382      CoinWorkDouble sumNegativeGap = 0.0;
    2383      CoinWorkDouble largeGap = 1.0e2 * solutionNorm_;
    2384      if (largeGap < 1.0e10) {
    2385           largeGap = 1.0e10;
    2386      }
    2387      largeGap = 1.0e30;
    2388      CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
    2389      CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance];
    2390      dualTolerance = dualTolerance / scaleFactor_;
    2391      for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
    2392           if (!fixedOrFree(iColumn)) {
    2393                numberComplementarityPairs++;
    2394                //can collapse as if no lower bound both zVec and deltaZ 0.0
    2395                if (lowerBound(iColumn)) {
    2396                     numberComplementarityItems++;
    2397                     CoinWorkDouble dualValue;
    2398                     CoinWorkDouble primalValue;
    2399                     if (!phase) {
    2400                          dualValue = zVec_[iColumn];
    2401                          primalValue = lowerSlack_[iColumn];
    2402                     } else {
    2403                          CoinWorkDouble change;
    2404                          change = solution_[iColumn] + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
    2405                          dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
    2406                          primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
    2407                     }
    2408                     //reduce primalValue
    2409                     if (primalValue > largeGap) {
    2410                          primalValue = largeGap;
    2411                     }
    2412                     CoinWorkDouble gapProduct = dualValue * primalValue;
    2413                     if (gapProduct < 0.0) {
    2414                          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
    2415                          //primalValue<<endl;
    2416                          numberNegativeGaps++;
    2417                          sumNegativeGap -= gapProduct;
    2418                          gapProduct = 0.0;
    2419                     }
    2420                     gap += gapProduct;
    2421                     //printf("l %d prim %g dual %g totalGap %g\n",
    2422                     //   iColumn,primalValue,dualValue,gap);
    2423                     if (gapProduct > largestGap) {
    2424                          largestGap = gapProduct;
    2425                     }
    2426                     smallestGap = CoinMin(smallestGap, gapProduct);
    2427                     if (dualValue > dualTolerance && primalValue > primalTolerance) {
    2428                          toleranceGap += dualValue * primalValue;
    2429                     }
    2430                }
    2431                if (upperBound(iColumn)) {
    2432                     numberComplementarityItems++;
    2433                     CoinWorkDouble dualValue;
    2434                     CoinWorkDouble primalValue;
    2435                     if (!phase) {
    2436                          dualValue = wVec_[iColumn];
    2437                          primalValue = upperSlack_[iColumn];
    2438                     } else {
    2439                          CoinWorkDouble change;
    2440                          change = upper_[iColumn] - solution_[iColumn] - deltaX_[iColumn] - upperSlack_[iColumn];
    2441                          dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
    2442                          primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
    2443                     }
    2444                     //reduce primalValue
    2445                     if (primalValue > largeGap) {
    2446                          primalValue = largeGap;
    2447                     }
    2448                     CoinWorkDouble gapProduct = dualValue * primalValue;
    2449                     if (gapProduct < 0.0) {
    2450                          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
    2451                          //primalValue<<endl;
    2452                          numberNegativeGaps++;
    2453                          sumNegativeGap -= gapProduct;
    2454                          gapProduct = 0.0;
    2455                     }
    2456                     gap += gapProduct;
    2457                     //printf("u %d prim %g dual %g totalGap %g\n",
    2458                     //   iColumn,primalValue,dualValue,gap);
    2459                     if (gapProduct > largestGap) {
    2460                          largestGap = gapProduct;
    2461                     }
    2462                     if (dualValue > dualTolerance && primalValue > primalTolerance) {
    2463                          toleranceGap += dualValue * primalValue;
    2464                     }
    2465                }
    2466           }
    2467      }
    2468      //if (numberIterations_>4)
    2469      //exit(9);
    2470      if (!phase && numberNegativeGaps) {
    2471           handler_->message(CLP_BARRIER_NEGATIVE_GAPS, messages_)
    2472                     << numberNegativeGaps << static_cast<double>(sumNegativeGap)
    2473                     << CoinMessageEol;
    2474      }
     2346  CoinWorkDouble gap = 0.0;
     2347  //seems to be same coding for phase = 1 or 2
     2348  numberComplementarityPairs = 0;
     2349  numberComplementarityItems = 0;
     2350  int numberTotal = numberRows_ + numberColumns_;
     2351  CoinWorkDouble toleranceGap = 0.0;
     2352  CoinWorkDouble largestGap = 0.0;
     2353  CoinWorkDouble smallestGap = COIN_DBL_MAX;
     2354  //seems to be same coding for phase = 1 or 2
     2355  int numberNegativeGaps = 0;
     2356  CoinWorkDouble sumNegativeGap = 0.0;
     2357  CoinWorkDouble largeGap = 1.0e2 * solutionNorm_;
     2358  if (largeGap < 1.0e10) {
     2359    largeGap = 1.0e10;
     2360  }
     2361  largeGap = 1.0e30;
     2362  CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
     2363  CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance];
     2364  dualTolerance = dualTolerance / scaleFactor_;
     2365  for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
     2366    if (!fixedOrFree(iColumn)) {
     2367      numberComplementarityPairs++;
     2368      //can collapse as if no lower bound both zVec and deltaZ 0.0
     2369      if (lowerBound(iColumn)) {
     2370        numberComplementarityItems++;
     2371        CoinWorkDouble dualValue;
     2372        CoinWorkDouble primalValue;
     2373        if (!phase) {
     2374          dualValue = zVec_[iColumn];
     2375          primalValue = lowerSlack_[iColumn];
     2376        } else {
     2377          CoinWorkDouble change;
     2378          change = solution_[iColumn] + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
     2379          dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
     2380          primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
     2381        }
     2382        //reduce primalValue
     2383        if (primalValue > largeGap) {
     2384          primalValue = largeGap;
     2385        }
     2386        CoinWorkDouble gapProduct = dualValue * primalValue;
     2387        if (gapProduct < 0.0) {
     2388          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
     2389          //primalValue<<endl;
     2390          numberNegativeGaps++;
     2391          sumNegativeGap -= gapProduct;
     2392          gapProduct = 0.0;
     2393        }
     2394        gap += gapProduct;
     2395        //printf("l %d prim %g dual %g totalGap %g\n",
     2396        //   iColumn,primalValue,dualValue,gap);
     2397        if (gapProduct > largestGap) {
     2398          largestGap = gapProduct;
     2399        }
     2400        smallestGap = CoinMin(smallestGap, gapProduct);
     2401        if (dualValue > dualTolerance && primalValue > primalTolerance) {
     2402          toleranceGap += dualValue * primalValue;
     2403        }
     2404      }
     2405      if (upperBound(iColumn)) {
     2406        numberComplementarityItems++;
     2407        CoinWorkDouble dualValue;
     2408        CoinWorkDouble primalValue;
     2409        if (!phase) {
     2410          dualValue = wVec_[iColumn];
     2411          primalValue = upperSlack_[iColumn];
     2412        } else {
     2413          CoinWorkDouble change;
     2414          change = upper_[iColumn] - solution_[iColumn] - deltaX_[iColumn] - upperSlack_[iColumn];
     2415          dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
     2416          primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
     2417        }
     2418        //reduce primalValue
     2419        if (primalValue > largeGap) {
     2420          primalValue = largeGap;
     2421        }
     2422        CoinWorkDouble gapProduct = dualValue * primalValue;
     2423        if (gapProduct < 0.0) {
     2424          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
     2425          //primalValue<<endl;
     2426          numberNegativeGaps++;
     2427          sumNegativeGap -= gapProduct;
     2428          gapProduct = 0.0;
     2429        }
     2430        gap += gapProduct;
     2431        //printf("u %d prim %g dual %g totalGap %g\n",
     2432        //   iColumn,primalValue,dualValue,gap);
     2433        if (gapProduct > largestGap) {
     2434          largestGap = gapProduct;
     2435        }
     2436        if (dualValue > dualTolerance && primalValue > primalTolerance) {
     2437          toleranceGap += dualValue * primalValue;
     2438        }
     2439      }
     2440    }
     2441  }
     2442  //if (numberIterations_>4)
     2443  //exit(9);
     2444  if (!phase && numberNegativeGaps) {
     2445    handler_->message(CLP_BARRIER_NEGATIVE_GAPS, messages_)
     2446      << numberNegativeGaps << static_cast< double >(sumNegativeGap)
     2447      << CoinMessageEol;
     2448  }
    24752449
    2476      //in case all free!
    2477      if (!numberComplementarityPairs) {
    2478           numberComplementarityPairs = 1;
    2479      }
     2450  //in case all free!
     2451  if (!numberComplementarityPairs) {
     2452    numberComplementarityPairs = 1;
     2453  }
    24802454#ifdef SOME_DEBUG
    2481      printf("with d,p steps %g,%g gap %g - smallest %g, largest %g, pairs %d\n",
    2482             actualDualStep_, actualPrimalStep_,
    2483             gap, smallestGap, largestGap, numberComplementarityPairs);
    2484 #endif
    2485      return gap;
     2455  printf("with d,p steps %g,%g gap %g - smallest %g, largest %g, pairs %d\n",
     2456    actualDualStep_, actualPrimalStep_,
     2457    gap, smallestGap, largestGap, numberComplementarityPairs);
     2458#endif
     2459  return gap;
    24862460}
    24872461// setupForSolve.
     
    24892463void ClpPredictorCorrector::setupForSolve(const int phase)
    24902464{
    2491      CoinWorkDouble extra = eExtra;
    2492      int numberTotal = numberRows_ + numberColumns_;
    2493      int iColumn;
     2465  CoinWorkDouble extra = eExtra;
     2466  int numberTotal = numberRows_ + numberColumns_;
     2467  int iColumn;
    24942468#ifdef SOME_DEBUG
    2495      printf("phase %d in setupForSolve, mu %.18g\n", phase, mu_);
    2496 #endif
    2497      CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
    2498      CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    2499      switch (phase) {
    2500      case 0:
    2501           CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
    2502           if (delta_ || dualR_) {
    2503                // add in regularization
    2504                CoinWorkDouble delta2 = delta_ * delta_;
    2505                for (int iRow = 0; iRow < numberRows_; iRow++) {
    2506                     rhsB_[iRow] -= delta2 * dualArray[iRow];
    2507                     if (dualR_)
    2508                          rhsB_[iRow] -= dualR_[iRow] * dualArray[iRow];
    2509                }
    2510           }
    2511           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    2512                rhsC_[iColumn] = 0.0;
    2513                rhsU_[iColumn] = 0.0;
    2514                rhsL_[iColumn] = 0.0;
    2515                rhsZ_[iColumn] = 0.0;
    2516                rhsW_[iColumn] = 0.0;
    2517                if (!flagged(iColumn)) {
    2518                     rhsC_[iColumn] = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn];
    2519                     rhsC_[iColumn] += gamma2 * solution_[iColumn];
    2520                     if (primalR_)
    2521                          rhsC_[iColumn] += primalR_[iColumn] * solution_[iColumn];
    2522                     if (lowerBound(iColumn)) {
    2523                          rhsZ_[iColumn] = -zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
    2524                          rhsL_[iColumn] = CoinMax(0.0, (lower_[iColumn] + lowerSlack_[iColumn]) - solution_[iColumn]);
    2525                     }
    2526                     if (upperBound(iColumn)) {
    2527                          rhsW_[iColumn] = -wVec_[iColumn] * (upperSlack_[iColumn] + extra);
    2528                          rhsU_[iColumn] = CoinMin(0.0, (upper_[iColumn] - upperSlack_[iColumn]) - solution_[iColumn]);
    2529                     }
    2530                }
    2531           }
     2469  printf("phase %d in setupForSolve, mu %.18g\n", phase, mu_);
     2470#endif
     2471  CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
     2472  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
     2473  switch (phase) {
     2474  case 0:
     2475    CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
     2476    if (delta_ || dualR_) {
     2477      // add in regularization
     2478      CoinWorkDouble delta2 = delta_ * delta_;
     2479      for (int iRow = 0; iRow < numberRows_; iRow++) {
     2480        rhsB_[iRow] -= delta2 * dualArray[iRow];
     2481        if (dualR_)
     2482          rhsB_[iRow] -= dualR_[iRow] * dualArray[iRow];
     2483      }
     2484    }
     2485    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2486      rhsC_[iColumn] = 0.0;
     2487      rhsU_[iColumn] = 0.0;
     2488      rhsL_[iColumn] = 0.0;
     2489      rhsZ_[iColumn] = 0.0;
     2490      rhsW_[iColumn] = 0.0;
     2491      if (!flagged(iColumn)) {
     2492        rhsC_[iColumn] = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn];
     2493        rhsC_[iColumn] += gamma2 * solution_[iColumn];
     2494        if (primalR_)
     2495          rhsC_[iColumn] += primalR_[iColumn] * solution_[iColumn];
     2496        if (lowerBound(iColumn)) {
     2497          rhsZ_[iColumn] = -zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
     2498          rhsL_[iColumn] = CoinMax(0.0, (lower_[iColumn] + lowerSlack_[iColumn]) - solution_[iColumn]);
     2499        }
     2500        if (upperBound(iColumn)) {
     2501          rhsW_[iColumn] = -wVec_[iColumn] * (upperSlack_[iColumn] + extra);
     2502          rhsU_[iColumn] = CoinMin(0.0, (upper_[iColumn] - upperSlack_[iColumn]) - solution_[iColumn]);
     2503        }
     2504      }
     2505    }
    25322506#if 0
    25332507          for (int i = 0; i < 3; i++) {
     
    25632537          }
    25642538#endif
    2565           break;
    2566      case 1:
    2567           // could be stored in delta2?
    2568           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    2569                rhsZ_[iColumn] = 0.0;
    2570                rhsW_[iColumn] = 0.0;
    2571                if (!flagged(iColumn)) {
    2572                     if (lowerBound(iColumn)) {
    2573                          rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra)
    2574                                           - deltaZ_[iColumn] * deltaX_[iColumn];
    2575                          // To bring in line with OSL
    2576                          rhsZ_[iColumn] += deltaZ_[iColumn] * rhsL_[iColumn];
    2577                     }
    2578                     if (upperBound(iColumn)) {
    2579                          rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra)
    2580                                           + deltaW_[iColumn] * deltaX_[iColumn];
    2581                          // To bring in line with OSL
    2582                          rhsW_[iColumn] -= deltaW_[iColumn] * rhsU_[iColumn];
    2583                     }
    2584                }
    2585           }
     2539    break;
     2540  case 1:
     2541    // could be stored in delta2?
     2542    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2543      rhsZ_[iColumn] = 0.0;
     2544      rhsW_[iColumn] = 0.0;
     2545      if (!flagged(iColumn)) {
     2546        if (lowerBound(iColumn)) {
     2547          rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra)
     2548            - deltaZ_[iColumn] * deltaX_[iColumn];
     2549          // To bring in line with OSL
     2550          rhsZ_[iColumn] += deltaZ_[iColumn] * rhsL_[iColumn];
     2551        }
     2552        if (upperBound(iColumn)) {
     2553          rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra)
     2554            + deltaW_[iColumn] * deltaX_[iColumn];
     2555          // To bring in line with OSL
     2556          rhsW_[iColumn] -= deltaW_[iColumn] * rhsU_[iColumn];
     2557        }
     2558      }
     2559    }
    25862560#if 0
    25872561          for (int i = 0; i < numberTotal; i++) {
     
    26182592          exit(66);
    26192593#endif
    2620           break;
    2621      case 2:
    2622           CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
    2623           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    2624                rhsZ_[iColumn] = 0.0;
    2625                rhsW_[iColumn] = 0.0;
    2626                if (!flagged(iColumn)) {
    2627                     if (lowerBound(iColumn)) {
    2628                          rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
    2629                     }
    2630                     if (upperBound(iColumn)) {
    2631                          rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra);
    2632                     }
    2633                }
    2634           }
    2635           break;
    2636      case 3: {
    2637           CoinWorkDouble minBeta = 0.1 * mu_;
    2638           CoinWorkDouble maxBeta = 10.0 * mu_;
    2639           CoinWorkDouble dualStep = CoinMin(1.0, actualDualStep_ + 0.1);
    2640           CoinWorkDouble primalStep = CoinMin(1.0, actualPrimalStep_ + 0.1);
     2594    break;
     2595  case 2:
     2596    CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
     2597    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2598      rhsZ_[iColumn] = 0.0;
     2599      rhsW_[iColumn] = 0.0;
     2600      if (!flagged(iColumn)) {
     2601        if (lowerBound(iColumn)) {
     2602          rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
     2603        }
     2604        if (upperBound(iColumn)) {
     2605          rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra);
     2606        }
     2607      }
     2608    }
     2609    break;
     2610  case 3: {
     2611    CoinWorkDouble minBeta = 0.1 * mu_;
     2612    CoinWorkDouble maxBeta = 10.0 * mu_;
     2613    CoinWorkDouble dualStep = CoinMin(1.0, actualDualStep_ + 0.1);
     2614    CoinWorkDouble primalStep = CoinMin(1.0, actualPrimalStep_ + 0.1);
    26412615#ifdef SOME_DEBUG
    2642           printf("good complementarity range %g to %g\n", minBeta, maxBeta);
    2643 #endif
    2644           //minBeta=0.0;
    2645           //maxBeta=COIN_DBL_MAX;
    2646           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    2647                if (!flagged(iColumn)) {
    2648                     if (lowerBound(iColumn)) {
    2649                          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
    2650                          CoinWorkDouble dualValue = zVec_[iColumn] + dualStep * deltaZ_[iColumn];
    2651                          CoinWorkDouble primalValue = lowerSlack_[iColumn] + primalStep * change;
    2652                          CoinWorkDouble gapProduct = dualValue * primalValue;
    2653                          if (gapProduct > 0.0 && dualValue < 0.0)
    2654                               gapProduct = - gapProduct;
     2616    printf("good complementarity range %g to %g\n", minBeta, maxBeta);
     2617#endif
     2618    //minBeta=0.0;
     2619    //maxBeta=COIN_DBL_MAX;
     2620    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2621      if (!flagged(iColumn)) {
     2622        if (lowerBound(iColumn)) {
     2623          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     2624          CoinWorkDouble dualValue = zVec_[iColumn] + dualStep * deltaZ_[iColumn];
     2625          CoinWorkDouble primalValue = lowerSlack_[iColumn] + primalStep * change;
     2626          CoinWorkDouble gapProduct = dualValue * primalValue;
     2627          if (gapProduct > 0.0 && dualValue < 0.0)
     2628            gapProduct = -gapProduct;
    26552629#ifdef FULL_DEBUG
    2656                          delta2Z_[iColumn] = gapProduct;
    2657                          if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
    2658                               printf("lower %d primal %g, dual %g, gap %g\n",
    2659                                      iColumn, primalValue, dualValue, gapProduct);
    2660 #endif
    2661                          CoinWorkDouble value = 0.0;
    2662                          if (gapProduct < minBeta) {
    2663                               value = 2.0 * (minBeta - gapProduct);
    2664                               value = (mu_ - gapProduct);
    2665                               value = (minBeta - gapProduct);
    2666                               assert (value > 0.0);
    2667                          } else if (gapProduct > maxBeta) {
    2668                               value = CoinMax(maxBeta - gapProduct, -maxBeta);
    2669                               assert (value < 0.0);
    2670                          }
    2671                          rhsZ_[iColumn] += value;
    2672                     }
    2673                     if (upperBound(iColumn)) {
    2674                          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
    2675                          CoinWorkDouble dualValue = wVec_[iColumn] + dualStep * deltaW_[iColumn];
    2676                          CoinWorkDouble primalValue = upperSlack_[iColumn] + primalStep * change;
    2677                          CoinWorkDouble gapProduct = dualValue * primalValue;
    2678                          if (gapProduct > 0.0 && dualValue < 0.0)
    2679                               gapProduct = - gapProduct;
     2630          delta2Z_[iColumn] = gapProduct;
     2631          if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
     2632            printf("lower %d primal %g, dual %g, gap %g\n",
     2633              iColumn, primalValue, dualValue, gapProduct);
     2634#endif
     2635          CoinWorkDouble value = 0.0;
     2636          if (gapProduct < minBeta) {
     2637            value = 2.0 * (minBeta - gapProduct);
     2638            value = (mu_ - gapProduct);
     2639            value = (minBeta - gapProduct);
     2640            assert(value > 0.0);
     2641          } else if (gapProduct > maxBeta) {
     2642            value = CoinMax(maxBeta - gapProduct, -maxBeta);
     2643            assert(value < 0.0);
     2644          }
     2645          rhsZ_[iColumn] += value;
     2646        }
     2647        if (upperBound(iColumn)) {
     2648          CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
     2649          CoinWorkDouble dualValue = wVec_[iColumn] + dualStep * deltaW_[iColumn];
     2650          CoinWorkDouble primalValue = upperSlack_[iColumn] + primalStep * change;
     2651          CoinWorkDouble gapProduct = dualValue * primalValue;
     2652          if (gapProduct > 0.0 && dualValue < 0.0)
     2653            gapProduct = -gapProduct;
    26802654#ifdef FULL_DEBUG
    2681                          delta2W_[iColumn] = gapProduct;
    2682                          if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
    2683                               printf("upper %d primal %g, dual %g, gap %g\n",
    2684                                      iColumn, primalValue, dualValue, gapProduct);
    2685 #endif
    2686                          CoinWorkDouble value = 0.0;
    2687                          if (gapProduct < minBeta) {
    2688                               value = (minBeta - gapProduct);
    2689                               assert (value > 0.0);
    2690                          } else if (gapProduct > maxBeta) {
    2691                               value = CoinMax(maxBeta - gapProduct, -maxBeta);
    2692                               assert (value < 0.0);
    2693                          }
    2694                          rhsW_[iColumn] += value;
    2695                     }
    2696                }
    2697           }
    2698      }
    2699      break;
    2700      } /* endswitch */
    2701      if (cholesky_->type() < 20) {
    2702           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    2703                CoinWorkDouble value = rhsC_[iColumn];
    2704                CoinWorkDouble zValue = rhsZ_[iColumn];
    2705                CoinWorkDouble wValue = rhsW_[iColumn];
     2655          delta2W_[iColumn] = gapProduct;
     2656          if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
     2657            printf("upper %d primal %g, dual %g, gap %g\n",
     2658              iColumn, primalValue, dualValue, gapProduct);
     2659#endif
     2660          CoinWorkDouble value = 0.0;
     2661          if (gapProduct < minBeta) {
     2662            value = (minBeta - gapProduct);
     2663            assert(value > 0.0);
     2664          } else if (gapProduct > maxBeta) {
     2665            value = CoinMax(maxBeta - gapProduct, -maxBeta);
     2666            assert(value < 0.0);
     2667          }
     2668          rhsW_[iColumn] += value;
     2669        }
     2670      }
     2671    }
     2672  } break;
     2673  } /* endswitch */
     2674  if (cholesky_->type() < 20) {
     2675    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2676      CoinWorkDouble value = rhsC_[iColumn];
     2677      CoinWorkDouble zValue = rhsZ_[iColumn];
     2678      CoinWorkDouble wValue = rhsW_[iColumn];
    27062679#if 0
    27072680#if 1
     
    27292702               }
    27302703#else
    2731                if (lowerBound(iColumn)) {
    2732                     CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
    2733                     value -= gHat / (lowerSlack_[iColumn] + extra);
    2734                }
    2735                if (upperBound(iColumn)) {
    2736                     CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
    2737                     value += hHat / (upperSlack_[iColumn] + extra);
    2738                }
    2739 #endif
    2740                workArray_[iColumn] = diagonal_[iColumn] * value;
    2741           }
     2704      if (lowerBound(iColumn)) {
     2705        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
     2706        value -= gHat / (lowerSlack_[iColumn] + extra);
     2707      }
     2708      if (upperBound(iColumn)) {
     2709        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
     2710        value += hHat / (upperSlack_[iColumn] + extra);
     2711      }
     2712#endif
     2713      workArray_[iColumn] = diagonal_[iColumn] * value;
     2714    }
    27422715#if 0
    27432716          if (solution_[0] > 0.0) {
     
    27502723          exit(66);
    27512724#endif
    2752      } else {
    2753           // KKT
    2754           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    2755                CoinWorkDouble value = rhsC_[iColumn];
    2756                CoinWorkDouble zValue = rhsZ_[iColumn];
    2757                CoinWorkDouble wValue = rhsW_[iColumn];
    2758                if (lowerBound(iColumn)) {
    2759                     CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
    2760                     value -= gHat / (lowerSlack_[iColumn] + extra);
    2761                }
    2762                if (upperBound(iColumn)) {
    2763                     CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
    2764                     value += hHat / (upperSlack_[iColumn] + extra);
    2765                }
    2766                workArray_[iColumn] = value;
    2767           }
    2768      }
     2725  } else {
     2726    // KKT
     2727    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     2728      CoinWorkDouble value = rhsC_[iColumn];
     2729      CoinWorkDouble zValue = rhsZ_[iColumn];
     2730      CoinWorkDouble wValue = rhsW_[iColumn];
     2731      if (lowerBound(iColumn)) {
     2732        CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
     2733        value -= gHat / (lowerSlack_[iColumn] + extra);
     2734      }
     2735      if (upperBound(iColumn)) {
     2736        CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
     2737        value += hHat / (upperSlack_[iColumn] + extra);
     2738      }
     2739      workArray_[iColumn] = value;
     2740    }
     2741  }
    27692742}
    27702743//method: sees if looks plausible change in complementarity
    27712744bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
    2772           CoinWorkDouble & bestNextGap,
    2773           bool allowIncreasingGap)
     2745  CoinWorkDouble &bestNextGap,
     2746  bool allowIncreasingGap)
    27742747{
    2775      const CoinWorkDouble beta3 = 0.99997;
    2776      bool goodMove = false;
    2777      int nextNumber;
    2778      int nextNumberItems;
    2779      int numberTotal = numberRows_ + numberColumns_;
    2780      CoinWorkDouble returnGap = bestNextGap;
    2781      CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
     2748  const CoinWorkDouble beta3 = 0.99997;
     2749  bool goodMove = false;
     2750  int nextNumber;
     2751  int nextNumberItems;
     2752  int numberTotal = numberRows_ + numberColumns_;
     2753  CoinWorkDouble returnGap = bestNextGap;
     2754  CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
    27822755#ifndef NO_RTTI
    2783      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     2756  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
    27842757#else
    2785      ClpQuadraticObjective * quadraticObj = NULL;
    2786      if (objective_->type() == 2)
    2787           quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    2788 #endif
    2789      if (nextGap > bestNextGap && nextGap > 0.9 * complementarityGap_ && doCorrector
    2790                && !quadraticObj && !allowIncreasingGap) {
     2758  ClpQuadraticObjective *quadraticObj = NULL;
     2759  if (objective_->type() == 2)
     2760    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     2761#endif
     2762  if (nextGap > bestNextGap && nextGap > 0.9 * complementarityGap_ && doCorrector
     2763    && !quadraticObj && !allowIncreasingGap) {
    27912764#ifdef SOME_DEBUG
    2792           printf("checkGood phase 1 next gap %.18g, phase 0 %.18g, old gap %.18g\n",
    2793                  nextGap, bestNextGap, complementarityGap_);
    2794 #endif
    2795           return false;
    2796      } else {
    2797           returnGap = nextGap;
    2798      }
    2799      CoinWorkDouble step;
    2800      if (actualDualStep_ > actualPrimalStep_) {
    2801           step = actualDualStep_;
    2802      } else {
    2803           step = actualPrimalStep_;
    2804      }
    2805      CoinWorkDouble testValue = 1.0 - step * (1.0 - beta3);
    2806      //testValue=0.0;
    2807      testValue *= complementarityGap_;
    2808      if (nextGap < testValue) {
    2809           //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
    2810           goodMove = true;
    2811      } else if(doCorrector) {
    2812           //if (actualDualStep_<actualPrimalStep_) {
    2813           //step=actualDualStep_;
    2814           //} else {
    2815           //step=actualPrimalStep_;
    2816           //}
    2817           CoinWorkDouble gap = bestNextGap;
    2818           goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
    2819           if (goodMove)
    2820                returnGap = gap;
    2821      } else {
    2822           goodMove = true;
    2823      }
    2824      if (goodMove)
    2825           goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
    2826      // Say good if small
    2827      //if (quadraticObj) {
    2828      if (CoinMax(actualDualStep_, actualPrimalStep_) < 1.0e-6)
    2829           goodMove = true;
    2830      if (!goodMove) {
    2831           //try smaller of two
    2832           if (actualDualStep_ < actualPrimalStep_) {
    2833                step = actualDualStep_;
    2834           } else {
    2835                step = actualPrimalStep_;
    2836           }
    2837           if (step > 1.0) {
    2838                step = 1.0;
    2839           }
    2840           actualPrimalStep_ = step;
    2841           //if (quadraticObj)
    2842           //actualPrimalStep_ *=0.5;
    2843           actualDualStep_ = step;
    2844           goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
    2845           int pass = 0;
    2846           while (!goodMove) {
    2847                pass++;
    2848                CoinWorkDouble gap = bestNextGap;
    2849                goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
    2850                if (goodMove || pass > 3) {
    2851                     returnGap = gap;
    2852                     break;
    2853                }
    2854                if (step < 1.0e-4) {
    2855                     break;
    2856                }
    2857                step *= 0.5;
    2858                actualPrimalStep_ = step;
    2859                //if (quadraticObj)
    2860                //actualPrimalStep_ *=0.5;
    2861                actualDualStep_ = step;
    2862           } /* endwhile */
    2863           if (doCorrector) {
    2864                //say bad move if both small
    2865                if (numberIterations_ & 1) {
    2866                     if (actualPrimalStep_ < 1.0e-2 && actualDualStep_ < 1.0e-2) {
    2867                          goodMove = false;
    2868                     }
    2869                } else {
    2870                     if (actualPrimalStep_ < 1.0e-5 && actualDualStep_ < 1.0e-5) {
    2871                          goodMove = false;
    2872                     }
    2873                     if (actualPrimalStep_ * actualDualStep_ < 1.0e-20) {
    2874                          goodMove = false;
    2875                     }
    2876                }
    2877           }
    2878      }
    2879      if (goodMove) {
    2880           //compute delta in objectives
    2881           CoinWorkDouble deltaObjectivePrimal = 0.0;
    2882           CoinWorkDouble deltaObjectiveDual =
    2883                innerProduct(deltaY_, numberRows_,
    2884                             rhsFixRegion_);
    2885           CoinWorkDouble error = 0.0;
    2886           CoinWorkDouble * workArray = workArray_;
    2887           CoinZeroN(workArray, numberColumns_);
    2888           CoinMemcpyN(deltaY_, numberRows_, workArray + numberColumns_);
    2889           matrix_->transposeTimes(-1.0, deltaY_, workArray);
    2890           //CoinWorkDouble sumPerturbCost=0.0;
    2891           for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
    2892                if (!flagged(iColumn)) {
    2893                     if (lowerBound(iColumn)) {
    2894                          //sumPerturbCost+=deltaX_[iColumn];
    2895                          deltaObjectiveDual += deltaZ_[iColumn] * lower_[iColumn];
    2896                     }
    2897                     if (upperBound(iColumn)) {
    2898                          //sumPerturbCost-=deltaX_[iColumn];
    2899                          deltaObjectiveDual -= deltaW_[iColumn] * upper_[iColumn];
    2900                     }
    2901                     CoinWorkDouble change = CoinAbs(workArray_[iColumn] - deltaZ_[iColumn] + deltaW_[iColumn]);
    2902                     error = CoinMax(change, error);
    2903                }
    2904                deltaObjectivePrimal += cost_[iColumn] * deltaX_[iColumn];
    2905           }
    2906           //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
    2907           CoinWorkDouble testValue;
    2908           if (error > 0.0) {
    2909                testValue = 1.0e1 * CoinMax(maximumDualError_, 1.0e-12) / error;
    2910           } else {
    2911                testValue = 1.0e1;
    2912           }
    2913           // If quadratic then primal step may compensate
    2914           if (testValue < actualDualStep_ && !quadraticObj) {
    2915                handler_->message(CLP_BARRIER_REDUCING, messages_)
    2916                          << "dual" << static_cast<double>(actualDualStep_)
    2917                          << static_cast<double>(testValue)
    2918                          << CoinMessageEol;
    2919                actualDualStep_ = testValue;
    2920           }
    2921      }
    2922      if (maximumRHSError_ < 1.0e1 * solutionNorm_ * primalTolerance()
    2923                && maximumRHSChange_ > 1.0e-16 * solutionNorm_) {
    2924           //check change in AX not too much
    2925           //??? could be dropped row going infeasible
    2926           CoinWorkDouble ratio = 1.0e1 * CoinMax(maximumRHSError_, 1.0e-12) / maximumRHSChange_;
    2927           if (ratio < actualPrimalStep_) {
    2928                handler_->message(CLP_BARRIER_REDUCING, messages_)
    2929                          << "primal" << static_cast<double>(actualPrimalStep_)
    2930                          << static_cast<double>(ratio)
    2931                          << CoinMessageEol;
    2932                if (ratio > 1.0e-6) {
    2933                     actualPrimalStep_ = ratio;
    2934                } else {
    2935                     actualPrimalStep_ = ratio;
    2936                     //std::cout <<"sign we should be stopping"<<std::endl;
    2937                }
    2938           }
    2939      }
    2940      if (goodMove)
    2941           bestNextGap = returnGap;
    2942      return goodMove;
     2765    printf("checkGood phase 1 next gap %.18g, phase 0 %.18g, old gap %.18g\n",
     2766      nextGap, bestNextGap, complementarityGap_);
     2767#endif
     2768    return false;
     2769  } else {
     2770    returnGap = nextGap;
     2771  }
     2772  CoinWorkDouble step;
     2773  if (actualDualStep_ > actualPrimalStep_) {
     2774    step = actualDualStep_;
     2775  } else {
     2776    step = actualPrimalStep_;
     2777  }
     2778  CoinWorkDouble testValue = 1.0 - step * (1.0 - beta3);
     2779  //testValue=0.0;
     2780  testValue *= complementarityGap_;
     2781  if (nextGap < testValue) {
     2782    //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
     2783    goodMove = true;
     2784  } else if (doCorrector) {
     2785    //if (actualDualStep_<actualPrimalStep_) {
     2786    //step=actualDualStep_;
     2787    //} else {
     2788    //step=actualPrimalStep_;
     2789    //}
     2790    CoinWorkDouble gap = bestNextGap;
     2791    goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
     2792    if (goodMove)
     2793      returnGap = gap;
     2794  } else {
     2795    goodMove = true;
     2796  }
     2797  if (goodMove)
     2798    goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
     2799  // Say good if small
     2800  //if (quadraticObj) {
     2801  if (CoinMax(actualDualStep_, actualPrimalStep_) < 1.0e-6)
     2802    goodMove = true;
     2803  if (!goodMove) {
     2804    //try smaller of two
     2805    if (actualDualStep_ < actualPrimalStep_) {
     2806      step = actualDualStep_;
     2807    } else {
     2808      step = actualPrimalStep_;
     2809    }
     2810    if (step > 1.0) {
     2811      step = 1.0;
     2812    }
     2813    actualPrimalStep_ = step;
     2814    //if (quadraticObj)
     2815    //actualPrimalStep_ *=0.5;
     2816    actualDualStep_ = step;
     2817    goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
     2818    int pass = 0;
     2819    while (!goodMove) {
     2820      pass++;
     2821      CoinWorkDouble gap = bestNextGap;
     2822      goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
     2823      if (goodMove || pass > 3) {
     2824        returnGap = gap;
     2825        break;
     2826      }
     2827      if (step < 1.0e-4) {
     2828        break;
     2829      }
     2830      step *= 0.5;
     2831      actualPrimalStep_ = step;
     2832      //if (quadraticObj)
     2833      //actualPrimalStep_ *=0.5;
     2834      actualDualStep_ = step;
     2835    } /* endwhile */
     2836    if (doCorrector) {
     2837      //say bad move if both small
     2838      if (numberIterations_ & 1) {
     2839        if (actualPrimalStep_ < 1.0e-2 && actualDualStep_ < 1.0e-2) {
     2840          goodMove = false;
     2841        }
     2842      } else {
     2843        if (actualPrimalStep_ < 1.0e-5 && actualDualStep_ < 1.0e-5) {
     2844          goodMove = false;
     2845        }
     2846        if (actualPrimalStep_ * actualDualStep_ < 1.0e-20) {
     2847          goodMove = false;
     2848        }
     2849      }
     2850    }
     2851  }
     2852  if (goodMove) {
     2853    //compute delta in objectives
     2854    CoinWorkDouble deltaObjectivePrimal = 0.0;
     2855    CoinWorkDouble deltaObjectiveDual = innerProduct(deltaY_, numberRows_,
     2856      rhsFixRegion_);
     2857    CoinWorkDouble error = 0.0;
     2858    CoinWorkDouble *workArray = workArray_;
     2859    CoinZeroN(workArray, numberColumns_);
     2860    CoinMemcpyN(deltaY_, numberRows_, workArray + numberColumns_);
     2861    matrix_->transposeTimes(-1.0, deltaY_, workArray);
     2862    //CoinWorkDouble sumPerturbCost=0.0;
     2863    for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
     2864      if (!flagged(iColumn)) {
     2865        if (lowerBound(iColumn)) {
     2866          //sumPerturbCost+=deltaX_[iColumn];
     2867          deltaObjectiveDual += deltaZ_[iColumn] * lower_[iColumn];
     2868        }
     2869        if (upperBound(iColumn)) {
     2870          //sumPerturbCost-=deltaX_[iColumn];
     2871          deltaObjectiveDual -= deltaW_[iColumn] * upper_[iColumn];
     2872        }
     2873        CoinWorkDouble change = CoinAbs(workArray_[iColumn] - deltaZ_[iColumn] + deltaW_[iColumn]);
     2874        error = CoinMax(change, error);
     2875      }
     2876      deltaObjectivePrimal += cost_[iColumn] * deltaX_[iColumn];
     2877    }
     2878    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
     2879    CoinWorkDouble testValue;
     2880    if (error > 0.0) {
     2881      testValue = 1.0e1 * CoinMax(maximumDualError_, 1.0e-12) / error;
     2882    } else {
     2883      testValue = 1.0e1;
     2884    }
     2885    // If quadratic then primal step may compensate
     2886    if (testValue < actualDualStep_ && !quadraticObj) {
     2887      handler_->message(CLP_BARRIER_REDUCING, messages_)
     2888        << "dual" << static_cast< double >(actualDualStep_)
     2889        << static_cast< double >(testValue)
     2890        << CoinMessageEol;
     2891      actualDualStep_ = testValue;
     2892    }
     2893  }
     2894  if (maximumRHSError_ < 1.0e1 * solutionNorm_ * primalTolerance()
     2895    && maximumRHSChange_ > 1.0e-16 * solutionNorm_) {
     2896    //check change in AX not too much
     2897    //??? could be dropped row going infeasible
     2898    CoinWorkDouble ratio = 1.0e1 * CoinMax(maximumRHSError_, 1.0e-12) / maximumRHSChange_;
     2899    if (ratio < actualPrimalStep_) {
     2900      handler_->message(CLP_BARRIER_REDUCING, messages_)
     2901        << "primal" << static_cast< double >(actualPrimalStep_)
     2902        << static_cast< double >(ratio)
     2903        << CoinMessageEol;
     2904      if (ratio > 1.0e-6) {
     2905        actualPrimalStep_ = ratio;
     2906      } else {
     2907        actualPrimalStep_ = ratio;
     2908        //std::cout <<"sign we should be stopping"<<std::endl;
     2909      }
     2910    }
     2911  }
     2912  if (goodMove)
     2913    bestNextGap = returnGap;
     2914  return goodMove;
    29432915}
    29442916//:  checks for one step size
    29452917bool ClpPredictorCorrector::checkGoodMove2(CoinWorkDouble move,
    2946           CoinWorkDouble & bestNextGap,
    2947           bool allowIncreasingGap)
     2918  CoinWorkDouble &bestNextGap,
     2919  bool allowIncreasingGap)
    29482920{
    2949      CoinWorkDouble complementarityMultiplier = 1.0 / numberComplementarityPairs_;
    2950      const CoinWorkDouble gamma = 1.0e-8;
    2951      const CoinWorkDouble gammap = 1.0e-8;
    2952      CoinWorkDouble gammad = 1.0e-8;
    2953      int nextNumber;
    2954      int nextNumberItems;
    2955      CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
    2956      if (nextGap > bestNextGap && !allowIncreasingGap)
    2957           return false;
    2958      CoinWorkDouble lowerBoundGap = gamma * nextGap * complementarityMultiplier;
    2959      bool goodMove = true;
    2960      int iColumn;
    2961      for ( iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
    2962           if (!flagged(iColumn)) {
    2963                if (lowerBound(iColumn)) {
    2964                     CoinWorkDouble part1 = lowerSlack_[iColumn] + actualPrimalStep_ * deltaSL_[iColumn];
    2965                     CoinWorkDouble part2 = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
    2966                     if (part1 * part2 < lowerBoundGap) {
    2967                          goodMove = false;
    2968                          break;
    2969                     }
    2970                }
    2971                if (upperBound(iColumn)) {
    2972                     CoinWorkDouble part1 = upperSlack_[iColumn] + actualPrimalStep_ * deltaSU_[iColumn];
    2973                     CoinWorkDouble part2 = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
    2974                     if (part1 * part2 < lowerBoundGap) {
    2975                          goodMove = false;
    2976                          break;
    2977                     }
    2978                }
    2979           }
    2980      }
    2981      CoinWorkDouble * nextDj = NULL;
    2982      CoinWorkDouble maximumDualError = maximumDualError_;
     2921  CoinWorkDouble complementarityMultiplier = 1.0 / numberComplementarityPairs_;
     2922  const CoinWorkDouble gamma = 1.0e-8;
     2923  const CoinWorkDouble gammap = 1.0e-8;
     2924  CoinWorkDouble gammad = 1.0e-8;
     2925  int nextNumber;
     2926  int nextNumberItems;
     2927  CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
     2928  if (nextGap > bestNextGap && !allowIncreasingGap)
     2929    return false;
     2930  CoinWorkDouble lowerBoundGap = gamma * nextGap * complementarityMultiplier;
     2931  bool goodMove = true;
     2932  int iColumn;
     2933  for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
     2934    if (!flagged(iColumn)) {
     2935      if (lowerBound(iColumn)) {
     2936        CoinWorkDouble part1 = lowerSlack_[iColumn] + actualPrimalStep_ * deltaSL_[iColumn];
     2937        CoinWorkDouble part2 = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
     2938        if (part1 * part2 < lowerBoundGap) {
     2939          goodMove = false;
     2940          break;
     2941        }
     2942      }
     2943      if (upperBound(iColumn)) {
     2944        CoinWorkDouble part1 = upperSlack_[iColumn] + actualPrimalStep_ * deltaSU_[iColumn];
     2945        CoinWorkDouble part2 = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
     2946        if (part1 * part2 < lowerBoundGap) {
     2947          goodMove = false;
     2948          break;
     2949        }
     2950      }
     2951    }
     2952  }
     2953  CoinWorkDouble *nextDj = NULL;
     2954  CoinWorkDouble maximumDualError = maximumDualError_;
    29832955#ifndef NO_RTTI
    2984      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     2956  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
    29852957#else
    2986      ClpQuadraticObjective * quadraticObj = NULL;
    2987      if (objective_->type() == 2)
    2988           quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    2989 #endif
    2990      CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    2991      if (quadraticObj) {
    2992           // change gammad
    2993           gammad = 1.0e-4;
    2994           CoinWorkDouble gamma2 = gamma_ * gamma_;
    2995           nextDj = new CoinWorkDouble [numberColumns_];
    2996           CoinWorkDouble * nextSolution = new CoinWorkDouble [numberColumns_];
    2997           // put next primal into nextSolution
    2998           for ( iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2999                if (!flagged(iColumn)) {
    3000                     nextSolution[iColumn] = solution_[iColumn] +
    3001                                             actualPrimalStep_ * deltaX_[iColumn];
    3002                } else {
    3003                     nextSolution[iColumn] = solution_[iColumn];
    3004                }
    3005           }
    3006           // do reduced costs
    3007           CoinMemcpyN(cost_, numberColumns_, nextDj);
    3008           matrix_->transposeTimes(-1.0, dualArray, nextDj);
    3009           matrix_->transposeTimes(-actualDualStep_, deltaY_, nextDj);
    3010           quadraticDjs(nextDj, nextSolution, 1.0);
    3011           delete [] nextSolution;
    3012           CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    3013           const int * columnQuadraticLength = quadratic->getVectorLengths();
    3014           for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    3015                if (!fixedOrFree(iColumn)) {
    3016                     CoinWorkDouble newZ = 0.0;
    3017                     CoinWorkDouble newW = 0.0;
    3018                     if (lowerBound(iColumn)) {
    3019                          newZ = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
    3020                     }
    3021                     if (upperBound(iColumn)) {
    3022                          newW = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
    3023                     }
    3024                     if (columnQuadraticLength[iColumn]) {
    3025                          CoinWorkDouble gammaTerm = gamma2;
    3026                          if (primalR_)
    3027                               gammaTerm += primalR_[iColumn];
    3028                          //CoinWorkDouble dualInfeasibility=
    3029                          //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
    3030                          //+gammaTerm*solution_[iColumn];
    3031                          CoinWorkDouble newInfeasibility =
    3032                               nextDj[iColumn] - newZ + newW
    3033                               + gammaTerm * (solution_[iColumn] + actualPrimalStep_ * deltaX_[iColumn]);
    3034                          maximumDualError = CoinMax(maximumDualError, newInfeasibility);
    3035                          //if (CoinAbs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
    3036                          //if (dualInfeasibility*newInfeasibility<0.0) {
    3037                          //  printf("%d current %g next %g\n",iColumn,dualInfeasibility,
    3038                          //       newInfeasibility);
    3039                          //  goodMove=false;
    3040                          //}
    3041                          //}
    3042                     }
    3043                }
    3044           }
    3045           delete [] nextDj;
    3046      }
    3047 //      Satisfy g_p(alpha)?
    3048      if (rhsNorm_ > solutionNorm_) {
    3049           solutionNorm_ = rhsNorm_;
    3050      }
    3051      CoinWorkDouble errorCheck = maximumRHSError_ / solutionNorm_;
    3052      if (errorCheck < maximumBoundInfeasibility_) {
    3053           errorCheck = maximumBoundInfeasibility_;
    3054      }
    3055      // scale back move
    3056      move = CoinMin(move, 0.95);
    3057      //scale
    3058      if ((1.0 - move)*errorCheck > primalTolerance()) {
    3059           if (nextGap < gammap*(1.0 - move)*errorCheck) {
    3060                goodMove = false;
    3061           }
    3062      }
    3063      //      Satisfy g_d(alpha)?
    3064      errorCheck = maximumDualError / objectiveNorm_;
    3065      if ((1.0 - move)*errorCheck > dualTolerance()) {
    3066           if (nextGap < gammad*(1.0 - move)*errorCheck) {
    3067                goodMove = false;
    3068           }
    3069      }
    3070      if (goodMove)
    3071           bestNextGap = nextGap;
    3072      return goodMove;
     2958  ClpQuadraticObjective *quadraticObj = NULL;
     2959  if (objective_->type() == 2)
     2960    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     2961#endif
     2962  CoinWorkDouble *dualArray = reinterpret_cast< CoinWorkDouble * >(dual_);
     2963  if (quadraticObj) {
     2964    // change gammad
     2965    gammad = 1.0e-4;
     2966    CoinWorkDouble gamma2 = gamma_ * gamma_;
     2967    nextDj = new CoinWorkDouble[numberColumns_];
     2968    CoinWorkDouble *nextSolution = new CoinWorkDouble[numberColumns_];
     2969    // put next primal into nextSolution
     2970    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2971      if (!flagged(iColumn)) {
     2972        nextSolution[iColumn] = solution_[iColumn] + actualPrimalStep_ * deltaX_[iColumn];
     2973      } else {
     2974        nextSolution[iColumn] = solution_[iColumn];
     2975      }
     2976    }
     2977    // do reduced costs
     2978    CoinMemcpyN(cost_, numberColumns_, nextDj);
     2979    matrix_->transposeTimes(-1.0, dualArray, nextDj);
     2980    matrix_->transposeTimes(-actualDualStep_, deltaY_, nextDj);
     2981    quadraticDjs(nextDj, nextSolution, 1.0);
     2982    delete[] nextSolution;
     2983    CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
     2984    const int *columnQuadraticLength = quadratic->getVectorLengths();
     2985    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2986      if (!fixedOrFree(iColumn)) {
     2987        CoinWorkDouble newZ = 0.0;
     2988        CoinWorkDouble newW = 0.0;
     2989        if (lowerBound(iColumn)) {
     2990          newZ = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
     2991        }
     2992        if (upperBound(iColumn)) {
     2993          newW = wVec_[iColumn] + actualDualStep_ * de