Changeset 212 for branches


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

lots of stuff

Location:
branches/pre
Files:
4 added
17 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpMessage.cpp

    r210 r212  
    5050  {CLP_PACKEDSCALE_FINAL,1003,2,"Final range of elements is %g to %g"},
    5151  {CLP_PACKEDSCALE_FORGET,1004,2,"Not bothering to scale as good enough"},
    52   {CLP_INITIALIZE_STEEP,1005,1,"Initializing steepest edge weights - old %g, new %g"},
     52  {CLP_INITIALIZE_STEEP,1005,2,"Initializing steepest edge weights - old %g, new %g"},
    5353  {CLP_UNABLE_OPEN,6001,0,"Unable to open file %s for reading"},
    5454  {CLP_BAD_BOUNDS,6002,0,"%d bad bound pairs or bad objectives were found - first at %c%d"},
     
    6464  {CLP_QUADRATIC_PRIMAL_DETAILS,109,32,"coeff %g, %g, %g - dj %g - deriv zero at %g, sj at %g"},
    6565  {CLP_IDIOT_ITERATION,30,1,"%d infeas %g, obj %g - mu %g, its %d, %d interior"},
    66   {CLP_INFEASIBLE,3003,1,"Analysis indicates model infeasible"},
    67   {CLP_MATRIX_CHANGE,31,1,"Matrix can not be converted into %s"},
    68   {CLP_TIMING,32,1,"%s objective %.9g - %d iterations time %.2f2%?, Presolve %.2f%?, Idiot %.2f%?"},
     66  {CLP_INFEASIBLE,3003,1,"Analysis indicates model infeasible or unbounded"},
     67  {CLP_MATRIX_CHANGE,31,2,"Matrix can not be converted into %s"},
     68  {CLP_TIMING,32,1,"%s objective %.10g - %d iterations time %.2f2%?, Presolve %.2f%?, Idiot %.2f%?"},
    6969  {CLP_INTERVAL_TIMING,33,2,"%s took %.2f seconds (total %.2f)"},
    7070  {CLP_SPRINT,34,1,"Pass %d took %d iterations, objective %g, dual infeasibilities %g( %d)"},
  • branches/pre/ClpPackedMatrix.cpp

    r208 r212  
    282282    else if (numberRows*2<numberColumns)
    283283      factor=0.2;
    284     if (model->numberIterations()%50==0)
    285       printf("%d nonzero\n",numberInRowArray);
     284    //if (model->numberIterations()%50==0)
     285    //printf("%d nonzero\n",numberInRowArray);
    286286  }
    287287  if (numberInRowArray>factor*numberRows||!rowCopy) {
     
    11261126  // mark empty and fixed columns
    11271127  // also see if worth scaling
    1128   assert (model->scalingFlag()==1); // dynamic not implemented
     1128  assert (model->scalingFlag()<4); // dynamic not implemented
    11291129  double largest=0.0;
    11301130  double smallest=1.0e50;
     
    11641164    return 1;
    11651165  } else {
     1166    int scalingMethod = model->scalingFlag();
     1167    if (scalingMethod==3) {
     1168      // Choose between 1 and 2
     1169      if (smallest<1.0e-5||smallest*largest<1.0)
     1170        scalingMethod=1;
     1171      else
     1172        scalingMethod=2;
     1173    }
    11661174    // and see if there any empty rows
    11671175    for (iRow=0;iRow<numberRows;iRow++) {
     
    11831191    double overallLargest=-1.0e-30;
    11841192    double overallSmallest=1.0e30;
    1185     if (model->scalingFlag()==1) {
     1193    if (scalingMethod==1) {
    11861194      // Maximum in each row
    11871195      for (iRow=0;iRow<numberRows;iRow++) {
     
    12021210      }
    12031211    } else {
    1204       assert(model->scalingFlag()==2);
     1212      assert(scalingMethod==2);
    12051213      int numberPass=3;
    12061214#ifdef USE_OBJECTIVE
     
    13141322      }
    13151323    }
     1324#define RANDOMIZE
     1325#ifdef RANDOMIZE
     1326    // randomize by up to 10%
     1327    for (iRow=0;iRow<numberRows;iRow++) {
     1328      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
     1329      rowScale[iRow] *= (1.0+0.1*value);
     1330    }
     1331#endif
    13161332    overallLargest=1.0;
    1317     if (overallSmallest<1.0e-3)
     1333    if (overallSmallest<1.0e-1)
    13181334      overallLargest = 1.0/sqrt(overallSmallest);
    13191335    overallLargest = min(1000.0,overallLargest);
     
    13341350        }
    13351351        columnScale[iColumn]=overallLargest/largest;
     1352#ifdef RANDOMIZE
     1353        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
     1354        columnScale[iColumn] *= (1.0+0.1*value);
     1355#endif
    13361356        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
    13371357      }
  • branches/pre/ClpPresolve.cpp

    r181 r212  
    894894void check_djs(CoinPostsolveMatrix *prob)
    895895{
    896   return;
     896  //return;
    897897  double *colels        = prob->colels_;
    898898  int *hrow             = prob->hrow_;
     
    935935        double coeff = colels[k];
    936936        k = link[k];
    937 
    938937        dj -= rowduals[row] * coeff;
    939938        rsol[row] += solutionValue*coeff;
    940939      }
    941       if (! (fabs(rcosts[colx] - dj) < 100*ZTOLDP))
     940      if (! (fabs(rcosts[colx] - dj) < 1.0e-4))
    942941        printf("BAD DJ:  %d %g %g\n",
    943942               colx, rcosts[colx], dj);
     
    965964        double dj = rowduals[rowx];
    966965        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
    967           if (dj <-1.0e-6)
     966          if (dj <-1.0e-5)
    968967            printf("neg rDJ:  %d %g  - %g %g %g\n",
    969968                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
    970969        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
    971           if (dj > 1.0e-6)
     970          if (dj > 1.0e-5)
    972971            printf("pos rDJ:  %d %g  - %g %g %g\n",
    973972                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
    974973        } else {
    975           if (fabs(dj) >1.0e-6)
     974          if (fabs(dj) >1.0e-5)
    976975            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
    977976                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
     
    10661065  return (! (i >= 0));
    10671066}
     1067#if     DEBUG_PRESOLVE
     1068static void matrix_bounds_ok(const double *lo, const double *up, int n)
     1069{
     1070  int i;
     1071  for (i=0; i<n; i++) {
     1072    PRESOLVEASSERT(lo[i] <= up[i]);
     1073    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
     1074    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
     1075  }
     1076}
     1077#endif
    10681078CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
    10691079                                     double maxmin_,
  • 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,
  • branches/pre/ClpSimplex.cpp

    r210 r212  
    9393  savedSolution_(NULL),
    9494  columnScale_(NULL),
    95   scalingFlag_(1),
     95  scalingFlag_(3),
    9696  numberTimesOptimal_(0),
    9797  changeMade_(1),
     
    196196  savedSolution_(NULL),
    197197  columnScale_(NULL),
    198   scalingFlag_(1),
     198  scalingFlag_(3),
    199199  numberTimesOptimal_(0),
    200200  changeMade_(1),
     
    692692
    693693  return status;
     694}
     695// Clean up status
     696void
     697ClpSimplex::cleanStatus()
     698{
     699  int iRow,iColumn;
     700  int numberBasic=0;
     701  // make row activities correct
     702  memset(rowActivityWork_,0,numberRows_*sizeof(double));
     703  times(1.0,columnActivityWork_,rowActivityWork_);
     704  if (!status_)
     705    createStatus();
     706  for (iRow=0;iRow<numberRows_;iRow++) {
     707    if (getRowStatus(iRow)==basic)
     708      numberBasic++;
     709    else {
     710      setRowStatus(iRow,superBasic);
     711      // but put to bound if close
     712      if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
     713          <=primalTolerance_) {
     714        rowActivityWork_[iRow]=rowLowerWork_[iRow];
     715        setRowStatus(iRow,atLowerBound);
     716      } else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
     717                 <=primalTolerance_) {
     718        rowActivityWork_[iRow]=rowUpperWork_[iRow];
     719        setRowStatus(iRow,atUpperBound);
     720      }
     721    }
     722  }
     723  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     724    if (getColumnStatus(iColumn)==basic) {
     725      if (numberBasic==numberRows_) {
     726        // take out of basis
     727        setColumnStatus(iColumn,superBasic);
     728        // but put to bound if close
     729        if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
     730            <=primalTolerance_) {
     731          columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     732          setColumnStatus(iColumn,atLowerBound);
     733        } else if (fabs(columnActivityWork_[iColumn]
     734                        -columnUpperWork_[iColumn])
     735                   <=primalTolerance_) {
     736          columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     737          setColumnStatus(iColumn,atUpperBound);
     738        }
     739      } else
     740        numberBasic++;
     741    } else {
     742      setColumnStatus(iColumn,superBasic);
     743      // but put to bound if close
     744      if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
     745          <=primalTolerance_) {
     746        columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     747        setColumnStatus(iColumn,atLowerBound);
     748      } else if (fabs(columnActivityWork_[iColumn]
     749                      -columnUpperWork_[iColumn])
     750                 <=primalTolerance_) {
     751        columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     752        setColumnStatus(iColumn,atUpperBound);
     753      }
     754    }
     755  }
    694756}
    695757
     
    9771039    } else {
    9781040      // values pass has less coding
    979       // make row activities correct
    980       memset(rowActivityWork_,0,numberRows_*sizeof(double));
    981       times(1.0,columnActivityWork_,rowActivityWork_);
     1041      // make row activities correct and clean basis a bit
     1042      cleanStatus();
    9821043      if (status_) {
    983         // set values from warm start (if sensible)
    9841044        int numberBasic=0;
    9851045        for (iRow=0;iRow<numberRows_;iRow++) {
    9861046          if (getRowStatus(iRow)==basic)
    9871047            numberBasic++;
    988           else {
    989             setRowStatus(iRow,superBasic);
    990             // but put to bound if close
    991             if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
    992                 <=primalTolerance_) {
    993               rowActivityWork_[iRow]=rowLowerWork_[iRow];
    994               setRowStatus(iRow,atLowerBound);
    995             } else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
    996                        <=primalTolerance_) {
    997               rowActivityWork_[iRow]=rowUpperWork_[iRow];
    998               setRowStatus(iRow,atUpperBound);
    999             }
    1000           }
    1001          
    10021048        }
    10031049        totalSlacks=numberBasic;
    10041050        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1005           if (getColumnStatus(iColumn)==basic) {
    1006             if (numberBasic==numberRows_) {
    1007               // take out of basis
    1008               setColumnStatus(iColumn,superBasic);
    1009               // but put to bound if close
    1010               if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    1011                   <=primalTolerance_) {
    1012                 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    1013                 setColumnStatus(iColumn,atLowerBound);
    1014               } else if (fabs(columnActivityWork_[iColumn]
    1015                             -columnUpperWork_[iColumn])
    1016                          <=primalTolerance_) {
    1017                 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    1018                 setColumnStatus(iColumn,atUpperBound);
    1019               }
    1020             } else
    1021               numberBasic++;
    1022           } else {
    1023             setColumnStatus(iColumn,superBasic);
    1024             // but put to bound if close
    1025             if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    1026                 <=primalTolerance_) {
    1027               columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    1028               setColumnStatus(iColumn,atLowerBound);
    1029             } else if (fabs(columnActivityWork_[iColumn]
    1030                           -columnUpperWork_[iColumn])
    1031                        <=primalTolerance_) {
    1032               columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    1033               setColumnStatus(iColumn,atUpperBound);
    1034             }
    1035           }
     1051          if (getColumnStatus(iColumn)==basic)
     1052            numberBasic++;
    10361053        }
    10371054      } else {
     
    12831300  savedSolution_(NULL),
    12841301  columnScale_(NULL),
    1285   scalingFlag_(1),
     1302  scalingFlag_(3),
    12861303  numberTimesOptimal_(0),
    12871304  changeMade_(1),
     
    13801397  savedSolution_(NULL),
    13811398  columnScale_(NULL),
    1382   scalingFlag_(1),
     1399  scalingFlag_(3),
    13831400  numberTimesOptimal_(0),
    13841401  changeMade_(1),
     
    22522269ClpSimplex::scaling(int mode)
    22532270{
    2254   if (mode>0&&mode<3) {
     2271  if (mode>0&&mode<4) {
    22552272    scalingFlag_=mode;
    22562273  } else if (!mode) {
     
    25942611      <<CoinMessageEol;
    25952612    // Set bounds slightly loose
     2613    double useTolerance = 1.0e-3;
    25962614    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2597       if (saveUpper[iColumn]>saveLower[iColumn]+tolerance) {
    2598         if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance+1.0e8) {
     2615      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
     2616        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e8) {
    25992617          // relax enough so will have correct dj
    26002618#if 1
    26012619          columnLower_[iColumn]=max(saveLower[iColumn],
    2602                                     columnLower_[iColumn]-100.0*tolerance);
     2620                                    columnLower_[iColumn]-100.0*useTolerance);
    26032621          columnUpper_[iColumn]=min(saveUpper[iColumn],
    2604                                     columnUpper_[iColumn]+100.0*tolerance);
     2622                                    columnUpper_[iColumn]+100.0*useTolerance);
    26052623#else
    26062624          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
    2607             if (columnUpper_[iColumn]- 100.0*tolerance>saveLower[iColumn]) {
    2608               columnLower_[iColumn]=columnUpper_[iColumn]-100.0*tolerance;
     2625            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
     2626              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
    26092627            } else {
    26102628              columnLower_[iColumn]=saveLower[iColumn];
    26112629              columnUpper_[iColumn]=min(saveUpper[iColumn],
    2612                                         saveLower[iColumn]+100.0*tolerance);
     2630                                        saveLower[iColumn]+100.0*useTolerance);
    26132631            }
    26142632          } else {
    2615             if (columnLower_[iColumn]+100.0*tolerance<saveUpper[iColumn]) {
    2616               columnUpper_[iColumn]=columnLower_[iColumn]+100.0*tolerance;
     2633            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
     2634              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
    26172635            } else {
    26182636              columnUpper_[iColumn]=saveUpper[iColumn];
    26192637              columnLower_[iColumn]=max(saveLower[iColumn],
    2620                                         saveUpper[iColumn]-100.0*tolerance);
     2638                                        saveUpper[iColumn]-100.0*useTolerance);
    26212639            }
    26222640          }
     
    26252643          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
    26262644            // relax a bit
    2627             columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
     2645            columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*useTolerance,
    26282646                                        saveUpper[iColumn]);
    26292647          }
    26302648          if (columnLower_[iColumn]>saveLower[iColumn]) {
    26312649            // relax a bit
    2632             columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
     2650            columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*useTolerance,
    26332651                                        saveLower[iColumn]);
    26342652          }
     
    26582676    int savePerturbation = perturbation_;
    26592677    perturbation_=100;
     2678    bool denseFactorization = initialDenseFactorization();
     2679    // It will be safe to allow dense
     2680    setInitialDenseFactorization(true);
    26602681    returnCode = ((ClpSimplexPrimal *) this)->primal(1);
     2682    setInitialDenseFactorization(denseFactorization);
    26612683    perturbation_=savePerturbation;
    26622684  }
     
    26712693    int savePerturbation = perturbation_;
    26722694    perturbation_=100;
     2695    bool denseFactorization = initialDenseFactorization();
     2696    // It will be safe to allow dense
     2697    setInitialDenseFactorization(true);
    26732698    returnCode = ((ClpSimplexDual *) this)->dual(0);
     2699    setInitialDenseFactorization(denseFactorization);
    26742700    perturbation_=savePerturbation;
    26752701  }
     
    39053931  }
    39063932}
    3907 int
    3908 ClpSimplex::nextSuperBasic()
    3909 {
    3910   if (firstFree_>=0) {
    3911     int returnValue=firstFree_;
    3912     int iColumn=firstFree_+1;
    3913     if (algorithm_>0) {
    3914       // primal
    3915       for (;iColumn<numberRows_+numberColumns_;iColumn++) {
    3916         if (getStatus(iColumn)==superBasic) {
    3917           if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
    3918             solution_[iColumn]=lower_[iColumn];
    3919             setStatus(iColumn,atLowerBound);
    3920           } else if (fabs(solution_[iColumn]-upper_[iColumn])
    3921                      <=primalTolerance_) {
    3922             solution_[iColumn]=upper_[iColumn];
    3923             setStatus(iColumn,atUpperBound);
    3924           } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
    3925             setStatus(iColumn,isFree);
    3926             break;
    3927           } else {
    3928             break;
    3929           }
    3930         }
    3931       }
    3932     } else {
    3933       // dual
    3934       for (;iColumn<numberRows_+numberColumns_;iColumn++) {
    3935         if (getStatus(iColumn)==isFree)
    3936           if (fabs(dj_[iColumn])>1.0e2*dualTolerance_)
    3937             break;
    3938       }
    3939     }
    3940     firstFree_ = iColumn;
    3941     if (firstFree_==numberRows_+numberColumns_)
    3942       firstFree_=-1;
    3943     return returnValue;
    3944   } else {
    3945     return -1;
    3946   }
    3947 }
    39483933/* Pivot in a variable and out a variable.  Returns 0 if okay,
    39493934   1 if inaccuracy forced re-factorization, -1 if would be singular.
  • branches/pre/ClpSimplexDual.cpp

    r210 r212  
    202202  ClpDataSave data = saveData();
    203203
     204  // Save so can see if doing after primal
     205  int initialStatus=problemStatus_;
     206
    204207  // If values pass then save given duals round check solution
    205208  double * saveDuals = NULL;
     
    208211    memcpy(saveDuals,dual_,numberRows_*sizeof(double));
    209212  }
    210 
     213 
    211214  // sanity check
    212215  // initialize - no values pass and algorithm_ is -1
     
    221224    // If values pass then scale pi
    222225    if (ifValuesPass) {
    223       if (!problemStatus_&&perturbation_<100)
     226      if (problemStatus_&&perturbation_<100)
    224227        perturb();
    225228      int i;
     
    311314      // If getting nowhere - why not give it a kick
    312315      // does not seem to work too well - do some more work
    313       if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)) {
     316      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
     317          &&initialStatus!=10) {
    314318        perturb();
    315319        // Can't get here if values pass
     
    29112915  if (perturbation_>100)
    29122916    return; //perturbed already
     2917  if (perturbation_==100)
     2918    perturbation_=50; // treat as normal
    29132919  bool modifyRowCosts=false;
    29142920  // dual perturbation
     
    29202926  int maxLength=0;
    29212927  int minLength=numberRows_;
    2922   if (perturbation_==50) {
    2923     int count[101];
    2924     int numberTotal=0;
    2925     memset(count,0,101*sizeof(int));
    2926     int iColumn;
    2927     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2928       if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
    2929         double value = fabs(objectiveWork_[iColumn]);
    2930         int k;
    2931         numberTotal++;
    2932         if (!value) {
    2933           k=0;
    2934         } else if (value>99.5) {
    2935           k=100;
    2936         } else {
    2937           double lo = floor(value+0.5);
    2938           if (fabs(lo-value)<1.0e-8)
    2939             k = (int) lo;
    2940           else
    2941             k=100;
    2942         }
    2943         count[k]++;
    2944         int length = lengths[iColumn];
    2945         if (length>2) {
    2946           maxLength = max(maxLength,length);
    2947           minLength = min(minLength,length);
    2948         }
    2949       }
    2950     }
    2951     int numberZero = count[0];
    2952     int numberNonZero = numberTotal-numberZero;
    2953     int numberOdd=count[100];
    2954     int numberInteger = numberNonZero-numberOdd;
    2955     if (3*numberOdd>numberTotal) {
    2956       // be more subtle
    2957       maximumFraction=1.0e-6;
    2958       constantPerturbation *= 0.1;
    2959     } else if (count[1]==numberInteger&&!numberOdd) {
    2960       // be less subtle
    2961       maximumFraction=2.0e-5;
    2962     }
    2963   } else {
    2964     int iColumn;
    2965     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2966       if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
    2967         int length = lengths[iColumn];
    2968         if (length>2) {
    2969           maxLength = max(maxLength,length);
    2970           minLength = min(minLength,length);
    2971         }
     2928  if (!numberIterations_&&perturbation_==50) {
     2929    // See if we need to perturb
     2930    double * sort = new double[numberColumns_];
     2931    int i;
     2932    for (i=0;i<numberColumns_;i++) {
     2933      double value = fabs(objectiveWork_[i]);
     2934      sort[i]=value;
     2935    }
     2936    std::sort(sort,sort+numberColumns_);
     2937    int number=1;
     2938    double last = sort[0];
     2939    for (i=1;i<numberColumns_;i++) {
     2940      if (last!=sort[i])
     2941        number++;
     2942      last=sort[i];
     2943    }
     2944    delete [] sort;
     2945    if (number*4>numberColumns_)
     2946      return; // good enough
     2947  }
     2948  int iColumn;
     2949  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     2950    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
     2951      int length = lengths[iColumn];
     2952      if (length>2) {
     2953        maxLength = max(maxLength,length);
     2954        minLength = min(minLength,length);
    29722955      }
    29732956    }
     
    29862969    maximumFraction = m[min(perturbation_-51,10)];
    29872970  }
    2988   int iRow,iColumn;
     2971  int iRow;
    29892972  double smallestNonZero=1.0e100;
    29902973  if (perturbation_>=50) {
    29912974    perturbation = 1.0e-8;
     2975    bool allSame=true;
     2976    double lastValue=0.0;
    29922977    for (iRow=0;iRow<numberRows_;iRow++) {
    2993       if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
     2978      double lo = rowLowerWork_[iRow];
     2979      double up = rowUpperWork_[iRow];
     2980      if (lo<up) {
    29942981        double value = fabs(rowObjectiveWork_[iRow]);
    29952982        perturbation = max(perturbation,value);
     
    29982985          smallestNonZero = min(smallestNonZero,value);
    29992986        }
     2987      }
     2988      if (lo&&lo>-1.0e10) {
     2989        lo=fabs(lo);
     2990        if (!lastValue)
     2991          lastValue=lo;
     2992        else if (fabs(lo-lastValue)>1.0e-7)
     2993          allSame=false;
     2994      }
     2995      if (up&&up<1.0e10) {
     2996        up=fabs(up);
     2997        if (!lastValue)
     2998          lastValue=up;
     2999        else if (fabs(up-lastValue)>1.0e-7)
     3000          allSame=false;
    30003001      }
    30013002    }
    30023003    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3003       if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
     3004      double lo = columnLowerWork_[iColumn];
     3005      double up = columnUpperWork_[iColumn];
     3006      if (lo<up) {
    30043007        double value =
    30053008          fabs(objectiveWork_[iColumn]);
     
    30093012        }
    30103013      }
     3014      if (lo&&lo>-1.0e10) {
     3015        lo=fabs(lo);
     3016        if (!lastValue)
     3017          lastValue=lo;
     3018        else if (fabs(lo-lastValue)>1.0e-7)
     3019          allSame=false;
     3020      }
     3021      if (up&&up<1.0e10) {
     3022        up=fabs(up);
     3023        if (!lastValue)
     3024          lastValue=up;
     3025        else if (fabs(up-lastValue)>1.0e-7)
     3026          allSame=false;
     3027      }
     3028    }
     3029    if (allSame) {
     3030      // Really hit perturbation
     3031      maximumFraction = max(1.0e-2*lastValue,maximumFraction);
    30113032    }
    30123033    perturbation = min(perturbation,smallestNonZero/maximumFraction);
     
    30393060  // modify costs
    30403061  bool printOut=(handler_->logLevel()==63);
     3062  printOut=false;
    30413063  if (modifyRowCosts) {
    30423064    for (iRow=0;iRow<numberRows_;iRow++) {
     
    39854007  delete [] scaled;
    39864008}
     4009int
     4010ClpSimplexDual::nextSuperBasic()
     4011{
     4012  if (firstFree_>=0) {
     4013    int returnValue=firstFree_;
     4014    int iColumn=firstFree_+1;
     4015    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
     4016      if (getStatus(iColumn)==isFree)
     4017        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_)
     4018          break;
     4019    }
     4020    firstFree_ = iColumn;
     4021    if (firstFree_==numberRows_+numberColumns_)
     4022      firstFree_=-1;
     4023    return returnValue;
     4024  } else {
     4025    return -1;
     4026  }
     4027}
  • branches/pre/ClpSimplexPrimal.cpp

    r210 r212  
    185185  // save data
    186186  ClpDataSave data = saveData();
    187 
    188 
    189     // initialize - maybe values pass and algorithm_ is +1
     187 
     188  // Save so can see if doing after dual
     189  int initialStatus=problemStatus_;
     190  // initialize - maybe values pass and algorithm_ is +1
    190191    if (!startup(ifValuesPass)) {
    191192
     
    242243      // If getting nowhere - why not give it a kick
    243244#if 1
    244       if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
     245      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
     246          &&initialStatus!=10)
    245247        perturb(1);
    246248#endif
     
    262264     
    263265      // Iterate
    264       whileIterating();
     266      whileIterating(ifValuesPass);
    265267    }
    266268  }
     
    292294*/
    293295int
    294 ClpSimplexPrimal::whileIterating()
     296ClpSimplexPrimal::whileIterating(int valuesOption)
    295297{
    296298  // Say if values pass
    297299  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
    298300  int returnCode=-1;
     301  int superBasicType=1;
     302  if (valuesOption>1)
     303    superBasicType=3;
    299304  // status stays at -1 while iterating, >=0 finished, -2 to invert
    300305  // status -3 to go to top without an invert
     
    352357    } else {
    353358      // in values pass
    354       int sequenceIn=nextSuperBasic();
     359      int sequenceIn=nextSuperBasic(superBasicType,columnArray_[0]);
     360      if (valuesOption>1)
     361        superBasicType=2;
    355362      if (sequenceIn<0) {
    356363        // end of values pass - initialize weights etc
     
    418425    }
    419426  }
     427  if (valuesOption>1)
     428    columnArray_[0]->setNumElements(0);
    420429  return returnCode;
    421430}
     
    604613          gutsOfSolution(NULL,NULL);
    605614          problemStatus_=-1; //continue
     615          goToDual=false;
    606616        } else {
    607617          // say infeasible
     
    13691379  if (perturbation_>100)
    13701380    return; //perturbed already
     1381  if (perturbation_==100)
     1382    perturbation_=50; // treat as normal
    13711383  int i;
     1384  if (!numberIterations_)
     1385    cleanStatus(); // make sure status okay
     1386  if (!numberIterations_&&perturbation_==50) {
     1387    // See if we need to perturb
     1388    double * sort = new double[numberRows_];
     1389    for (i=0;i<numberRows_;i++) {
     1390      double lo = fabs(lower_[i]);
     1391      double up = fabs(upper_[i]);
     1392      double value=0.0;
     1393      if (lo&&lo<1.0e20) {
     1394        if (up&&up<1.0e20)
     1395          value = 0.5*(lo+up);
     1396        else
     1397          value=lo;
     1398      } else {
     1399        if (up&&up<1.0e20)
     1400          value = up;
     1401      }
     1402      sort[i]=value;
     1403    }
     1404    std::sort(sort,sort+numberRows_);
     1405    int number=1;
     1406    double last = sort[0];
     1407    for (i=1;i<numberRows_;i++) {
     1408      if (last!=sort[i])
     1409        number++;
     1410      last=sort[i];
     1411    }
     1412    delete [] sort;
     1413    if (number*4>numberRows_)
     1414      return; // good enough
     1415  }
    13721416  // primal perturbation
    13731417  double perturbation=1.0e-20;
     1418  int numberNonZero=0;
    13741419  // maximum fraction of rhs/bounds to perturb
    1375   double maximumFraction = 1.0e-8;
     1420  double maximumFraction = 1.0e-7;
    13761421  if (perturbation_>=50) {
    13771422    perturbation = 1.0e-4;
     
    13871432        else
    13881433          upperValue=0.0;
    1389         double value = max(lowerValue,upperValue);
     1434        double value = max(fabs(lowerValue),fabs(upperValue));
    13901435        value = min(value,upper_[i]-lower_[i]);
     1436#if 1
     1437        if (value) {
     1438          perturbation += value;
     1439          numberNonZero++;
     1440        }
     1441#else
    13911442        perturbation = max(perturbation,value);
    1392       }
    1393     }
     1443#endif
     1444      }
     1445    }
     1446    if (numberNonZero)
     1447      perturbation /= (double) numberNonZero;
     1448    else
     1449      perturbation = 1.0e-1;
    13941450  } else if (perturbation_<100) {
    13951451    perturbation = pow(10.0,perturbation_);
     
    14011457  double largestPerCent=0.0;
    14021458  bool printOut=(handler_->logLevel()==63);
     1459  printOut=false; //off
    14031460  // Check if all slack
    14041461  int number=0;
     
    14801537                 i,lower_[i],lowerValue,upper_[i],upperValue);
    14811538      }
    1482       if (solution_[i]==lower_[i])
    1483         solution_[i]=lowerValue;
    1484       else if (solution_[i]==upper_[i])
    1485         solution_[i]=upperValue;
    14861539      lower_[i]=lowerValue;
    14871540      upper_[i]=upperValue;
     
    15181571        printf("row %d lower from %g to %g, upper from %g to %g\n",
    15191572               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
    1520       if (solution_[i]==lower_[i])
    1521         solution_[i]=lowerValue;
    1522       else if (solution_[i]==upper_[i])
    1523         solution_[i]=upperValue;
    15241573      lower_[i]=lowerValue;
    15251574      upper_[i]=upperValue;
     1575    }
     1576  }
     1577  // Clean up
     1578  for (i=0;i<numberColumns_+numberRows_;i++) {
     1579    switch(getStatus(i)) {
     1580     
     1581    case basic:
     1582      break;
     1583    case atUpperBound:
     1584      solution_[i]=upper_[i];
     1585      break;
     1586    case isFixed:
     1587    case atLowerBound:
     1588      solution_[i]=lower_[i];
     1589      break;
     1590    case isFree:
     1591      break;
     1592    case superBasic:
     1593      break;
    15261594    }
    15271595  }
     
    18101878      }
    18111879      // here do part of steepest - ready for next iteration
    1812       primalColumnPivot_->updateWeights(rowArray_[1]);
     1880      if (!ifValuesPass)
     1881        primalColumnPivot_->updateWeights(rowArray_[1]);
    18131882    } else {
    18141883      if (pivotRow_==-1) {
     
    19552024  }
    19562025}
     2026/* Get next superbasic -1 if none,
     2027   Normal type is 1
     2028   If type is 3 then initializes sorted list
     2029   if 2 uses list.
     2030*/
     2031int
     2032ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
     2033{
     2034  if (firstFree_>=0&&superBasicType) {
     2035    int returnValue=firstFree_;
     2036    int iColumn=firstFree_+1;
     2037    if (superBasicType>1) {
     2038      if (superBasicType>2) {
     2039        // Initialize list
     2040        // Wild guess that lower bound more natural than upper
     2041        int number=0;
     2042        double * work=columnArray->denseVector();
     2043        int * which=columnArray->getIndices();
     2044        for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     2045          if (getStatus(iColumn)==superBasic) {
     2046            if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
     2047              solution_[iColumn]=lower_[iColumn];
     2048              setStatus(iColumn,atLowerBound);
     2049            } else if (fabs(solution_[iColumn]-upper_[iColumn])
     2050                       <=primalTolerance_) {
     2051              solution_[iColumn]=upper_[iColumn];
     2052              setStatus(iColumn,atUpperBound);
     2053            } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
     2054              setStatus(iColumn,isFree);
     2055              break;
     2056            } else if (!flagged(iColumn)) {
     2057              // put ones near bounds at end after sorting
     2058              work[number]= - min(0.1*(solution_[iColumn]-lower_[iColumn]),
     2059                                  upper_[iColumn]-solution_[iColumn]);
     2060              which[number++] = iColumn;
     2061            }
     2062          }
     2063        }
     2064        CoinSort_2(work,work+number,which);
     2065        columnArray->setNumElements(number);
     2066        memset(work,0,number*sizeof(double));
     2067      }
     2068      int * which=columnArray->getIndices();
     2069      int number = columnArray->getNumElements();
     2070      if (!number) {
     2071        // finished
     2072        iColumn = numberRows_+numberColumns_;
     2073        returnValue=-1;
     2074      } else {
     2075        number--;
     2076        returnValue=which[number];
     2077        iColumn=returnValue;
     2078        columnArray->setNumElements(number);
     2079      }     
     2080    } else {
     2081      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
     2082        if (getStatus(iColumn)==superBasic) {
     2083          if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
     2084            solution_[iColumn]=lower_[iColumn];
     2085            setStatus(iColumn,atLowerBound);
     2086          } else if (fabs(solution_[iColumn]-upper_[iColumn])
     2087                     <=primalTolerance_) {
     2088            solution_[iColumn]=upper_[iColumn];
     2089            setStatus(iColumn,atUpperBound);
     2090          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
     2091            setStatus(iColumn,isFree);
     2092            break;
     2093          } else {
     2094            break;
     2095          }
     2096        }
     2097      }
     2098    }
     2099    firstFree_ = iColumn;
     2100    if (firstFree_==numberRows_+numberColumns_)
     2101      firstFree_=-1;
     2102    return returnValue;
     2103  } else {
     2104    return -1;
     2105  }
     2106}
     2107 
    19572108 
    19582109
  • branches/pre/ClpSolve.cpp

    r210 r212  
    103103    }
    104104    // We may be better off using original
    105     if (numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_) {
     105    if (presolve!=ClpSolve::presolveOff&&
     106        numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_) {
    106107      delete model2;
    107108      model2 = this;
     
    215216  if (model2->factorizationFrequency()==200) {
    216217    // User did not touch preset
    217     model2->setFactorizationFrequency(100+model2->numberRows()/100);
     218    model2->setFactorizationFrequency(100+model2->numberRows()/200);
    218219  }
    219220  if (method==ClpSolve::usePrimalorSprint) {
     
    241242      int nPasses=0;
    242243      Idiot info(*model2);
    243       if (numberRows>1000&&numberColumns>2*numberRows) {
     244      // Get ratio of average number of elements per column to 5
     245      double ratio  = ((double) numberElements/(double) numberColumns)/5.0;
     246      if (numberRows>200&&numberColumns>5000&&numberColumns>2*numberRows) {
    244247        if (plusMinus) {
    245248          // look at rhs
     
    265268            nPasses = max(nPasses,15);
    266269            if (numberElements<3*numberColumns)
    267               nPasses=0; // probably not worh it
     270              nPasses=10; // probably not worh it
     271            info.setLightweight(1); // say lightweight idiot
    268272          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
    269273            nPasses = 10+numberColumns/1000;
     
    284288          else
    285289            nPasses=5;
     290          info.setLightweight(1); // say lightweight idiot
    286291          if (numberElements<3*numberColumns)
    287             nPasses=0; // probably not worh it
     292            nPasses=2; // probably not worh it
    288293          else
    289294            nPasses = max(nPasses,5);
     295          nPasses = (int) (((double) nPasses)/ratio); // reduce if lots of elements per column
     296          if (nPasses<2)
     297            nPasses=0;
    290298          //info.setStartingWeight(1.0e-1);
    291299        }
     
    312320    if (doCrash)
    313321      model2->crash(1000,1);
    314     model2->primal(1);
     322    model2->primal(2);
    315323    time2 = CoinCpuTime();
    316324    timeCore = time2-timeX;
     
    479487    double originalOffset;
    480488    model2->getDblParam(ClpObjOffset,originalOffset);
     489    int totalIterations=0;
    481490    for (iPass=0;iPass<maxPass;iPass++) {
    482491      //printf("Bug until submodel new version\n");
     
    517526        currentModel = &small;
    518527      small.primal();
     528      totalIterations += small.numberIterations();
    519529      // move solution back
    520530      const double * solution = small.primalColumnSolution();
     
    594604      <<CoinMessageEol;
    595605    timeX=time2;
     606    model2->setNumberIterations(model2->numberIterations()+totalIterations);
    596607  } else {
    597608    assert (method!=ClpSolve::automatic); // later
  • branches/pre/Idiot.cpp

    r210 r212  
    122122  }
    123123  sum /= (double) (nnzero+1);
    124   // If mu not changed then compute
    125   if (mu_==1e-4)
    126     mu_= max(1.0e-3,sum*1.0e-5);
    127   maxIts_=20;
    128   maxIts2_=105;
     124  maxIts_=2;
    129125  if (numberPass<=0)
    130126    // Cast to double to avoid VACPP complaining
     
    132128  else
    133129    majorIterations_=numberPass;
     130  // If mu not changed then compute
     131  if (mu_==1e-4)
     132    mu_= max(1.0e-3,sum*1.0e-5);
     133  if (!lightWeight_) {
     134    maxIts2_=105;
     135  } else {
     136    mu_ *= 10.0;
     137    maxIts2_=23;
     138  }
    134139  //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
    135140  solve2(handler,messages);
  • branches/pre/Makefile.Clp

    r210 r212  
    77# between then specify the exact level you want, e.g., -O1 or -O2
    88OptLevel := -g
    9 OptLevel := -O1
     9#OptLevel := -O2
    1010
    1111LIBNAME := Clp
     
    5353ifeq ($(OptLevel),-g)
    5454#     CXXFLAGS += -DCLP_DEBUG
     55#CXXFLAGS += -DDEBUG_PRESOLVE
    5556endif
    5657ifeq ($(OptLevel),-O2)
  • branches/pre/Test/ClpMain.cpp

    r210 r212  
    2727#include "ClpDualRowSteepest.hpp"
    2828#include "ClpDualRowDantzig.hpp"
     29#include "ClpLinearObjective.hpp"
    2930#include "ClpPrimalColumnSteepest.hpp"
    3031#include "ClpPrimalColumnDantzig.hpp"
     
    6364  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
    6465  MAXIMIZE,MINIMIZE,EXIT,STDIN,UNITTEST,NETLIB_DUAL,NETLIB_PRIMAL,SOLUTION,
    65   TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,
     66  TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,
    6667
    6768  INVALID=1000
     
    857858    int allowImportErrors=0;
    858859    int keepImportNames=1;
    859     int doIdiot=-2;
     860    int doIdiot=-1;
    860861    int doCrash=0;
    861862    int doSprint=-1;
     
    880881    parameters[numberParameters++]=
    881882      ClpItem("allS!lack","Set basis back to all slack",
    882               ALLSLACK);
    883     parameters[numberParameters-1].setLonghelp
    884       (
    885        ""
     883              ALLSLACK,false);
     884    parameters[numberParameters-1].setLonghelp
     885      (
     886       "Useful for playing around"
    886887       );
    887888    parameters[numberParameters++]=
     
    916917    parameters[numberParameters-1].setLonghelp
    917918      (
    918        "This sets the directory which import, export, saveModel and restoreModel will use"
     919       "This sets the directory which import, export, saveModel and restoreModel will use.\
     920  It is initialized to './'"
    919921       );
    920922    parameters[numberParameters-1].setStringValue(directory);
     
    925927    parameters[numberParameters-1].setLonghelp
    926928      (
    927        ""
     929       "The dual algorithm in Clp is a single phase algorithm as opposed to a two phase\
     930 algorithm where you first get feasible then optimal.  If a problem has both upper and\
     931 lower bounds then it is trivial to get dual feasible by setting non basic variables\
     932 to correct bound.  If the gap between the upper and lower bounds of a variable is more\
     933 than the value of dualBound Clp introduces fake bounds so that it can make the problem\
     934 dual feasible.  This has the same effect as a composite objective function in the\
     935 primal algorithm.  Too high a value may mean more iterations, while too low a bound means\
     936 the code may go all the way and then have to increase the bounds.  OSL had a heuristic to\
     937 adjust bounds, maybe we need that here."
    928938       );
    929939    parameters[numberParameters-1].setDoubleValue(models->dualBound());
     
    936946    parameters[numberParameters-1].setLonghelp
    937947      (
    938        ""
     948       "Clp can use any pivot selection algorithm which the user codes as long as it\
     949 implements the features in the abstract pivot base class.  The Dantzig method is implemented\
     950 to show a simple method but its use is deprecated.  Steepest is the method of choice and there\
     951 are two variants which keep all weights updated but only scan a subset each iteration.\
     952 Partial switches this on while automatic decides at each iteration based on information\
     953 about the factorization."
    939954       );
    940955    parameters[numberParameters++]=
     
    943958    parameters[numberParameters-1].setLonghelp
    944959      (
    945        ""
     960       "This command solves the current model using the dual steepest algorithm.\
     961The time and iterations may be affected by settings such as presolve, scaling, crash\
     962 and also by dual pivot method, fake bound on variables and dual and primal tolerances."
    946963       );
    947964    parameters[numberParameters++]=
     
    951968    parameters[numberParameters-1].setLonghelp
    952969      (
    953        ""
     970       "Normally the default tolerance is fine, but you may want to increase it a\
     971 bit if a dual run seems to be having a hard time"
    954972       );
    955973    parameters[numberParameters-1].setDoubleValue(models->dualTolerance());
     
    959977    parameters[numberParameters-1].setLonghelp
    960978      (
    961        ""
     979       "This stops the execution of Clp, end, exit, quit and stop are synonyms"
    962980       );
    963981    parameters[numberParameters++]=
     
    966984    parameters[numberParameters-1].setLonghelp
    967985      (
    968        ""
     986       "The default is not to use any model which had errors when reading the mps file.\
     987  Setting this to 'on' will allow all errors from which the code can recover\
     988 by ignoring the error."
    969989       );
    970990    parameters[numberParameters++]=
     
    973993    parameters[numberParameters-1].setLonghelp
    974994      (
    975        ""
     995       "This stops the execution of Clp, end, exit, quit and stop are synonyms"
    976996       );
    977997    parameters[numberParameters++]=
     
    9801000    parameters[numberParameters-1].setLonghelp
    9811001      (
    982        ""
     1002       "This will write an MPS format file to the given file name.  It will use the default\
     1003 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
     1004 is initialized to 'default.mps'."
    9831005       );
    9841006    parameters[numberParameters-1].setStringValue(exportFile);
     
    9911013    parameters[numberParameters++]=
    9921014      ClpItem("idiot!Crash","Whether to try idiot crash",
    993               -2,200,IDIOT);
    994     parameters[numberParameters-1].setLonghelp
    995       (
    996        ""
     1015              -1,200,IDIOT);
     1016    parameters[numberParameters-1].setLonghelp
     1017      (
     1018       "This is a type of 'crash' which works well on some homogeneous problems.\
     1019 It works best on problems with unit elements and rhs but will do something to any model.  It should only be\
     1020 used before primal.  It can be set to -1 when the code decides for itself whether to use it,\
     1021 0 to switch off or n > 0 to do n passes."
    9971022       );
    9981023    parameters[numberParameters-1].setIntValue(doIdiot);
     
    10021027    parameters[numberParameters-1].setLonghelp
    10031028      (
    1004        ""
     1029       "This will read an MPS format file from the given file name.  It will use the default\
     1030 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
     1031 is initialized to '', i.e. it must be set.  If you have libgz then it can read compressed\
     1032 files 'xxxxxxxx.gz'.."
    10051033       );
    10061034    parameters[numberParameters-1].setStringValue(importFile);
     
    10111039    parameters[numberParameters-1].setLonghelp
    10121040      (
    1013        ""
     1041       "It saves space to get rid of names so if you need to you can set this to off."
    10141042       );
    10151043    parameters[numberParameters++]=
     
    10181046    parameters[numberParameters-1].setLonghelp
    10191047      (
    1020        ""
     1048       "If 0 then there should be no output in normal circumstances.  1 is probably the best\
     1049 value for most uses, while 2 and 3 give more information."
    10211050       );
    10221051    parameters[numberParameters-1].setIntValue(models->logLevel());
     
    10261055    parameters[numberParameters-1].setLonghelp
    10271056      (
    1028        ""
     1057       "The default is minimize - use 'maximize' for maximization.\n\
     1058You can also use the parameters 'direction maximize'."
    10291059       );
    10301060    parameters[numberParameters++]=
     
    10511081    parameters[numberParameters-1].setLonghelp
    10521082      (
    1053        ""
     1083       "The default is minimize - use 'maximize' for maximization.\n\
     1084This should only be necessary if you have previously set maximization \
     1085You can also use the parameters 'direction minimize'."
    10541086       );
    10551087    parameters[numberParameters++]=
     
    10581090    parameters[numberParameters-1].append("on");
    10591091    parameters[numberParameters-1].setLonghelp
    1060       (""
     1092      ("The default for the Clp library is to put out messages such as:\n\
     1093   Clp0005 2261  Objective 109.024 Primal infeas 944413 (758)\n\
     1094but this program turns this off to make it look more friendly.  It can be useful\
     1095 to turn them back on if you want to be able 'grep' for particular messages or if\
     1096 you intend to override the behavior of a particular message."
    10611097       );
    10621098    parameters[numberParameters++]=
     
    10711107    parameters[numberParameters-1].setLonghelp
    10721108      (
    1073        ""
     1109       "Clp will go faster if the matrix can be converted to a network.  The matrix\
     1110 operations may be a bit faster with more efficient storage, but the main advantage\
     1111 comes from using a network factorization.  It will probably not be as fast as a \
     1112specialized network code."
    10741113       );
    10751114    parameters[numberParameters++]=
     
    10821121    parameters[numberParameters-1].setLonghelp
    10831122      (
    1084        ""
     1123       "Normally Presolve does 5 passes but you may want to do less to make it\
     1124 more lightweight or do more if improvements are still being made.  As Presolve will return\
     1125 if nothing is being taken out, then you should not need to use this fine tuning."
    10851126       );
    10861127    parameters[numberParameters-1].setIntValue(preSolve);
     
    10951136    parameters[numberParameters-1].setLonghelp
    10961137      (
    1097        ""
     1138       "Perturbation helps to stop cycling, but Clp uses other measures for this.\
     1139  However large problems and especially ones with unit elements and unit rhs or costs\
     1140 benefit from perturbation.  Normally Clp tries to be intelligent, but you can switch this off.\
     1141  The Clp library has this off by default.  This program has it on."
    10981142       );
    10991143    parameters[numberParameters++]=
     
    11021146    parameters[numberParameters-1].setLonghelp
    11031147      (
    1104        ""
     1148       "Clp will go slightly faster if the matrix can be converted so that the elements are\
     1149 not stored and are known to be unit.  The main advantage is memory use.  Clp may automatically\
     1150 see if it can convert the problem so you should not need to use this."
    11051151       );
    11061152    parameters[numberParameters++]=
     
    11201166    parameters[numberParameters-1].append("part!ial");
    11211167    parameters[numberParameters-1].append("steep!est");
    1122     //parameters[numberParameters-1].append("change");
     1168    parameters[numberParameters-1].append("change");
    11231169    parameters[numberParameters-1].setLonghelp
    11241170      (
     
    11551201    parameters[numberParameters-1].setLonghelp
    11561202      (
    1157        ""
     1203       "This stops the execution of Clp, end, exit, quit and stop are synonyms"
    11581204       );
    11591205    parameters[numberParameters++]=
     
    11621208    parameters[numberParameters-1].setLonghelp
    11631209      (
     1210       "This will write an MPS format file to the given file name.  It will use the default\
     1211 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
     1212 is initialized to 'default.mps'."
    11641213       ""
    11651214       );
    11661215    parameters[numberParameters-1].setStringValue(restoreFile);
     1216    parameters[numberParameters++]=
     1217      ClpItem("reverse","Reverses sign of objective",
     1218              REVERSE,false);
     1219    parameters[numberParameters-1].setLonghelp
     1220      (
     1221       "Useful for testing if maximization works correctly"
     1222       );
    11671223    parameters[numberParameters++]=
    11681224      ClpItem("save!Model","Save model to binary file",
     
    11701226    parameters[numberParameters-1].setLonghelp
    11711227      (
     1228       "This will write an MPS format file to the given file name.  It will use the default\
     1229 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
     1230 is initialized to 'default.mps'."
    11721231       ""
    11731232       );
     
    11771236              "on",SCALING);
    11781237    parameters[numberParameters-1].append("off");
     1238    parameters[numberParameters-1].append("geo!metric");
     1239    parameters[numberParameters-1].append("auto!matic");
    11791240    parameters[numberParameters-1].setLonghelp
    11801241      (
     
    11941255    parameters[numberParameters-1].setLonghelp
    11951256      (
     1257       "This will write an MPS format file to the given file name.  It will use the default\
     1258 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
     1259 is initialized to 'default.mps'."
    11961260       ""
    11971261       );
     
    12171281    parameters[numberParameters-1].setLonghelp
    12181282      (
    1219        ""
     1283       "This stops the execution of Clp, end, exit, quit and stop are synonyms"
    12201284       );
    12211285    parameters[numberParameters++]=
     
    12301294    int numberGoodCommands=0;
    12311295    bool * goodModels = new bool[1];
    1232     int getNewMatrix=0;
    12331296   
    12341297   
     
    14261489                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    14271490              } else if (action==4) {
    1428                 ClpPrimalColumnSteepest steep(0);
     1491                ClpPrimalColumnSteepest steep(1);
    14291492                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    14301493              } else if (action==5) {
     
    14341497              break;
    14351498            case SCALING:
    1436               models[iModel].scaling(1-action);
     1499              if (action<2)
     1500                models[iModel].scaling(1-action);
     1501              else
     1502                models[iModel].scaling(action);
    14371503              break;
    14381504            case SPARSEFACTOR:
     
    14631529              break;
    14641530            case CRASH:
    1465               doCrash=-action;
     1531              doCrash=action;
    14661532              break;
    14671533            case MESSAGES:
     
    15381604            break;
    15391605          case PLUSMINUS:
    1540             getNewMatrix=1;
     1606            if (goodModels[iModel]) {
     1607              ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
     1608              ClpPackedMatrix* clpMatrix =
     1609                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
     1610              if (clpMatrix) {
     1611                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
     1612                if (newMatrix->getIndices()) {
     1613                  models[iModel].replaceMatrix(newMatrix);
     1614                  delete saveMatrix;
     1615                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
     1616                } else {
     1617                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
     1618                }
     1619              } else {
     1620                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
     1621              }
     1622            } else {
     1623              std::cout<<"** Current model not valid"<<std::endl;
     1624            }
    15411625            break;
    15421626          case NETWORK:
    1543             getNewMatrix=2;
     1627            if (goodModels[iModel]) {
     1628              ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
     1629              ClpPackedMatrix* clpMatrix =
     1630                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
     1631              if (clpMatrix) {
     1632                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
     1633                if (newMatrix->getIndices()) {
     1634                  models[iModel].replaceMatrix(newMatrix);
     1635                  delete saveMatrix;
     1636                  std::cout<<"Matrix converted to network matrix"<<std::endl;
     1637                } else {
     1638                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
     1639                }
     1640              } else {
     1641                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
     1642              }
     1643            } else {
     1644              std::cout<<"** Current model not valid"<<std::endl;
     1645            }
    15441646            break;
    15451647          case IMPORT:
     
    18381940          case ALLSLACK:
    18391941            models[iModel].createStatus();
     1942            {
     1943              // and do solution
     1944              int iColumn;
     1945              double * solution = models[iModel].primalColumnSolution();
     1946              const double * lower = models[iModel].columnLower();
     1947              const double * upper = models[iModel].columnUpper();
     1948              int numberColumns = models[iModel].numberColumns();
     1949              for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1950                if (lower[iColumn]>0.0) {
     1951                  solution[iColumn]=lower[iColumn];
     1952                } else if (upper[iColumn]<0.0) {
     1953                  solution[iColumn]=upper[iColumn];
     1954                } else {
     1955                  solution[iColumn]=0.0;
     1956                }
     1957              }
     1958            }
     1959            break;
     1960          case REVERSE:
     1961            if (goodModels[iModel]) {
     1962              int iColumn;
     1963              int numberColumns=models[iModel].numberColumns();
     1964              double * dualColumnSolution =
     1965                models[iModel].dualColumnSolution();
     1966              ClpObjective * obj = models[iModel].objectiveAsObject();
     1967              assert(dynamic_cast<ClpLinearObjective *> (obj));
     1968              double offset;
     1969              double * objective = obj->gradient(NULL,offset,true);
     1970              for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1971                dualColumnSolution[iColumn] = dualColumnSolution[iColumn];
     1972                objective[iColumn] = -objective[iColumn];
     1973              }
     1974              int iRow;
     1975              int numberRows=models[iModel].numberRows();
     1976              double * dualRowSolution =
     1977                models[iModel].dualRowSolution();
     1978              for (iRow=0;iRow<numberRows;iRow++)
     1979                dualRowSolution[iRow] = dualRowSolution[iRow];
     1980            }
    18401981            break;
    18411982          case DIRECTORY:
  • branches/pre/Test/unitTest.cpp

    r210 r212  
    717717        double value = result[iRow];
    718718        assert(eq(value,sol[iRow]));
    719         assert(value<upper[iRow]+1.0e-8);
    720         assert(value>lower[iRow]-1.0e-8);
     719        assert(value<upper[iRow]+2.0e-8);
     720        assert(value>lower[iRow]-2.0e-8);
    721721      }
    722722      CoinFillN ( result, numberRows,0.0);
  • branches/pre/include/ClpPrimalColumnSteepest.hpp

    r210 r212  
    2626 
    2727  /** Returns pivot column, -1 if none.
    28       updateArray has cost updates (also use pivotRow_ from last iteration)
     28      updateArray has cost updates (also use pivotRow_ from last iteration).
     29      Parts of operation split out into seperate functions for
     30      profiling and speed
    2931  */
    3032  virtual int pivotColumn(CoinIndexedVector * updates,
     
    3335                          CoinIndexedVector * spareColumn1,
    3436                          CoinIndexedVector * spareColumn2);
     37  /// For quadratic or funny nonlinearities
     38  int pivotColumnOldMethod(CoinIndexedVector * updates,
     39                          CoinIndexedVector * spareRow1,
     40                          CoinIndexedVector * spareRow2,
     41                          CoinIndexedVector * spareColumn1,
     42                          CoinIndexedVector * spareColumn2);
     43  /// Just update djs
     44  void justDjs(CoinIndexedVector * updates,
     45               CoinIndexedVector * spareRow1,
     46               CoinIndexedVector * spareRow2,
     47               CoinIndexedVector * spareColumn1,
     48               CoinIndexedVector * spareColumn2);
     49  /// Update djs, weights for Devex using djs
     50  void djsAndDevex(CoinIndexedVector * updates,
     51               CoinIndexedVector * spareRow1,
     52               CoinIndexedVector * spareRow2,
     53               CoinIndexedVector * spareColumn1,
     54               CoinIndexedVector * spareColumn2);
     55  /// Update djs, weights for Steepest using djs
     56  void djsAndSteepest(CoinIndexedVector * updates,
     57               CoinIndexedVector * spareRow1,
     58               CoinIndexedVector * spareRow2,
     59               CoinIndexedVector * spareColumn1,
     60               CoinIndexedVector * spareColumn2);
     61  /// Update djs, weights for Devex using pivot row
     62  void djsAndDevex2(CoinIndexedVector * updates,
     63               CoinIndexedVector * spareRow1,
     64               CoinIndexedVector * spareRow2,
     65               CoinIndexedVector * spareColumn1,
     66               CoinIndexedVector * spareColumn2);
     67  /// Update djs, weights for Steepest using pivot row
     68  void djsAndSteepest2(CoinIndexedVector * updates,
     69               CoinIndexedVector * spareRow1,
     70               CoinIndexedVector * spareRow2,
     71               CoinIndexedVector * spareColumn1,
     72               CoinIndexedVector * spareColumn2);
     73  /// Update weights for Devex
     74  void justDevex(CoinIndexedVector * updates,
     75               CoinIndexedVector * spareRow1,
     76               CoinIndexedVector * spareRow2,
     77               CoinIndexedVector * spareColumn1,
     78               CoinIndexedVector * spareColumn2);
     79  /// Update weights for Steepest
     80  void justSteepest(CoinIndexedVector * updates,
     81               CoinIndexedVector * spareRow1,
     82               CoinIndexedVector * spareRow2,
     83               CoinIndexedVector * spareColumn1,
     84               CoinIndexedVector * spareColumn2);
    3585
    3686  /// Updates weights - part 1 - also checks accuracy
     
    68118      0 is exact devex, 1 full steepest, 2 is partial exact devex
    69119      3 switches between 0 and 2 depending on factorization
    70       4 starts as partial dantzig but then may switch between 0 and 2.
     120      4 starts as partial dantzig/devex but then may switch between 0 and 2.
    71121      By partial exact devex is meant that the weights are updated as normal
    72122      but only part of the nonbasic variables are scanned. 
     
    130180      0 is exact devex, 1 full steepest, 2 is partial exact devex
    131181      3 switches between 0 and 2 depending on factorization
    132       4 starts as partial dantzig but then may switch between 0 and 2.
     182      4 starts as partial dantzig/devex but then may switch between 0 and 2.
    133183      By partial exact devex is meant that the weights are updated as normal
    134184      but only part of the nonbasic variables are scanned. 
  • branches/pre/include/ClpSimplex.hpp

    r210 r212  
    161161  /// Passes in factorization
    162162  void setFactorization( ClpFactorization & factorization);
    163   /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3 dynamic(later)
     163  /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
    164164  void scaling(int mode=1);
    165165  /// Gets scalingFlag
     
    360360    void restoreData(ClpDataSave saved);
    361361  int internalFactorize(int solveType);
     362  /// Clean up status
     363  void cleanStatus();
    362364  /// Factorizes using current basis. For external use
    363365  int factorize();
     
    545547  /// Sanity check on input rim data (after scaling) - returns true if okay
    546548  bool sanityCheck();
    547   /** Get next superbasic (primal) or next free (dual), -1 if none */
    548   int nextSuperBasic();
    549549  //@}
    550550  public:
     
    866866  /// Column scale factors
    867867  double * columnScale_;
    868   /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3 dynamic
     868  /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic
    869869  int scalingFlag_;
    870870  /// Number of times code has tentatively thought optimal
  • branches/pre/include/ClpSimplexDual.hpp

    r205 r212  
    265265  */
    266266  int pivotResult();
     267  /** Get next free , -1 if none */
     268  int nextSuperBasic();
    267269 
    268270  //@}
  • branches/pre/include/ClpSimplexPrimal.hpp

    r210 r212  
    138138      -5 Looks unbounded
    139139      +3 max iterations
     140     
     141      valuesOption has original value of valuesPass
    140142   */
    141   int whileIterating();
     143  int whileIterating(int valuesOption);
    142144
    143145  /** Do last half of an iteration.  This is split out so people can
     
    217219  /// Unflag all variables and return number unflagged
    218220  int unflag();
     221  /** Get next superbasic -1 if none,
     222      Normal type is 1
     223      If type is 3 then initializes sorted list
     224      if 2 uses list.
     225  */
     226  int nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray);
    219227
    220228  /// Create primal ray
  • branches/pre/include/Idiot.hpp

    r208 r212  
    139139  inline void setReduceIterations(int value)
    140140  { maxBigIts_ = value;};
    141   /// Amount of information - default of 1 should be okay */
     141  /// Amount of information - default of 1 should be okay
    142142  inline int getLogLevel() const
    143143  { return logLevel_;};
    144144  inline void setLogLevel(int value)
    145145  { logLevel_ = value;};
     146  /// How lightweight - 0 not, 1 yes
     147  inline int getLightweight() const
     148  { return lightWeight_;};
     149  inline void setLightweight(int value)
     150  { lightWeight_ = value;};
    146151  //@}
    147152
     
    211216                  2048 - keep lambda across mu change
    212217                  4096 - return best solution (not last found) */
     218  int lightWeight_; // 0 - normal, 1 lightweight
    213219};
    214220#endif
Note: See TracChangeset for help on using the changeset viewer.