- Timestamp:
- Jan 25, 2011 5:58:03 AM (10 years ago)
- Location:
- trunk/Clp/src
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Clp/src/ClpDynamicMatrix.cpp
r1677 r1678 177 177 startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1); 178 178 if (!numberGubColumns_) { 179 180 179 //startColumn_ = new CoinBigIndex [1]; 180 startColumn_[0] = 0; 181 181 } 182 182 CoinBigIndex numberElements = startColumn_[numberGubColumns_]; … … 235 235 numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements(); 236 236 matrix_ = originalMatrix; 237 //delete originalMatrixA; 237 238 flags_ &= ~1; 238 239 // resize model (matrix stays same) 240 // modify frequency 241 if (frequency>=50) 242 frequency = 50+(frequency-50)/2; 239 243 int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1; 240 244 model->resize(newRowSize, numberNeeded); … … 903 907 #endif 904 908 #endif 905 return 0; 909 if (numberStaticRows_+numberActiveSets_<model->numberRows()) 910 return 0; 911 else 912 return 1; 906 913 } 907 914 /* … … 2275 2282 } else { 2276 2283 // solo key 2277 bool needKey ;2284 bool needKey=false; 2278 2285 if (numberActive) { 2279 2286 if (whichKey<maximumGubColumns_) { -
trunk/Clp/src/ClpPackedMatrix.cpp
r1665 r1678 4013 4013 // Say no gaps 4014 4014 flags_ &= ~2; 4015 if (type_>=10) 4016 return true; // gub 4015 4017 if (check == 14 || check == 10) { 4016 4018 if (matrix_->getNumElements() < columnStart[numberColumns]) { -
trunk/Clp/src/ClpSimplex.cpp
r1677 r1678 1828 1828 1829 1829 // Update hidden stuff e.g. effective RHS and gub 1830 matrix_->updatePivot(this, oldIn, oldOut);1830 int invertNow=matrix_->updatePivot(this, oldIn, oldOut); 1831 1831 objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_); 1832 1832 if (handler_->logLevel() > 7) { … … 1985 1985 } 1986 1986 return 1; 1987 } else if (factorization_->timeToRefactorize() && !dontInvert) { 1987 } else if ((factorization_->timeToRefactorize() && !dontInvert) 1988 ||invertNow) { 1988 1989 //printf("ret after %d pivots\n",factorization_->pivots()); 1989 1990 return 1; -
trunk/Clp/src/ClpSimplexOther.cpp
r1677 r1678 4453 4453 // Sets basis from original 4454 4454 void 4455 ClpSimplexOther::setGubBasis( constClpSimplex &original,const int * whichRows,4455 ClpSimplexOther::setGubBasis(ClpSimplex &original,const int * whichRows, 4456 4456 const int * whichColumns) 4457 4457 { … … 4470 4470 //assert (firstOdd==numberNormal); 4471 4471 double * solution = primalColumnSolution(); 4472 constdouble * originalSolution = original.primalColumnSolution();4472 double * originalSolution = original.primalColumnSolution(); 4473 4473 const double * upperSet = gubMatrix->upperSet(); 4474 // Column copy of GUB part 4474 4475 int numberSets = gubMatrix->numberSets(); 4475 4476 const int * startSet = gubMatrix->startSets(); 4476 const CoinBigIndex * startColumn= gubMatrix->startColumn();4477 const CoinBigIndex * columnStart = gubMatrix->startColumn(); 4477 4478 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 4478 4724 for (int i=0;i<numberSets;i++) { 4479 4725 for (int j=startSet[i];j<startSet[i+1];j++) { 4480 4726 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound); 4481 4727 int iColumn = whichColumns[j+numberNormal]; 4482 if (iColumn<numberColumns) 4728 if (iColumn<numberColumns) { 4483 4729 columnIsGub[iColumn] = whichRows[numberNonGub+i]; 4730 } 4484 4731 } 4485 4732 } 4733 #endif 4486 4734 int * numberKey = new int [numberRows]; 4487 4735 memset(numberKey,0,numberRows*sizeof(int)); … … 4530 4778 int iSet = iOrig - numberColumns; 4531 4779 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 ) { 4533 4785 assert(numberKey[iRow]); 4534 4786 if (numberKey[iRow]==1) … … 4553 4805 int chosen=-1; 4554 4806 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]; 4556 4808 int iOrig = whichColumns[j+numberNormal]; 4557 4809 double value; 4558 4810 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 4560 4816 if (value>upper) 4561 4817 gubMatrix->setStatus(i,ClpSimplex::atLowerBound); -
trunk/Clp/src/ClpSimplexOther.hpp
r1676 r1678 197 197 int factorizationFrequency=50); 198 198 /// Sets basis from original 199 void setGubBasis( constClpSimplex &original,const int * whichRows,199 void setGubBasis(ClpSimplex &original,const int * whichRows, 200 200 const int * whichColumns); 201 201 /// Restores basis to original -
trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp
r1677 r1678 57 57 } 58 58 #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 } 62 71 } 63 72 bool deleteSolver; … … 787 796 return; 788 797 } 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 } 792 810 } 793 811 //void pclp(char *); … … 1286 1304 OsiHintStrength strength; 1287 1305 // Switch off printing if asked to 1288 bool gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength)); 1289 assert (gotHint); 1306 getHintParam(OsiDoReducePrint,takeHint,strength); 1290 1307 int saveMessageLevel=modelPtr_->logLevel(); 1291 1308 if (strength!=OsiHintIgnore&&takeHint) { … … 1296 1313 modelPtr_->messageHandler()->setLogLevel(0); 1297 1314 } 1315 //modelPtr_->messageHandler()->setLogLevel(1); 1298 1316 setBasis(basis_,modelPtr_); 1299 1317 // find gub … … 1304 1322 ClpSimplex * model2 = 1305 1323 static_cast<ClpSimplexOther *> (modelPtr_)->gubVersion(which,whichC, 1306 needed );1324 needed,100); 1307 1325 if (model2) { 1308 1326 // move in solution … … 1312 1330 ClpPrimalColumnSteepest steepest(5); 1313 1331 model2->setPrimalColumnPivotAlgorithm(steepest); 1314 double time1 = CoinCpuTime();1332 //double time1 = CoinCpuTime(); 1315 1333 model2->primal(); 1316 1334 //printf("Primal took %g seconds\n",CoinCpuTime()-time1); … … 1321 1339 modelPtr_->primal(1); 1322 1340 modelPtr_->setNumberIterations(totalIterations+modelPtr_->numberIterations()); 1341 delete model2; 1323 1342 } else { 1324 1343 modelPtr_->dual(); 1325 1344 } 1345 delete [] which; 1346 delete [] whichC; 1326 1347 basis_ = getBasis(modelPtr_); 1327 1348 modelPtr_->messageHandler()->setLogLevel(saveMessageLevel); … … 6333 6354 } 6334 6355 // set normal 6335 specialOptions_ &= (2047|3*8192|15*65536|2097152 );6356 specialOptions_ &= (2047|3*8192|15*65536|2097152|4194304); 6336 6357 if (otherInformation!=NULL) { 6337 6358 int * array = static_cast<int *> (otherInformation);
Note: See TracChangeset
for help on using the changeset viewer.