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

formatting

File:
1 edited

Legend:

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

    r2271 r2385  
    3232#endif
    3333// primal
    34 int ClpSimplexNonlinear::primal ()
     34int ClpSimplexNonlinear::primal()
    3535{
    3636
    37      int ifValuesPass = 1;
    38      algorithm_ = +3;
    39 
    40      // save data
    41      ClpDataSave data = saveData();
    42      matrix_->refresh(this); // make sure matrix okay
    43 
    44      // Save objective
    45      ClpObjective * saveObjective = NULL;
    46      if (objective_->type() > 1) {
    47           // expand to full if quadratic
     37  int ifValuesPass = 1;
     38  algorithm_ = +3;
     39
     40  // save data
     41  ClpDataSave data = saveData();
     42  matrix_->refresh(this); // make sure matrix okay
     43
     44  // Save objective
     45  ClpObjective *saveObjective = NULL;
     46  if (objective_->type() > 1) {
     47    // expand to full if quadratic
    4848#ifndef NO_RTTI
    49           ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     49    ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
    5050#else
    51           ClpQuadraticObjective * quadraticObj = NULL;
    52           if (objective_->type() == 2)
    53                quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    54 #endif
    55           // for moment only if no scaling
    56           // May be faster if switched off - but can't see why
    57           if (!quadraticObj->fullMatrix() && (!rowScale_ && !scalingFlag_) && objectiveScale_ == 1.0) {
    58                saveObjective = objective_;
    59                objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
    60           }
    61      }
    62      double bestObjectiveWhenFlagged = COIN_DBL_MAX;
    63      int pivotMode = 15;
    64      //pivotMode=20;
    65 
    66      // initialize - maybe values pass and algorithm_ is +1
    67      // true does something (? perturbs)
    68      if (!startup(true)) {
    69 
    70           // Set average theta
    71           nonLinearCost_->setAverageTheta(1.0e3);
    72           int lastCleaned = 0; // last time objective or bounds cleaned up
    73 
    74           // Say no pivot has occurred (for steepest edge and updates)
    75           pivotRow_ = -2;
    76 
    77           // This says whether to restore things etc
    78           int factorType = 0;
    79           // Start check for cycles
    80           progress_.startCheck();
    81           /*
     51    ClpQuadraticObjective *quadraticObj = NULL;
     52    if (objective_->type() == 2)
     53      quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     54#endif
     55    // for moment only if no scaling
     56    // May be faster if switched off - but can't see why
     57    if (!quadraticObj->fullMatrix() && (!rowScale_ && !scalingFlag_) && objectiveScale_ == 1.0) {
     58      saveObjective = objective_;
     59      objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
     60    }
     61  }
     62  double bestObjectiveWhenFlagged = COIN_DBL_MAX;
     63  int pivotMode = 15;
     64  //pivotMode=20;
     65
     66  // initialize - maybe values pass and algorithm_ is +1
     67  // true does something (? perturbs)
     68  if (!startup(true)) {
     69
     70    // Set average theta
     71    nonLinearCost_->setAverageTheta(1.0e3);
     72    int lastCleaned = 0; // last time objective or bounds cleaned up
     73
     74    // Say no pivot has occurred (for steepest edge and updates)
     75    pivotRow_ = -2;
     76
     77    // This says whether to restore things etc
     78    int factorType = 0;
     79    // Start check for cycles
     80    progress_.startCheck();
     81    /*
    8282            Status of problem:
    8383            0 - optimal
     
    9090            -5 - looks unbounded
    9191          */
    92           while (problemStatus_ < 0) {
    93                int iRow, iColumn;
    94                // clear
    95                for (iRow = 0; iRow < 4; iRow++) {
    96                     rowArray_[iRow]->clear();
    97                }
    98 
    99                for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
    100                     columnArray_[iColumn]->clear();
    101                }
    102 
    103                // give matrix (and model costs and bounds a chance to be
    104                // refreshed (normally null)
    105                matrix_->refresh(this);
    106                // If getting nowhere - why not give it a kick
    107                // If we have done no iterations - special
    108                if (lastGoodIteration_ == numberIterations_ && factorType)
    109                     factorType = 3;
    110 
    111                // may factorize, checks if problem finished
    112                if (objective_->type() > 1 && lastFlaggedIteration_ >= 0 &&
    113                          numberIterations_ > lastFlaggedIteration_ + 507) {
    114                     unflag();
    115                     lastFlaggedIteration_ = numberIterations_;
    116                     if (pivotMode >= 10) {
    117                          pivotMode--;
     92    while (problemStatus_ < 0) {
     93      int iRow, iColumn;
     94      // clear
     95      for (iRow = 0; iRow < 4; iRow++) {
     96        rowArray_[iRow]->clear();
     97      }
     98
     99      for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
     100        columnArray_[iColumn]->clear();
     101      }
     102
     103      // give matrix (and model costs and bounds a chance to be
     104      // refreshed (normally null)
     105      matrix_->refresh(this);
     106      // If getting nowhere - why not give it a kick
     107      // If we have done no iterations - special
     108      if (lastGoodIteration_ == numberIterations_ && factorType)
     109        factorType = 3;
     110
     111      // may factorize, checks if problem finished
     112      if (objective_->type() > 1 && lastFlaggedIteration_ >= 0 && numberIterations_ > lastFlaggedIteration_ + 507) {
     113        unflag();
     114        lastFlaggedIteration_ = numberIterations_;
     115        if (pivotMode >= 10) {
     116          pivotMode--;
    118117#ifdef CLP_DEBUG
    119                          if (handler_->logLevel() & 32)
    120                               printf("pivot mode now %d\n", pivotMode);
    121 #endif
    122                          if (pivotMode == 9)
    123                               pivotMode = 0;    // switch off fast attempt
    124                     }
    125                }
    126                statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true,
    127                                        bestObjectiveWhenFlagged);
    128 
    129                // Say good factorization
    130                factorType = 1;
    131 
    132                // Say no pivot has occurred (for steepest edge and updates)
    133                pivotRow_ = -2;
    134 
    135                // exit if victory declared
    136                if (problemStatus_ >= 0)
    137                     break;
    138 
    139                // test for maximum iterations
    140                if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
    141                     problemStatus_ = 3;
    142                     break;
    143                }
    144 
    145                if (firstFree_ < 0) {
    146                     if (ifValuesPass) {
    147                          // end of values pass
    148                          ifValuesPass = 0;
    149                          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
    150                          if (status >= 0) {
    151                               problemStatus_ = 5;
    152                               secondaryStatus_ = ClpEventHandler::endOfValuesPass;
    153                               break;
    154                          }
    155                     }
    156                }
    157                // Check event
    158                {
    159                     int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
    160                     if (status >= 0) {
    161                          problemStatus_ = 5;
    162                          secondaryStatus_ = ClpEventHandler::endOfFactorization;
    163                          break;
    164                     }
    165                }
    166                // Iterate
    167                whileIterating(pivotMode);
    168           }
    169      }
    170      // if infeasible get real values
    171      if (problemStatus_ == 1) {
    172           infeasibilityCost_ = 0.0;
    173           createRim(1 + 4);
    174           delete nonLinearCost_;
    175           nonLinearCost_ = new ClpNonLinearCost(this);
    176           nonLinearCost_->checkInfeasibilities(0.0);
    177           sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
    178           numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
    179           // and get good feasible duals
    180           computeDuals(NULL);
    181      }
    182      // correct objective value
    183      if (numberColumns_)
    184           objectiveValue_ = nonLinearCost_->feasibleCost() + objective_->nonlinearOffset();
    185      objectiveValue_ /= (objectiveScale_ * rhsScale_);
    186      // clean up
    187      unflag();
    188      finish();
    189      restoreData(data);
    190      // restore objective if full
    191      if (saveObjective) {
    192           delete objective_;
    193           objective_ = saveObjective;
    194      }
    195      return problemStatus_;
     118          if (handler_->logLevel() & 32)
     119            printf("pivot mode now %d\n", pivotMode);
     120#endif
     121          if (pivotMode == 9)
     122            pivotMode = 0; // switch off fast attempt
     123        }
     124      }
     125      statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true,
     126        bestObjectiveWhenFlagged);
     127
     128      // Say good factorization
     129      factorType = 1;
     130
     131      // Say no pivot has occurred (for steepest edge and updates)
     132      pivotRow_ = -2;
     133
     134      // exit if victory declared
     135      if (problemStatus_ >= 0)
     136        break;
     137
     138      // test for maximum iterations
     139      if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
     140        problemStatus_ = 3;
     141        break;
     142      }
     143
     144      if (firstFree_ < 0) {
     145        if (ifValuesPass) {
     146          // end of values pass
     147          ifValuesPass = 0;
     148          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
     149          if (status >= 0) {
     150            problemStatus_ = 5;
     151            secondaryStatus_ = ClpEventHandler::endOfValuesPass;
     152            break;
     153          }
     154        }
     155      }
     156      // Check event
     157      {
     158        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
     159        if (status >= 0) {
     160          problemStatus_ = 5;
     161          secondaryStatus_ = ClpEventHandler::endOfFactorization;
     162          break;
     163        }
     164      }
     165      // Iterate
     166      whileIterating(pivotMode);
     167    }
     168  }
     169  // if infeasible get real values
     170  if (problemStatus_ == 1) {
     171    infeasibilityCost_ = 0.0;
     172    createRim(1 + 4);
     173    delete nonLinearCost_;
     174    nonLinearCost_ = new ClpNonLinearCost(this);
     175    nonLinearCost_->checkInfeasibilities(0.0);
     176    sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
     177    numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
     178    // and get good feasible duals
     179    computeDuals(NULL);
     180  }
     181  // correct objective value
     182  if (numberColumns_)
     183    objectiveValue_ = nonLinearCost_->feasibleCost() + objective_->nonlinearOffset();
     184  objectiveValue_ /= (objectiveScale_ * rhsScale_);
     185  // clean up
     186  unflag();
     187  finish();
     188  restoreData(data);
     189  // restore objective if full
     190  if (saveObjective) {
     191    delete objective_;
     192    objective_ = saveObjective;
     193  }
     194  return problemStatus_;
    196195}
    197196/*  Refactorizes if necessary
     
    204203    - 2 restoring from saved
    205204*/
    206 void
    207 ClpSimplexNonlinear::statusOfProblemInPrimal(int & lastCleaned, int type,
    208           ClpSimplexProgress * progress,
    209           bool doFactorization,
    210           double & bestObjectiveWhenFlagged)
     205void ClpSimplexNonlinear::statusOfProblemInPrimal(int &lastCleaned, int type,
     206  ClpSimplexProgress *progress,
     207  bool doFactorization,
     208  double &bestObjectiveWhenFlagged)
    211209{
    212      int dummy; // for use in generalExpanded
    213      if (type == 2) {
    214           // trouble - restore solution
    215           CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_ );
    216           CoinMemcpyN(savedSolution_ + numberColumns_ , numberRows_, rowActivityWork_);
    217           CoinMemcpyN(savedSolution_ ,  numberColumns_, columnActivityWork_);
    218           // restore extra stuff
    219           matrix_->generalExpanded(this, 6, dummy);
    220           forceFactorization_ = 1; // a bit drastic but ..
    221           pivotRow_ = -1; // say no weights update
     210  int dummy; // for use in generalExpanded
     211  if (type == 2) {
     212    // trouble - restore solution
     213    CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_);
     214    CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_);
     215    CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_);
     216    // restore extra stuff
     217    matrix_->generalExpanded(this, 6, dummy);
     218    forceFactorization_ = 1; // a bit drastic but ..
     219    pivotRow_ = -1; // say no weights update
     220    changeMade_++; // say change made
     221  }
     222  int saveThreshold = factorization_->sparseThreshold();
     223  int tentativeStatus = problemStatus_;
     224  int numberThrownOut = 1; // to loop round on bad factorization in values pass
     225  while (numberThrownOut) {
     226    if (problemStatus_ > -3 || problemStatus_ == -4) {
     227      // factorize
     228      // later on we will need to recover from singularities
     229      // also we could skip if first time
     230      // do weights
     231      // This may save pivotRow_ for use
     232      if (doFactorization)
     233        primalColumnPivot_->saveWeights(this, 1);
     234
     235      if (type && doFactorization) {
     236        // is factorization okay?
     237        int factorStatus = internalFactorize(1);
     238        if (factorStatus) {
     239          if (type != 1 || largestPrimalError_ > 1.0e3
     240            || largestDualError_ > 1.0e3) {
     241            // was ||largestDualError_>1.0e3||objective_->type()>1) {
     242            // switch off dense
     243            int saveDense = factorization_->denseThreshold();
     244            factorization_->setDenseThreshold(0);
     245            // make sure will do safe factorization
     246            pivotVariable_[0] = -1;
     247            internalFactorize(2);
     248            factorization_->setDenseThreshold(saveDense);
     249            // Go to safe
     250            factorization_->pivotTolerance(0.99);
     251            // restore extra stuff
     252            matrix_->generalExpanded(this, 6, dummy);
     253          } else {
     254            // no - restore previous basis
     255            CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_);
     256            CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_);
     257            CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_);
     258            // restore extra stuff
     259            matrix_->generalExpanded(this, 6, dummy);
     260            matrix_->generalExpanded(this, 5, dummy);
     261            forceFactorization_ = 1; // a bit drastic but ..
     262            type = 2;
     263            // Go to safe
     264            factorization_->pivotTolerance(0.99);
     265            if (internalFactorize(1) != 0)
     266              largestPrimalError_ = 1.0e4; // force other type
     267          }
     268          // flag incoming
     269          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
     270            setFlagged(sequenceIn_);
     271            saveStatus_[sequenceIn_] = status_[sequenceIn_];
     272          }
    222273          changeMade_++; // say change made
    223      }
    224      int saveThreshold = factorization_->sparseThreshold();
    225      int tentativeStatus = problemStatus_;
    226      int numberThrownOut = 1; // to loop round on bad factorization in values pass
    227      while (numberThrownOut) {
    228           if (problemStatus_ > -3 || problemStatus_ == -4) {
    229                // factorize
    230                // later on we will need to recover from singularities
    231                // also we could skip if first time
    232                // do weights
    233                // This may save pivotRow_ for use
    234                if (doFactorization)
    235                     primalColumnPivot_->saveWeights(this, 1);
    236 
    237                if (type && doFactorization) {
    238                     // is factorization okay?
    239                     int factorStatus = internalFactorize(1);
    240                     if (factorStatus) {
    241                          if (type != 1 || largestPrimalError_ > 1.0e3
    242                                    || largestDualError_ > 1.0e3) {
    243                               // was ||largestDualError_>1.0e3||objective_->type()>1) {
    244                               // switch off dense
    245                               int saveDense = factorization_->denseThreshold();
    246                               factorization_->setDenseThreshold(0);
    247                               // make sure will do safe factorization
    248                               pivotVariable_[0] = -1;
    249                               internalFactorize(2);
    250                               factorization_->setDenseThreshold(saveDense);
    251                               // Go to safe
    252                               factorization_->pivotTolerance(0.99);
    253                               // restore extra stuff
    254                               matrix_->generalExpanded(this, 6, dummy);
    255                          } else {
    256                               // no - restore previous basis
    257                               CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_ );
    258                               CoinMemcpyN(savedSolution_ + numberColumns_ ,             numberRows_, rowActivityWork_);
    259                               CoinMemcpyN(savedSolution_ ,              numberColumns_, columnActivityWork_);
    260                               // restore extra stuff
    261                               matrix_->generalExpanded(this, 6, dummy);
    262                               matrix_->generalExpanded(this, 5, dummy);
    263                               forceFactorization_ = 1; // a bit drastic but ..
    264                               type = 2;
    265                               // Go to safe
    266                               factorization_->pivotTolerance(0.99);
    267                               if (internalFactorize(1) != 0)
    268                                    largestPrimalError_ = 1.0e4; // force other type
    269                          }
    270                          // flag incoming
    271                          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
    272                               setFlagged(sequenceIn_);
    273                               saveStatus_[sequenceIn_] = status_[sequenceIn_];
    274                          }
    275                          changeMade_++; // say change made
    276                     }
    277                }
    278                if (problemStatus_ != -4)
    279                     problemStatus_ = -3;
    280           }
    281           // at this stage status is -3 or -5 if looks unbounded
    282           // get primal and dual solutions
    283           // put back original costs and then check
    284           createRim(4);
    285           // May need to do more if column generation
    286           dummy = 4;
    287           matrix_->generalExpanded(this, 9, dummy);
    288           numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
    289           if (numberThrownOut) {
    290                problemStatus_ = tentativeStatus;
    291                doFactorization = true;
    292           }
    293      }
    294      // Double check reduced costs if no action
    295      if (progress->lastIterationNumber(0) == numberIterations_) {
    296           if (primalColumnPivot_->looksOptimal()) {
    297                numberDualInfeasibilities_ = 0;
    298                sumDualInfeasibilities_ = 0.0;
    299           }
    300      }
    301      // Check if looping
    302      int loop;
    303      if (type != 2)
    304           loop = progress->looping();
    305      else
    306           loop = -1;
    307      if (loop >= 0) {
    308           if (!problemStatus_) {
    309                // declaring victory
    310                numberPrimalInfeasibilities_ = 0;
    311                sumPrimalInfeasibilities_ = 0.0;
    312           } else {
    313                problemStatus_ = loop; //exit if in loop
    314                problemStatus_ = 10; // instead - try other algorithm
    315           }
    316           problemStatus_ = 10; // instead - try other algorithm
    317           return ;
    318      } else if (loop < -1) {
    319           // Is it time for drastic measures
    320           if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 &&
    321                     progress->oddState() < 10 && progress->oddState() >= 0) {
    322                progress->newOddState();
    323                nonLinearCost_->zapCosts();
    324           }
    325           // something may have changed
    326           gutsOfSolution(NULL, NULL, true);
    327      }
    328      // If progress then reset costs
    329      if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
    330           createRim(4, false); // costs back
    331           delete nonLinearCost_;
    332           nonLinearCost_ = new ClpNonLinearCost(this);
    333           progress->endOddState();
    334           gutsOfSolution(NULL, NULL, true);
    335      }
    336      // Flag to say whether to go to dual to clean up
    337      //bool goToDual = false;
    338      // really for free variables in
    339      //if((progressFlag_&2)!=0)
    340      //problemStatus_=-1;;
    341      progressFlag_ = 0; //reset progress flag
    342 
    343      handler_->message(CLP_SIMPLEX_STATUS, messages_)
    344                << numberIterations_ << nonLinearCost_->feasibleReportCost();
    345      handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
    346                << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
    347      handler_->printing(sumDualInfeasibilities_ > 0.0)
    348                << sumDualInfeasibilities_ << numberDualInfeasibilities_;
    349      handler_->printing(numberDualInfeasibilitiesWithoutFree_
    350                         < numberDualInfeasibilities_)
    351                << numberDualInfeasibilitiesWithoutFree_;
    352      handler_->message() << CoinMessageEol;
    353      if (!primalFeasible()) {
    354           nonLinearCost_->checkInfeasibilities(primalTolerance_);
    355           gutsOfSolution(NULL, NULL, true);
    356           nonLinearCost_->checkInfeasibilities(primalTolerance_);
    357      }
    358      double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
    359      if (trueInfeasibility > 1.0) {
    360           // If infeasibility going up may change weights
    361           double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
    362           if(progress->lastInfeasibility() < testValue) {
    363                if (infeasibilityCost_ < 1.0e14) {
    364                     infeasibilityCost_ *= 1.5;
    365                     if (handler_->logLevel() == 63)
    366                          printf("increasing weight to %g\n", infeasibilityCost_);
    367                     gutsOfSolution(NULL, NULL, true);
    368                }
    369           }
    370      }
    371      // we may wish to say it is optimal even if infeasible
    372      bool alwaysOptimal = (specialOptions_ & 1) != 0;
    373      // give code benefit of doubt
    374      if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
    375                sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
    376           // say optimal (with these bounds etc)
    377           numberDualInfeasibilities_ = 0;
    378           sumDualInfeasibilities_ = 0.0;
    379           numberPrimalInfeasibilities_ = 0;
    380           sumPrimalInfeasibilities_ = 0.0;
    381      }
    382      // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
    383      if (dualFeasible() || problemStatus_ == -4) {
    384           // see if extra helps
    385           if (nonLinearCost_->numberInfeasibilities() &&
    386                     (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
    387                     && !alwaysOptimal) {
    388                //may need infeasiblity cost changed
    389                // we can see if we can construct a ray
    390                // make up a new objective
    391                double saveWeight = infeasibilityCost_;
    392                // save nonlinear cost as we are going to switch off costs
    393                ClpNonLinearCost * nonLinear = nonLinearCost_;
    394                // do twice to make sure Primal solution has settled
    395                // put non-basics to bounds in case tolerance moved
    396                // put back original costs
    397                createRim(4);
    398                nonLinearCost_->checkInfeasibilities(0.0);
    399                gutsOfSolution(NULL, NULL, true);
    400 
    401                infeasibilityCost_ = 1.0e100;
    402                // put back original costs
    403                createRim(4);
    404                nonLinearCost_->checkInfeasibilities(primalTolerance_);
    405                // may have fixed infeasibilities - double check
    406                if (nonLinearCost_->numberInfeasibilities() == 0) {
    407                     // carry on
    408                     problemStatus_ = -1;
    409                     infeasibilityCost_ = saveWeight;
    410                     nonLinearCost_->checkInfeasibilities(primalTolerance_);
    411                } else {
    412                     nonLinearCost_ = NULL;
    413                     // scale
    414                     int i;
    415                     for (i = 0; i < numberRows_ + numberColumns_; i++)
    416                          cost_[i] *= 1.0e-95;
    417                     gutsOfSolution(NULL, NULL, false);
    418                     nonLinearCost_ = nonLinear;
    419                     infeasibilityCost_ = saveWeight;
    420                     if ((infeasibilityCost_ >= 1.0e18 ||
    421                               numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
    422                          /* goToDual = unPerturb(); // stop any further perturbation
     274        }
     275      }
     276      if (problemStatus_ != -4)
     277        problemStatus_ = -3;
     278    }
     279    // at this stage status is -3 or -5 if looks unbounded
     280    // get primal and dual solutions
     281    // put back original costs and then check
     282    createRim(4);
     283    // May need to do more if column generation
     284    dummy = 4;
     285    matrix_->generalExpanded(this, 9, dummy);
     286    numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
     287    if (numberThrownOut) {
     288      problemStatus_ = tentativeStatus;
     289      doFactorization = true;
     290    }
     291  }
     292  // Double check reduced costs if no action
     293  if (progress->lastIterationNumber(0) == numberIterations_) {
     294    if (primalColumnPivot_->looksOptimal()) {
     295      numberDualInfeasibilities_ = 0;
     296      sumDualInfeasibilities_ = 0.0;
     297    }
     298  }
     299  // Check if looping
     300  int loop;
     301  if (type != 2)
     302    loop = progress->looping();
     303  else
     304    loop = -1;
     305  if (loop >= 0) {
     306    if (!problemStatus_) {
     307      // declaring victory
     308      numberPrimalInfeasibilities_ = 0;
     309      sumPrimalInfeasibilities_ = 0.0;
     310    } else {
     311      problemStatus_ = loop; //exit if in loop
     312      problemStatus_ = 10; // instead - try other algorithm
     313    }
     314    problemStatus_ = 10; // instead - try other algorithm
     315    return;
     316  } else if (loop < -1) {
     317    // Is it time for drastic measures
     318    if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 && progress->oddState() < 10 && progress->oddState() >= 0) {
     319      progress->newOddState();
     320      nonLinearCost_->zapCosts();
     321    }
     322    // something may have changed
     323    gutsOfSolution(NULL, NULL, true);
     324  }
     325  // If progress then reset costs
     326  if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
     327    createRim(4, false); // costs back
     328    delete nonLinearCost_;
     329    nonLinearCost_ = new ClpNonLinearCost(this);
     330    progress->endOddState();
     331    gutsOfSolution(NULL, NULL, true);
     332  }
     333  // Flag to say whether to go to dual to clean up
     334  //bool goToDual = false;
     335  // really for free variables in
     336  //if((progressFlag_&2)!=0)
     337  //problemStatus_=-1;;
     338  progressFlag_ = 0; //reset progress flag
     339
     340  handler_->message(CLP_SIMPLEX_STATUS, messages_)
     341    << numberIterations_ << nonLinearCost_->feasibleReportCost();
     342  handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
     343    << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
     344  handler_->printing(sumDualInfeasibilities_ > 0.0)
     345    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
     346  handler_->printing(numberDualInfeasibilitiesWithoutFree_
     347    < numberDualInfeasibilities_)
     348    << numberDualInfeasibilitiesWithoutFree_;
     349  handler_->message() << CoinMessageEol;
     350  if (!primalFeasible()) {
     351    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     352    gutsOfSolution(NULL, NULL, true);
     353    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     354  }
     355  double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
     356  if (trueInfeasibility > 1.0) {
     357    // If infeasibility going up may change weights
     358    double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
     359    if (progress->lastInfeasibility() < testValue) {
     360      if (infeasibilityCost_ < 1.0e14) {
     361        infeasibilityCost_ *= 1.5;
     362        if (handler_->logLevel() == 63)
     363          printf("increasing weight to %g\n", infeasibilityCost_);
     364        gutsOfSolution(NULL, NULL, true);
     365      }
     366    }
     367  }
     368  // we may wish to say it is optimal even if infeasible
     369  bool alwaysOptimal = (specialOptions_ & 1) != 0;
     370  // give code benefit of doubt
     371  if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     372    // say optimal (with these bounds etc)
     373    numberDualInfeasibilities_ = 0;
     374    sumDualInfeasibilities_ = 0.0;
     375    numberPrimalInfeasibilities_ = 0;
     376    sumPrimalInfeasibilities_ = 0.0;
     377  }
     378  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
     379  if (dualFeasible() || problemStatus_ == -4) {
     380    // see if extra helps
     381    if (nonLinearCost_->numberInfeasibilities() && (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
     382      && !alwaysOptimal) {
     383      //may need infeasiblity cost changed
     384      // we can see if we can construct a ray
     385      // make up a new objective
     386      double saveWeight = infeasibilityCost_;
     387      // save nonlinear cost as we are going to switch off costs
     388      ClpNonLinearCost *nonLinear = nonLinearCost_;
     389      // do twice to make sure Primal solution has settled
     390      // put non-basics to bounds in case tolerance moved
     391      // put back original costs
     392      createRim(4);
     393      nonLinearCost_->checkInfeasibilities(0.0);
     394      gutsOfSolution(NULL, NULL, true);
     395
     396      infeasibilityCost_ = 1.0e100;
     397      // put back original costs
     398      createRim(4);
     399      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     400      // may have fixed infeasibilities - double check
     401      if (nonLinearCost_->numberInfeasibilities() == 0) {
     402        // carry on
     403        problemStatus_ = -1;
     404        infeasibilityCost_ = saveWeight;
     405        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     406      } else {
     407        nonLinearCost_ = NULL;
     408        // scale
     409        int i;
     410        for (i = 0; i < numberRows_ + numberColumns_; i++)
     411          cost_[i] *= 1.0e-95;
     412        gutsOfSolution(NULL, NULL, false);
     413        nonLinearCost_ = nonLinear;
     414        infeasibilityCost_ = saveWeight;
     415        if ((infeasibilityCost_ >= 1.0e18 || numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
     416          /* goToDual = unPerturb(); // stop any further perturbation
    423417                         if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
    424418                              goToDual = false;
    425419                         */
    426                          unPerturb();
    427                          nonLinearCost_->checkInfeasibilities(primalTolerance_);
    428                          numberDualInfeasibilities_ = 1; // carry on
    429                          problemStatus_ = -1;
    430                     }
    431                     if (infeasibilityCost_ >= 1.0e20 ||
    432                               numberDualInfeasibilities_ == 0) {
    433                          // we are infeasible - use as ray
    434                          delete [] ray_;
    435                          ray_ = new double [numberRows_];
    436                          CoinMemcpyN(dual_, numberRows_, ray_);
    437                          // and get feasible duals
    438                          infeasibilityCost_ = 0.0;
    439                          createRim(4);
    440                          nonLinearCost_->checkInfeasibilities(primalTolerance_);
    441                          gutsOfSolution(NULL, NULL, true);
    442                          // so will exit
    443                          infeasibilityCost_ = 1.0e30;
    444                          // reset infeasibilities
    445                          sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();;
    446                          numberPrimalInfeasibilities_ =
    447                               nonLinearCost_->numberInfeasibilities();
    448                     }
    449                     if (infeasibilityCost_ < 1.0e20) {
    450                          infeasibilityCost_ *= 5.0;
    451                          changeMade_++; // say change made
    452                          unflag();
    453                          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
    454                                    << infeasibilityCost_
    455                                    << CoinMessageEol;
    456                          // put back original costs and then check
    457                          createRim(4);
    458                          nonLinearCost_->checkInfeasibilities(0.0);
    459                          gutsOfSolution(NULL, NULL, true);
    460                          problemStatus_ = -1; //continue
    461                          //goToDual = false;
    462                     } else {
    463                          // say infeasible
    464                          problemStatus_ = 1;
    465                     }
    466                }
     420          unPerturb();
     421          nonLinearCost_->checkInfeasibilities(primalTolerance_);
     422          numberDualInfeasibilities_ = 1; // carry on
     423          problemStatus_ = -1;
     424        }
     425        if (infeasibilityCost_ >= 1.0e20 || numberDualInfeasibilities_ == 0) {
     426          // we are infeasible - use as ray
     427          delete[] ray_;
     428          ray_ = new double[numberRows_];
     429          CoinMemcpyN(dual_, numberRows_, ray_);
     430          // and get feasible duals
     431          infeasibilityCost_ = 0.0;
     432          createRim(4);
     433          nonLinearCost_->checkInfeasibilities(primalTolerance_);
     434          gutsOfSolution(NULL, NULL, true);
     435          // so will exit
     436          infeasibilityCost_ = 1.0e30;
     437          // reset infeasibilities
     438          sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
     439          ;
     440          numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
     441        }
     442        if (infeasibilityCost_ < 1.0e20) {
     443          infeasibilityCost_ *= 5.0;
     444          changeMade_++; // say change made
     445          unflag();
     446          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
     447            << infeasibilityCost_
     448            << CoinMessageEol;
     449          // put back original costs and then check
     450          createRim(4);
     451          nonLinearCost_->checkInfeasibilities(0.0);
     452          gutsOfSolution(NULL, NULL, true);
     453          problemStatus_ = -1; //continue
     454          //goToDual = false;
     455        } else {
     456          // say infeasible
     457          problemStatus_ = 1;
     458        }
     459      }
     460    } else {
     461      // may be optimal
     462      if (perturbation_ == 101) {
     463        /* goToDual =*/unPerturb(); // stop any further perturbation
     464        lastCleaned = -1; // carry on
     465      }
     466      bool unflagged = unflag() != 0;
     467      if (lastCleaned != numberIterations_ || unflagged) {
     468        handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
     469          << primalTolerance_
     470          << CoinMessageEol;
     471        double current = nonLinearCost_->feasibleReportCost();
     472        if (numberTimesOptimal_ < 4) {
     473          if (bestObjectiveWhenFlagged <= current) {
     474            numberTimesOptimal_++;
     475#ifdef CLP_DEBUG
     476            if (handler_->logLevel() & 32)
     477              printf("%d times optimal, current %g best %g\n", numberTimesOptimal_,
     478                current, bestObjectiveWhenFlagged);
     479#endif
    467480          } else {
    468                // may be optimal
    469                if (perturbation_ == 101) {
    470                     /* goToDual =*/ unPerturb(); // stop any further perturbation
    471                     lastCleaned = -1; // carry on
    472                }
    473                bool unflagged = unflag() != 0;
    474                if ( lastCleaned != numberIterations_ || unflagged) {
    475                     handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
    476                               << primalTolerance_
    477                               << CoinMessageEol;
    478                     double current = nonLinearCost_->feasibleReportCost();
    479                     if (numberTimesOptimal_ < 4) {
    480                          if (bestObjectiveWhenFlagged <= current) {
    481                               numberTimesOptimal_++;
    482 #ifdef CLP_DEBUG
    483                               if (handler_->logLevel() & 32)
    484                                    printf("%d times optimal, current %g best %g\n", numberTimesOptimal_,
    485                                           current, bestObjectiveWhenFlagged);
    486 #endif
    487                          } else {
    488                               bestObjectiveWhenFlagged = current;
    489                          }
    490                          changeMade_++; // say change made
    491                          if (numberTimesOptimal_ == 1) {
    492                               // better to have small tolerance even if slower
    493                               factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
    494                          }
    495                          lastCleaned = numberIterations_;
    496                          if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
    497                               handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
    498                                         << CoinMessageEol;
    499                          double oldTolerance = primalTolerance_;
    500                          primalTolerance_ = dblParam_[ClpPrimalTolerance];
    501                          // put back original costs and then check
    502                          createRim(4);
    503                          nonLinearCost_->checkInfeasibilities(oldTolerance);
    504                          gutsOfSolution(NULL, NULL, true);
    505                          if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
    506                                    sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
    507                               // say optimal (with these bounds etc)
    508                               numberDualInfeasibilities_ = 0;
    509                               sumDualInfeasibilities_ = 0.0;
    510                               numberPrimalInfeasibilities_ = 0;
    511                               sumPrimalInfeasibilities_ = 0.0;
    512                          }
    513                          if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
    514                               problemStatus_ = 0;
    515                          else
    516                               problemStatus_ = -1;
    517                     } else {
    518                          problemStatus_ = 0; // optimal
    519                          if (lastCleaned < numberIterations_) {
    520                               handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
    521                                         << CoinMessageEol;
    522                          }
    523                     }
    524                } else {
    525                     problemStatus_ = 0; // optimal
    526                }
    527           }
    528      } else {
    529           // see if looks unbounded
    530           if (problemStatus_ == -5) {
    531                if (nonLinearCost_->numberInfeasibilities()) {
    532                     if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) {
    533                          // back off weight
    534                          infeasibilityCost_ = 1.0e13;
    535                          unPerturb(); // stop any further perturbation
    536                     }
    537                     //we need infeasiblity cost changed
    538                     if (infeasibilityCost_ < 1.0e20) {
    539                          infeasibilityCost_ *= 5.0;
    540                          changeMade_++; // say change made
    541                          unflag();
    542                          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
    543                                    << infeasibilityCost_
    544                                    << CoinMessageEol;
    545                          // put back original costs and then check
    546                          createRim(4);
    547                          gutsOfSolution(NULL, NULL, true);
    548                          problemStatus_ = -1; //continue
    549                     } else {
    550                          // say unbounded
    551                          problemStatus_ = 2;
    552                     }
    553                } else {
    554                     // say unbounded
    555                     problemStatus_ = 2;
    556                }
    557           } else {
    558                if(type == 3 && problemStatus_ != -5)
    559                     unflag(); // odd
    560                // carry on
    561                problemStatus_ = -1;
    562           }
    563      }
    564      // save extra stuff
    565      matrix_->generalExpanded(this, 5, dummy);
    566      if (type == 0 || type == 1) {
    567           if (type != 1 || !saveStatus_) {
    568                // create save arrays
    569                delete [] saveStatus_;
    570                delete [] savedSolution_;
    571                saveStatus_ = new unsigned char [numberRows_+numberColumns_];
    572                savedSolution_ = new double [numberRows_+numberColumns_];
    573           }
    574           // save arrays
    575           CoinMemcpyN(status_, (numberColumns_ + numberRows_), saveStatus_);
    576           CoinMemcpyN(rowActivityWork_, numberRows_, savedSolution_ + numberColumns_ );
    577           CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_ );
    578      }
    579      if (doFactorization) {
    580           // restore weights (if saved) - also recompute infeasibility list
    581           if (tentativeStatus > -3)
    582                primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
     481            bestObjectiveWhenFlagged = current;
     482          }
     483          changeMade_++; // say change made
     484          if (numberTimesOptimal_ == 1) {
     485            // better to have small tolerance even if slower
     486            factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
     487          }
     488          lastCleaned = numberIterations_;
     489          if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
     490            handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
     491              << CoinMessageEol;
     492          double oldTolerance = primalTolerance_;
     493          primalTolerance_ = dblParam_[ClpPrimalTolerance];
     494          // put back original costs and then check
     495          createRim(4);
     496          nonLinearCost_->checkInfeasibilities(oldTolerance);
     497          gutsOfSolution(NULL, NULL, true);
     498          if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     499            // say optimal (with these bounds etc)
     500            numberDualInfeasibilities_ = 0;
     501            sumDualInfeasibilities_ = 0.0;
     502            numberPrimalInfeasibilities_ = 0;
     503            sumPrimalInfeasibilities_ = 0.0;
     504          }
     505          if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
     506            problemStatus_ = 0;
    583507          else
    584                primalColumnPivot_->saveWeights(this, 3);
    585           if (saveThreshold) {
    586                // use default at present
    587                factorization_->sparseThreshold(0);
    588                factorization_->goSparse();
    589           }
    590      }
    591      if (problemStatus_ < 0 && !changeMade_) {
    592           problemStatus_ = 4; // unknown
    593      }
    594      lastGoodIteration_ = numberIterations_;
    595      //if (goToDual)
    596      //problemStatus_=10; // try dual
    597      // Allow matrices to be sorted etc
    598      int fake = -999; // signal sort
    599      matrix_->correctSequence(this, fake, fake);
     508            problemStatus_ = -1;
     509        } else {
     510          problemStatus_ = 0; // optimal
     511          if (lastCleaned < numberIterations_) {
     512            handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
     513              << CoinMessageEol;
     514          }
     515        }
     516      } else {
     517        problemStatus_ = 0; // optimal
     518      }
     519    }
     520  } else {
     521    // see if looks unbounded
     522    if (problemStatus_ == -5) {
     523      if (nonLinearCost_->numberInfeasibilities()) {
     524        if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) {
     525          // back off weight
     526          infeasibilityCost_ = 1.0e13;
     527          unPerturb(); // stop any further perturbation
     528        }
     529        //we need infeasiblity cost changed
     530        if (infeasibilityCost_ < 1.0e20) {
     531          infeasibilityCost_ *= 5.0;
     532          changeMade_++; // say change made
     533          unflag();
     534          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
     535            << infeasibilityCost_
     536            << CoinMessageEol;
     537          // put back original costs and then check
     538          createRim(4);
     539          gutsOfSolution(NULL, NULL, true);
     540          problemStatus_ = -1; //continue
     541        } else {
     542          // say unbounded
     543          problemStatus_ = 2;
     544        }
     545      } else {
     546        // say unbounded
     547        problemStatus_ = 2;
     548      }
     549    } else {
     550      if (type == 3 && problemStatus_ != -5)
     551        unflag(); // odd
     552      // carry on
     553      problemStatus_ = -1;
     554    }
     555  }
     556  // save extra stuff
     557  matrix_->generalExpanded(this, 5, dummy);
     558  if (type == 0 || type == 1) {
     559    if (type != 1 || !saveStatus_) {
     560      // create save arrays
     561      delete[] saveStatus_;
     562      delete[] savedSolution_;
     563      saveStatus_ = new unsigned char[numberRows_ + numberColumns_];
     564      savedSolution_ = new double[numberRows_ + numberColumns_];
     565    }
     566    // save arrays
     567    CoinMemcpyN(status_, (numberColumns_ + numberRows_), saveStatus_);
     568    CoinMemcpyN(rowActivityWork_, numberRows_, savedSolution_ + numberColumns_);
     569    CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
     570  }
     571  if (doFactorization) {
     572    // restore weights (if saved) - also recompute infeasibility list
     573    if (tentativeStatus > -3)
     574      primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
     575    else
     576      primalColumnPivot_->saveWeights(this, 3);
     577    if (saveThreshold) {
     578      // use default at present
     579      factorization_->sparseThreshold(0);
     580      factorization_->goSparse();
     581    }
     582  }
     583  if (problemStatus_ < 0 && !changeMade_) {
     584    problemStatus_ = 4; // unknown
     585  }
     586  lastGoodIteration_ = numberIterations_;
     587  //if (goToDual)
     588  //problemStatus_=10; // try dual
     589  // Allow matrices to be sorted etc
     590  int fake = -999; // signal sort
     591  matrix_->correctSequence(this, fake, fake);
    600592}
    601593/*
     
    609601  +3 max iterations
    610602*/
    611 int
    612 ClpSimplexNonlinear::whileIterating(int & pivotMode)
     603int ClpSimplexNonlinear::whileIterating(int &pivotMode)
    613604{
    614      // Say if values pass
    615      //int ifValuesPass=(firstFree_>=0) ? 1 : 0;
    616      int ifValuesPass = 1;
    617      int returnCode = -1;
    618      // status stays at -1 while iterating, >=0 finished, -2 to invert
    619      // status -3 to go to top without an invert
    620      int numberInterior = 0;
    621      int nextUnflag = 10;
    622      int nextUnflagIteration = numberIterations_ + 10;
    623      // get two arrays
    624      double * array1 = new double[2*(numberRows_+numberColumns_)];
    625      double solutionError = -1.0;
    626      while (problemStatus_ == -1) {
    627           int result;
    628           rowArray_[1]->clear();
    629           //#define CLP_DEBUG
     605  // Say if values pass
     606  //int ifValuesPass=(firstFree_>=0) ? 1 : 0;
     607  int ifValuesPass = 1;
     608  int returnCode = -1;
     609  // status stays at -1 while iterating, >=0 finished, -2 to invert
     610  // status -3 to go to top without an invert
     611  int numberInterior = 0;
     612  int nextUnflag = 10;
     613  int nextUnflagIteration = numberIterations_ + 10;
     614  // get two arrays
     615  double *array1 = new double[2 * (numberRows_ + numberColumns_)];
     616  double solutionError = -1.0;
     617  while (problemStatus_ == -1) {
     618    int result;
     619    rowArray_[1]->clear();
     620    //#define CLP_DEBUG
    630621#if CLP_DEBUG > 1
    631           rowArray_[0]->checkClear();
    632           rowArray_[1]->checkClear();
    633           rowArray_[2]->checkClear();
    634           rowArray_[3]->checkClear();
    635           columnArray_[0]->checkClear();
    636 #endif
    637           if (numberInterior >= 5) {
    638                // this can go when we have better minimization
    639                if (pivotMode < 10)
    640                     pivotMode = 1;
    641                unflag();
     622    rowArray_[0]->checkClear();
     623    rowArray_[1]->checkClear();
     624    rowArray_[2]->checkClear();
     625    rowArray_[3]->checkClear();
     626    columnArray_[0]->checkClear();
     627#endif
     628    if (numberInterior >= 5) {
     629      // this can go when we have better minimization
     630      if (pivotMode < 10)
     631        pivotMode = 1;
     632      unflag();
    642633#ifdef CLP_DEBUG
    643                if (handler_->logLevel() & 32)
    644                     printf("interior unflag\n");
    645 #endif
    646                numberInterior = 0;
    647                nextUnflag = 10;
    648                nextUnflagIteration = numberIterations_ + 10;
    649           } else {
    650                if (numberInterior > nextUnflag &&
    651                          numberIterations_ > nextUnflagIteration) {
    652                     nextUnflagIteration = numberIterations_ + 10;
    653                     nextUnflag += 10;
    654                     unflag();
     634      if (handler_->logLevel() & 32)
     635        printf("interior unflag\n");
     636#endif
     637      numberInterior = 0;
     638      nextUnflag = 10;
     639      nextUnflagIteration = numberIterations_ + 10;
     640    } else {
     641      if (numberInterior > nextUnflag && numberIterations_ > nextUnflagIteration) {
     642        nextUnflagIteration = numberIterations_ + 10;
     643        nextUnflag += 10;
     644        unflag();
    655645#ifdef CLP_DEBUG
    656                     if (handler_->logLevel() & 32)
    657                          printf("unflagging as interior\n");
    658 #endif
    659                }
    660           }
    661           pivotRow_ = -1;
    662           result = pivotColumn(rowArray_[3], rowArray_[0],
    663                                columnArray_[0], rowArray_[1], pivotMode, solutionError,
    664                                array1);
    665           if (result) {
    666                if (result == 2 && sequenceIn_ < 0) {
    667                     // does not look good
    668                     double currentObj;
    669                     double thetaObj;
    670                     double predictedObj;
    671                     objective_->stepLength(this, solution_, solution_, 0.0,
    672                                            currentObj, thetaObj, predictedObj);
    673                     if (currentObj == predictedObj) {
     646        if (handler_->logLevel() & 32)
     647          printf("unflagging as interior\n");
     648#endif
     649      }
     650    }
     651    pivotRow_ = -1;
     652    result = pivotColumn(rowArray_[3], rowArray_[0],
     653      columnArray_[0], rowArray_[1], pivotMode, solutionError,
     654      array1);
     655    if (result) {
     656      if (result == 2 && sequenceIn_ < 0) {
     657        // does not look good
     658        double currentObj;
     659        double thetaObj;
     660        double predictedObj;
     661        objective_->stepLength(this, solution_, solution_, 0.0,
     662          currentObj, thetaObj, predictedObj);
     663        if (currentObj == predictedObj) {
    674664#ifdef CLP_INVESTIGATE
    675                          printf("looks bad - no change in obj %g\n", currentObj);
    676 #endif
    677                          if (factorization_->pivots())
    678                               result = 3;
    679                          else
    680                               problemStatus_ = 0;
    681                     }
    682                }
    683                if (result == 3)
    684                     break; // null vector not accurate
     665          printf("looks bad - no change in obj %g\n", currentObj);
     666#endif
     667          if (factorization_->pivots())
     668            result = 3;
     669          else
     670            problemStatus_ = 0;
     671        }
     672      }
     673      if (result == 3)
     674        break; // null vector not accurate
    685675#ifdef CLP_DEBUG
    686                if (handler_->logLevel() & 32) {
    687                     double currentObj;
    688                     double thetaObj;
    689                     double predictedObj;
    690                     objective_->stepLength(this, solution_, solution_, 0.0,
    691                                            currentObj, thetaObj, predictedObj);
    692                     printf("obj %g after interior move\n", currentObj);
    693                }
    694 #endif
    695                // just move and try again
    696                if (pivotMode < 10) {
    697                     pivotMode = result - 1;
    698                     numberInterior++;
    699                }
    700                continue;
    701           } else {
    702                if (pivotMode < 10) {
    703                     if (theta_ > 0.001)
    704                          pivotMode = 0;
    705                     else if (pivotMode == 2)
    706                          pivotMode = 1;
    707                }
    708                numberInterior = 0;
    709                nextUnflag = 10;
    710                nextUnflagIteration = numberIterations_ + 10;
    711           }
    712           sequenceOut_ = -1;
    713           rowArray_[1]->clear();
    714           if (sequenceIn_ >= 0) {
    715                // we found a pivot column
    716                assert (!flagged(sequenceIn_));
     676      if (handler_->logLevel() & 32) {
     677        double currentObj;
     678        double thetaObj;
     679        double predictedObj;
     680        objective_->stepLength(this, solution_, solution_, 0.0,
     681          currentObj, thetaObj, predictedObj);
     682        printf("obj %g after interior move\n", currentObj);
     683      }
     684#endif
     685      // just move and try again
     686      if (pivotMode < 10) {
     687        pivotMode = result - 1;
     688        numberInterior++;
     689      }
     690      continue;
     691    } else {
     692      if (pivotMode < 10) {
     693        if (theta_ > 0.001)
     694          pivotMode = 0;
     695        else if (pivotMode == 2)
     696          pivotMode = 1;
     697      }
     698      numberInterior = 0;
     699      nextUnflag = 10;
     700      nextUnflagIteration = numberIterations_ + 10;
     701    }
     702    sequenceOut_ = -1;
     703    rowArray_[1]->clear();
     704    if (sequenceIn_ >= 0) {
     705      // we found a pivot column
     706      assert(!flagged(sequenceIn_));
    717707#ifdef CLP_DEBUG
    718                if ((handler_->logLevel() & 32)) {
    719                     char x = isColumn(sequenceIn_) ? 'C' : 'R';
    720                     std::cout << "pivot column " <<
    721                               x << sequenceWithin(sequenceIn_) << std::endl;
    722                }
    723 #endif
    724                // do second half of iteration
    725                if (pivotRow_ < 0 && theta_ < 1.0e-8) {
    726                     assert (sequenceIn_ >= 0);
    727                     returnCode = pivotResult(ifValuesPass);
    728                } else {
    729                     // specialized code
    730                     returnCode = pivotNonlinearResult();
    731                     //printf("odd pivrow %d\n",sequenceOut_);
    732                     if (sequenceOut_ >= 0 && theta_ < 1.0e-5) {
    733                          if (getStatus(sequenceOut_) != isFixed) {
    734                               if (getStatus(sequenceOut_) == atUpperBound)
    735                                    solution_[sequenceOut_] = upper_[sequenceOut_];
    736                               else if (getStatus(sequenceOut_) == atLowerBound)
    737                                    solution_[sequenceOut_] = lower_[sequenceOut_];
    738                               setFlagged(sequenceOut_);
    739                          }
    740                     }
    741                }
    742                if (returnCode < -1 && returnCode > -5) {
    743                     problemStatus_ = -2; //
    744                } else if (returnCode == -5) {
    745                     // something flagged - continue;
    746                } else if (returnCode == 2) {
    747                     problemStatus_ = -5; // looks unbounded
    748                } else if (returnCode == 4) {
    749                     problemStatus_ = -2; // looks unbounded but has iterated
    750                } else if (returnCode != -1) {
    751                     assert(returnCode == 3);
    752                     problemStatus_ = 3;
    753                }
    754           } else {
    755                // no pivot column
     708      if ((handler_->logLevel() & 32)) {
     709        char x = isColumn(sequenceIn_) ? 'C' : 'R';
     710        std::cout << "pivot column " << x << sequenceWithin(sequenceIn_) << std::endl;
     711      }
     712#endif
     713      // do second half of iteration
     714      if (pivotRow_ < 0 && theta_ < 1.0e-8) {
     715        assert(sequenceIn_ >= 0);
     716        returnCode = pivotResult(ifValuesPass);
     717      } else {
     718        // specialized code
     719        returnCode = pivotNonlinearResult();
     720        //printf("odd pivrow %d\n",sequenceOut_);
     721        if (sequenceOut_ >= 0 && theta_ < 1.0e-5) {
     722          if (getStatus(sequenceOut_) != isFixed) {
     723            if (getStatus(sequenceOut_) == atUpperBound)
     724              solution_[sequenceOut_] = upper_[sequenceOut_];
     725            else if (getStatus(sequenceOut_) == atLowerBound)
     726              solution_[sequenceOut_] = lower_[sequenceOut_];
     727            setFlagged(sequenceOut_);
     728          }
     729        }
     730      }
     731      if (returnCode < -1 && returnCode > -5) {
     732        problemStatus_ = -2; //
     733      } else if (returnCode == -5) {
     734        // something flagged - continue;
     735      } else if (returnCode == 2) {
     736        problemStatus_ = -5; // looks unbounded
     737      } else if (returnCode == 4) {
     738        problemStatus_ = -2; // looks unbounded but has iterated
     739      } else if (returnCode != -1) {
     740        assert(returnCode == 3);
     741        problemStatus_ = 3;
     742      }
     743    } else {
     744      // no pivot column
    756745#ifdef CLP_DEBUG
    757                if (handler_->logLevel() & 32)
    758                     printf("** no column pivot\n");
    759 #endif
    760                if (pivotMode < 10) {
    761                     // looks optimal
    762                     primalColumnPivot_->setLooksOptimal(true);
    763                } else {
    764                     pivotMode--;
     746      if (handler_->logLevel() & 32)
     747        printf("** no column pivot\n");
     748#endif
     749      if (pivotMode < 10) {
     750        // looks optimal
     751        primalColumnPivot_->setLooksOptimal(true);
     752      } else {
     753        pivotMode--;
    765754#ifdef CLP_DEBUG
    766                     if (handler_->logLevel() & 32)
    767                          printf("pivot mode now %d\n", pivotMode);
    768 #endif
    769                     if (pivotMode == 9)
    770                          pivotMode = 0; // switch off fast attempt
    771                     unflag();
    772                }
    773                if (nonLinearCost_->numberInfeasibilities())
    774                     problemStatus_ = -4; // might be infeasible
    775                returnCode = 0;
    776                break;
    777           }
    778      }
    779      delete [] array1;
    780      return returnCode;
     755        if (handler_->logLevel() & 32)
     756          printf("pivot mode now %d\n", pivotMode);
     757#endif
     758        if (pivotMode == 9)
     759          pivotMode = 0; // switch off fast attempt
     760        unflag();
     761      }
     762      if (nonLinearCost_->numberInfeasibilities())
     763        problemStatus_ = -4; // might be infeasible
     764      returnCode = 0;
     765      break;
     766    }
     767  }
     768  delete[] array1;
     769  return returnCode;
    781770}
    782771// Creates direction vector
    783 void
    784 ClpSimplexNonlinear::directionVector (CoinIndexedVector * vectorArray,
    785                                       CoinIndexedVector * spare1, CoinIndexedVector * spare2,
    786                                       int pivotMode2,
    787                                       double & normFlagged, double & normUnflagged,
    788                                       int & numberNonBasic)
     772void ClpSimplexNonlinear::directionVector(CoinIndexedVector *vectorArray,
     773  CoinIndexedVector *spare1, CoinIndexedVector *spare2,
     774  int pivotMode2,
     775  double &normFlagged, double &normUnflagged,
     776  int &numberNonBasic)
    789777{
    790778#if CLP_DEBUG > 1
    791      vectorArray->checkClear();
    792      spare1->checkClear();
    793      spare2->checkClear();
    794 #endif
    795      double *array = vectorArray->denseVector();
    796      int * index = vectorArray->getIndices();
    797      int number = 0;
    798      sequenceIn_ = -1;
    799      normFlagged = 0.0;
    800      normUnflagged = 1.0;
    801      double dualTolerance2 = CoinMin(1.0e-8, 1.0e-2 * dualTolerance_);
    802      double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
    803      if (!numberNonBasic) {
    804           //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
    805           //printf("infeasible\n");
    806           if (!pivotMode2 || pivotMode2 >= 10) {
    807                normUnflagged = 0.0;
    808                double bestDj = 0.0;
    809                double bestSuper = 0.0;
    810                double sumSuper = 0.0;
    811                sequenceIn_ = -1;
    812                int nSuper = 0;
    813                for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
    814                     array[iSequence] = 0.0;
    815                     if (flagged(iSequence)) {
    816                          // accumulate  norm
    817                          switch(getStatus(iSequence)) {
    818 
    819                          case basic:
    820                          case ClpSimplex::isFixed:
    821                               break;
    822                          case atUpperBound:
    823                               if (dj_[iSequence] > dualTolerance3)
    824                                    normFlagged += dj_[iSequence] * dj_[iSequence];
    825                               break;
    826                          case atLowerBound:
    827                               if (dj_[iSequence] < -dualTolerance3)
    828                                    normFlagged += dj_[iSequence] * dj_[iSequence];
    829                               break;
    830                          case isFree:
    831                          case superBasic:
    832                               if (fabs(dj_[iSequence]) > dualTolerance3)
    833                                    normFlagged += dj_[iSequence] * dj_[iSequence];
    834                               break;
    835                          }
    836                          continue;
    837                     }
    838                     switch(getStatus(iSequence)) {
    839 
    840                     case basic:
    841                     case ClpSimplex::isFixed:
    842                          break;
    843                     case atUpperBound:
    844                          if (dj_[iSequence] > dualTolerance_) {
    845                               if (dj_[iSequence] > dualTolerance3)
    846                                    normUnflagged += dj_[iSequence] * dj_[iSequence];
    847                               if (pivotMode2 < 10) {
    848                                    array[iSequence] = -dj_[iSequence];
    849                                    index[number++] = iSequence;
    850                               } else {
    851                                    if (dj_[iSequence] > bestDj) {
    852                                         bestDj = dj_[iSequence];
    853                                         sequenceIn_ = iSequence;
    854                                    }
    855                               }
    856                          }
    857                          break;
    858                     case atLowerBound:
    859                          if (dj_[iSequence] < -dualTolerance_) {
    860                               if (dj_[iSequence] < -dualTolerance3)
    861                                    normUnflagged += dj_[iSequence] * dj_[iSequence];
    862                               if (pivotMode2 < 10) {
    863                                    array[iSequence] = -dj_[iSequence];
    864                                    index[number++] = iSequence;
    865                               } else {
    866                                    if (-dj_[iSequence] > bestDj) {
    867                                         bestDj = -dj_[iSequence];
    868                                         sequenceIn_ = iSequence;
    869                                    }
    870                               }
    871                          }
    872                          break;
    873                     case isFree:
    874                     case superBasic:
    875                          //#define ALLSUPER
     779  vectorArray->checkClear();
     780  spare1->checkClear();
     781  spare2->checkClear();
     782#endif
     783  double *array = vectorArray->denseVector();
     784  int *index = vectorArray->getIndices();
     785  int number = 0;
     786  sequenceIn_ = -1;
     787  normFlagged = 0.0;
     788  normUnflagged = 1.0;
     789  double dualTolerance2 = CoinMin(1.0e-8, 1.0e-2 * dualTolerance_);
     790  double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
     791  if (!numberNonBasic) {
     792    //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
     793    //printf("infeasible\n");
     794    if (!pivotMode2 || pivotMode2 >= 10) {
     795      normUnflagged = 0.0;
     796      double bestDj = 0.0;
     797      double bestSuper = 0.0;
     798      double sumSuper = 0.0;
     799      sequenceIn_ = -1;
     800      int nSuper = 0;
     801      for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
     802        array[iSequence] = 0.0;
     803        if (flagged(iSequence)) {
     804          // accumulate  norm
     805          switch (getStatus(iSequence)) {
     806
     807          case basic:
     808          case ClpSimplex::isFixed:
     809            break;
     810          case atUpperBound:
     811            if (dj_[iSequence] > dualTolerance3)
     812              normFlagged += dj_[iSequence] * dj_[iSequence];
     813            break;
     814          case atLowerBound:
     815            if (dj_[iSequence] < -dualTolerance3)
     816              normFlagged += dj_[iSequence] * dj_[iSequence];
     817            break;
     818          case isFree:
     819          case superBasic:
     820            if (fabs(dj_[iSequence]) > dualTolerance3)
     821              normFlagged += dj_[iSequence] * dj_[iSequence];
     822            break;
     823          }
     824          continue;
     825        }
     826        switch (getStatus(iSequence)) {
     827
     828        case basic:
     829        case ClpSimplex::isFixed:
     830          break;
     831        case atUpperBound:
     832          if (dj_[iSequence] > dualTolerance_) {
     833            if (dj_[iSequence] > dualTolerance3)
     834              normUnflagged += dj_[iSequence] * dj_[iSequence];
     835            if (pivotMode2 < 10) {
     836              array[iSequence] = -dj_[iSequence];
     837              index[number++] = iSequence;
     838            } else {
     839              if (dj_[iSequence] > bestDj) {
     840                bestDj = dj_[iSequence];
     841                sequenceIn_ = iSequence;
     842              }
     843            }
     844          }
     845          break;
     846        case atLowerBound:
     847          if (dj_[iSequence] < -dualTolerance_) {
     848            if (dj_[iSequence] < -dualTolerance3)
     849              normUnflagged += dj_[iSequence] * dj_[iSequence];
     850            if (pivotMode2 < 10) {
     851              array[iSequence] = -dj_[iSequence];
     852              index[number++] = iSequence;
     853            } else {
     854              if (-dj_[iSequence] > bestDj) {
     855                bestDj = -dj_[iSequence];
     856                sequenceIn_ = iSequence;
     857              }
     858            }
     859          }
     860          break;
     861        case isFree:
     862        case superBasic:
     863          //#define ALLSUPER
    876864#define NOSUPER
    877865#ifndef ALLSUPER
    878                          if (fabs(dj_[iSequence]) > dualTolerance_) {
    879                               if (fabs(dj_[iSequence]) > dualTolerance3)
    880                                    normUnflagged += dj_[iSequence] * dj_[iSequence];
    881                               nSuper++;
    882                               bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
    883                               sumSuper += fabs(dj_[iSequence]);
    884                          }
    885                          if (fabs(dj_[iSequence]) > dualTolerance2) {
    886                               array[iSequence] = -dj_[iSequence];
    887                               index[number++] = iSequence;
    888                          }
     866          if (fabs(dj_[iSequence]) > dualTolerance_) {
     867            if (fabs(dj_[iSequence]) > dualTolerance3)
     868              normUnflagged += dj_[iSequence] * dj_[iSequence];
     869            nSuper++;
     870            bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
     871            sumSuper += fabs(dj_[iSequence]);
     872          }
     873          if (fabs(dj_[iSequence]) > dualTolerance2) {
     874            array[iSequence] = -dj_[iSequence];
     875            index[number++] = iSequence;
     876          }
    889877#else
    890                          array[iSequence] = -dj_[iSequence];
    891                          index[number++] = iSequence;
    892                          if (pivotMode2 >= 10)
    893                               bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
    894 #endif
    895                          break;
    896                     }
    897                }
     878          array[iSequence] = -dj_[iSequence];
     879          index[number++] = iSequence;
     880          if (pivotMode2 >= 10)
     881            bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
     882#endif
     883          break;
     884        }
     885      }
    898886#ifdef NOSUPER
    899                // redo
    900                bestSuper = sumSuper;
    901                if(sequenceIn_ >= 0 && bestDj > bestSuper) {
    902                     int j;
    903                     // get rid of superbasics
    904                     for (j = 0; j < number; j++) {
    905                          int iSequence = index[j];
    906                          array[iSequence] = 0.0;
    907                     }
    908                     number = 0;
    909                     array[sequenceIn_] = -dj_[sequenceIn_];
    910                     index[number++] = sequenceIn_;
    911                } else {
    912                     sequenceIn_ = -1;
    913                }
     887      // redo
     888      bestSuper = sumSuper;
     889      if (sequenceIn_ >= 0 && bestDj > bestSuper) {
     890        int j;
     891        // get rid of superbasics
     892        for (j = 0; j < number; j++) {
     893          int iSequence = index[j];
     894          array[iSequence] = 0.0;
     895        }
     896        number = 0;
     897        array[sequenceIn_] = -dj_[sequenceIn_];
     898        index[number++] = sequenceIn_;
     899      } else {
     900        sequenceIn_ = -1;
     901      }
    914902#else
    915                if (pivotMode2 >= 10 || !nSuper) {
    916                     bool takeBest = true;
    917                     if (pivotMode2 == 100 && nSuper > 1)
    918                          takeBest = false;
    919                     if(sequenceIn_ >= 0 && takeBest) {
    920                          if (fabs(dj_[sequenceIn_]) > bestSuper) {
    921                               array[sequenceIn_] = -dj_[sequenceIn_];
    922                               index[number++] = sequenceIn_;
    923                          } else {
    924                               sequenceIn_ = -1;
    925                          }
    926                     } else {
    927                          sequenceIn_ = -1;
    928                     }
    929                }
     903      if (pivotMode2 >= 10 || !nSuper) {
     904        bool takeBest = true;
     905        if (pivotMode2 == 100 && nSuper > 1)
     906          takeBest = false;
     907        if (sequenceIn_ >= 0 && takeBest) {
     908          if (fabs(dj_[sequenceIn_]) > bestSuper) {
     909            array[sequenceIn_] = -dj_[sequenceIn_];
     910            index[number++] = sequenceIn_;
     911          } else {
     912            sequenceIn_ = -1;
     913          }
     914        } else {
     915          sequenceIn_ = -1;
     916        }
     917      }
    930918#endif
    931919#ifdef CLP_DEBUG
    932                if (handler_->logLevel() & 32) {
    933                     if (sequenceIn_ >= 0)
    934                          printf("%d superBasic, chosen %d - dj %g\n", nSuper, sequenceIn_,
    935                                 dj_[sequenceIn_]);
    936                     else
    937                          printf("%d superBasic - none chosen\n", nSuper);
    938                }
    939 #endif
    940           } else {
    941                double bestDj = 0.0;
    942                double saveDj = 0.0;
    943                if (sequenceOut_ >= 0) {
    944                     saveDj = dj_[sequenceOut_];
    945                     dj_[sequenceOut_] = 0.0;
    946                     switch(getStatus(sequenceOut_)) {
    947 
    948                     case basic:
    949                          sequenceOut_ = -1;
    950                     case ClpSimplex::isFixed:
    951                          break;
    952                     case atUpperBound:
    953                          if (dj_[sequenceOut_] > dualTolerance_) {
     920      if (handler_->logLevel() & 32) {
     921        if (sequenceIn_ >= 0)
     922          printf("%d superBasic, chosen %d - dj %g\n", nSuper, sequenceIn_,
     923            dj_[sequenceIn_]);
     924        else
     925          printf("%d superBasic - none chosen\n", nSuper);
     926      }
     927#endif
     928    } else {
     929      double bestDj = 0.0;
     930      double saveDj = 0.0;
     931      if (sequenceOut_ >= 0) {
     932        saveDj = dj_[sequenceOut_];
     933        dj_[sequenceOut_] = 0.0;
     934        switch (getStatus(sequenceOut_)) {
     935
     936        case basic:
     937          sequenceOut_ = -1;
     938        case ClpSimplex::isFixed:
     939          break;
     940        case atUpperBound:
     941          if (dj_[sequenceOut_] > dualTolerance_) {
    954942#ifdef CLP_DEBUG
    955                               if (handler_->logLevel() & 32)
    956                                    printf("after pivot out %d values %g %g %g, dj %g\n",
    957                                           sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
    958                                           upper_[sequenceOut_], dj_[sequenceOut_]);
    959 #endif
    960                          }
    961                          break;
    962                     case atLowerBound:
    963                          if (dj_[sequenceOut_] < -dualTolerance_) {
     943            if (handler_->logLevel() & 32)
     944              printf("after pivot out %d values %g %g %g, dj %g\n",
     945                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
     946                upper_[sequenceOut_], dj_[sequenceOut_]);
     947#endif
     948          }
     949          break;
     950        case atLowerBound:
     951          if (dj_[sequenceOut_] < -dualTolerance_) {
    964952#ifdef CLP_DEBUG
    965                               if (handler_->logLevel() & 32)
    966                                    printf("after pivot out %d values %g %g %g, dj %g\n",
    967                                           sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
    968                                           upper_[sequenceOut_], dj_[sequenceOut_]);
    969 #endif
    970                          }
    971                          break;
    972                     case isFree:
    973                     case superBasic:
    974                          if (dj_[sequenceOut_] > dualTolerance_) {
     953            if (handler_->logLevel() & 32)
     954              printf("after pivot out %d values %g %g %g, dj %g\n",
     955                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
     956                upper_[sequenceOut_], dj_[sequenceOut_]);
     957#endif
     958          }
     959          break;
     960        case isFree:
     961        case superBasic:
     962          if (dj_[sequenceOut_] > dualTolerance_) {
    975963#ifdef CLP_DEBUG
    976                               if (handler_->logLevel() & 32)
    977                                    printf("after pivot out %d values %g %g %g, dj %g\n",
    978                                           sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
    979                                           upper_[sequenceOut_], dj_[sequenceOut_]);
    980 #endif
    981                          } else if (dj_[sequenceOut_] < -dualTolerance_) {
     964            if (handler_->logLevel() & 32)
     965              printf("after pivot out %d values %g %g %g, dj %g\n",
     966                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
     967                upper_[sequenceOut_], dj_[sequenceOut_]);
     968#endif
     969          } else if (dj_[sequenceOut_] < -dualTolerance_) {
    982970#ifdef CLP_DEBUG
    983                               if (handler_->logLevel() & 32)
    984                                    printf("after pivot out %d values %g %g %g, dj %g\n",
    985                                           sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
    986                                           upper_[sequenceOut_], dj_[sequenceOut_]);
    987 #endif
    988                          }
    989                          break;
    990                     }
    991                }
    992                // Go for dj
    993                pivotMode2 = 3;
    994                for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
    995                     array[iSequence] = 0.0;
    996                     if (flagged(iSequence))
    997                          continue;
    998                     switch(getStatus(iSequence)) {
    999 
    1000                     case basic:
    1001                     case ClpSimplex::isFixed:
    1002                          break;
    1003                     case atUpperBound:
    1004                          if (dj_[iSequence] > dualTolerance_) {
    1005                               double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]);
    1006                               double merit = distance * dj_[iSequence];
    1007                               if (pivotMode2 == 1)
    1008                                    merit *= 1.0e-20; // discourage
    1009                               if (pivotMode2 == 3)
    1010                                    merit = fabs(dj_[iSequence]);
    1011                               if (merit > bestDj) {
    1012                                    sequenceIn_ = iSequence;
    1013                                    bestDj = merit;
    1014                               }
    1015                          }
    1016                          break;
    1017                     case atLowerBound:
    1018                          if (dj_[iSequence] < -dualTolerance_) {
    1019                               double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]);
    1020                               double merit = -distance * dj_[iSequence];
    1021                               if (pivotMode2 == 1)
    1022                                    merit *= 1.0e-20; // discourage
    1023                               if (pivotMode2 == 3)
    1024                                    merit = fabs(dj_[iSequence]);
    1025                               if (merit > bestDj) {
    1026                                    sequenceIn_ = iSequence;
    1027                                    bestDj = merit;
    1028                               }
    1029                          }
    1030                          break;
    1031                     case isFree:
    1032                     case superBasic:
    1033                          if (dj_[iSequence] > dualTolerance_) {
    1034                               double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]);
    1035                               distance = CoinMin(solution_[iSequence] - lower_[iSequence],
    1036                                                  upper_[iSequence] - solution_[iSequence]);
    1037                               double merit = distance * dj_[iSequence];
    1038                               if (pivotMode2 == 1)
    1039                                    merit = distance;
    1040                               if (pivotMode2 == 3)
    1041                                    merit = fabs(dj_[iSequence]);
    1042                               if (merit > bestDj) {
    1043                                    sequenceIn_ = iSequence;
    1044                                    bestDj = merit;
    1045                               }
    1046                          } else if (dj_[iSequence] < -dualTolerance_) {
    1047                               double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]);
    1048                               distance = CoinMin(solution_[iSequence] - lower_[iSequence],
    1049                                                  upper_[iSequence] - solution_[iSequence]);
    1050                               double merit = -distance * dj_[iSequence];
    1051                               if (pivotMode2 == 1)
    1052                                    merit = distance;
    1053                               if (pivotMode2 == 3)
    1054                                    merit = fabs(dj_[iSequence]);
    1055                               if (merit > bestDj) {
    1056                                    sequenceIn_ = iSequence;
    1057                                    bestDj = merit;
    1058                               }
    1059                          }
    1060                          break;
    1061                     }
    1062                }
    1063                if (sequenceOut_ >= 0) {
    1064                     dj_[sequenceOut_] = saveDj;
    1065                     sequenceOut_ = -1;
    1066                }
    1067                if (sequenceIn_ >= 0) {
    1068                     array[sequenceIn_] = -dj_[sequenceIn_];
    1069                     index[number++] = sequenceIn_;
    1070                }
    1071           }
    1072           numberNonBasic = number;
    1073      } else {
    1074           // compute norms
    1075           normUnflagged = 0.0;
    1076           for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
    1077                if (flagged(iSequence)) {
    1078                     // accumulate  norm
    1079                     switch(getStatus(iSequence)) {
    1080 
    1081                     case basic:
    1082                     case ClpSimplex::isFixed:
    1083                          break;
    1084                     case atUpperBound:
    1085                          if (dj_[iSequence] > dualTolerance_)
    1086                               normFlagged += dj_[iSequence] * dj_[iSequence];
    1087                          break;
    1088                     case atLowerBound:
    1089                          if (dj_[iSequence] < -dualTolerance_)
    1090                               normFlagged += dj_[iSequence] * dj_[iSequence];
    1091                          break;
    1092                     case isFree:
    1093                     case superBasic:
    1094                          if (fabs(dj_[iSequence]) > dualTolerance_)
    1095                               normFlagged += dj_[iSequence] * dj_[iSequence];
    1096                          break;
    1097                     }
    1098                }
    1099           }
    1100           // re-use list
    1101           number = 0;
    1102           int j;
    1103           for (j = 0; j < numberNonBasic; j++) {
    1104                int iSequence = index[j];
    1105                if (flagged(iSequence))
    1106                     continue;
    1107                switch(getStatus(iSequence)) {
    1108 
    1109                case basic:
    1110                case ClpSimplex::isFixed:
    1111                     continue; //abort();
    1112                     break;
    1113                case atUpperBound:
    1114                     if (dj_[iSequence] > dualTolerance_) {
    1115                          number++;
    1116                          normUnflagged += dj_[iSequence] * dj_[iSequence];
    1117                     }
    1118                     break;
    1119                case atLowerBound:
    1120                     if (dj_[iSequence] < -dualTolerance_) {
    1121                          number++;
    1122                          normUnflagged += dj_[iSequence] * dj_[iSequence];
    1123                     }
    1124                     break;
    1125                case isFree:
    1126                case superBasic:
    1127                     if (fabs(dj_[iSequence]) > dualTolerance_) {
    1128                          number++;
    1129                          normUnflagged += dj_[iSequence] * dj_[iSequence];
    1130                     }
    1131                     break;
    1132                }
    1133                array[iSequence] = -dj_[iSequence];
    1134           }
    1135           // switch to large
    1136           normUnflagged = 1.0;
    1137           if (!number) {
    1138                for (j = 0; j < numberNonBasic; j++) {
    1139                     int iSequence = index[j];
    1140                     array[iSequence] = 0.0;
    1141                }
    1142                numberNonBasic = 0;
    1143           }
    1144           number = numberNonBasic;
    1145      }
    1146      if (number) {
    1147           int j;
    1148           // Now do basic ones - how do I compensate for small basic infeasibilities?
    1149           int iRow;
    1150           for (iRow = 0; iRow < numberRows_; iRow++) {
    1151                int iPivot = pivotVariable_[iRow];
    1152                double value = 0.0;
    1153                if (solution_[iPivot] > upper_[iPivot]) {
    1154                     value = upper_[iPivot] - solution_[iPivot];
    1155                } else if (solution_[iPivot] < lower_[iPivot]) {
    1156                     value = lower_[iPivot] - solution_[iPivot];
    1157                }
    1158                //if (value)
    1159                //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
    1160                //     upper_[iPivot]);
    1161                //value=0.0;
    1162                value *= -1.0e0;
    1163                if (value) {
    1164                     array[iPivot] = value;
    1165                     index[number++] = iPivot;
    1166                }
    1167           }
    1168           double * array2 = spare1->denseVector();
    1169           int * index2 = spare1->getIndices();
    1170           int number2 = 0;
    1171           times(-1.0, array, array2);
    1172           array = array + numberColumns_;
    1173           for (iRow = 0; iRow < numberRows_; iRow++) {
    1174                double value = array2[iRow] + array[iRow];
    1175                if (value) {
    1176                     array2[iRow] = value;
    1177                     index2[number2++] = iRow;
    1178                } else {
    1179                     array2[iRow] = 0.0;
    1180                }
    1181           }
    1182           array -= numberColumns_;
    1183           spare1->setNumElements(number2);
    1184           // Ftran
    1185           factorization_->updateColumn(spare2, spare1);
    1186           number2 = spare1->getNumElements();
    1187           for (j = 0; j < number2; j++) {
    1188                int iSequence = index2[j];
    1189                double value = array2[iSequence];
    1190                array2[iSequence] = 0.0;
    1191                if (value) {
    1192                     int iPivot = pivotVariable_[iSequence];
    1193                     double oldValue = array[iPivot];
    1194                     if (!oldValue) {
    1195                          array[iPivot] = value;
    1196                          index[number++] = iPivot;
    1197                     } else {
    1198                          // something already there
    1199                          array[iPivot] = value + oldValue;
    1200                     }
    1201                }
    1202           }
    1203           spare1->setNumElements(0);
    1204      }
    1205      vectorArray->setNumElements(number);
     971            if (handler_->logLevel() & 32)
     972              printf("after pivot out %d values %g %g %g, dj %g\n",
     973                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
     974                upper_[sequenceOut_], dj_[sequenceOut_]);
     975#endif
     976          }
     977          break;
     978        }
     979      }
     980      // Go for dj
     981      pivotMode2 = 3;
     982      for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
     983        array[iSequence] = 0.0;
     984        if (flagged(iSequence))
     985          continue;
     986        switch (getStatus(iSequence)) {
     987
     988        case basic:
     989        case ClpSimplex::isFixed:
     990          break;
     991        case atUpperBound:
     992          if (dj_[iSequence] > dualTolerance_) {
     993            double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]);
     994            double merit = distance * dj_[iSequence];
     995            if (pivotMode2 == 1)
     996              merit *= 1.0e-20; // discourage
     997            if (pivotMode2 == 3)
     998              merit = fabs(dj_[iSequence]);
     999            if (merit > bestDj) {
     1000              sequenceIn_ = iSequence;
     1001              bestDj = merit;
     1002            }
     1003          }
     1004          break;
     1005        case atLowerBound:
     1006          if (dj_[iSequence] < -dualTolerance_) {
     1007            double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]);
     1008            double merit = -distance * dj_[iSequence];
     1009            if (pivotMode2 == 1)
     1010              merit *= 1.0e-20; // discourage
     1011            if (pivotMode2 == 3)
     1012              merit = fabs(dj_[iSequence]);
     1013            if (merit > bestDj) {
     1014              sequenceIn_ = iSequence;
     1015              bestDj = merit;
     1016            }
     1017          }
     1018          break;
     1019        case isFree:
     1020        case superBasic:
     1021          if (dj_[iSequence] > dualTolerance_) {
     1022            double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]);
     1023            distance = CoinMin(solution_[iSequence] - lower_[iSequence],
     1024              upper_[iSequence] - solution_[iSequence]);
     1025            double merit = distance * dj_[iSequence];
     1026            if (pivotMode2 == 1)
     1027              merit = distance;
     1028            if (pivotMode2 == 3)
     1029              merit = fabs(dj_[iSequence]);
     1030            if (merit > bestDj) {
     1031              sequenceIn_ = iSequence;
     1032              bestDj = merit;
     1033            }
     1034          } else if (dj_[iSequence] < -dualTolerance_) {
     1035            double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]);
     1036            distance = CoinMin(solution_[iSequence] - lower_[iSequence],
     1037              upper_[iSequence] - solution_[iSequence]);
     1038            double merit = -distance * dj_[iSequence];
     1039            if (pivotMode2 == 1)
     1040              merit = distance;
     1041            if (pivotMode2 == 3)
     1042              merit = fabs(dj_[iSequence]);
     1043            if (merit > bestDj) {
     1044              sequenceIn_ = iSequence;
     1045              bestDj = merit;
     1046            }
     1047          }
     1048          break;
     1049        }
     1050      }
     1051      if (sequenceOut_ >= 0) {
     1052        dj_[sequenceOut_] = saveDj;
     1053        sequenceOut_ = -1;
     1054      }
     1055      if (sequenceIn_ >= 0) {
     1056        array[sequenceIn_] = -dj_[sequenceIn_];
     1057        index[number++] = sequenceIn_;
     1058      }
     1059    }
     1060    numberNonBasic = number;
     1061  } else {
     1062    // compute norms
     1063    normUnflagged = 0.0;
     1064    for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
     1065      if (flagged(iSequence)) {
     1066        // accumulate  norm
     1067        switch (getStatus(iSequence)) {
     1068
     1069        case basic:
     1070        case ClpSimplex::isFixed:
     1071          break;
     1072        case atUpperBound:
     1073          if (dj_[iSequence] > dualTolerance_)
     1074            normFlagged += dj_[iSequence] * dj_[iSequence];
     1075          break;
     1076        case atLowerBound:
     1077          if (dj_[iSequence] < -dualTolerance_)
     1078            normFlagged += dj_[iSequence] * dj_[iSequence];
     1079          break;
     1080        case isFree:
     1081        case superBasic:
     1082          if (fabs(dj_[iSequence]) > dualTolerance_)
     1083            normFlagged += dj_[iSequence] * dj_[iSequence];
     1084          break;
     1085        }
     1086      }
     1087    }
     1088    // re-use list
     1089    number = 0;
     1090    int j;
     1091    for (j = 0; j < numberNonBasic; j++) {
     1092      int iSequence = index[j];
     1093      if (flagged(iSequence))
     1094        continue;
     1095      switch (getStatus(iSequence)) {
     1096
     1097      case basic:
     1098      case ClpSimplex::isFixed:
     1099        continue; //abort();
     1100        break;
     1101      case atUpperBound:
     1102        if (dj_[iSequence] > dualTolerance_) {
     1103          number++;
     1104          normUnflagged += dj_[iSequence] * dj_[iSequence];
     1105        }
     1106        break;
     1107      case atLowerBound:
     1108        if (dj_[iSequence] < -dualTolerance_) {
     1109          number++;
     1110          normUnflagged += dj_[iSequence] * dj_[iSequence];
     1111        }
     1112        break;
     1113      case isFree:
     1114      case superBasic:
     1115        if (fabs(dj_[iSequence]) > dualTolerance_) {
     1116          number++;
     1117          normUnflagged += dj_[iSequence] * dj_[iSequence];
     1118        }
     1119        break;
     1120      }
     1121      array[iSequence] = -dj_[iSequence];
     1122    }
     1123    // switch to large
     1124    normUnflagged = 1.0;
     1125    if (!number) {
     1126      for (j = 0; j < numberNonBasic; j++) {
     1127        int iSequence = index[j];
     1128        array[iSequence] = 0.0;
     1129      }
     1130      numberNonBasic = 0;
     1131    }
     1132    number = numberNonBasic;
     1133  }
     1134  if (number) {
     1135    int j;
     1136    // Now do basic ones - how do I compensate for small basic infeasibilities?
     1137    int iRow;
     1138    for (iRow = 0; iRow < numberRows_; iRow++) {
     1139      int iPivot = pivotVariable_[iRow];
     1140      double value = 0.0;
     1141      if (solution_[iPivot] > upper_[iPivot]) {
     1142        value = upper_[iPivot] - solution_[iPivot];
     1143      } else if (solution_[iPivot] < lower_[iPivot]) {
     1144        value = lower_[iPivot] - solution_[iPivot];
     1145      }
     1146      //if (value)
     1147      //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
     1148      //     upper_[iPivot]);
     1149      //value=0.0;
     1150      value *= -1.0e0;
     1151      if (value) {
     1152        array[iPivot] = value;
     1153        index[number++] = iPivot;
     1154      }
     1155    }
     1156    double *array2 = spare1->denseVector();
     1157    int *index2 = spare1->getIndices();
     1158    int number2 = 0;
     1159    times(-1.0, array, array2);
     1160    array = array + numberColumns_;
     1161    for (iRow = 0; iRow < numberRows_; iRow++) {
     1162      double value = array2[iRow] + array[iRow];
     1163      if (value) {
     1164        array2[iRow] = value;
     1165        index2[number2++] = iRow;
     1166      } else {
     1167        array2[iRow] = 0.0;
     1168      }
     1169    }
     1170    array -= numberColumns_;
     1171    spare1->setNumElements(number2);
     1172    // Ftran
     1173    factorization_->updateColumn(spare2, spare1);
     1174    number2 = spare1->getNumElements();
     1175    for (j = 0; j < number2; j++) {
     1176      int iSequence = index2[j];
     1177      double value = array2[iSequence];
     1178      array2[iSequence] = 0.0;
     1179      if (value) {
     1180        int iPivot = pivotVariable_[iSequence];
     1181        double oldValue = array[iPivot];
     1182        if (!oldValue) {
     1183          array[iPivot] = value;
     1184          index[number++] = iPivot;
     1185        } else {
     1186          // something already there
     1187          array[iPivot] = value + oldValue;
     1188        }
     1189      }
     1190    }
     1191    spare1->setNumElements(0);
     1192  }
     1193  vectorArray->setNumElements(number);
    12061194}
    12071195#define MINTYPE 1
    1208 #if MINTYPE==2
     1196#if MINTYPE == 2
    12091197static double
    1210 innerProductIndexed(const double * region1, int size, const double * region2, const int * which)
     1198innerProductIndexed(const double *region1, int size, const double *region2, const int *which)
    12111199{
    1212      int i;
    1213      double value = 0.0;
    1214      for (i = 0; i < size; i++) {
    1215           int j = which[i];
    1216           value += region1[j] * region2[j];
    1217      }
    1218      return value;
     1200  int i;
     1201  double value = 0.0;
     1202  for (i = 0; i < size; i++) {
     1203    int j = which[i];
     1204    value += region1[j] * region2[j];
     1205  }
     1206  return value;
    12191207}
    12201208#endif
     
    12241212   1 - no basis change
    12251213*/
    1226 int
    1227 ClpSimplexNonlinear::pivotColumn(CoinIndexedVector * longArray,
    1228                                  CoinIndexedVector * rowArray,
    1229                                  CoinIndexedVector * columnArray,
    1230                                  CoinIndexedVector * spare,
    1231                                  int & pivotMode,
    1232                                  double & solutionError,
    1233                                  double * dArray)
     1214int ClpSimplexNonlinear::pivotColumn(CoinIndexedVector *longArray,
     1215  CoinIndexedVector *rowArray,
     1216  CoinIndexedVector *columnArray,
     1217  CoinIndexedVector *spare,
     1218  int &pivotMode,
     1219  double &solutionError,
     1220  double *dArray)
    12341221{
    1235      // say not optimal
    1236      primalColumnPivot_->setLooksOptimal(false);
    1237      double acceptablePivot = 1.0e-10;
    1238      int lastSequenceIn = -1;
    1239      if (pivotMode && pivotMode < 10) {
    1240           acceptablePivot = 1.0e-6;
    1241           if (factorization_->pivots())
    1242                acceptablePivot = 1.0e-5; // if we have iterated be more strict
    1243      }
    1244      double acceptableBasic = 1.0e-7;
    1245 
    1246      int number = longArray->getNumElements();
    1247      int numberTotal = numberRows_ + numberColumns_;
    1248      int bestSequence = -1;
    1249      int bestBasicSequence = -1;
    1250      double eps = 1.0e-1;
    1251      eps = 1.0e-6;
    1252      double basicTheta = 1.0e30;
    1253      double objTheta = 0.0;
    1254      bool finished = false;
    1255      sequenceIn_ = -1;
    1256      int nPasses = 0;
    1257      int nTotalPasses = 0;
    1258      int nBigPasses = 0;
    1259      double djNorm0 = 0.0;
    1260      double djNorm = 0.0;
    1261      double normFlagged = 0.0;
    1262      double normUnflagged = 0.0;
    1263      int localPivotMode = pivotMode;
    1264      bool allFinished = false;
    1265      bool justOne = false;
    1266      int returnCode = 1;
    1267      double currentObj;
    1268      double predictedObj;
    1269      double thetaObj;
    1270      objective_->stepLength(this, solution_, solution_, 0.0,
    1271                             currentObj, predictedObj, thetaObj);
    1272      double saveObj = currentObj;
    1273 #if MINTYPE ==2
    1274      // try Shanno's method
    1275      //would be memory leak
    1276      //double * saveY=new double[numberTotal];
    1277      //double * saveS=new double[numberTotal];
    1278      //double * saveY2=new double[numberTotal];
    1279      //double * saveS2=new double[numberTotal];
    1280      double saveY[100];
    1281      double saveS[100];
    1282      double saveY2[100];
    1283      double saveS2[100];
    1284      double zz[10000];
    1285 #endif
    1286      double * dArray2 = dArray + numberTotal;
    1287      // big big loop
    1288      while (!allFinished) {
    1289           double * work = longArray->denseVector();
    1290           int * which = longArray->getIndices();
    1291           allFinished = true;
    1292           // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
    1293           //#define SMALLTHETA1 1.0e-25
    1294           //#define SMALLTHETA2 1.0e-25
     1222  // say not optimal
     1223  primalColumnPivot_->setLooksOptimal(false);
     1224  double acceptablePivot = 1.0e-10;
     1225  int lastSequenceIn = -1;
     1226  if (pivotMode && pivotMode < 10) {
     1227    acceptablePivot = 1.0e-6;
     1228    if (factorization_->pivots())
     1229      acceptablePivot = 1.0e-5; // if we have iterated be more strict
     1230  }
     1231  double acceptableBasic = 1.0e-7;
     1232
     1233  int number = longArray->getNumElements();
     1234  int numberTotal = numberRows_ + numberColumns_;
     1235  int bestSequence = -1;
     1236  int bestBasicSequence = -1;
     1237  double eps = 1.0e-1;
     1238  eps = 1.0e-6;
     1239  double basicTheta = 1.0e30;
     1240  double objTheta = 0.0;
     1241  bool finished = false;
     1242  sequenceIn_ = -1;
     1243  int nPasses = 0;
     1244  int nTotalPasses = 0;
     1245  int nBigPasses = 0;
     1246  double djNorm0 = 0.0;
     1247  double djNorm = 0.0;
     1248  double normFlagged = 0.0;
     1249  double normUnflagged = 0.0;
     1250  int localPivotMode = pivotMode;
     1251  bool allFinished = false;
     1252  bool justOne = false;
     1253  int returnCode = 1;
     1254  double currentObj;
     1255  double predictedObj;
     1256  double thetaObj;
     1257  objective_->stepLength(this, solution_, solution_, 0.0,
     1258    currentObj, predictedObj, thetaObj);
     1259  double saveObj = currentObj;
     1260#if MINTYPE == 2
     1261  // try Shanno's method
     1262  //would be memory leak
     1263  //double * saveY=new double[numberTotal];
     1264  //double * saveS=new double[numberTotal];
     1265  //double * saveY2=new double[numberTotal];
     1266  //double * saveS2=new double[numberTotal];
     1267  double saveY[100];
     1268  double saveS[100];
     1269  double saveY2[100];
     1270  double saveS2[100];
     1271  double zz[10000];
     1272#endif
     1273  double *dArray2 = dArray + numberTotal;
     1274  // big big loop
     1275  while (!allFinished) {
     1276    double *work = longArray->denseVector();
     1277    int *which = longArray->getIndices();
     1278    allFinished = true;
     1279    // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
     1280    //#define SMALLTHETA1 1.0e-25
     1281    //#define SMALLTHETA2 1.0e-25
    12951282#define SMALLTHETA1 1.0e-10
    12961283#define SMALLTHETA2 1.0e-10
    12971284#define CONJUGATE 2
    12981285#if CONJUGATE == 0
    1299           int conjugate = 0;
     1286    int conjugate = 0;
    13001287#elif CONJUGATE == 1
    1301           int conjugate = (pivotMode < 10) ? MINTYPE : 0;
     1288    int conjugate = (pivotMode < 10) ? MINTYPE : 0;
    13021289#elif CONJUGATE == 2
    1303           int conjugate = MINTYPE;
     1290    int conjugate = MINTYPE;
    13041291#else
    1305           int conjugate = MINTYPE;
    1306 #endif
    1307 
    1308           //if (pivotMode==1)
    1309           //localPivotMode=11;;
     1292    int conjugate = MINTYPE;
     1293#endif
     1294
     1295    //if (pivotMode==1)
     1296    //localPivotMode=11;;
    13101297#if CLP_DEBUG > 1
    1311           longArray->checkClear();
    1312 #endif
    1313           int numberNonBasic = 0;
    1314           int startLocalMode = -1;
    1315           while (!finished) {
    1316                double simpleObjective = COIN_DBL_MAX;
    1317                returnCode = 1;
    1318                int iSequence;
    1319                objective_->reducedGradient(this, dj_, false);
    1320                // get direction vector in longArray
    1321                longArray->clear();
    1322                // take out comment to try slightly different idea
    1323                if (nPasses + nTotalPasses > 3000 || nBigPasses > 100) {
    1324                     if (factorization_->pivots())
    1325                          returnCode = 3;
    1326                     break;
    1327                }
    1328                if (!nPasses) {
    1329                     numberNonBasic = 0;
    1330                     nBigPasses++;
    1331                }
    1332                // just do superbasic if in cleanup mode
    1333                int local = localPivotMode;
    1334                if (!local && pivotMode >= 10 && nBigPasses < 10) {
    1335                     local = 100;
    1336                } else if (local == -1 || nBigPasses >= 10) {
    1337                     local = 0;
    1338                     localPivotMode = 0;
    1339                }
    1340                if (justOne) {
    1341                     local = 2;
    1342                     //local=100;
    1343                     justOne = false;
    1344                }
    1345                if (!nPasses)
    1346                     startLocalMode = local;
    1347                directionVector(longArray, spare, rowArray, local,
    1348                                normFlagged, normUnflagged, numberNonBasic);
    1349                {
    1350                     // check null vector
    1351                     double * rhs = spare->denseVector();
     1298    longArray->checkClear();
     1299#endif
     1300    int numberNonBasic = 0;
     1301    int startLocalMode = -1;
     1302    while (!finished) {
     1303      double simpleObjective = COIN_DBL_MAX;
     1304      returnCode = 1;
     1305      int iSequence;
     1306      objective_->reducedGradient(this, dj_, false);
     1307      // get direction vector in longArray
     1308      longArray->clear();
     1309      // take out comment to try slightly different idea
     1310      if (nPasses + nTotalPasses > 3000 || nBigPasses > 100) {
     1311        if (factorization_->pivots())
     1312          returnCode = 3;
     1313        break;
     1314      }
     1315      if (!nPasses) {
     1316        numberNonBasic = 0;
     1317        nBigPasses++;
     1318      }
     1319      // just do superbasic if in cleanup mode
     1320      int local = localPivotMode;
     1321      if (!local && pivotMode >= 10 && nBigPasses < 10) {
     1322        local = 100;
     1323      } else if (local == -1 || nBigPasses >= 10) {
     1324        local = 0;
     1325        localPivotMode = 0;
     1326      }
     1327      if (justOne) {
     1328        local = 2;
     1329        //local=100;
     1330        justOne = false;
     1331      }
     1332      if (!nPasses)
     1333        startLocalMode = local;
     1334      directionVector(longArray, spare, rowArray, local,
     1335        normFlagged, normUnflagged, numberNonBasic);
     1336      {
     1337        // check null vector
     1338        double *rhs = spare->denseVector();
    13521339#if CLP_DEBUG > 1
    1353                     spare->checkClear();
    1354 #endif
    1355                     int iRow;
    1356                     multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, rhs, 0.0);
    1357                     matrix_->times(1.0, solution_, rhs, rowScale_, columnScale_);
    1358                     double largest = 0.0;
     1340        spare->checkClear();
     1341#endif
     1342        int iRow;
     1343        multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, rhs, 0.0);
     1344        matrix_->times(1.0, solution_, rhs, rowScale_, columnScale_);
     1345        double largest = 0.0;
    13591346#if CLP_DEBUG > 0
    1360                     int iLargest = -1;
    1361 #endif
    1362                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1363                          double value = fabs(rhs[iRow]);
    1364                          rhs[iRow] = 0.0;
    1365                          if (value > largest) {
    1366                               largest = value;
     1347        int iLargest = -1;
     1348#endif
     1349        for (iRow = 0; iRow < numberRows_; iRow++) {
     1350          double value = fabs(rhs[iRow]);
     1351          rhs[iRow] = 0.0;
     1352          if (value > largest) {
     1353            largest = value;
    13671354#if CLP_DEBUG > 0
    1368                               iLargest = iRow;
    1369 #endif
    1370                          }
    1371                     }
     1355            iLargest = iRow;
     1356#endif
     1357          }
     1358        }
    13721359#if CLP_DEBUG > 0
    1373                     if ((handler_->logLevel() & 32) && largest > 1.0e-8)
    1374                          printf("largest rhs error %g on row %d\n", largest, iLargest);
    1375 #endif
    1376                     if (solutionError < 0.0) {
    1377                          solutionError = largest;
    1378                     } else if (largest > CoinMax(1.0e-8, 1.0e2 * solutionError) &&
    1379                                factorization_->pivots()) {
    1380                          longArray->clear();
    1381                          pivotRow_ = -1;
    1382                          theta_ = 0.0;
    1383                          return 3;
    1384                     }
    1385                }
    1386                if (sequenceIn_ >= 0)
    1387                     lastSequenceIn = sequenceIn_;
    1388 #if MINTYPE!=2
    1389                double djNormSave = djNorm;
    1390 #endif
    1391                djNorm = 0.0;
    1392                int iIndex;
    1393                for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
    1394                     int iSequence = which[iIndex];
    1395                     double alpha = work[iSequence];
    1396                     djNorm += alpha * alpha;
    1397                }
    1398                // go to conjugate gradient if necessary
    1399                if (numberNonBasic && localPivotMode >= 10 && (nPasses > 4 || sequenceIn_ < 0)) {
    1400                     localPivotMode = 0;
    1401                     nTotalPasses += nPasses;
    1402                     nPasses = 0;
    1403                }
     1360        if ((handler_->logLevel() & 32) && largest > 1.0e-8)
     1361          printf("largest rhs error %g on row %d\n", largest, iLargest);
     1362#endif
     1363        if (solutionError < 0.0) {
     1364          solutionError = largest;
     1365        } else if (largest > CoinMax(1.0e-8, 1.0e2 * solutionError) && factorization_->pivots()) {
     1366          longArray->clear();
     1367          pivotRow_ = -1;
     1368          theta_ = 0.0;
     1369          return 3;
     1370        }
     1371      }
     1372      if (sequenceIn_ >= 0)
     1373        lastSequenceIn = sequenceIn_;
     1374#if MINTYPE != 2
     1375      double djNormSave = djNorm;
     1376#endif
     1377      djNorm = 0.0;
     1378      int iIndex;
     1379      for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
     1380        int iSequence = which[iIndex];
     1381        double alpha = work[iSequence];
     1382        djNorm += alpha * alpha;
     1383      }
     1384      // go to conjugate gradient if necessary
     1385      if (numberNonBasic && localPivotMode >= 10 && (nPasses > 4 || sequenceIn_ < 0)) {
     1386        localPivotMode = 0;
     1387        nTotalPasses += nPasses;
     1388        nPasses = 0;
     1389      }
    14041390#if CONJUGATE == 2
    1405                conjugate = (localPivotMode < 10) ? MINTYPE : 0;
    1406 #endif
    1407                //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
    1408                //     nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
    1409                if (!nPasses) {
    1410                     //printf("numberNon %d\n",numberNonBasic);
    1411 #if MINTYPE==2
    1412                     assert (numberNonBasic < 100);
    1413                     memset(zz, 0, numberNonBasic * numberNonBasic * sizeof(double));
    1414                     int put = 0;
    1415                     for (int iVariable = 0; iVariable < numberNonBasic; iVariable++) {
    1416                          zz[put] = 1.0;
    1417                          put += numberNonBasic + 1;
    1418                     }
    1419 #endif
    1420                     djNorm0 = CoinMax(djNorm, 1.0e-20);
    1421                     CoinMemcpyN(work, numberTotal, dArray);
    1422                     CoinMemcpyN(work, numberTotal, dArray2);
    1423                     if (sequenceIn_ >= 0 && numberNonBasic == 1) {
    1424                          // see if simple move
    1425                          double objTheta2 = objective_->stepLength(this, solution_, work, 1.0e30,
    1426                                             currentObj, predictedObj, thetaObj);
    1427                          rowArray->clear();
    1428 
    1429                          // update the incoming column
    1430                          unpackPacked(rowArray);
    1431                          factorization_->updateColumnFT(spare, rowArray);
    1432                          theta_ = 0.0;
    1433                          double * work2 = rowArray->denseVector();
    1434                          int number = rowArray->getNumElements();
    1435                          int * which2 = rowArray->getIndices();
    1436                          int iIndex;
    1437                          bool easyMove = false;
    1438                          double way;
    1439                          if (dj_[sequenceIn_] > 0.0)
    1440                               way = -1.0;
    1441                          else
    1442                               way = 1.0;
    1443                          double largest = COIN_DBL_MAX;
     1391      conjugate = (localPivotMode < 10) ? MINTYPE : 0;
     1392#endif
     1393      //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
     1394      //     nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
     1395      if (!nPasses) {
     1396        //printf("numberNon %d\n",numberNonBasic);
     1397#if MINTYPE == 2
     1398        assert(numberNonBasic < 100);
     1399        memset(zz, 0, numberNonBasic * numberNonBasic * sizeof(double));
     1400        int put = 0;
     1401        for (int iVariable = 0; iVariable < numberNonBasic; iVariable++) {
     1402          zz[put] = 1.0;
     1403          put += numberNonBasic + 1;
     1404        }
     1405#endif
     1406        djNorm0 = CoinMax(djNorm, 1.0e-20);
     1407        CoinMemcpyN(work, numberTotal, dArray);
     1408        CoinMemcpyN(work, numberTotal, dArray2);
     1409        if (sequenceIn_ >= 0 && numberNonBasic == 1) {
     1410          // see if simple move
     1411          double objTheta2 = objective_->stepLength(this, solution_, work, 1.0e30,
     1412            currentObj, predictedObj, thetaObj);
     1413          rowArray->clear();
     1414
     1415          // update the incoming column
     1416          unpackPacked(rowArray);
     1417          factorization_->updateColumnFT(spare, rowArray);
     1418          theta_ = 0.0;
     1419          double *work2 = rowArray->denseVector();
     1420          int number = rowArray->getNumElements();
     1421          int *which2 = rowArray->getIndices();
     1422          int iIndex;
     1423          bool easyMove = false;
     1424          double way;
     1425          if (dj_[sequenceIn_] > 0.0)
     1426            way = -1.0;
     1427          else
     1428            way = 1.0;
     1429          double largest = COIN_DBL_MAX;
    14441430#ifdef CLP_DEBUG
    1445                          int kPivot = -1;
    1446 #endif
    1447                          for (iIndex = 0; iIndex < number; iIndex++) {
    1448                               int iRow = which2[iIndex];
    1449                               double alpha = way * work2[iIndex];
    1450                               int iPivot = pivotVariable_[iRow];
    1451                               if (alpha < -1.0e-5) {
    1452                                    double distance = upper_[iPivot] - solution_[iPivot];
    1453                                    if (distance < -largest * alpha) {
     1431          int kPivot = -1;
     1432#endif
     1433          for (iIndex = 0; iIndex < number; iIndex++) {
     1434            int iRow = which2[iIndex];
     1435            double alpha = way * work2[iIndex];
     1436            int iPivot = pivotVariable_[iRow];
     1437            if (alpha < -1.0e-5) {
     1438              double distance = upper_[iPivot] - solution_[iPivot];
     1439              if (distance < -largest * alpha) {
    14541440#ifdef CLP_DEBUG
    1455                                         kPivot = iPivot;
    1456 #endif
    1457                                         largest = CoinMax(0.0, -distance / alpha);
    1458                                    }
    1459                                    if (distance < -1.0e-12 * alpha) {
    1460                                         easyMove = true;
    1461                                         break;
    1462                                    }
    1463                               } else if (alpha > 1.0e-5) {
    1464                                    double distance = solution_[iPivot] - lower_[iPivot];
    1465                                    if (distance < largest * alpha) {
     1441                kPivot = iPivot;
     1442#endif
     1443                largest = CoinMax(0.0, -distance / alpha);
     1444              }
     1445              if (distance < -1.0e-12 * alpha) {
     1446                easyMove = true;
     1447                break;
     1448              }
     1449            } else if (alpha > 1.0e-5) {
     1450              double distance = solution_[iPivot] - lower_[iPivot];
     1451              if (distance < largest * alpha) {
    14661452#ifdef CLP_DEBUG
    1467                                         kPivot = iPivot;
    1468 #endif
    1469                                         largest = CoinMax(0.0, distance / alpha);
    1470                                    }
    1471                                    if (distance < 1.0e-12 * alpha) {
    1472                                         easyMove = true;
    1473                                         break;
    1474                                    }
    1475                               }
    1476                          }
    1477                          // But largest has to match up with change
    1478                          assert (work[sequenceIn_]);
    1479                          largest /= fabs(work[sequenceIn_]);
    1480                          if (largest < objTheta2) {
    1481                               easyMove = true;
    1482                          } else if (!easyMove) {
    1483                               double objDrop = currentObj - predictedObj;
    1484                               double th = objective_->stepLength(this, solution_, work, largest,
    1485                                                                  currentObj, predictedObj, simpleObjective);
    1486                               simpleObjective = CoinMax(simpleObjective, predictedObj);
    1487                               double easyDrop = currentObj - simpleObjective;
    1488                               if (easyDrop > 1.0e-8 && easyDrop > 0.5 * objDrop) {
    1489                                    easyMove = true;
     1453                kPivot = iPivot;
     1454#endif
     1455                largest = CoinMax(0.0, distance / alpha);
     1456              }
     1457              if (distance < 1.0e-12 * alpha) {
     1458                easyMove = true;
     1459                break;
     1460              }
     1461            }
     1462          }
     1463          // But largest has to match up with change
     1464          assert(work[sequenceIn_]);
     1465          largest /= fabs(work[sequenceIn_]);
     1466          if (largest < objTheta2) {
     1467            easyMove = true;
     1468          } else if (!easyMove) {
     1469            double objDrop = currentObj - predictedObj;
     1470            double th = objective_->stepLength(this, solution_, work, largest,
     1471              currentObj, predictedObj, simpleObjective);
     1472            simpleObjective = CoinMax(simpleObjective, predictedObj);
     1473            double easyDrop = currentObj - simpleObjective;
     1474            if (easyDrop > 1.0e-8 && easyDrop > 0.5 * objDrop) {
     1475              easyMove = true;
    14901476#ifdef CLP_DEBUG
    1491                                    if (handler_->logLevel() & 32)
    1492                                         printf("easy - obj drop %g, easy drop %g\n", objDrop, easyDrop);
    1493 #endif
    1494                                    if (easyDrop > objDrop) {
    1495                                         // debug
    1496                                         printf("****** th %g simple %g\n", th, simpleObjective);
    1497                                         objective_->stepLength(this, solution_, work, 1.0e30,
    1498                                                                currentObj, predictedObj, simpleObjective);
    1499                                         objective_->stepLength(this, solution_, work, largest,
    1500                                                                currentObj, predictedObj, simpleObjective);
    1501                                    }
    1502                               }
    1503                          }
    1504                          rowArray->clear();
    1505 #ifdef CLP_DEBUG
    1506                          if (handler_->logLevel() & 32)
    1507                               printf("largest %g piv %d\n", largest, kPivot);
    1508 #endif
    1509                          if (easyMove) {
    1510                               valueIn_ = solution_[sequenceIn_];
    1511                               dualIn_ = dj_[sequenceIn_];
    1512                               lowerIn_ = lower_[sequenceIn_];
    1513                               upperIn_ = upper_[sequenceIn_];
    1514                               if (dualIn_ > 0.0)
    1515                                    directionIn_ = -1;
    1516                               else
    1517                                    directionIn_ = 1;
    1518                               longArray->clear();
    1519                               pivotRow_ = -1;
    1520                               theta_ = 0.0;
    1521                               return 0;
    1522                          }
    1523                     }
    1524                } else {
    1525 #if MINTYPE==1
    1526                     if (conjugate) {
    1527                          double djNorm2 = djNorm;
    1528                          if (numberNonBasic && false) {
    1529                               int iIndex;
    1530                               djNorm2 = 0.0;
    1531                               for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
    1532                                    int iSequence = which[iIndex];
    1533                                    double alpha = work[iSequence];
    1534                                    //djNorm2 += alpha*alpha;
    1535                                    double alpha2 = work[iSequence] - dArray2[iSequence];
    1536                                    djNorm2 += alpha * alpha2;
    1537                               }
    1538                               //printf("a %.18g b %.18g\n",djNorm,djNorm2);
    1539                          }
    1540                          djNorm = djNorm2;
    1541                          double beta = djNorm2 / djNormSave;
    1542                          // reset beta every so often
    1543                          //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1)
    1544                          //beta=0.0;
    1545                          //printf("current norm %g, old %g - beta %g\n",
    1546                          //     djNorm,djNormSave,beta);
    1547                          for (iSequence = 0; iSequence < numberTotal; iSequence++) {
    1548                               dArray[iSequence] = work[iSequence] + beta * dArray[iSequence];
    1549                               dArray2[iSequence] = work[iSequence];
    1550                          }
    1551                     } else {
    1552                          for (iSequence = 0; iSequence < numberTotal; iSequence++)
    1553                               dArray[iSequence] = work[iSequence];
    1554                     }
    1555 #else
    1556                     int number2 = numberNonBasic;
    1557                     if (number2) {
    1558                          // pack down into dArray
    1559                          int jLast = -1;
    1560                          for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
    1561                               int j = which[iSequence];
    1562                               assert(j > jLast);
    1563                               jLast = j;
    1564                               double value = work[j];
    1565                               dArray[iSequence] = -value;
    1566                          }
    1567                          // see whether to restart
    1568                          // check signs - gradient
    1569                          double g1 = innerProduct(dArray, number2, dArray);
    1570                          double g2 = innerProduct(dArray, number2, saveY2);
    1571                          // Get differences
    1572                          for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
    1573                               saveY2[iSequence] = dArray[iSequence] - saveY2[iSequence];
    1574                               saveS2[iSequence] = solution_[iSequence] - saveS2[iSequence];
    1575                          }
    1576                          double g3 = innerProduct(saveS2, number2, saveY2);
    1577                          printf("inner %g\n", g3);
    1578                          //assert(g3>0);
    1579                          double zzz[10000];
    1580                          int iVariable;
    1581                          g2 = 1.0e50; // temp
    1582                          if (fabs(g2) >= 0.2 * fabs(g1)) {
    1583                               // restart
    1584                               double delta = innerProduct(saveY2, number2, saveS2) /
    1585                                              innerProduct(saveY2, number2, saveY2);
    1586                               delta = 1.0; //temp
    1587                               memset(zz, 0, number2 * sizeof(double));
    1588                               int put = 0;
    1589                               for (iVariable = 0; iVariable < number2; iVariable++) {
    1590                                    zz[put] = delta;
    1591                                    put += number2 + 1;
    1592                               }
    1593                          } else {
    1594                          }
    1595                          CoinMemcpyN(zz, number2 * number2, zzz);
    1596                          double ww[100];
    1597                          // get sk -Hkyk
    1598                          for (iVariable = 0; iVariable < number2; iVariable++) {
    1599                               double value = 0.0;
    1600                               for (int jVariable = 0; jVariable < number2; jVariable++) {
    1601                                    value += saveY2[jVariable] * zzz[iVariable+number2*jVariable];
    1602                               }
    1603                               ww[iVariable] = saveS2[iVariable] - value;
    1604                          }
    1605                          double ys = innerProduct(saveY2, number2, saveS2);
    1606                          double multiplier1 = 1.0 / ys;
    1607                          double multiplier2 = innerProduct(saveY2, number2, ww) / (ys * ys);
    1608 #if 1
    1609                          // and second way
    1610                          // Hy
    1611                          double h[100];
    1612                          for (iVariable = 0; iVariable < number2; iVariable++) {
    1613                               double value = 0.0;
    1614                               for (int jVariable = 0; jVariable < number2; jVariable++) {
    1615                                    value += saveY2[jVariable] * zzz[iVariable+number2*jVariable];
    1616                               }
    1617                               h[iVariable] = value;
    1618                          }
    1619                          double hh[10000];
    1620                          double yhy1 = innerProduct(h, number2, saveY2) * multiplier1 + 1.0;
    1621                          yhy1 *= multiplier1;
    1622                          for (iVariable = 0; iVariable < number2; iVariable++) {
    1623                               for (int jVariable = 0; jVariable < number2; jVariable++) {
    1624                                    int put = iVariable + number2 * jVariable;
    1625                                    double value = zzz[put];
    1626                                    value += yhy1 * saveS2[iVariable] * saveS2[jVariable];
    1627                                    hh[put] = value;
    1628                               }
    1629                          }
    1630                          for (iVariable = 0; iVariable < number2; iVariable++) {
    1631                               for (int jVariable = 0; jVariable < number2; jVariable++) {
    1632                                    int put = iVariable + number2 * jVariable;
    1633                                    double value = hh[put];
    1634                                    value -= multiplier1 * (saveS2[iVariable] * h[jVariable]);
    1635                                    value -= multiplier1 * (saveS2[jVariable] * h[iVariable]);
    1636                                    hh[put] = value;
    1637                               }
    1638                          }
    1639 #endif
    1640                          // now update H
    1641                          for (iVariable = 0; iVariable < number2; iVariable++) {
    1642                               for (int jVariable = 0; jVariable < number2; jVariable++) {
    1643                                    int put = iVariable + number2 * jVariable;
    1644                                    double value = zzz[put];
    1645                                    value += multiplier1 * (ww[iVariable] * saveS2[jVariable]
    1646                                                            + ww[jVariable] * saveS2[iVariable]);
    1647                                    value -= multiplier2 * saveS2[iVariable] * saveS2[jVariable];
    1648                                    zzz[put] = value;
    1649                               }
    1650                          }
    1651                          //memcpy(zzz,hh,size*sizeof(double));
    1652                          // do search direction
    1653                          memset(dArray, 0, numberTotal * sizeof(double));
    1654                          for (iVariable = 0; iVariable < numberNonBasic; iVariable++) {
    1655                               double value = 0.0;
    1656                               for (int jVariable = 0; jVariable < number2; jVariable++) {
    1657                                    int k = which[jVariable];
    1658                                    value += work[k] * zzz[iVariable+number2*jVariable];
    1659                               }
    1660                               int i = which[iVariable];
    1661                               dArray[i] = value;
    1662                          }
    1663                          // Now fill out dArray
    1664                          {
    1665                               int j;
    1666                               // Now do basic ones
    1667                               int iRow;
    1668                               CoinIndexedVector * spare1 = spare;
    1669                               CoinIndexedVector * spare2 = rowArray;
    1670 #if CLP_DEBUG > 1
    1671                               spare1->checkClear();
    1672                               spare2->checkClear();
    1673 #endif
    1674                               double * array2 = spare1->denseVector();
    1675                               int * index2 = spare1->getIndices();
    1676                               int number2 = 0;
    1677                               times(-1.0, dArray, array2);
    1678                               dArray = dArray + numberColumns_;
    1679                               for (iRow = 0; iRow < numberRows_; iRow++) {
    1680                                    double value = array2[iRow] + dArray[iRow];
    1681                                    if (value) {
    1682                                         array2[iRow] = value;
    1683                                         index2[number2++] = iRow;
    1684                                    } else {
    1685                                         array2[iRow] = 0.0;
    1686                                    }
    1687                               }
    1688                               dArray -= numberColumns_;
    1689                               spare1->setNumElements(number2);
    1690                               // Ftran
    1691                               factorization_->updateColumn(spare2, spare1);
    1692                               number2 = spare1->getNumElements();
    1693                               for (j = 0; j < number2; j++) {
    1694                                    int iSequence = index2[j];
    1695                                    double value = array2[iSequence];
    1696                                    array2[iSequence] = 0.0;
    1697                                    if (value) {
    1698                                         int iPivot = pivotVariable_[iSequence];
    1699                                         double oldValue = dArray[iPivot];
    1700                                         dArray[iPivot] = value + oldValue;
    1701                                    }
    1702                               }
    1703                               spare1->setNumElements(0);
    1704                          }
    1705                          //assert (innerProductIndexed(dArray,number2,work,which)>0);
    1706                          //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
    1707                          printf ("innerP1 %g innerP2 %g\n", innerProduct(dArray, numberTotal, work),
    1708                                  innerProductIndexed(dArray, numberNonBasic, work, which));
    1709                          assert (innerProduct(dArray, numberTotal, work) > 0);
    1710 #if 1
    1711                          {
    1712                               // check null vector
    1713                               int iRow;
    1714                               double qq[10];
    1715                               memset(qq, 0, numberRows_ * sizeof(double));
    1716                               multiplyAdd(dArray + numberColumns_, numberRows_, -1.0, qq, 0.0);
    1717                               matrix_->times(1.0, dArray, qq, rowScale_, columnScale_);
    1718                               double largest = 0.0;
    1719                               int iLargest = -1;
    1720                               for (iRow = 0; iRow < numberRows_; iRow++) {
    1721                                    double value = fabs(qq[iRow]);
    1722                                    if (value > largest) {
    1723                                         largest = value;
    1724                                         iLargest = iRow;
    1725                                    }
    1726                               }
    1727                               printf("largest null error %g on row %d\n", largest, iLargest);
    1728                               for (iSequence = 0; iSequence < numberTotal; iSequence++) {
    1729                                    if (getStatus(iSequence) == basic)
    1730                                         assert (fabs(dj_[iSequence]) < 1.0e-3);
    1731                               }
    1732                          }
    1733 #endif
    1734                          CoinMemcpyN(saveY2, numberNonBasic, saveY);
    1735                          CoinMemcpyN(saveS2, numberNonBasic, saveS);
    1736                     }
    1737 #endif
    1738                }
    1739 #if MINTYPE==2
    1740                for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
    1741                     int j = which[iSequence];
    1742                     saveY2[iSequence] = -work[j];
    1743                     saveS2[iSequence] = solution_[j];
    1744                }
    1745 #endif
    1746                if (djNorm < eps * djNorm0 || (nPasses > 100 && djNorm < CoinMin(1.0e-1 * djNorm0, 1.0e-12))) {
    1747 #ifdef CLP_DEBUG
    1748                     if (handler_->logLevel() & 32)
    1749                          printf("dj norm reduced from %g to %g\n", djNorm0, djNorm);
    1750 #endif
    1751                     if (pivotMode < 10 || !numberNonBasic) {
    1752                          finished = true;
    1753                     } else {
    1754                          localPivotMode = pivotMode;
    1755                          nTotalPasses += nPasses;
    1756                          nPasses = 0;
    1757                          continue;
    1758                     }
    1759                }
    1760                //if (nPasses==100||nPasses==50)
    1761                //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
    1762                //        normFlagged);
    1763                if (nPasses > 100 && djNorm < 1.0e-2 * normFlagged && !startLocalMode) {
    1764 #ifdef CLP_DEBUG
    1765                     if (handler_->logLevel() & 32)
    1766                          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n", djNorm0, djNorm,
    1767                                 normFlagged);
    1768 #endif
    1769                     unflag();
    1770                     localPivotMode = 0;
    1771                     nTotalPasses += nPasses;
    1772                     nPasses = 0;
    1773                     continue;
    1774                }
    1775                if (djNorm > 0.99 * djNorm0 && nPasses > 1500) {
    1776                     finished = true;
    1777 #ifdef CLP_DEBUG
    1778                     if (handler_->logLevel() & 32)
    1779                          printf("dj norm NOT reduced from %g to %g\n", djNorm0, djNorm);
    1780 #endif
    1781                     djNorm = 1.2345e-20;
    1782                }
    1783                number = longArray->getNumElements();
    1784                if (!numberNonBasic) {
    1785                     // looks optimal
    1786                     // check if any dj look plausible
    1787                     int nSuper = 0;
    1788                     int nFlagged = 0;
    1789                     int nNormal = 0;
    1790                     for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
    1791                          if (flagged(iSequence)) {
    1792                               switch(getStatus(iSequence)) {
    1793 
    1794                               case basic:
    1795                               case ClpSimplex::isFixed:
    1796                                    break;
    1797                               case atUpperBound:
    1798                                    if (dj_[iSequence] > dualTolerance_)
    1799                                         nFlagged++;
    1800                                    break;
    1801                               case atLowerBound:
    1802                                    if (dj_[iSequence] < -dualTolerance_)
    1803                                         nFlagged++;
    1804                                    break;
    1805                               case isFree:
    1806                               case superBasic:
    1807                                    if (fabs(dj_[iSequence]) > dualTolerance_)
    1808                                         nFlagged++;
    1809                                    break;
    1810                               }
    1811                               continue;
    1812                          }
    1813                          switch(getStatus(iSequence)) {
    1814 
    1815                          case basic:
    1816                          case ClpSimplex::isFixed:
    1817                               break;
    1818                          case atUpperBound:
    1819                               if (dj_[iSequence] > dualTolerance_)
    1820                                    nNormal++;
    1821                               break;
    1822                          case atLowerBound:
    1823                               if (dj_[iSequence] < -dualTolerance_)
    1824                                    nNormal++;
    1825                               break;
    1826                          case isFree:
    1827                          case superBasic:
    1828                               if (fabs(dj_[iSequence]) > dualTolerance_)
    1829                                    nSuper++;
    1830                               break;
    1831                          }
    1832                     }
    1833 #ifdef CLP_DEBUG
    1834                     if (handler_->logLevel() & 32)
    1835                          printf("%d super, %d normal, %d flagged\n",
    1836                                 nSuper, nNormal, nFlagged);
    1837 #endif
    1838 
    1839                     int nFlagged2 = 1;
    1840                     if (lastSequenceIn < 0 && !nNormal && !nSuper) {
    1841                          nFlagged2 = unflag();
    1842                          if (pivotMode >= 10) {
    1843                               pivotMode--;
    1844 #ifdef CLP_DEBUG
    1845                               if (handler_->logLevel() & 32)
    1846                                    printf("pivot mode now %d\n", pivotMode);
    1847 #endif
    1848                               if (pivotMode == 9)
    1849                                    pivotMode = 0;       // switch off fast attempt
    1850                          }
    1851                     } else {
    1852                          lastSequenceIn = -1;
    1853                     }
    1854                     if (!nFlagged2 || (normFlagged + normUnflagged < 1.0e-8)) {
    1855                          primalColumnPivot_->setLooksOptimal(true);
    1856                          return 0;
    1857                     } else {
    1858                          localPivotMode = -1;
    1859                          nTotalPasses += nPasses;
    1860                          nPasses = 0;
    1861                          finished = false;
    1862                          continue;
    1863                     }
    1864                }
    1865                bestSequence = -1;
    1866                bestBasicSequence = -1;
    1867 
    1868                // temp
    1869                nPasses++;
    1870                if (nPasses > 2000)
    1871                     finished = true;
    1872                double theta = 1.0e30;
    1873                basicTheta = 1.0e30;
    1874                theta_ = 1.0e30;
    1875                double basicTolerance = 1.0e-4 * primalTolerance_;
    1876                for (iSequence = 0; iSequence < numberTotal; iSequence++) {
    1877                     //if (flagged(iSequence)
    1878                     //  continue;
    1879                     double alpha = dArray[iSequence];
    1880                     Status thisStatus = getStatus(iSequence);
    1881                     double oldValue = solution_[iSequence];
    1882                     if (thisStatus != basic) {
    1883                          if (fabs(alpha) >= acceptablePivot) {
    1884                               if (alpha < 0.0) {
    1885                                    // variable going towards lower bound
    1886                                    double bound = lower_[iSequence];
    1887                                    oldValue -= bound;
    1888                                    if (oldValue + theta * alpha < 0.0) {
    1889                                         bestSequence = iSequence;
    1890                                         theta = CoinMax(0.0, oldValue / (-alpha));
    1891                                    }
    1892                               } else {
    1893                                    // variable going towards upper bound
    1894                                    double bound = upper_[iSequence];
    1895                                    oldValue = bound - oldValue;
    1896                                    if (oldValue - theta * alpha < 0.0) {
    1897                                         bestSequence = iSequence;
    1898                                         theta = CoinMax(0.0, oldValue / alpha);
    1899                                    }
    1900                               }
    1901                          }
    1902                     } else {
    1903                          if (fabs(alpha) >= acceptableBasic) {
    1904                               if (alpha < 0.0) {
    1905                                    // variable going towards lower bound
    1906                                    double bound = lower_[iSequence];
    1907                                    oldValue -= bound;
    1908                                    if (oldValue + basicTheta * alpha < -basicTolerance) {
    1909                                         bestBasicSequence = iSequence;
    1910                                         basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / (-alpha));
    1911                                    }
    1912                               } else {
    1913                                    // variable going towards upper bound
    1914                                    double bound = upper_[iSequence];
    1915                                    oldValue = bound - oldValue;
    1916                                    if (oldValue - basicTheta * alpha < -basicTolerance) {
    1917                                         bestBasicSequence = iSequence;
    1918                                         basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / alpha);
    1919                                    }
    1920                               }
    1921                          }
    1922                     }
    1923                }
    1924                theta_ = CoinMin(theta, basicTheta);
    1925                // Now find minimum of function
    1926                double objTheta2 = objective_->stepLength(this, solution_, dArray, CoinMin(theta, basicTheta),
    1927                                   currentObj, predictedObj, thetaObj);
    1928 #ifdef CLP_DEBUG
    1929                if (handler_->logLevel() & 32)
    1930                     printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
    1931 #endif
    1932                objTheta2=CoinMin(objTheta2,1.0e29);
    1933 #if MINTYPE==1
    1934                if (conjugate) {
    1935                     double offset;
    1936                     const double * gradient = objective_->gradient(this,
    1937                                               dArray, offset,
    1938                                               true, 0);
    1939                     double product = 0.0;
    1940                     for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
    1941                          double alpha = dArray[iSequence];
    1942                          double value = alpha * gradient[iSequence];
    1943                          product += value;
    1944                     }
    1945                     //#define INCLUDESLACK
    1946 #ifdef INCLUDESLACK
    1947                     for (; iSequence < numberColumns_ + numberRows_; iSequence++) {
    1948                          double alpha = dArray[iSequence];
    1949                          double value = alpha * cost_[iSequence];
    1950                          product += value;
    1951                     }
    1952 #endif
    1953                     if (product > 0.0)
    1954                          objTheta = djNorm / product;
    1955                     else
    1956                          objTheta = COIN_DBL_MAX;
    1957 #ifndef NDEBUG
    1958                     if (product < -1.0e-8 && handler_->logLevel() > 1)
    1959                          printf("bad product %g\n", product);
    1960 #endif
    1961                     product = CoinMax(product, 0.0);
    1962                } else {
    1963                     objTheta = objTheta2;
    1964                }
    1965 #else
    1966                objTheta = objTheta2;
    1967 #endif
    1968                // if very small difference then take pivot (depends on djNorm?)
    1969                // distinguish between basic and non-basic
    1970                bool chooseObjTheta = objTheta < theta_;
    1971                if (chooseObjTheta) {
    1972                     if (thetaObj < currentObj - 1.0e-12 && objTheta + 1.0e-10 > theta_)
    1973                          chooseObjTheta = false;
    1974                     //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
    1975                     //chooseObjTheta=false;
    1976                }
    1977                //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
    1978                if (chooseObjTheta) {
    1979                     theta_ = objTheta;
    1980                } else {
    1981                     objTheta = CoinMax(objTheta, 1.00000001 * theta_ + 1.0e-12);
    1982                     //if (theta+1.0e-13>basicTheta) {
    1983                     //theta = CoinMax(theta,1.00000001*basicTheta);
    1984                     //theta_ = basicTheta;
    1985                     //}
    1986                }
    1987                // always out if one variable in and zero move
    1988                if (theta_ == basicTheta || (sequenceIn_ >= 0 && theta_ < 1.0e-10))
    1989                     finished = true; // come out
    1990 #ifdef CLP_DEBUG
    1991                if (handler_->logLevel() & 32) {
    1992                     printf("%d pass,", nPasses);
    1993                     if (sequenceIn_ >= 0)
    1994                          printf (" Sin %d,", sequenceIn_);
    1995                     if (basicTheta == theta_)
    1996                          printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
    1997                     else
    1998                          printf(" basicTheta %g", basicTheta);
    1999                     if (theta == theta_)
    2000                          printf(" X(%d) non-basicTheta %g", bestSequence, theta);
    2001                     else
    2002                          printf(" non-basicTheta %g", theta);
    2003                     printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
    2004                     printf(" djNorm %g\n", djNorm);
    2005                }
    2006 #endif
    2007                if (handler_->logLevel() > 3 && objTheta != theta_) {
    2008                     printf("%d pass obj %g,", nPasses, currentObj);
    2009                     if (sequenceIn_ >= 0)
    2010                          printf (" Sin %d,", sequenceIn_);
    2011                     if (basicTheta == theta_)
    2012                          printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
    2013                     else
    2014                          printf(" basicTheta %g", basicTheta);
    2015                     if (theta == theta_)
    2016                          printf(" X(%d) non-basicTheta %g", bestSequence, theta);
    2017                     else
    2018                          printf(" non-basicTheta %g", theta);
    2019                     printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
    2020                     printf(" djNorm %g\n", djNorm);
    2021                }
    2022                if (objTheta != theta_) {
    2023                     //printf("hit boundary after %d passes\n",nPasses);
    2024                     nTotalPasses += nPasses;
    2025                     nPasses = 0; // start again
    2026                }
    2027                if (localPivotMode < 10 || lastSequenceIn == bestSequence) {
    2028                     if (theta_ == theta && theta_ < basicTheta && theta_ < 1.0e-5)
    2029                          setFlagged(bestSequence); // out of active set
    2030                }
    2031                // Update solution
    2032                for (iSequence = 0; iSequence < numberTotal; iSequence++) {
    2033                     //for (iIndex=0;iIndex<number;iIndex++) {
    2034 
    2035                     //int iSequence = which[iIndex];
    2036                     double alpha = dArray[iSequence];
    2037                     if (alpha) {
    2038                          double value = solution_[iSequence] + theta_ * alpha;
    2039                          solution_[iSequence] = value;
    2040                          switch(getStatus(iSequence)) {
    2041 
    2042                          case basic:
    2043                          case isFixed:
    2044                          case isFree:
    2045                          case atUpperBound:
    2046                          case atLowerBound:
    2047                               nonLinearCost_->setOne(iSequence, value);
    2048                               break;
    2049                          case superBasic:
    2050                               // To get correct action
    2051                               setStatus(iSequence, isFixed);
    2052                               nonLinearCost_->setOne(iSequence, value);
    2053                               //assert (getStatus(iSequence)!=isFixed);
    2054                               break;
    2055                          }
    2056                     }
    2057                }
    2058                if (objTheta < SMALLTHETA2 && objTheta == theta_) {
    2059                     if (sequenceIn_ >= 0 && basicTheta < 1.0e-9) {
    2060                          // try just one
    2061                          localPivotMode = 0;
    2062                          sequenceIn_ = -1;
    2063                          nTotalPasses += nPasses;
    2064                          nPasses = 0;
    2065                          //finished=true;
    2066                          //objTheta=0.0;
    2067                          //theta_=0.0;
    2068                     } else if (sequenceIn_ < 0 && nTotalPasses > 10) {
    2069                          if (objTheta < 1.0e-10) {
    2070                               finished = true;
    2071                               //printf("zero move\n");
    2072                               break;
    2073                          }
    2074                     }
    2075                }
    2076 #ifdef CLP_DEBUG
    2077                if (handler_->logLevel() & 32) {
    2078                     objective_->stepLength(this, solution_, work, 0.0,
    2079                                            currentObj, predictedObj, thetaObj);
    2080                     printf("current obj %g after update - simple was %g\n", currentObj, simpleObjective);
    2081                }
    2082 #endif
    2083                if (sequenceIn_ >= 0 && !finished && objTheta > 1.0e-4) {
    2084                     // we made some progress - back to normal
    2085                     if (localPivotMode < 10) {
    2086                          localPivotMode = 0;
    2087                          sequenceIn_ = -1;
    2088                          nTotalPasses += nPasses;
    2089                          nPasses = 0;
    2090                     }
    2091 #ifdef CLP_DEBUG
    2092                     if (handler_->logLevel() & 32)
    2093                          printf("some progress?\n");
    2094 #endif
    2095                }
    2096 #if CLP_DEBUG > 1
    2097                longArray->checkClean();
    2098 #endif
    2099           }
     1477              if (handler_->logLevel() & 32)
     1478                printf("easy - obj drop %g, easy drop %g\n", objDrop, easyDrop);
     1479#endif
     1480              if (easyDrop > objDrop) {
     1481                // debug
     1482                printf("****** th %g simple %g\n", th, simpleObjective);
     1483                objective_->stepLength(this, solution_, work, 1.0e30,
     1484                  currentObj, predictedObj, simpleObjective);
     1485                objective_->stepLength(this, solution_, work, largest,
     1486                  currentObj, predictedObj, simpleObjective);
     1487              }
     1488            }
     1489          }
     1490          rowArray->clear();
    21001491#ifdef CLP_DEBUG
    21011492          if (handler_->logLevel() & 32)
    2102                printf("out of loop after %d (%d) passes\n", nPasses, nTotalPasses);
    2103 #endif
    2104           if (nTotalPasses >= 1000 || (nTotalPasses > 10 && sequenceIn_ < 0 && theta_ < 1.0e-10))
    2105                returnCode = 2;
    2106           bool ordinaryDj = false;
    2107           //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
    2108           //printf("could try ordinary iteration %g\n",theta_);
    2109           if(sequenceIn_ >= 0 && numberNonBasic == 1 && theta_ < 1.0e-15) {
    2110                //printf("trying ordinary iteration\n");
    2111                ordinaryDj = true;
    2112           }
    2113           if (!basicTheta && !ordinaryDj) {
    2114                //returnCode=2;
    2115                //objTheta=-1.0; // so we fall through
    2116           }
    2117           assert (theta_ < 1.0e30); // for now
    2118           // See if we need to pivot
    2119           if (theta_ == basicTheta || ordinaryDj) {
    2120                if (!ordinaryDj) {
    2121                     // find an incoming column which will force pivot
    2122                     int iRow;
    2123                     pivotRow_ = -1;
    2124                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2125                          if (pivotVariable_[iRow] == bestBasicSequence) {
    2126                               pivotRow_ = iRow;
    2127                               break;
    2128                          }
    2129                     }
    2130                     assert (pivotRow_ >= 0);
    2131                     // Get good size for pivot
    2132                     double acceptablePivot = 1.0e-7;
    2133                     if (factorization_->pivots() > 10)
    2134                          acceptablePivot = 1.0e-5; // if we have iterated be more strict
    2135                     else if (factorization_->pivots() > 5)
    2136                          acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
    2137                     // should be dArray but seems better this way!
    2138                     double direction = work[bestBasicSequence] > 0.0 ? -1.0 : 1.0;
    2139                     // create as packed
    2140                     rowArray->createPacked(1, &pivotRow_, &direction);
    2141                     factorization_->updateColumnTranspose(spare, rowArray);
    2142                     // put row of tableau in rowArray and columnArray
    2143                     matrix_->transposeTimes(this, -1.0,
    2144                                             rowArray, spare, columnArray);
    2145                     // choose one futhest away from bound which has reasonable pivot
    2146                     // If increasing we want negative alpha
    2147 
    2148                     double * work2;
    2149                     int iSection;
    2150 
    2151                     sequenceIn_ = -1;
    2152                     double bestValue = -1.0;
    2153                     double bestDirection = 0.0;
    2154                     // First pass we take correct direction and large pivots
    2155                     // then correct direction
    2156                     // then any
    2157                     double check[] = {1.0e-8, -1.0e-12, -1.0e30};
    2158                     double mult[] = {100.0, 1.0, 1.0};
    2159                     for (int iPass = 0; iPass < 3; iPass++) {
    2160                          //if (!bestValue&&iPass==2)
    2161                          //bestValue=-1.0;
    2162                          double acceptable = acceptablePivot * mult[iPass];
    2163                          double checkValue = check[iPass];
    2164                          for (iSection = 0; iSection < 2; iSection++) {
    2165 
    2166                               int addSequence;
    2167 
    2168                               if (!iSection) {
    2169                                    work2 = rowArray->denseVector();
    2170                                    number = rowArray->getNumElements();
    2171                                    which = rowArray->getIndices();
    2172                                    addSequence = numberColumns_;
    2173                               } else {
    2174                                    work2 = columnArray->denseVector();
    2175                                    number = columnArray->getNumElements();
    2176                                    which = columnArray->getIndices();
    2177                                    addSequence = 0;
    2178                               }
    2179                               int i;
    2180 
    2181                               for (i = 0; i < number; i++) {
    2182                                    int iSequence = which[i] + addSequence;
    2183                                    if (flagged(iSequence))
    2184                                         continue;
    2185                                    //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
    2186                                    //             upper_[iSequence]-solution_[iSequence]);
    2187                                    double alpha = work2[i];
    2188                                    // should be dArray but seems better this way!
    2189                                    double change = work[iSequence];
    2190                                    Status thisStatus = getStatus(iSequence);
    2191                                    double direction = 0;;
    2192                                    switch(thisStatus) {
    2193 
    2194                                    case basic:
    2195                                    case ClpSimplex::isFixed:
    2196                                         break;
    2197                                    case isFree:
    2198                                    case superBasic:
    2199                                         if (alpha < -acceptable && change > checkValue)
    2200                                              direction = 1.0;
    2201                                         else if (alpha > acceptable && change < -checkValue)
    2202                                              direction = -1.0;
    2203                                         break;
    2204                                    case atUpperBound:
    2205                                         if (alpha > acceptable && change < -checkValue)
    2206                                              direction = -1.0;
    2207                                         else if (iPass == 2 && alpha < -acceptable && change < -checkValue)
    2208                                              direction = 1.0;
    2209                                         break;
    2210                                    case atLowerBound:
    2211                                         if (alpha < -acceptable && change > checkValue)
    2212                                              direction = 1.0;
    2213                                         else if (iPass == 2 && alpha > acceptable && change > checkValue)
    2214                                              direction = -1.0;
    2215                                         break;
    2216                                    }
    2217                                    if (direction) {
    2218                                         if (sequenceIn_ != lastSequenceIn || localPivotMode < 10) {
    2219                                              if (CoinMin(solution_[iSequence] - lower_[iSequence],
    2220                                                          upper_[iSequence] - solution_[iSequence]) > bestValue) {
    2221                                                   bestValue = CoinMin(solution_[iSequence] - lower_[iSequence],
    2222                                                                       upper_[iSequence] - solution_[iSequence]);
    2223                                                   sequenceIn_ = iSequence;
    2224                                                   bestDirection = direction;
    2225                                              }
    2226                                         } else {
    2227                                              // choose
    2228                                              bestValue = COIN_DBL_MAX;
    2229                                              sequenceIn_ = iSequence;
    2230                                              bestDirection = direction;
    2231                                         }
    2232                                    }
    2233                               }
    2234                          }
    2235                          if (sequenceIn_ >= 0 && bestValue > 0.0)
    2236                               break;
    2237                     }
    2238                     if (sequenceIn_ >= 0) {
    2239                          valueIn_ = solution_[sequenceIn_];
    2240                          dualIn_ = dj_[sequenceIn_];
    2241                          if (bestDirection < 0.0) {
    2242                               // we want positive dj
    2243                               if (dualIn_ <= 0.0) {
    2244                                    //printf("bad dj - xx %g\n",dualIn_);
    2245                                    dualIn_ = 1.0;
    2246                               }
    2247                          } else {
    2248                               // we want negative dj
    2249                               if (dualIn_ >= 0.0) {
    2250                                    //printf("bad dj - xx %g\n",dualIn_);
    2251                                    dualIn_ = -1.0;
    2252                               }
    2253                          }
    2254                          lowerIn_ = lower_[sequenceIn_];
    2255                          upperIn_ = upper_[sequenceIn_];
    2256                          if (dualIn_ > 0.0)
    2257                               directionIn_ = -1;
    2258                          else
    2259                               directionIn_ = 1;
    2260                     } else {
    2261                          //ordinaryDj=true;
     1493            printf("largest %g piv %d\n", largest, kPivot);
     1494#endif
     1495          if (easyMove) {
     1496            valueIn_ = solution_[sequenceIn_];
     1497            dualIn_ = dj_[sequenceIn_];
     1498            lowerIn_ = lower_[sequenceIn_];
     1499            upperIn_ = upper_[sequenceIn_];
     1500            if (dualIn_ > 0.0)
     1501              directionIn_ = -1;
     1502            else
     1503              directionIn_ = 1;
     1504            longArray->clear();
     1505            pivotRow_ = -1;
     1506            theta_ = 0.0;
     1507            return 0;
     1508          }
     1509        }
     1510      } else {
     1511#if MINTYPE == 1
     1512        if (conjugate) {
     1513          double djNorm2 = djNorm;
     1514          if (numberNonBasic && false) {
     1515            int iIndex;
     1516            djNorm2 = 0.0;
     1517            for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
     1518              int iSequence = which[iIndex];
     1519              double alpha = work[iSequence];
     1520              //djNorm2 += alpha*alpha;
     1521              double alpha2 = work[iSequence] - dArray2[iSequence];
     1522              djNorm2 += alpha * alpha2;
     1523            }
     1524            //printf("a %.18g b %.18g\n",djNorm,djNorm2);
     1525          }
     1526          djNorm = djNorm2;
     1527          double beta = djNorm2 / djNormSave;
     1528          // reset beta every so often
     1529          //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1)
     1530          //beta=0.0;
     1531          //printf("current norm %g, old %g - beta %g\n",
     1532          //     djNorm,djNormSave,beta);
     1533          for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     1534            dArray[iSequence] = work[iSequence] + beta * dArray[iSequence];
     1535            dArray2[iSequence] = work[iSequence];
     1536          }
     1537        } else {
     1538          for (iSequence = 0; iSequence < numberTotal; iSequence++)
     1539            dArray[iSequence] = work[iSequence];
     1540        }
     1541#else
     1542        int number2 = numberNonBasic;
     1543        if (number2) {
     1544          // pack down into dArray
     1545          int jLast = -1;
     1546          for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
     1547            int j = which[iSequence];
     1548            assert(j > jLast);
     1549            jLast = j;
     1550            double value = work[j];
     1551            dArray[iSequence] = -value;
     1552          }
     1553          // see whether to restart
     1554          // check signs - gradient
     1555          double g1 = innerProduct(dArray, number2, dArray);
     1556          double g2 = innerProduct(dArray, number2, saveY2);
     1557          // Get differences
     1558          for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
     1559            saveY2[iSequence] = dArray[iSequence] - saveY2[iSequence];
     1560            saveS2[iSequence] = solution_[iSequence] - saveS2[iSequence];
     1561          }
     1562          double g3 = innerProduct(saveS2, number2, saveY2);
     1563          printf("inner %g\n", g3);
     1564          //assert(g3>0);
     1565          double zzz[10000];
     1566          int iVariable;
     1567          g2 = 1.0e50; // temp
     1568          if (fabs(g2) >= 0.2 * fabs(g1)) {
     1569            // restart
     1570            double delta = innerProduct(saveY2, number2, saveS2) / innerProduct(saveY2, number2, saveY2);
     1571            delta = 1.0; //temp
     1572            memset(zz, 0, number2 * sizeof(double));
     1573            int put = 0;
     1574            for (iVariable = 0; iVariable < number2; iVariable++) {
     1575              zz[put] = delta;
     1576              put += number2 + 1;
     1577            }
     1578          } else {
     1579          }
     1580          CoinMemcpyN(zz, number2 * number2, zzz);
     1581          double ww[100];
     1582          // get sk -Hkyk
     1583          for (iVariable = 0; iVariable < number2; iVariable++) {
     1584            double value = 0.0;
     1585            for (int jVariable = 0; jVariable < number2; jVariable++) {
     1586              value += saveY2[jVariable] * zzz[iVariable + number2 * jVariable];
     1587            }
     1588            ww[iVariable] = saveS2[iVariable] - value;
     1589          }
     1590          double ys = innerProduct(saveY2, number2, saveS2);
     1591          double multiplier1 = 1.0 / ys;
     1592          double multiplier2 = innerProduct(saveY2, number2, ww) / (ys * ys);
     1593#if 1
     1594          // and second way
     1595          // Hy
     1596          double h[100];
     1597          for (iVariable = 0; iVariable < number2; iVariable++) {
     1598            double value = 0.0;
     1599            for (int jVariable = 0; jVariable < number2; jVariable++) {
     1600              value += saveY2[jVariable] * zzz[iVariable + number2 * jVariable];
     1601            }
     1602            h[iVariable] = value;
     1603          }
     1604          double hh[10000];
     1605          double yhy1 = innerProduct(h, number2, saveY2) * multiplier1 + 1.0;
     1606          yhy1 *= multiplier1;
     1607          for (iVariable = 0; iVariable < number2; iVariable++) {
     1608            for (int jVariable = 0; jVariable < number2; jVariable++) {
     1609              int put = iVariable + number2 * jVariable;
     1610              double value = zzz[put];
     1611              value += yhy1 * saveS2[iVariable] * saveS2[jVariable];
     1612              hh[put] = value;
     1613            }
     1614          }
     1615          for (iVariable = 0; iVariable < number2; iVariable++) {
     1616            for (int jVariable = 0; jVariable < number2; jVariable++) {
     1617              int put = iVariable + number2 * jVariable;
     1618              double value = hh[put];
     1619              value -= multiplier1 * (saveS2[iVariable] * h[jVariable]);
     1620              value -= multiplier1 * (saveS2[jVariable] * h[iVariable]);
     1621              hh[put] = value;
     1622            }
     1623          }
     1624#endif
     1625          // now update H
     1626          for (iVariable = 0; iVariable < number2; iVariable++) {
     1627            for (int jVariable = 0; jVariable < number2; jVariable++) {
     1628              int put = iVariable + number2 * jVariable;
     1629              double value = zzz[put];
     1630              value += multiplier1 * (ww[iVariable] * saveS2[jVariable] + ww[jVariable] * saveS2[iVariable]);
     1631              value -= multiplier2 * saveS2[iVariable] * saveS2[jVariable];
     1632              zzz[put] = value;
     1633            }
     1634          }
     1635          //memcpy(zzz,hh,size*sizeof(double));
     1636          // do search direction
     1637          memset(dArray, 0, numberTotal * sizeof(double));
     1638          for (iVariable = 0; iVariable < numberNonBasic; iVariable++) {
     1639            double value = 0.0;
     1640            for (int jVariable = 0; jVariable < number2; jVariable++) {
     1641              int k = which[jVariable];
     1642              value += work[k] * zzz[iVariable + number2 * jVariable];
     1643            }
     1644            int i = which[iVariable];
     1645            dArray[i] = value;
     1646          }
     1647          // Now fill out dArray
     1648          {
     1649            int j;
     1650            // Now do basic ones
     1651            int iRow;
     1652            CoinIndexedVector *spare1 = spare;
     1653            CoinIndexedVector *spare2 = rowArray;
     1654#if CLP_DEBUG > 1
     1655            spare1->checkClear();
     1656            spare2->checkClear();
     1657#endif
     1658            double *array2 = spare1->denseVector();
     1659            int *index2 = spare1->getIndices();
     1660            int number2 = 0;
     1661            times(-1.0, dArray, array2);
     1662            dArray = dArray + numberColumns_;
     1663            for (iRow = 0; iRow < numberRows_; iRow++) {
     1664              double value = array2[iRow] + dArray[iRow];
     1665              if (value) {
     1666                array2[iRow] = value;
     1667                index2[number2++] = iRow;
     1668              } else {
     1669                array2[iRow] = 0.0;
     1670              }
     1671            }
     1672            dArray -= numberColumns_;
     1673            spare1->setNumElements(number2);
     1674            // Ftran
     1675            factorization_->updateColumn(spare2, spare1);
     1676            number2 = spare1->getNumElements();
     1677            for (j = 0; j < number2; j++) {
     1678              int iSequence = index2[j];
     1679              double value = array2[iSequence];
     1680              array2[iSequence] = 0.0;
     1681              if (value) {
     1682                int iPivot = pivotVariable_[iSequence];
     1683                double oldValue = dArray[iPivot];
     1684                dArray[iPivot] = value + oldValue;
     1685              }
     1686            }
     1687            spare1->setNumElements(0);
     1688          }
     1689          //assert (innerProductIndexed(dArray,number2,work,which)>0);
     1690          //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
     1691          printf("innerP1 %g innerP2 %g\n", innerProduct(dArray, numberTotal, work),
     1692            innerProductIndexed(dArray, numberNonBasic, work, which));
     1693          assert(innerProduct(dArray, numberTotal, work) > 0);
     1694#if 1
     1695          {
     1696            // check null vector
     1697            int iRow;
     1698            double qq[10];
     1699            memset(qq, 0, numberRows_ * sizeof(double));
     1700            multiplyAdd(dArray + numberColumns_, numberRows_, -1.0, qq, 0.0);
     1701            matrix_->times(1.0, dArray, qq, rowScale_, columnScale_);
     1702            double largest = 0.0;
     1703            int iLargest = -1;
     1704            for (iRow = 0; iRow < numberRows_; iRow++) {
     1705              double value = fabs(qq[iRow]);
     1706              if (value > largest) {
     1707                largest = value;
     1708                iLargest = iRow;
     1709              }
     1710            }
     1711            printf("largest null error %g on row %d\n", largest, iLargest);
     1712            for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     1713              if (getStatus(iSequence) == basic)
     1714                assert(fabs(dj_[iSequence]) < 1.0e-3);
     1715            }
     1716          }
     1717#endif
     1718          CoinMemcpyN(saveY2, numberNonBasic, saveY);
     1719          CoinMemcpyN(saveS2, numberNonBasic, saveS);
     1720        }
     1721#endif
     1722      }
     1723#if MINTYPE == 2
     1724      for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
     1725        int j = which[iSequence];
     1726        saveY2[iSequence] = -work[j];
     1727        saveS2[iSequence] = solution_[j];
     1728      }
     1729#endif
     1730      if (djNorm < eps * djNorm0 || (nPasses > 100 && djNorm < CoinMin(1.0e-1 * djNorm0, 1.0e-12))) {
    22621731#ifdef CLP_DEBUG
    2263                          if (handler_->logLevel() & 32) {
    2264                               printf("no easy pivot - norm %g mode %d\n", djNorm, localPivotMode);
    2265                               if (rowArray->getNumElements() + columnArray->getNumElements() < 12) {
    2266                                    for (iSection = 0; iSection < 2; iSection++) {
    2267 
    2268                                         int addSequence;
    2269 
    2270                                         if (!iSection) {
    2271                                              work2 = rowArray->denseVector();
    2272                                              number = rowArray->getNumElements();
    2273                                              which = rowArray->getIndices();
    2274                                              addSequence = numberColumns_;
    2275                                         } else {
    2276                                              work2 = columnArray->denseVector();
    2277                                              number = columnArray->getNumElements();
    2278                                              which = columnArray->getIndices();
    2279                                              addSequence = 0;
    2280                                         }
    2281                                         int i;
    2282                                         char section[] = {'R', 'C'};
    2283                                         for (i = 0; i < number; i++) {
    2284                                              int iSequence = which[i] + addSequence;
    2285                                              if (flagged(iSequence)) {
    2286                                                   printf("%c%d flagged\n", section[iSection], which[i]);
    2287                                                   continue;
    2288                                              } else {
    2289                                                   printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
    2290                                                          section[iSection], which[i], status_[iSequence],
    2291                                                          lower_[iSequence], solution_[iSequence], upper_[iSequence],
    2292                                                          work2[i], work[iSequence]);
    2293                                              }
    2294                                         }
    2295                                    }
    2296                               }
    2297                          }
    2298 #endif
    2299                          assert (sequenceIn_ < 0);
    2300                          justOne = true;
    2301                          allFinished = false; // Round again
    2302                          finished = false;
    2303                          nTotalPasses += nPasses;
    2304                          nPasses = 0;
    2305                          if (djNorm < 0.9 * djNorm0 && djNorm < 1.0e-3 && !localPivotMode) {
     1732        if (handler_->logLevel() & 32)
     1733          printf("dj norm reduced from %g to %g\n", djNorm0, djNorm);
     1734#endif
     1735        if (pivotMode < 10 || !numberNonBasic) {
     1736          finished = true;
     1737        } else {
     1738          localPivotMode = pivotMode;
     1739          nTotalPasses += nPasses;
     1740          nPasses = 0;
     1741          continue;
     1742        }
     1743      }
     1744      //if (nPasses==100||nPasses==50)
     1745      //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
     1746      //         normFlagged);
     1747      if (nPasses > 100 && djNorm < 1.0e-2 * normFlagged && !startLocalMode) {
    23061748#ifdef CLP_DEBUG
    2307                               if (handler_->logLevel() & 32)
    2308                                    printf("no pivot - mode %d norms %g %g - unflagging\n",
    2309                                           localPivotMode, djNorm0, djNorm);
    2310 #endif
    2311                               unflag(); //unflagging
    2312                               returnCode = 1;
    2313                          } else {
    2314                               returnCode = 2; // do single incoming
    2315                               returnCode = 1;
    2316                          }
    2317                     }
    2318                }
    2319                rowArray->clear();
    2320                columnArray->clear();
    2321                longArray->clear();
    2322                if (ordinaryDj) {
    2323                     valueIn_ = solution_[sequenceIn_];
    2324                     dualIn_ = dj_[sequenceIn_];
    2325                     lowerIn_ = lower_[sequenceIn_];
    2326                     upperIn_ = upper_[sequenceIn_];
    2327                     if (dualIn_ > 0.0)
    2328                          directionIn_ = -1;
    2329                     else
    2330                          directionIn_ = 1;
    2331                }
    2332                if (returnCode == 1)
    2333                     returnCode = 0;
     1749        if (handler_->logLevel() & 32)
     1750          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n", djNorm0, djNorm,
     1751            normFlagged);
     1752#endif
     1753        unflag();
     1754        localPivotMode = 0;
     1755        nTotalPasses += nPasses;
     1756        nPasses = 0;
     1757        continue;
     1758      }
     1759      if (djNorm > 0.99 * djNorm0 && nPasses > 1500) {
     1760        finished = true;
     1761#ifdef CLP_DEBUG
     1762        if (handler_->logLevel() & 32)
     1763          printf("dj norm NOT reduced from %g to %g\n", djNorm0, djNorm);
     1764#endif
     1765        djNorm = 1.2345e-20;
     1766      }
     1767      number = longArray->getNumElements();
     1768      if (!numberNonBasic) {
     1769        // looks optimal
     1770        // check if any dj look plausible
     1771        int nSuper = 0;
     1772        int nFlagged = 0;
     1773        int nNormal = 0;
     1774        for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
     1775          if (flagged(iSequence)) {
     1776            switch (getStatus(iSequence)) {
     1777
     1778            case basic:
     1779            case ClpSimplex::isFixed:
     1780              break;
     1781            case atUpperBound:
     1782              if (dj_[iSequence] > dualTolerance_)
     1783                nFlagged++;
     1784              break;
     1785            case atLowerBound:
     1786              if (dj_[iSequence] < -dualTolerance_)
     1787                nFlagged++;
     1788              break;
     1789            case isFree:
     1790            case superBasic:
     1791              if (fabs(dj_[iSequence]) > dualTolerance_)
     1792                nFlagged++;
     1793              break;
     1794            }
     1795            continue;
     1796          }
     1797          switch (getStatus(iSequence)) {
     1798
     1799          case basic:
     1800          case ClpSimplex::isFixed:
     1801            break;
     1802          case atUpperBound:
     1803            if (dj_[iSequence] > dualTolerance_)
     1804              nNormal++;
     1805            break;
     1806          case atLowerBound:
     1807            if (dj_[iSequence] < -dualTolerance_)
     1808              nNormal++;
     1809            break;
     1810          case isFree:
     1811          case superBasic:
     1812            if (fabs(dj_[iSequence]) > dualTolerance_)
     1813              nSuper++;
     1814            break;
     1815          }
     1816        }
     1817#ifdef CLP_DEBUG
     1818        if (handler_->logLevel() & 32)
     1819          printf("%d super, %d normal, %d flagged\n",
     1820            nSuper, nNormal, nFlagged);
     1821#endif
     1822
     1823        int nFlagged2 = 1;
     1824        if (lastSequenceIn < 0 && !nNormal && !nSuper) {
     1825          nFlagged2 = unflag();
     1826          if (pivotMode >= 10) {
     1827            pivotMode--;
     1828#ifdef CLP_DEBUG
     1829            if (handler_->logLevel() & 32)
     1830              printf("pivot mode now %d\n", pivotMode);
     1831#endif
     1832            if (pivotMode == 9)
     1833              pivotMode = 0; // switch off fast attempt
     1834          }
     1835        } else {
     1836          lastSequenceIn = -1;
     1837        }
     1838        if (!nFlagged2 || (normFlagged + normUnflagged < 1.0e-8)) {
     1839          primalColumnPivot_->setLooksOptimal(true);
     1840          return 0;
     1841        } else {
     1842          localPivotMode = -1;
     1843          nTotalPasses += nPasses;
     1844          nPasses = 0;
     1845          finished = false;
     1846          continue;
     1847        }
     1848      }
     1849      bestSequence = -1;
     1850      bestBasicSequence = -1;
     1851
     1852      // temp
     1853      nPasses++;
     1854      if (nPasses > 2000)
     1855        finished = true;
     1856      double theta = 1.0e30;
     1857      basicTheta = 1.0e30;
     1858      theta_ = 1.0e30;
     1859      double basicTolerance = 1.0e-4 * primalTolerance_;
     1860      for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     1861        //if (flagged(iSequence)
     1862        //  continue;
     1863        double alpha = dArray[iSequence];
     1864        Status thisStatus = getStatus(iSequence);
     1865        double oldValue = solution_[iSequence];
     1866        if (thisStatus != basic) {
     1867          if (fabs(alpha) >= acceptablePivot) {
     1868            if (alpha < 0.0) {
     1869              // variable going towards lower bound
     1870              double bound = lower_[iSequence];
     1871              oldValue -= bound;
     1872              if (oldValue + theta * alpha < 0.0) {
     1873                bestSequence = iSequence;
     1874                theta = CoinMax(0.0, oldValue / (-alpha));
     1875              }
     1876            } else {
     1877              // variable going towards upper bound
     1878              double bound = upper_[iSequence];
     1879              oldValue = bound - oldValue;
     1880              if (oldValue - theta * alpha < 0.0) {
     1881                bestSequence = iSequence;
     1882                theta = CoinMax(0.0, oldValue / alpha);
     1883              }
     1884            }
     1885          }
     1886        } else {
     1887          if (fabs(alpha) >= acceptableBasic) {
     1888            if (alpha < 0.0) {
     1889              // variable going towards lower bound
     1890              double bound = lower_[iSequence];
     1891              oldValue -= bound;
     1892              if (oldValue + basicTheta * alpha < -basicTolerance) {
     1893                bestBasicSequence = iSequence;
     1894                basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / (-alpha));
     1895              }
     1896            } else {
     1897              // variable going towards upper bound
     1898              double bound = upper_[iSequence];
     1899              oldValue = bound - oldValue;
     1900              if (oldValue - basicTheta * alpha < -basicTolerance) {
     1901                bestBasicSequence = iSequence;
     1902                basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / alpha);
     1903              }
     1904            }
     1905          }
     1906        }
     1907      }
     1908      theta_ = CoinMin(theta, basicTheta);
     1909      // Now find minimum of function
     1910      double objTheta2 = objective_->stepLength(this, solution_, dArray, CoinMin(theta, basicTheta),
     1911        currentObj, predictedObj, thetaObj);
     1912#ifdef CLP_DEBUG
     1913      if (handler_->logLevel() & 32)
     1914        printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
     1915#endif
     1916      objTheta2 = CoinMin(objTheta2, 1.0e29);
     1917#if MINTYPE == 1
     1918      if (conjugate) {
     1919        double offset;
     1920        const double *gradient = objective_->gradient(this,
     1921          dArray, offset,
     1922          true, 0);
     1923        double product = 0.0;
     1924        for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
     1925          double alpha = dArray[iSequence];
     1926          double value = alpha * gradient[iSequence];
     1927          product += value;
     1928        }
     1929        //#define INCLUDESLACK
     1930#ifdef INCLUDESLACK
     1931        for (; iSequence < numberColumns_ + numberRows_; iSequence++) {
     1932          double alpha = dArray[iSequence];
     1933          double value = alpha * cost_[iSequence];
     1934          product += value;
     1935        }
     1936#endif
     1937        if (product > 0.0)
     1938          objTheta = djNorm / product;
     1939        else
     1940          objTheta = COIN_DBL_MAX;
     1941#ifndef NDEBUG
     1942        if (product < -1.0e-8 && handler_->logLevel() > 1)
     1943          printf("bad product %g\n", product);
     1944#endif
     1945        product = CoinMax(product, 0.0);
     1946      } else {
     1947        objTheta = objTheta2;
     1948      }
     1949#else
     1950      objTheta = objTheta2;
     1951#endif
     1952      // if very small difference then take pivot (depends on djNorm?)
     1953      // distinguish between basic and non-basic
     1954      bool chooseObjTheta = objTheta < theta_;
     1955      if (chooseObjTheta) {
     1956        if (thetaObj < currentObj - 1.0e-12 && objTheta + 1.0e-10 > theta_)
     1957          chooseObjTheta = false;
     1958        //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
     1959        //chooseObjTheta=false;
     1960      }
     1961      //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
     1962      if (chooseObjTheta) {
     1963        theta_ = objTheta;
     1964      } else {
     1965        objTheta = CoinMax(objTheta, 1.00000001 * theta_ + 1.0e-12);
     1966        //if (theta+1.0e-13>basicTheta) {
     1967        //theta = CoinMax(theta,1.00000001*basicTheta);
     1968        //theta_ = basicTheta;
     1969        //}
     1970      }
     1971      // always out if one variable in and zero move
     1972      if (theta_ == basicTheta || (sequenceIn_ >= 0 && theta_ < 1.0e-10))
     1973        finished = true; // come out
     1974#ifdef CLP_DEBUG
     1975      if (handler_->logLevel() & 32) {
     1976        printf("%d pass,", nPasses);
     1977        if (sequenceIn_ >= 0)
     1978          printf(" Sin %d,", sequenceIn_);
     1979        if (basicTheta == theta_)
     1980          printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
     1981        else
     1982          printf(" basicTheta %g", basicTheta);
     1983        if (theta == theta_)
     1984          printf(" X(%d) non-basicTheta %g", bestSequence, theta);
     1985        else
     1986          printf(" non-basicTheta %g", theta);
     1987        printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
     1988        printf(" djNorm %g\n", djNorm);
     1989      }
     1990#endif
     1991      if (handler_->logLevel() > 3 && objTheta != theta_) {
     1992        printf("%d pass obj %g,", nPasses, currentObj);
     1993        if (sequenceIn_ >= 0)
     1994          printf(" Sin %d,", sequenceIn_);
     1995        if (basicTheta == theta_)
     1996          printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
     1997        else
     1998          printf(" basicTheta %g", basicTheta);
     1999        if (theta == theta_)
     2000          printf(" X(%d) non-basicTheta %g", bestSequence, theta);
     2001        else
     2002          printf(" non-basicTheta %g", theta);
     2003        printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
     2004        printf(" djNorm %g\n", djNorm);
     2005      }
     2006      if (objTheta != theta_) {
     2007        //printf("hit boundary after %d passes\n",nPasses);
     2008        nTotalPasses += nPasses;
     2009        nPasses = 0; // start again
     2010      }
     2011      if (localPivotMode < 10 || lastSequenceIn == bestSequence) {
     2012        if (theta_ == theta && theta_ < basicTheta && theta_ < 1.0e-5)
     2013          setFlagged(bestSequence); // out of active set
     2014      }
     2015      // Update solution
     2016      for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     2017        //for (iIndex=0;iIndex<number;iIndex++) {
     2018
     2019        //int iSequence = which[iIndex];
     2020        double alpha = dArray[iSequence];
     2021        if (alpha) {
     2022          double value = solution_[iSequence] + theta_ * alpha;
     2023          solution_[iSequence] = value;
     2024          switch (getStatus(iSequence)) {
     2025
     2026          case basic:
     2027          case isFixed:
     2028          case isFree:
     2029          case atUpperBound:
     2030          case atLowerBound:
     2031            nonLinearCost_->setOne(iSequence, value);
     2032            break;
     2033          case superBasic:
     2034            // To get correct action
     2035            setStatus(iSequence, isFixed);
     2036            nonLinearCost_->setOne(iSequence, value);
     2037            //assert (getStatus(iSequence)!=isFixed);
     2038            break;
     2039          }
     2040        }
     2041      }
     2042      if (objTheta < SMALLTHETA2 && objTheta == theta_) {
     2043        if (sequenceIn_ >= 0 && basicTheta < 1.0e-9) {
     2044          // try just one
     2045          localPivotMode = 0;
     2046          sequenceIn_ = -1;
     2047          nTotalPasses += nPasses;
     2048          nPasses = 0;
     2049          //finished=true;
     2050          //objTheta=0.0;
     2051          //theta_=0.0;
     2052        } else if (sequenceIn_ < 0 && nTotalPasses > 10) {
     2053          if (objTheta < 1.0e-10) {
     2054            finished = true;
     2055            //printf("zero move\n");
     2056            break;
     2057          }
     2058        }
     2059      }
     2060#ifdef CLP_DEBUG
     2061      if (handler_->logLevel() & 32) {
     2062        objective_->stepLength(this, solution_, work, 0.0,
     2063          currentObj, predictedObj, thetaObj);
     2064        printf("current obj %g after update - simple was %g\n", currentObj, simpleObjective);
     2065      }
     2066#endif
     2067      if (sequenceIn_ >= 0 && !finished && objTheta > 1.0e-4) {
     2068        // we made some progress - back to normal
     2069        if (localPivotMode < 10) {
     2070          localPivotMode = 0;
     2071          sequenceIn_ = -1;
     2072          nTotalPasses += nPasses;
     2073          nPasses = 0;
     2074        }
     2075#ifdef CLP_DEBUG
     2076        if (handler_->logLevel() & 32)
     2077          printf("some progress?\n");
     2078#endif
     2079      }
     2080#if CLP_DEBUG > 1
     2081      longArray->checkClean();
     2082#endif
     2083    }
     2084#ifdef CLP_DEBUG
     2085    if (handler_->logLevel() & 32)
     2086      printf("out of loop after %d (%d) passes\n", nPasses, nTotalPasses);
     2087#endif
     2088    if (nTotalPasses >= 1000 || (nTotalPasses > 10 && sequenceIn_ < 0 && theta_ < 1.0e-10))
     2089      returnCode = 2;
     2090    bool ordinaryDj = false;
     2091    //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
     2092    //printf("could try ordinary iteration %g\n",theta_);
     2093    if (sequenceIn_ >= 0 && numberNonBasic == 1 && theta_ < 1.0e-15) {
     2094      //printf("trying ordinary iteration\n");
     2095      ordinaryDj = true;
     2096    }
     2097    if (!basicTheta && !ordinaryDj) {
     2098      //returnCode=2;
     2099      //objTheta=-1.0; // so we fall through
     2100    }
     2101    assert(theta_ < 1.0e30); // for now
     2102    // See if we need to pivot
     2103    if (theta_ == basicTheta || ordinaryDj) {
     2104      if (!ordinaryDj) {
     2105        // find an incoming column which will force pivot
     2106        int iRow;
     2107        pivotRow_ = -1;
     2108        for (iRow = 0; iRow < numberRows_; iRow++) {
     2109          if (pivotVariable_[iRow] == bestBasicSequence) {
     2110            pivotRow_ = iRow;
     2111            break;
     2112          }
     2113        }
     2114        assert(pivotRow_ >= 0);
     2115        // Get good size for pivot
     2116        double acceptablePivot = 1.0e-7;
     2117        if (factorization_->pivots() > 10)
     2118          acceptablePivot = 1.0e-5; // if we have iterated be more strict
     2119        else if (factorization_->pivots() > 5)
     2120          acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
     2121        // should be dArray but seems better this way!
     2122        double direction = work[bestBasicSequence] > 0.0 ? -1.0 : 1.0;
     2123        // create as packed
     2124        rowArray->createPacked(1, &pivotRow_, &direction);
     2125        factorization_->updateColumnTranspose(spare, rowArray);
     2126        // put row of tableau in rowArray and columnArray
     2127        matrix_->transposeTimes(this, -1.0,
     2128          rowArray, spare, columnArray);
     2129        // choose one futhest away from bound which has reasonable pivot
     2130        // If increasing we want negative alpha
     2131
     2132        double *work2;
     2133        int iSection;
     2134
     2135        sequenceIn_ = -1;
     2136        double bestValue = -1.0;
     2137        double bestDirection = 0.0;
     2138        // First pass we take correct direction and large pivots
     2139        // then correct direction
     2140        // then any
     2141        double check[] = { 1.0e-8, -1.0e-12, -1.0e30 };
     2142        double mult[] = { 100.0, 1.0, 1.0 };
     2143        for (int iPass = 0; iPass < 3; iPass++) {
     2144          //if (!bestValue&&iPass==2)
     2145          //bestValue=-1.0;
     2146          double acceptable = acceptablePivot * mult[iPass];
     2147          double checkValue = check[iPass];
     2148          for (iSection = 0; iSection < 2; iSection++) {
     2149
     2150            int addSequence;
     2151
     2152            if (!iSection) {
     2153              work2 = rowArray->denseVector();
     2154              number = rowArray->getNumElements();
     2155              which = rowArray->getIndices();
     2156              addSequence = numberColumns_;
     2157            } else {
     2158              work2 = columnArray->denseVector();
     2159              number = columnArray->getNumElements();
     2160              which = columnArray->getIndices();
     2161              addSequence = 0;
     2162            }
     2163            int i;
     2164
     2165            for (i = 0; i < number; i++) {
     2166              int iSequence = which[i] + addSequence;
     2167              if (flagged(iSequence))
     2168                continue;
     2169              //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
     2170              //                  upper_[iSequence]-solution_[iSequence]);
     2171              double alpha = work2[i];
     2172              // should be dArray but seems better this way!
     2173              double change = work[iSequence];
     2174              Status thisStatus = getStatus(iSequence);
     2175              double direction = 0;
     2176              ;
     2177              switch (thisStatus) {
     2178
     2179              case basic:
     2180              case ClpSimplex::isFixed:
     2181                break;
     2182              case isFree:
     2183              case superBasic:
     2184                if (alpha < -acceptable && change > checkValue)
     2185                  direction = 1.0;
     2186                else if (alpha > acceptable && change < -checkValue)
     2187                  direction = -1.0;
     2188                break;
     2189              case atUpperBound:
     2190                if (alpha > acceptable && change < -checkValue)
     2191                  direction = -1.0;
     2192                else if (iPass == 2 && alpha < -acceptable && change < -checkValue)
     2193                  direction = 1.0;
     2194                break;
     2195              case atLowerBound:
     2196                if (alpha < -acceptable && change > checkValue)
     2197                  direction = 1.0;
     2198                else if (iPass == 2 && alpha > acceptable && change > checkValue)
     2199                  direction = -1.0;
     2200                break;
     2201              }
     2202              if (direction) {
     2203                if (sequenceIn_ != lastSequenceIn || localPivotMode < 10) {
     2204                  if (CoinMin(solution_[iSequence] - lower_[iSequence],
     2205                        upper_[iSequence] - solution_[iSequence])
     2206                    > bestValue) {
     2207                    bestValue = CoinMin(solution_[iSequence] - lower_[iSequence],
     2208                      upper_[iSequence] - solution_[iSequence]);
     2209                    sequenceIn_ = iSequence;
     2210                    bestDirection = direction;
     2211                  }
     2212                } else {
     2213                  // choose
     2214                  bestValue = COIN_DBL_MAX;
     2215                  sequenceIn_ = iSequence;
     2216                  bestDirection = direction;
     2217                }
     2218              }
     2219            }
     2220          }
     2221          if (sequenceIn_ >= 0 && bestValue > 0.0)
     2222            break;
     2223        }
     2224        if (sequenceIn_ >= 0) {
     2225          valueIn_ = solution_[sequenceIn_];
     2226          dualIn_ = dj_[sequenceIn_];
     2227          if (bestDirection < 0.0) {
     2228            // we want positive dj
     2229            if (dualIn_ <= 0.0) {
     2230              //printf("bad dj - xx %g\n",dualIn_);
     2231              dualIn_ = 1.0;
     2232            }
    23342233          } else {
    2335                // round again
    2336                longArray->clear();
    2337                if (djNorm < 1.0e-3 && !localPivotMode) {
    2338                     if (djNorm == 1.2345e-20 && djNorm0 > 1.0e-4) {
     2234            // we want negative dj
     2235            if (dualIn_ >= 0.0) {
     2236              //printf("bad dj - xx %g\n",dualIn_);
     2237              dualIn_ = -1.0;
     2238            }
     2239          }
     2240          lowerIn_ = lower_[sequenceIn_];
     2241          upperIn_ = upper_[sequenceIn_];
     2242          if (dualIn_ > 0.0)
     2243            directionIn_ = -1;
     2244          else
     2245            directionIn_ = 1;
     2246        } else {
     2247          //ordinaryDj=true;
    23392248#ifdef CLP_DEBUG
    2340                          if (handler_->logLevel() & 32)
    2341                               printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n", djNorm0, nPasses,
    2342                                      localPivotMode, returnCode);
    2343 #endif
    2344                          //if (!localPivotMode)
    2345                          //returnCode=2; // force singleton
    2346                     } else {
     2249          if (handler_->logLevel() & 32) {
     2250            printf("no easy pivot - norm %g mode %d\n", djNorm, localPivotMode);
     2251            if (rowArray->getNumElements() + columnArray->getNumElements() < 12) {
     2252              for (iSection = 0; iSection < 2; iSection++) {
     2253
     2254                int addSequence;
     2255
     2256                if (!iSection) {
     2257                  work2 = rowArray->denseVector();
     2258                  number = rowArray->getNumElements();
     2259                  which = rowArray->getIndices();
     2260                  addSequence = numberColumns_;
     2261                } else {
     2262                  work2 = columnArray->denseVector();
     2263                  number = columnArray->getNumElements();
     2264                  which = columnArray->getIndices();
     2265                  addSequence = 0;
     2266                }
     2267                int i;
     2268                char section[] = { 'R', 'C' };
     2269                for (i = 0; i < number; i++) {
     2270                  int iSequence = which[i] + addSequence;
     2271                  if (flagged(iSequence)) {
     2272                    printf("%c%d flagged\n", section[iSection], which[i]);
     2273                    continue;
     2274                  } else {
     2275                    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
     2276                      section[iSection], which[i], status_[iSequence],
     2277                      lower_[iSequence], solution_[iSequence], upper_[iSequence],
     2278                      work2[i], work[iSequence]);
     2279                  }
     2280                }
     2281              }
     2282            }
     2283          }
     2284#endif
     2285          assert(sequenceIn_ < 0);
     2286          justOne = true;
     2287          allFinished = false; // Round again
     2288          finished = false;
     2289          nTotalPasses += nPasses;
     2290          nPasses = 0;
     2291          if (djNorm < 0.9 * djNorm0 && djNorm < 1.0e-3 && !localPivotMode) {
    23472292#ifdef CLP_DEBUG
    2348                          if (handler_->logLevel() & 32)
    2349                               printf("unflagging as djNorm %g %g, %d passes\n", djNorm, djNorm0, nPasses);
    2350 #endif
    2351                          if (pivotMode >= 10) {
    2352                               pivotMode--;
    2353 #ifdef CLP_DEBUG
    2354                               if (handler_->logLevel() & 32)
    2355                                    printf("pivot mode now %d\n", pivotMode);
    2356 #endif
    2357                               if (pivotMode == 9)
    2358                                    pivotMode = 0;       // switch off fast attempt
    2359                          }
    2360                          bool unflagged = unflag() != 0;
    2361                          if (!unflagged && djNorm < 1.0e-10) {
    2362                               // ? declare victory
    2363                               sequenceIn_ = -1;
    2364                               returnCode = 0;
    2365                          }
    2366                     }
    2367                }
    2368           }
    2369      }
    2370      if (djNorm0 < 1.0e-12 * normFlagged) {
     2293            if (handler_->logLevel() & 32)
     2294              printf("no pivot - mode %d norms %g %g - unflagging\n",
     2295                localPivotMode, djNorm0, djNorm);
     2296#endif
     2297            unflag(); //unflagging
     2298            returnCode = 1;
     2299          } else {
     2300            returnCode = 2; // do single incoming
     2301            returnCode = 1;
     2302          }
     2303        }
     2304      }
     2305      rowArray->clear();
     2306      columnArray->clear();
     2307      longArray->clear();
     2308      if (ordinaryDj) {
     2309        valueIn_ = solution_[sequenceIn_];
     2310        dualIn_ = dj_[sequenceIn_];
     2311        lowerIn_ = lower_[sequenceIn_];
     2312        upperIn_ = upper_[sequenceIn_];
     2313        if (dualIn_ > 0.0)
     2314          directionIn_ = -1;
     2315        else
     2316          directionIn_ = 1;
     2317      }
     2318      if (returnCode == 1)
     2319        returnCode = 0;
     2320    } else {
     2321      // round again
     2322      longArray->clear();
     2323      if (djNorm < 1.0e-3 && !localPivotMode) {
     2324        if (djNorm == 1.2345e-20 && djNorm0 > 1.0e-4) {
    23712325#ifdef CLP_DEBUG
    23722326          if (handler_->logLevel() & 32)
    2373                printf("unflagging as djNorm %g %g and flagged norm %g\n", djNorm, djNorm0, normFlagged);
    2374 #endif
    2375           unflag();
    2376      }
    2377      if (saveObj - currentObj < 1.0e-5 && nTotalPasses > 2000) {
    2378           normUnflagged = 0.0;
    2379           double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
    2380           for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
    2381                switch(getStatus(iSequence)) {
    2382 
    2383                case basic:
    2384                case ClpSimplex::isFixed:
    2385                     break;
    2386                case atUpperBound:
    2387                     if (dj_[iSequence] > dualTolerance3)
    2388                          normFlagged += dj_[iSequence] * dj_[iSequence];
    2389                     break;
    2390                case atLowerBound:
    2391                     if (dj_[iSequence] < -dualTolerance3)
    2392                          normFlagged += dj_[iSequence] * dj_[iSequence];
    2393                     break;
    2394                case isFree:
    2395                case superBasic:
    2396                     if (fabs(dj_[iSequence]) > dualTolerance3)
    2397                          normFlagged += dj_[iSequence] * dj_[iSequence];
    2398                     break;
    2399                }
    2400           }
    2401           if (handler_->logLevel() > 2)
    2402                printf("possible optimal  %d %d %g %g\n",
    2403                       nBigPasses, nTotalPasses, saveObj - currentObj, normFlagged);
    2404           if (normFlagged < 1.0e-5) {
    2405                unflag();
    2406                primalColumnPivot_->setLooksOptimal(true);
    2407                returnCode = 0;
    2408           }
    2409      }
    2410      return returnCode;
     2327            printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n", djNorm0, nPasses,
     2328              localPivotMode, returnCode);
     2329#endif
     2330          //if (!localPivotMode)
     2331          //returnCode=2; // force singleton
     2332        } else {
     2333#ifdef CLP_DEBUG
     2334          if (handler_->logLevel() & 32)
     2335            printf("unflagging as djNorm %g %g, %d passes\n", djNorm, djNorm0, nPasses);
     2336#endif
     2337          if (pivotMode >= 10) {
     2338            pivotMode--;
     2339#ifdef CLP_DEBUG
     2340            if (handler_->logLevel() & 32)
     2341              printf("pivot mode now %d\n", pivotMode);
     2342#endif
     2343            if (pivotMode == 9)
     2344              pivotMode = 0; // switch off fast attempt
     2345          }
     2346          bool unflagged = unflag() != 0;
     2347          if (!unflagged && djNorm < 1.0e-10) {
     2348            // ? declare victory
     2349            sequenceIn_ = -1;
     2350            returnCode = 0;
     2351          }
     2352        }
     2353      }
     2354    }
     2355  }
     2356  if (djNorm0 < 1.0e-12 * normFlagged) {
     2357#ifdef CLP_DEBUG
     2358    if (handler_->logLevel() & 32)
     2359      printf("unflagging as djNorm %g %g and flagged norm %g\n", djNorm, djNorm0, normFlagged);
     2360#endif
     2361    unflag();
     2362  }
     2363  if (saveObj - currentObj < 1.0e-5 && nTotalPasses > 2000) {
     2364    normUnflagged = 0.0;
     2365    double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
     2366    for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
     2367      switch (getStatus(iSequence)) {
     2368
     2369      case basic:
     2370      case ClpSimplex::isFixed:
     2371        break;
     2372      case atUpperBound:
     2373        if (dj_[iSequence] > dualTolerance3)
     2374          normFlagged += dj_[iSequence] * dj_[iSequence];
     2375        break;
     2376      case atLowerBound:
     2377        if (dj_[iSequence] < -dualTolerance3)
     2378          normFlagged += dj_[iSequence] * dj_[iSequence];
     2379        break;
     2380      case isFree:
     2381      case superBasic:
     2382        if (fabs(dj_[iSequence]) > dualTolerance3)
     2383          normFlagged += dj_[iSequence] * dj_[iSequence];
     2384        break;
     2385      }
     2386    }
     2387    if (handler_->logLevel() > 2)
     2388      printf("possible optimal  %d %d %g %g\n",
     2389        nBigPasses, nTotalPasses, saveObj - currentObj, normFlagged);
     2390    if (normFlagged < 1.0e-5) {
     2391      unflag();
     2392      primalColumnPivot_->setLooksOptimal(true);
     2393      returnCode = 0;
     2394    }
     2395  }
     2396  return returnCode;
    24112397}
    24122398/* Do last half of an iteration.
     
    24222408
    24232409*/
    2424 int
    2425 ClpSimplexNonlinear::pivotNonlinearResult()
     2410int ClpSimplexNonlinear::pivotNonlinearResult()
    24262411{
    24272412
    2428      int returnCode = -1;
    2429 
    2430      rowArray_[1]->clear();
    2431 
    2432      // we found a pivot column
    2433      // update the incoming column
    2434      unpackPacked(rowArray_[1]);
    2435      factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
    2436      theta_ = 0.0;
    2437      double * work = rowArray_[1]->denseVector();
    2438      int number = rowArray_[1]->getNumElements();
    2439      int * which = rowArray_[1]->getIndices();
    2440      bool keepValue = false;
    2441      double saveValue = 0.0;
    2442      if (pivotRow_ >= 0) {
    2443           sequenceOut_ = pivotVariable_[pivotRow_];;
    2444           valueOut_ = solution(sequenceOut_);
    2445           keepValue = true;
    2446           saveValue = valueOut_;
    2447           lowerOut_ = lower_[sequenceOut_];
    2448           upperOut_ = upper_[sequenceOut_];
    2449           int iIndex;
    2450           for (iIndex = 0; iIndex < number; iIndex++) {
    2451 
    2452                int iRow = which[iIndex];
    2453                if (iRow == pivotRow_) {
    2454                     alpha_ = work[iIndex];
    2455                     break;
    2456                }
    2457           }
    2458      } else {
    2459           int iIndex;
    2460           double smallest = COIN_DBL_MAX;
    2461           for (iIndex = 0; iIndex < number; iIndex++) {
    2462                int iRow = which[iIndex];
    2463                double alpha = work[iIndex];
    2464                if (fabs(alpha) > 1.0e-6) {
    2465                     int iPivot = pivotVariable_[iRow];
    2466                     double distance = CoinMin(upper_[iPivot] - solution_[iPivot],
    2467                                               solution_[iPivot] - lower_[iPivot]);
    2468                     if (distance < smallest) {
    2469                          pivotRow_ = iRow;
    2470                          alpha_ = alpha;
    2471                          smallest = distance;
    2472                     }
    2473                }
    2474           }
    2475           if (smallest > primalTolerance_) {
    2476                smallest = COIN_DBL_MAX;
    2477                for (iIndex = 0; iIndex < number; iIndex++) {
    2478                     int iRow = which[iIndex];
    2479                     double alpha = work[iIndex];
    2480                     if (fabs(alpha) > 1.0e-6) {
    2481                          double distance = randomNumberGenerator_.randomDouble();
    2482                          if (distance < smallest) {
    2483                               pivotRow_ = iRow;
    2484                               alpha_ = alpha;
    2485                               smallest = distance;
    2486                          }
    2487                     }
    2488                }
    2489           }
    2490           assert (pivotRow_ >= 0);
    2491           sequenceOut_ = pivotVariable_[pivotRow_];;
    2492           valueOut_ = solution(sequenceOut_);
    2493           lowerOut_ = lower_[sequenceOut_];
    2494           upperOut_ = upper_[sequenceOut_];
    2495      }
    2496      double newValue = valueOut_ - theta_ * alpha_;
    2497      bool isSuperBasic = false;
    2498      if (valueOut_ >= upperOut_ - primalTolerance_) {
    2499           directionOut_ = -1;    // to upper bound
    2500           upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
    2501           upperOut_ = newValue;
    2502      } else if (valueOut_ <= lowerOut_ + primalTolerance_) {
    2503           directionOut_ = 1;    // to lower bound
    2504           lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
    2505      } else {
    2506           lowerOut_ = valueOut_;
    2507           upperOut_ = valueOut_;
    2508           isSuperBasic = true;
    2509           //printf("XX superbasic out\n");
    2510      }
    2511      dualOut_ = dj_[sequenceOut_];
    2512      // if stable replace in basis
    2513 
    2514      int updateStatus = factorization_->replaceColumn(this,
    2515                         rowArray_[2],
    2516                         rowArray_[1],
    2517                         pivotRow_,
    2518                         alpha_);
    2519 
    2520      // if no pivots, bad update but reasonable alpha - take and invert
    2521      if (updateStatus == 2 &&
    2522                lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
    2523           updateStatus = 4;
    2524      if (updateStatus == 1 || updateStatus == 4) {
    2525           // slight error
    2526           if (factorization_->pivots() > 5 || updateStatus == 4) {
    2527                returnCode = -3;
    2528           }
    2529      } else if (updateStatus == 2) {
    2530           // major error
    2531           // better to have small tolerance even if slower
    2532           factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
    2533           int maxFactor = factorization_->maximumPivots();
    2534           if (maxFactor > 10) {
    2535                if (forceFactorization_ < 0)
    2536                     forceFactorization_ = maxFactor;
    2537                forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
    2538           }
    2539           // later we may need to unwind more e.g. fake bounds
    2540           if(lastGoodIteration_ != numberIterations_) {
    2541                clearAll();
    2542                pivotRow_ = -1;
    2543                returnCode = -4;
    2544           } else {
    2545                // need to reject something
    2546                char x = isColumn(sequenceIn_) ? 'C' : 'R';
    2547                handler_->message(CLP_SIMPLEX_FLAG, messages_)
    2548                          << x << sequenceWithin(sequenceIn_)
    2549                          << CoinMessageEol;
    2550                setFlagged(sequenceIn_);
    2551                progress_.clearBadTimes();
    2552                lastBadIteration_ = numberIterations_; // say be more cautious
    2553                clearAll();
    2554                pivotRow_ = -1;
    2555                sequenceOut_ = -1;
    2556                returnCode = -5;
    2557 
    2558           }
    2559           return returnCode;
    2560      } else if (updateStatus == 3) {
    2561           // out of memory
    2562           // increase space if not many iterations
    2563           if (factorization_->pivots() <
    2564                     0.5 * factorization_->maximumPivots() &&
    2565                     factorization_->pivots() < 200)
    2566                factorization_->areaFactor(
    2567                     factorization_->areaFactor() * 1.1);
    2568           returnCode = -2; // factorize now
    2569      } else if (updateStatus == 5) {
    2570           problemStatus_ = -2; // factorize now
    2571      }
    2572 
    2573      // update primal solution
    2574 
    2575      double objectiveChange = 0.0;
    2576      // after this rowArray_[1] is not empty - used to update djs
    2577      // If pivot row >= numberRows then may be gub
    2578      updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, 1);
    2579 
    2580      double oldValue = valueIn_;
    2581      if (directionIn_ == -1) {
    2582           // as if from upper bound
    2583           if (sequenceIn_ != sequenceOut_) {
    2584                // variable becoming basic
    2585                valueIn_ -= fabs(theta_);
    2586           } else {
    2587                valueIn_ = lowerIn_;
    2588           }
    2589      } else {
    2590           // as if from lower bound
    2591           if (sequenceIn_ != sequenceOut_) {
    2592                // variable becoming basic
    2593                valueIn_ += fabs(theta_);
    2594           } else {
    2595                valueIn_ = upperIn_;
    2596           }
    2597      }
    2598      objectiveChange += dualIn_ * (valueIn_ - oldValue);
    2599      // outgoing
    2600      if (sequenceIn_ != sequenceOut_) {
    2601           if (directionOut_ > 0) {
    2602                valueOut_ = lowerOut_;
    2603           } else {
    2604                valueOut_ = upperOut_;
    2605           }
    2606           if(valueOut_ < lower_[sequenceOut_] - primalTolerance_)
    2607                valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
    2608           else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
    2609                valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
    2610           // may not be exactly at bound and bounds may have changed
    2611           // Make sure outgoing looks feasible
    2612           if (!isSuperBasic)
    2613                directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
    2614           solution_[sequenceOut_] = valueOut_;
    2615      }
    2616      // change cost and bounds on incoming if primal
    2617      nonLinearCost_->setOne(sequenceIn_, valueIn_);
    2618      int whatNext = housekeeping(objectiveChange);
    2619      if (keepValue)
    2620           solution_[sequenceOut_] = saveValue;
    2621      if (isSuperBasic)
    2622           setStatus(sequenceOut_, superBasic);
    2623      //#define CLP_DEBUG
     2413  int returnCode = -1;
     2414
     2415  rowArray_[1]->clear();
     2416
     2417  // we found a pivot column
     2418  // update the incoming column
     2419  unpackPacked(rowArray_[1]);
     2420  factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
     2421  theta_ = 0.0;
     2422  double *work = rowArray_[1]->denseVector();
     2423  int number = rowArray_[1]->getNumElements();
     2424  int *which = rowArray_[1]->getIndices();
     2425  bool keepValue = false;
     2426  double saveValue = 0.0;
     2427  if (pivotRow_ >= 0) {
     2428    sequenceOut_ = pivotVariable_[pivotRow_];
     2429    ;
     2430    valueOut_ = solution(sequenceOut_);
     2431    keepValue = true;
     2432    saveValue = valueOut_;
     2433    lowerOut_ = lower_[sequenceOut_];
     2434    upperOut_ = upper_[sequenceOut_];
     2435    int iIndex;
     2436    for (iIndex = 0; iIndex < number; iIndex++) {
     2437
     2438      int iRow = which[iIndex];
     2439      if (iRow == pivotRow_) {
     2440        alpha_ = work[iIndex];
     2441        break;
     2442      }
     2443    }
     2444  } else {
     2445    int iIndex;
     2446    double smallest = COIN_DBL_MAX;
     2447    for (iIndex = 0; iIndex < number; iIndex++) {
     2448      int iRow = which[iIndex];
     2449      double alpha = work[iIndex];
     2450      if (fabs(alpha) > 1.0e-6) {
     2451        int iPivot = pivotVariable_[iRow];
     2452        double distance = CoinMin(upper_[iPivot] - solution_[iPivot],
     2453          solution_[iPivot] - lower_[iPivot]);
     2454        if (distance < smallest) {
     2455          pivotRow_ = iRow;
     2456          alpha_ = alpha;
     2457          smallest = distance;
     2458        }
     2459      }
     2460    }
     2461    if (smallest > primalTolerance_) {
     2462      smallest = COIN_DBL_MAX;
     2463      for (iIndex = 0; iIndex < number; iIndex++) {
     2464        int iRow = which[iIndex];
     2465        double alpha = work[iIndex];
     2466        if (fabs(alpha) > 1.0e-6) {
     2467          double distance = randomNumberGenerator_.randomDouble();
     2468          if (distance < smallest) {
     2469            pivotRow_ = iRow;
     2470            alpha_ = alpha;
     2471            smallest = distance;
     2472          }
     2473        }
     2474      }
     2475    }
     2476    assert(pivotRow_ >= 0);
     2477    sequenceOut_ = pivotVariable_[pivotRow_];
     2478    ;
     2479    valueOut_ = solution(sequenceOut_);
     2480    lowerOut_ = lower_[sequenceOut_];
     2481    upperOut_ = upper_[sequenceOut_];
     2482  }
     2483  double newValue = valueOut_ - theta_ * alpha_;
     2484  bool isSuperBasic = false;
     2485  if (valueOut_ >= upperOut_ - primalTolerance_) {
     2486    directionOut_ = -1; // to upper bound
     2487    upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
     2488    upperOut_ = newValue;
     2489  } else if (valueOut_ <= lowerOut_ + primalTolerance_) {
     2490    directionOut_ = 1; // to lower bound
     2491    lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
     2492  } else {
     2493    lowerOut_ = valueOut_;
     2494    upperOut_ = valueOut_;
     2495    isSuperBasic = true;
     2496    //printf("XX superbasic out\n");
     2497  }
     2498  dualOut_ = dj_[sequenceOut_];
     2499  // if stable replace in basis
     2500
     2501  int updateStatus = factorization_->replaceColumn(this,
     2502    rowArray_[2],
     2503    rowArray_[1],
     2504    pivotRow_,
     2505    alpha_);
     2506
     2507  // if no pivots, bad update but reasonable alpha - take and invert
     2508  if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
     2509    updateStatus = 4;
     2510  if (updateStatus == 1 || updateStatus == 4) {
     2511    // slight error
     2512    if (factorization_->pivots() > 5 || updateStatus == 4) {
     2513      returnCode = -3;
     2514    }
     2515  } else if (updateStatus == 2) {
     2516    // major error
     2517    // better to have small tolerance even if slower
     2518    factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
     2519    int maxFactor = factorization_->maximumPivots();
     2520    if (maxFactor > 10) {
     2521      if (forceFactorization_ < 0)
     2522        forceFactorization_ = maxFactor;
     2523      forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
     2524    }
     2525    // later we may need to unwind more e.g. fake bounds
     2526    if (lastGoodIteration_ != numberIterations_) {
     2527      clearAll();
     2528      pivotRow_ = -1;
     2529      returnCode = -4;
     2530    } else {
     2531      // need to reject something
     2532      char x = isColumn(sequenceIn_) ? 'C' : 'R';
     2533      handler_->message(CLP_SIMPLEX_FLAG, messages_)
     2534        << x << sequenceWithin(sequenceIn_)
     2535        << CoinMessageEol;
     2536      setFlagged(sequenceIn_);
     2537      progress_.clearBadTimes();
     2538      lastBadIteration_ = numberIterations_; // say be more cautious
     2539      clearAll();
     2540      pivotRow_ = -1;
     2541      sequenceOut_ = -1;
     2542      returnCode = -5;
     2543    }
     2544    return returnCode;
     2545  } else if (updateStatus == 3) {
     2546    // out of memory
     2547    // increase space if not many iterations
     2548    if (factorization_->pivots() < 0.5 * factorization_->maximumPivots() && factorization_->pivots() < 200)
     2549      factorization_->areaFactor(
     2550        factorization_->areaFactor() * 1.1);
     2551    returnCode = -2; // factorize now
     2552  } else if (updateStatus == 5) {
     2553    problemStatus_ = -2; // factorize now
     2554  }
     2555
     2556  // update primal solution
     2557
     2558  double objectiveChange = 0.0;
     2559  // after this rowArray_[1] is not empty - used to update djs
     2560  // If pivot row >= numberRows then may be gub
     2561  updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, 1);
     2562
     2563  double oldValue = valueIn_;
     2564  if (directionIn_ == -1) {
     2565    // as if from upper bound
     2566    if (sequenceIn_ != sequenceOut_) {
     2567      // variable becoming basic
     2568      valueIn_ -= fabs(theta_);
     2569    } else {
     2570      valueIn_ = lowerIn_;
     2571    }
     2572  } else {
     2573    // as if from lower bound
     2574    if (sequenceIn_ != sequenceOut_) {
     2575      // variable becoming basic
     2576      valueIn_ += fabs(theta_);
     2577    } else {
     2578      valueIn_ = upperIn_;
     2579    }
     2580  }
     2581  objectiveChange += dualIn_ * (valueIn_ - oldValue);
     2582  // outgoing
     2583  if (sequenceIn_ != sequenceOut_) {
     2584    if (directionOut_ > 0) {
     2585      valueOut_ = lowerOut_;
     2586    } else {
     2587      valueOut_ = upperOut_;
     2588    }
     2589    if (valueOut_ < lower_[sequenceOut_] - primalTolerance_)
     2590      valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
     2591    else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
     2592      valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
     2593    // may not be exactly at bound and bounds may have changed
     2594    // Make sure outgoing looks feasible
     2595    if (!isSuperBasic)
     2596      directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
     2597    solution_[sequenceOut_] = valueOut_;
     2598  }
     2599  // change cost and bounds on incoming if primal
     2600  nonLinearCost_->setOne(sequenceIn_, valueIn_);
     2601  int whatNext = housekeeping(objectiveChange);
     2602  if (keepValue)
     2603    solution_[sequenceOut_] = saveValue;
     2604  if (isSuperBasic)
     2605    setStatus(sequenceOut_, superBasic);
     2606    //#define CLP_DEBUG
    26242607#if CLP_DEBUG > 1
    2625      {
    2626           int ninf = matrix_->checkFeasible(this);
    2627           if (ninf)
    2628                printf("infeas %d\n", ninf);
    2629      }
    2630 #endif
    2631      if (whatNext == 1) {
    2632           returnCode = -2; // refactorize
    2633      } else if (whatNext == 2) {
    2634           // maximum iterations or equivalent
    2635           returnCode = 3;
    2636      } else if(numberIterations_ == lastGoodIteration_
    2637                + 2 * factorization_->maximumPivots()) {
    2638           // done a lot of flips - be safe
    2639           returnCode = -2; // refactorize
    2640      }
    2641      // Check event
    2642      {
    2643           int status = eventHandler_->event(ClpEventHandler::endOfIteration);
    2644           if (status >= 0) {
    2645                problemStatus_ = 5;
    2646                secondaryStatus_ = ClpEventHandler::endOfIteration;
    2647                returnCode = 4;
    2648           }
    2649      }
    2650      return returnCode;
     2608  {
     2609    int ninf = matrix_->checkFeasible(this);
     2610    if (ninf)
     2611      printf("infeas %d\n", ninf);
     2612  }
     2613#endif
     2614  if (whatNext == 1) {
     2615    returnCode = -2; // refactorize
     2616  } else if (whatNext == 2) {
     2617    // maximum iterations or equivalent
     2618    returnCode = 3;
     2619  } else if (numberIterations_ == lastGoodIteration_ + 2 * factorization_->maximumPivots()) {
     2620    // done a lot of flips - be safe
     2621    returnCode = -2; // refactorize
     2622  }
     2623  // Check event
     2624  {
     2625    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
     2626    if (status >= 0) {
     2627      problemStatus_ = 5;
     2628      secondaryStatus_ = ClpEventHandler::endOfIteration;
     2629      returnCode = 4;
     2630    }
     2631  }
     2632  return returnCode;
    26512633}
    26522634// May use a cut approach for solving any LP
    2653 int
    2654 ClpSimplexNonlinear::primalDualCuts(char * rowsIn, int startUp,int algorithm)
     2635int ClpSimplexNonlinear::primalDualCuts(char *rowsIn, int startUp, int algorithm)
    26552636{
    26562637  if (!rowsIn) {
    2657     if (algorithm>0)
     2638    if (algorithm > 0)
    26582639      return ClpSimplex::primal(startUp);
    26592640    else
     
    26612642      return ClpSimplex::dual(startUp);
    26622643  } else {
    2663     int numberUsed=0;
    2664     int rowsThreshold=CoinMax(100,numberRows_/2);
     2644    int numberUsed = 0;
     2645    int rowsThreshold = CoinMax(100, numberRows_ / 2);
    26652646    //int rowsTry=CoinMax(50,numberRows_/4);
    26662647    // Just add this number of rows each time in small problem
     
    26682649    smallNumberRows = CoinMin(smallNumberRows, numberRows_ / 20);
    26692650    // We will need arrays to choose rows to add
    2670     double * weight = new double [numberRows_];
    2671     int * sort = new int [numberRows_+numberColumns_];
    2672     int * whichColumns = sort+numberRows_;
     2651    double *weight = new double[numberRows_];
     2652    int *sort = new int[numberRows_ + numberColumns_];
     2653    int *whichColumns = sort + numberRows_;
    26732654    int numberSort = 0;
    2674     for (int i=0;i<numberRows_;i++) {
     2655    for (int i = 0; i < numberRows_; i++) {
    26752656      if (rowsIn[i])
    2676         numberUsed++;
     2657        numberUsed++;
    26772658    }
    26782659    if (numberUsed) {
     
    26802661    } else {
    26812662      // If non slack use that information
    2682       int numberBinding=0;
    2683       numberPrimalInfeasibilities_=0;
    2684       sumPrimalInfeasibilities_=0.0;
     2663      int numberBinding = 0;
     2664      numberPrimalInfeasibilities_ = 0;
     2665      sumPrimalInfeasibilities_ = 0.0;
    26852666      memset(rowActivity_, 0, numberRows_ * sizeof(double));
    26862667      times(1.0, columnActivity_, rowActivity_);
    2687       for (int i=0;i<numberRows_;i++) {
    2688         double lowerDifference = rowActivity_[i]-rowLower_[i];
    2689         double upperDifference = rowActivity_[i]-rowUpper_[i];
    2690         if (lowerDifference<-10*primalTolerance_||upperDifference>10.0*primalTolerance_) {
    2691           numberPrimalInfeasibilities_++;
    2692           if (lowerDifference<0.0)
    2693             sumPrimalInfeasibilities_ -= lowerDifference;
    2694           else
    2695             sumPrimalInfeasibilities_ += upperDifference;
    2696           rowsIn[i]=1;
    2697         } else if (getRowStatus(i)!=basic) {
    2698           numberBinding++;
    2699           rowsIn[i]=1;
    2700         }
    2701       }
    2702       if (numberBinding<rowsThreshold) {
    2703         // try
     2668      for (int i = 0; i < numberRows_; i++) {
     2669        double lowerDifference = rowActivity_[i] - rowLower_[i];
     2670        double upperDifference = rowActivity_[i] - rowUpper_[i];
     2671        if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
     2672          numberPrimalInfeasibilities_++;
     2673          if (lowerDifference < 0.0)
     2674            sumPrimalInfeasibilities_ -= lowerDifference;
     2675          else
     2676            sumPrimalInfeasibilities_ += upperDifference;
     2677          rowsIn[i] = 1;
     2678        } else if (getRowStatus(i) != basic) {
     2679          numberBinding++;
     2680          rowsIn[i] = 1;
     2681        }
     2682      }
     2683      if (numberBinding < rowsThreshold) {
     2684        // try
    27042685      } else {
    2705         // random?? - or give up
    2706         // Set up initial list
    2707         numberSort = 0;
    2708         for (int iRow = 0; iRow < numberRows_; iRow++) {
     2686        // random?? - or give up
     2687        // Set up initial list
     2688        numberSort = 0;
     2689        for (int iRow = 0; iRow < numberRows_; iRow++) {
    27092690          weight[iRow] = 1.123e50;
    27102691          if (rowLower_[iRow] == rowUpper_[iRow]) {
    2711             sort[numberSort++] = iRow;
    2712             weight[iRow] = 0.0;
    2713           }
    2714         }
    2715         numberSort /= 2;
    2716         // and pad out with random rows
    2717         double ratio = ((double)(smallNumberRows - numberSort)) / ((double) numberRows_);
    2718         for (int iRow = 0; iRow < numberRows_; iRow++) {
     2692            sort[numberSort++] = iRow;
     2693            weight[iRow] = 0.0;
     2694          }
     2695        }
     2696        numberSort /= 2;
     2697        // and pad out with random rows
     2698        double ratio = ((double)(smallNumberRows - numberSort)) / ((double)numberRows_);
     2699        for (int iRow = 0; iRow < numberRows_; iRow++) {
    27192700          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
    2720             sort[numberSort++] = iRow;
    2721         }
    2722         // sort
    2723         CoinSort_2(weight, weight + numberRows_, sort);
    2724         numberSort = CoinMin(numberRows_, smallNumberRows);
    2725         memset(rowsIn, 0, numberRows_);
    2726         for (int iRow = 0; iRow < numberSort; iRow++)
    2727           rowsIn[sort[iRow]] = 1;
     2701            sort[numberSort++] = iRow;
     2702        }
     2703        // sort
     2704        CoinSort_2(weight, weight + numberRows_, sort);
     2705        numberSort = CoinMin(numberRows_, smallNumberRows);
     2706        memset(rowsIn, 0, numberRows_);
     2707        for (int iRow = 0; iRow < numberSort; iRow++)
     2708          rowsIn[sort[iRow]] = 1;
    27282709      }
    27292710    }
    27302711    // see if feasible
    2731     numberPrimalInfeasibilities_=0;
     2712    numberPrimalInfeasibilities_ = 0;
    27322713    memset(rowActivity_, 0, numberRows_ * sizeof(double));
    27332714    times(1.0, columnActivity_, rowActivity_);
    2734     for (int i=0;i<numberRows_;i++) {
    2735       double lowerDifference = rowActivity_[i]-rowLower_[i];
    2736       double upperDifference = rowActivity_[i]-rowUpper_[i];
    2737       if (lowerDifference<-10*primalTolerance_||upperDifference>10.0*primalTolerance_) {
    2738         if (lowerDifference<0.0)
    2739           sumPrimalInfeasibilities_ -= lowerDifference;
    2740         else
    2741           sumPrimalInfeasibilities_ += upperDifference;
    2742         numberPrimalInfeasibilities_++;
     2715    for (int i = 0; i < numberRows_; i++) {
     2716      double lowerDifference = rowActivity_[i] - rowLower_[i];
     2717      double upperDifference = rowActivity_[i] - rowUpper_[i];
     2718      if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
     2719        if (lowerDifference < 0.0)
     2720          sumPrimalInfeasibilities_ -= lowerDifference;
     2721        else
     2722          sumPrimalInfeasibilities_ += upperDifference;
     2723        numberPrimalInfeasibilities_++;
    27432724      }
    27442725    }
    27452726    printf("Initial infeasibilities - %g (%d)\n",
    2746            sumPrimalInfeasibilities_,numberPrimalInfeasibilities_);
     2727      sumPrimalInfeasibilities_, numberPrimalInfeasibilities_);
    27472728    // Just do this number of passes
    27482729    int maxPass = 50;
     
    27502731    int takeOutPass = 30;
    27512732    int iPass;
    2752    
    2753     const CoinBigIndex * start = this->clpMatrix()->getVectorStarts();
    2754     const int * length = this->clpMatrix()->getVectorLengths();
    2755     const int * row = this->clpMatrix()->getIndices();
    2756     problemStatus_=1;
     2733
     2734    const CoinBigIndex *start = this->clpMatrix()->getVectorStarts();
     2735    const int *length = this->clpMatrix()->getVectorLengths();
     2736    const int *row = this->clpMatrix()->getIndices();
     2737    problemStatus_ = 1;
    27572738    for (iPass = 0; iPass < maxPass; iPass++) {
    27582739      printf("Start of pass %d\n", iPass);
    27592740      int numberSmallColumns = 0;
    27602741      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2761         int n = 0;
    2762         for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
    2763           int iRow = row[j];
    2764           if (rowsIn[iRow])
    2765             n++;
    2766         }
    2767         if (n)
    2768           whichColumns[numberSmallColumns++] = iColumn;
    2769       }
    2770       numberSort=0;
    2771       for (int i=0;i<numberRows_;i++) {
    2772         if (rowsIn[i])
    2773           sort[numberSort++]=i;
     2742        int n = 0;
     2743        for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
     2744          int iRow = row[j];
     2745          if (rowsIn[iRow])
     2746            n++;
     2747        }
     2748        if (n)
     2749          whichColumns[numberSmallColumns++] = iColumn;
     2750      }
     2751      numberSort = 0;
     2752      for (int i = 0; i < numberRows_; i++) {
     2753        if (rowsIn[i])
     2754          sort[numberSort++] = i;
    27742755      }
    27752756      // Create small problem
    27762757      ClpSimplex small(this, numberSort, sort, numberSmallColumns, whichColumns);
    27772758      printf("Small model has %d rows, %d columns and %d elements\n",
    2778              small.numberRows(),small.numberColumns(),small.getNumElements());
     2759        small.numberRows(), small.numberColumns(), small.getNumElements());
    27792760      small.setFactorizationFrequency(100 + numberSort / 200);
    27802761      // Solve
    2781       small.setLogLevel(CoinMax(0,logLevel()-1));
    2782       if (iPass>20) {
    2783         if (sumPrimalInfeasibilities_>1.0e-1) {
    2784           small.dual();
    2785         } else {
    2786           small.primal(1);
    2787         }
     2762      small.setLogLevel(CoinMax(0, logLevel() - 1));
     2763      if (iPass > 20) {
     2764        if (sumPrimalInfeasibilities_ > 1.0e-1) {
     2765          small.dual();
     2766        } else {
     2767          small.primal(1);
     2768        }
    27882769      } else {
    2789         // presolve
     2770        // presolve
    27902771#if 0
    27912772        ClpSolve::SolveType method;
     
    28012782#else
    28022783#if 1
    2803         ClpPresolve * pinfo = new ClpPresolve();
    2804         ClpSimplex * small2 = pinfo->presolvedModel(small,1.0e-5);
    2805         assert (small2);
    2806         if (sumPrimalInfeasibilities_>1.0e-1) {
    2807           small2->dual();
    2808         } else {
    2809           small2->primal(1);
    2810         }
    2811         pinfo->postsolve(true);
    2812         delete pinfo;
     2784        ClpPresolve *pinfo = new ClpPresolve();
     2785        ClpSimplex *small2 = pinfo->presolvedModel(small, 1.0e-5);
     2786        assert(small2);
     2787        if (sumPrimalInfeasibilities_ > 1.0e-1) {
     2788          small2->dual();
     2789        } else {
     2790          small2->primal(1);
     2791        }
     2792        pinfo->postsolve(true);
     2793        delete pinfo;
    28132794#else
    2814         char * types = new char[small.numberRows()+small.numberColumns()];
    2815         memset(types,0,small.numberRows()+small.numberColumns());
    2816         if (sumPrimalInfeasibilities_>1.0e-1) {
    2817           small.miniSolve(types,types+small.numberRows(),
    2818                           -1,0);
    2819         } else {
    2820           small.miniSolve(types,types+small.numberRows(),
    2821                           1,1);
    2822         }
    2823         delete [] types;
    2824 #endif
    2825         if (small.sumPrimalInfeasibilities()>1.0)
    2826           small.primal(1);
     2795        char *types = new char[small.numberRows() + small.numberColumns()];
     2796        memset(types, 0, small.numberRows() + small.numberColumns());
     2797        if (sumPrimalInfeasibilities_ > 1.0e-1) {
     2798          small.miniSolve(types, types + small.numberRows(),
     2799            -1, 0);
     2800        } else {
     2801          small.miniSolve(types, types + small.numberRows(),
     2802            1, 1);
     2803        }
     2804        delete[] types;
     2805#endif
     2806        if (small.sumPrimalInfeasibilities() > 1.0)
     2807          small.primal(1);
    28272808#endif
    28282809      }
    28292810      bool dualInfeasible = (small.status() == 2);
    28302811      // move solution back
    2831       const double * smallSolution = small.primalColumnSolution();
     2812      const double *smallSolution = small.primalColumnSolution();
    28322813      for (int j = 0; j < numberSmallColumns; j++) {
    2833         int iColumn = whichColumns[j];
    2834         columnActivity_[iColumn] = smallSolution[j];
    2835         this->setColumnStatus(iColumn, small.getColumnStatus(j));
     2814        int iColumn = whichColumns[j];
     2815        columnActivity_[iColumn] = smallSolution[j];
     2816        this->setColumnStatus(iColumn, small.getColumnStatus(j));
    28362817      }
    28372818      for (int iRow = 0; iRow < numberSort; iRow++) {
    2838         int kRow = sort[iRow];
    2839         this->setRowStatus(kRow, small.getRowStatus(iRow));
     2819        i