Ignore:
Timestamp:
Sep 4, 2008 4:46:53 PM (11 years ago)
Author:
forrest
Message:

trying to make faster

File:
1 edited

Legend:

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

    r1243 r1266  
    295295    numberFake_ =0; // Number of variables at fake bounds
    296296    numberChanged_ =0; // Number of variables with changed costs
    297     changeBounds(true,NULL,objectiveChange);
     297    changeBounds(1,NULL,objectiveChange);
    298298   
    299299    if (!ifValuesPass) {
     
    357357  finish(startFinishOptions);
    358358}
     359#ifdef CLP_INVESTIGATE
     360static int z_reason[7]={0,0,0,0,0,0,0};
     361static int z_thinks=-1;
     362#endif
    359363void
    360364ClpSimplexDual::gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
    361365                           ClpDataSave & data)
    362366{
     367#ifdef CLP_INVESTIGATE
     368  z_reason[0]++;
     369  z_thinks=-1;
     370  int nPivots=9999;
     371#endif
    363372  int lastCleaned=0; // last time objective or bounds cleaned up
    364373 
     
    487496    // Do iterations
    488497    int returnCode=whileIterating(saveDuals,ifValuesPass);
     498#ifdef CLP_INVESTIGATE
     499    nPivots=factorization_->pivots();
     500#endif
    489501    if (returnCode==-2)
    490502      factorType=3;
    491503  }
     504#ifdef CLP_INVESTIGATE
     505  if (z_thinks!=-1) {
     506    assert (z_thinks<4);
     507    if ((!factorization_->pivots()&&nPivots<20)&&z_thinks>=0&&z_thinks<2)
     508      z_thinks += 4;
     509    z_reason[1+z_thinks]++;
     510  }
     511  if ((z_reason[0]%1000)==0) {
     512    printf("Reason");
     513    for (int i=0;i<7;i++)
     514      printf(" %d",z_reason[i]);
     515    printf("\n");
     516  }
     517#endif
    492518}
    493519int
     
    495521{
    496522  algorithm_ = -1;
    497  
     523  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
    498524  // save data
    499525  ClpDataSave data = saveData();
     
    618644    numberFake_ =0; // Number of variables at fake bounds
    619645    numberChanged_ =0; // Number of variables with changed costs
    620     changeBounds(true,NULL,objectiveChange);
     646    changeBounds(1,NULL,objectiveChange);
    621647   
    622648    int lastCleaned=0; // last time objective or bounds cleaned up
     
    772798ClpSimplexDual::whileIterating(double * & givenDuals,int ifValuesPass)
    773799{
     800#ifdef CLP_INVESTIGATE
     801  z_thinks=-1;
     802#endif
    774803#if 0
    775804  if (!numberIterations_&&auxiliaryModel_) {
     
    13471376                                                         rowArray_[1],
    13481377                                                         pivotRow_,
    1349                                                          alpha_);
     1378                                                         alpha_,
     1379                                                     (moreSpecialOptions_&16)!=0);
    13501380        // if no pivots, bad update but reasonable alpha - take and invert
    13511381        if (updateStatus==2&&
     
    13621392          dualRowPivot_->unrollWeights();
    13631393          // later we may need to unwind more e.g. fake bounds
    1364           if (factorization_->pivots()) {
     1394          if (factorization_->pivots()&&
     1395              ((moreSpecialOptions_&16)==0||factorization_->pivots()>4)) {
    13651396            problemStatus_=-2; // factorize now
    13661397            returnCode=-2;
     1398            moreSpecialOptions_ |= 16;
    13671399            break;
    13681400          } else {
     
    15471579        }
    15481580      } else {
     1581#ifdef CLP_INVESTIGATE
     1582        z_thinks=1;
     1583#endif
    15491584        // no incoming column is valid
    15501585        pivotRow_=-1;
     
    15531588          printf("** no column pivot\n");
    15541589#endif
    1555         if (factorization_->pivots()<=CoinMax(dontFactorizePivots_,5)&&acceptablePivot_<=1.0e-8) {
     1590        if (factorization_->pivots()<2&&acceptablePivot_<=1.0e-8) {
    15561591          //&&goodAccuracy()) {
    15571592          // If not in branch and bound etc save ray
     
    15651600          // If we have just factorized and infeasibility reasonable say infeas
    15661601          if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
    1567             if (valueOut_>upperOut_+1.0e-3||valueOut_<lowerOut_-1.0e-3
     1602            if (valueOut_>upperOut_+1.0e-4||valueOut_<lowerOut_-1.0e-4
    15681603                || (specialOptions_&64)==0) {
    15691604              // say infeasible
     
    16611696      }
    16621697    } else {
     1698#ifdef CLP_INVESTIGATE
     1699      z_thinks=0;
     1700#endif
    16631701      // no pivot row
    16641702#ifdef CLP_DEBUG
     
    16971735      returnCode=0;
    16981736      if (numberPivots<=CoinMax(dontFactorizePivots_,20)&&
    1699           (specialOptions_&2048)!=0&&!numberChanged_
    1700           &&dualBound_>1.0e8) {
     1737          (specialOptions_&2048)!=0&&(true||!numberChanged_||perturbation_==101)
     1738          &&dualBound_>=1.0e8) {
    17011739        specialCase=true;
    17021740        // as dual bound high - should be okay
     
    17391777            problemStatus_=-5;
    17401778          } else {
    1741             if (numberPivots) {
     1779            problemStatus_=0;
     1780            // May be perturbed
     1781            if (perturbation_==101||numberChanged_) {
     1782              perturbation_=102; // stop any perturbations
     1783              //double changeCost;
     1784              //changeBounds(1,NULL,changeCost);
     1785              createRim4(false);
     1786              // make sure duals are current
     1787              computeDuals(givenDuals);
     1788              checkDualSolution();
     1789              if (numberDualInfeasibilities_) {
     1790                problemStatus_=-3;
     1791              } else {
     1792                computeObjectiveValue(true);
     1793              }
     1794            } else if (numberPivots) {
     1795              computeObjectiveValue(true);
     1796            }
     1797            if (numberPivots<-1000) {
    17421798              // objective may be wrong
    17431799              objectiveValue_ = innerProduct(cost_,numberColumns_+numberRows_,solution_);
     
    17751831              }
    17761832            }
    1777             problemStatus_=0;
    17781833            sumPrimalInfeasibilities_=0.0;
    1779             if ((specialOptions_&(1024+16384))!=0) {
     1834            if ((specialOptions_&(1024+16384))!=0&&!problemStatus_) {
    17801835              CoinIndexedVector * arrayVector = rowArray_[1];
    17811836              arrayVector->clear();
     
    24112466// and cost of change vector
    24122467int
    2413 ClpSimplexDual::changeBounds(bool initialize,
     2468ClpSimplexDual::changeBounds(int initialize,
    24142469                                 CoinIndexedVector * outputArray,
    24152470                                 double & changeCost)
     
    25082563    }
    25092564    return numberInfeasibilities;
    2510   } else {
     2565  } else if (initialize==1) {
    25112566    int iSequence;
    25122567     
     
    25592614
    25602615    return 1;
     2616  } else {
     2617    // just reset changed ones
     2618    if (columnScale_) {
     2619      int iSequence;
     2620      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     2621        FakeBound fakeStatus = getFakeBound(iSequence);
     2622        if (fakeStatus!=noFake) {
     2623          if (((int) fakeStatus&1)!=0) {
     2624            // lower
     2625            double value = columnLower_[iSequence];
     2626            if (value>-1.0e30) {
     2627              double multiplier = rhsScale_/columnScale_[iSequence];
     2628              value *= multiplier;
     2629            }
     2630            columnLowerWork_[iSequence]=value;
     2631          }
     2632          if (((int) fakeStatus&2)!=0) {
     2633            // upper
     2634            double value = columnUpper_[iSequence];
     2635            if (value<1.0e30) {
     2636              double multiplier = rhsScale_/columnScale_[iSequence];
     2637              value *= multiplier;
     2638            }
     2639            columnUpperWork_[iSequence]=value;
     2640          }
     2641        }
     2642      }
     2643      for (iSequence=0;iSequence<numberRows_;iSequence++) {
     2644        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
     2645        if (fakeStatus!=noFake) {
     2646          if (((int) fakeStatus&1)!=0) {
     2647            // lower
     2648            double value = rowLower_[iSequence];
     2649            if (value>-1.0e30) {
     2650              double multiplier = rhsScale_*rowScale_[iSequence];
     2651              value *= multiplier;
     2652            }
     2653            rowLowerWork_[iSequence]=value;
     2654          }
     2655          if (((int) fakeStatus&2)!=0) {
     2656            // upper
     2657            double value = rowUpper_[iSequence];
     2658            if (value<1.0e30) {
     2659              double multiplier = rhsScale_*rowScale_[iSequence];
     2660              value *= multiplier;
     2661            }
     2662            rowUpperWork_[iSequence]=value;
     2663          }
     2664        }
     2665      }
     2666    } else {
     2667      int iSequence;
     2668      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     2669        FakeBound fakeStatus = getFakeBound(iSequence);
     2670        if (((int) fakeStatus&1)!=0) {
     2671          // lower
     2672          columnLowerWork_[iSequence]=columnLower_[iSequence];
     2673        }
     2674        if (((int) fakeStatus&2)!=0) {
     2675          // upper
     2676          columnUpperWork_[iSequence]=columnUpper_[iSequence];
     2677        }
     2678      }
     2679      for (iSequence=0;iSequence<numberRows_;iSequence++) {
     2680        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
     2681        if (((int) fakeStatus&1)!=0) {
     2682          // lower
     2683          rowLowerWork_[iSequence]=rowLower_[iSequence];
     2684        }
     2685        if (((int) fakeStatus&2)!=0) {
     2686          // upper
     2687          rowUpperWork_[iSequence]=rowUpper_[iSequence];
     2688        }
     2689      }
     2690    }
     2691    return 0;
    25612692  }
    25622693}
     
    34323563                                      int ifValuesPass)
    34333564{
     3565#ifdef CLP_INVESTIGATE
     3566  if (z_thinks>0&&z_thinks<2)
     3567    z_thinks+=2;
     3568#endif
    34343569  // If lots of iterations then adjust costs if large ones
    34353570  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
     
    38674002      progress_.modifyObjective(objectiveValue_
    38684003                               -sumDualInfeasibilities_*dualBound_);
     4004      // see if cutoff reached
     4005      double limit = 0.0;
     4006      getDblParam(ClpDualObjectiveLimit, limit);
     4007      if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
     4008         limit+1.0e-7+1.0e-8*fabs(limit)&&!numberAtFakeBound()) {
     4009        //looks infeasible on objective
     4010        if (perturbation_==101) {
     4011          perturbation_=102; // stop any perturbations
     4012          cleanDuals=1;
     4013          // make sure fake bounds are back
     4014          changeBounds(1,NULL,changeCost);
     4015          createRim4(false);
     4016          // make sure duals are current
     4017          computeDuals(givenDuals);
     4018          checkDualSolution();
     4019        }
     4020      }
    38694021      if (primalFeasible()&&!givenDuals) {
    38704022        // may be optimal - or may be bounds are wrong
     
    38754027        double * saveColumnSolution = NULL;
    38764028        double * saveRowSolution = NULL;
    3877         bool inCbc = (specialOptions_&0x01000000)!=0;
     4029        bool inCbc = (specialOptions_&(0x01000000|16384))!=0;
    38784030        if (!inCbc) {
    38794031          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
    38804032          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
    38814033        }
    3882         numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
     4034        numberChangedBounds=changeBounds(0,rowArray_[3],changeCost);
    38834035        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
    38844036          //looks optimal - do we need to reset tolerance
     
    38884040            // make sure fake bounds are back
    38894041            //computeObjectiveValue();
    3890             changeBounds(true,NULL,changeCost);
     4042            changeBounds(1,NULL,changeCost);
    38914043            //computeObjectiveValue();
    38924044            createRim4(false);
     
    38944046            computeDuals(givenDuals);
    38954047            checkDualSolution();
    3896             //if (numberDualInfeasibilities_)
     4048#define DUAL_TRY_FASTER
     4049#ifdef DUAL_TRY_FASTER
     4050            if (numberDualInfeasibilities_) {
     4051#endif
    38974052              numberChanged_=1; // force something to happen
    3898             //else
    3899             //computeObjectiveValue();
     4053#ifdef DUAL_TRY_FASTER
     4054            } else {
     4055              //double value = objectiveValue_;
     4056              computeObjectiveValue(true);
     4057              //printf("old %g new %g\n",value,objectiveValue_);
     4058              //numberChanged_=1;
     4059            }
     4060#endif
    39004061          }
    39014062          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
     
    40004161      if (problemStatus_==-4||problemStatus_==-5) {
    40014162        // may be infeasible - or may be bounds are wrong
    4002         numberChangedBounds=changeBounds(false,NULL,changeCost);
     4163        numberChangedBounds=changeBounds(0,NULL,changeCost);
    40034164        /* Should this be here as makes no difference to being feasible.
    40044165           But seems to make a difference to run times. */
     
    40084169          numberChangedBounds=1;
    40094170          // make sure fake bounds are back
    4010           changeBounds(true,NULL,changeCost);
     4171          changeBounds(1,NULL,changeCost);
    40114172          createRim4(false);
    40124173        }
     
    40924253          */
    40934254          int saveNumberFake = numberFake_;
    4094           changeBounds(true,NULL,changeCost);
     4255          changeBounds(1,NULL,changeCost);
    40954256          numberFake_ += saveNumberFake;
    40964257          cleanDuals=2;
     
    42054366    }
    42064367  }
    4207   if (type==0||type==1) {
    4208     if (!type) {
    4209       // create save arrays
    4210       delete [] saveStatus_;
    4211       delete [] savedSolution_;
    4212       saveStatus_ = new unsigned char [numberRows_+numberColumns_];
    4213       savedSolution_ = new double [numberRows_+numberColumns_];
    4214     }
    4215     // save arrays
    4216     CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
    4217     CoinMemcpyN(rowActivityWork_,
    4218            numberRows_,savedSolution_+numberColumns_);
    4219     CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
    4220     // save extra stuff
    4221     int dummy;
    4222     matrix_->generalExpanded(this,5,dummy);
    4223   }
    4224   if (weightsSaved) {
    4225     // restore weights (if saved) - also recompute infeasibility list
    4226     if (tentativeStatus>-3)
    4227       dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
    4228     else
    4229       dualRowPivot_->saveWeights(this,3);
    4230   }
    42314368  // unflag all variables (we may want to wait a bit?)
    42324369  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
     
    42664403#endif
    42674404      }
     4405    }
     4406  }
     4407  if (problemStatus_<0) {
     4408    if (type==0||type==1) {
     4409      if (!type) {
     4410        // create save arrays
     4411        delete [] saveStatus_;
     4412        delete [] savedSolution_;
     4413        saveStatus_ = new unsigned char [numberRows_+numberColumns_];
     4414        savedSolution_ = new double [numberRows_+numberColumns_];
     4415      }
     4416      // save arrays
     4417      CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
     4418      CoinMemcpyN(rowActivityWork_,
     4419                  numberRows_,savedSolution_+numberColumns_);
     4420      CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
     4421      // save extra stuff
     4422      int dummy;
     4423      matrix_->generalExpanded(this,5,dummy);
     4424    }
     4425    if (weightsSaved) {
     4426      // restore weights (if saved) - also recompute infeasibility list
     4427      if (tentativeStatus>-3)
     4428        dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
     4429      else
     4430        dualRowPivot_->saveWeights(this,3);
    42684431    }
    42694432  }
     
    43554518  }
    43564519  lastGoodIteration_ = numberIterations_;
     4520  if (numberIterations_>lastBadIteration_+100)
     4521    moreSpecialOptions_ &= ~16; // clear check accuracy flag
    43574522  if (problemStatus_<0) {
    43584523    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
     
    46014766  }
    46024767  int iColumn;
     4768  const int * columnLength = matrix_->getVectorLengths();
    46034769  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    46044770    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
    4605       int length = matrix_->getVectorLength(iColumn);
     4771      int length = columnLength[iColumn];
    46064772      if (length>2) {
    46074773        maxLength = CoinMax(maxLength,length);
     
    46174783  }
    46184784  bool uniformChange=false;
     4785  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
    46194786  if (perturbation_>50) {
    46204787    // Experiment
    46214788    // maximumFraction could be 1.0e-10 to 1.0
    46224789    double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
    4623     maximumFraction = m[CoinMin(perturbation_-51,10)];
     4790    int whichOne = perturbation_-51;
     4791    //if (inCbcOrOther&&whichOne>0)
     4792    //whichOne--;
     4793    maximumFraction = m[CoinMin(whichOne,10)];
     4794  } else if (inCbcOrOther) {
     4795    //maximumFraction = 1.0e-6;
    46244796  }
    46254797  int iRow;
     
    48114983      }
    48124984      if (value) {
    4813         int length = matrix_->getVectorLength(iColumn);
     4985        int length = columnLength[iColumn];
    48144986        if (length>3) {
    48154987          length = (int) ((double) length * factor);
     
    52055377  numberFake_ =0; // Number of variables at fake bounds
    52065378  numberChanged_ =0; // Number of variables with changed costs
    5207   changeBounds(true,NULL,objectiveChange);
     5379  changeBounds(1,NULL,objectiveChange);
    52085380
    52095381  problemStatus_ = -1;
     
    53305502  // Say not in fast dual
    53315503  specialOptions_ &= ~16384;
    5332   assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>1.0e8)
     5504  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>=1.0e8)
    53335505         ||returnCode||problemStatus_); // all bounds should be okay
     5506  if (numberFake_>0) {
     5507    // Set back
     5508    double dummy;
     5509    changeBounds(2,NULL,dummy);
     5510  }
    53345511  // Restore any saved stuff
    53355512  restoreData(data);
     
    59026079  // Get a row copy in standard format
    59036080  CoinPackedMatrix copy;
     6081  copy.setExtraGap(0.0);
     6082  copy.setExtraMajor(0.0);
    59046083  copy.reverseOrderedCopyOf(*columnCopy);
    59056084  // get matrix data pointers
     
    60636242  createRim1(false);
    60646243  double dummyChangeCost=0.0;
    6065   changeBounds(true,rowArray_[2],dummyChangeCost);
     6244  changeBounds(1,rowArray_[2],dummyChangeCost);
    60666245  // throw away change
    60676246  for (int i=0;i<4;i++)
Note: See TracChangeset for help on using the changeset viewer.