Changeset 1785 for trunk


Ignore:
Timestamp:
Aug 25, 2011 6:17:49 AM (8 years ago)
Author:
forrest
Message:

patches to try and make parametrics faster

Location:
trunk/Clp/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpSimplex.cpp

    r1769 r1785  
    669669          }
    670670     }
     671#if CAN_HAVE_ZERO_OBJ>1
     672     if ((specialOptions_&2097152)==0) {
     673#endif
    671674     computeDuals(givenDuals);
    672675     if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
     
    692695     checkBothSolutions();
    693696     objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
     697#if CAN_HAVE_ZERO_OBJ>1
     698     } else {
     699       checkPrimalSolution( rowActivityWork_, columnActivityWork_);
     700#ifndef COIN_REUSE_RANDOM
     701       memset(dj_,0,(numberRows_+numberColumns_)*sizeof(double));
     702#else
     703       for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     704         double value;
     705         switch (getStatus(iSequence)) {
     706         case atLowerBound:
     707           value=1.0e-9*(1.0+CoinDrand48());
     708           break;
     709         case atUpperBound:
     710           value=-1.0e-9*(1.0+CoinDrand48());
     711           break;
     712         default:
     713           value=0.0;
     714           break;
     715         }
     716       }
     717#endif
     718       objectiveValue_=0.0;
     719     }
     720#endif
    694721     if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 ||
    695722                                      largestDualError_ > 1.0e-2))
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1780 r1785  
    12231223                    bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
    12241224                                                   columnArray_[1], acceptablePivot, dubiousWeights);
     1225#if CAN_HAVE_ZERO_OBJ>1
     1226                    if ((specialOptions_&2097152)!=0)
     1227                      theta_=0.0;
     1228#endif
    12251229               } else {
    12261230                    // Make sure direction plausible
     
    13991403                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
    14001404                    if (candidate == -1) {
     1405#if CLP_CAN_HAVE_ZERO_OBJ>1
     1406                    if ((specialOptions_&2097152)==0) {
     1407#endif
    14011408                         // make sure incoming doesn't count
    14021409                         Status saveStatus = getStatus(sequenceIn_);
     
    14061413                                                      objectiveChange, false);
    14071414                         setStatus(sequenceIn_, saveStatus);
     1415#if CLP_CAN_HAVE_ZERO_OBJ>1
     1416                    } else {
     1417                      rowArray_[0]->clear();
     1418                      rowArray_[2]->clear();
     1419                      columnArray_[0]->clear();
     1420                    }
     1421#endif
    14081422                    } else {
    14091423                         updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
     
    37783792                              if (value > dualTolerance_) {
    37793793                                   thisIncrease = true;
     3794#if CLP_CAN_HAVE_ZERO_OBJ<2
    37803795#define MODIFYCOST 2
     3796#endif
    37813797#if MODIFYCOST
    37823798                                   // modify cost to hit new tolerance
     
    46664682                         if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
    46674683                                   (numberChanged_ || (specialOptions_ & 4096) == 0)) {
     4684#if CLP_CAN_HAVE_ZERO_OBJ
    46684685                           if ((specialOptions_&2097152)==0) {
     4686#endif
    46694687                              doOriginalTolerance = 2;
    46704688                              numberTimesOptimal_++;
     
    46814699                              }
    46824700                              cleanDuals = 2; // If nothing changed optimal else primal
     4701#if CLP_CAN_HAVE_ZERO_OBJ
    46834702                           } else {
    46844703                             // no cost - skip checks
    46854704                             problemStatus_=0;
    46864705                           }
     4706#endif
    46874707                         } else {
    46884708                              problemStatus_ = 0; // optimal
     
    49394959                    if (givenDuals)
    49404960                         dualTolerance_ = 1.0e50;
     4961#if CLP_CAN_HAVE_ZERO_OBJ>1
     4962                    if ((specialOptions_&2097152)==0) {
     4963#endif
    49414964                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
    49424965                                      0.0, objectiveChange, true);
     4966#if CLP_CAN_HAVE_ZERO_OBJ>1
     4967                    } else {
     4968                      rowArray_[0]->clear();
     4969                      rowArray_[1]->clear();
     4970                      columnArray_[0]->clear();
     4971                    }
     4972#endif
    49434973                    dualTolerance_ = saveTolerance;
    49444974#if 0
     
    49715001                    if (givenDuals)
    49725002                         dualTolerance_ = 1.0e50;
     5003#if CLP_CAN_HAVE_ZERO_OBJ>1
     5004                    if ((specialOptions_&2097152)==0) {
     5005#endif
    49735006                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
    49745007                                      0.0, objectiveChange, true);
     5008#if CLP_CAN_HAVE_ZERO_OBJ>1
     5009                    } else {
     5010                      rowArray_[0]->clear();
     5011                      rowArray_[1]->clear();
     5012                      columnArray_[0]->clear();
     5013                    }
     5014#endif
    49755015                    dualTolerance_ = saveTolerance;
    49765016                    if (!numberIterations_ && sumPrimalInfeasibilities_ >
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1784 r1785  
    942942// Restores solution from dualized problem
    943943int
    944 ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
     944ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem,
     945                                 bool checkAccuracy)
    945946{
    946947     int returnCode = 0;;
     
    11871188     // Below will go to ..DEBUG later
    11881189#if 1 //ndef NDEBUG
    1189      // Check if correct
    1190      double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
    1191      double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
    1192      double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
    1193      double * dual = CoinCopyOfArray(dual_, numberRows_);
    1194      this->dual(); //primal();
    1195      CoinRelFltEq eq(1.0e-5);
    1196      for (iRow = 0; iRow < numberRows_; iRow++) {
    1197           assert(eq(dual[iRow], dual_[iRow]));
     1190     if (checkAccuracy) {
     1191       // Check if correct
     1192       double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
     1193       double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
     1194       double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
     1195       double * dual = CoinCopyOfArray(dual_, numberRows_);
     1196       this->dual(); //primal();
     1197       CoinRelFltEq eq(1.0e-5);
     1198       for (iRow = 0; iRow < numberRows_; iRow++) {
     1199         assert(eq(dual[iRow], dual_[iRow]));
     1200       }
     1201       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1202         assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
     1203       }
     1204       for (iRow = 0; iRow < numberRows_; iRow++) {
     1205         assert(eq(rowActivity[iRow], rowActivity_[iRow]));
     1206       }
     1207       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1208         assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
     1209       }
     1210       delete [] columnActivity;
     1211       delete [] rowActivity;
     1212       delete [] reducedCost;
     1213       delete [] dual;
    11981214     }
    1199      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1200           assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
    1201      }
    1202      for (iRow = 0; iRow < numberRows_; iRow++) {
    1203           assert(eq(rowActivity[iRow], rowActivity_[iRow]));
    1204      }
    1205      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1206           assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
    1207      }
    1208      delete [] columnActivity;
    1209      delete [] rowActivity;
    1210      delete [] reducedCost;
    1211      delete [] dual;
    12121215#endif
    12131216     return returnCode;
     
    19921995               while (!returnCode) {
    19931996                    //assert (reportIncrement);
    1994                     returnCode = parametricsLoop(startingTheta, endingTheta, reportIncrement,
     1997                 parametricsData paramData;
     1998                 paramData.startingTheta=startingTheta;
     1999                 paramData.endingTheta=endingTheta;
     2000                 paramData.lowerChange = chgLower;
     2001                 paramData.upperChange = chgUpper;
     2002                    returnCode = parametricsLoop(paramData, reportIncrement,
    19952003                                                 chgLower, chgUpper, chgObjective, data,
    19962004                                                 canTryQuick);
     2005                 startingTheta=paramData.startingTheta;
     2006                 endingTheta=paramData.endingTheta;
    19972007                    if (!returnCode) {
    19982008                         //double change = endingTheta-startingTheta;
     
    26562666}
    26572667int
    2658 ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
     2668ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
    26592669                                 const double * lowerChange, const double * upperChange,
    26602670                                 const double * changeObjective, ClpDataSave & data,
    26612671                                 bool canTryQuick)
    26622672{
     2673  double startingTheta = paramData.startingTheta;
     2674  double & endingTheta = paramData.endingTheta;
    26632675     // stuff is already at starting
    26642676     // For this crude version just try and go to end
     
    27572769               reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
    27582770          } else {
    2759                whileIterating(startingTheta,  endingTheta, reportIncrement,
    2760                               lowerChange, upperChange,
     2771               whileIterating(paramData, reportIncrement,
    27612772                              changeObjective);
    27622773               startingTheta = endingTheta;
     
    28042815  assert (ratio==1||ratio==2);
    28052816  // allow for unscaled - even if not needed
    2806   int lengthArrays = 4*numberTotal+(numberTotal+2)*ratio;
     2817  int lengthArrays = 4*numberTotal+(2*numberTotal+2)*ratio;
    28072818  /*
    28082819    Save information and modify
     
    28192830  char * markDone = reinterpret_cast<char *>(lowerList+numberTotal);
    28202831  memset(markDone,0,numberTotal);
     2832  int * backwardBasic = lowerList+2*numberTotal;
     2833  parametricsData paramData;
     2834  paramData.lowerChange = lowerChange;
     2835  paramData.lowerList=lowerList;
     2836  paramData.upperChange = upperChange;
     2837  paramData.upperList=upperList;
     2838  paramData.markDone=markDone;
     2839  paramData.backwardBasic=backwardBasic;
    28212840  // Find theta when bounds will cross over and create arrays
    28222841  memset(lowerChange, 0, numberTotal * sizeof(double));
     
    30363055        bool canSkipFactorization=true;
    30373056        while (!returnCode) {
    3038           returnCode = parametricsLoop(startingTheta, endingTheta,
     3057                 paramData.startingTheta=startingTheta;
     3058                 paramData.endingTheta=endingTheta;
     3059                 returnCode = parametricsLoop(paramData,
    30393060                                       data,canSkipFactorization);
     3061                 startingTheta=paramData.startingTheta;
     3062                 endingTheta=paramData.endingTheta;
    30403063          canSkipFactorization=false;
    30413064          if (!returnCode) {
     
    30643087      ray_ = new double [numberColumns_];
    30653088    }
    3066     reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(1);
    30673089    if (swapped&&lower_) {
    30683090      double * temp=saveLower;
     
    30733095      upper_=temp;
    30743096    }
     3097    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
    30753098  }   
    30763099  if (!scalingFlag_) {
     
    31173140}
    31183141int
    3119 ClpSimplexOther::parametricsLoop(double & startingTheta, double & endingTheta,
     3142ClpSimplexOther::parametricsLoop(parametricsData & paramData,
    31203143                                 ClpDataSave & data,bool canSkipFactorization)
    31213144{
     3145  double & startingTheta = paramData.startingTheta;
     3146  double & endingTheta = paramData.endingTheta;
    31223147  int numberTotal = numberRows_+numberColumns_;
    3123   const double * lowerChange = lower_+numberTotal;
    3124   const double * upperChange = upper_+numberTotal;
    31253148  // stuff is already at starting
    3126   int * lowerList = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
    3127   int * upperList = (reinterpret_cast<int *>(upper_+4*numberTotal))+2;
     3149  const int * lowerList = paramData.lowerList;
     3150  const int * upperList = paramData.upperList;
    31283151  problemStatus_ = -1;
    31293152  //double saveEndingTheta=endingTheta;
     
    32593282    // Do iterations
    32603283    problemStatus_=-1;
    3261     whileIterating(startingTheta,  endingTheta, 0.0,
    3262                    lowerChange, upperChange,
     3284    whileIterating(paramData, 0.0,
    32633285                   NULL);
    32643286    //startingTheta = endingTheta;
     
    34173439*/
    34183440int
    3419 ClpSimplexOther::whileIterating(double & startingTheta, double & endingTheta, double /*reportIncrement*/,
    3420                                 const double * lowerChange, const double * upperChange,
     3441ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
    34213442                                const double * /*changeObjective*/)
    34223443{
     3444  double & startingTheta = paramData.startingTheta;
     3445  double & endingTheta = paramData.endingTheta;
     3446  const double * lowerChange = paramData.lowerChange;
     3447  const double * upperChange = paramData.upperChange;
    34233448  int numberTotal = numberColumns_ + numberRows_;
    3424   const int * lowerList = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
    3425   const int * upperList = (reinterpret_cast<int *>(upper_+4*numberTotal))+2;
     3449  const int * lowerList = paramData.lowerList;
     3450  const int * upperList = paramData.upperList;
     3451  // do basic pointers
     3452  int * backwardBasic = paramData.backwardBasic;
     3453  for (int i=0;i<numberTotal;i++)
     3454    backwardBasic[i]=-1;
     3455  for (int i=0;i<numberRows_;i++) {
     3456    int iPivot=pivotVariable_[i];
     3457    backwardBasic[iPivot]=i;
     3458  }
    34263459     {
    34273460          int i;
     
    34463479          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
    34473480          // Get theta for bounds - we know can't crossover
    3448           int pivotType = nextTheta(1, increaseTheta,
    3449                                     lowerChange, upperChange, NULL);
     3481          int pivotType = nextTheta(1, increaseTheta, paramData,
     3482                                    NULL);
    34503483          useTheta += theta_;
    34513484          double change = useTheta - lastTheta;
     
    35913624                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
    35923625                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
     3626#if CLP_CAN_HAVE_ZERO_OBJ
    35933627                    if ((specialOptions_&2097152)==0) {
     3628#endif
    35943629                      nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
    35953630                                                                                               rowArray_[2], theta_,
    35963631                                                                                               objectiveChange, false);
    35973632                      assert (!nswapped);
     3633#if CLP_CAN_HAVE_ZERO_OBJ
    35983634                    } else {
    35993635                      rowArray_[0]->clear();
     
    36013637                      columnArray_[0]->clear();
    36023638                    }
     3639#endif
    36033640                    // which will change basic solution
    36043641                    if (nswapped) {
     
    37103747                    }
    37113748                    // update primal solution
     3749#if CLP_CAN_HAVE_ZERO_OBJ
     3750                    if ((specialOptions_&2097152)!=0)
     3751                      theta_=0.0;
     3752#endif
    37123753                    if (theta_ < 0.0) {
    37133754#ifdef CLP_DEBUG
     
    37373778                    }
    37383779                    objectiveChange = 0.0;
     3780#if CLP_CAN_HAVE_ZERO_OBJ
    37393781                    if ((specialOptions_&2097152)==0) {
     3782#endif
    37403783                      for (int i=0;i<numberTotal;i++)
    37413784                        objectiveChange += solution_[i]*cost_[i];
    37423785                      objectiveChange -= objectiveValue_;
     3786#if CLP_CAN_HAVE_ZERO_OBJ
    37433787                    }
     3788#endif
    37443789                    // outgoing
    37453790                    originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
     
    37503795                         valueOut_ = lowerOut_;
    37513796                         dj_[sequenceOut_] = theta_;
     3797#if CLP_CAN_HAVE_ZERO_OBJ>1
     3798#ifdef COIN_REUSE_RANDOM
     3799                         if ((specialOptions_&2097152)!=0) {
     3800                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
     3801                         }
     3802#endif
     3803#endif
    37523804                    } else {
    37533805                         valueOut_ = upperOut_;
    37543806                         dj_[sequenceOut_] = -theta_;
     3807#if CLP_CAN_HAVE_ZERO_OBJ>1
     3808#ifdef COIN_REUSE_RANDOM
     3809                         if ((specialOptions_&2097152)!=0) {
     3810                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
     3811                         }
     3812#endif
     3813#endif
    37553814                    }
    37563815                    solution_[sequenceOut_] = valueOut_;
    37573816                    int whatNext = housekeeping(objectiveChange);
     3817                    assert (backwardBasic[sequenceOut_]==pivotRow_);
     3818                    backwardBasic[sequenceOut_]=-1;
     3819                    backwardBasic[sequenceIn_]=pivotRow_;
    37583820                    {
    37593821                      char in[200],out[200];
     
    38283890                   theta_ = 0.0;
    38293891                   //adjust [4] from handler - but
    3830                    rowArray_[4]->clear(); // temp
     3892                   //rowArray_[4]->clear(); // temp
    38313893                   if (action>=0&&action<10)
    38323894                     problemStatus_=-1; // carry on
     
    40634125// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
    40644126int
    4065 ClpSimplexOther::nextTheta(int type, double maxTheta,
    4066                            const double * lowerChange, const double * upperChange,
     4127ClpSimplexOther::nextTheta(int type, double maxTheta, parametricsData & paramData,
    40674128                           const double * /*changeObjective*/)
    40684129{
    4069   int numberTotal = numberColumns_ + numberRows_;
    4070   const int * lowerList = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
    4071   const int * upperList = (reinterpret_cast<int *>(upper_+4*numberTotal))+2;
     4130  const double * lowerChange = paramData.lowerChange;
     4131  const double * upperChange = paramData.upperChange;
     4132  const int * lowerList = paramData.lowerList;
     4133  const int * upperList = paramData.upperList;
    40724134  int iSequence;
    40734135  theta_ = maxTheta;
     
    40814143  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    40824144  const double * elementByColumn = matrix_->getElements();
     4145#if 0
     4146  double tempArray[5000];
     4147  bool checkIt=false;
     4148  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
     4149    memcpy(tempArray,array,numberRows_*sizeof(double));
     4150    checkIt=true;
     4151    needFullUpdate=true;
     4152  }
     4153#endif
    40834154  if (!factorization_->pivots()||needFullUpdate) {
    40844155    rowArray_[4]->clear();
     
    41874258    // ftran it
    41884259    factorization_->updateColumn(rowArray_[0], rowArray_[4]);
    4189   } else {
    4190     assert (sequenceIn_>=0);
     4260#if 0
     4261    if (checkIt) {
     4262      for (int i=0;i<numberRows_;i++) {
     4263        assert (fabs(tempArray[i]-array[i])<1.0e-8);
     4264      }
     4265    }
     4266#endif
     4267  } else if (sequenceIn_>=0) {
     4268    //assert (sequenceIn_>=0);
    41914269    assert (sequenceOut_>=0);
    41924270    assert (sequenceIn_!=sequenceOut_);
     
    42604338  const int * index = rowArray_[4]->getIndices();
    42614339  int number = rowArray_[4]->getNumElements();
    4262   int * lowerList2 = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
    4263   char * markDone = reinterpret_cast<char *>(lowerList2+numberTotal);
    4264 #if 0
    4265   for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
    4266     //int iPivot = index[iRow];
    4267     iSequence = pivotVariable_[iPivot];
    4268     // solution value will be sol - theta*alpha
    4269     // bounds will be bounds + change *theta
    4270     double currentSolution = solution_[iSequence];
    4271     double currentLower = lower_[iSequence];
    4272     double currentUpper = upper_[iSequence];
    4273     double alpha = array[iPivot];
    4274     assert (currentSolution >= currentLower - 100.0*primalTolerance_);
    4275     assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
    4276     double thetaCoefficient;
    4277     double hitsLower = COIN_DBL_MAX;
    4278     thetaCoefficient = lowerChange[iSequence] + alpha;
    4279     if (thetaCoefficient > 1.0e-8)
    4280       hitsLower = (currentSolution - currentLower) / thetaCoefficient;
    4281     //if (hitsLower < 0.0) {
    4282     // does not hit - but should we check further
    4283     //   hitsLower = COIN_DBL_MAX;
    4284     //}
    4285     double hitsUpper = COIN_DBL_MAX;
    4286     thetaCoefficient = upperChange[iSequence] + alpha;
    4287     if (thetaCoefficient < -1.0e-8)
    4288       hitsUpper = (currentSolution - currentUpper) / thetaCoefficient;
    4289     //if (hitsUpper < 0.0) {
    4290     // does not hit - but should we check further
    4291     //   hitsUpper = COIN_DBL_MAX;
    4292     //}
    4293     if (CoinMin(hitsLower, hitsUpper) < theta_) {
    4294       theta_ = CoinMin(hitsLower, hitsUpper);
    4295       toLower = hitsLower < hitsUpper;
    4296       pivotRow_ = iPivot;
    4297     }
    4298   }
    4299 #else
     4340  char * markDone = paramData.markDone;
     4341  const int * backwardBasic = paramData.backwardBasic;
    43004342  // first ones with alpha
    43014343  for (int i=0;i<number;i++) {
     
    43074349    // bounds will be bounds + change *theta
    43084350    double currentSolution = solution_[iSequence];
    4309     double currentLower = lower_[iSequence];
    4310     double currentUpper = upper_[iSequence];
    43114351    double alpha = array[iPivot];
    4312     assert (currentSolution >= currentLower - 100.0*primalTolerance_);
    4313     assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
    4314     double thetaCoefficient;
    4315     thetaCoefficient = lowerChange[iSequence] + alpha;
    4316     if (thetaCoefficient > 1.0e-8) {
     4352    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
     4353    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
     4354    if (thetaCoefficientLower > 1.0e-8) {
     4355      double currentLower = lower_[iSequence];
     4356      assert (currentSolution >= currentLower - 100.0*primalTolerance_);
    43174357      double gap=currentSolution-currentLower;
    4318       if (thetaCoefficient*theta_>gap) {
    4319         theta_ = gap/thetaCoefficient;
     4358      if (thetaCoefficientLower*theta_>gap) {
     4359        theta_ = gap/thetaCoefficientLower;
    43204360        toLower=true;
    43214361        pivotRow_=iPivot;
    43224362      }
    43234363    }
    4324     thetaCoefficient = upperChange[iSequence] + alpha;
    4325     if (thetaCoefficient < -1.0e-8) {
     4364    if (thetaCoefficientUpper < -1.0e-8) {
     4365      double currentUpper = upper_[iSequence];
     4366      assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
    43264367      double gap=currentSolution-currentUpper; //negative
    4327       if (thetaCoefficient*theta_<gap) {
    4328         theta_ = gap/thetaCoefficient;
     4368      if (thetaCoefficientUpper*theta_<gap) {
     4369        theta_ = gap/thetaCoefficientUpper;
    43294370        toLower=false;
    43304371        pivotRow_=iPivot;
     
    43464387          theta_ = gap/thetaCoefficient;
    43474388          toLower=true;
    4348           pivotRow_ = -2-iSequence;
     4389          pivotRow_ = backwardBasic[iSequence];
    43494390        }
    43504391      }
     
    43644405          theta_ = gap/thetaCoefficient;
    43654406          toLower=false;
    4366           pivotRow_ = -2-iSequence;
     4407          pivotRow_ = backwardBasic[iSequence];
    43674408        }
    43684409      }
    43694410    }
    43704411  }
    4371   if (pivotRow_<-1) {
    4372     // find
    4373     int iSequence = -pivotRow_-2;
    4374     for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
    4375       if (iSequence == pivotVariable_[iPivot]) {
    4376         pivotRow_=iPivot;
    4377         break;
    4378       }
    4379     }
    4380     assert (pivotRow_>=0);
    4381   }
    4382 #endif
    43834412  theta_ = CoinMax(theta_,0.0);
    43844413  if (theta_>1.0e-15) {
     
    67076736  double * weights = new double [numberLook+numberColumns_];
    67086737  double * columnWeights = weights+numberLook;
     6738#ifndef COIN_REUSE_RANDOM
    67096739  coin_init_random_vec(columnWeights,numberColumns_);
     6740#else
     6741  for (int i=0;i<numberColumns_;i++)
     6742    columnWeights[i]=CoinDrand48();
     6743#endif
    67106744  // get row copy
    67116745  CoinPackedMatrix rowCopy = *matrix();
     
    69616995        }
    69626996      }
    6963       if (sum==0.0||sum>=lastSum)
     6997      if (sum==0.0||sum>=lastSum-1.0e-8)
    69646998        break;
    69656999      lastSum=sum;
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1780 r1785  
    103103    /// Finds best possible pivot
    104104    double bestPivot(bool justColumns=false);
     105  typedef struct {
     106    double startingTheta;
     107    double endingTheta;
     108    double * lowerChange; // full array of lower bound changes
     109    int * lowerList; // list of lower bound changes
     110    double * upperChange; // full array of upper bound changes
     111    int * upperList; // list of upper bound changes
     112    char * markDone; // mark which ones looked at
     113    int * backwardBasic; // from sequence to pivot row
     114  } parametricsData;
    105115
    106116private:
     
    113123         if event handler exists it may do more
    114124     */
    115      int parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
     125     int parametricsLoop(parametricsData & paramData, double reportIncrement,
    116126                         const double * changeLower, const double * changeUpper,
    117127                         const double * changeObjective, ClpDataSave & data,
    118128                         bool canTryQuick);
    119      int parametricsLoop(double & startingTheta, double & endingTheta,
     129     int parametricsLoop(parametricsData & paramData,
    120130                         ClpDataSave & data,bool canSkipFactorization=false);
    121131     /**  Refactorizes if necessary
     
    137147         +3 max iterations
    138148      */
    139      int whileIterating(double & startingTheta, double & endingTheta, double reportIncrement,
    140                         const double * changeLower, const double * changeUpper,
     149     int whileIterating(parametricsData & paramData, double reportIncrement,
    141150                        const double * changeObjective);
    142151     /** Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none).
     
    144153         type 1 bounds, 2 objective, 3 both.
    145154     */
    146      int nextTheta(int type, double maxTheta,
    147                    const double * changeLower, const double * changeUpper,
     155     int nextTheta(int type, double maxTheta, parametricsData & paramData,
    148156                   const double * changeObjective);
    149157     /// Restores bound to original bound
     
    197205         non-zero return code indicates minor problems
    198206     */
    199      int restoreFromDual(const ClpSimplex * dualProblem);
     207  int restoreFromDual(const ClpSimplex * dualProblem,
     208                      bool checkAccuracy=false);
    200209     /** Does very cursory presolve.
    201210         rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
Note: See TracChangeset for help on using the changeset viewer.