Ignore:
Timestamp:
Oct 2, 2003 1:21:02 PM (17 years ago)
Author:
forrest
Message:

lots of stuff

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpPrimalColumnSteepest.cpp

    r210 r212  
    160160{
    161161  assert(model_);
     162  if (model_->nonLinearCost()->lookBothWays()||model_->algorithm()==2) {
     163    // Do old way
     164    return pivotColumnOldMethod(updates,spareRow1,spareRow2,
     165                                spareColumn1,spareColumn2);
     166  }
     167  int number=0;
     168  int * index;
     169  // dj could be very small (or even zero - take care)
     170  double dj = model_->dualIn();
     171  double tolerance=model_->currentDualTolerance();
     172  // we can't really trust infeasibilities if there is dual error
     173  // this coding has to mimic coding in checkDualSolution
     174  double error = min(1.0e-3,model_->largestDualError());
     175  // allow tolerance at least slightly bigger than standard
     176  tolerance = tolerance +  error;
     177  int pivotRow = model_->pivotRow();
     178  int anyUpdates;
     179  double * infeas = infeasible_->denseVector();
     180
     181  // Local copy of mode so can decide what to do
     182  int switchType;
     183  if (mode_==4)
     184    switchType = 5-numberSwitched_;
     185  else
     186    switchType = mode_;
     187  /* switchType -
     188     0 - all exact devex
     189     1 - all steepest
     190     2 - some exact devex
     191     3 - auto some exact devex
     192     4 - devex
     193     5 - dantzig
     194  */
     195  if (updates->getNumElements()) {
     196    // would have to have two goes for devex, three for steepest
     197    anyUpdates=2;
     198    // add in pivot contribution
     199    if (pivotRow>=0)
     200      updates->add(pivotRow,-dj);
     201  } else if (pivotRow>=0) {
     202    if (fabs(dj)>1.0e-15) {
     203      // some dj
     204      updates->insert(pivotRow,-dj);
     205      if (fabs(dj)>1.0e-6) {
     206        // reasonable size
     207        anyUpdates=1;
     208      } else {
     209        // too small
     210        anyUpdates=2;
     211      }
     212    } else {
     213      // zero dj
     214      anyUpdates=-1;
     215    }
     216  } else if (pivotSequence_>=0){
     217    // just after re-factorization
     218    anyUpdates=-1;
     219  } else {
     220    // sub flip - nothing to do
     221    anyUpdates=0;
     222  }
     223  int sequenceOut = model_->sequenceOut();
     224  if (switchType==5) {
     225    if (anyUpdates>0) {
     226      justDjs(updates,spareRow1,spareRow2,
     227        spareColumn1,spareColumn2);
     228    }
     229  } else if (anyUpdates==1) {
     230    if (switchType<4) {
     231      // exact etc when can use dj
     232      djsAndSteepest(updates,spareRow1,spareRow2,
     233        spareColumn1,spareColumn2);
     234    } else {
     235      // devex etc when can use dj
     236      djsAndDevex(updates,spareRow1,spareRow2,
     237        spareColumn1,spareColumn2);
     238    }
     239  } else if (anyUpdates==-1) {
     240    if (switchType<4) {
     241      // exact etc when djs okay
     242      justSteepest(updates,spareRow1,spareRow2,
     243        spareColumn1,spareColumn2);
     244    } else {
     245      // devex etc when djs okay
     246      justDevex(updates,spareRow1,spareRow2,
     247        spareColumn1,spareColumn2);
     248    }
     249  } else if (anyUpdates==2) {
     250    if (switchType<4) {
     251      // exact etc when have to use pivot
     252      djsAndSteepest2(updates,spareRow1,spareRow2,
     253        spareColumn1,spareColumn2);
     254    } else {
     255      // devex etc when have to use pivot
     256      djsAndDevex2(updates,spareRow1,spareRow2,
     257        spareColumn1,spareColumn2);
     258    }
     259  }
     260  // make sure outgoing from last iteration okay
     261  if (sequenceOut>=0) {
     262    ClpSimplex::Status status = model_->getStatus(sequenceOut);
     263    double value = model_->reducedCost(sequenceOut);
     264   
     265    switch(status) {
     266     
     267    case ClpSimplex::basic:
     268    case ClpSimplex::isFixed:
     269      break;
     270    case ClpSimplex::isFree:
     271    case ClpSimplex::superBasic:
     272      if (fabs(value)>FREE_ACCEPT*tolerance) {
     273        // we are going to bias towards free (but only if reasonable)
     274        value *= FREE_BIAS;
     275        // store square in list
     276        if (infeas[sequenceOut])
     277          infeas[sequenceOut] = value*value; // already there
     278        else
     279          infeasible_->quickAdd(sequenceOut,value*value);
     280      } else {
     281        infeasible_->zero(sequenceOut);
     282      }
     283      break;
     284    case ClpSimplex::atUpperBound:
     285      if (value>tolerance) {
     286        // store square in list
     287        if (infeas[sequenceOut])
     288          infeas[sequenceOut] = value*value; // already there
     289        else
     290          infeasible_->quickAdd(sequenceOut,value*value);
     291      } else {
     292        infeasible_->zero(sequenceOut);
     293      }
     294      break;
     295    case ClpSimplex::atLowerBound:
     296      if (value<-tolerance) {
     297        // store square in list
     298        if (infeas[sequenceOut])
     299          infeas[sequenceOut] = value*value; // already there
     300        else
     301          infeasible_->quickAdd(sequenceOut,value*value);
     302      } else {
     303        infeasible_->zero(sequenceOut);
     304      }
     305    }
     306  }
     307  // update of duals finished - now do pricing
     308  // See what sort of pricing
     309  int numberWanted=0;
     310  number = infeasible_->getNumElements();
     311  int numberColumns = model_->numberColumns();
     312  if (switchType==5) {
     313    pivotSequence_=-1;
     314    pivotRow=-1;
     315    // See if to switch
     316    int numberElements = model_->factorization()->numberElements();
     317    int numberRows = model_->numberRows();
     318    // ratio is done on number of columns here
     319    //double ratio = (double) numberElements/(double) numberColumns;
     320    double ratio = (double) numberElements/(double) numberRows;
     321    //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
     322    if (ratio<0.1) {
     323      numberWanted = max(100,number/200);
     324    } else if (ratio<0.15) {
     325      numberWanted = max(500,number/40);
     326    } else if (ratio<0.2) {
     327      numberWanted = max(2000,number/10);
     328      numberWanted = max(numberWanted,numberColumns/30);
     329    } else {
     330      switchType=4;
     331      // initialize
     332      numberSwitched_++;
     333      // Make sure will re-do
     334      delete [] weights_;
     335      weights_=NULL;
     336      saveWeights(model_,4);
     337      printf("switching to devex %d nel ratio %g\n",numberElements,ratio);
     338    }
     339    if (model_->numberIterations()%1000==0)
     340      printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     341  }
     342  if(switchType==4) {
     343    // Still in devex mode
     344    int numberElements = model_->factorization()->numberElements();
     345    int numberRows = model_->numberRows();
     346    // ratio is done on number of rows here
     347    double ratio = (double) numberElements/(double) numberRows;
     348    // Go to steepest if lot of iterations?
     349    if (ratio<1.0) {
     350      numberWanted = max(2000,number/20);
     351    } else if (ratio<2.0) {
     352      numberWanted = max(2000,number/10);
     353      numberWanted = max(numberWanted,numberColumns/20);
     354    } else {
     355      // we can zero out
     356      updates->clear();
     357      spareColumn1->clear();
     358      switchType=3;
     359      // initialize
     360      pivotSequence_=-1;
     361      pivotRow=-1;
     362      numberSwitched_++;
     363      // Make sure will re-do
     364      delete [] weights_;
     365      weights_=NULL;
     366      saveWeights(model_,4);
     367      printf("switching to exact %d nel ratio %g\n",numberElements,ratio);
     368      updates->clear();
     369    }
     370    if (model_->numberIterations()%1000==0)
     371      printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     372  }
     373  if (switchType<4) {
     374    if (switchType<2 ) {
     375      numberWanted = number+1;
     376    } else if (switchType==2) {
     377      numberWanted = max(2000,number/8);
     378    } else {
     379      int numberElements = model_->factorization()->numberElements();
     380      double ratio = (double) numberElements/(double) model_->numberRows();
     381      if (ratio<1.0) {
     382        numberWanted = max(2000,number/20);
     383      } else if (ratio<5.0) {
     384        numberWanted = max(2000,number/10);
     385        numberWanted = max(numberWanted,numberColumns/20);
     386      } else if (ratio<10.0) {
     387        numberWanted = max(2000,number/8);
     388        numberWanted = max(numberWanted,numberColumns/20);
     389      } else {
     390        ratio = number * (ratio/80.0);
     391        if (ratio>number) {
     392          numberWanted=number+1;
     393        } else {
     394          numberWanted = max(2000,(int) ratio);
     395          numberWanted = max(numberWanted,numberColumns/10);
     396        }
     397      }
     398    }
     399  }
     400
     401
     402  double bestDj = 1.0e-30;
     403  int bestSequence=-1;
     404
     405  int i,iSequence;
     406  index = infeasible_->getIndices();
     407  number = infeasible_->getNumElements();
     408  // Re-sort infeasible every 100 pivots
     409  if (0&&model_->factorization()->pivots()>0&&
     410      (model_->factorization()->pivots()%100)==0) {
     411    int nLook = model_->numberRows()+numberColumns;
     412    number=0;
     413    for (i=0;i<nLook;i++) {
     414      if (infeas[i]) {
     415        if (fabs(infeas[i])>COIN_INDEXED_TINY_ELEMENT)
     416          index[number++]=i;
     417        else
     418          infeas[i]=0.0;
     419      }
     420    }
     421    infeasible_->setNumElements(number);
     422  }
     423  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     424    // we can't really trust infeasibilities if there is dual error
     425    double checkTolerance = 1.0e-8;
     426    if (!model_->factorization()->pivots())
     427      checkTolerance = 1.0e-6;
     428    if (model_->largestDualError()>checkTolerance)
     429      tolerance *= model_->largestDualError()/checkTolerance;
     430    // But cap
     431    tolerance = min(1000.0,tolerance);
     432  }
     433#ifdef CLP_DEBUG
     434  if (model_->numberDualInfeasibilities()==1)
     435    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
     436           model_->largestDualError(),model_,model_->messageHandler(),
     437           number);
     438#endif
     439  // stop last one coming immediately
     440  double saveOutInfeasibility=0.0;
     441  if (sequenceOut>=0) {
     442    saveOutInfeasibility = infeas[sequenceOut];
     443    infeas[sequenceOut]=0.0;
     444  }
     445  tolerance *= tolerance; // as we are using squares
     446
     447  int iPass;
     448  // Setup two passes
     449  int start[4];
     450  start[1]=number;
     451  start[2]=0;
     452  double dstart = ((double) number) * CoinDrand48();
     453  start[0]=(int) dstart;
     454  start[3]=start[0];
     455  //double largestWeight=0.0;
     456  //double smallestWeight=1.0e100;
     457  for (iPass=0;iPass<2;iPass++) {
     458    int end = start[2*iPass+1];
     459    if (switchType<5) {
     460      for (i=start[2*iPass];i<end;i++) {
     461        iSequence = index[i];
     462        double value = infeas[iSequence];
     463        if (value>tolerance) {
     464          double weight = weights_[iSequence];
     465          //weight=1.0;
     466          if (value>bestDj*weight) {
     467            // check flagged variable and correct dj
     468            if (!model_->flagged(iSequence)) {
     469              bestDj=value/weight;
     470              bestSequence = iSequence;
     471            } else {
     472              // just to make sure we don't exit before got something
     473              numberWanted++;
     474            }
     475          }
     476        }
     477        numberWanted--;
     478        if (!numberWanted)
     479          break;
     480      }
     481    } else {
     482      // Dantzig
     483      for (i=start[2*iPass];i<end;i++) {
     484        iSequence = index[i];
     485        double value = infeas[iSequence];
     486        if (value>tolerance) {
     487          if (value>bestDj) {
     488            // check flagged variable and correct dj
     489            if (!model_->flagged(iSequence)) {
     490              bestDj=value;
     491              bestSequence = iSequence;
     492            } else {
     493              // just to make sure we don't exit before got something
     494              numberWanted++;
     495            }
     496          }
     497        }
     498        numberWanted--;
     499        if (!numberWanted)
     500          break;
     501      }
     502    }
     503    if (!numberWanted)
     504      break;
     505  }
     506  if (sequenceOut>=0) {
     507    infeas[sequenceOut]=saveOutInfeasibility;
     508  }
     509  /*if (model_->numberIterations()%100==0)
     510    printf("%d best %g\n",bestSequence,bestDj);*/
     511 
     512#ifdef CLP_DEBUG
     513  if (bestSequence>=0) {
     514    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
     515      assert(model_->reducedCost(bestSequence)<0.0);
     516    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
     517      assert(model_->reducedCost(bestSequence)>0.0);
     518  }
     519#endif
     520  return bestSequence;
     521}
     522// Just update djs
     523void
     524ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
     525               CoinIndexedVector * spareRow1,
     526               CoinIndexedVector * spareRow2,
     527               CoinIndexedVector * spareColumn1,
     528               CoinIndexedVector * spareColumn2)
     529{
    162530  int iSection,j;
    163   int number;
     531  int number=0;
     532  int * index;
     533  double * updateBy;
     534  double * reducedCost;
     535  double tolerance=model_->currentDualTolerance();
     536  // we can't really trust infeasibilities if there is dual error
     537  // this coding has to mimic coding in checkDualSolution
     538  double error = min(1.0e-3,model_->largestDualError());
     539  // allow tolerance at least slightly bigger than standard
     540  tolerance = tolerance +  error;
     541  int pivotRow = model_->pivotRow();
     542  double * infeas = infeasible_->denseVector();
     543  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     544 
     545  // put row of tableau in rowArray and columnArray
     546  model_->clpMatrix()->transposeTimes(model_,-1.0,
     547                                      updates,spareColumn2,spareColumn1);
     548  // normal
     549  for (iSection=0;iSection<2;iSection++) {
     550   
     551    reducedCost=model_->djRegion(iSection);
     552    int addSequence;
     553   
     554    if (!iSection) {
     555      number = updates->getNumElements();
     556      index = updates->getIndices();
     557      updateBy = updates->denseVector();
     558      addSequence = model_->numberColumns();
     559    } else {
     560      number = spareColumn1->getNumElements();
     561      index = spareColumn1->getIndices();
     562      updateBy = spareColumn1->denseVector();
     563      addSequence = 0;
     564    }
     565   
     566    for (j=0;j<number;j++) {
     567      int iSequence = index[j];
     568      double value = reducedCost[iSequence];
     569      value -= updateBy[iSequence];
     570      updateBy[iSequence]=0.0;
     571      reducedCost[iSequence] = value;
     572      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     573     
     574      switch(status) {
     575       
     576      case ClpSimplex::basic:
     577        infeasible_->zero(iSequence+addSequence);
     578      case ClpSimplex::isFixed:
     579        break;
     580      case ClpSimplex::isFree:
     581      case ClpSimplex::superBasic:
     582        if (fabs(value)>FREE_ACCEPT*tolerance) {
     583          // we are going to bias towards free (but only if reasonable)
     584          value *= FREE_BIAS;
     585          // store square in list
     586          if (infeas[iSequence+addSequence])
     587            infeas[iSequence+addSequence] = value*value; // already there
     588          else
     589            infeasible_->quickAdd(iSequence+addSequence,value*value);
     590        } else {
     591          infeasible_->zero(iSequence+addSequence);
     592        }
     593        break;
     594      case ClpSimplex::atUpperBound:
     595        if (value>tolerance) {
     596          // store square in list
     597          if (infeas[iSequence+addSequence])
     598            infeas[iSequence+addSequence] = value*value; // already there
     599          else
     600            infeasible_->quickAdd(iSequence+addSequence,value*value);
     601        } else {
     602          infeasible_->zero(iSequence+addSequence);
     603        }
     604        break;
     605      case ClpSimplex::atLowerBound:
     606        if (value<-tolerance) {
     607          // store square in list
     608          if (infeas[iSequence+addSequence])
     609            infeas[iSequence+addSequence] = value*value; // already there
     610          else
     611            infeasible_->quickAdd(iSequence+addSequence,value*value);
     612        } else {
     613          infeasible_->zero(iSequence+addSequence);
     614        }
     615      }
     616    }
     617  }
     618  updates->setNumElements(0);
     619  spareColumn1->setNumElements(0);
     620  if (pivotRow>=0) {
     621    // make sure infeasibility on incoming is 0.0
     622    int sequenceIn = model_->sequenceIn();
     623    infeasible_->zero(sequenceIn);
     624  }
     625}
     626// Update djs, weights for Devex
     627void
     628ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
     629               CoinIndexedVector * spareRow1,
     630               CoinIndexedVector * spareRow2,
     631               CoinIndexedVector * spareColumn1,
     632               CoinIndexedVector * spareColumn2)
     633{
     634  int j;
     635  int number=0;
     636  int * index;
     637  double * updateBy;
     638  double * reducedCost;
     639  // dj could be very small (or even zero - take care)
     640  double dj = model_->dualIn();
     641  double tolerance=model_->currentDualTolerance();
     642  // we can't really trust infeasibilities if there is dual error
     643  // this coding has to mimic coding in checkDualSolution
     644  double error = min(1.0e-3,model_->largestDualError());
     645  // allow tolerance at least slightly bigger than standard
     646  tolerance = tolerance +  error;
     647  // for weights update we use pivotSequence
     648  // unset in case sub flip
     649  assert (pivotSequence_>=0);
     650  pivotSequence_=-1;
     651  double * infeas = infeasible_->denseVector();
     652  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     653  // and we can see if reference
     654  double referenceIn=0.0;
     655  int sequenceIn = model_->sequenceIn();
     656  if (mode_!=1&&reference(sequenceIn))
     657    referenceIn=1.0;
     658  // save outgoing weight round update
     659  double outgoingWeight=0.0;
     660  int sequenceOut = model_->sequenceOut();
     661  if (sequenceOut>=0)
     662    outgoingWeight=weights_[sequenceOut];
     663   
     664  // put row of tableau in rowArray and columnArray
     665  model_->clpMatrix()->transposeTimes(model_,-1.0,
     666                                      updates,spareColumn2,spareColumn1);
     667  // update weights
     668  double * weight;
     669  int numberColumns = model_->numberColumns();
     670  double scaleFactor = -1.0/dj; // as formula is with 1.0
     671  // rows
     672  reducedCost=model_->djRegion(0);
     673  int addSequence=model_->numberColumns();;
     674   
     675  number = updates->getNumElements();
     676  index = updates->getIndices();
     677  updateBy = updates->denseVector();
     678  weight = weights_+numberColumns;
     679  // Devex
     680  for (j=0;j<number;j++) {
     681    double thisWeight;
     682    double pivot;
     683    double value3;
     684    int iSequence = index[j];
     685    double value = reducedCost[iSequence];
     686    double value2 = updateBy[iSequence];
     687    updateBy[iSequence]=0.0;
     688    value -= value2;
     689    reducedCost[iSequence] = value;
     690    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     691   
     692    switch(status) {
     693     
     694    case ClpSimplex::basic:
     695      infeasible_->zero(iSequence+addSequence);
     696    case ClpSimplex::isFixed:
     697      break;
     698    case ClpSimplex::isFree:
     699    case ClpSimplex::superBasic:
     700      thisWeight = weight[iSequence];
     701      // row has -1
     702      pivot = value2*scaleFactor;
     703      value3 = pivot * pivot*devex_;
     704      if (reference(iSequence+numberColumns))
     705        value3 += 1.0;
     706      weight[iSequence] = max(0.99*thisWeight,value3);
     707      if (fabs(value)>FREE_ACCEPT*tolerance) {
     708        // we are going to bias towards free (but only if reasonable)
     709        value *= FREE_BIAS;
     710        // store square in list
     711        if (infeas[iSequence+addSequence])
     712          infeas[iSequence+addSequence] = value*value; // already there
     713        else
     714          infeasible_->quickAdd(iSequence+addSequence,value*value);
     715      } else {
     716        infeasible_->zero(iSequence+addSequence);
     717      }
     718      break;
     719    case ClpSimplex::atUpperBound:
     720      thisWeight = weight[iSequence];
     721      // row has -1
     722      pivot = value2*scaleFactor;
     723      value3 = pivot * pivot*devex_;
     724      if (reference(iSequence+numberColumns))
     725        value3 += 1.0;
     726      weight[iSequence] = max(0.99*thisWeight,value3);
     727      if (value>tolerance) {
     728        // store square in list
     729        if (infeas[iSequence+addSequence])
     730          infeas[iSequence+addSequence] = value*value; // already there
     731        else
     732          infeasible_->quickAdd(iSequence+addSequence,value*value);
     733      } else {
     734        infeasible_->zero(iSequence+addSequence);
     735      }
     736      break;
     737    case ClpSimplex::atLowerBound:
     738      thisWeight = weight[iSequence];
     739      // row has -1
     740      pivot = value2*scaleFactor;
     741      value3 = pivot * pivot*devex_;
     742      if (reference(iSequence+numberColumns))
     743        value3 += 1.0;
     744      weight[iSequence] = max(0.99*thisWeight,value3);
     745      if (value<-tolerance) {
     746        // store square in list
     747        if (infeas[iSequence+addSequence])
     748          infeas[iSequence+addSequence] = value*value; // already there
     749        else
     750          infeasible_->quickAdd(iSequence+addSequence,value*value);
     751      } else {
     752        infeasible_->zero(iSequence+addSequence);
     753      }
     754    }
     755  }
     756 
     757  // columns
     758  weight = weights_;
     759 
     760  scaleFactor = -scaleFactor;
     761  reducedCost=model_->djRegion(1);
     762  number = spareColumn1->getNumElements();
     763  index = spareColumn1->getIndices();
     764  updateBy = spareColumn1->denseVector();
     765
     766  // Devex
     767 
     768  for (j=0;j<number;j++) {
     769    double thisWeight;
     770    double pivot;
     771    double value3;
     772    int iSequence = index[j];
     773    double value = reducedCost[iSequence];
     774    double value2 = updateBy[iSequence];
     775    value -= value2;
     776    updateBy[iSequence]=0.0;
     777    reducedCost[iSequence] = value;
     778    ClpSimplex::Status status = model_->getStatus(iSequence);
     779   
     780    switch(status) {
     781     
     782    case ClpSimplex::basic:
     783      infeasible_->zero(iSequence);
     784    case ClpSimplex::isFixed:
     785      break;
     786    case ClpSimplex::isFree:
     787    case ClpSimplex::superBasic:
     788      thisWeight = weight[iSequence];
     789      // row has -1
     790      pivot = value2*scaleFactor;
     791      value3 = pivot * pivot*devex_;
     792      if (reference(iSequence))
     793        value3 += 1.0;
     794      weight[iSequence] = max(0.99*thisWeight,value3);
     795      if (fabs(value)>FREE_ACCEPT*tolerance) {
     796        // we are going to bias towards free (but only if reasonable)
     797        value *= FREE_BIAS;
     798        // store square in list
     799        if (infeas[iSequence])
     800          infeas[iSequence] = value*value; // already there
     801        else
     802          infeasible_->quickAdd(iSequence,value*value);
     803      } else {
     804        infeasible_->zero(iSequence);
     805      }
     806      break;
     807    case ClpSimplex::atUpperBound:
     808      thisWeight = weight[iSequence];
     809      // row has -1
     810      pivot = value2*scaleFactor;
     811      value3 = pivot * pivot*devex_;
     812      if (reference(iSequence))
     813        value3 += 1.0;
     814      weight[iSequence] = max(0.99*thisWeight,value3);
     815      if (value>tolerance) {
     816        // store square in list
     817        if (infeas[iSequence])
     818          infeas[iSequence] = value*value; // already there
     819        else
     820          infeasible_->quickAdd(iSequence,value*value);
     821      } else {
     822        infeasible_->zero(iSequence);
     823      }
     824      break;
     825    case ClpSimplex::atLowerBound:
     826      thisWeight = weight[iSequence];
     827      // row has -1
     828      pivot = value2*scaleFactor;
     829      value3 = pivot * pivot*devex_;
     830      if (reference(iSequence))
     831        value3 += 1.0;
     832      weight[iSequence] = max(0.99*thisWeight,value3);
     833      if (value<-tolerance) {
     834        // store square in list
     835        if (infeas[iSequence])
     836          infeas[iSequence] = value*value; // already there
     837        else
     838          infeasible_->quickAdd(iSequence,value*value);
     839      } else {
     840        infeasible_->zero(iSequence);
     841      }
     842    }
     843  }
     844  // restore outgoing weight
     845  if (sequenceOut>=0)
     846    weights_[sequenceOut]=outgoingWeight;
     847  // make sure infeasibility on incoming is 0.0
     848  infeasible_->zero(sequenceIn);
     849  spareRow2->setNumElements(0);
     850  //#define SOME_DEBUG_1
     851#ifdef SOME_DEBUG_1
     852  // check for accuracy
     853  int iCheck=229;
     854  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     855  int numberRows = model_->numberRows();
     856  //int numberColumns = model_->numberColumns();
     857  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     858    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     859        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     860      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     861  }
     862#endif
     863  updates->setNumElements(0);
     864  spareColumn1->setNumElements(0);
     865}
     866// Update djs, weights for Steepest
     867void
     868ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates,
     869               CoinIndexedVector * spareRow1,
     870               CoinIndexedVector * spareRow2,
     871               CoinIndexedVector * spareColumn1,
     872               CoinIndexedVector * spareColumn2)
     873{
     874  int j;
     875  int number=0;
     876  int * index;
     877  double * updateBy;
     878  double * reducedCost;
     879  // dj could be very small (or even zero - take care)
     880  double dj = model_->dualIn();
     881  double tolerance=model_->currentDualTolerance();
     882  // we can't really trust infeasibilities if there is dual error
     883  // this coding has to mimic coding in checkDualSolution
     884  double error = min(1.0e-3,model_->largestDualError());
     885  // allow tolerance at least slightly bigger than standard
     886  tolerance = tolerance +  error;
     887  // for weights update we use pivotSequence
     888  // unset in case sub flip
     889  assert (pivotSequence_>=0);
     890  pivotSequence_=-1;
     891  double * infeas = infeasible_->denseVector();
     892  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     893
     894  model_->factorization()->updateColumnTranspose(spareRow2,
     895                                                 alternateWeights_);
     896  // and we can see if reference
     897  double referenceIn=0.0;
     898  int sequenceIn = model_->sequenceIn();
     899  if (mode_!=1&&reference(sequenceIn))
     900    referenceIn=1.0;
     901  // save outgoing weight round update
     902  double outgoingWeight=0.0;
     903  int sequenceOut = model_->sequenceOut();
     904  if (sequenceOut>=0)
     905    outgoingWeight=weights_[sequenceOut];
     906   
     907  // put row of tableau in rowArray and columnArray
     908  model_->clpMatrix()->transposeTimes(model_,-1.0,
     909                                      updates,spareColumn2,spareColumn1);
     910    // get subset which have nonzero tableau elements
     911    // Luckily spareRow2 is long enough (rowArray_[3])
     912    model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
     913                                              spareColumn1,
     914                                              spareRow2);
     915  // update weights
     916  double * weight;
     917  double * other = alternateWeights_->denseVector();
     918  int numberColumns = model_->numberColumns();
     919  double scaleFactor = -1.0/dj; // as formula is with 1.0
     920  // rows
     921  reducedCost=model_->djRegion(0);
     922  int addSequence=model_->numberColumns();;
     923   
     924  number = updates->getNumElements();
     925  index = updates->getIndices();
     926  updateBy = updates->denseVector();
     927  weight = weights_+numberColumns;
     928 
     929  for (j=0;j<number;j++) {
     930    double thisWeight;
     931    double pivot;
     932    double modification;
     933    double pivotSquared;
     934    int iSequence = index[j];
     935    double value = reducedCost[iSequence];
     936    double value2 = updateBy[iSequence];
     937    value -= value2;
     938    reducedCost[iSequence] = value;
     939    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     940    updateBy[iSequence]=0.0;
     941   
     942    switch(status) {
     943     
     944    case ClpSimplex::basic:
     945      infeasible_->zero(iSequence+addSequence);
     946    case ClpSimplex::isFixed:
     947      break;
     948    case ClpSimplex::isFree:
     949    case ClpSimplex::superBasic:
     950      thisWeight = weight[iSequence];
     951      // row has -1
     952      pivot = value2*scaleFactor;
     953      modification = other[iSequence];
     954      pivotSquared = pivot * pivot;
     955     
     956      thisWeight += pivotSquared * devex_ + pivot * modification;
     957      if (thisWeight<TRY_NORM) {
     958        if (mode_==1) {
     959          // steepest
     960          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     961        } else {
     962          // exact
     963          thisWeight = referenceIn*pivotSquared;
     964          if (reference(iSequence+numberColumns))
     965            thisWeight += 1.0;
     966          thisWeight = max(thisWeight,TRY_NORM);
     967        }
     968      }
     969      weight[iSequence] = thisWeight;
     970      if (fabs(value)>FREE_ACCEPT*tolerance) {
     971        // we are going to bias towards free (but only if reasonable)
     972        value *= FREE_BIAS;
     973        // store square in list
     974        if (infeas[iSequence+addSequence])
     975          infeas[iSequence+addSequence] = value*value; // already there
     976        else
     977          infeasible_->quickAdd(iSequence+addSequence,value*value);
     978      } else {
     979        infeasible_->zero(iSequence+addSequence);
     980      }
     981      break;
     982    case ClpSimplex::atUpperBound:
     983      thisWeight = weight[iSequence];
     984      // row has -1
     985      pivot = value2*scaleFactor;
     986      modification = other[iSequence];
     987      pivotSquared = pivot * pivot;
     988     
     989      thisWeight += pivotSquared * devex_ + pivot * modification;
     990      if (thisWeight<TRY_NORM) {
     991        if (mode_==1) {
     992          // steepest
     993          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     994        } else {
     995          // exact
     996          thisWeight = referenceIn*pivotSquared;
     997          if (reference(iSequence+numberColumns))
     998            thisWeight += 1.0;
     999          thisWeight = max(thisWeight,TRY_NORM);
     1000        }
     1001      }
     1002      weight[iSequence] = thisWeight;
     1003      if (value>tolerance) {
     1004        // store square in list
     1005        if (infeas[iSequence+addSequence])
     1006          infeas[iSequence+addSequence] = value*value; // already there
     1007        else
     1008          infeasible_->quickAdd(iSequence+addSequence,value*value);
     1009      } else {
     1010        infeasible_->zero(iSequence+addSequence);
     1011      }
     1012      break;
     1013    case ClpSimplex::atLowerBound:
     1014      thisWeight = weight[iSequence];
     1015      // row has -1
     1016      pivot = value2*scaleFactor;
     1017      modification = other[iSequence];
     1018      pivotSquared = pivot * pivot;
     1019     
     1020      thisWeight += pivotSquared * devex_ + pivot * modification;
     1021      if (thisWeight<TRY_NORM) {
     1022        if (mode_==1) {
     1023          // steepest
     1024          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1025        } else {
     1026          // exact
     1027          thisWeight = referenceIn*pivotSquared;
     1028          if (reference(iSequence+numberColumns))
     1029            thisWeight += 1.0;
     1030          thisWeight = max(thisWeight,TRY_NORM);
     1031        }
     1032      }
     1033      weight[iSequence] = thisWeight;
     1034      if (value<-tolerance) {
     1035        // store square in list
     1036        if (infeas[iSequence+addSequence])
     1037          infeas[iSequence+addSequence] = value*value; // already there
     1038        else
     1039          infeasible_->quickAdd(iSequence+addSequence,value*value);
     1040      } else {
     1041        infeasible_->zero(iSequence+addSequence);
     1042      }
     1043    }
     1044  }
     1045 
     1046  // columns
     1047  weight = weights_;
     1048 
     1049  scaleFactor = -scaleFactor;
     1050  reducedCost=model_->djRegion(1);
     1051  number = spareColumn1->getNumElements();
     1052  index = spareColumn1->getIndices();
     1053  updateBy = spareColumn1->denseVector();
     1054  double * updateBy2 = spareRow2->denseVector();
     1055
     1056  for (j=0;j<number;j++) {
     1057    double thisWeight;
     1058    double pivot;
     1059    double modification;
     1060    double pivotSquared;
     1061    int iSequence = index[j];
     1062    double value = reducedCost[iSequence];
     1063    double value2 = updateBy[iSequence];
     1064    updateBy[iSequence]=0.0;
     1065    value -= value2;
     1066    reducedCost[iSequence] = value;
     1067    ClpSimplex::Status status = model_->getStatus(iSequence);
     1068   
     1069    switch(status) {
     1070     
     1071    case ClpSimplex::basic:
     1072      infeasible_->zero(iSequence);
     1073    case ClpSimplex::isFixed:
     1074      updateBy2[iSequence]=0.0;
     1075      break;
     1076    case ClpSimplex::isFree:
     1077    case ClpSimplex::superBasic:
     1078      thisWeight = weight[iSequence];
     1079      pivot = value2*scaleFactor;
     1080      modification = updateBy2[iSequence];
     1081      updateBy2[iSequence]=0.0;
     1082      pivotSquared = pivot * pivot;
     1083     
     1084      thisWeight += pivotSquared * devex_ + pivot * modification;
     1085      if (thisWeight<TRY_NORM) {
     1086        if (mode_==1) {
     1087          // steepest
     1088          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1089        } else {
     1090          // exact
     1091          thisWeight = referenceIn*pivotSquared;
     1092          if (reference(iSequence))
     1093            thisWeight += 1.0;
     1094          thisWeight = max(thisWeight,TRY_NORM);
     1095        }
     1096      }
     1097      weight[iSequence] = thisWeight;
     1098      if (fabs(value)>FREE_ACCEPT*tolerance) {
     1099        // we are going to bias towards free (but only if reasonable)
     1100        value *= FREE_BIAS;
     1101        // store square in list
     1102        if (infeas[iSequence])
     1103          infeas[iSequence] = value*value; // already there
     1104        else
     1105          infeasible_->quickAdd(iSequence,value*value);
     1106      } else {
     1107        infeasible_->zero(iSequence);
     1108      }
     1109      break;
     1110    case ClpSimplex::atUpperBound:
     1111      thisWeight = weight[iSequence];
     1112      pivot = value2*scaleFactor;
     1113      modification = updateBy2[iSequence];
     1114      updateBy2[iSequence]=0.0;
     1115      pivotSquared = pivot * pivot;
     1116     
     1117      thisWeight += pivotSquared * devex_ + pivot * modification;
     1118      if (thisWeight<TRY_NORM) {
     1119        if (mode_==1) {
     1120          // steepest
     1121          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1122        } else {
     1123          // exact
     1124          thisWeight = referenceIn*pivotSquared;
     1125          if (reference(iSequence))
     1126            thisWeight += 1.0;
     1127          thisWeight = max(thisWeight,TRY_NORM);
     1128        }
     1129      }
     1130      weight[iSequence] = thisWeight;
     1131      if (value>tolerance) {
     1132        // store square in list
     1133        if (infeas[iSequence])
     1134          infeas[iSequence] = value*value; // already there
     1135        else
     1136          infeasible_->quickAdd(iSequence,value*value);
     1137      } else {
     1138        infeasible_->zero(iSequence);
     1139      }
     1140      break;
     1141    case ClpSimplex::atLowerBound:
     1142      thisWeight = weight[iSequence];
     1143      pivot = value2*scaleFactor;
     1144      modification = updateBy2[iSequence];
     1145      updateBy2[iSequence]=0.0;
     1146      pivotSquared = pivot * pivot;
     1147     
     1148      thisWeight += pivotSquared * devex_ + pivot * modification;
     1149      if (thisWeight<TRY_NORM) {
     1150        if (mode_==1) {
     1151          // steepest
     1152          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1153        } else {
     1154          // exact
     1155          thisWeight = referenceIn*pivotSquared;
     1156          if (reference(iSequence))
     1157            thisWeight += 1.0;
     1158          thisWeight = max(thisWeight,TRY_NORM);
     1159        }
     1160      }
     1161      weight[iSequence] = thisWeight;
     1162      if (value<-tolerance) {
     1163        // store square in list
     1164        if (infeas[iSequence])
     1165          infeas[iSequence] = value*value; // already there
     1166        else
     1167          infeasible_->quickAdd(iSequence,value*value);
     1168      } else {
     1169        infeasible_->zero(iSequence);
     1170      }
     1171    }
     1172  }
     1173  // restore outgoing weight
     1174  if (sequenceOut>=0)
     1175    weights_[sequenceOut]=outgoingWeight;
     1176  // make sure infeasibility on incoming is 0.0
     1177  infeasible_->zero(sequenceIn);
     1178  alternateWeights_->clear();
     1179  spareRow2->setNumElements(0);
     1180  //#define SOME_DEBUG_1
     1181#ifdef SOME_DEBUG_1
     1182  // check for accuracy
     1183  int iCheck=229;
     1184  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1185  int numberRows = model_->numberRows();
     1186  //int numberColumns = model_->numberColumns();
     1187  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     1188    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     1189        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     1190      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     1191  }
     1192#endif
     1193  updates->setNumElements(0);
     1194  spareColumn1->setNumElements(0);
     1195}
     1196// Update djs, weights for Devex
     1197void
     1198ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
     1199               CoinIndexedVector * spareRow1,
     1200               CoinIndexedVector * spareRow2,
     1201               CoinIndexedVector * spareColumn1,
     1202               CoinIndexedVector * spareColumn2)
     1203{
     1204  int iSection,j;
     1205  int number=0;
    1641206  int * index;
    1651207  double * updateBy;
     
    1741216  tolerance = tolerance +  error;
    1751217  int pivotRow = model_->pivotRow();
    176   // Local copy of mode so can decide what to do
    177   int switchType = mode_;
    178   if (mode_==4) {
    179     if (numberSwitched_) {
    180       switchType=3;
     1218  double * infeas = infeasible_->denseVector();
     1219  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1220 
     1221  // put row of tableau in rowArray and columnArray
     1222  model_->clpMatrix()->transposeTimes(model_,-1.0,
     1223                                      updates,spareColumn2,spareColumn1);
     1224  // normal
     1225  for (iSection=0;iSection<2;iSection++) {
     1226   
     1227    reducedCost=model_->djRegion(iSection);
     1228    int addSequence;
     1229   
     1230    if (!iSection) {
     1231      number = updates->getNumElements();
     1232      index = updates->getIndices();
     1233      updateBy = updates->denseVector();
     1234      addSequence = model_->numberColumns();
    1811235    } else {
    182       // See if to switch
    183       int numberElements = model_->factorization()->numberElements();
    184       int numberRows = model_->numberRows();
    185       // ratio is done on number of columns here
    186       int numberColumns = model_->numberColumns();
    187       double ratio = (double) numberElements/(double) numberColumns;
    188       //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
    189       int numberWanted;
    190       bool doPartial=true;
    191       int numberTotal = numberColumns+model_->numberRows();
    192       if (ratio<0.1) {
    193         numberWanted = max(100,numberTotal/200);
    194       } else if (ratio<1.0) {
    195         numberWanted = max(100,numberTotal/20);
    196       } else if (ratio<2.0) {
    197         numberWanted = max(200,numberTotal/10);
     1236      number = spareColumn1->getNumElements();
     1237      index = spareColumn1->getIndices();
     1238      updateBy = spareColumn1->denseVector();
     1239      addSequence = 0;
     1240    }
     1241   
     1242    for (j=0;j<number;j++) {
     1243      int iSequence = index[j];
     1244      double value = reducedCost[iSequence];
     1245      value -= updateBy[iSequence];
     1246      reducedCost[iSequence] = value;
     1247      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     1248     
     1249      switch(status) {
     1250       
     1251      case ClpSimplex::basic:
     1252        infeasible_->zero(iSequence+addSequence);
     1253      case ClpSimplex::isFixed:
     1254        break;
     1255      case ClpSimplex::isFree:
     1256      case ClpSimplex::superBasic:
     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        if (value>tolerance) {
     1271          // store square in list
     1272          if (infeas[iSequence+addSequence])
     1273            infeas[iSequence+addSequence] = value*value; // already there
     1274          else
     1275            infeasible_->quickAdd(iSequence+addSequence,value*value);
     1276        } else {
     1277          infeasible_->zero(iSequence+addSequence);
     1278        }
     1279        break;
     1280      case ClpSimplex::atLowerBound:
     1281        if (value<-tolerance) {
     1282          // store square in list
     1283          if (infeas[iSequence+addSequence])
     1284            infeas[iSequence+addSequence] = value*value; // already there
     1285          else
     1286            infeasible_->quickAdd(iSequence+addSequence,value*value);
     1287        } else {
     1288          infeasible_->zero(iSequence+addSequence);
     1289        }
     1290      }
     1291    }
     1292  }
     1293  // we can zero out as will have to get pivot row
     1294  updates->clear();
     1295  spareColumn1->clear();
     1296  // make sure infeasibility on incoming is 0.0
     1297  int sequenceIn = model_->sequenceIn();
     1298  infeasible_->zero(sequenceIn);
     1299  // for weights update we use pivotSequence
     1300  assert (pivotSequence_>=0);
     1301  pivotRow = pivotSequence_;
     1302  // unset in case sub flip
     1303  pivotSequence_=-1;
     1304  // make sure infeasibility on incoming is 0.0
     1305  const int * pivotVariable = model_->pivotVariable();
     1306  sequenceIn = pivotVariable[pivotRow];
     1307  infeasible_->zero(sequenceIn);
     1308  // and we can see if reference
     1309  double referenceIn=0.0;
     1310  if (mode_!=1&&reference(sequenceIn))
     1311    referenceIn=1.0;
     1312  // save outgoing weight round update
     1313  double outgoingWeight=0.0;
     1314  int sequenceOut = model_->sequenceOut();
     1315  if (sequenceOut>=0)
     1316    outgoingWeight=weights_[sequenceOut];
     1317  // update weights
     1318  updates->setNumElements(0);
     1319  spareColumn1->setNumElements(0);
     1320  // might as well set dj to 1
     1321  dj = 1.0;
     1322  updates->insert(pivotRow,-dj);
     1323  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1324  // put row of tableau in rowArray and columnArray
     1325  model_->clpMatrix()->transposeTimes(model_,-1.0,
     1326                                      updates,spareColumn2,spareColumn1);
     1327  double * weight;
     1328  int numberColumns = model_->numberColumns();
     1329  double scaleFactor = -1.0/dj; // as formula is with 1.0
     1330  // rows
     1331  number = updates->getNumElements();
     1332  index = updates->getIndices();
     1333  updateBy = updates->denseVector();
     1334  weight = weights_+numberColumns;
     1335 
     1336  assert (devex_>0.0);
     1337  for (j=0;j<number;j++) {
     1338    int iSequence = index[j];
     1339    double thisWeight = weight[iSequence];
     1340    // row has -1
     1341    double pivot = updateBy[iSequence]*scaleFactor;
     1342    updateBy[iSequence]=0.0;
     1343    double value = pivot * pivot*devex_;
     1344    if (reference(iSequence+numberColumns))
     1345      value += 1.0;
     1346    weight[iSequence] = max(0.99*thisWeight,value);
     1347  }
     1348 
     1349  // columns
     1350  weight = weights_;
     1351 
     1352  scaleFactor = -scaleFactor;
     1353 
     1354  number = spareColumn1->getNumElements();
     1355  index = spareColumn1->getIndices();
     1356  updateBy = spareColumn1->denseVector();
     1357  for (j=0;j<number;j++) {
     1358    int iSequence = index[j];
     1359    double thisWeight = weight[iSequence];
     1360    // row has -1
     1361    double pivot = updateBy[iSequence]*scaleFactor;
     1362    updateBy[iSequence]=0.0;
     1363    double value = pivot * pivot*devex_;
     1364    if (reference(iSequence))
     1365      value += 1.0;
     1366    weight[iSequence] = max(0.99*thisWeight,value);
     1367  }
     1368  // restore outgoing weight
     1369  if (sequenceOut>=0)
     1370    weights_[sequenceOut]=outgoingWeight;
     1371  spareColumn2->setNumElements(0);
     1372  //#define SOME_DEBUG_1
     1373#ifdef SOME_DEBUG_1
     1374  // check for accuracy
     1375  int iCheck=229;
     1376  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1377  int numberRows = model_->numberRows();
     1378  //int numberColumns = model_->numberColumns();
     1379  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     1380    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     1381        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     1382      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     1383  }
     1384#endif
     1385  updates->setNumElements(0);
     1386  spareColumn1->setNumElements(0);
     1387}
     1388// Update djs, weights for Steepest
     1389void
     1390ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates,
     1391               CoinIndexedVector * spareRow1,
     1392               CoinIndexedVector * spareRow2,
     1393               CoinIndexedVector * spareColumn1,
     1394               CoinIndexedVector * spareColumn2)
     1395{
     1396  int iSection,j;
     1397  int number=0;
     1398  int * index;
     1399  double * updateBy;
     1400  double * reducedCost;
     1401  // dj could be very small (or even zero - take care)
     1402  double dj = model_->dualIn();
     1403  double tolerance=model_->currentDualTolerance();
     1404  // we can't really trust infeasibilities if there is dual error
     1405  // this coding has to mimic coding in checkDualSolution
     1406  double error = min(1.0e-3,model_->largestDualError());
     1407  // allow tolerance at least slightly bigger than standard
     1408  tolerance = tolerance +  error;
     1409  int pivotRow = model_->pivotRow();
     1410  double * infeas = infeasible_->denseVector();
     1411  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1412 
     1413  // put row of tableau in rowArray and columnArray
     1414  model_->clpMatrix()->transposeTimes(model_,-1.0,
     1415                                      updates,spareColumn2,spareColumn1);
     1416  // normal
     1417  for (iSection=0;iSection<2;iSection++) {
     1418   
     1419    reducedCost=model_->djRegion(iSection);
     1420    int addSequence;
     1421   
     1422    if (!iSection) {
     1423      number = updates->getNumElements();
     1424      index = updates->getIndices();
     1425      updateBy = updates->denseVector();
     1426      addSequence = model_->numberColumns();
     1427    } else {
     1428      number = spareColumn1->getNumElements();
     1429      index = spareColumn1->getIndices();
     1430      updateBy = spareColumn1->denseVector();
     1431      addSequence = 0;
     1432    }
     1433   
     1434    for (j=0;j<number;j++) {
     1435      int iSequence = index[j];
     1436      double value = reducedCost[iSequence];
     1437      value -= updateBy[iSequence];
     1438      reducedCost[iSequence] = value;
     1439      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     1440     
     1441      switch(status) {
     1442       
     1443      case ClpSimplex::basic:
     1444        infeasible_->zero(iSequence+addSequence);
     1445      case ClpSimplex::isFixed:
     1446        break;
     1447      case ClpSimplex::isFree:
     1448      case ClpSimplex::superBasic:
     1449        if (fabs(value)>FREE_ACCEPT*tolerance) {
     1450          // we are going to bias towards free (but only if reasonable)
     1451          value *= FREE_BIAS;
     1452          // store square in list
     1453          if (infeas[iSequence+addSequence])
     1454            infeas[iSequence+addSequence] = value*value; // already there
     1455          else
     1456            infeasible_->quickAdd(iSequence+addSequence,value*value);
     1457        } else {
     1458          infeasible_->zero(iSequence+addSequence);
     1459        }
     1460        break;
     1461      case ClpSimplex::atUpperBound:
     1462        if (value>tolerance) {
     1463          // store square in list
     1464          if (infeas[iSequence+addSequence])
     1465            infeas[iSequence+addSequence] = value*value; // already there
     1466          else
     1467            infeasible_->quickAdd(iSequence+addSequence,value*value);
     1468        } else {
     1469          infeasible_->zero(iSequence+addSequence);
     1470        }
     1471        break;
     1472      case ClpSimplex::atLowerBound:
     1473        if (value<-tolerance) {
     1474          // store square in list
     1475          if (infeas[iSequence+addSequence])
     1476            infeas[iSequence+addSequence] = value*value; // already there
     1477          else
     1478            infeasible_->quickAdd(iSequence+addSequence,value*value);
     1479        } else {
     1480          infeasible_->zero(iSequence+addSequence);
     1481        }
     1482      }
     1483    }
     1484  }
     1485  // we can zero out as will have to get pivot row
     1486  // ***** move
     1487  updates->clear();
     1488  spareColumn1->clear();
     1489  if (pivotRow>=0) {
     1490    // make sure infeasibility on incoming is 0.0
     1491    int sequenceIn = model_->sequenceIn();
     1492    infeasible_->zero(sequenceIn);
     1493  }
     1494  // for weights update we use pivotSequence
     1495  pivotRow = pivotSequence_;
     1496  // unset in case sub flip
     1497  pivotSequence_=-1;
     1498  if (pivotRow>=0) {
     1499    // make sure infeasibility on incoming is 0.0
     1500    const int * pivotVariable = model_->pivotVariable();
     1501    int sequenceIn = pivotVariable[pivotRow];
     1502    infeasible_->zero(sequenceIn);
     1503    // and we can see if reference
     1504    double referenceIn=0.0;
     1505    if (mode_!=1&&reference(sequenceIn))
     1506      referenceIn=1.0;
     1507    // save outgoing weight round update
     1508    double outgoingWeight=0.0;
     1509    int sequenceOut = model_->sequenceOut();
     1510    if (sequenceOut>=0)
     1511      outgoingWeight=weights_[sequenceOut];
     1512    // update weights
     1513    updates->setNumElements(0);
     1514    spareColumn1->setNumElements(0);
     1515    // might as well set dj to 1
     1516    dj = 1.0;
     1517    updates->insert(pivotRow,-dj);
     1518    model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1519    // put row of tableau in rowArray and columnArray
     1520    model_->clpMatrix()->transposeTimes(model_,-1.0,
     1521                                        updates,spareColumn2,spareColumn1);
     1522    double * weight;
     1523    double * other = alternateWeights_->denseVector();
     1524    int numberColumns = model_->numberColumns();
     1525    double scaleFactor = -1.0/dj; // as formula is with 1.0
     1526    // rows
     1527    number = updates->getNumElements();
     1528    index = updates->getIndices();
     1529    updateBy = updates->denseVector();
     1530    weight = weights_+numberColumns;
     1531   
     1532    if (mode_<4||numberSwitched_>1) {
     1533      // Exact
     1534      // now update weight update array
     1535      model_->factorization()->updateColumnTranspose(spareRow2,
     1536                                                     alternateWeights_);
     1537      for (j=0;j<number;j++) {
     1538        int iSequence = index[j];
     1539        double thisWeight = weight[iSequence];
     1540        // row has -1
     1541        double pivot = updateBy[iSequence]*scaleFactor;
     1542        updateBy[iSequence]=0.0;
     1543        double modification = other[iSequence];
     1544        double pivotSquared = pivot * pivot;
     1545       
     1546        thisWeight += pivotSquared * devex_ + pivot * modification;
     1547        if (thisWeight<TRY_NORM) {
     1548          if (mode_==1) {
     1549            // steepest
     1550            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1551          } else {
     1552            // exact
     1553            thisWeight = referenceIn*pivotSquared;
     1554            if (reference(iSequence+numberColumns))
     1555              thisWeight += 1.0;
     1556            thisWeight = max(thisWeight,TRY_NORM);
     1557          }
     1558        }
     1559        weight[iSequence] = thisWeight;
     1560      }
     1561    } else if (mode_==4) {
     1562      // Devex
     1563      assert (devex_>0.0);
     1564      for (j=0;j<number;j++) {
     1565        int iSequence = index[j];
     1566        double thisWeight = weight[iSequence];
     1567        // row has -1
     1568        double pivot = updateBy[iSequence]*scaleFactor;
     1569        updateBy[iSequence]=0.0;
     1570        double value = pivot * pivot*devex_;
     1571        if (reference(iSequence+numberColumns))
     1572          value += 1.0;
     1573        weight[iSequence] = max(0.99*thisWeight,value);
     1574      }
     1575    }
     1576   
     1577    // columns
     1578    weight = weights_;
     1579   
     1580    scaleFactor = -scaleFactor;
     1581   
     1582    number = spareColumn1->getNumElements();
     1583    index = spareColumn1->getIndices();
     1584    updateBy = spareColumn1->denseVector();
     1585    if (mode_<4||numberSwitched_>1) {
     1586      // Exact
     1587      // get subset which have nonzero tableau elements
     1588      model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
     1589                                                spareColumn1,
     1590                                                spareColumn2);
     1591      double * updateBy2 = spareColumn2->denseVector();
     1592      for (j=0;j<number;j++) {
     1593        int iSequence = index[j];
     1594        double thisWeight = weight[iSequence];
     1595        double pivot = updateBy[iSequence]*scaleFactor;
     1596        updateBy[iSequence]=0.0;
     1597        double modification = updateBy2[iSequence];
     1598        updateBy2[iSequence]=0.0;
     1599        double pivotSquared = pivot * pivot;
     1600       
     1601        thisWeight += pivotSquared * devex_ + pivot * modification;
     1602        if (thisWeight<TRY_NORM) {
     1603          if (mode_==1) {
     1604            // steepest
     1605            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1606          } else {
     1607            // exact
     1608            thisWeight = referenceIn*pivotSquared;
     1609            if (reference(iSequence))
     1610              thisWeight += 1.0;
     1611            thisWeight = max(thisWeight,TRY_NORM);
     1612          }
     1613        }
     1614        weight[iSequence] = thisWeight;
     1615      }
     1616    } else if (mode_==4) {
     1617      // Devex
     1618      for (j=0;j<number;j++) {
     1619        int iSequence = index[j];
     1620        double thisWeight = weight[iSequence];
     1621        // row has -1
     1622        double pivot = updateBy[iSequence]*scaleFactor;
     1623        updateBy[iSequence]=0.0;
     1624        double value = pivot * pivot*devex_;
     1625        if (reference(iSequence))
     1626          value += 1.0;
     1627        weight[iSequence] = max(0.99*thisWeight,value);
     1628      }
     1629    }
     1630    // restore outgoing weight
     1631    if (sequenceOut>=0)
     1632      weights_[sequenceOut]=outgoingWeight;
     1633    alternateWeights_->clear();
     1634    spareColumn2->setNumElements(0);
     1635    //#define SOME_DEBUG_1
     1636#ifdef SOME_DEBUG_1
     1637    // check for accuracy
     1638    int iCheck=229;
     1639    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1640    int numberRows = model_->numberRows();
     1641    //int numberColumns = model_->numberColumns();
     1642    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     1643      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     1644          !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     1645        checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     1646    }
     1647#endif
     1648  }
     1649  updates->setNumElements(0);
     1650  spareColumn1->setNumElements(0);
     1651}
     1652// Update weights for Devex
     1653void
     1654ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
     1655               CoinIndexedVector * spareRow1,
     1656               CoinIndexedVector * spareRow2,
     1657               CoinIndexedVector * spareColumn1,
     1658               CoinIndexedVector * spareColumn2)
     1659{
     1660  int j;
     1661  int number=0;
     1662  int * index;
     1663  double * updateBy;
     1664  // dj could be very small (or even zero - take care)
     1665  double dj = model_->dualIn();
     1666  double tolerance=model_->currentDualTolerance();
     1667  // we can't really trust infeasibilities if there is dual error
     1668  // this coding has to mimic coding in checkDualSolution
     1669  double error = min(1.0e-3,model_->largestDualError());
     1670  // allow tolerance at least slightly bigger than standard
     1671  tolerance = tolerance +  error;
     1672  int pivotRow = model_->pivotRow();
     1673  // for weights update we use pivotSequence
     1674  pivotRow = pivotSequence_;
     1675  assert (pivotRow>=0);
     1676  // make sure infeasibility on incoming is 0.0
     1677  const int * pivotVariable = model_->pivotVariable();
     1678  int sequenceIn = pivotVariable[pivotRow];
     1679  infeasible_->zero(sequenceIn);
     1680  // and we can see if reference
     1681  double referenceIn=0.0;
     1682  if (mode_!=1&&reference(sequenceIn))
     1683    referenceIn=1.0;
     1684  // save outgoing weight round update
     1685  double outgoingWeight=0.0;
     1686  int sequenceOut = model_->sequenceOut();
     1687  if (sequenceOut>=0)
     1688    outgoingWeight=weights_[sequenceOut];
     1689  assert (!updates->getNumElements());
     1690  assert (!spareColumn1->getNumElements());
     1691  // unset in case sub flip
     1692  pivotSequence_=-1;
     1693  // might as well set dj to 1
     1694  dj = 1.0;
     1695  updates->insert(pivotRow,-dj);
     1696  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1697  // put row of tableau in rowArray and columnArray
     1698  model_->clpMatrix()->transposeTimes(model_,-1.0,
     1699                                      updates,spareColumn2,spareColumn1);
     1700  double * weight;
     1701  int numberColumns = model_->numberColumns();
     1702  double scaleFactor = -1.0/dj; // as formula is with 1.0
     1703  // rows
     1704  number = updates->getNumElements();
     1705  index = updates->getIndices();
     1706  updateBy = updates->denseVector();
     1707  weight = weights_+numberColumns;
     1708 
     1709  // Devex
     1710  assert (devex_>0.0);
     1711  for (j=0;j<number;j++) {
     1712    int iSequence = index[j];
     1713    double thisWeight = weight[iSequence];
     1714    // row has -1
     1715    double pivot = updateBy[iSequence]*scaleFactor;
     1716    updateBy[iSequence]=0.0;
     1717    double value = pivot * pivot*devex_;
     1718    if (reference(iSequence+numberColumns))
     1719      value += 1.0;
     1720    weight[iSequence] = max(0.99*thisWeight,value);
     1721  }
     1722 
     1723  // columns
     1724  weight = weights_;
     1725 
     1726  scaleFactor = -scaleFactor;
     1727 
     1728  number = spareColumn1->getNumElements();
     1729  index = spareColumn1->getIndices();
     1730  updateBy = spareColumn1->denseVector();
     1731  // Devex
     1732  for (j=0;j<number;j++) {
     1733    int iSequence = index[j];
     1734    double thisWeight = weight[iSequence];
     1735    // row has -1
     1736    double pivot = updateBy[iSequence]*scaleFactor;
     1737    updateBy[iSequence]=0.0;
     1738    double value = pivot * pivot*devex_;
     1739    if (reference(iSequence))
     1740      value += 1.0;
     1741    weight[iSequence] = max(0.99*thisWeight,value);
     1742  }
     1743  // restore outgoing weight
     1744  if (sequenceOut>=0)
     1745    weights_[sequenceOut]=outgoingWeight;
     1746  spareColumn2->setNumElements(0);
     1747  //#define SOME_DEBUG_1
     1748#ifdef SOME_DEBUG_1
     1749  // check for accuracy
     1750  int iCheck=229;
     1751  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1752  int numberRows = model_->numberRows();
     1753  //int numberColumns = model_->numberColumns();
     1754  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     1755    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     1756        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     1757      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     1758  }
     1759#endif
     1760  updates->setNumElements(0);
     1761  spareColumn1->setNumElements(0);
     1762}
     1763// Update weights for Steepest
     1764void
     1765ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates,
     1766               CoinIndexedVector * spareRow1,
     1767               CoinIndexedVector * spareRow2,
     1768               CoinIndexedVector * spareColumn1,
     1769               CoinIndexedVector * spareColumn2)
     1770{
     1771  int j;
     1772  int number=0;
     1773  int * index;
     1774  double * updateBy;
     1775  // dj could be very small (or even zero - take care)
     1776  double dj = model_->dualIn();
     1777  double tolerance=model_->currentDualTolerance();
     1778  // we can't really trust infeasibilities if there is dual error
     1779  // this coding has to mimic coding in checkDualSolution
     1780  double error = min(1.0e-3,model_->largestDualError());
     1781  // allow tolerance at least slightly bigger than standard
     1782  tolerance = tolerance +  error;
     1783  int pivotRow = model_->pivotRow();
     1784  // for weights update we use pivotSequence
     1785  pivotRow = pivotSequence_;
     1786  // unset in case sub flip
     1787  pivotSequence_=-1;
     1788  assert (pivotRow>=0);
     1789  // make sure infeasibility on incoming is 0.0
     1790  const int * pivotVariable = model_->pivotVariable();
     1791  int sequenceIn = pivotVariable[pivotRow];
     1792  infeasible_->zero(sequenceIn);
     1793  // and we can see if reference
     1794  double referenceIn=0.0;
     1795  if (mode_!=1&&reference(sequenceIn))
     1796    referenceIn=1.0;
     1797  // save outgoing weight round update
     1798  double outgoingWeight=0.0;
     1799  int sequenceOut = model_->sequenceOut();
     1800  if (sequenceOut>=0)
     1801    outgoingWeight=weights_[sequenceOut];
     1802  assert (!updates->getNumElements());
     1803  assert (!spareColumn1->getNumElements());
     1804  // update weights
     1805  //updates->setNumElements(0);
     1806  //spareColumn1->setNumElements(0);
     1807  // might as well set dj to 1
     1808  dj = 1.0;
     1809  updates->insert(pivotRow,-dj);
     1810  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1811  // put row of tableau in rowArray and columnArray
     1812  model_->clpMatrix()->transposeTimes(model_,-1.0,
     1813                                      updates,spareColumn2,spareColumn1);
     1814  double * weight;
     1815  double * other = alternateWeights_->denseVector();
     1816  int numberColumns = model_->numberColumns();
     1817  double scaleFactor = -1.0/dj; // as formula is with 1.0
     1818  // rows
     1819  number = updates->getNumElements();
     1820  index = updates->getIndices();
     1821  updateBy = updates->denseVector();
     1822  weight = weights_+numberColumns;
     1823 
     1824  // Exact
     1825  // now update weight update array
     1826  model_->factorization()->updateColumnTranspose(spareRow2,
     1827                                                 alternateWeights_);
     1828  for (j=0;j<number;j++) {
     1829    int iSequence = index[j];
     1830    double thisWeight = weight[iSequence];
     1831    // row has -1
     1832    double pivot = updateBy[iSequence]*scaleFactor;
     1833    updateBy[iSequence]=0.0;
     1834    double modification = other[iSequence];
     1835    double pivotSquared = pivot * pivot;
     1836   
     1837    thisWeight += pivotSquared * devex_ + pivot * modification;
     1838    if (thisWeight<TRY_NORM) {
     1839      if (mode_==1) {
     1840        // steepest
     1841        thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
    1981842      } else {
    199         doPartial=false; // switch off
    200         numberWanted=0; // just for compliler message
    201       }
    202       if (doPartial) {
    203         if (model_->numberIterations()%100==0)
    204           printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
    205         // No do partial
    206         // Always do slacks
    207         int anyUpdates;
    208         double * dual = model_->dualRowSolution();
    209         if (updates->getNumElements()) {
    210           // would have to have two goes for devex, three for steepest
    211           anyUpdates=2;
    212           // add in pivot contribution
    213           if (pivotRow>=0)
    214             updates->add(pivotRow,-dj);
    215         } else if (pivotRow>=0) {
    216           if (fabs(dj)>1.0e-15) {
    217             // some dj
    218             updates->insert(pivotRow,-dj);
    219             if (fabs(dj)>1.0e-6) {
    220               // reasonable size
    221               anyUpdates=1;
    222             } else {
    223               // too small
    224               anyUpdates=2;
    225             }
    226           } else {
    227             // zero dj
    228             anyUpdates=-1;
    229           }
    230         } else if (pivotSequence_>=0){
    231           // just after re-factorization
    232           anyUpdates=-1;
    233         } else {
    234           // sub flip - nothing to do
    235           anyUpdates=0;
    236         }
    237 
    238         // For now (as we can always switch off)
    239         assert (!model_->nonLinearCost()->lookBothWays());
    240         double * slackDj =model_->djRegion(0);
    241         double * dj =model_->djRegion(1);
    242         if (anyUpdates>0) {
    243           model_->factorization()->updateColumnTranspose(spareRow2,updates);
     1843        // exact
     1844        thisWeight = referenceIn*pivotSquared;
     1845        if (reference(iSequence+numberColumns))
     1846          thisWeight += 1.0;
     1847        thisWeight = max(thisWeight,TRY_NORM);
     1848      }
     1849    }
     1850    weight[iSequence] = thisWeight;
     1851  }
     1852 
     1853  // columns
     1854  weight = weights_;
     1855 
     1856  scaleFactor = -scaleFactor;
     1857 
     1858  number = spareColumn1->getNumElements();
     1859  index = spareColumn1->getIndices();
     1860  updateBy = spareColumn1->denseVector();
     1861  // Exact
     1862  // get subset which have nonzero tableau elements
     1863  model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
     1864                                            spareColumn1,
     1865                                            spareColumn2);
     1866  double * updateBy2 = spareColumn2->denseVector();
     1867  for (j=0;j<number;j++) {
     1868    int iSequence = index[j];
     1869    double thisWeight = weight[iSequence];
     1870    double pivot = updateBy[iSequence]*scaleFactor;
     1871    updateBy[iSequence]=0.0;
     1872    double modification = updateBy2[iSequence];
     1873    updateBy2[iSequence]=0.0;
     1874    double pivotSquared = pivot * pivot;
    2441875   
    245           number = updates->getNumElements();
    246           index = updates->getIndices();
    247           updateBy = updates->denseVector();
    248          
    249 
    250           for (j=0;j<number;j++) {
    251             int iSequence = index[j];
    252             double value = updateBy[iSequence];
    253             slackDj[iSequence] -= value;
    254             dual[iSequence] -= value;
    255             updateBy[iSequence]=0.0;
    256           }
    257           updates->setNumElements(0);
    258         }
    259        
    260         int iPass;
    261         // Setup two passes AND treat slacks differently
    262         int start[2][4];
    263         // slacks
    264         start[0][1]=numberRows;
    265         start[0][2]=0;
    266         double dstart = ((double) numberRows) * CoinDrand48();
    267         start[0][0]=(int) dstart;
    268         start[0][3]=start[0][0];
    269         for (iPass=0;iPass<4;iPass++)
    270           start[0][iPass] += numberColumns;
    271         start[1][1]=numberColumns;
    272         start[1][2]=0;
    273         dstart = ((double) numberColumns) * CoinDrand48();
    274         start[1][0]=(int) dstart;
    275         start[1][3]=start[1][0];
    276         int bestSequence=-1;
    277         double bestDj=tolerance;
    278         // We are going to fake it so pi looks sparse
    279         double * saveDense = updates->denseVector();
    280         updates->setDenseVector(dual);
    281         int kChunk = max(2000,numberWanted);
    282         for (iPass=0;iPass<2;iPass++) {
    283           // Slacks
    284           int end = start[0][2*iPass+1];
    285           for (int iSequence=start[0][2*iPass];iSequence<end;iSequence++) {
    286             double value;
    287             ClpSimplex::Status status = model_->getStatus(iSequence);
    288          
    289             switch(status) {
    290              
    291             case ClpSimplex::basic:
    292             case ClpSimplex::isFixed:
    293               break;
    294             case ClpSimplex::isFree:
    295             case ClpSimplex::superBasic:
    296               value = dj[iSequence];
    297               if (fabs(value)>FREE_ACCEPT*tolerance) {
    298                 numberWanted--;
    299                 // we are going to bias towards free (but only if reasonable)
    300                 value = fabs(value*FREE_BIAS);
    301                 if (value>bestDj) {
    302                   // check flagged variable and correct dj
    303                   if (!model_->flagged(iSequence)) {
    304                     bestDj=value;
    305                     bestSequence = iSequence;
    306                   } else {
    307                     // just to make sure we don't exit before got something
    308                     numberWanted++;
    309                   }
    310                 }
    311               }
    312               break;
    313             case ClpSimplex::atUpperBound:
    314               value = dj[iSequence];
    315               if (value>tolerance) {
    316                 numberWanted--;
    317                 if (value>bestDj) {
    318                   // check flagged variable and correct dj
    319                   if (!model_->flagged(iSequence)) {
    320                     bestDj=value;
    321                     bestSequence = iSequence;
    322                   } else {
    323                     // just to make sure we don't exit before got something
    324                     numberWanted++;
    325                   }
    326                 }
    327               }
    328               break;
    329             case ClpSimplex::atLowerBound:
    330               value = dj[iSequence];
    331               if (value<-tolerance) {
    332                 numberWanted--;
    333                 if (-value>bestDj) {
    334                   // check flagged variable and correct dj
    335                   if (!model_->flagged(iSequence)) {
    336                     bestDj=-value;
    337                     bestSequence = iSequence;
    338                   } else {
    339                     // just to make sure we don't exit before got something
    340                     numberWanted++;
    341                   }
    342                 }
    343               }
    344             }
    345             if (!numberWanted)
    346               break;
    347           }
    348           if (!numberWanted)
    349             break;
    350           // For Columns break up into manageable chunks
    351           end = start[1][2*iPass+1];
    352           index = spareColumn1->getIndices();
    353           updateBy = spareColumn1->denseVector();
    354           for (int jSequence=start[1][2*iPass];jSequence<end;jSequence+=kChunk) {
    355             int endThis = min(end,jSequence+kChunk);
    356             number=0;
    357             for (int iSequence=jSequence;iSequence<endThis;iSequence++) {
    358               ClpSimplex::Status status = model_->getStatus(iSequence);
    359          
    360               if(status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) {
    361                 index[number++]=iSequence;
    362               }
    363             }
    364             spareColumn1->setNumElements(number);
    365             // get subset which have nonzero djs
    366             // updates is actually slack djs
    367             model_->clpMatrix()->subsetTransposeTimes(model_,updates,
    368                                                       spareColumn1,
    369                                                       spareColumn2);
    370             double * updateBy2 = spareColumn2->denseVector();
    371             const double * cost = model_->costRegion();
    372             for (int j=0;j<number;j++) {
    373               int iSequence = index[j];
    374               double value = cost[iSequence]-updateBy2[iSequence];
    375               updateBy2[iSequence]=0.0;
    376               ClpSimplex::Status status = model_->getStatus(iSequence);
    377          
    378               switch(status) {
    379                
    380               case ClpSimplex::basic:
    381               case ClpSimplex::isFixed:
    382                 break;
    383               case ClpSimplex::isFree:
    384               case ClpSimplex::superBasic:
    385                 if (fabs(value)>FREE_ACCEPT*tolerance) {
    386                   numberWanted--;
    387                   // we are going to bias towards free (but only if reasonable)
    388                   if (fabs(value*FREE_BIAS)>bestDj) {
    389                     // check flagged variable and correct dj
    390                     if (!model_->flagged(iSequence)) {
    391                       bestDj=value*FREE_BIAS;
    392                       dj[iSequence]=value;
    393                       bestSequence = iSequence;
    394                     } else {
    395                       // just to make sure we don't exit before got something
    396                       numberWanted++;
    397                     }
    398                   }
    399                 }
    400                 break;
    401               case ClpSimplex::atUpperBound:
    402                 if (value>tolerance) {
    403                   numberWanted--;
    404                   if (value>bestDj) {
    405                     // check flagged variable and correct dj
    406                     if (!model_->flagged(iSequence)) {
    407                       bestDj=value;
    408                       dj[iSequence]=value;
    409                       bestSequence = iSequence;
    410                     } else {
    411                       // just to make sure we don't exit before got something
    412                       numberWanted++;
    413                     }
    414                   }
    415                 }
    416                 break;
    417               case ClpSimplex::atLowerBound:
    418                 if (value<-tolerance) {
    419                   numberWanted--;
    420                   if (-value>bestDj) {
    421                     // check flagged variable and correct dj
    422                     if (!model_->flagged(iSequence)) {
    423                       bestDj=-value;
    424                       dj[iSequence]=value;
    425                       bestSequence = iSequence;
    426                     } else {
    427                       // just to make sure we don't exit before got something
    428                       numberWanted++;
    429                     }
    430                   }
    431                 }
    432                 break;
    433               }
    434               if (!numberWanted) {
    435                 // erase rest
    436                 for (;j<number;j++) {
    437                   int iSequence = index[j];
    438                   updateBy2[iSequence]=0.0;
    439                 }
    440                 break;
    441               }
    442             }
    443             if (!numberWanted)
    444               break;
    445           }
    446           if (!numberWanted)
    447             break;
    448         }
    449         spareColumn1->setNumElements(0);
    450         spareColumn2->setNumElements(0);
    451         updates->setDenseVector(saveDense);
    452         /*if (model_->numberIterations()%100==0)
    453           printf("%d best %g\n",bestSequence,bestDj);*/
    454        
    455 #ifdef CLP_DEBUG
    456         if (bestSequence>=0) {
    457           if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
    458             assert(model_->reducedCost(bestSequence)<0.0);
    459           if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
    460             assert(model_->reducedCost(bestSequence)>0.0);
    461         }
     1876    thisWeight += pivotSquared * devex_ + pivot * modification;
     1877    if (thisWeight<TRY_NORM) {
     1878      if (mode_==1) {
     1879        // steepest
     1880        thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     1881      } else {
     1882        // exact
     1883        thisWeight = referenceIn*pivotSquared;
     1884        if (reference(iSequence))
     1885          thisWeight += 1.0;
     1886        thisWeight = max(thisWeight,TRY_NORM);
     1887      }
     1888    }
     1889    weight[iSequence] = thisWeight;
     1890  }
     1891  // restore outgoing weight
     1892  if (sequenceOut>=0)
     1893    weights_[sequenceOut]=outgoingWeight;
     1894  alternateWeights_->clear();
     1895  spareColumn2->setNumElements(0);
     1896  //#define SOME_DEBUG_1
     1897#ifdef SOME_DEBUG_1
     1898  // check for accuracy
     1899  int iCheck=229;
     1900  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1901  int numberRows = model_->numberRows();
     1902  //int numberColumns = model_->numberColumns();
     1903  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     1904    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     1905        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     1906      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     1907  }
    4621908#endif
    463         return bestSequence;
    464       } else {
    465         // Yes initialize
    466         pivotSequence_=-1;
    467         pivotRow=-1;
    468         numberSwitched_++;
    469         // Redo dualsolution
    470         model_->computeDuals(NULL);
    471         saveWeights(model_,4);
    472         printf("switching %d nel ratio %g\n",numberElements,ratio);
    473         updates->clear();
    474       }
    475     }
    476   }
     1909  updates->setNumElements(0);
     1910  spareColumn1->setNumElements(0);
     1911}
     1912// Returns pivot column, -1 if none
     1913int
     1914ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates,
     1915                                    CoinIndexedVector * spareRow1,
     1916                                    CoinIndexedVector * spareRow2,
     1917                                    CoinIndexedVector * spareColumn1,
     1918                                    CoinIndexedVector * spareColumn2)
     1919{
     1920  assert(model_);
     1921  int iSection,j;
     1922  int number=0;
     1923  int * index;
     1924  double * updateBy;
     1925  double * reducedCost;
     1926  // dj could be very small (or even zero - take care)
     1927  double dj = model_->dualIn();
     1928  double tolerance=model_->currentDualTolerance();
     1929  // we can't really trust infeasibilities if there is dual error
     1930  // this coding has to mimic coding in checkDualSolution
     1931  double error = min(1.0e-3,model_->largestDualError());
     1932  // allow tolerance at least slightly bigger than standard
     1933  tolerance = tolerance +  error;
     1934  int pivotRow = model_->pivotRow();
    4771935  int anyUpdates;
    4781936  double * infeas = infeasible_->denseVector();
    4791937
     1938  // Local copy of mode so can decide what to do
     1939  int switchType;
     1940  if (mode_==4)
     1941    switchType = 5-numberSwitched_;
     1942  else
     1943    switchType = mode_;
     1944  /* switchType -
     1945     0 - all exact devex
     1946     1 - all steepest
     1947     2 - some exact devex
     1948     3 - auto some exact devex
     1949     4 - devex
     1950     5 - dantzig
     1951  */
    4801952  if (updates->getNumElements()) {
    4811953    // would have to have two goes for devex, three for steepest
     
    8542326    }
    8552327  }
     2328  // See what sort of pricing
     2329  int numberWanted=0;
     2330  number = infeasible_->getNumElements();
     2331  int numberColumns = model_->numberColumns();
     2332  if (switchType==5) {
     2333    // we can zero out
     2334    updates->clear();
     2335    spareColumn1->clear();
     2336    pivotSequence_=-1;
     2337    pivotRow=-1;
     2338    // See if to switch
     2339    int numberElements = model_->factorization()->numberElements();
     2340    int numberRows = model_->numberRows();
     2341    // ratio is done on number of columns here
     2342    //double ratio = (double) numberElements/(double) numberColumns;
     2343    double ratio = (double) numberElements/(double) numberRows;
     2344    //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
     2345    if (ratio<0.1) {
     2346      numberWanted = max(100,number/200);
     2347    } else if (ratio<0.3) {
     2348      numberWanted = max(500,number/40);
     2349    } else if (ratio<0.5) {
     2350      numberWanted = max(2000,number/10);
     2351      numberWanted = max(numberWanted,numberColumns/30);
     2352    } else {
     2353      switchType=4;
     2354      // initialize
     2355      numberSwitched_++;
     2356      // Make sure will re-do
     2357      delete [] weights_;
     2358      weights_=NULL;
     2359      saveWeights(model_,4);
     2360      printf("switching to devex %d nel ratio %g\n",numberElements,ratio);
     2361      updates->clear();
     2362    }
     2363    if (model_->numberIterations()%1000==0)
     2364      printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     2365  }
     2366  if(switchType==4) {
     2367    // Still in devex mode
     2368    int numberElements = model_->factorization()->numberElements();
     2369    int numberRows = model_->numberRows();
     2370    // ratio is done on number of rows here
     2371    double ratio = (double) numberElements/(double) numberRows;
     2372    // Go to steepest if lot of iterations?
     2373    if (ratio<1.0) {
     2374      numberWanted = max(2000,number/20);
     2375    } else if (ratio<5.0) {
     2376      numberWanted = max(2000,number/10);
     2377      numberWanted = max(numberWanted,numberColumns/20);
     2378    } else {
     2379      // we can zero out
     2380      updates->clear();
     2381      spareColumn1->clear();
     2382      switchType=3;
     2383      // initialize
     2384      pivotSequence_=-1;
     2385      pivotRow=-1;
     2386      numberSwitched_++;
     2387      // Make sure will re-do
     2388      delete [] weights_;
     2389      weights_=NULL;
     2390      saveWeights(model_,4);
     2391      printf("switching to exact %d nel ratio %g\n",numberElements,ratio);
     2392      updates->clear();
     2393    }
     2394    if (model_->numberIterations()%1000==0)
     2395      printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     2396  }
     2397  if (switchType<4) {
     2398    if (switchType<2 ) {
     2399      numberWanted = number+1;
     2400    } else if (switchType==2) {
     2401      numberWanted = max(2000,number/8);
     2402    } else {
     2403      int numberElements = model_->factorization()->numberElements();
     2404      double ratio = (double) numberElements/(double) model_->numberRows();
     2405      if (ratio<1.0) {
     2406        numberWanted = max(2000,number/20);
     2407      } else if (ratio<5.0) {
     2408        numberWanted = max(2000,number/10);
     2409        numberWanted = max(numberWanted,numberColumns/20);
     2410      } else if (ratio<10.0) {
     2411        numberWanted = max(2000,number/8);
     2412        numberWanted = max(numberWanted,numberColumns/20);
     2413      } else {
     2414        ratio = number * (ratio/80.0);
     2415        if (ratio>number) {
     2416          numberWanted=number+1;
     2417        } else {
     2418          numberWanted = max(2000,(int) ratio);
     2419          numberWanted = max(numberWanted,numberColumns/10);
     2420        }
     2421      }
     2422    }
     2423  }
    8562424  // for weights update we use pivotSequence
    8572425  pivotRow = pivotSequence_;
     
    8832451                                          updates,spareColumn2,spareColumn1);
    8842452    }
    885     // now update weight update array
    886     model_->factorization()->updateColumnTranspose(spareRow2,
    887                                                    alternateWeights_);
     2453    double * weight;
    8882454    double * other = alternateWeights_->denseVector();
    889     double * weight;
    8902455    int numberColumns = model_->numberColumns();
    8912456    double scaleFactor = -1.0/dj; // as formula is with 1.0
     
    8952460    updateBy = updates->denseVector();
    8962461    weight = weights_+numberColumns;
    897    
    898     for (j=0;j<number;j++) {
    899       int iSequence = index[j];
    900       double thisWeight = weight[iSequence];
    901       // row has -1
    902       double pivot = updateBy[iSequence]*scaleFactor;
    903       updateBy[iSequence]=0.0;
    904       double modification = other[iSequence];
    905       double pivotSquared = pivot * pivot;
    9062462     
    907       thisWeight += pivotSquared * devex_ + pivot * modification;
    908       if (thisWeight<TRY_NORM) {
    909         if (switchType==1) {
    910           // steepest
    911           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
    912         } else {
    913           // exact
    914           thisWeight = referenceIn*pivotSquared;
    915           if (reference(iSequence+numberColumns))
    916             thisWeight += 1.0;
    917           thisWeight = max(thisWeight,TRY_NORM);
    918         }
    919       }
    920       weight[iSequence] = thisWeight;
     2463    if (switchType<4) {
     2464      // Exact
     2465      // now update weight update array
     2466      model_->factorization()->updateColumnTranspose(spareRow2,
     2467                                                     alternateWeights_);
     2468      for (j=0;j<number;j++) {
     2469        int iSequence = index[j];
     2470        double thisWeight = weight[iSequence];
     2471        // row has -1
     2472        double pivot = updateBy[iSequence]*scaleFactor;
     2473        updateBy[iSequence]=0.0;
     2474        double modification = other[iSequence];
     2475        double pivotSquared = pivot * pivot;
     2476       
     2477        thisWeight += pivotSquared * devex_ + pivot * modification;
     2478        if (thisWeight<TRY_NORM) {
     2479          if (switchType==1) {
     2480            // steepest
     2481            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     2482          } else {
     2483            // exact
     2484            thisWeight = referenceIn*pivotSquared;
     2485            if (reference(iSequence+numberColumns))
     2486              thisWeight += 1.0;
     2487            thisWeight = max(thisWeight,TRY_NORM);
     2488          }
     2489        }
     2490        weight[iSequence] = thisWeight;
     2491      }
     2492    } else if (switchType==4) {
     2493      // Devex
     2494      assert (devex_>0.0);
     2495      for (j=0;j<number;j++) {
     2496        int iSequence = index[j];
     2497        double thisWeight = weight[iSequence];
     2498        // row has -1
     2499        double pivot = updateBy[iSequence]*scaleFactor;
     2500        updateBy[iSequence]=0.0;
     2501        double value = pivot * pivot*devex_;
     2502        if (reference(iSequence+numberColumns))
     2503          value += 1.0;
     2504        weight[iSequence] = max(0.99*thisWeight,value);
     2505      }
    9212506    }
    9222507   
     
    9292514    index = spareColumn1->getIndices();
    9302515    updateBy = spareColumn1->denseVector();
    931     // get subset which have nonzero tableau elements
    932     model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
    933                                               spareColumn1,
    934                                               spareColumn2);
    935     double * updateBy2 = spareColumn2->denseVector();
    936     for (j=0;j<number;j++) {
    937       int iSequence = index[j];
    938       double thisWeight = weight[iSequence];
    939       double pivot = updateBy[iSequence]*scaleFactor;
    940       updateBy[iSequence]=0.0;
    941       double modification = updateBy2[iSequence];
    942       updateBy2[iSequence]=0.0;
    943       double pivotSquared = pivot * pivot;
    944      
    945       thisWeight += pivotSquared * devex_ + pivot * modification;
    946       if (thisWeight<TRY_NORM) {
    947         if (switchType==1) {
    948           // steepest
    949           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
    950         } else {
    951           // exact
    952           thisWeight = referenceIn*pivotSquared;
    953           if (reference(iSequence))
    954             thisWeight += 1.0;
    955           thisWeight = max(thisWeight,TRY_NORM);
    956         }
    957       }
    958       weight[iSequence] = thisWeight;
     2516    if (switchType<4) {
     2517      // Exact
     2518      // get subset which have nonzero tableau elements
     2519      model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
     2520                                                spareColumn1,
     2521                                                spareColumn2);
     2522      double * updateBy2 = spareColumn2->denseVector();
     2523      for (j=0;j<number;j++) {
     2524        int iSequence = index[j];
     2525        double thisWeight = weight[iSequence];
     2526        double pivot = updateBy[iSequence]*scaleFactor;
     2527        updateBy[iSequence]=0.0;
     2528        double modification = updateBy2[iSequence];
     2529        updateBy2[iSequence]=0.0;
     2530        double pivotSquared = pivot * pivot;
     2531       
     2532        thisWeight += pivotSquared * devex_ + pivot * modification;
     2533        if (thisWeight<TRY_NORM) {
     2534          if (switchType==1) {
     2535            // steepest
     2536            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     2537          } else {
     2538            // exact
     2539            thisWeight = referenceIn*pivotSquared;
     2540            if (reference(iSequence))
     2541              thisWeight += 1.0;
     2542            thisWeight = max(thisWeight,TRY_NORM);
     2543          }
     2544        }
     2545        weight[iSequence] = thisWeight;
     2546      }
     2547    } else if (switchType==4) {
     2548      // Devex
     2549      for (j=0;j<number;j++) {
     2550        int iSequence = index[j];
     2551        double thisWeight = weight[iSequence];
     2552        // row has -1
     2553        double pivot = updateBy[iSequence]*scaleFactor;
     2554        updateBy[iSequence]=0.0;
     2555        double value = pivot * pivot*devex_;
     2556        if (reference(iSequence))
     2557          value += 1.0;
     2558        weight[iSequence] = max(0.99*thisWeight,value);
     2559      }
    9592560    }
    9602561    // restore outgoing weight
     
    9632564    alternateWeights_->clear();
    9642565    spareColumn2->setNumElements(0);
     2566    //#define SOME_DEBUG_1
    9652567#ifdef SOME_DEBUG_1
    9662568    // check for accuracy
    9672569    int iCheck=229;
    968     printf("weight for iCheck is %g\n",weights_[iCheck]);
    969     if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
    970         !model_->getStatus(iCheck)!=isFixed)
    971       checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     2570    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     2571    int numberRows = model_->numberRows();
     2572    //int numberColumns = model_->numberColumns();
     2573    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
     2574      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
     2575          !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
     2576        checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
     2577    }
    9722578#endif
    9732579    updates->setNumElements(0);
     
    10082614  tolerance *= tolerance; // as we are using squares
    10092615
    1010   int numberWanted;
    1011   if (switchType<2 ) {
    1012     numberWanted = number+1;
    1013   } else if (switchType==2) {
    1014     numberWanted = max(2000,number/8);
    1015   } else {
    1016     int numberElements = model_->factorization()->numberElements();
    1017     double ratio = (double) numberElements/(double) model_->numberRows();
    1018     numberWanted = max(2000,number/8);
    1019     if (ratio<1.0) {
    1020       numberWanted = max(2000,number/20);
    1021     } else if (ratio>10.0) {
    1022       ratio = number * (ratio/80.0);
    1023       if (ratio>number)
    1024         numberWanted=number+1;
    1025       else
    1026         numberWanted = max(2000,(int) ratio);
    1027     }
    1028   }
    10292616  int iPass;
    10302617  // Setup two passes
     
    10392626  for (iPass=0;iPass<2;iPass++) {
    10402627    int end = start[2*iPass+1];
    1041     for (i=start[2*iPass];i<end;i++) {
    1042       iSequence = index[i];
    1043       double value = infeas[iSequence];
    1044       if (value>tolerance) {
    1045         double weight = weights_[iSequence];
    1046         //weight=1.0;
    1047         if (value>bestDj*weight) {
    1048           // check flagged variable and correct dj
    1049           if (!model_->flagged(iSequence)) {
    1050             bestDj=value/weight;
    1051             bestSequence = iSequence;
    1052           } else {
    1053             // just to make sure we don't exit before got something
    1054             numberWanted++;
     2628    if (switchType<5) {
     2629      for (i=start[2*iPass];i<end;i++) {
     2630        iSequence = index[i];
     2631        double value = infeas[iSequence];
     2632        if (value>tolerance) {
     2633          double weight = weights_[iSequence];
     2634          //weight=1.0;
     2635          if (value>bestDj*weight) {
     2636            // check flagged variable and correct dj
     2637            if (!model_->flagged(iSequence)) {
     2638              bestDj=value/weight;
     2639              bestSequence = iSequence;
     2640            } else {
     2641              // just to make sure we don't exit before got something
     2642              numberWanted++;
     2643            }
    10552644          }
    10562645        }
    1057       }
    1058       numberWanted--;
    1059       if (!numberWanted)
    1060         break;
     2646        numberWanted--;
     2647        if (!numberWanted)
     2648          break;
     2649      }
     2650    } else {
     2651      // Dantzig
     2652      for (i=start[2*iPass];i<end;i++) {
     2653        iSequence = index[i];
     2654        double value = infeas[iSequence];
     2655        if (value>tolerance) {
     2656          if (value>bestDj) {
     2657            // check flagged variable and correct dj
     2658            if (!model_->flagged(iSequence)) {
     2659              bestDj=value;
     2660              bestSequence = iSequence;
     2661            } else {
     2662              // just to make sure we don't exit before got something
     2663              numberWanted++;
     2664            }
     2665          }
     2666        }
     2667        numberWanted--;
     2668        if (!numberWanted)
     2669          break;
     2670      }
    10612671    }
    10622672    if (!numberWanted)
     
    10922702    if (mode==1&&!weights_)
    10932703      numberSwitched_=0; // Reset
    1094     if(!numberSwitched_)
    1095       return;
    10962704  }
    10972705  // alternateWeights_ is defined as indexed but is treated oddly
     
    13412949      alternateWeights_->setNumElements(number);
    13422950    } else {
    1343       for (i=0;i<number;i++) {
    1344         int iRow = which[i];
    1345         int iPivot = pivotVariable[iRow];
    1346         if (reference(iPivot)) {
    1347           devex_ += work[iRow]*work[iRow];
    1348           newWork[iRow] = -2.0*work[iRow];
    1349           newWhich[newNumber++]=iRow;
    1350         }
    1351       }
    1352       if (!newWork[pivotRow])
    1353         newWhich[newNumber++]=pivotRow; // add if not already in
    1354       newWork[pivotRow] = -2.0*max(devex_,0.0);
     2951      if (mode_!=4||numberSwitched_>1) {
     2952        for (i=0;i<number;i++) {
     2953          int iRow = which[i];
     2954          int iPivot = pivotVariable[iRow];
     2955          if (reference(iPivot)) {
     2956            devex_ += work[iRow]*work[iRow];
     2957            newWork[iRow] = -2.0*work[iRow];
     2958            newWhich[newNumber++]=iRow;
     2959          }
     2960        }
     2961        if (!newWork[pivotRow])
     2962          newWhich[newNumber++]=pivotRow; // add if not already in
     2963        newWork[pivotRow] = -2.0*max(devex_,0.0);
     2964      } else {
     2965        for (i=0;i<number;i++) {
     2966          int iRow = which[i];
     2967          int iPivot = pivotVariable[iRow];
     2968          if (reference(iPivot))
     2969            devex_ += work[iRow]*work[iRow];
     2970        }
     2971      }
    13552972      if (reference(sequenceIn)) {
    13562973        devex_+=1.0;
     
    13883005    printf("old weight %g, new %g\n",oldDevex,devex_);
    13893006#endif
    1390   double check = max(devex_,oldDevex);;
     3007  double check = max(devex_,oldDevex)+0.1;
    13913008  weights_[sequenceIn] = devex_;
    1392   if ( fabs ( devex_ - oldDevex ) > 1.0e-1 * check ) {
     3009  double testValue=0.1;
     3010  if (mode_==4&&numberSwitched_==1)
     3011    testValue = 0.5;
     3012  if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
    13933013#ifdef CLP_DEBUG
    13943014    if ((model_->messageHandler()->logLevel()&48)==16)
    13953015      printf("old weight %g, new %g\n",oldDevex,devex_);
    13963016#endif
    1397     if ( fabs ( devex_ - oldDevex ) > 1.0e5 * check ) {
     3017    //printf("old weight %g, new %g\n",oldDevex,devex_);
     3018    testValue=0.99;
     3019    if (mode_==1)
     3020      testValue=1.01e1; // make unlikely to do if steepest
     3021    else if (mode_==4&&numberSwitched_==1)
     3022      testValue=0.9;
     3023    double difference = fabs(devex_-oldDevex);
     3024    if ( difference > testValue * check ) {
    13983025      // need to redo
    13993026      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
Note: See TracChangeset for help on using the changeset viewer.