Changeset 1678 for trunk/Clp


Ignore:
Timestamp:
Jan 25, 2011 5:58:03 AM (9 years ago)
Author:
forrest
Message:

gub stuff

Location:
trunk/Clp/src
Files:
6 edited

Legend:

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

    r1677 r1678  
    177177     startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
    178178     if (!numberGubColumns_) {
    179           startColumn_ = new CoinBigIndex [1];
    180           startColumn_[0] = 0;
     179       //startColumn_ = new CoinBigIndex [1];
     180       startColumn_[0] = 0;
    181181     }
    182182     CoinBigIndex numberElements = startColumn_[numberGubColumns_];
     
    235235     numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
    236236     matrix_ = originalMatrix;
     237     //delete originalMatrixA;
    237238     flags_ &= ~1;
    238239     // resize model (matrix stays same)
     240     // modify frequency
     241     if (frequency>=50)
     242       frequency = 50+(frequency-50)/2;
    239243     int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1;
    240244     model->resize(newRowSize, numberNeeded);
     
    903907#endif
    904908#endif
    905      return 0;
     909     if (numberStaticRows_+numberActiveSets_<model->numberRows())
     910       return 0;
     911     else
     912       return 1;
    906913}
    907914/*
     
    22752282          } else {
    22762283            // solo key
    2277             bool needKey;
     2284            bool needKey=false;
    22782285            if (numberActive) {
    22792286              if (whichKey<maximumGubColumns_) {
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1665 r1678  
    40134013     // Say no gaps
    40144014     flags_ &= ~2;
     4015     if (type_>=10)
     4016       return true; // gub
    40154017     if (check == 14 || check == 10) {
    40164018          if (matrix_->getNumElements() < columnStart[numberColumns]) {
  • trunk/Clp/src/ClpSimplex.cpp

    r1677 r1678  
    18281828
    18291829     // Update hidden stuff e.g. effective RHS and gub
    1830      matrix_->updatePivot(this, oldIn, oldOut);
     1830     int invertNow=matrix_->updatePivot(this, oldIn, oldOut);
    18311831     objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
    18321832     if (handler_->logLevel() > 7) {
     
    19851985          }
    19861986          return 1;
    1987      } else if (factorization_->timeToRefactorize() && !dontInvert) {
     1987     } else if ((factorization_->timeToRefactorize() && !dontInvert)
     1988                ||invertNow) {
    19881989          //printf("ret after %d pivots\n",factorization_->pivots());
    19891990          return 1;
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1677 r1678  
    44534453// Sets basis from original
    44544454void
    4455 ClpSimplexOther::setGubBasis(const ClpSimplex &original,const int * whichRows,
     4455ClpSimplexOther::setGubBasis(ClpSimplex &original,const int * whichRows,
    44564456                             const int * whichColumns)
    44574457{
     
    44704470  //assert (firstOdd==numberNormal);
    44714471  double * solution = primalColumnSolution();
    4472   const double * originalSolution = original.primalColumnSolution();
     4472  double * originalSolution = original.primalColumnSolution();
    44734473  const double * upperSet = gubMatrix->upperSet();
     4474  // Column copy of GUB part
    44744475  int numberSets = gubMatrix->numberSets();
    44754476  const int * startSet = gubMatrix->startSets();
    4476   const CoinBigIndex * startColumn = gubMatrix->startColumn();
     4477  const CoinBigIndex * columnStart = gubMatrix->startColumn();
    44774478  const double * columnLower = gubMatrix->columnLower();
     4479#ifdef TRY_IMPROVE
     4480  const double * columnUpper = gubMatrix->columnUpper();
     4481  const double * lowerSet = gubMatrix->lowerSet();
     4482  const double * element = gubMatrix->element();
     4483  const int * row = gubMatrix->row();
     4484  bool allPositive=true;
     4485  double * rowActivity = new double[numberNonGub];
     4486  memset(rowActivity, 0, numberNonGub*sizeof(double));
     4487  {
     4488    // Non gub contribution
     4489    const double * element = matrix_->getElements();
     4490    const int * row = matrix_->getIndices();
     4491    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     4492    const int * columnLength = matrix_->getVectorLengths();
     4493    for (int i=0;i<numberNormal;i++) {
     4494      int iColumn = whichColumns[i];
     4495      double value = originalSolution[iColumn];
     4496      if (value) {
     4497        for (CoinBigIndex j = columnStart[i];
     4498             j < columnStart[i] + columnLength[i]; j++) {
     4499          int iRow = row[j];
     4500          rowActivity[iRow] += value*element[j];
     4501        }
     4502      }
     4503    }
     4504  }
     4505  double * newSolution = new double [numberGubColumns];
     4506  int * slacks = new int [numberSets];
     4507  for (int i=0;i<numberSets;i++) {
     4508    double sum=0.0;
     4509    int iSlack=-1;
     4510    for (int j=startSet[i];j<startSet[i+1];j++) {
     4511      gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
     4512      int iColumn = whichColumns[j+numberNormal];
     4513      if (iColumn<numberColumns) {
     4514        columnIsGub[iColumn] = whichRows[numberNonGub+i];
     4515        double value = originalSolution[iColumn];
     4516        sum += value;
     4517        newSolution[j]=value;
     4518        for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
     4519          int iRow = row[k];
     4520          rowActivity[iRow] += value*element[k];
     4521          if (element[k] < 0.0)
     4522            allPositive=false;
     4523        }
     4524        if (columnStart[j]==columnStart[j+1])
     4525          iSlack=j;
     4526      } else {
     4527        newSolution[j]=0.0;
     4528        iSlack=j;
     4529        allPositive=false; // for now
     4530      }
     4531    }
     4532    slacks[i]=iSlack;
     4533    if (sum>upperSet[i]+1.0e-8) {
     4534      double gap = sum-upperSet[i];
     4535      if (iSlack>=0) {
     4536        double value=newSolution[iSlack];
     4537        if (value>0.0) {
     4538          double down = CoinMin(gap,value);
     4539          gap -= down;
     4540          sum -= down;
     4541          newSolution[iSlack] = value-down;
     4542        }
     4543      }
     4544      if (gap>1.0e-8) {
     4545        for (int j=startSet[i];j<startSet[i+1];j++) {
     4546          int iColumn = whichColumns[j+numberNormal];
     4547          if (newSolution[j]>0.0&&iColumn<numberColumns) {
     4548            double value = newSolution[j];
     4549            double down = CoinMin(gap,value);
     4550            gap -= down;
     4551            sum -= down;
     4552            newSolution[iSlack] = value-down;
     4553            for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
     4554              int iRow = row[k];
     4555              rowActivity[iRow] -= down*element[k];
     4556            }
     4557          }
     4558        }
     4559      }
     4560      assert (gap<1.0e-8);
     4561    } else if (sum<lowerSet[i]-1.0e-8) {
     4562      double gap = lowerSet[i]-sum;
     4563      if (iSlack>=0) {
     4564        double value=newSolution[iSlack];
     4565        if (value<columnUpper[iSlack]) {
     4566          double up = CoinMin(gap,columnUpper[iSlack]-value);
     4567          gap -= up;
     4568          sum += up;
     4569          newSolution[iSlack] = value+up;
     4570        }
     4571      }
     4572      if (gap>1.0e-8) {
     4573        for (int j=startSet[i];j<startSet[i+1];j++) {
     4574          int iColumn = whichColumns[j+numberNormal];
     4575          if (newSolution[j]<columnUpper[j]&&iColumn<numberColumns) {
     4576            double value = newSolution[j];
     4577            double up = CoinMin(gap,columnUpper[j]-value);
     4578            gap -= up;
     4579            sum += up;
     4580            newSolution[iSlack] = value+up;
     4581            for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
     4582              int iRow = row[k];
     4583              rowActivity[iRow] += up*element[k];
     4584            }
     4585          }
     4586        }
     4587      }
     4588      assert (gap<1.0e-8);
     4589    }
     4590    if (fabs(sum-upperSet[i])>1.0e-7)
     4591      printf("Sum for set %d is %g - lower %g, upper %g\n",i,
     4592             sum,lowerSet[i],upperSet[i]);
     4593  }
     4594  if (allPositive) {
     4595    // See if we can improve solution
     4596    // first reduce if over
     4597    double * gaps = new double [numberNonGub];
     4598    double direction = optimizationDirection_;
     4599    const double * cost = gubMatrix->cost();
     4600    bool over=false;
     4601    for (int i=0;i<numberNonGub;i++) {
     4602      double activity = rowActivity[i];
     4603      gaps[i]=0.0;
     4604      if (activity>rowUpper_[i]+1.0e-6) {
     4605        gaps[i]=activity-rowUpper_[i];
     4606        over=true;
     4607      }
     4608    }
     4609    double * weights = new double [numberGubColumns];
     4610    int * which = new int [numberGubColumns];
     4611    int * whichSet = new int [numberGubColumns];
     4612    if (over) {
     4613      int n=0;
     4614      for (int i=0;i<numberSets;i++) {
     4615        int iSlack = slacks[i];
     4616        if (iSlack<0||newSolution[iSlack]>upperSet[i]-1.0e-8)
     4617          continue;
     4618        double slackCost = cost[iSlack]*direction;
     4619        for (int j=startSet[i];j<startSet[i+1];j++) {
     4620          whichSet[j]=i;
     4621          double value = newSolution[j];
     4622          double thisCost = cost[j]*direction;
     4623          if (value>columnLower[j]&&j!=iSlack) {
     4624            if(thisCost<slackCost) {
     4625              double sum = 1.0e-30;
     4626              for (CoinBigIndex k = columnStart[j];
     4627                   k < columnStart[j+1] ; k++) {
     4628                int iRow = row[k];
     4629                sum += gaps[iRow]*element[k];
     4630              }
     4631              which[n]=j;
     4632              // big drop and small difference in cost better
     4633              weights[n++]=(slackCost-thisCost)/sum;
     4634            } else {
     4635              // slack better anyway
     4636              double move = value-columnLower[j];
     4637              newSolution[iSlack]=CoinMin(upperSet[i],
     4638                                          newSolution[iSlack]+move);
     4639              newSolution[j]=columnLower[j];
     4640              for (CoinBigIndex k = columnStart[j];
     4641                   k < columnStart[j+1] ; k++) {
     4642                int iRow = row[k];
     4643                rowActivity[iRow] -= move*element[k];
     4644              }
     4645            }
     4646          }
     4647        }
     4648      }
     4649      // sort
     4650      CoinSort_2(weights,weights+n,which);
     4651      for (int i=0;i<n;i++) {
     4652        int j= which[i];
     4653        int iSet = whichSet[j];
     4654        int iSlack = slacks[iSet];
     4655        assert (iSlack>=0);
     4656        double move = 0.0;
     4657        for (CoinBigIndex k = columnStart[j];
     4658             k < columnStart[j+1] ; k++) {
     4659          int iRow = row[k];
     4660          if(rowActivity[iRow]-rowUpper_[iRow]>move*element[k]) {
     4661            move = (rowActivity[iRow]-rowUpper_[iRow])/element[k];
     4662          }
     4663        }
     4664        move=CoinMin(move,newSolution[j]-columnLower[j]);
     4665        if (move) {
     4666          newSolution[j] -= move;
     4667          newSolution[iSlack] += move;
     4668          for (CoinBigIndex k = columnStart[j];
     4669               k < columnStart[j+1] ; k++) {
     4670            int iRow = row[k];
     4671            rowActivity[iRow] -= move*element[k];
     4672          }
     4673        }
     4674      }
     4675    }
     4676    delete [] whichSet;
     4677    delete [] which;
     4678    delete [] weights;
     4679    delete [] gaps;
     4680    // redo original status!
     4681    for (int i=0;i<numberSets;i++) {
     4682      int numberBasic=0;
     4683      int numberNewBasic=0;
     4684      int j1=-1;
     4685      int j2=-1;
     4686      for (int j=startSet[i];j<startSet[i+1];j++) {
     4687        if (newSolution[j]>columnLower[j]) {
     4688          numberNewBasic++;
     4689          j2=j;
     4690        }
     4691        int iOrig = whichColumns[j+numberNormal];
     4692        if (iOrig<numberColumns) {
     4693          if (original.getColumnStatus(iOrig)!=ClpSimplex::atLowerBound) {
     4694            numberBasic++;
     4695            j1=j;
     4696          }
     4697        } else {
     4698          int iSet = iOrig - numberColumns;
     4699          int iRow = whichRows[iSet+numberNonGub];
     4700          if (original.getRowStatus(iRow)==ClpSimplex::basic) {
     4701            numberBasic++;
     4702            j1=j;
     4703            abort();
     4704          }
     4705        }
     4706      }
     4707      if (numberBasic==1&&numberNewBasic==1&&
     4708          j1!=j2) {
     4709        int iOrig1=whichColumns[j1+numberNormal];
     4710        int iOrig2=whichColumns[j2+numberNormal];
     4711        ClpSimplex::Status status1 = original.getColumnStatus(iOrig1);
     4712        ClpSimplex::Status status2 = original.getColumnStatus(iOrig2);
     4713        originalSolution[iOrig1] = newSolution[j1];
     4714        originalSolution[iOrig2] = newSolution[j2];
     4715        original.setColumnStatus(iOrig1,status2);
     4716        original.setColumnStatus(iOrig2,status1);
     4717      }
     4718    }
     4719  }
     4720  delete [] newSolution;
     4721  delete [] slacks;
     4722  delete [] rowActivity;
     4723#else
    44784724  for (int i=0;i<numberSets;i++) {
    44794725    for (int j=startSet[i];j<startSet[i+1];j++) {
    44804726      gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
    44814727      int iColumn = whichColumns[j+numberNormal];
    4482       if (iColumn<numberColumns)
     4728      if (iColumn<numberColumns) {
    44834729        columnIsGub[iColumn] = whichRows[numberNonGub+i];
     4730      }
    44844731    }
    44854732  }
     4733#endif
    44864734  int * numberKey = new int [numberRows];
    44874735  memset(numberKey,0,numberRows*sizeof(int));
     
    45304778      int iSet = iOrig - numberColumns;
    45314779      int iRow = whichRows[iSet+numberNonGub];
    4532       if (original.getRowStatus(iRow)==ClpSimplex::basic) {
     4780      if (original.getRowStatus(iRow)==ClpSimplex::basic
     4781#ifdef TRY_IMPROVE
     4782          ||newSolution[i]>columnLower[i]+1.0e-8
     4783#endif
     4784          ) {
    45334785        assert(numberKey[iRow]);
    45344786        if (numberKey[iRow]==1)
     
    45534805      int chosen=-1;
    45544806      for (int j=startSet[i];j<startSet[i+1];j++) {
    4555         int length=startColumn[j+1]-startColumn[j];
     4807        int length=columnStart[j+1]-columnStart[j];
    45564808        int iOrig = whichColumns[j+numberNormal];
    45574809        double value;
    45584810        if (iOrig<numberColumns) {
    4559           value=originalSolution[iOrig]-columnLower[j];
     4811#ifdef TRY_IMPROVE
     4812          value=newSolution[j]-columnLower[j];
     4813#else
     4814          value = originalSolution[iOrig]-columnLower[j];
     4815#endif
    45604816          if (value>upper)
    45614817            gubMatrix->setStatus(i,ClpSimplex::atLowerBound);
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1676 r1678  
    197197                             int factorizationFrequency=50);
    198198     /// Sets basis from original
    199      void setGubBasis(const ClpSimplex &original,const int * whichRows,
     199     void setGubBasis(ClpSimplex &original,const int * whichRows,
    200200                      const int * whichColumns);
    201201     /// Restores basis to original
  • trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp

    r1677 r1678  
    5757  }
    5858#endif
    59   if ((specialOptions_&2097152)!=0) {
    60     resolveGub((9*modelPtr_->numberRows())/10);
    61     return;
     59  if ((specialOptions_&2097152)!=0||(specialOptions_&4194304)!=0) {
     60    bool takeHint;
     61    OsiHintStrength strength;
     62    int algorithm = 0;
     63    getHintParam(OsiDoDualInInitial,takeHint,strength);
     64    if (strength!=OsiHintIgnore)
     65      algorithm = takeHint ? -1 : 1;
     66    if (algorithm>0||(specialOptions_&4194304)!=0) {
     67      // Gub
     68      resolveGub((9*modelPtr_->numberRows())/10);
     69      return;
     70    }
    6271  }
    6372  bool deleteSolver;
     
    787796    return;
    788797  }
    789   if ((specialOptions_&2097152)!=0) {
    790     resolveGub((9*modelPtr_->numberRows())/10);
    791     return;
     798  if ((specialOptions_&2097152)!=0||(specialOptions_&4194304)!=0) {
     799    bool takeHint;
     800    OsiHintStrength strength;
     801    int algorithm = 0;
     802    getHintParam(OsiDoDualInResolve,takeHint,strength);
     803    if (strength!=OsiHintIgnore)
     804      algorithm = takeHint ? -1 : 1;
     805    if (algorithm>0||(specialOptions_&4194304)!=0) {
     806      // Gub
     807      resolveGub((9*modelPtr_->numberRows())/10);
     808      return;
     809    }
    792810  }
    793811  //void pclp(char *);
     
    12861304  OsiHintStrength strength;
    12871305  // Switch off printing if asked to
    1288   bool gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength));
    1289   assert (gotHint);
     1306  getHintParam(OsiDoReducePrint,takeHint,strength);
    12901307  int saveMessageLevel=modelPtr_->logLevel();
    12911308  if (strength!=OsiHintIgnore&&takeHint) {
     
    12961313      modelPtr_->messageHandler()->setLogLevel(0);
    12971314  }
     1315  //modelPtr_->messageHandler()->setLogLevel(1);
    12981316  setBasis(basis_,modelPtr_);
    12991317  // find gub
     
    13041322  ClpSimplex * model2 =
    13051323    static_cast<ClpSimplexOther *> (modelPtr_)->gubVersion(which,whichC,
    1306                                                            needed);
     1324                                                           needed,100);
    13071325  if (model2) {
    13081326    // move in solution
     
    13121330    ClpPrimalColumnSteepest steepest(5);
    13131331    model2->setPrimalColumnPivotAlgorithm(steepest);
    1314     double time1 = CoinCpuTime();
     1332    //double time1 = CoinCpuTime();
    13151333    model2->primal();
    13161334    //printf("Primal took %g seconds\n",CoinCpuTime()-time1);
     
    13211339    modelPtr_->primal(1);
    13221340    modelPtr_->setNumberIterations(totalIterations+modelPtr_->numberIterations());
     1341    delete model2;
    13231342  } else {
    13241343    modelPtr_->dual();
    13251344  }
     1345  delete [] which;
     1346  delete [] whichC;
    13261347  basis_ = getBasis(modelPtr_);
    13271348  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
     
    63336354      }
    63346355      // set normal
    6335       specialOptions_ &= (2047|3*8192|15*65536|2097152);
     6356      specialOptions_ &= (2047|3*8192|15*65536|2097152|4194304);
    63366357      if (otherInformation!=NULL) {
    63376358        int * array = static_cast<int *> (otherInformation);
Note: See TracChangeset for help on using the changeset viewer.