[1370]  1  /* $Id: ClpPrimalColumnSteepest.cpp 2444 20190405 16:11:24Z stefan $ */ 

[2]  2  // Copyright (C) 2002, International Business Machines 

 3  // Corporation and others. All Rights Reserved. 

[1665]  4  // This code is licensed under the terms of the Eclipse Public License (EPL). 

[2]  5  

[63]  6  #include "CoinPragma.hpp" 

[2235]  7  #if ABOCA_LITE 

 8  // 1 is not owner of abcState_ 

 9  #define ABCSTATE_LITE 1 

 10  #endif 

 11  //#define FAKE_CILK 

[2]  12  

 13  #include "ClpSimplex.hpp" 

 14  #include "ClpPrimalColumnSteepest.hpp" 

[50]  15  #include "CoinIndexedVector.hpp" 

[2]  16  #include "ClpFactorization.hpp" 

[162]  17  #include "ClpNonLinearCost.hpp" 

[2]  18  #include "ClpMessage.hpp" 

[2235]  19  #include "ClpEventHandler.hpp" 

 20  #include "ClpPackedMatrix.hpp" 

[2]  21  #include "CoinHelperFunctions.hpp" 

 22  #include <stdio.h> 

[2235]  23  #if ABOCA_LITE 

 24  #undef ALT_UPDATE_WEIGHTS 

 25  #define ALT_UPDATE_WEIGHTS 2 

 26  #endif 

[2385]  27  #if ALT_UPDATE_WEIGHTS == 1 

 28  extern CoinIndexedVector *altVector[3]; 

[2235]  29  #endif 

[2385]  30  static void debug1(int iSequence, double value, double weight) 

[2235]  31  { 

 32  printf("xx %d inf %.20g wt %.20g\n", 

[2385]  33  iSequence, value, weight); 

[2235]  34  } 

[341]  35  //#define CLP_DEBUG 

[2]  36  //############################################################################# 

 37  // Constructors / Destructor / Assignment 

 38  //############################################################################# 

 39  

 40  // 

[1525]  41  // Default Constructor 

[2]  42  // 

[2385]  43  ClpPrimalColumnSteepest::ClpPrimalColumnSteepest(int mode) 

 44  : ClpPrimalColumnPivot() 

 45  , devex_(0.0) 

 46  , weights_(NULL) 

 47  , infeasible_(NULL) 

 48  , alternateWeights_(NULL) 

 49  , savedWeights_(NULL) 

 50  , reference_(NULL) 

 51  , state_(1) 

 52  , mode_(mode) 

 53  , infeasibilitiesState_(0) 

 54  , persistence_(normal) 

 55  , numberSwitched_(0) 

 56  , pivotSequence_(1) 

 57  , savedPivotSequence_(1) 

 58  , savedSequenceOut_(1) 

 59  , sizeFactorization_(0) 

[2]  60  { 

[2385]  61  type_ = 2 + 64 * mode; 

[2]  62  } 

 63  // 

[1525]  64  // Copy constructor 

[2]  65  // 

[2385]  66  ClpPrimalColumnSteepest::ClpPrimalColumnSteepest(const ClpPrimalColumnSteepest &rhs) 

 67  : ClpPrimalColumnPivot(rhs) 

[1525]  68  { 

[2385]  69  state_ = rhs.state_; 

 70  mode_ = rhs.mode_; 

 71  infeasibilitiesState_ = rhs.infeasibilitiesState_; 

 72  persistence_ = rhs.persistence_; 

 73  numberSwitched_ = rhs.numberSwitched_; 

 74  model_ = rhs.model_; 

 75  pivotSequence_ = rhs.pivotSequence_; 

 76  savedPivotSequence_ = rhs.savedPivotSequence_; 

 77  savedSequenceOut_ = rhs.savedSequenceOut_; 

 78  sizeFactorization_ = rhs.sizeFactorization_; 

 79  devex_ = rhs.devex_; 

 80  if ((model_ && model_>whatsChanged() & 1) != 0) { 

 81  if (rhs.infeasible_) { 

 82  infeasible_ = new CoinIndexedVector(rhs.infeasible_); 

 83  } else { 

 84  infeasible_ = NULL; 

 85  } 

 86  reference_ = NULL; 

 87  if (rhs.weights_) { 

 88  assert(model_); 

 89  int number = model_>numberRows() + model_>numberColumns(); 

 90  assert(number == rhs.model_>numberRows() + rhs.model_>numberColumns()); 

 91  weights_ = new double[number]; 

 92  CoinMemcpyN(rhs.weights_, number, weights_); 

 93  savedWeights_ = new double[number]; 

 94  CoinMemcpyN(rhs.savedWeights_, number, savedWeights_); 

 95  if (mode_ != 1) { 

 96  reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5); 

 97  } 

 98  } else { 

 99  weights_ = NULL; 

 100  savedWeights_ = NULL; 

 101  } 

 102  if (rhs.alternateWeights_) { 

 103  alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 

 104  } else { 

 105  alternateWeights_ = NULL; 

 106  } 

 107  } else { 

 108  infeasible_ = NULL; 

 109  reference_ = NULL; 

 110  weights_ = NULL; 

 111  savedWeights_ = NULL; 

 112  alternateWeights_ = NULL; 

 113  } 

[2]  114  } 

 115  

 116  // 

[1525]  117  // Destructor 

[2]  118  // 

[2385]  119  ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest() 

[2]  120  { 

[2385]  121  delete[] weights_; 

 122  delete infeasible_; 

 123  delete alternateWeights_; 

 124  delete[] savedWeights_; 

 125  delete[] reference_; 

[2]  126  } 

 127  

 128  // 

[1525]  129  // Assignment operator 

[2]  130  // 

 131  ClpPrimalColumnSteepest & 

[2385]  132  ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest &rhs) 

[2]  133  { 

[2385]  134  if (this != &rhs) { 

 135  ClpPrimalColumnPivot::operator=(rhs); 

 136  state_ = rhs.state_; 

 137  mode_ = rhs.mode_; 

 138  infeasibilitiesState_ = rhs.infeasibilitiesState_; 

 139  persistence_ = rhs.persistence_; 

 140  numberSwitched_ = rhs.numberSwitched_; 

 141  model_ = rhs.model_; 

 142  pivotSequence_ = rhs.pivotSequence_; 

 143  savedPivotSequence_ = rhs.savedPivotSequence_; 

 144  savedSequenceOut_ = rhs.savedSequenceOut_; 

 145  sizeFactorization_ = rhs.sizeFactorization_; 

 146  devex_ = rhs.devex_; 

 147  delete[] weights_; 

 148  delete[] reference_; 

 149  reference_ = NULL; 

 150  delete infeasible_; 

 151  delete alternateWeights_; 

 152  delete[] savedWeights_; 

 153  savedWeights_ = NULL; 

 154  if (rhs.infeasible_ != NULL) { 

 155  infeasible_ = new CoinIndexedVector(rhs.infeasible_); 

 156  } else { 

 157  infeasible_ = NULL; 

 158  } 

 159  if (rhs.weights_ != NULL) { 

 160  assert(model_); 

 161  int number = model_>numberRows() + model_>numberColumns(); 

 162  assert(number == rhs.model_>numberRows() + rhs.model_>numberColumns()); 

 163  weights_ = new double[number]; 

 164  CoinMemcpyN(rhs.weights_, number, weights_); 

 165  savedWeights_ = new double[number]; 

 166  CoinMemcpyN(rhs.savedWeights_, number, savedWeights_); 

 167  if (mode_ != 1) { 

 168  reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5); 

 169  } 

 170  } else { 

 171  weights_ = NULL; 

 172  } 

 173  if (rhs.alternateWeights_ != NULL) { 

 174  alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 

 175  } else { 

 176  alternateWeights_ = NULL; 

 177  } 

 178  } 

 179  return *this; 

[2]  180  } 

[517]  181  // These have to match ClpPackedMatrix version 

[2]  182  #define TRY_NORM 1.0e4 

 183  #define ADD_ONE 1.0 

[2235]  184  static void 

[2385]  185  pivotColumnBit(clpTempInfo &info) 

[2235]  186  { 

[2385]  187  double *COIN_RESTRICT weights = const_cast< double * >(info.lower); 

 188  const unsigned char *COIN_RESTRICT status = info.status; 

 189  double tolerance = info.tolerance; 

 190  double bestDj = info.primalRatio; 

 191  int bestSequence = 1; 

 192  double *COIN_RESTRICT infeas = const_cast< double * >(info.infeas); 

 193  const int *COIN_RESTRICT start = info.which; 

 194  const int *COIN_RESTRICT index = info.index; 

[2235]  195  for (int i = start[0]; i < start[1]; i++) { 

 196  int iSequence = index[i]; 

 197  double value = infeas[iSequence]; 

 198  double weight = weights[iSequence]; 

 199  if (value > tolerance) { 

 200  if (value > bestDj * weight) { 

[2385]  201  // check flagged variable and correct dj 

 202  if ((status[iSequence] & 64) == 0) { 

 203  bestDj = value / weight; 

 204  bestSequence = iSequence; 

 205  } 

[2235]  206  } 

 207  } 

 208  } 

[2385]  209  info.primalRatio = bestDj; 

 210  info.numberAdded = bestSequence; 

[2235]  211  } 

[2]  212  // Returns pivot column, 1 if none 

[259]  213  /* The Packed CoinIndexedVector updates has cost updates  for normal LP 

[1525]  214  that is just +weight where a feasibility changed. It also has 

[259]  215  reduced cost from last iteration in pivot row*/ 

[2385]  216  int ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector *updates, 

 217  CoinIndexedVector *spareRow1, 

 218  CoinIndexedVector *spareRow2, 

 219  CoinIndexedVector *spareColumn1, 

 220  CoinIndexedVector *spareColumn2) 

[2]  221  { 

[2385]  222  assert(model_); 

 223  if (model_>nonLinearCost()>lookBothWays()  model_>algorithm() == 2) { 

 224  // Do old way 

 225  updates>expand(); 

 226  return pivotColumnOldMethod(updates, spareRow1, spareRow2, 

 227  spareColumn1, spareColumn2); 

 228  } 

 229  int number = 0; 

 230  int *index; 

 231  double tolerance = model_>currentDualTolerance(); 

 232  // we can't really trust infeasibilities if there is dual error 

 233  // this coding has to mimic coding in checkDualSolution 

 234  double error = CoinMin(1.0e2, model_>largestDualError()); 

 235  // allow tolerance at least slightly bigger than standard 

 236  tolerance = tolerance + error; 

 237  int pivotRow = model_>pivotRow(); 

 238  int anyUpdates; 

 239  double *infeas = infeasible_>denseVector(); 

[225]  240  

[2385]  241  // Local copy of mode so can decide what to do 

 242  int switchType; 

 243  if (mode_ == 4) 

 244  switchType = 5  numberSwitched_; 

 245  else if (mode_ >= 10) 

 246  switchType = 3; 

 247  else 

 248  switchType = mode_; 

 249  /* switchType  

[1525]  250  0  all exact devex 

 251  1  all steepest 

 252  2  some exact devex 

 253  3  auto some exact devex 

 254  4  devex 

 255  5  dantzig 

 256  10  can go to minisprint 

 257  */ 

[2385]  258  // Look at gub 

[336]  259  #if 1 

[2385]  260  model_>clpMatrix()>dualExpanded(model_, updates, NULL, 4); 

[336]  261  #else 

[2385]  262  updates>clear(); 

 263  model_>computeDuals(NULL); 

[336]  264  #endif 

[2385]  265  if (updates>getNumElements() > 1) { 

 266  // would have to have two goes for devex, three for steepest 

 267  anyUpdates = 2; 

 268  } else if (updates>getNumElements()) { 

 269  if (updates>getIndices()[0] == pivotRow && fabs(updates>denseVector()[0]) > 1.0e6) { 

 270  // reasonable size 

 271  anyUpdates = 1; 

 272  //if (fabs(model_>dualIn())<1.0e4fabs(fabs(model_>dualIn())fabs(updates>denseVector()[0]))>1.0e5) 

 273  //printf("dualin %g pivot %g\n",model_>dualIn(),updates>denseVector()[0]); 

 274  } else { 

 275  // too small 

 276  anyUpdates = 2; 

 277  } 

 278  } else if (pivotSequence_ >= 0) { 

 279  // just after refactorization 

 280  anyUpdates = 1; 

 281  } else { 

 282  // sub flip  nothing to do 

 283  anyUpdates = 0; 

 284  } 

 285  int sequenceOut = model_>sequenceOut(); 

 286  if (switchType == 5) { 

 287  // If known matrix then we will do partial pricing 

 288  if (model_>clpMatrix()>canDoPartialPricing()) { 

 289  pivotSequence_ = 1; 

 290  pivotRow = 1; 

 291  // See if to switch 

 292  int numberRows = model_>numberRows(); 

 293  int numberWanted = 10; 

 294  int numberColumns = model_>numberColumns(); 

 295  int numberHiddenRows = model_>clpMatrix()>hiddenRows(); 

 296  double ratio = static_cast< double >(sizeFactorization_ + numberHiddenRows) / static_cast< double >(numberRows + 2 * numberHiddenRows); 

 297  // Number of dual infeasibilities at last invert 

 298  int numberDual = model_>numberDualInfeasibilities(); 

 299  int numberLook = CoinMin(numberDual, numberColumns / 10); 

 300  if (ratio < 1.0) { 

 301  numberWanted = 100; 

 302  numberLook /= 20; 

 303  numberWanted = CoinMax(numberWanted, numberLook); 

 304  } else if (ratio < 3.0) { 

 305  numberWanted = 500; 

 306  numberLook /= 15; 

 307  numberWanted = CoinMax(numberWanted, numberLook); 

 308  } else if (ratio < 4.0  mode_ == 5) { 

 309  numberWanted = 1000; 

 310  numberLook /= 10; 

 311  numberWanted = CoinMax(numberWanted, numberLook); 

 312  } else if (mode_ != 5) { 

 313  switchType = 4; 

 314  // initialize 

 315  numberSwitched_++; 

 316  // Make sure will redo 

 317  delete[] weights_; 

 318  weights_ = NULL; 

 319  model_>computeDuals(NULL); 

 320  saveWeights(model_, 4); 

 321  anyUpdates = 0; 

 322  COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio)); 

 323  } 

 324  if (switchType == 5) { 

 325  numberLook *= 5; // needs tuning for gub 

 326  if (model_>numberIterations() % 1000 == 0 && model_>logLevel() > 1) { 

 327  COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n", 

 328  sizeFactorization_, ratio, numberWanted, numberLook)); 

 329  } 

 330  // Update duals and row djs 

 331  // Do partial pricing 

 332  return partialPricing(updates, spareRow2, 

 333  numberWanted, numberLook); 

 334  } 

 335  } 

 336  } 

 337  int bestSequence = 1; 

 338  model_>spareIntArray_[3] = 3; 

 339  if (switchType == 5) { 

 340  if (anyUpdates > 0) { 

 341  justDjs(updates, spareRow2, 

 342  spareColumn1, spareColumn2); 

 343  } 

 344  } else if (anyUpdates == 1) { 

 345  if (switchType < 4) { 

 346  // exact etc when can use dj 

 347  djsAndSteepest(updates, spareRow2, 

 348  spareColumn1, spareColumn2); 

 349  if (model_>spareIntArray_[3] > 2) { 

 350  bestSequence = model_>spareIntArray_[3]; 

 351  infeasibilitiesState_ = 2; 

 352  } else if (model_>spareIntArray_[3] == 2) { 

 353  redoInfeasibilities(); 

 354  } 

 355  } else { 

 356  // devex etc when can use dj 

 357  djsAndDevex(updates, spareRow2, 

 358  spareColumn1, spareColumn2); 

 359  } 

 360  } else if (anyUpdates == 1) { 

 361  if (switchType < 4) { 

 362  // exact etc when djs okay 

 363  justSteepest(updates, spareRow2, 

 364  spareColumn1, spareColumn2); 

 365  } else { 

 366  // devex etc when djs okay 

 367  justDevex(updates, spareRow2, 

 368  spareColumn1, spareColumn2); 

 369  } 

 370  } else if (anyUpdates == 2) { 

 371  if (switchType < 4) { 

 372  // exact etc when have to use pivot 

 373  djsAndSteepest2(updates, spareRow2, 

 374  spareColumn1, spareColumn2); 

 375  } else { 

 376  // devex etc when have to use pivot 

 377  djsAndDevex2(updates, spareRow2, 

 378  spareColumn1, spareColumn2); 

 379  } 

 380  } 

 381  // everything may have been done if vector copy 

 382  if (infeasibilitiesState_ == 2) { 

 383  // all done 

 384  infeasibilitiesState_ = 1; 

 385  model_>clpMatrix()>setSavedBestSequence(bestSequence); 

 386  if (bestSequence >= 0) 

 387  model_>clpMatrix()>setSavedBestDj(model_>djRegion()[bestSequence]); 

 388  assert(sequenceOut != bestSequence); 

 389  return bestSequence; 

 390  } else if (infeasibilitiesState_ == 1) { 

 391  // need to redo 

 392  //infeasibilitiesState_ = 0; 

 393  redoInfeasibilities(); 

 394  } 

 395  

[451]  396  #ifdef CLP_DEBUG 

[2385]  397  alternateWeights_>checkClear(); 

[451]  398  #endif 

[2385]  399  // make sure outgoing from last iteration okay 

 400  if (sequenceOut >= 0) { 

 401  ClpSimplex::Status status = model_>getStatus(sequenceOut); 

 402  double value = model_>reducedCost(sequenceOut); 

[225]  403  

[2385]  404  switch (status) { 

[225]  405  

[2385]  406  case ClpSimplex::basic: 

 407  case ClpSimplex::isFixed: 

 408  break; 

 409  case ClpSimplex::isFree: 

 410  case ClpSimplex::superBasic: 

 411  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 412  // we are going to bias towards free (but only if reasonable) 

 413  value *= FREE_BIAS; 

 414  // store square in list 

 415  if (infeas[sequenceOut]) 

 416  infeas[sequenceOut] = value * value; // already there 

 417  else 

 418  infeasible_>quickAdd(sequenceOut, value * value); 

 419  } else { 

 420  infeasible_>zero(sequenceOut); 

 421  } 

 422  break; 

 423  case ClpSimplex::atUpperBound: 

 424  if (value > tolerance) { 

 425  // store square in list 

 426  if (infeas[sequenceOut]) 

 427  infeas[sequenceOut] = value * value; // already there 

 428  else 

 429  infeasible_>quickAdd(sequenceOut, value * value); 

 430  } else { 

 431  infeasible_>zero(sequenceOut); 

 432  } 

 433  break; 

 434  case ClpSimplex::atLowerBound: 

 435  if (value < tolerance) { 

 436  // store square in list 

 437  if (infeas[sequenceOut]) 

 438  infeas[sequenceOut] = value * value; // already there 

 439  else 

 440  infeasible_>quickAdd(sequenceOut, value * value); 

 441  } else { 

 442  infeasible_>zero(sequenceOut); 

 443  } 

 444  } 

 445  } 

 446  // update of duals finished  now do pricing 

 447  // See what sort of pricing 

 448  int numberWanted = 10; 

 449  number = infeasible_>getNumElements(); 

 450  int numberColumns = model_>numberColumns(); 

 451  if (switchType == 5) { 

 452  pivotSequence_ = 1; 

 453  pivotRow = 1; 

 454  // See if to switch 

 455  int numberRows = model_>numberRows(); 

 456  // ratio is done on number of columns here 

 457  //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns; 

 458  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows); 

 459  //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_>clpMatrix()>getNumElements(); 

 460  if (ratio < 1.0) { 

 461  numberWanted = CoinMax(100, number / 200); 

 462  } else if (ratio < 2.0  1.0) { 

 463  numberWanted = CoinMax(500, number / 40); 

 464  } else if (ratio < 4.0  3.0  mode_ == 5) { 

 465  numberWanted = CoinMax(2000, number / 10); 

 466  numberWanted = CoinMax(numberWanted, numberColumns / 30); 

 467  } else if (mode_ != 5) { 

 468  switchType = 4; 

 469  // initialize 

 470  numberSwitched_++; 

 471  // Make sure will redo 

 472  delete[] weights_; 

 473  weights_ = NULL; 

 474  saveWeights(model_, 4); 

 475  COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio)); 

 476  } 

 477  //if (model_>numberIterations()%1000==0) 

 478  //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted); 

 479  } 

 480  int numberRows = model_>numberRows(); 

 481  // ratio is done on number of rows here 

 482  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows); 

 483  if (switchType == 4) { 

 484  // Still in devex mode 

 485  // Go to steepest if lot of iterations? 

 486  if (ratio < 5.0) { 

 487  numberWanted = CoinMax(2000, number / 10); 

 488  numberWanted = CoinMax(numberWanted, numberColumns / 20); 

 489  } else if (ratio < 7.0) { 

 490  numberWanted = CoinMax(2000, number / 5); 

 491  numberWanted = CoinMax(numberWanted, numberColumns / 10); 

 492  } else { 

 493  // we can zero out 

 494  updates>clear(); 

 495  spareColumn1>clear(); 

 496  switchType = 3; 

 497  // initialize 

 498  pivotSequence_ = 1; 

 499  pivotRow = 1; 

 500  numberSwitched_++; 

 501  // Make sure will redo 

 502  delete[] weights_; 

 503  weights_ = NULL; 

 504  saveWeights(model_, 4); 

 505  COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio)); 

 506  updates>clear(); 

 507  } 

 508  if (model_>numberIterations() % 1000 == 0) 

 509  COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted)); 

 510  } 

 511  if (switchType < 4) { 

 512  if (switchType < 2) { 

 513  numberWanted = COIN_INT_MAX  1; 

 514  } else if (switchType == 2) { 

 515  numberWanted = CoinMax(2000, number / 8); 

 516  } else { 

 517  if (ratio < 1.0) { 

 518  numberWanted = CoinMax(2000, number / 20); 

 519  } else if (ratio < 5.0) { 

 520  numberWanted = CoinMax(2000, number / 10); 

 521  numberWanted = CoinMax(numberWanted, numberColumns / 40); 

 522  } else if (ratio < 10.0) { 

 523  numberWanted = CoinMax(2000, number / 8); 

 524  numberWanted = CoinMax(numberWanted, numberColumns / 20); 

 525  } else { 

 526  ratio = number * (ratio / 80.0); 

 527  if (ratio > number) { 

 528  numberWanted = number + 1; 

 529  } else { 

 530  numberWanted = CoinMax(2000, static_cast< int >(ratio)); 

 531  numberWanted = CoinMax(numberWanted, numberColumns / 10); 

 532  } 

 533  } 

 534  } 

 535  //if (model_>numberIterations()%1000==0) 

 536  //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted, 

 537  //switchType); 

 538  } 

[225]  539  

[2385]  540  double bestDj = 1.0e30; 

 541  bestSequence = 1; 

[2235]  542  #ifdef CLP_USER_DRIVEN 

[2385]  543  // could go parallel? 

 544  model_>eventHandler()>event(ClpEventHandler::beforeChooseIncoming); 

[2235]  545  #endif 

[1525]  546  

[2385]  547  int i, iSequence; 

 548  index = infeasible_>getIndices(); 

 549  number = infeasible_>getNumElements(); 

 550  // Resort infeasible every 100 pivots 

 551  // Not a good idea 

 552  if (0 && model_>factorization()>pivots() > 0 && (model_>factorization()>pivots() % 100) == 0) { 

 553  int nLook = model_>numberRows() + numberColumns; 

 554  number = 0; 

 555  for (i = 0; i < nLook; i++) { 

 556  if (infeas[i]) { 

 557  if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT) 

 558  index[number++] = i; 

 559  else 

 560  infeas[i] = 0.0; 

 561  } 

 562  } 

 563  infeasible_>setNumElements(number); 

 564  } 

 565  if (model_>numberIterations() < model_>lastBadIteration() + 200 && model_>factorization()>pivots() > 10) { 

 566  // we can't really trust infeasibilities if there is dual error 

 567  double checkTolerance = 1.0e8; 

 568  if (model_>largestDualError() > checkTolerance) 

 569  tolerance *= model_>largestDualError() / checkTolerance; 

 570  // But cap 

 571  tolerance = CoinMin(1000.0, tolerance); 

 572  } 

[225]  573  #ifdef CLP_DEBUG 

[2385]  574  if (model_>numberDualInfeasibilities() == 1) 

 575  printf("** %g %g %g %x %x %d\n", tolerance, model_>dualTolerance(), 

 576  model_>largestDualError(), model_, model_>messageHandler(), 

 577  number); 

[225]  578  #endif 

[2385]  579  // stop last one coming immediately 

 580  double saveOutInfeasibility = 0.0; 

 581  if (sequenceOut >= 0) { 

 582  saveOutInfeasibility = infeas[sequenceOut]; 

 583  infeas[sequenceOut] = 0.0; 

 584  } 

 585  if (model_>factorization()>pivots() && model_>numberPrimalInfeasibilities()) 

 586  tolerance = CoinMax(tolerance, 1.0e15 * model_>infeasibilityCost()); 

 587  tolerance *= tolerance; // as we are using squares 

[225]  588  

[2385]  589  int iPass; 

 590  // Setup two passes (unless all) 

 591  if (mode_ > 1) { 

 592  int start[4]; 

 593  start[1] = number; 

 594  start[2] = 0; 

 595  double dstart = static_cast< double >(number) * model_>randomNumberGenerator()>randomDouble(); 

 596  start[0] = static_cast< int >(dstart); 

 597  start[3] = start[0]; 

 598  for (iPass = 0; iPass < 2; iPass++) { 

 599  int end = start[2 * iPass + 1]; 

 600  if (switchType < 5) { 

 601  for (i = start[2 * iPass]; i < end; i++) { 

 602  iSequence = index[i]; 

 603  double value = infeas[iSequence]; 

 604  double weight = weights_[iSequence]; 

 605  //assert (weight); 

 606  if (value > tolerance) { 

 607  //weight=1.0; 

 608  if (value > bestDj * weight) { 

 609  // check flagged variable and correct dj 

 610  if (!model_>flagged(iSequence)) { 

 611  bestDj = value / weight; 

 612  bestSequence = iSequence; 

 613  } else { 

 614  // just to make sure we don't exit before got something 

 615  numberWanted++; 

 616  } 

 617  } 

 618  numberWanted; 

 619  } 

 620  if (!numberWanted) 

 621  break; 

 622  } 

 623  } else { 

 624  // Dantzig 

 625  for (i = start[2 * iPass]; i < end; i++) { 

 626  iSequence = index[i]; 

 627  double value = infeas[iSequence]; 

 628  if (value > tolerance) { 

 629  if (value > bestDj) { 

 630  // check flagged variable and correct dj 

 631  if (!model_>flagged(iSequence)) { 

 632  bestDj = value; 

 633  bestSequence = iSequence; 

 634  } else { 

 635  // just to make sure we don't exit before got something 

 636  numberWanted++; 

 637  } 

 638  } 

 639  numberWanted; 

 640  } 

 641  if (!numberWanted) 

 642  break; 

 643  } 

 644  } 

 645  if (!numberWanted) 

 646  break; 

 647  } 

 648  } else { 

[2235]  649  #if ABOCA_LITE 

[2385]  650  int numberThreads = CoinMax(abcState(), 1); 

[2235]  651  #define ABOCA_LITE_MAX ABOCA_LITE 

 652  #else 

[2385]  653  const int numberThreads = 1; 

[2235]  654  #define ABOCA_LITE_MAX 1 

 655  #endif 

[2385]  656  if (0) { 

 657  int iSequence; 

 658  double value; 

 659  double weight; 

 660  iSequence = 34841; 

 661  value = infeas[iSequence]; 

 662  weight = weights_[iSequence]; 

 663  debug1(iSequence, value, weight); 

 664  iSequence = 34845; 

 665  value = infeas[iSequence]; 

 666  weight = weights_[iSequence]; 

 667  debug1(iSequence, value, weight); 

 668  } 

 669  clpTempInfo info[ABOCA_LITE_MAX]; 

 670  int start_lite[2 * ABOCA_LITE_MAX]; 

 671  int chunk0 = (number + numberThreads  1) / numberThreads; 

 672  int n0 = 0; 

 673  for (int i = 0; i < numberThreads; i++) { 

 674  int *startX = start_lite + 2 * i; 

 675  info[i].primalRatio = bestDj; 

 676  info[i].lower = weights_; 

 677  info[i].infeas = infeas; 

 678  info[i].index = index; 

 679  info[i].status = const_cast< unsigned char * >(model_>statusArray()); 

 680  info[i].which = startX; 

 681  info[i].tolerance = tolerance; 

 682  startX[0] = n0; 

 683  startX[1] = CoinMin(n0 + chunk0, number); 

 684  n0 += chunk0; 

 685  } 

 686  if (numberThreads == 1) { 

 687  pivotColumnBit(info[0]); 

 688  } else { 

 689  for (int i = 0; i < numberThreads; i++) { 

 690  cilk_spawn pivotColumnBit(info[i]); 

 691  } 

 692  cilk_sync; 

 693  } 

 694  for (int i = 0; i < numberThreads; i++) { 

 695  double bestDjX = info[i].primalRatio; 

 696  if (bestDjX > bestDj) { 

 697  bestDj = bestDjX; 

 698  bestSequence = info[i].numberAdded; 

 699  } 

 700  } 

 701  } 

 702  //double largestWeight=0.0; 

 703  //double smallestWeight=1.0e100; 

 704  model_>clpMatrix()>setSavedBestSequence(bestSequence); 

 705  if (bestSequence >= 0) 

 706  model_>clpMatrix()>setSavedBestDj(model_>djRegion()[bestSequence]); 

 707  if (sequenceOut >= 0) { 

 708  infeas[sequenceOut] = saveOutInfeasibility; 

 709  } 

 710  /*if (model_>numberIterations()%100==0) 

[1525]  711  printf("%d best %g\n",bestSequence,bestDj);*/ 

 712  

[2235]  713  #ifdef CLP_USER_DRIVEN 

[2385]  714  // swap when working 

 715  struct { 

 716  int bestSequence; 

 717  double bestDj; 

 718  } tempInfo; 

 719  tempInfo.bestDj = bestDj; 

 720  tempInfo.bestSequence = bestSequence; 

 721  model_>eventHandler()>eventWithInfo(ClpEventHandler::afterChooseIncoming, 

 722  reinterpret_cast< void * >(&tempInfo)); 

[2235]  723  #endif 

[225]  724  #ifndef NDEBUG 

[2385]  725  if (bestSequence >= 0) { 

 726  if (model_>getStatus(bestSequence) == ClpSimplex::atLowerBound) 

 727  assert(model_>reducedCost(bestSequence) < 0.0); 

 728  if (model_>getStatus(bestSequence) == ClpSimplex::atUpperBound) { 

 729  assert(model_>reducedCost(bestSequence) > 0.0); 

 730  } 

 731  } 

[225]  732  #endif 

[2235]  733  #ifdef ALT_UPDATE_WEIGHTSz 

[2385]  734  printf("weights"); 

 735  for (int i = 0; i < model_>numberColumns() + model_>numberRows(); i++) { 

 736  if (model_>getColumnStatus(i) == ClpSimplex::isFixed  model_>getColumnStatus(i) == ClpSimplex::basic) 

 737  weights_[i] = 1.0; 

 738  printf(" %g", weights_[i]); 

 739  if ((i % 10) == 9) 

 740  printf("\n"); 

 741  } 

 742  printf("\n"); 

[2235]  743  #endif 

[1878]  744  #if 0 

 745  for (int i=0;i<numberRows;i++) 

 746  printf("row %d weight %g infeas %g\n",i,weights_[i+numberColumns],infeas[i+numberColumns]); 

 747  for (int i=0;i<numberColumns;i++) 

 748  printf("column %d weight %g infeas %g\n",i,weights_[i],infeas[i]); 

 749  #endif 

[2385]  750  return bestSequence; 

[225]  751  } 

 752  // Just update djs 

[2385]  753  void ClpPrimalColumnSteepest::justDjs(CoinIndexedVector *updates, 

 754  CoinIndexedVector *spareRow2, 

 755  CoinIndexedVector *spareColumn1, 

 756  CoinIndexedVector *spareColumn2) 

[225]  757  { 

[2385]  758  int iSection, j; 

 759  int number = 0; 

 760  int *index; 

 761  double *updateBy; 

 762  double *reducedCost; 

 763  double tolerance = model_>currentDualTolerance(); 

 764  // we can't really trust infeasibilities if there is dual error 

 765  // this coding has to mimic coding in checkDualSolution 

 766  double error = CoinMin(1.0e2, model_>largestDualError()); 

 767  // allow tolerance at least slightly bigger than standard 

 768  tolerance = tolerance + error; 

 769  int pivotRow = model_>pivotRow(); 

 770  double *infeas = infeasible_>denseVector(); 

 771  //updates>scanAndPack(); 

 772  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

[1197]  773  

[2385]  774  // put row of tableau in rowArray and columnArray (packed mode) 

 775  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 776  updates, spareColumn2, spareColumn1); 

 777  // normal 

 778  for (iSection = 0; iSection < 2; iSection++) { 

[1525]  779  

[2385]  780  reducedCost = model_>djRegion(iSection); 

 781  int addSequence; 

 782  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 783  double slack_multiplier; 

[1732]  784  #endif 

[1525]  785  

[2385]  786  if (!iSection) { 

 787  number = updates>getNumElements(); 

 788  index = updates>getIndices(); 

 789  updateBy = updates>denseVector(); 

 790  addSequence = model_>numberColumns(); 

 791  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 792  slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  793  #endif 

[2385]  794  } else { 

 795  number = spareColumn1>getNumElements(); 

 796  index = spareColumn1>getIndices(); 

 797  updateBy = spareColumn1>denseVector(); 

 798  addSequence = 0; 

 799  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 800  slack_multiplier = 1.0; 

[1732]  801  #endif 

[2385]  802  } 

[1525]  803  

[2385]  804  for (j = 0; j < number; j++) { 

 805  int iSequence = index[j]; 

 806  double value = reducedCost[iSequence]; 

 807  value = updateBy[j]; 

 808  updateBy[j] = 0.0; 

 809  reducedCost[iSequence] = value; 

 810  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  811  

[2385]  812  switch (status) { 

[1525]  813  

[2385]  814  case ClpSimplex::basic: 

 815  infeasible_>zero(iSequence + addSequence); 

 816  case ClpSimplex::isFixed: 

 817  break; 

 818  case ClpSimplex::isFree: 

 819  case ClpSimplex::superBasic: 

 820  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 821  // we are going to bias towards free (but only if reasonable) 

 822  value *= FREE_BIAS; 

 823  // store square in list 

 824  if (infeas[iSequence + addSequence]) 

 825  infeas[iSequence + addSequence] = value * value; // already there 

 826  else 

 827  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 828  } else { 

 829  infeasible_>zero(iSequence + addSequence); 

 830  } 

 831  break; 

 832  case ClpSimplex::atUpperBound: 

 833  iSequence += addSequence; 

 834  if (value > tolerance) { 

[1732]  835  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  836  value *= value * slack_multiplier; 

[1732]  837  #else 

[2385]  838  value *= value; 

[1732]  839  #endif 

[2385]  840  // store square in list 

 841  if (infeas[iSequence]) 

 842  infeas[iSequence] = value; // already there 

 843  else 

 844  infeasible_>quickAdd(iSequence, value); 

 845  } else { 

 846  infeasible_>zero(iSequence); 

 847  } 

 848  break; 

 849  case ClpSimplex::atLowerBound: 

 850  iSequence += addSequence; 

 851  if (value < tolerance) { 

[1732]  852  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  853  value *= value * slack_multiplier; 

[1732]  854  #else 

[2385]  855  value *= value; 

[1732]  856  #endif 

[2385]  857  // store square in list 

 858  if (infeas[iSequence]) 

 859  infeas[iSequence] = value; // already there 

 860  else 

 861  infeasible_>quickAdd(iSequence, value); 

 862  } else { 

 863  infeasible_>zero(iSequence); 

 864  } 

 865  } 

 866  } 

 867  } 

 868  updates>setNumElements(0); 

 869  spareColumn1>setNumElements(0); 

 870  if (pivotRow >= 0) { 

 871  // make sure infeasibility on incoming is 0.0 

 872  int sequenceIn = model_>sequenceIn(); 

 873  infeasible_>zero(sequenceIn); 

 874  } 

[225]  875  } 

 876  // Update djs, weights for Devex 

[2385]  877  void ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector *updates, 

 878  CoinIndexedVector *spareRow2, 

 879  CoinIndexedVector *spareColumn1, 

 880  CoinIndexedVector *spareColumn2) 

[225]  881  { 

[2385]  882  int j; 

 883  int number = 0; 

 884  int *index; 

 885  double *updateBy; 

 886  double *reducedCost; 

 887  double tolerance = model_>currentDualTolerance(); 

 888  // we can't really trust infeasibilities if there is dual error 

 889  // this coding has to mimic coding in checkDualSolution 

 890  double error = CoinMin(1.0e2, model_>largestDualError()); 

 891  // allow tolerance at least slightly bigger than standard 

 892  tolerance = tolerance + error; 

 893  // for weights update we use pivotSequence 

 894  // unset in case sub flip 

 895  assert(pivotSequence_ >= 0); 

 896  assert(model_>pivotVariable()[pivotSequence_] == model_>sequenceIn()); 

 897  pivotSequence_ = 1; 

 898  double *infeas = infeasible_>denseVector(); 

 899  //updates>scanAndPack(); 

 900  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 901  // and we can see if reference 

 902  //double referenceIn = 0.0; 

 903  int sequenceIn = model_>sequenceIn(); 

 904  //if (mode_ != 1 && reference(sequenceIn)) 

 905  // referenceIn = 1.0; 

 906  // save outgoing weight round update 

 907  double outgoingWeight = 0.0; 

 908  int sequenceOut = model_>sequenceOut(); 

 909  if (sequenceOut >= 0) 

 910  outgoingWeight = weights_[sequenceOut]; 

[225]  911  

[2385]  912  double scaleFactor = 1.0 / updates>denseVector()[0]; // as formula is with 1.0 

 913  // put row of tableau in rowArray and columnArray (packed mode) 

 914  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 915  updates, spareColumn2, spareColumn1); 

 916  // update weights 

 917  double *weight; 

 918  int numberColumns = model_>numberColumns(); 

 919  // rows 

 920  reducedCost = model_>djRegion(0); 

 921  int addSequence = model_>numberColumns(); 

 922  ; 

[1525]  923  

[2385]  924  number = updates>getNumElements(); 

 925  index = updates>getIndices(); 

 926  updateBy = updates>denseVector(); 

 927  weight = weights_ + numberColumns; 

 928  // Devex 

 929  for (j = 0; j < number; j++) { 

 930  double thisWeight; 

 931  double pivot; 

 932  double value3; 

 933  int iSequence = index[j]; 

 934  double value = reducedCost[iSequence]; 

 935  double value2 = updateBy[j]; 

 936  updateBy[j] = 0.0; 

 937  value = value2; 

 938  reducedCost[iSequence] = value; 

 939  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  940  

[2385]  941  switch (status) { 

[1525]  942  

[2385]  943  case ClpSimplex::basic: 

 944  infeasible_>zero(iSequence + addSequence); 

 945  case ClpSimplex::isFixed: 

 946  break; 

 947  case ClpSimplex::isFree: 

 948  case ClpSimplex::superBasic: 

 949  thisWeight = weight[iSequence]; 

 950  // row has 1 

 951  pivot = value2 * scaleFactor; 

 952  value3 = pivot * pivot * devex_; 

 953  if (reference(iSequence + numberColumns)) 

 954  value3 += 1.0; 

 955  weight[iSequence] = CoinMax(0.99 * thisWeight, value3); 

 956  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 957  // we are going to bias towards free (but only if reasonable) 

 958  value *= FREE_BIAS; 

 959  // store square in list 

 960  if (infeas[iSequence + addSequence]) 

 961  infeas[iSequence + addSequence] = value * value; // already there 

 962  else 

 963  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 964  } else { 

 965  infeasible_>zero(iSequence + addSequence); 

 966  } 

 967  break; 

 968  case ClpSimplex::atUpperBound: 

 969  thisWeight = weight[iSequence]; 

 970  // row has 1 

 971  pivot = value2 * scaleFactor; 

 972  value3 = pivot * pivot * devex_; 

 973  if (reference(iSequence + numberColumns)) 

 974  value3 += 1.0; 

 975  weight[iSequence] = CoinMax(0.99 * thisWeight, value3); 

 976  iSequence += addSequence; 

 977  if (value > tolerance) { 

 978  // store square in list 

[1732]  979  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  980  value *= value * CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  981  #else 

[2385]  982  value *= value; 

[1732]  983  #endif 

[2385]  984  if (infeas[iSequence]) 

 985  infeas[iSequence] = value; // already there 

 986  else 

 987  infeasible_>quickAdd(iSequence, value); 

 988  } else { 

 989  infeasible_>zero(iSequence); 

 990  } 

 991  break; 

 992  case ClpSimplex::atLowerBound: 

 993  thisWeight = weight[iSequence]; 

 994  // row has 1 

 995  pivot = value2 * scaleFactor; 

 996  value3 = pivot * pivot * devex_; 

 997  if (reference(iSequence + numberColumns)) 

 998  value3 += 1.0; 

 999  weight[iSequence] = CoinMax(0.99 * thisWeight, value3); 

 1000  iSequence += addSequence; 

 1001  if (value < tolerance) { 

 1002  // store square in list 

[1732]  1003  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1004  value *= value * CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  1005  #else 

[2385]  1006  value *= value; 

[1732]  1007  #endif 

[2385]  1008  if (infeas[iSequence]) 

 1009  infeas[iSequence] = value; // already there 

 1010  else 

 1011  infeasible_>quickAdd(iSequence, value); 

 1012  } else { 

 1013  infeasible_>zero(iSequence); 

 1014  } 

 1015  } 

 1016  } 

[1525]  1017  

[2385]  1018  // columns 

 1019  weight = weights_; 

[1525]  1020  

[2385]  1021  scaleFactor = scaleFactor; 

 1022  reducedCost = model_>djRegion(1); 

 1023  number = spareColumn1>getNumElements(); 

 1024  index = spareColumn1>getIndices(); 

 1025  updateBy = spareColumn1>denseVector(); 

[1525]  1026  

[2385]  1027  // Devex 

[1525]  1028  

[2385]  1029  for (j = 0; j < number; j++) { 

 1030  double thisWeight; 

 1031  double pivot; 

 1032  double value3; 

 1033  int iSequence = index[j]; 

 1034  double value = reducedCost[iSequence]; 

 1035  double value2 = updateBy[j]; 

 1036  value = value2; 

 1037  updateBy[j] = 0.0; 

 1038  reducedCost[iSequence] = value; 

 1039  ClpSimplex::Status status = model_>getStatus(iSequence); 

[1525]  1040  

[2385]  1041  switch (status) { 

[1525]  1042  

[2385]  1043  case ClpSimplex::basic: 

 1044  infeasible_>zero(iSequence); 

 1045  case ClpSimplex::isFixed: 

 1046  break; 

 1047  case ClpSimplex::isFree: 

 1048  case ClpSimplex::superBasic: 

 1049  thisWeight = weight[iSequence]; 

 1050  // row has 1 

 1051  pivot = value2 * scaleFactor; 

 1052  value3 = pivot * pivot * devex_; 

 1053  if (reference(iSequence)) 

 1054  value3 += 1.0; 

 1055  weight[iSequence] = CoinMax(0.99 * thisWeight, value3); 

 1056  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 1057  // we are going to bias towards free (but only if reasonable) 

 1058  value *= FREE_BIAS; 

 1059  // store square in list 

 1060  if (infeas[iSequence]) 

 1061  infeas[iSequence] = value * value; // already there 

 1062  else 

 1063  infeasible_>quickAdd(iSequence, value * value); 

 1064  } else { 

 1065  infeasible_>zero(iSequence); 

 1066  } 

 1067  break; 

 1068  case ClpSimplex::atUpperBound: 

 1069  thisWeight = weight[iSequence]; 

 1070  // row has 1 

 1071  pivot = value2 * scaleFactor; 

 1072  value3 = pivot * pivot * devex_; 

 1073  if (reference(iSequence)) 

 1074  value3 += 1.0; 

 1075  weight[iSequence] = CoinMax(0.99 * thisWeight, value3); 

 1076  if (value > tolerance) { 

 1077  // store square in list 

 1078  if (infeas[iSequence]) 

 1079  infeas[iSequence] = value * value; // already there 

 1080  else 

 1081  infeasible_>quickAdd(iSequence, value * value); 

 1082  } else { 

 1083  infeasible_>zero(iSequence); 

 1084  } 

 1085  break; 

 1086  case ClpSimplex::atLowerBound: 

 1087  thisWeight = weight[iSequence]; 

 1088  // row has 1 

 1089  pivot = value2 * scaleFactor; 

 1090  value3 = pivot * pivot * devex_; 

 1091  if (reference(iSequence)) 

 1092  value3 += 1.0; 

 1093  weight[iSequence] = CoinMax(0.99 * thisWeight, value3); 

 1094  if (value < tolerance) { 

 1095  // store square in list 

 1096  if (infeas[iSequence]) 

 1097  infeas[iSequence] = value * value; // already there 

 1098  else 

 1099  infeasible_>quickAdd(iSequence, value * value); 

 1100  } else { 

 1101  infeasible_>zero(iSequence); 

 1102  } 

 1103  } 

 1104  } 

 1105  // restore outgoing weight 

 1106  if (sequenceOut >= 0) 

 1107  weights_[sequenceOut] = outgoingWeight; 

 1108  // make sure infeasibility on incoming is 0.0 

 1109  infeasible_>zero(sequenceIn); 

 1110  spareRow2>setNumElements(0); 

 1111  //#define SOME_DEBUG_1 

[225]  1112  #ifdef SOME_DEBUG_1 

[2385]  1113  // check for accuracy 

 1114  int iCheck = 892; 

 1115  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 1116  int numberRows = model_>numberRows(); 

 1117  //int numberColumns = model_>numberColumns(); 

 1118  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 1119  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 1120  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 1121  } 

[225]  1122  #endif 

[2385]  1123  updates>setNumElements(0); 

 1124  spareColumn1>setNumElements(0); 

[225]  1125  } 

 1126  // Update djs, weights for Steepest 

[2385]  1127  void ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector *updates, 

 1128  CoinIndexedVector *spareRow2, 

 1129  CoinIndexedVector *spareColumn1, 

 1130  CoinIndexedVector *spareColumn2) 

[225]  1131  { 

[2385]  1132  int j; 

 1133  int number = 0; 

 1134  int *index; 

 1135  double *updateBy; 

 1136  double *reducedCost; 

 1137  double tolerance = model_>currentDualTolerance(); 

 1138  // we can't really trust infeasibilities if there is dual error 

 1139  // this coding has to mimic coding in checkDualSolution 

 1140  double error = CoinMin(1.0e2, model_>largestDualError()); 

 1141  // allow tolerance at least slightly bigger than standard 

 1142  tolerance = tolerance + error; 

 1143  // for weights update we use pivotSequence 

 1144  // unset in case sub flip 

 1145  assert(pivotSequence_ >= 0); 

 1146  assert(model_>pivotVariable()[pivotSequence_] == model_>sequenceIn()); 

 1147  double *infeas = infeasible_>denseVector(); 

 1148  double scaleFactor = 1.0 / updates>denseVector()[0]; // as formula is with 1.0 

 1149  assert(updates>getIndices()[0] == pivotSequence_); 

 1150  pivotSequence_ = 1; 

[2235]  1151  #if 1 

[2385]  1152  //updates>scanAndPack(); 

 1153  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 1154  //alternateWeights_>scanAndPack(); 

 1155  #if ALT_UPDATE_WEIGHTS != 2 

 1156  model_>factorization()>updateColumnTranspose(spareRow2, 

 1157  alternateWeights_); 

 1158  #elif ALT_UPDATE_WEIGHTS == 1 

 1159  if (altVector[1]) { 

 1160  int numberRows = model_>numberRows(); 

 1161  double *work1 = altVector[1]>denseVector(); 

 1162  double *worka = alternateWeights_>denseVector(); 

 1163  int iRow = 1; 

 1164  double diff = 1.0e8; 

 1165  for (int i = 0; i < numberRows; i++) { 

 1166  double dd = CoinMax(fabs(work1[i]), fabs(worka[i])); 

 1167  double d = fabs(work1[i]  worka[i]); 

 1168  if (dd > 1.0e6 && d > diff * dd) { 

 1169  diff = d / dd; 

 1170  iRow = i; 

 1171  } 

 1172  } 

 1173  if (iRow >= 0) 

 1174  printf("largest difference of %g (%g,%g) on row %d\n", 

 1175  diff, work1[iRow], worka[iRow], iRow); 

 1176  } 

[2235]  1177  #endif 

 1178  #else 

[2385]  1179  model_>factorization()>updateTwoColumnsTranspose(spareRow2, updates, alternateWeights_); 

[2235]  1180  #endif 

[2385]  1181  // and we can see if reference 

 1182  int sequenceIn = model_>sequenceIn(); 

 1183  double referenceIn; 

 1184  if (mode_ != 1) { 

 1185  if (reference(sequenceIn)) 

 1186  referenceIn = 1.0; 

 1187  else 

 1188  referenceIn = 0.0; 

 1189  } else { 

 1190  referenceIn = 1.0; 

 1191  } 

 1192  // save outgoing weight round update 

 1193  double outgoingWeight = 0.0; 

 1194  int sequenceOut = model_>sequenceOut(); 

 1195  if (sequenceOut >= 0) 

 1196  outgoingWeight = weights_[sequenceOut]; 

 1197  // update row weights here so we can scale alternateWeights_ 

 1198  // update weights 

 1199  double *weight; 

 1200  double *other = alternateWeights_>denseVector(); 

 1201  int numberColumns = model_>numberColumns(); 

 1202  // if if (model_>clpMatrix()>canCombine(model_, pi1)) no need to index 

 1203  // rows 

 1204  reducedCost = model_>djRegion(0); 

 1205  int addSequence = model_>numberColumns(); 

 1206  ; 

[225]  1207  

[2385]  1208  number = updates>getNumElements(); 

 1209  index = updates>getIndices(); 

 1210  updateBy = updates>denseVector(); 

 1211  weight = weights_ + numberColumns; 

[2235]  1212  #ifdef CLP_USER_DRIVEN 

[2385]  1213  model_>eventHandler()>eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates); 

[2235]  1214  #endif 

[1525]  1215  

[2385]  1216  for (j = 0; j < number; j++) { 

 1217  double thisWeight; 

 1218  double pivot; 

 1219  double modification; 

 1220  double pivotSquared; 

 1221  int iSequence = index[j]; 

 1222  double value2 = updateBy[j]; 

 1223  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

 1224  double value; 

[1525]  1225  

[2385]  1226  switch (status) { 

[1525]  1227  

[2385]  1228  case ClpSimplex::basic: 

 1229  infeasible_>zero(iSequence + addSequence); 

 1230  reducedCost[iSequence] = 0.0; 

 1231  case ClpSimplex::isFixed: 

 1232  break; 

 1233  case ClpSimplex::isFree: 

 1234  case ClpSimplex::superBasic: 

 1235  value = reducedCost[iSequence]  value2; 

 1236  modification = other[iSequence]; 

 1237  thisWeight = weight[iSequence]; 

 1238  // row has 1 

 1239  pivot = value2 * scaleFactor; 

 1240  pivotSquared = pivot * pivot; 

[1525]  1241  

[2385]  1242  thisWeight += pivotSquared * devex_ + pivot * modification; 

 1243  reducedCost[iSequence] = value; 

 1244  if (thisWeight < TRY_NORM) { 

 1245  if (mode_ == 1) { 

 1246  // steepest 

 1247  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 1248  } else { 

 1249  // exact 

 1250  thisWeight = referenceIn * pivotSquared; 

 1251  if (reference(iSequence + numberColumns)) 

 1252  thisWeight += 1.0; 

 1253  thisWeight = CoinMax(thisWeight, TRY_NORM); 

 1254  } 

 1255  } 

 1256  weight[iSequence] = thisWeight; 

 1257  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 1258  // we are going to bias towards free (but only if reasonable) 

 1259  value *= FREE_BIAS; 

 1260  // store square in list 

 1261  if (infeas[iSequence + addSequence]) 

 1262  infeas[iSequence + addSequence] = value * value; // already there 

 1263  else 

 1264  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 1265  } else { 

 1266  infeasible_>zero(iSequence + addSequence); 

 1267  } 

 1268  break; 

 1269  case ClpSimplex::atUpperBound: 

 1270  value = reducedCost[iSequence]  value2; 

 1271  modification = other[iSequence]; 

 1272  thisWeight = weight[iSequence]; 

 1273  // row has 1 

 1274  pivot = value2 * scaleFactor; 

 1275  pivotSquared = pivot * pivot; 

[1525]  1276  

[2385]  1277  thisWeight += pivotSquared * devex_ + pivot * modification; 

 1278  reducedCost[iSequence] = value; 

 1279  if (thisWeight < TRY_NORM) { 

 1280  if (mode_ == 1) { 

 1281  // steepest 

 1282  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 1283  } else { 

 1284  // exact 

 1285  thisWeight = referenceIn * pivotSquared; 

 1286  if (reference(iSequence + numberColumns)) 

 1287  thisWeight += 1.0; 

 1288  thisWeight = CoinMax(thisWeight, TRY_NORM); 

 1289  } 

 1290  } 

 1291  weight[iSequence] = thisWeight; 

 1292  iSequence += addSequence; 

 1293  if (value > tolerance) { 

 1294  // store square in list 

[1732]  1295  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1296  value *= value * CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  1297  #else 

[2385]  1298  value *= value; 

[1732]  1299  #endif 

[2385]  1300  if (infeas[iSequence]) 

 1301  infeas[iSequence] = value; // already there 

 1302  else 

 1303  infeasible_>quickAdd(iSequence, value); 

 1304  } else { 

 1305  infeasible_>zero(iSequence); 

 1306  } 

 1307  break; 

 1308  case ClpSimplex::atLowerBound: 

 1309  value = reducedCost[iSequence]  value2; 

 1310  modification = other[iSequence]; 

 1311  thisWeight = weight[iSequence]; 

 1312  // row has 1 

 1313  pivot = value2 * scaleFactor; 

 1314  pivotSquared = pivot * pivot; 

[1525]  1315  

[2385]  1316  thisWeight += pivotSquared * devex_ + pivot * modification; 

 1317  reducedCost[iSequence] = value; 

 1318  if (thisWeight < TRY_NORM) { 

 1319  if (mode_ == 1) { 

 1320  // steepest 

 1321  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 1322  } else { 

 1323  // exact 

 1324  thisWeight = referenceIn * pivotSquared; 

 1325  if (reference(iSequence + numberColumns)) 

 1326  thisWeight += 1.0; 

 1327  thisWeight = CoinMax(thisWeight, TRY_NORM); 

 1328  } 

 1329  } 

 1330  weight[iSequence] = thisWeight; 

 1331  iSequence += addSequence; 

 1332  if (value < tolerance) { 

 1333  // store square in list 

[1732]  1334  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1335  value *= value * CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  1336  #else 

[2385]  1337  value *= value; 

[1732]  1338  #endif 

[2385]  1339  if (infeas[iSequence]) 

 1340  infeas[iSequence] = value; // already there 

 1341  else 

 1342  infeasible_>quickAdd(iSequence, value); 

 1343  } else { 

 1344  infeasible_>zero(iSequence); 

 1345  } 

 1346  } 

 1347  } 

[2235]  1348  #ifdef CLP_USER_DRIVEN 

[2385]  1349  // could go parallel? 

 1350  //model_>eventHandler()>eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates); 

[2235]  1351  #endif 

[2385]  1352  // put row of tableau in rowArray and columnArray (packed) 

 1353  // get subset which have nonzero tableau elements 

 1354  int returnCode = transposeTimes2(updates, spareColumn1, 

 1355  alternateWeights_, spareColumn2, spareRow2, 

 1356  scaleFactor); 

 1357  // zero updateBy 

 1358  CoinZeroN(updateBy, number); 

 1359  alternateWeights_>clear(); 

 1360  // columns 

 1361  assert(scaleFactor); 

 1362  number = spareColumn1>getNumElements(); 

 1363  index = spareColumn1>getIndices(); 

 1364  updateBy = spareColumn1>denseVector(); 

 1365  if (returnCode != 2 && infeasibilitiesState_) { 

 1366  //spareColumn1>clear(); 

 1367  redoInfeasibilities(); 

 1368  } 

 1369  if (returnCode == 1) { 

 1370  // most work already done 

 1371  for (j = 0; j < number; j++) { 

 1372  int iSequence = index[j]; 

 1373  double value = updateBy[j]; 

 1374  if (value) { 

 1375  updateBy[j] = 0.0; 

 1376  infeasible_>quickAdd(iSequence, value); 

 1377  } else { 

 1378  infeasible_>zero(iSequence); 

 1379  } 

 1380  } 

 1381  } else if (returnCode == 0) { 

 1382  reducedCost = model_>djRegion(1); 

[1525]  1383  

[2385]  1384  for (j = 0; j < number; j++) { 

 1385  int iSequence = index[j]; 

 1386  double value = reducedCost[iSequence]; 

 1387  double value2 = updateBy[j]; 

 1388  updateBy[j] = 0.0; 

 1389  value = value2; 

 1390  reducedCost[iSequence] = value; 

 1391  ClpSimplex::Status status = model_>getStatus(iSequence); 

 1392  

 1393  switch (status) { 

 1394  

 1395  case ClpSimplex::basic: 

 1396  case ClpSimplex::isFixed: 

 1397  break; 

 1398  case ClpSimplex::isFree: 

 1399  case ClpSimplex::superBasic: 

 1400  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 1401  // we are going to bias towards free (but only if reasonable) 

 1402  value *= FREE_BIAS; 

 1403  // store square in list 

 1404  if (infeas[iSequence]) 

 1405  infeas[iSequence] = value * value; // already there 

 1406  else 

 1407  infeasible_>quickAdd(iSequence, value * value); 

 1408  } else { 

 1409  infeasible_>zero(iSequence); 

 1410  } 

 1411  break; 

 1412  case ClpSimplex::atUpperBound: 

 1413  if (value > tolerance) { 

 1414  // store square in list 

 1415  if (infeas[iSequence]) 

 1416  infeas[iSequence] = value * value; // already there 

 1417  else 

 1418  infeasible_>quickAdd(iSequence, value * value); 

 1419  } else { 

 1420  infeasible_>zero(iSequence); 

 1421  } 

 1422  break; 

 1423  case ClpSimplex::atLowerBound: 

 1424  if (value < tolerance) { 

 1425  // store square in list 

 1426  if (infeas[iSequence]) 

 1427  infeas[iSequence] = value * value; // already there 

 1428  else 

 1429  infeasible_>quickAdd(iSequence, value * value); 

 1430  } else { 

 1431  infeasible_>zero(iSequence); 

 1432  } 

 1433  } 

 1434  } 

 1435  } else { 

 1436  //assert(!number); 

 1437  } 

 1438  // restore outgoing weight 

 1439  if (sequenceOut >= 0) { 

 1440  //#define GROW_REFERENCE 

[2235]  1441  #ifdef GROW_REFERENCE 

[2385]  1442  if (!reference(sequenceOut)) { 

 1443  outgoingWeight += 1.0; 

 1444  setReference(sequenceOut, true); 

 1445  } 

[2235]  1446  #endif 

[2385]  1447  weights_[sequenceOut] = outgoingWeight; 

 1448  //if (model_>getStatus(sequenceOut) != ClpSimplex::basic && 

 1449  // model_>getStatus(sequenceOut) != ClpSimplex::isFixed) 

 1450  //checkAccuracy(sequenceOut, 1.0e1, updates, spareRow2); 

 1451  } 

 1452  // make sure infeasibility on incoming is 0.0 

 1453  infeasible_>zero(sequenceIn); 

 1454  spareColumn2>setNumElements(0); 

 1455  //#define SOME_DEBUG_1 

[225]  1456  #ifdef SOME_DEBUG_1 

[2385]  1457  // check for accuracy 

 1458  int iCheck = 892; 

 1459  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 1460  int numberRows = model_>numberRows(); 

 1461  //int numberColumns = model_>numberColumns(); 

 1462  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 1463  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 1464  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 1465  } 

[225]  1466  #endif 

[2385]  1467  updates>setNumElements(0); 

 1468  spareColumn1>setNumElements(0); 

[225]  1469  } 

 1470  // Update djs, weights for Devex 

[2385]  1471  void ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector *updates, 

 1472  CoinIndexedVector *spareRow2, 

 1473  CoinIndexedVector *spareColumn1, 

 1474  CoinIndexedVector *spareColumn2) 

[225]  1475  { 

[2385]  1476  int iSection, j; 

 1477  int number = 0; 

 1478  int *index; 

 1479  double *updateBy; 

 1480  double *reducedCost; 

 1481  // dj could be very small (or even zero  take care) 

 1482  double dj = model_>dualIn(); 

 1483  double tolerance = model_>currentDualTolerance(); 

 1484  // we can't really trust infeasibilities if there is dual error 

 1485  // this coding has to mimic coding in checkDualSolution 

 1486  double error = CoinMin(1.0e2, model_>largestDualError()); 

 1487  // allow tolerance at least slightly bigger than standard 

 1488  tolerance = tolerance + error; 

 1489  int pivotRow = model_>pivotRow(); 

 1490  double *infeas = infeasible_>denseVector(); 

 1491  //updates>scanAndPack(); 

 1492  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

[1197]  1493  

[2385]  1494  // put row of tableau in rowArray and columnArray 

 1495  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 1496  updates, spareColumn2, spareColumn1); 

 1497  // normal 

 1498  for (iSection = 0; iSection < 2; iSection++) { 

[1525]  1499  

[2385]  1500  reducedCost = model_>djRegion(iSection); 

 1501  int addSequence; 

 1502  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 1503  double slack_multiplier; 

[1732]  1504  #endif 

[1525]  1505  

[2385]  1506  if (!iSection) { 

 1507  number = updates>getNumElements(); 

 1508  index = updates>getIndices(); 

 1509  updateBy = updates>denseVector(); 

 1510  addSequence = model_>numberColumns(); 

 1511  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 1512  slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  1513  #endif 

[2385]  1514  } else { 

 1515  number = spareColumn1>getNumElements(); 

 1516  index = spareColumn1>getIndices(); 

 1517  updateBy = spareColumn1>denseVector(); 

 1518  addSequence = 0; 

 1519  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 1520  slack_multiplier = 1.0; 

[1732]  1521  #endif 

[2385]  1522  } 

[1525]  1523  

[2385]  1524  for (j = 0; j < number; j++) { 

 1525  int iSequence = index[j]; 

 1526  double value = reducedCost[iSequence]; 

 1527  value = updateBy[j]; 

 1528  updateBy[j] = 0.0; 

 1529  reducedCost[iSequence] = value; 

 1530  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  1531  

[2385]  1532  switch (status) { 

[1525]  1533  

[2385]  1534  case ClpSimplex::basic: 

 1535  infeasible_>zero(iSequence + addSequence); 

 1536  case ClpSimplex::isFixed: 

 1537  break; 

 1538  case ClpSimplex::isFree: 

 1539  case ClpSimplex::superBasic: 

 1540  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 1541  // we are going to bias towards free (but only if reasonable) 

 1542  value *= FREE_BIAS; 

 1543  // store square in list 

 1544  if (infeas[iSequence + addSequence]) 

 1545  infeas[iSequence + addSequence] = value * value; // already there 

 1546  else 

 1547  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 1548  } else { 

 1549  infeasible_>zero(iSequence + addSequence); 

 1550  } 

 1551  break; 

 1552  case ClpSimplex::atUpperBound: 

 1553  iSequence += addSequence; 

 1554  if (value > tolerance) { 

[1732]  1555  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1556  value *= value * slack_multiplier; 

[1732]  1557  #else 

[2385]  1558  value *= value; 

[1732]  1559  #endif 

[2385]  1560  // store square in list 

 1561  if (infeas[iSequence]) 

 1562  infeas[iSequence] = value; // already there 

 1563  else 

 1564  infeasible_>quickAdd(iSequence, value); 

 1565  } else { 

 1566  infeasible_>zero(iSequence); 

 1567  } 

 1568  break; 

 1569  case ClpSimplex::atLowerBound: 

 1570  iSequence += addSequence; 

 1571  if (value < tolerance) { 

[1732]  1572  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1573  value *= value * slack_multiplier; 

[1732]  1574  #else 

[2385]  1575  value *= value; 

[1732]  1576  #endif 

[2385]  1577  // store square in list 

 1578  if (infeas[iSequence]) 

 1579  infeas[iSequence] = value; // already there 

 1580  else 

 1581  infeasible_>quickAdd(iSequence, value); 

 1582  } else { 

 1583  infeasible_>zero(iSequence); 

 1584  } 

 1585  } 

 1586  } 

 1587  } 

 1588  // They are empty 

 1589  updates>setNumElements(0); 

 1590  spareColumn1>setNumElements(0); 

 1591  // make sure infeasibility on incoming is 0.0 

 1592  int sequenceIn = model_>sequenceIn(); 

 1593  infeasible_>zero(sequenceIn); 

 1594  // for weights update we use pivotSequence 

 1595  if (pivotSequence_ >= 0) { 

 1596  pivotRow = pivotSequence_; 

 1597  // unset in case sub flip 

 1598  pivotSequence_ = 1; 

 1599  // make sure infeasibility on incoming is 0.0 

 1600  const int *pivotVariable = model_>pivotVariable(); 

 1601  sequenceIn = pivotVariable[pivotRow]; 

 1602  infeasible_>zero(sequenceIn); 

 1603  // and we can see if reference 

 1604  //double referenceIn = 0.0; 

 1605  //if (mode_ != 1 && reference(sequenceIn)) 

 1606  // referenceIn = 1.0; 

 1607  // save outgoing weight round update 

 1608  double outgoingWeight = 0.0; 

 1609  int sequenceOut = model_>sequenceOut(); 

 1610  if (sequenceOut >= 0) 

 1611  outgoingWeight = weights_[sequenceOut]; 

 1612  // update weights 

 1613  updates>setNumElements(0); 

 1614  spareColumn1>setNumElements(0); 

 1615  // might as well set dj to 1 

 1616  dj = 1.0; 

 1617  updates>insert(pivotRow, dj); 

 1618  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 1619  // put row of tableau in rowArray and columnArray 

 1620  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 1621  updates, spareColumn2, spareColumn1); 

 1622  double *weight; 

 1623  int numberColumns = model_>numberColumns(); 

 1624  // rows 

 1625  number = updates>getNumElements(); 

 1626  index = updates>getIndices(); 

 1627  updateBy = updates>denseVector(); 

 1628  weight = weights_ + numberColumns; 

[1525]  1629  

[2385]  1630  assert(devex_ > 0.0); 

 1631  for (j = 0; j < number; j++) { 

 1632  int iSequence = index[j]; 

 1633  double thisWeight = weight[iSequence]; 

 1634  // row has 1 

 1635  double pivot = updateBy[iSequence]; 

 1636  updateBy[iSequence] = 0.0; 

 1637  double value = pivot * pivot * devex_; 

 1638  if (reference(iSequence + numberColumns)) 

 1639  value += 1.0; 

 1640  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 1641  } 

[1525]  1642  

[2385]  1643  // columns 

 1644  weight = weights_; 

[1525]  1645  

[2385]  1646  number = spareColumn1>getNumElements(); 

 1647  index = spareColumn1>getIndices(); 

 1648  updateBy = spareColumn1>denseVector(); 

 1649  for (j = 0; j < number; j++) { 

 1650  int iSequence = index[j]; 

 1651  double thisWeight = weight[iSequence]; 

 1652  // row has 1 

 1653  double pivot = updateBy[iSequence]; 

 1654  updateBy[iSequence] = 0.0; 

 1655  double value = pivot * pivot * devex_; 

 1656  if (reference(iSequence)) 

 1657  value += 1.0; 

 1658  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 1659  } 

 1660  // restore outgoing weight 

 1661  if (sequenceOut >= 0) 

 1662  weights_[sequenceOut] = outgoingWeight; 

 1663  spareColumn2>setNumElements(0); 

 1664  //#define SOME_DEBUG_1 

[225]  1665  #ifdef SOME_DEBUG_1 

[2385]  1666  // check for accuracy 

 1667  int iCheck = 892; 

 1668  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 1669  int numberRows = model_>numberRows(); 

 1670  //int numberColumns = model_>numberColumns(); 

 1671  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 1672  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 1673  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 1674  } 

[225]  1675  #endif 

[2385]  1676  updates>setNumElements(0); 

 1677  spareColumn1>setNumElements(0); 

 1678  } 

[225]  1679  } 

 1680  // Update djs, weights for Steepest 

[2385]  1681  void ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector *updates, 

 1682  CoinIndexedVector *spareRow2, 

 1683  CoinIndexedVector *spareColumn1, 

 1684  CoinIndexedVector *spareColumn2) 

[225]  1685  { 

[2385]  1686  int iSection, j; 

 1687  int number = 0; 

 1688  int *index; 

 1689  double *updateBy; 

 1690  double *reducedCost; 

 1691  // dj could be very small (or even zero  take care) 

 1692  double dj = model_>dualIn(); 

 1693  double tolerance = model_>currentDualTolerance(); 

 1694  // we can't really trust infeasibilities if there is dual error 

 1695  // this coding has to mimic coding in checkDualSolution 

 1696  double error = CoinMin(1.0e2, model_>largestDualError()); 

 1697  // allow tolerance at least slightly bigger than standard 

 1698  tolerance = tolerance + error; 

 1699  int pivotRow = model_>pivotRow(); 

 1700  double *infeas = infeasible_>denseVector(); 

 1701  //updates>scanAndPack(); 

 1702  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

[1197]  1703  

[2385]  1704  // put row of tableau in rowArray and columnArray 

 1705  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 1706  updates, spareColumn2, spareColumn1); 

[2235]  1707  #ifdef CLP_USER_DRIVEN 

[2385]  1708  model_>eventHandler()>eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates); 

[2235]  1709  #endif 

[2385]  1710  // normal 

 1711  for (iSection = 0; iSection < 2; iSection++) { 

[517]  1712  

[2385]  1713  reducedCost = model_>djRegion(iSection); 

 1714  int addSequence; 

 1715  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 1716  double slack_multiplier; 

[1732]  1717  #endif 

[1197]  1718  

[2385]  1719  if (!iSection) { 

 1720  number = updates>getNumElements(); 

 1721  index = updates>getIndices(); 

 1722  updateBy = updates>denseVector(); 

 1723  addSequence = model_>numberColumns(); 

 1724  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 1725  slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER; 

[1732]  1726  #endif 

[2385]  1727  } else { 

 1728  number = spareColumn1>getNumElements(); 

 1729  index = spareColumn1>getIndices(); 

 1730  updateBy = spareColumn1>denseVector(); 

 1731  addSequence = 0; 

 1732  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

 1733  slack_multiplier = 1.0; 

[1732]  1734  #endif 

[2385]  1735  } 

[1525]  1736  

[2385]  1737  for (j = 0; j < number; j++) { 

 1738  int iSequence = index[j]; 

 1739  double value = reducedCost[iSequence]; 

 1740  value = updateBy[j]; 

 1741  updateBy[j] = 0.0; 

 1742  reducedCost[iSequence] = value; 

 1743  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  1744  

[2385]  1745  switch (status) { 

[1525]  1746  

[2385]  1747  case ClpSimplex::basic: 

 1748  infeasible_>zero(iSequence + addSequence); 

 1749  case ClpSimplex::isFixed: 

 1750  break; 

 1751  case ClpSimplex::isFree: 

 1752  case ClpSimplex::superBasic: 

 1753  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 1754  // we are going to bias towards free (but only if reasonable) 

 1755  value *= FREE_BIAS; 

 1756  // store square in list 

 1757  if (infeas[iSequence + addSequence]) 

 1758  infeas[iSequence + addSequence] = value * value; // already there 

 1759  else 

 1760  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 1761  } else { 

 1762  infeasible_>zero(iSequence + addSequence); 

 1763  } 

 1764  break; 

 1765  case ClpSimplex::atUpperBound: 

 1766  iSequence += addSequence; 

 1767  if (value > tolerance) { 

[1732]  1768  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1769  value *= value * slack_multiplier; 

[1732]  1770  #else 

[2385]  1771  value *= value; 

[1732]  1772  #endif 

[2385]  1773  // store square in list 

 1774  if (infeas[iSequence]) 

 1775  infeas[iSequence] = value; // already there 

 1776  else 

 1777  infeasible_>quickAdd(iSequence, value); 

 1778  } else { 

 1779  infeasible_>zero(iSequence); 

 1780  } 

 1781  break; 

 1782  case ClpSimplex::atLowerBound: 

 1783  iSequence += addSequence; 

 1784  if (value < tolerance) { 

[1732]  1785  #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 

[2385]  1786  value *= value * slack_multiplier; 

[1732]  1787  #else 

[2385]  1788  value *= value; 

[1732]  1789  #endif 

[2385]  1790  // store square in list 

 1791  if (infeas[iSequence]) 

 1792  infeas[iSequence] = value; // already there 

 1793  else 

 1794  infeasible_>quickAdd(iSequence, value); 

 1795  } else { 

 1796  infeasible_>zero(iSequence); 

 1797  } 

 1798  } 

 1799  } 

 1800  } 

 1801  // we can zero out as will have to get pivot row 

 1802  // ***** move 

 1803  updates>setNumElements(0); 

 1804  spareColumn1>setNumElements(0); 

 1805  if (pivotRow >= 0) { 

 1806  // make sure infeasibility on incoming is 0.0 

 1807  int sequenceIn = model_>sequenceIn(); 

 1808  infeasible_>zero(sequenceIn); 

 1809  } 

 1810  // for weights update we use pivotSequence 

 1811  pivotRow = pivotSequence_; 

 1812  // unset in case sub flip 

 1813  pivotSequence_ = 1; 

 1814  if (pivotRow >= 0) { 

 1815  // make sure infeasibility on incoming is 0.0 

 1816  const int *pivotVariable = model_>pivotVariable(); 

 1817  int sequenceIn = pivotVariable[pivotRow]; 

 1818  assert(sequenceIn == model_>sequenceIn()); 

 1819  infeasible_>zero(sequenceIn); 

 1820  // and we can see if reference 

 1821  double referenceIn; 

 1822  if (mode_ != 1) { 

 1823  if (reference(sequenceIn)) 

 1824  referenceIn = 1.0; 

 1825  else 

 1826  referenceIn = 0.0; 

 1827  } else { 

 1828  referenceIn = 1.0; 

 1829  } 

 1830  // save outgoing weight round update 

 1831  double outgoingWeight = 0.0; 

 1832  int sequenceOut = model_>sequenceOut(); 

 1833  if (sequenceOut >= 0) 

 1834  outgoingWeight = weights_[sequenceOut]; 

 1835  // update weights 

 1836  updates>setNumElements(0); 

 1837  spareColumn1>setNumElements(0); 

 1838  // might as well set dj to 1 

 1839  dj = 1.0; 

 1840  updates>createPacked(1, &pivotRow, &dj); 

 1841  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 1842  bool needSubset = (mode_ < 4  numberSwitched_ > 1  mode_ >= 10); 

 1843  

 1844  double *weight; 

 1845  double *other = alternateWeights_>denseVector(); 

 1846  int numberColumns = model_>numberColumns(); 

 1847  // rows 

 1848  number = updates>getNumElements(); 

 1849  index = updates>getIndices(); 

 1850  updateBy = updates>denseVector(); 

 1851  weight = weights_ + numberColumns; 

 1852  if (needSubset) { 

 1853  #if ALT_UPDATE_WEIGHTS != 2 

 1854  model_>factorization()>updateColumnTranspose(spareRow2, 

 1855  alternateWeights_); 

 1856  #elif ALT_UPDATE_WEIGHTS == 1 

 1857  if (altVector[1]) { 

 1858  int numberRows = model_>numberRows(); 

 1859  double *work1 = altVector[1]>denseVector(); 

 1860  double *worka = alternateWeights_>denseVector(); 

 1861  int iRow = 1; 

 1862  double diff = 1.0e8; 

 1863  for (int i = 0; i < numberRows; i++) { 

 1864  double dd = CoinMax(fabs(work1[i]), fabs(worka[i])); 

 1865  double d = fabs(work1[i]  worka[i]); 

 1866  if (dd > 1.0e6 && d > diff * dd) { 

 1867  diff = d / dd; 

 1868  iRow = i; 

[1525]  1869  } 

[2385]  1870  } 

 1871  if (iRow >= 0) 

 1872  printf("largest2 difference of %g (%g,%g) on row %d\n", 

 1873  diff, work1[iRow], worka[iRow], iRow); 

 1874  } 

[2235]  1875  #endif 

[2385]  1876  // do alternateWeights_ here so can scale 

 1877  for (j = 0; j < number; j++) { 

 1878  int iSequence = index[j]; 

 1879  assert(iSequence >= 0 && iSequence < model_>numberRows()); 

 1880  double thisWeight = weight[iSequence]; 

 1881  // row has 1 

 1882  double pivot = updateBy[j]; 

 1883  double modification = other[iSequence]; 

 1884  double pivotSquared = pivot * pivot; 

[1525]  1885  

[2385]  1886  thisWeight += pivotSquared * devex_ + pivot * modification; 

 1887  if (thisWeight < TRY_NORM) { 

 1888  if (mode_ == 1) { 

 1889  // steepest 

 1890  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

[1525]  1891  } else { 

[2385]  1892  // exact 

 1893  thisWeight = referenceIn * pivotSquared; 

 1894  if (reference(iSequence + numberColumns)) 

 1895  thisWeight += 1.0; 

 1896  thisWeight = CoinMax(thisWeight, TRY_NORM); 

[1525]  1897  } 

[2385]  1898  } 

 1899  weight[iSequence] = thisWeight; 

 1900  } 

 1901  transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2, 0.0); 

 1902  } else { 

 1903  // put row of tableau in rowArray and columnArray 

 1904  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 1905  updates, spareColumn2, spareColumn1); 

 1906  } 

[1525]  1907  

[2385]  1908  if (needSubset) { 

 1909  CoinZeroN(updateBy, number); 

 1910  } else if (mode_ == 4) { 

 1911  // Devex 

 1912  assert(devex_ > 0.0); 

 1913  for (j = 0; j < number; j++) { 

 1914  int iSequence = index[j]; 

 1915  double thisWeight = weight[iSequence]; 

 1916  // row has 1 

 1917  double pivot = updateBy[j]; 

 1918  updateBy[j] = 0.0; 

 1919  double value = pivot * pivot * devex_; 

 1920  if (reference(iSequence + numberColumns)) 

 1921  value += 1.0; 

 1922  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 1923  } 

 1924  } 

[1525]  1925  

[2385]  1926  // columns 

 1927  weight = weights_; 

[1525]  1928  

[2385]  1929  number = spareColumn1>getNumElements(); 

 1930  index = spareColumn1>getIndices(); 

 1931  updateBy = spareColumn1>denseVector(); 

 1932  if (needSubset) { 

 1933  // Exact  already done 

 1934  } else if (mode_ == 4) { 

 1935  // Devex 

 1936  for (j = 0; j < number; j++) { 

 1937  int iSequence = index[j]; 

 1938  double thisWeight = weight[iSequence]; 

 1939  // row has 1 

 1940  double pivot = updateBy[j]; 

 1941  updateBy[j] = 0.0; 

 1942  double value = pivot * pivot * devex_; 

 1943  if (reference(iSequence)) 

 1944  value += 1.0; 

 1945  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 1946  } 

 1947  } 

 1948  // restore outgoing weight 

 1949  if (sequenceOut >= 0) 

 1950  weights_[sequenceOut] = outgoingWeight; 

 1951  alternateWeights_>clear(); 

 1952  spareColumn2>setNumElements(0); 

 1953  //#define SOME_DEBUG_1 

[225]  1954  #ifdef SOME_DEBUG_1 

[2385]  1955  // check for accuracy 

 1956  int iCheck = 892; 

 1957  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 1958  int numberRows = model_>numberRows(); 

 1959  //int numberColumns = model_>numberColumns(); 

 1960  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 1961  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 1962  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 1963  } 

[225]  1964  #endif 

[2385]  1965  } 

 1966  updates>setNumElements(0); 

 1967  spareColumn1>setNumElements(0); 

[225]  1968  } 

[517]  1969  // Updates two arrays for steepest 

[2385]  1970  int ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector *pi1, CoinIndexedVector *dj1, 

 1971  const CoinIndexedVector *pi2, CoinIndexedVector *dj2, 

 1972  CoinIndexedVector *spare, 

 1973  double scaleFactor) 

[517]  1974  { 

[2385]  1975  // see if reference 

 1976  int sequenceIn = model_>sequenceIn(); 

 1977  double referenceIn; 

 1978  if (mode_ != 1) { 

 1979  if (reference(sequenceIn)) 

 1980  referenceIn = 1.0; 

 1981  else 

 1982  referenceIn = 0.0; 

 1983  } else { 

 1984  referenceIn = 1.0; 

 1985  } 

 1986  int returnCode = 0; 

 1987  if (model_>clpMatrix()>canCombine(model_, pi1)) { 

 1988  double *infeas = scaleFactor ? infeasible_>denseVector() : NULL; 

 1989  // put row of tableau in rowArray and columnArray 

 1990  returnCode = model_>clpMatrix()>transposeTimes2(model_, pi1, 

 1991  dj1, pi2, spare, 

 1992  infeas, 

 1993  model_>djRegion(1), 

 1994  referenceIn, devex_, 

 1995  reference_, 

 1996  weights_, scaleFactor); 

 1997  if (model_>spareIntArray_[3] > 2) 

 1998  returnCode = 2; 

 1999  } else { 

 2000  // put row of tableau in rowArray and columnArray 

 2001  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 2002  pi1, dj2, dj1); 

 2003  // get subset which have nonzero tableau elements 

 2004  model_>clpMatrix()>subsetTransposeTimes(model_, pi2, dj1, dj2); 

 2005  bool killDjs = (scaleFactor == 0.0); 

 2006  if (!scaleFactor) 

 2007  scaleFactor = 1.0; 

 2008  // columns 

 2009  double *weight = weights_; 

[1525]  2010  

[2385]  2011  int number = dj1>getNumElements(); 

 2012  const int *index = dj1>getIndices(); 

 2013  double *updateBy = dj1>denseVector(); 

 2014  double *updateBy2 = dj2>denseVector(); 

[1525]  2015  

[2385]  2016  for (int j = 0; j < number; j++) { 

 2017  double thisWeight; 

 2018  double pivot; 

 2019  double pivotSquared; 

 2020  int iSequence = index[j]; 

 2021  double value2 = updateBy[j]; 

 2022  if (killDjs) 

 2023  updateBy[j] = 0.0; 

 2024  double modification = updateBy2[j]; 

 2025  updateBy2[j] = 0.0; 

 2026  ClpSimplex::Status status = model_>getStatus(iSequence); 

[1525]  2027  

[2385]  2028  if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) { 

 2029  thisWeight = weight[iSequence]; 

 2030  pivot = value2 * scaleFactor; 

 2031  pivotSquared = pivot * pivot; 

[1525]  2032  

[2385]  2033  thisWeight += pivotSquared * devex_ + pivot * modification; 

 2034  if (thisWeight < TRY_NORM) { 

 2035  if (referenceIn < 0.0) { 

 2036  // steepest 

 2037  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 2038  } else { 

 2039  // exact 

 2040  thisWeight = referenceIn * pivotSquared; 

 2041  if (reference(iSequence)) 

 2042  thisWeight += 1.0; 

 2043  thisWeight = CoinMax(thisWeight, TRY_NORM); 

[517]  2044  } 

[2385]  2045  } 

 2046  weight[iSequence] = thisWeight; 

 2047  } 

 2048  } 

 2049  } 

 2050  dj2>setNumElements(0); 

 2051  return returnCode; 

[517]  2052  } 

[225]  2053  // Update weights for Devex 

[2385]  2054  void ClpPrimalColumnSteepest::justDevex(CoinIndexedVector *updates, 

 2055  CoinIndexedVector *spareRow2, 

 2056  CoinIndexedVector *spareColumn1, 

 2057  CoinIndexedVector *spareColumn2) 

[225]  2058  { 

[2385]  2059  int j; 

 2060  int number = 0; 

 2061  int *index; 

 2062  double *updateBy; 

 2063  // dj could be very small (or even zero  take care) 

 2064  double dj = model_>dualIn(); 

 2065  double tolerance = model_>currentDualTolerance(); 

 2066  // we can't really trust infeasibilities if there is dual error 

 2067  // this coding has to mimic coding in checkDualSolution 

 2068  double error = CoinMin(1.0e2, model_>largestDualError()); 

 2069  // allow tolerance at least slightly bigger than standard 

 2070  tolerance = tolerance + error; 

 2071  int pivotRow = model_>pivotRow(); 

 2072  // for weights update we use pivotSequence 

 2073  pivotRow = pivotSequence_; 

 2074  assert(pivotRow >= 0); 

 2075  // make sure infeasibility on incoming is 0.0 

 2076  const int *pivotVariable = model_>pivotVariable(); 

 2077  int sequenceIn = pivotVariable[pivotRow]; 

 2078  infeasible_>zero(sequenceIn); 

 2079  // and we can see if reference 

 2080  //double referenceIn = 0.0; 

 2081  //if (mode_ != 1 && reference(sequenceIn)) 

 2082  // referenceIn = 1.0; 

 2083  // save outgoing weight round update 

 2084  double outgoingWeight = 0.0; 

 2085  int sequenceOut = model_>sequenceOut(); 

 2086  if (sequenceOut >= 0) 

 2087  outgoingWeight = weights_[sequenceOut]; 

 2088  assert(!updates>getNumElements()); 

 2089  assert(!spareColumn1>getNumElements()); 

 2090  // unset in case sub flip 

 2091  pivotSequence_ = 1; 

 2092  // might as well set dj to 1 

 2093  dj = 1.0; 

 2094  updates>createPacked(1, &pivotRow, &dj); 

 2095  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 2096  // put row of tableau in rowArray and columnArray 

 2097  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 2098  updates, spareColumn2, spareColumn1); 

 2099  double *weight; 

 2100  int numberColumns = model_>numberColumns(); 

 2101  // rows 

 2102  number = updates>getNumElements(); 

 2103  index = updates>getIndices(); 

 2104  updateBy = updates>denseVector(); 

 2105  weight = weights_ + numberColumns; 

[1525]  2106  

[2385]  2107  // Devex 

 2108  assert(devex_ > 0.0); 

 2109  for (j = 0; j < number; j++) { 

 2110  int iSequence = index[j]; 

 2111  double thisWeight = weight[iSequence]; 

 2112  // row has 1 

 2113  double pivot = updateBy[j]; 

 2114  updateBy[j] = 0.0; 

 2115  double value = pivot * pivot * devex_; 

 2116  if (reference(iSequence + numberColumns)) 

 2117  value += 1.0; 

 2118  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 2119  } 

[1525]  2120  

[2385]  2121  // columns 

 2122  weight = weights_; 

[1525]  2123  

[2385]  2124  number = spareColumn1>getNumElements(); 

 2125  index = spareColumn1>getIndices(); 

 2126  updateBy = spareColumn1>denseVector(); 

 2127  // Devex 

 2128  for (j = 0; j < number; j++) { 

 2129  int iSequence = index[j]; 

 2130  double thisWeight = weight[iSequence]; 

 2131  // row has 1 

 2132  double pivot = updateBy[j]; 

 2133  updateBy[j] = 0.0; 

 2134  double value = pivot * pivot * devex_; 

 2135  if (reference(iSequence)) 

 2136  value += 1.0; 

 2137  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 2138  } 

 2139  // restore outgoing weight 

 2140  if (sequenceOut >= 0) 

 2141  weights_[sequenceOut] = outgoingWeight; 

 2142  spareColumn2>setNumElements(0); 

 2143  //#define SOME_DEBUG_1 

[225]  2144  #ifdef SOME_DEBUG_1 

[2385]  2145  // check for accuracy 

 2146  int iCheck = 892; 

 2147  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 2148  int numberRows = model_>numberRows(); 

 2149  //int numberColumns = model_>numberColumns(); 

 2150  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 2151  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 2152  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 2153  } 

[225]  2154  #endif 

[2385]  2155  updates>setNumElements(0); 

 2156  spareColumn1>setNumElements(0); 

[225]  2157  } 

 2158  // Update weights for Steepest 

[2385]  2159  void ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector *updates, 

 2160  CoinIndexedVector *spareRow2, 

 2161  CoinIndexedVector *spareColumn1, 

 2162  CoinIndexedVector *spareColumn2) 

[225]  2163  { 

[2385]  2164  int j; 

 2165  int number = 0; 

 2166  int *index; 

 2167  double *updateBy; 

 2168  // dj could be very small (or even zero  take care) 

 2169  double dj = model_>dualIn(); 

 2170  double tolerance = model_>currentDualTolerance(); 

 2171  // we can't really trust infeasibilities if there is dual error 

 2172  // this coding has to mimic coding in checkDualSolution 

 2173  double error = CoinMin(1.0e2, model_>largestDualError()); 

 2174  // allow tolerance at least slightly bigger than standard 

 2175  tolerance = tolerance + error; 

 2176  int pivotRow = model_>pivotRow(); 

 2177  // for weights update we use pivotSequence 

 2178  pivotRow = pivotSequence_; 

 2179  // unset in case sub flip 

 2180  pivotSequence_ = 1; 

 2181  assert(pivotRow >= 0); 

 2182  // make sure infeasibility on incoming is 0.0 

 2183  const int *pivotVariable = model_>pivotVariable(); 

 2184  int sequenceIn = pivotVariable[pivotRow]; 

 2185  infeasible_>zero(sequenceIn); 

 2186  // and we can see if reference 

 2187  double referenceIn = 0.0; 

 2188  if (mode_ != 1 && reference(sequenceIn)) 

 2189  referenceIn = 1.0; 

 2190  // save outgoing weight round update 

 2191  double outgoingWeight = 0.0; 

 2192  int sequenceOut = model_>sequenceOut(); 

 2193  if (sequenceOut >= 0) 

 2194  outgoingWeight = weights_[sequenceOut]; 

 2195  assert(!updates>getNumElements()); 

 2196  assert(!spareColumn1>getNumElements()); 

 2197  // update weights 

 2198  //updates>setNumElements(0); 

 2199  //spareColumn1>setNumElements(0); 

 2200  // might as well set dj to 1 

 2201  dj = 1.0; 

 2202  updates>createPacked(1, &pivotRow, &dj); 

 2203  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 2204  // put row of tableau in rowArray and columnArray 

 2205  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 2206  updates, spareColumn2, spareColumn1); 

 2207  double *weight; 

 2208  double *other = alternateWeights_>denseVector(); 

 2209  int numberColumns = model_>numberColumns(); 

 2210  // rows 

 2211  number = updates>getNumElements(); 

 2212  index = updates>getIndices(); 

 2213  updateBy = updates>denseVector(); 

 2214  weight = weights_ + numberColumns; 

[1525]  2215  

[2385]  2216  // Exact 

 2217  // now update weight update array 

 2218  //alternateWeights_>scanAndPack(); 

 2219  #if ALT_UPDATE_WEIGHTS != 2 

 2220  model_>factorization()>updateColumnTranspose(spareRow2, 

 2221  alternateWeights_); 

 2222  #elif ALT_UPDATE_WEIGHTS == 1 

 2223  if (altVector[1]) { 

 2224  int numberRows = model_>numberRows(); 

 2225  double *work1 = altVector[1]>denseVector(); 

 2226  double *worka = alternateWeights_>denseVector(); 

 2227  int iRow = 1; 

 2228  double diff = 1.0e8; 

 2229  for (int i = 0; i < numberRows; i++) { 

 2230  double dd = CoinMax(fabs(work1[i]), fabs(worka[i])); 

 2231  double d = fabs(work1[i]  worka[i]); 

 2232  if (dd > 1.0e6 && d > diff * dd) { 

 2233  diff = d / dd; 

 2234  iRow = i; 

 2235  } 

 2236  } 

 2237  if (iRow >= 0) 

 2238  printf("largest3 difference of %g (%g,%g) on row %d\n", 

 2239  diff, work1[iRow], worka[iRow], iRow); 

 2240  } 

[2235]  2241  #endif 

[2385]  2242  // get subset which have nonzero tableau elements 

 2243  model_>clpMatrix()>subsetTransposeTimes(model_, alternateWeights_, 

 2244  spareColumn1, 

 2245  spareColumn2); 

 2246  for (j = 0; j < number; j++) { 

 2247  int iSequence = index[j]; 

 2248  double thisWeight = weight[iSequence]; 

 2249  // row has 1 

 2250  double pivot = updateBy[j]; 

 2251  updateBy[j] = 0.0; 

 2252  double modification = other[iSequence]; 

 2253  double pivotSquared = pivot * pivot; 

[1525]  2254  

[2385]  2255  thisWeight += pivotSquared * devex_ + pivot * modification; 

 2256  if (thisWeight < TRY_NORM) { 

 2257  if (mode_ == 1) { 

 2258  // steepest 

 2259  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 2260  } else { 

 2261  // exact 

 2262  thisWeight = referenceIn * pivotSquared; 

 2263  if (reference(iSequence + numberColumns)) 

 2264  thisWeight += 1.0; 

 2265  thisWeight = CoinMax(thisWeight, TRY_NORM); 

 2266  } 

 2267  } 

 2268  weight[iSequence] = thisWeight; 

 2269  } 

[1525]  2270  

[2385]  2271  // columns 

 2272  weight = weights_; 

[1525]  2273  

[2385]  2274  number = spareColumn1>getNumElements(); 

 2275  index = spareColumn1>getIndices(); 

 2276  updateBy = spareColumn1>denseVector(); 

 2277  // Exact 

 2278  double *updateBy2 = spareColumn2>denseVector(); 

 2279  for (j = 0; j < number; j++) { 

 2280  int iSequence = index[j]; 

 2281  double thisWeight = weight[iSequence]; 

 2282  double pivot = updateBy[j]; 

 2283  updateBy[j] = 0.0; 

 2284  double modification = updateBy2[j]; 

 2285  updateBy2[j] = 0.0; 

 2286  double pivotSquared = pivot * pivot; 

[1525]  2287  

[2385]  2288  thisWeight += pivotSquared * devex_ + pivot * modification; 

 2289  if (thisWeight < TRY_NORM) { 

 2290  if (mode_ == 1) { 

 2291  // steepest 

 2292  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 2293  } else { 

 2294  // exact 

 2295  thisWeight = referenceIn * pivotSquared; 

 2296  if (reference(iSequence)) 

 2297  thisWeight += 1.0; 

 2298  thisWeight = CoinMax(thisWeight, TRY_NORM); 

 2299  } 

 2300  } 

 2301  weight[iSequence] = thisWeight; 

 2302  } 

 2303  // restore outgoing weight 

 2304  if (sequenceOut >= 0) 

 2305  weights_[sequenceOut] = outgoingWeight; 

 2306  alternateWeights_>clear(); 

 2307  spareColumn2>setNumElements(0); 

 2308  //#define SOME_DEBUG_1 

[225]  2309  #ifdef SOME_DEBUG_1 

[2385]  2310  // check for accuracy 

 2311  int iCheck = 892; 

 2312  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 2313  int numberRows = model_>numberRows(); 

 2314  //int numberColumns = model_>numberColumns(); 

 2315  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 2316  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 2317  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 2318  } 

[225]  2319  #endif 

[2385]  2320  updates>setNumElements(0); 

 2321  spareColumn1>setNumElements(0); 

[225]  2322  } 

 2323  // Returns pivot column, 1 if none 

[2385]  2324  int ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector *updates, 

 2325  CoinIndexedVector *, 

 2326  CoinIndexedVector *spareRow2, 

 2327  CoinIndexedVector *spareColumn1, 

 2328  CoinIndexedVector *spareColumn2) 

[225]  2329  { 

[2385]  2330  assert(model_); 

 2331  int iSection, j; 

 2332  int number = 0; 

 2333  int *index; 

 2334  double *updateBy; 

 2335  double *reducedCost; 

 2336  // dj could be very small (or even zero  take care) 

 2337  double dj = model_>dualIn(); 

 2338  double tolerance = model_>currentDualTolerance(); 

 2339  // we can't really trust infeasibilities if there is dual error 

 2340  // this coding has to mimic coding in checkDualSolution 

 2341  double error = CoinMin(1.0e2, model_>largestDualError()); 

 2342  // allow tolerance at least slightly bigger than standard 

 2343  tolerance = tolerance + error; 

 2344  int pivotRow = model_>pivotRow(); 

 2345  int anyUpdates; 

 2346  double *infeas = infeasible_>denseVector(); 

[2]  2347  

[2385]  2348  // Local copy of mode so can decide what to do 

 2349  int switchType; 

 2350  if (mode_ == 4) 

 2351  switchType = 5  numberSwitched_; 

 2352  else if (mode_ >= 10) 

 2353  switchType = 3; 

 2354  else 

 2355  switchType = mode_; 

 2356  /* switchType  

[1525]  2357  0  all exact devex 

 2358  1  all steepest 

 2359  2  some exact devex 

 2360  3  auto some exact devex 

 2361  4  devex 

 2362  5  dantzig 

 2363  */ 

[2385]  2364  if (updates>getNumElements()) { 

 2365  // would have to have two goes for devex, three for steepest 

 2366  anyUpdates = 2; 

 2367  // add in pivot contribution 

 2368  if (pivotRow >= 0) 

 2369  updates>add(pivotRow, dj); 

 2370  } else if (pivotRow >= 0) { 

 2371  if (fabs(dj) > 1.0e15) { 

 2372  // some dj 

 2373  updates>insert(pivotRow, dj); 

 2374  if (fabs(dj) > 1.0e6) { 

 2375  // reasonable size 

 2376  anyUpdates = 1; 

 2377  } else { 

 2378  // too small 

 2379  anyUpdates = 2; 

 2380  } 

 2381  } else { 

 2382  // zero dj 

 2383  anyUpdates = 1; 

 2384  } 

 2385  } else if (pivotSequence_ >= 0) { 

 2386  // just after refactorization 

 2387  anyUpdates = 1; 

 2388  } else { 

 2389  // sub flip  nothing to do 

 2390  anyUpdates = 0; 

 2391  } 

[2]  2392  

[2385]  2393  if (anyUpdates > 0) { 

 2394  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

[225]  2395  

[2385]  2396  // put row of tableau in rowArray and columnArray 

 2397  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 2398  updates, spareColumn2, spareColumn1); 

 2399  // normal 

 2400  for (iSection = 0; iSection < 2; iSection++) { 

[2]  2401  

[2385]  2402  reducedCost = model_>djRegion(iSection); 

 2403  int addSequence; 

[1197]  2404  

[2385]  2405  if (!iSection) { 

 2406  number = updates>getNumElements(); 

 2407  index = updates>getIndices(); 

 2408  updateBy = updates>denseVector(); 

 2409  addSequence = model_>numberColumns(); 

 2410  } else { 

 2411  number = spareColumn1>getNumElements(); 

 2412  index = spareColumn1>getIndices(); 

 2413  updateBy = spareColumn1>denseVector(); 

 2414  addSequence = 0; 

 2415  } 

 2416  if (!model_>nonLinearCost()>lookBothWays()) { 

[1197]  2417  

[2385]  2418  for (j = 0; j < number; j++) { 

 2419  int iSequence = index[j]; 

 2420  double value = reducedCost[iSequence]; 

 2421  value = updateBy[iSequence]; 

 2422  reducedCost[iSequence] = value; 

 2423  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  2424  

[2385]  2425  switch (status) { 

[1525]  2426  

[2385]  2427  case ClpSimplex::basic: 

 2428  infeasible_>zero(iSequence + addSequence); 

 2429  case ClpSimplex::isFixed: 

 2430  break; 

 2431  case ClpSimplex::isFree: 

 2432  case ClpSimplex::superBasic: 

 2433  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 2434  // we are going to bias towards free (but only if reasonable) 

 2435  value *= FREE_BIAS; 

 2436  // store square in list 

 2437  if (infeas[iSequence + addSequence]) 

 2438  infeas[iSequence + addSequence] = value * value; // already there 

 2439  else 

 2440  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2441  } else { 

 2442  infeasible_>zero(iSequence + addSequence); 

 2443  } 

 2444  break; 

 2445  case ClpSimplex::atUpperBound: 

 2446  if (value > tolerance) { 

 2447  // store square in list 

 2448  if (infeas[iSequence + addSequence]) 

 2449  infeas[iSequence + addSequence] = value * value; // already there 

 2450  else 

 2451  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2452  } else { 

 2453  infeasible_>zero(iSequence + addSequence); 

 2454  } 

 2455  break; 

 2456  case ClpSimplex::atLowerBound: 

 2457  if (value < tolerance) { 

 2458  // store square in list 

 2459  if (infeas[iSequence + addSequence]) 

 2460  infeas[iSequence + addSequence] = value * value; // already there 

 2461  else 

 2462  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2463  } else { 

 2464  infeasible_>zero(iSequence + addSequence); 

 2465  } 

[1525]  2466  } 

[2385]  2467  } 

 2468  } else { 

 2469  ClpNonLinearCost *nonLinear = model_>nonLinearCost(); 

 2470  // We can go up OR down 

 2471  for (j = 0; j < number; j++) { 

 2472  int iSequence = index[j]; 

 2473  double value = reducedCost[iSequence]; 

 2474  value = updateBy[iSequence]; 

 2475  reducedCost[iSequence] = value; 

 2476  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  2477  

[2385]  2478  switch (status) { 

[1525]  2479  

 2480  case ClpSimplex::basic: 

[2385]  2481  infeasible_>zero(iSequence + addSequence); 

[1525]  2482  case ClpSimplex::isFixed: 

[2385]  2483  break; 

[1525]  2484  case ClpSimplex::isFree: 

 2485  case ClpSimplex::superBasic: 

[2385]  2486  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 2487  // we are going to bias towards free (but only if reasonable) 

 2488  value *= FREE_BIAS; 

 2489  // store square in list 

 2490  if (infeas[iSequence + addSequence]) 

 2491  infeas[iSequence + addSequence] = value * value; // already there 

 2492  else 

 2493  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2494  } else { 

 2495  infeasible_>zero(iSequence + addSequence); 

 2496  } 

 2497  break; 

[1525]  2498  case ClpSimplex::atUpperBound: 

[2385]  2499  if (value > tolerance) { 

 2500  // store square in list 

 2501  if (infeas[iSequence + addSequence]) 

 2502  infeas[iSequence + addSequence] = value * value; // already there 

 2503  else 

 2504  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2505  } else { 

 2506  // look other way  change up should be negative 

 2507  value = nonLinear>changeUpInCost(iSequence + addSequence); 

 2508  if (value > tolerance) { 

 2509  infeasible_>zero(iSequence + addSequence); 

 2510  } else { 

 2511  // store square in list 

 2512  if (infeas[iSequence + addSequence]) 

 2513  infeas[iSequence + addSequence] = value * value; // already there 

 2514  else 

 2515  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2516  } 

 2517  } 

 2518  break; 

[1525]  2519  case ClpSimplex::atLowerBound: 

[2385]  2520  if (value < tolerance) { 

 2521  // store square in list 

 2522  if (infeas[iSequence + addSequence]) 

 2523  infeas[iSequence + addSequence] = value * value; // already there 

 2524  else 

 2525  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2526  } else { 

 2527  // look other way  change down should be positive 

 2528  value = nonLinear>changeDownInCost(iSequence + addSequence); 

 2529  if (value < tolerance) { 

 2530  infeasible_>zero(iSequence + addSequence); 

 2531  } else { 

 2532  // store square in list 

 2533  if (infeas[iSequence + addSequence]) 

 2534  infeas[iSequence + addSequence] = value * value; // already there 

 2535  else 

 2536  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2537  } 

 2538  } 

[1525]  2539  } 

[2385]  2540  } 

 2541  } 

 2542  } 

 2543  if (anyUpdates == 2) { 

 2544  // we can zero out as will have to get pivot row 

 2545  updates>clear(); 

 2546  spareColumn1>clear(); 

 2547  } 

 2548  if (pivotRow >= 0) { 

 2549  // make sure infeasibility on incoming is 0.0 

 2550  int sequenceIn = model_>sequenceIn(); 

 2551  infeasible_>zero(sequenceIn); 

 2552  } 

 2553  } 

 2554  // make sure outgoing from last iteration okay 

 2555  int sequenceOut = model_>sequenceOut(); 

 2556  if (sequenceOut >= 0) { 

 2557  ClpSimplex::Status status = model_>getStatus(sequenceOut); 

 2558  double value = model_>reducedCost(sequenceOut); 

[1525]  2559  

[2385]  2560  switch (status) { 

[1525]  2561  

[2385]  2562  case ClpSimplex::basic: 

 2563  case ClpSimplex::isFixed: 

 2564  break; 

 2565  case ClpSimplex::isFree: 

 2566  case ClpSimplex::superBasic: 

 2567  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 2568  // we are going to bias towards free (but only if reasonable) 

 2569  value *= FREE_BIAS; 

 2570  // store square in list 

 2571  if (infeas[sequenceOut]) 

 2572  infeas[sequenceOut] = value * value; // already there 

 2573  else 

 2574  infeasible_>quickAdd(sequenceOut, value * value); 

 2575  } else { 

 2576  infeasible_>zero(sequenceOut); 

 2577  } 

 2578  break; 

 2579  case ClpSimplex::atUpperBound: 

 2580  if (value > tolerance) { 

 2581  // store square in list 

 2582  if (infeas[sequenceOut]) 

 2583  infeas[sequenceOut] = value * value; // already there 

 2584  else 

 2585  infeasible_>quickAdd(sequenceOut, value * value); 

 2586  } else { 

 2587  infeasible_>zero(sequenceOut); 

 2588  } 

 2589  break; 

 2590  case ClpSimplex::atLowerBound: 

 2591  if (value < tolerance) { 

 2592  // store square in list 

 2593  if (infeas[sequenceOut]) 

 2594  infeas[sequenceOut] = value * value; // already there 

 2595  else 

 2596  infeasible_>quickAdd(sequenceOut, value * value); 

 2597  } else { 

 2598  infeasible_>zero(sequenceOut); 

 2599  } 

 2600  } 

 2601  } 

[1525]  2602  

[2385]  2603  // If in quadratic recompute all 

 2604  if (model_>algorithm() == 2) { 

 2605  for (iSection = 0; iSection < 2; iSection++) { 

[1525]  2606  

[2385]  2607  reducedCost = model_>djRegion(iSection); 

 2608  int addSequence; 

 2609  int iSequence; 

[1525]  2610  

[2385]  2611  if (!iSection) { 

 2612  number = model_>numberRows(); 

 2613  addSequence = model_>numberColumns(); 

 2614  } else { 

 2615  number = model_>numberColumns(); 

 2616  addSequence = 0; 

 2617  } 

[1525]  2618  

[2385]  2619  if (!model_>nonLinearCost()>lookBothWays()) { 

 2620  for (iSequence = 0; iSequence < number; iSequence++) { 

 2621  double value = reducedCost[iSequence]; 

 2622  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

[1525]  2623  

[2385]  2624  switch (status) { 

[1525]  2625  

[2385]  2626  case ClpSimplex::basic: 

 2627  infeasible_>zero(iSequence + addSequence); 

 2628  case ClpSimplex::isFixed: 

 2629  break; 

 2630  case ClpSimplex::isFree: 

 2631  case ClpSimplex::superBasic: 

 2632  if (fabs(value) > tolerance) { 

 2633  // we are going to bias towards free (but only if reasonable) 

 2634  value *= FREE_BIAS; 

 2635  // store square in list 

 2636  if (infeas[iSequence + addSequence]) 

 2637  infeas[iSequence + addSequence] = value * value; // already there 

 2638  else 

 2639  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2640  } else { 

 2641  infeasible_>zero(iSequence + addSequence); 

 2642  } 

 2643  break; 

 2644  case ClpSimplex::atUpperBound: 

 2645  if (value > tolerance) { 

 2646  // store square in list 

 2647  if (infeas[iSequence + addSequence]) 

 2648  infeas[iSequence + addSequence] = value * value; // already there 

 2649  else 

 2650  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2651  } else { 

 2652  infeasible_>zero(iSequence + addSequence); 

 2653  } 

 2654  break; 

 2655  case ClpSimplex::atLowerBound: 

 2656  if (value < tolerance) { 

 2657  // store square in list 

 2658  if (infeas[iSequence + addSequence]) 

 2659  infeas[iSequence + addSequence] = value * value; // already there 

 2660  else 

 2661  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2662  } else { 

 2663  infeasible_>zero(iSequence + addSequence); 

 2664  } 

[1525]  2665  } 

[2385]  2666  } 

 2667  } else { 

 2668  // we can go both ways 

 2669  ClpNonLinearCost *nonLinear = model_>nonLinearCost(); 

 2670  for (iSequence = 0; iSequence < number; iSequence++) { 

 2671  double value = reducedCost[iSequence]; 

 2672  ClpSimplex::Status status = model_>getStatus(iSequence + addSequence); 

 2673  

 2674  switch (status) { 

 2675  

 2676  case ClpSimplex::basic: 

 2677  infeasible_>zero(iSequence + addSequence); 

 2678  case ClpSimplex::isFixed: 

 2679  break; 

 2680  case ClpSimplex::isFree: 

 2681  case ClpSimplex::superBasic: 

 2682  if (fabs(value) > tolerance) { 

 2683  // we are going to bias towards free (but only if reasonable) 

 2684  value *= FREE_BIAS; 

 2685  // store square in list 

 2686  if (infeas[iSequence + addSequence]) 

 2687  infeas[iSequence + addSequence] = value * value; // already there 

 2688  else 

 2689  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2690  } else { 

 2691  infeasible_>zero(iSequence + addSequence); 

 2692  } 

 2693  break; 

 2694  case ClpSimplex::atUpperBound: 

 2695  if (value > tolerance) { 

 2696  // store square in list 

 2697  if (infeas[iSequence + addSequence]) 

 2698  infeas[iSequence + addSequence] = value * value; // already there 

 2699  else 

 2700  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2701  } else { 

 2702  // look other way  change up should be negative 

 2703  value = nonLinear>changeUpInCost(iSequence + addSequence); 

 2704  if (value > tolerance) { 

 2705  infeasible_>zero(iSequence + addSequence); 

 2706  } else { 

 2707  // store square in list 

 2708  if (infeas[iSequence + addSequence]) 

 2709  infeas[iSequence + addSequence] = value * value; // already there 

 2710  else 

 2711  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2712  } 

 2713  } 

 2714  break; 

 2715  case ClpSimplex::atLowerBound: 

 2716  if (value < tolerance) { 

 2717  // store square in list 

 2718  if (infeas[iSequence + addSequence]) 

 2719  infeas[iSequence + addSequence] = value * value; // already there 

 2720  else 

 2721  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2722  } else { 

 2723  // look other way  change down should be positive 

 2724  value = nonLinear>changeDownInCost(iSequence + addSequence); 

 2725  if (value < tolerance) { 

 2726  infeasible_>zero(iSequence + addSequence); 

 2727  } else { 

 2728  // store square in list 

 2729  if (infeas[iSequence + addSequence]) 

 2730  infeas[iSequence + addSequence] = value * value; // already there 

 2731  else 

 2732  infeasible_>quickAdd(iSequence + addSequence, value * value); 

 2733  } 

 2734  } 

[1525]  2735  } 

[2385]  2736  } 

 2737  } 

 2738  } 

 2739  } 

 2740  // See what sort of pricing 

 2741  int numberWanted = 10; 

 2742  number = infeasible_>getNumElements(); 

 2743  int numberColumns = model_>numberColumns(); 

 2744  if (switchType == 5) { 

 2745  // we can zero out 

 2746  updates>clear(); 

 2747  spareColumn1>clear(); 

 2748  pivotSequence_ = 1; 

 2749  pivotRow = 1; 

 2750  // See if to switch 

 2751  int numberRows = model_>numberRows(); 

 2752  // ratio is done on number of columns here 

 2753  //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns; 

 2754  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows); 

 2755  //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_>clpMatrix()>getNumElements(); 

 2756  if (ratio < 0.1) { 

 2757  numberWanted = CoinMax(100, number / 200); 

 2758  } else if (ratio < 0.3) { 

 2759  numberWanted = CoinMax(500, number / 40); 

 2760  } else if (ratio < 0.5  mode_ == 5) { 

 2761  numberWanted = CoinMax(2000, number / 10); 

 2762  numberWanted = CoinMax(numberWanted, numberColumns / 30); 

 2763  } else if (mode_ != 5) { 

 2764  switchType = 4; 

 2765  // initialize 

 2766  numberSwitched_++; 

 2767  // Make sure will redo 

 2768  delete[] weights_; 

 2769  weights_ = NULL; 

 2770  saveWeights(model_, 4); 

 2771  COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio)); 

 2772  updates>clear(); 

 2773  } 

 2774  if (model_>numberIterations() % 1000 == 0) 

 2775  COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted)); 

 2776  } 

 2777  if (switchType == 4) { 

 2778  // Still in devex mode 

 2779  int numberRows = model_>numberRows(); 

 2780  // ratio is done on number of rows here 

 2781  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows); 

 2782  // Go to steepest if lot of iterations? 

 2783  if (ratio < 1.0) { 

 2784  numberWanted = CoinMax(2000, number / 20); 

 2785  } else if (ratio < 5.0) { 

 2786  numberWanted = CoinMax(2000, number / 10); 

 2787  numberWanted = CoinMax(numberWanted, numberColumns / 20); 

 2788  } else { 

 2789  // we can zero out 

 2790  updates>clear(); 

 2791  spareColumn1>clear(); 

 2792  switchType = 3; 

 2793  // initialize 

 2794  pivotSequence_ = 1; 

 2795  pivotRow = 1; 

 2796  numberSwitched_++; 

 2797  // Make sure will redo 

 2798  delete[] weights_; 

 2799  weights_ = NULL; 

 2800  saveWeights(model_, 4); 

 2801  COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio)); 

 2802  updates>clear(); 

 2803  } 

 2804  if (model_>numberIterations() % 1000 == 0) 

 2805  COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted)); 

 2806  } 

 2807  if (switchType < 4) { 

 2808  if (switchType < 2) { 

 2809  numberWanted = number + 1; 

 2810  } else if (switchType == 2) { 

 2811  numberWanted = CoinMax(2000, number / 8); 

 2812  } else { 

 2813  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(model_>numberRows()); 

 2814  if (ratio < 1.0) { 

 2815  numberWanted = CoinMax(2000, number / 20); 

 2816  } else if (ratio < 5.0) { 

 2817  numberWanted = CoinMax(2000, number / 10); 

 2818  numberWanted = CoinMax(numberWanted, numberColumns / 20); 

 2819  } else if (ratio < 10.0) { 

 2820  numberWanted = CoinMax(2000, number / 8); 

 2821  numberWanted = CoinMax(numberWanted, numberColumns / 20); 

 2822  } else { 

 2823  ratio = number * (ratio / 80.0); 

 2824  if (ratio > number) { 

 2825  numberWanted = number + 1; 

 2826  } else { 

 2827  numberWanted = CoinMax(2000, static_cast< int >(ratio)); 

 2828  numberWanted = CoinMax(numberWanted, numberColumns / 10); 

 2829  } 

 2830  } 

 2831  } 

 2832  } 

 2833  // for weights update we use pivotSequence 

 2834  pivotRow = pivotSequence_; 

 2835  // unset in case sub flip 

 2836  pivotSequence_ = 1; 

 2837  if (pivotRow >= 0) { 

 2838  // make sure infeasibility on incoming is 0.0 

 2839  const int *pivotVariable = model_>pivotVariable(); 

 2840  int sequenceIn = pivotVariable[pivotRow]; 

 2841  infeasible_>zero(sequenceIn); 

 2842  // and we can see if reference 

 2843  double referenceIn = 0.0; 

 2844  if (switchType != 1 && reference(sequenceIn)) 

 2845  referenceIn = 1.0; 

 2846  // save outgoing weight round update 

 2847  double outgoingWeight = 0.0; 

 2848  if (sequenceOut >= 0) 

 2849  outgoingWeight = weights_[sequenceOut]; 

 2850  // update weights 

 2851  if (anyUpdates != 1) { 

 2852  updates>setNumElements(0); 

 2853  spareColumn1>setNumElements(0); 

 2854  // might as well set dj to 1 

 2855  dj = 1.0; 

 2856  updates>insert(pivotRow, dj); 

 2857  model_>factorization()>updateColumnTranspose(spareRow2, updates); 

 2858  // put row of tableau in rowArray and columnArray 

 2859  model_>clpMatrix()>transposeTimes(model_, 1.0, 

 2860  updates, spareColumn2, spareColumn1); 

 2861  } 

 2862  double *weight; 

 2863  double *other = alternateWeights_>denseVector(); 

 2864  int numberColumns = model_>numberColumns(); 

 2865  double scaleFactor = 1.0 / dj; // as formula is with 1.0 

 2866  // rows 

 2867  number = updates>getNumElements(); 

 2868  index = updates>getIndices(); 

 2869  updateBy = updates>denseVector(); 

 2870  weight = weights_ + numberColumns; 

[1525]  2871  

[2385]  2872  if (switchType < 4) { 

 2873  // Exact 

 2874  // now update weight update array 

 2875  model_>factorization()>updateColumnTranspose(spareRow2, 

 2876  alternateWeights_); 

[2235]  2877  #ifdef ALT_UPDATE_WEIGHTS 

[2385]  2878  abort(); 

[2235]  2879  #endif 

[2385]  2880  for (j = 0; j < number; j++) { 

 2881  int iSequence = index[j]; 

 2882  double thisWeight = weight[iSequence]; 

 2883  // row has 1 

 2884  double pivot = updateBy[iSequence] * scaleFactor; 

 2885  updateBy[iSequence] = 0.0; 

 2886  double modification = other[iSequence]; 

 2887  double pivotSquared = pivot * pivot; 

[1525]  2888  

[2385]  2889  thisWeight += pivotSquared * devex_ + pivot * modification; 

 2890  if (thisWeight < TRY_NORM) { 

 2891  if (switchType == 1) { 

 2892  // steepest 

 2893  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 2894  } else { 

 2895  // exact 

 2896  thisWeight = referenceIn * pivotSquared; 

 2897  if (reference(iSequence + numberColumns)) 

 2898  thisWeight += 1.0; 

 2899  thisWeight = CoinMax(thisWeight, TRY_NORM); 

[1525]  2900  } 

[2385]  2901  } 

 2902  weight[iSequence] = thisWeight; 

 2903  } 

 2904  } else if (switchType == 4) { 

 2905  // Devex 

 2906  assert(devex_ > 0.0); 

 2907  for (j = 0; j < number; j++) { 

 2908  int iSequence = index[j]; 

 2909  double thisWeight = weight[iSequence]; 

 2910  // row has 1 

 2911  double pivot = updateBy[iSequence] * scaleFactor; 

 2912  updateBy[iSequence] = 0.0; 

 2913  double value = pivot * pivot * devex_; 

 2914  if (reference(iSequence + numberColumns)) 

 2915  value += 1.0; 

 2916  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 2917  } 

 2918  } 

[1525]  2919  

[2385]  2920  // columns 

 2921  weight = weights_; 

[1525]  2922  

[2385]  2923  scaleFactor = scaleFactor; 

[1525]  2924  

[2385]  2925  number = spareColumn1>getNumElements(); 

 2926  index = spareColumn1>getIndices(); 

 2927  updateBy = spareColumn1>denseVector(); 

 2928  if (switchType < 4) { 

 2929  // Exact 

 2930  // get subset which have nonzero tableau elements 

 2931  model_>clpMatrix()>subsetTransposeTimes(model_, alternateWeights_, 

 2932  spareColumn1, 

 2933  spareColumn2); 

 2934  double *updateBy2 = spareColumn2>denseVector(); 

 2935  for (j = 0; j < number; j++) { 

 2936  int iSequence = index[j]; 

 2937  double thisWeight = weight[iSequence]; 

 2938  double pivot = updateBy[iSequence] * scaleFactor; 

 2939  updateBy[iSequence] = 0.0; 

 2940  double modification = updateBy2[j]; 

 2941  updateBy2[j] = 0.0; 

 2942  double pivotSquared = pivot * pivot; 

[1525]  2943  

[2385]  2944  thisWeight += pivotSquared * devex_ + pivot * modification; 

 2945  if (thisWeight < TRY_NORM) { 

 2946  if (switchType == 1) { 

 2947  // steepest 

 2948  thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared); 

 2949  } else { 

 2950  // exact 

 2951  thisWeight = referenceIn * pivotSquared; 

 2952  if (reference(iSequence)) 

 2953  thisWeight += 1.0; 

 2954  thisWeight = CoinMax(thisWeight, TRY_NORM); 

[1525]  2955  } 

[2385]  2956  } 

 2957  weight[iSequence] = thisWeight; 

 2958  } 

 2959  } else if (switchType == 4) { 

 2960  // Devex 

 2961  for (j = 0; j < number; j++) { 

 2962  int iSequence = index[j]; 

 2963  double thisWeight = weight[iSequence]; 

 2964  // row has 1 

 2965  double pivot = updateBy[iSequence] * scaleFactor; 

 2966  updateBy[iSequence] = 0.0; 

 2967  double value = pivot * pivot * devex_; 

 2968  if (reference(iSequence)) 

 2969  value += 1.0; 

 2970  weight[iSequence] = CoinMax(0.99 * thisWeight, value); 

 2971  } 

 2972  } 

 2973  // restore outgoing weight 

 2974  if (sequenceOut >= 0) 

 2975  weights_[sequenceOut] = outgoingWeight; 

 2976  alternateWeights_>clear(); 

 2977  spareColumn2>setNumElements(0); 

 2978  //#define SOME_DEBUG_1 

[2]  2979  #ifdef SOME_DEBUG_1 

[2385]  2980  // check for accuracy 

 2981  int iCheck = 892; 

 2982  //printf("weight for iCheck is %g\n",weights_[iCheck]); 

 2983  int numberRows = model_>numberRows(); 

 2984  //int numberColumns = model_>numberColumns(); 

 2985  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) { 

 2986  if (model_>getStatus(iCheck) != ClpSimplex::basic && model_>getStatus(iCheck) != ClpSimplex::isFixed) 

 2987  checkAccuracy(iCheck, 1.0e1, updates, spareRow2); 

 2988  } 

[2]  2989  #endif 

[2385]  2990  updates>setNumElements(0); 

 2991  spareColumn1>setNumElements(0); 

 2992  } 

[2]  2993  

[2385]  2994  // update of duals finished  now do pricing 

[2]  2995  

[2385]  2996  double bestDj = 1.0e30; 

 2997  int bestSequence = 1; 

[2]  2998  

[2385]  2999  int i, iSequence; 

 3000  index = infeasible_>getIndices(); 

 3001  number = infeasible_>getNumElements(); 

 3002  if (model_>numberIterations() < model_>lastBadIteration() + 200) { 

 3003  // we can't really trust infeasibilities if there is dual error 

 3004  double checkTolerance = 1.0e8; 

 3005  if (!model_>factorization()>pivots()) 

 3006  checkTolerance = 1.0e6; 

 3007  if (model_>largestDualError() > checkTolerance) 

 3008  tolerance *= model_>largestDualError() / checkTolerance; 

 3009  // But cap 

 3010  tolerance = CoinMin(1000.0, tolerance); 

 3011  } 

[225]  3012  #ifdef CLP_DEBUG 

[2385]  3013  if (model_>numberDualInfeasibilities() == 1) 

 3014  printf("** %g %g %g %x %x %d\n", tolerance, model_>dualTolerance(), 

 3015  model_>largestDualError(), model_, model_>messageHandler(), 

 3016  number); 

[225]  3017  #endif 

[2385]  3018  // stop last one coming immediately 

 3019  double saveOutInfeasibility = 0.0; 

 3020  if (sequenceOut >= 0) { 

 3021  saveOutInfeasibility = infeas[sequenceOut]; 

 3022  infeas[sequenceOut] = 0.0; 

 3023  } 

 3024  tolerance *= tolerance; // as we are using squares 

[225]  3025  

[2385]  3026  int iPass; 

 3027  // Setup two passes 

 3028  int start[4]; 

 3029  start[1] = number; 

 3030  start[2] = 0; 

 3031  double dstart = static_cast< double >(number) * model_>randomNumberGenerator()>randomDouble(); 

 3032  start[0] = static_cast< int >(dstart); 

 3033  start[3] = start[0]; 

 3034  //double largestWeight=0.0; 

 3035  //double smallestWeight=1.0e100; 

 3036  for (iPass = 0; iPass < 2; iPass++) { 

 3037  int end = start[2 * iPass + 1]; 

 3038  if (switchType < 5) { 

 3039  for (i = start[2 * iPass]; i < end; i++) { 

 3040  iSequence = index[i]; 

 3041  double value = infeas[iSequence]; 

 3042  if (value > tolerance) { 

 3043  double weight = weights_[iSequence]; 

 3044  //weight=1.0; 

 3045  if (value > bestDj * weight) { 

 3046  // check flagged variable and correct dj 

 3047  if (!model_>flagged(iSequence)) { 

 3048  bestDj = value / weight; 

 3049  bestSequence = iSequence; 

 3050  } else { 

 3051  // just to make sure we don't exit before got something 

 3052  numberWanted++; 

 3053  } 

[1525]  3054  } 

[2385]  3055  } 

 3056  numberWanted; 

 3057  if (!numberWanted) 

 3058  break; 

 3059  } 

 3060  } else { 

 3061  // Dantzig 

 3062  for (i = start[2 * iPass]; i < end; i++) { 

 3063  iSequence = index[i]; 

 3064  double value = infeas[iSequence]; 

 3065  if (value > tolerance) { 

 3066  if (value > bestDj) { 

 3067  // check flagged variable and correct dj 

 3068  if (!model_>flagged(iSequence)) { 

 3069  bestDj = value; 

 3070  bestSequence = iSequence; 

 3071  } else { 

 3072  // just to make sure we don't exit before got something 

 3073  numberWanted++; 

 3074  } 

 3075  } 

 3076  } 

 3077  numberWanted; 

 3078  if (!numberWanted) 

 3079  break; 

 3080  } 

 3081  } 

 3082  if (!numberWanted) 

 3083  break; 

 3084  } 

 3085  if (sequenceOut >= 0) { 

 3086  infeas[sequenceOut] = saveOutInfeasibility; 

 3087  } 

 3088  /*if (model_>numberIterations()%100==0) 

[1525]  3089  printf("%d best %g\n",bestSequence,bestDj);*/ 

[2385]  3090  reducedCost = model_>djRegion(); 

 3091  model_>clpMatrix()>setSavedBestSequence(bestSequence); 

 3092  if (bestSequence >= 0) 

 3093  model_>clpMatrix()>setSavedBestDj(reducedCost[bestSequence]); 

[1525]  3094  

[2]  3095  #ifdef CLP_DEBUG 

[2385]  3096  if (bestSequence >= 0) { 

 3097  if (model_>getStatus(bestSequence) == ClpSimplex::atLowerBound) 

 3098  assert(model_>reducedCost(bestSequence) < 0.0); 

 3099  if (model_>getStatus(bestSequence) == ClpSimplex::atUpperBound) 

 3100  assert(model_>reducedCost(bestSequence) > 0.0); 

 3101  } 

[2]  3102  #endif 

[2385]  3103  return bestSequence; 

[2]  3104  } 

[618]  3105  // Called when maximum pivots changes 

[2385]  3106  void ClpPrimalColumnSteepest::maximumPivotsChanged() 

[618]  3107  { 

[2385]  3108  if (alternateWeights_ && alternateWeights_>capacity() != model_>numberRows() + model_>factorization()>maximumPivots()) { 

 3109  delete alternateWeights_; 

 3110  alternateWeights_ = new CoinIndexedVector(); 

 3111  // enough space so can use it for factorization 

 3112  alternateWeights_>reserve(model_>numberRows() + model_>factorization()>maximumPivots()); 

 3113  } 

[618]  3114  } 

[2385]  3115  void ClpPrimalColumnSteepest::redoInfeasibilities() 

[2235]  3116  { 

[2385]  3117  double *COIN_RESTRICT infeas = infeasible_>denseVector(); 

 3118  int *COIN_RESTRICT index = infeasible_>getIndices(); 

[2235]  3119  double tolerance = model_>currentDualTolerance(); 

 3120  // we can't really trust infeasibilities if there is dual error 

 3121  // this coding has to mimic coding in checkDualSolution 

 3122  double error = CoinMin(1.0e2, model_>largestDualError()); 

 3123  // allow tolerance at least slightly bigger than standard 

[2385]  3124  tolerance = tolerance + error; 

[2235]  3125  // reverse sign so test is cleaner 

[2385]  3126  tolerance = tolerance; 

[2235]  3127  int number = model_>numberRows() + model_>numberColumns(); 

[2385]  3128  int numberNonZero = 0; 

 3129  const double *COIN_RESTRICT reducedCost = model_>djRegion(); 

 3130  const unsigned char *COIN_RESTRICT status = model_>statusArray(); 

 3131  for (int iSequence = 0; iSequence < number; iSequence++) { 

 3132  unsigned char thisStatus = status[iSequence] & 7; 

[2235]  3133  double value = reducedCost[iSequence]; 

[2385]  3134  infeas[iSequence] = 0.0; 

 3135  if (thisStatus == 3) { 

 3136  } else if ((thisStatus & 1) != 0) { 

[2235]  3137  // basic or fixed 

[2385]  3138  value = 0.0; 

 3139  } else if (thisStatus == 2) { 

 3140  value = value; 

[2235]  3141  } else { 

 3142  // free or superbasic 

 3143  if (fabs(value) > FREE_ACCEPT * tolerance) { 

[2385]  3144  // we are going to bias towards free (but only if reasonable) 

 3145  value = fabs(value) * FREE_BIAS; 

[2235]  3146  } else { 

[2385]  3147  value = 0.0; 

[2235]  3148  } 

 3149  } 

 3150  if (value < tolerance) { 

 3151  // store square in list 

 3152  infeas[iSequence] = value * value; 

 3153  index[numberNonZero++] = iSequence; 

 3154  } 

 3155  } 

 3156  infeasible_>setNumElements(numberNonZero); 

 3157  infeasibilitiesState_ = 0; 

 3158  } 

[1525]  3159  /* 

[2]  3160  1) before factorization 

 3161  2) after factorization 

 3162  3) just redo infeasibilities 

 3163  4) restore weights 

[1367]  3164  5) at end of values pass (so need initialization) 

[2]  3165  */ 

[2385]  3166  void ClpPrimalColumnSteepest::saveWeights(ClpSimplex *model, int mode) 

[2]  3167  { 

[2385]  3168  model_ = model; 

 3169  if (mode == 6) { 

 3170  // If incoming weight is 1.0 then return else as 5 

 3171  int sequenceIn = model_>sequenceIn(); 

 3172  assert(sequenceIn >= 0 && sequenceIn < model_>numberRows() + model_>numberColumns()); 

[2444]  3173  // possible weights_ was never set up // assert(weights_); 

 3174  if (weights_ && weights_[sequenceIn] == (mode_ != 1) ? 1.0 : 1.0 + ADD_ONE) 

[2385]  3175  return; 

 3176  else 

 3177  mode = 5; 

 3178  } 

 3179  if (mode_ == 4  mode_ == 5) { 

 3180  if (mode == 1 && !weights_) 

 3181  numberSwitched_ = 0; // Reset 

 3182  } 

 3183  // alternateWeights_ is defined as indexed but is treated oddly 

 3184  // at times 

 3185  int numberRows = model_>numberRows(); 

 3186  int numberColumns = model_>numberColumns(); 

 3187  const int *pivotVariable = model_>pivotVariable(); 

 3188  bool doInfeasibilities = true; 

 3189  if (mode == 1) { 

 3190  if (!model_>numberIterations()) 

 3191  pivotSequence_ = 1; 

[2235]  3192  #if ABOCA_LITE 

[2385]  3193  int numberThreads = abcState(); 

 3194  if (numberThreads && mode_ > 1) 

 3195  mode_ = 0; // force exact 

[2235]  3196  #endif 

[2385]  3197  if (weights_) { 

 3198  // Check if size has changed 

 3199  if (infeasible_>capacity() == numberRows + numberColumns && alternateWeights_>capacity() == numberRows + model_>factorization()>maximumPivots()) { 

 3200  //alternateWeights_>clear(); 

 3201  if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) { 

 3202  #if ALT_UPDATE_WEIGHTS != 2 

 3203  // save pivot order 

 3204  CoinMemcpyN(pivotVariable, 

 3205  numberRows, alternateWeights_>getIndices()); 

[2235]  3206  #endif 

[2385]  3207  // change from pivot row number to sequence number 

 3208  pivotSequence_ = pivotVariable[pivotSequence_]; 

 3209  } else { 

 3210  pivotSequence_ = 1; 

 3211  } 

 3212  state_ = 1; 

 3213  } else { 

 3214  // size has changed  clear everything 

 3215  delete[] weights_; 

 3216  weights_ = NULL; 

 3217  delete infeasible_; 

 3218  infeasible_ = NULL; 

 3219  delete alternateWeights_; 

 3220  alternateWeights_ = NULL; 

 3221  delete[] savedWeights_; 

 3222  savedWeights_ = NULL; 

 3223  delete[] reference_; 

 3224  reference_ = NULL; 

 3225  state_ = 1; 

 3226  pivotSequence_ = 1; 

 3227  } 

 3228  } 

 3229  } else if (mode == 2  mode == 4  mode == 5) { 

 3230  // restore 

 3231  if (!weights_  state_ == 1  mode == 5) { 

 3232  // Partial is only allowed with certain types of matrix 

 3233  if ((mode_ != 4 && mode_ != 5)  numberSwitched_  !model_>clpMatrix()>canDoPartialPricing()) { 

 3234  // initialize weights 

 3235  delete[] weights_; 

 3236  delete alternateWeights_; 

 3237  weights_ = new double[numberRows + numberColumns]; 

 3238  alternateWeights_ = new CoinIndexedVector(); 

 3239  // enough space so can use it for factorization 

 3240  alternateWeights_>reserve(numberRows + model_>factorization()>maximumPivots()); 

 3241  initializeWeights(); 

 3242  // create saved weights 

 3243  delete[] savedWeights_; 

 3244  savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns); 

 3245  // just do initialization 

 3246  mode = 3; 

 3247  } else { 

 3248  // Partial pricing 

 3249  // use region as somewhere to save nonfixed slacks 

 3250  // set up infeasibilities 

 3251  if (!infeasible_) { 

 3252  infeasible_ = new CoinIndexedVector(); 

 3253  infeasible_>reserve(numberColumns + numberRows); 

 3254  } 

 3255  infeasible_>clear(); 

 3256  int number = model_>numberRows() + model_>numberColumns(); 

 3257  int iSequence; 

 3258  int numberLook = 0; 

 3259  int *which = infeasible_>getIndices(); 

 3260  for (iSequence = model_>numberColumns(); iSequence < number; iSequence++) { 

 3261  ClpSimplex::Status status = model_>getStatus(iSequence); 

 3262  if (status != ClpSimplex::isFixed) 

 3263  which[numberLook++] = iSequence; 

 3264  } 

 3265  infeasible_>setNumElements(numberLook); 

 3266  doInfeasibilities = false; 

 3267  } 

 3268  savedPivotSequence_ = 2; 

 3269  savedSequenceOut_ = 2; 

 3270  if (pivotSequence_ < 0  pivotSequence_ >= numberRows + numberColumns) 

 3271  pivotSequence_ = 1; 

[1525]  3272  

[2385]  3273  } else { 

 3274  if (mode != 4) { 

 3275  // save 

 3276  CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_); 

 3277  savedPivotSequence_ = pivotSequence_; 

 3278  savedSequenceOut_ = model_>sequenceOut(); 

 3279  } else { 

 3280  // restore 

 3281  CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_); 

 3282  // was  but I think should not be 

 3283  //pivotSequence_= savedPivotSequence_; 

 3284  //model_>setSequenceOut(savedSequenceOut_); 

 3285  pivotSequence_ = 1; 

 3286  model_>setSequenceOut(1); 

 3287  // indices are wrong so clear by hand 

 3288  //alternateWeights_>clear(); 

 3289  CoinZeroN(alternateWeights_>denseVector(), 

 3290  alternateWeights_>capacity()); 

 3291  alternateWeights_>setNumElements(0); 

 3292  } 

 3293  } 

 3294  state_ = 0; 

 3295  // set up infeasibilities 

 3296  if (!infeasible_) { 

 3297  infeasible_ = new CoinIndexedVector(); 

 3298  infeasible_>reserve(numberColumns + numberRows); 

 3299  } 

 3300  } 

 3301  if (mode >= 2 && mode != 5) { 

 3302  if (mode != 3) { 

 3303  if (pivotSequence_ >= 0) { 

 3304  #if ALT_UPDATE_WEIGHTS != 2 

 3305  // restore pivot row 

 3306  int iRow; 

 3307  // permute alternateWeights 

 3308  double *temp = model_>rowArray(3)>denseVector(); 

 3309  ; 

 3310  double *work = alternateWeights_>denseVector(); 

 3311  int *savePivotOrder = model_>rowArray(3)>getIndices(); 

 3312  int *oldPivotOrder = alternateWeights_>getIndices(); 

 3313  for (iRow = 0; iRow < numberRows; iRow++) { 

 3314  int iPivot = oldPivotOrder[iRow]; 

 3315  temp[iPivot] = work[iRow]; 

 3316  savePivotOrder[iRow] = iPivot; 

 3317  } 

 3318  int number = 0; 

 3319  int found = 1; 

 3320  int *which = oldPivotOrder; 

 3321  // find pivot row and recreate alternateWeights 

 3322  for (iRow = 0; iRow < numberRows; iRow++) { 

 3323  int iPivot = pivotVariable[iRow]; 

 3324  if (iPivot == pivotSequence_) 

 3325  found = iRow; 

 3326  work[iRow] = temp[iPivot]; 

 3327  if (work[iRow]) 

 3328  which[number++] = iRow; 

 3329  } 

 3330  alternateWeights_>setNumElements(number); 

[50]  3331  #ifdef CLP_DEBUG 

[2385]  3332  // Can happen but I should clean up 

 3333  assert(found >= 0); 

[50]  3334  #endif 

[2385]  3335  pivotSequence_ = found; 

 3336  for (iRow = 0; iRow < numberRows; iRow++) { 

 3337  int iPivot = savePivotOrder[iRow]; 

 3338  temp[iPivot] = 0.0; 

 3339  } 

[2235]  3340  #else 

[2385]  3341  for (int iRow = 0; iRow < numberRows; iRow++) { 

 3342  int iPivot = pivotVariable[iRow]; 

 3343  if (iPivot == pivotSequence_) { 

 3344  pivotSequence_ = iRow; 

 3345  break; 

 3346  } 

 3347  } 

[2235]  3348  #endif 

[2385]  3349  } else { 

 3350  // Just clean up 

 3351  /* If this happens when alternateWeights_ is 

[2281]  3352  in "save" mode then alternateWeights_>clear() 

 3353  is disastrous. 

 3354  As will be fairly dense anyway and this 

 3355  rarely happens just zero out */ 

[2385]  3356  if (alternateWeights_ && alternateWeights_>getNumElements()) { 

 3357  //alternateWeights_>clear(); 

 3358  CoinZeroN(alternateWeights_>denseVector(), 

 3359  alternateWeights_>capacity()); 

 3360  alternateWeights_>setNumElements(0); 

 3361  } 

 3362  } 

 3363  } 

 3364  // Save size of factorization 

 3365  if (!model_>factorization()>pivots()) 

 3366  sizeFactorization_ = model_>factorization()>numberElements(); 

 3367  if (!doInfeasibilities) 

 3368  return; // don't disturb infeasibilities 

 3369  double *COIN_RESTRICT infeas = infeasible_>denseVector(); 

 3370  int *COIN_RESTRICT index = infeasible_>getIndices(); 

 3371  int numberNonZero = 0; 

 3372  infeasibilitiesState_ = 0; 

 3373  double tolerance = model_>currentDualTolerance(); 

 3374  int number = model_>numberRows() + model_>numberColumns(); 

 3375  int iSequence; 

[225]  3376  

[2385]  3377  double *reducedCost = model_>djRegion(); 

 3378  const double *lower = model_>lowerRegion(); 

 3379  const double *upper = model_>upperRegion(); 

 3380  const double *solution = model_>solutionRegion(); 

 3381  double primalTolerance = model_>currentPrimalTolerance(); 

[225]  3382  

[2385]  3383  if (!model_>nonLinearCost()>lookBothWays()) { 

 3384  const unsigned char *COIN_RESTRICT status = model_>statusArray(); 

 3385  #ifndef CLP_PRIMAL_SLACK_MULTIPLIER 

 3386  for (iSequence = 0; iSequence < number; iSequence++) { 

 3387  double value = reducedCost[iSequence]; 

 3388  infeas[iSequence] = 0.0; 

 3389  unsigned char thisStatus = status[iSequence] & 7; 

 3390  if (thisStatus == 3) { 

 3391  } else if ((thisStatus & 1) != 0) { 

 3392  // basic or fixed 

 3393  value = 0.0; 

 3394  } else if (thisStatus == 2) { 

 3395  value = value; 

 3396  } else { 

 3397  // free or superbasic 

 3398  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 3399  // we are going to bias towards free (but only if reasonable) 

 3400  value = fabs(value) * FREE_BIAS; 

 3401  } else { 

 3402  value = 0.0; 

 3403  } 

 3404  } 

 3405  if (value < tolerance) { 

 3406  infeas[iSequence] = value * value; 

 3407  index[numberNonZero++] = iSequence; 

 3408  } 

 3409  } 

[1732]  3410  #else 

[2385]  3411  // Columns 

 3412  int numberColumns = model_>numberColumns(); 

 3413  for (iSequence = 0; iSequence < numberColumns; iSequence++) { 

 3414  infeas[iSequence] = 0.0; 

 3415  double value = reducedCost[iSequence]; 

 3416  unsigned char thisStatus = status[iSequence] & 7; 

 3417  if (thisStatus == 3) { 

 3418  } else if ((thisStatus & 1) != 0) { 

 3419  // basic or fixed 

 3420  value = 0.0; 

 3421  } else if (thisStatus == 2) { 

 3422  value = value; 

 3423  } else { 

 3424  // free or superbasic 

 3425  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 3426  // check hasn't slipped through 

 3427  if (solution[iSequence] < lower[iSequence] + primalTolerance) { 

 3428  model_>setStatus(iSequence, ClpSimplex::atLowerBound); 

 3429  } else if (solution[iSequence] > upper[iSequence]  primalTolerance) { 

 3430  model_>setStatus(iSequence, ClpSimplex::atUpperBound); 

 3431  value = value; 

 3432  } else { 

 3433  // we are going to bias towards free (but only if reasonable) 

 3434  value = fabs(value) * FREE_BIAS; 

 3435  } 

 3436  } else { 

 3437  value = 0.0; 

 3438  } 

 3439  } 

 3440  if (value < tolerance) { 

 3441  infeas[iSequence] = value * value; 

 3442  index[numberNonZero++] = iSequence; 

 3443  } 

 3444  } 

 3445  // Rows 

 3446  for (; iSequence < number; iSequence++) { 

 3447  double value = reducedCost[iSequence]; 

 3448  infeas[iSequence] = 0.0; 

 3449  unsigned char thisStatus = status[iSequence] & 7; 

 3450  if (thisStatus == 3) { 

 3451  } else if ((thisStatus & 1) != 0) { 

 3452  // basic or fixed 

 3453  value = 0.0; 

 3454  } else if (thisStatus == 2) { 

 3455  value = value; 

 3456  } else { 

 3457  // free or superbasic 

 3458  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 3459  // we are going to bias towards free (but only if reasonable) 

 3460  value = fabs(value) * FREE_BIAS; 

 3461  } else { 

 3462  value = 0.0; 

 3463  } 

 3464  } 

 3465  if (value < tolerance) { 

 3466  infeas[iSequence] = value * value; 

 3467  index[numberNonZero++] = iSequence; 

 3468  } 

 3469  } 

[1732]  3470  #endif 

[2385]  3471  infeasible_>setNumElements(numberNonZero); 

 3472  } else { 

 3473  ClpNonLinearCost *nonLinear = model_>nonLinearCost(); 

 3474  infeasible_>clear(); 

 3475  // can go both ways 

 3476  for (iSequence = 0; iSequence < number; iSequence++) { 

 3477  double value = reducedCost[iSequence]; 

 3478  ClpSimplex::Status status = model_>getStatus(iSequence); 

[1525]  3479  

[2385]  3480  switch (status) { 

[1525]  3481  

[2385]  3482  case ClpSimplex::basic: 

 3483  case ClpSimplex::isFixed: 

 3484  break; 

 3485  case ClpSimplex::isFree: 

 3486  case ClpSimplex::superBasic: 

 3487  if (fabs(value) > FREE_ACCEPT * tolerance) { 

 3488  // we are going to bias towards free (but only if reasonable) 

 3489  value *= FREE_BIAS; 

 3490  // store square in list 

 3491  infeasible_>quickAdd(iSequence, value * value); 

[1525]  3492  } 

[2385]  3493  break; 

 3494  case ClpSimplex::atUpperBound: 

 3495  if (value > tolerance) { 

 3496  infeasible_>quickAdd(iSequence, value * value); 

 3497  } else { 

 3498  // look other way  change up should be negative 

 3499  value = nonLinear>changeUpInCost(iSequence); 

 3500  if (value < tolerance) { 

 3501  // store square in list 

 3502  infeasible_>quickAdd(iSequence, value * value); 

 3503  } 

 3504  } 

 3505  break; 

 3506  case ClpSimplex::atLowerBound: 

 3507  if (value < tolerance) { 

 3508  infeasible_>quickAdd(iSequence, value * value); 

 3509  } else { 

 3510  // look other way  change down should be positive 

 3511  value = nonLinear>changeDownInCost(iSequence); 

 3512  if (value > tolerance) { 

 3513  // store square in list 

 3514  infeasible_>quickAdd(iSequence, value * value); 

 3515  } 

 3516  } 

 3517  } 

 3518  } 

