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

formatting

File:
1 edited

Legend:

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

    r2326 r2385  
    33// Corporation and others.  All Rights Reserved.
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    5 
    65
    76/* Notes on implementation of dual simplex algorithm.
     
    122121//#define FORCE_FOLLOW
    123122#ifdef FORCE_FOLLOW
    124 static FILE * fpFollow = NULL;
    125 static char * forceFile = "old.log";
     123static FILE *fpFollow = NULL;
     124static char *forceFile = "old.log";
    126125static int force_in = -1;
    127126static int force_out = -1;
     
    130129//#define VUB
    131130#ifdef VUB
    132 extern int * vub;
    133 extern int * toVub;
    134 extern int * nextDescendent;
     131extern int *vub;
     132extern int *toVub;
     133extern int *nextDescendent;
    135134#endif
    136135#ifdef NDEBUG
     
    229228#define CLEAN_FIXED 0
    230229// Startup part of dual (may be extended to other algorithms)
    231 int
    232 ClpSimplexDual::startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions)
     230int ClpSimplexDual::startupSolve(int ifValuesPass, double *saveDuals, int startFinishOptions)
    233231{
    234      // If values pass then save given duals round check solution
    235      // sanity check
    236      // initialize - no values pass and algorithm_ is -1
    237      // put in standard form (and make row copy)
    238      // create modifiable copies of model rim and do optional scaling
    239      // If problem looks okay
    240      // Do initial factorization
    241      // If user asked for perturbation - do it
    242      numberFake_ = 0; // Number of variables at fake bounds
    243      numberChanged_ = 0; // Number of variables with changed costs
    244      if (!startup(0/* ? fix valuesPass */, startFinishOptions)) {
    245           int usePrimal = 0;
    246           // looks okay
    247           // Superbasic variables not allowed
    248           // If values pass then scale pi
    249           if (ifValuesPass) {
    250                if (problemStatus_ && perturbation_ < 100)
    251                     usePrimal = perturb();
    252                int i;
    253                if (scalingFlag_ > 0) {
    254                     for (i = 0; i < numberRows_; i++) {
    255                          dual_[i] = saveDuals[i] * inverseRowScale_[i];
    256                     }
    257                } else {
    258                     CoinMemcpyN(saveDuals, numberRows_, dual_);
    259                }
    260                // now create my duals
    261                for (i = 0; i < numberRows_; i++) {
    262                     // slack
    263                     double value = dual_[i];
    264                     value += rowObjectiveWork_[i];
    265                     saveDuals[i+numberColumns_] = value;
    266                }
    267                CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
    268                transposeTimes(-1.0, dual_, saveDuals);
    269                // make reduced costs okay
    270                for (i = 0; i < numberColumns_; i++) {
    271                     if (getStatus(i) == atLowerBound) {
    272                          if (saveDuals[i] < 0.0) {
    273                               //if (saveDuals[i]<-1.0e-3)
    274                               //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
    275                               saveDuals[i] = 0.0;
    276                          }
    277                     } else if (getStatus(i) == atUpperBound) {
    278                          if (saveDuals[i] > 0.0) {
    279                               //if (saveDuals[i]>1.0e-3)
    280                               //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
    281                               saveDuals[i] = 0.0;
    282                          }
    283                     }
    284                }
    285                CoinMemcpyN(saveDuals, (numberColumns_ + numberRows_), dj_);
    286                // set up possible ones
    287                for (i = 0; i < numberRows_ + numberColumns_; i++)
    288                     clearPivoted(i);
    289                int iRow;
    290                for (iRow = 0; iRow < numberRows_; iRow++) {
    291                     int iPivot = pivotVariable_[iRow];
    292                     if (fabs(saveDuals[iPivot]) > dualTolerance_) {
    293                          if (getStatus(iPivot) != isFree)
    294                               setPivoted(iPivot);
    295                     }
    296                }
    297           } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
    298                // set up possible ones
    299                for (int i = 0; i < numberRows_ + numberColumns_; i++)
    300                     clearPivoted(i);
    301                int iRow;
    302                for (iRow = 0; iRow < numberRows_; iRow++) {
    303                     int iPivot = pivotVariable_[iRow];
    304                     if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
    305                          setPivoted(iPivot);
    306                     }
    307                }
    308           }
    309 
    310           double objectiveChange;
    311           assert (!numberFake_);
    312           assert (numberChanged_ == 0);
    313           if (!numberFake_) // if nonzero then adjust
    314                changeBounds(1, NULL, objectiveChange);
    315 
    316           if (!ifValuesPass) {
    317                // Check optimal
    318                if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
    319                     problemStatus_ = 0;
    320           }
    321           if (problemStatus_ < 0 && perturbation_ < 100) {
    322                bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
    323                if (!inCbcOrOther)
    324                     usePrimal = perturb();
    325                // Can't get here if values pass
    326                gutsOfSolution(NULL, NULL);
     232  // If values pass then save given duals round check solution
     233  // sanity check
     234  // initialize - no values pass and algorithm_ is -1
     235  // put in standard form (and make row copy)
     236  // create modifiable copies of model rim and do optional scaling
     237  // If problem looks okay
     238  // Do initial factorization
     239  // If user asked for perturbation - do it
     240  numberFake_ = 0; // Number of variables at fake bounds
     241  numberChanged_ = 0; // Number of variables with changed costs
     242  if (!startup(0 /* ? fix valuesPass */, startFinishOptions)) {
     243    int usePrimal = 0;
     244    // looks okay
     245    // Superbasic variables not allowed
     246    // If values pass then scale pi
     247    if (ifValuesPass) {
     248      if (problemStatus_ && perturbation_ < 100)
     249        usePrimal = perturb();
     250      int i;
     251      if (scalingFlag_ > 0) {
     252        for (i = 0; i < numberRows_; i++) {
     253          dual_[i] = saveDuals[i] * inverseRowScale_[i];
     254        }
     255      } else {
     256        CoinMemcpyN(saveDuals, numberRows_, dual_);
     257      }
     258      // now create my duals
     259      for (i = 0; i < numberRows_; i++) {
     260        // slack
     261        double value = dual_[i];
     262        value += rowObjectiveWork_[i];
     263        saveDuals[i + numberColumns_] = value;
     264      }
     265      CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
     266      transposeTimes(-1.0, dual_, saveDuals);
     267      // make reduced costs okay
     268      for (i = 0; i < numberColumns_; i++) {
     269        if (getStatus(i) == atLowerBound) {
     270          if (saveDuals[i] < 0.0) {
     271            //if (saveDuals[i]<-1.0e-3)
     272            //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
     273            saveDuals[i] = 0.0;
     274          }
     275        } else if (getStatus(i) == atUpperBound) {
     276          if (saveDuals[i] > 0.0) {
     277            //if (saveDuals[i]>1.0e-3)
     278            //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
     279            saveDuals[i] = 0.0;
     280          }
     281        }
     282      }
     283      CoinMemcpyN(saveDuals, (numberColumns_ + numberRows_), dj_);
     284      // set up possible ones
     285      for (i = 0; i < numberRows_ + numberColumns_; i++)
     286        clearPivoted(i);
     287      int iRow;
     288      for (iRow = 0; iRow < numberRows_; iRow++) {
     289        int iPivot = pivotVariable_[iRow];
     290        if (fabs(saveDuals[iPivot]) > dualTolerance_) {
     291          if (getStatus(iPivot) != isFree)
     292            setPivoted(iPivot);
     293        }
     294      }
     295    } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
     296      // set up possible ones
     297      for (int i = 0; i < numberRows_ + numberColumns_; i++)
     298        clearPivoted(i);
     299      int iRow;
     300      for (iRow = 0; iRow < numberRows_; iRow++) {
     301        int iPivot = pivotVariable_[iRow];
     302        if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
     303          setPivoted(iPivot);
     304        }
     305      }
     306    }
     307
     308    double objectiveChange;
     309    assert(!numberFake_);
     310    assert(numberChanged_ == 0);
     311    if (!numberFake_) // if nonzero then adjust
     312      changeBounds(1, NULL, objectiveChange);
     313
     314    if (!ifValuesPass) {
     315      // Check optimal
     316      if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
     317        problemStatus_ = 0;
     318    }
     319    if (problemStatus_ < 0 && perturbation_ < 100) {
     320      bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
     321      if (!inCbcOrOther)
     322        usePrimal = perturb();
     323      // Can't get here if values pass
     324      gutsOfSolution(NULL, NULL);
    327325#ifdef CLP_INVESTIGATE
    328                if (numberDualInfeasibilities_)
    329                     printf("ZZZ %d primal %d dual - sumdinf %g\n",
    330                            numberPrimalInfeasibilities_,
    331                            numberDualInfeasibilities_, sumDualInfeasibilities_);
    332 #endif
    333                if (handler_->logLevel() > 2) {
    334                     handler_->message(CLP_SIMPLEX_STATUS, messages_)
    335                               << numberIterations_ << objectiveValue();
    336                     handler_->printing(sumPrimalInfeasibilities_ > 0.0)
    337                               << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
    338                     handler_->printing(sumDualInfeasibilities_ > 0.0)
    339                               << sumDualInfeasibilities_ << numberDualInfeasibilities_;
    340                     handler_->printing(numberDualInfeasibilitiesWithoutFree_
    341                                        < numberDualInfeasibilities_)
    342                               << numberDualInfeasibilitiesWithoutFree_;
    343                     handler_->message() << CoinMessageEol;
    344                }
    345                if (inCbcOrOther) {
    346                     if (numberPrimalInfeasibilities_) {
    347                          usePrimal = perturb();
    348                          if (perturbation_ >= 101) {
    349                               computeDuals(NULL);
    350                               //gutsOfSolution(NULL,NULL);
    351                               checkDualSolution(); // recompute objective
    352                          }
    353                     } else if (numberDualInfeasibilities_) {
    354                          problemStatus_ = 10;
    355                          if ((moreSpecialOptions_ & 32) != 0 && false)
    356                               problemStatus_ = 0; // say optimal!!
    357 #if COIN_DEVELOP>2
    358 
    359                          printf("returning at %d\n", __LINE__);
    360 #endif
    361                          return 1; // to primal
    362                     }
    363                }
    364           } else if (!ifValuesPass) {
    365                gutsOfSolution(NULL, NULL);
    366                // double check
    367                if (numberDualInfeasibilities_ || numberPrimalInfeasibilities_)
    368                     problemStatus_ = -1;
    369           }
    370           if (usePrimal) {
    371                problemStatus_ = 10;
    372 #if COIN_DEVELOP>2
    373                printf("returning to use primal (no obj) at %d\n", __LINE__);
    374 #endif
    375           }
    376           return usePrimal;
    377      } else {
    378           return 1;
    379      }
     326      if (numberDualInfeasibilities_)
     327        printf("ZZZ %d primal %d dual - sumdinf %g\n",
     328          numberPrimalInfeasibilities_,
     329          numberDualInfeasibilities_, sumDualInfeasibilities_);
     330#endif
     331      if (handler_->logLevel() > 2) {
     332        handler_->message(CLP_SIMPLEX_STATUS, messages_)
     333          << numberIterations_ << objectiveValue();
     334        handler_->printing(sumPrimalInfeasibilities_ > 0.0)
     335          << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
     336        handler_->printing(sumDualInfeasibilities_ > 0.0)
     337          << sumDualInfeasibilities_ << numberDualInfeasibilities_;
     338        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     339          < numberDualInfeasibilities_)
     340          << numberDualInfeasibilitiesWithoutFree_;
     341        handler_->message() << CoinMessageEol;
     342      }
     343      if (inCbcOrOther) {
     344        if (numberPrimalInfeasibilities_) {
     345          usePrimal = perturb();
     346          if (perturbation_ >= 101) {
     347            computeDuals(NULL);
     348            //gutsOfSolution(NULL,NULL);
     349            checkDualSolution(); // recompute objective
     350          }
     351        } else if (numberDualInfeasibilities_) {
     352          problemStatus_ = 10;
     353          if ((moreSpecialOptions_ & 32) != 0 && false)
     354            problemStatus_ = 0; // say optimal!!
     355#if COIN_DEVELOP > 2
     356
     357          printf("returning at %d\n", __LINE__);
     358#endif
     359          return 1; // to primal
     360        }
     361      }
     362    } else if (!ifValuesPass) {
     363      gutsOfSolution(NULL, NULL);
     364      // double check
     365      if (numberDualInfeasibilities_ || numberPrimalInfeasibilities_)
     366        problemStatus_ = -1;
     367    }
     368    if (usePrimal) {
     369      problemStatus_ = 10;
     370#if COIN_DEVELOP > 2
     371      printf("returning to use primal (no obj) at %d\n", __LINE__);
     372#endif
     373    }
     374    return usePrimal;
     375  } else {
     376    return 1;
     377  }
    380378}
    381 void
    382 ClpSimplexDual::finishSolve(int startFinishOptions)
     379void ClpSimplexDual::finishSolve(int startFinishOptions)
    383380{
    384      assert (problemStatus_ || !sumPrimalInfeasibilities_);
    385 
    386      // clean up
    387      finish(startFinishOptions);
     381  assert(problemStatus_ || !sumPrimalInfeasibilities_);
     382
     383  // clean up
     384  finish(startFinishOptions);
    388385}
    389386//#define CLP_REPORT_PROGRESS
     
    393390#endif
    394391#ifdef CLP_INVESTIGATE_SERIAL
    395 static int z_reason[7] = {0, 0, 0, 0, 0, 0, 0};
     392static int z_reason[7] = { 0, 0, 0, 0, 0, 0, 0 };
    396393static int z_thinks = -1;
    397394#endif
    398 void
    399 ClpSimplexDual::gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
    400                            ClpDataSave & data)
     395void ClpSimplexDual::gutsOfDual(int ifValuesPass, double *&saveDuals, int initialStatus,
     396  ClpDataSave &data)
    401397{
    402398#ifdef CLP_INVESTIGATE_SERIAL
    403      z_reason[0]++;
    404      z_thinks = -1;
    405      int nPivots = 9999;
    406 #endif
    407      double largestPrimalError = 0.0;
    408      double largestDualError = 0.0;
    409      double smallestPrimalInfeasibility=COIN_DBL_MAX;
    410      int numberRayTries=0;
    411      // Start can skip some things in transposeTimes
    412      specialOptions_ |= 131072;
    413      int lastCleaned = 0; // last time objective or bounds cleaned up
    414 
    415      // This says whether to restore things etc
    416      // startup will have factorized so can skip
    417      int factorType = 0;
    418      // Start check for cycles
    419      progress_.startCheck();
    420      // Say change made on first iteration
    421      changeMade_ = 1;
    422      // Say last objective infinite
    423      //lastObjectiveValue_=-COIN_DBL_MAX;
    424      progressFlag_ = 0;
    425      /*
     399  z_reason[0]++;
     400  z_thinks = -1;
     401  int nPivots = 9999;
     402#endif
     403  double largestPrimalError = 0.0;
     404  double largestDualError = 0.0;
     405  double smallestPrimalInfeasibility = COIN_DBL_MAX;
     406  int numberRayTries = 0;
     407  // Start can skip some things in transposeTimes
     408  specialOptions_ |= 131072;
     409  int lastCleaned = 0; // last time objective or bounds cleaned up
     410
     411  // This says whether to restore things etc
     412  // startup will have factorized so can skip
     413  int factorType = 0;
     414  // Start check for cycles
     415  progress_.startCheck();
     416  // Say change made on first iteration
     417  changeMade_ = 1;
     418  // Say last objective infinite
     419  //lastObjectiveValue_=-COIN_DBL_MAX;
     420  progressFlag_ = 0;
     421  /*
    426422       Status of problem:
    427423       0 - optimal
     
    433429       -4 - looks infeasible
    434430     */
    435      while (problemStatus_ < 0) {
    436           int iRow, iColumn;
    437           // clear
    438           for (iRow = 0; iRow < 4; iRow++) {
    439                rowArray_[iRow]->clear();
    440           }
    441 
    442           for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
    443                columnArray_[iColumn]->clear();
    444           }
    445 
    446           // give matrix (and model costs and bounds a chance to be
    447           // refreshed (normally null)
    448           matrix_->refresh(this);
    449           // If getting nowhere - why not give it a kick
    450           // does not seem to work too well - do some more work
    451           if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (moreSpecialOptions_&1048576)==0
    452                     && initialStatus != 10) {
    453                perturb();
    454                // Can't get here if values pass
    455                gutsOfSolution(NULL, NULL);
    456                if (handler_->logLevel() > 2) {
    457                     handler_->message(CLP_SIMPLEX_STATUS, messages_)
    458                               << numberIterations_ << objectiveValue();
    459                     handler_->printing(sumPrimalInfeasibilities_ > 0.0)
    460                               << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
    461                     handler_->printing(sumDualInfeasibilities_ > 0.0)
    462                               << sumDualInfeasibilities_ << numberDualInfeasibilities_;
    463                     handler_->printing(numberDualInfeasibilitiesWithoutFree_
    464                                        < numberDualInfeasibilities_)
    465                               << numberDualInfeasibilitiesWithoutFree_;
    466                     handler_->message() << CoinMessageEol;
    467                }
    468           }
    469           // see if in Cbc etc
    470           bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
     431  while (problemStatus_ < 0) {
     432    int iRow, iColumn;
     433    // clear
     434    for (iRow = 0; iRow < 4; iRow++) {
     435      rowArray_[iRow]->clear();
     436    }
     437
     438    for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
     439      columnArray_[iColumn]->clear();
     440    }
     441
     442    // give matrix (and model costs and bounds a chance to be
     443    // refreshed (normally null)
     444    matrix_->refresh(this);
     445    // If getting nowhere - why not give it a kick
     446    // does not seem to work too well - do some more work
     447    if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (moreSpecialOptions_ & 1048576) == 0
     448      && initialStatus != 10) {
     449      perturb();
     450      // Can't get here if values pass
     451      gutsOfSolution(NULL, NULL);
     452      if (handler_->logLevel() > 2) {
     453        handler_->message(CLP_SIMPLEX_STATUS, messages_)
     454          << numberIterations_ << objectiveValue();
     455        handler_->printing(sumPrimalInfeasibilities_ > 0.0)
     456          << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
     457        handler_->printing(sumDualInfeasibilities_ > 0.0)
     458          << sumDualInfeasibilities_ << numberDualInfeasibilities_;
     459        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     460          < numberDualInfeasibilities_)
     461          << numberDualInfeasibilitiesWithoutFree_;
     462        handler_->message() << CoinMessageEol;
     463      }
     464    }
     465    // see if in Cbc etc
     466    bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
    471467#if 0
    472468          bool gotoPrimal = false;
     
    486482          }
    487483#endif
    488           bool disaster = false;
    489           if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
    490                disasterArea_->saveInfo();
    491                disaster = true;
    492           }
    493           // may factorize, checks if problem finished
    494           statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
    495                                 ifValuesPass);
    496           smallestPrimalInfeasibility=CoinMin(smallestPrimalInfeasibility,
    497                                               sumPrimalInfeasibilities_);
    498           if (sumPrimalInfeasibilities_>1.0e5 &&
    499               sumPrimalInfeasibilities_>1.0e5*smallestPrimalInfeasibility &&
    500               (moreSpecialOptions_&256)==0 &&
    501               ((progress_.lastObjective(0)<-1.0e10 &&
    502 -               progress_.lastObjective(1)>-1.0e5)||sumPrimalInfeasibilities_>1.0e10*smallestPrimalInfeasibility)&&problemStatus_<0) {
    503             // problems - try primal
    504             problemStatus_=10;
    505             // mark as large infeasibility cost wanted
    506             sumPrimalInfeasibilities_ = -123456789.0;
    507             //for (int i=0;i<numberRows_+numberColumns_;i++) {
    508             //if (fabs(cost_[i]*solution_[i])>1.0e4)
    509             //  printf("col %d cost %g sol %g bounds %g %g\n",
    510             //         i,cost_[i],solution_[i],lower_[i],upper_[i]);
    511             //}
    512           } else if ((specialOptions_&(32|2097152))!=0&&
    513                      problemStatus_==1&&!ray_&&
    514                      !numberRayTries && numberIterations_) {
    515             numberRayTries=1;
    516             problemStatus_=-1;
    517           }
    518           largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
    519           largestDualError = CoinMax(largestDualError, largestDualError_);
    520           if (disaster)
    521                problemStatus_ = 3;
    522           // If values pass then do easy ones on first time
    523           if (ifValuesPass &&
    524                     progress_.lastIterationNumber(0) < 0 && saveDuals) {
    525                doEasyOnesInValuesPass(saveDuals);
    526           }
    527 
    528           // Say good factorization
    529           factorType = 1;
    530           if (data.sparseThreshold_) {
    531                // use default at present
    532                factorization_->sparseThreshold(0);
    533                factorization_->goSparse();
    534           }
    535 
    536           // exit if victory declared
    537           if (problemStatus_ >= 0)
    538                break;
    539 
    540           // test for maximum iterations
    541           if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
    542                problemStatus_ = 3;
    543                break;
    544           }
    545           if (ifValuesPass && !saveDuals) {
    546                // end of values pass
    547                ifValuesPass = 0;
    548                int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
    549                if (status >= 0) {
    550                     problemStatus_ = 5;
    551                     secondaryStatus_ = ClpEventHandler::endOfValuesPass;
    552                     break;
    553                }
    554           }
    555           // Check event
    556           {
    557                int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
    558                if (status >= 0) {
    559                     problemStatus_ = 5;
    560                     secondaryStatus_ = ClpEventHandler::endOfFactorization;
    561                     break;
    562                }
    563           }
    564           // If looks odd try other way
    565           if ((moreSpecialOptions_&256)==0 &&
    566               fabs(objectiveValue_)>1.0e20&&sumDualInfeasibilities_>1.0
    567               &&problemStatus_<0) {
    568             problemStatus_=10;
    569             break;
    570           }
    571           // Do iterations
    572           int returnCode = whileIterating(saveDuals, ifValuesPass);
    573           if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
    574               fabs(objectiveValue_) > 1.0e10 )
    575             problemStatus_ = 10; // infeasible - but has looked feasible
     484    bool disaster = false;
     485    if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
     486      disasterArea_->saveInfo();
     487      disaster = true;
     488    }
     489    // may factorize, checks if problem finished
     490    statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
     491      ifValuesPass);
     492    smallestPrimalInfeasibility = CoinMin(smallestPrimalInfeasibility,
     493      sumPrimalInfeasibilities_);
     494    if (sumPrimalInfeasibilities_ > 1.0e5 && sumPrimalInfeasibilities_ > 1.0e5 * smallestPrimalInfeasibility && (moreSpecialOptions_ & 256) == 0 && ((progress_.lastObjective(0) < -1.0e10 && -progress_.lastObjective(1) > -1.0e5) || sumPrimalInfeasibilities_ > 1.0e10 * smallestPrimalInfeasibility) && problemStatus_ < 0) {
     495      // problems - try primal
     496      problemStatus_ = 10;
     497      // mark as large infeasibility cost wanted
     498      sumPrimalInfeasibilities_ = -123456789.0;
     499      //for (int i=0;i<numberRows_+numberColumns_;i++) {
     500      //if (fabs(cost_[i]*solution_[i])>1.0e4)
     501      //        printf("col %d cost %g sol %g bounds %g %g\n",
     502      //               i,cost_[i],solution_[i],lower_[i],upper_[i]);
     503      //}
     504    } else if ((specialOptions_ & (32 | 2097152)) != 0 && problemStatus_ == 1 && !ray_ && !numberRayTries && numberIterations_) {
     505      numberRayTries = 1;
     506      problemStatus_ = -1;
     507    }
     508    largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
     509    largestDualError = CoinMax(largestDualError, largestDualError_);
     510    if (disaster)
     511      problemStatus_ = 3;
     512    // If values pass then do easy ones on first time
     513    if (ifValuesPass && progress_.lastIterationNumber(0) < 0 && saveDuals) {
     514      doEasyOnesInValuesPass(saveDuals);
     515    }
     516
     517    // Say good factorization
     518    factorType = 1;
     519    if (data.sparseThreshold_) {
     520      // use default at present
     521      factorization_->sparseThreshold(0);
     522      factorization_->goSparse();
     523    }
     524
     525    // exit if victory declared
     526    if (problemStatus_ >= 0)
     527      break;
     528
     529    // test for maximum iterations
     530    if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
     531      problemStatus_ = 3;
     532      break;
     533    }
     534    if (ifValuesPass && !saveDuals) {
     535      // end of values pass
     536      ifValuesPass = 0;
     537      int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
     538      if (status >= 0) {
     539        problemStatus_ = 5;
     540        secondaryStatus_ = ClpEventHandler::endOfValuesPass;
     541        break;
     542      }
     543    }
     544    // Check event
     545    {
     546      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
     547      if (status >= 0) {
     548        problemStatus_ = 5;
     549        secondaryStatus_ = ClpEventHandler::endOfFactorization;
     550        break;
     551      }
     552    }
     553    // If looks odd try other way
     554    if ((moreSpecialOptions_ & 256) == 0 && fabs(objectiveValue_) > 1.0e20 && sumDualInfeasibilities_ > 1.0
     555      && problemStatus_ < 0) {
     556      problemStatus_ = 10;
     557      break;
     558    }
     559    // Do iterations
     560    int returnCode = whileIterating(saveDuals, ifValuesPass);
     561    if (problemStatus_ == 1 && (progressFlag_ & 8) != 0 && fabs(objectiveValue_) > 1.0e10)
     562      problemStatus_ = 10; // infeasible - but has looked feasible
    576563#ifdef CLP_INVESTIGATE_SERIAL
    577           nPivots = factorization_->pivots();
    578 #endif
    579           if (!problemStatus_ && factorization_->pivots())
    580                computeDuals(NULL); // need to compute duals
    581           if (returnCode == -2)
    582                factorType = 3;
    583      }
     564    nPivots = factorization_->pivots();
     565#endif
     566    if (!problemStatus_ && factorization_->pivots())
     567      computeDuals(NULL); // need to compute duals
     568    if (returnCode == -2)
     569      factorType = 3;
     570  }
    584571#ifdef CLP_INVESTIGATE_SERIAL
    585      // NOTE - can fail if parallel
    586      if (z_thinks != -1) {
    587           assert (z_thinks < 4);
    588           if ((!factorization_->pivots() && nPivots < 20) && z_thinks >= 0 && z_thinks < 2)
    589                z_thinks += 4;
    590           z_reason[1+z_thinks]++;
    591      }
    592      if ((z_reason[0] % 1000) == 0) {
    593           printf("Reason");
    594           for (int i = 0; i < 7; i++)
    595                printf(" %d", z_reason[i]);
    596           printf("\n");
    597      }
    598 #endif
    599      // Stop can skip some things in transposeTimes
    600      specialOptions_ &= ~131072;
    601      largestPrimalError_ = largestPrimalError;
    602      largestDualError_ = largestDualError;
     572  // NOTE - can fail if parallel
     573  if (z_thinks != -1) {
     574    assert(z_thinks < 4);
     575    if ((!factorization_->pivots() && nPivots < 20) && z_thinks >= 0 && z_thinks < 2)
     576      z_thinks += 4;
     577    z_reason[1 + z_thinks]++;
     578  }
     579  if ((z_reason[0] % 1000) == 0) {
     580    printf("Reason");
     581    for (int i = 0; i < 7; i++)
     582      printf(" %d", z_reason[i]);
     583    printf("\n");
     584  }
     585#endif
     586  // Stop can skip some things in transposeTimes
     587  specialOptions_ &= ~131072;
     588  largestPrimalError_ = largestPrimalError;
     589  largestDualError_ = largestDualError;
    603590}
    604 int
    605 ClpSimplexDual::dual(int ifValuesPass, int startFinishOptions)
     591int ClpSimplexDual::dual(int ifValuesPass, int startFinishOptions)
    606592{
    607593  //handler_->setLogLevel(63);
    608594  //yprintf("STARTing dual %d rows\n",numberRows_);
    609      bestObjectiveValue_ = -COIN_DBL_MAX;
    610      algorithm_ = -1;
    611      moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
    612      delete [] ray_;
    613      ray_ = NULL;
    614      // save data
    615      ClpDataSave data = saveData();
    616      double * saveDuals = NULL;
    617      int saveDont = dontFactorizePivots_;
    618      if ((specialOptions_ & 2048) == 0)
    619           dontFactorizePivots_ = 0;
    620      else if(!dontFactorizePivots_)
    621           dontFactorizePivots_ = 20;
    622      if (ifValuesPass) {
    623           saveDuals = new double [numberRows_+numberColumns_];
    624           CoinMemcpyN(dual_, numberRows_, saveDuals);
    625      }
    626      if (alphaAccuracy_ != -1.0)
    627           alphaAccuracy_ = 1.0;
    628      minimumPrimalTolerance_=primalTolerance();
    629      int returnCode = startupSolve(ifValuesPass, saveDuals, startFinishOptions);
    630      // Save so can see if doing after primal
    631      int initialStatus = problemStatus_;
    632      if (!returnCode && !numberDualInfeasibilities_ &&
    633                !numberPrimalInfeasibilities_ && perturbation_ < 101) {
    634           returnCode = 1; // to skip gutsOfDual
    635           problemStatus_ = 0;
    636      } else if (maximumIterations()==0) {
    637           returnCode = 1; // to skip gutsOfDual
    638           problemStatus_ = 3;
    639      }
    640 
    641      if (!returnCode)
    642           gutsOfDual(ifValuesPass, saveDuals, initialStatus, data);
    643      if (!problemStatus_) {
    644           // see if cutoff reached
    645           double limit = 0.0;
    646           getDblParam(ClpDualObjectiveLimit, limit);
    647           if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
    648                     limit + 1.0e-7 + 1.0e-8 * fabs(limit)) {
    649                // actually infeasible on objective
    650                problemStatus_ = 1;
    651                secondaryStatus_ = 1;
    652           }
    653      }
    654      // If infeasible but primal errors - try primal
    655      if (problemStatus_==1 && numberPrimalInfeasibilities_) {
    656        bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
    657        double factor = (!inCbcOrOther) ? 1.0 : 0.3;
    658        double averageInfeasibility = sumPrimalInfeasibilities_/
    659          static_cast<double>(numberPrimalInfeasibilities_);
    660        if (averageInfeasibility<factor*largestPrimalError_)
    661          problemStatus_= 10;
    662      }
    663        
    664      if (problemStatus_ == 10)
    665           startFinishOptions |= 1;
    666      finishSolve(startFinishOptions);
    667      delete [] saveDuals;
    668 
    669      // Restore any saved stuff
    670      restoreData(data);
    671      dontFactorizePivots_ = saveDont;
    672      if (problemStatus_ == 3)
    673           objectiveValue_ = CoinMax(bestObjectiveValue_, objectiveValue_ - bestPossibleImprovement_);
    674      return problemStatus_;
     595  bestObjectiveValue_ = -COIN_DBL_MAX;
     596  algorithm_ = -1;
     597  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
     598  delete[] ray_;
     599  ray_ = NULL;
     600  // save data
     601  ClpDataSave data = saveData();
     602  double *saveDuals = NULL;
     603  int saveDont = dontFactorizePivots_;
     604  if ((specialOptions_ & 2048) == 0)
     605    dontFactorizePivots_ = 0;
     606  else if (!dontFactorizePivots_)
     607    dontFactorizePivots_ = 20;
     608  if (ifValuesPass) {
     609    saveDuals = new double[numberRows_ + numberColumns_];
     610    CoinMemcpyN(dual_, numberRows_, saveDuals);
     611  }
     612  if (alphaAccuracy_ != -1.0)
     613    alphaAccuracy_ = 1.0;
     614  minimumPrimalTolerance_ = primalTolerance();
     615  int returnCode = startupSolve(ifValuesPass, saveDuals, startFinishOptions);
     616  // Save so can see if doing after primal
     617  int initialStatus = problemStatus_;
     618  if (!returnCode && !numberDualInfeasibilities_ && !numberPrimalInfeasibilities_ && perturbation_ < 101) {
     619    returnCode = 1; // to skip gutsOfDual
     620    problemStatus_ = 0;
     621  } else if (maximumIterations() == 0) {
     622    returnCode = 1; // to skip gutsOfDual
     623    problemStatus_ = 3;
     624  }
     625
     626  if (!returnCode)
     627    gutsOfDual(ifValuesPass, saveDuals, initialStatus, data);
     628  if (!problemStatus_) {
     629    // see if cutoff reached
     630    double limit = 0.0;
     631    getDblParam(ClpDualObjectiveLimit, limit);
     632    if (fabs(limit) < 1.0e30 && objectiveValue() * optimizationDirection_ > limit + 1.0e-7 + 1.0e-8 * fabs(limit)) {
     633      // actually infeasible on objective
     634      problemStatus_ = 1;
     635      secondaryStatus_ = 1;
     636    }
     637  }
     638  // If infeasible but primal errors - try primal
     639  if (problemStatus_ == 1 && numberPrimalInfeasibilities_) {
     640    bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
     641    double factor = (!inCbcOrOther) ? 1.0 : 0.3;
     642    double averageInfeasibility = sumPrimalInfeasibilities_ / static_cast< double >(numberPrimalInfeasibilities_);
     643    if (averageInfeasibility < factor * largestPrimalError_)
     644      problemStatus_ = 10;
     645  }
     646
     647  if (problemStatus_ == 10)
     648    startFinishOptions |= 1;
     649  finishSolve(startFinishOptions);
     650  delete[] saveDuals;
     651
     652  // Restore any saved stuff
     653  restoreData(data);
     654  dontFactorizePivots_ = saveDont;
     655  if (problemStatus_ == 3)
     656    objectiveValue_ = CoinMax(bestObjectiveValue_, objectiveValue_ - bestPossibleImprovement_);
     657  return problemStatus_;
    675658}
    676659// old way
     
    917900   +3 max iterations
    918901 */
    919 int
    920 ClpSimplexDual::whileIterating(double * & givenDuals, int ifValuesPass)
     902int ClpSimplexDual::whileIterating(double *&givenDuals, int ifValuesPass)
    921903{
    922904#ifdef CLP_INVESTIGATE_SERIAL
    923      z_thinks = -1;
     905  z_thinks = -1;
    924906#endif
    925907#ifdef CLP_DEBUG
    926      int debugIteration = -1;
    927 #endif
    928      {
    929           int i;
    930           for (i = 0; i < 4; i++) {
    931                rowArray_[i]->clear();
    932           }
    933           for (i = 0; i < SHORT_REGION; i++) {
    934                columnArray_[i]->clear();
    935           }
    936      }
     908  int debugIteration = -1;
     909#endif
     910  {
     911    int i;
     912    for (i = 0; i < 4; i++) {
     913      rowArray_[i]->clear();
     914    }
     915    for (i = 0; i < SHORT_REGION; i++) {
     916      columnArray_[i]->clear();
     917    }
     918  }
    937919#ifdef CLP_REPORT_PROGRESS
    938      double * savePSol = new double [numberRows_+numberColumns_];
    939      double * saveDj = new double [numberRows_+numberColumns_];
    940      double * saveCost = new double [numberRows_+numberColumns_];
    941      unsigned char * saveStat = new unsigned char [numberRows_+numberColumns_];
    942 #endif
    943      // if can't trust much and long way from optimal then relax
    944      if (largestPrimalError_ > 10.0)
    945           factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
    946      else
    947           factorization_->relaxAccuracyCheck(1.0);
    948      // status stays at -1 while iterating, >=0 finished, -2 to invert
    949      // status -3 to go to top without an invert
    950      int returnCode = -1;
    951      double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
     920  double *savePSol = new double[numberRows_ + numberColumns_];
     921  double *saveDj = new double[numberRows_ + numberColumns_];
     922  double *saveCost = new double[numberRows_ + numberColumns_];
     923  unsigned char *saveStat = new unsigned char[numberRows_ + numberColumns_];
     924#endif
     925  // if can't trust much and long way from optimal then relax
     926  if (largestPrimalError_ > 10.0)
     927    factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
     928  else
     929    factorization_->relaxAccuracyCheck(1.0);
     930  // status stays at -1 while iterating, >=0 finished, -2 to invert
     931  // status -3 to go to top without an invert
     932  int returnCode = -1;
     933  double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
    952934
    953935#if 0
     
    957939#endif
    958940
    959      // Get dubious weights
    960      CoinBigIndex * dubiousWeights = NULL;
     941  // Get dubious weights
     942  CoinBigIndex *dubiousWeights = NULL;
    961943#ifdef DUBIOUS_WEIGHTS
    962      factorization_->getWeights(rowArray_[0]->getIndices());
    963      dubiousWeights = matrix_->dubiousWeights(this, rowArray_[0]->getIndices());
    964 #endif
    965      // If values pass then get list of candidates
    966      int * candidateList = NULL;
    967      int numberCandidates = 0;
     944  factorization_->getWeights(rowArray_[0]->getIndices());
     945  dubiousWeights = matrix_->dubiousWeights(this, rowArray_[0]->getIndices());
     946#endif
     947  // If values pass then get list of candidates
     948  int *candidateList = NULL;
     949  int numberCandidates = 0;
    968950#ifdef CLP_DEBUG
    969      bool wasInValuesPass = (givenDuals != NULL);
    970 #endif
    971      int candidate = -1;
    972      if (givenDuals) {
    973           assert (ifValuesPass);
    974           ifValuesPass = 1;
    975           candidateList = new int[numberRows_];
    976           // move reduced costs across
    977           CoinMemcpyN(givenDuals, numberRows_ + numberColumns_, dj_);
    978           int iRow;
    979           for (iRow = 0; iRow < numberRows_; iRow++) {
    980                int iPivot = pivotVariable_[iRow];
    981                if (flagged(iPivot))
    982                     continue;
    983                if (fabs(dj_[iPivot]) > dualTolerance_) {
    984                     // for now safer to ignore free ones
    985                     if (lower_[iPivot] > -1.0e50 || upper_[iPivot] < 1.0e50)
    986                          if (pivoted(iPivot))
    987                               candidateList[numberCandidates++] = iRow;
    988                } else {
    989                     clearPivoted(iPivot);
    990                }
    991           }
    992           // and set first candidate
    993           if (!numberCandidates) {
    994                delete [] candidateList;
    995                delete [] givenDuals;
    996                givenDuals = NULL;
    997                candidateList = NULL;
    998                int iRow;
    999                for (iRow = 0; iRow < numberRows_; iRow++) {
    1000                     int iPivot = pivotVariable_[iRow];
    1001                     clearPivoted(iPivot);
    1002                }
    1003           }
    1004      } else {
    1005           assert (!ifValuesPass);
    1006      }
     951  bool wasInValuesPass = (givenDuals != NULL);
     952#endif
     953  int candidate = -1;
     954  if (givenDuals) {
     955    assert(ifValuesPass);
     956    ifValuesPass = 1;
     957    candidateList = new int[numberRows_];
     958    // move reduced costs across
     959    CoinMemcpyN(givenDuals, numberRows_ + numberColumns_, dj_);
     960    int iRow;
     961    for (iRow = 0; iRow < numberRows_; iRow++) {
     962      int iPivot = pivotVariable_[iRow];
     963      if (flagged(iPivot))
     964        continue;
     965      if (fabs(dj_[iPivot]) > dualTolerance_) {
     966        // for now safer to ignore free ones
     967        if (lower_[iPivot] > -1.0e50 || upper_[iPivot] < 1.0e50)
     968          if (pivoted(iPivot))
     969            candidateList[numberCandidates++] = iRow;
     970      } else {
     971        clearPivoted(iPivot);
     972      }
     973    }
     974    // and set first candidate
     975    if (!numberCandidates) {
     976      delete[] candidateList;
     977      delete[] givenDuals;
     978      givenDuals = NULL;
     979      candidateList = NULL;
     980      int iRow;
     981      for (iRow = 0; iRow < numberRows_; iRow++) {
     982        int iPivot = pivotVariable_[iRow];
     983        clearPivoted(iPivot);
     984      }
     985    }
     986  } else {
     987    assert(!ifValuesPass);
     988  }
    1007989#ifdef CHECK_ACCURACY
    1008      {
    1009           if (numberIterations_) {
    1010                int il = -1;
    1011                double largest = 1.0e-1;
    1012                int ilnb = -1;
    1013                double largestnb = 1.0e-8;
    1014                for (int i = 0; i < numberRows_ + numberColumns_; i++) {
    1015                     double diff = fabs(solution_[i] - zzzzzz[i]);
    1016                     if (diff > largest) {
    1017                          largest = diff;
    1018                          il = i;
    1019                     }
    1020                     if (getColumnStatus(i) != basic) {
    1021                          if (diff > largestnb) {
    1022                               largestnb = diff;
    1023                               ilnb = i;
    1024                          }
    1025                     }
    1026                }
    1027                if (il >= 0 && ilnb < 0)
    1028                     printf("largest diff of %g at %d, nonbasic %g at %d\n",
    1029                            largest, il, largestnb, ilnb);
    1030           }
    1031      }
    1032 #endif
    1033      while (problemStatus_ == -1) {
    1034           //if (numberIterations_>=101624)
    1035           //resetFakeBounds(-1);
     990  {
     991    if (numberIterations_) {
     992      int il = -1;
     993      double largest = 1.0e-1;
     994      int ilnb = -1;
     995      double largestnb = 1.0e-8;
     996      for (int i = 0; i < numberRows_ + numberColumns_; i++) {
     997        double diff = fabs(solution_[i] - zzzzzz[i]);
     998        if (diff > largest) {
     999          largest = diff;
     1000          il = i;
     1001        }
     1002        if (getColumnStatus(i) != basic) {
     1003          if (diff > largestnb) {
     1004            largestnb = diff;
     1005            ilnb = i;
     1006          }
     1007        }
     1008      }
     1009      if (il >= 0 && ilnb < 0)
     1010        printf("largest diff of %g at %d, nonbasic %g at %d\n",
     1011          largest, il, largestnb, ilnb);
     1012    }
     1013  }
     1014#endif
     1015  while (problemStatus_ == -1) {
     1016    //if (numberIterations_>=101624)
     1017    //resetFakeBounds(-1);
    10361018#ifdef CLP_DEBUG
    1037           if (givenDuals) {
    1038                double value5 = 0.0;
    1039                int i;
    1040                for (i = 0; i < numberRows_ + numberColumns_; i++) {
    1041                     if (dj_[i] < -1.0e-6)
    1042                          if (upper_[i] < 1.0e20)
    1043                               value5 += dj_[i] * upper_[i];
    1044                          else
    1045                               printf("bad dj %g on %d with large upper status %d\n",
    1046                                      dj_[i], i, status_[i] & 7);
    1047                     else if (dj_[i] > 1.0e-6)
    1048                          if (lower_[i] > -1.0e20)
    1049                               value5 += dj_[i] * lower_[i];
    1050                          else
    1051                               printf("bad dj %g on %d with large lower status %d\n",
    1052                                      dj_[i], i, status_[i] & 7);
    1053                }
    1054                printf("Values objective Value %g\n", value5);
    1055           }
    1056           if ((handler_->logLevel() & 32) && wasInValuesPass) {
    1057                double value5 = 0.0;
    1058                int i;
    1059                for (i = 0; i < numberRows_ + numberColumns_; i++) {
    1060                     if (dj_[i] < -1.0e-6)
    1061                          if (upper_[i] < 1.0e20)
    1062                               value5 += dj_[i] * upper_[i];
    1063                          else if (dj_[i] > 1.0e-6)
    1064                               if (lower_[i] > -1.0e20)
    1065                                    value5 += dj_[i] * lower_[i];
    1066                }
    1067                printf("Values objective Value %g\n", value5);
    1068                {
    1069                     int i;
    1070                     for (i = 0; i < numberRows_ + numberColumns_; i++) {
    1071                          int iSequence = i;
    1072                          double oldValue;
    1073 
    1074                          switch(getStatus(iSequence)) {
    1075 
    1076                          case basic:
    1077                          case ClpSimplex::isFixed:
    1078                               break;
    1079                          case isFree:
    1080                          case superBasic:
    1081                               abort();
    1082                               break;
    1083                          case atUpperBound:
    1084                               oldValue = dj_[iSequence];
    1085                               //assert (oldValue<=tolerance);
    1086                               assert (fabs(solution_[iSequence] - upper_[iSequence]) < 1.0e-7);
    1087                               break;
    1088                          case atLowerBound:
    1089                               oldValue = dj_[iSequence];
    1090                               //assert (oldValue>=-tolerance);
    1091                               assert (fabs(solution_[iSequence] - lower_[iSequence]) < 1.0e-7);
    1092                               break;
    1093                          }
    1094                     }
    1095                }
    1096           }
     1019    if (givenDuals) {
     1020      double value5 = 0.0;
     1021      int i;
     1022      for (i = 0; i < numberRows_ + numberColumns_; i++) {
     1023        if (dj_[i] < -1.0e-6)
     1024          if (upper_[i] < 1.0e20)
     1025            value5 += dj_[i] * upper_[i];
     1026          else
     1027            printf("bad dj %g on %d with large upper status %d\n",
     1028              dj_[i], i, status_[i] & 7);
     1029        else if (dj_[i] > 1.0e-6)
     1030          if (lower_[i] > -1.0e20)
     1031            value5 += dj_[i] * lower_[i];
     1032          else
     1033            printf("bad dj %g on %d with large lower status %d\n",
     1034              dj_[i], i, status_[i] & 7);
     1035      }
     1036      printf("Values objective Value %g\n", value5);
     1037    }
     1038    if ((handler_->logLevel() & 32) && wasInValuesPass) {
     1039      double value5 = 0.0;
     1040      int i;
     1041      for (i = 0; i < numberRows_ + numberColumns_; i++) {
     1042        if (dj_[i] < -1.0e-6)
     1043          if (upper_[i] < 1.0e20)
     1044            value5 += dj_[i] * upper_[i];
     1045          else if (dj_[i] > 1.0e-6)
     1046            if (lower_[i] > -1.0e20)
     1047              value5 += dj_[i] * lower_[i];
     1048      }
     1049      printf("Values objective Value %g\n", value5);
     1050      {
     1051        int i;
     1052        for (i = 0; i < numberRows_ + numberColumns_; i++) {
     1053          int iSequence = i;
     1054          double oldValue;
     1055
     1056          switch (getStatus(iSequence)) {
     1057
     1058          case basic:
     1059          case ClpSimplex::isFixed:
     1060            break;
     1061          case isFree:
     1062          case superBasic:
     1063            abort();
     1064            break;
     1065          case atUpperBound:
     1066            oldValue = dj_[iSequence];
     1067            //assert (oldValue<=tolerance);
     1068            assert(fabs(solution_[iSequence] - upper_[iSequence]) < 1.0e-7);
     1069            break;
     1070          case atLowerBound:
     1071            oldValue = dj_[iSequence];
     1072            //assert (oldValue>=-tolerance);
     1073            assert(fabs(solution_[iSequence] - lower_[iSequence]) < 1.0e-7);
     1074            break;
     1075          }
     1076        }
     1077      }
     1078    }
    10971079#endif
    10981080#ifdef CLP_DEBUG
    1099           {
    1100                int i;
    1101                for (i = 0; i < 4; i++) {
    1102                     rowArray_[i]->checkClear();
    1103                }
    1104                for (i = 0; i < 2; i++) {
    1105                     columnArray_[i]->checkClear();
    1106                }
    1107           }
    1108 #endif
    1109 #if CLP_DEBUG>2
    1110           // very expensive
    1111           if (numberIterations_ > 3063 && numberIterations_ < 30700) {
    1112                //handler_->setLogLevel(63);
    1113                double saveValue = objectiveValue_;
    1114                double * saveRow1 = new double[numberRows_];
    1115                double * saveRow2 = new double[numberRows_];
    1116                CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
    1117                CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
    1118                double * saveColumn1 = new double[numberColumns_];
    1119                double * saveColumn2 = new double[numberColumns_];
    1120                CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
    1121                CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
    1122                gutsOfSolution(NULL, NULL);
    1123                printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
    1124                       numberIterations_,
    1125                       saveValue, objectiveValue_, sumDualInfeasibilities_);
    1126                if (saveValue > objectiveValue_ + 1.0e-2)
    1127                     printf("**bad**\n");
    1128                CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
    1129                CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
    1130                CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
    1131                CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
    1132                delete [] saveRow1;
    1133                delete [] saveRow2;
    1134                delete [] saveColumn1;
    1135                delete [] saveColumn2;
    1136                objectiveValue_ = saveValue;
    1137           }
     1081    {
     1082      int i;
     1083      for (i = 0; i < 4; i++) {
     1084        rowArray_[i]->checkClear();
     1085      }
     1086      for (i = 0; i < 2; i++) {
     1087        columnArray_[i]->checkClear();
     1088      }
     1089    }
     1090#endif
     1091#if CLP_DEBUG > 2
     1092    // very expensive
     1093    if (numberIterations_ > 3063 && numberIterations_ < 30700) {
     1094      //handler_->setLogLevel(63);
     1095      double saveValue = objectiveValue_;
     1096      double *saveRow1 = new double[numberRows_];
     1097      double *saveRow2 = new double[numberRows_];
     1098      CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
     1099      CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
     1100      double *saveColumn1 = new double[numberColumns_];
     1101      double *saveColumn2 = new double[numberColumns_];
     1102      CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
     1103      CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
     1104      gutsOfSolution(NULL, NULL);
     1105      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
     1106        numberIterations_,
     1107        saveValue, objectiveValue_, sumDualInfeasibilities_);
     1108      if (saveValue > objectiveValue_ + 1.0e-2)
     1109        printf("**bad**\n");
     1110      CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
     1111      CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
     1112      CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
     1113      CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
     1114      delete[] saveRow1;
     1115      delete[] saveRow2;
     1116      delete[] saveColumn1;
     1117      delete[] saveColumn2;
     1118      objectiveValue_ = saveValue;
     1119    }
    11381120#endif
    11391121#if 0
     
    11561138#endif
    11571139#ifdef CLP_DEBUG
    1158           {
    1159                int iSequence, number = numberRows_ + numberColumns_;
    1160                for (iSequence = 0; iSequence < number; iSequence++) {
    1161                     double lowerValue = lower_[iSequence];
    1162                     double upperValue = upper_[iSequence];
    1163                     double value = solution_[iSequence];
    1164                     if(getStatus(iSequence) != basic && getStatus(iSequence) != isFree) {
    1165                          assert(lowerValue > -1.0e20);
    1166                          assert(upperValue < 1.0e20);
    1167                     }
    1168                     switch(getStatus(iSequence)) {
    1169 
    1170                     case basic:
    1171                          break;
    1172                     case isFree:
    1173                     case superBasic:
    1174                          break;
    1175                     case atUpperBound:
    1176                          assert (fabs(value - upperValue) <= primalTolerance_) ;
    1177                          break;
    1178                     case atLowerBound:
    1179                     case ClpSimplex::isFixed:
    1180                          assert (fabs(value - lowerValue) <= primalTolerance_) ;
    1181                          break;
    1182                     }
    1183                }
    1184           }
    1185           if(numberIterations_ == debugIteration) {
    1186                printf("dodgy iteration coming up\n");
    1187           }
     1140    {
     1141      int iSequence, number = numberRows_ + numberColumns_;
     1142      for (iSequence = 0; iSequence < number; iSequence++) {
     1143        double lowerValue = lower_[iSequence];
     1144        double upperValue = upper_[iSequence];
     1145        double value = solution_[iSequence];
     1146        if (getStatus(iSequence) != basic && getStatus(iSequence) != isFree) {
     1147          assert(lowerValue > -1.0e20);
     1148          assert(upperValue < 1.0e20);
     1149        }
     1150        switch (getStatus(iSequence)) {
     1151
     1152        case basic:
     1153          break;
     1154        case isFree:
     1155        case superBasic:
     1156          break;
     1157        case atUpperBound:
     1158          assert(fabs(value - upperValue) <= primalTolerance_);
     1159          break;
     1160        case atLowerBound:
     1161        case ClpSimplex::isFixed:
     1162          assert(fabs(value - lowerValue) <= primalTolerance_);
     1163          break;
     1164        }
     1165      }
     1166    }
     1167    if (numberIterations_ == debugIteration) {
     1168      printf("dodgy iteration coming up\n");
     1169    }
    11881170#endif
    11891171#if 0
     
    11961178          }
    11971179#endif
    1198           // choose row to go out
    1199           // dualRow will go to virtual row pivot choice algorithm
    1200           // make sure values pass off if it should be
    1201           if (numberCandidates)
    1202                candidate = candidateList[--numberCandidates];
    1203           else
    1204                candidate = -1;
    1205           dualRow(candidate);
    1206           if (pivotRow_ >= 0) {
     1180    // choose row to go out
     1181    // dualRow will go to virtual row pivot choice algorithm
     1182    // make sure values pass off if it should be
     1183    if (numberCandidates)
     1184      candidate = candidateList[--numberCandidates];
     1185    else
     1186      candidate = -1;
     1187    dualRow(candidate);
     1188    if (pivotRow_ >= 0) {
    12071189#if ABOCA_LITE_FACTORIZATION
    1208             int numberThreads=abcState();
    1209             if (numberThreads)
    1210               cilk_spawn factorization_->replaceColumn1(columnArray_[1],
    1211                                                       pivotRow_);
    1212 #endif
    1213                // we found a pivot row
    1214                if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
    1215                     handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
    1216                               << pivotRow_
    1217                               << CoinMessageEol;
    1218                }
    1219                // check accuracy of weights
    1220                dualRowPivot_->checkAccuracy();
    1221                // Get good size for pivot
    1222                // Allow first few iterations to take tiny
    1223                double acceptablePivot = 1.0e-1 * acceptablePivot_;
    1224                if (numberIterations_ > 100)
    1225                     acceptablePivot = acceptablePivot_;
    1226                if (factorization_->pivots() > 10 ||
    1227                          (factorization_->pivots() && saveSumDual))
    1228                     acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
    1229                else if (factorization_->pivots() > 5)
    1230                     acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
    1231                else if (factorization_->pivots())
    1232                     acceptablePivot = acceptablePivot_; // relax
    1233                // But factorizations complain if <1.0e-8
    1234                //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
    1235                double bestPossiblePivot = 1.0;
    1236                // get sign for finding row of tableau
    1237                if (candidate < 0) {
    1238                     // normal iteration
    1239                     // create as packed
    1240                     double direction = directionOut_;
    1241                     rowArray_[0]->createPacked(1, &pivotRow_, &direction);
    1242                     factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
    1243                     // Allow to do dualColumn0
    1244                     if (numberThreads_ < -1)
    1245                          spareIntArray_[0] = 1;
    1246                     spareDoubleArray_[0] = acceptablePivot;
    1247                     rowArray_[3]->clear();
    1248                     sequenceIn_ = -1;
    1249                     // put row of tableau in rowArray[0] and columnArray[0]
    1250                     assert (!rowArray_[1]->getNumElements());
    1251                     if (!scaledMatrix_) {
    1252                          if ((moreSpecialOptions_ & 8) != 0 && !rowScale_)
    1253                               spareIntArray_[0] = 1;
    1254                          matrix_->transposeTimes(this, -1.0,
    1255                                                  rowArray_[0], rowArray_[1], columnArray_[0]);
    1256                     } else {
    1257                          double * saveR = rowScale_;
    1258                          double * saveC = columnScale_;
    1259                          rowScale_ = NULL;
    1260                          columnScale_ = NULL;
    1261                          if ((moreSpecialOptions_ & 8) != 0)
    1262                               spareIntArray_[0] = 1;
    1263                          scaledMatrix_->transposeTimes(this, -1.0,
    1264                                                        rowArray_[0], rowArray_[1], columnArray_[0]);
    1265                          rowScale_ = saveR;
    1266                          columnScale_ = saveC;
    1267                     }
     1190      int numberThreads = abcState();
     1191      if (numberThreads)
     1192        cilk_spawn factorization_->replaceColumn1(columnArray_[1],
     1193          pivotRow_);
     1194#endif
     1195      // we found a pivot row
     1196      if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
     1197        handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
     1198          << pivotRow_
     1199          << CoinMessageEol;
     1200      }
     1201      // check accuracy of weights
     1202      dualRowPivot_->checkAccuracy();
     1203      // Get good size for pivot
     1204      // Allow first few iterations to take tiny
     1205      double acceptablePivot = 1.0e-1 * acceptablePivot_;
     1206      if (numberIterations_ > 100)
     1207        acceptablePivot = acceptablePivot_;
     1208      if (factorization_->pivots() > 10 || (factorization_->pivots() && saveSumDual))
     1209        acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
     1210      else if (factorization_->pivots() > 5)
     1211        acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
     1212      else if (factorization_->pivots())
     1213        acceptablePivot = acceptablePivot_; // relax
     1214      // But factorizations complain if <1.0e-8
     1215      //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
     1216      double bestPossiblePivot = 1.0;
     1217      // get sign for finding row of tableau
     1218      if (candidate < 0) {
     1219        // normal iteration
     1220        // create as packed
     1221        double direction = directionOut_;
     1222        rowArray_[0]->createPacked(1, &pivotRow_, &direction);
     1223        factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
     1224        // Allow to do dualColumn0
     1225        if (numberThreads_ < -1)
     1226          spareIntArray_[0] = 1;
     1227        spareDoubleArray_[0] = acceptablePivot;
     1228        rowArray_[3]->clear();
     1229        sequenceIn_ = -1;
     1230        // put row of tableau in rowArray[0] and columnArray[0]
     1231        assert(!rowArray_[1]->getNumElements());
     1232        if (!scaledMatrix_) {
     1233          if ((moreSpecialOptions_ & 8) != 0 && !rowScale_)
     1234            spareIntArray_[0] = 1;
     1235          matrix_->transposeTimes(this, -1.0,
     1236            rowArray_[0], rowArray_[1], columnArray_[0]);
     1237        } else {
     1238          double *saveR = rowScale_;
     1239          double *saveC = columnScale_;
     1240          rowScale_ = NULL;
     1241          columnScale_ = NULL;
     1242          if ((moreSpecialOptions_ & 8) != 0)
     1243            spareIntArray_[0] = 1;
     1244          scaledMatrix_->transposeTimes(this, -1.0,
     1245            rowArray_[0], rowArray_[1], columnArray_[0]);
     1246          rowScale_ = saveR;
     1247          columnScale_ = saveC;
     1248        }
    12681249#ifdef CLP_REPORT_PROGRESS
    1269                     memcpy(savePSol, solution_, (numberColumns_ + numberRows_)*sizeof(double));
    1270                     memcpy(saveDj, dj_, (numberColumns_ + numberRows_)*sizeof(double));
    1271                     memcpy(saveCost, cost_, (numberColumns_ + numberRows_)*sizeof(double));
    1272                     memcpy(saveStat, status_, (numberColumns_ + numberRows_)*sizeof(char));
    1273 #endif
    1274                     // do ratio test for normal iteration
    1275                     bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
    1276 #ifdef LONG_REGION_2 
    1277                                                    rowArray_[2],
    1278 #else 
    1279                                                    columnArray_[1],
    1280 #endif 
    1281                                                    acceptablePivot, dubiousWeights);
    1282                     if (sequenceIn_<0&&acceptablePivot>acceptablePivot_)
    1283                       acceptablePivot_ = - fabs(acceptablePivot_); // stop early exit
    1284 #if CAN_HAVE_ZERO_OBJ>1
    1285                     if ((specialOptions_&16777216)!=0)
    1286                       theta_=0.0;
    1287 #endif
    1288                } else {
    1289                     // Make sure direction plausible
    1290                     CoinAssert (upperOut_ < 1.0e50 || lowerOut_ > -1.0e50);
    1291                     // If in integer cleanup do direction using duals
    1292                     // may be wrong way round
    1293                     if(ifValuesPass == 2) {
    1294                          if (dual_[pivotRow_] > 0.0) {
    1295                               // this will give a -1 in pivot row (as slacks are -1.0)
    1296                               directionOut_ = 1;
    1297                          } else {
    1298                               directionOut_ = -1;
    1299                          }
    1300                     }
    1301                     if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) {
    1302                          if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_))
    1303                               directionOut_ = 1;
    1304                     } else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) {
    1305                          if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_))
    1306                               directionOut_ = -1;
    1307                     }
    1308                     double direction = directionOut_;
    1309                     rowArray_[0]->createPacked(1, &pivotRow_, &direction);
    1310                     factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
    1311                     // put row of tableau in rowArray[0] and columnArray[0]
    1312                     if (!scaledMatrix_) {
    1313                          matrix_->transposeTimes(this, -1.0,
    1314                                                  rowArray_[0], rowArray_[3], columnArray_[0]);
    1315                     } else {
    1316                          double * saveR = rowScale_;
    1317                          double * saveC = columnScale_;
    1318                          rowScale_ = NULL;
    1319                          columnScale_ = NULL;
    1320                          scaledMatrix_->transposeTimes(this, -1.0,
    1321                                                        rowArray_[0], rowArray_[3], columnArray_[0]);
    1322                          rowScale_ = saveR;
    1323                          columnScale_ = saveC;
    1324                     }
    1325                     acceptablePivot *= 10.0;
    1326                     // do ratio test
    1327                     if (ifValuesPass == 1) {
    1328                          checkPossibleValuesMove(rowArray_[0], columnArray_[0],
    1329                                                  acceptablePivot);
    1330                     } else {
    1331                          checkPossibleCleanup(rowArray_[0], columnArray_[0],
    1332                                               acceptablePivot);
    1333                          if (sequenceIn_ < 0) {
    1334                               rowArray_[0]->clear();
    1335                               columnArray_[0]->clear();
    1336                               continue; // can't do anything
    1337                          }
    1338                     }
    1339 
    1340                     // recompute true dualOut_
    1341                     if (directionOut_ < 0) {
    1342                          dualOut_ = valueOut_ - upperOut_;
    1343                     } else {
    1344                          dualOut_ = lowerOut_ - valueOut_;
    1345                     }
    1346                     // check what happened if was values pass
    1347                     // may want to move part way i.e. movement
    1348                     bool normalIteration = (sequenceIn_ != sequenceOut_);
    1349 
    1350                     clearPivoted(sequenceOut_); // make sure won't be done again
    1351                     // see if end of values pass
    1352                     if (!numberCandidates) {
    1353                          int iRow;
    1354                          delete [] candidateList;
    1355                          delete [] givenDuals;
    1356                          candidate = -2; // -2 signals end
    1357                          givenDuals = NULL;
    1358                          candidateList = NULL;
    1359                          ifValuesPass = 1;
    1360                          for (iRow = 0; iRow < numberRows_; iRow++) {
    1361                               int iPivot = pivotVariable_[iRow];
    1362                               //assert (fabs(dj_[iPivot]),1.0e-5);
    1363                               clearPivoted(iPivot);
    1364                          }
    1365                     }
    1366                     if (!normalIteration) {
    1367                          //rowArray_[0]->cleanAndPackSafe(1.0e-60);
    1368                          //columnArray_[0]->cleanAndPackSafe(1.0e-60);
    1369                          updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
    1370                          if (candidate == -2)
    1371                               problemStatus_ = -2;
    1372                          continue; // skip rest of iteration
    1373                     } else {
    1374                          // recompute dualOut_
    1375                          if (directionOut_ < 0) {
    1376                               dualOut_ = valueOut_ - upperOut_;
    1377                          } else {
    1378                               dualOut_ = lowerOut_ - valueOut_;
    1379                          }
    1380                     }
    1381                }
    1382                if (sequenceIn_ >= 0) {
    1383                     // normal iteration
    1384                     // update the incoming column
    1385                     double btranAlpha = -alpha_ * directionOut_; // for check
    1386                     unpackPacked(rowArray_[1]);
    1387                     // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
    1388                     // and update dual weights (can do in parallel - with extra array)
    1389                     alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
    1390                                                           rowArray_[2],
    1391                                                           rowArray_[3],
    1392                                                           rowArray_[1]);
    1393                     // see if update stable
     1250        memcpy(savePSol, solution_, (numberColumns_ + numberRows_) * sizeof(double));
     1251        memcpy(saveDj, dj_, (numberColumns_ + numberRows_) * sizeof(double));
     1252        memcpy(saveCost, cost_, (numberColumns_ + numberRows_) * sizeof(double));
     1253        memcpy(saveStat, status_, (numberColumns_ + numberRows_) * sizeof(char));
     1254#endif
     1255        // do ratio test for normal iteration
     1256        bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
     1257#ifdef LONG_REGION_2
     1258          rowArray_[2],
     1259#else
     1260          columnArray_[1],
     1261#endif
     1262          acceptablePivot, dubiousWeights);
     1263        if (sequenceIn_ < 0 && acceptablePivot > acceptablePivot_)
     1264          acceptablePivot_ = -fabs(acceptablePivot_); // stop early exit
     1265#if CAN_HAVE_ZERO_OBJ > 1
     1266        if ((specialOptions_ & 16777216) != 0)
     1267          theta_ = 0.0;
     1268#endif
     1269      } else {
     1270        // Make sure direction plausible
     1271        CoinAssert(upperOut_ < 1.0e50 || lowerOut_ > -1.0e50);
     1272        // If in integer cleanup do direction using duals
     1273        // may be wrong way round
     1274        if (ifValuesPass == 2) {
     1275          if (dual_[pivotRow_] > 0.0) {
     1276            // this will give a -1 in pivot row (as slacks are -1.0)
     1277            directionOut_ = 1;
     1278          } else {
     1279            directionOut_ = -1;
     1280          }
     1281        }
     1282        if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) {
     1283          if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_))
     1284            directionOut_ = 1;
     1285        } else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) {
     1286          if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_))
     1287            directionOut_ = -1;
     1288        }
     1289        double direction = directionOut_;
     1290        rowArray_[0]->createPacked(1, &pivotRow_, &direction);
     1291        factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
     1292        // put row of tableau in rowArray[0] and columnArray[0]
     1293        if (!scaledMatrix_) {
     1294          matrix_->transposeTimes(this, -1.0,
     1295            rowArray_[0], rowArray_[3], columnArray_[0]);
     1296        } else {
     1297          double *saveR = rowScale_;
     1298          double *saveC = columnScale_;
     1299          rowScale_ = NULL;
     1300          columnScale_ = NULL;
     1301          scaledMatrix_->transposeTimes(this, -1.0,
     1302            rowArray_[0], rowArray_[3], columnArray_[0]);
     1303          rowScale_ = saveR;
     1304          columnScale_ = saveC;
     1305        }
     1306        acceptablePivot *= 10.0;
     1307        // do ratio test
     1308        if (ifValuesPass == 1) {
     1309          checkPossibleValuesMove(rowArray_[0], columnArray_[0],
     1310            acceptablePivot);
     1311        } else {
     1312          checkPossibleCleanup(rowArray_[0], columnArray_[0],
     1313            acceptablePivot);
     1314          if (sequenceIn_ < 0) {
     1315            rowArray_[0]->clear();
     1316            columnArray_[0]->clear();
     1317            continue; // can't do anything
     1318          }
     1319        }
     1320
     1321        // recompute true dualOut_
     1322        if (directionOut_ < 0) {
     1323          dualOut_ = valueOut_ - upperOut_;
     1324        } else {
     1325          dualOut_ = lowerOut_ - valueOut_;
     1326        }
     1327        // check what happened if was values pass
     1328        // may want to move part way i.e. movement
     1329        bool normalIteration = (sequenceIn_ != sequenceOut_);
     1330
     1331        clearPivoted(sequenceOut_); // make sure won't be done again
     1332        // see if end of values pass
     1333        if (!numberCandidates) {
     1334          int iRow;
     1335          delete[] candidateList;
     1336          delete[] givenDuals;
     1337          candidate = -2; // -2 signals end
     1338          givenDuals = NULL;
     1339          candidateList = NULL;
     1340          ifValuesPass = 1;
     1341          for (iRow = 0; iRow < numberRows_; iRow++) {
     1342            int iPivot = pivotVariable_[iRow];
     1343            //assert (fabs(dj_[iPivot]),1.0e-5);
     1344            clearPivoted(iPivot);
     1345          }
     1346        }
     1347        if (!normalIteration) {
     1348          //rowArray_[0]->cleanAndPackSafe(1.0e-60);
     1349          //columnArray_[0]->cleanAndPackSafe(1.0e-60);
     1350          updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
     1351          if (candidate == -2)
     1352            problemStatus_ = -2;
     1353          continue; // skip rest of iteration
     1354        } else {
     1355          // recompute dualOut_
     1356          if (directionOut_ < 0) {
     1357            dualOut_ = valueOut_ - upperOut_;
     1358          } else {
     1359            dualOut_ = lowerOut_ - valueOut_;
     1360          }
     1361        }
     1362      }
     1363      if (sequenceIn_ >= 0) {
     1364        // normal iteration
     1365        // update the incoming column
     1366        double btranAlpha = -alpha_ * directionOut_; // for check
     1367        unpackPacked(rowArray_[1]);
     1368        // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
     1369        // and update dual weights (can do in parallel - with extra array)
     1370        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
     1371          rowArray_[2],
     1372          rowArray_[3],
     1373          rowArray_[1]);
     1374        // see if update stable
    13941375#ifdef CLP_DEBUG
    1395                     if ((handler_->logLevel() & 32))
    1396                          printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
    1397 #endif
    1398                     double checkValue = 1.0e-7;
    1399                     // if can't trust much and long way from optimal then relax
    1400                     if (largestPrimalError_ > 10.0)
    1401                          checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
    1402                     if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
    1403                               fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
    1404                          handler_->message(CLP_DUAL_CHECK, messages_)
    1405                                    << btranAlpha
    1406                                    << alpha_
    1407                                    << CoinMessageEol;
    1408                          if (factorization_->pivots()) {
    1409                               dualRowPivot_->unrollWeights();
    1410                               problemStatus_ = -2; // factorize now
    1411                               rowArray_[0]->clear();
    1412                               rowArray_[1]->clear();
    1413                               columnArray_[0]->clear();
    1414                               returnCode = -2;
    1415                               break;
    1416                          } else {
    1417                               // take on more relaxed criterion
    1418                               double test;
    1419                               if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
    1420                                    test = 1.0e-1 * fabs(alpha_);
    1421                               else
    1422                                    test = 1.0e-4 * (1.0 + fabs(alpha_));
    1423                               if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
    1424                                         fabs(btranAlpha - alpha_) > test) {
    1425                                    dualRowPivot_->unrollWeights();
    1426                                    // need to reject something
    1427                                    char x = isColumn(sequenceOut_) ? 'C' : 'R';
    1428                                    handler_->message(CLP_SIMPLEX_FLAG, messages_)
    1429                                              << x << sequenceWithin(sequenceOut_)
    1430                                              << CoinMessageEol;
     1376        if ((handler_->logLevel() & 32))
     1377          printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
     1378#endif
     1379        double checkValue = 1.0e-7;
     1380        // if can't trust much and long way from optimal then relax
     1381        if (largestPrimalError_ > 10.0)
     1382          checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
     1383        if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha - alpha_) > checkValue * (1.0 + fabs(alpha_))) {
     1384          handler_->message(CLP_DUAL_CHECK, messages_)
     1385            << btranAlpha
     1386            << alpha_
     1387            << CoinMessageEol;
     1388          if (factorization_->pivots()) {
     1389            dualRowPivot_->unrollWeights();
     1390            problemStatus_ = -2; // factorize now
     1391            rowArray_[0]->clear();
     1392            rowArray_[1]->clear();
     1393            columnArray_[0]->clear();
     1394            returnCode = -2;
     1395            break;
     1396          } else {
     1397            // take on more relaxed criterion
     1398            double test;
     1399            if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
     1400              test = 1.0e-1 * fabs(alpha_);
     1401            else
     1402              test = 1.0e-4 * (1.0 + fabs(alpha_));
     1403            if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha - alpha_) > test) {
     1404              dualRowPivot_->unrollWeights();
     1405              // need to reject something
     1406              char x = isColumn(sequenceOut_) ? 'C' : 'R';
     1407              handler_->message(CLP_SIMPLEX_FLAG, messages_)
     1408                << x << sequenceWithin(sequenceOut_)
     1409                << CoinMessageEol;
    14311410#ifdef COIN_DEVELOP
    1432                                    printf("flag a %g %g\n", btranAlpha, alpha_);
    1433 #endif
    1434                                    //#define FEB_TRY
    1435 #if 1 //def FEB_TRY
    1436                                    // Make safer?
    1437                                    factorization_->saferTolerances (-0.99, -1.03);
    1438 #endif
    1439                                    setFlagged(sequenceOut_);
    1440                                    progress_.clearBadTimes();
    1441                                    lastBadIteration_ = numberIterations_; // say be more cautious
    1442                                    rowArray_[0]->clear();
    1443                                    rowArray_[1]->clear();
    1444                                    columnArray_[0]->clear();
    1445                                    if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
    1446                                         //printf("I think should declare infeasible\n");
    1447                                         problemStatus_ = 1;
    1448                                         returnCode = 1;
    1449                                         break;
    1450                                    }
    1451                                    continue;
    1452                               }
    1453                          }
    1454                     }
    1455                     // update duals BEFORE replaceColumn so can do updateColumn
    1456                     double objectiveChange = 0.0;
    1457                     // do duals first as variables may flip bounds
    1458                     // rowArray_[0] and columnArray_[0] may have flips
    1459                     // so use rowArray_[3] for work array from here on
    1460                     int nswapped = 0;
    1461                     //rowArray_[0]->cleanAndPackSafe(1.0e-60);
    1462                     //columnArray_[0]->cleanAndPackSafe(1.0e-60);
    1463                     if (candidate == -1) {
    1464 #if CLP_CAN_HAVE_ZERO_OBJ>1
    1465                     if ((specialOptions_&16777216)==0) {
    1466 #endif
    1467                          // make sure incoming doesn't count
    1468                          Status saveStatus = getStatus(sequenceIn_);
    1469                          setStatus(sequenceIn_, basic);
    1470                          nswapped = updateDualsInDual(rowArray_[0], columnArray_[0],
    1471                                                       rowArray_[2], theta_,
    1472                                                       objectiveChange, false);
    1473                          setStatus(sequenceIn_, saveStatus);
    1474 #if CLP_CAN_HAVE_ZERO_OBJ>1
    1475                     } else {
    1476                       rowArray_[0]->clear();
    1477                       rowArray_[2]->clear();
    1478                       columnArray_[0]->clear();
    1479                     }
    1480 #endif
    1481                     } else {
    1482                          updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
    1483                     }
    1484                     double oldDualOut = dualOut_;
    1485                     // which will change basic solution
    1486                     if (nswapped) {
    1487                          if (rowArray_[2]->getNumElements()) {
    1488                               factorization_->updateColumn(rowArray_[3], rowArray_[2]);
    1489                               dualRowPivot_->updatePrimalSolution(rowArray_[2],
    1490                                                                   1.0, objectiveChange);
    1491                          }
    1492                          // recompute dualOut_
    1493                          valueOut_ = solution_[sequenceOut_];
    1494                          if (directionOut_ < 0) {
    1495                               dualOut_ = valueOut_ - upperOut_;
    1496                          } else {
    1497                               dualOut_ = lowerOut_ - valueOut_;
    1498                          }
     1411              printf("flag a %g %g\n", btranAlpha, alpha_);
     1412#endif
     1413              //#define FEB_TRY
     1414#if 1 //def FEB_TRY \
     1415  // Make safer?
     1416              factorization_->saferTolerances(-0.99, -1.03);
     1417#endif
     1418              setFlagged(sequenceOut_);
     1419              progress_.clearBadTimes();
     1420              lastBadIteration_ = numberIterations_; // say be more cautious
     1421              rowArray_[0]->clear();
     1422              rowArray_[1]->clear();
     1423              columnArray_[0]->clear();
     1424              if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
     1425                //printf("I think should declare infeasible\n");
     1426                problemStatus_ = 1;
     1427                returnCode = 1;
     1428                break;
     1429              }
     1430              continue;
     1431            }
     1432          }
     1433        }
     1434        // update duals BEFORE replaceColumn so can do updateColumn
     1435        double objectiveChange = 0.0;
     1436        // do duals first as variables may flip bounds
     1437        // rowArray_[0] and columnArray_[0] may have flips
     1438        // so use rowArray_[3] for work array from here on
     1439        int nswapped = 0;
     1440        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
     1441        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
     1442        if (candidate == -1) {
     1443#if CLP_CAN_HAVE_ZERO_OBJ > 1
     1444          if ((specialOptions_ & 16777216) == 0) {
     1445#endif
     1446            // make sure incoming doesn't count
     1447            Status saveStatus = getStatus(sequenceIn_);
     1448            setStatus(sequenceIn_, basic);
     1449            nswapped = updateDualsInDual(rowArray_[0], columnArray_[0],
     1450              rowArray_[2], theta_,
     1451              objectiveChange, false);
     1452            setStatus(sequenceIn_, saveStatus);
     1453#if CLP_CAN_HAVE_ZERO_OBJ > 1
     1454          } else {
     1455            rowArray_[0]->clear();
     1456            rowArray_[2]->clear();
     1457            columnArray_[0]->clear();
     1458          }
     1459#endif
     1460        } else {
     1461          updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
     1462        }
     1463        double oldDualOut = dualOut_;
     1464        // which will change basic solution
     1465        if (nswapped) {
     1466          if (rowArray_[2]->getNumElements()) {
     1467            factorization_->updateColumn(rowArray_[3], rowArray_[2]);
     1468            dualRowPivot_->updatePrimalSolution(rowArray_[2],
     1469              1.0, objectiveChange);
     1470          }
     1471          // recompute dualOut_
     1472          valueOut_ = solution_[sequenceOut_];
     1473          if (directionOut_ < 0) {
     1474            dualOut_ = valueOut_ - upperOut_;
     1475          } else {
     1476            dualOut_ = lowerOut_ - valueOut_;
     1477          }
    14991478#if 0
    15001479                         if (dualOut_ < 0.0) {
     
    15191498                         }
    15201499#endif
    1521                     }
    1522                     // amount primal will move
    1523                     double movement = -dualOut_ * directionOut_ / alpha_;
    1524                     double movementOld = oldDualOut * directionOut_ / alpha_;
    1525                     // so objective should increase by fabs(dj)*movement
    1526                     // but we already have objective change - so check will be good
    1527                     if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) {
     1500        }
     1501        // amount primal will move
     1502        double movement = -dualOut_ * directionOut_ / alpha_;
     1503        double movementOld = oldDualOut * directionOut_ / alpha_;
     1504        // so objective should increase by fabs(dj)*movement
     1505        // but we already have objective change - so check will be good
     1506        if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) {
    15281507#ifdef CLP_DEBUG
    1529                          if (handler_->logLevel() & 32)
    1530                               printf("movement %g, swap change %g, rest %g  * %g\n",
    1531                                      objectiveChange + fabs(movement * dualIn_),
    1532                                      objectiveChange, movement, dualIn_);
    1533 #endif
    1534                          if(factorization_->pivots()) {
    1535                               // going backwards - factorize
    1536                               dualRowPivot_->unrollWeights();
    1537                               problemStatus_ = -2; // factorize now
    1538                               returnCode = -2;
    1539                               break;
    1540                          }
    1541                     }
    1542                     // if stable replace in basis
    1543                     int updateStatus=123456789;
     1508          if (handler_->logLevel() & 32)
     1509            printf("movement %g, swap change %g, rest %g  * %g\n",
     1510              objectiveChange + fabs(movement * dualIn_),
     1511              objectiveChange, movement, dualIn_);
     1512#endif
     1513          if (factorization_->pivots()) {
     1514            // going backwards - factorize
     1515            dualRowPivot_->unrollWeights();
     1516            problemStatus_ = -2; // factorize now
     1517            returnCode = -2;
     1518            break;
     1519          }
     1520        }
     1521        // if stable replace in basis
     1522        int updateStatus = 123456789;
    15441523#if ABOCA_LITE_FACTORIZATION
    1545                     if (numberThreads)
    1546                       cilk_sync;
    1547                     if (columnArray_[1]->getNumElements())
    1548                       updateStatus = factorization_->replaceColumn2(columnArray_[1],
    1549                                                                     pivotRow_,alpha_);
    1550                     if (updateStatus==123456789)
    1551 #endif
    1552                      updateStatus = factorization_->replaceColumn(this,
    1553                                        rowArray_[2],
    1554                                        rowArray_[1],
    1555                                        pivotRow_,
    1556                                        alpha_,
    1557                                        (moreSpecialOptions_ & 16) != 0,
    1558                                        acceptablePivot);
    1559                     // If looks like bad pivot - refactorize
    1560                     if (fabs(dualOut_) > 1.0e50)
    1561                          updateStatus = 2;
    1562                     // if no pivots, bad update but reasonable alpha - take and invert
    1563                     if (updateStatus == 2 &&
    1564                               !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
    1565                          updateStatus = 4;
    1566                     if (updateStatus == 1 || updateStatus == 4) {
    1567                          // slight error
    1568                          if (factorization_->pivots() > 5 || updateStatus == 4) {
    1569                               problemStatus_ = -2; // factorize now
    1570                               returnCode = -3;
    1571                          }
    1572                     } else if (updateStatus == 2) {
    1573                          // major error
    1574                          dualRowPivot_->unrollWeights();
    1575                          // later we may need to unwind more e.g. fake bounds
    1576                          if (factorization_->pivots() &&
    1577                                    ((moreSpecialOptions_ & 16) == 0 || factorization_->pivots() > 4)) {
    1578                               problemStatus_ = -2; // factorize now
    1579                               returnCode = -2;
    1580                               moreSpecialOptions_ |= 16;
    1581                               double pivotTolerance = factorization_->pivotTolerance();
    1582                               if (pivotTolerance<0.4&&factorization_->pivots()<100) {
    1583                                 factorization_->pivotTolerance(1.05*pivotTolerance);
     1524        if (numberThreads)
     1525          cilk_sync;
     1526        if (columnArray_[1]->getNumElements())
     1527          updateStatus = factorization_->replaceColumn2(columnArray_[1],
     1528            pivotRow_, alpha_);
     1529        if (updateStatus == 123456789)
     1530#endif
     1531          updateStatus = factorization_->replaceColumn(this,
     1532            rowArray_[2],
     1533            rowArray_[1],
     1534            pivotRow_,
     1535            alpha_,
     1536            (moreSpecialOptions_ & 16) != 0,
     1537            acceptablePivot);
     1538        // If looks like bad pivot - refactorize
     1539        if (fabs(dualOut_) > 1.0e50)
     1540          updateStatus = 2;
     1541        // if no pivots, bad update but reasonable alpha - take and invert
     1542        if (updateStatus == 2 && !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
     1543          updateStatus = 4;
     1544        if (updateStatus == 1 || updateStatus == 4) {
     1545          // slight error
     1546          if (factorization_->pivots() > 5 || updateStatus == 4) {
     1547            problemStatus_ = -2; // factorize now
     1548            returnCode = -3;
     1549          }
     1550        } else if (updateStatus == 2) {
     1551          // major error
     1552          dualRowPivot_->unrollWeights();
     1553          // later we may need to unwind more e.g. fake bounds
     1554          if (factorization_->pivots() && ((moreSpecialOptions_ & 16) == 0 || factorization_->pivots() > 4)) {
     1555            problemStatus_ = -2; // factorize now
     1556            returnCode = -2;
     1557            moreSpecialOptions_ |= 16;
     1558            double pivotTolerance = factorization_->pivotTolerance();
     1559            if (pivotTolerance < 0.4 && factorization_->pivots() < 100) {
     1560              factorization_->pivotTolerance(1.05 * pivotTolerance);
    15841561#ifdef CLP_USEFUL_PRINTOUT
    1585                                 printf("Changing pivot tolerance from %g to %g as ftran/btran error %g/%g\n",
    1586                                        pivotTolerance,factorization_->pivotTolerance(),
    1587                                        alpha_,btranAlpha);
    1588 #endif
    1589                               }
    1590                               break;
    1591                          } else {
    1592                               // need to reject something
    1593                               char x = isColumn(sequenceOut_) ? 'C' : 'R';
    1594                               handler_->message(CLP_SIMPLEX_FLAG, messages_)
    1595                                         << x << sequenceWithin(sequenceOut_)
    1596                                         << CoinMessageEol;
     1562              printf("Changing pivot tolerance from %g to %g as ftran/btran error %g/%g\n",
     1563                pivotTolerance, factorization_->pivotTolerance(),
     1564                alpha_, btranAlpha);
     1565#endif
     1566            }
     1567            break;
     1568          } else {
     1569            // need to reject something
     1570            char x = isColumn(sequenceOut_) ? 'C' : 'R';
     1571            handler_->message(CLP_SIMPLEX_FLAG, messages_)
     1572              << x << sequenceWithin(sequenceOut_)
     1573              << CoinMessageEol;
    15971574#ifdef COIN_DEVELOP
    1598                               printf("flag b %g\n", alpha_);
    1599 #endif
    1600                               setFlagged(sequenceOut_);
    1601                               progress_.clearBadTimes();
    1602                               lastBadIteration_ = numberIterations_; // say be more cautious
    1603                               rowArray_[0]->clear();
    1604                               rowArray_[1]->clear();
    1605                               columnArray_[0]->clear();
    1606                               // make sure dual feasible
    1607                               // look at all rows and columns
    1608                               double objectiveChange = 0.0;
    1609                               updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
    1610                                                 0.0, objectiveChange, true);
    1611                               rowArray_[1]->clear();
    1612                               columnArray_[0]->clear();
    1613                               continue;
    1614                          }
    1615                     } else if (updateStatus == 3) {
    1616                          // out of memory
    1617                          // increase space if not many iterations
    1618                          if (factorization_->pivots() <
    1619                                    0.5 * factorization_->maximumPivots() &&
    1620                                    factorization_->pivots() < 200)
    1621                               factorization_->areaFactor(
    1622                                    factorization_->areaFactor() * 1.1);
    1623                          problemStatus_ = -2; // factorize now
    1624                     } else if (updateStatus == 5) {
    1625                          problemStatus_ = -2; // factorize now
    1626                     }
    1627                     // update primal solution
    1628                     if (theta_ < 0.0 && candidate == -1) {
     1575            printf("flag b %g\n", alpha_);
     1576#endif
     1577            setFlagged(sequenceOut_);
     1578            progress_.clearBadTimes();
     1579            lastBadIteration_ = numberIterations_; // say be more cautious
     1580            rowArray_[0]->clear();
     1581            rowArray_[1]->clear();
     1582            columnArray_[0]->clear();
     1583            // make sure dual feasible
     1584            // look at all rows and columns
     1585            double objectiveChange = 0.0;
     1586            updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
     1587              0.0, objectiveChange, true);
     1588            rowArray_[1]->clear();
     1589            columnArray_[0]->clear();
     1590            continue;
     1591          }
     1592        } else if (updateStatus == 3) {
     1593          // out of memory
     1594          // increase space if not many iterations
     1595          if (factorization_->pivots() < 0.5 * factorization_->maximumPivots() && factorization_->pivots() < 200)
     1596            factorization_->areaFactor(
     1597              factorization_->areaFactor() * 1.1);
     1598          problemStatus_ = -2; // factorize now
     1599        } else if (updateStatus == 5) {
     1600          problemStatus_ = -2; // factorize now
     1601        }
     1602        // update primal solution
     1603        if (theta_ < 0.0 && candidate == -1) {
    16291604#ifdef CLP_DEBUG
    1630                          if (handler_->logLevel() & 32)
    1631                               printf("negative theta %g\n", theta_);
    1632 #endif
    1633                          theta_ = 0.0;
    1634                     }
    1635                     // do actual flips
    1636                     flipBounds(rowArray_[0], columnArray_[0]);
    1637                     //rowArray_[1]->expand();
    1638                     dualRowPivot_->updatePrimalSolution(rowArray_[1],
    1639                                                         movement,
    1640                                                         objectiveChange);
     1605          if (handler_->logLevel() & 32)
     1606            printf("negative theta %g\n", theta_);
     1607#endif
     1608          theta_ = 0.0;
     1609        }
     1610        // do actual flips
     1611        flipBounds(rowArray_[0], columnArray_[0]);
     1612        //rowArray_[1]->expand();
     1613        dualRowPivot_->updatePrimalSolution(rowArray_[1],
     1614          movement,
     1615          objectiveChange);
    16411616#ifdef CLP_DEBUG
    1642                     double oldobj = objectiveValue_;
    1643 #endif
    1644                     // modify dualout
    1645                     dualOut_ /= alpha_;
    1646                     dualOut_ *= -directionOut_;
    1647                     //setStatus(sequenceIn_,basic);
    1648                     dj_[sequenceIn_] = 0.0;
    1649                     double oldValue = valueIn_;
    1650                     if (directionIn_ == -1) {
    1651                          // as if from upper bound
    1652                          valueIn_ = upperIn_ + dualOut_;
    1653                     } else {
    1654                          // as if from lower bound
    1655                          valueIn_ = lowerIn_ + dualOut_;
    1656                     }
    1657                     objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
    1658                     // outgoing
    1659                     // set dj to zero unless values pass
    1660                     if (directionOut_ > 0) {
    1661                          valueOut_ = lowerOut_;
    1662                          if (candidate == -1)
    1663                               dj_[sequenceOut_] = theta_;
    1664                     } else {
    1665                          valueOut_ = upperOut_;
    1666                          if (candidate == -1)
    1667                               dj_[sequenceOut_] = -theta_;
    1668                     }
    1669                     solution_[sequenceOut_] = valueOut_;
    1670                     int whatNext = housekeeping(objectiveChange);
     1617        double oldobj = objectiveValue_;
     1618#endif
     1619        // modify dualout
     1620        dualOut_ /= alpha_;
     1621        dualOut_ *= -directionOut_;
     1622        //setStatus(sequenceIn_,basic);
     1623        dj_[sequenceIn_] = 0.0;
     1624        double oldValue = valueIn_;
     1625        if (directionIn_ == -1) {
     1626          // as if from upper bound
     1627          valueIn_ = upperIn_ + dualOut_;
     1628        } else {
     1629          // as if from lower bound
     1630          valueIn_ = lowerIn_ + dualOut_;
     1631        }
     1632        objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
     1633        // outgoing
     1634        // set dj to zero unless values pass
     1635        if (directionOut_ > 0) {
     1636          valueOut_ = lowerOut_;
     1637          if (candidate == -1)
     1638            dj_[sequenceOut_] = theta_;
     1639        } else {
     1640          valueOut_ = upperOut_;
     1641          if (candidate == -1)
     1642            dj_[sequenceOut_] = -theta_;
     1643        }
     1644        solution_[sequenceOut_] = valueOut_;
     1645        int whatNext = housekeeping(objectiveChange);
    16711646#if 0
    16721647                    for (int i=0;i<numberRows_+numberColumns_;i++) {
     
    16811656#endif
    16821657#ifdef CLP_REPORT_PROGRESS
    1683                     if (ixxxxxx > ixxyyyy - 5) {
    1684                          handler_->setLogLevel(63);
    1685                          int nTotal = numberColumns_ + numberRows_;
    1686                          double oldObj = 0.0;
    1687                          double newObj = 0.0;
    1688                          for (int i = 0; i < nTotal; i++) {
    1689                               if (savePSol[i])
    1690                                    oldObj += savePSol[i] * saveCost[i];
    1691                               if (solution_[i])
    1692                                    newObj += solution_[i] * cost_[i];
    1693                               bool printIt = false;
    1694                               if (cost_[i] != saveCost[i])
    1695                                    printIt = true;
    1696                               if (status_[i] != saveStat[i])
    1697                                    printIt = true;
    1698                               if (printIt)
    1699                                    printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
    1700                                           i, saveStat[i], saveCost[i], savePSol[i],
    1701                                           status_[i], cost_[i], solution_[i]);
    1702                               // difference
    1703                               savePSol[i] = solution_[i] - savePSol[i];
    1704                          }
    1705                          printf("pivots %d, old obj %g new %g\n",
    1706                                 factorization_->pivots(),
    1707                                 oldObj, newObj);
    1708                          memset(saveDj, 0, numberRows_ * sizeof(double));
    1709                          times(1.0, savePSol, saveDj);
    1710                          double largest = 1.0e-6;
    1711                          int k = -1;
    1712                          for (int i = 0; i < numberRows_; i++) {
    1713                               saveDj[i] -= savePSol[i+numberColumns_];
    1714                               if (fabs(saveDj[i]) > largest) {
    1715                                    largest = fabs(saveDj[i]);
    1716                                    k = i;
    1717                               }
    1718                          }
    1719                          if (k >= 0)
    1720                               printf("Not null %d %g\n", k, largest);
    1721                     }
     1658        if (ixxxxxx > ixxyyyy - 5) {
     1659          handler_->setLogLevel(63);
     1660          int nTotal = numberColumns_ + numberRows_;
     1661          double oldObj = 0.0;
     1662          double newObj = 0.0;
     1663          for (int i = 0; i < nTotal; i++) {
     1664            if (savePSol[i])
     1665              oldObj += savePSol[i] * saveCost[i];
     1666            if (solution_[i])
     1667              newObj += solution_[i] * cost_[i];
     1668            bool printIt = false;
     1669            if (cost_[i] != saveCost[i])
     1670              printIt = true;
     1671            if (status_[i] != saveStat[i])
     1672              printIt = true;
     1673            if (printIt)
     1674              printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
     1675                i, saveStat[i], saveCost[i], savePSol[i],
     1676                status_[i], cost_[i], solution_[i]);
     1677            // difference
     1678            savePSol[i] = solution_[i] - savePSol[i];
     1679          }
     1680          printf("pivots %d, old obj %g new %g\n",
     1681            factorization_->pivots(),
     1682            oldObj, newObj);
     1683          memset(saveDj, 0, numberRows_ * sizeof(double));
     1684          times(1.0, savePSol, saveDj);
     1685          double largest = 1.0e-6;
     1686          int k = -1;
     1687          for (int i = 0; i < numberRows_; i++) {
     1688            saveDj[i] -= savePSol[i + numberColumns_];
     1689            if (fabs(saveDj[i]) > largest) {
     1690              largest = fabs(saveDj[i]);
     1691              k = i;
     1692            }
     1693          }
     1694          if (k >= 0)
     1695            printf("Not null %d %g\n", k, largest);
     1696        }
    17221697#endif
    17231698#ifdef VUB
    1724                     {
    1725                          if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) || toVub[sequenceIn_] >= 0 ||
    1726                                    (sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) || toVub[sequenceOut_] >= 0) {
    1727                               int inSequence = sequenceIn_;
    1728                               int inVub = -1;
    1729                               if (sequenceIn_ < numberColumns_)
    1730                                    inVub = vub[sequenceIn_];
    1731                               int inBack = toVub[inSequence];
    1732                               int inSlack = -1;
    1733                               if (inSequence >= numberColumns_ && inBack >= 0) {
    1734                                    inSlack = inSequence - numberColumns_;
    1735                                    inSequence = inBack;
    1736                                    inBack = toVub[inSequence];
    1737                               }
    1738                               if (inVub >= 0)
    1739                                    printf("Vub %d in ", inSequence);
    1740                               if (inBack >= 0 && inSlack < 0)
    1741                                    printf("%d (descendent of %d) in ", inSequence, inBack);
    1742                               if (inSlack >= 0)
    1743                                    printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack);
    1744                               int outSequence = sequenceOut_;
    1745                               int outVub = -1;
    1746                               if (sequenceOut_ < numberColumns_)
    1747                                    outVub = vub[sequenceOut_];
    1748                               int outBack = toVub[outSequence];
    1749                               int outSlack = -1;
    1750                               if (outSequence >= numberColumns_ && outBack >= 0) {
    1751                                    outSlack = outSequence - numberColumns_;
    1752                                    outSequence = outBack;
    1753                                    outBack = toVub[outSequence];
    1754                               }
    1755                               if (outVub >= 0)
    1756                                    printf("Vub %d out ", outSequence);
    1757                               if (outBack >= 0 && outSlack < 0)
    1758                                    printf("%d (descendent of %d) out ", outSequence, outBack);
    1759                               if (outSlack >= 0)
    1760                                    printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack);
    1761                               printf("\n");
    1762                          }
    1763                     }
     1699        {
     1700          if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) || toVub[sequenceIn_] >= 0 || (sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) || toVub[sequenceOut_] >= 0) {
     1701            int inSequence = sequenceIn_;
     1702            int inVub = -1;
     1703            if (sequenceIn_ < numberColumns_)
     1704              inVub = vub[sequenceIn_];
     1705            int inBack = toVub[inSequence];
     1706            int inSlack = -1;
     1707            if (inSequence >= numberColumns_ && inBack >= 0) {
     1708              inSlack = inSequence - numberColumns_;
     1709              inSequence = inBack;
     1710              inBack = toVub[inSequence];
     1711            }
     1712            if (inVub >= 0)
     1713              printf("Vub %d in ", inSequence);
     1714            if (inBack >= 0 && inSlack < 0)
     1715              printf("%d (descendent of %d) in ", inSequence, inBack);
     1716            if (inSlack >= 0)
     1717              printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack);
     1718            int outSequence = sequenceOut_;
     1719            int outVub = -1;
     1720            if (sequenceOut_ < numberColumns_)
     1721              outVub = vub[sequenceOut_];
     1722            int outBack = toVub[outSequence];
     1723            int outSlack = -1;
     1724            if (outSequence >= numberColumns_ && outBack >= 0) {
     1725              outSlack = outSequence - numberColumns_;
     1726              outSequence = outBack;
     1727              outBack = toVub[outSequence];
     1728            }
     1729            if (outVub >= 0)
     1730              printf("Vub %d out ", outSequence);
     1731            if (outBack >= 0 && outSlack < 0)
     1732              printf("%d (descendent of %d) out ", outSequence, outBack);
     1733            if (outSlack >= 0)
     1734              printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack);
     1735            printf("\n");
     1736          }
     1737        }
    17641738#endif
    17651739#if 0
     
    17691743                         exit(77);
    17701744#endif
    1771                     if (!givenDuals && ifValuesPass && ifValuesPass != 2) {
    1772                          handler_->message(CLP_END_VALUES_PASS, messages_)
    1773                                    << numberIterations_;
    1774                          whatNext = 1;
    1775                     }
     1745        if (!givenDuals && ifValuesPass && ifValuesPass != 2) {
     1746          handler_->message(CLP_END_VALUES_PASS, messages_)
     1747            << numberIterations_;
     1748          whatNext = 1;
     1749        }
    17761750#ifdef CHECK_ACCURACY
    1777                     if (whatNext) {
    1778                          CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz);
    1779                     }
    1780 #endif
    1781                     //if (numberIterations_==1890)
    1782                     //whatNext=1;
    1783                     //if (numberIterations_>2000)
    1784                     //exit(77);
    1785                     // and set bounds correctly
    1786                     originalBound(sequenceIn_);
    1787                     changeBound(sequenceOut_);
     1751        if (whatNext) {
     1752          CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz);
     1753        }
     1754#endif
     1755        //if (numberIterations_==1890)
     1756        //whatNext=1;
     1757        //if (numberIterations_>2000)
     1758        //exit(77);
     1759        // and set bounds correctly
     1760        originalBound(sequenceIn_);
     1761        changeBound(sequenceOut_);
    17881762#ifdef CLP_DEBUG
    1789                     if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
    1790                          printf("obj backwards %g %g\n", objectiveValue_, oldobj);
     1763        if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
     1764          printf("obj backwards %g %g\n", objectiveValue_, oldobj);
    17911765#endif
    17921766#if 0
     
    18051779                    }
    18061780#endif
    1807                     if (whatNext == 1 || candidate == -2) {
    1808                          problemStatus_ = -2; // refactorize
    1809                     } else if (whatNext == 2) {
    1810                          // maximum iterations or equivalent
    1811                          problemStatus_ = 3;
    1812                          returnCode = 3;
    1813                          break;
     1781        if (whatNext == 1 || candidate == -2) {
     1782          problemStatus_ = -2; // refactorize
     1783        } else if (whatNext == 2) {
     1784          // maximum iterations or equivalent
     1785          problemStatus_ = 3;
     1786          returnCode = 3;
     1787          break;
     1788        }
     1789        // Check event
     1790        {
     1791          int status = eventHandler_->event(ClpEventHandler::endOfIteration);
     1792          if (status >= 0) {
     1793            problemStatus_ = 5;
     1794            secondaryStatus_ = ClpEventHandler::endOfIteration;
     1795            returnCode = 4;
     1796            break;
     1797          }
     1798        }
     1799      } else {
     1800#ifdef CLP_INVESTIGATE_SERIAL
     1801        z_thinks = 1;
     1802#endif
     1803        // no incoming column is valid
     1804        spareIntArray_[3] = pivotRow_;
     1805        pivotRow_ = -1;
     1806#ifdef CLP_DEBUG
     1807        if (handler_->logLevel() & 32)
     1808          printf("** no column pivot\n");
     1809#endif
     1810        delete[] ray_;
     1811        ray_ = NULL;
     1812        if ((factorization_->pivots() < 2
     1813              || ((specialOptions_ & 2097152) != 0 && factorization_->pivots() < 50))
     1814          && acceptablePivot_ <= 1.0e-8 && acceptablePivot_ > 0.0) {
     1815          //&&goodAccuracy()) {
     1816          // If not in branch and bound etc save ray
     1817          if ((specialOptions_ & (1024 | 4096)) == 0 || (specialOptions_ & (32 | 2097152)) != 0) {
     1818            // create ray anyway
     1819            ray_ = new double[numberRows_];
     1820            rowArray_[0]->expand(); // in case packed
     1821            const double *array = rowArray_[0]->denseVector();
     1822            for (int i = 0; i < numberRows_; i++)
     1823              ray_[i] = array[i];
     1824#ifdef PRINT_RAY_METHOD
     1825            {
     1826              double *farkas = new double[2 * numberColumns_ + numberRows_];
     1827              int nBasic = 0;
     1828              int nPlusLower = 0;
     1829              int nPlusFixedLower = 0;
     1830              int nMinusLower = 0;
     1831              int nMinusFixedLower = 0;
     1832              int nPlusUpper = 0;
     1833              int nPlusFixedUpper = 0;
     1834              int nMinusUpper = 0;
     1835              int nMinusFixedUpper = 0;
     1836              memset(farkas, 0, (2 * numberColumns_ + numberRows_) * sizeof(double));
     1837              transposeTimes(-1.0, ray_, farkas);
     1838              for (int i = 0; i < numberRows_; i++) {
     1839                if (fabs(ray_[i]) > 1.0e-7) {
     1840                  if (getRowStatus(i) == basic) {
     1841                    nBasic++;
     1842                  } else if (getRowStatus(i) == atLowerBound) {
     1843                    if (ray_[i] > 0.0)
     1844                      nPlusLower++;
     1845                    else
     1846                      nMinusLower++;
     1847                  } else if (getRowStatus(i) == atUpperBound) {
     1848                    if (ray_[i] > 0.0)
     1849                      nPlusUpper++;
     1850                    else
     1851                      nMinusUpper++;
     1852                  } else {
     1853                    // fixed slack
     1854                  }
     1855                }
     1856              }
     1857              printf("Slacks %d basic lower +,- %d,%d upper +,- %d,%d\n",
     1858                nBasic, nPlusLower, nMinusLower, nPlusUpper, nMinusLower);
     1859              for (int i = 0; i < numberColumns_; i++) {
     1860                if (fabs(farkas[i]) > 1.0e-7) {
     1861                  if (getColumnStatus(i) == basic) {
     1862                    nBasic++;
     1863                  } else if (getColumnStatus(i) == atLowerBound) {
     1864                    if (farkas[i] > 0.0)
     1865                      nPlusLower++;
     1866                    else
     1867                      nMinusLower++;
     1868                  } else if (getColumnStatus(i) == atUpperBound) {
     1869                    if (farkas[i] > 0.0)
     1870                      nPlusUpper++;
     1871                    else
     1872                      nMinusUpper++;
     1873                  } else {
     1874                    if (!lower_[i]) {
     1875                      if (farkas[i] > 0.0) {
     1876                        nPlusFixedLower++;
     1877                      } else {
     1878                        nMinusFixedLower++;
     1879                      }
     1880                    } else {
     1881                      if (farkas[i] > 0.0) {
     1882                        nPlusFixedUpper++;
     1883                      } else {
     1884                        nMinusFixedUpper++;
     1885                      }
    18141886                    }
    1815                     // Check event
    1816                     {
    1817                          int status = eventHandler_->event(ClpEventHandler::endOfIteration);
    1818                          if (status >= 0) {
    1819                               problemStatus_ = 5;
    1820                               secondaryStatus_ = ClpEventHandler::endOfIteration;
    1821                               returnCode = 4;
    1822                               break;
    1823                          }
    1824                     }
    1825                } else {
     1887                  }
     1888                }
     1889              }
     1890              printf("End %d basic lower +,- %d,%d upper +,- %d,%d fixed %d,%d %d,%d\n",
     1891                nBasic, nPlusLower, nMinusLower, nPlusUpper, nMinusUpper,
     1892                nPlusFixedLower, nMinusFixedLower, nPlusFixedUpper, nMinusFixedUpper);
     1893              printf("Dual creating infeasibility ray direction out %d - pivRow %d seqOut %d lower %g,val %g,upper %g\n",
     1894                directionOut_, spareIntArray_[3], sequenceOut_, lowerOut_, valueOut_, upperOut_);
     1895              delete[] farkas;
     1896            }
     1897#endif
     1898          } else {
     1899            ray_ = NULL;
     1900          }
     1901          // If we have just factorized and infeasibility reasonable say infeas
     1902          double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13;
     1903          // but if none at fake bounds
     1904          if (!checkFakeBounds())
     1905            dualTest = 0.0;
     1906          if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) {
     1907            double testValue = 1.0e-4;
     1908            if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1)
     1909              testValue = 1.0e-6;
     1910            if (valueOut_ > upperOut_ + testValue || valueOut_ < lowerOut_ - testValue
     1911              || (specialOptions_ & 64) == 0) {
     1912              // say infeasible
     1913              problemStatus_ = 1;
     1914              // unless primal feasible!!!!
     1915              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
     1916              //   numberDualInfeasibilities_,sumDualInfeasibilities_);
     1917              //#define TEST_CLP_NODE
     1918#ifndef TEST_CLP_NODE
     1919              // Should be correct - but ...
     1920              int numberFake = numberAtFakeBound();
     1921              double sumPrimal = (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
     1922              if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-5 || (sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) {
     1923                if (sumPrimal > 50.0 && factorization_->pivots() > 2) {
     1924                  problemStatus_ = -4;
     1925#ifdef COIN_DEVELOP
     1926                  printf("status to -4 at %d - primalinf %g pivots %d\n",
     1927                    __LINE__, sumPrimalInfeasibilities_,
     1928                    factorization_->pivots());
     1929#endif
     1930                } else {
     1931                  problemStatus_ = 10;
     1932#if COIN_DEVELOP > 1
     1933                  printf("returning at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d - options (1024-16384) %d %d %d %d %d\n",
     1934                    __LINE__, numberPrimalInfeasibilities_,
     1935                    sumPrimalInfeasibilities_,
     1936                    numberDualInfeasibilities_, sumDualInfeasibilities_,
     1937                    numberFake_, dualBound_, factorization_->pivots(),
     1938                    (specialOptions_ & 1024) != 0 ? 1 : 0,
     1939                    (specialOptions_ & 2048) != 0 ? 1 : 0,
     1940                    (specialOptions_ & 4096) != 0 ? 1 : 0,
     1941                    (specialOptions_ & 8192) != 0 ? 1 : 0,
     1942                    (specialOptions_ & 16384) != 0 ? 1 : 0);
     1943#endif
     1944                  // Get rid of objective
     1945                  if ((specialOptions_ & 16384) == 0)
     1946                    objective_ = new ClpLinearObjective(NULL, numberColumns_);
     1947                }
     1948              }
     1949#else
     1950              if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) {
     1951#ifdef COIN_DEVELOP
     1952                printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
     1953                  __LINE__, numberPrimalInfeasibilities_,
     1954                  sumPrimalInfeasibilities_,
     1955                  numberDualInfeasibilities_, sumDualInfeasibilities_,
     1956                  numberFake_, dualBound_, factorization_->pivots());
     1957#endif
     1958                if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
     1959                  problemStatus_ = 10;
     1960#if COIN_DEVELOP > 1
     1961                  printf("returning at %d\n", __LINE__);
     1962#endif
     1963                  // Get rid of objective
     1964                  if ((specialOptions_ & 16384) == 0)
     1965                    objective_ = new ClpLinearObjective(NULL, numberColumns_);
     1966                }
     1967              }
     1968#endif
     1969              rowArray_[0]->clear();
     1970              columnArray_[0]->clear();
     1971              returnCode = 1;
     1972              break;
     1973            }
     1974          }
     1975          // If special option set - put off as long as possible
     1976          if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) {
     1977            if (factorization_->pivots() == 0)
     1978              problemStatus_ = -4; //say looks infeasible
     1979          } else {
     1980            // flag
     1981            char x = isColumn(sequenceOut_) ? 'C' : 'R';
     1982            handler_->message(CLP_SIMPLEX_FLAG, messages_)
     1983              << x << sequenceWithin(sequenceOut_)
     1984              << CoinMessageEol;
     1985#ifdef COIN_DEVELOP
     1986            printf("flag c\n");
     1987#endif
     1988            setFlagged(sequenceOut_);
     1989            if (!factorization_->pivots()) {
     1990              rowArray_[0]->clear();
     1991              columnArray_[0]->clear();
     1992              continue;
     1993            }
     1994          }
     1995        }
     1996        acceptablePivot_ = fabs(acceptablePivot_);
     1997        if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
     1998          acceptablePivot_ = 1.0e-8;
     1999        rowArray_[0]->clear();
     2000        columnArray_[0]->clear();
     2001        returnCode = 1;
     2002        break;
     2003      }
     2004    } else {
    18262005#ifdef CLP_INVESTIGATE_SERIAL
    1827                     z_thinks = 1;
    1828 #endif
    1829                     // no incoming column is valid
    1830                     spareIntArray_[3]=pivotRow_;
    1831                     pivotRow_ = -1;
     2006      z_thinks = 0;
     2007#endif
     2008      // no pivot row
    18322009#ifdef CLP_DEBUG
    1833                     if (handler_->logLevel() & 32)
    1834                          printf("** no column pivot\n");
    1835 #endif
    1836                     delete [] ray_;
    1837                     ray_ = NULL;
    1838                     if ((factorization_->pivots() < 2
    1839                          ||((specialOptions_&2097152)!=0&&factorization_->pivots()<50))
    1840                         && acceptablePivot_ <= 1.0e-8 && acceptablePivot_ > 0.0) {
    1841                          //&&goodAccuracy()) {
    1842                          // If not in branch and bound etc save ray
    1843                          if ((specialOptions_&(1024 | 4096)) == 0 || (specialOptions_ & (32|2097152)) != 0) {
    1844                               // create ray anyway
    1845                               ray_ = new double [ numberRows_];
    1846                               rowArray_[0]->expand(); // in case packed
    1847                               const double * array = rowArray_[0]->denseVector();
    1848                               for (int i=0;i<numberRows_;i++)
    1849                                 ray_[i] = array[i];
    1850 #ifdef PRINT_RAY_METHOD
    1851                               {
    1852                                 double * farkas = new double [2*numberColumns_+numberRows_];
    1853                                 int nBasic=0;
    1854                                 int nPlusLower=0;
    1855                                 int nPlusFixedLower=0;
    1856                                 int nMinusLower=0;
    1857                                 int nMinusFixedLower=0;
    1858                                 int nPlusUpper=0;
    1859                                 int nPlusFixedUpper=0;
    1860                                 int nMinusUpper=0;
    1861                                 int nMinusFixedUpper=0;
    1862                                 memset(farkas,0,(2*numberColumns_+numberRows_)*sizeof(double));
    1863                                 transposeTimes(-1.0,ray_,farkas);
    1864                                 for (int i=0;i<numberRows_;i++) {
    1865                                   if (fabs(ray_[i])>1.0e-7) {
    1866                                     if (getRowStatus(i)==basic) {
    1867                                       nBasic++;
    1868                                     } else if (getRowStatus(i)==atLowerBound) {
    1869                                       if (ray_[i]>0.0)
    1870                                         nPlusLower++;
    1871                                       else
    1872                                         nMinusLower++;
    1873                                     } else if (getRowStatus(i)==atUpperBound) {
    1874                                       if (ray_[i]>0.0)
    1875                                         nPlusUpper++;
    1876                                       else
    1877                                         nMinusUpper++;
    1878                                     } else {
    1879                                       // fixed slack
    1880                                     }
    1881                                   }
    1882                                 }
    1883                                 printf("Slacks %d basic lower +,- %d,%d upper +,- %d,%d\n",
    1884                                        nBasic,nPlusLower,nMinusLower,nPlusUpper,nMinusLower);
    1885                                 for (int i=0;i<numberColumns_;i++) {
    1886                                   if (fabs(farkas[i])>1.0e-7) {
    1887                                     if (getColumnStatus(i)==basic) {
    1888                                       nBasic++;
    1889                                     } else if (getColumnStatus(i)==atLowerBound) {
    1890                                       if (farkas[i]>0.0)
    1891                                         nPlusLower++;
    1892                                       else
    1893                                         nMinusLower++;
    1894                                     } else if (getColumnStatus(i)==atUpperBound) {
    1895                                       if (farkas[i]>0.0)
    1896                                         nPlusUpper++;
    1897                                       else
    1898                                         nMinusUpper++;
    1899                                     } else {
    1900                                       if (!lower_[i]) {
    1901                                         if (farkas[i]>0.0) {
    1902                                           nPlusFixedLower++;
    1903                                         } else {
    1904                                           nMinusFixedLower++;
    1905                                         }
    1906                                       } else {
    1907                                         if (farkas[i]>0.0) {
    1908                                           nPlusFixedUpper++;
    1909                                         } else {
    1910                                           nMinusFixedUpper++;
    1911                                         }
    1912                                       }
    1913                                     }
    1914                                   }
    1915                                 }
    1916                                 printf("End %d basic lower +,- %d,%d upper +,- %d,%d fixed %d,%d %d,%d\n",
    1917                                        nBasic,nPlusLower,nMinusLower,nPlusUpper,nMinusUpper,
    1918                                        nPlusFixedLower,nMinusFixedLower,nPlusFixedUpper,nMinusFixedUpper);
    1919                                 printf("Dual creating infeasibility ray direction out %d - pivRow %d seqOut %d lower %g,val %g,upper %g\n",
    1920                                        directionOut_,spareIntArray_[3],sequenceOut_,lowerOut_,valueOut_,upperOut_);
    1921                                 delete [] farkas;
    1922                               }
    1923 #endif
    1924                          } else {
    1925                               ray_ = NULL;
    1926                          }
    1927                          // If we have just factorized and infeasibility reasonable say infeas
    1928                          double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13;
    1929                          // but if none at fake bounds
    1930                          if (!checkFakeBounds())
    1931                            dualTest=0.0;
    1932                          if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) {
    1933                               double testValue = 1.0e-4;
    1934                               if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1)
    1935                                    testValue = 1.0e-6;
    1936                               if (valueOut_ > upperOut_ + testValue || valueOut_ < lowerOut_ - testValue
    1937                                         || (specialOptions_ & 64) == 0) {
    1938                                    // say infeasible
    1939                                    problemStatus_ = 1;
    1940                                    // unless primal feasible!!!!
    1941                                    //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
    1942                                    //   numberDualInfeasibilities_,sumDualInfeasibilities_);
    1943                                    //#define TEST_CLP_NODE
    1944 #ifndef TEST_CLP_NODE
    1945                                    // Should be correct - but ...
    1946                                    int numberFake = numberAtFakeBound();
    1947                                    double sumPrimal =  (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
    1948                                    if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-5 ||
    1949                                              (sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) {
    1950                                         if (sumPrimal > 50.0 && factorization_->pivots() > 2) {
    1951                                              problemStatus_ = -4;
     2010      if (handler_->logLevel() & 32)
     2011        printf("** no row pivot\n");
     2012#endif
     2013      // If in branch and bound try and get rid of fixed variables
     2014      if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
     2015        assert(!candidateList);
     2016        candidateList = new int[numberRows_];
     2017        int iRow;
     2018        for (iRow = 0; iRow < numberRows_; iRow++) {
     2019          int iPivot = pivotVariable_[iRow];
     2020          if (flagged(iPivot) || !pivoted(iPivot))
     2021            continue;
     2022          assert(iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
     2023          candidateList[numberCandidates++] = iRow;
     2024        }
     2025        // and set first candidate
     2026        if (!numberCandidates) {
     2027          delete[] candidateList;
     2028          candidateList = NULL;
     2029          int iRow;
     2030          for (iRow = 0; iRow < numberRows_; iRow++) {
     2031            int iPivot = pivotVariable_[iRow];
     2032            clearPivoted(iPivot);
     2033          }
     2034        } else {
     2035          ifValuesPass = 2;
     2036          continue;
     2037        }
     2038      }
     2039      int numberPivots = factorization_->pivots();
     2040      bool specialCase;
     2041      int useNumberFake;
     2042      returnCode = 0;
     2043      if (numberPivots <= CoinMax(dontFactorizePivots_, 20) && (specialOptions_ & 2048) != 0 && (true || !numberChanged_ || perturbation_ == 101)
     2044        && dualBound_ >= 1.0e8) {
     2045        specialCase = true;
     2046        // as dual bound high - should be okay
     2047        useNumberFake = 0;
     2048      } else {
     2049        specialCase = false;
     2050        useNumberFake = numberFake_;
     2051      }
     2052      if (!numberPivots || specialCase) {
     2053        if (numberPrimalInfeasibilities_ && problemStatus_ == -1)
     2054          problemStatus_ = -4;
     2055        // may have crept through - so may be optimal
     2056        // check any flagged variables
     2057        int iRow;
     2058        for (iRow = 0; iRow < numberRows_; iRow++) {
     2059          int iPivot = pivotVariable_[iRow];
     2060          if (flagged(iPivot))
     2061            break;
     2062        }
     2063        if (iRow < numberRows_ && numberPivots) {
     2064          // try factorization
     2065          returnCode = -2;
     2066        }
     2067
     2068        if (useNumberFake || numberDualInfeasibilities_) {
     2069          // may be dual infeasible
     2070          if ((specialOptions_ & 1024) == 0)
     2071            problemStatus_ = -5;
     2072          else if (!useNumberFake && numberPrimalInfeasibilities_
     2073            && !numberPivots)
     2074            problemStatus_ = 1;
     2075        } else {
     2076          if (iRow < numberRows_) {
    19522077#ifdef COIN_DEVELOP
    1953                                              printf("status to -4 at %d - primalinf %g pivots %d\n",
    1954                                                     __LINE__, sumPrimalInfeasibilities_,
    1955                                                     factorization_->pivots());
    1956 #endif
    1957                                         } else {
    1958                                              problemStatus_ = 10;
    1959 #if COIN_DEVELOP>1
    1960                                              printf("returning at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d - options (1024-16384) %d %d %d %d %d\n",
    1961                                                     __LINE__, numberPrimalInfeasibilities_,
    1962                                                     sumPrimalInfeasibilities_,
    1963                                                     numberDualInfeasibilities_, sumDualInfeasibilities_,
    1964                                                     numberFake_, dualBound_, factorization_->pivots(),
    1965                                                     (specialOptions_ & 1024) != 0 ? 1 : 0,
    1966                                                     (specialOptions_ & 2048) != 0 ? 1 : 0,
    1967                                                     (specialOptions_ & 4096) != 0 ? 1 : 0,
    1968                                                     (specialOptions_ & 8192) != 0 ? 1 : 0,
    1969                                                     (specialOptions_ & 16384) != 0 ? 1 : 0
    1970                                                    );
    1971 #endif
    1972                                              // Get rid of objective
    1973                                              if ((specialOptions_ & 16384) == 0)
    1974                                                   objective_ = new ClpLinearObjective(NULL, numberColumns_);
    1975                                         }
    1976                                    }
    1977 #else
    1978                                    if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) {
    1979 #ifdef COIN_DEVELOP
    1980                                         printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
    1981                                                __LINE__, numberPrimalInfeasibilities_,
    1982                                                sumPrimalInfeasibilities_,
    1983                                                numberDualInfeasibilities_, sumDualInfeasibilities_,
    1984                                                numberFake_, dualBound_, factorization_->pivots());
    1985 #endif
    1986                                         if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
    1987                                              problemStatus_ = 10;
    1988 #if COIN_DEVELOP>1
    1989                                              printf("returning at %d\n", __LINE__);
    1990 #endif
    1991                                              // Get rid of objective
    1992                                              if ((specialOptions_ & 16384) == 0)
    1993                                                   objective_ = new ClpLinearObjective(NULL, numberColumns_);
    1994                                         }
    1995                                    }
    1996 #endif
    1997                                    rowArray_[0]->clear();
    1998                                    columnArray_[0]->clear();
    1999                                    returnCode = 1;
    2000                                    break;
    2001                               }
    2002                          }
    2003                          // If special option set - put off as long as possible
    2004                          if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) {
    2005                               if (factorization_->pivots() == 0)
    2006                                    problemStatus_ = -4; //say looks infeasible
    2007                          } else {
    2008                               // flag
    2009                               char x = isColumn(sequenceOut_) ? 'C' : 'R';
    2010                               handler_->message(CLP_SIMPLEX_FLAG, messages_)
    2011                                         << x << sequenceWithin(sequenceOut_)
    2012                                         << CoinMessageEol;
    2013 #ifdef COIN_DEVELOP
    2014                               printf("flag c\n");
    2015 #endif
    2016                               setFlagged(sequenceOut_);
    2017                               if (!factorization_->pivots()) {
    2018                                    rowArray_[0]->clear();
    2019                                    columnArray_[0]->clear();
    2020                                    continue;
    2021                               }
    2022                          }
    2023                     }
    2024                     acceptablePivot_ = fabs(acceptablePivot_);
    2025                     if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
    2026                          acceptablePivot_ = 1.0e-8;
    2027                     rowArray_[0]->clear();
    2028                     columnArray_[0]->clear();
    2029                     returnCode = 1;
    2030                     break;
    2031                }
     2078            std::cout << "Flagged variables at end - infeasible?" << std::endl;
     2079            printf("Probably infeasible - pivot was %g\n", alpha_);
     2080#endif
     2081            //if (fabs(alpha_)<1.0e-4) {
     2082            //problemStatus_=1;
     2083            //} else {
     2084#ifdef CLP_DEBUG
     2085            abort();
     2086#endif
     2087            //}
     2088            problemStatus_ = -5;
    20322089          } else {
    2033 #ifdef CLP_INVESTIGATE_SERIAL
    2034                z_thinks = 0;
    2035 #endif
    2036                // no pivot row
    2037 #ifdef CLP_DEBUG
    2038                if (handler_->logLevel() & 32)
    2039                     printf("** no row pivot\n");
    2040 #endif
    2041                // If in branch and bound try and get rid of fixed variables
    2042                if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
    2043                     assert (!candidateList);
    2044                     candidateList = new int[numberRows_];
    2045                     int iRow;
    2046                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2047                          int iPivot = pivotVariable_[iRow];
    2048                          if (flagged(iPivot) || !pivoted(iPivot))
    2049                               continue;
    2050                          assert (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
    2051                          candidateList[numberCandidates++] = iRow;
    2052                     }
    2053                     // and set first candidate
    2054                     if (!numberCandidates) {
    2055                          delete [] candidateList;
    2056                          candidateList = NULL;
    2057                          int iRow;
    2058                          for (iRow = 0; iRow < numberRows_; iRow++) {
    2059                               int iPivot = pivotVariable_[iRow];
    2060                               clearPivoted(iPivot);
    2061                          }
    2062                     } else {
    2063                          ifValuesPass = 2;
    2064                          continue;
    2065                     }
    2066                }
    2067                int numberPivots = factorization_->pivots();
    2068                bool specialCase;
    2069                int useNumberFake;
    2070                returnCode = 0;
    2071                if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
    2072                          (specialOptions_ & 2048) != 0 && (true || !numberChanged_ || perturbation_ == 101)
    2073                          && dualBound_ >= 1.0e8) {
    2074                     specialCase = true;
    2075                     // as dual bound high - should be okay
    2076                     useNumberFake = 0;
    2077                } else {
    2078                     specialCase = false;
    2079                     useNumberFake = numberFake_;
    2080                }
    2081                if (!numberPivots || specialCase) {
    2082                  if (numberPrimalInfeasibilities_ && problemStatus_==-1)
    2083                    problemStatus_=-4;
    2084                     // may have crept through - so may be optimal
    2085                     // check any flagged variables
    2086                     int iRow;
    2087                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2088                          int iPivot = pivotVariable_[iRow];
    2089                          if (flagged(iPivot))
    2090                               break;
    2091                     }
    2092                     if (iRow < numberRows_ && numberPivots) {
    2093                          // try factorization
    2094                          returnCode = -2;
    2095                     }
    2096 
    2097                     if (useNumberFake || numberDualInfeasibilities_) {
    2098                          // may be dual infeasible
    2099                          if ((specialOptions_ & 1024) == 0)
    2100                               problemStatus_ = -5;
    2101                          else if (!useNumberFake && numberPrimalInfeasibilities_
    2102                                    && !numberPivots)
    2103                               problemStatus_ = 1;
    2104                     } else {
    2105                          if (iRow < numberRows_) {
    2106 #ifdef COIN_DEVELOP
    2107                               std::cout << "Flagged variables at end - infeasible?" << std::endl;
    2108                               printf("Probably infeasible - pivot was %g\n", alpha_);
    2109 #endif
    2110                               //if (fabs(alpha_)<1.0e-4) {
    2111                               //problemStatus_=1;
    2112                               //} else {
    2113 #ifdef CLP_DEBUG
    2114                               abort();
    2115 #endif
    2116                               //}
    2117                               problemStatus_ = -5;
    2118                          } else {
    2119                               problemStatus_ = 0;
     2090            problemStatus_ = 0;
    21202091#ifndef CLP_CHECK_NUMBER_PIVOTS
    21212092#define CLP_CHECK_NUMBER_PIVOTS 10
    21222093#endif
    21232094#if CLP_CHECK_NUMBER_PIVOTS < 20
    2124                               if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
     2095            if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
    21252096#ifndef NDEBUG_CLP
    2126                                    int nTotal = numberRows_ + numberColumns_;
    2127                                    double * comp = CoinCopyOfArray(solution_, nTotal);
    2128 #endif
    2129                                    computePrimals(rowActivityWork_, columnActivityWork_);
     2097              int nTotal = numberRows_ + numberColumns_;
     2098              double *comp = CoinCopyOfArray(solution_, nTotal);
     2099#endif
     2100              computePrimals(rowActivityWork_, columnActivityWork_);
    21302101#ifndef NDEBUG_CLP
    2131                                    double largest = 1.0e-5;
    2132                                    int bad = -1;
    2133                                    for (int i = 0; i < nTotal; i++) {
    2134                                         double value = solution_[i];
    2135                                         double larger = CoinMax(fabs(value), fabs(comp[i]));
    2136                                         double tol = 1.0e-5 + 1.0e-5 * larger;
    2137                                         double diff = fabs(value - comp[i]);
    2138                                         if (diff - tol > largest) {
    2139                                              bad = i;
    2140                                              largest = diff - tol;
    2141                                         }
    2142                                    }
    2143                                    if (bad >= 0)
    2144                                      COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad]));
    2145 #endif
    2146                                    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
    2147                                    if (numberPrimalInfeasibilities_) {
     2102              double largest = 1.0e-5;
     2103              int bad = -1;
     2104              for (int i = 0; i < nTotal; i++) {
     2105                double value = solution_[i];
     2106                double larger = CoinMax(fabs(value), fabs(comp[i]));
     2107                double tol = 1.0e-5 + 1.0e-5 * larger;
     2108                double diff = fabs(value - comp[i]);
     2109                if (diff - tol > largest) {
     2110                  bad = i;
     2111                  largest = diff - tol;
     2112                }
     2113              }
     2114              if (bad >= 0)
     2115                COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad]));
     2116#endif
     2117              checkPrimalSolution(rowActivityWork_, columnActivityWork_);
     2118              if (numberPrimalInfeasibilities_) {
    21482119#ifdef CLP_INVESTIGATE
    2149                                         printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
    2150                                                sumPrimalInfeasibilities_);
    2151 #endif
    2152                                         problemStatus_ = -1;
    2153                                         returnCode = -2;
    2154                                    }
     2120                printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
     2121                  sumPrimalInfeasibilities_);
     2122#endif
     2123                problemStatus_ = -1;
     2124                returnCode = -2;
     2125              }
    21552126#ifndef NDEBUG_CLP
    2156                                    memcpy(solution_, comp, nTotal * sizeof(double));
    2157                                    delete [] comp;
    2158 #endif
    2159                               }
    2160 #endif
    2161                               if (!problemStatus_) {
    2162                                    // make it look OK
    2163                                    numberPrimalInfeasibilities_ = 0;
    2164                                    sumPrimalInfeasibilities_ = 0.0;
    2165                                    numberDualInfeasibilities_ = 0;
    2166                                    sumDualInfeasibilities_ = 0.0;
    2167                                    // May be perturbed
    2168                                    if (perturbation_ == 101 || numberChanged_) {
    2169                                         numberChanged_ = 0; // Number of variables with changed costs
    2170                                         perturbation_ = 102; // stop any perturbations
    2171                                         //double changeCost;
    2172                                         //changeBounds(1,NULL,changeCost);
    2173                                         createRim4(false);
    2174                                         // make sure duals are current
    2175                                         computeDuals(givenDuals);
    2176                                         checkDualSolution();
    2177                                         progress_.modifyObjective(-COIN_DBL_MAX);
    2178                                         if (numberDualInfeasibilities_) {
    2179                                              problemStatus_ = 10; // was -3;
    2180                                         } else {
    2181                                              computeObjectiveValue(true);
    2182                                         }
    2183                                    } else if (numberPivots) {
    2184                                         computeObjectiveValue(true);
    2185                                    }
    2186                                    if (numberPivots < -1000) {
    2187                                         // objective may be wrong
    2188                                         objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
    2189                                         objectiveValue_ += objective_->nonlinearOffset();
    2190                                         objectiveValue_ /= (objectiveScale_ * rhsScale_);
    2191                                         if ((specialOptions_ & 16384) == 0) {
    2192                                              // and dual_ may be wrong (i.e. for fixed or basic)
    2193                                              CoinIndexedVector * arrayVector = rowArray_[1];
    2194                                              arrayVector->clear();
    2195                                              int iRow;
    2196                                              double * array = arrayVector->denseVector();
    2197                                              /* Use dual_ instead of array
     2127              memcpy(solution_, comp, nTotal * sizeof(double));
     2128              delete[] comp;
     2129#endif
     2130            }
     2131#endif
     2132            if (!problemStatus_) {
     2133              // make it look OK
     2134              numberPrimalInfeasibilities_ = 0;
     2135              sumPrimalInfeasibilities_ = 0.0;
     2136              numberDualInfeasibilities_ = 0;
     2137              sumDualInfeasibilities_ = 0.0;
     2138              // May be perturbed
     2139              if (perturbation_ == 101 || numberChanged_) {
     2140                numberChanged_ = 0; // Number of variables with changed costs
     2141                perturbation_ = 102; // stop any perturbations
     2142                //double changeCost;
     2143                //changeBounds(1,NULL,changeCost);
     2144                createRim4(false);
     2145                // make sure duals are current
     2146                computeDuals(givenDuals);
     2147                checkDualSolution();
     2148                progress_.modifyObjective(-COIN_DBL_MAX);
     2149                if (numberDualInfeasibilities_) {
     2150                  problemStatus_ = 10; // was -3;
     2151                } else {
     2152                  computeObjectiveValue(true);
     2153                }
     2154              } else if (numberPivots) {
     2155                computeObjectiveValue(true);
     2156              }
     2157              if (numberPivots < -1000) {
     2158                // objective may be wrong
     2159                objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
     2160                objectiveValue_ += objective_->nonlinearOffset();
     2161                objectiveValue_ /= (objectiveScale_ * rhsScale_);
     2162                if ((specialOptions_ & 16384) == 0) {
     2163                  // and dual_ may be wrong (i.e. for fixed or basic)
     2164                  CoinIndexedVector *arrayVector = rowArray_[1];
     2165                  arrayVector->clear();
     2166                  int iRow;
     2167                  double *array = arrayVector->denseVector();
     2168                  /* Use dual_ instead of array
    21982169                                                Even though dual_ is only numberRows_ long this is
    21992170                                                okay as gets permuted to longer rowArray_[2]
    22002171                                             */
    2201                                              arrayVector->setDenseVector(dual_);
    2202                                              int * index = arrayVector->getIndices();
    2203                                              int number = 0;
    2204                                              for (iRow = 0; iRow < numberRows_; iRow++) {
    2205                                                   int iPivot = pivotVariable_[iRow];
    2206                                                   double value = cost_[iPivot];
    2207                                                   dual_[iRow] = value;
    2208                                                   if (value) {
    2209                                                        index[number++] = iRow;
    2210                                                   }
    2211                                              }
    2212                                              arrayVector->setNumElements(number);
    2213                                              // Extended duals before "updateTranspose"
    2214                                              matrix_->dualExpanded(this, arrayVector, NULL, 0);
    2215                                              // Btran basic costs
    2216                                              rowArray_[2]->clear();
    2217                                              factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
    2218                                              // and return vector
    2219                                              arrayVector->setDenseVector(array);
    2220                                         }
    2221                                    }
    2222                                    sumPrimalInfeasibilities_ = 0.0;
    2223                               }
    2224                               if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) {
    2225                                    CoinIndexedVector * arrayVector = rowArray_[1];
    2226                                    arrayVector->clear();
    2227                                    double * rhs = arrayVector->denseVector();
    2228                                    times(1.0, solution_, rhs);
     2172                  arrayVector->setDenseVector(dual_);
     2173                  int *index = arrayVector->getIndices();
     2174                  int number = 0;
     2175                  for (iRow = 0; iRow < numberRows_; iRow++) {
     2176                    int iPivot = pivotVariable_[iRow];
     2177                    double value = cost_[iPivot];
     2178                    dual_[iRow] = value;
     2179                    if (value) {
     2180                      index[number++] = iRow;
     2181                    }
     2182                  }
     2183                  arrayVector->setNumElements(number);
     2184                  // Extended duals before "updateTranspose"
     2185                  matrix_->dualExpanded(this, arrayVector, NULL, 0);
     2186                  // Btran basic costs
     2187                  rowArray_[2]->clear();
     2188                  factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
     2189                  // and return vector
     2190                  arrayVector->setDenseVector(array);
     2191                }
     2192              }
     2193              sumPrimalInfeasibilities_ = 0.0;
     2194            }
     2195            if ((specialOptions_ & (1024 + 16384)) != 0 && !problemStatus_) {
     2196              CoinIndexedVector *arrayVector = rowArray_[1];
     2197              arrayVector->clear();
     2198              double *rhs = arrayVector->denseVector();
     2199              times(1.0, solution_, rhs);
    22292200#ifdef CHECK_ACCURACY
    2230                                    bool bad = false;
    2231 #endif
    2232                                    bool bad2 = false;
    2233                                    int i;
    2234                                    for ( i = 0; i < numberRows_; i++) {
    2235                                         if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
    2236                                                   rhs[i] > rowUpperWork_[i] + primalTolerance_) {
    2237                                              bad2 = true;
     2201              bool bad = false;
     2202#endif
     2203              bool bad2 = false;
     2204              int i;
     2205              for (i = 0; i < numberRows_; i++) {
     2206                if (rhs[i] < rowLowerWork_[i] - primalTolerance_ || rhs[i] > rowUpperWork_[i] + primalTolerance_) {
     2207                  bad2 = true;
    22382208#ifdef CHECK_ACCURACY
    2239                                              printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
    2240                                                     rowLowerWork_[i], rowUpperWork_[i],
    2241                                                     rhs[i], rowActivityWork_[i]);
    2242 #endif
    2243                                         } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
     2209                  printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
     2210                    rowLowerWork_[i], rowUpperWork_[i],
     2211                    rhs[i], rowActivityWork_[i]);
     2212#endif
     2213                } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
    22442214#ifdef CHECK_ACCURACY
    2245                                              bad = true;
    2246                                              printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
    2247 #endif
    2248                                         }
    2249                                         rhs[i] = 0.0;
    2250                                    }
    2251                                    for ( i = 0; i < numberColumns_; i++) {
    2252                                         if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
    2253                                                   solution_[i] > columnUpperWork_[i] + primalTolerance_) {
    2254                                              bad2 = true;
     2215                  bad = true;
     2216                  printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
     2217#endif
     2218                }
     2219                rhs[i] = 0.0;
     2220              }
     2221              for (i = 0; i < numberColumns_; i++) {
     2222                if (solution_[i] < columnLowerWork_[i] - primalTolerance_ || solution_[i] > columnUpperWork_[i] + primalTolerance_) {
     2223                  bad2 = true;
    22552224#ifdef CHECK_ACCURACY
    2256                                              printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
    2257                                                     columnLowerWork_[i], columnUpperWork_[i],
    2258                                                     solution_[i], columnActivityWork_[i]);
    2259 #endif
    2260                                         }
    2261                                    }
    2262                                    if (bad2) {
    2263                                         problemStatus_ = -3;
    2264                                         returnCode = -2;
    2265                                         // Force to re-factorize early next time
    2266                                         int numberPivots = factorization_->pivots();
    2267                                         forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
    2268                                    }
    2269                               }
    2270                          }
    2271                     }
    2272                } else {
    2273                     problemStatus_ = -3;
    2274                     returnCode = -2;
    2275                     // Force to re-factorize early next time
    2276                     int numberPivots = factorization_->pivots();
    2277                     forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
    2278                }
    2279                break;
    2280           }
    2281      }
    2282      if (givenDuals) {
    2283           CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
    2284           // get rid of any values pass array
    2285           delete [] candidateList;
    2286      }
    2287      delete [] dubiousWeights;
     2225                  printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
     2226                    columnLowerWork_[i], columnUpperWork_[i],
     2227                    solution_[i], columnActivityWork_[i]);
     2228#endif
     2229                }
     2230              }
     2231              if (bad2) {
     2232                problemStatus_ = -3;
     2233                returnCode = -2;
     2234                // Force to re-factorize early next time
     2235                int numberPivots = factorization_->pivots();
     2236                forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
     2237              }
     2238            }
     2239          }
     2240        }
     2241      } else {
     2242        problemStatus_ = -3;
     2243        returnCode = -2;
     2244        // Force to re-factorize early next time
     2245        int numberPivots = factorization_->pivots();
     2246        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
     2247      }
     2248      break;
     2249    }
     2250  }
     2251  if (givenDuals) {
     2252    CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
     2253    // get rid of any values pass array
     2254    delete[] candidateList;
     2255  }
     2256  delete[] dubiousWeights;
    22882257#ifdef CLP_REPORT_PROGRESS
    2289      if (ixxxxxx > ixxyyyy - 5) {
    2290           int nTotal = numberColumns_ + numberRows_;
    2291           double oldObj = 0.0;
    2292           double newObj = 0.0;
    2293           for (int i = 0; i < nTotal; i++) {
    2294                if (savePSol[i])
    2295                     oldObj += savePSol[i] * saveCost[i];
    2296                if (solution_[i])
    2297                     newObj += solution_[i] * cost_[i];
    2298                bool printIt = false;
    2299                if (cost_[i] != saveCost[i])
    2300                     printIt = true;
    2301                if (status_[i] != saveStat[i])
    2302                     printIt = true;
    2303                if (printIt)
    2304                     printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
    2305                            i, saveStat[i], saveCost[i], savePSol[i],
    2306                            status_[i], cost_[i], solution_[i]);
    2307                // difference
    2308                savePSol[i] = solution_[i] - savePSol[i];
    2309           }
    2310           printf("exit pivots %d, old obj %g new %g\n",
    2311                  factorization_->pivots(),
    2312                  oldObj, newObj);
    2313           memset(saveDj, 0, numberRows_ * sizeof(double));
    2314           times(1.0, savePSol, saveDj);
    2315           double largest = 1.0e-6;
    2316           int k = -1;
    2317           for (int i = 0; i < numberRows_; i++) {
    2318                saveDj[i] -= savePSol[i+numberColumns_];
    2319                if (fabs(saveDj[i]) > largest) {
    2320                     largest = fabs(saveDj[i]);
    2321                     k = i;
    2322                }
    2323           }
    2324           if (k >= 0)
    2325                printf("Not null %d %g\n", k, largest);
    2326      }
    2327      delete [] savePSol ;
    2328      delete [] saveDj ;
    2329      delete [] saveCost ;
    2330      delete [] saveStat ;
    2331 #endif
    2332      return returnCode;
     2258  if (ixxxxxx > ixxyyyy - 5) {
     2259    int nTotal = numberColumns_ + numberRows_;
     2260    double oldObj = 0.0;
     2261    double newObj = 0.0;
     2262    for (int i = 0; i < nTotal; i++) {
     2263      if (savePSol[i])
     2264        oldObj += savePSol[i] * saveCost[i];
     2265      if (solution_[i])
     2266        newObj += solution_[i] * cost_[i];
     2267      bool printIt = false;
     2268      if (cost_[i] != saveCost[i])
     2269        printIt = true;
     2270      if (status_[i] != saveStat[i])
     2271        printIt = true;
     2272      if (printIt)
     2273        printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
     2274          i, saveStat[i], saveCost[i], savePSol[i],
     2275          status_[i], cost_[i], solution_[i]);
     2276      // difference
     2277      savePSol[i] = solution_[i] - savePSol[i];
     2278    }
     2279    printf("exit pivots %d, old obj %g new %g\n",
     2280      factorization_->pivots(),
     2281      oldObj, newObj);
     2282    memset(saveDj, 0, numberRows_ * sizeof(double));
     2283    times(1.0, savePSol, saveDj);
     2284    double largest = 1.0e-6;
     2285    int k = -1;
     2286    for (int i = 0; i < numberRows_; i++) {
     2287      saveDj[i] -= savePSol[i + numberColumns_];
     2288      if (fabs(saveDj[i]) > largest) {
     2289        largest = fabs(saveDj[i]);
     2290        k = i;
     2291      }
     2292    }
     2293    if (k >= 0)
     2294      printf("Not null %d %g\n", k, largest);
     2295  }
     2296  delete[] savePSol;
     2297  delete[] saveDj;
     2298  delete[] saveCost;
     2299  delete[] saveStat;
     2300#endif
     2301  return returnCode;
    23332302}
    23342303#if ABOCA_LITE
    23352304static void
    2336 updateDualBit(clpTempInfo & info)
     2305updateDualBit(clpTempInfo &info)
    23372306{
    23382307  int numberInfeasibilities = 0;
    23392308  double tolerance = info.tolerance;
    23402309  double theta = info.theta;
    2341   double * COIN_RESTRICT reducedCost = info.reducedCost;
    2342   const double * COIN_RESTRICT lower = info.lower;
    2343   const double * COIN_RESTRICT upper = info.upper;
    2344   double * COIN_RESTRICT work = info.work;
     2310  double *COIN_RESTRICT reducedCost = info.reducedCost;
     2311  const double *COIN_RESTRICT lower = info.lower;
     2312  const double *COIN_RESTRICT upper = info.upper;
     2313  double *COIN_RESTRICT work = info.work;
    23452314  int number = info.numberToDo;
    2346   int * COIN_RESTRICT which = info.which;
    2347   const unsigned char * COIN_RESTRICT statusArray = info.status;
    2348   double multiplier[] = { -1.0, 1.0};
     2315  int *COIN_RESTRICT which = info.which;
     2316  const unsigned char *COIN_RESTRICT statusArray = info.status;
     2317  double multiplier[] = { -1.0, 1.0 };
    23492318  for (int i = 0; i < number; i++) {
    23502319    int iSequence = which[i];
    23512320    double alphaI = work[i];
    23522321    work[i] = 0.0;
    2353    
     2322
    23542323    int iStatus = (statusArray[iSequence] & 3) - 1;
    23552324    if (iStatus) {
    23562325      double value = reducedCost[iSequence] - theta * alphaI;
    2357       reducedCost[iSequence] = value; 
     2326      reducedCost[iSequence] = value;
    23582327      //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
    2359       double mult = multiplier[iStatus-1];
     2328      double mult = multiplier[iStatus - 1];
    23602329      value *= mult;
    23612330      // skip if free
    2362       if (value < -tolerance&&iStatus > 0) {
    2363         // flipping bounds
    2364         double movement = mult * (upper[iSequence] - lower[iSequence]);
    2365         work[numberInfeasibilities] = movement;
    2366         which[numberInfeasibilities++] = iSequence;
    2367       }
    2368     }
    2369   }
    2370   info.numberInfeasibilities=numberInfeasibilities;
     2331      if (value < -tolerance && iStatus > 0) {
     2332        // flipping bounds
     2333        double movement = mult * (upper[iSequence] - lower[iSequence]);
     2334        work[numberInfeasibilities] = movement;
     2335        which[numberInfeasibilities++] = iSequence;
     2336      }
     2337    }
     2338  }
     2339  info.numberInfeasibilities = numberInfeasibilities;
    23712340}
    23722341#endif
     
    23752344   rowArray and columnarray will have flipped
    23762345   The output vector has movement (row length array) */
    2377 int
    2378 ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
    2379                                   CoinIndexedVector * columnArray,
    2380                                   CoinIndexedVector * outputArray,
    2381                                   double theta,
    2382                                   double & objectiveChange,
    2383                                   bool fullRecompute)
     2346int ClpSimplexDual::updateDualsInDual(CoinIndexedVector *rowArray,
     2347  CoinIndexedVector *columnArray,
     2348  CoinIndexedVector *outputArray,
     2349  double theta,
     2350  double &objectiveChange,
     2351  bool fullRecompute)
    23842352{
    23852353
    2386      outputArray->clear();
    2387 
    2388 
    2389      int numberInfeasibilities = 0;
    2390      int numberRowInfeasibilities = 0;
    2391 
    2392      // get a tolerance
    2393      double tolerance = dualTolerance_;
    2394      // we can't really trust infeasibilities if there is dual error
    2395      double error = CoinMin(1.0e-2, largestDualError_);
    2396      // allow tolerance at least slightly bigger than standard
    2397      tolerance = tolerance +  error;
    2398 
    2399      double changeObj = 0.0;
    2400 
    2401      // Coding is very similar but we can save a bit by splitting
    2402      // Do rows
    2403      if (!fullRecompute) {
    2404           int i;
    2405           double * COIN_RESTRICT reducedCost = djRegion(0);
    2406           const double * COIN_RESTRICT lower = lowerRegion(0);
    2407           const double * COIN_RESTRICT upper = upperRegion(0);
    2408           const double * COIN_RESTRICT cost = costRegion(0);
    2409           double * COIN_RESTRICT work;
    2410           int number;
    2411           int * COIN_RESTRICT which;
    2412           const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_;
    2413           assert(rowArray->packedMode());
    2414           work = rowArray->denseVector();
    2415           number = rowArray->getNumElements();
    2416           which = rowArray->getIndices();
    2417           double multiplier[] = { -1.0, 1.0};
    2418           for (i = 0; i < number; i++) {
    2419                int iSequence = which[i];
    2420                double alphaI = work[i];
    2421                work[i] = 0.0;
    2422                int iStatus = (statusArray[iSequence] & 3) - 1;
    2423                if (iStatus) {
    2424                     double value = reducedCost[iSequence] - theta * alphaI;
    2425                     // NO - can have free assert (iStatus>0);
    2426                     reducedCost[iSequence] = value;
    2427                     double mult = multiplier[iStatus-1];
    2428                     value *= mult;
    2429                     // skip if free
    2430                     if (value < -tolerance&&iStatus > 0) {
    2431                          // flipping bounds
    2432                          double movement = mult * (lower[iSequence] - upper[iSequence]);
    2433                          which[numberInfeasibilities++] = iSequence;
     2354  outputArray->clear();
     2355
     2356  int numberInfeasibilities = 0;
     2357  int numberRowInfeasibilities = 0;
     2358
     2359  // get a tolerance
     2360  double tolerance = dualTolerance_;
     2361  // we can't really trust infeasibilities if there is dual error
     2362  double error = CoinMin(1.0e-2, largestDualError_);
     2363  // allow tolerance at least slightly bigger than standard
     2364  tolerance = tolerance + error;
     2365
     2366  double changeObj = 0.0;
     2367
     2368  // Coding is very similar but we can save a bit by splitting
     2369  // Do rows
     2370  if (!fullRecompute) {
     2371    int i;
     2372    double *COIN_RESTRICT reducedCost = djRegion(0);
     2373    const double *COIN_RESTRICT lower = lowerRegion(0);
     2374    const double *COIN_RESTRICT upper = upperRegion(0);
     2375    const double *COIN_RESTRICT cost = costRegion(0);
     2376    double *COIN_RESTRICT work;
     2377    int number;
     2378    int *COIN_RESTRICT which;
     2379    const unsigned char *COIN_RESTRICT statusArray = status_ + numberColumns_;
     2380    assert(rowArray->packedMode());
     2381    work = rowArray->denseVector();
     2382    number = rowArray->getNumElements();
     2383    which = rowArray->getIndices();
     2384    double multiplier[] = { -1.0, 1.0 };
     2385    for (i = 0; i < number; i++) {
     2386      int iSequence = which[i];
     2387      double alphaI = work[i];
     2388      work[i] = 0.0;
     2389      int iStatus = (statusArray[iSequence] & 3) - 1;
     2390      if (iStatus) {
     2391        double value = reducedCost[iSequence] - theta * alphaI;
     2392        // NO - can have free assert (iStatus>0);
     2393        reducedCost[iSequence] = value;
     2394        double mult = multiplier[iStatus - 1];
     2395        value *= mult;
     2396        // skip if free
     2397        if (value < -tolerance && iStatus > 0) {
     2398          // flipping bounds
     2399          double movement = mult * (lower[iSequence] - upper[iSequence]);
     2400          which[numberInfeasibilities++] = iSequence;
    24342401#ifndef NDEBUG
    2435                          if (fabs(movement) >= 1.0e30)
    2436                               resetFakeBounds(-1000 - iSequence);
     2402          if (fabs(movement) >= 1.0e30)
     2403            resetFakeBounds(-1000 - iSequence);
    24372404#endif
    24382405#ifdef CLP_DEBUG
    2439                          if ((handler_->logLevel() & 32))
    2440                               printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2441                                      0, iSequence, value, alphaI, movement);
    2442 #endif
    2443                          changeObj -= movement * cost[iSequence];
    2444                          outputArray->quickAdd(iSequence, movement);
    2445                     }
    2446                }
    2447           }
    2448           // Do columns
    2449           reducedCost = djRegion(1);
    2450           lower = lowerRegion(1);
    2451           upper = upperRegion(1);
    2452           cost = costRegion(1);
    2453           // set number of infeasibilities in row array
    2454           numberRowInfeasibilities = numberInfeasibilities;
    2455           rowArray->setNumElements(numberInfeasibilities);
    2456           numberInfeasibilities = 0;
    2457           work = columnArray->denseVector();
    2458           number = columnArray->getNumElements();
    2459           which = columnArray->getIndices();
    2460           if ((moreSpecialOptions_ & 8) != 0) {
    2461                const unsigned char * COIN_RESTRICT statusArray = status_;
     2406          if ((handler_->logLevel() & 32))
     2407            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2408              0, iSequence, value, alphaI, movement);
     2409#endif
     2410          changeObj -= movement * cost[iSequence];
     2411          outputArray->quickAdd(iSequence, movement);
     2412        }
     2413      }
     2414    }
     2415    // Do columns
     2416    reducedCost = djRegion(1);
     2417    lower = lowerRegion(1);
     2418    upper = upperRegion(1);
     2419    cost = costRegion(1);
     2420    // set number of infeasibilities in row array
     2421    numberRowInfeasibilities = numberInfeasibilities;
     2422    rowArray->setNumElements(numberInfeasibilities);
     2423    numberInfeasibilities = 0;
     2424    work = columnArray->denseVector();
     2425    number = columnArray->getNumElements();
     2426    which = columnArray->getIndices();
     2427    if ((moreSpecialOptions_ & 8) != 0) {
     2428      const unsigned char *COIN_RESTRICT statusArray = status_;
    24622429#if ABOCA_LITE
    2463                int numberThreads=abcState();
    2464                if (numberThreads) {
    2465                clpTempInfo info[ABOCA_LITE];
    2466                int chunk = (number+numberThreads-1)/numberThreads;
    2467                int n=0;
    2468                int * whichX = which;
    2469                for (i=0;i<numberThreads;i++) {
    2470                  info[i].theta=theta;
    2471                  info[i].tolerance=tolerance;
    2472                 info[i].reducedCost = reducedCost;
    2473                 info[i].lower = lower;
    2474                 info[i].upper = upper;
    2475                  info[i].status=statusArray;
    2476                  info[i].which=which+n;
    2477                  info[i].work=work+n;
    2478                  info[i].numberToDo=CoinMin(chunk,number-n);
    2479                 n += chunk;
    2480                }
    2481                for (i=0;i<numberThreads;i++) {
    2482                 cilk_spawn updateDualBit(info[i]);
    2483                }
    2484                cilk_sync;
    2485                for (i=0;i<numberThreads;i++) {
    2486                 int n = info[i].numberInfeasibilities;
    2487                  double * workV = info[i].work;
    2488                  int * whichV = info[i].which;
    2489                 for (int j = 0; j < n; j++) {
    2490                    int iSequence = whichV[j];
    2491                    double movement = workV[j];
    2492                    workV[j] = 0.0;
    2493                    whichX[numberInfeasibilities++]=iSequence;
     2430      int numberThreads = abcState();
     2431      if (numberThreads) {
     2432        clpTempInfo info[ABOCA_LITE];
     2433        int chunk = (number + numberThreads - 1) / numberThreads;
     2434        int n = 0;
     2435        int *whichX = which;
     2436        for (i = 0; i < numberThreads; i++) {
     2437          info[i].theta = theta;
     2438          info[i].tolerance = tolerance;
     2439          info[i].reducedCost = reducedCost;
     2440          info[i].lower = lower;
     2441          info[i].upper = upper;
     2442          info[i].status = statusArray;
     2443          info[i].which = which + n;
     2444          info[i].work = work + n;
     2445          info[i].numberToDo = CoinMin(chunk, number - n);
     2446          n += chunk;
     2447        }
     2448        for (i = 0; i < numberThreads; i++) {
     2449          cilk_spawn updateDualBit(info[i]);
     2450        }
     2451        cilk_sync;
     2452        for (i = 0; i < numberThreads; i++) {
     2453          int n = info[i].numberInfeasibilities;
     2454          double *workV = info[i].work;
     2455          int *whichV = info[i].which;
     2456          for (int j = 0; j < n; j++) {
     2457            int iSequence = whichV[j];
     2458            double movement = workV[j];
     2459            workV[j] = 0.0;
     2460            whichX[numberInfeasibilities++] = iSequence;
    24942461#ifndef NDEBUG
    2495                    if (fabs(movement) >= 1.0e30)
    2496                      resetFakeBounds(-1000 - iSequence);
    2497 #endif
    2498                     changeObj += movement * cost[iSequence];
    2499                     matrix_->add(this, outputArray, iSequence, movement);
    2500                 }
    2501                }
    2502                } else {
    2503 #endif
    2504                for (i = 0; i < number; i++) {
    2505                     int iSequence = which[i];
    2506                     double alphaI = work[i];
    2507                     work[i] = 0.0;
    2508 
    2509                     int iStatus = (statusArray[iSequence] & 3) - 1;
    2510                     if (iStatus) {
    2511                          double value = reducedCost[iSequence] - theta * alphaI;
    2512                          assert (iStatus>0);
    2513                          reducedCost[iSequence] = value;
    2514                         //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
    2515                          double mult = multiplier[iStatus-1];
    2516                          value *= mult;
    2517                         // skip if free
    2518                          if (value < -tolerance&&iStatus > 0) {
    2519                               // flipping bounds
    2520                               double movement = mult * (upper[iSequence] - lower[iSequence]);
    2521                               which[numberInfeasibilities++] = iSequence;
     2462            if (fabs(movement) >= 1.0e30)
     2463              resetFakeBounds(-1000 - iSequence);
     2464#endif
     2465            changeObj += movement * cost[iSequence];
     2466            matrix_->add(this, outputArray, iSequence, movement);
     2467          }
     2468        }
     2469      } else {
     2470#endif
     2471        for (i = 0; i < number; i++) {
     2472          int iSequence = which[i];
     2473          double alphaI = work[i];
     2474          work[i] = 0.0;
     2475
     2476          int iStatus = (statusArray[iSequence] & 3) - 1;
     2477          if (iStatus) {
     2478            double value = reducedCost[iSequence] - theta * alphaI;
     2479            assert(iStatus > 0);
     2480            reducedCost[iSequence] = value;
     2481            //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
     2482            double mult = multiplier[iStatus - 1];
     2483            value *= mult;
     2484            // skip if free
     2485            if (value < -tolerance && iStatus > 0) {
     2486              // flipping bounds
     2487              double movement = mult * (upper[iSequence] - lower[iSequence]);
     2488              which[numberInfeasibilities++] = iSequence;
    25222489#ifndef NDEBUG
    2523                               if (fabs(movement) >= 1.0e30)
    2524                                    resetFakeBounds(-1000 - iSequence);
     2490              if (fabs(movement) >= 1.0e30)
     2491                resetFakeBounds(-1000 - iSequence);
    25252492#endif
    25262493#ifdef CLP_DEBUG
    2527                               if ((handler_->logLevel() & 32))
    2528                                    printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2529                                           1, iSequence, value, alphaI, movement);
    2530 #endif
    2531                               changeObj += movement * cost[iSequence];
    2532                               matrix_->add(this, outputArray, iSequence, movement);
    2533                          }
    2534                     }
    2535                }
     2494              if ((handler_->logLevel() & 32))
     2495                printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2496                  1, iSequence, value, alphaI, movement);
     2497#endif
     2498              changeObj += movement * cost[iSequence];
     2499              matrix_->add(this, outputArray, iSequence, movement);
     2500            }
     2501          }
     2502        }
    25362503#if ABOCA_LITE
    2537                }
    2538 #endif
    2539           } else {
    2540                for (i = 0; i < number; i++) {
    2541                     int iSequence = which[i];
    2542                     double alphaI = work[i];
    2543                     work[i] = 0.0;
    2544 
    2545                     Status status = getStatus(iSequence);
    2546                     if (status == atLowerBound) {
    2547                          double value = reducedCost[iSequence] - theta * alphaI;
    2548                          reducedCost[iSequence] = value;
    2549                          double movement = 0.0;
    2550 
    2551                          if (value < -tolerance) {
    2552                               // to upper bound
    2553                               which[numberInfeasibilities++] = iSequence;
    2554                               movement = upper[iSequence] - lower[iSequence];
     2504      }
     2505#endif
     2506    } else {
     2507      for (i = 0; i < number; i++) {
     2508        int iSequence = which[i];
     2509        double alphaI = work[i];
     2510        work[i] = 0.0;
     2511
     2512        Status status = getStatus(iSequence);
     2513        if (status == atLowerBound) {
     2514          double value = reducedCost[iSequence] - theta * alphaI;
     2515          reducedCost[iSequence] = value;
     2516          double movement = 0.0;
     2517
     2518          if (value < -tolerance) {
     2519            // to upper bound
     2520            which[numberInfeasibilities++] = iSequence;
     2521            movement = upper[iSequence] - lower[iSequence];
    25552522#ifndef NDEBUG
    2556                               if (fabs(movement) >= 1.0e30)
    2557                                    resetFakeBounds(-1000 - iSequence);
     2523            if (fabs(movement) >= 1.0e30)
     2524              resetFakeBounds(-1000 - iSequence);
    25582525#endif
    25592526#ifdef CLP_DEBUG
    2560                               if ((handler_->logLevel() & 32))
    2561                                    printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2562                                           1, iSequence, value, alphaI, movement);
    2563 #endif
    2564                               changeObj += movement * cost[iSequence];
    2565                               matrix_->add(this, outputArray, iSequence, movement);
    2566                          }
    2567                     } else if (status == atUpperBound) {
    2568                          double value = reducedCost[iSequence] - theta * alphaI;
    2569                          reducedCost[iSequence] = value;
    2570                          double movement = 0.0;
    2571 
    2572                          if (value > tolerance) {
    2573                               // to lower bound (if swap)
    2574                               which[numberInfeasibilities++] = iSequence;
    2575                               movement = lower[iSequence] - upper[iSequence];
     2527            if ((handler_->logLevel() & 32))
     2528              printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2529                1, iSequence, value, alphaI, movement);
     2530#endif
     2531            changeObj += movement * cost[iSequence];
     2532            matrix_->add(this, outputArray, iSequence, movement);
     2533          }
     2534        } else if (status == atUpperBound) {
     2535          double value = reducedCost[iSequence] - theta * alphaI;
     2536          reducedCost[iSequence] = value;
     2537          double movement = 0.0;
     2538
     2539          if (value > tolerance) {
     2540            // to lower bound (if swap)
     2541            which[numberInfeasibilities++] = iSequence;
     2542            movement = lower[iSequence] - upper[iSequence];
    25762543#ifndef NDEBUG
    2577                               if (fabs(movement) >= 1.0e30)
    2578                                    resetFakeBounds(-1000 - iSequence);
     2544            if (fabs(movement) >= 1.0e30)
     2545              resetFakeBounds(-1000 - iSequence);
    25792546#endif
    25802547#ifdef CLP_DEBUG
    2581                               if ((handler_->logLevel() & 32))
    2582                                    printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2583                                           1, iSequence, value, alphaI, movement);
    2584 #endif
    2585                               changeObj += movement * cost[iSequence];
    2586                               matrix_->add(this, outputArray, iSequence, movement);
    2587                          }
    2588                     } else if (status == isFree) {
    2589                          double value = reducedCost[iSequence] - theta * alphaI;
    2590                          reducedCost[iSequence] = value;
    2591                     }
    2592                }
    2593           }
    2594      } else {
    2595           double * COIN_RESTRICT solution = solutionRegion(0);
    2596           double * COIN_RESTRICT reducedCost = djRegion(0);
    2597           double * COIN_RESTRICT lower = lowerRegion(0);
    2598           double * COIN_RESTRICT upper = upperRegion(0);
    2599           const double * COIN_RESTRICT cost = costRegion(0);
    2600           int * COIN_RESTRICT which;
    2601           which = rowArray->getIndices();
    2602           int iSequence;
    2603           for (iSequence = 0; iSequence < numberRows_; iSequence++) {
    2604                double value = reducedCost[iSequence];
    2605 
    2606                Status status = getStatus(iSequence + numberColumns_);
    2607                // more likely to be at upper bound ?
    2608                if (status == atUpperBound) {
    2609                     double movement = 0.0;
    2610                     //#define NO_SWAP7
    2611                     if (value > tolerance) {
    2612                          // to lower bound (if swap)
    2613                          // put back alpha
    2614                          which[numberInfeasibilities++] = iSequence;
    2615                          movement = lower[iSequence] - upper[iSequence];
     2548            if ((handler_->logLevel() & 32))
     2549              printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2550                1, iSequence, value, alphaI, movement);
     2551#endif
     2552            changeObj += movement * cost[iSequence];
     2553            matrix_->add(this, outputArray, iSequence, movement);
     2554          }
     2555        } else if (status == isFree) {
     2556          double value = reducedCost[iSequence] - theta * alphaI;
     2557          reducedCost[iSequence] = value;
     2558        }
     2559      }
     2560    }
     2561  } else {
     2562    double *COIN_RESTRICT solution = solutionRegion(0);
     2563    double *COIN_RESTRICT reducedCost = djRegion(0);
     2564    double *COIN_RESTRICT lower = lowerRegion(0);
     2565    double *COIN_RESTRICT upper = upperRegion(0);
     2566    const double *COIN_RESTRICT cost = costRegion(0);
     2567    int *COIN_RESTRICT which;
     2568    which = rowArray->getIndices();
     2569    int iSequence;
     2570    for (iSequence = 0; iSequence < numberRows_; iSequence++) {
     2571      double value = reducedCost[iSequence];
     2572
     2573      Status status = getStatus(iSequence + numberColumns_);
     2574      // more likely to be at upper bound ?
     2575      if (status == atUpperBound) {
     2576        double movement = 0.0;
     2577        //#define NO_SWAP7
     2578        if (value > tolerance) {
     2579          // to lower bound (if swap)
     2580          // put back alpha
     2581          which[numberInfeasibilities++] = iSequence;
     2582          movement = lower[iSequence] - upper[iSequence];
    26162583#define TRY_SET_FAKE
    26172584#ifdef TRY_SET_FAKE
    2618                         if (fabs(movement) > dualBound_) {
    2619                            FakeBound bound = getFakeBound(iSequence + numberColumns_);
    2620                            if (bound == ClpSimplexDual::noFake) {
    2621                              setFakeBound(iSequence + numberColumns_,
    2622                                           ClpSimplexDual::lowerFake);
    2623                              lower[iSequence] = upper[iSequence] - dualBound_;
    2624                              assert (fabs(lower[iSequence])<1.0e30);
    2625                              movement = lower[iSequence] - upper[iSequence];
    2626                              numberFake_++;
     2585          if (fabs(movement) > dualBound_) {
     2586            FakeBound bound = getFakeBound(iSequence + numberColumns_);
     2587            if (bound == ClpSimplexDual::noFake) {
     2588              setFakeBound(iSequence + numberColumns_,
     2589                ClpSimplexDual::lowerFake);
     2590              lower[iSequence] = upper[iSequence] - dualBound_;
     2591              assert(fabs(lower[iSequence]) < 1.0e30);
     2592              movement = lower[iSequence] - upper[iSequence];
     2593              numberFake_++;
    26272594#ifndef NDEBUG
    2628                            } else {
    2629                              if (fabs(movement) >= 1.0e30)
    2630                                resetFakeBounds(-1000 - iSequence);
    2631 #endif
    2632                            }
    2633                         }
    2634 #endif
    2635                          changeObj += movement * cost[iSequence];
    2636                          outputArray->quickAdd(iSequence, -movement);
     2595            } else {
     2596              if (fabs(movement) >= 1.0e30)
     2597                resetFakeBounds(-1000 - iSequence);
     2598#endif
     2599            }
     2600          }
     2601#endif
     2602          changeObj += movement * cost[iSequence];
     2603          outputArray->quickAdd(iSequence, -movement);
    26372604#ifndef NO_SWAP7
    2638                     } else if (value > -tolerance) {
    2639                          // at correct bound but may swap
    2640                          FakeBound bound = getFakeBound(iSequence + numberColumns_);
    2641                          if (bound == ClpSimplexDual::upperFake) {
    2642                               movement = lower[iSequence] - upper[iSequence];
     2605        } else if (value > -tolerance) {
     2606          // at correct bound but may swap
     2607          FakeBound bound = getFakeBound(iSequence + numberColumns_);
     2608          if (bound == ClpSimplexDual::upperFake) {
     2609            movement = lower[iSequence] - upper[iSequence];
    26432610#ifndef NDEBUG
    2644                               if (fabs(movement) >= 1.0e30)
    2645                                 resetFakeBounds(-1000 - iSequence);
    2646 #endif
    2647                               setStatus(iSequence + numberColumns_, atLowerBound);
    2648                               solution[iSequence] = lower[iSequence];
    2649                               changeObj += movement * cost[iSequence];
    2650                               //numberFake_--;
    2651                               //setFakeBound(iSequence+numberColumns_,noFake);
    2652                          }
    2653 #endif
    2654                     }
    2655                } else if (status == atLowerBound) {
    2656                     double movement = 0.0;
    2657 
    2658                     if (value < -tolerance) {
    2659                          // to upper bound
    2660                          // put back alpha
    2661                          which[numberInfeasibilities++] = iSequence;
    2662                          movement = upper[iSequence] - lower[iSequence];
     2611            if (fabs(movement) >= 1.0e30)
     2612              resetFakeBounds(-1000 - iSequence);
     2613#endif
     2614            setStatus(iSequence + numberColumns_, atLowerBound);
     2615            solution[iSequence] = lower[iSequence];
     2616            changeObj += movement * cost[iSequence];
     2617            //numberFake_--;
     2618            //setFakeBound(iSequence+numberColumns_,noFake);
     2619          }
     2620#endif
     2621        }
     2622      } else if (status == atLowerBound) {
     2623        double movement = 0.0;
     2624
     2625        if (value < -tolerance) {
     2626          // to upper bound
     2627          // put back alpha
     2628          which[numberInfeasibilities++] = iSequence;
     2629          movement = upper[iSequence] - lower[iSequence];
    26632630#ifdef TRY_SET_FAKE
    2664                         if (fabs(movement) > dualBound_) {
    2665                            FakeBound bound = getFakeBound(iSequence + numberColumns_);
    2666                            if (bound == ClpSimplexDual::noFake) {
    2667                              setFakeBound(iSequence + numberColumns_,
    2668                                           ClpSimplexDual::upperFake);
    2669                              upper[iSequence] = lower[iSequence] + dualBound_;
    2670                              assert (fabs(upper[iSequence])<1.0e30);
    2671                              movement = upper[iSequence] - lower[iSequence];
    2672                              numberFake_++;
     2631          if (fabs(movement) > dualBound_) {
     2632            FakeBound bound = getFakeBound(iSequence + numberColumns_);
     2633            if (bound == ClpSimplexDual::noFake) {
     2634              setFakeBound(iSequence + numberColumns_,
     2635                ClpSimplexDual::upperFake);
     2636              upper[iSequence] = lower[iSequence] + dualBound_;
     2637              assert(fabs(upper[iSequence]) < 1.0e30);
     2638              movement = upper[iSequence] - lower[iSequence];
     2639              numberFake_++;
    26732640#ifndef NDEBUG
    2674                            } else {
    2675                              if (fabs(movement) >= 1.0e30)
    2676                                resetFakeBounds(-1000 - iSequence);
    2677 #endif
    2678                            }
    2679                         }
    2680 #endif
    2681                          changeObj += movement * cost[iSequence];
    2682                          outputArray->quickAdd(iSequence, -movement);
     2641            } else {
     2642              if (fabs(movement) >= 1.0e30)
     2643                resetFakeBounds(-1000 - iSequence);
     2644#endif
     2645            }
     2646          }
     2647#endif
     2648          changeObj += movement * cost[iSequence];
     2649          outputArray->quickAdd(iSequence, -movement);
    26832650#ifndef NO_SWAP7
    2684                     } else if (value < tolerance) {
    2685                          // at correct bound but may swap
    2686                          FakeBound bound = getFakeBound(iSequence + numberColumns_);
    2687                          if (bound == ClpSimplexDual::lowerFake) {
    2688                               movement = upper[iSequence] - lower[iSequence];
     2651        } else if (value < tolerance) {
     2652          // at correct bound but may swap
     2653          FakeBound bound = getFakeBound(iSequence + numberColumns_);
     2654          if (bound == ClpSimplexDual::lowerFake) {
     2655            movement = upper[iSequence] - lower[iSequence];
    26892656#ifndef NDEBUG
    2690                               if (fabs(movement) >= 1.0e30)
    2691                                 resetFakeBounds(-1000 - iSequence);
    2692 #endif
    2693                               setStatus(iSequence + numberColumns_, atUpperBound);
    2694                               solution[iSequence] = upper[iSequence];
    2695                               changeObj += movement * cost[iSequence];
    2696                               //numberFake_--;
    2697                               //setFakeBound(iSequence+numberColumns_,noFake);
    2698                          }
    2699 #endif
    2700                     }
    2701                }
    2702           }
    2703           // Do columns
    2704           solution = solutionRegion(1);
    2705           reducedCost = djRegion(1);
    2706           lower = lowerRegion(1);
    2707           upper = upperRegion(1);
    2708           cost = costRegion(1);
    2709           // set number of infeasibilities in row array
    2710           numberRowInfeasibilities = numberInfeasibilities;
    2711           rowArray->setNumElements(numberInfeasibilities);
    2712           numberInfeasibilities = 0;
    2713           which = columnArray->getIndices();
    2714           for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
    2715                double value = reducedCost[iSequence];
    2716 
    2717                Status status = getStatus(iSequence);
    2718                if (status == atLowerBound) {
    2719                     double movement = 0.0;
    2720 
    2721                     if (value < -tolerance) {
    2722                          // to upper bound
    2723                          // put back alpha
    2724                          which[numberInfeasibilities++] = iSequence;
    2725                          movement = upper[iSequence] - lower[iSequence];
     2657            if (fabs(movement) >= 1.0e30)
     2658              resetFakeBounds(-1000 - iSequence);
     2659#endif
     2660            setStatus(iSequence + numberColumns_, atUpperBound);
     2661            solution[iSequence] = upper[iSequence];
     2662            changeObj += movement * cost[iSequence];
     2663            //numberFake_--;
     2664            //setFakeBound(iSequence+numberColumns_,noFake);
     2665          }
     2666#endif
     2667        }
     2668      }
     2669    }
     2670    // Do columns
     2671    solution = solutionRegion(1);
     2672    reducedCost = djRegion(1);
     2673    lower = lowerRegion(1);
     2674    upper = upperRegion(1);
     2675    cost = costRegion(1);
     2676    // set number of infeasibilities in row array
     2677    numberRowInfeasibilities = numberInfeasibilities;
     2678    rowArray->setNumElements(numberInfeasibilities);
     2679    numberInfeasibilities = 0;
     2680    which = columnArray->getIndices();
     2681    for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
     2682      double value = reducedCost[iSequence];
     2683
     2684      Status status = getStatus(iSequence);
     2685      if (status == atLowerBound) {
     2686        double movement = 0.0;
     2687
     2688        if (value < -tolerance) {
     2689          // to upper bound
     2690          // put back alpha
     2691          which[numberInfeasibilities++] = iSequence;
     2692          movement = upper[iSequence] - lower[iSequence];
    27262693#ifdef TRY_SET_FAKE
    2727                         if (fabs(movement) > dualBound_) {
    2728                            FakeBound bound = getFakeBound(iSequence);
    2729                            if (bound == ClpSimplexDual::noFake) {
    2730                              setFakeBound(iSequence,
    2731                                           ClpSimplexDual::upperFake);
    2732                              upper[iSequence] = lower[iSequence] + dualBound_;
    2733                              assert (fabs(upper[iSequence])<1.0e30);
    2734                              movement = upper[iSequence] - lower[iSequence];
    2735                              numberFake_++;
     2694          if (fabs(movement) > dualBound_) {
     2695            FakeBound bound = getFakeBound(iSequence);
     2696            if (bound == ClpSimplexDual::noFake) {
     2697              setFakeBound(iSequence,
     2698                ClpSimplexDual::upperFake);
     2699              upper[iSequence] = lower[iSequence] + dualBound_;
     2700              assert(fabs(upper[iSequence]) < 1.0e30);
     2701              movement = upper[iSequence] - lower[iSequence];
     2702              numberFake_++;
    27362703#ifndef NDEBUG
    2737                            } else {
    2738                              if (fabs(movement) >= 1.0e30)
    2739                                resetFakeBounds(-1000 - iSequence);
    2740 #endif
    2741                            }
    2742                         }
    2743 #endif
    2744                          changeObj += movement * cost[iSequence];
    2745                          matrix_->add(this, outputArray, iSequence, movement);
     2704            } else {
     2705              if (fabs(movement) >= 1.0e30)
     2706                resetFakeBounds(-1000 - iSequence);
     2707#endif
     2708            }
     2709          }
     2710#endif
     2711          changeObj += movement * cost[iSequence];
     2712          matrix_->add(this, outputArray, iSequence, movement);
    27462713#ifndef NO_SWAP7
    2747                     } else if (value < tolerance) {
    2748                          // at correct bound but may swap
    2749                          FakeBound bound = getFakeBound(iSequence);
    2750                          if (bound == ClpSimplexDual::lowerFake) {
    2751                               movement = upper[iSequence] - lower[iSequence];
     2714        } else if (value < tolerance) {
     2715          // at correct bound but may swap
     2716          FakeBound bound = getFakeBound(iSequence);
     2717          if (bound == ClpSimplexDual::lowerFake) {
     2718            movement = upper[iSequence] - lower[iSequence];
    27522719#ifndef NDEBUG
    2753                               if (fabs(movement) >= 1.0e30)
    2754                                 resetFakeBounds(-1000 - iSequence);
    2755 #endif
    2756                               setStatus(iSequence, atUpperBound);
    2757                               solution[iSequence] = upper[iSequence];
    2758                               changeObj += movement * cost[iSequence];
    2759                               //numberFake_--;
    2760                               //setFakeBound(iSequence,noFake);
    2761                          }
    2762 #endif
    2763                     }
    2764                } else if (status == atUpperBound) {
    2765                     double movement = 0.0;
    2766 
    2767                     if (value > tolerance) {
    2768                          // to lower bound (if swap)
    2769                          // put back alpha
    2770                          which[numberInfeasibilities++] = iSequence;
    2771                          movement = lower[iSequence] - upper[iSequence];
     2720            if (fabs(movement) >= 1.0e30)
     2721              resetFakeBounds(-1000 - iSequence);
     2722#endif
     2723            setStatus(iSequence, atUpperBound);
     2724            solution[iSequence] = upper[iSequence];
     2725            changeObj += movement * cost[iSequence];
     2726            //numberFake_--;
     2727            //setFakeBound(iSequence,noFake);
     2728          }
     2729#endif
     2730        }
     2731      } else if (status == atUpperBound) {
     2732        double movement = 0.0;
     2733
     2734        if (value > tolerance) {
     2735          // to lower bound (if swap)
     2736          // put back alpha
     2737          which[numberInfeasibilities++] = iSequence;
     2738          movement = lower[iSequence] - upper[iSequence];
    27722739#ifdef TRY_SET_FAKE
    2773                         if (fabs(movement) > dualBound_) {
    2774                            FakeBound bound = getFakeBound(iSequence);
    2775                            if (bound == ClpSimplexDual::noFake) {
    2776                              setFakeBound(iSequence,
    2777                                           ClpSimplexDual::lowerFake);
    2778                              lower[iSequence] = upper[iSequence] - dualBound_;
    2779                              assert (fabs(lower[iSequence])<1.0e30);
    2780                              movement = lower[iSequence] - upper[iSequence];
    2781                              numberFake_++;
     2740          if (fabs(movement) > dualBound_) {
     2741            FakeBound bound = getFakeBound(iSequence);
     2742            if (bound == ClpSimplexDual::noFake) {
     2743              setFakeBound(iSequence,
     2744                ClpSimplexDual::lowerFake);
     2745              lower[iSequence] = upper[iSequence] - dualBound_;
     2746              assert(fabs(lower[iSequence]) < 1.0e30);
     2747              movement = lower[iSequence] - upper[iSequence];
     2748              numberFake_++;
    27822749#ifndef NDEBUG
    2783                            } else {
    2784                              if (fabs(movement) >= 1.0e30)
    2785                                resetFakeBounds(-1000 - iSequence);
    2786 #endif
    2787                            }
    2788                         }
    2789 #endif
    2790                          changeObj += movement * cost[iSequence];
    2791                          matrix_->add(this, outputArray, iSequence, movement);
     2750            } else {
     2751              if (fabs(movement) >= 1.0e30)
     2752                resetFakeBounds(-1000 - iSequence);
     2753#endif
     2754            }
     2755          }
     2756#endif
     2757          changeObj += movement * cost[iSequence];
     2758          matrix_->add(this, outputArray, iSequence, movement);
    27922759#ifndef NO_SWAP7
    2793                     } else if (value > -tolerance) {
    2794                          // at correct bound but may swap
    2795                          FakeBound bound = getFakeBound(iSequence);
    2796                          if (bound == ClpSimplexDual::upperFake) {
    2797                               movement = lower[iSequence] - upper[iSequence];
     2760        } else if (value > -tolerance) {
     2761          // at correct bound but may swap
     2762          FakeBound bound = getFakeBound(iSequence);
     2763          if (bound == ClpSimplexDual::upperFake) {
     2764            movement = lower[iSequence] - upper[iSequence];
    27982765#ifndef NDEBUG
    2799                               if (fabs(movement) >= 1.0e30)
    2800                                 resetFakeBounds(-1000 - iSequence);
    2801 #endif
    2802                               setStatus(iSequence, atLowerBound);
    2803                               solution[iSequence] = lower[iSequence];
    2804                               changeObj += movement * cost[iSequence];
    2805                               //numberFake_--;
    2806                               //setFakeBound(iSequence,noFake);
    2807                          }
    2808 #endif
    2809                     }
    2810                }
    2811           }
    2812      }
     2766            if (fabs(movement) >= 1.0e30)
     2767              resetFakeBounds(-1000 - iSequence);
     2768#endif
     2769            setStatus(iSequence, atLowerBound);
     2770            solution[iSequence] = lower[iSequence];
     2771            changeObj += movement * cost[iSequence];
     2772            //numberFake_--;
     2773            //setFakeBound(iSequence,noFake);
     2774          }
     2775#endif
     2776        }
     2777      }
     2778    }
     2779  }
    28132780
    28142781#ifdef CLP_DEBUG
    2815      if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
    2816           printf("%d fake after full update\n", numberFake_);
    2817 #endif
    2818      // set number of infeasibilities
    2819      columnArray->setNumElements(numberInfeasibilities);
    2820      numberInfeasibilities += numberRowInfeasibilities;
    2821      if (fullRecompute) {
    2822           // do actual flips
    2823           flipBounds(rowArray, columnArray);
    2824      }
    2825      objectiveChange += changeObj;
    2826      return numberInfeasibilities;
     2782  if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
     2783    printf("%d fake after full update\n", numberFake_);
     2784#endif
     2785  // set number of infeasibilities
     2786  columnArray->setNumElements(numberInfeasibilities);
     2787  numberInfeasibilities += numberRowInfeasibilities;
     2788  if (fullRecompute) {
     2789    // do actual flips
     2790    flipBounds(rowArray, columnArray);
     2791  }
     2792  objectiveChange += changeObj;
     2793  return numberInfeasibilities;
    28272794}
    2828 void
    2829 ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
    2830                                         CoinIndexedVector * columnArray,
    2831                                         double theta)
     2795void ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector *rowArray,
     2796  CoinIndexedVector *columnArray,
     2797  double theta)
    28322798{
    28332799
    2834      // use a tighter tolerance except for all being okay
    2835      double tolerance = dualTolerance_;
    2836 
    2837      // Coding is very similar but we can save a bit by splitting
    2838      // Do rows
    2839      {
    2840           int i;
    2841           double * reducedCost = djRegion(0);
    2842           double * work;
    2843           int number;
    2844           int * which;
    2845           work = rowArray->denseVector();
    2846           number = rowArray->getNumElements();
    2847           which = rowArray->getIndices();
    2848           for (i = 0; i < number; i++) {
    2849                int iSequence = which[i];
    2850                double alphaI = work[i];
    2851                double value = reducedCost[iSequence] - theta * alphaI;
    2852                work[i] = 0.0;
    2853                reducedCost[iSequence] = value;
    2854 
    2855                Status status = getStatus(iSequence + numberColumns_);
    2856                // more likely to be at upper bound ?
    2857                if (status == atUpperBound) {
    2858 
    2859                     if (value > tolerance)
    2860                          reducedCost[iSequence] = 0.0;
    2861                } else if (status == atLowerBound) {
    2862 
    2863                     if (value < -tolerance) {
    2864                          reducedCost[iSequence] = 0.0;
    2865                     }
    2866                }
    2867           }
    2868      }
    2869      rowArray->setNumElements(0);
    2870 
    2871      // Do columns
    2872      {
    2873           int i;
    2874           double * reducedCost = djRegion(1);
    2875           double * work;
    2876           int number;
    2877           int * which;
    2878           work = columnArray->denseVector();
    2879           number = columnArray->getNumElements();
    2880           which = columnArray->getIndices();
    2881 
    2882           for (i = 0; i < number; i++) {
    2883                int iSequence = which[i];
    2884                double alphaI = work[i];
    2885                double value = reducedCost[iSequence] - theta * alphaI;
    2886                work[i] = 0.0;
    2887                reducedCost[iSequence] = value;
    2888 
    2889                Status status = getStatus(iSequence);
    2890                if (status == atLowerBound) {
    2891                     if (value < -tolerance)
    2892                          reducedCost[iSequence] = 0.0;
    2893                } else if (status == atUpperBound) {
    2894                     if (value > tolerance)
    2895                          reducedCost[iSequence] = 0.0;
    2896                }
    2897           }
    2898      }
    2899      columnArray->setNumElements(0);
     2800  // use a tighter tolerance except for all being okay
     2801  double tolerance = dualTolerance_;
     2802
     2803  // Coding is very similar but we can save a bit by splitting
     2804  // Do rows
     2805  {
     2806    int i;
     2807    double *reducedCost = djRegion(0);
     2808    double *work;
     2809    int number;
     2810    int *which;
     2811    work = rowArray->denseVector();
     2812    number = rowArray->getNumElements();
     2813    which = rowArray->getIndices();
     2814    for (i = 0; i < number; i++) {
     2815      int iSequence = which[i];
     2816      double alphaI = work[i];
     2817      double value = reducedCost[iSequence] - theta * alphaI;
     2818      work[i] = 0.0;
     2819      reducedCost[iSequence] = value;
     2820
     2821      Status status = getStatus(iSequence + numberColumns_);
     2822      // more likely to be at upper bound ?
     2823      if (status == atUpperBound) {
     2824
     2825        if (value > tolerance)
     2826          reducedCost[iSequence] = 0.0;
     2827      } else if (status == atLowerBound) {
     2828
     2829        if (value < -tolerance) {
     2830          reducedCost[iSequence] = 0.0;
     2831        }
     2832      }
     2833    }
     2834  }
     2835  rowArray->setNumElements(0);
     2836
     2837  // Do columns
     2838  {
     2839    int i;
     2840    double *reducedCost = djRegion(1);
     2841    double *work;
     2842    int number;
     2843    int *which;
     2844    work = columnArray->denseVector();
     2845    number = columnArray->getNumElements();
     2846    which = columnArray->getIndices();
     2847
     2848    for (i = 0; i < number; i++) {
     2849      int iSequence = which[i];
     2850      double alphaI = work[i];
     2851      double value = reducedCost[iSequence] - theta * alphaI;
     2852      work[i] = 0.0;
     2853      reducedCost[iSequence] = value;
     2854
     2855      Status status = getStatus(iSequence);
     2856      if (status == atLowerBound) {
     2857        if (value < -tolerance)
     2858          reducedCost[iSequence] = 0.0;
     2859      } else if (status == atUpperBound) {
     2860        if (value > tolerance)
     2861          reducedCost[iSequence] = 0.0;
     2862      }
     2863    }
     2864  }
     2865  columnArray->setNumElements(0);
    29002866}
    29012867/*
     
    29052871   For easy problems we can just choose one of the first rows we look at
    29062872*/
    2907 void
    2908 ClpSimplexDual::dualRow(int alreadyChosen)
     2873void ClpSimplexDual::dualRow(int alreadyChosen)
    29092874{
    2910      // get pivot row using whichever method it is
    2911      int chosenRow = -1;
     2875  // get pivot row using whichever method it is
     2876  int chosenRow = -1;
    29122877#ifdef FORCE_FOLLOW
    2913      bool forceThis = false;
    2914      if (!fpFollow && strlen(forceFile)) {
    2915           fpFollow = fopen(forceFile, "r");
    2916           assert (fpFollow);
    2917      }
    2918      if (fpFollow) {
    2919           if (numberIterations_ <= force_iteration) {
    2920                // read to next Clp0102
    2921                char temp[300];
    2922                while (fgets(temp, 250, fpFollow)) {
    2923                     if (strncmp(temp, "Clp0102", 7))
    2924                          continue;
    2925                     char cin, cout;
    2926                     sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
    2927                            &force_iteration, &cin, &force_in, &cout, &force_out);
    2928                     if (cin == 'R')
    2929                          force_in += numberColumns_;
    2930                     if (cout == 'R')
    2931                          force_out += numberColumns_;
    2932                     forceThis = true;
    2933                     assert (numberIterations_ == force_iteration - 1);
    2934                     printf("Iteration %d will force %d out and %d in\n",
    2935                            force_iteration, force_out, force_in);
    2936                     alreadyChosen = force_out;
    2937                     break;
    2938                }
    2939           } else {
    2940                // use old
    2941                forceThis = true;
    2942           }
    2943           if (!forceThis) {
    2944                fclose(fpFollow);
    2945                fpFollow = NULL;
    2946                forceFile = "";
    2947           }
    2948      }
    2949 #endif
    2950      //double freeAlpha = 0.0;
    2951      if (alreadyChosen < 0) {
    2952           // first see if any free variables and put them in basis
    2953           int nextFree = nextSuperBasic();
    2954           //nextFree=-1; //off
    2955           if (nextFree >= 0) {
    2956                // unpack vector and find a good pivot
    2957                unpack(rowArray_[1], nextFree);
    2958                factorization_->updateColumn(rowArray_[2], rowArray_[1]);
    2959 
    2960                double * work = rowArray_[1]->denseVector();
    2961                int number = rowArray_[1]->getNumElements();
    2962                int * which = rowArray_[1]->getIndices();
    2963                double bestFeasibleAlpha = 0.0;
    2964                int bestFeasibleRow = -1;
    2965                double bestInfeasibleAlpha = 0.0;
    2966                int bestInfeasibleRow = -1;
    2967                int i;
    2968 
    2969                for (i = 0; i < number; i++) {
    2970                     int iRow = which[i];
    2971                     double alpha = fabs(work[iRow]);
    2972                     if (alpha > 1.0e-3) {
    2973                          int iSequence = pivotVariable_[iRow];
    2974                          double value = solution_[iSequence];
    2975                          double lower = lower_[iSequence];
    2976                          double upper = upper_[iSequence];
    2977                          double infeasibility = 0.0;
    2978                          if (value > upper)
    2979                               infeasibility = value - upper;
    2980                          else if (value < lower)
    2981                               infeasibility = lower - value;
    2982                          if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
    2983                               if (!flagged(iSequence)) {
    2984                                    bestInfeasibleAlpha = infeasibility * alpha;
    2985                                    bestInfeasibleRow = iRow;
    2986                               }
    2987                          }
    2988                          if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
    2989                               bestFeasibleAlpha = alpha;
    2990                               bestFeasibleRow = iRow;
    2991                          }
    2992                     }
    2993                }
    2994                if (bestInfeasibleRow >= 0)
    2995                     chosenRow = bestInfeasibleRow;
    2996                else if (bestFeasibleAlpha > 1.0e-2)
    2997                     chosenRow = bestFeasibleRow;
    2998                if (chosenRow >= 0) {
    2999                     pivotRow_ = chosenRow;
    3000                     //freeAlpha = work[chosenRow];
    3001                }
    3002                rowArray_[1]->clear();
    3003           }
    3004      } else {
    3005           // in values pass
    3006           chosenRow = alreadyChosen;
     2878  bool forceThis = false;
     2879  if (!fpFollow && strlen(forceFile)) {
     2880    fpFollow = fopen(forceFile, "r");
     2881    assert(fpFollow);
     2882  }
     2883  if (fpFollow) {
     2884    if (numberIterations_ <= force_iteration) {
     2885      // read to next Clp0102
     2886      char temp[300];
     2887      while (fgets(temp, 250, fpFollow)) {
     2888        if (strncmp(temp, "Clp0102", 7))
     2889          continue;
     2890        char cin, cout;
     2891        sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
     2892          &force_iteration, &cin, &force_in, &cout, &force_out);
     2893        if (cin == 'R')
     2894          force_in += numberColumns_;
     2895        if (cout == 'R')
     2896          force_out += numberColumns_;
     2897        forceThis = true;
     2898        assert(numberIterations_ == force_iteration - 1);
     2899        printf("Iteration %d will force %d out and %d in\n",
     2900          force_iteration, force_out, force_in);
     2901        alreadyChosen = force_out;
     2902        break;
     2903      }
     2904    } else {
     2905      // use old
     2906      forceThis = true;
     2907    }
     2908    if (!forceThis) {
     2909      fclose(fpFollow);
     2910      fpFollow = NULL;
     2911      forceFile = "";
     2912    }
     2913  }
     2914#endif
     2915  //double freeAlpha = 0.0;
     2916  if (alreadyChosen < 0) {
     2917    // first see if any free variables and put them in basis
     2918    int nextFree = nextSuperBasic();
     2919    //nextFree=-1; //off
     2920    if (nextFree >= 0) {
     2921      // unpack vector and find a good pivot
     2922      unpack(rowArray_[1], nextFree);
     2923      factorization_->updateColumn(rowArray_[2], rowArray_[1]);
     2924
     2925      double *work = rowArray_[1]->denseVector();
     2926      int number = rowArray_[1]->getNumElements();
     2927      int *which = rowArray_[1]->getIndices();
     2928      double bestFeasibleAlpha = 0.0;
     2929      int bestFeasibleRow = -1;
     2930      double bestInfeasibleAlpha = 0.0;
     2931      int bestInfeasibleRow = -1;
     2932      int i;
     2933
     2934      for (i = 0; i < number; i++) {
     2935        int iRow = which[i];
     2936        double alpha = fabs(work[iRow]);
     2937        if (alpha > 1.0e-3) {
     2938          int iSequence = pivotVariable_[iRow];
     2939          double value = solution_[iSequence];
     2940          double lower = lower_[iSequence];
     2941          double upper = upper_[iSequence];
     2942          double infeasibility = 0.0;
     2943          if (value > upper)
     2944            infeasibility = value - upper;
     2945          else if (value < lower)
     2946            infeasibility = lower - value;
     2947          if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
     2948            if (!flagged(iSequence)) {
     2949              bestInfeasibleAlpha = infeasibility * alpha;
     2950              bestInfeasibleRow = iRow;
     2951            }
     2952          }
     2953          if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
     2954            bestFeasibleAlpha = alpha;
     2955            bestFeasibleRow = iRow;
     2956          }
     2957        }
     2958      }
     2959      if (bestInfeasibleRow >= 0)
     2960        chosenRow = bestInfeasibleRow;
     2961      else if (bestFeasibleAlpha > 1.0e-2)
     2962        chosenRow = bestFeasibleRow;
     2963      if (chosenRow >= 0) {
     2964        pivotRow_ = chosenRow;
     2965        //freeAlpha = work[chosenRow];
     2966      }
     2967      rowArray_[1]->clear();
     2968    }
     2969  } else {
     2970    // in values pass
     2971    chosenRow = alreadyChosen;
    30072972#ifdef FORCE_FOLLOW
    3008           if(forceThis) {
    3009                alreadyChosen = -1;
    3010                chosenRow = -1;
    3011                for (int i = 0; i < numberRows_; i++) {
    3012                     if (pivotVariable_[i] == force_out) {
    3013                          chosenRow = i;
    3014                          break;
    3015                     }
    3016                }
    3017                assert (chosenRow >= 0);
    3018           }
    3019 #endif
    3020           pivotRow_ = chosenRow;
    3021      }
    3022      if (chosenRow < 0)
    3023           pivotRow_ = dualRowPivot_->pivotRow();
    3024 
    3025      if (pivotRow_ >= 0) {
    3026           sequenceOut_ = pivotVariable_[pivotRow_];
    3027           valueOut_ = solution_[sequenceOut_];
    3028           lowerOut_ = lower_[sequenceOut_];
    3029           upperOut_ = upper_[sequenceOut_];
    3030           if (alreadyChosen < 0) {
    3031                // if we have problems we could try other way and hope we get a
    3032                // zero pivot?
    3033                if (valueOut_ > upperOut_) {
    3034                     directionOut_ = -1;
    3035                     dualOut_ = valueOut_ - upperOut_;
    3036                } else if (valueOut_ < lowerOut_) {
    3037                     directionOut_ = 1;
    3038                     dualOut_ = lowerOut_ - valueOut_;
    3039                } else {
     2973    if (forceThis) {
     2974      alreadyChosen = -1;
     2975      chosenRow = -1;
     2976      for (int i = 0; i < numberRows_; i++) {
     2977        if (pivotVariable_[i] == force_out) {
     2978          chosenRow = i;
     2979          break;
     2980        }
     2981      }
     2982      assert(chosenRow >= 0);
     2983    }
     2984#endif
     2985    pivotRow_ = chosenRow;
     2986  }
     2987  if (chosenRow < 0)
     2988    pivotRow_ = dualRowPivot_->pivotRow();
     2989
     2990  if (pivotRow_ >= 0) {
     2991    sequenceOut_ = pivotVariable_[pivotRow_];
     2992    valueOut_ = solution_[sequenceOut_];
     2993    lowerOut_ = lower_[sequenceOut_];
     2994    upperOut_ = upper_[sequenceOut_];
     2995    if (alreadyChosen < 0) {
     2996      // if we have problems we could try other way and hope we get a
     2997      // zero pivot?
     2998      if (valueOut_ > upperOut_) {
     2999        directionOut_ = -1;
     3000        dualOut_ = valueOut_ - upperOut_;
     3001      } else if (valueOut_ < lowerOut_) {
     3002        directionOut_ = 1;
     3003        dualOut_ = lowerOut_ - valueOut_;
     3004      } else {
    30403005#if 1
    3041                     // odd (could be free) - it's feasible - go to nearest
    3042                     if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
    3043                          directionOut_ = 1;
    3044                          dualOut_ = lowerOut_ - valueOut_;
    3045                     } else {
    3046                          directionOut_ = -1;
    3047                          dualOut_ = valueOut_ - upperOut_;
    3048                     }
     3006        // odd (could be free) - it's feasible - go to nearest
     3007        if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
     3008          directionOut_ = 1;
     3009          dualOut_ = lowerOut_ - valueOut_;
     3010        } else {
     3011          directionOut_ = -1;
     3012          dualOut_ = valueOut_ - upperOut_;
     3013        }
    30493014#else
    3050                     // odd (could be free) - it's feasible - improve obj
    3051                     printf("direction from alpha of %g is %d\n",
    3052                            freeAlpha, freeAlpha > 0.0 ? 1 : -1);
    3053                     if (valueOut_ - lowerOut_ > 1.0e20)
    3054                          freeAlpha = 1.0;
    3055                     else if(upperOut_ - valueOut_ > 1.0e20)
    3056                          freeAlpha = -1.0;
    3057                     //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
    3058                     if (freeAlpha < 0.0) {
    3059                          directionOut_ = 1;
    3060                          dualOut_ = lowerOut_ - valueOut_;
    3061                     } else {
    3062                          directionOut_ = -1;
    3063                          dualOut_ = valueOut_ - upperOut_;
    3064                     }
    3065                     printf("direction taken %d - bounds %g %g %g\n",
    3066                            directionOut_, lowerOut_, valueOut_, upperOut_);
    3067 #endif
    3068                }
     3015        // odd (could be free) - it's feasible - improve obj
     3016        printf("direction from alpha of %g is %d\n",
     3017          freeAlpha, freeAlpha > 0.0 ? 1 : -1);
     3018        if (valueOut_ - lowerOut_ > 1.0e20)
     3019          freeAlpha = 1.0;
     3020        else if (upperOut_ - valueOut_ > 1.0e20)
     3021          freeAlpha = -1.0;
     3022        //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
     3023        if (freeAlpha < 0.0) {
     3024          directionOut_ = 1;
     3025          dualOut_ = lowerOut_ - valueOut_;
     3026        } else {
     3027          directionOut_ = -1;
     3028          dualOut_ = valueOut_ - upperOut_;
     3029        }
     3030        printf("direction taken %d - bounds %g %g %g\n",
     3031          directionOut_, lowerOut_, valueOut_, upperOut_);
     3032#endif
     3033      }
    30693034#ifdef CLP_DEBUG
    3070                assert(dualOut_ >= 0.0);
    3071 #endif
    3072           } else {
    3073                // in values pass so just use sign of dj
    3074                // We don't want to go through any barriers so set dualOut low
    3075                // free variables will never be here
    3076                dualOut_ = 1.0e-6;
    3077                if (dj_[sequenceOut_] > 0.0) {
    3078                     // this will give a -1 in pivot row (as slacks are -1.0)
    3079                     directionOut_ = 1;
    3080                } else {
    3081                     directionOut_ = -1;
    3082                }
    3083           }
    3084      }
    3085      return ;
     3035      assert(dualOut_ >= 0.0);
     3036#endif
     3037    } else {
     3038      // in values pass so just use sign of dj
     3039      // We don't want to go through any barriers so set dualOut low
     3040      // free variables will never be here
     3041      dualOut_ = 1.0e-6;
     3042      if (dj_[sequenceOut_] > 0.0) {
     3043        // this will give a -1 in pivot row (as slacks are -1.0)
     3044        directionOut_ = 1;
     3045      } else {
     3046        directionOut_ = -1;
     3047      }
     3048    }
     3049  }
     3050  return;
    30863051}
    30873052// Checks if any fake bounds active - if so returns number and modifies
     
    30923057// Fills in changeVector which can be used to see if unbounded
    30933058// and cost of change vector
    3094 int
    3095 ClpSimplexDual::changeBounds(int initialize,
    3096                              CoinIndexedVector * outputArray,
    3097                              double & changeCost)
     3059int ClpSimplexDual::changeBounds(int initialize,
     3060  CoinIndexedVector *outputArray,
     3061  double &changeCost)
    30983062{
    3099      numberFake_ = 0;
    3100      if (!initialize) {
    3101           int numberInfeasibilities;
    3102           double newBound;
    3103           newBound = 5.0 * dualBound_;
    3104           numberInfeasibilities = 0;
    3105           changeCost = 0.0;
    3106           // put back original bounds and then check
    3107           createRim1(false);
    3108           int iSequence;
    3109           // bounds will get bigger - just look at ones at bounds
    3110           for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
    3111                double lowerValue = lower_[iSequence];
    3112                double upperValue = upper_[iSequence];
    3113                double value = solution_[iSequence];
    3114                setFakeBound(iSequence, ClpSimplexDual::noFake);
    3115                switch(getStatus(iSequence)) {
    3116 
    3117                case basic:
    3118                case ClpSimplex::isFixed:
    3119                     break;
    3120                case isFree:
    3121                case superBasic:
    3122                     break;
    3123                case atUpperBound:
    3124                     if (fabs(value - upperValue) > primalTolerance_) {
    3125                       if(fabs(dj_[iSequence])>1.0e-9) {
    3126                          numberInfeasibilities++;
    3127                       } else {
    3128                         setStatus(iSequence,superBasic);
    3129                         moreSpecialOptions_ &= ~8;
    3130                       }
    3131                     }
    3132                     break;
    3133                case atLowerBound:
    3134                     if (fabs(value - lowerValue) > primalTolerance_) {
    3135                       if(fabs(dj_[iSequence])>1.0e-9) {
    3136                          numberInfeasibilities++;
    3137                       } else {
    3138                         setStatus(iSequence,superBasic);
    3139                         moreSpecialOptions_ &= ~8;
    3140                       }
    3141                     }
    3142                     break;
    3143                }
    3144           }
    3145           // If dual infeasible then carry on
    3146           if (numberInfeasibilities) {
    3147                handler_->message(CLP_DUAL_CHECKB, messages_)
    3148                          << newBound
    3149                          << CoinMessageEol;
    3150                int iSequence;
    3151                for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
    3152                     double lowerValue = lower_[iSequence];
    3153                     double upperValue = upper_[iSequence];
    3154                     double newLowerValue;
    3155                     double newUpperValue;
    3156                     Status status = getStatus(iSequence);
    3157                     if (status == atUpperBound ||
    3158                               status == atLowerBound) {
    3159                          double value = solution_[iSequence];
    3160                          if (value - lowerValue <= upperValue - value) {
    3161                               newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
    3162                               newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
    3163                          } else {
    3164                               newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
    3165                               newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
    3166                          }
    3167                          if (newLowerValue > lowerValue) {
    3168                               if (newUpperValue < upperValue) {
    3169                                    setFakeBound(iSequence, ClpSimplexDual::bothFake);
    3170                                    // redo
    3171                                    if (status == atLowerBound) {
    3172                                      newLowerValue = value;
    3173                                      newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
    3174                                    } else {
    3175                                      newUpperValue = value;
    3176                                      newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
    3177                                    }
    3178                                    numberFake_++;
    3179                               } else {
    3180                                    setFakeBound(iSequence, ClpSimplexDual::lowerFake);
    3181                                    numberFake_++;
    3182                               }
    3183                          } else {
    3184                               if (newUpperValue < upperValue) {
    3185                                    setFakeBound(iSequence, ClpSimplexDual::upperFake);
    3186                                    numberFake_++;
    3187                               }
    3188                          }
    3189                          lower_[iSequence] = newLowerValue;
    3190                          upper_[iSequence] = newUpperValue;
    3191                          if (status == atUpperBound)
    3192                               solution_[iSequence] = newUpperValue;
    3193                          else
    3194                               solution_[iSequence] = newLowerValue;
    3195                          double movement = solution_[iSequence] - value;
    3196                          if (movement && outputArray) {
    3197                               if (iSequence >= numberColumns_) {
    3198                                    outputArray->quickAdd(iSequence, -movement);
    3199                                    changeCost += movement * cost_[iSequence];
    3200                               } else {
    3201                                    matrix_->add(this, outputArray, iSequence, movement);
    3202                                    changeCost += movement * cost_[iSequence];
    3203                               }
    3204                          }
    3205                     }
    3206                }
    3207                dualBound_ = newBound;
     3063  numberFake_ = 0;
     3064  if (!initialize) {
     3065    int numberInfeasibilities;
     3066    double newBound;
     3067    newBound = 5.0 * dualBound_;
     3068    numberInfeasibilities = 0;
     3069    changeCost = 0.0;
     3070    // put back original bounds and then check
     3071    createRim1(false);
     3072    int iSequence;
     3073    // bounds will get bigger - just look at ones at bounds
     3074    for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
     3075      double lowerValue = lower_[iSequence];
     3076      double upperValue = upper_[iSequence];
     3077      double value = solution_[iSequence];
     3078      setFakeBound(iSequence, ClpSimplexDual::noFake);
     3079      switch (getStatus(iSequence)) {
     3080
     3081      case basic:
     3082      case ClpSimplex::isFixed:
     3083        break;
     3084      case isFree:
     3085      case superBasic:
     3086        break;
     3087      case atUpperBound:
     3088        if (fabs(value - upperValue) > primalTolerance_) {
     3089          if (fabs(dj_[iSequence]) > 1.0e-9) {
     3090            numberInfeasibilities++;
    32083091          } else {
    3209                numberInfeasibilities = -1;
    3210           }
    3211           return numberInfeasibilities;
    3212      } else if (initialize == 1 || initialize == 3) {
    3213           int iSequence;
    3214           if (initialize == 3) {
    3215             if (columnScale_) {
    3216               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
    3217                 if (getFakeBound(iSequence) != ClpSimplexDual::noFake) {
    3218                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
    3219                   // lower
    3220                   double value = columnLower_[iSequence];
    3221                   if (value > -1.0e30) {
    3222                     value *= multiplier;
    3223                   }
    3224                   lower_[iSequence] = value;
    3225                   // upper
    3226                   value = columnUpper_[iSequence];
    3227                   if (value < 1.0e30) {
    3228                     value *= multiplier;
    3229                   }
    3230                   upper_[iSequence] = value;
    3231                   setFakeBound(iSequence, ClpSimplexDual::noFake);
    3232                 }
    3233               }
    3234               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
    3235                 // lower
    3236                 double multiplier = rhsScale_ * rowScale_[iSequence];
    3237                 double value = rowLower_[iSequence];
    3238                 if (value > -1.0e30) {
    3239                   value *= multiplier;
    3240                 }
    3241                 lower_[iSequence+numberColumns_] = value;
    3242                 // upper
    3243                 value = rowUpper_[iSequence];
    3244                 if (value < 1.0e30) {
    3245                   value *= multiplier;
    3246                 }
    3247                 upper_[iSequence+numberColumns_] = value;
    3248                 setFakeBound(iSequence+numberColumns_, ClpSimplexDual::noFake);
    3249               }
    3250             } else {
    3251               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
    3252                 if (getFakeBound(iSequence) != ClpSimplexDual::noFake) {
    3253                   lower_[iSequence] = columnLower_[iSequence];
    3254                   upper_[iSequence] = columnUpper_[iSequence];
    3255                   setFakeBound(iSequence, ClpSimplexDual::noFake);
    3256                 }
    3257               }
    3258               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
    3259                 if (getFakeBound(iSequence+numberColumns_) != ClpSimplexDual::noFake) {
    3260                   lower_[iSequence+numberColumns_] = rowLower_[iSequence];
    3261                   upper_[iSequence+numberColumns_] = rowUpper_[iSequence];
    3262                   setFakeBound(iSequence+numberColumns_, ClpSimplexDual::noFake);
    3263                 }
    3264               }
    3265             }
    3266           }
    3267           double testBound = 0.999999 * dualBound_;
    3268           for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
    3269                Status status = getStatus(iSequence);
    3270                if (status == atUpperBound ||
    3271                          status == atLowerBound) {
    3272                     double lowerValue = lower_[iSequence];
    3273                     double upperValue = upper_[iSequence];
    3274                     double value = solution_[iSequence];
    3275                     if (lowerValue > -largeValue_ || upperValue < largeValue_) {
    3276                          if (true || lowerValue - value > -0.5 * dualBound_ ||
    3277                                    upperValue - value < 0.5 * dualBound_) {
    3278                               if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
    3279                                    if (upperValue > lowerValue + testBound) {
    3280                                         if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
    3281                                              numberFake_++;
    3282                                         upper_[iSequence] = lowerValue + dualBound_;
    3283                                         setFakeBound(iSequence, ClpSimplexDual::upperFake);
    3284                                    }
    3285                               } else {
    3286                                    if (lowerValue < upperValue - testBound) {
    3287                                         if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
    3288                                              numberFake_++;
    3289                                         lower_[iSequence] = upperValue - dualBound_;
    3290                                         setFakeBound(iSequence, ClpSimplexDual::lowerFake);
    3291                                    }
    3292                               }
    3293                          } else {
    3294                               if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
    3295                                    numberFake_++;
    3296                               lower_[iSequence] = -0.5 * dualBound_;
    3297                               upper_[iSequence] = 0.5 * dualBound_;
    3298                               setFakeBound(iSequence, ClpSimplexDual::bothFake);