Changeset 336
- Timestamp:
- Mar 22, 2004 11:30:03 AM (16 years ago)
- Location:
- trunk
- Files:
-
- 34 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/ClpDualRowSteepest.cpp
r266 r336 9 9 #include "CoinHelperFunctions.hpp" 10 10 #include <cstdio> 11 12 11 //############################################################################# 13 12 // Constructors / Destructor / Assignment -
trunk/ClpDynamicExampleMatrix.cpp
r335 r336 352 352 for (int i=0;i<numberIds;i++) { 353 353 int sequence = ids[i]; 354 int iSet = set[sequence]; 354 355 CoinBigIndex start = startColumnGen_[sequence]; 355 356 addColumn(startColumnGen_[sequence+1]-start, … … 359 360 columnLowerGen_ ? columnLowerGen_[sequence] : 0, 360 361 columnUpperGen_ ? columnUpperGen_[sequence] : 1.0e30, 361 set[sequence],getDynamicStatus(i));362 idGen_[i Set]=sequence; // say which one in362 iSet,getDynamicStatus(i)); 363 idGen_[i]=sequence; // say which one in 363 364 setDynamicStatusGen(sequence,inSmall); 364 365 } … … 427 428 // Partial pricing 428 429 void 429 ClpDynamicExampleMatrix::partialPricing(ClpSimplex * model, int start, int end,430 ClpDynamicExampleMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 430 431 int & bestSequence, int & numberWanted) 431 432 { 433 numberWanted=currentWanted_; 432 434 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_) { 440 436 // 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 { 445 439 // 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); 451 442 endG2 = min(endG2,numberSets_); 452 443 //printf("gub price - set start %d end %d\n", … … 460 451 int structuralOffset = slackOffset+numberSets_; 461 452 int structuralOffset2 = structuralOffset+maximumGubColumns_; 462 int saveWanted=numberWanted;463 453 // If nothing found yet can go all the way to end 464 454 int endAll = endG2; … … 479 469 // startG2,endG2,numberWanted); 480 470 int bestSet=-1; 471 int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; 472 int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_; 481 473 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) { 484 475 // give up 485 476 numberWanted=0; 486 477 break; 487 } else if (iSet==endG2 -1&&bestSequence>=0) {478 } else if (iSet==endG2&&bestSequence>=0) { 488 479 break; 489 480 } … … 618 609 savedBestSet_ = bestSet; 619 610 } 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; 620 622 // See if may be finished 621 623 if (!startG2&&bestSequence<0) 622 624 infeasibilityWeight_=model_->infeasibilityCost(); 623 else 625 else if (bestSequence>=0) 624 626 infeasibilityWeight_=-1.0; 627 currentWanted_=numberWanted; 625 628 } 626 629 } -
trunk/ClpDynamicMatrix.cpp
r335 r336 32 32 sumOfRelaxedPrimalInfeasibilities_(0.0), 33 33 savedBestGubDual_(0.0), 34 savedBestDj_(0.0),35 savedBestSequence_(-1),36 34 savedBestSet_(0), 37 35 backToPivotRow_(NULL), … … 103 101 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; 104 102 savedBestGubDual_=rhs.savedBestGubDual_; 105 savedBestDj_ = rhs.savedBestDj_;106 savedBestSequence_=rhs.savedBestSequence_;107 103 savedBestSet_=rhs.savedBestSet_; 108 104 noCheck_ = rhs.noCheck_; … … 163 159 numberStaticRows_ = numberRows; 164 160 savedBestGubDual_=0.0; 165 savedBestSequence_=-1;166 savedBestDj_ = 0.0;167 161 savedBestSet_=0; 168 162 // Number of columns needed … … 365 359 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; 366 360 savedBestGubDual_=rhs.savedBestGubDual_; 367 savedBestDj_ = rhs.savedBestDj_;368 savedBestSequence_=rhs.savedBestSequence_;369 361 savedBestSet_=rhs.savedBestSet_; 370 362 noCheck_ = rhs.noCheck_; … … 396 388 // Partial pricing 397 389 void 398 ClpDynamicMatrix::partialPricing(ClpSimplex * model, int start, int end,390 ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 399 391 int & bestSequence, int & numberWanted) 400 392 { 393 numberWanted=currentWanted_; 401 394 assert(!model->rowScale()); 402 395 if (numberSets_) { … … 404 397 // always??? 405 398 //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); 408 400 } else { 409 401 // no gub 410 ClpPackedMatrix::partialPricing(model,start ,end,bestSequence,numberWanted);402 ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted); 411 403 return; 412 404 } 413 405 if (numberWanted>0) { 414 406 // 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); 420 409 endG2 = min(endG2,numberSets_); 421 410 //printf("gub price - set start %d end %d\n", … … 428 417 int slackOffset = lastDynamic_+numberRows; 429 418 int structuralOffset = slackOffset+numberSets_; 430 int saveWanted=numberWanted;431 419 // If nothing found yet can go all the way to end 432 420 int endAll = endG2; … … 458 446 } 459 447 #endif 448 int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; 449 int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_; 460 450 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) { 463 452 // give up 464 453 numberWanted=0; 465 454 break; 466 } else if (iSet==endG2 -1&&bestSequence>=0) {455 } else if (iSet==endG2&&bestSequence>=0) { 467 456 break; 468 457 } … … 568 557 if (!startG2&&bestSequence<0) 569 558 infeasibilityWeight_=model_->infeasibilityCost(); 570 else 559 else if (bestSequence>=0) 571 560 infeasibilityWeight_=-1.0; 572 561 } 562 currentWanted_=numberWanted; 573 563 } 574 564 /* Returns effective RHS if it is being used. This is used for long problems … … 1711 1701 modifyOffset(key,valueOfKey); 1712 1702 rhsOffset_[newRow] = -shift; // sign? 1703 #ifdef CLP_DEBUG 1704 model->rowArray(1)->checkClear(); 1705 #endif 1713 1706 // 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]; 1717 1710 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(); 1721 1716 if (updateStatus) { 1722 1717 if (updateStatus==3) { -
trunk/ClpFactorization.cpp
r328 r336 171 171 numberColumnBasic, 172 172 NULL,NULL,NULL); 173 173 // and recompute as network side say different 174 if (model->numberIterations()) 175 numberRowBasic = numberBasic - numberColumnBasic; 174 176 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, 176 180 2 * numberElements ); 177 181 //fill … … 214 218 // If we get here status is 0 or -1 215 219 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); 216 225 int * permuteBack = permuteBack_; 217 226 int * back = pivotColumnBack_; 218 227 //int * pivotTemp = pivotColumn_; 219 //ClpDisjointCopyN ( pivotVariable, numberRows _, pivotTemp );228 //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp ); 220 229 // Redo pivot order 221 230 for (i=0;i<numberRowBasic;i++) { … … 224 233 pivotVariable[permuteBack[back[i]]]=k+numberColumns; 225 234 } 226 for (;i< numberRows;i++) {235 for (;i<useNumberRows;i++) { 227 236 int k = pivotTemp[i]; 228 237 // so rowIsBasic[k] would be permuteBack[back[i]] … … 232 241 // these arrays start off as copies of permute 233 242 // (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_ ); 236 245 if (networkMatrix) { 237 246 maximumPivots(saveMaximumPivots); … … 249 258 if (!doCheck) { 250 259 gutsOfDestructor(); 260 status_=0; 251 261 #if 0 252 262 // but put back permute arrays so odd things will work … … 460 470 partial update already in U */ 461 471 int 462 ClpFactorization::replaceColumn ( CoinIndexedVector * regionSparse, 463 int pivotRow, 464 double pivotCheck , 465 bool checkBeforeModifying) 472 ClpFactorization::replaceColumn ( const ClpSimplex * model, 473 CoinIndexedVector * regionSparse, 474 CoinIndexedVector * tableauColumn, 475 int pivotRow, 476 double pivotCheck , 477 bool checkBeforeModifying) 466 478 { 467 479 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 472 490 } else { 473 491 if (doCheck) { -
trunk/ClpGubDynamicMatrix.cpp
r335 r336 16 16 #include "ClpGubDynamicMatrix.hpp" 17 17 #include "ClpMessage.hpp" 18 #define CLP_DEBUG18 //#define CLP_DEBUG 19 19 //#define CLP_DEBUG_PRINT 20 20 //############################################################################# … … 39 39 lowerSet_(NULL), 40 40 upperSet_(NULL), 41 model_(NULL),42 41 numberGubColumns_(0), 43 42 firstAvailable_(0), … … 75 74 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_); 76 75 upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_); 77 model_ = rhs.model_;78 76 } 79 77 … … 188 186 backward_ = new int[numberNeeded]; 189 187 backToPivotRow_ = new int[numberNeeded]; 188 // We know a bit better 189 delete [] changeCost_; 190 190 changeCost_ = new double [numberRows+numberSets_]; 191 191 keyVariable_ = new int[numberSets_]; … … 272 272 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_); 273 273 upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_); 274 model_ = rhs.model_;275 274 } 276 275 return *this; … … 285 284 // Partial pricing 286 285 void 287 ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, int start, int end,286 ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 288 287 int & bestSequence, int & numberWanted) 289 288 { 290 289 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; 298 295 } else { 299 // no gub300 ClpPackedMatrix::partialPricing(model,start,end,bestSequence,numberWanted);301 return;302 }303 if (numberWanted>0) {304 296 // 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); 310 299 endG2 = min(endG2,numberSets_); 311 300 //printf("gub price - set start %d end %d\n", … … 318 307 int numberRows = model->numberRows(); 319 308 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_; 321 313 if (bestSequence>=0) 322 314 bestDj = fabs(reducedCost[bestSequence]); … … 346 338 } 347 339 #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) { 350 344 // give up 351 345 numberWanted=0; 346 break; 347 } else if (iSet==endG2&&bestSequence>=0) { 352 348 break; 353 349 } … … 442 438 break; 443 439 } 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; 444 455 } 445 456 if (bestSequence!=saveSequence||bestType>=0) { … … 531 542 upperColumn[bestSequence],0.0); 532 543 } 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; 535 554 } 536 555 // This is local to Gub to allow synchronization when status is good … … 1040 1059 for (iSet=0;iSet<numberSets_;iSet++) 1041 1060 toIndex_[iSet]=-1; 1042 fromIndex_ = new int [ max(getNumRows()+1,numberSets_)];1061 fromIndex_ = new int [numberRows+numberSets_]; 1043 1062 double tolerance = model->primalTolerance(); 1044 1063 double * element = matrix_->getMutableElements(); … … 1791 1810 bool print=false; 1792 1811 int iSet; 1812 int trueIn=-1; 1813 int trueOut=-1; 1814 int numberRows = model->numberRows(); 1815 int numberColumns = model->numberColumns(); 1793 1816 if (sequenceIn==firstAvailable_) { 1794 1817 if (doPrinting) … … 1811 1834 if (iSet>=0) { 1812 1835 int bigSequence = id_[sequenceIn-firstDynamic_]; 1836 trueIn=bigSequence+numberRows+numberColumns+numberSets_; 1813 1837 if (doPrinting) 1814 1838 printf(" incoming set %d big seq %d",iSet,bigSequence); 1815 1839 print=true; 1816 1840 } 1841 } else if (sequenceIn>=numberRows+numberColumns) { 1842 trueIn = numberRows+numberColumns+gubSlackIn_; 1817 1843 } 1818 1844 if (sequenceOut<lastDynamic_) { … … 1820 1846 if (iSet>=0) { 1821 1847 int bigSequence = id_[sequenceOut-firstDynamic_]; 1848 trueOut=bigSequence+firstDynamic_; 1822 1849 if (model->getStatus(sequenceOut)==ClpSimplex::atUpperBound) 1823 1850 setDynamicStatus(bigSequence,atUpperBound); … … 1835 1862 printf("\n"); 1836 1863 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; 1837 1869 if (doPrinting&&0) { 1838 1870 for (int i=0;i<numberSets_;i++) { … … 1927 1959 } 1928 1960 } 1961 /* Just for debug - may be extended to other matrix types later. 1962 Returns number of primal infeasibilities. 1963 */ 1964 int 1965 ClpGubDynamicMatrix::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 31 31 sumOfRelaxedDualInfeasibilities_(0.0), 32 32 sumOfRelaxedPrimalInfeasibilities_(0.0), 33 infeasibilityWeight_(0.0), 33 34 start_(NULL), 34 35 end_(NULL), … … 45 46 toIndex_(NULL), 46 47 fromIndex_(NULL), 48 model_(NULL), 47 49 numberDualInfeasibilities_(0), 48 50 numberPrimalInfeasibilities_(0), 51 noCheck_(-1), 49 52 numberSets_(0), 50 53 saveNumber_(0), … … 79 82 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,numberColumns); 80 83 changeCost_ = ClpCopyOfArray(rhs.changeCost_,getNumRows()+numberSets_); 84 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+numberSets_+1); 81 85 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_); 82 86 // find longest set … … 94 98 next_ = ClpCopyOfArray(rhs.next_,numberColumns+numberSets_+2*length); 95 99 toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_); 96 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1);97 100 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; 98 101 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; 99 102 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; 100 103 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; 104 infeasibilityWeight_=rhs.infeasibilityWeight_; 101 105 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; 102 106 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; 107 noCheck_ = rhs.noCheck_; 103 108 firstGub_ = rhs.firstGub_; 104 109 lastGub_ = rhs.lastGub_; 105 110 gubType_ = rhs.gubType_; 111 model_ = rhs.model_; 106 112 } 107 113 … … 115 121 sumOfRelaxedDualInfeasibilities_(0.0), 116 122 sumOfRelaxedPrimalInfeasibilities_(0.0), 123 infeasibilityWeight_(0.0), 117 124 start_(NULL), 118 125 end_(NULL), … … 129 136 toIndex_(NULL), 130 137 fromIndex_(NULL), 138 model_(NULL), 131 139 numberDualInfeasibilities_(0), 132 140 numberPrimalInfeasibilities_(0), 141 noCheck_(-1), 133 142 numberSets_(0), 134 143 saveNumber_(0), … … 159 168 gubSlackIn_(-1) 160 169 { 170 model_=NULL; 161 171 numberSets_ = numberSets; 162 172 start_ = ClpCopyOfArray(start,numberSets_); … … 230 240 savedKeyVariable_= new int [numberSets_]; 231 241 memset(savedKeyVariable_,0,numberSets_*sizeof(int)); 242 noCheck_ = -1; 243 infeasibilityWeight_=0.0; 232 244 } 233 245 … … 238 250 sumOfRelaxedDualInfeasibilities_(0.0), 239 251 sumOfRelaxedPrimalInfeasibilities_(0.0), 252 infeasibilityWeight_(0.0), 240 253 start_(NULL), 241 254 end_(NULL), … … 252 265 toIndex_(NULL), 253 266 fromIndex_(NULL), 267 model_(NULL), 254 268 numberDualInfeasibilities_(0), 255 269 numberPrimalInfeasibilities_(0), 270 noCheck_(-1), 256 271 numberSets_(0), 257 272 saveNumber_(0), … … 324 339 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,numberColumns); 325 340 changeCost_ = ClpCopyOfArray(rhs.changeCost_,getNumRows()+numberSets_); 341 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+numberSets_+1); 326 342 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_); 327 343 // find longest set … … 339 355 next_ = ClpCopyOfArray(rhs.next_,numberColumns+numberSets_+2*length); 340 356 toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_); 341 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1);342 357 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; 343 358 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; 344 359 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; 345 360 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; 361 infeasibilityWeight_=rhs.infeasibilityWeight_; 346 362 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; 347 363 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; 364 noCheck_ = rhs.noCheck_; 348 365 firstGub_ = rhs.firstGub_; 349 366 lastGub_ = rhs.lastGub_; 350 367 gubType_ = rhs.gubType_; 368 model_ = rhs.model_; 351 369 } 352 370 return *this; … … 530 548 // reduce for gub 531 549 factor *= 0.5; 550 assert (!y->getNumElements()); 532 551 if (numberInRowArray>factor*numberRows||!rowCopy) { 533 552 // do by column … … 542 561 int iSet = -1; 543 562 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]; 575 590 } 576 591 } else { 577 592 djMod = 0.0; 578 593 } 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; 581 667 CoinBigIndex j; 582 668 for (j=columnStart[iColumn]; … … 586 672 } 587 673 if (fabs(value)>zeroTolerance) { 588 array[numberNonZero]=value;589 674 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; 590 690 } 591 691 } 592 692 } else { 593 // scaled594 // modify pi so can collapse to one loop595 for (i=0;i<numberInRowArray;i++) {596 int iRow = whichRow[i];597 pi[iRow]=scalar*piOld[i]*rowScale[iRow];598 }599 693 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; 623 713 CoinBigIndex j; 624 714 const double * columnScale = model->columnScale(); … … 626 716 j<columnStart[iColumn]+columnLength[iColumn];j++) { 627 717 int iRow = row[j]; 628 value += pi[iRow]*elementByColumn[j] ;718 value += pi[iRow]*elementByColumn[j]*rowScale[iRow]; 629 719 } 630 720 value *= columnScale[iColumn]; 631 721 if (fabs(value)>zeroTolerance) { 632 array[numberNonZero]=value;633 722 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; 688 740 } 689 741 } 690 742 } 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 } 782 757 } 783 758 } … … 1378 1353 // Partial pricing 1379 1354 void 1380 ClpGubMatrix::partialPricing(ClpSimplex * model, int start, int end,1355 ClpGubMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 1381 1356 int & bestSequence, int & numberWanted) 1382 1357 { 1358 numberWanted=currentWanted_; 1383 1359 if (numberSets_) { 1384 1360 // 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) { 1387 1366 // do gub 1388 1367 const double * element =matrix_->getElements(); … … 1402 1381 int numberRows = model->numberRows(); 1403 1382 if (bestSequence>=0) 1404 bestDj = fabs( reducedCost[bestSequence]);1383 bestDj = fabs(this->reducedCost(model,bestSequence)); 1405 1384 else 1406 1385 bestDj=tolerance; 1407 1386 int sequenceOut = model->sequenceOut(); 1408 1387 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; 1411 1398 int iSet = -1; 1412 1399 double djMod=0.0; … … 1415 1402 double bestDjMod=0.0; 1416 1403 // 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 } 1418 1412 if (backward_[iSequence]!=iSet) { 1419 1413 // get pi on gub row 1420 1414 iSet = backward_[iSequence]; 1421 1415 if (iSet>=0) { 1416 nSets++; 1422 1417 int iBasic = keyVariable_[iSet]; 1423 1418 if (iBasic>=numberColumns) { … … 1592 1587 model->costRegion()[bestSequence] = 0.0; 1593 1588 } 1589 savedBestSequence_ = bestSequence; 1590 savedBestDj_ = reducedCost[savedBestSequence_]; 1594 1591 } 1595 1592 } else { … … 1598 1595 // startG,endG,numberWanted); 1599 1596 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 } 1600 1604 if (backward_[iSequence]!=iSet) { 1601 1605 // get pi on gub row 1602 1606 iSet = backward_[iSequence]; 1603 1607 if (iSet>=0) { 1608 nSets++; 1604 1609 int iBasic = keyVariable_[iSet]; 1605 1610 if (iBasic>=numberColumns) { … … 1772 1777 } 1773 1778 } 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; 1774 1784 } 1775 1785 if (numberWanted) { 1776 1786 // 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); 1778 1792 } 1779 1793 } else { 1780 1794 // 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; 1783 1800 } 1784 1801 /* expands an updated column to allow for extra rows which the main … … 2046 2063 if (number>saveNumber_) { 2047 2064 // clear 2065 double theta = model->theta(); 2066 double * solution = model->solutionRegion(); 2048 2067 for (i=saveNumber_;i<number;i++) { 2049 #ifdef CLP_DEBUG_PRINT2050 2068 int iRow = index[i]; 2051 2069 int iColumn = pivotVariable[iRow]; 2070 #ifdef CLP_DEBUG_PRINT 2052 2071 printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n", 2053 2072 iColumn,fromIndex_[i-saveNumber_],lower[iColumn],upper[iColumn],array[i], 2054 2073 solution[iColumn],solution[iColumn]-model->theta()*array[i],model->theta()); 2055 2074 #endif 2075 double value = array[i]; 2056 2076 array[i]=0.0; 2057 2077 int iSet=fromIndex_[i-saveNumber_]; 2058 2078 toIndex_[iSet]=-1; 2079 if (iSet==iSetIn&&iColumn<numberColumns) { 2080 // update as may need value 2081 solution[iColumn] -= theta*value; 2082 } 2059 2083 } 2060 2084 } … … 2295 2319 backToPivotRow_[iPivot]=i; 2296 2320 } 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 } 2297 2330 double * dj = model->djRegion(); 2298 2331 double * dual = model->dualRowSolution(); … … 2409 2442 // and get statistics for column generation 2410 2443 synchronize(model,4); 2444 infeasibilityWeight_=-1.0; 2411 2445 } 2412 2446 break; … … 2605 2639 } 2606 2640 number = numberBasic; 2641 if (model->numberIterations()) 2642 assert (number==model->numberRows()); 2607 2643 } 2608 2644 break; … … 2720 2756 //if (!alpha) 2721 2757 //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), 2723 2761 iRow, alpha); 2724 2762 returnCode = max (updateStatus, returnCode); … … 2741 2779 //if (!alpha) 2742 2780 //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), 2744 2784 pivotRow, alpha); 2745 2785 returnCode = max (updateStatus, returnCode); … … 2754 2794 printf("TTTTTTry 4 %d\n",possiblePivotKey_); 2755 2795 #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), 2757 2799 possiblePivotKey_, 2758 bestAlpha);2800 bestAlpha); 2759 2801 returnCode = max (updateStatus, returnCode); 2760 2802 incomingColumn = pivotVariable[possiblePivotKey_]; … … 2801 2843 //if (!alpha) 2802 2844 //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), 2804 2848 iRow, alpha); 2805 2849 returnCode = max (updateStatus, returnCode); … … 2848 2892 //if (!alpha) 2849 2893 //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), 2851 2897 iRow, alpha); 2852 2898 returnCode = max (updateStatus, returnCode); … … 2868 2914 //if (!alpha) 2869 2915 //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), 2871 2919 pivotRow, alpha); 2872 2920 returnCode = max (updateStatus, returnCode); … … 2877 2925 } else { 2878 2926 // 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), 2880 2930 model->pivotRow(), 2881 2931 model->alpha()); … … 2999 3049 returnCode=synchronize(model,8); 3000 3050 } 3051 break; 3052 default: 3001 3053 break; 3002 3054 } … … 3603 3655 int pivotRow = model->pivotRow(); 3604 3656 int iSetIn; 3605 if (sequenceIn<numberColumns) 3657 // Correct sequence in 3658 trueSequenceIn_=sequenceIn; 3659 if (sequenceIn<numberColumns) { 3606 3660 iSetIn = backward_[sequenceIn]; 3607 else if (sequenceIn<numberColumns+numberRows)3661 } else if (sequenceIn<numberColumns+numberRows) { 3608 3662 iSetIn = -1; 3609 else3663 } else { 3610 3664 iSetIn = gubSlackIn_; 3665 trueSequenceIn_ = numberColumns+numberRows+iSetIn; 3666 } 3611 3667 int iSetOut=-1; 3668 trueSequenceOut_=sequenceOut; 3612 3669 if (sequenceOut<numberColumns) { 3613 3670 iSetOut = backward_[sequenceOut]; … … 3620 3677 else 3621 3678 assert(iSetOut == fromIndex_[iExtra]); 3679 trueSequenceOut_ = numberColumns+numberRows+iSetOut; 3622 3680 } 3623 3681 if (rhsOffset_) { … … 4017 4075 return 0; 4018 4076 } 4019 4020 4077 // Switches off dj checking each factorization (for BIG models) 4078 void 4079 ClpGubMatrix::switchOffCheck() 4080 { 4081 noCheck_=0; 4082 infeasibilityWeight_=0.0; 4083 } 4084 // Correct sequence in and out to give true value 4085 void 4086 ClpGubMatrix::correctSequence(int & sequenceIn, int & sequenceOut) const 4087 { 4088 sequenceIn = trueSequenceIn_; 4089 sequenceOut = trueSequenceOut_; 4090 } -
trunk/ClpMatrixBase.cpp
r328 r336 20 20 ClpMatrixBase::ClpMatrixBase () : 21 21 rhsOffset_(NULL), 22 startFraction_(0.0), 23 endFraction_(1.0), 24 savedBestDj_(0.0), 25 originalWanted_(0), 26 currentWanted_(0), 27 savedBestSequence_(-1), 22 28 type_(-1), 23 29 lastRefresh_(-1), 24 30 refreshFrequency_(0), 31 minimumObjectsScan_(-1), 32 minimumGoodReducedCosts_(-1), 33 trueSequenceIn_(-1), 34 trueSequenceOut_(-1), 25 35 skipDualCheck_(false) 26 36 { … … 35 45 skipDualCheck_(rhs.skipDualCheck_) 36 46 { 47 startFraction_ = rhs.startFraction_; 48 endFraction_ = rhs.endFraction_; 49 savedBestDj_ = rhs.savedBestDj_; 50 originalWanted_ = rhs.originalWanted_; 51 currentWanted_ = rhs.currentWanted_; 52 savedBestSequence_ = rhs.savedBestSequence_; 37 53 lastRefresh_ = rhs.lastRefresh_; 38 54 refreshFrequency_ = rhs.refreshFrequency_; 55 minimumObjectsScan_ = rhs.minimumObjectsScan_; 56 minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_; 57 trueSequenceIn_ = rhs.trueSequenceIn_; 58 trueSequenceOut_ = rhs.trueSequenceOut_; 39 59 skipDualCheck_ = rhs.skipDualCheck_; 40 60 int numberRows = rhs.getNumRows(); … … 69 89 rhsOffset_=NULL; 70 90 } 91 startFraction_ = rhs.startFraction_; 92 endFraction_ = rhs.endFraction_; 93 savedBestDj_ = rhs.savedBestDj_; 94 originalWanted_ = rhs.originalWanted_; 95 currentWanted_ = rhs.currentWanted_; 96 savedBestSequence_ = rhs.savedBestSequence_; 71 97 lastRefresh_ = rhs.lastRefresh_; 72 98 refreshFrequency_ = rhs.refreshFrequency_; 99 minimumObjectsScan_ = rhs.minimumObjectsScan_; 100 minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_; 101 trueSequenceIn_ = rhs.trueSequenceIn_; 102 trueSequenceOut_ = rhs.trueSequenceOut_; 73 103 skipDualCheck_ = rhs.skipDualCheck_; 74 104 } … … 167 197 // Partial pricing 168 198 void 169 ClpMatrixBase::partialPricing(ClpSimplex * model, int start, intend,199 ClpMatrixBase::partialPricing(ClpSimplex * model, double start, double end, 170 200 int & bestSequence, int & numberWanted) 171 201 { … … 355 385 ClpMatrixBase::reducedCost(ClpSimplex * model,int sequence) const 356 386 { 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 */ 397 int 398 ClpMatrixBase::checkFeasible() const 399 { 400 return 0; 401 } 402 // Correct sequence in and out to give true value 403 void 404 ClpMatrixBase::correctSequence(int & sequenceIn, int & sequenceOut) const 405 { 406 } 407 408 -
trunk/ClpMessage.cpp
r325 r336 94 94 {CLP_BARRIER_FEASIBLE,57,2,"Infeasibilities - bound %g , primal %g ,dual %g"}, 95 95 {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"}, 96 97 {CLP_DUMMY_END,999999,0,""} 97 98 }; -
trunk/ClpNetworkMatrix.cpp
r309 r336 772 772 // Partial pricing 773 773 void 774 ClpNetworkMatrix::partialPricing(ClpSimplex * model, int start, int end,774 ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 775 775 int & bestSequence, int & numberWanted) 776 776 { 777 numberWanted=currentWanted_; 777 778 int j; 779 int start = (int) (startFraction*numberColumns_); 780 int end = min((int) (endFraction*numberColumns_+1),numberColumns_); 778 781 double tolerance=model->currentDualTolerance(); 779 782 double * reducedCost = model->djRegion(); … … 895 898 value -= duals[iRowP]; 896 899 reducedCost[bestSequence] = value; 900 savedBestSequence_ = bestSequence; 901 savedBestDj_ = reducedCost[savedBestSequence_]; 897 902 } 898 903 } else { … … 992 997 value -= duals[iRowP]; 993 998 reducedCost[bestSequence] = value; 994 } 995 } 999 savedBestSequence_ = bestSequence; 1000 savedBestDj_ = reducedCost[savedBestSequence_]; 1001 } 1002 } 1003 currentWanted_=numberWanted; 996 1004 } 997 1005 // Allow any parts of a created CoinMatrix to be deleted -
trunk/ClpNonLinearCost.cpp
r328 r336 480 480 break; 481 481 case ClpSimplex::isFree: 482 if (toNearest)483 482 //if (toNearest) 483 //solution[iSequence] = 0.0; 484 484 break; 485 485 case ClpSimplex::atUpperBound: … … 577 577 //assert (iRange==whichRange_[iSequence]); 578 578 } 579 //feasibleCost_ /= (model_->rhsScale()*model_->objScale()); 579 580 } 580 581 /* Goes through one bound for each variable. … … 984 985 double value; 985 986 model_->getDblParam(ClpObjOffset,value); 986 return feasibleCost_*model_->optimizationDirection()-value; 987 return feasibleCost_*model_->optimizationDirection()/ 988 (model_->objectiveScale()*model_->rhsScale())-value; 987 989 } 988 990 // Get rid of real costs (just for moment) -
trunk/ClpPackedMatrix.cpp
r328 r336 26 26 : ClpMatrixBase(), 27 27 matrix_(NULL), 28 zeroElements_(false) 28 numberActiveColumns_(0), 29 zeroElements_(false), 30 hasGaps_(true) 29 31 { 30 32 setType(1); … … 38 40 { 39 41 matrix_ = new CoinPackedMatrix(*(rhs.matrix_)); 42 numberActiveColumns_ = rhs.numberActiveColumns_; 40 43 zeroElements_ = rhs.zeroElements_; 44 hasGaps_ = rhs.hasGaps_; 41 45 int numberRows = getNumRows(); 42 46 if (rhs.rhsOffset_&&numberRows) { … … 56 60 matrix_ = rhs; 57 61 zeroElements_ = false; 62 hasGaps_ = true; 63 numberActiveColumns_ = matrix_->getNumCols(); 58 64 setType(1); 59 65 … … 64 70 { 65 71 matrix_ = new CoinPackedMatrix(rhs); 72 numberActiveColumns_ = matrix_->getNumCols(); 66 73 zeroElements_ = false; 74 hasGaps_ = true; 67 75 setType(1); 68 76 … … 87 95 delete matrix_; 88 96 matrix_ = new CoinPackedMatrix(*(rhs.matrix_)); 97 numberActiveColumns_ = rhs.numberActiveColumns_; 89 98 zeroElements_ = rhs.zeroElements_; 99 hasGaps_ = rhs.hasGaps_; 90 100 } 91 101 return *this; … … 118 128 matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows, 119 129 numberColumns,whichColumns); 130 numberActiveColumns_ = matrix_->getNumCols(); 120 131 zeroElements_ = rhs.zeroElements_; 132 hasGaps_ = rhs.hasGaps_; 121 133 } 122 134 ClpPackedMatrix::ClpPackedMatrix ( … … 128 140 matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows, 129 141 numberColumns,whichColumns); 142 numberActiveColumns_ = matrix_->getNumCols(); 130 143 zeroElements_ = false; 144 hasGaps_=true; 131 145 setType(1); 132 146 } … … 140 154 copy->matrix_->reverseOrderedCopyOf(*matrix_); 141 155 copy->matrix_->removeGaps(); 156 copy->numberActiveColumns_ = copy->matrix_->getNumCols(); 157 copy->hasGaps_=false; 142 158 return copy; 143 159 } … … 153 169 const int * columnLength = matrix_->getVectorLengths(); 154 170 const double * elementByColumn = matrix_->getElements(); 155 int numberColumns = matrix_->getNumCols();156 171 //memset(y,0,matrix_->getNumRows()*sizeof(double)); 157 for (iColumn=0;iColumn<number Columns;iColumn++) {172 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 158 173 CoinBigIndex j; 159 174 double value = scalar*x[iColumn]; … … 177 192 const int * columnLength = matrix_->getVectorLengths(); 178 193 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++) { 182 196 CoinBigIndex j; 183 197 double value=0.0; … … 203 217 const int * columnLength = matrix_->getVectorLengths(); 204 218 const double * elementByColumn = matrix_->getElements(); 205 int numberColumns = matrix_->getNumCols();206 219 //memset(y,0,matrix_->getNumRows()*sizeof(double)); 207 for (iColumn=0;iColumn<number Columns;iColumn++) {220 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 208 221 CoinBigIndex j; 209 222 double value = x[iColumn]; … … 239 252 const int * columnLength = matrix_->getVectorLengths(); 240 253 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++) { 244 256 CoinBigIndex j; 245 257 double value=0.0; … … 278 290 double factor = 0.3; 279 291 // We may not want to do by row if there may be cache problems 280 int numberColumns = model->numberColumns();281 292 // It would be nice to find L2 cache size - for moment 512K 282 293 // Be slightly optimistic 283 if (number Columns*sizeof(double)>1000000) {284 if (numberRows*10<number Columns)294 if (numberActiveColumns_*sizeof(double)>1000000) { 295 if (numberRows*10<numberActiveColumns_) 285 296 factor=0.1; 286 else if (numberRows*4<number Columns)297 else if (numberRows*4<numberActiveColumns_) 287 298 factor=0.15; 288 else if (numberRows*2<number Columns)299 else if (numberRows*2<numberActiveColumns_) 289 300 factor=0.2; 290 301 //if (model->numberIterations()%50==0) 291 302 //printf("%d nonzero\n",numberInRowArray); 292 303 } 304 assert (!y->getNumElements()); 293 305 if (numberInRowArray>factor*numberRows||!rowCopy) { 294 306 // 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 } 295 313 int iColumn; 296 314 // get matrix data pointers … … 300 318 const double * elementByColumn = matrix_->getElements(); 301 319 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++) { 318 378 double value = 0.0; 319 379 CoinBigIndex j; … … 324 384 } 325 385 if (fabs(value)>zeroTolerance) { 326 array[numberNonZero]=value;327 386 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; 328 402 } 329 403 } 330 404 } 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++) { 338 424 double value = 0.0; 339 425 CoinBigIndex j; … … 342 428 j<columnStart[iColumn]+columnLength[iColumn];j++) { 343 429 int iRow = row[j]; 344 value += pi[iRow]*elementByColumn[j] ;430 value += pi[iRow]*elementByColumn[j]*rowScale[iRow]; 345 431 } 346 432 value *= columnScale[iColumn]; 347 433 if (fabs(value)>zeroTolerance) { 348 array[numberNonZero]=value;349 434 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; 402 452 } 403 453 } 404 454 } 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 } 494 469 } 495 470 } … … 521 496 } 522 497 } 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 */ 502 void 503 ClpPackedMatrix::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 */ 788 void 789 ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar, 790 const CoinIndexedVector * rowArray, 791 CoinIndexedVector * columnArray) const 792 { 793 } 523 794 /* Return <code>x * A + y</code> in <code>z</code>. 524 795 Squashes small elements and knows about ClpSimplex */ … … 542 813 const int * whichRow = rowArray->getIndices(); 543 814 bool packed = rowArray->packedMode(); 544 if (numberInRowArray>2 ||y->getNumElements()) {815 if (numberInRowArray>2) { 545 816 // do by rows 546 817 // ** Row copy is already scaled 547 818 int iRow; 548 int * mark = y->getIndices();549 int numberOriginal=y->getNumElements();550 819 int i; 820 int numberOriginal = 0; 551 821 if (packed) { 552 assert(!numberOriginal);553 822 numberNonZero=0; 554 823 // and set up mark as char array … … 556 825 double * array2 = y->denseVector(); 557 826 #ifdef CLP_DEBUG 558 int numberColumns = model->numberColumns(); 559 for (i=0;i<numberColumns;i++) { 827 for (i=0;i<numberActiveColumns_;i++) { 560 828 assert(!marked[i]); 561 829 assert(!array2[i]); … … 591 859 } 592 860 } 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; 601 863 // and set up mark as char array 602 864 char * marked = (char *) markVector; … … 788 1050 assert (!rowArray->packedMode()); 789 1051 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]; 793 1056 double value = 0.0; 794 1057 CoinBigIndex j; 795 1058 for (j=columnStart[iColumn]; 796 j<columnStart[iColumn ]+columnLength[iColumn];j++) {1059 j<columnStart[iColumn+1];j++) { 797 1060 int iRow = row[j]; 798 1061 value += pi[iRow]*elementByColumn[j]; 799 1062 } 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 } 800 1074 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]; 806 1079 double value = 0.0; 1080 double scale=columnScale[iColumn]; 807 1081 CoinBigIndex j; 808 const double * columnScale = model->columnScale();809 1082 for (j=columnStart[iColumn]; 810 j<columnStart[iColumn ]+columnLength[iColumn];j++) {1083 j<columnStart[iColumn+1];j++) { 811 1084 int iRow = row[j]; 812 1085 value += pi[iRow]*elementByColumn[j]*rowScale[iRow]; 813 1086 } 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; 815 1101 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 } 816 1132 } 817 1133 } … … 921 1237 { 922 1238 int numberRows = model->numberRows(); 923 int numberColumns = m odel->numberColumns();1239 int numberColumns = matrix_->getNumCols(); 924 1240 // If empty - return as sanityCheck will trap 925 1241 if (!numberRows||!numberColumns) … … 1170 1486 } 1171 1487 } 1172 #define RANDOMIZE1488 //#define RANDOMIZE 1173 1489 #ifdef RANDOMIZE 1174 1490 // randomize by up to 10% … … 1181 1497 if (overallSmallest<1.0e-1) 1182 1498 overallLargest = 1.0/sqrt(overallSmallest); 1183 overallLargest = min(100 0.0,overallLargest);1499 overallLargest = min(100.0,overallLargest); 1184 1500 overallSmallest=1.0e50; 1185 1501 for (iColumn=0;iColumn<numberColumns;iColumn++) { … … 1445 1761 const int * columnLength = matrix_->getVectorLengths(); 1446 1762 const double * elementByColumn = matrix_->getElements(); 1763 int numberRows = matrix_->getNumRows(); 1447 1764 int numberColumns = matrix_->getNumCols(); 1448 int numberRows = matrix_->getNumRows();1449 1765 int * mark = new int [numberRows]; 1450 1766 int i; 1767 // Say no gaps 1768 hasGaps_=false; 1451 1769 for (i=0;i<numberRows;i++) 1452 1770 mark[i]=-1; 1453 1771 for (iColumn=0;iColumn<numberColumns;iColumn++) { 1454 1772 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++) { 1457 1778 double value = fabs(elementByColumn[j]); 1458 1779 int iRow = row[j]; … … 1527 1848 { 1528 1849 int numberRows = model->numberRows(); 1529 int numberColumns = model->numberColumns();1850 int numberColumns = matrix_->getNumCols(); 1530 1851 int number = numberRows+numberColumns; 1531 1852 CoinBigIndex * weights = new CoinBigIndex[number]; … … 1560 1881 smallestPositive=COIN_DBL_MAX; 1561 1882 largestPositive=0.0; 1562 int numberColumns = matrix_->getNumCols();1563 1883 // get matrix data pointers 1564 1884 const double * elementByColumn = matrix_->getElements(); 1565 1885 const CoinBigIndex * columnStart = matrix_->getVectorStarts(); 1566 1886 const int * columnLength = matrix_->getVectorLengths(); 1887 int numberColumns = matrix_->getNumCols(); 1567 1888 int i; 1568 1889 for (i=0;i<numberColumns;i++) { … … 1588 1909 // Partial pricing 1589 1910 void 1590 ClpPackedMatrix::partialPricing(ClpSimplex * model, int start, int end,1911 ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 1591 1912 int & bestSequence, int & numberWanted) 1592 1913 { 1914 numberWanted=currentWanted_; 1915 int start = (int) (startFraction*numberActiveColumns_); 1916 int end = min((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_); 1593 1917 const double * element =matrix_->getElements(); 1594 1918 const int * row = matrix_->getIndices(); … … 1610 1934 int sequenceOut = model->sequenceOut(); 1611 1935 int saveSequence = bestSequence; 1936 int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 1937 int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_; 1612 1938 if (rowScale) { 1613 1939 // scaled … … 1696 2022 } 1697 2023 } 2024 if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) { 2025 // give up 2026 break; 2027 } 1698 2028 if (!numberWanted) 1699 2029 break; … … 1709 2039 } 1710 2040 reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence]; 2041 savedBestSequence_ = bestSequence; 2042 savedBestDj_ = reducedCost[savedBestSequence_]; 1711 2043 } 1712 2044 } else { … … 1793 2125 } 1794 2126 } 2127 if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) { 2128 // give up 2129 break; 2130 } 1795 2131 if (!numberWanted) 1796 2132 break; … … 1805 2141 } 1806 2142 reducedCost[bestSequence] = value; 1807 } 1808 } 2143 savedBestSequence_ = bestSequence; 2144 savedBestDj_ = reducedCost[savedBestSequence_]; 2145 } 2146 } 2147 currentWanted_=numberWanted; 1809 2148 } 1810 2149 // Sets up an effective RHS … … 1817 2156 rhsOffset(model,true); 1818 2157 } 2158 // makes sure active columns correct 2159 int 2160 ClpPackedMatrix::refresh(ClpSimplex * model) 2161 { 2162 numberActiveColumns_ = matrix_->getNumCols(); 2163 return 0; 2164 } 1819 2165 -
trunk/ClpPlusMinusOneMatrix.cpp
r309 r336 571 571 const int * whichRow = rowArray->getIndices(); 572 572 bool packed = rowArray->packedMode(); 573 if (numberInRowArray>2 ||y->getNumElements()) {573 if (numberInRowArray>2) { 574 574 // do by rows 575 575 int iRow; 576 576 double * markVector = y->denseVector(); // probably empty .. but 577 int * mark = y->getIndices(); 578 int numberOriginal=y->getNumElements(); 577 int numberOriginal=0; 579 578 int i; 580 579 if (packed) { 581 assert(!numberOriginal);582 580 numberNonZero=0; 583 581 // and set up mark as char array … … 628 626 } 629 627 } 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; 637 629 // and set up mark as char array 638 630 char * marked = (char *) markVector; … … 1423 1415 // Partial pricing 1424 1416 void 1425 ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, int start, int end,1417 ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 1426 1418 int & bestSequence, int & numberWanted) 1427 1419 { 1420 numberWanted=currentWanted_; 1421 int start = (int) (startFraction*numberColumns_); 1422 int end = min((int) (endFraction*numberColumns_+1),numberColumns_); 1428 1423 CoinBigIndex j; 1429 1424 double tolerance=model->currentDualTolerance(); … … 1547 1542 } 1548 1543 reducedCost[bestSequence] = value; 1549 } 1544 savedBestSequence_ = bestSequence; 1545 savedBestDj_ = reducedCost[savedBestSequence_]; 1546 } 1547 currentWanted_=numberWanted; 1550 1548 } 1551 1549 // Allow any parts of a created CoinMatrix to be deleted -
trunk/ClpPredictorCorrector.cpp
r328 r336 9 9 10 10 #include "CoinPragma.hpp" 11 11 static int xxxxxx=-1; 12 12 #include <math.h> 13 13 … … 99 99 double * savePrimal=NULL; 100 100 int numberTotal = numberRows_+numberColumns_; 101 bool tryJustPredictor=false; 101 102 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 102 157 complementarityGap_=complementarityGap(numberComplementarityPairs_,0); 103 158 handler_->message(CLP_BARRIER_ITERATION,messages_) … … 187 242 <<CoinMessageEol; 188 243 } 244 //if (complementarityGap_>=0.98*lastComplementarityGap) { 245 //tryJustPredictor=true; 246 //printf("trying just predictor\n"); 247 //} 189 248 if (complementarityGap_>=1.05*lastComplementarityGap) { 190 249 handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_) … … 340 399 double part2=nextGap/complementarityGap_; 341 400 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 342 408 } else { 343 409 double phi; … … 364 430 //std::cout<<"gap part "<< 365 431 //complementarityGap_*(beta2-tau)<<" after adding product = "<<xx<<std::endl; 366 //xx=-1.0;432 xx=-1.0; 367 433 if (xx>0.0) { 368 434 double saveMu = mu_; … … 391 457 } 392 458 //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) { 394 461 doCorrector=false; 395 462 bestNextGap=COIN_DBL_MAX; … … 405 472 <<CoinMessageEol; 406 473 phase=2; 474 tryJustPredictor=false; 407 475 } else { 408 476 phase=1; … … 464 532 } /* endwhile */ 465 533 } /* 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 } 466 573 numberFixed=updateSolution(); 467 574 numberFixedTotal+=numberFixed; … … 695 802 } 696 803 } else { 697 // iColumnust primal dual804 //just primal dual 698 805 for (int iColumn=0;iColumn<numberTotal;iColumn++) { 699 806 if (!flagged(iColumn)) { … … 901 1008 #ifndef KKT 902 1009 double maximumRHS = maximumAbsElement(updateRegion_,numberRows_); 1010 double saveMaximum = maximumRHS; 903 1011 double scale=1.0; 904 1012 double unscale=1.0; … … 985 1093 } 986 1094 relativeError = maximumRHSError/solutionNorm_; 1095 relativeError = maximumRHSError/saveMaximum; 987 1096 if (relativeError>tryError) 988 1097 relativeError=tryError; … … 1157 1266 } 1158 1267 // modify fixed RHS 1159 multiplyAdd(dj_+numberColumns_,numberRows_, -1.0,rhsFixRegion_,0.0);1268 multiplyAdd(dj_+numberColumns_,numberRows_,1.0,rhsFixRegion_,0.0); 1160 1269 // create plausible RHS? 1161 1270 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_); 1164 1273 rhsNorm_=maximumAbsElement(errorRegion_,numberRows_); 1165 1274 if (rhsNorm_<1.0) { … … 1229 1338 <<initialValue<<objectiveNorm_ 1230 1339 <<CoinMessageEol; 1231 int strategy= 1;1340 int strategy=0; 1232 1341 double extra=1.0e-10; 1233 1342 double largeGap=1.0e15; 1234 1343 double safeObjectiveValue=2.0*objectiveNorm_; 1344 //safeObjectiveValue=1.0+objectiveNorm_; 1345 //printf("temp safe obj value of %g\n",safeObjectiveValue); 1235 1346 double safeFree=1.0e-1*initialValue; 1236 1347 safeFree=1.0; … … 1245 1356 double low; 1246 1357 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; 1249 1360 double setPrimal=initialValue; 1250 1361 if (strategy==0) { … … 2253 2364 <<largestDiagonal<<smallestDiagonal 2254 2365 <<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 2255 2376 // update rhs region 2256 2377 multiplyAdd(work1+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0); … … 2274 2395 solutionNorm_ = norm; 2275 2396 //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); 2277 2398 matrix_->times(1.0,solution_,errorRegion_); 2278 2399 maximumDualError_=maximumDualError; -
trunk/ClpPrimalColumnSteepest.cpp
r328 r336 12 12 #include "CoinHelperFunctions.hpp" 13 13 #include <stdio.h> 14 15 14 16 15 //############################################################################# … … 198 197 */ 199 198 // Look at gub 199 #if 1 200 200 model_->clpMatrix()->dualExpanded(model_,updates,NULL,4); 201 #else 202 updates->clear(); 203 model_->computeDuals(NULL); 204 #endif 201 205 if (updates->getNumElements()>1) { 202 206 // would have to have two goes for devex, three for steepest … … 514 518 iSequence = index[i]; 515 519 double value = infeas[iSequence]; 520 double weight = weights_[iSequence]; 516 521 if (value>tolerance) { 517 double weight = weights_[iSequence];518 522 //weight=1.0; 519 523 if (value>bestDj*weight) { … … 527 531 } 528 532 } 529 }530 numberWanted--;533 numberWanted--; 534 } 531 535 if (!numberWanted) 532 536 break; … … 548 552 } 549 553 } 550 }551 numberWanted--;554 numberWanted--; 555 } 552 556 if (!numberWanted) 553 557 break; … … 557 561 break; 558 562 } 563 model_->clpMatrix()->setSavedBestSequence(bestSequence); 564 if (bestSequence>=0) 565 model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]); 559 566 if (sequenceOut>=0) { 560 567 infeas[sequenceOut]=saveOutInfeasibility; … … 985 992 double pivotSquared; 986 993 int iSequence = index[j]; 987 double value = reducedCost[iSequence];988 994 double value2 = updateBy[j]; 989 value -= value2; 990 reducedCost[iSequence] = value; 995 updateBy[j]=0.0; 991 996 ClpSimplex::Status status = model_->getStatus(iSequence+addSequence); 992 modification = other[iSequence]; 993 updateBy[j]=0.0; 997 double value; 994 998 995 999 switch(status) { … … 997 1001 case ClpSimplex::basic: 998 1002 infeasible_->zero(iSequence+addSequence); 1003 reducedCost[iSequence] = 0.0; 999 1004 case ClpSimplex::isFixed: 1000 1005 break; 1001 1006 case ClpSimplex::isFree: 1002 1007 case ClpSimplex::superBasic: 1008 value = reducedCost[iSequence] - value2; 1009 modification = other[iSequence]; 1003 1010 thisWeight = weight[iSequence]; 1004 1011 // row has -1 … … 1007 1014 1008 1015 thisWeight += pivotSquared * devex_ + pivot * modification; 1016 reducedCost[iSequence] = value; 1009 1017 if (thisWeight<TRY_NORM) { 1010 1018 if (mode_==1) { … … 1033 1041 break; 1034 1042 case ClpSimplex::atUpperBound: 1043 value = reducedCost[iSequence] - value2; 1044 modification = other[iSequence]; 1035 1045 thisWeight = weight[iSequence]; 1036 1046 // row has -1 … … 1039 1049 1040 1050 thisWeight += pivotSquared * devex_ + pivot * modification; 1051 reducedCost[iSequence] = value; 1041 1052 if (thisWeight<TRY_NORM) { 1042 1053 if (mode_==1) { … … 1063 1074 break; 1064 1075 case ClpSimplex::atLowerBound: 1076 value = reducedCost[iSequence] - value2; 1077 modification = other[iSequence]; 1065 1078 thisWeight = weight[iSequence]; 1066 1079 // row has -1 … … 1069 1082 1070 1083 thisWeight += pivotSquared * devex_ + pivot * modification; 1084 reducedCost[iSequence] = value; 1071 1085 if (thisWeight<TRY_NORM) { 1072 1086 if (mode_==1) { … … 2723 2737 /*if (model_->numberIterations()%100==0) 2724 2738 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]); 2725 2743 2726 2744 #ifdef CLP_DEBUG … … 3521 3539 const double * cost = model_->costRegion(1); 3522 3540 3541 model_->clpMatrix()->setOriginalWanted(numberWanted); 3542 model_->clpMatrix()->setCurrentWanted(numberWanted); 3523 3543 int iPassR=0,iPassC=0; 3524 3544 // Setup two passes … … 3534 3554 startR[0]=(int) dstart; 3535 3555 startR[3]=startR[0]; 3536 intstartC[4];3537 startC[1]= numberColumns;3556 double startC[4]; 3557 startC[1]=1.0; 3538 3558 startC[2]=0; 3539 3559 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; 3543 3562 reducedCost = model_->djRegion(1); 3544 3563 int sequenceOut = model_->sequenceOut(); … … 3549 3568 int i; 3550 3569 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]); 3552 3571 printf("\n"); 3553 3572 } … … 3555 3574 bool finishedR=false,finishedC=false; 3556 3575 bool doingR = randomR>randomC; 3557 doingR=false;3576 //doingR=false; 3558 3577 int saveNumberWanted = numberWanted; 3559 3578 while (!finishedR||!finishedC) { … … 3639 3658 reducedCost[bestSequence] = cost[bestSequence] +duals[bestSequence-numberColumns]; 3640 3659 bestDj=fabs(reducedCost[bestSequence]); 3641 } 3660 model_->clpMatrix()->setSavedBestSequence(bestSequence); 3661 model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]); 3662 } 3663 model_->clpMatrix()->setCurrentWanted(numberWanted); 3642 3664 if (!numberWanted) 3643 3665 break; … … 3657 3679 int saveSequence = bestSequence; 3658 3680 // 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 3661 3689 model_->clpMatrix()->partialPricing(model_,start,end,bestSequence,numberWanted); 3662 numberLook -= (end-start); 3690 numberWanted=model_->clpMatrix()->currentWanted(); 3691 numberLook -= (int) ((end-start)*numberColumns); 3663 3692 if (numberLook<0&&(10*(saveNumberWanted-numberWanted)>saveNumberWanted)) 3664 3693 numberWanted=0; // give up … … 3671 3700 doingR=true; 3672 3701 // 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) { 3676 3704 if (iPassC) 3677 3705 finishedC=true; -
trunk/ClpSimplex.cpp
r335 r336 43 43 remainingDualInfeasibility_(0.0), 44 44 largeValue_(1.0e15), 45 objectiveScale_(1.0), 46 rhsScale_(1.0), 45 47 largestPrimalError_(0.0), 46 48 largestDualError_(0.0), … … 110 112 incomingInfeasibility_(1.0), 111 113 allowedInfeasibility_(10.0), 114 automaticScale_(0), 112 115 progress_(NULL) 113 116 { … … 148 151 remainingDualInfeasibility_(0.0), 149 152 largeValue_(1.0e15), 153 objectiveScale_(1.0), 154 rhsScale_(1.0), 150 155 largestPrimalError_(0.0), 151 156 largestDualError_(0.0), … … 215 220 incomingInfeasibility_(1.0), 216 221 allowedInfeasibility_(10.0), 222 automaticScale_(0), 217 223 progress_(NULL) 218 224 { … … 348 354 // now check solutions 349 355 checkPrimalSolution( rowActivityWork_, columnActivityWork_); 350 objectiveValue_ += objectiveModification ;356 objectiveValue_ += objectiveModification/(objectiveScale_*rhsScale_); 351 357 checkDualSolution(); 352 358 if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2|| … … 423 429 } 424 430 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 426 447 // Ftran adjusted RHS and iterate to improve accuracy 427 448 double lastError=COIN_DBL_MAX; … … 431 452 CoinIndexedVector * lastVector = &previousVector; 432 453 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 433 463 bool goodSolution=true; 434 464 for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) { … … 525 555 } 526 556 } 557 #ifdef CLP_DEBUG 558 if (numberIterations_==-3840) { 559 exit(77); 560 } 561 #endif 527 562 } 528 563 // now dual side … … 1247 1282 // Update hidden stuff e.g. effective RHS and gub 1248 1283 matrix_->updatePivot(this,oldIn,oldOut); 1249 objectiveValue_ += objectiveChange ;1284 objectiveValue_ += objectiveChange/(objectiveScale_*rhsScale_); 1250 1285 handler_->message(CLP_SIMPLEX_HOUSE2,messages_) 1251 1286 <<numberIterations_<<objectiveValue() … … 1263 1298 //exit(77); 1264 1299 // 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, 1266 1304 directionIn_,directionOut_); 1267 1305 if (cycle>0) { … … 1323 1361 remainingDualInfeasibility_(0.0), 1324 1362 largeValue_(1.0e15), 1363 objectiveScale_(1.0), 1364 rhsScale_(1.0), 1325 1365 largestPrimalError_(0.0), 1326 1366 largestDualError_(0.0), … … 1390 1430 incomingInfeasibility_(1.0), 1391 1431 allowedInfeasibility_(10.0), 1432 automaticScale_(0), 1392 1433 progress_(NULL) 1393 1434 { … … 1422 1463 remainingDualInfeasibility_(0.0), 1423 1464 largeValue_(1.0e15), 1465 objectiveScale_(1.0), 1466 rhsScale_(1.0), 1424 1467 largestPrimalError_(0.0), 1425 1468 largestDualError_(0.0), … … 1489 1532 incomingInfeasibility_(1.0), 1490 1533 allowedInfeasibility_(10.0), 1534 automaticScale_(0), 1491 1535 progress_(NULL) 1492 1536 { … … 1585 1629 remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_; 1586 1630 largeValue_ = rhs.largeValue_; 1631 objectiveScale_ = rhs.objectiveScale_; 1632 rhsScale_ = rhs.rhsScale_; 1587 1633 largestPrimalError_ = rhs.largestPrimalError_; 1588 1634 largestDualError_ = rhs.largestDualError_; … … 1630 1676 incomingInfeasibility_ = rhs.incomingInfeasibility_; 1631 1677 allowedInfeasibility_ = rhs.allowedInfeasibility_; 1678 automaticScale_ = rhs.automaticScale_; 1632 1679 if (rhs.progress_) 1633 1680 progress_ = new ClpSimplexProgress(*rhs.progress_); … … 1814 1861 } 1815 1862 } 1863 objectiveValue_ /= (objectiveScale_*rhsScale_); 1816 1864 } 1817 1865 void … … 2231 2279 if (matrix_->scale(this)) 2232 2280 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 } 2233 2389 } 2234 2390 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_; 2239 2392 // but also scale by scale factors 2240 2393 // not really sure about row scaling … … 2253 2406 } 2254 2407 } 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 } 2281 2460 } 2282 2461 } … … 2366 2545 ray_=NULL; 2367 2546 } 2547 // set upperOut_ to furthest away from bound so can use in dual for dualBound_ 2548 upperOut_=0.0; 2368 2549 // ray may be null if in branch and bound 2369 2550 if (rowScale_) { … … 2373 2554 int numberDualScaled=0; 2374 2555 int numberDualUnscaled=0; 2556 double scaleC = 1.0/objectiveScale_; 2557 double scaleR = 1.0/rhsScale_; 2375 2558 for (i=0;i<numberColumns_;i++) { 2376 2559 double scaleFactor = columnScale_[i]; 2377 2560 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; 2383 2574 double value = columnActivity_[i]; 2384 2575 if (value<columnLower_[i]-primalTolerance_) … … 2391 2582 if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_) 2392 2583 numberDualScaled++; 2393 reducedCost_[i] = valueScaledDual/scaleFactor;2584 reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor; 2394 2585 double valueDual = reducedCost_[i]; 2395 2586 if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_) … … 2401 2592 double scaleFactor = rowScale_[i]; 2402 2593 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; 2408 2607 double value = rowActivity_[i]; 2409 2608 if (value<rowLower_[i]-primalTolerance_) … … 2416 2615 if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_) 2417 2616 numberDualScaled++; 2418 dual_[i] = valueScaledDual*scaleFactor ;2617 dual_[i] = valueScaledDual*scaleFactor*scaleC; 2419 2618 double valueDual = dual_[i]; 2420 2619 if (rowObjective_) … … 2446 2645 } 2447 2646 } 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 } 2448 2733 } else { 2449 2734 if (columnActivityWork_) { 2450 2735 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 } 2451 2745 columnActivity_[i] = columnActivityWork_[i]; 2452 2746 reducedCost_[i] = reducedCostWork_[i]; 2453 2747 } 2454 2748 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 } 2455 2758 rowActivity_[i] = rowActivityWork_[i]; 2456 2759 } 2457 2760 } 2761 } 2762 // switch off scalefactor if auto 2763 if (automaticScale_) { 2764 rhsScale_=1.0; 2765 objectiveScale_=1.0; 2458 2766 } 2459 2767 if (optimizationDirection_!=1.0) { … … 2468 2776 if(getRidOfFactorizationData>=0) 2469 2777 gutsOfDelete(getRidOfFactorizationData+1); 2778 // get rid of data 2779 matrix_->generalExpanded(this,13,scalingFlag_); 2470 2780 } 2471 2781 void … … 3120 3430 // check which algorithms allowed 3121 3431 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_; 3123 3437 returnCode = ((ClpSimplexPrimal *) this)->dual(0); 3124 else 3438 dualBound_=saveBound; 3439 } else { 3125 3440 returnCode = ((ClpSimplexDual *) this)->primal(0); 3441 } 3126 3442 setInitialDenseFactorization(denseFactorization); 3127 3443 perturbation_=savePerturbation; … … 3191 3507 perturbation_ = otherModel.perturbation_; 3192 3508 specialOptions_ = otherModel.specialOptions_; 3509 automaticScale_ = otherModel.automaticScale_; 3510 objectiveScale_ = otherModel.objectiveScale_; 3511 rhsScale_ = otherModel.rhsScale_; 3193 3512 } 3194 3513 typedef struct { … … 4493 4812 4494 4813 // 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_, 4497 4818 alpha_); 4498 4819 bool takePivot=true; -
trunk/ClpSimplexDual.cpp
r333 r336 573 573 #endif 574 574 #if 0 575 if (!factorization_->pivots()){ 575 // if (factorization_->pivots()){ 576 { 576 577 int iPivot; 577 578 double * array = rowArray_[3]->denseVector(); … … 596 597 double upperValue=upper_[iSequence]; 597 598 double value=solution_[iSequence]; 598 if(getStatus(iSequence)!=basic ) {599 if(getStatus(iSequence)!=basic&&getStatus(iSequence)!=isFree) { 599 600 assert(lowerValue>-1.0e20); 600 601 assert(upperValue<1.0e20); … … 843 844 assert(fabs(dualOut_)<1.0e50); 844 845 // if stable replace in basis 845 int updateStatus = factorization_->replaceColumn(rowArray_[2], 846 int updateStatus = factorization_->replaceColumn(this, 847 rowArray_[2], 848 rowArray_[1], 846 849 pivotRow_, 847 850 alpha_); … … 2958 2961 if (rowScale_) { 2959 2962 if (rowLowerWork_[iRow]>-1.0e50) 2960 rowLowerWork_[iRow] *= rowScale_[iRow] ;2963 rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_; 2961 2964 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_; 2963 2971 } 2964 2972 } else { … … 2969 2977 double multiplier = 1.0/columnScale_[iSequence]; 2970 2978 if (columnLowerWork_[iSequence]>-1.0e50) 2971 columnLowerWork_[iSequence] *= multiplier ;2979 columnLowerWork_[iSequence] *= multiplier*rhsScale_; 2972 2980 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_; 2974 2987 } 2975 2988 } … … 3033 3046 int minLength=numberRows_; 3034 3047 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; 3035 3058 int numberNonZero=0; 3036 3059 if (!numberIterations_&&perturbation_==50) { … … 3068 3091 //exit(0); 3069 3092 #endif 3093 //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_), 3094 // elementRatio); 3070 3095 //number=0; 3071 if (number*4>numberColumns_ ) {3096 if (number*4>numberColumns_||elementRatio>1.0e12) { 3072 3097 perturbation_=100; 3073 3098 return; // good enough … … 3421 3446 columnUpper_[iColumn] = newUpper[i]; 3422 3447 if (scalingFlag_<=0) 3423 upper_[iColumn] = newUpper[i] ;3448 upper_[iColumn] = newUpper[i]*rhsScale_; 3424 3449 else 3425 upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale3450 upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale 3426 3451 // Start of fast iterations 3427 3452 int status = fastDual(alwaysFinish); … … 3482 3507 columnLower_[iColumn] = newLower[i]; 3483 3508 if (scalingFlag_<=0) 3484 lower_[iColumn] = newLower[i] ;3509 lower_[iColumn] = newLower[i]*rhsScale_; 3485 3510 else 3486 lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale3511 lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale 3487 3512 // Start of fast iterations 3488 3513 status = fastDual(alwaysFinish); -
trunk/ClpSimplexPrimal.cpp
r335 r336 1 1 // Copyright (C) 2002, International Business Machines 2 2 // Corporation and others. All Rights Reserved. 3 4 3 5 4 /* Notes on implementation of primal simplex algorithm. … … 745 744 } 746 745 // 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) { 749 749 //may need infeasiblity cost changed 750 750 // we can see if we can construct a ray … … 1071 1071 // we also need somewhere for effective rhs 1072 1072 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(); 1074 1076 //int numberFlip=0; // Those which may change if flips 1075 1077 … … 1112 1114 bool gotList=false; 1113 1115 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 1114 1133 while (!gotList) { 1115 1134 pivotOne=-1; … … 1128 1147 double oldValue = solution_[iPivot]; 1129 1148 // get where in bound sequence 1149 // note that after this alpha is actually fabs(alpha) 1130 1150 if (alpha>0.0) { 1131 1151 // basic variable going towards lower bound … … 1136 1156 double bound = upper_[iPivot]; 1137 1157 oldValue = bound-oldValue; 1158 alpha = - alpha; 1138 1159 } 1139 1160 1140 double value = oldValue-tentativeTheta* fabs(alpha);1161 double value = oldValue-tentativeTheta*alpha; 1141 1162 assert (oldValue>=-tolerance); 1142 1163 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; 1146 1167 pivotOne=numberRemaining; 1147 1168 } … … 1149 1170 spare[numberRemaining]=alpha; 1150 1171 rhs[numberRemaining]=oldValue; 1172 indexPoint[numberRemaining]=iIndex; 1151 1173 index[numberRemaining++]=iRow; 1152 totalThru += fabs(alpha);1174 totalThru += alpha; 1153 1175 setActive(iRow); 1154 1176 //} else if (value<primalTolerance_*1.002) { … … 1168 1190 } 1169 1191 totalThru=0.0; 1170 // we need to keep where rhs non-zeros are1171 int numberInRhs=numberRemaining;1172 memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));1173 rhsArray->setNumElements(numberInRhs);1174 rhsArray->setPacked();1175 1192 1176 1193 theta_=maximumMovement; … … 1184 1201 // First check if previously chosen one will work 1185 1202 if (pivotOne>=0&&0) { 1186 double thruCost = infeasibilityCost_* fabs(spare[pivotOne]);1203 double thruCost = infeasibilityCost_*spare[pivotOne]; 1187 1204 if (thruCost>=0.99*fabs(dualIn_)) 1188 1205 printf("Could pivot on %d as change %g dj %g\n", … … 1190 1207 double alpha = spare[pivotOne]; 1191 1208 double oldValue = rhs[pivotOne]; 1192 theta_ = oldValue/ fabs(alpha);1209 theta_ = oldValue/alpha; 1193 1210 pivotRow_=pivotOne; 1194 1211 // Stop loop … … 1203 1220 for (iIndex=0;iIndex<numberRemaining;iIndex++) { 1204 1221 1205 double alpha = fabs(spare[iIndex]);1222 double alpha = spare[iIndex]; 1206 1223 double oldValue = rhs[iIndex]; 1207 1224 double value = oldValue-upperTheta*alpha; … … 1221 1238 double alpha = spare[iIndex]; 1222 1239 double oldValue = rhs[iIndex]; 1223 double value = oldValue-upperTheta* fabs(alpha);1240 double value = oldValue-upperTheta*alpha; 1224 1241 1225 1242 if (value<=0||iBest==iIndex) { 1226 1243 // 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]); 1228 1246 setActive(iRow); 1229 if ( fabs(alpha)>bestPivot) {1230 bestPivot= fabs(alpha);1247 if (alpha>bestPivot) { 1248 bestPivot=alpha; 1231 1249 theta_ = oldValue/bestPivot; 1232 1250 pivotRow_=iIndex; … … 1254 1272 if (bestPivot>bestEverPivot) 1255 1273 bestEverPivot=bestPivot; 1256 } 1257 } 1274 } } 1258 1275 // can get here without pivotRow_ set but with lastPivotRow 1259 1276 if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) { … … 1274 1291 int position=pivotRow_; // position in list 1275 1292 pivotRow_=index[position]; 1276 //alpha_ = work[pivotRow_]; 1277 alpha_ = way*spare[position]; 1293 alpha_=work[indexPoint[position]]; 1278 1294 // translate to sequence 1279 1295 sequenceOut_ = pivotVariable_[pivotRow_]; … … 1291 1307 // will we need to increase tolerance 1292 1308 //#define CLP_DEBUG 1293 #ifdef CLP_DEBUG1294 bool found=false;1295 #endif1296 1309 double largestInfeasibility = primalTolerance_; 1297 1310 if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) { … … 1299 1312 for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) { 1300 1313 largestInfeasibility = max (largestInfeasibility, 1301 -(rhs[iIndex]- fabs(spare[iIndex])*theta_));1314 -(rhs[iIndex]-spare[iIndex]*theta_)); 1302 1315 } 1303 1316 #define CLP_DEBUG … … 1365 1378 1366 1379 // put back original bounds etc 1380 memcpy(rhsArray->getIndices(),index,numberRemaining*sizeof(int)); 1381 rhsArray->setNumElements(numberRemaining); 1382 rhsArray->setPacked(); 1367 1383 nonLinearCost_->goBackAll(rhsArray); 1368 1369 1384 rhsArray->clear(); 1370 1385 … … 1613 1628 if (!numberIterations_) 1614 1629 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; 1615 1640 if (!numberIterations_&&perturbation_==50) { 1616 1641 // See if we need to perturb … … 1640 1665 } 1641 1666 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) { 1643 1670 perturbation_=100; 1644 1671 return; // good enough … … 1695 1722 if (getRowStatus(iSequence)==basic) 1696 1723 number++; 1724 } 1725 if (rhsScale_>100.0) { 1726 // tone down perturbation 1727 maximumFraction *= 0.1; 1697 1728 } 1698 1729 if (number!=numberRows_) … … 2166 2197 int updateStatus = matrix_->generalExpanded(this,3,updateType); 2167 2198 if (updateType>=0) 2168 updateStatus = factorization_->replaceColumn(rowArray_[2], 2199 updateStatus = factorization_->replaceColumn(this, 2200 rowArray_[2], 2201 rowArray_[1], 2169 2202 pivotRow_, 2170 2203 alpha_); … … 2317 2350 nonLinearCost_->setOne(sequenceIn_,valueIn_); 2318 2351 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 } 2325 2358 #endif 2326 2359 if (whatNext==1) { -
trunk/ClpSimplexPrimalQuadratic.cpp
r300 r336 1051 1051 if (cleanupIteration) 1052 1052 factorization_->relaxAccuracyCheck(1.0e3*saveCheck); 1053 updateStatus=factorization_->replaceColumn(rowArray_[2], 1053 updateStatus=factorization_->replaceColumn(this, 1054 rowArray_[2], 1055 rowArray_[1], 1054 1056 pivotRow_, 1055 1057 alpha_); -
trunk/Samples/decompose.cpp
r225 r336 32 32 // ken-18 33 33 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++) 34 38 rowBlock[iRow]=-1; 35 39 } else if (numberRows==810) { … … 142 146 sub[iBlock]=ClpSimplex(&model,numberRow2,whichRow, 143 147 numberColumn2,whichColumn); 144 #if 1148 #if 0 145 149 // temp 146 150 double * upper = sub[iBlock].columnUpper(); … … 219 223 master.scaling(false); 220 224 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"); 229 226 } 230 227 231 228 master.primal(); 229 int problemStatus = master.status(); // do here as can change (delcols) 232 230 if (master.numberIterations()==0&&iPass) 233 231 break; // finished … … 263 261 const double * dual=NULL; 264 262 bool deleteDual=false; 265 if ( master.isProvenOptimal()) {263 if (problemStatus==0) { 266 264 dual = master.dualRowSolution(); 267 } else if ( master.isProvenPrimalInfeasible()) {265 } else if (problemStatus==1) { 268 266 // could do composite objective 269 267 dual = master.infeasibilityRay(); … … 290 288 top[iBlock].transposeTimes(dual,objective2); 291 289 int i; 292 if ( master.isProvenOptimal()) {290 if (problemStatus==0) { 293 291 for (i=0;i<numberColumns2;i++) 294 292 objective2[i] = saveObj[i]-objective2[i]; … … 312 310 int number=start; 313 311 double dj = objValue; 314 if ( !master.isProvenOptimal())312 if (problemStatus) 315 313 dj=0.0; 316 314 double smallest=1.0e100; … … 377 375 } 378 376 } 377 delete [] saveObj; 379 378 } 380 379 if (deleteDual) … … 395 394 const double * fullLower = model.columnLower(); 396 395 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 } 397 405 int kColumn=0; 398 for (iRow=0;iRow<numberRows;iRow++)399 model.setRowStatus(iRow,ClpSimplex::superBasic);400 406 for (iColumn=0;iColumn<numberColumns;iColumn++) { 401 407 if (columnBlock[iColumn]==-1) { … … 413 419 if (whichBlock[i-numberMasterColumns]==iBlock) 414 420 sol[i] = solution[i]; 415 memset(lower,0, numberMasterRows*sizeof(double));421 memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double)); 416 422 master.times(1.0,sol,lower); 417 423 for (iRow=0;iRow<numberMasterRows;iRow++) { … … 433 439 sub[iBlock].primal(); 434 440 const double * subSolution = sub[iBlock].primalColumnSolution(); 441 const double * subRowSolution = sub[iBlock].primalRowSolution(); 435 442 // move solution 436 443 kColumn=0; … … 442 449 } 443 450 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); 444 459 } 445 460 for (iColumn=0;iColumn<numberColumns;iColumn++) { 446 461 if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&& 447 462 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); 451 474 model.primal(1); 452 475 delete [] sol; -
trunk/Samples/driver.cpp
r238 r336 20 20 exit(1); 21 21 } 22 22 #ifdef STYLE1 23 23 if (argc<3 ||!strstr(argv[2],"primal")) { 24 24 // Use the dual algorithm unless user said "primal" … … 27 27 model.initialPrimalSolve(); 28 28 } 29 #else 30 model.setLogLevel(8); 31 ClpSolve solvectl; 29 32 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 30 47 std::string modelName; 31 48 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 */ 3 11 4 12 #include "ClpSimplex.hpp" -
trunk/Samples/testGub.cpp
r332 r336 38 38 char * status = new char [numberColumns]; 39 39 char * rowStatus = new char[numberRows]; 40 int i;41 40 int n; 42 41 n = fread(solution,sizeof(double),numberColumns,fp); … … 49 48 assert (n==numberRows); 50 49 #if 0 50 int i; 51 51 for (i=0;i<numberColumns;i++) { 52 52 if (status[i]==0) -
trunk/Test/ClpMain.cpp
r327 r336 56 56 57 57 DUALTOLERANCE=1,PRIMALTOLERANCE,DUALBOUND,PRIMALWEIGHT,MAXTIME,OBJSCALE, 58 RHSSCALE, 58 59 59 60 LOGLEVEL=101,MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT, 60 61 61 62 DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR, 62 PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES, 63 PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE, 63 64 64 65 DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX, … … 567 568 break; 568 569 case OBJSCALE: 569 model->setOptimizationDirection(value); 570 model->setObjectiveScale(value); 571 break; 572 case RHSSCALE: 573 model->setRhsScale(value); 570 574 break; 571 575 default: … … 596 600 break; 597 601 case OBJSCALE: 598 value=model->optimizationDirection(); 602 value=model->objectiveScale(); 603 break; 604 case RHSSCALE: 605 value=model->rhsScale(); 599 606 break; 600 607 default: … … 891 898 ); 892 899 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++]= 893 910 ClpItem("barr!ier","Solve using primal dual predictor corrector algorithm", 894 911 BARRIER); … … 1135 1152 ClpItem("objective!Scale","Scale factor to apply to objective", 1136 1153 -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); 1138 1160 parameters[numberParameters++]= 1139 1161 ClpItem("passP!resolve","How many passes in presolve", … … 1258 1280 "Useful for testing if maximization works correctly" 1259 1281 ); 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); 1260 1291 parameters[numberParameters++]= 1261 1292 ClpItem("save!Model","Save model to binary file", … … 1548 1579 models[iModel].scaling(action); 1549 1580 break; 1581 case AUTOSCALE: 1582 models[iModel].setAutomaticScaling(action!=0); 1583 break; 1550 1584 case SPARSEFACTOR: 1551 1585 models[iModel].setSparseFactorization((1-action)!=0); -
trunk/include/ClpDynamicExampleMatrix.hpp
r335 r336 32 32 public: 33 33 /// Partial pricing 34 virtual void partialPricing(ClpSimplex * model, int start, intend,34 virtual void partialPricing(ClpSimplex * model, double start, double end, 35 35 int & bestSequence, int & numberWanted); 36 36 -
trunk/include/ClpDynamicMatrix.hpp
r335 r336 26 26 }; 27 27 /// Partial pricing 28 virtual void partialPricing(ClpSimplex * model, int start, intend,28 virtual void partialPricing(ClpSimplex * model, double start, double end, 29 29 int & bestSequence, int & numberWanted); 30 30 … … 71 71 mode=10 - return 1 if there may be changing bounds on variable (column generation) 72 72 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 73 75 */ 74 76 virtual int generalExpanded(ClpSimplex * model,int mode,int & number); … … 251 253 /// Saved best dual on gub row in pricing 252 254 double savedBestGubDual_; 253 /// Saved best dj in pricing254 double savedBestDj_;255 /// Saved best sequence in pricing256 int savedBestSequence_;257 255 /// Saved best set in pricing 258 256 int savedBestSet_; -
trunk/include/ClpFactorization.hpp
r225 r336 64 64 after factorization and thereafter re-factorize 65 65 partial update already in U */ 66 int replaceColumn ( CoinIndexedVector * regionSparse, 66 int replaceColumn ( const ClpSimplex * model, 67 CoinIndexedVector * regionSparse, 68 CoinIndexedVector * tableauColumn, 67 69 int pivotRow, 68 70 double pivotCheck , -
trunk/include/ClpGubDynamicMatrix.hpp
r335 r336 19 19 public: 20 20 /// Partial pricing 21 virtual void partialPricing(ClpSimplex * model, int start, intend,21 virtual void partialPricing(ClpSimplex * model, double start, double end, 22 22 int & bestSequence, int & numberWanted); 23 23 /** This is local to Gub to allow synchronization: … … 50 50 virtual void times(double scalar, 51 51 const double * x, double * y) const; 52 /** Just for debug 53 Returns number of primal infeasibilities. Recomputes keys 54 */ 55 virtual int checkFeasible() const; 52 56 //@} 53 57 … … 156 160 inline int numberElements() const 157 161 { return numberElements_;}; 158 /// Switches off dj checking each factorization (for BIG models)159 void switchOffCheck();160 162 /// Status region for gub slacks 161 163 inline unsigned char * gubRowStatus() const … … 197 199 /// Optional true upper bounds on sets 198 200 float * upperSet_; 199 /// Pointer back to model200 ClpSimplex * model_;201 201 /// size 202 202 int numberGubColumns_; -
trunk/include/ClpGubMatrix.hpp
r311 r336 49 49 int column, double multiplier) const; 50 50 /// Partial pricing 51 virtual void partialPricing(ClpSimplex * model, int start, intend,51 virtual void partialPricing(ClpSimplex * model, double start, double end, 52 52 int & bestSequence, int & numberWanted); 53 53 /// Returns number of hidden rows e.g. gub … … 119 119 mode=10 - return 1 if there may be changing bounds on variable (column generation) 120 120 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 121 123 */ 122 124 virtual int generalExpanded(ClpSimplex * model,int mode,int & number); … … 145 147 */ 146 148 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; 147 151 //@} 148 152 … … 263 267 inline int numberSets() const 264 268 { return numberSets_;}; 269 /// Switches off dj checking each factorization (for BIG models) 270 void ClpGubMatrix::switchOffCheck(); 265 271 //@} 266 272 … … 278 284 /// Sum of Primal infeasibilities using tolerance based on error in primals 279 285 double sumOfRelaxedPrimalInfeasibilities_; 286 /// Infeasibility weight when last full pass done 287 double infeasibilityWeight_; 280 288 /// Starts 281 289 int * start_; … … 308 316 // Reverse pointer from index to set 309 317 int * fromIndex_; 318 /// Pointer back to model 319 ClpSimplex * model_; 310 320 /// Number of dual infeasibilities 311 321 int numberDualInfeasibilities_; 312 322 /// Number of primal infeasibilities 313 323 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_; 314 330 /// Number of sets (gub rows) 315 331 int numberSets_; -
trunk/include/ClpMatrixBase.hpp
r328 r336 133 133 virtual int hiddenRows() const; 134 134 /// Partial pricing 135 virtual void partialPricing(ClpSimplex * model, int start, intend,135 virtual void partialPricing(ClpSimplex * model, double start, double end, 136 136 int & bestSequence, int & numberWanted); 137 137 /** expands an updated column to allow for extra rows which the main … … 178 178 mode=10 - return 1 if there may be changing bounds on variable (column generation) 179 179 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 180 182 */ 181 183 virtual int generalExpanded(ClpSimplex * model,int mode,int & number); … … 188 190 */ 189 191 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 ; 190 195 /// 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; 192 199 //@} 193 200 … … 276 283 inline void setSkipDualCheck(bool yes) 277 284 { 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;}; 278 326 //@} 279 327 … … 306 354 expensive */ 307 355 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_; 308 368 /// type (may be useful) 309 369 int type_; … … 312 372 /// If rhsOffset used this is refresh frequency (0==off) 313 373 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_; 314 382 /// whether to skip dual checks most of time 315 383 bool skipDualCheck_; -
trunk/include/ClpMessage.hpp
r263 r336 93 93 CLP_BARRIER_FEASIBLE, 94 94 CLP_BARRIER_STEP, 95 CLP_RIM_SCALE, 95 96 CLP_DUMMY_END 96 97 }; -
trunk/include/ClpNetworkMatrix.hpp
r286 r336 95 95 virtual bool canDoPartialPricing() const; 96 96 /// Partial pricing 97 virtual void partialPricing(ClpSimplex * model, int start, intend,97 virtual void partialPricing(ClpSimplex * model, double start, double end, 98 98 int & bestSequence, int & numberWanted); 99 99 //@} -
trunk/include/ClpPackedMatrix.hpp
r299 r336 54 54 /** Delete the columns whose indices are listed in <code>indDel</code>. */ 55 55 virtual void deleteCols(const int numDel, const int * indDel) 56 { matrix_->deleteCols(numDel,indDel); };56 { matrix_->deleteCols(numDel,indDel);numberActiveColumns_ = matrix_->getNumCols();}; 57 57 /** Delete the rows whose indices are listed in <code>indDel</code>. */ 58 58 virtual void deleteRows(const int numDel, const int * indDel) 59 { matrix_->deleteRows(numDel,indDel); };59 { matrix_->deleteRows(numDel,indDel);numberActiveColumns_ = matrix_->getNumCols();}; 60 60 /// Append Columns 61 61 virtual void appendCols(int number, const CoinPackedVectorBase * const * columns) 62 { matrix_->appendCols(number,columns); };62 { matrix_->appendCols(number,columns);numberActiveColumns_ = matrix_->getNumCols();}; 63 63 /// Append Rows 64 64 virtual void appendRows(int number, const CoinPackedVectorBase * const * rows) 65 { matrix_->appendRows(number,rows); };65 { matrix_->appendRows(number,rows);numberActiveColumns_ = matrix_->getNumCols();}; 66 66 /** Replace the elements of a vector. The indices remain the same. 67 67 This is only needed if scaling and a row copy is used. … … 125 125 virtual bool canDoPartialPricing() const; 126 126 /// Partial pricing 127 virtual void partialPricing(ClpSimplex * model, int start, intend,127 virtual void partialPricing(ClpSimplex * model, double start, double end, 128 128 int & bestSequence, int & numberWanted); 129 /// makes sure active columns correct 130 virtual int refresh(ClpSimplex * model); 129 131 //@} 130 132 … … 160 162 CoinIndexedVector * z) const; 161 163 /** 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>. 162 179 Can use y as temporary array (will be empty at end) 163 180 Note - If x packed mode - then z packed mode … … 235 252 /// Data 236 253 CoinPackedMatrix * matrix_; 254 /// number of active columns (normally same as number of columns) 255 int numberActiveColumns_; 237 256 /// Zero element flag - set true if any zero elements 238 257 bool zeroElements_; 258 /// Gaps flag - set true if column start and length don't say contiguous 259 bool hasGaps_; 239 260 //@} 240 261 }; -
trunk/include/ClpPlusMinusOneMatrix.hpp
r286 r336 191 191 virtual bool canDoPartialPricing() const; 192 192 /// Partial pricing 193 virtual void partialPricing(ClpSimplex * model, int start, intend,193 virtual void partialPricing(ClpSimplex * model, double start, double end, 194 194 int & bestSequence, int & numberWanted); 195 195 //@} -
trunk/include/ClpSimplex.hpp
r335 r336 527 527 inline int * pivotVariable() const 528 528 { return pivotVariable_;}; 529 /// Scaling of objective 12345.0 (auto), 0.0 (off), other user529 /// Scaling of objective 530 530 inline double objectiveScale() const 531 531 { return objectiveScale_;} ; 532 532 inline void setObjectiveScale(double value) 533 533 { 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;}; 534 544 /// Current dual tolerance 535 545 inline double currentDualTolerance() const … … 828 838 /// Scaling of objective 829 839 double objectiveScale_; 840 /// Scaling of rhs and bounds 841 double rhsScale_; 830 842 /// Largest error on Ax-b 831 843 double largestPrimalError_; … … 1000 1012 float incomingInfeasibility_; 1001 1013 float allowedInfeasibility_; 1014 /// Automatic scaling of objective and rhs and bounds 1015 int automaticScale_; 1002 1016 /// For dealing with all issues of cycling etc 1003 1017 ClpSimplexProgress * progress_;
Note: See TracChangeset
for help on using the changeset viewer.