00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 #include "CoinPragma.hpp"
00091
00092 #include <math.h>
00093
00094 #include "CoinHelperFunctions.hpp"
00095 #include "ClpSimplexDual.hpp"
00096 #include "ClpFactorization.hpp"
00097 #include "CoinPackedMatrix.hpp"
00098 #include "CoinIndexedVector.hpp"
00099 #include "CoinWarmStartBasis.hpp"
00100 #include "ClpDualRowDantzig.hpp"
00101 #include "ClpPlusMinusOneMatrix.hpp"
00102 #include "ClpMessage.hpp"
00103 #include <cfloat>
00104 #include <cassert>
00105 #include <string>
00106 #include <stdio.h>
00107 #include <iostream>
00108
00109 #define CHECK_DJ
00110
00111 int ClpSimplexDual::dual (int ifValuesPass )
00112 {
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 algorithm_ = -1;
00200
00201
00202 ClpDataSave data = saveData();
00203
00204
00205 int initialStatus=problemStatus_;
00206
00207
00208 double * saveDuals = NULL;
00209 if (ifValuesPass) {
00210 saveDuals = new double [numberRows_+numberColumns_];
00211 memcpy(saveDuals,dual_,numberRows_*sizeof(double));
00212 }
00213
00214
00215
00216
00217
00218
00219
00220
00221 if (!startup(0)) {
00222
00223
00224
00225 if (ifValuesPass) {
00226 if (problemStatus_&&perturbation_<100)
00227 perturb();
00228 int i;
00229 if (scalingFlag_>0) {
00230 for (i=0;i<numberRows_;i++) {
00231 dual_[i] = saveDuals[i]/rowScale_[i];
00232 }
00233 } else {
00234 memcpy(dual_,saveDuals,numberRows_*sizeof(double));
00235 }
00236
00237 for (i=0;i<numberRows_;i++) {
00238
00239 double value = dual_[i];
00240 value += rowObjectiveWork_[i];
00241 saveDuals[i+numberColumns_]=value;
00242 }
00243 ClpDisjointCopyN(objectiveWork_,numberColumns_,saveDuals);
00244 transposeTimes(-1.0,dual_,saveDuals);
00245 memcpy(dj_,saveDuals,(numberColumns_+numberRows_)*sizeof(double));
00246
00247 for (i=0;i<numberRows_+numberColumns_;i++)
00248 clearPivoted(i);
00249 int iRow;
00250 for (iRow=0;iRow<numberRows_;iRow++) {
00251 int iPivot=pivotVariable_[iRow];
00252 if (fabs(saveDuals[iPivot])>dualTolerance_) {
00253 if (getStatus(iPivot)!=isFree)
00254 setPivoted(iPivot);
00255 }
00256 }
00257 }
00258
00259 double objectiveChange;
00260 numberFake_ =0;
00261 changeBounds(true,NULL,objectiveChange);
00262
00263 int lastCleaned=0;
00264
00265 if (!ifValuesPass) {
00266
00267 if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
00268 problemStatus_=0;
00269 }
00270 if (problemStatus_<0&&perturbation_<100) {
00271 perturb();
00272
00273 gutsOfSolution(NULL,NULL);
00274 if (handler_->logLevel()>2) {
00275 handler_->message(CLP_SIMPLEX_STATUS,messages_)
00276 <<numberIterations_<<objectiveValue();
00277 handler_->printing(sumPrimalInfeasibilities_>0.0)
00278 <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
00279 handler_->printing(sumDualInfeasibilities_>0.0)
00280 <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
00281 handler_->printing(numberDualInfeasibilitiesWithoutFree_
00282 <numberDualInfeasibilities_)
00283 <<numberDualInfeasibilitiesWithoutFree_;
00284 handler_->message()<<CoinMessageEol;
00285 }
00286 }
00287
00288
00289 int factorType=0;
00290
00291 progress_->startCheck();
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 while (problemStatus_<0) {
00303 int iRow, iColumn;
00304
00305 for (iRow=0;iRow<4;iRow++) {
00306 rowArray_[iRow]->clear();
00307 }
00308
00309 for (iColumn=0;iColumn<2;iColumn++) {
00310 columnArray_[iColumn]->clear();
00311 }
00312
00313
00314
00315 matrix_->refresh(this);
00316
00317
00318 if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
00319 &&initialStatus!=10) {
00320 perturb();
00321
00322 gutsOfSolution(NULL,NULL);
00323 if (handler_->logLevel()>2) {
00324 handler_->message(CLP_SIMPLEX_STATUS,messages_)
00325 <<numberIterations_<<objectiveValue();
00326 handler_->printing(sumPrimalInfeasibilities_>0.0)
00327 <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
00328 handler_->printing(sumDualInfeasibilities_>0.0)
00329 <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
00330 handler_->printing(numberDualInfeasibilitiesWithoutFree_
00331 <numberDualInfeasibilities_)
00332 <<numberDualInfeasibilitiesWithoutFree_;
00333 handler_->message()<<CoinMessageEol;
00334 }
00335 }
00336
00337 statusOfProblemInDual(lastCleaned,factorType,progress_,saveDuals);
00338
00339 if (ifValuesPass&&
00340 progress_->lastIterationNumber(0)<0) {
00341 doEasyOnesInValuesPass(saveDuals);
00342 }
00343
00344
00345 factorType=1;
00346 if (data.sparseThreshold_) {
00347
00348 factorization_->sparseThreshold(0);
00349 factorization_->goSparse();
00350 }
00351
00352
00353 if (problemStatus_>=0)
00354 break;
00355
00356
00357 whileIterating(saveDuals);
00358 }
00359 }
00360
00361 assert (problemStatus_||!sumPrimalInfeasibilities_);
00362
00363
00364 finish();
00365 delete [] saveDuals;
00366
00367
00368 restoreData(data);
00369 return problemStatus_;
00370 }
00371
00372
00373
00374
00375
00376
00377
00378
00379 int
00380 ClpSimplexDual::whileIterating(double * & givenDuals)
00381 {
00382 #ifdef CLP_DEBUG
00383 int debugIteration=-1;
00384 #endif
00385 {
00386 int i;
00387 for (i=0;i<4;i++) {
00388 rowArray_[i]->clear();
00389 }
00390 for (i=0;i<2;i++) {
00391 columnArray_[i]->clear();
00392 }
00393 }
00394
00395 if (largestPrimalError_>10.0)
00396 factorization_->relaxAccuracyCheck(min(1.0e2,largestPrimalError_/10.0));
00397 else
00398 factorization_->relaxAccuracyCheck(1.0);
00399
00400
00401 int returnCode = -1;
00402
00403 #if 0
00404
00405 double averagePrimalInfeasibility = sumPrimalInfeasibilities_/
00406 ((double ) (numberPrimalInfeasibilities_+1));
00407 #endif
00408
00409
00410 factorization_->getWeights(rowArray_[0]->getIndices());
00411 CoinBigIndex * dubiousWeights = matrix_->dubiousWeights(this,rowArray_[0]->getIndices());
00412
00413 int * candidateList = NULL;
00414 int numberCandidates = 0;
00415 #ifdef CLP_DEBUG
00416 bool wasInValuesPass= (givenDuals!=NULL);
00417 #endif
00418 int candidate=-1;
00419 if (givenDuals) {
00420 candidateList = new int[numberRows_];
00421
00422 memcpy(dj_,givenDuals,(numberRows_+numberColumns_)*sizeof(double));
00423 int iRow;
00424 for (iRow=0;iRow<numberRows_;iRow++) {
00425 int iPivot=pivotVariable_[iRow];
00426 if (flagged(iPivot))
00427 continue;
00428 if (fabs(dj_[iPivot])>dualTolerance_) {
00429 if (pivoted(iPivot))
00430 candidateList[numberCandidates++]=iRow;
00431 } else {
00432 clearPivoted(iPivot);
00433 }
00434 }
00435
00436 if (!numberCandidates) {
00437 delete [] candidateList;
00438 delete [] givenDuals;
00439 givenDuals=NULL;
00440 candidateList=NULL;
00441 int iRow;
00442 for (iRow=0;iRow<numberRows_;iRow++) {
00443 int iPivot=pivotVariable_[iRow];
00444 clearPivoted(iPivot);
00445 }
00446 }
00447 }
00448 while (problemStatus_==-1) {
00449 #ifdef CLP_DEBUG
00450 if (givenDuals) {
00451 double value5=0.0;
00452 int i;
00453 for (i=0;i<numberRows_+numberColumns_;i++) {
00454 if (dj_[i]<-1.0e-6)
00455 if (upper_[i]<1.0e20)
00456 value5 += dj_[i]*upper_[i];
00457 else
00458 printf("bad dj %g on %d with large upper status %d\n",
00459 dj_[i],i,status_[i]&7);
00460 else if (dj_[i] >1.0e-6)
00461 if (lower_[i]>-1.0e20)
00462 value5 += dj_[i]*lower_[i];
00463 else
00464 printf("bad dj %g on %d with large lower status %d\n",
00465 dj_[i],i,status_[i]&7);
00466 }
00467 printf("Values objective Value %g\n",value5);
00468 }
00469 if ((handler_->logLevel()&32)&&wasInValuesPass) {
00470 double value5=0.0;
00471 int i;
00472 for (i=0;i<numberRows_+numberColumns_;i++) {
00473 if (dj_[i]<-1.0e-6)
00474 if (upper_[i]<1.0e20)
00475 value5 += dj_[i]*upper_[i];
00476 else if (dj_[i] >1.0e-6)
00477 if (lower_[i]>-1.0e20)
00478 value5 += dj_[i]*lower_[i];
00479 }
00480 printf("Values objective Value %g\n",value5);
00481 {
00482 int i;
00483 for (i=0;i<numberRows_+numberColumns_;i++) {
00484 int iSequence = i;
00485 double oldValue;
00486
00487 switch(getStatus(iSequence)) {
00488
00489 case basic:
00490 case ClpSimplex::isFixed:
00491 break;
00492 case isFree:
00493 case superBasic:
00494 abort();
00495 break;
00496 case atUpperBound:
00497 oldValue = dj_[iSequence];
00498
00499 assert (fabs(solution_[iSequence]-upper_[iSequence])<1.0e-7);
00500 break;
00501 case atLowerBound:
00502 oldValue = dj_[iSequence];
00503
00504 assert (fabs(solution_[iSequence]-lower_[iSequence])<1.0e-7);
00505 break;
00506 }
00507 }
00508 }
00509 }
00510 #endif
00511 #ifdef CLP_DEBUG
00512 {
00513 int i;
00514 for (i=0;i<4;i++) {
00515 rowArray_[i]->checkClear();
00516 }
00517 for (i=0;i<2;i++) {
00518 columnArray_[i]->checkClear();
00519 }
00520 }
00521 #endif
00522 #if CLP_DEBUG>2
00523
00524 if (numberIterations_>0&&numberIterations_<-801) {
00525 handler_->setLogLevel(63);
00526 double saveValue = objectiveValue_;
00527 double * saveRow1 = new double[numberRows_];
00528 double * saveRow2 = new double[numberRows_];
00529 memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
00530 memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
00531 double * saveColumn1 = new double[numberColumns_];
00532 double * saveColumn2 = new double[numberColumns_];
00533 memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
00534 memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
00535 gutsOfSolution(NULL,NULL);
00536 printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
00537 numberIterations_,
00538 saveValue,objectiveValue_,sumDualInfeasibilities_);
00539 if (saveValue>objectiveValue_+1.0e-2)
00540 printf("**bad**\n");
00541 memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
00542 memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
00543 memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
00544 memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
00545 delete [] saveRow1;
00546 delete [] saveRow2;
00547 delete [] saveColumn1;
00548 delete [] saveColumn2;
00549 objectiveValue_=saveValue;
00550 }
00551 #endif
00552 #if 0
00553 if (!factorization_->pivots()){
00554 int iPivot;
00555 double * array = rowArray_[3]->denseVector();
00556 int i;
00557 for (iPivot=0;iPivot<numberRows_;iPivot++) {
00558 int iSequence = pivotVariable_[iPivot];
00559 unpack(rowArray_[3],iSequence);
00560 factorization_->updateColumn(rowArray_[2],rowArray_[3]);
00561 assert (fabs(array[iPivot]-1.0)<1.0e-4);
00562 array[iPivot]=0.0;
00563 for (i=0;i<numberRows_;i++)
00564 assert (fabs(array[i])<1.0e-4);
00565 rowArray_[3]->clear();
00566 }
00567 }
00568 #endif
00569 #ifdef CLP_DEBUG
00570 {
00571 int iSequence, number=numberRows_+numberColumns_;
00572 for (iSequence=0;iSequence<number;iSequence++) {
00573 double lowerValue=lower_[iSequence];
00574 double upperValue=upper_[iSequence];
00575 double value=solution_[iSequence];
00576 if(getStatus(iSequence)!=basic) {
00577 assert(lowerValue>-1.0e20);
00578 assert(upperValue<1.0e20);
00579 }
00580 switch(getStatus(iSequence)) {
00581
00582 case basic:
00583 break;
00584 case isFree:
00585 case superBasic:
00586 break;
00587 case atUpperBound:
00588 assert (fabs(value-upperValue)<=primalTolerance_) ;
00589 break;
00590 case atLowerBound:
00591 case ClpSimplex::isFixed:
00592 assert (fabs(value-lowerValue)<=primalTolerance_) ;
00593 break;
00594 }
00595 }
00596 }
00597 if(numberIterations_==debugIteration) {
00598 printf("dodgy iteration coming up\n");
00599 }
00600 #endif
00601
00602
00603
00604 if (numberCandidates)
00605 candidate = candidateList[--numberCandidates];
00606 else
00607 candidate=-1;
00608 dualRow(candidate);
00609 if (pivotRow_>=0) {
00610
00611 handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
00612 <<pivotRow_
00613 <<CoinMessageEol;
00614
00615 dualRowPivot_->checkAccuracy();
00616
00617 double acceptablePivot=1.0e-7;
00618 if (factorization_->pivots()>10)
00619 acceptablePivot=1.0e-5;
00620 else if (factorization_->pivots()>5)
00621 acceptablePivot=1.0e-6;
00622
00623 if (candidate<0) {
00624
00625
00626 double direction=directionOut_;
00627 rowArray_[0]->createPacked(1,&pivotRow_,&direction);
00628 factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
00629
00630 matrix_->transposeTimes(this,-1.0,
00631 rowArray_[0],rowArray_[3],columnArray_[0]);
00632
00633 dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
00634 rowArray_[3],acceptablePivot,dubiousWeights);
00635 } else {
00636
00637 assert (upperOut_<1.0e50||lowerOut_>-1.0e50);
00638 if (directionOut_<0&&fabs(valueOut_-upperOut_)>dualBound_+primalTolerance_) {
00639 if (fabs(valueOut_-upperOut_)>fabs(valueOut_-lowerOut_))
00640 directionOut_=1;
00641 } else if (directionOut_>0&&fabs(valueOut_-upperOut_)<dualBound_+primalTolerance_) {
00642 if (fabs(valueOut_-upperOut_)>fabs(valueOut_-lowerOut_))
00643 directionOut_=-1;
00644 }
00645 double direction=directionOut_;
00646 rowArray_[0]->createPacked(1,&pivotRow_,&direction);
00647 factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
00648
00649 matrix_->transposeTimes(this,-1.0,
00650 rowArray_[0],rowArray_[3],columnArray_[0]);
00651 acceptablePivot *= 10.0;
00652
00653 checkPossibleValuesMove(rowArray_[0],columnArray_[0],
00654 acceptablePivot,NULL);
00655
00656
00657 if (directionOut_<0) {
00658 dualOut_ = valueOut_ - upperOut_;
00659 } else {
00660 dualOut_ = lowerOut_ - valueOut_;
00661 }
00662
00663
00664 bool normalIteration = (sequenceIn_!=sequenceOut_);
00665
00666 clearPivoted(sequenceOut_);
00667
00668 if (!numberCandidates) {
00669 int iRow;
00670 delete [] candidateList;
00671 delete [] givenDuals;
00672 candidate=-2;
00673 givenDuals=NULL;
00674 handler_->message(CLP_END_VALUES_PASS,messages_)
00675 <<numberIterations_;
00676 candidateList=NULL;
00677 for (iRow=0;iRow<numberRows_;iRow++) {
00678 int iPivot=pivotVariable_[iRow];
00679
00680 clearPivoted(iPivot);
00681 }
00682 }
00683 if (!normalIteration) {
00684
00685
00686 updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
00687 if (candidate==-2)
00688 problemStatus_=-2;
00689 continue;
00690 } else {
00691
00692 if (directionOut_<0) {
00693 dualOut_ = valueOut_ - upperOut_;
00694 } else {
00695 dualOut_ = lowerOut_ - valueOut_;
00696 }
00697 }
00698 }
00699 if (sequenceIn_>=0) {
00700
00701
00702 double btranAlpha = -alpha_*directionOut_;
00703 unpackPacked(rowArray_[1]);
00704 factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
00705
00706 alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
00707 rowArray_[2],
00708 rowArray_[1]);
00709
00710 #ifdef CLP_DEBUG
00711 if ((handler_->logLevel()&32))
00712 printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
00713 #endif
00714 double checkValue=1.0e-7;
00715
00716 if (largestPrimalError_>10.0)
00717 checkValue = min(1.0e-4,1.0e-8*largestPrimalError_);
00718 if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
00719 fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
00720 handler_->message(CLP_DUAL_CHECK,messages_)
00721 <<btranAlpha
00722 <<alpha_
00723 <<CoinMessageEol;
00724 if (factorization_->pivots()) {
00725 dualRowPivot_->unrollWeights();
00726 problemStatus_=-2;
00727 rowArray_[0]->clear();
00728 rowArray_[1]->clear();
00729 columnArray_[0]->clear();
00730 returnCode=-2;
00731 break;
00732 } else {
00733
00734 if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
00735 fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
00736 dualRowPivot_->unrollWeights();
00737
00738 char x = isColumn(sequenceOut_) ? 'C' :'R';
00739 handler_->message(CLP_SIMPLEX_FLAG,messages_)
00740 <<x<<sequenceWithin(sequenceOut_)
00741 <<CoinMessageEol;
00742 setFlagged(sequenceOut_);
00743 progress_->clearBadTimes();
00744 lastBadIteration_ = numberIterations_;
00745 rowArray_[0]->clear();
00746 rowArray_[1]->clear();
00747 columnArray_[0]->clear();
00748 continue;
00749 }
00750 }
00751 }
00752
00753 double objectiveChange=0.0;
00754
00755
00756
00757 int nswapped = 0;
00758
00759
00760 if (candidate==-1)
00761 nswapped = updateDualsInDual(rowArray_[0],columnArray_[0],
00762 rowArray_[2],theta_,
00763 objectiveChange,false);
00764 else
00765 updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
00766
00767
00768 if (nswapped) {
00769 factorization_->updateColumn(rowArray_[3],rowArray_[2]);
00770 dualRowPivot_->updatePrimalSolution(rowArray_[2],
00771 1.0,objectiveChange);
00772 #ifdef CLP_DEBUG
00773 double saveDualOut = dualOut_;
00774 #endif
00775
00776 valueOut_ = solution_[sequenceOut_];
00777 if (directionOut_<0) {
00778 dualOut_ = valueOut_ - upperOut_;
00779 } else {
00780 dualOut_ = lowerOut_ - valueOut_;
00781 }
00782 #if 0
00783 if (dualOut_<0.0) {
00784 #ifdef CLP_DEBUG
00785 if (handler_->logLevel()&32) {
00786 printf(" dualOut_ %g %g save %g\n",dualOut_,averagePrimalInfeasibility,saveDualOut);
00787 printf("values %g %g %g %g %g %g %g\n",lowerOut_,valueOut_,upperOut_,
00788 objectiveChange,);
00789 }
00790 #endif
00791 if (upperOut_==lowerOut_)
00792 dualOut_=0.0;
00793 }
00794 if(dualOut_<-max(1.0e-12*averagePrimalInfeasibility,1.0e-8)
00795 &&factorization_->pivots()>100&&
00796 getStatus(sequenceIn_)!=isFree) {
00797
00798 dualRowPivot_->unrollWeights();
00799 problemStatus_=-2;
00800 returnCode=-2;
00801 break;
00802 }
00803 #endif
00804 }
00805
00806 double movement = -dualOut_*directionOut_/alpha_;
00807
00808
00809 if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
00810 #ifdef CLP_DEBUG
00811 if (handler_->logLevel()&32)
00812 printf("movement %g, swap change %g, rest %g * %g\n",
00813 objectiveChange+fabs(movement*dualIn_),
00814 objectiveChange,movement,dualIn_);
00815 #endif
00816 if(factorization_->pivots()) {
00817
00818 dualRowPivot_->unrollWeights();
00819 problemStatus_=-2;
00820 returnCode=-2;
00821 break;
00822 }
00823 }
00824 assert(fabs(dualOut_)<1.0e50);
00825
00826 int updateStatus = factorization_->replaceColumn(rowArray_[2],
00827 pivotRow_,
00828 alpha_);
00829
00830 if (updateStatus==2&&
00831 !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
00832 updateStatus=4;
00833 if (updateStatus==1||updateStatus==4) {
00834
00835 if (factorization_->pivots()>5||updateStatus==4) {
00836 problemStatus_=-2;
00837 returnCode=-3;
00838 }
00839 } else if (updateStatus==2) {
00840
00841 dualRowPivot_->unrollWeights();
00842
00843 if (factorization_->pivots()) {
00844 problemStatus_=-2;
00845 returnCode=-2;
00846 break;
00847 } else {
00848
00849 char x = isColumn(sequenceOut_) ? 'C' :'R';
00850 handler_->message(CLP_SIMPLEX_FLAG,messages_)
00851 <<x<<sequenceWithin(sequenceOut_)
00852 <<CoinMessageEol;
00853 setFlagged(sequenceOut_);
00854 progress_->clearBadTimes();
00855 lastBadIteration_ = numberIterations_;
00856 rowArray_[0]->clear();
00857 rowArray_[1]->clear();
00858 columnArray_[0]->clear();
00859
00860
00861 double objectiveChange=0.0;
00862 updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
00863 0.0,objectiveChange,true);
00864 continue;
00865 }
00866 } else if (updateStatus==3) {
00867
00868
00869 if (factorization_->pivots()<
00870 0.5*factorization_->maximumPivots()&&
00871 factorization_->pivots()<200)
00872 factorization_->areaFactor(
00873 factorization_->areaFactor() * 1.1);
00874 problemStatus_=-2;
00875 } else if (updateStatus==5) {
00876 problemStatus_=-2;
00877 }
00878
00879 if (theta_<0.0&&candidate==-1) {
00880 #ifdef CLP_DEBUG
00881 if (handler_->logLevel()&32)
00882 printf("negative theta %g\n",theta_);
00883 #endif
00884 theta_=0.0;
00885 }
00886
00887 flipBounds(rowArray_[0],columnArray_[0],theta_);
00888
00889 dualRowPivot_->updatePrimalSolution(rowArray_[1],
00890 movement,
00891 objectiveChange);
00892 #ifdef CLP_DEBUG
00893 double oldobj=objectiveValue_;
00894 #endif
00895
00896 dualOut_ /= alpha_;
00897 dualOut_ *= -directionOut_;
00898
00899 dj_[sequenceIn_]=0.0;
00900 double oldValue=valueIn_;
00901 if (directionIn_==-1) {
00902
00903 valueIn_ = upperIn_+dualOut_;
00904 } else {
00905
00906 valueIn_ = lowerIn_+dualOut_;
00907 }
00908 objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
00909
00910
00911 if (directionOut_>0) {
00912 valueOut_ = lowerOut_;
00913 if (candidate==-1)
00914 dj_[sequenceOut_] = theta_;
00915 } else {
00916 valueOut_ = upperOut_;
00917 if (candidate==-1)
00918 dj_[sequenceOut_] = -theta_;
00919 }
00920 solution_[sequenceOut_]=valueOut_;
00921 int whatNext=housekeeping(objectiveChange);
00922
00923 originalBound(sequenceIn_);
00924 changeBound(sequenceOut_);
00925 #ifdef CLP_DEBUG
00926 if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
00927 printf("obj backwards %g %g\n",objectiveValue_,oldobj);
00928 #endif
00929 if (whatNext==1||candidate==-2) {
00930 problemStatus_ =-2;
00931 } else if (whatNext==2) {
00932
00933 problemStatus_= 3;
00934 returnCode=3;
00935 break;
00936 }
00937 } else {
00938
00939 #ifdef CLP_DEBUG
00940 if (handler_->logLevel()&32)
00941 printf("** no column pivot\n");
00942 #endif
00943 if (factorization_->pivots()<5) {
00944 problemStatus_=-4;
00945
00946 delete [] ray_;
00947 ray_ = new double [ numberRows_];
00948 rowArray_[0]->expand();
00949 ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
00950 }
00951 rowArray_[0]->clear();
00952 columnArray_[0]->clear();
00953 returnCode=1;
00954 break;
00955 }
00956 } else {
00957
00958 #ifdef CLP_DEBUG
00959 if (handler_->logLevel()&32)
00960 printf("** no row pivot\n");
00961 #endif
00962 if (!factorization_->pivots()) {
00963
00964
00965 int iRow;
00966 for (iRow=0;iRow<numberRows_;iRow++) {
00967 int iPivot=pivotVariable_[iRow];
00968 if (flagged(iPivot))
00969 break;
00970 }
00971 if (numberFake_||numberDualInfeasibilities_) {
00972
00973 problemStatus_=-5;
00974 } else {
00975 if (iRow<numberRows_) {
00976 #ifdef CLP_DEBUG
00977 std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
00978 printf("Probably infeasible - pivot was %g\n",alpha_);
00979 #endif
00980 if (fabs(alpha_)<1.0e-4) {
00981 problemStatus_=1;
00982 } else {
00983 #ifdef CLP_DEBUG
00984 abort();
00985 #endif
00986 }
00987 } else {
00988 problemStatus_=0;
00989 }
00990 }
00991 } else {
00992 problemStatus_=-3;
00993 }
00994 returnCode=0;
00995 break;
00996 }
00997 }
00998 if (givenDuals) {
00999 memcpy(givenDuals,dj_,(numberRows_+numberColumns_)*sizeof(double));
01000
01001 delete [] candidateList;
01002 }
01003 delete [] dubiousWeights;
01004 return returnCode;
01005 }
01006
01007
01008
01009
01010 int
01011 ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
01012 CoinIndexedVector * columnArray,
01013 CoinIndexedVector * outputArray,
01014 double theta,
01015 double & objectiveChange,
01016 bool fullRecompute)
01017 {
01018
01019 outputArray->clear();
01020
01021
01022 int numberInfeasibilities=0;
01023 int numberRowInfeasibilities=0;
01024
01025
01026 double tolerance = dualTolerance_;
01027
01028 double changeObj=0.0;
01029
01030
01031
01032 if (!fullRecompute) {
01033 int i;
01034 double * reducedCost = djRegion(0);
01035 const double * lower = lowerRegion(0);
01036 const double * upper = upperRegion(0);
01037 const double * cost = costRegion(0);
01038 double * work;
01039 int number;
01040 int * which;
01041 assert(rowArray->packedMode());
01042 work = rowArray->denseVector();
01043 number = rowArray->getNumElements();
01044 which = rowArray->getIndices();
01045 for (i=0;i<number;i++) {
01046 int iSequence = which[i];
01047 double alphaI = work[i];
01048 work[i]=0.0;
01049
01050 Status status = getStatus(iSequence+numberColumns_);
01051
01052 if (status==atUpperBound) {
01053
01054 double value = reducedCost[iSequence]-theta*alphaI;
01055 reducedCost[iSequence]=value;
01056
01057 if (value>tolerance) {
01058
01059 which[numberInfeasibilities++]=iSequence;
01060 double movement = lower[iSequence]-upper[iSequence];
01061 assert (fabs(movement)<1.0e30);
01062 #ifdef CLP_DEBUG
01063 if ((handler_->logLevel()&32))
01064 printf("%d %d, new dj %g, alpha %g, movement %g\n",
01065 0,iSequence,value,alphaI,movement);
01066 #endif
01067 changeObj += movement*cost[iSequence];
01068 outputArray->quickAdd(iSequence,-movement);
01069 }
01070 } else if (status==atLowerBound) {
01071
01072 double value = reducedCost[iSequence]-theta*alphaI;
01073 reducedCost[iSequence]=value;
01074
01075 if (value<-tolerance) {
01076
01077 which[numberInfeasibilities++]=iSequence;
01078 double movement = upper[iSequence] - lower[iSequence];
01079 assert (fabs(movement)<1.0e30);
01080 #ifdef CLP_DEBUG
01081 if ((handler_->logLevel()&32))
01082 printf("%d %d, new dj %g, alpha %g, movement %g\n",
01083 0,iSequence,value,alphaI,movement);
01084 #endif
01085 changeObj += movement*cost[iSequence];
01086 outputArray->quickAdd(iSequence,-movement);
01087 }
01088 }
01089 }
01090 } else {
01091 double * solution = solutionRegion(0);
01092 double * reducedCost = djRegion(0);
01093 const double * lower = lowerRegion(0);
01094 const double * upper = upperRegion(0);
01095 const double * cost = costRegion(0);
01096 int * which;
01097 which = rowArray->getIndices();
01098 int iSequence;
01099 for (iSequence=0;iSequence<numberRows_;iSequence++) {
01100 double value = reducedCost[iSequence];
01101
01102 Status status = getStatus(iSequence+numberColumns_);
01103
01104 if (status==atUpperBound) {
01105 double movement=0.0;
01106
01107 if (value>tolerance) {
01108
01109
01110 which[numberInfeasibilities++]=iSequence;
01111 movement = lower[iSequence]-upper[iSequence];
01112 changeObj += movement*cost[iSequence];
01113 outputArray->quickAdd(iSequence,-movement);
01114 assert (fabs(movement)<1.0e30);
01115 } else if (value>-tolerance) {
01116
01117 FakeBound bound = getFakeBound(iSequence+numberColumns_);
01118 if (bound==ClpSimplexDual::upperFake) {
01119 movement = lower[iSequence]-upper[iSequence];
01120 assert (fabs(movement)<1.0e30);
01121 setStatus(iSequence+numberColumns_,atLowerBound);
01122 solution[iSequence] = lower[iSequence];
01123 changeObj += movement*cost[iSequence];
01124 numberFake_--;
01125 setFakeBound(iSequence+numberColumns_,noFake);
01126 }
01127 }
01128 } else if (status==atLowerBound) {
01129 double movement=0.0;
01130
01131 if (value<-tolerance) {
01132
01133
01134 which[numberInfeasibilities++]=iSequence;
01135 movement = upper[iSequence] - lower[iSequence];
01136 assert (fabs(movement)<1.0e30);
01137 changeObj += movement*cost[iSequence];
01138 outputArray->quickAdd(iSequence,-movement);
01139 } else if (value<tolerance) {
01140
01141 FakeBound bound = getFakeBound(iSequence+numberColumns_);
01142 if (bound==ClpSimplexDual::lowerFake) {
01143 movement = upper[iSequence]-lower[iSequence];
01144 assert (fabs(movement)<1.0e30);
01145 setStatus(iSequence+numberColumns_,atUpperBound);
01146 solution[iSequence] = upper[iSequence];
01147 changeObj += movement*cost[iSequence];
01148 numberFake_--;
01149 setFakeBound(iSequence+numberColumns_,noFake);
01150 }
01151 }
01152 }
01153 }
01154 }
01155
01156
01157 if (!fullRecompute) {
01158 int i;
01159 double * reducedCost = djRegion(1);
01160 const double * lower = lowerRegion(1);
01161 const double * upper = upperRegion(1);
01162 const double * cost = costRegion(1);
01163 double * work;
01164 int number;
01165 int * which;
01166
01167 numberRowInfeasibilities=numberInfeasibilities;
01168 rowArray->setNumElements(numberInfeasibilities);
01169 numberInfeasibilities=0;
01170 work = columnArray->denseVector();
01171 number = columnArray->getNumElements();
01172 which = columnArray->getIndices();
01173
01174 for (i=0;i<number;i++) {
01175 int iSequence = which[i];
01176 double alphaI = work[i];
01177 work[i]=0.0;
01178
01179 Status status = getStatus(iSequence);
01180 if (status==atLowerBound) {
01181 double value = reducedCost[iSequence]-theta*alphaI;
01182 reducedCost[iSequence]=value;
01183 double movement=0.0;
01184
01185 if (value<-tolerance) {
01186
01187 which[numberInfeasibilities++]=iSequence;
01188 movement = upper[iSequence] - lower[iSequence];
01189 assert (fabs(movement)<1.0e30);
01190 #ifdef CLP_DEBUG
01191 if ((handler_->logLevel()&32))
01192 printf("%d %d, new dj %g, alpha %g, movement %g\n",
01193 1,iSequence,value,alphaI,movement);
01194 #endif
01195 changeObj += movement*cost[iSequence];
01196 matrix_->add(this,outputArray,iSequence,movement);
01197 }
01198 } else if (status==atUpperBound) {
01199 double value = reducedCost[iSequence]-theta*alphaI;
01200 reducedCost[iSequence]=value;
01201 double movement=0.0;
01202
01203 if (value>tolerance) {
01204
01205 which[numberInfeasibilities++]=iSequence;
01206 movement = lower[iSequence]-upper[iSequence];
01207 assert (fabs(movement)<1.0e30);
01208 #ifdef CLP_DEBUG
01209 if ((handler_->logLevel()&32))
01210 printf("%d %d, new dj %g, alpha %g, movement %g\n",
01211 1,iSequence,value,alphaI,movement);
01212 #endif
01213 changeObj += movement*cost[iSequence];
01214 matrix_->add(this,outputArray,iSequence,movement);
01215 }
01216 } else if (status==isFree) {
01217 double value = reducedCost[iSequence]-theta*alphaI;
01218 reducedCost[iSequence]=value;
01219 }
01220 }
01221 } else {
01222 double * solution = solutionRegion(1);
01223 double * reducedCost = djRegion(1);
01224 const double * lower = lowerRegion(1);
01225 const double * upper = upperRegion(1);
01226 const double * cost = costRegion(1);
01227 int * which;
01228
01229 numberRowInfeasibilities=numberInfeasibilities;
01230 rowArray->setNumElements(numberInfeasibilities);
01231 numberInfeasibilities=0;
01232 which = columnArray->getIndices();
01233 int iSequence;
01234 for (iSequence=0;iSequence<numberColumns_;iSequence++) {
01235 double value = reducedCost[iSequence];
01236
01237 Status status = getStatus(iSequence);
01238 if (status==atLowerBound) {
01239 double movement=0.0;
01240
01241 if (value<-tolerance) {
01242
01243
01244 which[numberInfeasibilities++]=iSequence;
01245 movement = upper[iSequence] - lower[iSequence];
01246 assert (fabs(movement)<1.0e30);
01247 changeObj += movement*cost[iSequence];
01248 matrix_->add(this,outputArray,iSequence,movement);
01249 } else if (value<tolerance) {
01250
01251 FakeBound bound = getFakeBound(iSequence);
01252 if (bound==ClpSimplexDual::lowerFake) {
01253 movement = upper[iSequence]-lower[iSequence];
01254 assert (fabs(movement)<1.0e30);
01255 setStatus(iSequence,atUpperBound);
01256 solution[iSequence] = upper[iSequence];
01257 changeObj += movement*cost[iSequence];
01258 numberFake_--;
01259 setFakeBound(iSequence,noFake);
01260 }
01261 }
01262 } else if (status==atUpperBound) {
01263 double movement=0.0;
01264
01265 if (value>tolerance) {
01266
01267
01268 which[numberInfeasibilities++]=iSequence;
01269 movement = lower[iSequence]-upper[iSequence];
01270 assert (fabs(movement)<1.0e30);
01271 changeObj += movement*cost[iSequence];
01272 matrix_->add(this,outputArray,iSequence,movement);
01273 } else if (value>-tolerance) {
01274
01275 FakeBound bound = getFakeBound(iSequence);
01276 if (bound==ClpSimplexDual::upperFake) {
01277 movement = lower[iSequence]-upper[iSequence];
01278 assert (fabs(movement)<1.0e30);
01279 setStatus(iSequence,atLowerBound);
01280 solution[iSequence] = lower[iSequence];
01281 changeObj += movement*cost[iSequence];
01282 numberFake_--;
01283 setFakeBound(iSequence,noFake);
01284 }
01285 }
01286 }
01287 }
01288 }
01289 #ifdef CLP_DEBUG
01290 if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0)
01291 printf("%d fake after full update\n",numberFake_);
01292 #endif
01293
01294 columnArray->setNumElements(numberInfeasibilities);
01295 numberInfeasibilities += numberRowInfeasibilities;
01296 if (fullRecompute) {
01297
01298 flipBounds(rowArray,columnArray,0.0);
01299 }
01300 objectiveChange += changeObj;
01301 return numberInfeasibilities;
01302 }
01303 void
01304 ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
01305 CoinIndexedVector * columnArray,
01306 double theta)
01307 {
01308
01309
01310 double tolerance = dualTolerance_;
01311
01312
01313
01314 {
01315 int i;
01316 double * reducedCost = djRegion(0);
01317 double * work;
01318 int number;
01319 int * which;
01320 work = rowArray->denseVector();
01321 number = rowArray->getNumElements();
01322 which = rowArray->getIndices();
01323 for (i=0;i<number;i++) {
01324 int iSequence = which[i];
01325 double alphaI = work[i];
01326 double value = reducedCost[iSequence]-theta*alphaI;
01327 work[i]=0.0;
01328 reducedCost[iSequence]=value;
01329
01330 Status status = getStatus(iSequence+numberColumns_);
01331
01332 if (status==atUpperBound) {
01333
01334 if (value>tolerance)
01335 reducedCost[iSequence]=0.0;
01336 } else if (status==atLowerBound) {
01337
01338 if (value<-tolerance) {
01339 reducedCost[iSequence]=0.0;
01340 }
01341 }
01342 }
01343 }
01344 rowArray->setNumElements(0);
01345
01346
01347 {
01348 int i;
01349 double * reducedCost = djRegion(1);
01350 double * work;
01351 int number;
01352 int * which;
01353 work = columnArray->denseVector();
01354 number = columnArray->getNumElements();
01355 which = columnArray->getIndices();
01356
01357 for (i=0;i<number;i++) {
01358 int iSequence = which[i];
01359 double alphaI = work[i];
01360 double value = reducedCost[iSequence]-theta*alphaI;
01361 work[i]=0.0;
01362 reducedCost[iSequence]=value;
01363
01364 Status status = getStatus(iSequence);
01365 if (status==atLowerBound) {
01366 if (value<-tolerance)
01367 reducedCost[iSequence]=0.0;
01368 } else if (status==atUpperBound) {
01369 if (value>tolerance)
01370 reducedCost[iSequence]=0.0;
01371 }
01372 }
01373 }
01374 columnArray->setNumElements(0);
01375 }
01376
01377
01378
01379
01380
01381
01382 void
01383 ClpSimplexDual::dualRow(int alreadyChosen)
01384 {
01385
01386 int chosenRow=-1;
01387 if (alreadyChosen<0) {
01388
01389 int nextFree = nextSuperBasic();
01390
01391 if (nextFree>=0) {
01392
01393 unpack(rowArray_[1],nextFree);
01394 factorization_->updateColumn(rowArray_[2],rowArray_[1]);
01395
01396 double * work=rowArray_[1]->denseVector();
01397 int number=rowArray_[1]->getNumElements();
01398 int * which=rowArray_[1]->getIndices();
01399 double bestFeasibleAlpha=0.0;
01400 int bestFeasibleRow=-1;
01401 double bestInfeasibleAlpha=0.0;
01402 int bestInfeasibleRow=-1;
01403 int i;
01404
01405 for (i=0;i<number;i++) {
01406 int iRow = which[i];
01407 double alpha = fabs(work[iRow]);
01408 if (alpha>1.0e-3) {
01409 int iSequence=pivotVariable_[iRow];
01410 double value = solution_[iSequence];
01411 double lower = lower_[iSequence];
01412 double upper = upper_[iSequence];
01413 double infeasibility=0.0;
01414 if (value>upper)
01415 infeasibility = value-upper;
01416 else if (value<lower)
01417 infeasibility = lower-value;
01418 if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
01419 if (!flagged(iSequence)) {
01420 bestInfeasibleAlpha = infeasibility*alpha;
01421 bestInfeasibleRow=iRow;
01422 }
01423 }
01424 if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
01425 bestFeasibleAlpha = alpha;
01426 bestFeasibleRow=iRow;
01427 }
01428 }
01429 }
01430 if (bestInfeasibleRow>=0)
01431 chosenRow = bestInfeasibleRow;
01432 else if (bestFeasibleAlpha>1.0e-2)
01433 chosenRow = bestFeasibleRow;
01434 if (chosenRow>=0)
01435 pivotRow_=chosenRow;
01436 rowArray_[1]->clear();
01437 }
01438 } else {
01439
01440 chosenRow=alreadyChosen;
01441 pivotRow_=chosenRow;
01442 }
01443 if (chosenRow<0)
01444 pivotRow_=dualRowPivot_->pivotRow();
01445
01446 if (pivotRow_>=0) {
01447 sequenceOut_ = pivotVariable_[pivotRow_];
01448 valueOut_ = solution_[sequenceOut_];
01449 lowerOut_ = lower_[sequenceOut_];
01450 upperOut_ = upper_[sequenceOut_];
01451 if (alreadyChosen<0) {
01452
01453
01454 if (valueOut_>upperOut_) {
01455 directionOut_ = -1;
01456 dualOut_ = valueOut_ - upperOut_;
01457 } else if (valueOut_<lowerOut_) {
01458 directionOut_ = 1;
01459 dualOut_ = lowerOut_ - valueOut_;
01460 } else {
01461
01462 if (valueOut_-lowerOut_<upperOut_-valueOut_) {
01463 directionOut_ = 1;
01464 dualOut_ = lowerOut_ - valueOut_;
01465 } else {
01466 directionOut_ = -1;
01467 dualOut_ = valueOut_ - upperOut_;
01468 }
01469 }
01470 #ifdef CLP_DEBUG
01471 assert(dualOut_>=0.0);
01472 #endif
01473 } else {
01474
01475
01476
01477 dualOut_ = 1.0e-6;
01478 if (dj_[sequenceOut_]>0.0) {
01479
01480 directionOut_ = 1;
01481 } else {
01482 directionOut_ = -1;
01483 }
01484 }
01485 }
01486 return ;
01487 }
01488
01489
01490
01491
01492
01493
01494
01495 int
01496 ClpSimplexDual::changeBounds(bool initialize,
01497 CoinIndexedVector * outputArray,
01498 double & changeCost)
01499 {
01500 numberFake_ = 0;
01501 if (!initialize) {
01502 int numberInfeasibilities;
01503 double newBound;
01504 newBound = 5.0*dualBound_;
01505 numberInfeasibilities=0;
01506 changeCost=0.0;
01507
01508 createRim(3);
01509 int iSequence;
01510
01511 for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
01512 double lowerValue=lower_[iSequence];
01513 double upperValue=upper_[iSequence];
01514 double value=solution_[iSequence];
01515 setFakeBound(iSequence,ClpSimplexDual::noFake);
01516 switch(getStatus(iSequence)) {
01517
01518 case basic:
01519 case ClpSimplex::isFixed:
01520 break;
01521 case isFree:
01522 case superBasic:
01523 break;
01524 case atUpperBound:
01525 if (fabs(value-upperValue)>primalTolerance_)
01526 numberInfeasibilities++;
01527 break;
01528 case atLowerBound:
01529 if (fabs(value-lowerValue)>primalTolerance_)
01530 numberInfeasibilities++;
01531 break;
01532 }
01533 }
01534
01535 if (numberInfeasibilities) {
01536 handler_->message(CLP_DUAL_CHECKB,messages_)
01537 <<newBound
01538 <<CoinMessageEol;
01539 int iSequence;
01540 for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
01541 double lowerValue=lower_[iSequence];
01542 double upperValue=upper_[iSequence];
01543 double newLowerValue;
01544 double newUpperValue;
01545 Status status = getStatus(iSequence);
01546 if (status==atUpperBound||
01547 status==atLowerBound) {
01548 double value = solution_[iSequence];
01549 if (value-lowerValue<=upperValue-value) {
01550 newLowerValue = max(lowerValue,value-0.666667*newBound);
01551 newUpperValue = min(upperValue,newLowerValue+newBound);
01552 } else {
01553 newUpperValue = min(upperValue,value+0.666667*newBound);
01554 newLowerValue = max(lowerValue,newUpperValue-newBound);
01555 }
01556 lower_[iSequence]=newLowerValue;
01557 upper_[iSequence]=newUpperValue;
01558 if (newLowerValue > lowerValue) {
01559 if (newUpperValue < upperValue) {
01560 setFakeBound(iSequence,ClpSimplexDual::bothFake);
01561 numberFake_++;
01562 } else {
01563 setFakeBound(iSequence,ClpSimplexDual::lowerFake);
01564 numberFake_++;
01565 }
01566 } else {
01567 if (newUpperValue < upperValue) {
01568 setFakeBound(iSequence,ClpSimplexDual::upperFake);
01569 numberFake_++;
01570 }
01571 }
01572 if (status==atUpperBound)
01573 solution_[iSequence] = newUpperValue;
01574 else
01575 solution_[iSequence] = newLowerValue;
01576 double movement = solution_[iSequence] - value;
01577 if (movement&&outputArray) {
01578 if (iSequence>=numberColumns_) {
01579 outputArray->quickAdd(iSequence,-movement);
01580 changeCost += movement*cost_[iSequence];
01581 } else {
01582 matrix_->add(this,outputArray,iSequence,movement);
01583 changeCost += movement*cost_[iSequence];
01584 }
01585 }
01586 }
01587 }
01588 dualBound_ = newBound;
01589 } else {
01590 numberInfeasibilities=-1;
01591 }
01592 return numberInfeasibilities;
01593 } else {
01594 int iSequence;
01595
01596 for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
01597 Status status = getStatus(iSequence);
01598 if (status==atUpperBound||
01599 status==atLowerBound) {
01600 double lowerValue=lower_[iSequence];
01601 double upperValue=upper_[iSequence];
01602 double value = solution_[iSequence];
01603 if (lowerValue>-largeValue_||upperValue<largeValue_) {
01604 if (lowerValue-value>-0.5*dualBound_||
01605 upperValue-value<0.5*dualBound_) {
01606 if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
01607 if (upperValue > lowerValue + dualBound_) {
01608 upper_[iSequence]=lowerValue+dualBound_;
01609 setFakeBound(iSequence,ClpSimplexDual::upperFake);
01610 numberFake_++;
01611 }
01612 } else {
01613 if (lowerValue < upperValue - dualBound_) {
01614 lower_[iSequence]=upperValue-dualBound_;
01615 setFakeBound(iSequence,ClpSimplexDual::lowerFake);
01616 numberFake_++;
01617 }
01618 }
01619 } else {
01620 lower_[iSequence]=-0.5*dualBound_;
01621 upper_[iSequence]= 0.5*dualBound_;
01622 setFakeBound(iSequence,ClpSimplexDual::bothFake);
01623 numberFake_++;
01624 }
01625 } else {
01626
01627
01628 assert(!("should not be here"));
01629 lower_[iSequence]=-0.5*dualBound_;
01630 upper_[iSequence]= 0.5*dualBound_;
01631 setFakeBound(iSequence,ClpSimplexDual::bothFake);
01632 numberFake_++;
01633 setStatus(iSequence,atUpperBound);
01634 solution_[iSequence]=0.5*dualBound_;
01635 }
01636 }
01637 }
01638
01639 return 1;
01640 }
01641 }
01642
01643
01644
01645
01646
01647
01648
01649
01650 void
01651 ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
01652 CoinIndexedVector * columnArray,
01653 CoinIndexedVector * spareArray,
01654 CoinIndexedVector * spareArray2,
01655 double acceptablePivot,
01656 CoinBigIndex * dubiousWeights)
01657 {
01658 double * work;
01659 int number;
01660 int * which;
01661 double * reducedCost;
01662 int iSection;
01663
01664 sequenceIn_=-1;
01665 int numberPossiblySwapped=0;
01666 int numberRemaining=0;
01667
01668 double totalThru=0.0;
01669
01670 double bestEverPivot=acceptablePivot;
01671 int lastSequence = -1;
01672 double lastPivot=0.0;
01673 double upperTheta;
01674 double newTolerance = dualTolerance_;
01675
01676
01677 bool thisIncrease=false;
01678
01679 bool modifyCosts=false;
01680
01681 double increaseInObjective=0.0;
01682
01683
01684
01685 int iFlip = 0;
01686
01687 int interesting[2];
01688
01689 int swapped[2];
01690
01691 int marker[2][2];
01692
01693 double * array[2], * spare, * spare2;
01694
01695 int * indices[2], * index, * index2;
01696 spareArray->clear();
01697 spareArray2->clear();
01698 array[0] = spareArray->denseVector();
01699 indices[0] = spareArray->getIndices();
01700 spare = array[0];
01701 index = indices[0];
01702 array[1] = spareArray2->denseVector();
01703 indices[1] = spareArray2->getIndices();
01704 int i;
01705 const double * lower;
01706 const double * upper;
01707
01708
01709 for (i=0;i<2;i++) {
01710 interesting[i]=0;
01711 swapped[i]=numberColumns_;
01712 marker[i][0]=0;
01713 marker[i][1]=numberColumns_;
01714 }
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736 double tentativeTheta = 1.0e22;
01737 upperTheta = 1.0e31;
01738 double freePivot = acceptablePivot;
01739
01740 for (iSection=0;iSection<2;iSection++) {
01741
01742 int addSequence;
01743
01744 if (!iSection) {
01745 lower = rowLowerWork_;
01746 upper = rowUpperWork_;
01747 work = rowArray->denseVector();
01748 number = rowArray->getNumElements();
01749 which = rowArray->getIndices();
01750 reducedCost = rowReducedCost_;
01751 addSequence = numberColumns_;
01752 } else {
01753 lower = columnLowerWork_;
01754 upper = columnUpperWork_;
01755 work = columnArray->denseVector();
01756 number = columnArray->getNumElements();
01757 which = columnArray->getIndices();
01758 reducedCost = reducedCostWork_;
01759 addSequence = 0;
01760 }
01761
01762 for (i=0;i<number;i++) {
01763 int iSequence = which[i];
01764 double alpha;
01765 double oldValue;
01766 double value;
01767 bool keep;
01768
01769 switch(getStatus(iSequence+addSequence)) {
01770
01771 case basic:
01772 case ClpSimplex::isFixed:
01773 break;
01774 case isFree:
01775 case superBasic:
01776 alpha = work[i];
01777 oldValue = reducedCost[iSequence];
01778 if (oldValue>dualTolerance_) {
01779 keep = true;
01780 } else if (oldValue<-dualTolerance_) {
01781 keep = true;
01782 } else {
01783 if (fabs(alpha)>10.0*acceptablePivot)
01784 keep = true;
01785 else
01786 keep=false;
01787 }
01788 if (keep) {
01789
01790 if (fabs(alpha)>freePivot) {
01791 freePivot=fabs(alpha);
01792 sequenceIn_ = iSequence + addSequence;
01793 theta_=oldValue/alpha;
01794 alpha_=alpha;
01795 }
01796 }
01797 break;
01798 case atUpperBound:
01799 alpha = work[i];
01800 oldValue = reducedCost[iSequence];
01801 value = oldValue-tentativeTheta*alpha;
01802 assert (oldValue<=dualTolerance_*1.0001);
01803 if (value>newTolerance) {
01804 value = oldValue-upperTheta*alpha;
01805 if (value>newTolerance && -alpha>=acceptablePivot)
01806 upperTheta = (oldValue-newTolerance)/alpha;
01807
01808 spare[numberRemaining]=alpha;
01809 index[numberRemaining++]=iSequence+addSequence;
01810 }
01811 break;
01812 case atLowerBound:
01813 alpha = work[i];
01814 oldValue = reducedCost[iSequence];
01815 value = oldValue-tentativeTheta*alpha;
01816 assert (oldValue>=-dualTolerance_*1.0001);
01817 if (value<-newTolerance) {
01818 value = oldValue-upperTheta*alpha;
01819 if (value<-newTolerance && alpha>=acceptablePivot)
01820 upperTheta = (oldValue+newTolerance)/alpha;
01821
01822 spare[numberRemaining]=alpha;
01823 index[numberRemaining++]=iSequence+addSequence;
01824 }
01825 break;
01826 }
01827 }
01828 }
01829 interesting[0]=numberRemaining;
01830 marker[0][0] = numberRemaining;
01831
01832 if (!numberRemaining&&sequenceIn_<0)
01833 return ;
01834
01835 if (sequenceIn_>=0) {
01836
01837 } else {
01838
01839 theta_=1.0e50;
01840
01841 tentativeTheta = max(10.0*upperTheta,1.0e-7);
01842
01843
01844
01845 while (tentativeTheta < 1.0e22) {
01846 double thruThis = 0.0;
01847
01848 double bestPivot=acceptablePivot;
01849 int bestSequence=-1;
01850
01851 numberPossiblySwapped = numberColumns_;
01852 numberRemaining = 0;
01853
01854 upperTheta = 1.0e50;
01855
01856 spare = array[iFlip];
01857 index = indices[iFlip];
01858 spare2 = array[1-iFlip];
01859 index2 = indices[1-iFlip];
01860
01861
01862
01863
01864
01865 #define TRYBIAS 1
01866
01867
01868 double increaseInThis=0.0;
01869
01870 for (i=0;i<interesting[iFlip];i++) {
01871 int iSequence = index[i];
01872 double alpha = spare[i];
01873 double oldValue = dj_[iSequence];
01874 double value = oldValue-tentativeTheta*alpha;
01875
01876 if (alpha < 0.0) {
01877
01878 if (value>newTolerance) {
01879 double range = upper_[iSequence] - lower_[iSequence];
01880 thruThis -= range*alpha;
01881 #if TRYBIAS==1
01882 if (oldValue>0.0)
01883 increaseInThis -= oldValue*range;
01884 #elif TRYBIAS==2
01885 increaseInThis -= oldValue*range;
01886 #else
01887 increaseInThis -= (oldValue+dualTolerance_)*range;
01888 #endif
01889
01890 spare2[--numberPossiblySwapped]=alpha;
01891 index2[numberPossiblySwapped]=iSequence;
01892 if (fabs(alpha)>bestPivot) {
01893 bestPivot=fabs(alpha);
01894 bestSequence=numberPossiblySwapped;
01895 }
01896 } else {
01897 value = oldValue-upperTheta*alpha;
01898 if (value>newTolerance && -alpha>=acceptablePivot)
01899 upperTheta = (oldValue-newTolerance)/alpha;
01900 spare2[numberRemaining]=alpha;
01901 index2[numberRemaining++]=iSequence;
01902 }
01903 } else {
01904
01905 if (value<-newTolerance) {
01906 double range = upper_[iSequence] - lower_[iSequence];
01907 thruThis += range*alpha;
01908
01909 #if TRYBIAS==1
01910 if (oldValue<0.0)
01911 increaseInThis += oldValue*range;
01912 #elif TRYBIAS==2
01913 increaseInThis += oldValue*range;
01914 #else
01915 increaseInThis += (oldValue-dualTolerance_)*range;
01916 #endif
01917
01918 spare2[--numberPossiblySwapped]=alpha;
01919 index2[numberPossiblySwapped]=iSequence;
01920 if (fabs(alpha)>bestPivot) {
01921 bestPivot=fabs(alpha);
01922 bestSequence=numberPossiblySwapped;
01923 }
01924 } else {
01925 value = oldValue-upperTheta*alpha;
01926 if (value<-newTolerance && alpha>=acceptablePivot)
01927 upperTheta = (oldValue+newTolerance)/alpha;
01928 spare2[numberRemaining]=alpha;
01929 index2[numberRemaining++]=iSequence;
01930 }
01931 }
01932 }
01933 swapped[1-iFlip]=numberPossiblySwapped;
01934 interesting[1-iFlip]=numberRemaining;
01935 marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
01936 marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
01937
01938 if (totalThru+thruThis>=fabs(dualOut_)||
01939 increaseInObjective+increaseInThis<0.0) {
01940
01941
01942 numberRemaining=0;
01943 for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
01944 spare[numberRemaining]=spare2[i];
01945 index[numberRemaining++]=index2[i];
01946 }
01947 interesting[iFlip]=numberRemaining;
01948 int iTry;
01949 #define MAXTRY 100
01950
01951 for (iTry=0;iTry<MAXTRY;iTry++) {
01952
01953 upperTheta=1.0e50;
01954 numberPossiblySwapped = numberColumns_;
01955 numberRemaining = 0;
01956
01957 increaseInThis=0.0;
01958
01959 thruThis=0.0;
01960
01961 spare = array[iFlip];
01962 index = indices[iFlip];
01963 spare2 = array[1-iFlip];
01964 index2 = indices[1-iFlip];
01965
01966 for (i=0;i<interesting[iFlip];i++) {
01967 int iSequence=index[i];
01968 double alpha=spare[i];
01969 double oldValue = dj_[iSequence];
01970 double value = oldValue-upperTheta*alpha;
01971
01972 if (alpha < 0.0) {
01973
01974 if (value>newTolerance) {
01975 if (-alpha>=acceptablePivot) {
01976 upperTheta = (oldValue-newTolerance)/alpha;
01977
01978 value = oldValue-upperTheta*alpha;
01979 if (value<0)
01980 upperTheta *= 1.0 +1.0e-11;
01981 }
01982 }
01983 } else {
01984
01985 if (value<-newTolerance) {
01986 if (alpha>=acceptablePivot) {
01987 upperTheta = (oldValue+newTolerance)/alpha;
01988
01989 value = oldValue-upperTheta*alpha;
01990 if (value>0)
01991 upperTheta *= 1.0 +1.0e-11;
01992 }
01993 }
01994 }
01995 }
01996 bestPivot=acceptablePivot;
01997 sequenceIn_=-1;
01998 double bestWeight=COIN_DBL_MAX;
01999 double largestPivot=acceptablePivot;
02000
02001
02002 for (i=0;i<interesting[iFlip];i++) {
02003 int iSequence=index[i];
02004 double alpha=spare[i];
02005 double value = dj_[iSequence]-upperTheta*alpha;
02006 double badDj=0.0;
02007
02008 bool addToSwapped=false;
02009
02010 if (alpha < 0.0) {
02011
02012 if (value>=0.0) {
02013 addToSwapped=true;
02014 #if TRYBIAS==1
02015 badDj = -max(dj_[iSequence],0.0);
02016 #elif TRYBIAS==2
02017 badDj = -dj_[iSequence];
02018 #else
02019 badDj = -dj_[iSequence]-dualTolerance_;
02020 #endif
02021 }
02022 } else {
02023
02024 if (value<=0.0) {
02025 addToSwapped=true;
02026 #if TRYBIAS==1
02027 badDj = min(dj_[iSequence],0.0);
02028 #elif TRYBIAS==2
02029 badDj = dj_[iSequence];
02030 #else
02031 badDj = dj_[iSequence]-dualTolerance_;
02032 #endif
02033 }
02034 }
02035 if (!addToSwapped) {
02036
02037 spare2[numberRemaining]=alpha;
02038 index2[numberRemaining++]=iSequence;
02039 } else {
02040
02041 spare2[--numberPossiblySwapped]=alpha;
02042 index2[numberPossiblySwapped]=iSequence;
02043
02044 bool take=false;
02045 double absAlpha = fabs(alpha);
02046
02047 double weight;
02048 if (dubiousWeights)
02049 weight=dubiousWeights[iSequence];
02050 else
02051 weight=1.0;
02052 weight += CoinDrand48()*1.0e-2;
02053 if (absAlpha>2.0*bestPivot) {
02054 take=true;
02055 } else if (absAlpha>0.5*largestPivot) {
02056
02057 if (weight*bestPivot<bestWeight*absAlpha)
02058 take=true;
02059 }
02060 if (take) {
02061 sequenceIn_ = numberPossiblySwapped;
02062 bestPivot = absAlpha;
02063 theta_ = dj_[iSequence]/alpha;
02064 largestPivot = max(largestPivot,bestPivot);
02065 bestWeight = weight;
02066
02067
02068 } else {
02069
02070
02071 }
02072 double range = upper[iSequence] - lower[iSequence];
02073 thruThis += range*fabs(alpha);
02074 increaseInThis += badDj*range;
02075 }
02076 }
02077 swapped[1-iFlip]=numberPossiblySwapped;
02078 interesting[1-iFlip]=numberRemaining;
02079 marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
02080 marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
02081
02082 double increase = (fabs(dualOut_)-totalThru)*theta_;
02083 increase += increaseInObjective;
02084 if (theta_<0.0)
02085 thruThis += fabs(dualOut_);
02086 if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
02087
02088
02089
02090 bestPivot=0.0;
02091 } else {
02092
02093 totalThru += thruThis;
02094 increaseInObjective += increaseInThis;
02095 }
02096 if (bestPivot<0.1*bestEverPivot&&
02097 bestEverPivot>1.0e-6&&
02098 (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
02099
02100 sequenceIn_=lastSequence;
02101
02102 iFlip = 1-iFlip;
02103 break;
02104 } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
02105 if (lastPivot>acceptablePivot) {
02106
02107 sequenceIn_=lastSequence;
02108
02109 iFlip = 1-iFlip;
02110 } else {
02111
02112 }
02113 break;
02114 } else if (totalThru>=fabs(dualOut_)) {
02115 modifyCosts=true;
02116 break;
02117 } else {
02118 lastSequence=sequenceIn_;
02119 if (bestPivot>bestEverPivot)
02120 bestEverPivot=bestPivot;
02121 iFlip = 1 -iFlip;
02122 modifyCosts=true;
02123 }
02124 }
02125 if (iTry==MAXTRY)
02126 iFlip = 1-iFlip;
02127 break;
02128 } else {
02129
02130 if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
02131 bestEverPivot=bestPivot;
02132 lastSequence=bestSequence;
02133 } else {
02134
02135 memcpy(array[1-iFlip]+swapped[iFlip],
02136 array[iFlip]+swapped[iFlip],
02137 (numberColumns_-swapped[iFlip])*sizeof(double));
02138 memcpy(indices[1-iFlip]+swapped[iFlip],
02139 indices[iFlip]+swapped[iFlip],
02140 (numberColumns_-swapped[iFlip])*sizeof(int));
02141 marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
02142 swapped[1-iFlip]=swapped[iFlip];
02143 }
02144 increaseInObjective += increaseInThis;
02145 iFlip = 1 - iFlip;
02146 tentativeTheta = 2.0*upperTheta;
02147 totalThru += thruThis;
02148 }
02149 }
02150
02151
02152 if (sequenceIn_<0&&lastSequence>=0) {
02153
02154 sequenceIn_=lastSequence;
02155
02156 iFlip = 1-iFlip;
02157 }
02158
02159 #define MINIMUMTHETA 1.0e-12
02160
02161
02162 double minimumTheta;
02163 if (upperOut_>lowerOut_)
02164 minimumTheta=MINIMUMTHETA;
02165 else
02166 minimumTheta=0.0;
02167 if (sequenceIn_>=0) {
02168
02169
02170 iFlip = 1 -iFlip;
02171 spare = array[iFlip];
02172 index = indices[iFlip];
02173 double oldValue;
02174 alpha_ = spare[sequenceIn_];
02175 sequenceIn_ = indices[iFlip][sequenceIn_];
02176 oldValue = dj_[sequenceIn_];
02177 theta_ = oldValue/alpha_;
02178 if (theta_<minimumTheta) {
02179
02180 if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
02181 theta_=minimumTheta;
02182 } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
02183 theta_=minimumTheta;
02184 thisIncrease=true;
02185 } else {
02186 theta_=(oldValue+newTolerance)/alpha_;
02187 assert(theta_>=0.0);
02188 thisIncrease=true;
02189 }
02190 }
02191
02192
02193 if (modifyCosts) {
02194 int i;
02195 for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
02196 int iSequence=index[i];
02197 double alpha=spare[i];
02198 double value = dj_[iSequence]-theta_*alpha;
02199
02200
02201
02202 if (alpha < 0.0) {
02203
02204 if (value>dualTolerance_) {
02205 thisIncrease=true;
02206 #define MODIFYCOST 2
02207 #if MODIFYCOST
02208
02209 double modification = alpha*theta_-dj_[iSequence]
02210 +newTolerance;
02211
02212
02213 dj_[iSequence] += modification;
02214 cost_[iSequence] += modification;
02215
02216 #endif
02217 }
02218 } else {
02219
02220 if (-value>dualTolerance_) {
02221 thisIncrease=true;
02222 #if MODIFYCOST
02223
02224 double modification = alpha*theta_-dj_[iSequence]
02225 -newTolerance;
02226
02227
02228 dj_[iSequence] += modification;
02229 cost_[iSequence] += modification;
02230
02231 #endif
02232 }
02233 }
02234 }
02235 }
02236 }
02237 }
02238
02239 if (sequenceIn_>=0) {
02240 lowerIn_ = lower_[sequenceIn_];
02241 upperIn_ = upper_[sequenceIn_];
02242 valueIn_ = solution_[sequenceIn_];
02243 dualIn_ = dj_[sequenceIn_];
02244
02245 if (numberTimesOptimal_) {
02246
02247
02248 }
02249 #if MODIFYCOST>1
02250
02251
02252 double modification = theta_*alpha_-dualIn_;
02253 dualIn_ += modification;
02254 dj_[sequenceIn_]=dualIn_;
02255 cost_[sequenceIn_] += modification;
02256
02257
02258
02259 #ifdef CLP_DEBUG
02260 if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
02261 printf("exact %d new cost %g, change %g\n",sequenceIn_,
02262 cost_[sequenceIn_],modification);
02263 #endif
02264 #endif
02265
02266 if (alpha_<0.0) {
02267
02268 directionIn_=-1;
02269 upperIn_=valueIn_;
02270 } else {
02271
02272 directionIn_=1;
02273 lowerIn_=valueIn_;
02274 }
02275 }
02276
02277
02278
02279
02280
02281 for (i=0;i<2;i++) {
02282 memset(array[i],0,marker[i][0]*sizeof(double));
02283 memset(array[i]+marker[i][1],0,
02284 (numberColumns_-marker[i][1])*sizeof(double));
02285 }
02286 }
02287
02288
02289 int
02290 ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
02291 CoinIndexedVector * spare,
02292 double changeCost)
02293 {
02294 int status=2;
02295 factorization_->updateColumn(spare,ray);
02296
02297 int i;
02298 int number=ray->getNumElements();
02299 int * index = ray->getIndices();
02300 double * array = ray->denseVector();
02301 for (i=0;i<number;i++) {
02302 int iRow=index[i];
02303 int iPivot=pivotVariable_[iRow];
02304 changeCost -= cost(iPivot)*array[iRow];
02305 }
02306 double way;
02307 if (changeCost>0.0) {
02308
02309 way=1.0;
02310 } else if (changeCost<0.0) {
02311
02312 way=-1.0;
02313 } else {
02314 #ifdef CLP_DEBUG
02315 printf("can't decide on up or down\n");
02316 #endif
02317 way=0.0;
02318 status=-3;
02319 }
02320 double movement=1.0e10*way;
02321 double zeroTolerance = 1.0e-14*dualBound_;
02322 for (i=0;i<number;i++) {
02323 int iRow=index[i];
02324 int iPivot=pivotVariable_[iRow];
02325 double arrayValue = array[iRow];
02326 if (fabs(arrayValue)<zeroTolerance)
02327 arrayValue=0.0;
02328 double newValue=solution(iPivot)+movement*arrayValue;
02329 if (newValue>upper(iPivot)+primalTolerance_||
02330 newValue<lower(iPivot)-primalTolerance_)
02331 status=-3;
02332 }
02333 if (status==2) {
02334
02335 delete [] ray_;
02336 ray_ = new double [numberColumns_];
02337 ClpFillN(ray_,numberColumns_,0.0);
02338 for (i=0;i<number;i++) {
02339 int iRow=index[i];
02340 int iPivot=pivotVariable_[iRow];
02341 double arrayValue = array[iRow];
02342 if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
02343 ray_[iPivot] = way* array[iRow];
02344 }
02345 }
02346 ray->clear();
02347 return status;
02348 }
02349
02350 void
02351 ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
02352 ClpSimplexProgress * progress,
02353 double * givenDuals)
02354 {
02355 bool normalType=true;
02356 if (type==2) {
02357
02358 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
02359 memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
02360 numberRows_*sizeof(double));
02361 memcpy(columnActivityWork_,savedSolution_ ,
02362 numberColumns_*sizeof(double));
02363 forceFactorization_=1;
02364 changeMade_++;
02365 }
02366 int tentativeStatus = problemStatus_;
02367 double changeCost;
02368 bool unflagVariables = true;
02369 if (problemStatus_>-3||factorization_->pivots()) {
02370
02371
02372
02373
02374 dualRowPivot_->saveWeights(this,1);
02375 if (type) {
02376
02377 if (internalFactorize(1)) {
02378
02379 unflagVariables = false;
02380 assert (type==1);
02381 changeMade_++;
02382 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
02383 memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
02384 numberRows_*sizeof(double));
02385 memcpy(columnActivityWork_,savedSolution_ ,
02386 numberColumns_*sizeof(double));
02387
02388 double dummyChangeCost=0.0;
02389 changeBounds(true,rowArray_[2],dummyChangeCost);
02390
02391 int i;
02392 for (i=0;i<4;i++)
02393 rowArray_[i]->clear();
02394
02395 char x = isColumn(sequenceOut_) ? 'C' :'R';
02396 handler_->message(CLP_SIMPLEX_FLAG,messages_)
02397 <<x<<sequenceWithin(sequenceOut_)
02398 <<CoinMessageEol;
02399 setFlagged(sequenceOut_);
02400 progress_->clearBadTimes();
02401
02402 forceFactorization_=1;
02403 type = 2;
02404
02405 if (internalFactorize(1)) {
02406 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
02407 memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
02408 numberRows_*sizeof(double));
02409 memcpy(columnActivityWork_,savedSolution_ ,
02410 numberColumns_*sizeof(double));
02411
02412 #ifndef NDEBUG
02413 int returnCode = internalFactorize(1);
02414 assert (returnCode==0);
02415 #else
02416 internalFactorize(1);
02417 #endif
02418 }
02419 }
02420 }
02421 if (problemStatus_!=-4||factorization_->pivots()>10)
02422 problemStatus_=-3;
02423 }
02424
02425
02426 gutsOfSolution(givenDuals,NULL);
02427
02428 if (progress->lastIterationNumber(0)==numberIterations_) {
02429 if (dualRowPivot_->looksOptimal()) {
02430 numberPrimalInfeasibilities_ = 0;
02431 sumPrimalInfeasibilities_ = 0.0;
02432 }
02433 }
02434
02435 int loop;
02436 if (!givenDuals&&type!=2)
02437 loop = progress->looping();
02438 else
02439 loop=-1;
02440 int situationChanged=0;
02441 if (loop>=0) {
02442 problemStatus_ = loop;
02443 if (!problemStatus_) {
02444
02445 numberPrimalInfeasibilities_ = 0;
02446 sumPrimalInfeasibilities_ = 0.0;
02447 } else {
02448 problemStatus_ = 10;
02449 }
02450 return;
02451 } else if (loop<-1) {
02452
02453 gutsOfSolution(NULL,NULL);
02454 situationChanged=1;
02455 }
02456
02457 if((progressFlag_&2)!=0) {
02458 situationChanged=2;
02459 }
02460 progressFlag_ = 0;
02461 handler_->message(CLP_SIMPLEX_STATUS,messages_)
02462 <<numberIterations_<<objectiveValue();
02463 handler_->printing(sumPrimalInfeasibilities_>0.0)
02464 <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
02465 handler_->printing(sumDualInfeasibilities_>0.0)
02466 <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
02467 handler_->printing(numberDualInfeasibilitiesWithoutFree_
02468 <numberDualInfeasibilities_)
02469 <<numberDualInfeasibilitiesWithoutFree_;
02470 handler_->message()<<CoinMessageEol;
02471
02472
02473
02474 while (problemStatus_<=-3) {
02475 int cleanDuals=0;
02476 if (situationChanged!=0)
02477 cleanDuals=1;
02478 int numberChangedBounds=0;
02479 int doOriginalTolerance=0;
02480 if ( lastCleaned==numberIterations_)
02481 doOriginalTolerance=1;
02482
02483
02484 if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
02485 sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
02486
02487 numberDualInfeasibilities_ = 0;
02488 sumDualInfeasibilities_ = 0.0;
02489 numberPrimalInfeasibilities_ = 0;
02490 sumPrimalInfeasibilities_ = 0.0;
02491 }
02492
02493 if (dualFeasible()||problemStatus_==-4) {
02494 progress->modifyObjective(objectiveValue_
02495 -sumDualInfeasibilities_*dualBound_);
02496 if (primalFeasible()) {
02497 normalType=false;
02498
02499 handler_->message(CLP_DUAL_BOUNDS,messages_)
02500 <<dualBound_
02501 <<CoinMessageEol;
02502
02503 ClpDisjointCopyN(columnActivityWork_,numberColumns_,
02504 columnArray_[0]->denseVector());
02505 ClpDisjointCopyN(rowActivityWork_,numberRows_,
02506 rowArray_[2]->denseVector());
02507 numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
02508 if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
02509
02510 if (perturbation_==101) {
02511 perturbation_=102;
02512 cleanDuals=1;
02513 createRim(4);
02514
02515 changeBounds(true,NULL,changeCost);
02516 }
02517 if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
02518 doOriginalTolerance=2;
02519 numberTimesOptimal_++;
02520 changeMade_++;
02521 if (numberTimesOptimal_==1) {
02522 dualTolerance_ = dblParam_[ClpDualTolerance];
02523
02524 factorization_->zeroTolerance(1.0e-15);
02525 } else {
02526 dualTolerance_ = dblParam_[ClpDualTolerance];
02527 dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
02528 }
02529 cleanDuals=2;
02530 } else {
02531 problemStatus_=0;
02532 if (lastCleaned<numberIterations_) {
02533 handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
02534 <<CoinMessageEol;
02535 }
02536 }
02537 } else {
02538 cleanDuals=1;
02539 if (doOriginalTolerance==1) {
02540
02541
02542 int iSequence;
02543 int iChosen=-1;
02544 double largest = 100.0*primalTolerance_;
02545 for (iSequence=0;iSequence<numberRows_+numberColumns_;
02546 iSequence++) {
02547 double djValue = dj_[iSequence];
02548 double originalLo = originalLower(iSequence);
02549 double originalUp = originalUpper(iSequence);
02550 if (fabs(djValue)>fabs(largest)) {
02551 if (getStatus(iSequence)!=basic) {
02552 if (djValue>0&&originalLo<-1.0e20) {
02553 if (djValue>fabs(largest)) {
02554 largest=djValue;
02555 iChosen=iSequence;
02556 }
02557 } else if (djValue<0&&originalUp>1.0e20) {
02558 if (-djValue>fabs(largest)) {
02559 largest=djValue;
02560 iChosen=iSequence;
02561 }
02562 }
02563 }
02564 }
02565 }
02566 if (iChosen>=0) {
02567 int iSave=sequenceIn_;
02568 sequenceIn_=iChosen;
02569 unpack(rowArray_[1]);
02570 sequenceIn_ = iSave;
02571
02572 if (numberDualInfeasibilities_) {
02573 if (fabs(changeCost)>1.0e-5)
02574 printf("Odd free/unbounded combo\n");
02575 changeCost += cost_[iChosen];
02576 }
02577 problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
02578 changeCost);
02579 rowArray_[1]->clear();
02580 } else {
02581 problemStatus_=-3;
02582 }
02583 if (problemStatus_==2&&perturbation_==101) {
02584 perturbation_=102;
02585 cleanDuals=1;
02586 createRim(4);
02587 problemStatus_=-1;
02588 }
02589 if (problemStatus_==2) {
02590
02591
02592 int iColumn;
02593 double * original = columnArray_[0]->denseVector();
02594 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
02595 if(getColumnStatus(iColumn)!= basic)
02596 ray_[iColumn] +=
02597 columnActivityWork_[iColumn]-original[iColumn];
02598 columnActivityWork_[iColumn] = original[iColumn];
02599 }
02600 ClpDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
02601 rowActivityWork_);
02602 }
02603 } else {
02604 doOriginalTolerance=2;
02605 rowArray_[0]->clear();
02606 }
02607 }
02608 ClpFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
02609 ClpFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
02610 }
02611 if (problemStatus_==-4||problemStatus_==5) {
02612
02613 numberChangedBounds=changeBounds(false,NULL,changeCost);
02614
02615
02616 if (perturbation_==101&&0) {
02617 perturbation_=102;
02618 cleanDuals=1;
02619 numberChangedBounds=1;
02620
02621 changeBounds(true,NULL,changeCost);
02622 createRim(4);
02623 }
02624 if (numberChangedBounds<=0||dualBound_>1.0e20||
02625 (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
02626 problemStatus_=1;
02627 if (perturbation_==101) {
02628 perturbation_=102;
02629
02630
02631
02632 }
02633 } else {
02634 normalType=false;
02635 problemStatus_=-1;
02636 cleanDuals=1;
02637 if (numberChangedBounds<=0)
02638 doOriginalTolerance=2;
02639
02640 delete [] ray_;
02641 ray_ = NULL;
02642 }
02643
02644 }
02645 } else {
02646 cleanDuals=1;
02647 }
02648 if (problemStatus_<0) {
02649 if (doOriginalTolerance==2) {
02650
02651 lastCleaned=numberIterations_;
02652 handler_->message(CLP_DUAL_ORIGINAL,messages_)
02653 <<CoinMessageEol;
02654 perturbation_=102;
02655 #if 0
02656 double * xcost = new double[numberRows_+numberColumns_];
02657 double * xlower = new double[numberRows_+numberColumns_];
02658 double * xupper = new double[numberRows_+numberColumns_];
02659 double * xdj = new double[numberRows_+numberColumns_];
02660 double * xsolution = new double[numberRows_+numberColumns_];
02661 memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
02662 memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
02663 memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
02664 memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
02665 memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
02666 #endif
02667 createRim(4);
02668
02669 computeDuals(givenDuals);
02670 checkDualSolution();
02671 #if 0
02672 int i;
02673 for (i=0;i<numberRows_+numberColumns_;i++) {
02674 if (cost_[i]!=xcost[i])
02675 printf("** %d old cost %g new %g sol %g\n",
02676 i,xcost[i],cost_[i],solution_[i]);
02677 if (lower_[i]!=xlower[i])
02678 printf("** %d old lower %g new %g sol %g\n",
02679 i,xlower[i],lower_[i],solution_[i]);
02680 if (upper_[i]!=xupper[i])
02681 printf("** %d old upper %g new %g sol %g\n",
02682 i,xupper[i],upper_[i],solution_[i]);
02683 if (dj_[i]!=xdj[i])
02684 printf("** %d old dj %g new %g sol %g\n",
02685 i,xdj[i],dj_[i],solution_[i]);
02686 if (solution_[i]!=xsolution[i])
02687 printf("** %d old solution %g new %g sol %g\n",
02688 i,xsolution[i],solution_[i],solution_[i]);
02689 }
02690
02691
02692
02693
02694
02695 #endif
02696
02697 if (doOriginalTolerance==2) {
02698 changeMade_++;
02699 changeBounds(true,NULL,changeCost);
02700 cleanDuals=2;
02701
02702 }
02703 #if 0
02704
02705 for (i=0;i<numberRows_+numberColumns_;i++) {
02706 if (cost_[i]!=xcost[i])
02707 printf("** %d old cost %g new %g sol %g\n",
02708 i,xcost[i],cost_[i],solution_[i]);
02709 if (lower_[i]!=xlower[i])
02710 printf("** %d old lower %g new %g sol %g\n",
02711 i,xlower[i],lower_[i],solution_[i]);
02712 if (upper_[i]!=xupper[i])
02713 printf("** %d old upper %g new %g sol %g\n",
02714 i,xupper[i],upper_[i],solution_[i]);
02715 if (dj_[i]!=xdj[i])
02716 printf("** %d old dj %g new %g sol %g\n",
02717 i,xdj[i],dj_[i],solution_[i]);
02718 if (solution_[i]!=xsolution[i])
02719 printf("** %d old solution %g new %g sol %g\n",
02720 i,xsolution[i],solution_[i],solution_[i]);
02721 }
02722 delete [] xcost;
02723 delete [] xupper;
02724 delete [] xlower;
02725 delete [] xdj;
02726 delete [] xsolution;
02727 #endif
02728 }
02729 if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
02730
02731
02732 rowArray_[0]->clear();
02733 columnArray_[0]->clear();
02734 double objectiveChange=0.0;
02735 #if 0
02736 double * xcost = new double[numberRows_+numberColumns_];
02737 double * xlower = new double[numberRows_+numberColumns_];
02738 double * xupper = new double[numberRows_+numberColumns_];
02739 double * xdj = new double[numberRows_+numberColumns_];
02740 double * xsolution = new double[numberRows_+numberColumns_];
02741 memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
02742 memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
02743 memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
02744 memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
02745 memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
02746 #endif
02747 updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
02748 0.0,objectiveChange,true);
02749 #if 0
02750 int i;
02751 for (i=0;i<numberRows_+numberColumns_;i++) {
02752 if (cost_[i]!=xcost[i])
02753 printf("** %d old cost %g new %g sol %g\n",
02754 i,xcost[i],cost_[i],solution_[i]);
02755 if (lower_[i]!=xlower[i])
02756 printf("** %d old lower %g new %g sol %g\n",
02757 i,xlower[i],lower_[i],solution_[i]);
02758 if (upper_[i]!=xupper[i])
02759 printf("** %d old upper %g new %g sol %g\n",
02760 i,xupper[i],upper_[i],solution_[i]);
02761 if (dj_[i]!=xdj[i])
02762 printf("** %d old dj %g new %g sol %g\n",
02763 i,xdj[i],dj_[i],solution_[i]);
02764 if (solution_[i]!=xsolution[i])
02765 printf("** %d old solution %g new %g sol %g\n",
02766 i,xsolution[i],solution_[i],solution_[i]);
02767 }
02768 delete [] xcost;
02769 delete [] xupper;
02770 delete [] xlower;
02771 delete [] xdj;
02772 delete [] xsolution;
02773 #endif
02774
02775 gutsOfSolution(NULL,NULL);
02776 updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
02777 0.0,objectiveChange,true);
02778
02779
02780 if (numberDualInfeasibilities_||situationChanged==2)
02781 problemStatus_=-1;
02782 situationChanged=0;
02783 } else {
02784
02785 if (cleanDuals!=2)
02786 problemStatus_=-1;
02787 else
02788 problemStatus_ = 10;
02789 }
02790 }
02791 }
02792 if (type==0||type==1) {
02793 if (!type) {
02794
02795 delete [] saveStatus_;
02796 delete [] savedSolution_;
02797 saveStatus_ = new unsigned char [numberRows_+numberColumns_];
02798 savedSolution_ = new double [numberRows_+numberColumns_];
02799 }
02800
02801 memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
02802 memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
02803 numberRows_*sizeof(double));
02804 memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
02805 }
02806
02807
02808 if (tentativeStatus>-3)
02809 dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
02810 else
02811 dualRowPivot_->saveWeights(this,3);
02812
02813 if (tentativeStatus!=-2&&unflagVariables) {
02814 int iRow;
02815 for (iRow=0;iRow<numberRows_;iRow++) {
02816 int iPivot=pivotVariable_[iRow];
02817 clearFlagged(iPivot);
02818 }
02819 }
02820
02821 double limit = 0.0;
02822 getDblParam(ClpDualObjectiveLimit, limit);
02823 if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
02824 limit&&
02825 !numberAtFakeBound()&&!numberDualInfeasibilities_) {
02826 problemStatus_=1;
02827 secondaryStatus_ = 1;
02828 }
02829 if (problemStatus_<0&&!changeMade_) {
02830 problemStatus_=4;
02831 }
02832 lastGoodIteration_ = numberIterations_;
02833 #if 0
02834 double thisObj = progress->lastObjective(0);
02835 double lastObj = progress->lastObjective(1);
02836 if (lastObj>thisObj+1.0e-6*max(fabs(thisObj),fabs(lastObj))+1.0e-8
02837 &&givenDuals==NULL) {
02838 int maxFactor = factorization_->maximumPivots();
02839 if (maxFactor>10) {
02840 if (forceFactorization_<0)
02841 forceFactorization_= maxFactor;
02842 forceFactorization_ = max (1,(forceFactorization_>>1));
02843 printf("Reducing factorization frequency\n");
02844 }
02845 }
02846 #endif
02847 }
02848
02849
02850
02851
02852 void
02853 ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
02854 CoinIndexedVector * columnArray,
02855 double change)
02856 {
02857 int number;
02858 int * which;
02859
02860 int iSection;
02861
02862 for (iSection=0;iSection<2;iSection++) {
02863 int i;
02864 double * solution = solutionRegion(iSection);
02865 double * lower = lowerRegion(iSection);
02866 double * upper = upperRegion(iSection);
02867 int addSequence;
02868 if (!iSection) {
02869 number = rowArray->getNumElements();
02870 which = rowArray->getIndices();
02871 addSequence = numberColumns_;
02872 } else {
02873 number = columnArray->getNumElements();
02874 which = columnArray->getIndices();
02875 addSequence = 0;
02876 }
02877
02878 for (i=0;i<number;i++) {
02879 int iSequence = which[i];
02880 Status status = getStatus(iSequence+addSequence);
02881
02882 switch(status) {
02883
02884 case basic:
02885 case isFree:
02886 case superBasic:
02887 case ClpSimplex::isFixed:
02888 break;
02889 case atUpperBound:
02890
02891 setStatus(iSequence+addSequence,atLowerBound);
02892 solution[iSequence] = lower[iSequence];
02893 break;
02894 case atLowerBound:
02895
02896 setStatus(iSequence+addSequence,atUpperBound);
02897 solution[iSequence] = upper[iSequence];
02898 break;
02899 }
02900 }
02901 }
02902 rowArray->setNumElements(0);
02903 columnArray->setNumElements(0);
02904 }
02905
02906 void
02907 ClpSimplexDual::originalBound( int iSequence)
02908 {
02909 if (getFakeBound(iSequence)!=noFake)
02910 numberFake_--;;
02911 if (iSequence>=numberColumns_) {
02912
02913 int iRow = iSequence-numberColumns_;
02914 rowLowerWork_[iRow]=rowLower_[iRow];
02915 rowUpperWork_[iRow]=rowUpper_[iRow];
02916 if (rowScale_) {
02917 if (rowLowerWork_[iRow]>-1.0e50)
02918 rowLowerWork_[iRow] *= rowScale_[iRow];
02919 if (rowUpperWork_[iRow]<1.0e50)
02920 rowUpperWork_[iRow] *= rowScale_[iRow];
02921 }
02922 } else {
02923
02924 columnLowerWork_[iSequence]=columnLower_[iSequence];
02925 columnUpperWork_[iSequence]=columnUpper_[iSequence];
02926 if (rowScale_) {
02927 double multiplier = 1.0/columnScale_[iSequence];
02928 if (columnLowerWork_[iSequence]>-1.0e50)
02929 columnLowerWork_[iSequence] *= multiplier;
02930 if (columnUpperWork_[iSequence]<1.0e50)
02931 columnUpperWork_[iSequence] *= multiplier;
02932 }
02933 }
02934 setFakeBound(iSequence,noFake);
02935 }
02936
02937
02938 bool
02939 ClpSimplexDual::changeBound( int iSequence)
02940 {
02941
02942 double oldLower=lower_[iSequence];
02943 double oldUpper=upper_[iSequence];
02944 double value=solution_[iSequence];
02945 bool modified=false;
02946 originalBound(iSequence);
02947
02948 double lowerValue=lower_[iSequence];
02949 double upperValue=upper_[iSequence];
02950
02951 lower_[iSequence] = oldLower;
02952 upper_[iSequence] = oldUpper;
02953 if (getFakeBound(iSequence)!=noFake)
02954 numberFake_--;;
02955 if (value==oldLower) {
02956 if (upperValue > oldLower + dualBound_) {
02957 upper_[iSequence]=oldLower+dualBound_;
02958 setFakeBound(iSequence,upperFake);
02959 modified=true;
02960 numberFake_++;
02961 }
02962 } else if (value==oldUpper) {
02963 if (lowerValue < oldUpper - dualBound_) {
02964 lower_[iSequence]=oldUpper-dualBound_;
02965 setFakeBound(iSequence,lowerFake);
02966 modified=true;
02967 numberFake_++;
02968 }
02969 } else {
02970 assert(value==oldLower||value==oldUpper);
02971 }
02972 return modified;
02973 }
02974
02975 void
02976 ClpSimplexDual::perturb()
02977 {
02978 if (perturbation_>100)
02979 return;
02980 if (perturbation_==100)
02981 perturbation_=50;
02982 int savePerturbation = perturbation_;
02983 bool modifyRowCosts=false;
02984
02985 double perturbation=1.0e-20;
02986
02987 double maximumFraction = 1.0e-5;
02988 double constantPerturbation = 100.0*dualTolerance_;
02989 const int * lengths = matrix_->getVectorLengths();
02990 int maxLength=0;
02991 int minLength=numberRows_;
02992 double averageCost = 0.0;
02993 int numberNonZero=0;
02994 if (!numberIterations_&&perturbation_==50) {
02995
02996 double * sort = new double[numberColumns_];
02997
02998 const double * obj = objective();
02999 int i;
03000 for (i=0;i<numberColumns_;i++) {
03001 double value = fabs(obj[i]);
03002 sort[i]=value;
03003 averageCost += value;
03004 if (value)
03005 numberNonZero++;
03006 }
03007 if (numberNonZero)
03008 averageCost /= (double) numberNonZero;
03009 else
03010 averageCost = 1.0;
03011 std::sort(sort,sort+numberColumns_);
03012 int number=1;
03013 double last = sort[0];
03014 for (i=1;i<numberColumns_;i++) {
03015 if (last!=sort[i])
03016 number++;
03017 last=sort[i];
03018 }
03019 delete [] sort;
03020 #if 0
03021 printf("nnz %d percent %d",number,(number*100)/numberColumns_);
03022 if (number*4>numberColumns_)
03023 printf(" - Would not perturb\n");
03024 else
03025 printf(" - Would perturb\n");
03026 exit(0);
03027 #endif
03028 if (number*4>numberColumns_) {
03029 perturbation_=100;
03030 return;
03031 }
03032 }
03033 int iColumn;
03034 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
03035 if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
03036 int length = lengths[iColumn];
03037 if (length>2) {
03038 maxLength = max(maxLength,length);
03039 minLength = min(minLength,length);
03040 }
03041 }
03042 }
03043
03044 if (perturbation_>=70) {
03045 modifyRowCosts=true;
03046 perturbation_ -= 20;
03047 printf("Row costs modified, ");
03048 }
03049 bool uniformChange=false;
03050 if (perturbation_>50) {
03051
03052
03053 double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
03054 maximumFraction = m[min(perturbation_-51,10)];
03055 }
03056 int iRow;
03057 double smallestNonZero=1.0e100;
03058 numberNonZero=0;
03059 if (perturbation_>=50) {
03060 perturbation = 1.0e-8;
03061 bool allSame=true;
03062 double lastValue=0.0;
03063 for (iRow=0;iRow<numberRows_;iRow++) {
03064 double lo = rowLowerWork_[iRow];
03065 double up = rowUpperWork_[iRow];
03066 if (lo<up) {
03067 double value = fabs(rowObjectiveWork_[iRow]);
03068 perturbation = max(perturbation,value);
03069 if (value) {
03070 modifyRowCosts=true;
03071 smallestNonZero = min(smallestNonZero,value);
03072 }
03073 }
03074 if (lo&&lo>-1.0e10) {
03075 numberNonZero++;
03076 lo=fabs(lo);
03077 if (!lastValue)
03078 lastValue=lo;
03079 else if (fabs(lo-lastValue)>1.0e-7)
03080 allSame=false;
03081 }
03082 if (up&&up<1.0e10) {
03083 numberNonZero++;
03084 up=fabs(up);
03085 if (!lastValue)
03086 lastValue=up;
03087 else if (fabs(up-lastValue)>1.0e-7)
03088 allSame=false;
03089 }
03090 }
03091 double lastValue2=0.0;
03092 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
03093 double lo = columnLowerWork_[iColumn];
03094 double up = columnUpperWork_[iColumn];
03095 if (lo<up) {
03096 double value =
03097 fabs(objectiveWork_[iColumn]);
03098 perturbation = max(perturbation,value);
03099 if (value) {
03100 smallestNonZero = min(smallestNonZero,value);
03101 }
03102 }
03103 if (lo&&lo>-1.0e10) {
03104
03105 lo=fabs(lo);
03106 if (!lastValue2)
03107 lastValue2=lo;
03108 else if (fabs(lo-lastValue2)>1.0e-7)
03109 allSame=false;
03110 }
03111 if (up&&up<1.0e10) {
03112
03113 up=fabs(up);
03114 if (!lastValue2)
03115 lastValue2=up;
03116 else if (fabs(up-lastValue2)>1.0e-7)
03117 allSame=false;
03118 }
03119 }
03120 if (allSame) {
03121
03122 double smallestNegative;
03123 double largestNegative;
03124 double smallestPositive;
03125 double largestPositive;
03126 matrix_->rangeOfElements(smallestNegative,largestNegative,
03127 smallestPositive,largestPositive);
03128 if (smallestNegative==largestNegative&&
03129 smallestPositive==largestPositive) {
03130
03131 double adjust = min(100.0*maximumFraction,1.0e-3*max(lastValue,lastValue2));
03132 maximumFraction = max(adjust,maximumFraction);
03133 }
03134 }
03135 perturbation = min(perturbation,smallestNonZero/maximumFraction);
03136 } else {
03137
03138 maximumFraction = 1.0e-1;
03139
03140 if (perturbation_<=-900) {
03141 modifyRowCosts=true;
03142 perturbation_ += 1000;
03143 printf("Row costs modified, ");
03144 }
03145 if (perturbation_<=-10) {
03146 perturbation_ += 10;
03147 maximumFraction = 1.0;
03148 if ((-perturbation_)%100>=10) {
03149 uniformChange=true;
03150 perturbation_+=20;
03151 }
03152 while (perturbation_<-10) {
03153 perturbation_ += 100;
03154 maximumFraction *= 1.0e-1;
03155 }
03156 }
03157 perturbation = pow(10.0,perturbation_);
03158 }
03159 double largestZero=0.0;
03160 double largest=0.0;
03161 double largestPerCent=0.0;
03162
03163 bool printOut=(handler_->logLevel()==63);
03164 printOut=false;
03165 if (modifyRowCosts) {
03166 for (iRow=0;iRow<numberRows_;iRow++) {
03167 if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
03168 double value = perturbation;
03169 double currentValue = rowObjectiveWork_[iRow];
03170 value = min(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
03171 if (rowLowerWork_[iRow]>-largeValue_) {
03172 if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow]))
03173 value *= CoinDrand48();
03174 else
03175 value *= -CoinDrand48();
03176 } else if (rowUpperWork_[iRow]<largeValue_) {
03177 value *= -CoinDrand48();
03178 } else {
03179 value=0.0;
03180 }
03181 if (currentValue) {
03182 largest = max(largest,fabs(value));
03183 if (fabs(value)>fabs(currentValue)*largestPerCent)
03184 largestPerCent=fabs(value/currentValue);
03185 } else {
03186 largestZero=max(largestZero,fabs(value));
03187 }
03188 if (printOut)
03189 printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
03190 rowObjectiveWork_[iRow] += value;
03191 }
03192 }
03193 }
03194
03195
03196 double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
03197
03198 double extraWeight=10.0;
03199
03200 double weight2[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0};
03201 if (constantPerturbation<99.0*dualTolerance_) {
03202 perturbation *= 0.1;
03203 extraWeight=0.5;
03204 memcpy(weight,weight2,sizeof(weight2));
03205 }
03206
03207 double factor=1.0;
03208 if (maxLength) {
03209 factor = 3.0/(double) minLength;
03210 }
03211
03212 const double m1 = 0.5;
03213 double smallestAllowed = min(1.0e-2*dualTolerance_,maximumFraction);
03214 double largestAllowed = max(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
03215 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
03216 if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
03217 double value = perturbation;
03218 double currentValue = objectiveWork_[iColumn];
03219 value = min(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
03220
03221 double value2 = constantPerturbation+1.0e-1*smallestNonZero;
03222 if (uniformChange) {
03223 value = maximumFraction;
03224 value2=maximumFraction;
03225 }
03226 if (columnLowerWork_[iColumn]>-largeValue_) {
03227 if (fabs(columnLowerWork_[iColumn])<
03228 fabs(columnUpperWork_[iColumn])) {
03229 value *= (1.0-m1 +m1*CoinDrand48());
03230 value2 *= (1.0-m1 +m1*CoinDrand48());
03231 } else {
03232 value *= -(1.0-m1+m1*CoinDrand48());
03233 value2 *= -(1.0-m1+m1*CoinDrand48());
03234 }
03235 } else if (columnUpperWork_[iColumn]<largeValue_) {
03236 value *= -(1.0-m1+m1*CoinDrand48());
03237 value2 *= -(1.0-m1+m1*CoinDrand48());
03238 } else {
03239 value=0.0;
03240 }
03241 if (value) {
03242 int length = lengths[iColumn];
03243 if (length>3) {
03244 length = (int) ((double) length * factor);
03245 length = max(3,length);
03246 }
03247 double multiplier;
03248 #if 1
03249 if (length<10)
03250 multiplier=weight[length];
03251 else
03252 multiplier = weight[10];
03253 #else
03254 if (length<10)
03255 multiplier=weight[length];
03256 else
03257 multiplier = weight[10]+extraWeight*(length-10);
03258 multiplier *= 0.5;
03259 #endif
03260 value *= multiplier;
03261 value = min (value,value2);
03262 if (savePerturbation<50||savePerturbation>60) {
03263 if (fabs(value)<=dualTolerance_)
03264 value=0.0;
03265 } else if (value) {
03266
03267 if (fabs(value)<=smallestAllowed) {
03268 value *= 10.0;
03269 while (fabs(value)<=smallestAllowed)
03270 value *= 10.0;
03271 } else if (fabs(value)>largestAllowed) {
03272 value *= 0.1;
03273 while (fabs(value)>largestAllowed)
03274 value *= 0.1;
03275 }
03276 }
03277 if (currentValue) {
03278 largest = max(largest,fabs(value));
03279 if (fabs(value)>fabs(currentValue)*largestPerCent)
03280 largestPerCent=fabs(value/currentValue);
03281 } else {
03282 largestZero=max(largestZero,fabs(value));
03283 }
03284 if (printOut)
03285 printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
03286 objectiveWork_[iColumn] += value;
03287 }
03288 }
03289 }
03290 handler_->message(CLP_SIMPLEX_PERTURB,messages_)
03291 <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
03292 <<CoinMessageEol;
03293
03294
03295
03296
03297 perturbation_=101;
03298
03299 }
03300
03301
03302
03303
03304
03305
03306 int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
03307 double * newLower, double * newUpper,
03308 double ** outputSolution,
03309 int * outputStatus, int * outputIterations,
03310 bool stopOnFirstInfeasible,
03311 bool alwaysFinish)
03312 {
03313 int i;
03314 int returnCode=0;
03315 double saveObjectiveValue = objectiveValue_;
03316 #if 1
03317 algorithm_ = -1;
03318
03319
03320
03321
03322
03323 createRim(7+8+16,true);
03324
03325
03326
03327
03328
03329
03330
03331
03332 factorization_->increasingRows(2);
03333
03334 factorization_->slackValue(-1.0);
03335 factorization_->zeroTolerance(1.0e-13);
03336
03337 int factorizationStatus = internalFactorize(0);
03338 if (factorizationStatus<0)
03339 return 1;
03340 else if (factorizationStatus)
03341 handler_->message(CLP_SINGULARITIES,messages_)
03342 <<factorizationStatus
03343 <<CoinMessageEol;
03344
03345
03346 ClpFactorization saveFactorization(*factorization_);
03347
03348 double * saveSolution = new double[numberRows_+numberColumns_];
03349 memcpy(saveSolution,solution_,
03350 (numberRows_+numberColumns_)*sizeof(double));
03351 unsigned char * saveStatus =
03352 new unsigned char [numberRows_+numberColumns_];
03353 memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
03354
03355 double * saveLower = new double[numberRows_+numberColumns_];
03356 memcpy(saveLower,lower_,
03357 (numberRows_+numberColumns_)*sizeof(double));
03358 double * saveUpper = new double[numberRows_+numberColumns_];
03359 memcpy(saveUpper,upper_,
03360 (numberRows_+numberColumns_)*sizeof(double));
03361 double * saveObjective = new double[numberRows_+numberColumns_];
03362 memcpy(saveObjective,cost_,
03363 (numberRows_+numberColumns_)*sizeof(double));
03364 int * savePivot = new int [numberRows_];
03365 memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
03366
03367
03368 int iSolution = 0;
03369 for (i=0;i<numberVariables;i++) {
03370 int iColumn = variables[i];
03371 double objectiveChange;
03372 double saveBound;
03373
03374
03375
03376 saveBound = columnUpper_[iColumn];
03377
03378 columnUpper_[iColumn] = newUpper[i];
03379 if (scalingFlag_<=0)
03380 upper_[iColumn] = newUpper[i];
03381 else
03382 upper_[iColumn] = newUpper[i]/columnScale_[iColumn];
03383
03384 int status = fastDual(alwaysFinish);
03385 if (status) {
03386
03387 checkPrimalSolution(rowActivityWork_,columnActivityWork_);
03388 double limit = 0.0;
03389 getDblParam(ClpDualObjectiveLimit, limit);
03390 if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
03391 limit) {
03392 problemStatus_=0;
03393 }
03394 status=problemStatus_;
03395 }
03396 if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
03397 objectiveChange = objectiveValue_-saveObjectiveValue;
03398 } else {
03399 objectiveChange = 1.0e100;
03400 status=1;
03401 }
03402
03403 if (scalingFlag_<=0) {
03404 memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
03405 } else {
03406 int j;
03407 double * sol = outputSolution[iSolution];
03408 for (j=0;j<numberColumns_;j++)
03409 sol[j] = solution_[j]*columnScale_[j];
03410 }
03411 outputStatus[iSolution]=status;
03412 outputIterations[iSolution]=numberIterations_;
03413 iSolution++;
03414
03415 memcpy(solution_,saveSolution,
03416 (numberRows_+numberColumns_)*sizeof(double));
03417 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
03418 memcpy(lower_,saveLower,
03419 (numberRows_+numberColumns_)*sizeof(double));
03420 memcpy(upper_,saveUpper,
03421 (numberRows_+numberColumns_)*sizeof(double));
03422 memcpy(cost_,saveObjective,
03423 (numberRows_+numberColumns_)*sizeof(double));
03424 columnUpper_[iColumn] = saveBound;
03425 memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
03426 delete factorization_;
03427 factorization_ = new ClpFactorization(saveFactorization);
03428
03429 newUpper[i]=objectiveChange;
03430 #ifdef CLP_DEBUG
03431 printf("down on %d costs %g\n",iColumn,objectiveChange);
03432 #endif
03433
03434
03435
03436 saveBound = columnLower_[iColumn];
03437
03438 columnLower_[iColumn] = newLower[i];
03439 if (scalingFlag_<=0)
03440 lower_[iColumn] = newLower[i];
03441 else
03442 lower_[iColumn] = newLower[i]/columnScale_[iColumn];
03443
03444 status = fastDual(alwaysFinish);
03445 if (status) {
03446
03447 checkPrimalSolution(rowActivityWork_,columnActivityWork_);
03448 double limit = 0.0;
03449 getDblParam(ClpDualObjectiveLimit, limit);
03450 if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
03451 limit) {
03452 problemStatus_=0;
03453 }
03454 status=problemStatus_;
03455 }
03456 if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
03457 objectiveChange = objectiveValue_-saveObjectiveValue;
03458 } else {
03459 objectiveChange = 1.0e100;
03460 status=1;
03461 }
03462 if (scalingFlag_<=0) {
03463 memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
03464 } else {
03465 int j;
03466 double * sol = outputSolution[iSolution];
03467 for (j=0;j<numberColumns_;j++)
03468 sol[j] = solution_[j]*columnScale_[j];
03469 }
03470 outputStatus[iSolution]=status;
03471 outputIterations[iSolution]=numberIterations_;
03472 iSolution++;
03473
03474
03475 memcpy(solution_,saveSolution,
03476 (numberRows_+numberColumns_)*sizeof(double));
03477 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
03478 memcpy(lower_,saveLower,
03479 (numberRows_+numberColumns_)*sizeof(double));
03480 memcpy(upper_,saveUpper,
03481 (numberRows_+numberColumns_)*sizeof(double));
03482 memcpy(cost_,saveObjective,
03483 (numberRows_+numberColumns_)*sizeof(double));
03484 columnLower_[iColumn] = saveBound;
03485 memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
03486 delete factorization_;
03487 factorization_ = new ClpFactorization(saveFactorization);
03488
03489 newLower[i]=objectiveChange;
03490 #ifdef CLP_DEBUG
03491 printf("up on %d costs %g\n",iColumn,objectiveChange);
03492 #endif
03493
03494
03495
03496
03497
03498
03499 if (newUpper[i]<1.0e100) {
03500 if(newLower[i]<1.0e100) {
03501
03502 } else {
03503
03504 returnCode=1;
03505 if (stopOnFirstInfeasible)
03506 break;
03507 }
03508 } else {
03509 if(newLower[i]<1.0e100) {
03510
03511 returnCode=1;
03512 if (stopOnFirstInfeasible)
03513 break;
03514 } else {
03515
03516 returnCode=-1;
03517 break;
03518 }
03519 }
03520 }
03521 delete [] saveSolution;
03522 delete [] saveLower;
03523 delete [] saveUpper;
03524 delete [] saveObjective;
03525 delete [] saveStatus;
03526 delete [] savePivot;
03527
03528
03529 deleteRim();
03530 #else
03531
03532 double * rowSolution = new double[numberRows_];
03533 memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
03534 double * columnSolution = new double[numberColumns_];
03535 memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
03536 unsigned char * saveStatus =
03537 new unsigned char [numberRows_+numberColumns_];
03538 if (!status_) {
03539
03540 status_ = new unsigned char [numberColumns_+numberRows_];
03541 memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
03542 }
03543 memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
03544 int iSolution =0;
03545 for (i=0;i<numberVariables;i++) {
03546 int iColumn = variables[i];
03547 double objectiveChange;
03548
03549
03550
03551 double saveUpper = columnUpper_[iColumn];
03552 columnUpper_[iColumn] = newUpper[i];
03553 int status=dual(0);
03554 memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
03555 outputStatus[iSolution]=status;
03556 outputIterations[iSolution]=numberIterations_;
03557 iSolution++;
03558
03559
03560 columnUpper_[iColumn] = saveUpper;
03561 memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
03562 memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
03563 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
03564
03565 if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
03566 objectiveChange = objectiveValue_-saveObjectiveValue;
03567 } else {
03568 objectiveChange = 1.0e100;
03569 }
03570 newUpper[i]=objectiveChange;
03571 #ifdef CLP_DEBUG
03572 printf("down on %d costs %g\n",iColumn,objectiveChange);
03573 #endif
03574
03575
03576
03577 double saveLower = columnLower_[iColumn];
03578 columnLower_[iColumn] = newLower[i];
03579 status=dual(0);
03580 memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
03581 outputStatus[iSolution]=status;
03582 outputIterations[iSolution]=numberIterations_;
03583 iSolution++;
03584
03585
03586 columnLower_[iColumn] = saveLower;
03587 memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
03588 memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
03589 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
03590
03591 if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
03592 objectiveChange = objectiveValue_-saveObjectiveValue;
03593 } else {
03594 objectiveChange = 1.0e100;
03595 }
03596 newLower[i]=objectiveChange;
03597 #ifdef CLP_DEBUG
03598 printf("up on %d costs %g\n",iColumn,objectiveChange);
03599 #endif
03600
03601
03602
03603
03604
03605
03606 if (newUpper[i]<1.0e100) {
03607 if(newLower[i]<1.0e100) {
03608
03609 } else {
03610
03611 returnCode=1;
03612 if (stopOnFirstInfeasible)
03613 break;
03614 }
03615 } else {
03616 if(newLower[i]<1.0e100) {
03617
03618 returnCode=1;
03619 if (stopOnFirstInfeasible)
03620 break;
03621 } else {
03622
03623 returnCode=-1;
03624 break;
03625 }
03626 }
03627 }
03628 delete [] rowSolution;
03629 delete [] columnSolution;
03630 delete [] saveStatus;
03631 #endif
03632 objectiveValue_ = saveObjectiveValue;
03633 return returnCode;
03634 }
03635
03636 int ClpSimplexDual::fastDual(bool alwaysFinish)
03637 {
03638 algorithm_ = -1;
03639
03640 ClpDataSave data = saveData();
03641 dualTolerance_=dblParam_[ClpDualTolerance];
03642 primalTolerance_=dblParam_[ClpPrimalTolerance];
03643
03644
03645 double saveDualBound = dualBound_;
03646
03647 double objectiveChange;
03648
03649
03650 gutsOfSolution(NULL,NULL);
03651 numberFake_ =0;
03652 changeBounds(true,NULL,objectiveChange);
03653
03654 problemStatus_ = -1;
03655 numberIterations_=0;
03656 factorization_->sparseThreshold(0);
03657 factorization_->goSparse();
03658
03659 int lastCleaned=0;
03660
03661
03662 numberTimesOptimal_=0;
03663
03664
03665 int factorType=0;
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675
03676
03677
03678
03679
03680
03681
03682
03683
03684
03685
03686
03687 int returnCode = 0;
03688
03689 int iRow,iColumn;
03690 while (problemStatus_<0) {
03691
03692 for (iRow=0;iRow<4;iRow++) {
03693 rowArray_[iRow]->clear();
03694 }
03695
03696 for (iColumn=0;iColumn<2;iColumn++) {
03697 columnArray_[iColumn]->clear();
03698 }
03699
03700
03701
03702 matrix_->refresh(this);
03703
03704
03705 statusOfProblemInDual(lastCleaned,factorType,progress_,NULL);
03706
03707
03708 factorType=1;
03709
03710
03711 if (problemStatus_<0) {
03712 double * givenPi=NULL;
03713 returnCode = whileIterating(givenPi);
03714 if (!alwaysFinish&&(returnCode<1||returnCode==3)) {
03715 double limit = 0.0;
03716 getDblParam(ClpDualObjectiveLimit, limit);
03717 if(fabs(limit)>1.0e30||objectiveValue()*optimizationDirection_<
03718 limit||
03719 numberAtFakeBound()) {
03720 returnCode=1;
03721 secondaryStatus_ = 1;
03722
03723 #ifdef CLP_DEBUG
03724 printf("returning from fastDual after %d iterations with code %d\n",
03725 numberIterations_,returnCode);
03726 #endif
03727 break;
03728 }
03729 }
03730 returnCode=0;
03731 }
03732 }
03733
03734
03735 for (iRow=0;iRow<4;iRow++) {
03736 rowArray_[iRow]->clear();
03737 }
03738
03739 for (iColumn=0;iColumn<2;iColumn++) {
03740 columnArray_[iColumn]->clear();
03741 }
03742 assert(!numberFake_||returnCode||problemStatus_);
03743
03744 restoreData(data);
03745 dualBound_ = saveDualBound;
03746 return returnCode;
03747 }
03748
03749
03750 int
03751 ClpSimplexDual::numberAtFakeBound()
03752 {
03753 int iSequence;
03754 int numberFake=0;
03755
03756 for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
03757 FakeBound bound = getFakeBound(iSequence);
03758 switch(getStatus(iSequence)) {
03759
03760 case basic:
03761 break;
03762 case isFree:
03763 case superBasic:
03764 case ClpSimplex::isFixed:
03765 assert (bound==noFake);
03766 break;
03767 case atUpperBound:
03768 if (bound==upperFake||bound==bothFake)
03769 numberFake++;
03770 break;
03771 case atLowerBound:
03772 if (bound==lowerFake||bound==bothFake)
03773 numberFake++;
03774 break;
03775 }
03776 }
03777 numberFake_ = numberFake;
03778 return numberFake;
03779 }
03780
03781
03782
03783
03784
03785
03786 int
03787 ClpSimplexDual::pivotResult()
03788 {
03789 abort();
03790 return 0;
03791 }
03792
03793
03794
03795
03796
03797 int
03798 ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
03799 CoinIndexedVector * columnArray,
03800 double acceptablePivot,
03801 CoinBigIndex * dubiousWeights)
03802 {
03803 double * work;
03804 int number;
03805 int * which;
03806 int iSection;
03807
03808 double tolerance = dualTolerance_*1.001;
03809
03810 double thetaDown = 1.0e31;
03811 double changeDown ;
03812 double thetaUp = 1.0e31;
03813 double bestAlphaDown = acceptablePivot*0.99999;
03814 double bestAlphaUp = acceptablePivot*0.99999;
03815 int sequenceDown =-1;
03816 int sequenceUp = sequenceOut_;
03817
03818 double djBasic = dj_[sequenceOut_];
03819 if (djBasic>0.0) {
03820
03821
03822 thetaUp = djBasic;
03823 changeDown = -lower_[sequenceOut_];
03824 } else {
03825
03826
03827 thetaUp = -djBasic;
03828 changeDown = upper_[sequenceOut_];
03829 }
03830 bestAlphaUp = 1.0;
03831 int addSequence;
03832
03833 double alphaUp=0.0;
03834 double alphaDown=0.0;
03835
03836 for (iSection=0;iSection<2;iSection++) {
03837
03838 int i;
03839 if (!iSection) {
03840 work = rowArray->denseVector();
03841 number = rowArray->getNumElements();
03842 which = rowArray->getIndices();
03843 addSequence = numberColumns_;
03844 } else {
03845 work = columnArray->denseVector();
03846 number = columnArray->getNumElements();
03847 which = columnArray->getIndices();
03848 addSequence = 0;
03849 }
03850
03851 for (i=0;i<number;i++) {
03852 int iSequence = which[i];
03853 int iSequence2 = iSequence + addSequence;
03854 double alpha;
03855 double oldValue;
03856 double value;
03857
03858 switch(getStatus(iSequence2)) {
03859
03860 case basic:
03861 break;
03862 case ClpSimplex::isFixed:
03863 alpha = work[i];
03864 changeDown += alpha*upper_[iSequence2];
03865 break;
03866 case isFree:
03867 case superBasic:
03868 alpha = work[i];
03869
03870 if (fabs(alpha)>bestAlphaUp) {
03871 thetaDown = 0.0;
03872 thetaUp = 0.0;
03873 bestAlphaDown = fabs(alpha);
03874 bestAlphaUp = bestAlphaUp;
03875 sequenceDown =iSequence2;
03876 sequenceUp = sequenceDown;
03877 alphaUp = alpha;
03878 alphaDown = alpha;
03879 }
03880 break;
03881 case atUpperBound:
03882 alpha = work[i];
03883 oldValue = dj_[iSequence2];
03884 changeDown += alpha*upper_[iSequence2];
03885 if (alpha>=acceptablePivot) {
03886
03887 value = oldValue+thetaUp*alpha;
03888 if (value>-tolerance) {
03889 if (value>tolerance||fabs(alpha)>bestAlphaUp) {
03890 thetaUp = -oldValue/alpha;
03891 bestAlphaUp = fabs(alpha);
03892 sequenceUp = iSequence2;
03893 alphaUp = alpha;
03894 }
03895 }
03896 } else if (alpha<=-acceptablePivot) {
03897
03898 value = oldValue-thetaDown*alpha;
03899 if (value>-tolerance) {
03900 if (value>tolerance||fabs(alpha)>bestAlphaDown) {
03901 thetaDown = oldValue/alpha;
03902 bestAlphaDown = fabs(alpha);
03903 sequenceDown = iSequence2;
03904 alphaDown = alpha;
03905 }
03906 }
03907 }
03908 break;
03909 case atLowerBound:
03910 alpha = work[i];
03911 oldValue = dj_[iSequence2];
03912 changeDown += alpha*lower_[iSequence2];
03913 if (alpha<=-acceptablePivot) {
03914
03915 value = oldValue+thetaUp*alpha;
03916 if (value<tolerance) {
03917 if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
03918 thetaUp = -oldValue/alpha;
03919 bestAlphaUp = fabs(alpha);
03920 sequenceUp = iSequence2;
03921 alphaUp = alpha;
03922 }
03923 }
03924 } else if (alpha>=acceptablePivot) {
03925
03926 value = oldValue-thetaDown*alpha;
03927 if (value<tolerance) {
03928 if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
03929 thetaDown = oldValue/alpha;
03930 bestAlphaDown = fabs(alpha);
03931 sequenceDown = iSequence2;
03932 alphaDown = alpha;
03933 }
03934 }
03935 }
03936 break;
03937 }
03938 }
03939 }
03940 int returnCode;
03941 thetaUp *= -1.0;
03942 double changeUp = -thetaUp * changeDown;
03943 changeDown = -thetaDown*changeDown;
03944 if (max(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
03945
03946 if (fabs(alphaDown)<fabs(alphaUp)) {
03947 sequenceDown =-1;
03948 }
03949 }
03950
03951 if (changeDown>changeUp&&sequenceDown>=0) {
03952 theta_ = thetaDown;
03953 sequenceIn_ = sequenceDown;
03954 alpha_ = alphaDown;
03955 #ifdef CLP_DEBUG
03956 if ((handler_->logLevel()&32))
03957 printf("predicted way - dirout %d, change %g,%g theta %g\n",
03958 directionOut_,changeDown,changeUp,theta_);
03959 #endif
03960 returnCode = 0;
03961 } else {
03962 theta_ = thetaUp;
03963 sequenceIn_ = sequenceUp;
03964 alpha_ = alphaUp;
03965 if (sequenceIn_!=sequenceOut_) {
03966 #ifdef CLP_DEBUG
03967 if ((handler_->logLevel()&32))
03968 printf("opposite way - dirout %d, change %g,%g theta %g\n",
03969 directionOut_,changeDown,changeUp,theta_);
03970 #endif
03971 returnCode = 0;
03972 } else {
03973 #ifdef CLP_DEBUG
03974 if ((handler_->logLevel()&32))
03975 printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
03976 directionOut_,changeDown,changeUp,theta_);
03977 #endif
03978 returnCode = 1;
03979 }
03980 }
03981 return returnCode;
03982 }
03983
03984
03985
03986
03987 void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
03988 {
03989
03990 CoinPackedMatrix * columnCopy = matrix();
03991
03992 CoinPackedMatrix copy;
03993 copy.reverseOrderedCopyOf(*columnCopy);
03994
03995 const int * column = copy.getIndices();
03996 const CoinBigIndex * rowStart = copy.getVectorStarts();
03997 const int * rowLength = copy.getVectorLengths();
03998 const double * elementByRow = copy.getElements();
03999 double tolerance = dualTolerance_*1.001;
04000
04001 int iRow;
04002 #ifdef CLP_DEBUG
04003 {
04004 double value5=0.0;
04005 int i;
04006 for (i=0;i<numberRows_+numberColumns_;i++) {
04007 if (dj[i]<-1.0e-6)
04008 value5 += dj[i]*upper_[i];
04009 else if (dj[i] >1.0e-6)
04010 value5 += dj[i]*lower_[i];
04011 }
04012 printf("Values objective Value before %g\n",value5);
04013 }
04014 #endif
04015
04016 double * scaled=NULL;
04017 if (rowScale_)
04018 scaled = new double[numberColumns_];
04019 for (iRow=0;iRow<numberRows_;iRow++) {
04020
04021 int iSequence = iRow + numberColumns_;
04022 double djBasic = dj[iSequence];
04023 if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
04024
04025 double changeUp ;
04026
04027 if (djBasic>0.0) {
04028
04029 changeUp = -lower_[iSequence];
04030 } else {
04031
04032 changeUp = upper_[iSequence];
04033 }
04034 bool canMove=true;
04035 int i;
04036 const double * thisElements = elementByRow + rowStart[iRow];
04037 const int * thisIndices = column+rowStart[iRow];
04038 if (rowScale_) {
04039
04040 double scale = rowScale_[iRow];
04041 for (i = 0;i<rowLength[iRow];i++) {
04042 int iColumn = thisIndices[i];
04043 double alpha = thisElements[i];
04044 scaled[i] = scale*alpha*columnScale_[iColumn];
04045 }
04046 thisElements = scaled;
04047 }
04048 for (i = 0;i<rowLength[iRow];i++) {
04049 int iColumn = thisIndices[i];
04050 double alpha = thisElements[i];
04051 double oldValue = dj[iColumn];;
04052 double value;
04053
04054 switch(getStatus(iColumn)) {
04055
04056 case basic:
04057 if (dj[iColumn]<-tolerance&&
04058 fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
04059
04060 changeUp += alpha*upper_[iColumn];
04061
04062 value = oldValue+djBasic*alpha;
04063 if (value>tolerance)
04064 canMove=false;
04065 } else if (dj[iColumn]>tolerance&&
04066 fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
04067 changeUp += alpha*lower_[iColumn];
04068
04069 value = oldValue+djBasic*alpha;
04070 if (value<-tolerance)
04071 canMove=false;
04072 } else {
04073 canMove=false;
04074 }
04075 break;
04076 case ClpSimplex::isFixed:
04077 changeUp += alpha*upper_[iColumn];
04078 break;
04079 case isFree:
04080 case superBasic:
04081 canMove=false;
04082 break;
04083 case atUpperBound:
04084 changeUp += alpha*upper_[iColumn];
04085
04086 value = oldValue+djBasic*alpha;
04087 if (value>tolerance)
04088 canMove=false;
04089 break;
04090 case atLowerBound:
04091 changeUp += alpha*lower_[iColumn];
04092
04093 value = oldValue+djBasic*alpha;
04094 if (value<-tolerance)
04095 canMove=false;
04096 break;
04097 }
04098 }
04099 if (canMove) {
04100 if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
04101
04102 for (i = 0;i<rowLength[iRow];i++) {
04103 int iColumn = thisIndices[i];
04104 double alpha = thisElements[i];
04105 dj[iColumn] += djBasic * alpha;
04106 }
04107 dj[iSequence]=0.0;
04108 #ifdef CLP_DEBUG
04109 {
04110 double value5=0.0;
04111 int i;
04112 for (i=0;i<numberRows_+numberColumns_;i++) {
04113 if (dj[i]<-1.0e-6)
04114 value5 += dj[i]*upper_[i];
04115 else if (dj[i] >1.0e-6)
04116 value5 += dj[i]*lower_[i];
04117 }
04118 printf("Values objective Value after row %d old dj %g %g\n",
04119 iRow,djBasic,value5);
04120 }
04121 #endif
04122 }
04123 }
04124 }
04125 }
04126 delete [] scaled;
04127 }
04128 int
04129 ClpSimplexDual::nextSuperBasic()
04130 {
04131 if (firstFree_>=0) {
04132 int returnValue=firstFree_;
04133 int iColumn=firstFree_+1;
04134 for (;iColumn<numberRows_+numberColumns_;iColumn++) {
04135 if (getStatus(iColumn)==isFree)
04136 if (fabs(dj_[iColumn])>1.0e2*dualTolerance_)
04137 break;
04138 }
04139 firstFree_ = iColumn;
04140 if (firstFree_==numberRows_+numberColumns_)
04141 firstFree_=-1;
04142 return returnValue;
04143 } else {
04144 return -1;
04145 }
04146 }