Changeset 259


Ignore:
Timestamp:
Nov 20, 2003 10:48:15 AM (16 years ago)
Author:
forrest
Message:

Hopefully improving primal and presolve

Location:
trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpMessage.cpp

    r256 r259  
    1919  {CLP_SIMPLEX_ERROR,4,1,"Stopped due to errors - objective value %g"},
    2020  //{CLP_SIMPLEX_STATUS,5,1,"%d  Objective %g%? Primal infeas %g (%d)%? Dual infeas %g (%d)%? without free dual infeas (%d)"},
    21   {CLP_SIMPLEX_STATUS,5,1,"%d  Obj %g%? Primal inf %g (%d)%? Dual inf %g (%d)%? w.u. free dual inf (%d)"},
     21  {CLP_SIMPLEX_STATUS,5,1,"%d  Obj %g%? Primal inf %g (%d)%? Dual inf %g (%d)%? w.o. free dual inf (%d)"},
    2222  {CLP_DUAL_BOUNDS,25,3,"Looking optimal checking bounds with %g"},
    2323  {CLP_SIMPLEX_ACCURACY,6,3,"Primal error %g, dual error %g"},
  • trunk/ClpPresolve.cpp

    r247 r259  
    439439    // Check number rows dropped
    440440    int lastDropped=0;
     441    prob->pass_=0;
    441442    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
    442443#ifdef PRESOLVE_SUMMARY
     
    458459            paction_ = slack_doubleton_action::presolve(prob, paction_,
    459460                                                        notFinished);
     461          if (prob->status_)
     462            break;
     463        }
     464        if (dual&&whichPass==1) {
     465          // this can also make E rows so do one bit here
     466          paction_ = remove_dual_action::presolve(prob, paction_);
    460467          if (prob->status_)
    461468            break;
  • trunk/ClpPrimalColumnDantzig.cpp

    r225 r259  
    6969  double * updateBy;
    7070  double * reducedCost;
    71   double dj = model_->dualIn();
    7271
    7372  bool anyUpdates;
     
    7574  if (updates->getNumElements()) {
    7675    anyUpdates=true;
    77     // add in pivot contribution
    78     if (model_->pivotRow()>=0)
    79       updates->add(model_->pivotRow(),-dj);
    80   } else if (model_->pivotRow()>=0) {
    81     updates->insert(model_->pivotRow(),-dj);
    82     anyUpdates=true;
    8376  } else {
    8477    // sub flip - nothing to do
    8578    anyUpdates=false;
    8679  }
    87   if (!updates->getNumElements())
    88     anyUpdates=false; // in case dj 0.0 (for quadratic)
    8980  if (anyUpdates) {
    9081    model_->factorization()->updateColumnTranspose(spareRow2,updates);
     
    109100        int iSequence = index[j];
    110101        double value = reducedCost[iSequence];
    111         value -= updateBy[iSequence];
    112         updateBy[iSequence]=0.0;
     102        value -= updateBy[j];
     103        updateBy[j]=0.0;
    113104        reducedCost[iSequence] = value;
    114105      }
  • trunk/ClpPrimalColumnSteepest.cpp

    r256 r259  
    1414
    1515
    16 // bias for free variables
    17 #define FREE_BIAS 1.0e1
    18 // Acceptance criteria for free variables
    19 #define FREE_ACCEPT 1.0e2
    2016//#############################################################################
    2117// Constructors / Destructor / Assignment
     
    3834    pivotSequence_(-1),
    3935    savedPivotSequence_(-1),
    40     savedSequenceOut_(-1)
     36    savedSequenceOut_(-1),
     37    sizeFactorization_(0)
    4138{
    4239  type_=2+64*mode;
     
    5653  savedPivotSequence_ = rhs.savedPivotSequence_;
    5754  savedSequenceOut_ = rhs.savedSequenceOut_;
     55  sizeFactorization_ = rhs.sizeFactorization_;
    5856  devex_ = rhs.devex_;
    5957  if (rhs.infeasible_) {
     
    112110    savedPivotSequence_ = rhs.savedPivotSequence_;
    113111    savedSequenceOut_ = rhs.savedSequenceOut_;
     112    sizeFactorization_ = rhs.sizeFactorization_;
    114113    devex_ = rhs.devex_;
    115114    delete [] weights_;
     
    152151#define ADD_ONE 1.0
    153152// Returns pivot column, -1 if none
     153/*      The Packed CoinIndexedVector updates has cost updates - for normal LP
     154        that is just +-weight where a feasibility changed.  It also has
     155        reduced cost from last iteration in pivot row*/
    154156int
    155157ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
     
    162164  if (model_->nonLinearCost()->lookBothWays()||model_->algorithm()==2) {
    163165    // Do old way
     166    updates->expand();
    164167    return pivotColumnOldMethod(updates,spareRow1,spareRow2,
    165168                                spareColumn1,spareColumn2);
     
    167170  int number=0;
    168171  int * index;
    169   // dj could be very small (or even zero - take care)
    170   double dj = model_->dualIn();
    171172  double tolerance=model_->currentDualTolerance();
    172173  // we can't really trust infeasibilities if there is dual error
     
    196197     10 - can go to mini-sprint
    197198  */
    198   if (updates->getNumElements()) {
     199  if (updates->getNumElements()>1) {
    199200    // would have to have two goes for devex, three for steepest
    200201    anyUpdates=2;
    201     // add in pivot contribution
    202     if (pivotRow>=0)
    203       updates->add(pivotRow,-dj);
    204   } else if (pivotRow>=0) {
    205     if (fabs(dj)>1.0e-15) {
    206       // some dj
    207       updates->insert(pivotRow,-dj);
    208       if (fabs(dj)>1.0e-6) {
    209         // reasonable size
    210         anyUpdates=1;
    211       } else {
    212         // too small
    213         anyUpdates=2;
    214       }
     202  } else if (updates->getNumElements()) {
     203    if (updates->getIndices()[0]==pivotRow&&fabs(updates->denseVector()[0])>1.0e-6) {
     204      // reasonable size
     205      anyUpdates=1;
     206      //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
     207      //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
    215208    } else {
    216       // zero dj
    217       anyUpdates=-1;
     209      // too small
     210      anyUpdates=2;
    218211    }
    219212  } else if (pivotSequence_>=0){
     
    225218  }
    226219  int sequenceOut = model_->sequenceOut();
     220  if (switchType==5) {
     221    // If known matrix then we will do partial pricing
     222    if (model_->clpMatrix()->canDoPartialPricing()) {
     223      pivotSequence_=-1;
     224      pivotRow=-1;
     225      // See if to switch
     226      int numberRows = model_->numberRows();
     227      int numberWanted=0;
     228      int numberColumns = model_->numberColumns();
     229      double ratio = (double) sizeFactorization_/(double) numberRows;
     230      // Number of dual infeasibilities at last invert
     231      int numberDual = model_->numberDualInfeasibilities();
     232      int numberLook = min(numberDual,numberColumns/10);
     233      if (ratio<1.0) {
     234        numberWanted = 100;
     235        numberWanted = max(numberWanted,numberLook/20);
     236      } else if (ratio<3.0) {
     237        numberWanted = 500;
     238        numberWanted = max(numberWanted,numberLook/15);
     239      } else if (ratio<4.0) {
     240        numberWanted = 1000;
     241        numberWanted = max(numberWanted,numberLook/10);
     242      } else {
     243        switchType=4;
     244        // initialize
     245        numberSwitched_++;
     246        // Make sure will re-do
     247        delete [] weights_;
     248        weights_=NULL;
     249        model_->computeDuals(NULL);
     250        saveWeights(model_,4);
     251        anyUpdates=0;
     252        printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
     253      }
     254      if (switchType==5) {
     255        if (model_->numberIterations()%1000==0)
     256          printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
     257        // Update duals and row djs
     258        // Do partial pricing
     259        return partialPricing(updates,spareRow2,
     260                              numberWanted);
     261      }
     262    }
     263  }
    227264  if (switchType==5) {
    228265    if (anyUpdates>0) {
     
    261298    }
    262299  }
     300  alternateWeights_->checkClear();
    263301  // make sure outgoing from last iteration okay
    264302  if (sequenceOut>=0) {
     
    317355    pivotRow=-1;
    318356    // See if to switch
    319     int numberElements = model_->factorization()->numberElements();
    320357    int numberRows = model_->numberRows();
    321358    // ratio is done on number of columns here
    322     //double ratio = (double) numberElements/(double) numberColumns;
    323     double ratio = (double) numberElements/(double) numberRows;
    324     //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
    325     if (ratio<0.1) {
     359    //double ratio = (double) sizeFactorization_/(double) numberColumns;
     360    double ratio = (double) sizeFactorization_/(double) numberRows;
     361    //double ratio = (double) sizeFactorization_/(double) model_->clpMatrix()->getNumElements();
     362    if (ratio<1.0) {
    326363      numberWanted = max(100,number/200);
    327     } else if (ratio<0.15) {
     364    } else if (ratio<2.0-1.0) {
    328365      numberWanted = max(500,number/40);
    329     } else if (ratio<0.2) {
     366    } else if (ratio<4.0-3.0) {
    330367      numberWanted = max(2000,number/10);
    331368      numberWanted = max(numberWanted,numberColumns/30);
     
    338375      weights_=NULL;
    339376      saveWeights(model_,4);
    340       printf("switching to devex %d nel ratio %g\n",numberElements,ratio);
     377      printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
    341378    }
    342379    //if (model_->numberIterations()%1000==0)
    343     //printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
    344   }
    345   int numberElements = model_->factorization()->numberElements();
     380    //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
     381  }
    346382  int numberRows = model_->numberRows();
    347383  // ratio is done on number of rows here
    348   double ratio = (double) numberElements/(double) numberRows;
     384  double ratio = (double) sizeFactorization_/(double) numberRows;
    349385  if(switchType==4) {
    350386    // Still in devex mode
    351387    // Go to steepest if lot of iterations?
    352     if (ratio<1.0) {
    353       numberWanted = max(2000,number/20);
    354     } else if (ratio<4.0) {
     388    if (ratio<5.0) {
    355389      numberWanted = max(2000,number/10);
    356390      numberWanted = max(numberWanted,numberColumns/20);
     391    } else if (ratio<7.0) {
     392      numberWanted = max(2000,number/5);
     393      numberWanted = max(numberWanted,numberColumns/10);
    357394    } else {
    358395      // we can zero out
     
    368405      weights_=NULL;
    369406      saveWeights(model_,4);
    370       printf("switching to exact %d nel ratio %g\n",numberElements,ratio);
     407      printf("switching to exact %d nel ratio %g\n",sizeFactorization_,ratio);
    371408      updates->clear();
    372409    }
    373410    if (model_->numberIterations()%1000==0)
    374     printf("numels %d ratio %g wanted %d type x\n",numberElements,ratio,numberWanted);
     411    printf("numels %d ratio %g wanted %d type x\n",sizeFactorization_,ratio,numberWanted);
    375412  }
    376413  if (switchType<4) {
     
    399436    }
    400437    //if (model_->numberIterations()%1000==0)
    401     //printf("numels %d ratio %g wanted %d type %d\n",numberElements,ratio,numberWanted,
     438    //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
    402439    //switchType);
    403440  }
     
    411448  number = infeasible_->getNumElements();
    412449  // Re-sort infeasible every 100 pivots
     450  // Not a good idea
    413451  if (0&&model_->factorization()->pivots()>0&&
    414452      (model_->factorization()->pivots()%100)==0) {
     
    547585  int pivotRow = model_->pivotRow();
    548586  double * infeas = infeasible_->denseVector();
     587  //updates->scanAndPack();
    549588  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    550589 
    551   // put row of tableau in rowArray and columnArray
     590  // put row of tableau in rowArray and columnArray (packed mode)
    552591  model_->clpMatrix()->transposeTimes(model_,-1.0,
    553592                                      updates,spareColumn2,spareColumn1);
     
    573612      int iSequence = index[j];
    574613      double value = reducedCost[iSequence];
    575       value -= updateBy[iSequence];
    576       updateBy[iSequence]=0.0;
     614      value -= updateBy[j];
     615      updateBy[j]=0.0;
    577616      reducedCost[iSequence] = value;
    578617      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     
    643682  double * updateBy;
    644683  double * reducedCost;
    645   // dj could be very small (or even zero - take care)
    646   double dj = model_->dualIn();
    647684  double tolerance=model_->currentDualTolerance();
    648685  // we can't really trust infeasibilities if there is dual error
     
    656693  pivotSequence_=-1;
    657694  double * infeas = infeasible_->denseVector();
     695  //updates->scanAndPack();
    658696  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    659697  // and we can see if reference
     
    668706    outgoingWeight=weights_[sequenceOut];
    669707   
    670   // put row of tableau in rowArray and columnArray
     708  double scaleFactor = 1.0/updates->denseVector()[0]; // as formula is with 1.0
     709  // put row of tableau in rowArray and columnArray (packed mode)
    671710  model_->clpMatrix()->transposeTimes(model_,-1.0,
    672711                                      updates,spareColumn2,spareColumn1);
     
    674713  double * weight;
    675714  int numberColumns = model_->numberColumns();
    676   double scaleFactor = -1.0/dj; // as formula is with 1.0
    677715  // rows
    678716  reducedCost=model_->djRegion(0);
     
    690728    int iSequence = index[j];
    691729    double value = reducedCost[iSequence];
    692     double value2 = updateBy[iSequence];
    693     updateBy[iSequence]=0.0;
     730    double value2 = updateBy[j];
     731    updateBy[j]=0.0;
    694732    value -= value2;
    695733    reducedCost[iSequence] = value;
     
    778816    int iSequence = index[j];
    779817    double value = reducedCost[iSequence];
    780     double value2 = updateBy[iSequence];
     818    double value2 = updateBy[j];
    781819    value -= value2;
    782     updateBy[iSequence]=0.0;
     820    updateBy[j]=0.0;
    783821    reducedCost[iSequence] = value;
    784822    ClpSimplex::Status status = model_->getStatus(iSequence);
     
    883921  double * updateBy;
    884922  double * reducedCost;
    885   // dj could be very small (or even zero - take care)
    886   double dj = model_->dualIn();
    887923  double tolerance=model_->currentDualTolerance();
    888924  // we can't really trust infeasibilities if there is dual error
     
    896932  pivotSequence_=-1;
    897933  double * infeas = infeasible_->denseVector();
     934  double scaleFactor = 1.0/updates->denseVector()[0]; // as formula is with 1.0
     935  //updates->scanAndPack();
    898936  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    899 
     937  //alternateWeights_->scanAndPack();
    900938  model_->factorization()->updateColumnTranspose(spareRow2,
    901939                                                 alternateWeights_);
     
    911949    outgoingWeight=weights_[sequenceOut];
    912950   
    913   // put row of tableau in rowArray and columnArray
     951  // put row of tableau in rowArray and columnArray (packed)
    914952  model_->clpMatrix()->transposeTimes(model_,-1.0,
    915953                                      updates,spareColumn2,spareColumn1);
     
    919957                                              spareColumn1,
    920958                                              spareRow2);
    921   assert (spareRow2->packedMode());
    922959  // update weights
    923960  double * weight;
    924961  double * other = alternateWeights_->denseVector();
    925962  int numberColumns = model_->numberColumns();
    926   double scaleFactor = -1.0/dj; // as formula is with 1.0
    927963  // rows
    928964  reducedCost=model_->djRegion(0);
     
    941977    int iSequence = index[j];
    942978    double value = reducedCost[iSequence];
    943     double value2 = updateBy[iSequence];
     979    double value2 = updateBy[j];
    944980    value -= value2;
    945981    reducedCost[iSequence] = value;
    946982    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    947     updateBy[iSequence]=0.0;
     983    modification = other[iSequence];
     984    updateBy[j]=0.0;
    948985   
    949986    switch(status) {
     
    958995      // row has -1
    959996      pivot = value2*scaleFactor;
    960       modification = other[iSequence];
    961997      pivotSquared = pivot * pivot;
    962998     
     
    9911027      // row has -1
    9921028      pivot = value2*scaleFactor;
    993       modification = other[iSequence];
    9941029      pivotSquared = pivot * pivot;
    9951030     
     
    10221057      // row has -1
    10231058      pivot = value2*scaleFactor;
    1024       modification = other[iSequence];
    10251059      pivotSquared = pivot * pivot;
    10261060     
     
    10501084    }
    10511085  }
    1052  
     1086  alternateWeights_->clear();
    10531087  // columns
    10541088  weight = weights_;
     
    10681102    int iSequence = index[j];
    10691103    double value = reducedCost[iSequence];
    1070     double value2 = updateBy[iSequence];
    1071     updateBy[iSequence]=0.0;
     1104    double value2 = updateBy[j];
     1105    updateBy[j]=0.0;
    10721106    value -= value2;
    10731107    reducedCost[iSequence] = value;
     
    11841218  // make sure infeasibility on incoming is 0.0
    11851219  infeasible_->zero(sequenceIn);
    1186   alternateWeights_->clear();
    11871220  spareRow2->setNumElements(0);
    11881221  //#define SOME_DEBUG_1
     
    12251258  int pivotRow = model_->pivotRow();
    12261259  double * infeas = infeasible_->denseVector();
     1260  //updates->scanAndPack();
    12271261  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    12281262 
     
    12511285      int iSequence = index[j];
    12521286      double value = reducedCost[iSequence];
    1253       value -= updateBy[iSequence];
     1287      value -= updateBy[j];
     1288      updateBy[j]=0.0;
    12541289      reducedCost[iSequence] = value;
    12551290      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     
    12991334    }
    13001335  }
    1301   // we can zero out as will have to get pivot row
    1302   updates->clear();
    1303   spareColumn1->clear();
     1336  // They are empty
     1337  updates->setNumElements(0);
     1338  spareColumn1->setNumElements(0);
    13041339  // make sure infeasibility on incoming is 0.0
    13051340  int sequenceIn = model_->sequenceIn();
     
    13351370    double * weight;
    13361371    int numberColumns = model_->numberColumns();
    1337     double scaleFactor = -1.0/dj; // as formula is with 1.0
    13381372    // rows
    13391373    number = updates->getNumElements();
     
    13471381      double thisWeight = weight[iSequence];
    13481382      // row has -1
    1349       double pivot = updateBy[iSequence]*scaleFactor;
     1383      double pivot = - updateBy[iSequence];
    13501384      updateBy[iSequence]=0.0;
    13511385      double value = pivot * pivot*devex_;
     
    13571391    // columns
    13581392    weight = weights_;
    1359    
    1360     scaleFactor = -scaleFactor;
    13611393   
    13621394    number = spareColumn1->getNumElements();
     
    13671399      double thisWeight = weight[iSequence];
    13681400      // row has -1
    1369       double pivot = updateBy[iSequence]*scaleFactor;
     1401      double pivot = updateBy[iSequence];
    13701402      updateBy[iSequence]=0.0;
    13711403      double value = pivot * pivot*devex_;
     
    14181450  int pivotRow = model_->pivotRow();
    14191451  double * infeas = infeasible_->denseVector();
     1452  //updates->scanAndPack();
    14201453  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    14211454 
     
    14441477      int iSequence = index[j];
    14451478      double value = reducedCost[iSequence];
    1446       value -= updateBy[iSequence];
     1479      value -= updateBy[j];
     1480      updateBy[j]=0.0;
    14471481      reducedCost[iSequence] = value;
    14481482      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     
    14941528  // we can zero out as will have to get pivot row
    14951529  // ***** move
    1496   updates->clear();
    1497   spareColumn1->clear();
     1530  updates->setNumElements(0);
     1531  spareColumn1->setNumElements(0);
    14981532  if (pivotRow>=0) {
    14991533    // make sure infeasibility on incoming is 0.0
     
    15231557    spareColumn1->setNumElements(0);
    15241558    // might as well set dj to 1
    1525     dj = 1.0;
    1526     updates->insert(pivotRow,-dj);
     1559    dj = -1.0;
     1560    updates->createPacked(1,&pivotRow,&dj);
    15271561    model_->factorization()->updateColumnTranspose(spareRow2,updates);
    15281562    // put row of tableau in rowArray and columnArray
     
    15321566    double * other = alternateWeights_->denseVector();
    15331567    int numberColumns = model_->numberColumns();
    1534     double scaleFactor = -1.0/dj; // as formula is with 1.0
    15351568    // rows
    15361569    number = updates->getNumElements();
     
    15421575      // Exact
    15431576      // now update weight update array
     1577      //alternateWeights_->scanAndPack();
    15441578      model_->factorization()->updateColumnTranspose(spareRow2,
    15451579                                                     alternateWeights_);
     1580      // Exact
     1581      // get subset which have nonzero tableau elements
     1582      model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
     1583                                                spareColumn1,
     1584                                                spareColumn2);
    15461585      for (j=0;j<number;j++) {
    15471586        int iSequence = index[j];
    15481587        double thisWeight = weight[iSequence];
    15491588        // row has -1
    1550         double pivot = updateBy[iSequence]*scaleFactor;
    1551         updateBy[iSequence]=0.0;
     1589        double pivot = - updateBy[j];
     1590        updateBy[j]=0.0;
    15521591        double modification = other[iSequence];
    15531592        double pivotSquared = pivot * pivot;
     
    15751614        double thisWeight = weight[iSequence];
    15761615        // row has -1
    1577         double pivot = updateBy[iSequence]*scaleFactor;
    1578         updateBy[iSequence]=0.0;
     1616        double pivot = -updateBy[j];
     1617        updateBy[j]=0.0;
    15791618        double value = pivot * pivot*devex_;
    15801619        if (reference(iSequence+numberColumns))
     
    15861625    // columns
    15871626    weight = weights_;
    1588    
    1589     scaleFactor = -scaleFactor;
    15901627   
    15911628    number = spareColumn1->getNumElements();
     
    15941631    if (mode_<4||numberSwitched_>1||mode_>=10) {
    15951632      // Exact
    1596       // get subset which have nonzero tableau elements
    1597       model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
    1598                                                 spareColumn1,
    1599                                                 spareColumn2);
    1600       assert (spareColumn2->packedMode());
    16011633      double * updateBy2 = spareColumn2->denseVector();
    16021634      for (j=0;j<number;j++) {
    16031635        int iSequence = index[j];
    16041636        double thisWeight = weight[iSequence];
    1605         double pivot = updateBy[iSequence]*scaleFactor;
    1606         updateBy[iSequence]=0.0;
     1637        double pivot = updateBy[j];
     1638        updateBy[j]=0.0;
    16071639        double modification = updateBy2[j];
    16081640        updateBy2[j]=0.0;
     
    16301662        double thisWeight = weight[iSequence];
    16311663        // row has -1
    1632         double pivot = updateBy[iSequence]*scaleFactor;
    1633         updateBy[iSequence]=0.0;
     1664        double pivot = updateBy[j];
     1665        updateBy[j]=0.0;
    16341666        double value = pivot * pivot*devex_;
    16351667        if (reference(iSequence))
     
    17021734  pivotSequence_=-1;
    17031735  // might as well set dj to 1
    1704   dj = 1.0;
    1705   updates->insert(pivotRow,-dj);
     1736  dj = -1.0;
     1737  updates->createPacked(1,&pivotRow,&dj);
    17061738  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    17071739  // put row of tableau in rowArray and columnArray
     
    17101742  double * weight;
    17111743  int numberColumns = model_->numberColumns();
    1712   double scaleFactor = -1.0/dj; // as formula is with 1.0
    17131744  // rows
    17141745  number = updates->getNumElements();
     
    17231754    double thisWeight = weight[iSequence];
    17241755    // row has -1
    1725     double pivot = updateBy[iSequence]*scaleFactor;
    1726     updateBy[iSequence]=0.0;
     1756    double pivot = - updateBy[j];
     1757    updateBy[j]=0.0;
    17271758    double value = pivot * pivot*devex_;
    17281759    if (reference(iSequence+numberColumns))
     
    17331764  // columns
    17341765  weight = weights_;
    1735  
    1736   scaleFactor = -scaleFactor;
    17371766 
    17381767  number = spareColumn1->getNumElements();
     
    17441773    double thisWeight = weight[iSequence];
    17451774    // row has -1
    1746     double pivot = updateBy[iSequence]*scaleFactor;
    1747     updateBy[iSequence]=0.0;
     1775    double pivot = updateBy[j];
     1776    updateBy[j]=0.0;
    17481777    double value = pivot * pivot*devex_;
    17491778    if (reference(iSequence))
     
    18161845  //spareColumn1->setNumElements(0);
    18171846  // might as well set dj to 1
    1818   dj = 1.0;
    1819   updates->insert(pivotRow,-dj);
     1847  dj = -1.0;
     1848  updates->createPacked(1,&pivotRow,&dj);
    18201849  model_->factorization()->updateColumnTranspose(spareRow2,updates);
    18211850  // put row of tableau in rowArray and columnArray
     
    18251854  double * other = alternateWeights_->denseVector();
    18261855  int numberColumns = model_->numberColumns();
    1827   double scaleFactor = -1.0/dj; // as formula is with 1.0
    18281856  // rows
    18291857  number = updates->getNumElements();
     
    18341862  // Exact
    18351863  // now update weight update array
     1864  //alternateWeights_->scanAndPack();
    18361865  model_->factorization()->updateColumnTranspose(spareRow2,
    18371866                                                 alternateWeights_);
     1867  // get subset which have nonzero tableau elements
     1868  model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
     1869                                            spareColumn1,
     1870                                            spareColumn2);
    18381871  for (j=0;j<number;j++) {
    18391872    int iSequence = index[j];
    18401873    double thisWeight = weight[iSequence];
    18411874    // row has -1
    1842     double pivot = updateBy[iSequence]*scaleFactor;
    1843     updateBy[iSequence]=0.0;
     1875    double pivot = -updateBy[j];
     1876    updateBy[j]=0.0;
    18441877    double modification = other[iSequence];
    18451878    double pivotSquared = pivot * pivot;
     
    18641897  weight = weights_;
    18651898 
    1866   scaleFactor = -scaleFactor;
    1867  
    18681899  number = spareColumn1->getNumElements();
    18691900  index = spareColumn1->getIndices();
    18701901  updateBy = spareColumn1->denseVector();
    18711902  // Exact
    1872   // get subset which have nonzero tableau elements
    1873   model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
    1874                                             spareColumn1,
    1875                                             spareColumn2);
    1876   assert (spareColumn2->packedMode());
    18771903  double * updateBy2 = spareColumn2->denseVector();
    18781904  for (j=0;j<number;j++) {
    18791905    int iSequence = index[j];
    18801906    double thisWeight = weight[iSequence];
    1881     double pivot = updateBy[iSequence]*scaleFactor;
    1882     updateBy[iSequence]=0.0;
     1907    double pivot = updateBy[j];
     1908    updateBy[j]=0.0;
    18831909    double modification = updateBy2[j];
    18841910    updateBy2[j]=0.0;
     
    23502376    pivotRow=-1;
    23512377    // See if to switch
    2352     int numberElements = model_->factorization()->numberElements();
    23532378    int numberRows = model_->numberRows();
    23542379    // ratio is done on number of columns here
    2355     //double ratio = (double) numberElements/(double) numberColumns;
    2356     double ratio = (double) numberElements/(double) numberRows;
    2357     //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
     2380    //double ratio = (double) sizeFactorization_/(double) numberColumns;
     2381    double ratio = (double) sizeFactorization_/(double) numberRows;
     2382    //double ratio = (double) sizeFactorization_/(double) model_->clpMatrix()->getNumElements();
    23582383    if (ratio<0.1) {
    23592384      numberWanted = max(100,number/200);
     
    23712396      weights_=NULL;
    23722397      saveWeights(model_,4);
    2373       printf("switching to devex %d nel ratio %g\n",numberElements,ratio);
     2398      printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
    23742399      updates->clear();
    23752400    }
    23762401    if (model_->numberIterations()%1000==0)
    2377       printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     2402      printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
    23782403  }
    23792404  if(switchType==4) {
    23802405    // Still in devex mode
    2381     int numberElements = model_->factorization()->numberElements();
    23822406    int numberRows = model_->numberRows();
    23832407    // ratio is done on number of rows here
    2384     double ratio = (double) numberElements/(double) numberRows;
     2408    double ratio = (double) sizeFactorization_/(double) numberRows;
    23852409    // Go to steepest if lot of iterations?
    23862410    if (ratio<1.0) {
     
    24022426      weights_=NULL;
    24032427      saveWeights(model_,4);
    2404       printf("switching to exact %d nel ratio %g\n",numberElements,ratio);
     2428      printf("switching to exact %d nel ratio %g\n",sizeFactorization_,ratio);
    24052429      updates->clear();
    24062430    }
    24072431    if (model_->numberIterations()%1000==0)
    2408       printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     2432      printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
    24092433  }
    24102434  if (switchType<4) {
     
    24142438      numberWanted = max(2000,number/8);
    24152439    } else {
    2416       int numberElements = model_->factorization()->numberElements();
    2417       double ratio = (double) numberElements/(double) model_->numberRows();
     2440      double ratio = (double) sizeFactorization_/(double) model_->numberRows();
    24182441      if (ratio<1.0) {
    24192442        numberWanted = max(2000,number/20);
     
    25332556                                                spareColumn1,
    25342557                                                spareColumn2);
    2535       assert (spareColumn2->packedMode());
    25362558      double * updateBy2 = spareColumn2->denseVector();
    25372559      for (j=0;j<number;j++) {
     
    27222744  int numberColumns = model_->numberColumns();
    27232745  const int * pivotVariable = model_->pivotVariable();
     2746  bool doInfeasibilities=true;
    27242747  if (mode==1&&weights_) {
    27252748    if (pivotSequence_>=0) {
     
    27362759    // restore
    27372760    if (!weights_||state_==-1||mode==5) {
    2738       // initialize weights
    2739       delete [] weights_;
    2740       delete alternateWeights_;
    2741       weights_ = new double[numberRows+numberColumns];
    2742       alternateWeights_ = new CoinIndexedVector();
    2743       // enough space so can use it for factorization
    2744       alternateWeights_->reserve(numberRows+
    2745                                  model_->factorization()->maximumPivots());
    2746       initializeWeights();
    2747       // create saved weights
    2748       delete [] savedWeights_;
    2749       savedWeights_ = new double[numberRows+numberColumns];
    2750       memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
    2751              sizeof(double));
     2761      // Partial is only allowed with certain types of matrix
     2762      if (mode_!=4||numberSwitched_||!model_->clpMatrix()->canDoPartialPricing()) {
     2763        // initialize weights
     2764        delete [] weights_;
     2765        delete alternateWeights_;
     2766        weights_ = new double[numberRows+numberColumns];
     2767        alternateWeights_ = new CoinIndexedVector();
     2768        // enough space so can use it for factorization
     2769        alternateWeights_->reserve(numberRows+
     2770                                   model_->factorization()->maximumPivots());
     2771        initializeWeights();
     2772        // create saved weights
     2773        delete [] savedWeights_;
     2774        savedWeights_ = new double[numberRows+numberColumns];
     2775        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
     2776               sizeof(double));
     2777      } else {
     2778        // Partial pricing
     2779        // use region as somewhere to save non-fixed slacks
     2780        // set up infeasibilities
     2781        if (!infeasible_) {
     2782          infeasible_ = new CoinIndexedVector();
     2783          infeasible_->reserve(numberColumns+numberRows);
     2784        }
     2785        infeasible_->clear();
     2786        int number = model_->numberRows() + model_->numberColumns();
     2787        int iSequence;
     2788        int numberLook=0;
     2789        int * which = infeasible_->getIndices();
     2790        for (iSequence=model_->numberColumns();iSequence<number;iSequence++) {
     2791          ClpSimplex::Status status = model_->getStatus(iSequence);
     2792          if (status!=ClpSimplex::isFixed)
     2793            which[numberLook++]=iSequence;
     2794        }
     2795        infeasible_->setNumElements(numberLook);
     2796        doInfeasibilities=false;
     2797      }
    27522798      savedPivotSequence_=-2;
    27532799      savedSequenceOut_ = -2;
     
    28102856      delete [] temp;
    28112857    }
     2858    // Save size of factorization
     2859    if (!model->factorization()->pivots())
     2860      sizeFactorization_ = model_->factorization()->numberElements();
     2861    if(!doInfeasibilities)
     2862      return; // don't disturb infeasibilities
    28122863    infeasible_->clear();
    28132864    double tolerance=model_->currentDualTolerance();
     
    29493000
    29503001  devex_ =0.0;
    2951 
    2952   if (pivotRow>=0) {
    2953     if (switchType==1) {
    2954       for (i=0;i<number;i++) {
    2955         int iRow = which[i];
    2956         devex_ += work[iRow]*work[iRow];
    2957         newWork[iRow] = -2.0*work[iRow];
    2958       }
    2959       newWork[pivotRow] = -2.0*max(devex_,0.0);
    2960       devex_+=ADD_ONE;
    2961       weights_[sequenceOut]=1.0+ADD_ONE;
    2962       memcpy(newWhich,which,number*sizeof(int));
    2963       alternateWeights_->setNumElements(number);
     3002  // Can't create alternateWeights_ as packed as needed unpacked
     3003  if (!input->packedMode()) {
     3004    if (pivotRow>=0) {
     3005      if (switchType==1) {
     3006        for (i=0;i<number;i++) {
     3007          int iRow = which[i];
     3008          devex_ += work[iRow]*work[iRow];
     3009          newWork[iRow] = -2.0*work[iRow];
     3010        }
     3011        newWork[pivotRow] = -2.0*max(devex_,0.0);
     3012        devex_+=ADD_ONE;
     3013        weights_[sequenceOut]=1.0+ADD_ONE;
     3014        memcpy(newWhich,which,number*sizeof(int));
     3015        alternateWeights_->setNumElements(number);
     3016      } else {
     3017        if (mode_!=4||numberSwitched_>1) {
     3018          for (i=0;i<number;i++) {
     3019            int iRow = which[i];
     3020            int iPivot = pivotVariable[iRow];
     3021            if (reference(iPivot)) {
     3022              devex_ += work[iRow]*work[iRow];
     3023              newWork[iRow] = -2.0*work[iRow];
     3024              newWhich[newNumber++]=iRow;
     3025            }
     3026          }
     3027          if (!newWork[pivotRow]&&devex_>0.0)
     3028            newWhich[newNumber++]=pivotRow; // add if not already in
     3029          newWork[pivotRow] = -2.0*max(devex_,0.0);
     3030        } else {
     3031          for (i=0;i<number;i++) {
     3032            int iRow = which[i];
     3033            int iPivot = pivotVariable[iRow];
     3034            if (reference(iPivot))
     3035              devex_ += work[iRow]*work[iRow];
     3036          }
     3037        }
     3038        if (reference(sequenceIn)) {
     3039          devex_+=1.0;
     3040        } else {
     3041        }
     3042        if (reference(sequenceOut)) {
     3043          weights_[sequenceOut]=1.0+1.0;
     3044        } else {
     3045          weights_[sequenceOut]=1.0;
     3046        }
     3047        alternateWeights_->setNumElements(newNumber);
     3048      }
    29643049    } else {
    2965       if (mode_!=4||numberSwitched_>1) {
     3050      if (switchType==1) {
     3051        for (i=0;i<number;i++) {
     3052          int iRow = which[i];
     3053          devex_ += work[iRow]*work[iRow];
     3054        }
     3055        devex_ += ADD_ONE;
     3056      } else {
    29663057        for (i=0;i<number;i++) {
    29673058          int iRow = which[i];
     
    29693060          if (reference(iPivot)) {
    29703061            devex_ += work[iRow]*work[iRow];
    2971             newWork[iRow] = -2.0*work[iRow];
    2972             newWhich[newNumber++]=iRow;
    29733062          }
    29743063        }
    2975         if (!newWork[pivotRow])
    2976           newWhich[newNumber++]=pivotRow; // add if not already in
     3064        if (reference(sequenceIn))
     3065          devex_+=1.0;
     3066      }
     3067    }
     3068  } else {
     3069    // packed input
     3070    if (pivotRow>=0) {
     3071      if (switchType==1) {
     3072        for (i=0;i<number;i++) {
     3073          int iRow = which[i];
     3074          devex_ += work[i]*work[i];
     3075          newWork[iRow] = -2.0*work[i];
     3076        }
    29773077        newWork[pivotRow] = -2.0*max(devex_,0.0);
     3078        devex_+=ADD_ONE;
     3079        weights_[sequenceOut]=1.0+ADD_ONE;
     3080        memcpy(newWhich,which,number*sizeof(int));
     3081        alternateWeights_->setNumElements(number);
     3082      } else {
     3083        if (mode_!=4||numberSwitched_>1) {
     3084          for (i=0;i<number;i++) {
     3085            int iRow = which[i];
     3086            int iPivot = pivotVariable[iRow];
     3087            if (reference(iPivot)) {
     3088              devex_ += work[i]*work[i];
     3089              newWork[iRow] = -2.0*work[i];
     3090              newWhich[newNumber++]=iRow;
     3091            }
     3092          }
     3093          if (!newWork[pivotRow]&&devex_>0.0)
     3094            newWhich[newNumber++]=pivotRow; // add if not already in
     3095          newWork[pivotRow] = -2.0*max(devex_,0.0);
     3096        } else {
     3097          for (i=0;i<number;i++) {
     3098            int iRow = which[i];
     3099            int iPivot = pivotVariable[iRow];
     3100            if (reference(iPivot))
     3101              devex_ += work[i]*work[i];
     3102          }
     3103        }
     3104        if (reference(sequenceIn)) {
     3105          devex_+=1.0;
     3106        } else {
     3107        }
     3108        if (reference(sequenceOut)) {
     3109          weights_[sequenceOut]=1.0+1.0;
     3110        } else {
     3111          weights_[sequenceOut]=1.0;
     3112        }
     3113        alternateWeights_->setNumElements(newNumber);
     3114      }
     3115    } else {
     3116      if (switchType==1) {
     3117        for (i=0;i<number;i++) {
     3118          devex_ += work[i]*work[i];
     3119        }
     3120        devex_ += ADD_ONE;
    29783121      } else {
    29793122        for (i=0;i<number;i++) {
    29803123          int iRow = which[i];
    29813124          int iPivot = pivotVariable[iRow];
    2982           if (reference(iPivot))
    2983             devex_ += work[iRow]*work[iRow];
    2984         }
    2985       }
    2986       if (reference(sequenceIn)) {
    2987         devex_+=1.0;
    2988       } else {
    2989       }
    2990       if (reference(sequenceOut)) {
    2991         weights_[sequenceOut]=1.0+1.0;
    2992       } else {
    2993         weights_[sequenceOut]=1.0;
    2994       }
    2995       alternateWeights_->setNumElements(newNumber);
    2996     }
    2997   } else {
    2998     if (switchType==1) {
    2999       for (i=0;i<number;i++) {
    3000         int iRow = which[i];
    3001         devex_ += work[iRow]*work[iRow];
    3002       }
    3003       devex_ += ADD_ONE;
    3004     } else {
    3005       for (i=0;i<number;i++) {
    3006         int iRow = which[i];
    3007         int iPivot = pivotVariable[iRow];
    3008         if (reference(iPivot)) {
    3009           devex_ += work[iRow]*work[iRow];
    3010         }
    3011       }
    3012       if (reference(sequenceIn))
    3013         devex_+=1.0;
     3125          if (reference(iPivot)) {
     3126            devex_ += work[i]*work[i];
     3127          }
     3128        }
     3129        if (reference(sequenceIn))
     3130          devex_+=1.0;
     3131      }
    30143132    }
    30153133  }
     
    32933411  numberSwitched_=10;
    32943412}
     3413// Update djs doing partial pricing (dantzig)
     3414int
     3415ClpPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
     3416                                        CoinIndexedVector * spareRow2,
     3417                                        int numberWanted)
     3418{
     3419  int number=0;
     3420  int * index;
     3421  double * updateBy;
     3422  double * reducedCost;
     3423  double saveTolerance = model_->currentDualTolerance();
     3424  double tolerance=model_->currentDualTolerance();
     3425  // we can't really trust infeasibilities if there is dual error
     3426  // this coding has to mimic coding in checkDualSolution
     3427  double error = min(1.0e-3,model_->largestDualError());
     3428  // allow tolerance at least slightly bigger than standard
     3429  tolerance = tolerance +  error;
     3430  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     3431    // we can't really trust infeasibilities if there is dual error
     3432    double checkTolerance = 1.0e-8;
     3433    if (!model_->factorization()->pivots())
     3434      checkTolerance = 1.0e-6;
     3435    if (model_->largestDualError()>checkTolerance)
     3436      tolerance *= model_->largestDualError()/checkTolerance;
     3437    // But cap
     3438    tolerance = min(1000.0,tolerance);
     3439  }
     3440  if (model_->factorization()->pivots()&&model_->numberPrimalInfeasibilities())
     3441    tolerance = max(tolerance,1.0e-10*model_->infeasibilityCost());
     3442  // So partial pricing can use
     3443  model_->setCurrentDualTolerance(tolerance);
     3444  model_->factorization()->updateColumnTranspose(spareRow2,updates);
     3445  int numberColumns = model_->numberColumns();
     3446
     3447  // Rows
     3448  reducedCost=model_->djRegion(0);
     3449  int addSequence;
     3450   
     3451  number = updates->getNumElements();
     3452  index = updates->getIndices();
     3453  updateBy = updates->denseVector();
     3454  addSequence = numberColumns;
     3455  int j; 
     3456  double * duals = model_->dualRowSolution();
     3457  for (j=0;j<number;j++) {
     3458    int iSequence = index[j];
     3459    double value = duals[iSequence];
     3460    value -= updateBy[j];
     3461    updateBy[j]=0.0;
     3462    duals[iSequence] = value;
     3463  }
     3464  double bestDj = tolerance;
     3465  int bestSequence=-1;
     3466
     3467  const double * cost = model_->costRegion(1);
     3468
     3469  int iPassR=0,iPassC=0;
     3470  // Setup two passes
     3471  // This biases towards picking row variables
     3472  // This probably should be fixed
     3473  int startR[4];
     3474  const int * which = infeasible_->getIndices();
     3475  int nSlacks = infeasible_->getNumElements();
     3476  startR[1]=nSlacks;
     3477  startR[2]=0;
     3478  double randomR = CoinDrand48();
     3479  double dstart = ((double) nSlacks) * randomR;
     3480  startR[0]=(int) dstart;
     3481  startR[3]=startR[0];
     3482  int startC[4];
     3483  startC[1]=numberColumns;
     3484  startC[2]=0;
     3485  double randomC = CoinDrand48();
     3486  dstart = ((double) numberColumns) * randomC;
     3487  startC[0]=(int) dstart;
     3488  startC[3]=startC[0];
     3489  reducedCost = model_->djRegion(1);
     3490  int sequenceOut = model_->sequenceOut();
     3491  double * duals2 = duals-numberColumns;
     3492  int chunk = min(1024,(numberColumns+nSlacks)/32);
     3493  if (model_->numberIterations()%1000==0) {
     3494    printf("%d wanted, nSlacks %d, chunk %d\n",numberWanted,nSlacks,chunk);
     3495    int i;
     3496    for (i=0;i<4;i++)
     3497      printf("start R %d C %d ",startR[i],startC[i]);
     3498    printf("\n");
     3499  }
     3500  chunk = max(chunk,256);
     3501  bool finishedR=false,finishedC=false;
     3502  bool doingR = randomR>randomC;
     3503  doingR=false;
     3504  while (!finishedR||!finishedC) {
     3505    if (finishedR)
     3506      doingR=false;
     3507    if (doingR) {
     3508      int saveSequence = bestSequence;
     3509      int start = startR[iPassR];
     3510      int end = min(startR[iPassR+1],start+chunk/2);
     3511      int jSequence;
     3512      for (jSequence=start;jSequence<end;jSequence++) {
     3513        int iSequence=which[jSequence];
     3514        if (iSequence!=sequenceOut) {
     3515          double value;
     3516          ClpSimplex::Status status = model_->getStatus(iSequence);
     3517         
     3518          switch(status) {
     3519           
     3520          case ClpSimplex::basic:
     3521          case ClpSimplex::isFixed:
     3522            break;
     3523          case ClpSimplex::isFree:
     3524          case ClpSimplex::superBasic:
     3525            value = fabs(cost[iSequence]+duals2[iSequence]);
     3526            if (value>FREE_ACCEPT*tolerance) {
     3527              numberWanted--;
     3528              // we are going to bias towards free (but only if reasonable)
     3529              value *= FREE_BIAS;
     3530              if (value>bestDj) {
     3531                // check flagged variable and correct dj
     3532                if (!model_->flagged(iSequence)) {
     3533                  bestDj=value;
     3534                  bestSequence = iSequence;
     3535                } else {
     3536                  // just to make sure we don't exit before got something
     3537                  numberWanted++;
     3538                }
     3539              }
     3540            }
     3541            break;
     3542          case ClpSimplex::atUpperBound:
     3543            value = cost[iSequence]+duals2[iSequence];
     3544            if (value>tolerance) {
     3545              numberWanted--;
     3546              if (value>bestDj) {
     3547                // check flagged variable and correct dj
     3548                if (!model_->flagged(iSequence)) {
     3549                  bestDj=value;
     3550                  bestSequence = iSequence;
     3551                } else {
     3552                  // just to make sure we don't exit before got something
     3553                  numberWanted++;
     3554                }
     3555              }
     3556            }
     3557            break;
     3558          case ClpSimplex::atLowerBound:
     3559            value = -(cost[iSequence]+duals2[iSequence]);
     3560            if (value>tolerance) {
     3561              numberWanted--;
     3562              if (value>bestDj) {
     3563                // check flagged variable and correct dj
     3564                if (!model_->flagged(iSequence)) {
     3565                  bestDj=value;
     3566                  bestSequence = iSequence;
     3567                } else {
     3568                  // just to make sure we don't exit before got something
     3569                  numberWanted++;
     3570                }
     3571              }
     3572            }
     3573            break;
     3574          }
     3575        }
     3576        if (!numberWanted)
     3577          break;
     3578      }
     3579      if (saveSequence!=bestSequence) {
     3580        // dj
     3581        reducedCost[bestSequence] = cost[bestSequence] +duals[bestSequence-numberColumns];
     3582        bestDj=reducedCost[bestSequence];
     3583      }
     3584      if (!numberWanted)
     3585        break;
     3586      doingR=false;
     3587      // update start
     3588      startR[iPassR]=jSequence;
     3589      if (jSequence>=startR[iPassR+1]) {
     3590        if (iPassR)
     3591          finishedR=true;
     3592        else
     3593          iPassR=2;
     3594      }
     3595    }
     3596    if (finishedC)
     3597      doingR=true;
     3598    if (!doingR) {
     3599      int saveSequence = bestSequence;
     3600      // Columns
     3601      int start = startC[iPassC];
     3602      int end = min(startC[iPassC+1],start+chunk);;
     3603      model_->clpMatrix()->partialPricing(model_,start,end,bestSequence,numberWanted);
     3604      if (saveSequence!=bestSequence) {
     3605        // dj
     3606        bestDj=reducedCost[bestSequence];
     3607      }
     3608      if (!numberWanted)
     3609        break;
     3610      doingR=true;
     3611      // update start
     3612      int iSequence = end;
     3613      startC[iPassC]=iSequence;
     3614      if (iSequence>=startC[iPassC+1]) {
     3615        if (iPassC)
     3616          finishedC=true;
     3617        else
     3618          iPassC=2;
     3619      }
     3620    }
     3621  }
     3622  updates->setNumElements(0);
     3623
     3624  // Restore tolerance
     3625  model_->setCurrentDualTolerance(saveTolerance);
     3626#ifndef NDEBUG
     3627  if (bestSequence>=0) {
     3628    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
     3629      assert(model_->reducedCost(bestSequence)<0.0);
     3630    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
     3631      assert(model_->reducedCost(bestSequence)>0.0);
     3632  }
     3633#endif
     3634  return bestSequence;
     3635}
  • trunk/ClpSimplexDual.cpp

    r256 r259  
    29612961  int maxLength=0;
    29622962  int minLength=numberRows_;
     2963  double averageCost = 0.0;
     2964  int numberNonZero=0;
    29632965  if (!numberIterations_&&perturbation_==50) {
    29642966    // See if we need to perturb
     
    29702972      double value = fabs(obj[i]);
    29712973      sort[i]=value;
    2972     }
     2974      averageCost += value;
     2975      if (value)
     2976        numberNonZero++;
     2977    }
     2978    if (numberNonZero)
     2979      averageCost /= (double) numberNonZero;
     2980    else
     2981      averageCost = 1.0;
    29732982    std::sort(sort,sort+numberColumns_);
    29742983    int number=1;
     
    30183027  int iRow;
    30193028  double smallestNonZero=1.0e100;
    3020   int numberNonZero=0;
     3029  numberNonZero=0;
    30213030  if (perturbation_>=50) {
    30223031    perturbation = 1.0e-8;
     
    30913100          smallestPositive==largestPositive) {
    30923101        // Really hit perturbation
    3093         maximumFraction = max(1.0e-3*max(lastValue,lastValue2),maximumFraction);
     3102        double adjust = min(100.0*maximumFraction,1.0e-3*max(lastValue,lastValue2));
     3103        maximumFraction = max(adjust,maximumFraction);
    30943104      }
    30953105    }
     
    31723182  // Make variables with more elements more expensive
    31733183  const double m1 = 0.5;
     3184  double smallestAllowed = min(1.0e-2*dualTolerance_,maximumFraction);
     3185  double largestAllowed = max(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
    31743186  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    31753187    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
     
    32193231        value *= multiplier;
    32203232        value = min (value,value2);
    3221         if (savePerturbation!=50) {
     3233        if (savePerturbation<50||savePerturbation>60) {
    32223234          if (fabs(value)<=dualTolerance_)
    32233235            value=0.0;
    32243236        } else if (value) {
    32253237          // get in range
    3226           if (fabs(value)<=dualTolerance_) {
     3238          if (fabs(value)<=smallestAllowed) {
    32273239            value *= 10.0;
    3228             while (fabs(value)<=dualTolerance_)
     3240            while (fabs(value)<=smallestAllowed)
    32293241              value *= 10.0;
    3230           } else if (fabs(value)>0.1) {
     3242          } else if (fabs(value)>largestAllowed) {
    32313243            value *= 0.1;
    3232             while (fabs(value)>0.1)
     3244            while (fabs(value)>largestAllowed)
    32333245              value *= 0.1;
    32343246          }
  • trunk/ClpSimplexPrimal.cpp

    r256 r259  
    191191  if (!startup(ifValuesPass)) {
    192192   
     193    // Set average theta
     194    nonLinearCost_->setAverageTheta(1.0e3);
    193195    int lastCleaned=0; // last time objective or bounds cleaned up
    194196   
     
    415417  // status -3 to go to top without an invert
    416418  while (problemStatus_==-1) {
     419    //#define CLP_DEBUG 1
    417420#ifdef CLP_DEBUG
    418421    {
     
    940943                            int valuesPass)
    941944{
     945  //rowArray->scanAndPack();
    942946  if (valuesPass) {
    943947    dualIn_ = cost_[sequenceIn_];
     
    951955     
    952956      int iRow = which[iIndex];
    953       double alpha = work[iRow];
     957      double alpha = work[iIndex];
    954958      int iPivot=pivotVariable_[iRow];
    955       dualIn_ -= alpha*cost(iPivot);
     959      dualIn_ -= alpha*cost_[iPivot];
    956960    }
    957961    // determine direction here
     
    974978  // sequence stays as row number until end
    975979  pivotRow_=-1;
    976   int numberSwapped=0;
    977980  int numberRemaining=0;
    978981
    979   int numberThru =0; // number gone thru a barrier
    980   int lastThru =0; // number gone thru a barrier on last time
    981  
    982982  double totalThru=0.0; // for when variables flip
    983983  double acceptablePivot=1.0e-7;
     
    988988  double lastPivot=0.0;
    989989  double lastTheta=1.0e50;
    990   int lastNumberSwapped=0;
    991990
    992991  // use spareArrays to put ones looked at in
     
    996995  // with a backup list saved in double * part of indexed vector
    997996
    998   // for zeroing out arrays after
    999   int maximumSwapped=0;
    1000997  // pivot elements
    1001998  double * spare;
    1002999  // indices
    1003   int * index, * indexSwapped;
    1004   int * saveSwapped;
     1000  int * index;
    10051001  spareArray->clear();
    1006   spareArray2->clear();
    10071002  spare = spareArray->denseVector();
    10081003  index = spareArray->getIndices();
    1009   saveSwapped = (int *) spareArray2->denseVector();
    1010   indexSwapped = spareArray2->getIndices();
    10111004
    10121005  // we also need somewhere for effective rhs
    10131006  double * rhs=rhsArray->denseVector();
     1007  //int * indexRhs = rhsArray->getIndices();
     1008  //int numberFlip=0; // Those which may change if flips
    10141009
    10151010  /*
     
    10251020  // do first pass to get possibles
    10261021  // We can also see if unbounded
    1027   // We also re-compute reduced cost
    1028 
    1029   dualIn_ = cost_[sequenceIn_];
    10301022
    10311023  double * work=rowArray->denseVector();
     
    10411033    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
    10421034
    1043   double tentativeTheta = maximumMovement;
     1035  double averageTheta = nonLinearCost_->averageTheta();
     1036  double tentativeTheta = min(10.0*averageTheta,maximumMovement);
    10441037  double upperTheta = maximumMovement;
     1038  if (tentativeTheta>0.5*maximumMovement)
     1039    tentativeTheta=maximumMovement;
     1040
     1041  double dualCheck = fabs(dualIn_);
     1042  // but make a bit more pessimistic
     1043  dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
    10451044
    10461045  int iIndex;
    1047 #if 0
    1048   if (numberIterations_<=39)
    1049     handler_->setLogLevel(63);
    1050   else
    1051     handler_->setLogLevel(2);
    1052   if (numberIterations_==38)
    1053     printf("trouble\n");
    1054   assert (solution_[29176]>-1.0e20);
    1055 #endif
    1056   for (iIndex=0;iIndex<number;iIndex++) {
    1057 
    1058     int iRow = which[iIndex];
    1059     double alpha = work[iRow];
    1060     int iPivot=pivotVariable_[iRow];
    1061     dualIn_ -= alpha*cost(iPivot);
    1062     alpha *= way;
    1063     double oldValue = solution(iPivot);
    1064     // get where in bound sequence
    1065     if (alpha>0.0) {
    1066       // basic variable going towards lower bound
    1067       double bound = lower(iPivot);
    1068       oldValue -= bound;
     1046  bool gotList=false;
     1047  int pivotOne=-1;
     1048  while (!gotList) {
     1049    pivotOne=-1;
     1050    totalThru=0.0;
     1051    // We also re-compute reduced cost
     1052    numberRemaining=0;
     1053    dualIn_ = cost_[sequenceIn_];
     1054    double tolerance = primalTolerance_*1.002;
     1055    for (iIndex=0;iIndex<number;iIndex++) {
     1056
     1057      int iRow = which[iIndex];
     1058      double alpha = work[iIndex];
     1059      int iPivot=pivotVariable_[iRow];
     1060      dualIn_ -= alpha*cost_[iPivot];
     1061      alpha *= way;
     1062      double oldValue = solution_[iPivot];
     1063      // get where in bound sequence
     1064      if (alpha>0.0) {
     1065        // basic variable going towards lower bound
     1066        double bound = lower_[iPivot];
     1067        oldValue -= bound;
     1068      } else {
     1069        // basic variable going towards upper bound
     1070        double bound = upper_[iPivot];
     1071        oldValue = bound-oldValue;
     1072      }
     1073     
     1074      double value = oldValue-tentativeTheta*fabs(alpha);
     1075      assert (oldValue>=-tolerance);
     1076      if (value<=tolerance) {
     1077        value=oldValue-upperTheta*fabs(alpha);
     1078        if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot) {
     1079          upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
     1080          pivotOne=numberRemaining;
     1081        }
     1082        // add to list
     1083        spare[numberRemaining]=alpha;
     1084        rhs[numberRemaining]=oldValue;
     1085        index[numberRemaining++]=iRow;
     1086        totalThru += fabs(alpha);
     1087        setActive(iRow);
     1088        //} else if (value<primalTolerance_*1.002) {
     1089        // May change if is a flip
     1090        //indexRhs[numberFlip++]=iRow;
     1091      }
     1092    }
     1093    if (upperTheta<maximumMovement&&totalThru*infeasibilityCost_>=1.0001*dualCheck) {
     1094      // Can pivot here
     1095      gotList=true;
     1096    } else if (tentativeTheta<maximumMovement) {
     1097      //printf("Going round with average theta of %g\n",averageTheta);
     1098      tentativeTheta=maximumMovement;
    10691099    } else {
    1070       // basic variable going towards upper bound
    1071       double bound = upper(iPivot);
    1072       oldValue = bound-oldValue;
    1073     }
    1074      
    1075     double value = oldValue-tentativeTheta*fabs(alpha);
    1076     assert (oldValue>=-primalTolerance_*1.002);
    1077     if (value<-primalTolerance_) {
    1078       // add to list
    1079       spare[numberRemaining]=alpha;
    1080       rhs[iRow]=oldValue;
    1081       index[numberRemaining++]=iRow;
    1082       double value=oldValue-upperTheta*fabs(alpha);
    1083       if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
    1084         upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
    1085     }
    1086   }
    1087 #if 0
    1088   if (numberIterations_>17701)
    1089     handler_->setLogLevel(63);
    1090   if (!valuesPass&&fabs(dualIn_-saveDj)>1.0e-1*(1.0+fabs(saveDj))) {
    1091     double d=0.0;
    1092     for (iIndex=0;iIndex<number;iIndex++) {
    1093      
    1094       int iRow = which[iIndex];
    1095       double alpha = work[iRow];
    1096       int iPivot=pivotVariable_[iRow];
    1097       double value = alpha*cost(iPivot);
    1098       d -= value;
    1099       if (value>1.0e7)
    1100         printf("%d %g\n",iRow,value);
    1101     }
    1102   }
    1103 #endif
     1100      gotList=true;
     1101    }
     1102  }
     1103  totalThru=0.0;
    11041104  // we need to keep where rhs non-zeros are
    11051105  int numberInRhs=numberRemaining;
    11061106  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
    11071107  rhsArray->setNumElements(numberInRhs);
     1108  rhsArray->setPacked();
    11081109
    11091110  theta_=maximumMovement;
    11101111
    1111   double dualCheck = fabs(dualIn_);
    1112   // but make a bit more pessimistic
    1113   dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
    1114 
    11151112  bool goBackOne = false;
     1113
     1114  //printf("%d remain out of %d\n",numberRemaining,number);
    11161115  int iTry=0;
    1117   int ibigTry=0;
    1118   if (numberRemaining) {
    1119 
    1120     // looks like pivoting
    1121     // now try until reasonable theta
    1122     tentativeTheta = max(10.0*upperTheta,1.0e-7);
    1123     tentativeTheta = min(tentativeTheta,maximumMovement);
    1124    
    1125     // loops increasing tentative theta until can't go through
    1126    
    1127     while (tentativeTheta <= maximumMovement) {
    1128       ibigTry++;
    1129       double thruThis = 0.0;
     1116#define MAXTRY 1000
     1117  if (numberRemaining&&upperTheta<maximumMovement) {
     1118    // First check if previously chosen one will work
     1119    if (pivotOne>=0&&0) {
     1120      double thruCost = infeasibilityCost_*fabs(spare[pivotOne]);
     1121      if (thruCost>=0.99*fabs(dualIn_))
     1122        printf("Could pivot on %d as change %g dj %g\n",
     1123               index[pivotOne],thruCost,dualIn_);
     1124      double alpha = spare[pivotOne];
     1125      double oldValue = rhs[pivotOne];
     1126      theta_ = oldValue/fabs(alpha);
     1127      pivotRow_=pivotOne;
     1128      // Stop loop
     1129      iTry=MAXTRY;
     1130    }
     1131
     1132    // first get ratio with tolerance
     1133    for ( ;iTry<MAXTRY;iTry++) {
    11301134     
     1135      upperTheta=maximumMovement;
     1136      int iBest=-1;
     1137      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
     1138       
     1139        double alpha = fabs(spare[iIndex]);
     1140        double oldValue = rhs[iIndex];
     1141        double value = oldValue-upperTheta*alpha;
     1142       
     1143        if (value<-primalTolerance_ && alpha>=acceptablePivot) {
     1144          upperTheta = (oldValue+primalTolerance_)/alpha;
     1145          iBest=iIndex; // just in case weird numbers
     1146        }
     1147      }
     1148     
     1149      // now look at best in this lot
    11311150      double bestPivot=acceptablePivot;
    1132       pivotRow_ = -1;
    1133      
    1134       numberSwapped = 0;
    1135      
    1136       upperTheta = maximumMovement;
    1137      
     1151      pivotRow_=-1;
    11381152      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
    1139 
     1153       
    11401154        int iRow = index[iIndex];
    11411155        double alpha = spare[iIndex];
    1142         double oldValue = rhs[iRow];
    1143         double value = oldValue-tentativeTheta*fabs(alpha);
    1144 
    1145         if (value<-primalTolerance_) {
    1146           // how much would it cost to go thru
    1147           thruThis += alpha*
    1148             nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
    1149           // goes on swapped list (also means candidates if too many)
    1150           indexSwapped[numberSwapped++]=iRow;
     1156        double oldValue = rhs[iIndex];
     1157        double value = oldValue-upperTheta*fabs(alpha);
     1158       
     1159        if (value<=0||iBest==iIndex) {
     1160          // how much would it cost to go thru and modify bound
     1161          totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha,rhs[iIndex]);
     1162          setActive(iRow);
    11511163          if (fabs(alpha)>bestPivot) {
    11521164            bestPivot=fabs(alpha);
    1153             pivotRow_ = iRow;
    11541165            theta_ = oldValue/bestPivot;
    1155           }
     1166            pivotRow_=iIndex;
     1167          }
     1168        }
     1169      }
     1170      if (bestPivot<0.1*bestEverPivot&&
     1171          bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
     1172        // back to previous one
     1173        goBackOne = true;
     1174        break;
     1175      } else if (pivotRow_==-1&&upperTheta>largeValue_) {
     1176        if (lastPivot>acceptablePivot) {
     1177          // back to previous one
     1178          goBackOne = true;
    11561179        } else {
    1157           value = oldValue-upperTheta*fabs(alpha);
    1158           if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
    1159             upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
    1160         }
    1161       }
    1162      
    1163       maximumSwapped = max(maximumSwapped,numberSwapped);
    1164 
    1165       if (totalThru+thruThis>=dualCheck) {
    1166         // We should be pivoting in this batch
    1167         // so compress down to this lot
    1168 
    1169         int saveNumber = numberRemaining;
    1170         numberRemaining=0;
    1171         for (iIndex=0;iIndex<numberSwapped;iIndex++) {
    1172           int iRow = indexSwapped[iIndex];
    1173           spare[numberRemaining]=way*work[iRow];
    1174           index[numberRemaining++]=iRow;
    1175         }
    1176         memset(spare+numberRemaining,0,
    1177                (saveNumber-numberRemaining)*sizeof(double));
    1178         //int iTry;
    1179 #define MAXTRY 100
    1180         // first get ratio with tolerance
    1181         for (iTry=0;iTry<MAXTRY;iTry++) {
    1182          
    1183           upperTheta=maximumMovement;
    1184           numberSwapped = 0;
    1185           int iBest=-1;
    1186           for (iIndex=0;iIndex<numberRemaining;iIndex++) {
    1187            
    1188             int iRow = index[iIndex];
    1189             double alpha = fabs(spare[iIndex]);
    1190             double oldValue = rhs[iRow];
    1191             double value = oldValue-upperTheta*alpha;
    1192            
    1193             if (value<-primalTolerance_ && alpha>=acceptablePivot) {
    1194               upperTheta = (oldValue+primalTolerance_)/alpha;
    1195               iBest=iRow; // just in case weird numbers
    1196             }
    1197           }
    1198          
    1199           // now look at best in this lot
    1200           bestPivot=acceptablePivot;
    1201           pivotRow_=-1;
    1202           for (iIndex=0;iIndex<numberRemaining;iIndex++) {
    1203            
    1204             int iRow = index[iIndex];
    1205             double alpha = spare[iIndex];
    1206             double oldValue = rhs[iRow];
    1207             double value = oldValue-upperTheta*fabs(alpha);
    1208            
    1209             if (value<=0||iRow==iBest) {
    1210               // how much would it cost to go thru
    1211               totalThru += alpha*
    1212                 nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
    1213               // goes on swapped list (also means candidates if too many)
    1214               indexSwapped[numberSwapped++]=iRow;
    1215               if (fabs(alpha)>bestPivot) {
    1216                 bestPivot=fabs(alpha);
    1217                 theta_ = oldValue/bestPivot;
    1218                 pivotRow_=iRow;
    1219               }
    1220             } else {
    1221               value = oldValue-upperTheta*fabs(alpha);
    1222               if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
    1223                 upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
    1224             }
    1225           }
    1226          
    1227           maximumSwapped = max(maximumSwapped,numberSwapped);
    1228           // had (&&bestPivot<1.0e-3||totalThru>0.1*dualCheck) as well
    1229           if (bestPivot<0.1*bestEverPivot&&
    1230               bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
    1231             // back to previous one
    1232             goBackOne = true;
    1233             break;
    1234           } else if (pivotRow_==-1&&upperTheta>largeValue_) {
    1235             if (lastPivot>acceptablePivot) {
    1236               // back to previous one
    1237               goBackOne = true;
    1238               //break;
    1239             } else {
    1240               // can only get here if all pivots so far too small
    1241             }
    1242             break;
    1243           } else if (totalThru>=dualCheck) {
    1244             break; // no point trying another loop
    1245           } else {
    1246             // skip this lot
    1247             nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
    1248             lastPivotRow=pivotRow_;
    1249             lastTheta = theta_;
    1250             lastThru = numberThru;
    1251             numberThru += numberSwapped;
    1252             lastNumberSwapped = numberSwapped;
    1253             memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
    1254             if (bestPivot>bestEverPivot)
    1255               bestEverPivot=bestPivot;
    1256           }
     1180          // can only get here if all pivots so far too small
    12571181        }
    12581182        break;
     1183      } else if (totalThru>=dualCheck) {
     1184        break; // no point trying another loop
    12591185      } else {
    1260         // skip this lot
    1261         nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
    12621186        lastPivotRow=pivotRow_;
    12631187        lastTheta = theta_;
    1264         lastThru = numberThru;
    1265         numberThru += numberSwapped;
    1266         lastNumberSwapped = numberSwapped;
    1267         memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
    12681188        if (bestPivot>bestEverPivot)
    12691189          bestEverPivot=bestPivot;
    1270         totalThru += thruThis;
    1271         tentativeTheta = 2.0*upperTheta;
    12721190      }
    12731191    }
     
    12771195      pivotRow_=lastPivotRow;
    12781196      theta_ = lastTheta;
    1279             // undo this lot
    1280       nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
    1281       memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
    1282       numberSwapped = lastNumberSwapped;
    1283     }
    1284   }
    1285 
     1197    }
     1198  } else {
     1199    // mark ones which may move
     1200    //for (int i=0;i<numberFlip;i++) {
     1201    //int iRow= indexRhs[i];
     1202    //setActive(iRow);
     1203    //}
     1204  }
     1205  //if (iTry>50)
     1206  //printf("** %d tries\n",iTry);
    12861207  if (pivotRow_>=0) {
    1287    
    1288     alpha_ = work[pivotRow_];
     1208    int position=pivotRow_; // position in list
     1209    pivotRow_=index[position];
     1210    //alpha_ = work[pivotRow_];
     1211    alpha_ = way*spare[position];
    12891212    // translate to sequence
    12901213    sequenceOut_ = pivotVariable_[pivotRow_];
     
    13081231    if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
    13091232      theta_=minimumTheta;
    1310       for (iIndex=0;iIndex<numberSwapped;iIndex++) {
    1311         int iRow = indexSwapped[iIndex];
     1233      for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) {
     1234        largestInfeasibility = max (largestInfeasibility,
     1235                                    -(rhs[iIndex]-fabs(spare[iIndex])*theta_));
     1236      }
     1237#define CLP_DEBUG
    13121238#ifdef CLP_DEBUG
    1313         if (iRow==pivotRow_)
    1314           found=true;
    1315 #endif
    1316         largestInfeasibility = max (largestInfeasibility,
    1317                                     -(rhs[iRow]-fabs(work[iRow])*theta_));
    1318       }
    1319 #ifdef CLP_DEBUG
    1320       assert(found);
    13211239      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
    13221240        printf("Primal tolerance increased from %g to %g\n",
     
    13261244      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
    13271245    }
    1328 
     1246    // Need to look at all in some cases
     1247    if (theta_>tentativeTheta) {
     1248      for (iIndex=0;iIndex<number;iIndex++)
     1249        setActive(which[iIndex]);
     1250    }
    13291251    if (way<0.0)
    13301252      theta_ = - theta_;
     
    13651287  }
    13661288
     1289  double theta1 = max(theta_,1.0e-12);
     1290  double theta2 = numberIterations_*nonLinearCost_->averageTheta();
     1291  // Set average theta
     1292  nonLinearCost_->setAverageTheta((theta1+theta2)/((double) (numberIterations_+1)));
     1293  //if (numberIterations_%1000==0)
     1294  //printf("average theta is %g\n",nonLinearCost_->averageTheta());
     1295
    13671296  // clear arrays
    13681297
    13691298  memset(spare,0,numberRemaining*sizeof(double));
    1370   memset(saveSwapped,0,maximumSwapped*sizeof(int));
    13711299
    13721300  // put back original bounds etc
     
    14451373int
    14461374ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
    1447                   double theta,
    1448                   double & objectiveChange,
     1375                                        double theta,
     1376                                        double & objectiveChange,
    14491377                                        int valuesPass)
    14501378{
     1379  // Cost on pivot row may change - may need to change dualIn
     1380  double oldCost=0.0;
     1381  if (pivotRow_>=0)
     1382    oldCost = cost_[sequenceOut_];
     1383  //rowArray->scanAndPack();
    14511384  double * work=rowArray->denseVector();
    14521385  int number=rowArray->getNumElements();
     
    14541387
    14551388  int newNumber = 0;
    1456 
     1389  int pivotPosition = -1;
    14571390  nonLinearCost_->setChangeInCost(0.0);
    14581391  int iIndex;
    1459 
    1460   for (iIndex=0;iIndex<number;iIndex++) {
    1461 
    1462     int iRow = which[iIndex];
    1463     double alpha = work[iRow];
    1464     int iPivot=pivotVariable_[iRow];
    1465     double & value = solutionAddress(iPivot);
    1466     double change = theta*alpha;
    1467     value -= change;
    1468     // But make sure one going out is feasible
    1469     if (change>0.0) {
    1470       // going down
    1471       if (value<lower(iPivot)+primalTolerance_) {
    1472         if (iPivot==sequenceOut_&&value>lower(iPivot)-1.001*primalTolerance_)
    1473           value=lower(iPivot);
    1474         double difference = nonLinearCost_->setOne(iPivot,value);
    1475         work[iRow] = difference;
    1476         if (difference) {
    1477           //change reduced cost on this
    1478           reducedCostAddress(iPivot) = -difference;
    1479           which[newNumber++]=iRow;
     1392  if (!valuesPass) {
     1393    for (iIndex=0;iIndex<number;iIndex++) {
     1394     
     1395      int iRow = which[iIndex];
     1396      double alpha = work[iIndex];
     1397      work[iIndex]=0.0;
     1398      int iPivot=pivotVariable_[iRow];
     1399      double change = theta*alpha;
     1400      double value = solution_[iPivot] - change;
     1401      solution_[iPivot]=value;
     1402#ifndef NDEBUG
     1403      // check if not active then okay
     1404      if (!active(iRow)) {
     1405        // But make sure one going out is feasible
     1406        if (change>0.0) {
     1407          // going down
     1408          if (value<lower_[iPivot]+primalTolerance_) {
     1409            if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
     1410              value=lower_[iPivot];
     1411            double difference = nonLinearCost_->setOne(iPivot,value);
     1412            assert (!difference);
     1413          }
     1414        } else {
     1415          // going up
     1416          if (value>upper_[iPivot]-primalTolerance_) {
     1417            if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
     1418              value=upper_[iPivot];
     1419            double difference = nonLinearCost_->setOne(iPivot,value);
     1420            assert (!difference);
     1421          }
     1422        }
     1423      }
     1424#endif   
     1425      if (active(iRow)) {
     1426        clearActive(iRow);
     1427        // But make sure one going out is feasible
     1428        if (change>0.0) {
     1429          // going down
     1430          if (value<lower_[iPivot]+primalTolerance_) {
     1431            if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
     1432              value=lower_[iPivot];
     1433            double difference = nonLinearCost_->setOne(iPivot,value);
     1434            if (difference) {
     1435              if (iRow==pivotRow_)
     1436                pivotPosition=newNumber;
     1437              work[newNumber] = difference;
     1438              //change reduced cost on this
     1439              dj_[iPivot] = -difference;
     1440              which[newNumber++]=iRow;
     1441            }
     1442          }
     1443        } else {
     1444          // going up
     1445          if (value>upper_[iPivot]-primalTolerance_) {
     1446            if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
     1447              value=upper_[iPivot];
     1448            double difference = nonLinearCost_->setOne(iPivot,value);
     1449            if (difference) {
     1450              if (iRow==pivotRow_)
     1451                pivotPosition=newNumber;
     1452              work[newNumber] = difference;
     1453              //change reduced cost on this
     1454              dj_[iPivot] = -difference;
     1455              which[newNumber++]=iRow;
     1456            }
     1457          }
     1458        }
     1459      }
     1460    }
     1461  } else {
     1462    // values pass so look at all
     1463    for (iIndex=0;iIndex<number;iIndex++) {
     1464     
     1465      int iRow = which[iIndex];
     1466      double alpha = work[iIndex];
     1467      work[iIndex]=0.0;
     1468      int iPivot=pivotVariable_[iRow];
     1469      double change = theta*alpha;
     1470      double value = solution_[iPivot] - change;
     1471      solution_[iPivot]=value;
     1472      clearActive(iRow);
     1473      // But make sure one going out is feasible
     1474      if (change>0.0) {
     1475        // going down
     1476        if (value<lower_[iPivot]+primalTolerance_) {
     1477          if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
     1478            value=lower_[iPivot];
     1479          double difference = nonLinearCost_->setOne(iPivot,value);
     1480          if (difference) {
     1481            if (iRow==pivotRow_)
     1482              pivotPosition=newNumber;
     1483            work[newNumber] = difference;
     1484            //change reduced cost on this
     1485            dj_[iPivot] = -difference;
     1486            which[newNumber++]=iRow;
     1487          }
    14801488        }
    14811489      } else {
    1482         work[iRow]=0.0;
    1483       }
     1490        // going up
     1491        if (value>upper_[iPivot]-primalTolerance_) {
     1492          if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
     1493            value=upper_[iPivot];
     1494          double difference = nonLinearCost_->setOne(iPivot,value);
     1495          if (difference) {
     1496            if (iRow==pivotRow_)
     1497              pivotPosition=newNumber;
     1498            work[newNumber] = difference;
     1499            //change reduced cost on this
     1500            dj_[iPivot] = -difference;
     1501            which[newNumber++]=iRow;
     1502          }
     1503        }
     1504      }
     1505    }
     1506  }
     1507  objectiveChange += nonLinearCost_->changeInCost();
     1508  rowArray->setPacked();
     1509#if 0
     1510  rowArray->setNumElements(newNumber);
     1511  rowArray->expand();
     1512  if (pivotRow_>=0) {
     1513    dualIn_ += (oldCost-cost_[sequenceOut_]);
     1514    // update change vector to include pivot
     1515    rowArray->add(pivotRow_,-dualIn_);
     1516    // and convert to packed
     1517    rowArray->scanAndPack();
     1518  } else {
     1519    // and convert to packed
     1520    rowArray->scanAndPack();
     1521  }
     1522#else
     1523  if (pivotRow_>=0) {
     1524    dualIn_ += (oldCost-cost_[sequenceOut_]);
     1525    // update change vector to include pivot
     1526    if (pivotPosition>=0) {
     1527      work[pivotPosition] -= dualIn_;
    14841528    } else {
    1485       // going up
    1486       if (value>upper(iPivot)-primalTolerance_) {
    1487         if (iPivot==sequenceOut_&&value<upper(iPivot)+1.001*primalTolerance_)
    1488           value=upper(iPivot);
    1489         double difference = nonLinearCost_->setOne(iPivot,value);
    1490         work[iRow] = difference;
    1491         if (difference) {
    1492           //change reduced cost on this
    1493           reducedCostAddress(iPivot) = -difference;
    1494           which[newNumber++]=iRow;
    1495         }
    1496       } else {
    1497         work[iRow]=0.0;
    1498       }
    1499     }
    1500   }
    1501   objectiveChange += nonLinearCost_->changeInCost();
     1529      work[newNumber]=-dualIn_;
     1530      which[newNumber++]=pivotRow_;
     1531    }
     1532  }
    15021533  rowArray->setNumElements(newNumber);
    1503 #if 0
    1504   if (newNumber>10) {
    1505     printf("in %d out %d, row %d alpha %g\n",
    1506            sequenceIn_,sequenceOut_,pivotRow_,alpha_);
    1507     for (iIndex=0;iIndex<newNumber;iIndex++) {
    1508       int iRow = which[iIndex];
    1509       double alpha = work[iRow];
    1510       printf("%d %g\n",iRow,alpha);
    1511     }
    1512   }
    15131534#endif
    1514 
    15151535  return 0;
    15161536}
     
    15231543  if (perturbation_==100)
    15241544    perturbation_=50; // treat as normal
     1545  int savePerturbation = perturbation_;
    15251546  int i;
    15261547  if (!numberIterations_)
     
    15621583  int numberNonZero=0;
    15631584  // maximum fraction of rhs/bounds to perturb
    1564   double maximumFraction = 1.0e-7;
     1585  double maximumFraction = 1.0e-5;
    15651586  if (perturbation_>=50) {
    15661587    perturbation = 1.0e-4;
     
    16141635  // Change so at least 1.0e-5 and no more than 0.1
    16151636  // For now just no more than 0.1
     1637  // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
    16161638  if (type==1) {
    1617     double multiplier = perturbation*maximumFraction;
     1639    //double multiplier = perturbation*maximumFraction;
    16181640    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    16191641      if (getStatus(iSequence)==basic) {
     
    16241646        difference = min(difference,perturbation);
    16251647        difference = min(difference,fabs(solutionValue)+1.0);
    1626         double value = multiplier*(difference+1.0);
     1648        double value = maximumFraction*(difference+1.0);
    16271649        value = min(value,0.1);
    16281650        value *= CoinDrand48();
     
    16611683    }
    16621684  } else {
     1685    double tolerance = 100.0*primalTolerance_;
    16631686    for (i=0;i<numberColumns_;i++) {
    16641687      double lowerValue=lower_[i], upperValue=upper_[i];
     
    16671690        value = min(value,0.1);
    16681691        value *= CoinDrand48();
    1669         if (lowerValue>-1.0e20&&lowerValue)
    1670           lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1671         if (upperValue<1.0e20&&upperValue)
    1672           upperValue += value * (max(1.0,1.0e-5*fabs(upperValue)));
     1692        if (savePerturbation!=50) {
     1693          if (fabs(value)<=primalTolerance_)
     1694            value=0.0;
     1695          if (lowerValue>-1.0e20&&lowerValue)
     1696            lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
     1697          if (upperValue<1.0e20&&upperValue)
     1698            upperValue += value * (max(1.0e-2,1.0e-5*fabs(upperValue)));
     1699        } else if (value) {
     1700          double valueL =value *(max(1.0e-2,1.0e-5*fabs(lowerValue)));
     1701          // get in range
     1702          if (valueL<=tolerance) {
     1703            valueL *= 10.0;
     1704            while (valueL<=tolerance)
     1705              valueL *= 10.0;
     1706          } else if (valueL>1.0) {
     1707            valueL *= 0.1;
     1708            while (valueL>1.0)
     1709              valueL *= 0.1;
     1710          }
     1711          if (lowerValue>-1.0e20&&lowerValue)
     1712            lowerValue -= valueL;
     1713          double valueU =value *(max(1.0e-2,1.0e-5*fabs(upperValue)));
     1714          // get in range
     1715          if (valueU<=tolerance) {
     1716            valueU *= 10.0;
     1717            while (valueU<=tolerance)
     1718              valueU *= 10.0;
     1719          } else if (valueU>1.0) {
     1720            valueU *= 0.1;
     1721            while (valueU>1.0)
     1722              valueU *= 0.1;
     1723          }
     1724          if (upperValue<1.0e20&&upperValue)
     1725            upperValue += valueU;
     1726        }
    16731727        if (lowerValue!=lower_[i]) {
    16741728          double difference = fabs(lowerValue-lower_[i]);
     
    16951749      value = min(value,0.1);
    16961750      value *= CoinDrand48();
    1697       if (upperValue>lowerValue+primalTolerance_) {
    1698         if (lowerValue>-1.0e20&&lowerValue)
    1699           lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1700         if (upperValue<1.0e20&&upperValue)
    1701           upperValue += value * (max(1.0,1.0e-5*fabs(upperValue)));
     1751      if (upperValue>lowerValue+tolerance) {
     1752        if (savePerturbation!=50) {
     1753          if (fabs(value)<=primalTolerance_)
     1754            value=0.0;
     1755          if (lowerValue>-1.0e20&&lowerValue)
     1756            lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
     1757          if (upperValue<1.0e20&&upperValue)
     1758            upperValue += value * (max(1.0e-2,1.0e-5*fabs(upperValue)));
     1759        } else if (value) {
     1760          double valueL =value *(max(1.0e-2,1.0e-5*fabs(lowerValue)));
     1761          // get in range
     1762          if (valueL<=tolerance) {
     1763            valueL *= 10.0;
     1764            while (valueL<=tolerance)
     1765              valueL *= 10.0;
     1766          } else if (valueL>1.0) {
     1767            valueL *= 0.1;
     1768            while (valueL>1.0)
     1769              valueL *= 0.1;
     1770          }
     1771          if (lowerValue>-1.0e20&&lowerValue)
     1772            lowerValue -= valueL;
     1773          double valueU =value *(max(1.0e-2,1.0e-5*fabs(upperValue)));
     1774          // get in range
     1775          if (valueU<=tolerance) {
     1776            valueU *= 10.0;
     1777            while (valueU<=tolerance)
     1778              valueU *= 10.0;
     1779          } else if (valueU>1.0) {
     1780            valueU *= 0.1;
     1781            while (valueU>1.0)
     1782              valueU *= 0.1;
     1783          }
     1784          if (upperValue<1.0e20&&upperValue)
     1785            upperValue += valueU;
     1786        }
    17021787      } else if (upperValue>0.0) {
    1703         upperValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1704         lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1788        upperValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
     1789        lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
    17051790      } else if (upperValue<0.0) {
    1706         upperValue += value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1707         lowerValue += value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1791        upperValue += value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
     1792        lowerValue += value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
    17081793      } else {
    17091794      }
     
    18381923  bool roundAgain=true;
    18391924  int returnCode=-1;
     1925  double yyyy[20000];
     1926  if (numberIterations_==5141)
     1927    memcpy(yyyy,solution_,(numberColumns_+numberRows_)*sizeof(double));
    18401928
    18411929  // loop round if user setting and doing refactorization
     
    18461934    sequenceOut_=-1;
    18471935    rowArray_[1]->clear();
     1936#if 0
     1937    {
     1938      int seq[]={612,643};
     1939      int k;
     1940      for (k=0;k<sizeof(seq)/sizeof(int);k++) {
     1941        int iSeq=seq[k];
     1942        if (getColumnStatus(iSeq)!=basic) {
     1943          double djval;
     1944          double * work;
     1945          int number;
     1946          int * which;
     1947         
     1948          int iIndex;
     1949          unpack(rowArray_[1],iSeq);
     1950          factorization_->updateColumn(rowArray_[2],rowArray_[1]);
     1951          djval = cost_[iSeq];
     1952          work=rowArray_[1]->denseVector();
     1953          number=rowArray_[1]->getNumElements();
     1954          which=rowArray_[1]->getIndices();
     1955         
     1956          for (iIndex=0;iIndex<number;iIndex++) {
     1957           
     1958            int iRow = which[iIndex];
     1959            double alpha = work[iRow];
     1960            int iPivot=pivotVariable_[iRow];
     1961            djval -= alpha*cost_[iPivot];
     1962          }
     1963          double comp = 1.0e-8 + 1.0e-7*(max(fabs(dj_[iSeq]),fabs(djval)));
     1964          if (fabs(djval-dj_[iSeq])>comp)
     1965            printf("Bad dj %g for %d - true is %g\n",
     1966                   dj_[iSeq],iSeq,djval);
     1967          assert (fabs(djval)<1.0e-3||djval*dj_[iSeq]>0.0);
     1968          rowArray_[1]->clear();
     1969        }
     1970      }
     1971    }
     1972#endif
     1973       
    18481974    // we found a pivot column
    18491975    // update the incoming column
    1850     unpack(rowArray_[1]);
     1976    unpackPacked(rowArray_[1]);
    18511977    // save reduced cost
    18521978    double saveDj = dualIn_;
     
    18742000            progress_->clearBadTimes();
    18752001            lastBadIteration_ = numberIterations_; // say be more cautious
    1876             rowArray_[1]->clear();
     2002            clearAll();
    18772003            pivotRow_=-1;
    18782004          }
     
    18912017        <<CoinMessageEol;
    18922018      if(lastGoodIteration_ != numberIterations_) {
    1893         rowArray_[1]->clear();
     2019        clearAll();
    18942020        pivotRow_=-1; // say no weights update
    18952021        returnCode=-4;
     
    19132039          progress_->clearBadTimes();
    19142040          lastBadIteration_ = numberIterations_; // say be more cautious
    1915           rowArray_[1]->clear();
     2041          clearAll();
    19162042          pivotRow_=-1;
    19172043          returnCode=-5;
     
    19282054        // update duals
    19292055        if (pivotRow_>=0) {
    1930           alpha_ = rowArray_[1]->denseVector()[pivotRow_];
     2056          // as packed need to find pivot row
     2057          //assert (rowArray_[1]->packedMode());
     2058          //int i;
     2059         
     2060          //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
    19312061          assert (fabs(alpha_)>1.0e-8);
    19322062          double multiplier = dualIn_/alpha_;
     
    19872117        // later we may need to unwind more e.g. fake bounds
    19882118        if(lastGoodIteration_ != numberIterations_) {
    1989           rowArray_[1]->clear();
     2119          clearAll();
    19902120          pivotRow_=-1;
    19912121          if (solveType_==1) {
     
    20142144          }
    20152145          lastBadIteration_ = numberIterations_; // say be more cautious
    2016           rowArray_[1]->clear();
     2146          clearAll();
    20172147          pivotRow_=-1;
    20182148          sequenceOut_=-1;
     
    20672197   
    20682198    double objectiveChange=0.0;
    2069     // Cost on pivot row may change - may need to change dualIn
    2070     double oldCost=0.0;
    2071     if (pivotRow_>=0)
    2072       oldCost = cost(pivotVariable_[pivotRow_]);
    2073     // rowArray_[1] is not empty - used to update djs
     2199    // after this rowArray_[1] is not empty - used to update djs
    20742200    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
    2075     if (pivotRow_>=0)
    2076       dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
     2201   
    20772202    double oldValue = valueIn_;
    20782203    if (directionIn_==-1) {
     
    21012226        valueOut_ = upperOut_;
    21022227      }
    2103       assert(valueOut_>=lower_[sequenceOut_]-primalTolerance_&&
    2104              valueOut_<=upper_[sequenceOut_]+primalTolerance_);
     2228      if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
     2229        valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
     2230      else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
     2231        valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
    21052232      // may not be exactly at bound and bounds may have changed
    21062233      // Make sure outgoing looks feasible
     
    21112238    nonLinearCost_->setOne(sequenceIn_,valueIn_);
    21122239    int whatNext=housekeeping(objectiveChange);
     2240
    21132241#if 0
    21142242    if (numberIterations_==1148)
     
    21722300  if (sequenceIn_<numberColumns_)
    21732301    ray_[sequenceIn_]=directionIn_;
    2174   for (i=0;i<number;i++) {
    2175     int iRow=index[i];
    2176     int iPivot=pivotVariable_[iRow];
    2177     double arrayValue = array[iRow];
    2178     if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
    2179       ray_[iPivot] = way* array[iRow];
     2302  if (!rowArray->packedMode()) {
     2303    for (i=0;i<number;i++) {
     2304      int iRow=index[i];
     2305      int iPivot=pivotVariable_[iRow];
     2306      double arrayValue = array[iRow];
     2307      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
     2308        ray_[iPivot] = way* arrayValue;
     2309    }
     2310  } else {
     2311    for (i=0;i<number;i++) {
     2312      int iRow=index[i];
     2313      int iPivot=pivotVariable_[iRow];
     2314      double arrayValue = array[i];
     2315      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
     2316        ray_[iPivot] = way* arrayValue;
     2317    }
    21802318  }
    21812319}
     
    22612399  }
    22622400}
    2263  
     2401void
     2402ClpSimplexPrimal::clearAll()
     2403{
     2404  int number=rowArray_[1]->getNumElements();
     2405  int * which=rowArray_[1]->getIndices();
    22642406 
    2265 
     2407  int iIndex;
     2408  for (iIndex=0;iIndex<number;iIndex++) {
     2409   
     2410    int iRow = which[iIndex];
     2411    clearActive(iRow);
     2412  }
     2413  rowArray_[1]->clear();
     2414}
     2415
  • trunk/ClpSolve.cpp

    r256 r259  
    508508          info.setStartingWeight(1.0e3);
    509509          info.setReduceIterations(6);
     510          // also be more lenient on infeasibilities
     511          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
     512          info.setDropEnoughWeighted(-2.0);
    510513        } else if (nPasses>=50) {
    511514          info.setStartingWeight(1.0e3);
  • trunk/Test/ClpMain.cpp

    r248 r259  
    11// Copyright (C) 2002, International Business Machines
    22// Corporation and others.  All Rights Reserved.
    3 
     3   
    44#include "CoinPragma.hpp"
    55
     
    1414#include "CoinPragma.hpp"
    1515#include "CoinHelperFunctions.hpp"
    16 #define CLPVERSION "0.98.6"
     16#define CLPVERSION "0.99.0"
    1717
    1818#include "CoinMpsIO.hpp"
  • trunk/include/Idiot.hpp

    r256 r259  
    158158  inline void setStrategy(int value)
    159159  { strategy_ = value;};
     160  /// Fine tuning - okay if feasibility drop this factor
     161  inline double getDropEnoughFeasibility() const
     162  { return dropEnoughFeasibility_;};
     163  inline void setDropEnoughFeasibility(double value)
     164  { dropEnoughFeasibility_=value;};
     165  /// Fine tuning - okay if weighted obj drop this factor
     166  inline double getDropEnoughWeighted() const
     167  { return dropEnoughWeighted_;};
     168  inline void setDropEnoughWeighted(double value)
     169  { dropEnoughWeighted_=value;};
    160170  //@}
    161171
Note: See TracChangeset for help on using the changeset viewer.