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/ClpSimplexPrimal.cpp

    r2373 r2385  
    7878 */
    7979
    80 
    8180#include "CoinPragma.hpp"
    8281
     
    102101#define ALT_UPDATE_WEIGHTS 2
    103102#endif
    104 #if ALT_UPDATE_WEIGHTS==1
    105 CoinIndexedVector * altVector[3]={NULL,NULL,NULL};
     103#if ALT_UPDATE_WEIGHTS == 1
     104CoinIndexedVector *altVector[3] = { NULL, NULL, NULL };
    106105#endif
    107106#ifdef CLP_USER_DRIVEN1
     
    116115   Variable will change theta if currentValue - currentTheta*alpha < 0.0
    117116*/
    118 bool userChoiceValid1(const ClpSimplex * model,
    119                       int sequenceOut,
    120                       double currentValue,
    121                       double currentTheta,
    122                       double alpha,
    123                       double realAlpha);
     117bool userChoiceValid1(const ClpSimplex *model,
     118  int sequenceOut,
     119  double currentValue,
     120  double currentTheta,
     121  double alpha,
     122  double realAlpha);
    124123/* This returns true if chosen in/out pair valid.
    125124   The main thing to check would be variable flipping bounds may be
     
    127126   If you return false sequenceIn_ will be flagged as ineligible.
    128127*/
    129 bool userChoiceValid2(const ClpSimplex * model);
     128bool userChoiceValid2(const ClpSimplex *model);
    130129/* If a good pivot then you may wish to unflag some variables.
    131130 */
    132 void userChoiceWasGood(ClpSimplex * model);
     131void userChoiceWasGood(ClpSimplex *model);
    133132#endif
    134133// primal
    135 int ClpSimplexPrimal::primal (int ifValuesPass , int startFinishOptions)
     134int ClpSimplexPrimal::primal(int ifValuesPass, int startFinishOptions)
    136135{
    137136
    138      /*
     137  /*
    139138         Method
    140139
     
    216215     */
    217216
    218      algorithm_ = +1;
    219      moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
    220 
    221      // save data
    222      ClpDataSave data = saveData();
    223      if (problemStatus_==10&&sumPrimalInfeasibilities_==-123456789.0) {
    224        // large infeasibility cost wanted
    225        infeasibilityCost_ = CoinMax(infeasibilityCost_,1.0e13);
    226      }
    227      matrix_->refresh(this); // make sure matrix okay
    228 
    229      // Save so can see if doing after dual
    230      int initialStatus = problemStatus_;
    231      int initialIterations = numberIterations_;
    232      int initialNegDjs = -1;
    233      // initialize - maybe values pass and algorithm_ is +1
     217  algorithm_ = +1;
     218  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
     219
     220  // save data
     221  ClpDataSave data = saveData();
     222  if (problemStatus_ == 10 && sumPrimalInfeasibilities_ == -123456789.0) {
     223    // large infeasibility cost wanted
     224    infeasibilityCost_ = CoinMax(infeasibilityCost_, 1.0e13);
     225  }
     226  matrix_->refresh(this); // make sure matrix okay
     227
     228  // Save so can see if doing after dual
     229  int initialStatus = problemStatus_;
     230  int initialIterations = numberIterations_;
     231  int initialNegDjs = -1;
     232  // initialize - maybe values pass and algorithm_ is +1
    234233#if 0
    235234     // if so - put in any superbasic costed slacks
     
    262261     }
    263262#endif
    264      // Start can skip some things in transposeTimes
    265      specialOptions_ |= 131072;
    266      if (!startup(ifValuesPass, startFinishOptions)) {
    267          // See if better to use all slack basis
    268          if (nonLinearCost_->sumInfeasibilities()>1.0e15) {
    269            // If not all slack then make it
    270            int numberRowBasic = 0;
    271            for (int i=0;i<numberRows_;i++) {
    272              if (getRowStatus(i)==basic)
    273                numberRowBasic++;
    274            }
    275            if (numberRowBasic<numberRows_) {
    276              allSlackBasis(true);
    277              int lastCleaned=-10000;
    278              statusOfProblemInPrimal(lastCleaned, 1, &progress_, true, ifValuesPass, NULL);
    279            }
    280         }
    281 
    282           // Set average theta
    283           nonLinearCost_->setAverageTheta(1.0e3);
    284           int lastCleaned = 0; // last time objective or bounds cleaned up
    285 
    286           // Say no pivot has occurred (for steepest edge and updates)
    287           pivotRow_ = -2;
    288 
    289           // This says whether to restore things etc
    290           int factorType = 0;
    291           if (problemStatus_ < 0 && perturbation_ < 100 && !ifValuesPass) {
    292                perturb(0);
    293                // Can't get here if values pass
    294                assert (!ifValuesPass);
    295                gutsOfSolution(NULL, NULL);
    296                if (handler_->logLevel() > 2) {
    297                     handler_->message(CLP_SIMPLEX_STATUS, messages_)
    298                               << numberIterations_ << objectiveValue();
    299                     handler_->printing(sumPrimalInfeasibilities_ > 0.0)
    300                               << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
    301                     handler_->printing(sumDualInfeasibilities_ > 0.0)
    302                               << sumDualInfeasibilities_ << numberDualInfeasibilities_;
    303                     handler_->printing(numberDualInfeasibilitiesWithoutFree_
    304                                        < numberDualInfeasibilities_)
    305                               << numberDualInfeasibilitiesWithoutFree_;
    306                     handler_->message() << CoinMessageEol;
    307                }
    308           }
    309           ClpSimplex * saveModel = NULL;
    310           int stopSprint = -1;
    311           int sprintPass = 0;
    312           int reasonableSprintIteration = 0;
    313           int lastSprintIteration = 0;
    314           double lastObjectiveValue = COIN_DBL_MAX;
    315           // Start check for cycles
    316           progress_.fillFromModel(this);
    317           progress_.startCheck();
    318           /*
     263  // Start can skip some things in transposeTimes
     264  specialOptions_ |= 131072;
     265  if (!startup(ifValuesPass, startFinishOptions)) {
     266    // See if better to use all slack basis
     267    if (nonLinearCost_->sumInfeasibilities() > 1.0e15) {
     268      // If not all slack then make it
     269      int numberRowBasic = 0;
     270      for (int i = 0; i < numberRows_; i++) {
     271        if (getRowStatus(i) == basic)
     272          numberRowBasic++;
     273      }
     274      if (numberRowBasic < numberRows_) {
     275        allSlackBasis(true);
     276        int lastCleaned = -10000;
     277        statusOfProblemInPrimal(lastCleaned, 1, &progress_, true, ifValuesPass, NULL);
     278      }
     279    }
     280
     281    // Set average theta
     282    nonLinearCost_->setAverageTheta(1.0e3);
     283    int lastCleaned = 0; // last time objective or bounds cleaned up
     284
     285    // Say no pivot has occurred (for steepest edge and updates)
     286    pivotRow_ = -2;
     287
     288    // This says whether to restore things etc
     289    int factorType = 0;
     290    if (problemStatus_ < 0 && perturbation_ < 100 && !ifValuesPass) {
     291      perturb(0);
     292      // Can't get here if values pass
     293      assert(!ifValuesPass);
     294      gutsOfSolution(NULL, NULL);
     295      if (handler_->logLevel() > 2) {
     296        handler_->message(CLP_SIMPLEX_STATUS, messages_)
     297          << numberIterations_ << objectiveValue();
     298        handler_->printing(sumPrimalInfeasibilities_ > 0.0)
     299          << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
     300        handler_->printing(sumDualInfeasibilities_ > 0.0)
     301          << sumDualInfeasibilities_ << numberDualInfeasibilities_;
     302        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     303          < numberDualInfeasibilities_)
     304          << numberDualInfeasibilitiesWithoutFree_;
     305        handler_->message() << CoinMessageEol;
     306      }
     307    }
     308    ClpSimplex *saveModel = NULL;
     309    int stopSprint = -1;
     310    int sprintPass = 0;
     311    int reasonableSprintIteration = 0;
     312    int lastSprintIteration = 0;
     313    double lastObjectiveValue = COIN_DBL_MAX;
     314    // Start check for cycles
     315    progress_.fillFromModel(this);
     316    progress_.startCheck();
     317    /*
    319318            Status of problem:
    320319            0 - optimal
     
    327326            -5 - looks unbounded
    328327          */
    329           while (problemStatus_ < 0) {
    330                int iRow, iColumn;
    331                // clear
    332                for (iRow = 0; iRow < 4; iRow++) {
    333                     rowArray_[iRow]->clear();
    334                }
    335 
    336                for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
    337                     columnArray_[iColumn]->clear();
    338                }
    339 
    340                // give matrix (and model costs and bounds a chance to be
    341                // refreshed (normally null)
    342                matrix_->refresh(this);
    343                // If getting nowhere - why not give it a kick
     328    while (problemStatus_ < 0) {
     329      int iRow, iColumn;
     330      // clear
     331      for (iRow = 0; iRow < 4; iRow++) {
     332        rowArray_[iRow]->clear();
     333      }
     334
     335      for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
     336        columnArray_[iColumn]->clear();
     337      }
     338
     339      // give matrix (and model costs and bounds a chance to be
     340      // refreshed (normally null)
     341      matrix_->refresh(this);
     342      // If getting nowhere - why not give it a kick
    344343#if 1
    345                if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
    346                          && initialStatus != 10) {
    347                     perturb(1);
    348                     matrix_->rhsOffset(this, true, false);
    349                }
    350 #endif
    351                // If we have done no iterations - special
    352                if (lastGoodIteration_ == numberIterations_ && factorType)
    353                     factorType = 3;
    354                if (saveModel) {
    355                     // Doing sprint
    356                     if (sequenceIn_ < 0 || numberIterations_ >= stopSprint) {
    357                          problemStatus_ = -1;
    358                          originalModel(saveModel);
    359                          saveModel = NULL;
    360                          if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration &&
    361                                    sprintPass > 100)
    362                               primalColumnPivot_->switchOffSprint();
    363                          //lastSprintIteration=numberIterations_;
    364                          COIN_DETAIL_PRINT(printf("End small model\n"));
    365                     }
    366                }
    367 
    368                // may factorize, checks if problem finished
    369                statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
    370                if (initialStatus == 10) {
    371                     // cleanup phase
    372                     if(initialIterations != numberIterations_) {
    373                          if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
    374                               // getting worse - try perturbing
    375                               if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
    376                                    perturb(1);
    377                                    matrix_->rhsOffset(this, true, false);
    378                                    statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
    379                               }
    380                          }
    381                     } else {
    382                          // save number of negative djs
    383                          if (!numberPrimalInfeasibilities_)
    384                               initialNegDjs = numberDualInfeasibilities_;
    385                          // make sure weight won't be changed
    386                          if (infeasibilityCost_ == 1.0e10)
    387                               infeasibilityCost_ = 1.000001e10;
    388                     }
    389                }
    390                // See if sprint says redo because of problems
    391                if (numberDualInfeasibilities_ == -776) {
    392                     // Need new set of variables
    393                     problemStatus_ = -1;
    394                     originalModel(saveModel);
    395                     saveModel = NULL;
    396                     //lastSprintIteration=numberIterations_;
    397                     COIN_DETAIL_PRINT(printf("End small model after\n"));
    398                     statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
    399                }
    400                int numberSprintIterations = 0;
    401                int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
    402                if (problemStatus_ == 777) {
    403                     // problems so do one pass with normal
    404                     problemStatus_ = -1;
    405                     originalModel(saveModel);
    406                     saveModel = NULL;
    407                     // Skip factorization
    408                     //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
    409                     statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
    410                } else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
    411                     int numberSort = 0;
    412                     int numberFixed = 0;
    413                     int numberBasic = 0;
    414                     reasonableSprintIteration = numberIterations_ + 100;
    415                     int * whichColumns = new int[numberColumns_];
    416                     double * weight = new double[numberColumns_];
    417                     int numberNegative = 0;
    418                     double sumNegative = 0.0;
    419                     // now massage weight so all basic in plus good djs
    420                     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    421                          double dj = dj_[iColumn];
    422                          switch(getColumnStatus(iColumn)) {
    423 
    424                          case basic:
    425                               dj = -1.0e50;
    426                               numberBasic++;
    427                               break;
    428                          case atUpperBound:
    429                               dj = -dj;
    430                               break;
    431                          case isFixed:
    432                               dj = 1.0e50;
    433                               numberFixed++;
    434                               break;
    435                          case atLowerBound:
    436                               dj = dj;
    437                               break;
    438                          case isFree:
    439                               dj = -100.0 * fabs(dj);
    440                               break;
    441                          case superBasic:
    442                               dj = -100.0 * fabs(dj);
    443                               break;
    444                          }
    445                          if (dj < -dualTolerance_ && dj > -1.0e50) {
    446                               numberNegative++;
    447                               sumNegative -= dj;
    448                          }
    449                          weight[iColumn] = dj;
    450                          whichColumns[iColumn] = iColumn;
    451                     }
    452                     handler_->message(CLP_SPRINT, messages_)
    453                               << sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
    454                               << numberNegative
    455                               << CoinMessageEol;
    456                     sprintPass++;
    457                     lastSprintIteration = numberIterations_;
    458                     if (objectiveValue()*optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
    459                          // switch off
    460                       COIN_DETAIL_PRINT(printf("Switching off sprint\n"));
    461                          primalColumnPivot_->switchOffSprint();
    462                     } else {
    463                          lastObjectiveValue = objectiveValue() * optimizationDirection_;
    464                          // sort
    465                          CoinSort_2(weight, weight + numberColumns_, whichColumns);
    466                          numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
    467                          // Sort to make consistent ?
    468                          std::sort(whichColumns, whichColumns + numberSort);
    469                          saveModel = new ClpSimplex(this, numberSort, whichColumns);
    470                          delete [] whichColumns;
    471                          delete [] weight;
    472                          // Skip factorization
    473                          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
    474                          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
    475                          stopSprint = numberIterations_ + numberSprintIterations;
    476                          COIN_DETAIL_PRINT(printf("Sprint with %d columns for %d iterations\n",
    477                                                   numberSprintColumns, numberSprintIterations));
    478                     }
    479                }
    480 
    481                // Say good factorization
    482                factorType = 1;
    483 
    484                // Say no pivot has occurred (for steepest edge and updates)
    485                pivotRow_ = -2;
    486                // exit if feasible and done enough iterations
    487                if ((moreSpecialOptions_&1048576)!=0) {
    488                  int maxIterations=maximumIterations()-1000000;
    489                  if (maxIterations>0&&maxIterations<200000) {
    490                    if (!nonLinearCost_->numberInfeasibilities()&&
    491                        numberIterations_>=maxIterations) {
    492                      problemStatus_ = 3;
    493                      secondaryStatus_ = 10;
    494                    }
    495                  }
    496                }
    497                // exit if victory declared
    498                if (problemStatus_ >= 0) {
     344      if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
     345        && initialStatus != 10) {
     346        perturb(1);
     347        matrix_->rhsOffset(this, true, false);
     348      }
     349#endif
     350      // If we have done no iterations - special
     351      if (lastGoodIteration_ == numberIterations_ && factorType)
     352        factorType = 3;
     353      if (saveModel) {
     354        // Doing sprint
     355        if (sequenceIn_ < 0 || numberIterations_ >= stopSprint) {
     356          problemStatus_ = -1;
     357          originalModel(saveModel);
     358          saveModel = NULL;
     359          if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration && sprintPass > 100)
     360            primalColumnPivot_->switchOffSprint();
     361          //lastSprintIteration=numberIterations_;
     362          COIN_DETAIL_PRINT(printf("End small model\n"));
     363        }
     364      }
     365
     366      // may factorize, checks if problem finished
     367      statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
     368      if (initialStatus == 10) {
     369        // cleanup phase
     370        if (initialIterations != numberIterations_) {
     371          if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
     372            // getting worse - try perturbing
     373            if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
     374              perturb(1);
     375              matrix_->rhsOffset(this, true, false);
     376              statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
     377            }
     378          }
     379        } else {
     380          // save number of negative djs
     381          if (!numberPrimalInfeasibilities_)
     382            initialNegDjs = numberDualInfeasibilities_;
     383          // make sure weight won't be changed
     384          if (infeasibilityCost_ == 1.0e10)
     385            infeasibilityCost_ = 1.000001e10;
     386        }
     387      }
     388      // See if sprint says redo because of problems
     389      if (numberDualInfeasibilities_ == -776) {
     390        // Need new set of variables
     391        problemStatus_ = -1;
     392        originalModel(saveModel);
     393        saveModel = NULL;
     394        //lastSprintIteration=numberIterations_;
     395        COIN_DETAIL_PRINT(printf("End small model after\n"));
     396        statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
     397      }
     398      int numberSprintIterations = 0;
     399      int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
     400      if (problemStatus_ == 777) {
     401        // problems so do one pass with normal
     402        problemStatus_ = -1;
     403        originalModel(saveModel);
     404        saveModel = NULL;
     405        // Skip factorization
     406        //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
     407        statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
     408      } else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
     409        int numberSort = 0;
     410        int numberFixed = 0;
     411        int numberBasic = 0;
     412        reasonableSprintIteration = numberIterations_ + 100;
     413        int *whichColumns = new int[numberColumns_];
     414        double *weight = new double[numberColumns_];
     415        int numberNegative = 0;
     416        double sumNegative = 0.0;
     417        // now massage weight so all basic in plus good djs
     418        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     419          double dj = dj_[iColumn];
     420          switch (getColumnStatus(iColumn)) {
     421
     422          case basic:
     423            dj = -1.0e50;
     424            numberBasic++;
     425            break;
     426          case atUpperBound:
     427            dj = -dj;
     428            break;
     429          case isFixed:
     430            dj = 1.0e50;
     431            numberFixed++;
     432            break;
     433          case atLowerBound:
     434            dj = dj;
     435            break;
     436          case isFree:
     437            dj = -100.0 * fabs(dj);
     438            break;
     439          case superBasic:
     440            dj = -100.0 * fabs(dj);
     441            break;
     442          }
     443          if (dj < -dualTolerance_ && dj > -1.0e50) {
     444            numberNegative++;
     445            sumNegative -= dj;
     446          }
     447          weight[iColumn] = dj;
     448          whichColumns[iColumn] = iColumn;
     449        }
     450        handler_->message(CLP_SPRINT, messages_)
     451          << sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
     452          << numberNegative
     453          << CoinMessageEol;
     454        sprintPass++;
     455        lastSprintIteration = numberIterations_;
     456        if (objectiveValue() * optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
     457          // switch off
     458          COIN_DETAIL_PRINT(printf("Switching off sprint\n"));
     459          primalColumnPivot_->switchOffSprint();
     460        } else {
     461          lastObjectiveValue = objectiveValue() * optimizationDirection_;
     462          // sort
     463          CoinSort_2(weight, weight + numberColumns_, whichColumns);
     464          numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
     465          // Sort to make consistent ?
     466          std::sort(whichColumns, whichColumns + numberSort);
     467          saveModel = new ClpSimplex(this, numberSort, whichColumns);
     468          delete[] whichColumns;
     469          delete[] weight;
     470          // Skip factorization
     471          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
     472          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
     473          stopSprint = numberIterations_ + numberSprintIterations;
     474          COIN_DETAIL_PRINT(printf("Sprint with %d columns for %d iterations\n",
     475            numberSprintColumns, numberSprintIterations));
     476        }
     477      }
     478
     479      // Say good factorization
     480      factorType = 1;
     481
     482      // Say no pivot has occurred (for steepest edge and updates)
     483      pivotRow_ = -2;
     484      // exit if feasible and done enough iterations
     485      if ((moreSpecialOptions_ & 1048576) != 0) {
     486        int maxIterations = maximumIterations() - 1000000;
     487        if (maxIterations > 0 && maxIterations < 200000) {
     488          if (!nonLinearCost_->numberInfeasibilities() && numberIterations_ >= maxIterations) {
     489            problemStatus_ = 3;
     490            secondaryStatus_ = 10;
     491          }
     492        }
     493      }
     494      // exit if victory declared
     495      if (problemStatus_ >= 0) {
    499496#ifdef CLP_USER_DRIVEN
    500                  int status =
    501                    eventHandler_->event(ClpEventHandler::endInPrimal);
    502                  if (status>=0&&status<10) {
    503                    // carry on
    504                    problemStatus_=-1;
    505                    if (status==0)
    506                      continue; // re-factorize
    507                  } else if (status>=10) {
    508                    problemStatus_=status-10;
    509                    break;
    510                  } else {
    511                    break;
    512                  }
     497        int status = eventHandler_->event(ClpEventHandler::endInPrimal);
     498        if (status >= 0 && status < 10) {
     499          // carry on
     500          problemStatus_ = -1;
     501          if (status == 0)
     502            continue; // re-factorize
     503        } else if (status >= 10) {
     504          problemStatus_ = status - 10;
     505          break;
     506        } else {
     507          break;
     508        }
    513509#else
    514                  break;
    515 #endif
    516                }
    517 
    518                // test for maximum iterations
    519                if (hitMaximumIterations() ||
    520                    (ifValuesPass == 2 && firstFree_ < 0)) {
    521                     problemStatus_ = 3;
    522                     break;
    523                } else if ((moreSpecialOptions_&524288)!=0&&
    524                           !nonLinearCost_->numberInfeasibilities()&&
    525                           fabs(dblParam_[ClpDualObjectiveLimit])>1.0e30) {
    526                     problemStatus_ = 3;
    527                     secondaryStatus_ = 10;
    528                     break;
    529                }
    530 
    531                if (firstFree_ < 0) {
    532                     if (ifValuesPass) {
    533                          // end of values pass
    534                          ifValuesPass = 0;
    535                          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
    536                          if (status >= 0) {
    537                               problemStatus_ = 5;
    538                               secondaryStatus_ = ClpEventHandler::endOfValuesPass;
    539                               break;
    540                          }
    541                          //#define FEB_TRY
     510        break;
     511#endif
     512      }
     513
     514      // test for maximum iterations
     515      if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
     516        problemStatus_ = 3;
     517        break;
     518      } else if ((moreSpecialOptions_ & 524288) != 0 && !nonLinearCost_->numberInfeasibilities() && fabs(dblParam_[ClpDualObjectiveLimit]) > 1.0e30) {
     519        problemStatus_ = 3;
     520        secondaryStatus_ = 10;
     521        break;
     522      }
     523
     524      if (firstFree_ < 0) {
     525        if (ifValuesPass) {
     526          // end of values pass
     527          ifValuesPass = 0;
     528          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
     529          if (status >= 0) {
     530            problemStatus_ = 5;
     531            secondaryStatus_ = ClpEventHandler::endOfValuesPass;
     532            break;
     533          }
     534          //#define FEB_TRY
    542535#if 1 //def FEB_TRY
    543                          if (perturbation_ < 100)
    544                               perturb(0);
    545 #endif
    546                     }
    547                }
    548                // Check event
    549                {
    550                     int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
    551                     if (status >= 0) {
    552                         // if >=100 - then special e.g. unperturb
    553                         if (status!=101) {
    554                           problemStatus_ = 5;
    555                           secondaryStatus_ = ClpEventHandler::endOfFactorization;
    556                           break;
    557                         } else {
    558                           unPerturb();
    559                           continue;
    560                         }
    561                     }
    562                }
    563                // Iterate
    564                whileIterating(ifValuesPass ? 1 : 0);
    565                if (sequenceIn_<0&&ifValuesPass==2)
    566                  problemStatus_=3; // user wants to exit
    567           }
    568      }
    569      // if infeasible get real values
    570      //printf("XXXXY final cost %g\n",infeasibilityCost_);
    571      progress_.initialWeight_ = 0.0;
    572      if (problemStatus_ == 1 && secondaryStatus_ != 6) {
    573           double saveWeight = infeasibilityCost_;
     536          if (perturbation_ < 100)
     537            perturb(0);
     538#endif
     539        }
     540      }
     541      // Check event
     542      {
     543        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
     544        if (status >= 0) {
     545          // if >=100 - then special e.g. unperturb
     546          if (status != 101) {
     547            problemStatus_ = 5;
     548            secondaryStatus_ = ClpEventHandler::endOfFactorization;
     549            break;
     550          } else {
     551            unPerturb();
     552            continue;
     553          }
     554        }
     555      }
     556      // Iterate
     557      whileIterating(ifValuesPass ? 1 : 0);
     558      if (sequenceIn_ < 0 && ifValuesPass == 2)
     559        problemStatus_ = 3; // user wants to exit
     560    }
     561  }
     562  // if infeasible get real values
     563  //printf("XXXXY final cost %g\n",infeasibilityCost_);
     564  progress_.initialWeight_ = 0.0;
     565  if (problemStatus_ == 1 && secondaryStatus_ != 6) {
     566    double saveWeight = infeasibilityCost_;
    574567#ifndef WANT_INFEASIBLE_DUALS
    575           infeasibilityCost_ = 0.0;
    576           createRim(1 + 4);
     568    infeasibilityCost_ = 0.0;
     569    createRim(1 + 4);
    577570#else
    578           infeasibilityCost_ = 1.0;
    579           createRim(1);
    580           memset(cost_,0,(numberRows_+numberColumns_)*sizeof(double));
    581 #endif
    582           delete nonLinearCost_;
    583           nonLinearCost_ = new ClpNonLinearCost(this);
    584           nonLinearCost_->checkInfeasibilities(0.0);
    585           sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
    586           numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
    587           // and get good feasible duals
    588           computeDuals(NULL);
    589           infeasibilityCost_=saveWeight;
    590      }
    591      // Stop can skip some things in transposeTimes
    592      specialOptions_ &= ~131072;
    593      // clean up
    594      unflag();
    595      finish(startFinishOptions);
    596      restoreData(data);
    597      return problemStatus_;
     571    infeasibilityCost_ = 1.0;
     572    createRim(1);
     573    memset(cost_, 0, (numberRows_ + numberColumns_) * sizeof(double));
     574#endif
     575    delete nonLinearCost_;
     576    nonLinearCost_ = new ClpNonLinearCost(this);
     577    nonLinearCost_->checkInfeasibilities(0.0);
     578    sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
     579    numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
     580    // and get good feasible duals
     581    computeDuals(NULL);
     582    infeasibilityCost_ = saveWeight;
     583  }
     584  // Stop can skip some things in transposeTimes
     585  specialOptions_ &= ~131072;
     586  // clean up
     587  unflag();
     588  finish(startFinishOptions);
     589  restoreData(data);
     590  return problemStatus_;
    598591}
    599592/*
     
    607600  +3 max iterations
    608601*/
    609 int
    610 ClpSimplexPrimal::whileIterating(int valuesOption)
     602int ClpSimplexPrimal::whileIterating(int valuesOption)
    611603{
    612      // Say if values pass
    613      int ifValuesPass = (firstFree_ >= 0) ? 1 : 0;
    614      int returnCode = -1;
    615      int superBasicType = 1;
    616      if (valuesOption > 1)
    617           superBasicType = 3;
    618      // delete any rays
    619      delete [] ray_;
    620      ray_ = NULL;
    621      // status stays at -1 while iterating, >=0 finished, -2 to invert
    622      // status -3 to go to top without an invert
    623      while (problemStatus_ == -1) {
    624           //#define CLP_DEBUG 1
     604  // Say if values pass
     605  int ifValuesPass = (firstFree_ >= 0) ? 1 : 0;
     606  int returnCode = -1;
     607  int superBasicType = 1;
     608  if (valuesOption > 1)
     609    superBasicType = 3;
     610  // delete any rays
     611  delete[] ray_;
     612  ray_ = NULL;
     613  // status stays at -1 while iterating, >=0 finished, -2 to invert
     614  // status -3 to go to top without an invert
     615  while (problemStatus_ == -1) {
     616    //#define CLP_DEBUG 1
    625617#ifdef CLP_DEBUG
    626           {
    627                int i;
    628                // not [1] as has information
    629                for (i = 0; i < 4; i++) {
    630                     if (i != 1)
    631                          rowArray_[i]->checkClear();
    632                }
    633                for (i = 0; i < 2; i++) {
    634                     columnArray_[i]->checkClear();
    635                }
    636           }
     618    {
     619      int i;
     620      // not [1] as has information
     621      for (i = 0; i < 4; i++) {
     622        if (i != 1)
     623          rowArray_[i]->checkClear();
     624      }
     625      for (i = 0; i < 2; i++) {
     626        columnArray_[i]->checkClear();
     627      }
     628    }
    637629#endif
    638630#if 0
     
    663655                 nonLinearCost_->numberInfeasibilities());
    664656#endif
    665 #if CLP_DEBUG>2
    666           // very expensive
    667           if (numberIterations_ > 0 && numberIterations_ < 100 && !ifValuesPass) {
    668                handler_->setLogLevel(63);
    669                double saveValue = objectiveValue_;
    670                double * saveRow1 = new double[numberRows_];
    671                double * saveRow2 = new double[numberRows_];
    672                CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
    673                CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
    674                double * saveColumn1 = new double[numberColumns_];
    675                double * saveColumn2 = new double[numberColumns_];
    676                CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
    677                CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
    678                gutsOfSolution(NULL, NULL, false);
    679                printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
    680                       numberIterations_,
    681                       saveValue, objectiveValue_, sumPrimalInfeasibilities_);
    682                CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
    683                CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
    684                CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
    685                CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
    686                delete [] saveRow1;
    687                delete [] saveRow2;
    688                delete [] saveColumn1;
    689                delete [] saveColumn2;
    690                objectiveValue_ = saveValue;
    691           }
    692 #endif
    693           if (!ifValuesPass) {
    694                // choose column to come in
    695                // can use pivotRow_ to update weights
    696                // pass in list of cost changes so can do row updates (rowArray_[1])
    697                // NOTE rowArray_[0] is used by computeDuals which is a
    698                // slow way of getting duals but might be used
    699 #ifdef LONG_REGION_2 
    700                primalColumn(rowArray_[1], rowArray_[0], rowArray_[3],
    701                             columnArray_[0], rowArray_[2]);
    702 #else 
    703                primalColumn(rowArray_[1], rowArray_[2], rowArray_[3],
    704                             columnArray_[0], columnArray_[1]);
    705 #endif 
    706           } else {
    707                // in values pass
    708                int sequenceIn = nextSuperBasic(superBasicType, columnArray_[0]);
    709                if (valuesOption > 1)
    710                     superBasicType = 2;
    711                if (sequenceIn < 0) {
    712                     // end of values pass - initialize weights etc
    713                     handler_->message(CLP_END_VALUES_PASS, messages_)
    714                               << numberIterations_;
    715                     primalColumnPivot_->saveWeights(this, 5);
    716                     problemStatus_ = -2; // factorize now
    717                     pivotRow_ = -1; // say no weights update
    718                     returnCode = -4;
    719                     // Clean up
    720                     int i;
    721                     for (i = 0; i < numberRows_ + numberColumns_; i++) {
    722                          if (getColumnStatus(i) == atLowerBound || getColumnStatus(i) == isFixed)
    723                               solution_[i] = lower_[i];
    724                          else if (getColumnStatus(i) == atUpperBound)
    725                               solution_[i] = upper_[i];
    726                     }
    727                     break;
    728                } else {
    729                     // normal
    730                     sequenceIn_ = sequenceIn;
    731                     valueIn_ = solution_[sequenceIn_];
    732                     lowerIn_ = lower_[sequenceIn_];
    733                     upperIn_ = upper_[sequenceIn_];
    734                     dualIn_ = dj_[sequenceIn_];
    735                }
    736           }
    737           pivotRow_ = -1;
    738           sequenceOut_ = -1;
    739           rowArray_[1]->clear();
    740           if (sequenceIn_ >= 0) {
    741                // we found a pivot column
    742                assert (!flagged(sequenceIn_));
     657#if CLP_DEBUG > 2
     658    // very expensive
     659    if (numberIterations_ > 0 && numberIterations_ < 100 && !ifValuesPass) {
     660      handler_->setLogLevel(63);
     661      double saveValue = objectiveValue_;
     662      double *saveRow1 = new double[numberRows_];
     663      double *saveRow2 = new double[numberRows_];
     664      CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
     665      CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
     666      double *saveColumn1 = new double[numberColumns_];
     667      double *saveColumn2 = new double[numberColumns_];
     668      CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
     669      CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
     670      gutsOfSolution(NULL, NULL, false);
     671      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
     672        numberIterations_,
     673        saveValue, objectiveValue_, sumPrimalInfeasibilities_);
     674      CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
     675      CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
     676      CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
     677      CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
     678      delete[] saveRow1;
     679      delete[] saveRow2;
     680      delete[] saveColumn1;
     681      delete[] saveColumn2;
     682      objectiveValue_ = saveValue;
     683    }
     684#endif
     685    if (!ifValuesPass) {
     686      // choose column to come in
     687      // can use pivotRow_ to update weights
     688      // pass in list of cost changes so can do row updates (rowArray_[1])
     689      // NOTE rowArray_[0] is used by computeDuals which is a
     690      // slow way of getting duals but might be used
     691#ifdef LONG_REGION_2
     692      primalColumn(rowArray_[1], rowArray_[0], rowArray_[3],
     693        columnArray_[0], rowArray_[2]);
     694#else
     695      primalColumn(rowArray_[1], rowArray_[2], rowArray_[3],
     696        columnArray_[0], columnArray_[1]);
     697#endif
     698    } else {
     699      // in values pass
     700      int sequenceIn = nextSuperBasic(superBasicType, columnArray_[0]);
     701      if (valuesOption > 1)
     702        superBasicType = 2;
     703      if (sequenceIn < 0) {
     704        // end of values pass - initialize weights etc
     705        handler_->message(CLP_END_VALUES_PASS, messages_)
     706          << numberIterations_;
     707        primalColumnPivot_->saveWeights(this, 5);
     708        problemStatus_ = -2; // factorize now
     709        pivotRow_ = -1; // say no weights update
     710        returnCode = -4;
     711        // Clean up
     712        int i;
     713        for (i = 0; i < numberRows_ + numberColumns_; i++) {
     714          if (getColumnStatus(i) == atLowerBound || getColumnStatus(i) == isFixed)
     715            solution_[i] = lower_[i];
     716          else if (getColumnStatus(i) == atUpperBound)
     717            solution_[i] = upper_[i];
     718        }
     719        break;
     720      } else {
     721        // normal
     722        sequenceIn_ = sequenceIn;
     723        valueIn_ = solution_[sequenceIn_];
     724        lowerIn_ = lower_[sequenceIn_];
     725        upperIn_ = upper_[sequenceIn_];
     726        dualIn_ = dj_[sequenceIn_];
     727      }
     728    }
     729    pivotRow_ = -1;
     730    sequenceOut_ = -1;
     731    rowArray_[1]->clear();
     732    if (sequenceIn_ >= 0) {
     733      // we found a pivot column
     734      assert(!flagged(sequenceIn_));
    743735#ifdef CLP_DEBUG
    744                if ((handler_->logLevel() & 32)) {
    745                     char x = isColumn(sequenceIn_) ? 'C' : 'R';
    746                     std::cout << "pivot column " <<
    747                               x << sequenceWithin(sequenceIn_) << std::endl;
    748                }
     736      if ((handler_->logLevel() & 32)) {
     737        char x = isColumn(sequenceIn_) ? 'C' : 'R';
     738        std::cout << "pivot column " << x << sequenceWithin(sequenceIn_) << std::endl;
     739      }
    749740#endif
    750741#ifdef CLP_DEBUG
    751                {
    752                     int checkSequence = -2077;
    753                     if (checkSequence >= 0 && checkSequence < numberRows_ + numberColumns_ && !ifValuesPass) {
    754                          rowArray_[2]->checkClear();
    755                          rowArray_[3]->checkClear();
    756                          double * array = rowArray_[3]->denseVector();
    757                          int * index = rowArray_[3]->getIndices();
    758                          unpackPacked(rowArray_[3], checkSequence);
    759                          factorization_->updateColumnForDebug(rowArray_[2], rowArray_[3]);
    760                          int number = rowArray_[3]->getNumElements();
    761                          double dualIn = cost_[checkSequence];
    762                          int i;
    763                          for (i = 0; i < number; i++) {
    764                               int iRow = index[i];
    765                               int iPivot = pivotVariable_[iRow];
    766                               double alpha = array[i];
    767                               dualIn -= alpha * cost_[iPivot];
    768                          }
    769                          printf("old dj for %d was %g, recomputed %g\n", checkSequence,
    770                                 dj_[checkSequence], dualIn);
    771                          rowArray_[3]->clear();
    772                          if (numberIterations_ > 2000)
    773                               exit(1);
    774                     }
    775                }
    776 #endif
    777                // do second half of iteration
    778                returnCode = pivotResult(ifValuesPass);
    779                if (returnCode < -1 && returnCode > -5) {
    780                     problemStatus_ = -2; //
    781                } else if (returnCode == -5) {
    782                     if ((moreSpecialOptions_ & 16) == 0 && factorization_->pivots()) {
    783                          moreSpecialOptions_ |= 16;
    784                          problemStatus_ = -2;
    785                     }
    786                     // otherwise something flagged - continue;
    787                } else if (returnCode == 2) {
    788                     problemStatus_ = -5; // looks unbounded
    789                } else if (returnCode == 4) {
    790                     problemStatus_ = -2; // looks unbounded but has iterated
    791                } else if (returnCode != -1) {
    792                     assert(returnCode == 3);
    793                     if (problemStatus_ != 5)
    794                          problemStatus_ = 3;
    795                }
    796           } else {
    797                // no pivot column
     742      {
     743        int checkSequence = -2077;
     744        if (checkSequence >= 0 && checkSequence < numberRows_ + numberColumns_ && !ifValuesPass) {
     745          rowArray_[2]->checkClear();
     746          rowArray_[3]->checkClear();
     747          double *array = rowArray_[3]->denseVector();
     748          int *index = rowArray_[3]->getIndices();
     749          unpackPacked(rowArray_[3], checkSequence);
     750          factorization_->updateColumnForDebug(rowArray_[2], rowArray_[3]);
     751          int number = rowArray_[3]->getNumElements();
     752          double dualIn = cost_[checkSequence];
     753          int i;
     754          for (i = 0; i < number; i++) {
     755            int iRow = index[i];
     756            int iPivot = pivotVariable_[iRow];
     757            double alpha = array[i];
     758            dualIn -= alpha * cost_[iPivot];
     759          }
     760          printf("old dj for %d was %g, recomputed %g\n", checkSequence,
     761            dj_[checkSequence], dualIn);
     762          rowArray_[3]->clear();
     763          if (numberIterations_ > 2000)
     764            exit(1);
     765        }
     766      }
     767#endif
     768      // do second half of iteration
     769      returnCode = pivotResult(ifValuesPass);
     770      if (returnCode < -1 && returnCode > -5) {
     771        problemStatus_ = -2; //
     772      } else if (returnCode == -5) {
     773        if ((moreSpecialOptions_ & 16) == 0 && factorization_->pivots()) {
     774          moreSpecialOptions_ |= 16;
     775          problemStatus_ = -2;
     776        }
     777        // otherwise something flagged - continue;
     778      } else if (returnCode == 2) {
     779        problemStatus_ = -5; // looks unbounded
     780      } else if (returnCode == 4) {
     781        problemStatus_ = -2; // looks unbounded but has iterated
     782      } else if (returnCode != -1) {
     783        assert(returnCode == 3);
     784        if (problemStatus_ != 5)
     785          problemStatus_ = 3;
     786      }
     787    } else {
     788      // no pivot column
    798789#ifdef CLP_DEBUG
    799                if (handler_->logLevel() & 32)
    800                     printf("** no column pivot\n");
    801 #endif
    802                if (nonLinearCost_->numberInfeasibilities())
    803                     problemStatus_ = -4; // might be infeasible
    804                // Force to re-factorize early next time
    805                int numberPivots = factorization_->pivots();
    806                returnCode = 0;
     790      if (handler_->logLevel() & 32)
     791        printf("** no column pivot\n");
     792#endif
     793      if (nonLinearCost_->numberInfeasibilities())
     794        problemStatus_ = -4; // might be infeasible
     795      // Force to re-factorize early next time
     796      int numberPivots = factorization_->pivots();
     797      returnCode = 0;
    807798#ifdef CLP_USER_DRIVEN
    808                // If large number of pivots trap later?
    809                if (problemStatus_==-1 && numberPivots<1000000) {
    810                 int status = eventHandler_->event(ClpEventHandler::noCandidateInPrimal);
    811                  if (status>=0&&status<10) {
    812                    // carry on
    813                    problemStatus_=-1;
    814                    if (status==0)
    815                      break;
    816                  } else if (status>=10) {
    817                      problemStatus_=status-10;
    818                      break;
    819                 } else {
    820                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
    821                    break;
    822                 }
    823                }
     799      // If large number of pivots trap later?
     800      if (problemStatus_ == -1 && numberPivots < 1000000) {
     801        int status = eventHandler_->event(ClpEventHandler::noCandidateInPrimal);
     802        if (status >= 0 && status < 10) {
     803          // carry on
     804          problemStatus_ = -1;
     805          if (status == 0)
     806            break;
     807        } else if (status >= 10) {
     808          problemStatus_ = status - 10;
     809          break;
     810        } else {
     811          forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
     812          break;
     813        }
     814      }
    824815#else
    825                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
    826                    break;
    827 #endif
    828           }
    829      }
    830      if (valuesOption > 1)
    831           columnArray_[0]->setNumElements(0);
    832      return returnCode;
     816      forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
     817      break;
     818#endif
     819    }
     820  }
     821  if (valuesOption > 1)
     822    columnArray_[0]->setNumElements(0);
     823  return returnCode;
    833824}
    834825/* Checks if finished.  Updates status */
    835 void
    836 ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned, int type,
    837           ClpSimplexProgress * progress,
    838           bool doFactorization,
    839           int ifValuesPass,
    840           ClpSimplex * originalModel)
     826void ClpSimplexPrimal::statusOfProblemInPrimal(int &lastCleaned, int type,
     827  ClpSimplexProgress *progress,
     828  bool doFactorization,
     829  int ifValuesPass,
     830  ClpSimplex *originalModel)
    841831{
    842      int dummy; // for use in generalExpanded
    843      int saveFirstFree = firstFree_;
    844      double saveOriginalWeight = infeasibilityCost_;
    845      // number of pivots done
    846      int numberPivots = factorization_->pivots();
    847      if (type == 2) {
    848           // trouble - restore solution
    849           CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
    850           CoinMemcpyN(savedSolution_ + numberColumns_ ,
    851                       numberRows_, rowActivityWork_);
    852           CoinMemcpyN(savedSolution_ ,
    853                       numberColumns_, columnActivityWork_);
    854           // restore extra stuff
    855           matrix_->generalExpanded(this, 6, dummy);
    856           forceFactorization_ = 1; // a bit drastic but ..
    857           pivotRow_ = -1; // say no weights update
     832  int dummy; // for use in generalExpanded
     833  int saveFirstFree = firstFree_;
     834  double saveOriginalWeight = infeasibilityCost_;
     835  // number of pivots done
     836  int numberPivots = factorization_->pivots();
     837  if (type == 2) {
     838    // trouble - restore solution
     839    CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
     840    CoinMemcpyN(savedSolution_ + numberColumns_,
     841      numberRows_, rowActivityWork_);
     842    CoinMemcpyN(savedSolution_,
     843      numberColumns_, columnActivityWork_);
     844    // restore extra stuff
     845    matrix_->generalExpanded(this, 6, dummy);
     846    forceFactorization_ = 1; // a bit drastic but ..
     847    pivotRow_ = -1; // say no weights update
     848    changeMade_++; // say change made
     849  }
     850  int saveThreshold = factorization_->sparseThreshold();
     851  int tentativeStatus = problemStatus_;
     852  int numberThrownOut = 1; // to loop round on bad factorization in values pass
     853  double lastSumInfeasibility = COIN_DBL_MAX;
     854  if (numberIterations_)
     855    lastSumInfeasibility = nonLinearCost_->sumInfeasibilities();
     856  int nPass = 0;
     857  while (numberThrownOut) {
     858    int nSlackBasic = 0;
     859    if (nPass) {
     860      for (int i = 0; i < numberRows_; i++) {
     861        if (getRowStatus(i) == basic)
     862          nSlackBasic++;
     863      }
     864    }
     865    nPass++;
     866    if (problemStatus_ > -3 || problemStatus_ == -4) {
     867      // factorize
     868      // later on we will need to recover from singularities
     869      // also we could skip if first time
     870      // do weights
     871      // This may save pivotRow_ for use
     872      if (doFactorization)
     873        primalColumnPivot_->saveWeights(this, 1);
     874
     875      if ((type && doFactorization) || nSlackBasic == numberRows_) {
     876        // is factorization okay?
     877        int factorStatus = internalFactorize(1);
     878        if (factorStatus) {
     879          if (solveType_ == 2 + 8) {
     880            // say odd
     881            problemStatus_ = 5;
     882            return;
     883          }
     884          if (type != 1 || largestPrimalError_ > 1.0e3
     885            || largestDualError_ > 1.0e3) {
     886            // switch off dense
     887            int saveDense = factorization_->denseThreshold();
     888            factorization_->setDenseThreshold(0);
     889            // Go to safe
     890            factorization_->pivotTolerance(0.99);
     891            // make sure will do safe factorization
     892            pivotVariable_[0] = -1;
     893            internalFactorize(2);
     894            factorization_->setDenseThreshold(saveDense);
     895            // restore extra stuff
     896            matrix_->generalExpanded(this, 6, dummy);
     897          } else {
     898            // no - restore previous basis
     899            // Keep any flagged variables
     900            int i;
     901            for (i = 0; i < numberRows_ + numberColumns_; i++) {
     902              if (flagged(i))
     903                saveStatus_[i] |= 64; //say flagged
     904            }
     905            CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
     906            if (numberPivots <= 1) {
     907              // throw out something
     908              if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
     909                setFlagged(sequenceIn_);
     910              } else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
     911                setFlagged(sequenceOut_);
     912              }
     913              double newTolerance = CoinMax(0.5 + 0.499 * randomNumberGenerator_.randomDouble(), factorization_->pivotTolerance());
     914              factorization_->pivotTolerance(newTolerance);
     915            } else {
     916              // Go to safe
     917              factorization_->pivotTolerance(0.99);
     918            }
     919            CoinMemcpyN(savedSolution_ + numberColumns_,
     920              numberRows_, rowActivityWork_);
     921            CoinMemcpyN(savedSolution_,
     922              numberColumns_, columnActivityWork_);
     923            // restore extra stuff
     924            matrix_->generalExpanded(this, 6, dummy);
     925            matrix_->generalExpanded(this, 5, dummy);
     926            forceFactorization_ = 1; // a bit drastic but ..
     927            type = 2;
     928            if (internalFactorize(2) != 0) {
     929              largestPrimalError_ = 1.0e4; // force other type
     930            }
     931          }
    858932          changeMade_++; // say change made
    859      }
    860      int saveThreshold = factorization_->sparseThreshold();
    861      int tentativeStatus = problemStatus_;
    862      int numberThrownOut = 1; // to loop round on bad factorization in values pass
    863      double lastSumInfeasibility = COIN_DBL_MAX;
    864      if (numberIterations_)
    865           lastSumInfeasibility = nonLinearCost_->sumInfeasibilities();
    866      int nPass = 0;
    867      while (numberThrownOut) {
    868           int nSlackBasic = 0;
    869           if (nPass) {
    870                for (int i = 0; i < numberRows_; i++) {
    871                     if (getRowStatus(i) == basic)
    872                          nSlackBasic++;
    873                }
    874           }
    875           nPass++;
    876           if (problemStatus_ > -3 || problemStatus_ == -4) {
    877                // factorize
    878                // later on we will need to recover from singularities
    879                // also we could skip if first time
    880                // do weights
    881                // This may save pivotRow_ for use
    882                if (doFactorization)
    883                     primalColumnPivot_->saveWeights(this, 1);
    884 
    885                if ((type && doFactorization) || nSlackBasic == numberRows_) {
    886                     // is factorization okay?
    887                     int factorStatus = internalFactorize(1);
    888                     if (factorStatus) {
    889                          if (solveType_ == 2 + 8) {
    890                               // say odd
    891                               problemStatus_ = 5;
    892                               return;
    893                          }
    894                          if (type != 1 || largestPrimalError_ > 1.0e3
    895                                    || largestDualError_ > 1.0e3) {
    896                               // switch off dense
    897                               int saveDense = factorization_->denseThreshold();
    898                               factorization_->setDenseThreshold(0);
    899                               // Go to safe
    900                               factorization_->pivotTolerance(0.99);
    901                               // make sure will do safe factorization
    902                               pivotVariable_[0] = -1;
    903                               internalFactorize(2);
    904                               factorization_->setDenseThreshold(saveDense);
    905                               // restore extra stuff
    906                               matrix_->generalExpanded(this, 6, dummy);
    907                          } else {
    908                               // no - restore previous basis
    909                               // Keep any flagged variables
    910                               int i;
    911                               for (i = 0; i < numberRows_ + numberColumns_; i++) {
    912                                    if (flagged(i))
    913                                         saveStatus_[i] |= 64; //say flagged
    914                               }
    915                               CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
    916                               if (numberPivots <= 1) {
    917                                    // throw out something
    918                                    if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
    919                                         setFlagged(sequenceIn_);
    920                                    } else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
    921                                         setFlagged(sequenceOut_);
    922                                    }
    923                                    double newTolerance = CoinMax(0.5 + 0.499 * randomNumberGenerator_.randomDouble(), factorization_->pivotTolerance());
    924                                    factorization_->pivotTolerance(newTolerance);
    925                               } else {
    926                                    // Go to safe
    927                                    factorization_->pivotTolerance(0.99);
    928                               }
    929                               CoinMemcpyN(savedSolution_ + numberColumns_ ,
    930                                           numberRows_, rowActivityWork_);
    931                               CoinMemcpyN(savedSolution_ ,
    932                                           numberColumns_, columnActivityWork_);
    933                               // restore extra stuff
    934                               matrix_->generalExpanded(this, 6, dummy);
    935                               matrix_->generalExpanded(this, 5, dummy);
    936                               forceFactorization_ = 1; // a bit drastic but ..
    937                               type = 2;
    938                               if (internalFactorize(2) != 0) {
    939                                    largestPrimalError_ = 1.0e4; // force other type
    940                               }
    941                          }
    942                          changeMade_++; // say change made
    943                     }
    944                }
    945                if (problemStatus_ != -4)
    946                     problemStatus_ = -3;
    947           }
    948           if(progress->realInfeasibility_[0]<1.0e-1 &&
    949              primalTolerance_==1.0e-7&&progress->iterationNumber_[0]>0&&
    950              progress->iterationNumber_[CLP_PROGRESS-1]-progress->iterationNumber_[0]>25) {
    951             // default - so user did not set
    952             int iP;
    953             double minAverage=COIN_DBL_MAX;
    954             double maxAverage=0.0;
    955             for (iP=0;iP<CLP_PROGRESS;iP++) {
    956               int n=progress->numberInfeasibilities_[iP];
    957               if (!n) {
    958                 break;
    959               } else {
    960                 double average=progress->realInfeasibility_[iP];
    961                 if (average>0.1)
    962                   break;
    963                 average /= static_cast<double>(n);
    964                 minAverage=CoinMin(minAverage,average);
    965                 maxAverage=CoinMax(maxAverage,average);
    966               }
    967             }
    968             if (iP==CLP_PROGRESS&&minAverage<1.0e-5&&maxAverage<1.0e-3) {
    969               // change tolerance
    970 #if CBC_USEFUL_PRINTING>0
    971               printf("CCchanging tolerance\n");
    972 #endif
    973               primalTolerance_=1.0e-6;
    974               dblParam_[ClpPrimalTolerance]=1.0e-6;
    975               moreSpecialOptions_ |= 4194304;
    976             }
    977           }
    978           // at this stage status is -3 or -5 if looks unbounded
    979           // get primal and dual solutions
    980           // put back original costs and then check
    981           // createRim(4); // costs do not change
    982           // May need to do more if column generation
    983           dummy = 4;
    984           matrix_->generalExpanded(this, 9, dummy);
     933        }
     934      }
     935      if (problemStatus_ != -4)
     936        problemStatus_ = -3;
     937    }
     938    if (progress->realInfeasibility_[0] < 1.0e-1 && primalTolerance_ == 1.0e-7 && progress->iterationNumber_[0] > 0 && progress->iterationNumber_[CLP_PROGRESS - 1] - progress->iterationNumber_[0] > 25) {
     939      // default - so user did not set
     940      int iP;
     941      double minAverage = COIN_DBL_MAX;
     942      double maxAverage = 0.0;
     943      for (iP = 0; iP < CLP_PROGRESS; iP++) {
     944        int n = progress->numberInfeasibilities_[iP];
     945        if (!n) {
     946          break;
     947        } else {
     948          double average = progress->realInfeasibility_[iP];
     949          if (average > 0.1)
     950            break;
     951          average /= static_cast< double >(n);
     952          minAverage = CoinMin(minAverage, average);
     953          maxAverage = CoinMax(maxAverage, average);
     954        }
     955      }
     956      if (iP == CLP_PROGRESS && minAverage < 1.0e-5 && maxAverage < 1.0e-3) {
     957        // change tolerance
     958#if CBC_USEFUL_PRINTING > 0
     959        printf("CCchanging tolerance\n");
     960#endif
     961        primalTolerance_ = 1.0e-6;
     962        dblParam_[ClpPrimalTolerance] = 1.0e-6;
     963        moreSpecialOptions_ |= 4194304;
     964      }
     965    }
     966    // at this stage status is -3 or -5 if looks unbounded
     967    // get primal and dual solutions
     968    // put back original costs and then check
     969    // createRim(4); // costs do not change
     970    // May need to do more if column generation
     971    dummy = 4;
     972    matrix_->generalExpanded(this, 9, dummy);
    985973#ifndef CLP_CAUTION
    986974#define CLP_CAUTION 1
    987975#endif
    988976#if CLP_CAUTION
    989           double lastAverageInfeasibility = sumDualInfeasibilities_ /
    990                                             static_cast<double>(numberDualInfeasibilities_ + 10);
     977    double lastAverageInfeasibility = sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_ + 10);
    991978#endif
    992979#ifdef CLP_USER_DRIVEN
    993           int status = eventHandler_->event(ClpEventHandler::goodFactorization);
    994           if (status >= 0) {
    995             lastSumInfeasibility = COIN_DBL_MAX;
    996           }
    997 #endif
    998           if (ifValuesPass && firstFree_ <0) {
    999             largestPrimalError_=1.0e7;
    1000             largestDualError_=1.0e7;
    1001           }
    1002           bool sayValuesPass = (firstFree_>=0);
    1003           //if (ifValuesPass&&numberIterations_==baseIteration_)
    1004           //sayValuesPass=true;
    1005           numberThrownOut = gutsOfSolution(NULL, NULL, sayValuesPass);
    1006           double sumInfeasibility = nonLinearCost_->sumInfeasibilities();
    1007           int reason2 = 0;
     980    int status = eventHandler_->event(ClpEventHandler::goodFactorization);
     981    if (status >= 0) {
     982      lastSumInfeasibility = COIN_DBL_MAX;
     983    }
     984#endif
     985    if (ifValuesPass && firstFree_ < 0) {
     986      largestPrimalError_ = 1.0e7;
     987      largestDualError_ = 1.0e7;
     988    }
     989    bool sayValuesPass = (firstFree_ >= 0);
     990    //if (ifValuesPass&&numberIterations_==baseIteration_)
     991    //sayValuesPass=true;
     992    numberThrownOut = gutsOfSolution(NULL, NULL, sayValuesPass);
     993    double sumInfeasibility = nonLinearCost_->sumInfeasibilities();
     994    int reason2 = 0;
    1008995#if CLP_CAUTION
    1009 #if CLP_CAUTION==2
    1010           double test2 = 1.0e5;
     996#if CLP_CAUTION == 2
     997    double test2 = 1.0e5;
    1011998#else
    1012           double test2 = 1.0e-1;
    1013 #endif
    1014           if (!lastSumInfeasibility && sumInfeasibility &&
    1015                     lastAverageInfeasibility < test2 && numberPivots > 10)
    1016                reason2 = 3;
    1017           if (lastSumInfeasibility < 1.0e-6 && sumInfeasibility > 1.0e-3 &&
    1018                     numberPivots > 10)
    1019                reason2 = 4;
    1020 #endif
    1021           if (numberThrownOut)
    1022                reason2 = 1;
    1023           if ((sumInfeasibility > 1.0e7 && sumInfeasibility > 100.0 * lastSumInfeasibility
    1024                     && factorization_->pivotTolerance() < 0.11) ||
    1025                     (largestPrimalError_ > 1.0e10 && largestDualError_ > 1.0e10))
    1026                reason2 = 2;
    1027           if (reason2) {
    1028                problemStatus_ = tentativeStatus;
    1029                doFactorization = true;
    1030                if (numberPivots||(numberThrownOut==-123456789&&saveStatus_)) {
    1031                     // go back
    1032                     // trouble - restore solution
    1033                     CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
    1034                     CoinMemcpyN(savedSolution_ + numberColumns_ ,
    1035                                 numberRows_, rowActivityWork_);
    1036                     CoinMemcpyN(savedSolution_ ,
    1037                                 numberColumns_, columnActivityWork_);
    1038                     // restore extra stuff
    1039                     matrix_->generalExpanded(this, 6, dummy);
    1040                     if (reason2 < 3) {
    1041                          // Go to safe
    1042                          factorization_->pivotTolerance(CoinMin(0.99, 1.01 * factorization_->pivotTolerance()));
    1043                          forceFactorization_ = 1; // a bit drastic but ..
    1044                     } else if (forceFactorization_ < 0) {
    1045                          forceFactorization_ = CoinMin(numberPivots / 2, 100);
    1046                     } else {
    1047                          forceFactorization_ = CoinMin(forceFactorization_,
    1048                                                        CoinMax(3, numberPivots / 2));
    1049                     }
    1050                     pivotRow_ = -1; // say no weights update
    1051                     changeMade_++; // say change made
    1052                     if (numberPivots == 1) {
    1053                          // throw out something
    1054                          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
    1055                               setFlagged(sequenceIn_);
    1056                          } else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
    1057                               setFlagged(sequenceOut_);
    1058                          }
    1059                     }
    1060                     type = 2; // so will restore weights
    1061                     if (internalFactorize(2) != 0) {
    1062                          largestPrimalError_ = 1.0e4; // force other type
    1063                     }
    1064                     numberPivots = 0;
    1065                     numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
    1066                     //assert (!numberThrownOut);
     999    double test2 = 1.0e-1;
     1000#endif
     1001    if (!lastSumInfeasibility && sumInfeasibility && lastAverageInfeasibility < test2 && numberPivots > 10)
     1002      reason2 = 3;
     1003    if (lastSumInfeasibility < 1.0e-6 && sumInfeasibility > 1.0e-3 && numberPivots > 10)
     1004      reason2 = 4;
     1005#endif
     1006    if (numberThrownOut)
     1007      reason2 = 1;
     1008    if ((sumInfeasibility > 1.0e7 && sumInfeasibility > 100.0 * lastSumInfeasibility
     1009          && factorization_->pivotTolerance() < 0.11)
     1010      || (largestPrimalError_ > 1.0e10 && largestDualError_ > 1.0e10))
     1011      reason2 = 2;
     1012    if (reason2) {
     1013      problemStatus_ = tentativeStatus;
     1014      doFactorization = true;
     1015      if (numberPivots || (numberThrownOut == -123456789 && saveStatus_)) {
     1016        // go back
     1017        // trouble - restore solution
     1018        CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
     1019        CoinMemcpyN(savedSolution_ + numberColumns_,
     1020          numberRows_, rowActivityWork_);
     1021        CoinMemcpyN(savedSolution_,
     1022          numberColumns_, columnActivityWork_);
     1023        // restore extra stuff
     1024        matrix_->generalExpanded(this, 6, dummy);
     1025        if (reason2 < 3) {
     1026          // Go to safe
     1027          factorization_->pivotTolerance(CoinMin(0.99, 1.01 * factorization_->pivotTolerance()));
     1028          forceFactorization_ = 1; // a bit drastic but ..
     1029        } else if (forceFactorization_ < 0) {
     1030          forceFactorization_ = CoinMin(numberPivots / 2, 100);
     1031        } else {
     1032          forceFactorization_ = CoinMin(forceFactorization_,
     1033            CoinMax(3, numberPivots / 2));
     1034        }
     1035        pivotRow_ = -1; // say no weights update
     1036        changeMade_++; // say change made
     1037        if (numberPivots == 1) {
     1038          // throw out something
     1039          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
     1040            setFlagged(sequenceIn_);
     1041          } else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
     1042            setFlagged(sequenceOut_);
     1043          }
     1044        }
     1045        type = 2; // so will restore weights
     1046        if (internalFactorize(2) != 0) {
     1047          largestPrimalError_ = 1.0e4; // force other type
     1048        }
     1049        numberPivots = 0;
     1050        numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
     1051        //assert (!numberThrownOut);
    10671052#ifdef CLP_USEFUL_PRINTOUT
    1068                     if (numberThrownOut)
    1069                       printf("OUCH! - %d thrownout at %s line %d\n",
    1070                              numberThrownOut,__FILE__,__LINE__);
    1071 #endif
    1072                     sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
    1073                }
    1074           }
    1075      }
    1076      // Double check reduced costs if no action
    1077      if (progress->lastIterationNumber(0) == numberIterations_) {
    1078           if (primalColumnPivot_->looksOptimal()) {
    1079                numberDualInfeasibilities_ = 0;
    1080                sumDualInfeasibilities_ = 0.0;
    1081           }
    1082      }
    1083      // If in primal and small dj give up
    1084      if ((specialOptions_ & 1024) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) {
    1085           double average = sumDualInfeasibilities_ / (static_cast<double> (numberDualInfeasibilities_));
    1086           if (numberIterations_ > 300 && average < 1.0e-4) {
    1087                numberDualInfeasibilities_ = 0;
    1088                sumDualInfeasibilities_ = 0.0;
    1089           }
    1090      }
    1091      // Check if looping
    1092      int loop;
    1093      if (type != 2 && !ifValuesPass)
    1094           loop = progress->looping();
    1095      else
    1096           loop = -1;
    1097      if (loop >= 0) {
    1098           if (!problemStatus_) {
    1099                // declaring victory
    1100                numberPrimalInfeasibilities_ = 0;
    1101                sumPrimalInfeasibilities_ = 0.0;
     1053        if (numberThrownOut)
     1054          printf("OUCH! - %d thrownout at %s line %d\n",
     1055            numberThrownOut, __FILE__, __LINE__);
     1056#endif
     1057        sumInfeasibility = nonLinearCost_->sumInfeasibilities();
     1058      }
     1059    }
     1060  }
     1061  // Double check reduced costs if no action
     1062  if (progress->lastIterationNumber(0) == numberIterations_) {
     1063    if (primalColumnPivot_->looksOptimal()) {
     1064      numberDualInfeasibilities_ = 0;
     1065      sumDualInfeasibilities_ = 0.0;
     1066    }
     1067  }
     1068  // If in primal and small dj give up
     1069  if ((specialOptions_ & 1024) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) {
     1070    double average = sumDualInfeasibilities_ / (static_cast< double >(numberDualInfeasibilities_));
     1071    if (numberIterations_ > 300 && average < 1.0e-4) {
     1072      numberDualInfeasibilities_ = 0;
     1073      sumDualInfeasibilities_ = 0.0;
     1074    }
     1075  }
     1076  // Check if looping
     1077  int loop;
     1078  if (type != 2 && !ifValuesPass)
     1079    loop = progress->looping();
     1080  else
     1081    loop = -1;
     1082  if (loop >= 0) {
     1083    if (!problemStatus_) {
     1084      // declaring victory
     1085      numberPrimalInfeasibilities_ = 0;
     1086      sumPrimalInfeasibilities_ = 0.0;
     1087    } else {
     1088      problemStatus_ = loop; //exit if in loop
     1089      problemStatus_ = 10; // instead - try other algorithm
     1090      numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
     1091    }
     1092    problemStatus_ = 10; // instead - try other algorithm
     1093    return;
     1094  } else if (loop < -1) {
     1095    // Is it time for drastic measures
     1096    if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 && progress->oddState() < 10 && progress->oddState() >= 0) {
     1097      progress->newOddState();
     1098      nonLinearCost_->zapCosts();
     1099    }
     1100    // something may have changed
     1101    gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1102  }
     1103  // If progress then reset costs
     1104  if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
     1105    createRim(4, false); // costs back
     1106    delete nonLinearCost_;
     1107    nonLinearCost_ = new ClpNonLinearCost(this);
     1108    progress->endOddState();
     1109    gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1110  }
     1111  // Flag to say whether to go to dual to clean up
     1112  bool goToDual = false;
     1113  // really for free variables in
     1114  //if((progressFlag_&2)!=0)
     1115  //problemStatus_=-1;;
     1116  progressFlag_ = 0; //reset progress flag
     1117
     1118  handler_->message(CLP_SIMPLEX_STATUS, messages_)
     1119    << numberIterations_ << nonLinearCost_->feasibleReportCost();
     1120  handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
     1121    << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
     1122  handler_->printing(sumDualInfeasibilities_ > 0.0)
     1123    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
     1124  handler_->printing(numberDualInfeasibilitiesWithoutFree_
     1125    < numberDualInfeasibilities_)
     1126    << numberDualInfeasibilitiesWithoutFree_;
     1127  handler_->message() << CoinMessageEol;
     1128  if (!primalFeasible()) {
     1129    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1130    gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1131    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1132  }
     1133  if (nonLinearCost_->numberInfeasibilities() > 0 && !progress->initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10) {
     1134    // first time infeasible - start up weight computation
     1135    double *oldDj = dj_;
     1136    double *oldCost = cost_;
     1137    int numberRows2 = numberRows_ + numberExtraRows_;
     1138    int numberTotal = numberRows2 + numberColumns_;
     1139    dj_ = new double[numberTotal];
     1140    cost_ = new double[numberTotal];
     1141    reducedCostWork_ = dj_;
     1142    rowReducedCost_ = dj_ + numberColumns_;
     1143    objectiveWork_ = cost_;
     1144    rowObjectiveWork_ = cost_ + numberColumns_;
     1145    double direction = optimizationDirection_ * objectiveScale_;
     1146    const double *obj = objective();
     1147    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
     1148    int iSequence;
     1149    if (columnScale_)
     1150      for (iSequence = 0; iSequence < numberColumns_; iSequence++)
     1151        cost_[iSequence] = obj[iSequence] * direction * columnScale_[iSequence];
     1152    else
     1153      for (iSequence = 0; iSequence < numberColumns_; iSequence++)
     1154        cost_[iSequence] = obj[iSequence] * direction;
     1155    computeDuals(NULL);
     1156    int numberSame = 0;
     1157    int numberDifferent = 0;
     1158    int numberZero = 0;
     1159    int numberFreeSame = 0;
     1160    int numberFreeDifferent = 0;
     1161    int numberFreeZero = 0;
     1162    int n = 0;
     1163    for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     1164      if (getStatus(iSequence) != basic && !flagged(iSequence)) {
     1165        // not basic
     1166        double distanceUp = upper_[iSequence] - solution_[iSequence];
     1167        double distanceDown = solution_[iSequence] - lower_[iSequence];
     1168        double feasibleDj = dj_[iSequence];
     1169        double infeasibleDj = oldDj[iSequence] - feasibleDj;
     1170        double value = feasibleDj * infeasibleDj;
     1171        if (distanceUp > primalTolerance_) {
     1172          // Check if "free"
     1173          if (distanceDown > primalTolerance_) {
     1174            // free
     1175            if (value > dualTolerance_) {
     1176              numberFreeSame++;
     1177            } else if (value < -dualTolerance_) {
     1178              numberFreeDifferent++;
     1179              dj_[n++] = feasibleDj / infeasibleDj;
     1180            } else {
     1181              numberFreeZero++;
     1182            }
    11021183          } else {
    1103                problemStatus_ = loop; //exit if in loop
    1104                problemStatus_ = 10; // instead - try other algorithm
    1105                numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
    1106           }
    1107           problemStatus_ = 10; // instead - try other algorithm
    1108           return ;
    1109      } else if (loop < -1) {
    1110           // Is it time for drastic measures
    1111           if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 &&
    1112                     progress->oddState() < 10 && progress->oddState() >= 0) {
    1113                progress->newOddState();
    1114                nonLinearCost_->zapCosts();
    1115           }
    1116           // something may have changed
    1117           gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1118      }
    1119      // If progress then reset costs
    1120      if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
    1121           createRim(4, false); // costs back
    1122           delete nonLinearCost_;
    1123           nonLinearCost_ = new ClpNonLinearCost(this);
    1124           progress->endOddState();
    1125           gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1126      }
    1127      // Flag to say whether to go to dual to clean up
    1128      bool goToDual = false;
    1129      // really for free variables in
    1130      //if((progressFlag_&2)!=0)
    1131      //problemStatus_=-1;;
    1132      progressFlag_ = 0; //reset progress flag
    1133 
    1134      handler_->message(CLP_SIMPLEX_STATUS, messages_)
    1135                << numberIterations_ << nonLinearCost_->feasibleReportCost();
    1136      handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
    1137                << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
    1138      handler_->printing(sumDualInfeasibilities_ > 0.0)
    1139                << sumDualInfeasibilities_ << numberDualInfeasibilities_;
    1140      handler_->printing(numberDualInfeasibilitiesWithoutFree_
    1141                         < numberDualInfeasibilities_)
    1142                << numberDualInfeasibilitiesWithoutFree_;
    1143      handler_->message() << CoinMessageEol;
    1144      if (!primalFeasible()) {
    1145           nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1146           gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1147           nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1148      }
    1149      if (nonLinearCost_->numberInfeasibilities() > 0 && !progress->initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10) {
    1150           // first time infeasible - start up weight computation
    1151           double * oldDj = dj_;
    1152           double * oldCost = cost_;
    1153           int numberRows2 = numberRows_ + numberExtraRows_;
    1154           int numberTotal = numberRows2 + numberColumns_;
    1155           dj_ = new double[numberTotal];
    1156           cost_ = new double[numberTotal];
    1157           reducedCostWork_ = dj_;
    1158           rowReducedCost_ = dj_ + numberColumns_;
    1159           objectiveWork_ = cost_;
    1160           rowObjectiveWork_ = cost_ + numberColumns_;
    1161           double direction = optimizationDirection_ * objectiveScale_;
    1162           const double * obj = objective();
    1163           memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
    1164           int iSequence;
    1165           if (columnScale_)
    1166                for (iSequence = 0; iSequence < numberColumns_; iSequence++)
    1167                     cost_[iSequence] = obj[iSequence] * direction * columnScale_[iSequence];
    1168           else
    1169                for (iSequence = 0; iSequence < numberColumns_; iSequence++)
    1170                     cost_[iSequence] = obj[iSequence] * direction;
    1171           computeDuals(NULL);
    1172           int numberSame = 0;
    1173           int numberDifferent = 0;
    1174           int numberZero = 0;
    1175           int numberFreeSame = 0;
    1176           int numberFreeDifferent = 0;
    1177           int numberFreeZero = 0;
    1178           int n = 0;
    1179           for (iSequence = 0; iSequence < numberTotal; iSequence++) {
    1180                if (getStatus(iSequence) != basic && !flagged(iSequence)) {
    1181                     // not basic
    1182                     double distanceUp = upper_[iSequence] - solution_[iSequence];
    1183                     double distanceDown = solution_[iSequence] - lower_[iSequence];
    1184                     double feasibleDj = dj_[iSequence];
    1185                     double infeasibleDj = oldDj[iSequence] - feasibleDj;
    1186                     double value = feasibleDj * infeasibleDj;
    1187                     if (distanceUp > primalTolerance_) {
    1188                          // Check if "free"
    1189                          if (distanceDown > primalTolerance_) {
    1190                               // free
    1191                               if (value > dualTolerance_) {
    1192                                    numberFreeSame++;
    1193                               } else if(value < -dualTolerance_) {
    1194                                    numberFreeDifferent++;
    1195                                    dj_[n++] = feasibleDj / infeasibleDj;
    1196                               } else {
    1197                                    numberFreeZero++;
    1198                               }
    1199                          } else {
    1200                               // should not be negative
    1201                               if (value > dualTolerance_) {
    1202                                    numberSame++;
    1203                               } else if(value < -dualTolerance_) {
    1204                                    numberDifferent++;
    1205                                    dj_[n++] = feasibleDj / infeasibleDj;
    1206                               } else {
    1207                                    numberZero++;
    1208                               }
    1209                          }
    1210                     } else if (distanceDown > primalTolerance_) {
    1211                          // should not be positive
    1212                          if (value > dualTolerance_) {
    1213                               numberSame++;
    1214                          } else if(value < -dualTolerance_) {
    1215                               numberDifferent++;
    1216                               dj_[n++] = feasibleDj / infeasibleDj;
    1217                          } else {
    1218                               numberZero++;
    1219                          }
    1220                     }
    1221                }
    1222                progress->initialWeight_ = -1.0;
    1223           }
     1184            // should not be negative
     1185            if (value > dualTolerance_) {
     1186              numberSame++;
     1187            } else if (value < -dualTolerance_) {
     1188              numberDifferent++;
     1189              dj_[n++] = feasibleDj / infeasibleDj;
     1190            } else {
     1191              numberZero++;
     1192            }
     1193          }
     1194        } else if (distanceDown > primalTolerance_) {
     1195          // should not be positive
     1196          if (value > dualTolerance_) {
     1197            numberSame++;
     1198          } else if (value < -dualTolerance_) {
     1199            numberDifferent++;
     1200            dj_[n++] = feasibleDj / infeasibleDj;
     1201          } else {
     1202            numberZero++;
     1203          }
     1204        }
     1205      }
     1206      progress->initialWeight_ = -1.0;
     1207    }
    12241208#ifdef CLP_USEFUL_PRINTOUT
    1225           printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
    1226              numberSame,numberDifferent,numberZero,
    1227              numberFreeSame,numberFreeDifferent,numberFreeZero);
    1228 #endif
    1229           // we want most to be same
    1230           if (n) {
    1231                std::sort(dj_, dj_ + n);
     1209    printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
     1210      numberSame, numberDifferent, numberZero,
     1211      numberFreeSame, numberFreeDifferent, numberFreeZero);
     1212#endif
     1213    // we want most to be same
     1214    if (n) {
     1215      std::sort(dj_, dj_ + n);
    12321216#ifdef CLP_USEFUL_PRINTOUT
    1233                double most = 0.95;
    1234                int which = static_cast<int> ((1.0 - most) * static_cast<double> (n));
    1235                double take2 = -dj_[which] * infeasibilityCost_;
    1236                printf("XXXX inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take2,-dj_[0]*infeasibilityCost_,-dj_[n-1]*infeasibilityCost_);
    1237 #endif
    1238                double take = -dj_[0] * infeasibilityCost_;
    1239                // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);
    1240                infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10);
     1217      double most = 0.95;
     1218      int which = static_cast< int >((1.0 - most) * static_cast< double >(n));
     1219      double take2 = -dj_[which] * infeasibilityCost_;
     1220      printf("XXXX inf cost %g take %g (range %g %g)\n", infeasibilityCost_, take2, -dj_[0] * infeasibilityCost_, -dj_[n - 1] * infeasibilityCost_);
     1221#endif
     1222      double take = -dj_[0] * infeasibilityCost_;
     1223      // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);
     1224      infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10);
    12411225#ifdef CLP_USEFUL_PRINTOUT
    1242                printf("XXXX changing weight to %g\n",infeasibilityCost_);
    1243 #endif
    1244           }
    1245           delete [] dj_;
    1246           delete [] cost_;
    1247           dj_ = oldDj;
    1248           cost_ = oldCost;
    1249           reducedCostWork_ = dj_;
    1250           rowReducedCost_ = dj_ + numberColumns_;
    1251           objectiveWork_ = cost_;
    1252           rowObjectiveWork_ = cost_ + numberColumns_;
    1253           if (n||matrix_->type()>=15)
    1254                gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1255      }
    1256      double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
    1257      if (!nonLinearCost_->numberInfeasibilities() && infeasibilityCost_ == 1.0e10 && !ifValuesPass && true) {
    1258           // relax if default
    1259           infeasibilityCost_ = CoinMin(CoinMax(100.0 * sumDualInfeasibilities_, 1.0e8), 1.00000001e10);
    1260           // reset looping criterion
    1261           progress->reset();
    1262           trueInfeasibility = 1.123456e10;
    1263      }
    1264      if (trueInfeasibility > 1.0) {
    1265           // If infeasibility going up may change weights
    1266           double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
    1267           double lastInf = progress->lastInfeasibility(1);
    1268           double lastInf3 = progress->lastInfeasibility(3);
    1269           double thisObj = progress->lastObjective(0);
    1270           double thisInf = progress->lastInfeasibility(0);
    1271           thisObj += infeasibilityCost_ * 2.0 * thisInf;
    1272           double lastObj = progress->lastObjective(1);
    1273           lastObj += infeasibilityCost_ * 2.0 * lastInf;
    1274           double lastObj3 = progress->lastObjective(3);
    1275           lastObj3 += infeasibilityCost_ * 2.0 * lastInf3;
    1276           assert (thisObj<=COIN_DBL_MAX);
    1277           assert (lastObj<=COIN_DBL_MAX);
    1278           assert (lastObj3<=COIN_DBL_MAX);
    1279           if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-7
    1280                     && firstFree_ < 0 && thisInf >= lastInf) {
    1281                if (handler_->logLevel() == 63)
    1282                     printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_);
    1283                int maxFactor = factorization_->maximumPivots();
    1284                if (maxFactor > 10) {
    1285                     if (forceFactorization_ < 0)
    1286                          forceFactorization_ = maxFactor;
    1287                     forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2));
    1288                     if (handler_->logLevel() == 63)
    1289                          printf("Reducing factorization frequency to %d\n", forceFactorization_);
    1290                }
    1291           } else if (lastObj3 < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj3)) - 1.0e-7
    1292                      && firstFree_ < 0 && thisInf >= lastInf) {
    1293                if (handler_->logLevel() == 63)
    1294                     printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_);
    1295                int maxFactor = factorization_->maximumPivots();
    1296                if (maxFactor > 10) {
    1297                     if (forceFactorization_ < 0)
    1298                          forceFactorization_ = maxFactor;
    1299                     forceFactorization_ = CoinMax(1, (forceFactorization_ * 2) / 3);
    1300                     if (handler_->logLevel() == 63)
    1301                          printf("Reducing factorization frequency to %d\n", forceFactorization_);
    1302                }
    1303           } else if(lastInf < testValue || trueInfeasibility == 1.123456e10) {
    1304                if (infeasibilityCost_ < 1.0e14) {
    1305                     infeasibilityCost_ *= 1.5;
    1306                     // reset looping criterion
    1307                     progress->reset();
    1308                     if (handler_->logLevel() == 63)
    1309                          printf("increasing weight to %g\n", infeasibilityCost_);
    1310                     gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1311                }
    1312           }
    1313      }
    1314      // we may wish to say it is optimal even if infeasible
    1315      bool alwaysOptimal = (specialOptions_ & 1) != 0;
    1316      // give code benefit of doubt
    1317      if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
    1318                sumOfRelaxedPrimalInfeasibilities_ == 0.0 &&
    1319          progress->objective_[CLP_PROGRESS-1]>=
    1320          progress->objective_[CLP_PROGRESS-2]-1.0e-9*(10.0+fabs(objectiveValue_))) {
    1321           // say optimal (with these bounds etc)
    1322           numberDualInfeasibilities_ = 0;
    1323           sumDualInfeasibilities_ = 0.0;
    1324           numberPrimalInfeasibilities_ = 0;
    1325           sumPrimalInfeasibilities_ = 0.0;
    1326           // But check if in sprint
    1327           if (originalModel) {
    1328                // Carry on and re-do
    1329                numberDualInfeasibilities_ = -776;
    1330           }
    1331           // But if real primal infeasibilities nonzero carry on
    1332           if (nonLinearCost_->numberInfeasibilities()) {
    1333                // most likely to happen if infeasible
    1334                double relaxedToleranceP = primalTolerance_;
    1335                // we can't really trust infeasibilities if there is primal error
    1336                double error = CoinMin(1.0e-2, largestPrimalError_);
    1337                // allow tolerance at least slightly bigger than standard
    1338                relaxedToleranceP = relaxedToleranceP +  error;
    1339                int ninfeas = nonLinearCost_->numberInfeasibilities();
    1340                double sum = nonLinearCost_->sumInfeasibilities();
    1341                double average = sum / static_cast<double> (ninfeas);
     1226      printf("XXXX changing weight to %g\n", infeasibilityCost_);
     1227#endif
     1228    }
     1229    delete[] dj_;
     1230    delete[] cost_;
     1231    dj_ = oldDj;
     1232    cost_ = oldCost;
     1233    reducedCostWork_ = dj_;
     1234    rowReducedCost_ = dj_ + numberColumns_;
     1235    objectiveWork_ = cost_;
     1236    rowObjectiveWork_ = cost_ + numberColumns_;
     1237    if (n || matrix_->type() >= 15)
     1238      gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1239  }
     1240  double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
     1241  if (!nonLinearCost_->numberInfeasibilities() && infeasibilityCost_ == 1.0e10 && !ifValuesPass && true) {
     1242    // relax if default
     1243    infeasibilityCost_ = CoinMin(CoinMax(100.0 * sumDualInfeasibilities_, 1.0e8), 1.00000001e10);
     1244    // reset looping criterion
     1245    progress->reset();
     1246    trueInfeasibility = 1.123456e10;
     1247  }
     1248  if (trueInfeasibility > 1.0) {
     1249    // If infeasibility going up may change weights
     1250    double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
     1251    double lastInf = progress->lastInfeasibility(1);
     1252    double lastInf3 = progress->lastInfeasibility(3);
     1253    double thisObj = progress->lastObjective(0);
     1254    double thisInf = progress->lastInfeasibility(0);
     1255    thisObj += infeasibilityCost_ * 2.0 * thisInf;
     1256    double lastObj = progress->lastObjective(1);
     1257    lastObj += infeasibilityCost_ * 2.0 * lastInf;
     1258    double lastObj3 = progress->lastObjective(3);
     1259    lastObj3 += infeasibilityCost_ * 2.0 * lastInf3;
     1260    assert(thisObj <= COIN_DBL_MAX);
     1261    assert(lastObj <= COIN_DBL_MAX);
     1262    assert(lastObj3 <= COIN_DBL_MAX);
     1263    if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-7
     1264      && firstFree_ < 0 && thisInf >= lastInf) {
     1265      if (handler_->logLevel() == 63)
     1266        printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_);
     1267      int maxFactor = factorization_->maximumPivots();
     1268      if (maxFactor > 10) {
     1269        if (forceFactorization_ < 0)
     1270          forceFactorization_ = maxFactor;
     1271        forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2));
     1272        if (handler_->logLevel() == 63)
     1273          printf("Reducing factorization frequency to %d\n", forceFactorization_);
     1274      }
     1275    } else if (lastObj3 < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj3)) - 1.0e-7
     1276      && firstFree_ < 0 && thisInf >= lastInf) {
     1277      if (handler_->logLevel() == 63)
     1278        printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_);
     1279      int maxFactor = factorization_->maximumPivots();
     1280      if (maxFactor > 10) {
     1281        if (forceFactorization_ < 0)
     1282          forceFactorization_ = maxFactor;
     1283        forceFactorization_ = CoinMax(1, (forceFactorization_ * 2) / 3);
     1284        if (handler_->logLevel() == 63)
     1285          printf("Reducing factorization frequency to %d\n", forceFactorization_);
     1286      }
     1287    } else if (lastInf < testValue || trueInfeasibility == 1.123456e10) {
     1288      if (infeasibilityCost_ < 1.0e14) {
     1289        infeasibilityCost_ *= 1.5;
     1290        // reset looping criterion
     1291        progress->reset();
     1292        if (handler_->logLevel() == 63)
     1293          printf("increasing weight to %g\n", infeasibilityCost_);
     1294        gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1295      }
     1296    }
     1297  }
     1298  // we may wish to say it is optimal even if infeasible
     1299  bool alwaysOptimal = (specialOptions_ & 1) != 0;
     1300  // give code benefit of doubt
     1301  if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0 && progress->objective_[CLP_PROGRESS - 1] >= progress->objective_[CLP_PROGRESS - 2] - 1.0e-9 * (10.0 + fabs(objectiveValue_))) {
     1302    // say optimal (with these bounds etc)
     1303    numberDualInfeasibilities_ = 0;
     1304    sumDualInfeasibilities_ = 0.0;
     1305    numberPrimalInfeasibilities_ = 0;
     1306    sumPrimalInfeasibilities_ = 0.0;
     1307    // But check if in sprint
     1308    if (originalModel) {
     1309      // Carry on and re-do
     1310      numberDualInfeasibilities_ = -776;
     1311    }
     1312    // But if real primal infeasibilities nonzero carry on
     1313    if (nonLinearCost_->numberInfeasibilities()) {
     1314      // most likely to happen if infeasible
     1315      double relaxedToleranceP = primalTolerance_;
     1316      // we can't really trust infeasibilities if there is primal error
     1317      double error = CoinMin(1.0e-2, largestPrimalError_);
     1318      // allow tolerance at least slightly bigger than standard
     1319      relaxedToleranceP = relaxedToleranceP + error;
     1320      int ninfeas = nonLinearCost_->numberInfeasibilities();
     1321      double sum = nonLinearCost_->sumInfeasibilities();
     1322      double average = sum / static_cast< double >(ninfeas);
    13421323#ifdef COIN_DEVELOP
    1343                if (handler_->logLevel() > 0)
    1344                     printf("nonLinearCost says infeasible %d summing to %g\n",
    1345                            ninfeas, sum);
    1346 #endif
    1347                if (average > relaxedToleranceP) {
    1348                     sumOfRelaxedPrimalInfeasibilities_ = sum;
    1349                     numberPrimalInfeasibilities_ = ninfeas;
    1350                     sumPrimalInfeasibilities_ = sum;
     1324      if (handler_->logLevel() > 0)
     1325        printf("nonLinearCost says infeasible %d summing to %g\n",
     1326          ninfeas, sum);
     1327#endif
     1328      if (average > relaxedToleranceP) {
     1329        sumOfRelaxedPrimalInfeasibilities_ = sum;
     1330        numberPrimalInfeasibilities_ = ninfeas;
     1331        sumPrimalInfeasibilities_ = sum;
    13511332#ifdef COIN_DEVELOP
    1352                     bool unflagged =
    1353 #endif
    1354                          unflag();
     1333        bool unflagged =
     1334#endif
     1335          unflag();
    13551336#ifdef COIN_DEVELOP
    1356                     if (unflagged && handler_->logLevel() > 0)
    1357                          printf(" - but flagged variables\n");
    1358 #endif
    1359                }
    1360           }
    1361      }
    1362      bool looksOptimal = (!numberDualInfeasibilities_&&!nonLinearCost_->sumInfeasibilities());
    1363      // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
    1364      if ((dualFeasible() || problemStatus_ == -4) && (!ifValuesPass||looksOptimal||firstFree_<0)) {
    1365           // see if extra helps
    1366           if (nonLinearCost_->numberInfeasibilities() &&
    1367                     (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
    1368               /*&& alwaysOptimal*/) {
    1369                //may need infeasiblity cost changed
    1370                // we can see if we can construct a ray
    1371                // make up a new objective
    1372                double saveWeight = infeasibilityCost_;
    1373                // save nonlinear cost as we are going to switch off costs
    1374                //ClpNonLinearCost * nonLinear = nonLinearCost_;
    1375                // do twice to make sure Primal solution has settled
    1376                // put non-basics to bounds in case tolerance moved
    1377                // put back original costs
    1378                createRim(4);
    1379                nonLinearCost_->checkInfeasibilities(0.0);
    1380                gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1381 
    1382                //infeasibilityCost_ = 1.0e100;
    1383                // put back original costs
    1384                createRim(4);
    1385                nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1386                // may have fixed infeasibilities - double check
    1387                if (nonLinearCost_->numberInfeasibilities() == 0) {
    1388                     // carry on
    1389                     problemStatus_ = -1;
    1390                     infeasibilityCost_ = saveWeight;
    1391                     nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1392                } else {
     1337        if (unflagged && handler_->logLevel() > 0)
     1338          printf(" - but flagged variables\n");
     1339#endif
     1340      }
     1341    }
     1342  }
     1343  bool looksOptimal = (!numberDualInfeasibilities_ && !nonLinearCost_->sumInfeasibilities());
     1344  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
     1345  if ((dualFeasible() || problemStatus_ == -4) && (!ifValuesPass || looksOptimal || firstFree_ < 0)) {
     1346    // see if extra helps
     1347    if (nonLinearCost_->numberInfeasibilities() && (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
     1348      /*&& alwaysOptimal*/) {
     1349      //may need infeasiblity cost changed
     1350      // we can see if we can construct a ray
     1351      // make up a new objective
     1352      double saveWeight = infeasibilityCost_;
     1353      // save nonlinear cost as we are going to switch off costs
     1354      //ClpNonLinearCost * nonLinear = nonLinearCost_;
     1355      // do twice to make sure Primal solution has settled
     1356      // put non-basics to bounds in case tolerance moved
     1357      // put back original costs
     1358      createRim(4);
     1359      nonLinearCost_->checkInfeasibilities(0.0);
     1360      gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1361
     1362      //infeasibilityCost_ = 1.0e100;
     1363      // put back original costs
     1364      createRim(4);
     1365      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1366      // may have fixed infeasibilities - double check
     1367      if (nonLinearCost_->numberInfeasibilities() == 0) {
     1368        // carry on
     1369        problemStatus_ = -1;
     1370        infeasibilityCost_ = saveWeight;
     1371        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1372      } else {
    13931373#ifndef MAX_INFEASIBILITY_COST
    13941374#define MAX_INFEASIBILITY_COST 1.0e18
    13951375#endif
    1396                     infeasibilityCost_ = 1.0e30;
    1397                     gutsOfSolution(NULL, NULL, ifValuesPass != 0 && firstFree_>=0);
    1398                     infeasibilityCost_ = CoinMax(saveWeight,saveOriginalWeight);
    1399                     if ((infeasibilityCost_ >= MAX_INFEASIBILITY_COST ||
    1400                               numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
    1401                          goToDual = unPerturb(); // stop any further perturbation
    1402                          if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
    1403                            goToDual = false;
    1404                          nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1405                          numberDualInfeasibilities_ = 1; // carry on
    1406                          problemStatus_ = -1;
    1407                     } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e-2 &&
    1408                                (moreSpecialOptions_ & (256|8192)) == 0) {
    1409                          goToDual = true;
    1410                          factorization_->pivotTolerance(CoinMax(0.9, factorization_->pivotTolerance()));
    1411                     }
    1412                     if (!goToDual) {
    1413                          if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST ||
    1414                                    numberDualInfeasibilities_ == 0) {
    1415                               // we are infeasible - use as ray
    1416                               delete [] ray_;
    1417                               //if ((specialOptions_&(32|0x01000000))!=0x01000000) {
    1418                               if ((specialOptions_&0x03000000)==0) {
    1419                                 ray_ = new double [numberRows_];
    1420                                 double * saveObjective = CoinCopyOfArray(objective(),
    1421                                                                          numberColumns_);
    1422                                 memset(objective(),0,numberColumns_*sizeof(double));
    1423                                 infeasibilityCost_ = 1.0;
    1424                                 createRim(4);
    1425                                 nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1426                                 gutsOfSolution(NULL, NULL, false);
    1427                                 memcpy(objective(),saveObjective,numberColumns_*sizeof(double));
    1428                                 delete [] saveObjective;
    1429                                 // swap sign
    1430                                 for (int i=0;i<numberRows_;i++)
    1431                                   ray_[i] = -dual_[i];
    1432                               } else {
    1433                                 ray_ = NULL;
    1434                               }
     1376        infeasibilityCost_ = 1.0e30;
     1377        gutsOfSolution(NULL, NULL, ifValuesPass != 0 && firstFree_ >= 0);
     1378        infeasibilityCost_ = CoinMax(saveWeight, saveOriginalWeight);
     1379        if ((infeasibilityCost_ >= MAX_INFEASIBILITY_COST || numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
     1380          goToDual = unPerturb(); // stop any further perturbation
     1381          if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
     1382            goToDual = false;
     1383          nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1384          numberDualInfeasibilities_ = 1; // carry on
     1385          problemStatus_ = -1;
     1386        } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e-2 && (moreSpecialOptions_ & (256 | 8192)) == 0) {
     1387          goToDual = true;
     1388          factorization_->pivotTolerance(CoinMax(0.9, factorization_->pivotTolerance()));
     1389        }
     1390        if (!goToDual) {
     1391          if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST || numberDualInfeasibilities_ == 0) {
     1392            // we are infeasible - use as ray
     1393            delete[] ray_;
     1394            //if ((specialOptions_&(32|0x01000000))!=0x01000000) {
     1395            if ((specialOptions_ & 0x03000000) == 0) {
     1396              ray_ = new double[numberRows_];
     1397              double *saveObjective = CoinCopyOfArray(objective(),
     1398                numberColumns_);
     1399              memset(objective(), 0, numberColumns_ * sizeof(double));
     1400              infeasibilityCost_ = 1.0;
     1401              createRim(4);
     1402              nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1403              gutsOfSolution(NULL, NULL, false);
     1404              memcpy(objective(), saveObjective, numberColumns_ * sizeof(double));
     1405              delete[] saveObjective;
     1406              // swap sign
     1407              for (int i = 0; i < numberRows_; i++)
     1408                ray_[i] = -dual_[i];
     1409            } else {
     1410              ray_ = NULL;
     1411            }
    14351412#ifdef PRINT_RAY_METHOD
    1436                               printf("Primal creating infeasibility ray\n");
    1437 #endif
    1438                               // and get feasible duals
    1439                               infeasibilityCost_ = 0.0;
    1440                               createRim(4);
    1441                               nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1442                               gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1443                               // so will exit
    1444                               infeasibilityCost_ = 1.0e30;
    1445                               // reset infeasibilities
    1446                               sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();;
    1447                               numberPrimalInfeasibilities_ =
    1448                                    nonLinearCost_->numberInfeasibilities();
    1449                          }
    1450                          if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
    1451                            double testValue = MAX_INFEASIBILITY_COST;
    1452                            if (testValue == 1.0e18) {
    1453                              // Check it is not just noise
    1454                              const double * obj = objective();
    1455                              double largestCost=0.0;
    1456                              if (columnScale_) {
    1457                                for (int i=0; i<numberColumns_;i++) {
    1458                                  largestCost =
    1459                                    CoinMax(largestCost,fabs(obj[i]*columnScale_[i]));
    1460                                }
    1461                              } else {
    1462                                for (int i=0; i<numberColumns_;i++) {
    1463                                  largestCost =
    1464                                    CoinMax(largestCost,fabs(obj[i]));
    1465                                }
    1466                              }
    1467                              testValue = 1.0e12*(largestCost+1.0e-6);
    1468                              if (numberDualInfeasibilities_) {
    1469                                double average = sumDualInfeasibilities_/numberDualInfeasibilities_;
    1470                                testValue = CoinMax(testValue,average);
    1471                              }
    1472                              testValue =
    1473                                CoinMin(testValue,MAX_INFEASIBILITY_COST);
    1474                            }
    1475                            if (infeasibilityCost_<testValue) {
    1476                               infeasibilityCost_ *= 5.0;
    1477                               // reset looping criterion
    1478                               progress->reset();
    1479                               changeMade_++; // say change made
    1480                               handler_->message(CLP_PRIMAL_WEIGHT, messages_)
    1481                                         << infeasibilityCost_
    1482                                         << CoinMessageEol;
    1483                               // put back original costs and then check
    1484                               createRim(4);
    1485                               nonLinearCost_->checkInfeasibilities(0.0);
    1486                               gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1487                               problemStatus_ = -1; //continue
    1488                               goToDual = false;
    1489                            } else {
    1490                              // say infeasible
    1491                              problemStatus_ = 1;
    1492                            }
    1493                          } else {
    1494                               // say infeasible
    1495                               problemStatus_ = 1;
    1496                          }
    1497                     }
    1498                }
     1413            printf("Primal creating infeasibility ray\n");
     1414#endif
     1415            // and get feasible duals
     1416            infeasibilityCost_ = 0.0;
     1417            createRim(4);
     1418            nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1419            gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1420            // so will exit
     1421            infeasibilityCost_ = 1.0e30;
     1422            // reset infeasibilities
     1423            sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
     1424            ;
     1425            numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
     1426          }
     1427          if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
     1428            double testValue = MAX_INFEASIBILITY_COST;
     1429            if (testValue == 1.0e18) {
     1430              // Check it is not just noise
     1431              const double *obj = objective();
     1432              double largestCost = 0.0;
     1433              if (columnScale_) {
     1434                for (int i = 0; i < numberColumns_; i++) {
     1435                  largestCost = CoinMax(largestCost, fabs(obj[i] * columnScale_[i]));
     1436                }
     1437              } else {
     1438                for (int i = 0; i < numberColumns_; i++) {
     1439                  largestCost = CoinMax(largestCost, fabs(obj[i]));
     1440                }
     1441              }
     1442              testValue = 1.0e12 * (largestCost + 1.0e-6);
     1443              if (numberDualInfeasibilities_) {
     1444                double average = sumDualInfeasibilities_ / numberDualInfeasibilities_;
     1445                testValue = CoinMax(testValue, average);
     1446              }
     1447              testValue = CoinMin(testValue, MAX_INFEASIBILITY_COST);
     1448            }
     1449            if (infeasibilityCost_ < testValue) {
     1450              infeasibilityCost_ *= 5.0;
     1451              // reset looping criterion
     1452              progress->reset();
     1453              changeMade_++; // say change made
     1454              handler_->message(CLP_PRIMAL_WEIGHT, messages_)
     1455                << infeasibilityCost_
     1456                << CoinMessageEol;
     1457              // put back original costs and then check
     1458              createRim(4);
     1459              nonLinearCost_->checkInfeasibilities(0.0);
     1460              gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1461              problemStatus_ = -1; //continue
     1462              goToDual = false;
     1463            } else {
     1464              // say infeasible
     1465              problemStatus_ = 1;
     1466            }
    14991467          } else {
    1500                // may be optimal
    1501                if (perturbation_ == 101) {
    1502                     goToDual = unPerturb(); // stop any further perturbation
    1503                     if ((numberRows_ > 20000 || numberDualInfeasibilities_) && !numberTimesOptimal_)
    1504                          goToDual = false; // Better to carry on a bit longer
    1505                     lastCleaned = -1; // carry on
    1506                }
    1507                bool unflagged = (unflag() != 0);
    1508                if ( lastCleaned != numberIterations_ || unflagged) {
    1509                     handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
    1510                               << primalTolerance_
    1511                               << CoinMessageEol;
    1512                     if (numberTimesOptimal_ < 4) {
    1513                          numberTimesOptimal_++;
    1514                          changeMade_++; // say change made
    1515                          if (numberTimesOptimal_ == 1) {
    1516                               // better to have small tolerance even if slower
    1517                               factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
    1518                          }
    1519                          lastCleaned = numberIterations_;
    1520                          if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
    1521                               handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
    1522                                         << CoinMessageEol;
    1523                          double oldTolerance = primalTolerance_;
    1524                          primalTolerance_ = dblParam_[ClpPrimalTolerance];
     1468            // say infeasible
     1469            problemStatus_ = 1;
     1470          }
     1471        }
     1472      }
     1473    } else {
     1474      // may be optimal
     1475      if (perturbation_ == 101) {
     1476        goToDual = unPerturb(); // stop any further perturbation
     1477        if ((numberRows_ > 20000 || numberDualInfeasibilities_) && !numberTimesOptimal_)
     1478          goToDual = false; // Better to carry on a bit longer
     1479        lastCleaned = -1; // carry on
     1480      }
     1481      bool unflagged = (unflag() != 0);
     1482      if (lastCleaned != numberIterations_ || unflagged) {
     1483        handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
     1484          << primalTolerance_
     1485          << CoinMessageEol;
     1486        if (numberTimesOptimal_ < 4) {
     1487          numberTimesOptimal_++;
     1488          changeMade_++; // say change made
     1489          if (numberTimesOptimal_ == 1) {
     1490            // better to have small tolerance even if slower
     1491            factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
     1492          }
     1493          lastCleaned = numberIterations_;
     1494          if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
     1495            handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
     1496              << CoinMessageEol;
     1497          double oldTolerance = primalTolerance_;
     1498          primalTolerance_ = dblParam_[ClpPrimalTolerance];
    15251499#if 0
    15261500                         double * xcost = new double[numberRows_+numberColumns_];
     
    15351509                         CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
    15361510#endif
    1537                          // put back original costs and then check
    1538                          createRim(4);
    1539                          nonLinearCost_->checkInfeasibilities(oldTolerance);
     1511          // put back original costs and then check
     1512          createRim(4);
     1513          nonLinearCost_->checkInfeasibilities(oldTolerance);
    15401514#if 0
    15411515                         int i;
     
    15631537                         delete [] xsolution;
    15641538#endif
    1565                          gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1566                          if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
    1567                                    sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
    1568                               // say optimal (with these bounds etc)
    1569                               numberDualInfeasibilities_ = 0;
    1570                               sumDualInfeasibilities_ = 0.0;
    1571                               numberPrimalInfeasibilities_ = 0;
    1572                               sumPrimalInfeasibilities_ = 0.0;
    1573                          }
    1574                          if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
    1575                               problemStatus_ = 0;
    1576                          else
    1577                               problemStatus_ = -1;
    1578                     } else {
    1579                          problemStatus_ = 0; // optimal
    1580                          if (lastCleaned < numberIterations_) {
    1581                               handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
    1582                                         << CoinMessageEol;
    1583                          }
    1584                     }
    1585                } else {
    1586                  /* Previous code here mostly works but
     1539          gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1540          if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     1541            // say optimal (with these bounds etc)
     1542            numberDualInfeasibilities_ = 0;
     1543            sumDualInfeasibilities_ = 0.0;
     1544            numberPrimalInfeasibilities_ = 0;
     1545            sumPrimalInfeasibilities_ = 0.0;
     1546          }
     1547          if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
     1548            problemStatus_ = 0;
     1549          else
     1550            problemStatus_ = -1;
     1551        } else {
     1552          problemStatus_ = 0; // optimal
     1553          if (lastCleaned < numberIterations_) {
     1554            handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
     1555              << CoinMessageEol;
     1556          }
     1557        }
     1558      } else {
     1559        /* Previous code here mostly works but
    15871560                    sumOfRelaxed is rubbish in primal
    15881561                 - so give benefit of doubt still */
    1589                  double error = CoinMin(1.0e-4, largestPrimalError_);
    1590                  // allow bigger tolerance than standard
    1591                  double saveTolerance = primalTolerance_;
    1592                  primalTolerance_ = 2.0*primalTolerance_ + error;
    1593                  nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1594                  double relaxedSum = nonLinearCost_->sumInfeasibilities();
    1595                  // back
    1596                  primalTolerance_=saveTolerance;
    1597                  nonLinearCost_->checkInfeasibilities(primalTolerance_);
    1598                  if (alwaysOptimal || !relaxedSum)
    1599                    problemStatus_ = 0; // optimal
    1600                  else
    1601                    problemStatus_ = 1; // infeasible
    1602                }
    1603           }
    1604      } else {
    1605           // see if looks unbounded
    1606           if (problemStatus_ == -5) {
    1607                if (nonLinearCost_->numberInfeasibilities()) {
    1608                     if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST && perturbation_ == 101) {
    1609                          // back off weight
    1610                          infeasibilityCost_ = 1.0e13;
    1611                          // reset looping criterion
    1612                          progress->reset();
    1613                          unPerturb(); // stop any further perturbation
    1614                     }
    1615                     //we need infeasiblity cost changed
    1616                     if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
    1617                         infeasibilityCost_ *= 5.0;
    1618                         // reset looping criterion
    1619                          progress->reset();
    1620                          changeMade_++; // say change made
    1621                          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
    1622                                    << infeasibilityCost_
    1623                                    << CoinMessageEol;
    1624                          // put back original costs and then check
    1625                          createRim(4);
    1626                          gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    1627                          problemStatus_ = -1; //continue
    1628                     } else {
    1629                          // say infeasible
    1630                          problemStatus_ = 1;
    1631                          // we are infeasible - use as ray
    1632                          delete [] ray_;
    1633                          ray_ = new double [numberRows_];
    1634                          CoinMemcpyN(dual_, numberRows_, ray_);
    1635                     }
    1636                } else {
    1637                     // say unbounded
    1638                     problemStatus_ = 2;
    1639                }
    1640           } else {
    1641                // carry on
    1642                problemStatus_ = -1;
    1643                if(type == 3 && !ifValuesPass) {
    1644                     //bool unflagged =
    1645                     unflag();
    1646                     if (sumDualInfeasibilities_ < 1.0e-3 ||
    1647                               (sumDualInfeasibilities_ / static_cast<double> (numberDualInfeasibilities_)) < 1.0e-5 ||
    1648                               progress->lastIterationNumber(0) == numberIterations_) {
    1649                          if (!numberPrimalInfeasibilities_) {
    1650                               if (numberTimesOptimal_ < 4) {
    1651                                    numberTimesOptimal_++;
    1652                                    changeMade_++; // say change made
    1653                               } else {
    1654                                    problemStatus_ = 0;
    1655                                    secondaryStatus_ = 5;
    1656                               }
    1657                          }
    1658                     }
    1659                }
    1660           }
    1661      }
    1662      if (problemStatus_ == 0) {
    1663           double objVal = (nonLinearCost_->feasibleCost()
    1664                         + objective_->nonlinearOffset());
    1665           objVal /= (objectiveScale_ * rhsScale_);
    1666           double tol = 1.0e-10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e-8;
    1667           if (fabs(objVal - objectiveValue_) > tol) {
     1562        double error = CoinMin(1.0e-4, largestPrimalError_);
     1563        // allow bigger tolerance than standard
     1564        double saveTolerance = primalTolerance_;
     1565        primalTolerance_ = 2.0 * primalTolerance_ + error;
     1566        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1567        double relaxedSum = nonLinearCost_->sumInfeasibilities();
     1568        // back
     1569        primalTolerance_ = saveTolerance;
     1570        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1571        if (alwaysOptimal || !relaxedSum)
     1572          problemStatus_ = 0; // optimal
     1573        else
     1574          problemStatus_ = 1; // infeasible
     1575      }
     1576    }
     1577  } else {
     1578    // see if looks unbounded
     1579    if (problemStatus_ == -5) {
     1580      if (nonLinearCost_->numberInfeasibilities()) {
     1581        if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST && perturbation_ == 101) {
     1582          // back off weight
     1583          infeasibilityCost_ = 1.0e13;
     1584          // reset looping criterion
     1585          progress->reset();
     1586          unPerturb(); // stop any further perturbation
     1587        }
     1588        //we need infeasiblity cost changed
     1589        if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
     1590          infeasibilityCost_ *= 5.0;
     1591          // reset looping criterion
     1592          progress->reset();
     1593          changeMade_++; // say change made
     1594          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
     1595            << infeasibilityCost_
     1596            << CoinMessageEol;
     1597          // put back original costs and then check
     1598          createRim(4);
     1599          gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     1600          problemStatus_ = -1; //continue
     1601        } else {
     1602          // say infeasible
     1603          problemStatus_ = 1;
     1604          // we are infeasible - use as ray
     1605          delete[] ray_;
     1606          ray_ = new double[numberRows_];
     1607          CoinMemcpyN(dual_, numberRows_, ray_);
     1608        }
     1609      } else {
     1610        // say unbounded
     1611        problemStatus_ = 2;
     1612      }
     1613    } else {
     1614      // carry on
     1615      problemStatus_ = -1;
     1616      if (type == 3 && !ifValuesPass) {
     1617        //bool unflagged =
     1618        unflag();
     1619        if (sumDualInfeasibilities_ < 1.0e-3 || (sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_)) < 1.0e-5 || progress->lastIterationNumber(0) == numberIterations_) {
     1620          if (!numberPrimalInfeasibilities_) {
     1621            if (numberTimesOptimal_ < 4) {
     1622              numberTimesOptimal_++;
     1623              changeMade_++; // say change made
     1624            } else {
     1625              problemStatus_ = 0;
     1626              secondaryStatus_ = 5;
     1627            }
     1628          }
     1629        }
     1630      }
     1631    }
     1632  }
     1633  if (problemStatus_ == 0) {
     1634    double objVal = (nonLinearCost_->feasibleCost()
     1635      + objective_->nonlinearOffset());
     1636    objVal /= (objectiveScale_ * rhsScale_);
     1637    double tol = 1.0e-10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e-8;
     1638    if (fabs(objVal - objectiveValue_) > tol) {
    16681639#ifdef COIN_DEVELOP
    1669                if (handler_->logLevel() > 0)
    1670                     printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
    1671                            objVal, objectiveValue_);
    1672 #endif
    1673                objectiveValue_ = objVal;
    1674           }
    1675      }
    1676      // save extra stuff
    1677      matrix_->generalExpanded(this, 5, dummy);
    1678      if (type == 0 || type == 1) {
    1679           if (type != 1 || !saveStatus_) {
    1680                // create save arrays
    1681                delete [] saveStatus_;
    1682                delete [] savedSolution_;
    1683                saveStatus_ = new unsigned char [numberRows_+numberColumns_];
    1684                savedSolution_ = new double [numberRows_+numberColumns_];
    1685           }
    1686           // save arrays
    1687           CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
    1688           CoinMemcpyN(rowActivityWork_,
    1689                       numberRows_, savedSolution_ + numberColumns_);
    1690           CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
    1691      }
    1692      // see if in Cbc etc
    1693      bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
    1694      bool disaster = false;
    1695      if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
    1696           disasterArea_->saveInfo();
    1697           disaster = true;
    1698      }
    1699      if (disaster)
    1700           problemStatus_ = 3;
    1701      if (problemStatus_ < 0 && !changeMade_) {
    1702           problemStatus_ = 4; // unknown
    1703      }
    1704      lastGoodIteration_ = numberIterations_;
    1705      if (numberIterations_ > lastBadIteration_ + 100)
    1706           moreSpecialOptions_ &= ~16; // clear check accuracy flag
    1707      if ((moreSpecialOptions_ & 256) != 0)
    1708        goToDual=false;
    1709      if (goToDual || (numberIterations_ > 1000 && largestPrimalError_ > 1.0e6
    1710                       && largestDualError_ > 1.0e6)) {
    1711           problemStatus_ = 10; // try dual
    1712           // See if second call
    1713           if ((moreSpecialOptions_ & 256) != 0&&nonLinearCost_->sumInfeasibilities()>1.0e2) {
    1714                numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
    1715                sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
    1716                // say infeasible
    1717                if (numberPrimalInfeasibilities_ && largestPrimalError_<1.0e-1)
    1718                     problemStatus_ = 1;
    1719           }
    1720      }
    1721      // make sure first free monotonic
    1722      if (firstFree_ >= 0 && saveFirstFree >= 0) {
    1723        if (numberIterations_) {
    1724          firstFree_=saveFirstFree;
    1725        } else {
    1726          firstFree_ = -1;
    1727          nextSuperBasic(1, NULL);
    1728        }
    1729      }
    1730      if (doFactorization) {
    1731           // restore weights (if saved) - also recompute infeasibility list
    1732           if (tentativeStatus > -3)
    1733                primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
    1734           else
    1735                primalColumnPivot_->saveWeights(this, 3);
    1736           if (saveThreshold) {
    1737                // use default at present
    1738                factorization_->sparseThreshold(0);
    1739                factorization_->goSparse();
    1740           }
    1741      }
    1742      if (problemStatus_==1) {
    1743        // compute true objective value
    1744        computeObjectiveValue(true);
    1745      }
    1746      // Allow matrices to be sorted etc
    1747      int fake = -999; // signal sort
    1748      matrix_->correctSequence(this, fake, fake);
     1640      if (handler_->logLevel() > 0)
     1641        printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
     1642          objVal, objectiveValue_);
     1643#endif
     1644      objectiveValue_ = objVal;
     1645    }
     1646  }
     1647  // save extra stuff
     1648  matrix_->generalExpanded(this, 5, dummy);
     1649  if (type == 0 || type == 1) {
     1650    if (type != 1 || !saveStatus_) {
     1651      // create save arrays
     1652      delete[] saveStatus_;
     1653      delete[] savedSolution_;
     1654      saveStatus_ = new unsigned char[numberRows_ + numberColumns_];
     1655      savedSolution_ = new double[numberRows_ + numberColumns_];
     1656    }
     1657    // save arrays
     1658    CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
     1659    CoinMemcpyN(rowActivityWork_,
     1660      numberRows_, savedSolution_ + numberColumns_);
     1661    CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
     1662  }
     1663  // see if in Cbc etc
     1664  bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
     1665  bool disaster = false;
     1666  if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
     1667    disasterArea_->saveInfo();
     1668    disaster = true;
     1669  }
     1670  if (disaster)
     1671    problemStatus_ = 3;
     1672  if (problemStatus_ < 0 && !changeMade_) {
     1673    problemStatus_ = 4; // unknown
     1674  }
     1675  lastGoodIteration_ = numberIterations_;
     1676  if (numberIterations_ > lastBadIteration_ + 100)
     1677    moreSpecialOptions_ &= ~16; // clear check accuracy flag
     1678  if ((moreSpecialOptions_ & 256) != 0)
     1679    goToDual = false;
     1680  if (goToDual || (numberIterations_ > 1000 && largestPrimalError_ > 1.0e6 && largestDualError_ > 1.0e6)) {
     1681    problemStatus_ = 10; // try dual
     1682    // See if second call
     1683    if ((moreSpecialOptions_ & 256) != 0 && nonLinearCost_->sumInfeasibilities() > 1.0e2) {
     1684      numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
     1685      sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
     1686      // say infeasible
     1687      if (numberPrimalInfeasibilities_ && largestPrimalError_ < 1.0e-1)
     1688        problemStatus_ = 1;
     1689    }
     1690  }
     1691  // make sure first free monotonic
     1692  if (firstFree_ >= 0 && saveFirstFree >= 0) {
     1693    if (numberIterations_) {
     1694      firstFree_ = saveFirstFree;
     1695    } else {
     1696      firstFree_ = -1;
     1697      nextSuperBasic(1, NULL);
     1698    }
     1699  }
     1700  if (doFactorization) {
     1701    // restore weights (if saved) - also recompute infeasibility list
     1702    if (tentativeStatus > -3)
     1703      primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
     1704    else
     1705      primalColumnPivot_->saveWeights(this, 3);
     1706    if (saveThreshold) {
     1707      // use default at present
     1708      factorization_->sparseThreshold(0);
     1709      factorization_->goSparse();
     1710    }
     1711  }
     1712  if (problemStatus_ == 1) {
     1713    // compute true objective value
     1714    computeObjectiveValue(true);
     1715  }
     1716  // Allow matrices to be sorted etc
     1717  int fake = -999; // signal sort
     1718  matrix_->correctSequence(this, fake, fake);
    17491719}
    17501720/*
     
    17551725   On exit rhsArray will have changes in costs of basic variables
    17561726*/
    1757 void
    1758 ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
    1759                             CoinIndexedVector * rhsArray,
    1760                             CoinIndexedVector * spareArray,
    1761                             int valuesPass)
     1727void ClpSimplexPrimal::primalRow(CoinIndexedVector *rowArray,
     1728  CoinIndexedVector *rhsArray,
     1729  CoinIndexedVector *spareArray,
     1730  int valuesPass)
    17621731{
    1763      double saveDj = dualIn_;
    1764      if (valuesPass && objective_->type() < 2) {
    1765           dualIn_ = cost_[sequenceIn_];
    1766 
    1767           double * work = rowArray->denseVector();
    1768           int number = rowArray->getNumElements();
    1769           int * which = rowArray->getIndices();
    1770 
    1771           int iIndex;
    1772           for (iIndex = 0; iIndex < number; iIndex++) {
    1773 
    1774                int iRow = which[iIndex];
    1775                double alpha = work[iIndex];
    1776                int iPivot = pivotVariable_[iRow];
    1777                dualIn_ -= alpha * cost_[iPivot];
    1778           }
    1779           // determine direction here
    1780           if (dualIn_ < -dualTolerance_) {
    1781                directionIn_ = 1;
    1782           } else if (dualIn_ > dualTolerance_) {
    1783                directionIn_ = -1;
    1784           } else {
    1785                // towards nearest bound
    1786                if (valueIn_ - lowerIn_ < upperIn_ - valueIn_) {
    1787                     directionIn_ = -1;
    1788                     dualIn_ = dualTolerance_;
    1789                } else {
    1790                     directionIn_ = 1;
    1791                     dualIn_ = -dualTolerance_;
    1792                }
    1793           }
    1794      }
    1795 
    1796      // sequence stays as row number until end
    1797      pivotRow_ = -1;
    1798      int numberRemaining = 0;
    1799 
    1800      double totalThru = 0.0; // for when variables flip
    1801      // Allow first few iterations to take tiny
    1802      double acceptablePivot = 1.0e-1 * acceptablePivot_;
    1803      if (numberIterations_ > 100)
    1804           acceptablePivot = acceptablePivot_;
    1805      if (factorization_->pivots() > 10)
    1806           acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
    1807      else if (factorization_->pivots() > 5)
    1808           acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
    1809      else if (factorization_->pivots())
    1810           acceptablePivot = acceptablePivot_; // relax
    1811      double bestEverPivot = acceptablePivot;
    1812      int lastPivotRow = -1;
    1813      double lastPivot = 0.0;
    1814      double lastTheta = 1.0e50;
    1815 
    1816      // use spareArrays to put ones looked at in
    1817      // First one is list of candidates
    1818      // We could compress if we really know we won't need any more
    1819      // Second array has current set of pivot candidates
    1820      // with a backup list saved in double * part of indexed vector
    1821 
    1822      // pivot elements
    1823      double * spare;
    1824      // indices
    1825      int * index;
    1826      spareArray->clear();
    1827      spare = spareArray->denseVector();
    1828      index = spareArray->getIndices();
    1829 
    1830      // we also need somewhere for effective rhs
    1831      double * rhs = rhsArray->denseVector();
    1832      // and we can use indices to point to alpha
    1833      // that way we can store fabs(alpha)
    1834      int * indexPoint = rhsArray->getIndices();
    1835      //int numberFlip=0; // Those which may change if flips
    1836 
    1837      /*
     1732  double saveDj = dualIn_;
     1733  if (valuesPass && objective_->type() < 2) {
     1734    dualIn_ = cost_[sequenceIn_];
     1735
     1736    double *work = rowArray->denseVector();
     1737    int number = rowArray->getNumElements();
     1738    int *which = rowArray->getIndices();
     1739
     1740    int iIndex;
     1741    for (iIndex = 0; iIndex < number; iIndex++) {
     1742
     1743      int iRow = which[iIndex];
     1744      double alpha = work[iIndex];
     1745      int iPivot = pivotVariable_[iRow];
     1746      dualIn_ -= alpha * cost_[iPivot];
     1747    }
     1748    // determine direction here
     1749    if (dualIn_ < -dualTolerance_) {
     1750      directionIn_ = 1;
     1751    } else if (dualIn_ > dualTolerance_) {
     1752      directionIn_ = -1;
     1753    } else {
     1754      // towards nearest bound
     1755      if (valueIn_ - lowerIn_ < upperIn_ - valueIn_) {
     1756        directionIn_ = -1;
     1757        dualIn_ = dualTolerance_;
     1758      } else {
     1759        directionIn_ = 1;
     1760        dualIn_ = -dualTolerance_;
     1761      }
     1762    }
     1763  }
     1764
     1765  // sequence stays as row number until end
     1766  pivotRow_ = -1;
     1767  int numberRemaining = 0;
     1768
     1769  double totalThru = 0.0; // for when variables flip
     1770  // Allow first few iterations to take tiny
     1771  double acceptablePivot = 1.0e-1 * acceptablePivot_;
     1772  if (numberIterations_ > 100)
     1773    acceptablePivot = acceptablePivot_;
     1774  if (factorization_->pivots() > 10)
     1775    acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
     1776  else if (factorization_->pivots() > 5)
     1777    acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
     1778  else if (factorization_->pivots())
     1779    acceptablePivot = acceptablePivot_; // relax
     1780  double bestEverPivot = acceptablePivot;
     1781  int lastPivotRow = -1;
     1782  double lastPivot = 0.0;
     1783  double lastTheta = 1.0e50;
     1784
     1785  // use spareArrays to put ones looked at in
     1786  // First one is list of candidates
     1787  // We could compress if we really know we won't need any more
     1788  // Second array has current set of pivot candidates
     1789  // with a backup list saved in double * part of indexed vector
     1790
     1791  // pivot elements
     1792  double *spare;
     1793  // indices
     1794  int *index;
     1795  spareArray->clear();
     1796  spare = spareArray->denseVector();
     1797  index = spareArray->getIndices();
     1798
     1799  // we also need somewhere for effective rhs
     1800  double *rhs = rhsArray->denseVector();
     1801  // and we can use indices to point to alpha
     1802  // that way we can store fabs(alpha)
     1803  int *indexPoint = rhsArray->getIndices();
     1804  //int numberFlip=0; // Those which may change if flips
     1805
     1806  /*
    18381807       First we get a list of possible pivots.  We can also see if the
    18391808       problem looks unbounded.
     
    18451814      */
    18461815
    1847      // do first pass to get possibles
    1848      // We can also see if unbounded
    1849 
    1850      double * work = rowArray->denseVector();
    1851      int number = rowArray->getNumElements();
    1852      int * which = rowArray->getIndices();
    1853 
    1854      // we need to swap sign if coming in from ub
    1855      double way = directionIn_;
    1856      double maximumMovement;
    1857      if (way > 0.0)
    1858           maximumMovement = CoinMin(1.0e30, upperIn_ - valueIn_);
    1859      else
    1860           maximumMovement = CoinMin(1.0e30, valueIn_ - lowerIn_);
    1861 
    1862      double averageTheta = nonLinearCost_->averageTheta();
    1863      double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
    1864      double upperTheta = maximumMovement;
    1865      if (tentativeTheta > 0.5 * maximumMovement)
    1866           tentativeTheta = maximumMovement;
    1867      bool thetaAtMaximum = tentativeTheta == maximumMovement;
    1868      // In case tiny bounds increase
    1869      if (maximumMovement < 1.0)
    1870           tentativeTheta *= 1.1;
    1871      double dualCheck = fabs(dualIn_);
    1872      // but make a bit more pessimistic
    1873      dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
    1874 
    1875      int iIndex;
    1876      int pivotOne = -1;
    1877      //#define CLP_DEBUG
     1816  // do first pass to get possibles
     1817  // We can also see if unbounded
     1818
     1819  double *work = rowArray->denseVector();
     1820  int number = rowArray->getNumElements();
     1821  int *which = rowArray->getIndices();
     1822
     1823  // we need to swap sign if coming in from ub
     1824  double way = directionIn_;
     1825  double maximumMovement;
     1826  if (way > 0.0)
     1827    maximumMovement = CoinMin(1.0e30, upperIn_ - valueIn_);
     1828  else
     1829    maximumMovement = CoinMin(1.0e30, valueIn_ - lowerIn_);
     1830
     1831  double averageTheta = nonLinearCost_->averageTheta();
     1832  double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
     1833  double upperTheta = maximumMovement;
     1834  if (tentativeTheta > 0.5 * maximumMovement)
     1835    tentativeTheta = maximumMovement;
     1836  bool thetaAtMaximum = tentativeTheta == maximumMovement;
     1837  // In case tiny bounds increase
     1838  if (maximumMovement < 1.0)
     1839    tentativeTheta *= 1.1;
     1840  double dualCheck = fabs(dualIn_);
     1841  // but make a bit more pessimistic
     1842  dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
     1843
     1844  int iIndex;
     1845  int pivotOne = -1;
     1846  //#define CLP_DEBUG
    18781847#ifdef CLP_DEBUG
    1879      if (numberIterations_ == -3839 || numberIterations_ == -3840) {
    1880           double dj = cost_[sequenceIn_];
    1881           printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_);
    1882           for (iIndex = 0; iIndex < number; iIndex++) {
    1883 
    1884                int iRow = which[iIndex];
    1885                double alpha = work[iIndex];
    1886                int iPivot = pivotVariable_[iRow];
    1887                dj -= alpha * cost_[iPivot];
    1888                printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
    1889                       iRow, iPivot, lower_[iPivot], solution_[iPivot], upper_[iPivot],
    1890                       alpha, solution_[iPivot] - 1.0e9 * alpha, cost_[iPivot], dj);
    1891           }
    1892      }
    1893 #endif
    1894      while (true) {
    1895           pivotOne = -1;
    1896           totalThru = 0.0;
    1897           //double totalThruFake=0.0;
    1898           //int nThru=0;
    1899           // We also re-compute reduced cost
    1900           numberRemaining = 0;
    1901           dualIn_ = cost_[sequenceIn_];
     1848  if (numberIterations_ == -3839 || numberIterations_ == -3840) {
     1849    double dj = cost_[sequenceIn_];
     1850    printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_);
     1851    for (iIndex = 0; iIndex < number; iIndex++) {
     1852
     1853      int iRow = which[iIndex];
     1854      double alpha = work[iIndex];
     1855      int iPivot = pivotVariable_[iRow];
     1856      dj -= alpha * cost_[iPivot];
     1857      printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
     1858        iRow, iPivot, lower_[iPivot], solution_[iPivot], upper_[iPivot],
     1859        alpha, solution_[iPivot] - 1.0e9 * alpha, cost_[iPivot], dj);
     1860    }
     1861  }
     1862#endif
     1863  while (true) {
     1864    pivotOne = -1;
     1865    totalThru = 0.0;
     1866    //double totalThruFake=0.0;
     1867    //int nThru=0;
     1868    // We also re-compute reduced cost
     1869    numberRemaining = 0;
     1870    dualIn_ = cost_[sequenceIn_];
    19021871#ifndef NDEBUG
    1903           //double tolerance = primalTolerance_ * 1.002;
    1904 #endif
    1905           for (iIndex = 0; iIndex < number; iIndex++) {
    1906 
    1907                int iRow = which[iIndex];
    1908                double alpha = work[iIndex];
    1909                int iPivot = pivotVariable_[iRow];
    1910                // faster without if (cost_[iPivot])
    1911                dualIn_ -= alpha * cost_[iPivot];
    1912                alpha *= way;
    1913                double oldValue = solution_[iPivot];
    1914                // get where in bound sequence
    1915                // note that after this alpha is actually fabs(alpha)
    1916                bool possible;
    1917                // do computation same way as later on in primal
    1918                if (alpha > 0.0) {
    1919                     // basic variable going towards lower bound
    1920                     double bound = lower_[iPivot];
    1921                     // must be exactly same as when used
    1922                     double change = tentativeTheta * alpha;
    1923                     possible = (oldValue - change) <= bound + primalTolerance_;
    1924                     oldValue -= bound;
    1925                } else {
    1926                     // basic variable going towards upper bound
    1927                     double bound = upper_[iPivot];
    1928                     // must be exactly same as when used
    1929                     double change = tentativeTheta * alpha;
    1930                     possible = (oldValue - change) >= bound - primalTolerance_;
    1931                     oldValue = bound - oldValue;
    1932                     alpha = - alpha;
    1933                }
    1934                double value;
    1935                //assert (oldValue >= -10.0*tolerance);
    1936                if (possible) {
    1937                     value = oldValue - upperTheta * alpha;
     1872    //double tolerance = primalTolerance_ * 1.002;
     1873#endif
     1874    for (iIndex = 0; iIndex < number; iIndex++) {
     1875
     1876      int iRow = which[iIndex];
     1877      double alpha = work[iIndex];
     1878      int iPivot = pivotVariable_[iRow];
     1879      // faster without if (cost_[iPivot])
     1880      dualIn_ -= alpha * cost_[iPivot];
     1881      alpha *= way;
     1882      double oldValue = solution_[iPivot];
     1883      // get where in bound sequence
     1884      // note that after this alpha is actually fabs(alpha)
     1885      bool possible;
     1886      // do computation same way as later on in primal
     1887      if (alpha > 0.0) {
     1888        // basic variable going towards lower bound
     1889        double bound = lower_[iPivot];
     1890        // must be exactly same as when used
     1891        double change = tentativeTheta * alpha;
     1892        possible = (oldValue - change) <= bound + primalTolerance_;
     1893        oldValue -= bound;
     1894      } else {
     1895        // basic variable going towards upper bound
     1896        double bound = upper_[iPivot];
     1897        // must be exactly same as when used
     1898        double change = tentativeTheta * alpha;
     1899        possible = (oldValue - change) >= bound - primalTolerance_;
     1900        oldValue = bound - oldValue;
     1901        alpha = -alpha;
     1902      }
     1903      double value;
     1904      //assert (oldValue >= -10.0*tolerance);
     1905      if (possible) {
     1906        value = oldValue - upperTheta * alpha;
    19381907#ifdef CLP_USER_DRIVEN1
    1939                     if(!userChoiceValid1(this,iPivot,oldValue,
    1940                                          upperTheta,alpha,work[iIndex]*way))
    1941                       value =0.0; // say can't use
    1942 #endif
    1943                     if (value < -primalTolerance_ && alpha >= acceptablePivot) {
    1944                          upperTheta = (oldValue + primalTolerance_) / alpha;
    1945                          pivotOne = numberRemaining;
    1946                     }
    1947                     // add to list
    1948                     spare[numberRemaining] = alpha;
    1949                     rhs[numberRemaining] = oldValue;
    1950                     indexPoint[numberRemaining] = iIndex;
    1951                     index[numberRemaining++] = iRow;
    1952                     totalThru += alpha; // need more if dubious pricing
     1908        if (!userChoiceValid1(this, iPivot, oldValue,
     1909              upperTheta, alpha, work[iIndex] * way))
     1910          value = 0.0; // say can't use
     1911#endif
     1912        if (value < -primalTolerance_ && alpha >= acceptablePivot) {
     1913          upperTheta = (oldValue + primalTolerance_) / alpha;
     1914          pivotOne = numberRemaining;
     1915        }
     1916        // add to list
     1917        spare[numberRemaining] = alpha;
     1918        rhs[numberRemaining] = oldValue;
     1919        indexPoint[numberRemaining] = iIndex;
     1920        index[numberRemaining++] = iRow;
     1921        totalThru += alpha; // need more if dubious pricing
    19531922#if 0 // wasdef CLP_USER_DRIVEN
    19541923                    clpUserStruct info;
     
    19671936                      printf("%d thru - fake %g real %g\n",nThru,totalThruFake,totalThru);
    19681937                    nThru++;
    1969 #endif             
    1970                     setActive(iRow);
    1971                     //} else if (value<primalTolerance_*1.002) {
    1972                     // May change if is a flip
    1973                     //indexRhs[numberFlip++]=iRow;
    1974                }
    1975           }
     1938#endif
     1939        setActive(iRow);
     1940        //} else if (value<primalTolerance_*1.002) {
     1941        // May change if is a flip
     1942        //indexRhs[numberFlip++]=iRow;
     1943      }
     1944    }
    19761945#if 0 // was def CLP_USER_DRIVEN
    19771946          clpUserStruct info;
     
    19901959            }
    19911960          }
    1992 #endif             
    1993           if (upperTheta < maximumMovement && totalThru*infeasibilityCost_ >= 1.0001 * dualCheck) {
    1994                // Can pivot here
    1995                break;
    1996           } else if (!thetaAtMaximum) {
    1997                //printf("Going round with average theta of %g\n",averageTheta);
    1998                tentativeTheta = maximumMovement;
    1999                thetaAtMaximum = true; // seems to be odd compiler error
    2000           } else {
    2001                break;
    2002           }
    2003      }
    2004      totalThru = 0.0;
    2005 
    2006      theta_ = maximumMovement;
    2007 
    2008      bool goBackOne = false;
    2009      if (objective_->type() > 1)
    2010           dualIn_ = saveDj;
    2011 
    2012      //printf("%d remain out of %d\n",numberRemaining,number);
    2013      int iTry = 0;
     1961#endif
     1962    if (upperTheta < maximumMovement && totalThru * infeasibilityCost_ >= 1.0001 * dualCheck) {
     1963      // Can pivot here
     1964      break;
     1965    } else if (!thetaAtMaximum) {
     1966      //printf("Going round with average theta of %g\n",averageTheta);
     1967      tentativeTheta = maximumMovement;
     1968      thetaAtMaximum = true; // seems to be odd compiler error
     1969    } else {
     1970      break;
     1971    }
     1972  }
     1973  totalThru = 0.0;
     1974
     1975  theta_ = maximumMovement;
     1976
     1977  bool goBackOne = false;
     1978  if (objective_->type() > 1)
     1979    dualIn_ = saveDj;
     1980
     1981  //printf("%d remain out of %d\n",numberRemaining,number);
     1982  int iTry = 0;
    20141983#define MAXTRY 1000
    2015      if (numberRemaining && upperTheta < maximumMovement) {
    2016           // First check if previously chosen one will work
    2017           if (pivotOne >= 0 && 0) {
    2018                double thruCost = infeasibilityCost_ * spare[pivotOne];
    2019                if (thruCost >= 0.99 * fabs(dualIn_))
    2020                     COIN_DETAIL_PRINT(printf("Could pivot on %d as change %g dj %g\n",
    2021                                              index[pivotOne], thruCost, dualIn_));
    2022                double alpha = spare[pivotOne];
    2023                double oldValue = rhs[pivotOne];
    2024                theta_ = oldValue / alpha;
    2025                pivotRow_ = pivotOne;
    2026                // Stop loop
    2027                iTry = MAXTRY;
    2028           }
    2029 
    2030           // first get ratio with tolerance
    2031           for ( ; iTry < MAXTRY; iTry++) {
    2032 
    2033                upperTheta = maximumMovement;
    2034                int iBest = -1;
    2035                for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
    2036 
    2037                     double alpha = spare[iIndex];
    2038                     double oldValue = rhs[iIndex];
    2039                     double value = oldValue - upperTheta * alpha;
     1984  if (numberRemaining && upperTheta < maximumMovement) {
     1985    // First check if previously chosen one will work
     1986    if (pivotOne >= 0 && 0) {
     1987      double thruCost = infeasibilityCost_ * spare[pivotOne];
     1988      if (thruCost >= 0.99 * fabs(dualIn_))
     1989        COIN_DETAIL_PRINT(printf("Could pivot on %d as change %g dj %g\n",
     1990          index[pivotOne], thruCost, dualIn_));
     1991      double alpha = spare[pivotOne];
     1992      double oldValue = rhs[pivotOne];
     1993      theta_ = oldValue / alpha;
     1994      pivotRow_ = pivotOne;
     1995      // Stop loop
     1996      iTry = MAXTRY;
     1997    }
     1998
     1999    // first get ratio with tolerance
     2000    for (; iTry < MAXTRY; iTry++) {
     2001
     2002      upperTheta = maximumMovement;
     2003      int iBest = -1;
     2004      for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
     2005
     2006        double alpha = spare[iIndex];
     2007        double oldValue = rhs[iIndex];
     2008        double value = oldValue - upperTheta * alpha;
    20402009
    20412010#ifdef CLP_USER_DRIVEN1
    2042                     int sequenceOut=pivotVariable_[index[iIndex]];
    2043                     if(!userChoiceValid1(this,sequenceOut,oldValue,
    2044                                          upperTheta,alpha, 0.0))
    2045                       value =0.0; // say can't use
    2046 #endif
    2047                     if (value < -primalTolerance_ && alpha >= acceptablePivot) {
    2048                          upperTheta = (oldValue + primalTolerance_) / alpha;
    2049                          iBest = iIndex; // just in case weird numbers
    2050                     }
    2051                }
    2052                // now look at best in this lot
    2053                // But also see how infeasible small pivots will make
    2054                double sumInfeasibilities = 0.0;
    2055                double bestPivot = acceptablePivot;
    2056                pivotRow_ = -1;
    2057                for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
    2058 
    2059                     int iRow = index[iIndex];
    2060                     double alpha = spare[iIndex];
    2061                     double oldValue = rhs[iIndex];
    2062                     double value = oldValue - upperTheta * alpha;
    2063 
    2064                     if (value <= 0 || iBest == iIndex) {
    2065                          // how much would it cost to go thru and modify bound
    2066                          double trueAlpha = way * work[indexPoint[iIndex]];
    2067                          int jSequence=pivotVariable_[iRow];
    2068                         //if (printing)
    2069                         //jSequence += 1000000;
    2070                          totalThru += nonLinearCost_->changeInCost(jSequence, trueAlpha, rhs[iIndex]);
     2011        int sequenceOut = pivotVariable_[index[iIndex]];
     2012        if (!userChoiceValid1(this, sequenceOut, oldValue,
     2013              upperTheta, alpha, 0.0))
     2014          value = 0.0; // say can't use
     2015#endif
     2016        if (value < -primalTolerance_ && alpha >= acceptablePivot) {
     2017          upperTheta = (oldValue + primalTolerance_) / alpha;
     2018          iBest = iIndex; // just in case weird numbers
     2019        }
     2020      }
     2021      // now look at best in this lot
     2022      // But also see how infeasible small pivots will make
     2023      double sumInfeasibilities = 0.0;
     2024      double bestPivot = acceptablePivot;
     2025      pivotRow_ = -1;
     2026      for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
     2027
     2028        int iRow = index[iIndex];
     2029        double alpha = spare[iIndex];
     2030        double oldValue = rhs[iIndex];
     2031        double value = oldValue - upperTheta * alpha;
     2032
     2033        if (value <= 0 || iBest == iIndex) {
     2034          // how much would it cost to go thru and modify bound
     2035          double trueAlpha = way * work[indexPoint[iIndex]];
     2036          int jSequence = pivotVariable_[iRow];
     2037          //if (printing)
     2038          //jSequence += 1000000;
     2039          totalThru += nonLinearCost_->changeInCost(jSequence, trueAlpha, rhs[iIndex]);
    20712040#if 0 //was def CLP_USER_DRIVEN
    20722041                         clpUserStruct info;
     
    20822051                         totalThru=info.totalThru;
    20832052                         rhs[iIndex]=info.rhs;
    2084 #endif             
    2085                          setActive(iRow);
    2086                          if (alpha > bestPivot) {
    2087                               bestPivot = alpha;
    2088                               theta_ = oldValue / bestPivot;
    2089                               pivotRow_ = iIndex;
    2090                          } else if (alpha < acceptablePivot
     2053#endif
     2054          setActive(iRow);
     2055          if (alpha > bestPivot) {
     2056            bestPivot = alpha;
     2057            theta_ = oldValue / bestPivot;
     2058            pivotRow_ = iIndex;
     2059          } else if (alpha < acceptablePivot
    20912060#ifdef CLP_USER_DRIVEN1
    2092                       ||!userChoiceValid1(this,pivotVariable_[index[iIndex]],
    2093                                           oldValue,upperTheta,alpha,0.0)
    2094 #endif
    2095                                     ) {
    2096                               if (value < -primalTolerance_)
    2097                                    sumInfeasibilities += -value - primalTolerance_;
    2098                          }
    2099                     }
    2100                }
     2061            || !userChoiceValid1(this, pivotVariable_[index[iIndex]],
     2062                 oldValue, upperTheta, alpha, 0.0)
     2063#endif
     2064          ) {
     2065            if (value < -primalTolerance_)
     2066              sumInfeasibilities += -value - primalTolerance_;
     2067          }
     2068        }
     2069      }
    21012070#if 0 // was def CLP_USER_DRIVEN
    21022071               clpUserStruct info;
     
    21052074               eventHandler_->eventWithInfo(ClpEventHandler::pivotRow,
    21062075                                       &info);
    2107 #endif             
    2108                if (bestPivot < 0.1 * bestEverPivot &&
    2109                          bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
    2110                     // back to previous one
    2111                     goBackOne = true;
    2112                     break;
    2113                } else if (pivotRow_ == -1 && upperTheta > largeValue_) {
    2114                     if (lastPivot > acceptablePivot) {
    2115                          // back to previous one
    2116                          goBackOne = true;
    2117                     } else {
    2118                          // can only get here if all pivots so far too small
    2119                     }
    2120                     break;
    2121                } else if (totalThru >= dualCheck) {
    2122                     if (sumInfeasibilities > primalTolerance_ && !nonLinearCost_->numberInfeasibilities()) {
    2123                          // Looks a bad choice
    2124                          if (lastPivot > acceptablePivot) {
    2125                               goBackOne = true;
    2126                          } else {
    2127                               // say no good
    2128                               dualIn_ = 0.0;
    2129                          }
    2130                     }
    2131                     break; // no point trying another loop
    2132                } else {
    2133                     lastPivotRow = pivotRow_;
    2134                     lastTheta = theta_;
    2135                     if (bestPivot > bestEverPivot)
    2136                          bestEverPivot = bestPivot;
    2137                }
    2138           }
    2139           // can get here without pivotRow_ set but with lastPivotRow
    2140           if (goBackOne || (pivotRow_ < 0 && lastPivotRow >= 0)) {
    2141                // back to previous one
    2142                pivotRow_ = lastPivotRow;
    2143                theta_ = lastTheta;
    2144           }
    2145      } else if (pivotRow_ < 0 && maximumMovement > 1.0e20) {
    2146           // looks unbounded
    2147           valueOut_ = COIN_DBL_MAX; // say odd
    2148           if (nonLinearCost_->numberInfeasibilities()) {
    2149                // but infeasible??
    2150                // move variable but don't pivot
    2151                tentativeTheta = 1.0e50;
    2152                for (iIndex = 0; iIndex < number; iIndex++) {
    2153                     int iRow = which[iIndex];
    2154                     double alpha = work[iIndex];
    2155                     int iPivot = pivotVariable_[iRow];
    2156                     alpha *= way;
    2157                     double oldValue = solution_[iPivot];
    2158                     // get where in bound sequence
    2159                     // note that after this alpha is actually fabs(alpha)
    2160                     if (alpha > 0.0) {
    2161                          // basic variable going towards lower bound
    2162                          double bound = lower_[iPivot];
    2163                          oldValue -= bound;
    2164                     } else {
    2165                          // basic variable going towards upper bound
    2166                          double bound = upper_[iPivot];
    2167                          oldValue = bound - oldValue;
    2168                          alpha = - alpha;
    2169                     }
    2170                     if (oldValue - tentativeTheta * alpha < 0.0) {
    2171                          tentativeTheta = oldValue / alpha;
    2172                     }
    2173                }
    2174                // If free in then see if we can get to 0.0
    2175                if (lowerIn_ < -1.0e20 && upperIn_ > 1.0e20) {
    2176                     if (dualIn_ * valueIn_ > 0.0) {
    2177                          if (fabs(valueIn_) < 1.0e-2 && (tentativeTheta < fabs(valueIn_) || tentativeTheta > 1.0e20)) {
    2178                               tentativeTheta = fabs(valueIn_);
    2179                          }
    2180                     }
    2181                }
    2182                if (tentativeTheta < 1.0e10)
    2183                     valueOut_ = valueIn_ + way * tentativeTheta;
    2184           }
    2185      }
    2186      //if (iTry>50)
    2187      //printf("** %d tries\n",iTry);
    2188      if (pivotRow_ >= 0) {
    2189           int position = pivotRow_; // position in list
    2190           pivotRow_ = index[position];
    2191           alpha_ = work[indexPoint[position]];
    2192           // translate to sequence
    2193           sequenceOut_ = pivotVariable_[pivotRow_];
    2194           valueOut_ = solution(sequenceOut_);
    2195           lowerOut_ = lower_[sequenceOut_];
    2196           upperOut_ = upper_[sequenceOut_];
     2076#endif
     2077      if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
     2078        // back to previous one
     2079        goBackOne = true;
     2080        break;
     2081      } else if (pivotRow_ == -1 && upperTheta > largeValue_) {
     2082        if (lastPivot > acceptablePivot) {
     2083          // back to previous one
     2084          goBackOne = true;
     2085        } else {
     2086          // can only get here if all pivots so far too small
     2087        }
     2088        break;
     2089      } else if (totalThru >= dualCheck) {
     2090        if (sumInfeasibilities > primalTolerance_ && !nonLinearCost_->numberInfeasibilities()) {
     2091          // Looks a bad choice
     2092          if (lastPivot > acceptablePivot) {
     2093            goBackOne = true;
     2094          } else {
     2095            // say no good
     2096            dualIn_ = 0.0;
     2097          }
     2098        }
     2099        break; // no point trying another loop
     2100      } else {
     2101        lastPivotRow = pivotRow_;
     2102        lastTheta = theta_;
     2103        if (bestPivot > bestEverPivot)
     2104          bestEverPivot = bestPivot;
     2105      }
     2106    }
     2107    // can get here without pivotRow_ set but with lastPivotRow
     2108    if (goBackOne || (pivotRow_ < 0 && lastPivotRow >= 0)) {
     2109      // back to previous one
     2110      pivotRow_ = lastPivotRow;
     2111      theta_ = lastTheta;
     2112    }
     2113  } else if (pivotRow_ < 0 && maximumMovement > 1.0e20) {
     2114    // looks unbounded
     2115    valueOut_ = COIN_DBL_MAX; // say odd
     2116    if (nonLinearCost_->numberInfeasibilities()) {
     2117      // but infeasible??
     2118      // move variable but don't pivot
     2119      tentativeTheta = 1.0e50;
     2120      for (iIndex = 0; iIndex < number; iIndex++) {
     2121        int iRow = which[iIndex];
     2122        double alpha = work[iIndex];
     2123        int iPivot = pivotVariable_[iRow];
     2124        alpha *= way;
     2125        double oldValue = solution_[iPivot];
     2126        // get where in bound sequence
     2127        // note that after this alpha is actually fabs(alpha)
     2128        if (alpha > 0.0) {
     2129          // basic variable going towards lower bound
     2130          double bound = lower_[iPivot];
     2131          oldValue -= bound;
     2132        } else {
     2133          // basic variable going towards upper bound
     2134          double bound = upper_[iPivot];
     2135          oldValue = bound - oldValue;
     2136          alpha = -alpha;
     2137        }
     2138        if (oldValue - tentativeTheta * alpha < 0.0) {
     2139          tentativeTheta = oldValue / alpha;
     2140        }
     2141      }
     2142      // If free in then see if we can get to 0.0
     2143      if (lowerIn_ < -1.0e20 && upperIn_ > 1.0e20) {
     2144        if (dualIn_ * valueIn_ > 0.0) {
     2145          if (fabs(valueIn_) < 1.0e-2 && (tentativeTheta < fabs(valueIn_) || tentativeTheta > 1.0e20)) {
     2146            tentativeTheta = fabs(valueIn_);
     2147          }
     2148        }
     2149      }
     2150      if (tentativeTheta < 1.0e10)
     2151        valueOut_ = valueIn_ + way * tentativeTheta;
     2152    }
     2153  }
     2154  //if (iTry>50)
     2155  //printf("** %d tries\n",iTry);
     2156  if (pivotRow_ >= 0) {
     2157    int position = pivotRow_; // position in list
     2158    pivotRow_ = index[position];
     2159    alpha_ = work[indexPoint[position]];
     2160    // translate to sequence
     2161    sequenceOut_ = pivotVariable_[pivotRow_];
     2162    valueOut_ = solution(sequenceOut_);
     2163    lowerOut_ = lower_[sequenceOut_];
     2164    upperOut_ = upper_[sequenceOut_];
    21972165#define MINIMUMTHETA 1.0e-12
    2198           // Movement should be minimum for anti-degeneracy - unless
    2199           // fixed variable out
    2200           double minimumTheta;
    2201           if (upperOut_ > lowerOut_)
    2202                minimumTheta = MINIMUMTHETA;
    2203           else
    2204                minimumTheta = 0.0;
    2205           // But can't go infeasible
    2206           double distance;
    2207           if (alpha_ * way > 0.0)
    2208                distance = valueOut_ - lowerOut_;
    2209           else
    2210                distance = upperOut_ - valueOut_;
    2211           if (distance - minimumTheta * fabs(alpha_) < -primalTolerance_)
    2212                minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha_));
    2213           // will we need to increase tolerance
    2214           //#define CLP_DEBUG
    2215           double largestInfeasibility = primalTolerance_;
    2216           if (theta_ < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
    2217                theta_ = minimumTheta;
    2218                for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
    2219                     largestInfeasibility = CoinMax(largestInfeasibility,
    2220                                                    -(rhs[iIndex] - spare[iIndex] * theta_));
    2221                }
     2166    // Movement should be minimum for anti-degeneracy - unless
     2167    // fixed variable out
     2168    double minimumTheta;
     2169    if (upperOut_ > lowerOut_)
     2170      minimumTheta = MINIMUMTHETA;
     2171    else
     2172      minimumTheta = 0.0;
     2173    // But can't go infeasible
     2174    double distance;
     2175    if (alpha_ * way > 0.0)
     2176      distance = valueOut_ - lowerOut_;
     2177    else
     2178      distance = upperOut_ - valueOut_;
     2179    if (distance - minimumTheta * fabs(alpha_) < -primalTolerance_)
     2180      minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha_));
     2181    // will we need to increase tolerance
     2182    //#define CLP_DEBUG
     2183    double largestInfeasibility = primalTolerance_;
     2184    if (theta_ < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
     2185      theta_ = minimumTheta;
     2186      for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
     2187        largestInfeasibility = CoinMax(largestInfeasibility,
     2188          -(rhs[iIndex] - spare[iIndex] * theta_));
     2189      }
    22222190//#define CLP_DEBUG
    22232191#ifdef CLP_DEBUG
    2224                if (largestInfeasibility > primalTolerance_ && (handler_->logLevel() & 32) > -1)
    2225                     printf("Primal tolerance increased from %g to %g\n",
    2226                            primalTolerance_, largestInfeasibility);
    2227 #endif
    2228 //#undef CLP_DEBUG
    2229                primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
    2230           }
    2231           // Need to look at all in some cases
    2232           if (theta_ > tentativeTheta) {
    2233                for (iIndex = 0; iIndex < number; iIndex++)
    2234                     setActive(which[iIndex]);
    2235           }
    2236           if (way < 0.0)
    2237                theta_ = - theta_;
    2238           double newValue = valueOut_ - theta_ * alpha_;
    2239           // If  4 bit set - Force outgoing variables to exact bound (primal)
    2240           if (alpha_ * way < 0.0) {
    2241                directionOut_ = -1;    // to upper bound
    2242                if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
    2243                     upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
    2244                } else {
    2245                     upperOut_ = newValue;
    2246                }
    2247           } else {
    2248                directionOut_ = 1;    // to lower bound
    2249                if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
    2250                     lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
    2251                } else {
    2252                     lowerOut_ = newValue;
    2253                }
    2254           }
    2255           dualOut_ = reducedCost(sequenceOut_);
    2256      } else if (maximumMovement < 1.0e20) {
    2257           // flip
    2258           pivotRow_ = -2; // so we can tell its a flip
    2259           sequenceOut_ = sequenceIn_;
    2260           valueOut_ = valueIn_;
    2261           dualOut_ = dualIn_;
    2262           lowerOut_ = lowerIn_;
    2263           upperOut_ = upperIn_;
    2264           alpha_ = 0.0;
    2265           if (way < 0.0) {
    2266                directionOut_ = 1;    // to lower bound
    2267                theta_ = lowerOut_ - valueOut_;
    2268           } else {
    2269                directionOut_ = -1;    // to upper bound
    2270                theta_ = upperOut_ - valueOut_;
    2271           }
    2272      }
    2273 
    2274      double theta1 = CoinMax(theta_, 1.0e-12);
    2275      double theta2 = numberIterations_ * nonLinearCost_->averageTheta();
    2276      // Set average theta
    2277      nonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast<double> (numberIterations_ + 1)));
    2278      //if (numberIterations_%1000==0)
    2279      //printf("average theta is %g\n",nonLinearCost_->averageTheta());
    2280 
    2281      // clear arrays
    2282 
    2283      CoinZeroN(spare, numberRemaining);
    2284 
    2285      // put back original bounds etc
    2286      CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
    2287      rhsArray->setNumElements(numberRemaining);
    2288      rhsArray->setPacked();
    2289      nonLinearCost_->goBackAll(rhsArray);
    2290      rhsArray->clear();
    2291 
     2192      if (largestInfeasibility > primalTolerance_ && (handler_->logLevel() & 32) > -1)
     2193        printf("Primal tolerance increased from %g to %g\n",
     2194          primalTolerance_, largestInfeasibility);
     2195#endif
     2196      //#undef CLP_DEBUG
     2197      primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
     2198    }
     2199    // Need to look at all in some cases
     2200    if (theta_ > tentativeTheta) {
     2201      for (iIndex = 0; iIndex < number; iIndex++)
     2202        setActive(which[iIndex]);
     2203    }
     2204    if (way < 0.0)
     2205      theta_ = -theta_;
     2206    double newValue = valueOut_ - theta_ * alpha_;
     2207    // If  4 bit set - Force outgoing variables to exact bound (primal)
     2208    if (alpha_ * way < 0.0) {
     2209      directionOut_ = -1; // to upper bound
     2210      if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
     2211        upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
     2212      } else {
     2213        upperOut_ = newValue;
     2214      }
     2215    } else {
     2216      directionOut_ = 1; // to lower bound
     2217      if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
     2218        lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
     2219      } else {
     2220        lowerOut_ = newValue;
     2221      }
     2222    }
     2223    dualOut_ = reducedCost(sequenceOut_);
     2224  } else if (maximumMovement < 1.0e20) {
     2225    // flip
     2226    pivotRow_ = -2; // so we can tell its a flip
     2227    sequenceOut_ = sequenceIn_;
     2228    valueOut_ = valueIn_;
     2229    dualOut_ = dualIn_;
     2230    lowerOut_ = lowerIn_;
     2231    upperOut_ = upperIn_;
     2232    alpha_ = 0.0;
     2233    if (way < 0.0) {
     2234      directionOut_ = 1; // to lower bound
     2235      theta_ = lowerOut_ - valueOut_;
     2236    } else {
     2237      directionOut_ = -1; // to upper bound
     2238      theta_ = upperOut_ - valueOut_;
     2239    }
     2240  }
     2241
     2242  double theta1 = CoinMax(theta_, 1.0e-12);
     2243  double theta2 = numberIterations_ * nonLinearCost_->averageTheta();
     2244  // Set average theta
     2245  nonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast< double >(numberIterations_ + 1)));
     2246  //if (numberIterations_%1000==0)
     2247  //printf("average theta is %g\n",nonLinearCost_->averageTheta());
     2248
     2249  // clear arrays
     2250
     2251  CoinZeroN(spare, numberRemaining);
     2252
     2253  // put back original bounds etc
     2254  CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
     2255  rhsArray->setNumElements(numberRemaining);
     2256  rhsArray->setPacked();
     2257  nonLinearCost_->goBackAll(rhsArray);
     2258  rhsArray->clear();
    22922259}
    22932260/*
     
    22982265   For easy problems we can just choose one of the first columns we look at
    22992266*/
    2300 void
    2301 ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
    2302                                CoinIndexedVector * spareRow1,
    2303                                CoinIndexedVector * spareRow2,
    2304                                CoinIndexedVector * spareColumn1,
    2305                                CoinIndexedVector * spareColumn2)
     2267void ClpSimplexPrimal::primalColumn(CoinIndexedVector *updates,
     2268  CoinIndexedVector *spareRow1,
     2269  CoinIndexedVector *spareRow2,
     2270  CoinIndexedVector *spareColumn1,
     2271  CoinIndexedVector *spareColumn2)
    23062272{
    23072273
    2308      ClpMatrixBase * saveMatrix = matrix_;
    2309      double * saveRowScale = rowScale_;
    2310      if (scaledMatrix_) {
    2311           rowScale_ = NULL;
    2312           matrix_ = scaledMatrix_;
    2313      }
    2314      sequenceIn_ = primalColumnPivot_->pivotColumn(updates, spareRow1,
    2315                    spareRow2, spareColumn1,
    2316                    spareColumn2);
    2317      if (scaledMatrix_) {
    2318           matrix_ = saveMatrix;
    2319           rowScale_ = saveRowScale;
    2320      }
    2321      if (sequenceIn_ >= 0) {
    2322           valueIn_ = solution_[sequenceIn_];
    2323           dualIn_ = dj_[sequenceIn_];
    2324           if (nonLinearCost_->lookBothWays()) {
    2325                // double check
    2326                ClpSimplex::Status status = getStatus(sequenceIn_);
    2327 
    2328                switch(status) {
    2329                case ClpSimplex::atUpperBound:
    2330                     if (dualIn_ < 0.0) {
    2331                          // move to other side
    2332                          COIN_DETAIL_PRINT(printf("For %d U (%g, %g, %g) dj changed from %g",
    2333                                 sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
    2334                                                   upper_[sequenceIn_], dualIn_));
    2335                          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
    2336                          COIN_DETAIL_PRINT(printf(" to %g\n", dualIn_));
    2337                          nonLinearCost_->setOne(sequenceIn_, upper_[sequenceIn_] + 2.0 * currentPrimalTolerance());
    2338                          setStatus(sequenceIn_, ClpSimplex::atLowerBound);
    2339                     }
    2340                     break;
    2341                case ClpSimplex::atLowerBound:
    2342                     if (dualIn_ > 0.0) {
    2343                          // move to other side
    2344                          COIN_DETAIL_PRINT(printf("For %d L (%g, %g, %g) dj changed from %g",
    2345                                 sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
    2346                                                   upper_[sequenceIn_], dualIn_));
    2347                          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
    2348                          COIN_DETAIL_PRINT(printf(" to %g\n", dualIn_));
    2349                          nonLinearCost_->setOne(sequenceIn_, lower_[sequenceIn_] - 2.0 * currentPrimalTolerance());
    2350                          setStatus(sequenceIn_, ClpSimplex::atUpperBound);
    2351                     }
    2352                     break;
    2353                default:
    2354                     break;
    2355                }
    2356           }
    2357           lowerIn_ = lower_[sequenceIn_];
    2358           upperIn_ = upper_[sequenceIn_];
    2359           if (dualIn_ > 0.0)
    2360                directionIn_ = -1;
    2361           else
    2362                directionIn_ = 1;
    2363      } else {
    2364           sequenceIn_ = -1;
    2365      }
     2274  ClpMatrixBase *saveMatrix = matrix_;
     2275  double *saveRowScale = rowScale_;
     2276  if (scaledMatrix_) {
     2277    rowScale_ = NULL;
     2278    matrix_ = scaledMatrix_;
     2279  }
     2280  sequenceIn_ = primalColumnPivot_->pivotColumn(updates, spareRow1,
     2281    spareRow2, spareColumn1,
     2282    spareColumn2);
     2283  if (scaledMatrix_) {
     2284    matrix_ = saveMatrix;
     2285    rowScale_ = saveRowScale;
     2286  }
     2287  if (sequenceIn_ >= 0) {
     2288    valueIn_ = solution_[sequenceIn_];
     2289    dualIn_ = dj_[sequenceIn_];
     2290    if (nonLinearCost_->lookBothWays()) {
     2291      // double check
     2292      ClpSimplex::Status status = getStatus(sequenceIn_);
     2293
     2294      switch (status) {
     2295      case ClpSimplex::atUpperBound:
     2296        if (dualIn_ < 0.0) {
     2297          // move to other side
     2298          COIN_DETAIL_PRINT(printf("For %d U (%g, %g, %g) dj changed from %g",
     2299            sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
     2300            upper_[sequenceIn_], dualIn_));
     2301          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
     2302          COIN_DETAIL_PRINT(printf(" to %g\n", dualIn_));
     2303          nonLinearCost_->setOne(sequenceIn_, upper_[sequenceIn_] + 2.0 * currentPrimalTolerance());
     2304          setStatus(sequenceIn_, ClpSimplex::atLowerBound);
     2305        }
     2306        break;
     2307      case ClpSimplex::atLowerBound:
     2308        if (dualIn_ > 0.0) {
     2309          // move to other side
     2310          COIN_DETAIL_PRINT(printf("For %d L (%g, %g, %g) dj changed from %g",
     2311            sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
     2312            upper_[sequenceIn_], dualIn_));
     2313          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
     2314          COIN_DETAIL_PRINT(printf(" to %g\n", dualIn_));
     2315          nonLinearCost_->setOne(sequenceIn_, lower_[sequenceIn_] - 2.0 * currentPrimalTolerance());
     2316          setStatus(sequenceIn_, ClpSimplex::atUpperBound);
     2317        }
     2318        break;
     2319      default:
     2320        break;
     2321      }
     2322    }
     2323    lowerIn_ = lower_[sequenceIn_];
     2324    upperIn_ = upper_[sequenceIn_];
     2325    if (dualIn_ > 0.0)
     2326      directionIn_ = -1;
     2327    else
     2328      directionIn_ = 1;
     2329  } else {
     2330    sequenceIn_ = -1;
     2331  }
    23662332}
    23672333/* The primals are updated by the given array.
     
    23692335   After rowArray will have list of cost changes
    23702336*/
    2371 int
    2372 ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
    2373                                         double theta,
    2374                                         double & objectiveChange,
    2375                                         int valuesPass)
     2337int ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector *rowArray,
     2338  double theta,
     2339  double &objectiveChange,
     2340  int valuesPass)
    23762341{
    2377      // Cost on pivot row may change - may need to change dualIn
    2378      double oldCost = 0.0;
    2379      if (pivotRow_ >= 0)
    2380           oldCost = cost_[sequenceOut_];
    2381      //rowArray->scanAndPack();
    2382      double * work = rowArray->denseVector();
    2383      int number = rowArray->getNumElements();
    2384      int * which = rowArray->getIndices();
    2385 
    2386      int newNumber = 0;
    2387      int pivotPosition = -1;
    2388      nonLinearCost_->setChangeInCost(0.0);
    2389      //printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
    2390      //   solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
    2391      // allow for case where bound+tolerance == bound
    2392      //double tolerance = 0.999999*primalTolerance_;
    2393      double relaxedTolerance = 1.001 * primalTolerance_;
    2394      int iIndex;
    2395      if (!valuesPass) {
    2396           for (iIndex = 0; iIndex < number; iIndex++) {
    2397 
    2398                int iRow = which[iIndex];
    2399                double alpha = work[iIndex];
    2400                work[iIndex] = 0.0;
    2401                int iPivot = pivotVariable_[iRow];
    2402                double change = theta * alpha;
    2403                double value = solution_[iPivot] - change;
    2404                solution_[iPivot] = value;
     2342  // Cost on pivot row may change - may need to change dualIn
     2343  double oldCost = 0.0;
     2344  if (pivotRow_ >= 0)
     2345    oldCost = cost_[sequenceOut_];
     2346  //rowArray->scanAndPack();
     2347  double *work = rowArray->denseVector();
     2348  int number = rowArray->getNumElements();
     2349  int *which = rowArray->getIndices();
     2350
     2351  int newNumber = 0;
     2352  int pivotPosition = -1;
     2353  nonLinearCost_->setChangeInCost(0.0);
     2354  //printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
     2355  //   solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
     2356  // allow for case where bound+tolerance == bound
     2357  //double tolerance = 0.999999*primalTolerance_;
     2358  double relaxedTolerance = 1.001 * primalTolerance_;
     2359  int iIndex;
     2360  if (!valuesPass) {
     2361    for (iIndex = 0; iIndex < number; iIndex++) {
     2362
     2363      int iRow = which[iIndex];
     2364      double alpha = work[iIndex];
     2365      work[iIndex] = 0.0;
     2366      int iPivot = pivotVariable_[iRow];
     2367      double change = theta * alpha;
     2368      double value = solution_[iPivot] - change;
     2369      solution_[iPivot] = value;
    24052370#ifndef NDEBUG
    2406                // check if not active then okay
    2407                if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != -1) {
    2408                     // But make sure one going out is feasible
    2409                     if (change > 0.0) {
    2410                          // going down
    2411                          if (value <= lower_[iPivot] + primalTolerance_) {
    2412                               if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
    2413                                    value = lower_[iPivot];
    2414                               //double difference = nonLinearCost_->setOne(iPivot, value);
    2415                               //assert (!difference || fabs(change) > 1.0e9);
    2416                          }
    2417                     } else {
    2418                          // going up
    2419                          if (value >= upper_[iPivot] - primalTolerance_) {
    2420                               if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
    2421                                    value = upper_[iPivot];
    2422                               //double difference = nonLinearCost_->setOne(iPivot, value);
    2423                               //assert (!difference || fabs(change) > 1.0e9);
    2424                          }
    2425                     }
    2426                }
    2427 #endif
    2428                if (active(iRow) || theta_ < 0.0) {
    2429                     clearActive(iRow);
    2430                     // But make sure one going out is feasible
    2431                     if (change > 0.0) {
    2432                          // going down
    2433                          if (value <= lower_[iPivot] + primalTolerance_) {
    2434                               if (iPivot == sequenceOut_ && value >= lower_[iPivot] - relaxedTolerance)
    2435                                    value = lower_[iPivot];
    2436                               double difference = nonLinearCost_->setOne(iPivot, value);
    2437                               if (difference) {
    2438                                    if (iRow == pivotRow_)
    2439                                         pivotPosition = newNumber;
    2440                                    work[newNumber] = difference;
    2441                                    //change reduced cost on this
    2442                                    dj_[iPivot] = -difference;
    2443                                    which[newNumber++] = iRow;
    2444                               }
    2445                          }
    2446                     } else {
    2447                          // going up
    2448                          if (value >= upper_[iPivot] - primalTolerance_) {
    2449                               if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
    2450                                    value = upper_[iPivot];
    2451                               double difference = nonLinearCost_->setOne(iPivot, value);
    2452                               if (difference) {
    2453                                    if (iRow == pivotRow_)
    2454                                         pivotPosition = newNumber;
    2455                                    work[newNumber] = difference;
    2456                                    //change reduced cost on this
    2457                                    dj_[iPivot] = -difference;
    2458                                    which[newNumber++] = iRow;
    2459                               }
    2460                          }
    2461                     }
    2462                }
    2463           }
    2464      } else {
    2465           // values pass so look at all
    2466           for (iIndex = 0; iIndex < number; iIndex++) {
    2467 
    2468                int iRow = which[iIndex];
    2469                double alpha = work[iIndex];
    2470                work[iIndex] = 0.0;
    2471                int iPivot = pivotVariable_[iRow];
    2472                double change = theta * alpha;
    2473                double value = solution_[iPivot] - change;
    2474                solution_[iPivot] = value;
    2475                clearActive(iRow);
    2476                // But make sure one going out is feasible
    2477                if (change > 0.0) {
    2478                     // going down
    2479                     if (value <= lower_[iPivot] + primalTolerance_) {
    2480                          if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
    2481                               value = lower_[iPivot];
    2482                          double difference = nonLinearCost_->setOne(iPivot, value);
    2483                          if (difference) {
    2484                               if (iRow == pivotRow_)
    2485                                    pivotPosition = newNumber;
    2486                               work[newNumber] = difference;
    2487                               //change reduced cost on this
    2488                               dj_[iPivot] = -difference;
    2489                               which[newNumber++] = iRow;
    2490                          }
    2491                     }
    2492                } else {
    2493                     // going up
    2494                     if (value >= upper_[iPivot] - primalTolerance_) {
    2495                          if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
    2496                               value = upper_[iPivot];
    2497                          double difference = nonLinearCost_->setOne(iPivot, value);
    2498                          if (difference) {
    2499                               if (iRow == pivotRow_)
    2500                                    pivotPosition = newNumber;
    2501                               work[newNumber] = difference;
    2502                               //change reduced cost on this
    2503                               dj_[iPivot] = -difference;
    2504                               which[newNumber++] = iRow;
    2505                          }
    2506                     }
    2507                }
    2508           }
    2509      }
    2510      objectiveChange += nonLinearCost_->changeInCost();
    2511      rowArray->setPacked();
     2371      // check if not active then okay
     2372      if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != -1) {
     2373        // But make sure one going out is feasible
     2374        if (change > 0.0) {
     2375          // going down
     2376          if (value <= lower_[iPivot] + primalTolerance_) {
     2377            if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
     2378              value = lower_[iPivot];
     2379            //double difference = nonLinearCost_->setOne(iPivot, value);
     2380            //assert (!difference || fabs(change) > 1.0e9);
     2381          }
     2382        } else {
     2383          // going up
     2384          if (value >= upper_[iPivot] - primalTolerance_) {
     2385            if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
     2386              value = upper_[iPivot];
     2387            //double difference = nonLinearCost_->setOne(iPivot, value);
     2388            //assert (!difference || fabs(change) > 1.0e9);
     2389          }
     2390        }
     2391      }
     2392#endif
     2393      if (active(iRow) || theta_ < 0.0) {
     2394        clearActive(iRow);
     2395        // But make sure one going out is feasible
     2396        if (change > 0.0) {
     2397          // going down
     2398          if (value <= lower_[iPivot] + primalTolerance_) {
     2399            if (iPivot == sequenceOut_ && value >= lower_[iPivot] - relaxedTolerance)
     2400              value = lower_[iPivot];
     2401            double difference = nonLinearCost_->setOne(iPivot, value);
     2402            if (difference) {
     2403              if (iRow == pivotRow_)
     2404                pivotPosition = newNumber;
     2405              work[newNumber] = difference;
     2406              //change reduced cost on this
     2407              dj_[iPivot] = -difference;
     2408              which[newNumber++] = iRow;
     2409            }
     2410          }
     2411        } else {
     2412          // going up
     2413          if (value >= upper_[iPivot] - primalTolerance_) {
     2414            if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
     2415              value = upper_[iPivot];
     2416            double difference = nonLinearCost_->setOne(iPivot, value);
     2417            if (difference) {
     2418              if (iRow == pivotRow_)
     2419                pivotPosition = newNumber;
     2420              work[newNumber] = difference;
     2421              //change reduced cost on this
     2422              dj_[iPivot] = -difference;
     2423              which[newNumber++] = iRow;
     2424            }
     2425          }
     2426        }
     2427      }
     2428    }
     2429  } else {
     2430    // values pass so look at all
     2431    for (iIndex = 0; iIndex < number; iIndex++) {
     2432
     2433      int iRow = which[iIndex];
     2434      double alpha = work[iIndex];
     2435      work[iIndex] = 0.0;
     2436      int iPivot = pivotVariable_[iRow];
     2437      double change = theta * alpha;
     2438      double value = solution_[iPivot] - change;
     2439      solution_[iPivot] = value;
     2440      clearActive(iRow);
     2441      // But make sure one going out is feasible
     2442      if (change > 0.0) {
     2443        // going down
     2444        if (value <= lower_[iPivot] + primalTolerance_) {
     2445          if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
     2446            value = lower_[iPivot];
     2447          double difference = nonLinearCost_->setOne(iPivot, value);
     2448          if (difference) {
     2449            if (iRow == pivotRow_)
     2450              pivotPosition = newNumber;
     2451            work[newNumber] = difference;
     2452            //change reduced cost on this
     2453            dj_[iPivot] = -difference;
     2454            which[newNumber++] = iRow;
     2455          }
     2456        }
     2457      } else {
     2458        // going up
     2459        if (value >= upper_[iPivot] - primalTolerance_) {
     2460          if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
     2461            value = upper_[iPivot];
     2462          double difference = nonLinearCost_->setOne(iPivot, value);
     2463          if (difference) {
     2464            if (iRow == pivotRow_)
     2465              pivotPosition = newNumber;
     2466            work[newNumber] = difference;
     2467            //change reduced cost on this
     2468            dj_[iPivot] = -difference;
     2469            which[newNumber++] = iRow;
     2470          }
     2471        }
     2472      }
     2473    }
     2474  }
     2475  objectiveChange += nonLinearCost_->changeInCost();
     2476  rowArray->setPacked();
    25122477#if 0
    25132478     rowArray->setNumElements(newNumber);
     
    25242489     }
    25252490#else
    2526      if (pivotRow_ >= 0) {
    2527           double dualIn = dualIn_ + (oldCost - cost_[sequenceOut_]);
    2528           // update change vector to include pivot
    2529           if (pivotPosition >= 0) {
    2530                work[pivotPosition] -= dualIn;
    2531           } else {
    2532                work[newNumber] = -dualIn;
    2533                which[newNumber++] = pivotRow_;
    2534           }
    2535      }
    2536      rowArray->setNumElements(newNumber);
    2537 #endif
    2538      return 0;
     2491  if (pivotRow_ >= 0) {
     2492    double dualIn = dualIn_ + (oldCost - cost_[sequenceOut_]);
     2493    // update change vector to include pivot
     2494    if (pivotPosition >= 0) {
     2495      work[pivotPosition] -= dualIn;
     2496    } else {
     2497      work[newNumber] = -dualIn;
     2498      which[newNumber++] = pivotRow_;
     2499    }
     2500  }
     2501  rowArray->setNumElements(newNumber);
     2502#endif
     2503  return 0;
    25392504}
    25402505// Perturbs problem
    2541 void
    2542 ClpSimplexPrimal::perturb(int type)
     2506void ClpSimplexPrimal::perturb(int type)
    25432507{
    2544      if (perturbation_ > 100)
    2545           return; //perturbed already
    2546      if (perturbation_ == 100)
    2547           perturbation_ = 50; // treat as normal
    2548      int savePerturbation = perturbation_;
    2549      int i;
    2550      if (!numberIterations_)
    2551           cleanStatus(); // make sure status okay
    2552      // Make sure feasible bounds
    2553      if (nonLinearCost_) {
    2554        nonLinearCost_->checkInfeasibilities();
    2555        //nonLinearCost_->feasibleBounds();
    2556      }
    2557      // look at element range
    2558      double smallestNegative;
    2559      double largestNegative;
    2560      double smallestPositive;
    2561      double largestPositive;
    2562      matrix_->rangeOfElements(smallestNegative, largestNegative,
    2563                               smallestPositive, largestPositive);
    2564      smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
    2565      largestPositive = CoinMax(fabs(largestNegative), largestPositive);
    2566      double elementRatio = largestPositive / smallestPositive;
    2567      if (!numberIterations_ && perturbation_ == 50) {
    2568           // See if we need to perturb
    2569           int numberTotal = CoinMax(numberRows_, numberColumns_);
    2570           double * sort = new double[numberTotal];
    2571           int nFixed = 0;
    2572           for (i = 0; i < numberRows_; i++) {
    2573                double lo = fabs(rowLower_[i]);
    2574                double up = fabs(rowUpper_[i]);
    2575                double value = 0.0;
    2576                if (lo && lo < 1.0e20) {
    2577                     if (up && up < 1.0e20) {
    2578                          value = 0.5 * (lo + up);
    2579                          if (lo == up)
    2580                               nFixed++;
    2581                     } else {
    2582                          value = lo;
    2583                     }
    2584                } else {
    2585                     if (up && up < 1.0e20)
    2586                          value = up;
    2587                }
    2588                sort[i] = value;
    2589           }
    2590           std::sort(sort, sort + numberRows_);
    2591           int number = 1;
    2592           double last = sort[0];
    2593           for (i = 1; i < numberRows_; i++) {
    2594                if (last != sort[i])
    2595                     number++;
    2596                last = sort[i];
    2597           }
     2508  if (perturbation_ > 100)
     2509    return; //perturbed already
     2510  if (perturbation_ == 100)
     2511    perturbation_ = 50; // treat as normal
     2512  int savePerturbation = perturbation_;
     2513  int i;
     2514  if (!numberIterations_)
     2515    cleanStatus(); // make sure status okay
     2516  // Make sure feasible bounds
     2517  if (nonLinearCost_) {
     2518    nonLinearCost_->checkInfeasibilities();
     2519    //nonLinearCost_->feasibleBounds();
     2520  }
     2521  // look at element range
     2522  double smallestNegative;
     2523  double largestNegative;
     2524  double smallestPositive;
     2525  double largestPositive;
     2526  matrix_->rangeOfElements(smallestNegative, largestNegative,
     2527    smallestPositive, largestPositive);
     2528  smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
     2529  largestPositive = CoinMax(fabs(largestNegative), largestPositive);
     2530  double elementRatio = largestPositive / smallestPositive;
     2531  if (!numberIterations_ && perturbation_ == 50) {
     2532    // See if we need to perturb
     2533    int numberTotal = CoinMax(numberRows_, numberColumns_);
     2534    double *sort = new double[numberTotal];
     2535    int nFixed = 0;
     2536    for (i = 0; i < numberRows_; i++) {
     2537      double lo = fabs(rowLower_[i]);
     2538      double up = fabs(rowUpper_[i]);
     2539      double value = 0.0;
     2540      if (lo && lo < 1.0e20) {
     2541        if (up && up < 1.0e20) {
     2542          value = 0.5 * (lo + up);
     2543          if (lo == up)
     2544            nFixed++;
     2545        } else {
     2546          value = lo;
     2547        }
     2548      } else {
     2549        if (up && up < 1.0e20)
     2550          value = up;
     2551      }
     2552      sort[i] = value;
     2553    }
     2554    std::sort(sort, sort + numberRows_);
     2555    int number = 1;
     2556    double last = sort[0];
     2557    for (i = 1; i < numberRows_; i++) {
     2558      if (last != sort[i])
     2559        number++;
     2560      last = sort[i];
     2561    }
    25982562#ifdef KEEP_GOING_IF_FIXED
    2599           //printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
    2600           //   numberRows_,nFixed,elementRatio);
    2601 #endif
    2602           if (number * 4 > numberRows_ || elementRatio > 1.0e12) {
    2603                perturbation_ = 100;
    2604                delete [] sort;
    2605                return; // good enough
    2606           }
    2607           number = 0;
     2563    //printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
     2564    //   numberRows_,nFixed,elementRatio);
     2565#endif
     2566    if (number * 4 > numberRows_ || elementRatio > 1.0e12) {
     2567      perturbation_ = 100;
     2568      delete[] sort;
     2569      return; // good enough
     2570    }
     2571    number = 0;
    26082572#ifdef KEEP_GOING_IF_FIXED
    2609           if (!integerType_) {
    2610                // look at columns
    2611                nFixed = 0;
    2612                for (i = 0; i < numberColumns_; i++) {
    2613                     double lo = fabs(columnLower_[i]);
    2614                     double up = fabs(columnUpper_[i]);
    2615                     double value = 0.0;
    2616                     if (lo && lo < 1.0e20) {
    2617                          if (up && up < 1.0e20) {
    2618                               value = 0.5 * (lo + up);
    2619                               if (lo == up)
    2620                                    nFixed++;
    2621                          } else {
    2622                               value = lo;
    2623                          }
    2624                     } else {
    2625                          if (up && up < 1.0e20)
    2626                               value = up;
    2627                     }
    2628                     sort[i] = value;
    2629                }
    2630                std::sort(sort, sort + numberColumns_);
    2631                number = 1;
    2632                last = sort[0];
    2633                for (i = 1; i < numberColumns_; i++) {
    2634                     if (last != sort[i])
    2635                          number++;
    2636                     last = sort[i];
    2637                }
    2638                //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
    2639                //     numberColumns_,nFixed);
    2640           }
    2641 #endif
    2642           delete [] sort;
    2643           if (number * 4 > numberColumns_) {
    2644                perturbation_ = 100;
    2645                return; // good enough
    2646           }
    2647      }
    2648      // primal perturbation
    2649      double perturbation = 1.0e-20;
    2650      double bias = 1.0;
    2651      int numberNonZero = 0;
    2652      // maximum fraction of rhs/bounds to perturb
    2653      double maximumFraction = 1.0e-5;
    2654      if (perturbation_ >= 50) {
    2655           perturbation = 1.0e-4;
    2656           for (i = 0; i < numberColumns_ + numberRows_; i++) {
    2657                if (upper_[i] > lower_[i] + primalTolerance_) {
    2658                     double lowerValue, upperValue;
    2659                     if (lower_[i] > -1.0e20)
    2660                          lowerValue = fabs(lower_[i]);
    2661                     else
    2662                          lowerValue = 0.0;
    2663                     if (upper_[i] < 1.0e20)
    2664                          upperValue = fabs(upper_[i]);
    2665                     else
    2666                          upperValue = 0.0;
    2667                     double value = CoinMax(fabs(lowerValue), fabs(upperValue));
    2668                     value = CoinMin(value, upper_[i] - lower_[i]);
     2573    if (!integerType_) {
     2574      // look at columns
     2575      nFixed = 0;
     2576      for (i = 0; i < numberColumns_; i++) {
     2577        double lo = fabs(columnLower_[i]);
     2578        double up = fabs(columnUpper_[i]);
     2579        double value = 0.0;
     2580        if (lo && lo < 1.0e20) {
     2581          if (up && up < 1.0e20) {
     2582            value = 0.5 * (lo + up);
     2583            if (lo == up)
     2584              nFixed++;
     2585          } else {
     2586            value = lo;
     2587          }
     2588        } else {
     2589          if (up && up < 1.0e20)
     2590            value = up;
     2591        }
     2592        sort[i] = value;
     2593      }
     2594      std::sort(sort, sort + numberColumns_);
     2595      number = 1;
     2596      last = sort[0];
     2597      for (i = 1; i < numberColumns_; i++) {
     2598        if (last != sort[i])
     2599          number++;
     2600        last = sort[i];
     2601      }
     2602      //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
     2603      //     numberColumns_,nFixed);
     2604    }
     2605#endif
     2606    delete[] sort;
     2607    if (number * 4 > numberColumns_) {
     2608      perturbation_ = 100;
     2609      return; // good enough
     2610    }
     2611  }
     2612  // primal perturbation
     2613  double perturbation = 1.0e-20;
     2614  double bias = 1.0;
     2615  int numberNonZero = 0;
     2616  // maximum fraction of rhs/bounds to perturb
     2617  double maximumFraction = 1.0e-5;
     2618  if (perturbation_ >= 50) {
     2619    perturbation = 1.0e-4;
     2620    for (i = 0; i < numberColumns_ + numberRows_; i++) {
     2621      if (upper_[i] > lower_[i] + primalTolerance_) {
     2622        double lowerValue, upperValue;
     2623        if (lower_[i] > -1.0e20)
     2624          lowerValue = fabs(lower_[i]);
     2625        else
     2626          lowerValue = 0.0;
     2627        if (upper_[i] < 1.0e20)
     2628          upperValue = fabs(upper_[i]);
     2629        else
     2630          upperValue = 0.0;
     2631        double value = CoinMax(fabs(lowerValue), fabs(upperValue));
     2632        value = CoinMin(value, upper_[i] - lower_[i]);
    26692633#if 1
    2670                     if (value) {
    2671                          perturbation += value;
    2672                          numberNonZero++;
    2673                     }
     2634        if (value) {
     2635          perturbation += value;
     2636          numberNonZero++;
     2637        }
    26742638#else
    2675                     perturbation = CoinMax(perturbation, value);
    2676 #endif
    2677                }
    2678           }
    2679           if (numberNonZero)
    2680                perturbation /= static_cast<double> (numberNonZero);
    2681           else
    2682                perturbation = 1.0e-1;
    2683           if (perturbation_ > 50 && perturbation_ < 55) {
    2684                // reduce
    2685                while (perturbation_ < 55) {
    2686                     perturbation_++;
    2687                     perturbation *= 0.25;
    2688                     bias *= 0.25;
    2689                }
    2690           } else if (perturbation_ >= 55 && perturbation_ < 60) {
    2691                // increase
    2692                while (perturbation_ > 55) {
    2693                     perturbation_--;
    2694                     perturbation *= 4.0;
    2695                }
    2696                perturbation_ = 50;
    2697           }
    2698      } else if (perturbation_ < 100) {
    2699           perturbation = pow(10.0, perturbation_);
    2700           // user is in charge
    2701           maximumFraction = 1.0;
    2702      }
    2703      double largestZero = 0.0;
    2704      double largest = 0.0;
    2705      double largestPerCent = 0.0;
    2706      bool printOut = (handler_->logLevel() == 63);
    2707      printOut = false; //off
    2708      // Check if all slack
    2709      int number = 0;
    2710      int iSequence;
    2711      for (iSequence = 0; iSequence < numberRows_; iSequence++) {
    2712           if (getRowStatus(iSequence) == basic)
    2713                number++;
    2714      }
    2715      if (rhsScale_ > 100.0) {
    2716           // tone down perturbation
    2717           maximumFraction *= 0.1;
    2718      }
    2719      if (savePerturbation==51) {
    2720        perturbation = CoinMin(0.1,perturbation);
    2721        maximumFraction *=0.1;
    2722      }
    2723      if (number != numberRows_)
    2724           type = 1;
    2725      // modify bounds
    2726      // Change so at least 1.0e-5 and no more than 0.1
    2727      // For now just no more than 0.1
    2728      // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
    2729      // seems much slower???#define SAVE_PERT
     2639        perturbation = CoinMax(perturbation, value);
     2640#endif
     2641      }
     2642    }
     2643    if (numberNonZero)
     2644      perturbation /= static_cast< double >(numberNonZero);
     2645    else
     2646      perturbation = 1.0e-1;
     2647    if (perturbation_ > 50 && perturbation_ < 55) {
     2648      // reduce
     2649      while (perturbation_ < 55) {
     2650        perturbation_++;
     2651        perturbation *= 0.25;
     2652        bias *= 0.25;
     2653      }
     2654    } else if (perturbation_ >= 55 && perturbation_ < 60) {
     2655      // increase
     2656      while (perturbation_ > 55) {
     2657        perturbation_--;
     2658        perturbation *= 4.0;
     2659      }
     2660      perturbation_ = 50;
     2661    }
     2662  } else if (perturbation_ < 100) {
     2663    perturbation = pow(10.0, perturbation_);
     2664    // user is in charge
     2665    maximumFraction = 1.0;
     2666  }
     2667  double largestZero = 0.0;
     2668  double largest = 0.0;
     2669  double largestPerCent = 0.0;
     2670  bool printOut = (handler_->logLevel() == 63);
     2671  printOut = false; //off
     2672  // Check if all slack
     2673  int number = 0;
     2674  int iSequence;
     2675  for (iSequence = 0; iSequence < numberRows_; iSequence++) {
     2676    if (getRowStatus(iSequence) == basic)
     2677      number++;
     2678  }
     2679  if (rhsScale_ > 100.0) {
     2680    // tone down perturbation
     2681    maximumFraction *= 0.1;
     2682  }
     2683  if (savePerturbation == 51) {
     2684    perturbation = CoinMin(0.1, perturbation);
     2685    maximumFraction *= 0.1;
     2686  }
     2687  if (number != numberRows_)
     2688    type = 1;
     2689    // modify bounds
     2690    // Change so at least 1.0e-5 and no more than 0.1
     2691    // For now just no more than 0.1
     2692    // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
     2693    // seems much slower???#define SAVE_PERT
    27302694#ifdef SAVE_PERT
    2731      if (2 * numberColumns_ > maximumPerturbationSize_) {
    2732           delete [] perturbationArray_;
    2733           maximumPerturbationSize_ = 2 * numberColumns_;
    2734           perturbationArray_ = new double [maximumPerturbationSize_];
    2735           for (int iColumn = 0; iColumn < maximumPerturbationSize_; iColumn++) {
    2736                perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
    2737           }
    2738      }
    2739 #endif
    2740      if (type == 1) {
    2741           double tolerance = 100.0 * primalTolerance_;
    2742           tolerance = 10.0 * primalTolerance_; // try smaller
    2743           //double multiplier = perturbation*maximumFraction;
    2744           for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
    2745                if (getStatus(iSequence) == basic) {
    2746                     double lowerValue = lower_[iSequence];
    2747                     double upperValue = upper_[iSequence];
    2748                     if (upperValue > lowerValue + tolerance) {
    2749                          double solutionValue = solution_[iSequence];
    2750                          double difference = upperValue - lowerValue;
    2751                          difference = CoinMin(difference, perturbation);
    2752                          difference = CoinMin(difference, fabs(solutionValue) + 1.0);
    2753                          double value = maximumFraction * (difference + bias);
    2754                          value = CoinMin(value, 0.1);
    2755                          value = CoinMax(value,primalTolerance_);
     2695  if (2 * numberColumns_ > maximumPerturbationSize_) {
     2696    delete[] perturbationArray_;
     2697    maximumPerturbationSize_ = 2 * numberColumns_;
     2698    perturbationArray_ = new double[maximumPerturbationSize_];
     2699    for (int iColumn = 0; iColumn < maximumPerturbationSize_; iColumn++) {
     2700      perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
     2701    }
     2702  }
     2703#endif
     2704  if (type == 1) {
     2705    double tolerance = 100.0 * primalTolerance_;
     2706    tolerance = 10.0 * primalTolerance_; // try smaller
     2707    //double multiplier = perturbation*maximumFraction;
     2708    for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
     2709      if (getStatus(iSequence) == basic) {
     2710        double lowerValue = lower_[iSequence];
     2711        double upperValue = upper_[iSequence];
     2712        if (upperValue > lowerValue + tolerance) {
     2713          double solutionValue = solution_[iSequence];
     2714          double difference = upperValue - lowerValue;
     2715          difference = CoinMin(difference, perturbation);
     2716          difference = CoinMin(difference, fabs(solutionValue) + 1.0);
     2717          double value = maximumFraction * (difference + bias);
     2718          value = CoinMin(value, 0.1);
     2719          value = CoinMax(value, primalTolerance_);
    27562720#ifndef SAVE_PERT
    2757                          value *= randomNumberGenerator_.randomDouble();
     2721          value *= randomNumberGenerator_.randomDouble();
    27582722#else
    2759                          value *= perturbationArray_[2*iSequence];
    2760 #endif
    2761                         if (value) {
    2762                            while (value<tolerance)
    2763                              value *= 3.0;
    2764                         }
    2765                          if (solutionValue - lowerValue <= primalTolerance_) {
    2766                               lower_[iSequence] -= value;
    2767                          } else if (upperValue - solutionValue <= primalTolerance_) {
    2768                               upper_[iSequence] += value;
    2769                          } else {
     2723          value *= perturbationArray_[2 * iSequence];
     2724#endif
     2725          if (value) {
     2726            while (value < tolerance)
     2727              value *= 3.0;
     2728          }
     2729          if (solutionValue - lowerValue <= primalTolerance_) {
     2730            lower_[iSequence] -= value;
     2731          } else if (upperValue - solutionValue <= primalTolerance_) {
     2732            upper_[iSequence] += value;
     2733          } else {
    27702734#if 0
    27712735                              if (iSequence >= numberColumns_) {
     
    27792743                              }
    27802744#else
    2781                               value = 0.0;
    2782 #endif
    2783                          }
    2784                          if (value) {
    2785                               if (printOut)
    2786                                    printf("col %d lower from %g to %g, upper from %g to %g\n",
    2787                                           iSequence, lowerValue, lower_[iSequence], upperValue , upper_[iSequence]);
    2788                               if (solutionValue) {
    2789                                    largest = CoinMax(largest, value);
    2790                                    if (value > (fabs(solutionValue) + 1.0)*largestPerCent)
    2791                                         largestPerCent = value / (fabs(solutionValue) + 1.0);
    2792                               } else {
    2793                                    largestZero = CoinMax(largestZero, value);
    2794                               }
    2795                          }
    2796                     }
    2797                }
    2798           }
    2799      } else {
    2800           double tolerance = 100.0 * primalTolerance_;
    2801           tolerance = 10.0 * primalTolerance_; // try smaller
    2802           for (i = 0; i < numberColumns_; i++) {
    2803                double lowerValue = lower_[i], upperValue = upper_[i];
    2804                if (upperValue > lowerValue + primalTolerance_) {
    2805                     double value = perturbation * maximumFraction;
    2806                     value = CoinMin(value, 0.1);
     2745            value = 0.0;
     2746#endif
     2747          }
     2748          if (value) {
     2749            if (printOut)
     2750              printf("col %d lower from %g to %g, upper from %g to %g\n",
     2751                iSequence, lowerValue, lower_[iSequence], upperValue, upper_[iSequence]);
     2752            if (solutionValue) {
     2753              largest = CoinMax(largest, value);
     2754              if (value > (fabs(solutionValue) + 1.0) * largestPerCent)
     2755                largestPerCent = value / (fabs(solutionValue) + 1.0);
     2756            } else {
     2757              largestZero = CoinMax(largestZero, value);
     2758            }
     2759          }
     2760        }
     2761      }
     2762    }
     2763  } else {
     2764    double tolerance = 100.0 * primalTolerance_;
     2765    tolerance = 10.0 * primalTolerance_; // try smaller
     2766    for (i = 0; i < numberColumns_; i++) {
     2767      double lowerValue = lower_[i], upperValue = upper_[i];
     2768      if (upperValue > lowerValue + primalTolerance_) {
     2769        double value = perturbation * maximumFraction;
     2770        value = CoinMin(value, 0.1);
    28072771#ifndef SAVE_PERT
    2808                     value *= randomNumberGenerator_.randomDouble();
     2772        value *= randomNumberGenerator_.randomDouble();
    28092773#else
    2810                     value *= perturbationArray_[2*i+1];
    2811 #endif
    2812                     value *= randomNumberGenerator_.randomDouble();
    2813                     if (savePerturbation != 50) {
    2814                       if (fabs(value) <= primalTolerance_)
    2815                         value = 0.0;
    2816                     }
    2817                     if (value) {
    2818                          double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2819                          // get in range
    2820                          if (valueL <= tolerance) {
    2821                               valueL *= 10.0;
    2822                               while (valueL <= tolerance)
    2823                                    valueL *= 10.0;
    2824                          } else if (valueL > 1.0e-3) {
    2825                               valueL *= 0.1;
    2826                               while (valueL > 1.0e-3)
    2827                                    valueL *= 0.1;
    2828                          }
    2829                          if (lowerValue > -1.0e20 && lowerValue)
    2830                               lowerValue -= valueL;
    2831                          double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
    2832                          // get in range
    2833                          if (valueU <= tolerance) {
    2834                               valueU *= 10.0;
    2835                               while (valueU <= tolerance)
    2836                                    valueU *= 10.0;
    2837                          } else if (valueU > 1.0e-3) {
    2838                               valueU *= 0.1;
    2839                               while (valueU > 1.0e-3)
    2840                                    valueU *= 0.1;
    2841                          }
    2842                          if (upperValue < 1.0e20 && upperValue)
    2843                               upperValue += valueU;
    2844                     }
    2845                     if (lowerValue != lower_[i]) {
    2846                          double difference = fabs(lowerValue - lower_[i]);
    2847                          largest = CoinMax(largest, difference);
    2848                          if (difference > fabs(lower_[i])*largestPerCent)
    2849                               largestPerCent = fabs(difference / lower_[i]);
    2850                     }
    2851                     if (upperValue != upper_[i]) {
    2852                          double difference = fabs(upperValue - upper_[i]);
    2853                          largest = CoinMax(largest, difference);
    2854                          if (difference > fabs(upper_[i])*largestPerCent)
    2855                               largestPerCent = fabs(difference / upper_[i]);
    2856                     }
    2857                     if (printOut)
    2858                          printf("col %d lower from %g to %g, upper from %g to %g\n",
    2859                                 i, lower_[i], lowerValue, upper_[i], upperValue);
    2860                }
    2861                lower_[i] = lowerValue;
    2862                upper_[i] = upperValue;
    2863           }
    2864           // biased
    2865           const double * rowLower = rowLower_-numberColumns_;
    2866           const double * rowUpper = rowUpper_-numberColumns_;
    2867           for (; i < numberColumns_ + numberRows_; i++) {
    2868                double lowerValue = lower_[i], upperValue = upper_[i];
    2869                double value = perturbation * maximumFraction;
    2870                value = CoinMin(value, 0.1);
    2871                value *= randomNumberGenerator_.randomDouble();
    2872                if (rowLower[i]!=rowUpper[i] &&
    2873                    upperValue > lowerValue + tolerance) {
    2874                     if (savePerturbation != 50) {
    2875                          if (fabs(value) <= primalTolerance_)
    2876                               value = 0.0;
    2877                          if (lowerValue > -1.0e20 && lowerValue)
    2878                               lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2879                          if (upperValue < 1.0e20 && upperValue)
    2880                               upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
    2881                     } else if (value) {
    2882                          double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2883                          // get in range
    2884                          if (valueL <= tolerance) {
    2885                               valueL *= 10.0;
    2886                               while (valueL <= tolerance)
    2887                                    valueL *= 10.0;
    2888                          } else if (valueL > 1.0) {
    2889                               valueL *= 0.1;
    2890                               while (valueL > 1.0)
    2891                                    valueL *= 0.1;
    2892                          }
    2893                          if (lowerValue > -1.0e20 && lowerValue)
    2894                               lowerValue -= valueL;
    2895                          double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
    2896                          // get in range
    2897                          if (valueU <= tolerance) {
    2898                               valueU *= 10.0;
    2899                               while (valueU <= tolerance)
    2900                                    valueU *= 10.0;
    2901                          } else if (valueU > 1.0) {
    2902                               valueU *= 0.1;
    2903                               while (valueU > 1.0)
    2904                                    valueU *= 0.1;
    2905                          }
    2906                          if (upperValue < 1.0e20 && upperValue)
    2907                               upperValue += valueU;
    2908                     }
    2909                } else if (upperValue > 0.0) {
    2910                  //upperValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2911                  // lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2912                } else if (upperValue < 0.0) {
    2913                  // upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2914                  // lowerValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2915                } else {
    2916                }
    2917                if (lowerValue != lower_[i]) {
    2918                     double difference = fabs(lowerValue - lower_[i]);
    2919                     largest = CoinMax(largest, difference);
    2920                     if (difference > fabs(lower_[i])*largestPerCent)
    2921                          largestPerCent = fabs(difference / lower_[i]);
    2922                }
    2923                if (upperValue != upper_[i]) {
    2924                     double difference = fabs(upperValue - upper_[i]);
    2925                     largest = CoinMax(largest, difference);
    2926                     if (difference > fabs(upper_[i])*largestPerCent)
    2927                          largestPerCent = fabs(difference / upper_[i]);
    2928                }
    2929                if (printOut)
    2930                     printf("row %d lower from %g to %g, upper from %g to %g\n",
    2931                            i - numberColumns_, lower_[i], lowerValue, upper_[i], upperValue);
    2932                lower_[i] = lowerValue;
    2933                upper_[i] = upperValue;
    2934           }
    2935      }
    2936      // Clean up
    2937      for (i = 0; i < numberColumns_ + numberRows_; i++) {
    2938           switch(getStatus(i)) {
    2939 
    2940           case basic:
    2941                break;
    2942           case atUpperBound:
    2943                solution_[i] = upper_[i];
    2944                break;
    2945           case isFixed:
    2946           case atLowerBound:
    2947                solution_[i] = lower_[i];
    2948                break;
    2949           case isFree:
    2950                break;
    2951           case superBasic:
    2952                break;
    2953           }
    2954      }
    2955      if (largest || largestZero) {
    2956        handler_->message(CLP_SIMPLEX_PERTURB, messages_)
    2957                << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
    2958                << CoinMessageEol;
    2959        // say perturbed
    2960        perturbation_ = 101;
    2961      } else {
    2962        // say no need
    2963        perturbation_ = 102;
    2964      }
     2774        value *= perturbationArray_[2 * i + 1];
     2775#endif
     2776        value *= randomNumberGenerator_.randomDouble();
     2777        if (savePerturbation != 50) {
     2778          if (fabs(value) <= primalTolerance_)
     2779            value = 0.0;
     2780        }
     2781        if (value) {
     2782          double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2783          // get in range
     2784          if (valueL <= tolerance) {
     2785            valueL *= 10.0;
     2786            while (valueL <= tolerance)
     2787              valueL *= 10.0;
     2788          } else if (valueL > 1.0e-3) {
     2789            valueL *= 0.1;
     2790            while (valueL > 1.0e-3)
     2791              valueL *= 0.1;
     2792          }
     2793          if (lowerValue > -1.0e20 && lowerValue)
     2794            lowerValue -= valueL;
     2795          double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
     2796          // get in range
     2797          if (valueU <= tolerance) {
     2798            valueU *= 10.0;
     2799            while (valueU <= tolerance)
     2800              valueU *= 10.0;
     2801          } else if (valueU > 1.0e-3) {
     2802            valueU *= 0.1;
     2803            while (valueU > 1.0e-3)
     2804              valueU *= 0.1;
     2805          }
     2806          if (upperValue < 1.0e20 && upperValue)
     2807            upperValue += valueU;
     2808        }
     2809        if (lowerValue != lower_[i]) {
     2810          double difference = fabs(lowerValue - lower_[i]);
     2811          largest = CoinMax(largest, difference);
     2812          if (difference > fabs(lower_[i]) * largestPerCent)
     2813            largestPerCent = fabs(difference / lower_[i]);
     2814        }
     2815        if (upperValue != upper_[i]) {
     2816          double difference = fabs(upperValue - upper_[i]);
     2817          largest = CoinMax(largest, difference);
     2818          if (difference > fabs(upper_[i]) * largestPerCent)
     2819            largestPerCent = fabs(difference / upper_[i]);
     2820        }
     2821        if (printOut)
     2822          printf("col %d lower from %g to %g, upper from %g to %g\n",
     2823            i, lower_[i], lowerValue, upper_[i], upperValue);
     2824      }
     2825      lower_[i] = lowerValue;
     2826      upper_[i] = upperValue;
     2827    }
     2828    // biased
     2829    const double *rowLower = rowLower_ - numberColumns_;
     2830    const double *rowUpper = rowUpper_ - numberColumns_;
     2831    for (; i < numberColumns_ + numberRows_; i++) {
     2832      double lowerValue = lower_[i], upperValue = upper_[i];
     2833      double value = perturbation * maximumFraction;
     2834      value = CoinMin(value, 0.1);
     2835      value *= randomNumberGenerator_.randomDouble();
     2836      if (rowLower[i] != rowUpper[i] && upperValue > lowerValue + tolerance) {
     2837        if (savePerturbation != 50) {
     2838          if (fabs(value) <= primalTolerance_)
     2839            value = 0.0;
     2840          if (lowerValue > -1.0e20 && lowerValue)
     2841            lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2842          if (upperValue < 1.0e20 && upperValue)
     2843            upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
     2844        } else if (value) {
     2845          double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2846          // get in range
     2847          if (valueL <= tolerance) {
     2848            valueL *= 10.0;
     2849            while (valueL <= tolerance)
     2850              valueL *= 10.0;
     2851          } else if (valueL > 1.0) {
     2852            valueL *= 0.1;
     2853            while (valueL > 1.0)
     2854              valueL *= 0.1;
     2855          }
     2856          if (lowerValue > -1.0e20 && lowerValue)
     2857            lowerValue -= valueL;
     2858          double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
     2859          // get in range
     2860          if (valueU <= tolerance) {
     2861            valueU *= 10.0;
     2862            while (valueU <= tolerance)
     2863              valueU *= 10.0;
     2864          } else if (valueU > 1.0) {
     2865            valueU *= 0.1;
     2866            while (valueU > 1.0)
     2867              valueU *= 0.1;
     2868          }
     2869          if (upperValue < 1.0e20 && upperValue)
     2870            upperValue += valueU;
     2871        }
     2872      } else if (upperValue > 0.0) {
     2873        //upperValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2874        // lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2875      } else if (upperValue < 0.0) {
     2876        // upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2877        // lowerValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
     2878      } else {
     2879      }
     2880      if (lowerValue != lower_[i]) {
     2881        double difference = fabs(lowerValue - lower_[i]);
     2882        largest = CoinMax(largest, difference);
     2883        if (difference > fabs(lower_[i]) * largestPerCent)
     2884          largestPerCent = fabs(difference / lower_[i]);
     2885      }
     2886      if (upperValue != upper_[i]) {
     2887        double difference = fabs(upperValue - upper_[i]);
     2888        largest = CoinMax(largest, difference);
     2889        if (difference > fabs(upper_[i]) * largestPerCent)
     2890          largestPerCent = fabs(difference / upper_[i]);
     2891      }
     2892      if (printOut)
     2893        printf("row %d lower from %g to %g, upper from %g to %g\n",
     2894          i - numberColumns_, lower_[i], lowerValue, upper_[i], upperValue);
     2895      lower_[i] = lowerValue;
     2896      upper_[i] = upperValue;
     2897    }
     2898  }
     2899  // Clean up
     2900  for (i = 0; i < numberColumns_ + numberRows_; i++) {
     2901    switch (getStatus(i)) {
     2902
     2903    case basic:
     2904      break;
     2905    case atUpperBound:
     2906      solution_[i] = upper_[i];
     2907      break;
     2908    case isFixed:
     2909    case atLowerBound:
     2910      solution_[i] = lower_[i];
     2911      break;
     2912    case isFree:
     2913      break;
     2914    case superBasic:
     2915      break;
     2916    }
     2917  }
     2918  if (largest || largestZero) {
     2919    handler_->message(CLP_SIMPLEX_PERTURB, messages_)
     2920      << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
     2921      << CoinMessageEol;
     2922    // say perturbed
     2923    perturbation_ = 101;
     2924  } else {
     2925    // say no need
     2926    perturbation_ = 102;
     2927  }
    29652928}
    29662929// un perturb
    2967 bool
    2968 ClpSimplexPrimal::unPerturb()
     2930bool ClpSimplexPrimal::unPerturb()
    29692931{
    2970      if (perturbation_ != 101)
    2971           return false;
    2972      // put back original bounds and costs
    2973      createRim(1 + 4);
    2974      sanityCheck();
    2975      // unflag
    2976      unflag();
    2977      // get a valid nonlinear cost function
    2978      delete nonLinearCost_;
    2979      nonLinearCost_ = new ClpNonLinearCost(this);
    2980      perturbation_ = 102; // stop any further perturbation
    2981      // move non basic variables to new bounds
    2982      nonLinearCost_->checkInfeasibilities(0.0);
     2932  if (perturbation_ != 101)
     2933    return false;
     2934  // put back original bounds and costs
     2935  createRim(1 + 4);
     2936  sanityCheck();
     2937  // unflag
     2938  unflag();
     2939  // get a valid nonlinear cost function
     2940  delete nonLinearCost_;
     2941  nonLinearCost_ = new ClpNonLinearCost(this);
     2942  perturbation_ = 102; // stop any further perturbation
     2943  // move non basic variables to new bounds
     2944  nonLinearCost_->checkInfeasibilities(0.0);
    29832945#if 1
    2984      // Try using dual
    2985      return true;
     2946  // Try using dual
     2947  return true;
    29862948#else
    2987      gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    2988      return false;
    2989 #endif
    2990 
     2949  gutsOfSolution(NULL, NULL, ifValuesPass != 0);
     2950  return false;
     2951#endif
    29912952}
    29922953// Unflag all variables and return number unflagged
    2993 int
    2994 ClpSimplexPrimal::unflag()
     2954int ClpSimplexPrimal::unflag()
    29952955{
    2996      int i;
    2997      int number = numberRows_ + numberColumns_;
    2998      int numberFlagged = 0;
    2999      // we can't really trust infeasibilities if there is dual error
    3000      // allow tolerance bigger than standard to check on duals
    3001      double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
    3002      for (i = 0; i < number; i++) {
    3003           if (flagged(i)) {
    3004                clearFlagged(i);
    3005                // only say if reasonable dj
    3006                if (fabs(dj_[i]) > relaxedToleranceD)
    3007                     numberFlagged++;
    3008           }
    3009      }
    3010      numberFlagged += matrix_->generalExpanded(this, 8, i);
    3011      if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
    3012           printf("%d unflagged\n", numberFlagged);
    3013      return numberFlagged;
     2956  int i;
     2957  int number = numberRows_ + numberColumns_;
     2958  int numberFlagged = 0;
     2959  // we can't really trust infeasibilities if there is dual error
     2960  // allow tolerance bigger than standard to check on duals
     2961  double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
     2962  for (i = 0; i < number; i++) {
     2963    if (flagged(i)) {
     2964      clearFlagged(i);
     2965      // only say if reasonable dj
     2966      if (fabs(dj_[i]) > relaxedToleranceD)
     2967        numberFlagged++;
     2968    }
     2969  }
     2970  numberFlagged += matrix_->generalExpanded(this, 8, i);
     2971  if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
     2972    printf("%d unflagged\n", numberFlagged);
     2973  return numberFlagged;
    30142974}
    30152975// Do not change infeasibility cost and always say optimal
    3016 void
    3017 ClpSimplexPrimal::alwaysOptimal(bool onOff)
     2976void ClpSimplexPrimal::alwaysOptimal(bool onOff)
    30182977{
    3019      if (onOff)
    3020           specialOptions_ |= 1;
    3021      else
    3022           specialOptions_ &= ~1;
     2978  if (onOff)
     2979    specialOptions_ |= 1;
     2980  else
     2981    specialOptions_ &= ~1;
    30232982}
    3024 bool
    3025 ClpSimplexPrimal::alwaysOptimal() const
     2983bool ClpSimplexPrimal::alwaysOptimal() const
    30262984{
    3027      return (specialOptions_ & 1) != 0;
     2985  return (specialOptions_ & 1) != 0;
    30282986}
    30292987// Flatten outgoing variables i.e. - always to exact bound
    3030 void
    3031 ClpSimplexPrimal::exactOutgoing(bool onOff)
     2988void ClpSimplexPrimal::exactOutgoing(bool onOff)
    30322989{
    3033      if (onOff)
    3034           specialOptions_ |= 4;
    3035      else
    3036           specialOptions_ &= ~4;
     2990  if (onOff)
     2991    specialOptions_ |= 4;
     2992  else
     2993    specialOptions_ &= ~4;
    30372994}
    3038 bool
    3039 ClpSimplexPrimal::exactOutgoing() const
     2995bool ClpSimplexPrimal::exactOutgoing() const
    30402996{
    3041      return (specialOptions_ & 4) != 0;
     2997  return (specialOptions_ & 4) != 0;
    30422998}
    30432999#if ALT_UPDATE_WEIGHTS
    3044 static void doAlternate(ClpSimplex * model)
     3000static void doAlternate(ClpSimplex *model)
    30453001{
    30463002  // copy rowArray_[1]
    3047   ClpPrimalColumnSteepest * steep =
    3048     dynamic_cast<ClpPrimalColumnSteepest *>(model->primalColumnPivot());
     3003  ClpPrimalColumnSteepest *steep = dynamic_cast< ClpPrimalColumnSteepest * >(model->primalColumnPivot());
    30493004  if (steep) {
    3050 #if ALT_UPDATE_WEIGHTS==1
     3005#if ALT_UPDATE_WEIGHTS == 1
    30513006    if (!altVector[0]) {
    3052       altVector[1]=new CoinIndexedVector(2000);
    3053     }
    3054     CoinIndexedVector * alternateWeights2 = altVector[1];
     3007      altVector[1] = new CoinIndexedVector(2000);
     3008    }
     3009    CoinIndexedVector *alternateWeights2 = altVector[1];
    30553010#else
    3056     CoinIndexedVector * alternateWeights2 = steep->alternateWeights();
    3057 #endif
    3058     double * other2 = alternateWeights2->denseVector();
    3059     int * index2 = alternateWeights2->getIndices();
     3011    CoinIndexedVector *alternateWeights2 = steep->alternateWeights();
     3012#endif
     3013    double *other2 = alternateWeights2->denseVector();
     3014    int *index2 = alternateWeights2->getIndices();
    30603015    //rowArray packed
    3061     double * array = model->rowArray(1)->denseVector();
    3062     int * index = model->rowArray(1)->getIndices();
    3063     int n=model->rowArray(1)->getNumElements();
    3064     const int * pivotVariable = model->pivotVariable();
    3065     int n2=0;
     3016    double *array = model->rowArray(1)->denseVector();
     3017    int *index = model->rowArray(1)->getIndices();
     3018    int n = model->rowArray(1)->getNumElements();
     3019    const int *pivotVariable = model->pivotVariable();
     3020    int n2 = 0;
    30663021    alternateWeights2->clear();
    3067     if (steep->mode()!=1) {
    3068       for (int i=0;i<n;i++) {
    3069         int iRow = index[i];
    3070         int iPivot = pivotVariable[iRow];
    3071         if (steep->reference(iPivot)) {
    3072           other2[iRow] = -2.0*array[i];
    3073           index2[n2++] = iRow;
    3074         }
     3022    if (steep->mode() != 1) {
     3023      for (int i = 0; i < n; i++) {
     3024        int iRow = index[i];
     3025        int iPivot = pivotVariable[iRow];
     3026        if (steep->reference(iPivot)) {
     3027          other2[iRow] = -2.0 * array[i];
     3028          index2[n2++] = iRow;
     3029        }
    30753030      }
    30763031    } else {
    3077       for (int i=0;i<n;i++) {
    3078         int iRow = index[i];
    3079         other2[iRow] = -2.0*array[i];
    3080         index2[n2++] = iRow;
     3032      for (int i = 0; i < n; i++) {
     3033        int iRow = index[i];
     3034        other2[iRow] = -2.0 * array[i];
     3035        index2[n2++] = iRow;
    30813036      }
    30823037    }
    30833038    alternateWeights2->setNumElements(n2);
    30843039    model->factorization()->updateColumnTranspose(model->rowArray(0),
    3085                                            alternateWeights2);
     3040      alternateWeights2);
    30863041  }
    30873042}
     
    30973052  +3 max iterations (iteration done)
    30983053*/
    3099 int
    3100 ClpSimplexPrimal::pivotResult(int ifValuesPass)
     3054int ClpSimplexPrimal::pivotResult(int ifValuesPass)
    31013055{
    31023056
    3103      bool roundAgain = true;
    3104      int returnCode = -1;
    3105 
    3106      // loop round if user setting and doing refactorization
    3107      while (roundAgain) {
    3108           roundAgain = false;
    3109           returnCode = -1;
    3110           pivotRow_ = -1;
    3111           sequenceOut_ = -1;
    3112           rowArray_[1]->clear();
     3057  bool roundAgain = true;
     3058  int returnCode = -1;
     3059
     3060  // loop round if user setting and doing refactorization
     3061  while (roundAgain) {
     3062    roundAgain = false;
     3063    returnCode = -1;
     3064    pivotRow_ = -1;
     3065    sequenceOut_ = -1;
     3066    rowArray_[1]->clear();
    31133067#if 0
    31143068          {
     
    31493103#endif
    31503104
    3151           // we found a pivot column
    3152           // update the incoming column
    3153           unpackPacked(rowArray_[1]);
    3154           // save reduced cost
    3155           double saveDj = dualIn_;
    3156           factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
    3157           // Get extra rows
    3158           matrix_->extendUpdated(this, rowArray_[1], 0);
     3105    // we found a pivot column
     3106    // update the incoming column
     3107    unpackPacked(rowArray_[1]);
     3108    // save reduced cost
     3109    double saveDj = dualIn_;
     3110    factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
     3111    // Get extra rows
     3112    matrix_->extendUpdated(this, rowArray_[1], 0);
    31593113#ifdef ALT_UPDATE_WEIGHTS
    3160           cilk_spawn doAlternate(this);
    3161 #endif
    3162           // do ratio test and re-compute dj
     3114    cilk_spawn doAlternate(this);
     3115#endif
     3116    // do ratio test and re-compute dj
    31633117#ifdef CLP_USER_DRIVEN
    3164           if (solveType_ != 2 || (moreSpecialOptions_ & 512) == 0) {
    3165 #endif
    3166                primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
    3167                          ifValuesPass);
     3118    if (solveType_ != 2 || (moreSpecialOptions_ & 512) == 0) {
     3119#endif
     3120      primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
     3121        ifValuesPass);
    31683122#ifdef CLP_USER_DRIVEN
    3169                // user can tell which use it is
    3170                int status = eventHandler_->event(ClpEventHandler::pivotRow);
    3171                if (status >= 0) {
    3172                     problemStatus_ = 5;
    3173                     secondaryStatus_ = ClpEventHandler::pivotRow;
    3174                     break;
    3175                }
    3176           } else {
    3177                int status = eventHandler_->event(ClpEventHandler::pivotRow);
    3178                if (status >= 0) {
    3179                     problemStatus_ = 5;
    3180                     secondaryStatus_ = ClpEventHandler::pivotRow;
    3181                     break;
    3182                }
    3183           }
    3184 #endif
    3185           if (ifValuesPass) {
    3186                saveDj = dualIn_;
    3187                //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
    3188                if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
    3189                     if(fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {
    3190                          // try other way
    3191                          directionIn_ = -directionIn_;
    3192                          primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
    3193                                    0);
    3194                     }
    3195                     if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
    3196                          if (solveType_ == 1) {
    3197                               // reject it
    3198                               char x = isColumn(sequenceIn_) ? 'C' : 'R';
    3199                               handler_->message(CLP_SIMPLEX_FLAG, messages_)
    3200                                         << x << sequenceWithin(sequenceIn_)
    3201                                         << CoinMessageEol;
    3202                               setFlagged(sequenceIn_);
    3203                               progress_.clearBadTimes();
    3204                               lastBadIteration_ = numberIterations_; // say be more cautious
    3205                               clearAll();
    3206                               pivotRow_ = -1;
    3207                          }
    3208                          returnCode = -5;
    3209                          break;
    3210                     }
    3211                }
    3212           }
     3123      // user can tell which use it is
     3124      int status = eventHandler_->event(ClpEventHandler::pivotRow);
     3125      if (status >= 0) {
     3126        problemStatus_ = 5;
     3127        secondaryStatus_ = ClpEventHandler::pivotRow;
     3128        break;
     3129      }
     3130    } else {
     3131      int status = eventHandler_->event(ClpEventHandler::pivotRow);
     3132      if (status >= 0) {
     3133        problemStatus_ = 5;
     3134        secondaryStatus_ = ClpEventHandler::pivotRow;
     3135        break;
     3136      }
     3137    }
     3138#endif
     3139    if (ifValuesPass) {
     3140      saveDj = dualIn_;
     3141      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
     3142      if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
     3143        if (fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {