Changeset 336


Ignore:
Timestamp:
Mar 22, 2004 11:30:03 AM (17 years ago)
Author:
forrest
Message:

CPutting back changesVS: ----------------------------------------------------------------------

Location:
trunk
Files:
34 edited

Legend:

Unmodified
Added
Removed

  • trunk/ClpDualRowSteepest.cpp

    r266 r336  
    99#include "CoinHelperFunctions.hpp"
    1010#include <cstdio>
    11 
    1211//#############################################################################
    1312// Constructors / Destructor / Assignment

  • trunk/ClpDynamicExampleMatrix.cpp

    r335 r336  
    352352    for (int i=0;i<numberIds;i++) {
    353353      int sequence = ids[i];
     354      int iSet = set[sequence];
    354355      CoinBigIndex start = startColumnGen_[sequence];
    355356      addColumn(startColumnGen_[sequence+1]-start,
    … …  
    359360                columnLowerGen_ ? columnLowerGen_[sequence] : 0,
    360361                columnUpperGen_ ? columnUpperGen_[sequence] : 1.0e30,
    361                 set[sequence],getDynamicStatus(i));
    362       idGen_[iSet]=sequence; // say which one in
     362                iSet,getDynamicStatus(i));
     363      idGen_[i]=sequence; // say which one in
    363364      setDynamicStatusGen(sequence,inSmall);
    364365    }
    … …  
    427428// Partial pricing
    428429void
    429 ClpDynamicExampleMatrix::partialPricing(ClpSimplex * model, int start, int end,
     430ClpDynamicExampleMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    430431                              int & bestSequence, int & numberWanted)
    431432{
     433  numberWanted=currentWanted_;
    432434  assert(!model->rowScale());
    433   if (numberSets_) {
    434     // Do packed part before gub
    435     // always???
    436     //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
    437     //   start,min(end,firstAvailable_),end,firstAvailable_);
    438     ClpPackedMatrix::partialPricing(model,0,firstDynamic_,bestSequence,numberWanted);
    439   } else {
     435  if (!numberSets_) {
    440436    // no gub
    441     ClpPackedMatrix::partialPricing(model,start,end,bestSequence,numberWanted);
    442     return;
    443   }
    444   if (numberWanted>0) {
     437    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
     438  } else {
    445439    // and do some proportion of full set
    446     double ratio = ((double) numberSets_)/((double) lastDynamic_-firstDynamic_);
    447     int startG = max (start,firstDynamic_);
    448     int endG = min(lastDynamic_,end);
    449     int startG2 = (int) (ratio* (startG-firstDynamic_));
    450     int endG2 = ((int) (startG2 + ratio *(endG-startG)))+1;
     440    int startG2 = (int) (startFraction*numberSets_);
     441    int endG2 = (int) (endFraction*numberSets_+0.1);
    451442    endG2 = min(endG2,numberSets_);
    452443    //printf("gub price - set start %d end %d\n",
    … …  
    460451    int structuralOffset = slackOffset+numberSets_;
    461452    int structuralOffset2 = structuralOffset+maximumGubColumns_;
    462     int saveWanted=numberWanted;
    463453    // If nothing found yet can go all the way to end
    464454    int endAll = endG2;
    … …  
    479469    //     startG2,endG2,numberWanted);
    480470    int bestSet=-1;
     471    int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
     472    int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
    481473    for (int iSet = startG2;iSet<endAll;iSet++) {
    482       //if (numberWanted+50<saveWanted&&iSet>startG2+5) {
    483       if (numberWanted+5<saveWanted&&iSet>startG2+5) {
     474      if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
    484475        // give up
    485476        numberWanted=0;
    486477        break;
    487       } else if (iSet==endG2-1&&bestSequence>=0) {
     478      } else if (iSet==endG2&&bestSequence>=0) {
    488479        break;
    489480      }
    … …  
    618609      savedBestSet_ = bestSet;
    619610    }
     611    // Do packed part before gub
     612    // always???
     613    // Resize so just do to gub
     614    numberActiveColumns_=firstDynamic_;
     615    int saveMinNeg=minimumGoodReducedCosts_;
     616    if (bestSequence>=0)
     617      minimumGoodReducedCosts_=-2;
     618    currentWanted_=numberWanted;
     619    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
     620    numberActiveColumns_=matrix_->getNumCols();
     621    minimumGoodReducedCosts_=saveMinNeg;
    620622    // See if may be finished
    621623    if (!startG2&&bestSequence<0)
    622624      infeasibilityWeight_=model_->infeasibilityCost();
    623     else
     625    else if (bestSequence>=0)
    624626      infeasibilityWeight_=-1.0;
     627    currentWanted_=numberWanted;
    625628  }
    626629}

  • trunk/ClpDynamicMatrix.cpp

    r335 r336  
    3232    sumOfRelaxedPrimalInfeasibilities_(0.0),
    3333    savedBestGubDual_(0.0),
    34     savedBestDj_(0.0),
    35     savedBestSequence_(-1),
    3634    savedBestSet_(0),
    3735    backToPivotRow_(NULL),
    … …  
    103101  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    104102  savedBestGubDual_=rhs.savedBestGubDual_;
    105   savedBestDj_ = rhs.savedBestDj_;
    106   savedBestSequence_=rhs.savedBestSequence_;
    107103  savedBestSet_=rhs.savedBestSet_;
    108104  noCheck_ = rhs.noCheck_;
    … …  
    163159  numberStaticRows_ = numberRows;
    164160  savedBestGubDual_=0.0;
    165   savedBestSequence_=-1;
    166   savedBestDj_ = 0.0;
    167161  savedBestSet_=0;
    168162  // Number of columns needed
    … …  
    365359    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    366360    savedBestGubDual_=rhs.savedBestGubDual_;
    367     savedBestDj_ = rhs.savedBestDj_;
    368     savedBestSequence_=rhs.savedBestSequence_;
    369361    savedBestSet_=rhs.savedBestSet_;
    370362    noCheck_ = rhs.noCheck_;
    … …  
    396388// Partial pricing
    397389void
    398 ClpDynamicMatrix::partialPricing(ClpSimplex * model, int start, int end,
     390ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    399391                              int & bestSequence, int & numberWanted)
    400392{
     393  numberWanted=currentWanted_;
    401394  assert(!model->rowScale());
    402395  if (numberSets_) {
    … …  
    404397    // always???
    405398    //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
    406     //   start,min(end,firstAvailable_),end,firstAvailable_);
    407     ClpPackedMatrix::partialPricing(model,0,firstDynamic_,bestSequence,numberWanted);
     399    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
    408400  } else {
    409401    // no gub
    410     ClpPackedMatrix::partialPricing(model,start,end,bestSequence,numberWanted);
     402    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
    411403    return;
    412404  }
    413405  if (numberWanted>0) {
    414406    // and do some proportion of full set
    415     double ratio = ((double) numberSets_)/((double) lastDynamic_-firstDynamic_);
    416     int startG = max (start,firstDynamic_);
    417     int endG = min(lastDynamic_,end);
    418     int startG2 = (int) (ratio* (startG-firstDynamic_));
    419     int endG2 = ((int) (startG2 + ratio *(endG-startG)))+1;
     407    int startG2 = (int) (startFraction*numberSets_);
     408    int endG2 = (int) (endFraction*numberSets_+0.1);
    420409    endG2 = min(endG2,numberSets_);
    421410    //printf("gub price - set start %d end %d\n",
    … …  
    428417    int slackOffset = lastDynamic_+numberRows;
    429418    int structuralOffset = slackOffset+numberSets_;
    430     int saveWanted=numberWanted;
    431419    // If nothing found yet can go all the way to end
    432420    int endAll = endG2;
    … …  
    458446    }
    459447#endif
     448    int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
     449    int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
    460450    for (int iSet = startG2;iSet<endAll;iSet++) {
    461       //if (numberWanted+50<saveWanted&&iSet>startG2+5) {
    462       if (numberWanted+5<saveWanted&&iSet>startG2+5) {
     451      if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
    463452        // give up
    464453        numberWanted=0;
    465454        break;
    466       } else if (iSet==endG2-1&&bestSequence>=0) {
     455      } else if (iSet==endG2&&bestSequence>=0) {
    467456        break;
    468457      }
    … …  
    568557    if (!startG2&&bestSequence<0)
    569558      infeasibilityWeight_=model_->infeasibilityCost();
    570     else
     559    else if (bestSequence>=0)
    571560      infeasibilityWeight_=-1.0;
    572561  }
     562  currentWanted_=numberWanted;
    573563}
    574564/* Returns effective RHS if it is being used.  This is used for long problems
    … …  
    17111701        modifyOffset(key,valueOfKey);
    17121702        rhsOffset_[newRow] = -shift; // sign?
     1703#ifdef CLP_DEBUG
     1704        model->rowArray(1)->checkClear();
     1705#endif
    17131706        // now pivot in
    1714         unpack(model,model->rowArray(3),firstAvailable_);
    1715         model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
    1716         double alpha = model->rowArray(3)->denseVector()[newRow];
     1707        unpack(model,model->rowArray(1),firstAvailable_);
     1708        model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(1));
     1709        double alpha = model->rowArray(1)->denseVector()[newRow];
    17171710        int updateStatus =
    1718           model->factorization()->replaceColumn(model->rowArray(2),
    1719                                                                  newRow, alpha);
    1720         model->rowArray(3)->clear();
     1711          model->factorization()->replaceColumn(model,
     1712                                                model->rowArray(2),
     1713                                                model->rowArray(1),
     1714                                                newRow, alpha);
     1715        model->rowArray(1)->clear();
    17211716        if (updateStatus) {
    17221717          if (updateStatus==3) {

  • trunk/ClpFactorization.cpp

    r328 r336  
    171171                                           numberColumnBasic,
    172172                                           NULL,NULL,NULL);
    173 
     173        // and recompute as network side say different
     174        if (model->numberIterations())
     175          numberRowBasic = numberBasic - numberColumnBasic;
    174176        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
    175         getAreas ( numberRows, numberRowBasic+numberColumnBasic, numberElements,
     177        // If iteration not zero then may be compressed
     178        getAreas ( !model->numberIterations() ? numberRows : numberBasic,
     179                   numberRowBasic+numberColumnBasic, numberElements,
    176180                   2 * numberElements );
    177181        //fill
    … …  
    214218      // If we get here status is 0 or -1
    215219      if (status_ == 0) {
     220        // We may need to tamper with order and redo - e.g. network with side
     221        int useNumberRows = numberRows;
     222        // **** we will also need to add test in dual steepest to do
     223        // as we do for network
     224        matrix->generalExpanded(model,12,useNumberRows);
    216225        int * permuteBack = permuteBack_;
    217226        int * back = pivotColumnBack_;
    218227        //int * pivotTemp = pivotColumn_;
    219         //ClpDisjointCopyN ( pivotVariable, numberRows_ , pivotTemp  );
     228        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
    220229        // Redo pivot order
    221230        for (i=0;i<numberRowBasic;i++) {
    … …  
    224233          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
    225234        }
    226         for (;i<numberRows;i++) {
     235        for (;i<useNumberRows;i++) {
    227236          int k = pivotTemp[i];
    228237          // so rowIsBasic[k] would be permuteBack[back[i]]
    … …  
    232241        // these arrays start off as copies of permute
    233242        // (and we could use permute_ instead of pivotColumn (not back though))
    234         ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
    235         ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
     243        ClpDisjointCopyN ( permute_, useNumberRows , pivotColumn_  );
     244        ClpDisjointCopyN ( permuteBack_, useNumberRows , pivotColumnBack_  );
    236245        if (networkMatrix) {
    237246          maximumPivots(saveMaximumPivots);
    … …  
    249258          if (!doCheck) {
    250259            gutsOfDestructor();
     260            status_=0;
    251261#if 0
    252262            // but put back permute arrays so odd things will work
    … …  
    460470   partial update already in U */
    461471int
    462 ClpFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
    463                       int pivotRow,
    464                       double pivotCheck ,
    465                       bool checkBeforeModifying)
     472ClpFactorization::replaceColumn ( const ClpSimplex * model,
     473                                  CoinIndexedVector * regionSparse,
     474                                  CoinIndexedVector * tableauColumn,
     475                                  int pivotRow,
     476                                  double pivotCheck ,
     477                                  bool checkBeforeModifying)
    466478{
    467479  if (!networkBasis_) {
    468     return CoinFactorization::replaceColumn(regionSparse,
    469                                             pivotRow,
    470                                             pivotCheck,
    471                                             checkBeforeModifying);
     480    // see if FT
     481    if (doForrestTomlin_)
     482      return CoinFactorization::replaceColumn(regionSparse,
     483                                              pivotRow,
     484                                              pivotCheck,
     485                                              checkBeforeModifying);
     486    else
     487      return CoinFactorization::replaceColumnPFI(tableauColumn,
     488                                              pivotRow,pivotCheck); // Note array
     489
    472490  } else {
    473491    if (doCheck) {

  • trunk/ClpGubDynamicMatrix.cpp

    r335 r336  
    1616#include "ClpGubDynamicMatrix.hpp"
    1717#include "ClpMessage.hpp"
    18 #define CLP_DEBUG
     18//#define CLP_DEBUG
    1919//#define CLP_DEBUG_PRINT
    2020//#############################################################################
    … …  
    3939    lowerSet_(NULL),
    4040    upperSet_(NULL),
    41     model_(NULL),
    4241    numberGubColumns_(0),
    4342    firstAvailable_(0),
    … …  
    7574  lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
    7675  upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
    77   model_ = rhs.model_;
    7876}
    7977
    … …  
    188186  backward_ = new int[numberNeeded];
    189187  backToPivotRow_ = new int[numberNeeded];
     188  // We know a bit better
     189  delete [] changeCost_;
    190190  changeCost_ = new double [numberRows+numberSets_];
    191191  keyVariable_ = new int[numberSets_];
    … …  
    272272    lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
    273273    upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
    274     model_ = rhs.model_;
    275274  }
    276275  return *this;
    … …  
    285284// Partial pricing
    286285void
    287 ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, int start, int end,
     286ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    288287                              int & bestSequence, int & numberWanted)
    289288{
    290289  assert(!model->rowScale());
    291   if (numberSets_) {
    292     // Do packed part before gub
    293     // always???
    294     //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
    295     //   start,min(end,firstAvailable_),end,firstAvailable_);
    296     ClpPackedMatrix::partialPricing(model,0,firstDynamic_,bestSequence,numberWanted);
    297     ClpGubMatrix::partialPricing(model,start,min(firstAvailable_,end),bestSequence,numberWanted);
     290  numberWanted=currentWanted_;
     291  if (!numberSets_) {
     292    // no gub
     293    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
     294    return;
    298295  } else {
    299     // no gub
    300     ClpPackedMatrix::partialPricing(model,start,end,bestSequence,numberWanted);
    301     return;
    302   }
    303   if (numberWanted>0) {
    304296    // and do some proportion of full set
    305     double ratio = ((double) numberSets_)/((double) lastDynamic_-firstDynamic_);
    306     int startG = max (start,firstDynamic_);
    307     int endG = min(lastDynamic_,end);
    308     int startG2 = (int) (ratio* (startG-firstDynamic_));
    309     int endG2 = ((int) (startG2 + ratio *(endG-startG)))+1;
     297    int startG2 = (int) (startFraction*numberSets_);
     298    int endG2 = (int) (endFraction*numberSets_+0.1);
    310299    endG2 = min(endG2,numberSets_);
    311300    //printf("gub price - set start %d end %d\n",
    … …  
    318307    int numberRows = model->numberRows();
    319308    int numberColumns = lastDynamic_;
    320     int saveWanted=numberWanted;
     309    // If nothing found yet can go all the way to end
     310    int endAll = endG2;
     311    if (bestSequence<0&&!startG2)
     312      endAll = numberSets_;
    321313    if (bestSequence>=0)
    322314      bestDj = fabs(reducedCost[bestSequence]);
    … …  
    346338    }
    347339#endif
    348     for (int iSet = startG2;iSet<endG2;iSet++) {
    349       if (numberWanted+5<saveWanted&&iSet>startG2+5) {
     340    int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
     341    int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
     342    for (int iSet = startG2;iSet<endAll;iSet++) {
     343      if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
    350344        // give up
    351345        numberWanted=0;
     346        break;
     347      } else if (iSet==endG2&&bestSequence>=0) {
    352348        break;
    353349      }
    … …  
    442438        break;
    443439      }
     440    }
     441    // Do packed part before gub and small gub - but lightly
     442    int saveMinNeg=minimumGoodReducedCosts_;
     443    int saveSequence2 = bestSequence;
     444    if (bestSequence>=0)
     445      minimumGoodReducedCosts_=-2;
     446    int saveLast=lastGub_;
     447    lastGub_=firstAvailable_;
     448    currentWanted_=numberWanted;
     449    ClpGubMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
     450    minimumGoodReducedCosts_=saveMinNeg;
     451    lastGub_=saveLast;
     452    if (bestSequence!=saveSequence2) {
     453      bestType=-1; // in normal or small gub part
     454      saveSequence=bestSequence;
    444455    }
    445456    if (bestSequence!=saveSequence||bestType>=0) {
    … …  
    531542                                       upperColumn[bestSequence],0.0);
    532543      }
    533     }
    534   }
     544      savedBestSequence_ = bestSequence;
     545      savedBestDj_ = reducedCost[savedBestSequence_];
     546    }
     547    // See if may be finished
     548    if (!startG2&&bestSequence<0)
     549      infeasibilityWeight_=model_->infeasibilityCost();
     550    else if (bestSequence>=0)
     551      infeasibilityWeight_=-1.0;
     552  }
     553  currentWanted_=numberWanted;
    535554}
    536555// This is local to Gub to allow synchronization when status is good
    … …  
    10401059  for (iSet=0;iSet<numberSets_;iSet++)
    10411060    toIndex_[iSet]=-1;
    1042   fromIndex_ = new int [max(getNumRows()+1,numberSets_)];
     1061  fromIndex_ = new int [numberRows+numberSets_];
    10431062  double tolerance = model->primalTolerance();
    10441063  double * element =  matrix_->getMutableElements();
    … …  
    17911810  bool print=false;
    17921811  int iSet;
     1812  int trueIn=-1;
     1813  int trueOut=-1;
     1814  int numberRows = model->numberRows();
     1815  int numberColumns = model->numberColumns();
    17931816  if (sequenceIn==firstAvailable_) {
    17941817    if (doPrinting)
    … …  
    18111834    if (iSet>=0) {
    18121835      int bigSequence = id_[sequenceIn-firstDynamic_];
     1836      trueIn=bigSequence+numberRows+numberColumns+numberSets_;
    18131837      if (doPrinting)
    18141838        printf(" incoming set %d big seq %d",iSet,bigSequence);
    18151839      print=true;
    18161840    }
     1841  } else if (sequenceIn>=numberRows+numberColumns) {
     1842    trueIn = numberRows+numberColumns+gubSlackIn_;
    18171843  }
    18181844  if (sequenceOut<lastDynamic_) {
    … …  
    18201846    if (iSet>=0) {
    18211847      int bigSequence = id_[sequenceOut-firstDynamic_];
     1848      trueOut=bigSequence+firstDynamic_;
    18221849      if (model->getStatus(sequenceOut)==ClpSimplex::atUpperBound)
    18231850        setDynamicStatus(bigSequence,atUpperBound);
    … …  
    18351862    printf("\n");
    18361863  ClpGubMatrix::updatePivot(model,oldInValue,oldOutValue);
     1864  // Redo true in and out
     1865  if (trueIn>=0)
     1866    trueSequenceIn_=trueIn;
     1867  if (trueOut>=0)
     1868    trueSequenceOut_=trueOut;
    18371869  if (doPrinting&&0) {
    18381870    for (int i=0;i<numberSets_;i++) {
    … …  
    19271959  }
    19281960}
     1961/* Just for debug - may be extended to other matrix types later.
     1962   Returns number of primal infeasibilities.
     1963*/
     1964int
     1965ClpGubDynamicMatrix::checkFeasible() const
     1966{
     1967  int numberRows = model_->numberRows();
     1968  double * rhs = new double[numberRows];
     1969  int numberColumns = model_->numberColumns();
     1970  int iRow;
     1971  CoinZeroN(rhs,numberRows);
     1972  // do ones at bounds before gub
     1973  const double * smallSolution = model_->solutionRegion();
     1974  const double * element =matrix_->getElements();
     1975  const int * row = matrix_->getIndices();
     1976  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
     1977  const int * length = matrix_->getVectorLengths();
     1978  int iColumn;
     1979  int numberInfeasible=0;
     1980  const double * rowLower = model_->rowLower();
     1981  const double * rowUpper = model_->rowUpper();
     1982  for (iRow=0;iRow<numberRows;iRow++) {
     1983    double value=smallSolution[numberColumns+iRow];
     1984    if (value<rowLower[iRow]-1.0e-5||
     1985        value>rowUpper[iRow]+1.0e-5) {
     1986      //printf("row %d %g %g %g\n",
     1987      //     iRow,rowLower[iRow],value,rowUpper[iRow]);
     1988      numberInfeasible++;
     1989    }
     1990    rhs[iRow]=value;
     1991  }
     1992  const double * columnLower = model_->columnLower();
     1993  const double * columnUpper = model_->columnUpper();
     1994  for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
     1995    double value=smallSolution[iColumn];
     1996    if (value<columnLower[iColumn]-1.0e-5||
     1997        value>columnUpper[iColumn]+1.0e-5) {
     1998      //printf("column %d %g %g %g\n",
     1999      //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
     2000      numberInfeasible++;
     2001    }
     2002    for (CoinBigIndex j=startColumn[iColumn];
     2003         j<startColumn[iColumn]+length[iColumn];j++) {
     2004      int jRow=row[j];
     2005      rhs[jRow] -= value*element[j];
     2006    }
     2007  }
     2008  double * solution = new double [numberGubColumns_];
     2009  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
     2010    double value=0.0;
     2011    if(getDynamicStatus(iColumn)==atUpperBound)
     2012      value = upperColumn_[iColumn];
     2013    else if (lowerColumn_)
     2014      value = lowerColumn_[iColumn];
     2015    solution[iColumn]=value;
     2016  }
     2017  // ones in small and gub
     2018  for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
     2019    int jFull = id_[iColumn-firstDynamic_];
     2020    solution[jFull]=smallSolution[iColumn];
     2021  }
     2022  // fill in all basic in small model
     2023  int * pivotVariable = model_->pivotVariable();
     2024  for (iRow=0;iRow<numberRows;iRow++) {
     2025    int iColumn = pivotVariable[iRow];
     2026    if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
     2027      int iSequence = id_[iColumn-firstDynamic_];
     2028      solution[iSequence]=smallSolution[iColumn];
     2029    }
     2030  }
     2031  // and now compute value to use for key
     2032  ClpSimplex::Status iStatus;
     2033  for (int iSet=0;iSet<numberSets_;iSet++) {
     2034    iColumn = keyVariable_[iSet];
     2035    if (iColumn<numberColumns) {
     2036      int iSequence = id_[iColumn-firstDynamic_];
     2037      solution[iSequence]=0.0;
     2038      double b=0.0;
     2039      // key is structural - where is slack
     2040      iStatus = getStatus(iSet);
     2041      assert (iStatus!=ClpSimplex::basic);
     2042      if (iStatus==ClpSimplex::atLowerBound)
     2043        b=lower_[iSet];
     2044      else
     2045        b=upper_[iSet];
     2046      // subtract out others at bounds
     2047      for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++)
     2048        b -= solution[j];
     2049      solution[iSequence]=b;
     2050    }
     2051  }
     2052  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
     2053    double value = solution[iColumn];
     2054    if ((lowerColumn_&&value<lowerColumn_[iColumn]-1.0e-5)||
     2055        (!lowerColumn_&&value<-1.0e-5)||
     2056        (upperColumn_&&value>upperColumn_[iColumn]+1.0e-5)) {
     2057      //printf("column %d %g %g %g\n",
     2058      //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
     2059      numberInfeasible++;
     2060    }
     2061    if (value) {
     2062      for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
     2063        int iRow = row_[j];
     2064        rhs[iRow] -= element_[j]*value;
     2065      }
     2066    }
     2067  }
     2068  for (iRow=0;iRow<numberRows;iRow++) {
     2069    if (fabs(rhs[iRow])>1.0e-5)
     2070      printf("rhs mismatch %d %g\n",iRow,rhs[iRow]);
     2071  }
     2072  delete [] solution;
     2073  delete [] rhs;
     2074  return numberInfeasible;
     2075}

  • trunk/ClpGubMatrix.cpp

    r335 r336  
    3131    sumOfRelaxedDualInfeasibilities_(0.0),
    3232    sumOfRelaxedPrimalInfeasibilities_(0.0),
     33    infeasibilityWeight_(0.0),
    3334    start_(NULL),
    3435    end_(NULL),
    … …  
    4546    toIndex_(NULL),
    4647    fromIndex_(NULL),
     48    model_(NULL),
    4749    numberDualInfeasibilities_(0),
    4850    numberPrimalInfeasibilities_(0),
     51    noCheck_(-1),
    4952    numberSets_(0),
    5053    saveNumber_(0),
    … …  
    7982  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,numberColumns);
    8083  changeCost_ = ClpCopyOfArray(rhs.changeCost_,getNumRows()+numberSets_);
     84  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+numberSets_+1);
    8185  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
    8286  // find longest set
    … …  
    9498  next_ = ClpCopyOfArray(rhs.next_,numberColumns+numberSets_+2*length);
    9599  toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
    96   fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1);
    97100  sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
    98101  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    99102  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    100103  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     104  infeasibilityWeight_=rhs.infeasibilityWeight_;
    101105  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    102106  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
     107  noCheck_ = rhs.noCheck_;
    103108  firstGub_ = rhs.firstGub_;
    104109  lastGub_ = rhs.lastGub_;
    105110  gubType_ = rhs.gubType_;
     111  model_ = rhs.model_;
    106112}
    107113
    … …  
    115121    sumOfRelaxedDualInfeasibilities_(0.0),
    116122    sumOfRelaxedPrimalInfeasibilities_(0.0),
     123    infeasibilityWeight_(0.0),
    117124    start_(NULL),
    118125    end_(NULL),
    … …  
    129136    toIndex_(NULL),
    130137    fromIndex_(NULL),
     138    model_(NULL),
    131139    numberDualInfeasibilities_(0),
    132140    numberPrimalInfeasibilities_(0),
     141    noCheck_(-1),
    133142    numberSets_(0),
    134143    saveNumber_(0),
    … …  
    159168    gubSlackIn_(-1)
    160169{
     170  model_=NULL;
    161171  numberSets_ = numberSets;
    162172  start_ = ClpCopyOfArray(start,numberSets_);
    … …  
    230240  savedKeyVariable_= new int [numberSets_];
    231241  memset(savedKeyVariable_,0,numberSets_*sizeof(int));
     242  noCheck_ = -1;
     243  infeasibilityWeight_=0.0;
    232244}
    233245
    … …  
    238250    sumOfRelaxedDualInfeasibilities_(0.0),
    239251    sumOfRelaxedPrimalInfeasibilities_(0.0),
     252    infeasibilityWeight_(0.0),
    240253    start_(NULL),
    241254    end_(NULL),
    … …  
    252265    toIndex_(NULL),
    253266    fromIndex_(NULL),
     267    model_(NULL),
    254268    numberDualInfeasibilities_(0),
    255269    numberPrimalInfeasibilities_(0),
     270    noCheck_(-1),
    256271    numberSets_(0),
    257272    saveNumber_(0),
    … …  
    324339    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,numberColumns);
    325340    changeCost_ = ClpCopyOfArray(rhs.changeCost_,getNumRows()+numberSets_);
     341    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+numberSets_+1);
    326342    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
    327343    // find longest set
    … …  
    339355    next_ = ClpCopyOfArray(rhs.next_,numberColumns+numberSets_+2*length);
    340356    toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
    341     fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1);
    342357    sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
    343358    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    344359    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    345360    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     361    infeasibilityWeight_=rhs.infeasibilityWeight_;
    346362    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    347363    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
     364    noCheck_ = rhs.noCheck_;
    348365    firstGub_ = rhs.firstGub_;
    349366    lastGub_ = rhs.lastGub_;
    350367    gubType_ = rhs.gubType_;
     368    model_ = rhs.model_;
    351369  }
    352370  return *this;
    … …  
    530548  // reduce for gub
    531549  factor *= 0.5;
     550  assert (!y->getNumElements());
    532551  if (numberInRowArray>factor*numberRows||!rowCopy) {
    533552    // do by column
    … …  
    542561    int iSet = -1;
    543562    double djMod=0.0;
    544     if (!y->getNumElements()) {
    545       if (packed) {
    546         // need to expand pi into y
    547         assert(y->capacity()>=numberRows);
    548         double * piOld = pi;
    549         pi = y->denseVector();
    550         const int * whichRow = rowArray->getIndices();
    551         int i;
    552         if (!rowScale) {
    553           // modify pi so can collapse to one loop
    554           for (i=0;i<numberInRowArray;i++) {
    555             int iRow = whichRow[i];
    556             pi[iRow]=scalar*piOld[i];
    557           }
    558           for (iColumn=0;iColumn<numberColumns;iColumn++) {
    559             if (backward_[iColumn]!=iSet) {
    560               // get pi on gub row
    561               iSet = backward_[iColumn];
    562               if (iSet>=0) {
    563                 int iBasic = keyVariable_[iSet];
    564                 if (iBasic<numberColumns) {
    565                   // get dj without
    566                   assert (model->getStatus(iBasic)==ClpSimplex::basic);
    567                   djMod=0.0;
    568                   for (CoinBigIndex j=columnStart[iBasic];
    569                        j<columnStart[iBasic]+columnLength[iBasic];j++) {
    570                     int jRow=row[j];
    571                     djMod -= pi[jRow]*elementByColumn[j];
    572                   }
    573                 } else {
    574                   djMod = 0.0;
     563    if (packed) {
     564      // need to expand pi into y
     565      assert(y->capacity()>=numberRows);
     566      double * piOld = pi;
     567      pi = y->denseVector();
     568      const int * whichRow = rowArray->getIndices();
     569      int i;
     570      if (!rowScale) {
     571        // modify pi so can collapse to one loop
     572        for (i=0;i<numberInRowArray;i++) {
     573          int iRow = whichRow[i];
     574          pi[iRow]=scalar*piOld[i];
     575        }
     576        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     577          if (backward_[iColumn]!=iSet) {
     578            // get pi on gub row
     579            iSet = backward_[iColumn];
     580            if (iSet>=0) {
     581              int iBasic = keyVariable_[iSet];
     582              if (iBasic<numberColumns) {
     583                // get dj without
     584                assert (model->getStatus(iBasic)==ClpSimplex::basic);
     585                djMod=0.0;
     586                for (CoinBigIndex j=columnStart[iBasic];
     587                     j<columnStart[iBasic]+columnLength[iBasic];j++) {
     588                  int jRow=row[j];
     589                  djMod -= pi[jRow]*elementByColumn[j];
    575590                }
    576591              } else {
    577592                djMod = 0.0;
    578593              }
    579             }
    580             double value = -djMod;
     594            } else {
     595              djMod = 0.0;
     596            }
     597          }
     598          double value = -djMod;
     599          CoinBigIndex j;
     600          for (j=columnStart[iColumn];
     601               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     602            int iRow = row[j];
     603            value += pi[iRow]*elementByColumn[j];
     604          }
     605          if (fabs(value)>zeroTolerance) {
     606            array[numberNonZero]=value;
     607            index[numberNonZero++]=iColumn;
     608          }
     609        }
     610      } else {
     611        // scaled
     612        // modify pi so can collapse to one loop
     613        for (i=0;i<numberInRowArray;i++) {
     614          int iRow = whichRow[i];
     615          pi[iRow]=scalar*piOld[i]*rowScale[iRow];
     616        }
     617        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     618          if (backward_[iColumn]!=iSet) {
     619            // get pi on gub row
     620            iSet = backward_[iColumn];
     621            if (iSet>=0) {
     622              int iBasic = keyVariable_[iSet];
     623              if (iBasic<numberColumns) {
     624                // get dj without
     625                assert (model->getStatus(iBasic)==ClpSimplex::basic);
     626                djMod=0.0;
     627                // scaled
     628                for (CoinBigIndex j=columnStart[iBasic];
     629                     j<columnStart[iBasic]+columnLength[iBasic];j++) {
     630                  int jRow=row[j];
     631                  djMod -= pi[jRow]*elementByColumn[j]*rowScale[jRow];
     632                }
     633              } else {
     634                djMod = 0.0;
     635              }
     636            } else {
     637              djMod = 0.0;
     638            }
     639          }
     640          double value = -djMod;
     641          CoinBigIndex j;
     642          const double * columnScale = model->columnScale();
     643          for (j=columnStart[iColumn];
     644               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     645            int iRow = row[j];
     646            value += pi[iRow]*elementByColumn[j];
     647          }
     648          value *= columnScale[iColumn];
     649          if (fabs(value)>zeroTolerance) {
     650            array[numberNonZero]=value;
     651            index[numberNonZero++]=iColumn;
     652          }
     653        }
     654      }
     655      // zero out
     656      for (i=0;i<numberInRowArray;i++) {
     657        int iRow = whichRow[i];
     658        pi[iRow]=0.0;
     659      }
     660    } else {
     661      // code later
     662      assert (packed);
     663      if (!rowScale) {
     664        if (scalar==-1.0) {
     665          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     666            double value = 0.0;
    581667            CoinBigIndex j;
    582668            for (j=columnStart[iColumn];
    … …  
    586672            }
    587673            if (fabs(value)>zeroTolerance) {
    588               array[numberNonZero]=value;
    589674              index[numberNonZero++]=iColumn;
     675              array[iColumn]=-value;
     676            }
     677          }
     678        } else if (scalar==1.0) {
     679          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     680            double value = 0.0;
     681            CoinBigIndex j;
     682            for (j=columnStart[iColumn];
     683                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     684              int iRow = row[j];
     685              value += pi[iRow]*elementByColumn[j];
     686            }
     687            if (fabs(value)>zeroTolerance) {
     688              index[numberNonZero++]=iColumn;
     689              array[iColumn]=value;
    590690            }
    591691          }
    592692        } else {
    593           // scaled
    594           // modify pi so can collapse to one loop
    595           for (i=0;i<numberInRowArray;i++) {
    596             int iRow = whichRow[i];
    597             pi[iRow]=scalar*piOld[i]*rowScale[iRow];
    598           }
    599693          for (iColumn=0;iColumn<numberColumns;iColumn++) {
    600             if (backward_[iColumn]!=iSet) {
    601               // get pi on gub row
    602               iSet = backward_[iColumn];
    603               if (iSet>=0) {
    604                 int iBasic = keyVariable_[iSet];
    605                 if (iBasic<numberColumns) {
    606                   // get dj without
    607                   assert (model->getStatus(iBasic)==ClpSimplex::basic);
    608                   djMod=0.0;
    609                   // scaled
    610                   for (CoinBigIndex j=columnStart[iBasic];
    611                        j<columnStart[iBasic]+columnLength[iBasic];j++) {
    612                     int jRow=row[j];
    613                     djMod -= pi[jRow]*elementByColumn[j]*rowScale[jRow];
    614                   }
    615                 } else {
    616                   djMod = 0.0;
    617                 }
    618               } else {
    619                 djMod = 0.0;
    620               }
    621             }
    622             double value = -djMod;
     694            double value = 0.0;
     695            CoinBigIndex j;
     696            for (j=columnStart[iColumn];
     697                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     698              int iRow = row[j];
     699              value += pi[iRow]*elementByColumn[j];
     700            }
     701            value *= scalar;
     702            if (fabs(value)>zeroTolerance) {
     703              index[numberNonZero++]=iColumn;
     704              array[iColumn]=value;
     705            }
     706          }
     707        }
     708      } else {
     709        // scaled
     710        if (scalar==-1.0) {
     711          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     712            double value = 0.0;
    623713            CoinBigIndex j;
    624714            const double * columnScale = model->columnScale();
    … …  
    626716                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
    627717              int iRow = row[j];
    628               value += pi[iRow]*elementByColumn[j];
     718              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    629719            }
    630720            value *= columnScale[iColumn];
    631721            if (fabs(value)>zeroTolerance) {
    632               array[numberNonZero]=value;
    633722              index[numberNonZero++]=iColumn;
    634             }
    635           }
    636         }
    637         // zero out
    638         for (i=0;i<numberInRowArray;i++) {
    639           int iRow = whichRow[i];
    640           pi[iRow]=0.0;
    641         }
    642       } else {
    643         // code later
    644         assert (packed);
    645         if (!rowScale) {
    646           if (scalar==-1.0) {
    647             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    648               double value = 0.0;
    649               CoinBigIndex j;
    650               for (j=columnStart[iColumn];
    651                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    652                 int iRow = row[j];
    653                 value += pi[iRow]*elementByColumn[j];
    654               }
    655               if (fabs(value)>zeroTolerance) {
    656                 index[numberNonZero++]=iColumn;
    657                 array[iColumn]=-value;
    658               }
    659             }
    660           } else if (scalar==1.0) {
    661             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    662               double value = 0.0;
    663               CoinBigIndex j;
    664               for (j=columnStart[iColumn];
    665                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    666                 int iRow = row[j];
    667                 value += pi[iRow]*elementByColumn[j];
    668               }
    669               if (fabs(value)>zeroTolerance) {
    670                 index[numberNonZero++]=iColumn;
    671                 array[iColumn]=value;
    672               }
    673             }
    674           } else {
    675             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    676               double value = 0.0;
    677               CoinBigIndex j;
    678               for (j=columnStart[iColumn];
    679                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    680                 int iRow = row[j];
    681                 value += pi[iRow]*elementByColumn[j];
    682               }
    683               value *= scalar;
    684               if (fabs(value)>zeroTolerance) {
    685                 index[numberNonZero++]=iColumn;
    686                 array[iColumn]=value;
    687               }
     723              array[iColumn]=-value;
     724            }
     725          }
     726        } else if (scalar==1.0) {
     727          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     728            double value = 0.0;
     729            CoinBigIndex j;
     730            const double * columnScale = model->columnScale();
     731            for (j=columnStart[iColumn];
     732                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     733              int iRow = row[j];
     734              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     735            }
     736            value *= columnScale[iColumn];
     737            if (fabs(value)>zeroTolerance) {
     738              index[numberNonZero++]=iColumn;
     739              array[iColumn]=value;
    688740            }
    689741          }
    690742        } else {
    691           // scaled
    692           if (scalar==-1.0) {
    693             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    694               double value = 0.0;
    695               CoinBigIndex j;
    696               const double * columnScale = model->columnScale();
    697               for (j=columnStart[iColumn];
    698                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    699                 int iRow = row[j];
    700                 value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    701               }
    702               value *= columnScale[iColumn];
    703               if (fabs(value)>zeroTolerance) {
    704                 index[numberNonZero++]=iColumn;
    705                 array[iColumn]=-value;
    706               }
    707             }
    708           } else if (scalar==1.0) {
    709             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    710               double value = 0.0;
    711               CoinBigIndex j;
    712               const double * columnScale = model->columnScale();
    713               for (j=columnStart[iColumn];
    714                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    715                 int iRow = row[j];
    716                 value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    717               }
    718               value *= columnScale[iColumn];
    719               if (fabs(value)>zeroTolerance) {
    720                 index[numberNonZero++]=iColumn;
    721                 array[iColumn]=value;
    722               }
    723             }
    724           } else {
    725             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    726               double value = 0.0;
    727               CoinBigIndex j;
    728               const double * columnScale = model->columnScale();
    729               for (j=columnStart[iColumn];
    730                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    731                 int iRow = row[j];
    732                 value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    733               }
    734               value *= scalar*columnScale[iColumn];
    735               if (fabs(value)>zeroTolerance) {
    736                 index[numberNonZero++]=iColumn;
    737                 array[iColumn]=value;
    738               }
    739             }
    740           }
    741         }
    742       }
    743     } else {
    744       assert(!packed);
    745       // code later
    746       assert (!y->getNumElements());
    747       double * markVector = y->denseVector(); // not empty
    748       if (!rowScale) {
    749         for (iColumn=0;iColumn<numberColumns;iColumn++) {
    750           double value = markVector[iColumn];
    751           markVector[iColumn]=0.0;
    752           double value2 = 0.0;
    753           CoinBigIndex j;
    754           for (j=columnStart[iColumn];
    755                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    756             int iRow = row[j];
    757             value2 += pi[iRow]*elementByColumn[j];
    758           }
    759           value += value2*scalar;
    760           if (fabs(value)>zeroTolerance) {
    761             index[numberNonZero++]=iColumn;
    762             array[iColumn]=value;
    763           }
    764         }
    765       } else {
    766         // scaled
    767         for (iColumn=0;iColumn<numberColumns;iColumn++) {
    768           double value = markVector[iColumn];
    769           markVector[iColumn]=0.0;
    770           CoinBigIndex j;
    771           const double * columnScale = model->columnScale();
    772           double value2 = 0.0;
    773           for (j=columnStart[iColumn];
    774                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    775             int iRow = row[j];
    776             value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    777           }
    778           value += value2*scalar*columnScale[iColumn];
    779           if (fabs(value)>zeroTolerance) {
    780             index[numberNonZero++]=iColumn;
    781             array[iColumn]=value;
     743          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     744            double value = 0.0;
     745            CoinBigIndex j;
     746            const double * columnScale = model->columnScale();
     747            for (j=columnStart[iColumn];
     748                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     749              int iRow = row[j];
     750              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     751            }
     752            value *= scalar*columnScale[iColumn];
     753            if (fabs(value)>zeroTolerance) {
     754              index[numberNonZero++]=iColumn;
     755              array[iColumn]=value;
     756            }
    782757          }
    783758        }
    … …  
    13781353// Partial pricing
    13791354void
    1380 ClpGubMatrix::partialPricing(ClpSimplex * model, int start, int end,
     1355ClpGubMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    13811356                              int & bestSequence, int & numberWanted)
    13821357{
     1358  numberWanted=currentWanted_;
    13831359  if (numberSets_) {
    13841360    // Do packed part before gub
    1385     ClpPackedMatrix::partialPricing(model,start,firstGub_,bestSequence,numberWanted);
    1386     if (numberWanted) {
     1361    int numberColumns = matrix_->getNumCols();
     1362    double ratio = ((double) firstGub_)/((double) numberColumns);
     1363    ClpPackedMatrix::partialPricing(model,startFraction*ratio,
     1364                                    endFraction*ratio,bestSequence,numberWanted);
     1365    if (numberWanted||minimumGoodReducedCosts_<-1) {
    13871366      // do gub
    13881367      const double * element =matrix_->getElements();
    … …  
    14021381      int numberRows = model->numberRows();
    14031382      if (bestSequence>=0)
    1404         bestDj = fabs(reducedCost[bestSequence]);
     1383        bestDj = fabs(this->reducedCost(model,bestSequence));
    14051384      else
    14061385        bestDj=tolerance;
    14071386      int sequenceOut = model->sequenceOut();
    14081387      int saveSequence = bestSequence;
    1409       int startG = max (start,firstGub_);
    1410       int endG = min(lastGub_,end);
     1388      int startG = firstGub_+ (int) (startFraction* (lastGub_-firstGub_));
     1389      int endG = firstGub_+ (int) (endFraction* (lastGub_-firstGub_));
     1390      endG = min(lastGub_,endG+1);
     1391      // If nothing found yet can go all the way to end
     1392      int endAll = endG;
     1393      if (bestSequence<0&&!startG)
     1394        endAll = lastGub_;
     1395      int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
     1396      int minNeg = minimumGoodReducedCosts_==-1 ? 5 : minimumGoodReducedCosts_;
     1397      int nSets=0;
    14111398      int iSet = -1;
    14121399      double djMod=0.0;
    … …  
    14151402        double bestDjMod=0.0;
    14161403        // scaled
    1417         for (iSequence=startG;iSequence<endG;iSequence++) {
     1404        for (iSequence=startG;iSequence<endAll;iSequence++) {
     1405          if (numberWanted+minNeg<originalWanted_&&nSets>minSet) {
     1406            // give up
     1407            numberWanted=0;
     1408            break;
     1409          } else if (iSequence==endG&&bestSequence>=0) {
     1410            break;
     1411          }
    14181412          if (backward_[iSequence]!=iSet) {
    14191413            // get pi on gub row
    14201414            iSet = backward_[iSequence];
    14211415            if (iSet>=0) {
     1416              nSets++;
    14221417              int iBasic = keyVariable_[iSet];
    14231418              if (iBasic>=numberColumns) {
    … …  
    15921587            model->costRegion()[bestSequence] = 0.0;
    15931588          }
     1589          savedBestSequence_ = bestSequence;
     1590          savedBestDj_ = reducedCost[savedBestSequence_];
    15941591        }
    15951592      } else {
    … …  
    15981595        //     startG,endG,numberWanted);
    15991596        for (iSequence=startG;iSequence<endG;iSequence++) {
     1597          if (numberWanted+minNeg<originalWanted_&&nSets>minSet) {
     1598            // give up
     1599            numberWanted=0;
     1600            break;
     1601          } else if (iSequence==endG&&bestSequence>=0) {
     1602            break;
     1603          }
    16001604          if (backward_[iSequence]!=iSet) {
    16011605            // get pi on gub row
    16021606            iSet = backward_[iSequence];
    16031607            if (iSet>=0) {
     1608              nSets++;
    16041609              int iBasic = keyVariable_[iSet];
    16051610              if (iBasic>=numberColumns) {
    … …  
    17721777        }
    17731778      }
     1779      // See if may be finished
     1780      if (startG==firstGub_&&bestSequence<0)
     1781        infeasibilityWeight_=model_->infeasibilityCost();
     1782      else if (bestSequence>=0)
     1783        infeasibilityWeight_=-1.0;
    17741784    }
    17751785    if (numberWanted) {
    17761786      // Do packed part after gub
    1777       ClpPackedMatrix::partialPricing(model,max(lastGub_,start),end,bestSequence,numberWanted);
     1787      double offset = ((double) lastGub_)/((double) numberColumns);
     1788      double ratio = ((double) numberColumns)/((double) numberColumns)-offset;
     1789      double start2 = offset + ratio*startFraction;
     1790      double end2 = min(1.0,offset + ratio*endFraction+1.0e-6);
     1791      ClpPackedMatrix::partialPricing(model,start2,end2,bestSequence,numberWanted);
    17781792    }
    17791793  } else {
    17801794    // no gub
    1781     ClpPackedMatrix::partialPricing(model,start,end,bestSequence,numberWanted);
    1782   }
     1795    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
     1796  }
     1797  if (bestSequence>=0)
     1798    infeasibilityWeight_=-1.0; // not optimal
     1799  currentWanted_=numberWanted;
    17831800}
    17841801/* expands an updated column to allow for extra rows which the main
    … …  
    20462063    if (number>saveNumber_) {
    20472064      // clear
     2065      double theta = model->theta();
     2066      double * solution = model->solutionRegion();
    20482067      for (i=saveNumber_;i<number;i++) {
    2049 #ifdef CLP_DEBUG_PRINT
    20502068        int iRow = index[i];
    20512069        int iColumn = pivotVariable[iRow];
     2070#ifdef CLP_DEBUG_PRINT
    20522071        printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
    20532072               iColumn,fromIndex_[i-saveNumber_],lower[iColumn],upper[iColumn],array[i],
    20542073               solution[iColumn],solution[iColumn]-model->theta()*array[i],model->theta());
    20552074#endif
     2075        double value = array[i];
    20562076        array[i]=0.0;
    20572077        int iSet=fromIndex_[i-saveNumber_];
    20582078        toIndex_[iSet]=-1;
     2079        if (iSet==iSetIn&&iColumn<numberColumns) {
     2080          // update as may need value
     2081          solution[iColumn] -= theta*value;
     2082        }
    20592083      }
    20602084    }
    … …  
    22952319          backToPivotRow_[iPivot]=i;
    22962320      }
     2321      if (noCheck_>=0) {
     2322        if (infeasibilityWeight_!=model->infeasibilityCost()) {
     2323          // don't bother checking
     2324          sumDualInfeasibilities_=100.0;
     2325          numberDualInfeasibilities_=1;
     2326          sumOfRelaxedDualInfeasibilities_ =100.0;
     2327          return;
     2328        }
     2329      }
    22972330      double * dj = model->djRegion();
    22982331      double * dual = model->dualRowSolution();
    … …  
    24092442      // and get statistics for column generation
    24102443      synchronize(model,4);
     2444      infeasibilityWeight_=-1.0;
    24112445    }
    24122446    break;
    … …  
    26052639      }
    26062640      number = numberBasic;
     2641      if (model->numberIterations())
     2642        assert (number==model->numberRows());
    26072643    }
    26082644    break;
    … …  
    27202756                //if (!alpha)
    27212757                //printf("zero alpha a\n");
    2722                 int updateStatus = model->factorization()->replaceColumn(model->rowArray(2),
     2758                int updateStatus = model->factorization()->replaceColumn(model,
     2759                                                                         model->rowArray(2),
     2760                                                                         model->rowArray(3),
    27232761                                                                         iRow, alpha);
    27242762                returnCode = max (updateStatus, returnCode);
    … …  
    27412779              //if (!alpha)
    27422780              //printf("zero alpha b\n");
    2743               int updateStatus = model->factorization()->replaceColumn(model->rowArray(2),
     2781              int updateStatus = model->factorization()->replaceColumn(model,
     2782                                                                       model->rowArray(2),
     2783                                                                       model->rowArray(3),
    27442784                                                                       pivotRow, alpha);
    27452785              returnCode = max (updateStatus, returnCode);
    … …  
    27542794            printf("TTTTTTry 4 %d\n",possiblePivotKey_);
    27552795#endif
    2756             int updateStatus = model->factorization()->replaceColumn(model->rowArray(2),
     2796            int updateStatus = model->factorization()->replaceColumn(model,
     2797                                                                     model->rowArray(2),
     2798                                                                     model->rowArray(1),
    27572799                                                                     possiblePivotKey_,
    2758                                                                    bestAlpha);
     2800                                                                     bestAlpha);
    27592801            returnCode = max (updateStatus, returnCode);
    27602802            incomingColumn = pivotVariable[possiblePivotKey_];
    … …  
    28012843            //if (!alpha)
    28022844            //printf("zero alpha d\n");
    2803             int updateStatus = model->factorization()->replaceColumn(model->rowArray(2),
     2845            int updateStatus = model->factorization()->replaceColumn(model,
     2846                                                                     model->rowArray(2),
     2847                                                                     model->rowArray(3),
    28042848                                                                     iRow, alpha);
    28052849            returnCode = max (updateStatus, returnCode);
    … …  
    28482892            //if (!alpha)
    28492893            //printf("zero alpha e\n");
    2850             int updateStatus = model->factorization()->replaceColumn(model->rowArray(2),
     2894            int updateStatus = model->factorization()->replaceColumn(model,
     2895                                                                     model->rowArray(2),
     2896                                                                     model->rowArray(3),
    28512897                                                                     iRow, alpha);
    28522898            returnCode = max (updateStatus, returnCode);
    … …  
    28682914          //if (!alpha)
    28692915          //printf("zero alpha f\n");
    2870           int updateStatus = model->factorization()->replaceColumn(model->rowArray(2),
     2916          int updateStatus = model->factorization()->replaceColumn(model,
     2917                                                                   model->rowArray(2),
     2918                                                                   model->rowArray(3),
    28712919                                                                   pivotRow, alpha);
    28722920          returnCode = max (updateStatus, returnCode);
    … …  
    28772925      } else {
    28782926        // normal - but might as well do here
    2879         returnCode = model->factorization()->replaceColumn(model->rowArray(2),
     2927        returnCode = model->factorization()->replaceColumn(model,
     2928                                                           model->rowArray(2),
     2929                                                           model->rowArray(1),
    28802930                                                           model->pivotRow(),
    28812931                                                           model->alpha());
    … …  
    29993049      returnCode=synchronize(model,8);
    30003050    }
     3051    break;
     3052  default:
    30013053    break;
    30023054  }
    … …  
    36033655  int pivotRow = model->pivotRow();
    36043656  int iSetIn;
    3605   if (sequenceIn<numberColumns)
     3657  // Correct sequence in
     3658  trueSequenceIn_=sequenceIn;
     3659  if (sequenceIn<numberColumns) {
    36063660    iSetIn = backward_[sequenceIn];
    3607   else if (sequenceIn<numberColumns+numberRows)
     3661  } else if (sequenceIn<numberColumns+numberRows) {
    36083662    iSetIn = -1;
    3609   else
     3663  } else {
    36103664    iSetIn = gubSlackIn_;
     3665    trueSequenceIn_ = numberColumns+numberRows+iSetIn;
     3666  }
    36113667  int iSetOut=-1;
     3668  trueSequenceOut_=sequenceOut;
    36123669  if (sequenceOut<numberColumns) {
    36133670    iSetOut = backward_[sequenceOut];
    … …  
    36203677    else
    36213678      assert(iSetOut == fromIndex_[iExtra]);
     3679    trueSequenceOut_ = numberColumns+numberRows+iSetOut;
    36223680  }
    36233681  if (rhsOffset_) {
    … …  
    40174075  return 0;
    40184076}
    4019 
    4020 
     4077// Switches off dj checking each factorization (for BIG models)
     4078void
     4079ClpGubMatrix::switchOffCheck()
     4080{
     4081  noCheck_=0;
     4082  infeasibilityWeight_=0.0;
     4083}
     4084// Correct sequence in and out to give true value
     4085void
     4086ClpGubMatrix::correctSequence(int & sequenceIn, int & sequenceOut) const
     4087{
     4088  sequenceIn = trueSequenceIn_;
     4089  sequenceOut = trueSequenceOut_;
     4090}

  • trunk/ClpMatrixBase.cpp

    r328 r336  
    2020ClpMatrixBase::ClpMatrixBase () :
    2121  rhsOffset_(NULL),
     22  startFraction_(0.0),
     23  endFraction_(1.0),
     24  savedBestDj_(0.0),
     25  originalWanted_(0),
     26  currentWanted_(0),
     27  savedBestSequence_(-1),
    2228  type_(-1),
    2329  lastRefresh_(-1),
    2430  refreshFrequency_(0),
     31  minimumObjectsScan_(-1),
     32  minimumGoodReducedCosts_(-1),
     33  trueSequenceIn_(-1),
     34  trueSequenceOut_(-1),
    2535  skipDualCheck_(false)
    2636{
    … …  
    3545  skipDualCheck_(rhs.skipDualCheck_)
    3646{ 
     47  startFraction_ = rhs.startFraction_;
     48  endFraction_ = rhs.endFraction_;
     49  savedBestDj_ = rhs.savedBestDj_;
     50  originalWanted_ = rhs.originalWanted_;
     51  currentWanted_ = rhs.currentWanted_;
     52  savedBestSequence_ = rhs.savedBestSequence_;
    3753  lastRefresh_ = rhs.lastRefresh_;
    3854  refreshFrequency_ = rhs.refreshFrequency_;
     55  minimumObjectsScan_ = rhs.minimumObjectsScan_;
     56  minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
     57  trueSequenceIn_ = rhs.trueSequenceIn_;
     58  trueSequenceOut_ = rhs.trueSequenceOut_;
    3959  skipDualCheck_ = rhs.skipDualCheck_;
    4060  int numberRows = rhs.getNumRows();
    … …  
    6989      rhsOffset_=NULL;
    7090    }
     91    startFraction_ = rhs.startFraction_;
     92    endFraction_ = rhs.endFraction_;
     93    savedBestDj_ = rhs.savedBestDj_;
     94    originalWanted_ = rhs.originalWanted_;
     95    currentWanted_ = rhs.currentWanted_;
     96    savedBestSequence_ = rhs.savedBestSequence_;
    7197    lastRefresh_ = rhs.lastRefresh_;
    7298    refreshFrequency_ = rhs.refreshFrequency_;
     99    minimumObjectsScan_ = rhs.minimumObjectsScan_;
     100    minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
     101    trueSequenceIn_ = rhs.trueSequenceIn_;
     102    trueSequenceOut_ = rhs.trueSequenceOut_;
    73103    skipDualCheck_ = rhs.skipDualCheck_;
    74104  }
    … …  
    167197// Partial pricing
    168198void
    169 ClpMatrixBase::partialPricing(ClpSimplex * model, int start, int end,
     199ClpMatrixBase::partialPricing(ClpSimplex * model, double start, double end,
    170200                              int & bestSequence, int & numberWanted)
    171201{
    … …  
    355385ClpMatrixBase::reducedCost(ClpSimplex * model,int sequence) const
    356386{
    357   return model->djRegion()[sequence];
    358 }
    359 
    360 
     387  int numberRows = model->numberRows();
     388  int numberColumns = model->numberColumns();
     389  if (sequence<numberRows+numberColumns)
     390    return model->djRegion()[sequence];
     391  else
     392    return savedBestDj_;
     393}
     394/* Just for debug if odd type matrix.
     395   Returns number of primal infeasibilities.
     396*/
     397int
     398ClpMatrixBase::checkFeasible() const
     399{
     400  return 0;
     401}
     402// Correct sequence in and out to give true value
     403void
     404ClpMatrixBase::correctSequence(int & sequenceIn, int & sequenceOut) const
     405{
     406}
     407
     408

  • trunk/ClpMessage.cpp

    r325 r336  
    9494  {CLP_BARRIER_FEASIBLE,57,2,"Infeasibilities - bound %g , primal %g ,dual %g"},
    9595  {CLP_BARRIER_STEP,58,2,"Steps - primal %g ,dual %g , mu %g"},
     96  {CLP_RIM_SCALE,59,1,"Automatic rim scaling gives objective scale of %g and rhs/bounds scale of %g"},
    9697  {CLP_DUMMY_END,999999,0,""}
    9798};

  • trunk/ClpNetworkMatrix.cpp

    r309 r336  
    772772// Partial pricing
    773773void
    774 ClpNetworkMatrix::partialPricing(ClpSimplex * model, int start, int end,
     774ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    775775                              int & bestSequence, int & numberWanted)
    776776{
     777  numberWanted=currentWanted_;
    777778  int j;
     779  int start = (int) (startFraction*numberColumns_);
     780  int end = min((int) (endFraction*numberColumns_+1),numberColumns_);
    778781  double tolerance=model->currentDualTolerance();
    779782  double * reducedCost = model->djRegion();
    … …  
    895898        value -= duals[iRowP];
    896899      reducedCost[bestSequence] = value;
     900      savedBestSequence_ = bestSequence;
     901      savedBestDj_ = reducedCost[savedBestSequence_];
    897902    }
    898903  } else {
    … …  
    992997      value -= duals[iRowP];
    993998      reducedCost[bestSequence] = value;
    994     }
    995   }
     999      savedBestSequence_ = bestSequence;
     1000      savedBestDj_ = reducedCost[savedBestSequence_];
     1001    }
     1002  }
     1003  currentWanted_=numberWanted;
    9961004}
    9971005// Allow any parts of a created CoinMatrix to be deleted

  • trunk/ClpNonLinearCost.cpp

    r328 r336  
    480480      break;
    481481    case ClpSimplex::isFree:
    482       if (toNearest)
    483         solution[iSequence] = 0.0;
     482      //if (toNearest)
     483      //solution[iSequence] = 0.0;
    484484      break;
    485485    case ClpSimplex::atUpperBound:
    … …  
    577577    //assert (iRange==whichRange_[iSequence]);
    578578  }
     579  //feasibleCost_ /= (model_->rhsScale()*model_->objScale());
    579580}
    580581/* Goes through one bound for each variable.
    … …  
    984985  double value;
    985986  model_->getDblParam(ClpObjOffset,value);
    986   return feasibleCost_*model_->optimizationDirection()-value;
     987  return feasibleCost_*model_->optimizationDirection()/
     988    (model_->objectiveScale()*model_->rhsScale())-value;
    987989}
    988990// Get rid of real costs (just for moment)

  • trunk/ClpPackedMatrix.cpp

    r328 r336  
    2626  : ClpMatrixBase(),
    2727    matrix_(NULL),
    28     zeroElements_(false)
     28    numberActiveColumns_(0),
     29    zeroElements_(false),
     30    hasGaps_(true)
    2931{
    3032  setType(1);
    … …  
    3840{ 
    3941  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
     42  numberActiveColumns_ = rhs.numberActiveColumns_;
    4043  zeroElements_ = rhs.zeroElements_;
     44  hasGaps_ = rhs.hasGaps_;
    4145  int numberRows = getNumRows();
    4246  if (rhs.rhsOffset_&&numberRows) {
    … …  
    5660  matrix_ = rhs;
    5761  zeroElements_ = false;
     62  hasGaps_ = true;
     63  numberActiveColumns_ = matrix_->getNumCols();
    5864  setType(1);
    5965 
    … …  
    6470{ 
    6571  matrix_ = new CoinPackedMatrix(rhs);
     72  numberActiveColumns_ = matrix_->getNumCols();
    6673  zeroElements_ = false;
     74  hasGaps_ = true;
    6775  setType(1);
    6876 
    … …  
    8795    delete matrix_;
    8896    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
     97    numberActiveColumns_ = rhs.numberActiveColumns_;
    8998    zeroElements_ = rhs.zeroElements_;
     99    hasGaps_ = rhs.hasGaps_;
    90100  }
    91101  return *this;
    … …  
    118128  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
    119129                                 numberColumns,whichColumns);
     130  numberActiveColumns_ = matrix_->getNumCols();
    120131  zeroElements_ = rhs.zeroElements_;
     132  hasGaps_ = rhs.hasGaps_;
    121133}
    122134ClpPackedMatrix::ClpPackedMatrix (
    … …  
    128140  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
    129141                                 numberColumns,whichColumns);
     142  numberActiveColumns_ = matrix_->getNumCols();
    130143  zeroElements_ = false;
     144  hasGaps_=true;
    131145  setType(1);
    132146}
    … …  
    140154  copy->matrix_->reverseOrderedCopyOf(*matrix_);
    141155  copy->matrix_->removeGaps();
     156  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
     157  copy->hasGaps_=false;
    142158  return copy;
    143159}
    … …  
    153169  const int * columnLength = matrix_->getVectorLengths();
    154170  const double * elementByColumn = matrix_->getElements();
    155   int numberColumns = matrix_->getNumCols();
    156171  //memset(y,0,matrix_->getNumRows()*sizeof(double));
    157   for (iColumn=0;iColumn<numberColumns;iColumn++) {
     172  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    158173    CoinBigIndex j;
    159174    double value = scalar*x[iColumn];
    … …  
    177192  const int * columnLength = matrix_->getVectorLengths();
    178193  const double * elementByColumn = matrix_->getElements();
    179   int numberColumns = matrix_->getNumCols();
    180   //memset(y,0,numberColumns*sizeof(double));
    181   for (iColumn=0;iColumn<numberColumns;iColumn++) {
     194  //memset(y,0,numberActiveColumns_*sizeof(double));
     195  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    182196    CoinBigIndex j;
    183197    double value=0.0;
    … …  
    203217    const int * columnLength = matrix_->getVectorLengths();
    204218    const double * elementByColumn = matrix_->getElements();
    205     int numberColumns = matrix_->getNumCols();
    206219    //memset(y,0,matrix_->getNumRows()*sizeof(double));
    207     for (iColumn=0;iColumn<numberColumns;iColumn++) {
     220    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    208221      CoinBigIndex j;
    209222      double value = x[iColumn];
    … …  
    239252    const int * columnLength = matrix_->getVectorLengths();
    240253    const double * elementByColumn = matrix_->getElements();
    241     int numberColumns = matrix_->getNumCols();
    242     //memset(y,0,numberColumns*sizeof(double));
    243     for (iColumn=0;iColumn<numberColumns;iColumn++) {
     254    //memset(y,0,numberActiveColumns_*sizeof(double));
     255    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    244256      CoinBigIndex j;
    245257      double value=0.0;
    … …  
    278290  double factor = 0.3;
    279291  // We may not want to do by row if there may be cache problems
    280   int numberColumns = model->numberColumns();
    281292  // It would be nice to find L2 cache size - for moment 512K
    282293  // Be slightly optimistic
    283   if (numberColumns*sizeof(double)>1000000) {
    284     if (numberRows*10<numberColumns)
     294  if (numberActiveColumns_*sizeof(double)>1000000) {
     295    if (numberRows*10<numberActiveColumns_)
    285296      factor=0.1;
    286     else if (numberRows*4<numberColumns)
     297    else if (numberRows*4<numberActiveColumns_)
    287298      factor=0.15;
    288     else if (numberRows*2<numberColumns)
     299    else if (numberRows*2<numberActiveColumns_)
    289300      factor=0.2;
    290301    //if (model->numberIterations()%50==0)
    291302    //printf("%d nonzero\n",numberInRowArray);
    292303  }
     304  assert (!y->getNumElements());
    293305  if (numberInRowArray>factor*numberRows||!rowCopy) {
    294306    // do by column
     307    // If no gaps - can do a bit faster
     308    if (!hasGaps_&&true) {
     309      transposeTimesByColumn( model,  scalar,
     310                              rowArray, y, columnArray);
     311      return;
     312    }
    295313    int iColumn;
    296314    // get matrix data pointers
    … …  
    300318    const double * elementByColumn = matrix_->getElements();
    301319    const double * rowScale = model->rowScale();
    302     int numberColumns = model->numberColumns();
    303     if (!y->getNumElements()) {
    304       if (packed) {
    305         // need to expand pi into y
    306         assert(y->capacity()>=numberRows);
    307         double * piOld = pi;
    308         pi = y->denseVector();
    309         const int * whichRow = rowArray->getIndices();
    310         int i;
    311         if (!rowScale) {
    312           // modify pi so can collapse to one loop
    313           for (i=0;i<numberInRowArray;i++) {
    314             int iRow = whichRow[i];
    315             pi[iRow]=scalar*piOld[i];
    316           }
    317           for (iColumn=0;iColumn<numberColumns;iColumn++) {
     320    if (packed) {
     321      // need to expand pi into y
     322      assert(y->capacity()>=numberRows);
     323      double * piOld = pi;
     324      pi = y->denseVector();
     325      const int * whichRow = rowArray->getIndices();
     326      int i;
     327      if (!rowScale) {
     328        // modify pi so can collapse to one loop
     329        for (i=0;i<numberInRowArray;i++) {
     330          int iRow = whichRow[i];
     331          pi[iRow]=scalar*piOld[i];
     332        }
     333        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     334          double value = 0.0;
     335          CoinBigIndex j;
     336          for (j=columnStart[iColumn];
     337               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     338            int iRow = row[j];
     339            value += pi[iRow]*elementByColumn[j];
     340          }
     341          if (fabs(value)>zeroTolerance) {
     342            array[numberNonZero]=value;
     343            index[numberNonZero++]=iColumn;
     344          }
     345        }
     346      } else {
     347        // scaled
     348        // modify pi so can collapse to one loop
     349        for (i=0;i<numberInRowArray;i++) {
     350          int iRow = whichRow[i];
     351          pi[iRow]=scalar*piOld[i]*rowScale[iRow];
     352        }
     353        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     354          double value = 0.0;
     355          CoinBigIndex j;
     356          const double * columnScale = model->columnScale();
     357          for (j=columnStart[iColumn];
     358               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     359            int iRow = row[j];
     360            value += pi[iRow]*elementByColumn[j];
     361          }
     362          value *= columnScale[iColumn];
     363          if (fabs(value)>zeroTolerance) {
     364            array[numberNonZero]=value;
     365            index[numberNonZero++]=iColumn;
     366          }
     367        }
     368      }
     369      // zero out
     370      for (i=0;i<numberInRowArray;i++) {
     371        int iRow = whichRow[i];
     372        pi[iRow]=0.0;
     373      }
     374    } else {
     375      if (!rowScale) {
     376        if (scalar==-1.0) {
     377          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    318378            double value = 0.0;
    319379            CoinBigIndex j;
    … …  
    324384            }
    325385            if (fabs(value)>zeroTolerance) {
    326               array[numberNonZero]=value;
    327386              index[numberNonZero++]=iColumn;
     387              array[iColumn]=-value;
     388            }
     389          }
     390        } else if (scalar==1.0) {
     391          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     392            double value = 0.0;
     393            CoinBigIndex j;
     394            for (j=columnStart[iColumn];
     395                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     396              int iRow = row[j];
     397              value += pi[iRow]*elementByColumn[j];
     398            }
     399            if (fabs(value)>zeroTolerance) {
     400              index[numberNonZero++]=iColumn;
     401              array[iColumn]=value;
    328402            }
    329403          }
    330404        } else {
    331           // scaled
    332           // modify pi so can collapse to one loop
    333           for (i=0;i<numberInRowArray;i++) {
    334             int iRow = whichRow[i];
    335             pi[iRow]=scalar*piOld[i]*rowScale[iRow];
    336           }
    337           for (iColumn=0;iColumn<numberColumns;iColumn++) {
     405          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     406            double value = 0.0;
     407            CoinBigIndex j;
     408            for (j=columnStart[iColumn];
     409                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     410              int iRow = row[j];
     411              value += pi[iRow]*elementByColumn[j];
     412            }
     413            value *= scalar;
     414            if (fabs(value)>zeroTolerance) {
     415              index[numberNonZero++]=iColumn;
     416              array[iColumn]=value;
     417            }
     418          }
     419        }
     420      } else {
     421        // scaled
     422        if (scalar==-1.0) {
     423          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    338424            double value = 0.0;
    339425            CoinBigIndex j;
    … …  
    342428                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
    343429              int iRow = row[j];
    344               value += pi[iRow]*elementByColumn[j];
     430              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    345431            }
    346432            value *= columnScale[iColumn];
    347433            if (fabs(value)>zeroTolerance) {
    348               array[numberNonZero]=value;
    349434              index[numberNonZero++]=iColumn;
    350             }
    351           }
    352         }
    353         // zero out
    354         for (i=0;i<numberInRowArray;i++) {
    355           int iRow = whichRow[i];
    356           pi[iRow]=0.0;
    357         }
    358       } else {
    359         if (!rowScale) {
    360           if (scalar==-1.0) {
    361             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    362               double value = 0.0;
    363               CoinBigIndex j;
    364               for (j=columnStart[iColumn];
    365                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    366                 int iRow = row[j];
    367                 value += pi[iRow]*elementByColumn[j];
    368               }
    369               if (fabs(value)>zeroTolerance) {
    370                 index[numberNonZero++]=iColumn;
    371                 array[iColumn]=-value;
    372               }
    373             }
    374           } else if (scalar==1.0) {
    375             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    376               double value = 0.0;
    377               CoinBigIndex j;
    378               for (j=columnStart[iColumn];
    379                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    380                 int iRow = row[j];
    381                 value += pi[iRow]*elementByColumn[j];
    382               }
    383               if (fabs(value)>zeroTolerance) {
    384                 index[numberNonZero++]=iColumn;
    385                 array[iColumn]=value;
    386               }
    387             }
    388           } else {
    389             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    390               double value = 0.0;
    391               CoinBigIndex j;
    392               for (j=columnStart[iColumn];
    393                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    394                 int iRow = row[j];
    395                 value += pi[iRow]*elementByColumn[j];
    396               }
    397               value *= scalar;
    398               if (fabs(value)>zeroTolerance) {
    399                 index[numberNonZero++]=iColumn;
    400                 array[iColumn]=value;
    401               }
     435              array[iColumn]=-value;
     436            }
     437          }
     438        } else if (scalar==1.0) {
     439          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     440            double value = 0.0;
     441            CoinBigIndex j;
     442            const double * columnScale = model->columnScale();
     443            for (j=columnStart[iColumn];
     444                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     445              int iRow = row[j];
     446              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     447            }
     448            value *= columnScale[iColumn];
     449            if (fabs(value)>zeroTolerance) {
     450              index[numberNonZero++]=iColumn;
     451              array[iColumn]=value;
    402452            }
    403453          }
    404454        } else {
    405           // scaled
    406           if (scalar==-1.0) {
    407             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    408               double value = 0.0;
    409               CoinBigIndex j;
    410               const double * columnScale = model->columnScale();
    411               for (j=columnStart[iColumn];
    412                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    413                 int iRow = row[j];
    414                 value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    415               }
    416               value *= columnScale[iColumn];
    417               if (fabs(value)>zeroTolerance) {
    418                 index[numberNonZero++]=iColumn;
    419                 array[iColumn]=-value;
    420               }
    421             }
    422           } else if (scalar==1.0) {
    423             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    424               double value = 0.0;
    425               CoinBigIndex j;
    426               const double * columnScale = model->columnScale();
    427               for (j=columnStart[iColumn];
    428                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    429                 int iRow = row[j];
    430                 value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    431               }
    432               value *= columnScale[iColumn];
    433               if (fabs(value)>zeroTolerance) {
    434                 index[numberNonZero++]=iColumn;
    435                 array[iColumn]=value;
    436               }
    437             }
    438           } else {
    439             for (iColumn=0;iColumn<numberColumns;iColumn++) {
    440               double value = 0.0;
    441               CoinBigIndex j;
    442               const double * columnScale = model->columnScale();
    443               for (j=columnStart[iColumn];
    444                    j<columnStart[iColumn]+columnLength[iColumn];j++) {
    445                 int iRow = row[j];
    446                 value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    447               }
    448               value *= scalar*columnScale[iColumn];
    449               if (fabs(value)>zeroTolerance) {
    450                 index[numberNonZero++]=iColumn;
    451                 array[iColumn]=value;
    452               }
    453             }
    454           }
    455         }
    456       }
    457     } else {
    458       assert(!packed);
    459       double * markVector = y->denseVector(); // not empty
    460       if (!rowScale) {
    461         for (iColumn=0;iColumn<numberColumns;iColumn++) {
    462           double value = markVector[iColumn];
    463           markVector[iColumn]=0.0;
    464           double value2 = 0.0;
    465           CoinBigIndex j;
    466           for (j=columnStart[iColumn];
    467                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    468             int iRow = row[j];
    469             value2 += pi[iRow]*elementByColumn[j];
    470           }
    471           value += value2*scalar;
    472           if (fabs(value)>zeroTolerance) {
    473             index[numberNonZero++]=iColumn;
    474             array[iColumn]=value;
    475           }
    476         }
    477       } else {
    478         // scaled
    479         for (iColumn=0;iColumn<numberColumns;iColumn++) {
    480           double value = markVector[iColumn];
    481           markVector[iColumn]=0.0;
    482           CoinBigIndex j;
    483           const double * columnScale = model->columnScale();
    484           double value2 = 0.0;
    485           for (j=columnStart[iColumn];
    486                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    487             int iRow = row[j];
    488             value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    489           }
    490           value += value2*scalar*columnScale[iColumn];
    491           if (fabs(value)>zeroTolerance) {
    492             index[numberNonZero++]=iColumn;
    493             array[iColumn]=value;
     455          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     456            double value = 0.0;
     457            CoinBigIndex j;
     458            const double * columnScale = model->columnScale();
     459            for (j=columnStart[iColumn];
     460                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     461              int iRow = row[j];
     462              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     463            }
     464            value *= scalar*columnScale[iColumn];
     465            if (fabs(value)>zeroTolerance) {
     466              index[numberNonZero++]=iColumn;
     467              array[iColumn]=value;
     468            }
    494469          }
    495470        }
    … …  
    521496  }
    522497}
     498/* Return <code>x * scalar * A + y</code> in <code>z</code>.
     499   Note - If x packed mode - then z packed mode
     500   This does by column and knows no gaps
     501   Squashes small elements and knows about ClpSimplex */
     502void
     503ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
     504                                        const CoinIndexedVector * rowArray,
     505                                        CoinIndexedVector * y,
     506                                        CoinIndexedVector * columnArray) const
     507{
     508  double * pi = rowArray->denseVector();
     509  int numberNonZero=0;
     510  int * index = columnArray->getIndices();
     511  double * array = columnArray->denseVector();
     512  int numberInRowArray = rowArray->getNumElements();
     513  // maybe I need one in OsiSimplex
     514  double zeroTolerance = model->factorization()->zeroTolerance();
     515  bool packed = rowArray->packedMode();
     516  // do by column
     517  int iColumn;
     518  // get matrix data pointers
     519  const int * row = matrix_->getIndices();
     520  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     521  const double * elementByColumn = matrix_->getElements();
     522  const double * rowScale = model->rowScale();
     523  assert (!y->getNumElements());
     524  assert (numberActiveColumns_>0);
     525  if (packed) {
     526    // need to expand pi into y
     527    assert(y->capacity()>=model->numberRows());
     528    double * piOld = pi;
     529    pi = y->denseVector();
     530    const int * whichRow = rowArray->getIndices();
     531    int i;
     532    if (!rowScale) {
     533      // modify pi so can collapse to one loop
     534      for (i=0;i<numberInRowArray;i++) {
     535        int iRow = whichRow[i];
     536        pi[iRow]=scalar*piOld[i];
     537      }
     538      double value = 0.0;
     539      CoinBigIndex j;
     540      CoinBigIndex end = columnStart[1];
     541      for (j=columnStart[0]; j<end;j++) {
     542        int iRow = row[j];
     543        value += pi[iRow]*elementByColumn[j];
     544      }
     545      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     546        CoinBigIndex start = end;
     547        end = columnStart[iColumn+2];
     548        if (fabs(value)>zeroTolerance) {
     549          array[numberNonZero]=value;
     550          index[numberNonZero++]=iColumn;
     551        }
     552        value = 0.0;
     553        for (j=start; j<end;j++) {
     554          int iRow = row[j];
     555          value += pi[iRow]*elementByColumn[j];
     556        }
     557      }
     558      if (fabs(value)>zeroTolerance) {
     559        array[numberNonZero]=value;
     560        index[numberNonZero++]=iColumn;
     561      }
     562    } else {
     563      // scaled
     564      // modify pi so can collapse to one loop
     565      for (i=0;i<numberInRowArray;i++) {
     566        int iRow = whichRow[i];
     567        pi[iRow]=scalar*piOld[i]*rowScale[iRow];
     568      }
     569      const double * columnScale = model->columnScale();
     570      double value = 0.0;
     571      double scale=columnScale[0];
     572      CoinBigIndex j;
     573      CoinBigIndex end = columnStart[1];
     574      for (j=columnStart[0]; j<end;j++) {
     575        int iRow = row[j];
     576        value += pi[iRow]*elementByColumn[j];
     577      }
     578      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     579        value *= scale;
     580        CoinBigIndex start = end;
     581        scale = columnScale[iColumn+1];
     582        end = columnStart[iColumn+2];
     583        if (fabs(value)>zeroTolerance) {
     584          array[numberNonZero]=value;
     585          index[numberNonZero++]=iColumn;
     586        }
     587        value = 0.0;
     588        for (j=start; j<end;j++) {
     589          int iRow = row[j];
     590          value += pi[iRow]*elementByColumn[j];
     591        }
     592      }
     593      value *= scale;
     594      if (fabs(value)>zeroTolerance) {
     595        array[numberNonZero]=value;
     596        index[numberNonZero++]=iColumn;
     597      }
     598    }
     599    // zero out
     600    for (i=0;i<numberInRowArray;i++) {
     601      int iRow = whichRow[i];
     602      pi[iRow]=0.0;
     603    }
     604  } else {
     605    if (!rowScale) {
     606      if (scalar==-1.0) {
     607        double value = 0.0;
     608        CoinBigIndex j;
     609        CoinBigIndex end = columnStart[1];
     610        for (j=columnStart[0]; j<end;j++) {
     611          int iRow = row[j];
     612          value += pi[iRow]*elementByColumn[j];
     613        }
     614        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     615          CoinBigIndex start = end;
     616          end = columnStart[iColumn+2];
     617          if (fabs(value)>zeroTolerance) {
     618            array[numberNonZero]=-value;
     619            index[numberNonZero++]=iColumn;
     620          }
     621          value = 0.0;
     622          for (j=start; j<end;j++) {
     623            int iRow = row[j];
     624            value += pi[iRow]*elementByColumn[j];
     625          }
     626        }
     627        if (fabs(value)>zeroTolerance) {
     628          array[numberNonZero]=-value;
     629          index[numberNonZero++]=iColumn;
     630        }
     631      } else if (scalar==1.0) {
     632        double value = 0.0;
     633        CoinBigIndex j;
     634        CoinBigIndex end = columnStart[1];
     635        for (j=columnStart[0]; j<end;j++) {
     636          int iRow = row[j];
     637          value += pi[iRow]*elementByColumn[j];
     638        }
     639        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     640          CoinBigIndex start = end;
     641          end = columnStart[iColumn+2];
     642          if (fabs(value)>zeroTolerance) {
     643            array[numberNonZero]=value;
     644            index[numberNonZero++]=iColumn;
     645          }
     646          value = 0.0;
     647          for (j=start; j<end;j++) {
     648            int iRow = row[j];
     649            value += pi[iRow]*elementByColumn[j];
     650          }
     651        }
     652        if (fabs(value)>zeroTolerance) {
     653          array[numberNonZero]=value;
     654          index[numberNonZero++]=iColumn;
     655        }
     656      } else {
     657        double value = 0.0;
     658        CoinBigIndex j;
     659        CoinBigIndex end = columnStart[1];
     660        for (j=columnStart[0]; j<end;j++) {
     661          int iRow = row[j];
     662          value += pi[iRow]*elementByColumn[j];
     663        }
     664        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     665          value *= scalar;
     666          CoinBigIndex start = end;
     667          end = columnStart[iColumn+2];
     668          if (fabs(value)>zeroTolerance) {
     669            array[numberNonZero]=value;
     670            index[numberNonZero++]=iColumn;
     671          }
     672          value = 0.0;
     673          for (j=start; j<end;j++) {
     674            int iRow = row[j];
     675            value += pi[iRow]*elementByColumn[j];
     676          }
     677        }
     678        value *= scalar;
     679        if (fabs(value)>zeroTolerance) {
     680          array[numberNonZero]=value;
     681          index[numberNonZero++]=iColumn;
     682        }
     683      }
     684    } else {
     685      // scaled
     686      if (scalar==-1.0) {
     687        const double * columnScale = model->columnScale();
     688        double value = 0.0;
     689        double scale=columnScale[0];
     690        CoinBigIndex j;
     691        CoinBigIndex end = columnStart[1];
     692        for (j=columnStart[0]; j<end;j++) {
     693          int iRow = row[j];
     694          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     695        }
     696        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     697          value *= scale;
     698          CoinBigIndex start = end;
     699          end = columnStart[iColumn+2];
     700          scale = columnScale[iColumn+1];
     701          if (fabs(value)>zeroTolerance) {
     702            array[numberNonZero]=-value;
     703            index[numberNonZero++]=iColumn;
     704          }
     705          value = 0.0;
     706          for (j=start; j<end;j++) {
     707            int iRow = row[j];
     708            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     709          }
     710        }
     711        value *= scale;
     712        if (fabs(value)>zeroTolerance) {
     713          array[numberNonZero]=-value;
     714          index[numberNonZero++]=iColumn;
     715        }
     716      } else if (scalar==1.0) {
     717        const double * columnScale = model->columnScale();
     718        double value = 0.0;
     719        double scale=columnScale[0];
     720        CoinBigIndex j;
     721        CoinBigIndex end = columnStart[1];
     722        for (j=columnStart[0]; j<end;j++) {
     723          int iRow = row[j];
     724          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     725        }
     726        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     727          value *= scale;
     728          CoinBigIndex start = end;
     729          end = columnStart[iColumn+2];
     730          scale = columnScale[iColumn+1];
     731          if (fabs(value)>zeroTolerance) {
     732            array[numberNonZero]=value;
     733            index[numberNonZero++]=iColumn;
     734          }
     735          value = 0.0;
     736          for (j=start; j<end;j++) {
     737            int iRow = row[j];
     738            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     739          }
     740        }
     741        value *= scale;
     742        if (fabs(value)>zeroTolerance) {
     743          array[numberNonZero]=value;
     744          index[numberNonZero++]=iColumn;
     745        }
     746      } else {
     747        const double * columnScale = model->columnScale();
     748        double value = 0.0;
     749        double scale=columnScale[0]*scalar;
     750        CoinBigIndex j;
     751        CoinBigIndex end = columnStart[1];
     752        for (j=columnStart[0]; j<end;j++) {
     753          int iRow = row[j];
     754          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     755        }
     756        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     757          value *= scale;
     758          CoinBigIndex start = end;
     759          end = columnStart[iColumn+2];
     760          scale = columnScale[iColumn+1]*scalar;
     761          if (fabs(value)>zeroTolerance) {
     762            array[numberNonZero]=value;
     763            index[numberNonZero++]=iColumn;
     764          }
     765          value = 0.0;
     766          for (j=start; j<end;j++) {
     767            int iRow = row[j];
     768            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     769          }
     770        }
     771        value *= scale;
     772        if (fabs(value)>zeroTolerance) {
     773          array[numberNonZero]=value;
     774          index[numberNonZero++]=iColumn;
     775        }
     776      }
     777    }
     778  }
     779  columnArray->setNumElements(numberNonZero);
     780  y->setNumElements(0);
     781  if (packed)
     782    columnArray->setPackedMode(true);
     783}
     784/* Return <code>x * scalar * A in <code>z</code>.
     785   Note - this version when x packed mode and so will be z
     786   This does by column and knows no gaps and knows y empty
     787   Squashes small elements and knows about ClpSimplex */
     788void
     789ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
     790                                        const CoinIndexedVector * rowArray,
     791                                        CoinIndexedVector * columnArray) const
     792{
     793}
    523794/* Return <code>x * A + y</code> in <code>z</code>.
    524795        Squashes small elements and knows about ClpSimplex */
    … …  
    542813  const int * whichRow = rowArray->getIndices();
    543814  bool packed = rowArray->packedMode();
    544   if (numberInRowArray>2||y->getNumElements()) {
     815  if (numberInRowArray>2) {
    545816    // do by rows
    546817    // ** Row copy is already scaled
    547818    int iRow;
    548     int * mark = y->getIndices();
    549     int numberOriginal=y->getNumElements();
    550819    int i;
     820    int numberOriginal = 0;
    551821    if (packed) {
    552       assert(!numberOriginal);
    553822      numberNonZero=0;
    554823      // and set up mark as char array
    … …  
    556825      double * array2 = y->denseVector();
    557826#ifdef CLP_DEBUG
    558       int numberColumns = model->numberColumns();
    559       for (i=0;i<numberColumns;i++) {
     827      for (i=0;i<numberActiveColumns_;i++) {
    560828        assert(!marked[i]);
    561829        assert(!array2[i]);
    … …  
    591859      }
    592860    } else {
    593       double * markVector = y->denseVector(); // probably empty .. but
    594       for (i=0;i<numberOriginal;i++) {
    595         int iColumn = mark[i];
    596         index[i]=iColumn;
    597         array[iColumn]=markVector[iColumn];
    598         markVector[iColumn]=0.0;
    599       }
    600       numberNonZero=numberOriginal;
     861      double * markVector = y->denseVector();
     862      numberNonZero=0;
    601863      // and set up mark as char array
    602864      char * marked = (char *) markVector;
    … …  
    7881050  assert (!rowArray->packedMode());
    7891051  columnArray->setPacked();
    790   if (!rowScale) {
    791     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    792       int iColumn = which[jColumn];
     1052  if (!hasGaps_&&numberToDo>5) {
     1053    // no gaps and a reasonable number
     1054    if (!rowScale) {
     1055      int iColumn = which[0];
    7931056      double value = 0.0;
    7941057      CoinBigIndex j;
    7951058      for (j=columnStart[iColumn];
    796            j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1059           j<columnStart[iColumn+1];j++) {
    7971060        int iRow = row[j];
    7981061        value += pi[iRow]*elementByColumn[j];
    7991062      }
     1063      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
     1064        int iColumn = which[jColumn+1];
     1065        CoinBigIndex start = columnStart[iColumn];
     1066        CoinBigIndex end = columnStart[iColumn+1];
     1067        array[jColumn]=value;
     1068        value = 0.0;
     1069        for (j=start; j<end;j++) {
     1070          int iRow = row[j];
     1071          value += pi[iRow]*elementByColumn[j];
     1072        }
     1073      }
    8001074      array[jColumn]=value;
    801     }
    802   } else {
    803     // scaled
    804     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    805       int iColumn = which[jColumn];
     1075    } else {
     1076      // scaled
     1077      const double * columnScale = model->columnScale();
     1078      int iColumn = which[0];
    8061079      double value = 0.0;
     1080      double scale=columnScale[iColumn];
    8071081      CoinBigIndex j;
    808       const double * columnScale = model->columnScale();
    8091082      for (j=columnStart[iColumn];
    810            j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1083           j<columnStart[iColumn+1];j++) {
    8111084        int iRow = row[j];
    8121085        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    8131086      }
    814       value *= columnScale[iColumn];
     1087      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
     1088        int iColumn = which[jColumn+1];
     1089        value *= scale;
     1090        scale = columnScale[iColumn];
     1091        CoinBigIndex start = columnStart[iColumn];
     1092        CoinBigIndex end = columnStart[iColumn+1];
     1093        array[jColumn]=value;
     1094        value = 0.0;
     1095        for (j=start; j<end;j++) {
     1096          int iRow = row[j];
     1097          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     1098        }
     1099      }
     1100      value *= scale;
    8151101      array[jColumn]=value;
     1102    }
     1103  } else {
     1104    // gaps
     1105    if (!rowScale) {
     1106      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     1107        int iColumn = which[jColumn];
     1108        double value = 0.0;
     1109        CoinBigIndex j;
     1110        for (j=columnStart[iColumn];
     1111             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1112          int iRow = row[j];
     1113          value += pi[iRow]*elementByColumn[j];
     1114        }
     1115        array[jColumn]=value;
     1116      }
     1117    } else {
     1118      // scaled
     1119      const double * columnScale = model->columnScale();
     1120      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     1121        int iColumn = which[jColumn];
     1122        double value = 0.0;
     1123        CoinBigIndex j;
     1124        for (j=columnStart[iColumn];
     1125             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1126          int iRow = row[j];
     1127          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     1128        }
     1129        value *= columnScale[iColumn];
     1130        array[jColumn]=value;
     1131      }
    8161132    }
    8171133  }
    … …  
    9211237{
    9221238  int numberRows = model->numberRows();
    923   int numberColumns = model->numberColumns();
     1239  int numberColumns = matrix_->getNumCols();
    9241240  // If empty - return as sanityCheck will trap
    9251241  if (!numberRows||!numberColumns)
    … …  
    11701486      }
    11711487    }
    1172 #define RANDOMIZE
     1488    //#define RANDOMIZE
    11731489#ifdef RANDOMIZE
    11741490    // randomize by up to 10%
    … …  
    11811497    if (overallSmallest<1.0e-1)
    11821498      overallLargest = 1.0/sqrt(overallSmallest);
    1183     overallLargest = min(1000.0,overallLargest);
     1499    overallLargest = min(100.0,overallLargest);
    11841500    overallSmallest=1.0e50;
    11851501    for (iColumn=0;iColumn<numberColumns;iColumn++) {
    … …  
    14451761  const int * columnLength = matrix_->getVectorLengths();
    14461762  const double * elementByColumn = matrix_->getElements();
     1763  int numberRows = matrix_->getNumRows();
    14471764  int numberColumns = matrix_->getNumCols();
    1448   int numberRows = matrix_->getNumRows();
    14491765  int * mark = new int [numberRows];
    14501766  int i;
     1767  // Say no gaps
     1768  hasGaps_=false;
    14511769  for (i=0;i<numberRows;i++)
    14521770    mark[i]=-1;
    14531771  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    14541772    CoinBigIndex j;
    1455     for (j=columnStart[iColumn];
    1456          j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1773    CoinBigIndex start = columnStart[iColumn];
     1774    CoinBigIndex end = start + columnLength[iColumn];
     1775    if (end!=columnStart[iColumn+1])
     1776      hasGaps_=true;
     1777    for (j=start;j<end;j++) {
    14571778      double value = fabs(elementByColumn[j]);
    14581779      int iRow = row[j];
    … …  
    15271848{
    15281849  int numberRows = model->numberRows();
    1529   int numberColumns =model->numberColumns();
     1850  int numberColumns = matrix_->getNumCols();
    15301851  int number = numberRows+numberColumns;
    15311852  CoinBigIndex * weights = new CoinBigIndex[number];
    … …  
    15601881  smallestPositive=COIN_DBL_MAX;
    15611882  largestPositive=0.0;
    1562   int numberColumns = matrix_->getNumCols();
    15631883  // get matrix data pointers
    15641884  const double * elementByColumn = matrix_->getElements();
    15651885  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    15661886  const int * columnLength = matrix_->getVectorLengths();
     1887  int numberColumns = matrix_->getNumCols();
    15671888  int i;
    15681889  for (i=0;i<numberColumns;i++) {
    … …  
    15881909// Partial pricing
    15891910void
    1590 ClpPackedMatrix::partialPricing(ClpSimplex * model, int start, int end,
     1911ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    15911912                              int & bestSequence, int & numberWanted)
    15921913{
     1914  numberWanted=currentWanted_;
     1915  int start = (int) (startFraction*numberActiveColumns_);
     1916  int end = min((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
    15931917  const double * element =matrix_->getElements();
    15941918  const int * row = matrix_->getIndices();
    … …  
    16101934  int sequenceOut = model->sequenceOut();
    16111935  int saveSequence = bestSequence;
     1936  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_;
     1937  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
    16121938  if (rowScale) {
    16131939    // scaled
    … …  
    16962022        }
    16972023      }
     2024      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
     2025        // give up
     2026        break;
     2027      }
    16982028      if (!numberWanted)
    16992029        break;
    … …  
    17092039      }
    17102040      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
     2041      savedBestSequence_ = bestSequence;
     2042      savedBestDj_ = reducedCost[savedBestSequence_];
    17112043    }
    17122044  } else {
    … …  
    17932125        }
    17942126      }
     2127      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
     2128        // give up
     2129        break;
     2130      }
    17952131      if (!numberWanted)
    17962132        break;
    … …  
    18052141      }
    18062142      reducedCost[bestSequence] = value;
    1807     }
    1808   }
     2143      savedBestSequence_ = bestSequence;
     2144      savedBestDj_ = reducedCost[savedBestSequence_];
     2145    }
     2146  }
     2147  currentWanted_=numberWanted;
    18092148}
    18102149// Sets up an effective RHS
    … …  
    18172156  rhsOffset(model,true);
    18182157}
     2158// makes sure active columns correct
     2159int
     2160ClpPackedMatrix::refresh(ClpSimplex * model)
     2161{
     2162  numberActiveColumns_ = matrix_->getNumCols();
     2163  return 0;
     2164}
    18192165

  • trunk/ClpPlusMinusOneMatrix.cpp

    r309 r336  
    571571  const int * whichRow = rowArray->getIndices();
    572572  bool packed = rowArray->packedMode();
    573   if (numberInRowArray>2||y->getNumElements()) {
     573  if (numberInRowArray>2) {
    574574    // do by rows
    575575    int iRow;
    576576    double * markVector = y->denseVector(); // probably empty .. but
    577     int * mark = y->getIndices();
    578     int numberOriginal=y->getNumElements();
     577    int numberOriginal=0;
    579578    int i;
    580579    if (packed) {
    581       assert(!numberOriginal);
    582580      numberNonZero=0;
    583581      // and set up mark as char array
    … …  
    628626      }
    629627    } else {     
    630       for (i=0;i<numberOriginal;i++) {
    631         int iColumn = mark[i];
    632         index[i]=iColumn;
    633         array[iColumn]=markVector[iColumn];
    634         markVector[iColumn]=0.0;
    635       }
    636       numberNonZero=numberOriginal;
     628      numberNonZero=0;
    637629      // and set up mark as char array
    638630      char * marked = (char *) markVector;
    … …  
    14231415// Partial pricing
    14241416void
    1425 ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, int start, int end,
     1417ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    14261418                              int & bestSequence, int & numberWanted)
    14271419{
     1420  numberWanted=currentWanted_;
     1421  int start = (int) (startFraction*numberColumns_);
     1422  int end = min((int) (endFraction*numberColumns_+1),numberColumns_);
    14281423  CoinBigIndex j;
    14291424  double tolerance=model->currentDualTolerance();
    … …  
    15471542    }
    15481543    reducedCost[bestSequence] = value;
    1549   }
     1544    savedBestSequence_ = bestSequence;
     1545    savedBestDj_ = reducedCost[savedBestSequence_];
     1546  }
     1547  currentWanted_=numberWanted;
    15501548}
    15511549// Allow any parts of a created CoinMatrix to be deleted

  • trunk/ClpPredictorCorrector.cpp

    r328 r336  
    99
    1010#include "CoinPragma.hpp"
    11 
     11static int xxxxxx=-1;
    1212#include <math.h>
    1313
    … …  
    9999  double * savePrimal=NULL;
    100100  int numberTotal = numberRows_+numberColumns_;
     101  bool tryJustPredictor=false;
    101102  while (problemStatus_<0) {
     103    if (numberIterations_==xxxxxx) {
     104      FILE *fp = fopen("yy.yy","wb");
     105      if (fp) {
     106        int nrow=numberRows_;
     107        int ncol=numberColumns_;
     108        // and flip
     109        int i;
     110        for (i=0;i<nrow;i++) {
     111          solution_[ncol+i]= -solution_[ncol+i];
     112          double value;
     113          value=lowerSlack_[ncol+i];
     114          lowerSlack_[ncol+i]=upperSlack_[ncol+i];
     115          upperSlack_[ncol+i]=value;
     116          value=zVec_[ncol+i];
     117          zVec_[ncol+i]=wVec_[ncol+i];
     118          wVec_[ncol+i]=value;
     119        }
     120        fwrite(dual_,sizeof(double),nrow,fp);
     121        fwrite(errorRegion_,sizeof(double),nrow,fp);
     122        fwrite(rhsFixRegion_,sizeof(double),nrow,fp);
     123        fwrite(solution_+ncol,sizeof(double),nrow,fp);
     124        fwrite(diagonal_+ncol,sizeof(double),nrow,fp);
     125        fwrite(wVec_+ncol,sizeof(double),nrow,fp);
     126        fwrite(zVec_+ncol,sizeof(double),nrow,fp);
     127        fwrite(upperSlack_+ncol,sizeof(double),nrow,fp);
     128        fwrite(lowerSlack_+ncol,sizeof(double),nrow,fp);
     129        fwrite(solution_,sizeof(double),ncol,fp);
     130        fwrite(diagonal_,sizeof(double),ncol,fp);
     131        fwrite(wVec_,sizeof(double),ncol,fp);
     132        fwrite(zVec_,sizeof(double),ncol,fp);
     133        fwrite(upperSlack_,sizeof(double),ncol,fp);
     134        fwrite(lowerSlack_,sizeof(double),ncol,fp);
     135        fclose(fp);
     136        actualDualStep_=0.0;
     137        actualPrimalStep_=0.0;
     138      } else {
     139        printf("no file\n");
     140      }
     141    }
     142    //#define FULL_DEBUG
     143#ifdef FULL_DEBUG
     144    {
     145      int i;
     146      printf("row    pi          artvec       rhsfx\n");
     147      for (i=0;i<numberRows_;i++) {
     148        printf("%d %g %g %g\n",i,dual_[i],errorRegion_[i],rhsFixRegion_[i]);
     149      }
     150      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
     151      for (i=0;i<numberColumns_+numberRows_;i++) {
     152        printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],wVec_[i],
     153               zVec_[i],upperSlack_[i],lowerSlack_[i]);
     154      }
     155    }
     156#endif
    102157    complementarityGap_=complementarityGap(numberComplementarityPairs_,0);
    103158      handler_->message(CLP_BARRIER_ITERATION,messages_)
    … …  
    187242        <<CoinMessageEol;
    188243    }
     244    //if (complementarityGap_>=0.98*lastComplementarityGap) {
     245    //tryJustPredictor=true;
     246    //printf("trying just predictor\n");
     247    //}
    189248    if (complementarityGap_>=1.05*lastComplementarityGap) {
    190249      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
    … …  
    340399            double part2=nextGap/complementarityGap_;
    341400            mu_=part1*part2*part2;
     401#if 0
     402            double papermu =complementarityGap_/numberComplementarityPairs_;
     403            double affmu = nextGap/nextNumber;
     404            double sigma = pow(affmu/papermu,3);
     405            printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
     406                   mu_,papermu,affmu,sigma,sigma*papermu);
     407#endif             
    342408          } else {
    343409            double phi;
    … …  
    364430          //std::cout<<"gap part "<<
    365431          //complementarityGap_*(beta2-tau)<<" after adding product = "<<xx<<std::endl;
    366           //xx=-1.0;
     432          xx=-1.0;
    367433          if (xx>0.0) {
    368434            double saveMu = mu_;
    … …  
    391457          }
    392458          //printf("product %g mu %g\n",product,mu_);
    393           if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) {
     459          if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0||
     460              tryJustPredictor) {
    394461            doCorrector=false;
    395462            bestNextGap=COIN_DBL_MAX;
    … …  
    405472              <<CoinMessageEol;
    406473            phase=2;
     474            tryJustPredictor=false;
    407475          } else {
    408476            phase=1;
    … …  
    464532      } /* endwhile */
    465533    } /* endwhile */
     534    if (numberIterations_==1) {
     535      FILE *fp = fopen("xx.xx","rb");
     536      if (fp) {
     537        int nrow=numberRows_;
     538        int ncol=numberColumns_;
     539        fread(dual_,sizeof(double),nrow,fp);
     540        fread(errorRegion_,sizeof(double),nrow,fp);
     541        fread(rhsFixRegion_,sizeof(double),nrow,fp);
     542        fread(solution_+ncol,sizeof(double),nrow,fp);
     543        fread(diagonal_+ncol,sizeof(double),nrow,fp);
     544        fread(wVec_+ncol,sizeof(double),nrow,fp);
     545        fread(zVec_+ncol,sizeof(double),nrow,fp);
     546        fread(upperSlack_+ncol,sizeof(double),nrow,fp);
     547        fread(lowerSlack_+ncol,sizeof(double),nrow,fp);
     548        fread(solution_,sizeof(double),ncol,fp);
     549        fread(diagonal_,sizeof(double),ncol,fp);
     550        fread(wVec_,sizeof(double),ncol,fp);
     551        fread(zVec_,sizeof(double),ncol,fp);
     552        fread(upperSlack_,sizeof(double),ncol,fp);
     553        fread(lowerSlack_,sizeof(double),ncol,fp);
     554        fclose(fp);
     555        // and flip
     556        int i;
     557        for (i=0;i<nrow;i++) {
     558          solution_[ncol+i]= -solution_[ncol+i];
     559          double value;
     560          value=lowerSlack_[ncol+i];
     561          lowerSlack_[ncol+i]=upperSlack_[ncol+i];
     562          upperSlack_[ncol+i]=value;
     563          value=zVec_[ncol+i];
     564          zVec_[ncol+i]=wVec_[ncol+i];
     565          wVec_[ncol+i]=value;
     566        }
     567        actualDualStep_=0.0;
     568        actualPrimalStep_=0.0;
     569      } else {
     570        printf("no file\n");
     571      }
     572    }
    466573    numberFixed=updateSolution();
    467574    numberFixedTotal+=numberFixed;
    … …  
    695802    }
    696803  } else {
    697     //iColumnust primal dual
     804    //just primal dual
    698805    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
    699806      if (!flagged(iColumn)) {
    … …  
    9011008#ifndef KKT
    9021009    double maximumRHS = maximumAbsElement(updateRegion_,numberRows_);
     1010    double saveMaximum = maximumRHS;
    9031011    double scale=1.0;
    9041012    double unscale=1.0;
    … …  
    9851093    }
    9861094    relativeError = maximumRHSError/solutionNorm_;
     1095    relativeError = maximumRHSError/saveMaximum;
    9871096    if (relativeError>tryError)
    9881097      relativeError=tryError;
    … …  
    11571266  }
    11581267  //   modify fixed RHS
    1159   multiplyAdd(dj_+numberColumns_,numberRows_,-1.0,rhsFixRegion_,0.0);
     1268  multiplyAdd(dj_+numberColumns_,numberRows_,1.0,rhsFixRegion_,0.0);
    11601269  //   create plausible RHS?
    11611270  matrix_->times(-1.0,dj_,rhsFixRegion_);
    1162   multiplyAdd(solution_+numberColumns_,numberRows_,1.0,errorRegion_,0.0);
    1163   matrix_->times(-1.0,solution_,errorRegion_);
     1271  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
     1272  matrix_->times(1.0,solution_,errorRegion_);
    11641273  rhsNorm_=maximumAbsElement(errorRegion_,numberRows_);
    11651274  if (rhsNorm_<1.0) {
    … …  
    12291338    <<initialValue<<objectiveNorm_
    12301339    <<CoinMessageEol;
    1231   int strategy=1;
     1340  int strategy=0;
    12321341  double extra=1.0e-10;
    12331342  double largeGap=1.0e15;
    12341343  double safeObjectiveValue=2.0*objectiveNorm_;
     1344  //safeObjectiveValue=1.0+objectiveNorm_;
     1345  //printf("temp safe obj value of %g\n",safeObjectiveValue);
    12351346  double safeFree=1.0e-1*initialValue;
    12361347  safeFree=1.0;
    … …  
    12451356      double low;
    12461357      double high;
    1247       double randomZ =0.5*CoinDrand48()+0.5;
    1248       double randomW =0.5*CoinDrand48()+0.5;
     1358      double randomZ =0.1*CoinDrand48()+0.9;
     1359      double randomW =0.1*CoinDrand48()+0.9;
    12491360      double setPrimal=initialValue;
    12501361      if (strategy==0) {
    … …  
    22532364    <<largestDiagonal<<smallestDiagonal
    22542365    <<CoinMessageEol;
     2366#if 0
     2367  // If diagonal wild - kill some
     2368  if (largestDiagonal>1.0e17*smallestDiagonal) {
     2369    double killValue =largestDiagonal*1.0e-17;
     2370    for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     2371      if (fabs(diagonal_[iColumn])<killValue)
     2372        diagonal_[iolumn]=0.0;
     2373    }
     2374  }
     2375#endif
    22552376  // update rhs region
    22562377  multiplyAdd(work1+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
    … …  
    22742395  solutionNorm_ = norm;
    22752396  //compute error and fixed RHS
    2276   multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,1.0);
     2397  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
    22772398  matrix_->times(1.0,solution_,errorRegion_);
    22782399  maximumDualError_=maximumDualError;

  • trunk/ClpPrimalColumnSteepest.cpp

    r328 r336  
    1212#include "CoinHelperFunctions.hpp"
    1313#include <stdio.h>
    14 
    1514
    1615//#############################################################################
    … …  
    198197  */
    199198  // Look at gub
     199#if 1
    200200  model_->clpMatrix()->dualExpanded(model_,updates,NULL,4);
     201#else
     202  updates->clear();
     203  model_->computeDuals(NULL);
     204#endif
    201205  if (updates->getNumElements()>1) {
    202206    // would have to have two goes for devex, three for steepest
    … …  
    514518        iSequence = index[i];
    515519        double value = infeas[iSequence];
     520        double weight = weights_[iSequence];
    516521        if (value>tolerance) {
    517           double weight = weights_[iSequence];
    518522          //weight=1.0;
    519523          if (value>bestDj*weight) {
    … …  
    527531            }
    528532          }
    529         }
    530         numberWanted--;
     533          numberWanted--;
     534        }
    531535        if (!numberWanted)
    532536          break;
    … …  
    548552            }
    549553          }
    550         }
    551         numberWanted--;
     554          numberWanted--;
     555        }
    552556        if (!numberWanted)
    553557          break;
    … …  
    557561      break;
    558562  }
     563  model_->clpMatrix()->setSavedBestSequence(bestSequence);
     564  if (bestSequence>=0)
     565    model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
    559566  if (sequenceOut>=0) {
    560567    infeas[sequenceOut]=saveOutInfeasibility;
    … …  
    985992    double pivotSquared;
    986993    int iSequence = index[j];
    987     double value = reducedCost[iSequence];
    988994    double value2 = updateBy[j];
    989     value -= value2;
    990     reducedCost[iSequence] = value;
     995    updateBy[j]=0.0;
    991996    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    992     modification = other[iSequence];
    993     updateBy[j]=0.0;
     997    double value;
    994998   
    995999    switch(status) {
    … …  
    9971001    case ClpSimplex::basic:
    9981002      infeasible_->zero(iSequence+addSequence);
     1003      reducedCost[iSequence] = 0.0;
    9991004    case ClpSimplex::isFixed:
    10001005      break;
    10011006    case ClpSimplex::isFree:
    10021007    case ClpSimplex::superBasic:
     1008      value = reducedCost[iSequence] - value2;
     1009      modification = other[iSequence];
    10031010      thisWeight = weight[iSequence];
    10041011      // row has -1
    … …  
    10071014     
    10081015      thisWeight += pivotSquared * devex_ + pivot * modification;
     1016      reducedCost[iSequence] = value;
    10091017      if (thisWeight<TRY_NORM) {
    10101018        if (mode_==1) {
    … …  
    10331041      break;
    10341042    case ClpSimplex::atUpperBound:
     1043      value = reducedCost[iSequence] - value2;
     1044      modification = other[iSequence];
    10351045      thisWeight = weight[iSequence];
    10361046      // row has -1
    … …  
    10391049     
    10401050      thisWeight += pivotSquared * devex_ + pivot * modification;
     1051      reducedCost[iSequence] = value;
    10411052      if (thisWeight<TRY_NORM) {
    10421053        if (mode_==1) {
    … …  
    10631074      break;
    10641075    case ClpSimplex::atLowerBound:
     1076      value = reducedCost[iSequence] - value2;
     1077      modification = other[iSequence];
    10651078      thisWeight = weight[iSequence];
    10661079      // row has -1
    … …  
    10691082     
    10701083      thisWeight += pivotSquared * devex_ + pivot * modification;
     1084      reducedCost[iSequence] = value;
    10711085      if (thisWeight<TRY_NORM) {
    10721086        if (mode_==1) {
    … …  
    27232737  /*if (model_->numberIterations()%100==0)
    27242738    printf("%d best %g\n",bestSequence,bestDj);*/
     2739  reducedCost=model_->djRegion();
     2740  model_->clpMatrix()->setSavedBestSequence(bestSequence);
     2741  if (bestSequence>=0)
     2742    model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
    27252743 
    27262744#ifdef CLP_DEBUG
    … …  
    35213539  const double * cost = model_->costRegion(1);
    35223540
     3541  model_->clpMatrix()->setOriginalWanted(numberWanted);
     3542  model_->clpMatrix()->setCurrentWanted(numberWanted);
    35233543  int iPassR=0,iPassC=0;
    35243544  // Setup two passes
    … …  
    35343554  startR[0]=(int) dstart;
    35353555  startR[3]=startR[0];
    3536   int startC[4];
    3537   startC[1]=numberColumns;
     3556  double startC[4];
     3557  startC[1]=1.0;
    35383558  startC[2]=0;
    35393559  double randomC = CoinDrand48();
    3540   dstart = ((double) numberColumns) * randomC;
    3541   startC[0]=(int) dstart;
    3542   startC[3]=startC[0];
     3560  startC[0]=randomC;
     3561  startC[3]=randomC;
    35433562  reducedCost = model_->djRegion(1);
    35443563  int sequenceOut = model_->sequenceOut();
    … …  
    35493568    int i;
    35503569    for (i=0;i<4;i++)
    3551       printf("start R %d C %d ",startR[i],startC[i]);
     3570      printf("start R %d C %g ",startR[i],startC[i]);
    35523571    printf("\n");
    35533572  }
    … …  
    35553574  bool finishedR=false,finishedC=false;
    35563575  bool doingR = randomR>randomC;
    3557   doingR=false;
     3576  //doingR=false;
    35583577  int saveNumberWanted = numberWanted;
    35593578  while (!finishedR||!finishedC) {
    … …  
    36393658        reducedCost[bestSequence] = cost[bestSequence] +duals[bestSequence-numberColumns];
    36403659        bestDj=fabs(reducedCost[bestSequence]);
    3641       }
     3660        model_->clpMatrix()->setSavedBestSequence(bestSequence);
     3661        model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
     3662      }
     3663      model_->clpMatrix()->setCurrentWanted(numberWanted);
    36423664      if (!numberWanted)
    36433665        break;
    … …  
    36573679      int saveSequence = bestSequence;
    36583680      // Columns
    3659       int start = startC[iPassC];
    3660       int end = min(startC[iPassC+1],start+chunk);;
     3681      double start = startC[iPassC];
     3682      // If we put this idea back then each function needs to update endFraction **
     3683#if 0
     3684      double dchunk = ((double) chunk)/((double) numberColumns);
     3685      double end = min(startC[iPassC+1],start+dchunk);;
     3686#else
     3687      double end=startC[iPassC+1]; // force end
     3688#endif
    36613689      model_->clpMatrix()->partialPricing(model_,start,end,bestSequence,numberWanted);
    3662       numberLook -= (end-start);
     3690      numberWanted=model_->clpMatrix()->currentWanted();
     3691      numberLook -= (int) ((end-start)*numberColumns);
    36633692      if (numberLook<0&&(10*(saveNumberWanted-numberWanted)>saveNumberWanted))
    36643693        numberWanted=0; // give up
    … …  
    36713700      doingR=true;
    36723701      // update start
    3673       int iSequence = end;
    3674       startC[iPassC]=iSequence;
    3675       if (iSequence>=startC[iPassC+1]) {
     3702      startC[iPassC]=end;
     3703      if (end>=startC[iPassC+1]-1.0e-8) {
    36763704        if (iPassC)
    36773705          finishedC=true;

  • trunk/ClpSimplex.cpp

    r335 r336  
    4343  remainingDualInfeasibility_(0.0),
    4444  largeValue_(1.0e15),
     45  objectiveScale_(1.0),
     46  rhsScale_(1.0),
    4547  largestPrimalError_(0.0),
    4648  largestDualError_(0.0),
    … …  
    110112  incomingInfeasibility_(1.0),
    111113  allowedInfeasibility_(10.0),
     114  automaticScale_(0),
    112115  progress_(NULL)
    113116{
    … …  
    148151  remainingDualInfeasibility_(0.0),
    149152  largeValue_(1.0e15),
     153  objectiveScale_(1.0),
     154  rhsScale_(1.0),
    150155  largestPrimalError_(0.0),
    151156  largestDualError_(0.0),
    … …  
    215220  incomingInfeasibility_(1.0),
    216221  allowedInfeasibility_(10.0),
     222  automaticScale_(0),
    217223  progress_(NULL)
    218224{
    … …  
    348354  // now check solutions
    349355  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    350   objectiveValue_ += objectiveModification;
     356  objectiveValue_ += objectiveModification/(objectiveScale_*rhsScale_);
    351357  checkDualSolution();
    352358  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
    … …  
    423429  }
    424430  arrayVector.setNumElements(number);
    425 
     431#ifdef CLP_DEBUG
     432  if (numberIterations_==-3840) {
     433    int i;
     434    for (i=0;i<numberRows_+numberColumns_;i++)
     435      printf("%d status %d\n",i,status_[i]);
     436    printf("xxxxx1\n");
     437    for (i=0;i<numberRows_;i++)
     438      if (array[i])
     439        printf("%d rhs %g\n",i,array[i]);
     440    printf("xxxxx2\n");
     441    for (i=0;i<numberRows_+numberColumns_;i++)
     442      if (getStatus(i)!=basic)
     443        printf("%d non basic %g %g %g\n",i,lower_[i],solution_[i],upper_[i]);
     444    printf("xxxxx3\n");
     445  }
     446#endif
    426447  // Ftran adjusted RHS and iterate to improve accuracy
    427448  double lastError=COIN_DBL_MAX;
    … …  
    431452  CoinIndexedVector * lastVector = &previousVector;
    432453  factorization_->updateColumn(workSpace,thisVector);
     454#ifdef CLP_DEBUG
     455  if (numberIterations_==-3840) {
     456    int i;
     457    for (i=0;i<numberRows_;i++)
     458      if (array[i])
     459        printf("%d after rhs %g\n",i,array[i]);
     460    printf("xxxxx4\n");
     461  }
     462#endif
    433463  bool goodSolution=true;
    434464  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    … …  
    525555    }
    526556  }
     557#ifdef CLP_DEBUG
     558  if (numberIterations_==-3840) {
     559    exit(77);
     560  }
     561#endif
    527562}
    528563// now dual side
    … …  
    12471282  // Update hidden stuff e.g. effective RHS and gub
    12481283  matrix_->updatePivot(this,oldIn,oldOut);
    1249   objectiveValue_ += objectiveChange;
     1284  objectiveValue_ += objectiveChange/(objectiveScale_*rhsScale_);
    12501285  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
    12511286    <<numberIterations_<<objectiveValue()
    … …  
    12631298  //exit(77);
    12641299  // check for small cycles
    1265   int cycle=progress_->cycle(sequenceIn_,sequenceOut_,
     1300  int in = sequenceIn_;
     1301  int out = sequenceOut_;
     1302  matrix_->correctSequence(in,out);
     1303  int cycle=progress_->cycle(in,out,
    12661304                            directionIn_,directionOut_);
    12671305  if (cycle>0) {
    … …  
    13231361  remainingDualInfeasibility_(0.0),
    13241362  largeValue_(1.0e15),
     1363  objectiveScale_(1.0),
     1364  rhsScale_(1.0),
    13251365  largestPrimalError_(0.0),
    13261366  largestDualError_(0.0),
    … …  
    13901430  incomingInfeasibility_(1.0),
    13911431  allowedInfeasibility_(10.0),
     1432  automaticScale_(0),
    13921433  progress_(NULL)
    13931434{
    … …  
    14221463  remainingDualInfeasibility_(0.0),
    14231464  largeValue_(1.0e15),
     1465  objectiveScale_(1.0),
     1466  rhsScale_(1.0),
    14241467  largestPrimalError_(0.0),
    14251468  largestDualError_(0.0),
    … …  
    14891532  incomingInfeasibility_(1.0),
    14901533  allowedInfeasibility_(10.0),
     1534  automaticScale_(0),
    14911535  progress_(NULL)
    14921536{
    … …  
    15851629  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
    15861630  largeValue_ = rhs.largeValue_;
     1631  objectiveScale_ = rhs.objectiveScale_;
     1632  rhsScale_ = rhs.rhsScale_;
    15871633  largestPrimalError_ = rhs.largestPrimalError_;
    15881634  largestDualError_ = rhs.largestDualError_;
    … …  
    16301676  incomingInfeasibility_ = rhs.incomingInfeasibility_;
    16311677  allowedInfeasibility_ = rhs.allowedInfeasibility_;
     1678  automaticScale_ = rhs.automaticScale_;
    16321679  if (rhs.progress_)
    16331680    progress_ = new ClpSimplexProgress(*rhs.progress_);
    … …  
    18141861    }
    18151862  }
     1863  objectiveValue_ /= (objectiveScale_*rhsScale_);
    18161864}
    18171865void
    … …  
    22312279    if (matrix_->scale(this))
    22322280      scalingFlag_=-scalingFlag_; // not scaled after all
     2281    if (rowScale_&&automaticScale_) {
     2282      // try automatic scaling
     2283      double smallestObj=1.0e100;
     2284      double largestObj=0.0;
     2285      double largestRhs=0.0;
     2286      for (i=0;i<numberColumns_+numberRows_;i++) {
     2287        double value = fabs(cost_[i]);
     2288        if (i<numberColumns_)
     2289          value *= columnScale_[i];
     2290        if (value&&lower_[i]!=upper_[i]) {
     2291          smallestObj = min(smallestObj,value);
     2292          largestObj = max(largestObj,value);
     2293        }
     2294        if (lower_[i]>0.0||upper_[i]<0.0) {
     2295          double scale;
     2296          if (i<numberColumns_)
     2297            scale = 1.0/columnScale_[i];
     2298          else
     2299            scale = rowScale_[i-numberColumns_];
     2300          //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
     2301          if (lower_[i]>0)
     2302            largestRhs=max(largestRhs,lower_[i]*scale);
     2303          if (upper_[i]<0.0)
     2304            largestRhs=max(largestRhs,-upper_[i]*scale);
     2305        }
     2306      }
     2307      printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
     2308      bool scalingDone=false;
     2309      // look at element range
     2310      double smallestNegative;
     2311      double largestNegative;
     2312      double smallestPositive;
     2313      double largestPositive;
     2314      matrix_->rangeOfElements(smallestNegative, largestNegative,
     2315                               smallestPositive, largestPositive);
     2316      smallestPositive = min(fabs(smallestNegative),smallestPositive);
     2317      largestPositive = max(fabs(largestNegative),largestPositive);
     2318      if (largestObj) {
     2319        double ratio = largestObj/smallestObj;
     2320        double scale=1.0;
     2321        if (ratio<1.0e8) {
     2322          // reasonable
     2323          if (smallestObj<1.0e-4) {
     2324            // may as well scale up
     2325            scalingDone=true;
     2326            scale=1.0e-3/smallestObj;
     2327          } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
     2328            //done=true;
     2329          } else {
     2330            scalingDone=true;
     2331            if (algorithm_<0) {
     2332              scale = 1.0e6/largestObj;
     2333            } else {
     2334              scale = max(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
     2335            }
     2336          }
     2337        } else if (ratio<1.0e12) {
     2338          // not so good
     2339          if (smallestObj<1.0e-7) {
     2340            // may as well scale up
     2341            scalingDone=true;
     2342            scale=1.0e-6/smallestObj;
     2343          } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
     2344            //done=true;
     2345          } else {
     2346            scalingDone=true;
     2347            if (algorithm_<0) {
     2348              scale = 1.0e7/largestObj;
     2349            } else {
     2350              scale = max(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
     2351            }
     2352          }
     2353        } else {
     2354          // Really nasty problem
     2355          if (smallestObj<1.0e-8) {
     2356            // may as well scale up
     2357            scalingDone=true;
     2358            scale=1.0e-7/smallestObj;
     2359            largestObj *= scale;
     2360          }
     2361          if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
     2362            //done=true;
     2363          } else {
     2364            scalingDone=true;
     2365            if (algorithm_<0) {
     2366              scale = 1.0e7/largestObj;
     2367            } else {
     2368              scale = max(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
     2369            }
     2370          }
     2371        }
     2372        objectiveScale_=scale;
     2373      }
     2374      if (largestRhs>1.0e12) {
     2375        scalingDone=true;
     2376        rhsScale_=1.0e9/largestRhs;
     2377      } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
     2378        scalingDone=true;
     2379        rhsScale_=1.0e6/largestRhs;
     2380      } else {
     2381        rhsScale_=1.0;
     2382      }
     2383      if (scalingDone) {
     2384        handler_->message(CLP_RIM_SCALE,messages_)
     2385          <<objectiveScale_<<rhsScale_
     2386          <<CoinMessageEol;
     2387      }
     2388    }
    22332389  }
    22342390  if ((what&4)!=0) {
    2235     double direction = optimizationDirection_;
    2236     // direction is actually scale out not scale in
    2237     if (direction)
    2238       direction = 1.0/direction;
     2391    double direction = optimizationDirection_*objectiveScale_;
    22392392    // but also scale by scale factors
    22402393    // not really sure about row scaling
    … …  
    22532406    }
    22542407  }
    2255   if ((what&(1+32))!=0&&rowScale_) {
    2256     for (i=0;i<numberColumns_;i++) {
    2257       double multiplier = 1.0/columnScale_[i];
    2258       if (columnLowerWork_[i]>-1.0e50)
    2259         columnLowerWork_[i] *= multiplier;
    2260       if (columnUpperWork_[i]<1.0e50)
    2261         columnUpperWork_[i] *= multiplier;
    2262      
    2263     }
    2264     for (i=0;i<numberRows_;i++) {
    2265       if (rowLowerWork_[i]>-1.0e50)
    2266         rowLowerWork_[i] *= rowScale_[i];
    2267       if (rowUpperWork_[i]<1.0e50)
    2268         rowUpperWork_[i] *= rowScale_[i];
    2269     }
    2270   }
    2271   if ((what&(8+32))!=0&&rowScale_) {
    2272     // on entry
    2273     for (i=0;i<numberColumns_;i++) {
    2274       columnActivityWork_[i] /= columnScale_[i];
    2275       reducedCostWork_[i] *= columnScale_[i];
    2276     }
    2277     for (i=0;i<numberRows_;i++) {
    2278       rowActivityWork_[i] *= rowScale_[i];
    2279       dual_[i] /= rowScale_[i];
    2280       rowReducedCost_[i] = dual_[i];
     2408  if ((what&(1+32))!=0) {
     2409    if(rowScale_) {
     2410      for (i=0;i<numberColumns_;i++) {
     2411        double multiplier = rhsScale_/columnScale_[i];
     2412        if (columnLowerWork_[i]>-1.0e50)
     2413          columnLowerWork_[i] *= multiplier;
     2414        if (columnUpperWork_[i]<1.0e50)
     2415          columnUpperWork_[i] *= multiplier;
     2416       
     2417      }
     2418      for (i=0;i<numberRows_;i++) {
     2419        double multiplier = rhsScale_*rowScale_[i];
     2420        if (rowLowerWork_[i]>-1.0e50)
     2421          rowLowerWork_[i] *= multiplier;
     2422        if (rowUpperWork_[i]<1.0e50)
     2423          rowUpperWork_[i] *= multiplier;
     2424      }
     2425    } else if (rhsScale_!=1.0) {
     2426      for (i=0;i<numberColumns_+numberRows_;i++) {
     2427        if (lower_[i]>-1.0e50)
     2428          lower_[i] *= rhsScale_;
     2429        if (upper_[i]<1.0e50)
     2430          upper_[i] *= rhsScale_;
     2431       
     2432      }
     2433    }
     2434  }
     2435  if ((what&(8+32))!=0) {
     2436    if (rowScale_) {
     2437      // on entry
     2438      for (i=0;i<numberColumns_;i++) {
     2439        columnActivityWork_[i] /= columnScale_[i];
     2440        columnActivityWork_[i] *= rhsScale_;
     2441        reducedCostWork_[i] *= objectiveScale_*columnScale_[i];
     2442      }
     2443      for (i=0;i<numberRows_;i++) {
     2444        rowActivityWork_[i] *= rowScale_[i]*rhsScale_;
     2445        dual_[i] /= rowScale_[i];
     2446        dual_[i] *= objectiveScale_;
     2447        rowReducedCost_[i] = dual_[i];
     2448      }
     2449    } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
     2450      // on entry
     2451      for (i=0;i<numberColumns_;i++) {
     2452        columnActivityWork_[i] *= rhsScale_;
     2453        reducedCostWork_[i] *= objectiveScale_;
     2454      }
     2455      for (i=0;i<numberRows_;i++) {
     2456        rowActivityWork_[i] *= rhsScale_;
     2457        dual_[i] *= objectiveScale_;
     2458        rowReducedCost_[i] = dual_[i];
     2459      }
    22812460    }
    22822461  }
    … …  
    23662545    ray_=NULL;
    23672546  }
     2547  // set upperOut_ to furthest away from bound so can use in dual for dualBound_
     2548  upperOut_=0.0;
    23682549  // ray may be null if in branch and bound
    23692550  if (rowScale_) {
    … …  
    23732554    int numberDualScaled=0;
    23742555    int numberDualUnscaled=0;
     2556    double scaleC = 1.0/objectiveScale_;
     2557    double scaleR = 1.0/rhsScale_;
    23752558    for (i=0;i<numberColumns_;i++) {
    23762559      double scaleFactor = columnScale_[i];
    23772560      double valueScaled = columnActivityWork_[i];
    2378       if (valueScaled<columnLowerWork_[i]-primalTolerance_)
    2379         numberPrimalScaled++;
    2380       else if (valueScaled>columnUpperWork_[i]+primalTolerance_)
    2381         numberPrimalScaled++;
    2382       columnActivity_[i] = valueScaled*scaleFactor;
     2561      double lowerScaled = columnLowerWork_[i];
     2562      double upperScaled = columnUpperWork_[i];
     2563      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     2564        if (valueScaled<lowerScaled-primalTolerance_)
     2565          numberPrimalScaled++;
     2566        else
     2567          upperOut_ = max(upperOut_,valueScaled-lowerScaled);
     2568        if (valueScaled>upperScaled+primalTolerance_)
     2569          numberPrimalScaled++;
     2570        else
     2571          upperOut_ = max(upperOut_,upperScaled-valueScaled);
     2572      }
     2573      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
    23832574      double value = columnActivity_[i];
    23842575      if (value<columnLower_[i]-primalTolerance_)
    … …  
    23912582      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    23922583        numberDualScaled++;
    2393       reducedCost_[i] = valueScaledDual/scaleFactor;
     2584      reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
    23942585      double valueDual = reducedCost_[i];
    23952586      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    … …  
    24012592      double scaleFactor = rowScale_[i];
    24022593      double valueScaled = rowActivityWork_[i];
    2403       if (valueScaled<rowLowerWork_[i]-primalTolerance_)
    2404         numberPrimalScaled++;
    2405       else if (valueScaled>rowUpperWork_[i]+primalTolerance_)
    2406         numberPrimalScaled++;
    2407       rowActivity_[i] = valueScaled/scaleFactor;
     2594      double lowerScaled = rowLowerWork_[i];
     2595      double upperScaled = rowUpperWork_[i];
     2596      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     2597        if (valueScaled<lowerScaled-primalTolerance_)
     2598          numberPrimalScaled++;
     2599        else
     2600          upperOut_ = max(upperOut_,valueScaled-lowerScaled);
     2601        if (valueScaled>upperScaled+primalTolerance_)
     2602          numberPrimalScaled++;
     2603        else
     2604          upperOut_ = max(upperOut_,upperScaled-valueScaled);
     2605      }
     2606      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
    24082607      double value = rowActivity_[i];
    24092608      if (value<rowLower_[i]-primalTolerance_)
    … …  
    24162615      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    24172616        numberDualScaled++;
    2418       dual_[i] = valueScaledDual*scaleFactor;
     2617      dual_[i] = valueScaledDual*scaleFactor*scaleC;
    24192618      double valueDual = dual_[i];
    24202619      if (rowObjective_)
    … …  
    24462645      }
    24472646    }
     2647  } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
     2648    // Collect infeasibilities
     2649    int numberPrimalScaled=0;
     2650    int numberPrimalUnscaled=0;
     2651    int numberDualScaled=0;
     2652    int numberDualUnscaled=0;
     2653    double scaleC = 1.0/objectiveScale_;
     2654    double scaleR = 1.0/rhsScale_;
     2655    for (i=0;i<numberColumns_;i++) {
     2656      double valueScaled = columnActivityWork_[i];
     2657      double lowerScaled = columnLowerWork_[i];
     2658      double upperScaled = columnUpperWork_[i];
     2659      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     2660        if (valueScaled<lowerScaled-primalTolerance_)
     2661          numberPrimalScaled++;
     2662        else
     2663          upperOut_ = max(upperOut_,valueScaled-lowerScaled);
     2664        if (valueScaled>upperScaled+primalTolerance_)
     2665          numberPrimalScaled++;
     2666        else
     2667          upperOut_ = max(upperOut_,upperScaled-valueScaled);
     2668      }
     2669      columnActivity_[i] = valueScaled*scaleR;
     2670      double value = columnActivity_[i];
     2671      if (value<columnLower_[i]-primalTolerance_)
     2672        numberPrimalUnscaled++;
     2673      else if (value>columnUpper_[i]+primalTolerance_)
     2674        numberPrimalUnscaled++;
     2675      double valueScaledDual = reducedCostWork_[i];
     2676      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     2677        numberDualScaled++;
     2678      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     2679        numberDualScaled++;
     2680      reducedCost_[i] = valueScaledDual*scaleC;
     2681      double valueDual = reducedCost_[i];
     2682      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     2683        numberDualUnscaled++;
     2684      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     2685        numberDualUnscaled++;
     2686    }
     2687    for (i=0;i<numberRows_;i++) {
     2688      double valueScaled = rowActivityWork_[i];
     2689      double lowerScaled = rowLowerWork_[i];
     2690      double upperScaled = rowUpperWork_[i];
     2691      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     2692        if (valueScaled<lowerScaled-primalTolerance_)
     2693          numberPrimalScaled++;
     2694        else
     2695          upperOut_ = max(upperOut_,valueScaled-lowerScaled);
     2696        if (valueScaled>upperScaled+primalTolerance_)
     2697          numberPrimalScaled++;
     2698        else
     2699          upperOut_ = max(upperOut_,upperScaled-valueScaled);
     2700      }
     2701      rowActivity_[i] = valueScaled*scaleR;
     2702      double value = rowActivity_[i];
     2703      if (value<rowLower_[i]-primalTolerance_)
     2704        numberPrimalUnscaled++;
     2705      else if (value>rowUpper_[i]+primalTolerance_)
     2706        numberPrimalUnscaled++;
     2707      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
     2708      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     2709        numberDualScaled++;
     2710      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     2711        numberDualScaled++;
     2712      dual_[i] = valueScaledDual*scaleC;
     2713      double valueDual = dual_[i];
     2714      if (rowObjective_)
     2715        valueDual += rowObjective_[i];
     2716      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     2717        numberDualUnscaled++;
     2718      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     2719        numberDualUnscaled++;
     2720    }
     2721    if (!problemStatus_&&!secondaryStatus_) {
     2722      // See if we need to set secondary status
     2723      if (numberPrimalUnscaled) {
     2724        if (numberDualUnscaled)
     2725          secondaryStatus_=4;
     2726        else
     2727          secondaryStatus_=2;
     2728      } else {
     2729        if (numberDualUnscaled)
     2730          secondaryStatus_=3;
     2731      }
     2732    }
    24482733  } else {
    24492734    if (columnActivityWork_) {
    24502735      for (i=0;i<numberColumns_;i++) {
     2736        double value = columnActivityWork_[i];
     2737        double lower = columnLowerWork_[i];
     2738        double upper = columnUpperWork_[i];
     2739        if (lower>-1.0e20||upper<1.0e20) {
     2740          if (value>lower)
     2741            upperOut_ = max(upperOut_,value-lower);
     2742          if (value<upper)
     2743            upperOut_ = max(upperOut_,upper-value);
     2744        }
    24512745        columnActivity_[i] = columnActivityWork_[i];
    24522746        reducedCost_[i] = reducedCostWork_[i];
    24532747      }
    24542748      for (i=0;i<numberRows_;i++) {
     2749        double value = rowActivityWork_[i];
     2750        double lower = rowLowerWork_[i];
     2751        double upper = rowUpperWork_[i];
     2752        if (lower>-1.0e20||upper<1.0e20) {
     2753          if (value>lower)
     2754            upperOut_ = max(upperOut_,value-lower);
     2755          if (value<upper)
     2756            upperOut_ = max(upperOut_,upper-value);
     2757        }
    24552758        rowActivity_[i] = rowActivityWork_[i];
    24562759      }
    24572760    }
     2761  }
     2762  // switch off scalefactor if auto
     2763  if (automaticScale_) {
     2764    rhsScale_=1.0;
     2765    objectiveScale_=1.0;
    24582766  }
    24592767  if (optimizationDirection_!=1.0) {
    … …  
    24682776  if(getRidOfFactorizationData>=0)
    24692777    gutsOfDelete(getRidOfFactorizationData+1);
     2778  // get rid of data
     2779  matrix_->generalExpanded(this,13,scalingFlag_);
    24702780}
    24712781void
    … …  
    31203430    // check which algorithms allowed
    31213431    int dummy;
    3122     if ((matrix_->generalExpanded(this,4,dummy)&2)!=0)
     3432    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
     3433      // upperOut_ has largest away from bound
     3434      double saveBound=dualBound_;
     3435      if (upperOut_>0.0)
     3436        dualBound_=2.0*upperOut_;
    31233437      returnCode = ((ClpSimplexPrimal *) this)->dual(0);
    3124     else
     3438      dualBound_=saveBound;
     3439    } else {
    31253440      returnCode = ((ClpSimplexDual *) this)->primal(0);
     3441    }
    31263442    setInitialDenseFactorization(denseFactorization);
    31273443    perturbation_=savePerturbation;
    … …  
    31913507  perturbation_ = otherModel.perturbation_;
    31923508  specialOptions_ = otherModel.specialOptions_;
     3509  automaticScale_ = otherModel.automaticScale_;
     3510  objectiveScale_ = otherModel.objectiveScale_;
     3511  rhsScale_ = otherModel.rhsScale_;
    31933512}
    31943513typedef struct {
    … …  
    44934812   
    44944813    // if stable replace in basis
    4495     int updateStatus = factorization_->replaceColumn(rowArray_[2],
    4496                                                    pivotRow_,
     4814    int updateStatus = factorization_->replaceColumn(this,
     4815                                                     rowArray_[2],
     4816                                                     rowArray_[1],
     4817                                                     pivotRow_,
    44974818                                                     alpha_);
    44984819    bool takePivot=true;

  • trunk/ClpSimplexDual.cpp

    r333 r336  
    573573#endif
    574574#if 0
    575     if (!factorization_->pivots()){
     575    //    if (factorization_->pivots()){
     576    {
    576577      int iPivot;
    577578      double * array = rowArray_[3]->denseVector();
    … …  
    596597        double upperValue=upper_[iSequence];
    597598        double value=solution_[iSequence];
    598         if(getStatus(iSequence)!=basic) {
     599        if(getStatus(iSequence)!=basic&&getStatus(iSequence)!=isFree) {
    599600          assert(lowerValue>-1.0e20);
    600601          assert(upperValue<1.0e20);
    … …  
    843844        assert(fabs(dualOut_)<1.0e50);
    844845        // if stable replace in basis
    845         int updateStatus = factorization_->replaceColumn(rowArray_[2],
     846        int updateStatus = factorization_->replaceColumn(this,
     847                                                         rowArray_[2],
     848                                                         rowArray_[1],
    846849                                                         pivotRow_,
    847850                                                         alpha_);
    … …  
    29582961    if (rowScale_) {
    29592962      if (rowLowerWork_[iRow]>-1.0e50)
    2960         rowLowerWork_[iRow] *= rowScale_[iRow];
     2963        rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
    29612964      if (rowUpperWork_[iRow]<1.0e50)
    2962         rowUpperWork_[iRow] *= rowScale_[iRow];
     2965        rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
     2966    } else if (rhsScale_!=1.0) {
     2967      if (rowLowerWork_[iRow]>-1.0e50)
     2968        rowLowerWork_[iRow] *= rhsScale_;
     2969      if (rowUpperWork_[iRow]<1.0e50)
     2970        rowUpperWork_[iRow] *= rhsScale_;
    29632971    }
    29642972  } else {
    … …  
    29692977      double multiplier = 1.0/columnScale_[iSequence];
    29702978      if (columnLowerWork_[iSequence]>-1.0e50)
    2971         columnLowerWork_[iSequence] *= multiplier;
     2979        columnLowerWork_[iSequence] *= multiplier*rhsScale_;
    29722980      if (columnUpperWork_[iSequence]<1.0e50)
    2973         columnUpperWork_[iSequence] *= multiplier;
     2981        columnUpperWork_[iSequence] *= multiplier*rhsScale_;
     2982    } else if (rhsScale_!=1.0) {
     2983      if (columnLowerWork_[iSequence]>-1.0e50)
     2984        columnLowerWork_[iSequence] *= rhsScale_;
     2985      if (columnUpperWork_[iSequence]<1.0e50)
     2986        columnUpperWork_[iSequence] *= rhsScale_;
    29742987    }
    29752988  }
    … …  
    30333046  int minLength=numberRows_;
    30343047  double averageCost = 0.0;
     3048  // look at element range
     3049  double smallestNegative;
     3050  double largestNegative;
     3051  double smallestPositive;
     3052  double largestPositive;
     3053  matrix_->rangeOfElements(smallestNegative, largestNegative,
     3054                           smallestPositive, largestPositive);
     3055  smallestPositive = min(fabs(smallestNegative),smallestPositive);
     3056  largestPositive = max(fabs(largestNegative),largestPositive);
     3057  double elementRatio = largestPositive/smallestPositive;
    30353058  int numberNonZero=0;
    30363059  if (!numberIterations_&&perturbation_==50) {
    … …  
    30683091    //exit(0);
    30693092#endif
     3093    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
     3094    //                                                                elementRatio);
    30703095    //number=0;
    3071     if (number*4>numberColumns_) {
     3096    if (number*4>numberColumns_||elementRatio>1.0e12) {
    30723097      perturbation_=100;
    30733098      return; // good enough
    … …  
    34213446    columnUpper_[iColumn] = newUpper[i];
    34223447    if (scalingFlag_<=0)
    3423       upper_[iColumn] = newUpper[i];
     3448      upper_[iColumn] = newUpper[i]*rhsScale_;
    34243449    else
    3425       upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
     3450      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
    34263451    // Start of fast iterations
    34273452    int status = fastDual(alwaysFinish);
    … …  
    34823507    columnLower_[iColumn] = newLower[i];
    34833508    if (scalingFlag_<=0)
    3484       lower_[iColumn] = newLower[i];
     3509      lower_[iColumn] = newLower[i]*rhsScale_;
    34853510    else
    3486       lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
     3511      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
    34873512    // Start of fast iterations
    34883513    status = fastDual(alwaysFinish);

  • trunk/ClpSimplexPrimal.cpp

    r335 r336  
    11// Copyright (C) 2002, International Business Machines
    22// Corporation and others.  All Rights Reserved.
    3 
    43
    54/* Notes on implementation of primal simplex algorithm.
    … …  
    745744  }
    746745  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
    747   if (dualFeasible()||problemStatus_==-4) {
    748     if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
     746  if ((dualFeasible()||problemStatus_==-4)&&!ifValuesPass) {
     747    if (nonLinearCost_->numberInfeasibilities()
     748        &&!alwaysOptimal) {
    749749      //may need infeasiblity cost changed
    750750      // we can see if we can construct a ray
    … …  
    10711071  // we also need somewhere for effective rhs
    10721072  double * rhs=rhsArray->denseVector();
    1073   //int * indexRhs = rhsArray->getIndices();
     1073  // and we can use indices to point to alpha
     1074  // that way we can store fabs(alpha)
     1075  int * indexPoint = rhsArray->getIndices();
    10741076  //int numberFlip=0; // Those which may change if flips
    10751077
    … …  
    11121114  bool gotList=false;
    11131115  int pivotOne=-1;
     1116  //#define CLP_DEBUG
     1117#ifdef CLP_DEBUG
     1118  if (numberIterations_==3839||numberIterations_==3840) {
     1119    double dj=cost_[sequenceIn_];
     1120    printf("cost in on %d is %g, dual in %g\n",sequenceIn_,dj,dualIn_);
     1121    for (iIndex=0;iIndex<number;iIndex++) {
     1122
     1123      int iRow = which[iIndex];
     1124      double alpha = work[iIndex];
     1125      int iPivot=pivotVariable_[iRow];
     1126      dj -= alpha*cost_[iPivot];
     1127      printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
     1128             iRow,iPivot,lower_[iPivot],solution_[iPivot],upper_[iPivot],
     1129             alpha, solution_[iPivot]-1.0e9*alpha,cost_[iPivot],dj);
     1130    }
     1131  }
     1132#endif
    11141133  while (!gotList) {
    11151134    pivotOne=-1;
    … …  
    11281147      double oldValue = solution_[iPivot];
    11291148      // get where in bound sequence
     1149      // note that after this alpha is actually fabs(alpha)
    11301150      if (alpha>0.0) {
    11311151        // basic variable going towards lower bound
    … …  
    11361156        double bound = upper_[iPivot];
    11371157        oldValue = bound-oldValue;
     1158        alpha = - alpha;
    11381159      }
    11391160     
    1140       double value = oldValue-tentativeTheta*fabs(alpha);
     1161      double value = oldValue-tentativeTheta*alpha;
    11411162      assert (oldValue>=-tolerance);
    11421163      if (value<=tolerance) {
    1143         value=oldValue-upperTheta*fabs(alpha);
    1144         if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot) {
    1145           upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
     1164        value=oldValue-upperTheta*alpha;
     1165        if (value<-primalTolerance_&&alpha>=acceptablePivot) {
     1166          upperTheta = (oldValue+primalTolerance_)/alpha;
    11461167          pivotOne=numberRemaining;
    11471168        }
    … …  
    11491170        spare[numberRemaining]=alpha;
    11501171        rhs[numberRemaining]=oldValue;
     1172        indexPoint[numberRemaining]=iIndex;
    11511173        index[numberRemaining++]=iRow;
    1152         totalThru += fabs(alpha);
     1174        totalThru += alpha;
    11531175        setActive(iRow);
    11541176        //} else if (value<primalTolerance_*1.002) {
    … …  
    11681190  }
    11691191  totalThru=0.0;
    1170   // we need to keep where rhs non-zeros are
    1171   int numberInRhs=numberRemaining;
    1172   memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
    1173   rhsArray->setNumElements(numberInRhs);
    1174   rhsArray->setPacked();
    11751192
    11761193  theta_=maximumMovement;
    … …  
    11841201    // First check if previously chosen one will work
    11851202    if (pivotOne>=0&&0) {
    1186       double thruCost = infeasibilityCost_*fabs(spare[pivotOne]);
     1203      double thruCost = infeasibilityCost_*spare[pivotOne];
    11871204      if (thruCost>=0.99*fabs(dualIn_))
    11881205        printf("Could pivot on %d as change %g dj %g\n",
    … …  
    11901207      double alpha = spare[pivotOne];
    11911208      double oldValue = rhs[pivotOne];
    1192       theta_ = oldValue/fabs(alpha);
     1209      theta_ = oldValue/alpha;
    11931210      pivotRow_=pivotOne;
    11941211      // Stop loop
    … …  
    12031220      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
    12041221       
    1205         double alpha = fabs(spare[iIndex]);
     1222        double alpha = spare[iIndex];
    12061223        double oldValue = rhs[iIndex];
    12071224        double value = oldValue-upperTheta*alpha;
    … …  
    12211238        double alpha = spare[iIndex];
    12221239        double oldValue = rhs[iIndex];
    1223         double value = oldValue-upperTheta*fabs(alpha);
     1240        double value = oldValue-upperTheta*alpha;
    12241241       
    12251242        if (value<=0||iBest==iIndex) {
    12261243          // how much would it cost to go thru and modify bound
    1227           totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha,rhs[iIndex]);
     1244          double trueAlpha=way*work[indexPoint[iIndex]];
     1245          totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],trueAlpha,rhs[iIndex]);
    12281246          setActive(iRow);
    1229           if (fabs(alpha)>bestPivot) {
    1230             bestPivot=fabs(alpha);
     1247          if (alpha>bestPivot) {
     1248            bestPivot=alpha;
    12311249            theta_ = oldValue/bestPivot;
    12321250            pivotRow_=iIndex;
    … …  
    12541272        if (bestPivot>bestEverPivot)
    12551273          bestEverPivot=bestPivot;
    1256       }
    1257     }
     1274      }    }
    12581275    // can get here without pivotRow_ set but with lastPivotRow
    12591276    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
    … …  
    12741291    int position=pivotRow_; // position in list
    12751292    pivotRow_=index[position];
    1276     //alpha_ = work[pivotRow_];
    1277     alpha_ = way*spare[position];
     1293    alpha_=work[indexPoint[position]];
    12781294    // translate to sequence
    12791295    sequenceOut_ = pivotVariable_[pivotRow_];
    … …  
    12911307    // will we need to increase tolerance
    12921308    //#define CLP_DEBUG
    1293 #ifdef CLP_DEBUG
    1294     bool found=false;
    1295 #endif
    12961309    double largestInfeasibility = primalTolerance_;
    12971310    if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
    … …  
    12991312      for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) {
    13001313        largestInfeasibility = max (largestInfeasibility,
    1301                                     -(rhs[iIndex]-fabs(spare[iIndex])*theta_));
     1314                                    -(rhs[iIndex]-spare[iIndex]*theta_));
    13021315      }
    13031316#define CLP_DEBUG
    … …  
    13651378
    13661379  // put back original bounds etc
     1380  memcpy(rhsArray->getIndices(),index,numberRemaining*sizeof(int));
     1381  rhsArray->setNumElements(numberRemaining);
     1382  rhsArray->setPacked();
    13671383  nonLinearCost_->goBackAll(rhsArray);
    1368 
    13691384  rhsArray->clear();
    13701385
    … …  
    16131628  if (!numberIterations_)
    16141629    cleanStatus(); // make sure status okay
     1630  // look at element range
     1631  double smallestNegative;
     1632  double largestNegative;
     1633  double smallestPositive;
     1634  double largestPositive;
     1635  matrix_->rangeOfElements(smallestNegative, largestNegative,
     1636                           smallestPositive, largestPositive);
     1637  smallestPositive = min(fabs(smallestNegative),smallestPositive);
     1638  largestPositive = max(fabs(largestNegative),largestPositive);
     1639  double elementRatio = largestPositive/smallestPositive;
    16151640  if (!numberIterations_&&perturbation_==50) {
    16161641    // See if we need to perturb
    … …  
    16401665    }
    16411666    delete [] sort;
    1642     if (number*4>numberRows_) {
     1667    //printf("ratio number diff rhs %g, element ratio %g\n",((double)number)/((double) numberRows_),
     1668    //                                                                elementRatio);
     1669    if (number*4>numberRows_||elementRatio>1.0e12) {
    16431670      perturbation_=100;
    16441671      return; // good enough
    … …  
    16951722    if (getRowStatus(iSequence)==basic)
    16961723      number++;
     1724  }
     1725  if (rhsScale_>100.0) {
     1726    // tone down perturbation
     1727    maximumFraction *= 0.1;
    16971728  }
    16981729  if (number!=numberRows_)
    … …  
    21662197      int updateStatus = matrix_->generalExpanded(this,3,updateType);
    21672198      if (updateType>=0)
    2168         updateStatus = factorization_->replaceColumn(rowArray_[2],
     2199        updateStatus = factorization_->replaceColumn(this,
     2200                                                     rowArray_[2],
     2201                                                     rowArray_[1],
    21692202                                                     pivotRow_,
    21702203                                                     alpha_);
    … …  
    23172350    nonLinearCost_->setOne(sequenceIn_,valueIn_);
    23182351    int whatNext=housekeeping(objectiveChange);
    2319 
    2320 #if 0
    2321     if (numberIterations_==1148)
    2322       whatNext=1;
    2323     if (numberIterations_>1200)
    2324     exit(0);
     2352#ifdef CLP_DEBUG
     2353    {
     2354      int ninf= matrix_->checkFeasible();
     2355      if (ninf)
     2356        printf("infeas %d\n",ninf);
     2357    }
    23252358#endif
    23262359    if (whatNext==1) {

  • trunk/ClpSimplexPrimalQuadratic.cpp

    r300 r336  
    10511051            if (cleanupIteration)
    10521052              factorization_->relaxAccuracyCheck(1.0e3*saveCheck);
    1053             updateStatus=factorization_->replaceColumn(rowArray_[2],
     1053            updateStatus=factorization_->replaceColumn(this,
     1054                                                       rowArray_[2],
     1055                                                       rowArray_[1],
    10541056                                                       pivotRow_,
    10551057                                                       alpha_);

  • trunk/Samples/decompose.cpp

    r225 r336  
    3232    // ken-18
    3333    for (iRow=104976;iRow<numberRows;iRow++)
     34      rowBlock[iRow]=-1;
     35  } else if (numberRows==2426) {
     36    // ken-7
     37    for (iRow=2401;iRow<numberRows;iRow++)
    3438      rowBlock[iRow]=-1;
    3539  } else if (numberRows==810) {
    … …  
    142146    sub[iBlock]=ClpSimplex(&model,numberRow2,whichRow,
    143147                           numberColumn2,whichColumn);
    144 #if 1
     148#if 0
    145149    // temp
    146150    double * upper = sub[iBlock].columnUpper();
    … …  
    219223    master.scaling(false);
    220224    if (0) {
    221       CoinMpsIO writer;
    222       writer.setMpsData(*master.matrix(), COIN_DBL_MAX,
    223                         master.getColLower(), master.getColUpper(),
    224                         master.getObjCoefficients(),
    225                         (const char*) 0 /*integrality*/,
    226                         master.getRowLower(), master.getRowUpper(),
    227                         NULL,NULL);
    228       writer.writeMps("yy.mps");
     225      master.writeMps("yy.mps");
    229226    }
    230227
    231228    master.primal();
     229    int problemStatus = master.status(); // do here as can change (delcols)
    232230    if (master.numberIterations()==0&&iPass)
    233231      break; // finished
    … …  
    263261    const double * dual=NULL;
    264262    bool deleteDual=false;
    265     if (master.isProvenOptimal()) {
     263    if (problemStatus==0) {
    266264      dual = master.dualRowSolution();
    267     } else if (master.isProvenPrimalInfeasible()) {
     265    } else if (problemStatus==1) {
    268266      // could do composite objective
    269267      dual = master.infeasibilityRay();
    … …  
    290288      top[iBlock].transposeTimes(dual,objective2);
    291289      int i;
    292       if (master.isProvenOptimal()) {
     290      if (problemStatus==0) {
    293291        for (i=0;i<numberColumns2;i++)
    294292          objective2[i] = saveObj[i]-objective2[i];
    … …  
    312310          int number=start;
    313311          double dj = objValue;
    314           if (!master.isProvenOptimal())
     312          if (problemStatus)
    315313            dj=0.0;
    316314          double smallest=1.0e100;
    … …  
    377375        }
    378376      }
     377      delete [] saveObj;
    379378    }
    380379    if (deleteDual)
    … …  
    395394  const double * fullLower = model.columnLower();
    396395  const double * fullUpper = model.columnUpper();
     396  const double * rowSolution = master.primalRowSolution();
     397  double * fullRowSolution = model.primalRowSolution();
     398  int kRow=0;
     399  for (iRow=0;iRow<numberRows;iRow++) {
     400    if (rowBlock[iRow]==-1) {
     401      model.setRowStatus(iRow,master.getRowStatus(kRow));
     402      fullRowSolution[iRow]=rowSolution[kRow++];
     403    }
     404  }
    397405  int kColumn=0;
    398   for (iRow=0;iRow<numberRows;iRow++)
    399     model.setRowStatus(iRow,ClpSimplex::superBasic);
    400406  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    401407    if (columnBlock[iColumn]==-1) {
    … …  
    413419      if (whichBlock[i-numberMasterColumns]==iBlock)
    414420        sol[i] = solution[i];
    415     memset(lower,0,numberMasterRows*sizeof(double));
     421    memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));
    416422    master.times(1.0,sol,lower);
    417423    for (iRow=0;iRow<numberMasterRows;iRow++) {
    … …  
    433439    sub[iBlock].primal();
    434440    const double * subSolution = sub[iBlock].primalColumnSolution();
     441    const double * subRowSolution = sub[iBlock].primalRowSolution();
    435442    // move solution
    436443    kColumn=0;
    … …  
    442449    }
    443450    assert(kColumn==sub[iBlock].numberColumns());
     451    kRow=0;
     452    for (iRow=0;iRow<numberRows;iRow++) {
     453      if (rowBlock[iRow]==iBlock) {
     454        model.setRowStatus(iRow,sub[iBlock].getRowStatus(kRow));
     455        fullRowSolution[iRow]=subRowSolution[kRow++];
     456      }
     457    }
     458    assert(kRow == sub[iBlock].numberRows()-numberMasterRows);
    444459  }
    445460  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    446461    if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
    447462        fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
    448       model.setStatus(iColumn,ClpSimplex::basic);
    449     }
    450   }
     463      assert(model.getStatus(iColumn)==ClpSimplex::basic);
     464    } else if (fullSolution[iColumn]>=fullUpper[iColumn]-1.0e-8) {
     465      // may help to make rest non basic
     466      model.setStatus(iColumn,ClpSimplex::atUpperBound);
     467    } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {
     468      // may help to make rest non basic
     469      model.setStatus(iColumn,ClpSimplex::atLowerBound);
     470    }
     471  }
     472  for (iRow=0;iRow<numberRows;iRow++)
     473    model.setRowStatus(iRow,ClpSimplex::superBasic);
    451474  model.primal(1);
    452475  delete [] sol;

  • trunk/Samples/driver.cpp

    r238 r336  
    2020    exit(1);
    2121  }
    22 
     22#ifdef STYLE1
    2323  if (argc<3 ||!strstr(argv[2],"primal")) {
    2424    // Use the dual algorithm unless user said "primal"
    … …  
    2727    model.initialPrimalSolve();
    2828  }
     29#else
     30  model.setLogLevel(8);
     31  ClpSolve solvectl;
    2932
     33
     34  if (argc<3 ||!strstr(argv[2],"primal")) {
     35    // Use the dual algorithm unless user said "primal"
     36    std::cout << std::endl << " Solve using Dual: " << std::endl;
     37    solvectl.setSolveType(ClpSolve::useDual);
     38    solvectl.setPresolveType(ClpSolve::presolveOn);
     39    model.initialSolve(solvectl);
     40  } else {
     41    std::cout << std::endl << " Solve using Primal: " << std::endl;
     42    solvectl.setSolveType(ClpSolve::usePrimal);
     43    solvectl.setPresolveType(ClpSolve::presolveOn);
     44    model.initialSolve(solvectl);
     45  }
     46#endif
    3047  std::string modelName;
    3148  model.getStrParam(ClpProbName,modelName);

  • trunk/Samples/dualCuts.cpp

    r256 r336  
    1 // Copyright (C) 2003, International Business Machines
    2 // Corporation and others.  All Rights Reserved.
     1/* Copyright (C) 2003, International Business Machines Corporation
     2   and others.  All Rights Reserved.
     3
     4   This sample program is designed to illustrate programming
     5   techniques using CoinLP, has not been thoroughly tested
     6   and comes without any warranty whatsoever.
     7
     8   You may copy, modify and distribute this sample program without
     9   any restrictions whatsoever and without any payment to anyone.
     10*/
    311
    412#include "ClpSimplex.hpp"

  • trunk/Samples/testGub.cpp

    r332 r336  
    3838    char * status = new char [numberColumns];
    3939    char * rowStatus = new char[numberRows];
    40     int i;
    4140    int n;
    4241    n = fread(solution,sizeof(double),numberColumns,fp);
    … …  
    4948    assert (n==numberRows);
    5049#if 0
     50    int i;
    5151    for (i=0;i<numberColumns;i++) {
    5252      if (status[i]==0)

  • trunk/Test/ClpMain.cpp

    r327 r336  
    5656 
    5757  DUALTOLERANCE=1,PRIMALTOLERANCE,DUALBOUND,PRIMALWEIGHT,MAXTIME,OBJSCALE,
     58  RHSSCALE,
    5859
    5960  LOGLEVEL=101,MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    6061 
    6162  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    62   PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,
     63  PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
    6364 
    6465  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
    … …  
    567568      break;
    568569    case OBJSCALE:
    569       model->setOptimizationDirection(value);
     570      model->setObjectiveScale(value);
     571      break;
     572    case RHSSCALE:
     573      model->setRhsScale(value);
    570574      break;
    571575    default:
    … …  
    596600    break;
    597601  case OBJSCALE:
    598     value=model->optimizationDirection();
     602    value=model->objectiveScale();
     603    break;
     604  case RHSSCALE:
     605    value=model->rhsScale();
    599606    break;
    600607  default:
    … …  
    891898       );
    892899    parameters[numberParameters++]=
     900      ClpItem("auto!Scale","Whether to scale objective, rhs and bounds of problem if they look odd",
     901              "off",AUTOSCALE,false);
     902    parameters[numberParameters-1].append("on");
     903    parameters[numberParameters-1].setLonghelp
     904      (
     905       "If you think you may get odd objective values or large equality rows etc then\
     906 it may be worth setting this true.  It is still experimental and you may prefer\
     907 to use objective!Scale and rhs!Scale"
     908       );
     909    parameters[numberParameters++]=
    893910      ClpItem("barr!ier","Solve using primal dual predictor corrector algorithm",
    894911              BARRIER);
    … …  
    11351152      ClpItem("objective!Scale","Scale factor to apply to objective",
    11361153              -1.0e20,1.0e20,OBJSCALE,false);
    1137     parameters[numberParameters-1].setDoubleValue(models->optimizationDirection());
     1154    parameters[numberParameters-1].setLonghelp
     1155      (
     1156       "If the objective function has some very large values, you may wish to scale them\
     1157 internally by this amount.  It can also be set by autoscale.  It is applied after scaling"
     1158       );
     1159    parameters[numberParameters-1].setDoubleValue(1.0);
    11381160    parameters[numberParameters++]=
    11391161      ClpItem("passP!resolve","How many passes in presolve",
    … …  
    12581280       "Useful for testing if maximization works correctly"
    12591281       );
     1282    parameters[numberParameters++]=
     1283      ClpItem("rhs!Scale","Scale factor to apply to rhs and bounds",
     1284              -1.0e20,1.0e20,RHSSCALE,false);
     1285    parameters[numberParameters-1].setLonghelp
     1286      (
     1287       "If the rhs or bounds have some very large meaningful values, you may wish to scale them\
     1288 internally by this amount.  It can also be set by autoscale"
     1289       );
     1290    parameters[numberParameters-1].setDoubleValue(1.0);
    12601291    parameters[numberParameters++]=
    12611292      ClpItem("save!Model","Save model to binary file",
    … …  
    15481579              models[iModel].scaling(action);
    15491580              break;
     1581            case AUTOSCALE:
     1582              models[iModel].setAutomaticScaling(action!=0);
     1583              break;
    15501584            case SPARSEFACTOR:
    15511585              models[iModel].setSparseFactorization((1-action)!=0);

  • trunk/include/ClpDynamicExampleMatrix.hpp

    r335 r336  
    3232public:
    3333  /// Partial pricing
    34   virtual void partialPricing(ClpSimplex * model, int start, int end,
     34  virtual void partialPricing(ClpSimplex * model, double start, double end,
    3535                      int & bestSequence, int & numberWanted);
    3636 

  • trunk/include/ClpDynamicMatrix.hpp

    r335 r336  
    2626  };
    2727  /// Partial pricing
    28   virtual void partialPricing(ClpSimplex * model, int start, int end,
     28  virtual void partialPricing(ClpSimplex * model, double start, double end,
    2929                      int & bestSequence, int & numberWanted);
    3030
    … …  
    7171      mode=10  - return 1 if there may be changing bounds on variable (column generation)
    7272      mode=11  - make sure set is clean (used when a variable rejected - but not flagged)
     73      mode=12  - after factorize but before permute stuff
     74      mode=13  - at end of simplex to delete stuff
    7375  */
    7476  virtual int generalExpanded(ClpSimplex * model,int mode,int & number);
    … …  
    251253  /// Saved best dual on gub row in pricing
    252254  double savedBestGubDual_;
    253   /// Saved best dj in pricing
    254   double savedBestDj_;
    255   /// Saved best sequence in pricing
    256   int savedBestSequence_;
    257255  /// Saved best set in pricing
    258256  int savedBestSet_;

  • trunk/include/ClpFactorization.hpp

    r225 r336  
    6464      after factorization and thereafter re-factorize
    6565   partial update already in U */
    66   int replaceColumn ( CoinIndexedVector * regionSparse,
     66  int replaceColumn ( const ClpSimplex * model,
     67                      CoinIndexedVector * regionSparse,
     68                      CoinIndexedVector * tableauColumn,
    6769                      int pivotRow,
    6870                      double pivotCheck ,

  • trunk/include/ClpGubDynamicMatrix.hpp

    r335 r336  
    1919public:
    2020  /// Partial pricing
    21   virtual void partialPricing(ClpSimplex * model, int start, int end,
     21  virtual void partialPricing(ClpSimplex * model, double start, double end,
    2222                              int & bestSequence, int & numberWanted);
    2323  /** This is local to Gub to allow synchronization:
    … …  
    5050  virtual void times(double scalar,
    5151                       const double * x, double * y) const;
     52  /** Just for debug
     53      Returns number of primal infeasibilities. Recomputes keys
     54  */
     55  virtual int checkFeasible() const;
    5256  //@}
    5357
    … …  
    156160  inline int numberElements() const
    157161  { return numberElements_;};
    158   /// Switches off dj checking each factorization (for BIG models)
    159   void switchOffCheck();
    160162  /// Status region for gub slacks
    161163  inline unsigned char * gubRowStatus() const
    … …  
    197199  /// Optional true upper bounds on sets
    198200  float * upperSet_;
    199   /// Pointer back to model
    200   ClpSimplex * model_;
    201201  /// size
    202202  int numberGubColumns_;

  • trunk/include/ClpGubMatrix.hpp

    r311 r336  
    4949                   int column, double multiplier) const;
    5050  /// Partial pricing
    51   virtual void partialPricing(ClpSimplex * model, int start, int end,
     51  virtual void partialPricing(ClpSimplex * model, double start, double end,
    5252                      int & bestSequence, int & numberWanted);
    5353  /// Returns number of hidden rows e.g. gub
    … …  
    119119      mode=10  - return 1 if there may be changing bounds on variable (column generation)
    120120      mode=11  - make sure set is clean (used when a variable rejected - but not flagged)
     121      mode=12  - after factorize but before permute stuff
     122      mode=13  - at end of simplex to delete stuff
    121123  */
    122124  virtual int generalExpanded(ClpSimplex * model,int mode,int & number);
    … …  
    145147  */
    146148  virtual int synchronize(ClpSimplex * model,int mode);
     149  /// Correct sequence in and out to give true value
     150  virtual void correctSequence(int & sequenceIn, int & sequenceOut) const;
    147151  //@}
    148152
    … …  
    263267  inline int numberSets() const
    264268  { return numberSets_;};
     269  /// Switches off dj checking each factorization (for BIG models)
     270  void ClpGubMatrix::switchOffCheck();
    265271   //@}
    266272   
    … …  
    278284  /// Sum of Primal infeasibilities using tolerance based on error in primals
    279285  double sumOfRelaxedPrimalInfeasibilities_;
     286  /// Infeasibility weight when last full pass done
     287  double infeasibilityWeight_;
    280288  /// Starts
    281289  int * start_;
    … …  
    308316  // Reverse pointer from index to set
    309317  int * fromIndex_;
     318  /// Pointer back to model
     319  ClpSimplex * model_;
    310320  /// Number of dual infeasibilities
    311321  int numberDualInfeasibilities_;
    312322  /// Number of primal infeasibilities
    313323  int numberPrimalInfeasibilities_;
     324  /** If pricing will declare victory (i.e. no check every factorization).
     325      -1 - always check
     326      0  - don't check
     327      1  - in don't check mode but looks optimal
     328  */
     329  int noCheck_;
    314330  /// Number of sets (gub rows)
    315331  int numberSets_;

  • trunk/include/ClpMatrixBase.hpp

    r328 r336  
    133133  virtual int hiddenRows() const;
    134134  /// Partial pricing
    135   virtual void partialPricing(ClpSimplex * model, int start, int end,
     135  virtual void partialPricing(ClpSimplex * model, double start, double end,
    136136                              int & bestSequence, int & numberWanted);
    137137  /** expands an updated column to allow for extra rows which the main
    … …  
    178178      mode=10  - return 1 if there may be changing bounds on variable (column generation)
    179179      mode=11  - make sure set is clean (used when a variable rejected - but not flagged)
     180      mode=12  - after factorize but before permute stuff
     181      mode=13  - at end of simplex to delete stuff
    180182  */
    181183  virtual int generalExpanded(ClpSimplex * model,int mode,int & number);
    … …  
    188190  */
    189191  virtual void createVariable(ClpSimplex * model, int & bestSequence);
     192  /** Just for debug if odd type matrix.
     193      Returns number of primal infeasibilities. */
     194  virtual int checkFeasible() const ;
    190195  /// Returns reduced cost of a variable
    191   virtual double reducedCost(ClpSimplex * model,int sequence) const;
     196  double reducedCost(ClpSimplex * model,int sequence) const;
     197  /// Correct sequence in and out to give true value
     198  virtual void correctSequence(int & sequenceIn, int & sequenceOut) const;
    192199  //@}
    193200 
    … …  
    276283  inline void setSkipDualCheck(bool yes)
    277284  { skipDualCheck_=yes;};
     285  /** Partial pricing tuning parameter - minimum number of "objects" to scan.
     286      e.g. number of Gub sets but could be number of variables */
     287  inline int minimumObjectsScan() const
     288  { return minimumObjectsScan_;};
     289  inline void setMinimumObjectsScan(int value)
     290  { minimumObjectsScan_=value;};
     291  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
     292  inline int minimumGoodReducedCosts() const
     293  { return minimumGoodReducedCosts_;};
     294  inline void setMinimumGoodReducedCosts(int value)
     295  { minimumGoodReducedCosts_=value;};
     296  /// Current start of search space in matrix (as fraction)
     297  inline double startFraction() const
     298  { return startFraction_;};
     299  inline void setStartFraction(double value)
     300  { startFraction_ = value;};
     301  /// Current end of search space in matrix (as fraction)
     302  inline double endFraction() const
     303  { return endFraction_;};
     304  inline void setEndFraction(double value)
     305  { endFraction_ = value;};
     306  /// Current best reduced cost
     307  inline double savedBestDj() const
     308  { return savedBestDj_;};
     309  inline void setSavedBestDj(double value)
     310  { savedBestDj_ = value;};
     311  /// Initial number of negative reduced costs wanted
     312  inline int originalWanted() const
     313  { return originalWanted_;};
     314  inline void setOriginalWanted(int value)
     315  { originalWanted_ = value;};
     316  /// Current number of negative reduced costs which we still need
     317  inline int currentWanted() const
     318  { return currentWanted_;};
     319  inline void setCurrentWanted(int value)
     320  { currentWanted_ = value;};
     321  /// Current best sequence
     322  inline int savedBestSequence() const
     323  { return savedBestSequence_;};
     324  inline void setSavedBestSequence(int value)
     325  { savedBestSequence_ = value;};
    278326  //@}
    279327 
    … …  
    306354      expensive */
    307355  double * rhsOffset_;
     356  /// Current start of search space in matrix (as fraction)
     357  double startFraction_;
     358  /// Current end of search space in matrix (as fraction)
     359  double endFraction_;
     360  /// Best reduced cost so far
     361  double savedBestDj_;
     362  /// Initial number of negative reduced costs wanted
     363  int originalWanted_;
     364  /// Current number of negative reduced costs which we still need
     365  int currentWanted_;
     366  /// Saved best sequence in pricing
     367  int savedBestSequence_;
    308368  /// type (may be useful)
    309369  int type_;
    … …  
    312372  /// If rhsOffset used this is refresh frequency (0==off)
    313373  int refreshFrequency_;
     374  /// Partial pricing tuning parameter - minimum number of "objects" to scan
     375  int minimumObjectsScan_;
     376  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
     377  int minimumGoodReducedCosts_;
     378  /// True sequence in (i.e. from larger problem)
     379  int trueSequenceIn_;
     380  /// True sequence out (i.e. from larger problem)
     381  int trueSequenceOut_;
    314382  /// whether to skip dual checks most of time
    315383  bool skipDualCheck_;

  • trunk/include/ClpMessage.hpp

    r263 r336  
    9393  CLP_BARRIER_FEASIBLE,
    9494  CLP_BARRIER_STEP,
     95  CLP_RIM_SCALE,
    9596  CLP_DUMMY_END
    9697};

  • trunk/include/ClpNetworkMatrix.hpp

    r286 r336  
    9595  virtual bool canDoPartialPricing() const;
    9696  /// Partial pricing
    97   virtual void partialPricing(ClpSimplex * model, int start, int end,
     97  virtual void partialPricing(ClpSimplex * model, double start, double end,
    9898                      int & bestSequence, int & numberWanted);
    9999   //@}

  • trunk/include/ClpPackedMatrix.hpp

    r299 r336  
    5454    /** Delete the columns whose indices are listed in <code>indDel</code>. */
    5555    virtual void deleteCols(const int numDel, const int * indDel)
    56   { matrix_->deleteCols(numDel,indDel);};
     56  { matrix_->deleteCols(numDel,indDel);numberActiveColumns_ = matrix_->getNumCols();};
    5757    /** Delete the rows whose indices are listed in <code>indDel</code>. */
    5858    virtual void deleteRows(const int numDel, const int * indDel)
    59   { matrix_->deleteRows(numDel,indDel);};
     59  { matrix_->deleteRows(numDel,indDel);numberActiveColumns_ = matrix_->getNumCols();};
    6060  /// Append Columns
    6161  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns)
    62   { matrix_->appendCols(number,columns);};
     62  { matrix_->appendCols(number,columns);numberActiveColumns_ = matrix_->getNumCols();};
    6363  /// Append Rows
    6464  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows)
    65   { matrix_->appendRows(number,rows);};
     65  { matrix_->appendRows(number,rows);numberActiveColumns_ = matrix_->getNumCols();};
    6666  /** Replace the elements of a vector.  The indices remain the same.
    6767      This is only needed if scaling and a row copy is used.
    … …  
    125125  virtual bool canDoPartialPricing() const;
    126126  /// Partial pricing
    127   virtual void partialPricing(ClpSimplex * model, int start, int end,
     127  virtual void partialPricing(ClpSimplex * model, double start, double end,
    128128                      int & bestSequence, int & numberWanted);
     129  /// makes sure active columns correct
     130  virtual int refresh(ClpSimplex * model);
    129131   //@}
    130132
    … …  
    160162                              CoinIndexedVector * z) const;
    161163    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
     164        Note - If x packed mode - then z packed mode
     165        This does by column and knows no gaps
     166        Squashes small elements and knows about ClpSimplex */
     167  void transposeTimesByColumn(const ClpSimplex * model, double scalar,
     168                              const CoinIndexedVector * x,
     169                              CoinIndexedVector * y,
     170                              CoinIndexedVector * z) const;
     171    /** Return <code>x * scalar * A in <code>z</code>.
     172        Note - this version when x packed mode and so will be z
     173        This does by column and knows no gaps and knows y empty
     174        Squashes small elements and knows about ClpSimplex */
     175  void transposeTimesByColumn(const ClpSimplex * model, double scalar,
     176                              const CoinIndexedVector * x,
     177                              CoinIndexedVector * z) const;
     178    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
    162179        Can use y as temporary array (will be empty at end)
    163180        Note - If x packed mode - then z packed mode
    … …  
    235252  /// Data
    236253  CoinPackedMatrix * matrix_;
     254  /// number of active columns (normally same as number of columns)
     255  int numberActiveColumns_;
    237256  /// Zero element flag - set true if any zero elements
    238257  bool zeroElements_;
     258  /// Gaps flag - set true if column start and length don't say contiguous
     259  bool hasGaps_;
    239260   //@}
    240261};

  • trunk/include/ClpPlusMinusOneMatrix.hpp

    r286 r336  
    191191  virtual bool canDoPartialPricing() const;
    192192  /// Partial pricing
    193   virtual void partialPricing(ClpSimplex * model, int start, int end,
     193  virtual void partialPricing(ClpSimplex * model, double start, double end,
    194194                      int & bestSequence, int & numberWanted);
    195195   //@}

  • trunk/include/ClpSimplex.hpp

    r335 r336  
    527527  inline int * pivotVariable() const
    528528          { return pivotVariable_;};
    529   /// Scaling of objective 12345.0 (auto), 0.0 (off), other user
     529  /// Scaling of objective
    530530  inline double objectiveScale() const
    531531          { return objectiveScale_;} ;
    532532  inline void setObjectiveScale(double value)
    533533          { objectiveScale_ = value;} ;
     534  /// Scaling of rhs and bounds
     535  inline double rhsScale() const
     536          { return rhsScale_;} ;
     537  inline void setRhsScale(double value)
     538          { rhsScale_ = value;} ;
     539  /// If automatic scaling on
     540  inline bool automaticScaling() const
     541  { return automaticScale_!=0;};
     542  inline void setAutomaticScaling(bool onOff)
     543  { automaticScale_ = onOff ? 1: 0;};
    534544  /// Current dual tolerance
    535545  inline double currentDualTolerance() const
    … …  
    828838  /// Scaling of objective
    829839  double objectiveScale_;
     840  /// Scaling of rhs and bounds
     841  double rhsScale_;
    830842  /// Largest error on Ax-b
    831843  double largestPrimalError_;
    … …  
    10001012  float incomingInfeasibility_;
    10011013  float allowedInfeasibility_;
     1014  /// Automatic scaling of objective and rhs and bounds
     1015  int automaticScale_;
    10021016  /// For dealing with all issues of cycling etc
    10031017  ClpSimplexProgress * progress_;
Note: See TracChangeset for help on using the changeset viewer.