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 #include "CoinPragma.hpp"
00081
00082 #include <math.h>
00083
00084 #include "CoinHelperFunctions.hpp"
00085 #include "ClpSimplexPrimal.hpp"
00086 #include "ClpFactorization.hpp"
00087 #include "ClpNonLinearCost.hpp"
00088 #include "CoinPackedMatrix.hpp"
00089 #include "CoinIndexedVector.hpp"
00090 #include "CoinWarmStartBasis.hpp"
00091 #include "ClpPrimalColumnPivot.hpp"
00092 #include "ClpMessage.hpp"
00093 #include <cfloat>
00094 #include <cassert>
00095 #include <string>
00096 #include <stdio.h>
00097 #include <iostream>
00098
00099 int ClpSimplexPrimal::primal (int ifValuesPass )
00100 {
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
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 algorithm_ = +1;
00183
00184
00185
00186 ClpDataSave data = saveData();
00187
00188
00189 int initialStatus=problemStatus_;
00190
00191 if (!startup(ifValuesPass)) {
00192
00193
00194 nonLinearCost_->setAverageTheta(1.0e3);
00195 int lastCleaned=0;
00196
00197
00198 pivotRow_=-2;
00199
00200
00201 int factorType=0;
00202 if (problemStatus_<0&&perturbation_<100) {
00203 perturb(0);
00204
00205 assert (!ifValuesPass);
00206 gutsOfSolution(NULL,NULL);
00207 if (handler_->logLevel()>2) {
00208 handler_->message(CLP_SIMPLEX_STATUS,messages_)
00209 <<numberIterations_<<objectiveValue();
00210 handler_->printing(sumPrimalInfeasibilities_>0.0)
00211 <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
00212 handler_->printing(sumDualInfeasibilities_>0.0)
00213 <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
00214 handler_->printing(numberDualInfeasibilitiesWithoutFree_
00215 <numberDualInfeasibilities_)
00216 <<numberDualInfeasibilitiesWithoutFree_;
00217 handler_->message()<<CoinMessageEol;
00218 }
00219 }
00220 ClpSimplex * saveModel=NULL;
00221 int stopSprint=-1;
00222 int sprintPass=0;
00223 int reasonableSprintIteration=0;
00224 int lastSprintIteration=0;
00225 double lastObjectiveValue=COIN_DBL_MAX;
00226
00227 progress_->startCheck();
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 while (problemStatus_<0) {
00240 int iRow,iColumn;
00241
00242 for (iRow=0;iRow<4;iRow++) {
00243 rowArray_[iRow]->clear();
00244 }
00245
00246 for (iColumn=0;iColumn<2;iColumn++) {
00247 columnArray_[iColumn]->clear();
00248 }
00249
00250
00251
00252 matrix_->refresh(this);
00253
00254 #if 1
00255 if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
00256 &&initialStatus!=10)
00257 perturb(1);
00258 #endif
00259
00260 if (lastGoodIteration_==numberIterations_&&factorType)
00261 factorType=3;
00262 if (saveModel) {
00263
00264 if (sequenceIn_<0||numberIterations_>=stopSprint) {
00265 problemStatus_=-1;
00266 originalModel(saveModel);
00267 saveModel=NULL;
00268 if (sequenceIn_<0&&numberIterations_<reasonableSprintIteration&&
00269 sprintPass>100)
00270 primalColumnPivot_->switchOffSprint();
00271
00272 printf("End small model\n");
00273 }
00274 }
00275
00276
00277 statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
00278
00279 if (numberDualInfeasibilities_==-776) {
00280
00281 problemStatus_=-1;
00282 originalModel(saveModel);
00283 saveModel=NULL;
00284
00285 printf("End small model after\n");
00286 statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
00287 }
00288 int numberSprintIterations=0;
00289 int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
00290 if (problemStatus_==777) {
00291
00292 problemStatus_=-1;
00293 originalModel(saveModel);
00294 saveModel=NULL;
00295
00296
00297 statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
00298 } else if (problemStatus_<0&&!saveModel&&numberSprintColumns&&firstFree_<0) {
00299 int numberSort=0;
00300 int numberFixed=0;
00301 int numberBasic=0;
00302 reasonableSprintIteration = numberIterations_ + 100;
00303 int * whichColumns = new int[numberColumns_];
00304 double * weight = new double[numberColumns_];
00305 int numberNegative=0;
00306 double sumNegative = 0.0;
00307
00308 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
00309 double dj = dj_[iColumn];
00310 switch(getColumnStatus(iColumn)) {
00311
00312 case basic:
00313 dj = -1.0e50;
00314 numberBasic++;
00315 break;
00316 case atUpperBound:
00317 dj = -dj;
00318 break;
00319 case isFixed:
00320 dj=1.0e50;
00321 numberFixed++;
00322 break;
00323 case atLowerBound:
00324 dj = dj;
00325 break;
00326 case isFree:
00327 dj = -100.0*fabs(dj);
00328 break;
00329 case superBasic:
00330 dj = -100.0*fabs(dj);
00331 break;
00332 }
00333 if (dj<-dualTolerance_&&dj>-1.0e50) {
00334 numberNegative++;
00335 sumNegative -= dj;
00336 }
00337 weight[iColumn]=dj;
00338 whichColumns[iColumn] = iColumn;
00339 }
00340 handler_->message(CLP_SPRINT,messages_)
00341 <<sprintPass<<numberIterations_-lastSprintIteration<<objectiveValue()<<sumNegative
00342 <<numberNegative
00343 <<CoinMessageEol;
00344 sprintPass++;
00345 lastSprintIteration=numberIterations_;
00346 if (objectiveValue()*optimizationDirection_>lastObjectiveValue-1.0e-7&&sprintPass>5) {
00347
00348 printf("Switching off sprint\n");
00349 primalColumnPivot_->switchOffSprint();
00350 } else {
00351 lastObjectiveValue = objectiveValue()*optimizationDirection_;
00352
00353 CoinSort_2(weight,weight+numberColumns_,whichColumns);
00354 numberSort = min(numberColumns_-numberFixed,numberBasic+numberSprintColumns);
00355
00356 std::sort(whichColumns,whichColumns+numberSort);
00357 saveModel = new ClpSimplex(this,numberSort,whichColumns);
00358 delete [] whichColumns;
00359 delete [] weight;
00360
00361
00362
00363 stopSprint = numberIterations_+numberSprintIterations;
00364 printf("Sprint with %d columns for %d iterations\n",
00365 numberSprintColumns,numberSprintIterations);
00366 }
00367 }
00368
00369
00370 factorType=1;
00371
00372
00373 pivotRow_=-2;
00374
00375
00376 if (problemStatus_>=0)
00377 break;
00378
00379
00380 whileIterating(ifValuesPass);
00381 }
00382 }
00383
00384 if (problemStatus_==1) {
00385 infeasibilityCost_=0.0;
00386 createRim(7);
00387 nonLinearCost_->checkInfeasibilities(0.0);
00388 sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
00389 numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
00390
00391 computeDuals(NULL);
00392 }
00393
00394 unflag();
00395 finish();
00396 restoreData(data);
00397 return problemStatus_;
00398 }
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 int
00410 ClpSimplexPrimal::whileIterating(int valuesOption)
00411 {
00412
00413 int ifValuesPass=(firstFree_>=0) ? 1 : 0;
00414 int returnCode=-1;
00415 int superBasicType=1;
00416 if (valuesOption>1)
00417 superBasicType=3;
00418
00419
00420 while (problemStatus_==-1) {
00421
00422 #ifdef CLP_DEBUG
00423 {
00424 int i;
00425
00426 for (i=0;i<4;i++) {
00427 if (i!=1)
00428 rowArray_[i]->checkClear();
00429 }
00430 for (i=0;i<2;i++) {
00431 columnArray_[i]->checkClear();
00432 }
00433 }
00434 #endif
00435 #if CLP_DEBUG>2
00436
00437 if (numberIterations_>0&&numberIterations_<-2534) {
00438 handler_->setLogLevel(63);
00439 double saveValue = objectiveValue_;
00440 double * saveRow1 = new double[numberRows_];
00441 double * saveRow2 = new double[numberRows_];
00442 memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
00443 memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
00444 double * saveColumn1 = new double[numberColumns_];
00445 double * saveColumn2 = new double[numberColumns_];
00446 memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
00447 memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
00448 createRim(7);
00449 gutsOfSolution(NULL,NULL);
00450 printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
00451 numberIterations_,
00452 saveValue,objectiveValue_,sumPrimalInfeasibilities_);
00453 memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
00454 memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
00455 memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
00456 memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
00457 delete [] saveRow1;
00458 delete [] saveRow2;
00459 delete [] saveColumn1;
00460 delete [] saveColumn2;
00461 objectiveValue_=saveValue;
00462 }
00463 #endif
00464 if (!ifValuesPass) {
00465
00466
00467
00468
00469
00470 primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
00471 columnArray_[0],columnArray_[1]);
00472 } else {
00473
00474 int sequenceIn=nextSuperBasic(superBasicType,columnArray_[0]);
00475 if (valuesOption>1)
00476 superBasicType=2;
00477 if (sequenceIn<0) {
00478
00479 handler_->message(CLP_END_VALUES_PASS,messages_)
00480 <<numberIterations_;
00481 primalColumnPivot_->saveWeights(this,5);
00482 problemStatus_=-2;
00483 pivotRow_=-1;
00484 returnCode=-4;
00485
00486 int i;
00487 for (i=0;i<numberRows_+numberColumns_;i++) {
00488 if (getColumnStatus(i)==atLowerBound||getColumnStatus(i)==isFixed)
00489 solution_[i]=lower_[i];
00490 else if (getColumnStatus(i)==atUpperBound)
00491 solution_[i]=upper_[i];
00492 }
00493 break;
00494 } else {
00495
00496 sequenceIn_ = sequenceIn;
00497 valueIn_=solution_[sequenceIn_];
00498 lowerIn_=lower_[sequenceIn_];
00499 upperIn_=upper_[sequenceIn_];
00500 dualIn_=dj_[sequenceIn_];
00501 }
00502 }
00503 pivotRow_=-1;
00504 sequenceOut_=-1;
00505 rowArray_[1]->clear();
00506 if (sequenceIn_>=0) {
00507
00508 assert (!flagged(sequenceIn_));
00509 #ifdef CLP_DEBUG
00510 if ((handler_->logLevel()&32)) {
00511 char x = isColumn(sequenceIn_) ? 'C' :'R';
00512 std::cout<<"pivot column "<<
00513 x<<sequenceWithin(sequenceIn_)<<std::endl;
00514 }
00515 #endif
00516
00517 returnCode = pivotResult(ifValuesPass);
00518 if (returnCode<-1&&returnCode>-5) {
00519 problemStatus_=-2;
00520 } else if (returnCode==-5) {
00521
00522 } else if (returnCode==2) {
00523 problemStatus_=-5;
00524 } else if (returnCode==4) {
00525 problemStatus_=-2;
00526 } else if (returnCode!=-1) {
00527 assert(returnCode==3);
00528 problemStatus_=3;
00529 }
00530 } else {
00531
00532 #ifdef CLP_DEBUG
00533 if (handler_->logLevel()&32)
00534 printf("** no column pivot\n");
00535 #endif
00536 if (nonLinearCost_->numberInfeasibilities())
00537 problemStatus_=-4;
00538 returnCode=0;
00539 break;
00540 }
00541 }
00542 if (valuesOption>1)
00543 columnArray_[0]->setNumElements(0);
00544 return returnCode;
00545 }
00546
00547 void
00548 ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type,
00549 ClpSimplexProgress * progress,
00550 bool doFactorization,
00551 ClpSimplex * originalModel)
00552 {
00553 if (type==2) {
00554
00555 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
00556 memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
00557 numberRows_*sizeof(double));
00558 memcpy(columnActivityWork_,savedSolution_ ,
00559 numberColumns_*sizeof(double));
00560 forceFactorization_=1;
00561 pivotRow_=-1;
00562 changeMade_++;
00563 }
00564 int saveThreshold = factorization_->sparseThreshold();
00565 int tentativeStatus = problemStatus_;
00566 if (problemStatus_>-3||problemStatus_==-4) {
00567
00568
00569
00570
00571
00572 if (doFactorization)
00573 primalColumnPivot_->saveWeights(this,1);
00574
00575 if (type&&doFactorization) {
00576
00577 if (internalFactorize(1)) {
00578 if (solveType_==2) {
00579
00580 problemStatus_=5;
00581 return;
00582 }
00583 #if 1
00584
00585 int saveDense = factorization_->denseThreshold();
00586 factorization_->setDenseThreshold(0);
00587 internalFactorize(2);
00588 factorization_->setDenseThreshold(saveDense);
00589 #else
00590
00591
00592 assert (type==1);
00593 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
00594 memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
00595 numberRows_*sizeof(double));
00596 memcpy(columnActivityWork_,savedSolution_ ,
00597 numberColumns_*sizeof(double));
00598 forceFactorization_=1;
00599 type = 2;
00600 assert (internalFactorize(1)==0);
00601 #endif
00602 changeMade_++;
00603 }
00604 }
00605 if (problemStatus_!=-4)
00606 problemStatus_=-3;
00607 }
00608
00609
00610
00611 createRim(4);
00612 gutsOfSolution(NULL,NULL,(firstFree_>=0));
00613
00614 if (progress->lastIterationNumber(0)==numberIterations_) {
00615 if (primalColumnPivot_->looksOptimal()) {
00616 numberDualInfeasibilities_ = 0;
00617 sumDualInfeasibilities_ = 0.0;
00618 }
00619 }
00620
00621 int loop;
00622 if (type!=2)
00623 loop = progress->looping();
00624 else
00625 loop=-1;
00626 if (loop>=0) {
00627 if (!problemStatus_) {
00628
00629 numberPrimalInfeasibilities_ = 0;
00630 sumPrimalInfeasibilities_ = 0.0;
00631 } else {
00632 problemStatus_ = loop;
00633 problemStatus_ = 10;
00634 }
00635 problemStatus_ = 10;
00636 return ;
00637 } else if (loop<-1) {
00638
00639 if (nonLinearCost_->numberInfeasibilities()&&progress->badTimes()>5&&
00640 progress->oddState()<10&&progress->oddState()>=0) {
00641 progress->newOddState();
00642 nonLinearCost_->zapCosts();
00643 }
00644
00645 gutsOfSolution(NULL,NULL);
00646 }
00647
00648 if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
00649 createRim(4,false);
00650 delete nonLinearCost_;
00651 nonLinearCost_ = new ClpNonLinearCost(this);
00652 progress->endOddState();
00653 gutsOfSolution(NULL,NULL);
00654 }
00655
00656 bool goToDual=false;
00657
00658
00659
00660 progressFlag_ = 0;
00661
00662 handler_->message(CLP_SIMPLEX_STATUS,messages_)
00663 <<numberIterations_<<nonLinearCost_->feasibleReportCost();
00664 handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
00665 <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
00666 handler_->printing(sumDualInfeasibilities_>0.0)
00667 <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
00668 handler_->printing(numberDualInfeasibilitiesWithoutFree_
00669 <numberDualInfeasibilities_)
00670 <<numberDualInfeasibilitiesWithoutFree_;
00671 handler_->message()<<CoinMessageEol;
00672 if (!primalFeasible()) {
00673 nonLinearCost_->checkInfeasibilities(primalTolerance_);
00674 gutsOfSolution(NULL,NULL);
00675 }
00676
00677 bool alwaysOptimal = (specialOptions_&1)!=0;
00678
00679 if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
00680 sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
00681
00682 numberDualInfeasibilities_ = 0;
00683 sumDualInfeasibilities_ = 0.0;
00684 numberPrimalInfeasibilities_ = 0;
00685 sumPrimalInfeasibilities_ = 0.0;
00686
00687 if (originalModel) {
00688
00689 numberDualInfeasibilities_ = -776;
00690 }
00691 }
00692
00693 if (dualFeasible()||problemStatus_==-4) {
00694 if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
00695
00696
00697
00698 double saveWeight = infeasibilityCost_;
00699
00700 ClpNonLinearCost * nonLinear = nonLinearCost_;
00701
00702
00703
00704 createRim(4);
00705 nonLinearCost_->checkInfeasibilities(0.0);
00706 gutsOfSolution(NULL,NULL);
00707
00708 infeasibilityCost_=1.0e100;
00709
00710 createRim(4);
00711 nonLinearCost_->checkInfeasibilities(primalTolerance_);
00712
00713 if (nonLinearCost_->numberInfeasibilities()==0) {
00714
00715 problemStatus_ = -1;
00716 infeasibilityCost_=saveWeight;
00717 nonLinearCost_->checkInfeasibilities(primalTolerance_);
00718 } else {
00719 nonLinearCost_=NULL;
00720
00721 int i;
00722 for (i=0;i<numberRows_+numberColumns_;i++)
00723 cost_[i] *= 1.0e-95;
00724 gutsOfSolution(NULL,NULL);
00725 nonLinearCost_=nonLinear;
00726 infeasibilityCost_=saveWeight;
00727 if ((infeasibilityCost_>=1.0e18||
00728 numberDualInfeasibilities_==0)&&perturbation_==101) {
00729 goToDual=unPerturb();
00730 if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
00731 goToDual=false;
00732 nonLinearCost_->checkInfeasibilities(primalTolerance_);
00733 numberDualInfeasibilities_=1;
00734 problemStatus_=-1;
00735 }
00736 if (infeasibilityCost_>=1.0e20||
00737 numberDualInfeasibilities_==0) {
00738
00739 delete [] ray_;
00740 ray_ = new double [numberRows_];
00741 memcpy(ray_,dual_,numberRows_*sizeof(double));
00742
00743 infeasibilityCost_=0.0;
00744 createRim(4);
00745 nonLinearCost_->checkInfeasibilities(primalTolerance_);
00746 gutsOfSolution(NULL,NULL);
00747
00748 infeasibilityCost_=1.0e30;
00749
00750 sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
00751 numberPrimalInfeasibilities_=
00752 nonLinearCost_->numberInfeasibilities();
00753 }
00754 if (infeasibilityCost_<1.0e20) {
00755 infeasibilityCost_ *= 5.0;
00756 changeMade_++;
00757 handler_->message(CLP_PRIMAL_WEIGHT,messages_)
00758 <<infeasibilityCost_
00759 <<CoinMessageEol;
00760
00761 createRim(4);
00762 nonLinearCost_->checkInfeasibilities(0.0);
00763 gutsOfSolution(NULL,NULL);
00764 problemStatus_=-1;
00765 goToDual=false;
00766 } else {
00767
00768 problemStatus_ = 1;
00769 }
00770 }
00771 } else {
00772
00773 if (perturbation_==101) {
00774 goToDual=unPerturb();
00775 lastCleaned=-1;
00776 }
00777 bool unflagged = unflag();
00778 if ( lastCleaned!=numberIterations_||unflagged) {
00779 handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
00780 <<primalTolerance_
00781 <<CoinMessageEol;
00782 if (numberTimesOptimal_<4) {
00783 numberTimesOptimal_++;
00784 changeMade_++;
00785 if (numberTimesOptimal_==1) {
00786
00787 factorization_->zeroTolerance(1.0e-15);
00788 }
00789 lastCleaned=numberIterations_;
00790 if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
00791 handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
00792 <<CoinMessageEol;
00793 double oldTolerance = primalTolerance_;
00794 primalTolerance_=dblParam_[ClpPrimalTolerance];
00795 #if 0
00796 double * xcost = new double[numberRows_+numberColumns_];
00797 double * xlower = new double[numberRows_+numberColumns_];
00798 double * xupper = new double[numberRows_+numberColumns_];
00799 double * xdj = new double[numberRows_+numberColumns_];
00800 double * xsolution = new double[numberRows_+numberColumns_];
00801 memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
00802 memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
00803 memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
00804 memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
00805 memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
00806 #endif
00807
00808 createRim(4);
00809 nonLinearCost_->checkInfeasibilities(oldTolerance);
00810 #if 0
00811 int i;
00812 for (i=0;i<numberRows_+numberColumns_;i++) {
00813 if (cost_[i]!=xcost[i])
00814 printf("** %d old cost %g new %g sol %g\n",
00815 i,xcost[i],cost_[i],solution_[i]);
00816 if (lower_[i]!=xlower[i])
00817 printf("** %d old lower %g new %g sol %g\n",
00818 i,xlower[i],lower_[i],solution_[i]);
00819 if (upper_[i]!=xupper[i])
00820 printf("** %d old upper %g new %g sol %g\n",
00821 i,xupper[i],upper_[i],solution_[i]);
00822 if (dj_[i]!=xdj[i])
00823 printf("** %d old dj %g new %g sol %g\n",
00824 i,xdj[i],dj_[i],solution_[i]);
00825 if (solution_[i]!=xsolution[i])
00826 printf("** %d old solution %g new %g sol %g\n",
00827 i,xsolution[i],solution_[i],solution_[i]);
00828 }
00829 delete [] xcost;
00830 delete [] xupper;
00831 delete [] xlower;
00832 delete [] xdj;
00833 delete [] xsolution;
00834 #endif
00835 gutsOfSolution(NULL,NULL);
00836 if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
00837 sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
00838
00839 numberDualInfeasibilities_ = 0;
00840 sumDualInfeasibilities_ = 0.0;
00841 numberPrimalInfeasibilities_ = 0;
00842 sumPrimalInfeasibilities_ = 0.0;
00843 }
00844 if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
00845 problemStatus_=0;
00846 else
00847 problemStatus_ = -1;
00848 } else {
00849 problemStatus_=0;
00850 if (lastCleaned<numberIterations_) {
00851 handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
00852 <<CoinMessageEol;
00853 }
00854 }
00855 } else {
00856 problemStatus_=0;
00857 }
00858 }
00859 } else {
00860
00861 if (problemStatus_==-5) {
00862 if (nonLinearCost_->numberInfeasibilities()) {
00863 if (infeasibilityCost_>1.0e18&&perturbation_==101) {
00864
00865 infeasibilityCost_ = 1.0e13;
00866 goToDual=unPerturb();
00867 }
00868
00869 if (infeasibilityCost_<1.0e20) {
00870 infeasibilityCost_ *= 5.0;
00871 changeMade_++;
00872 handler_->message(CLP_PRIMAL_WEIGHT,messages_)
00873 <<infeasibilityCost_
00874 <<CoinMessageEol;
00875
00876 createRim(4);
00877 gutsOfSolution(NULL,NULL);
00878 problemStatus_=-1;
00879 } else {
00880
00881 problemStatus_ = 2;
00882 }
00883 } else {
00884
00885 problemStatus_ = 2;
00886 }
00887 } else {
00888 if(type==3&&problemStatus_!=-5)
00889 unflag();
00890
00891 problemStatus_ = -1;
00892 }
00893 }
00894 if (type==0||type==1) {
00895 if (type!=1||!saveStatus_) {
00896
00897 delete [] saveStatus_;
00898 delete [] savedSolution_;
00899 saveStatus_ = new unsigned char [numberRows_+numberColumns_];
00900 savedSolution_ = new double [numberRows_+numberColumns_];
00901 }
00902
00903 memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
00904 memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
00905 numberRows_*sizeof(double));
00906 memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
00907 }
00908 if (doFactorization) {
00909
00910 if (tentativeStatus>-3)
00911 primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
00912 else
00913 primalColumnPivot_->saveWeights(this,3);
00914 if (saveThreshold) {
00915
00916 factorization_->sparseThreshold(0);
00917 factorization_->goSparse();
00918 }
00919 }
00920 if (problemStatus_<0&&!changeMade_) {
00921 problemStatus_=4;
00922 }
00923 lastGoodIteration_ = numberIterations_;
00924 if (goToDual)
00925 problemStatus_=10;
00926 #if 0
00927 double thisObj = progress->lastObjective(0);
00928 double lastObj = progress->lastObjective(1);
00929 if (lastObj<thisObj-1.0e-7*max(fabs(thisObj),fabs(lastObj))-1.0e-8
00930 &&firstFree_<0) {
00931 int maxFactor = factorization_->maximumPivots();
00932 if (maxFactor>10) {
00933 if (forceFactorization_<0)
00934 forceFactorization_= maxFactor;
00935 forceFactorization_ = max (1,(forceFactorization_>>1));
00936 printf("Reducing factorization frequency\n");
00937 }
00938 }
00939 #endif
00940 }
00941
00942
00943
00944
00945
00946
00947
00948 void
00949 ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
00950 CoinIndexedVector * rhsArray,
00951 CoinIndexedVector * spareArray,
00952 CoinIndexedVector * spareArray2,
00953 int valuesPass)
00954 {
00955
00956 if (valuesPass) {
00957 dualIn_ = cost_[sequenceIn_];
00958
00959 double * work=rowArray->denseVector();
00960 int number=rowArray->getNumElements();
00961 int * which=rowArray->getIndices();
00962
00963 int iIndex;
00964 for (iIndex=0;iIndex<number;iIndex++) {
00965
00966 int iRow = which[iIndex];
00967 double alpha = work[iIndex];
00968 int iPivot=pivotVariable_[iRow];
00969 dualIn_ -= alpha*cost_[iPivot];
00970 }
00971
00972 if (dualIn_<-dualTolerance_) {
00973 directionIn_=1;
00974 } else if (dualIn_>dualTolerance_) {
00975 directionIn_=-1;
00976 } else {
00977
00978 if (valueIn_-lowerIn_<upperIn_-valueIn_) {
00979 directionIn_=-1;
00980 dualIn_=dualTolerance_;
00981 } else {
00982 directionIn_=1;
00983 dualIn_=-dualTolerance_;
00984 }
00985 }
00986 }
00987
00988
00989 pivotRow_=-1;
00990 int numberRemaining=0;
00991
00992 double totalThru=0.0;
00993 double acceptablePivot=1.0e-7;
00994 if (factorization_->pivots())
00995 acceptablePivot=1.0e-5;
00996 double bestEverPivot=acceptablePivot;
00997 int lastPivotRow = -1;
00998 double lastPivot=0.0;
00999 double lastTheta=1.0e50;
01000
01001
01002
01003
01004
01005
01006
01007
01008 double * spare;
01009
01010 int * index;
01011 spareArray->clear();
01012 spare = spareArray->denseVector();
01013 index = spareArray->getIndices();
01014
01015
01016 double * rhs=rhsArray->denseVector();
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033 double * work=rowArray->denseVector();
01034 int number=rowArray->getNumElements();
01035 int * which=rowArray->getIndices();
01036
01037
01038 double way = directionIn_;
01039 double maximumMovement;
01040 if (way>0.0)
01041 maximumMovement = min(1.0e30,upperIn_-valueIn_);
01042 else
01043 maximumMovement = min(1.0e30,valueIn_-lowerIn_);
01044
01045 double averageTheta = nonLinearCost_->averageTheta();
01046 double tentativeTheta = min(10.0*averageTheta,maximumMovement);
01047 double upperTheta = maximumMovement;
01048 if (tentativeTheta>0.5*maximumMovement)
01049 tentativeTheta=maximumMovement;
01050
01051 double dualCheck = fabs(dualIn_);
01052
01053 dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
01054
01055 int iIndex;
01056 bool gotList=false;
01057 int pivotOne=-1;
01058 while (!gotList) {
01059 pivotOne=-1;
01060 totalThru=0.0;
01061
01062 numberRemaining=0;
01063 dualIn_ = cost_[sequenceIn_];
01064 double tolerance = primalTolerance_*1.002;
01065 for (iIndex=0;iIndex<number;iIndex++) {
01066
01067 int iRow = which[iIndex];
01068 double alpha = work[iIndex];
01069 int iPivot=pivotVariable_[iRow];
01070 dualIn_ -= alpha*cost_[iPivot];
01071 alpha *= way;
01072 double oldValue = solution_[iPivot];
01073
01074 if (alpha>0.0) {
01075
01076 double bound = lower_[iPivot];
01077 oldValue -= bound;
01078 } else {
01079
01080 double bound = upper_[iPivot];
01081 oldValue = bound-oldValue;
01082 }
01083
01084 double value = oldValue-tentativeTheta*fabs(alpha);
01085 assert (oldValue>=-tolerance);
01086 if (value<=tolerance) {
01087 value=oldValue-upperTheta*fabs(alpha);
01088 if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot) {
01089 upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
01090 pivotOne=numberRemaining;
01091 }
01092
01093 spare[numberRemaining]=alpha;
01094 rhs[numberRemaining]=oldValue;
01095 index[numberRemaining++]=iRow;
01096 totalThru += fabs(alpha);
01097 setActive(iRow);
01098
01099
01100
01101 }
01102 }
01103 if (upperTheta<maximumMovement&&totalThru*infeasibilityCost_>=1.0001*dualCheck) {
01104
01105 gotList=true;
01106 } else if (tentativeTheta<maximumMovement) {
01107
01108 tentativeTheta=maximumMovement;
01109 } else {
01110 gotList=true;
01111 }
01112 }
01113 totalThru=0.0;
01114
01115 int numberInRhs=numberRemaining;
01116 memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
01117 rhsArray->setNumElements(numberInRhs);
01118 rhsArray->setPacked();
01119
01120 theta_=maximumMovement;
01121
01122 bool goBackOne = false;
01123
01124
01125 int iTry=0;
01126 #define MAXTRY 1000
01127 if (numberRemaining&&upperTheta<maximumMovement) {
01128
01129 if (pivotOne>=0&&0) {
01130 double thruCost = infeasibilityCost_*fabs(spare[pivotOne]);
01131 if (thruCost>=0.99*fabs(dualIn_))
01132 printf("Could pivot on %d as change %g dj %g\n",
01133 index[pivotOne],thruCost,dualIn_);
01134 double alpha = spare[pivotOne];
01135 double oldValue = rhs[pivotOne];
01136 theta_ = oldValue/fabs(alpha);
01137 pivotRow_=pivotOne;
01138
01139 iTry=MAXTRY;
01140 }
01141
01142
01143 for ( ;iTry<MAXTRY;iTry++) {
01144
01145 upperTheta=maximumMovement;
01146 int iBest=-1;
01147 for (iIndex=0;iIndex<numberRemaining;iIndex++) {
01148
01149 double alpha = fabs(spare[iIndex]);
01150 double oldValue = rhs[iIndex];
01151 double value = oldValue-upperTheta*alpha;
01152
01153 if (value<-primalTolerance_ && alpha>=acceptablePivot) {
01154 upperTheta = (oldValue+primalTolerance_)/alpha;
01155 iBest=iIndex;
01156 }
01157 }
01158
01159
01160 double bestPivot=acceptablePivot;
01161 pivotRow_=-1;
01162 for (iIndex=0;iIndex<numberRemaining;iIndex++) {
01163
01164 int iRow = index[iIndex];
01165 double alpha = spare[iIndex];
01166 double oldValue = rhs[iIndex];
01167 double value = oldValue-upperTheta*fabs(alpha);
01168
01169 if (value<=0||iBest==iIndex) {
01170
01171 totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha,rhs[iIndex]);
01172 setActive(iRow);
01173 if (fabs(alpha)>bestPivot) {
01174 bestPivot=fabs(alpha);
01175 theta_ = oldValue/bestPivot;
01176 pivotRow_=iIndex;
01177 }
01178 }
01179 }
01180 if (bestPivot<0.1*bestEverPivot&&
01181 bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
01182
01183 goBackOne = true;
01184 break;
01185 } else if (pivotRow_==-1&&upperTheta>largeValue_) {
01186 if (lastPivot>acceptablePivot) {
01187
01188 goBackOne = true;
01189 } else {
01190
01191 }
01192 break;
01193 } else if (totalThru>=dualCheck) {
01194 break;
01195 } else {
01196 lastPivotRow=pivotRow_;
01197 lastTheta = theta_;
01198 if (bestPivot>bestEverPivot)
01199 bestEverPivot=bestPivot;
01200 }
01201 }
01202
01203 if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
01204
01205 pivotRow_=lastPivotRow;
01206 theta_ = lastTheta;
01207 }
01208 } else {
01209
01210
01211
01212
01213
01214 }
01215
01216
01217 if (pivotRow_>=0) {
01218 int position=pivotRow_;
01219 pivotRow_=index[position];
01220
01221 alpha_ = way*spare[position];
01222
01223 sequenceOut_ = pivotVariable_[pivotRow_];
01224 valueOut_ = solution(sequenceOut_);
01225 lowerOut_=lower_[sequenceOut_];
01226 upperOut_=upper_[sequenceOut_];
01227 #define MINIMUMTHETA 1.0e-12
01228
01229
01230 double minimumTheta;
01231 if (upperOut_>lowerOut_)
01232 minimumTheta=MINIMUMTHETA;
01233 else
01234 minimumTheta=0.0;
01235
01236
01237 #ifdef CLP_DEBUG
01238 bool found=false;
01239 #endif
01240 double largestInfeasibility = primalTolerance_;
01241 if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
01242 theta_=minimumTheta;
01243 for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) {
01244 largestInfeasibility = max (largestInfeasibility,
01245 -(rhs[iIndex]-fabs(spare[iIndex])*theta_));
01246 }
01247 #define CLP_DEBUG
01248 #ifdef CLP_DEBUG
01249 if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
01250 printf("Primal tolerance increased from %g to %g\n",
01251 primalTolerance_,largestInfeasibility);
01252 #endif
01253 #undef CLP_DEBUG
01254 primalTolerance_ = max(primalTolerance_,largestInfeasibility);
01255 }
01256
01257 if (theta_>tentativeTheta) {
01258 for (iIndex=0;iIndex<number;iIndex++)
01259 setActive(which[iIndex]);
01260 }
01261 if (way<0.0)
01262 theta_ = - theta_;
01263 double newValue = valueOut_ - theta_*alpha_;
01264
01265 if (alpha_*way<0.0) {
01266 directionOut_=-1;
01267 if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
01268 upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
01269 } else {
01270 upperOut_ = newValue;
01271 }
01272 } else {
01273 directionOut_=1;
01274 if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
01275 lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
01276 } else {
01277 lowerOut_ = newValue;
01278 }
01279 }
01280 dualOut_ = reducedCost(sequenceOut_);
01281 } else if (maximumMovement<1.0e20) {
01282
01283 pivotRow_ = -2;
01284 sequenceOut_ = sequenceIn_;
01285 valueOut_ = valueIn_;
01286 dualOut_ = dualIn_;
01287 lowerOut_ = lowerIn_;
01288 upperOut_ = upperIn_;
01289 alpha_ = 0.0;
01290 if (way<0.0) {
01291 directionOut_=1;
01292 theta_ = lowerOut_ - valueOut_;
01293 } else {
01294 directionOut_=-1;
01295 theta_ = upperOut_ - valueOut_;
01296 }
01297 }
01298
01299 double theta1 = max(theta_,1.0e-12);
01300 double theta2 = numberIterations_*nonLinearCost_->averageTheta();
01301
01302 nonLinearCost_->setAverageTheta((theta1+theta2)/((double) (numberIterations_+1)));
01303
01304
01305
01306
01307
01308 memset(spare,0,numberRemaining*sizeof(double));
01309
01310
01311 nonLinearCost_->goBackAll(rhsArray);
01312
01313 rhsArray->clear();
01314
01315 }
01316
01317
01318
01319
01320
01321
01322
01323 void
01324 ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
01325 CoinIndexedVector * spareRow1,
01326 CoinIndexedVector * spareRow2,
01327 CoinIndexedVector * spareColumn1,
01328 CoinIndexedVector * spareColumn2)
01329 {
01330 sequenceIn_ = primalColumnPivot_->pivotColumn(updates,spareRow1,
01331 spareRow2,spareColumn1,
01332 spareColumn2);
01333 if (sequenceIn_>=0) {
01334 valueIn_=solution_[sequenceIn_];
01335 dualIn_=dj_[sequenceIn_];
01336 if (nonLinearCost_->lookBothWays()) {
01337
01338 ClpSimplex::Status status = getStatus(sequenceIn_);
01339
01340 switch(status) {
01341 case ClpSimplex::atUpperBound:
01342 if (dualIn_<0.0) {
01343
01344 printf("For %d U (%g, %g, %g) dj changed from %g",
01345 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
01346 upper_[sequenceIn_],dualIn_);
01347 dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
01348 printf(" to %g\n",dualIn_);
01349 nonLinearCost_->setOne(sequenceIn_,upper_[sequenceIn_]+2.0*currentPrimalTolerance());
01350 setStatus(sequenceIn_,ClpSimplex::atLowerBound);
01351 }
01352 break;
01353 case ClpSimplex::atLowerBound:
01354 if (dualIn_>0.0) {
01355
01356 printf("For %d L (%g, %g, %g) dj changed from %g",
01357 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
01358 upper_[sequenceIn_],dualIn_);
01359 dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
01360 printf(" to %g\n",dualIn_);
01361 nonLinearCost_->setOne(sequenceIn_,lower_[sequenceIn_]-2.0*currentPrimalTolerance());
01362 setStatus(sequenceIn_,ClpSimplex::atUpperBound);
01363 }
01364 break;
01365 default:
01366 break;
01367 }
01368 }
01369 lowerIn_=lower_[sequenceIn_];
01370 upperIn_=upper_[sequenceIn_];
01371 if (dualIn_>0.0)
01372 directionIn_ = -1;
01373 else
01374 directionIn_ = 1;
01375 } else {
01376 sequenceIn_ = -1;
01377 }
01378 }
01379
01380
01381
01382
01383 int
01384 ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
01385 double theta,
01386 double & objectiveChange,
01387 int valuesPass)
01388 {
01389
01390 double oldCost=0.0;
01391 if (pivotRow_>=0)
01392 oldCost = cost_[sequenceOut_];
01393
01394 double * work=rowArray->denseVector();
01395 int number=rowArray->getNumElements();
01396 int * which=rowArray->getIndices();
01397
01398 int newNumber = 0;
01399 int pivotPosition = -1;
01400 nonLinearCost_->setChangeInCost(0.0);
01401 int iIndex;
01402 if (!valuesPass) {
01403 for (iIndex=0;iIndex<number;iIndex++) {
01404
01405 int iRow = which[iIndex];
01406 double alpha = work[iIndex];
01407 work[iIndex]=0.0;
01408 int iPivot=pivotVariable_[iRow];
01409 double change = theta*alpha;
01410 double value = solution_[iPivot] - change;
01411 solution_[iPivot]=value;
01412 #ifndef NDEBUG
01413
01414 if (!active(iRow)) {
01415
01416 if (change>0.0) {
01417
01418 if (value<lower_[iPivot]+primalTolerance_) {
01419 if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
01420 value=lower_[iPivot];
01421 double difference = nonLinearCost_->setOne(iPivot,value);
01422 assert (!difference);
01423 }
01424 } else {
01425
01426 if (value>upper_[iPivot]-primalTolerance_) {
01427 if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
01428 value=upper_[iPivot];
01429 double difference = nonLinearCost_->setOne(iPivot,value);
01430 assert (!difference);
01431 }
01432 }
01433 }
01434 #endif
01435 if (active(iRow)) {
01436 clearActive(iRow);
01437
01438 if (change>0.0) {
01439
01440 if (value<lower_[iPivot]+primalTolerance_) {
01441 if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
01442 value=lower_[iPivot];
01443 double difference = nonLinearCost_->setOne(iPivot,value);
01444 if (difference) {
01445 if (iRow==pivotRow_)
01446 pivotPosition=newNumber;
01447 work[newNumber] = difference;
01448
01449 dj_[iPivot] = -difference;
01450 which[newNumber++]=iRow;
01451 }
01452 }
01453 } else {
01454
01455 if (value>upper_[iPivot]-primalTolerance_) {
01456 if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
01457 value=upper_[iPivot];
01458 double difference = nonLinearCost_->setOne(iPivot,value);
01459 if (difference) {
01460 if (iRow==pivotRow_)
01461 pivotPosition=newNumber;
01462 work[newNumber] = difference;
01463
01464 dj_[iPivot] = -difference;
01465 which[newNumber++]=iRow;
01466 }
01467 }
01468 }
01469 }
01470 }
01471 } else {
01472
01473 for (iIndex=0;iIndex<number;iIndex++) {
01474
01475 int iRow = which[iIndex];
01476 double alpha = work[iIndex];
01477 work[iIndex]=0.0;
01478 int iPivot=pivotVariable_[iRow];
01479 double change = theta*alpha;
01480 double value = solution_[iPivot] - change;
01481 solution_[iPivot]=value;
01482 clearActive(iRow);
01483
01484 if (change>0.0) {
01485
01486 if (value<lower_[iPivot]+primalTolerance_) {
01487 if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
01488 value=lower_[iPivot];
01489 double difference = nonLinearCost_->setOne(iPivot,value);
01490 if (difference) {
01491 if (iRow==pivotRow_)
01492 pivotPosition=newNumber;
01493 work[newNumber] = difference;
01494
01495 dj_[iPivot] = -difference;
01496 which[newNumber++]=iRow;
01497 }
01498 }
01499 } else {
01500
01501 if (value>upper_[iPivot]-primalTolerance_) {
01502 if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
01503 value=upper_[iPivot];
01504 double difference = nonLinearCost_->setOne(iPivot,value);
01505 if (difference) {
01506 if (iRow==pivotRow_)
01507 pivotPosition=newNumber;
01508 work[newNumber] = difference;
01509
01510 dj_[iPivot] = -difference;
01511 which[newNumber++]=iRow;
01512 }
01513 }
01514 }
01515 }
01516 }
01517 objectiveChange += nonLinearCost_->changeInCost();
01518 rowArray->setPacked();
01519 #if 0
01520 rowArray->setNumElements(newNumber);
01521 rowArray->expand();
01522 if (pivotRow_>=0) {
01523 dualIn_ += (oldCost-cost_[sequenceOut_]);
01524
01525 rowArray->add(pivotRow_,-dualIn_);
01526
01527 rowArray->scanAndPack();
01528 } else {
01529
01530 rowArray->scanAndPack();
01531 }
01532 #else
01533 if (pivotRow_>=0) {
01534 dualIn_ += (oldCost-cost_[sequenceOut_]);
01535
01536 if (pivotPosition>=0) {
01537 work[pivotPosition] -= dualIn_;
01538 } else {
01539 work[newNumber]=-dualIn_;
01540 which[newNumber++]=pivotRow_;
01541 }
01542 }
01543 rowArray->setNumElements(newNumber);
01544 #endif
01545 return 0;
01546 }
01547
01548 void
01549 ClpSimplexPrimal::perturb(int type)
01550 {
01551 if (perturbation_>100)
01552 return;
01553 if (perturbation_==100)
01554 perturbation_=50;
01555 int savePerturbation = perturbation_;
01556 int i;
01557 if (!numberIterations_)
01558 cleanStatus();
01559 if (!numberIterations_&&perturbation_==50) {
01560
01561 double * sort = new double[numberRows_];
01562 for (i=0;i<numberRows_;i++) {
01563 double lo = fabs(lower_[i]);
01564 double up = fabs(upper_[i]);
01565 double value=0.0;
01566 if (lo&&lo<1.0e20) {
01567 if (up&&up<1.0e20)
01568 value = 0.5*(lo+up);
01569 else
01570 value=lo;
01571 } else {
01572 if (up&&up<1.0e20)
01573 value = up;
01574 }
01575 sort[i]=value;
01576 }
01577 std::sort(sort,sort+numberRows_);
01578 int number=1;
01579 double last = sort[0];
01580 for (i=1;i<numberRows_;i++) {
01581 if (last!=sort[i])
01582 number++;
01583 last=sort[i];
01584 }
01585 delete [] sort;
01586 if (number*4>numberRows_) {
01587 perturbation_=100;
01588 return;
01589 }
01590 }
01591
01592 double perturbation=1.0e-20;
01593 int numberNonZero=0;
01594
01595 double maximumFraction = 1.0e-5;
01596 if (perturbation_>=50) {
01597 perturbation = 1.0e-4;
01598 for (i=0;i<numberColumns_+numberRows_;i++) {
01599 if (upper_[i]>lower_[i]+primalTolerance_) {
01600 double lowerValue, upperValue;
01601 if (lower_[i]>-1.0e20)
01602 lowerValue = fabs(lower_[i]);
01603 else
01604 lowerValue=0.0;
01605 if (upper_[i]<1.0e20)
01606 upperValue = fabs(upper_[i]);
01607 else
01608 upperValue=0.0;
01609 double value = max(fabs(lowerValue),fabs(upperValue));
01610 value = min(value,upper_[i]-lower_[i]);
01611 #if 1
01612 if (value) {
01613 perturbation += value;
01614 numberNonZero++;
01615 }
01616 #else
01617 perturbation = max(perturbation,value);
01618 #endif
01619 }
01620 }
01621 if (numberNonZero)
01622 perturbation /= (double) numberNonZero;
01623 else
01624 perturbation = 1.0e-1;
01625 } else if (perturbation_<100) {
01626 perturbation = pow(10.0,perturbation_);
01627
01628 maximumFraction = 1.0;
01629 }
01630 double largestZero=0.0;
01631 double largest=0.0;
01632 double largestPerCent=0.0;
01633 bool printOut=(handler_->logLevel()==63);
01634 printOut=false;
01635
01636 int number=0;
01637 int iSequence;
01638 for (iSequence=0;iSequence<numberRows_;iSequence++) {
01639 if (getRowStatus(iSequence)==basic)
01640 number++;
01641 }
01642 if (number!=numberRows_)
01643 type=1;
01644
01645
01646
01647
01648 if (type==1) {
01649
01650 for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
01651 if (getStatus(iSequence)==basic) {
01652 double solutionValue = solution_[iSequence];
01653 double lowerValue = lower_[iSequence];
01654 double upperValue = upper_[iSequence];
01655 double difference = upperValue-lowerValue;
01656 difference = min(difference,perturbation);
01657 difference = min(difference,fabs(solutionValue)+1.0);
01658 double value = maximumFraction*(difference+1.0);
01659 value = min(value,0.1);
01660 value *= CoinDrand48();
01661 if (solutionValue-lowerValue<=primalTolerance_) {
01662 lower_[iSequence] -= value;
01663 } else if (upperValue-solutionValue<=primalTolerance_) {
01664 upper_[iSequence] += value;
01665 } else {
01666 #if 0
01667 if (iSequence>=numberColumns_) {
01668
01669 if (upperValue>1.0e30&&lowerValue<-1.0e30)
01670 value=0.0;
01671 else
01672 value = - value;
01673 } else {
01674 value = 0.0;
01675 }
01676 #else
01677 value=0.0;
01678 #endif
01679 }
01680 if (value) {
01681 if (printOut)
01682 printf("col %d lower from %g to %g, upper from %g to %g\n",
01683 iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
01684 if (solutionValue) {
01685 largest = max(largest,value);
01686 if (value>(fabs(solutionValue)+1.0)*largestPerCent)
01687 largestPerCent=value/(fabs(solutionValue)+1.0);
01688 } else {
01689 largestZero = max(largestZero,value);
01690 }
01691 }
01692 }
01693 }
01694 } else {
01695 double tolerance = 100.0*primalTolerance_;
01696 for (i=0;i<numberColumns_;i++) {
01697 double lowerValue=lower_[i], upperValue=upper_[i];
01698 if (upperValue>lowerValue+primalTolerance_) {
01699 double value = perturbation*maximumFraction;
01700 value = min(value,0.1);
01701 value *= CoinDrand48();
01702 if (savePerturbation!=50) {
01703 if (fabs(value)<=primalTolerance_)
01704 value=0.0;
01705 if (lowerValue>-1.0e20&&lowerValue)
01706 lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
01707 if (upperValue<1.0e20&&upperValue)
01708 upperValue += value * (max(1.0e-2,1.0e-5*fabs(upperValue)));
01709 } else if (value) {
01710 double valueL =value *(max(1.0e-2,1.0e-5*fabs(lowerValue)));
01711
01712 if (valueL<=tolerance) {
01713 valueL *= 10.0;
01714 while (valueL<=tolerance)
01715 valueL *= 10.0;
01716 } else if (valueL>1.0) {
01717 valueL *= 0.1;
01718 while (valueL>1.0)
01719 valueL *= 0.1;
01720 }
01721 if (lowerValue>-1.0e20&&lowerValue)
01722 lowerValue -= valueL;
01723 double valueU =value *(max(1.0e-2,1.0e-5*fabs(upperValue)));
01724
01725 if (valueU<=tolerance) {
01726 valueU *= 10.0;
01727 while (valueU<=tolerance)
01728 valueU *= 10.0;
01729 } else if (valueU>1.0) {
01730 valueU *= 0.1;
01731 while (valueU>1.0)
01732 valueU *= 0.1;
01733 }
01734 if (upperValue<1.0e20&&upperValue)
01735 upperValue += valueU;
01736 }
01737 if (lowerValue!=lower_[i]) {
01738 double difference = fabs(lowerValue-lower_[i]);
01739 largest = max(largest,difference);
01740 if (difference>fabs(lower_[i])*largestPerCent)
01741 largestPerCent=fabs(difference/lower_[i]);
01742 }
01743 if (upperValue!=upper_[i]) {
01744 double difference = fabs(upperValue-upper_[i]);
01745 largest = max(largest,difference);
01746 if (difference>fabs(upper_[i])*largestPerCent)
01747 largestPerCent=fabs(difference/upper_[i]);
01748 }
01749 if (printOut)
01750 printf("col %d lower from %g to %g, upper from %g to %g\n",
01751 i,lower_[i],lowerValue,upper_[i],upperValue);
01752 }
01753 lower_[i]=lowerValue;
01754 upper_[i]=upperValue;
01755 }
01756 for (;i<numberColumns_+numberRows_;i++) {
01757 double lowerValue=lower_[i], upperValue=upper_[i];
01758 double value = perturbation*maximumFraction;
01759 value = min(value,0.1);
01760 value *= CoinDrand48();
01761 if (upperValue>lowerValue+tolerance) {
01762 if (savePerturbation!=50) {
01763 if (fabs(value)<=primalTolerance_)
01764 value=0.0;
01765 if (lowerValue>-1.0e20&&lowerValue)
01766 lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
01767 if (upperValue<1.0e20&&upperValue)
01768 upperValue += value * (max(1.0e-2,1.0e-5*fabs(upperValue)));
01769 } else if (value) {
01770 double valueL =value *(max(1.0e-2,1.0e-5*fabs(lowerValue)));
01771
01772 if (valueL<=tolerance) {
01773 valueL *= 10.0;
01774 while (valueL<=tolerance)
01775 valueL *= 10.0;
01776 } else if (valueL>1.0) {
01777 valueL *= 0.1;
01778 while (valueL>1.0)
01779 valueL *= 0.1;
01780 }
01781 if (lowerValue>-1.0e20&&lowerValue)
01782 lowerValue -= valueL;
01783 double valueU =value *(max(1.0e-2,1.0e-5*fabs(upperValue)));
01784
01785 if (valueU<=tolerance) {
01786 valueU *= 10.0;
01787 while (valueU<=tolerance)
01788 valueU *= 10.0;
01789 } else if (valueU>1.0) {
01790 valueU *= 0.1;
01791 while (valueU>1.0)
01792 valueU *= 0.1;
01793 }
01794 if (upperValue<1.0e20&&upperValue)
01795 upperValue += valueU;
01796 }
01797 } else if (upperValue>0.0) {
01798 upperValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
01799 lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
01800 } else if (upperValue<0.0) {
01801 upperValue += value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
01802 lowerValue += value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
01803 } else {
01804 }
01805 if (lowerValue!=lower_[i]) {
01806 double difference = fabs(lowerValue-lower_[i]);
01807 largest = max(largest,difference);
01808 if (difference>fabs(lower_[i])*largestPerCent)
01809 largestPerCent=fabs(difference/lower_[i]);
01810 }
01811 if (upperValue!=upper_[i]) {
01812 double difference = fabs(upperValue-upper_[i]);
01813 largest = max(largest,difference);
01814 if (difference>fabs(upper_[i])*largestPerCent)
01815 largestPerCent=fabs(difference/upper_[i]);
01816 }
01817 if (printOut)
01818 printf("row %d lower from %g to %g, upper from %g to %g\n",
01819 i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
01820 lower_[i]=lowerValue;
01821 upper_[i]=upperValue;
01822 }
01823 }
01824
01825 for (i=0;i<numberColumns_+numberRows_;i++) {
01826 switch(getStatus(i)) {
01827
01828 case basic:
01829 break;
01830 case atUpperBound:
01831 solution_[i]=upper_[i];
01832 break;
01833 case isFixed:
01834 case atLowerBound:
01835 solution_[i]=lower_[i];
01836 break;
01837 case isFree:
01838 break;
01839 case superBasic:
01840 break;
01841 }
01842 }
01843 handler_->message(CLP_SIMPLEX_PERTURB,messages_)
01844 <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
01845 <<CoinMessageEol;
01846
01847
01848 perturbation_=101;
01849 }
01850
01851 bool
01852 ClpSimplexPrimal::unPerturb()
01853 {
01854 if (perturbation_!=101)
01855 return false;
01856
01857 createRim(7);
01858 sanityCheck();
01859
01860 unflag();
01861
01862 delete nonLinearCost_;
01863 nonLinearCost_= new ClpNonLinearCost(this);
01864 perturbation_ = 102;
01865
01866 nonLinearCost_->checkInfeasibilities(0.0);
01867 #if 1
01868
01869 return true;
01870 #else
01871 gutsOfSolution(NULL,NULL);
01872 return false;
01873 #endif
01874
01875 }
01876
01877 int
01878 ClpSimplexPrimal::unflag()
01879 {
01880 int i;
01881 int number = numberRows_+numberColumns_;
01882 int numberFlagged=0;
01883 for (i=0;i<number;i++) {
01884 if (flagged(i)) {
01885 clearFlagged(i);
01886 numberFlagged++;
01887 }
01888 }
01889 return numberFlagged;
01890 }
01891
01892 void
01893 ClpSimplexPrimal::alwaysOptimal(bool onOff)
01894 {
01895 if (onOff)
01896 specialOptions_ |= 1;
01897 else
01898 specialOptions_ &= ~1;
01899 }
01900 bool
01901 ClpSimplexPrimal::alwaysOptimal() const
01902 {
01903 return (specialOptions_&1)!=0;
01904 }
01905
01906 void
01907 ClpSimplexPrimal::exactOutgoing(bool onOff)
01908 {
01909 if (onOff)
01910 specialOptions_ |= 4;
01911 else
01912 specialOptions_ &= ~4;
01913 }
01914 bool
01915 ClpSimplexPrimal::exactOutgoing() const
01916 {
01917 return (specialOptions_&4)!=0;
01918 }
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929 int
01930 ClpSimplexPrimal::pivotResult(int ifValuesPass)
01931 {
01932
01933 bool roundAgain=true;
01934 int returnCode=-1;
01935
01936
01937 while (roundAgain) {
01938 roundAgain=false;
01939 returnCode=-1;
01940 pivotRow_=-1;
01941 sequenceOut_=-1;
01942 rowArray_[1]->clear();
01943 #if 0
01944 {
01945 int seq[]={612,643};
01946 int k;
01947 for (k=0;k<sizeof(seq)/sizeof(int);k++) {
01948 int iSeq=seq[k];
01949 if (getColumnStatus(iSeq)!=basic) {
01950 double djval;
01951 double * work;
01952 int number;
01953 int * which;
01954
01955 int iIndex;
01956 unpack(rowArray_[1],iSeq);
01957 factorization_->updateColumn(rowArray_[2],rowArray_[1]);
01958 djval = cost_[iSeq];
01959 work=rowArray_[1]->denseVector();
01960 number=rowArray_[1]->getNumElements();
01961 which=rowArray_[1]->getIndices();
01962
01963 for (iIndex=0;iIndex<number;iIndex++) {
01964
01965 int iRow = which[iIndex];
01966 double alpha = work[iRow];
01967 int iPivot=pivotVariable_[iRow];
01968 djval -= alpha*cost_[iPivot];
01969 }
01970 double comp = 1.0e-8 + 1.0e-7*(max(fabs(dj_[iSeq]),fabs(djval)));
01971 if (fabs(djval-dj_[iSeq])>comp)
01972 printf("Bad dj %g for %d - true is %g\n",
01973 dj_[iSeq],iSeq,djval);
01974 assert (fabs(djval)<1.0e-3||djval*dj_[iSeq]>0.0);
01975 rowArray_[1]->clear();
01976 }
01977 }
01978 }
01979 #endif
01980
01981
01982
01983 unpackPacked(rowArray_[1]);
01984
01985 double saveDj = dualIn_;
01986 factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
01987
01988 primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
01989 ifValuesPass);
01990 if (ifValuesPass) {
01991 saveDj=dualIn_;
01992 if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
01993 if(fabs(dualIn_)<1.0e2*dualTolerance_) {
01994
01995 directionIn_=-directionIn_;
01996 primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
01997 0);
01998 }
01999 if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
02000 if (solveType_==1) {
02001
02002 char x = isColumn(sequenceIn_) ? 'C' :'R';
02003 handler_->message(CLP_SIMPLEX_FLAG,messages_)
02004 <<x<<sequenceWithin(sequenceIn_)
02005 <<CoinMessageEol;
02006 setFlagged(sequenceIn_);
02007 progress_->clearBadTimes();
02008 lastBadIteration_ = numberIterations_;
02009 clearAll();
02010 pivotRow_=-1;
02011 }
02012 returnCode=-5;
02013 break;
02014 }
02015 }
02016 }
02017 double checkValue=1.0e-2;
02018 if (largestDualError_>1.0e-5)
02019 checkValue=1.0e-1;
02020 if (solveType_==1&&((saveDj*dualIn_<1.0e-20&&!ifValuesPass)||
02021 fabs(saveDj-dualIn_)>checkValue*(1.0+fabs(saveDj)))) {
02022 handler_->message(CLP_PRIMAL_DJ,messages_)
02023 <<saveDj<<dualIn_
02024 <<CoinMessageEol;
02025 if(lastGoodIteration_ != numberIterations_) {
02026 clearAll();
02027 pivotRow_=-1;
02028 returnCode=-4;
02029 if(lastGoodIteration_+1 == numberIterations_) {
02030
02031
02032 nonLinearCost_->checkInfeasibilities(0.0);
02033 }
02034 sequenceOut_=-1;
02035 break;
02036 } else {
02037
02038 if (saveDj*dualIn_<1.0e-20||
02039 fabs(saveDj-dualIn_)>2.0e-1*(1.0+fabs(dualIn_))) {
02040
02041 char x = isColumn(sequenceIn_) ? 'C' :'R';
02042 handler_->message(CLP_SIMPLEX_FLAG,messages_)
02043 <<x<<sequenceWithin(sequenceIn_)
02044 <<CoinMessageEol;
02045 setFlagged(sequenceIn_);
02046 progress_->clearBadTimes();
02047 lastBadIteration_ = numberIterations_;
02048 clearAll();
02049 pivotRow_=-1;
02050 returnCode=-5;
02051 sequenceOut_=-1;
02052 break;
02053 }
02054 }
02055 }
02056 if (pivotRow_>=0) {
02057 if (solveType_==2) {
02058
02059
02060 primalRay(rowArray_[1]);
02061
02062 if (pivotRow_>=0) {
02063
02064
02065
02066
02067
02068 assert (fabs(alpha_)>1.0e-8);
02069 double multiplier = dualIn_/alpha_;
02070 rowArray_[0]->insert(pivotRow_,multiplier);
02071 factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
02072
02073 matrix_->transposeTimes(this,-1.0,
02074 rowArray_[0],columnArray_[1],columnArray_[0]);
02075
02076 int i;
02077 int * index = columnArray_[0]->getIndices();
02078 int number = columnArray_[0]->getNumElements();
02079 double * element = columnArray_[0]->denseVector();
02080 for (i=0;i<number;i++) {
02081 int ii = index[i];
02082 dj_[ii] += element[ii];
02083 element[ii]=0.0;
02084 }
02085 columnArray_[0]->setNumElements(0);
02086
02087 index = rowArray_[0]->getIndices();
02088 number = rowArray_[0]->getNumElements();
02089 element = rowArray_[0]->denseVector();
02090 for (i=0;i<number;i++) {
02091 int ii = index[i];
02092 dj_[ii+numberColumns_] += element[ii];
02093 dual_[ii] = dj_[ii+numberColumns_];
02094 element[ii]=0.0;
02095 }
02096 rowArray_[0]->setNumElements(0);
02097
02098 assert (fabs(dj_[sequenceIn_])<1.0e-6);
02099 }
02100 }
02101
02102 int updateStatus = factorization_->replaceColumn(rowArray_[2],
02103 pivotRow_,
02104 alpha_);
02105
02106 if (updateStatus==2&&
02107 lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
02108 updateStatus=4;
02109 if (updateStatus==1||updateStatus==4) {
02110
02111 if (factorization_->pivots()>5||updateStatus==4) {
02112 returnCode=-3;
02113 }
02114 } else if (updateStatus==2) {
02115
02116
02117 factorization_->zeroTolerance(1.0e-15);
02118 int maxFactor = factorization_->maximumPivots();
02119 if (maxFactor>10) {
02120 if (forceFactorization_<0)
02121 forceFactorization_= maxFactor;
02122 forceFactorization_ = max (1,(forceFactorization_>>1));
02123 }
02124
02125 if(lastGoodIteration_ != numberIterations_) {
02126 clearAll();
02127 pivotRow_=-1;
02128 if (solveType_==1) {
02129 returnCode=-4;
02130 break;
02131 } else {
02132
02133 int lastCleaned;
02134 ClpSimplexProgress dummyProgress;
02135 if (saveStatus_)
02136 statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true);
02137 else
02138 statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true);
02139 roundAgain=true;
02140 continue;
02141 }
02142 } else {
02143
02144 if (solveType_==1) {
02145 char x = isColumn(sequenceIn_) ? 'C' :'R';
02146 handler_->message(CLP_SIMPLEX_FLAG,messages_)
02147 <<x<<sequenceWithin(sequenceIn_)
02148 <<CoinMessageEol;
02149 setFlagged(sequenceIn_);
02150 progress_->clearBadTimes();
02151 }
02152 lastBadIteration_ = numberIterations_;
02153 clearAll();
02154 pivotRow_=-1;
02155 sequenceOut_=-1;
02156 returnCode = -5;
02157 break;
02158
02159 }
02160 } else if (updateStatus==3) {
02161
02162
02163 if (factorization_->pivots()<
02164 0.5*factorization_->maximumPivots()&&
02165 factorization_->pivots()<200)
02166 factorization_->areaFactor(
02167 factorization_->areaFactor() * 1.1);
02168 returnCode =-2;
02169 } else if (updateStatus==5) {
02170 problemStatus_=-2;
02171 }
02172
02173 if (!ifValuesPass)
02174 primalColumnPivot_->updateWeights(rowArray_[1]);
02175 } else {
02176 if (pivotRow_==-1) {
02177
02178 rowArray_[0]->clear();
02179 if (!factorization_->pivots()) {
02180 returnCode = 2;
02181
02182 primalRay(rowArray_[1]);
02183 } else if (solveType_==2) {
02184
02185 int lastCleaned;
02186 ClpSimplexProgress dummyProgress;
02187 if (saveStatus_)
02188 statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true);
02189 else
02190 statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true);
02191 roundAgain=true;
02192 continue;
02193 } else {
02194 returnCode = 4;
02195 }
02196 break;
02197 } else {
02198
02199 }
02200 }
02201
02202
02203
02204
02205 double objectiveChange=0.0;
02206
02207 updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
02208
02209 double oldValue = valueIn_;
02210 if (directionIn_==-1) {
02211
02212 if (sequenceIn_!=sequenceOut_) {
02213
02214 valueIn_ -= fabs(theta_);
02215 } else {
02216 valueIn_=lowerIn_;
02217 }
02218 } else {
02219
02220 if (sequenceIn_!=sequenceOut_) {
02221
02222 valueIn_ += fabs(theta_);
02223 } else {
02224 valueIn_=upperIn_;
02225 }
02226 }
02227 objectiveChange += dualIn_*(valueIn_-oldValue);
02228
02229 if (sequenceIn_!=sequenceOut_) {
02230 if (directionOut_>0) {
02231 valueOut_ = lowerOut_;
02232 } else {
02233 valueOut_ = upperOut_;
02234 }
02235 if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
02236 valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
02237 else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
02238 valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
02239
02240
02241 directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
02242 solution_[sequenceOut_]=valueOut_;
02243 }
02244
02245 nonLinearCost_->setOne(sequenceIn_,valueIn_);
02246 int whatNext=housekeeping(objectiveChange);
02247
02248 #if 0
02249 if (numberIterations_==1148)
02250 whatNext=1;
02251 if (numberIterations_>1200)
02252 exit(0);
02253 #endif
02254 if (whatNext==1) {
02255 returnCode =-2;
02256 } else if (whatNext==2) {
02257
02258 returnCode=3;
02259 } else if(numberIterations_ == lastGoodIteration_
02260 + 2 * factorization_->maximumPivots()) {
02261
02262 returnCode =-2;
02263 }
02264 }
02265 if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
02266
02267 int lastCleaned;
02268 ClpSimplexProgress dummyProgress;
02269 if (saveStatus_)
02270 statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true);
02271 else
02272 statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true);
02273 if (problemStatus_==5) {
02274 printf("Singular basis\n");
02275 problemStatus_=-1;
02276 returnCode=5;
02277 }
02278 }
02279 #ifdef CLP_DEBUG
02280 {
02281 int i;
02282
02283 for (i=0;i<4;i++) {
02284 if (i!=1)
02285 rowArray_[i]->checkClear();
02286 }
02287 for (i=0;i<2;i++) {
02288 columnArray_[i]->checkClear();
02289 }
02290 }
02291 #endif
02292 return returnCode;
02293 }
02294
02295 void
02296 ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
02297 {
02298 delete [] ray_;
02299 ray_ = new double [numberColumns_];
02300 ClpFillN(ray_,numberColumns_,0.0);
02301 int number=rowArray->getNumElements();
02302 int * index = rowArray->getIndices();
02303 double * array = rowArray->denseVector();
02304 double way=-directionIn_;
02305 int i;
02306 double zeroTolerance=1.0e-12;
02307 if (sequenceIn_<numberColumns_)
02308 ray_[sequenceIn_]=directionIn_;
02309 if (!rowArray->packedMode()) {
02310 for (i=0;i<number;i++) {
02311 int iRow=index[i];
02312 int iPivot=pivotVariable_[iRow];
02313 double arrayValue = array[iRow];
02314 if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
02315 ray_[iPivot] = way* arrayValue;
02316 }
02317 } else {
02318 for (i=0;i<number;i++) {
02319 int iRow=index[i];
02320 int iPivot=pivotVariable_[iRow];
02321 double arrayValue = array[i];
02322 if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
02323 ray_[iPivot] = way* arrayValue;
02324 }
02325 }
02326 }
02327
02328
02329
02330
02331
02332 int
02333 ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
02334 {
02335 if (firstFree_>=0&&superBasicType) {
02336 int returnValue=firstFree_;
02337 int iColumn=firstFree_+1;
02338 if (superBasicType>1) {
02339 if (superBasicType>2) {
02340
02341
02342 int number=0;
02343 double * work=columnArray->denseVector();
02344 int * which=columnArray->getIndices();
02345 for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
02346 if (getStatus(iColumn)==superBasic) {
02347 if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
02348 solution_[iColumn]=lower_[iColumn];
02349 setStatus(iColumn,atLowerBound);
02350 } else if (fabs(solution_[iColumn]-upper_[iColumn])
02351 <=primalTolerance_) {
02352 solution_[iColumn]=upper_[iColumn];
02353 setStatus(iColumn,atUpperBound);
02354 } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
02355 setStatus(iColumn,isFree);
02356 break;
02357 } else if (!flagged(iColumn)) {
02358
02359 work[number]= - min(0.1*(solution_[iColumn]-lower_[iColumn]),
02360 upper_[iColumn]-solution_[iColumn]);
02361 which[number++] = iColumn;
02362 }
02363 }
02364 }
02365 CoinSort_2(work,work+number,which);
02366 columnArray->setNumElements(number);
02367 memset(work,0,number*sizeof(double));
02368 }
02369 int * which=columnArray->getIndices();
02370 int number = columnArray->getNumElements();
02371 if (!number) {
02372
02373 iColumn = numberRows_+numberColumns_;
02374 returnValue=-1;
02375 } else {
02376 number--;
02377 returnValue=which[number];
02378 iColumn=returnValue;
02379 columnArray->setNumElements(number);
02380 }
02381 } else {
02382 for (;iColumn<numberRows_+numberColumns_;iColumn++) {
02383 if (getStatus(iColumn)==superBasic) {
02384 if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
02385 solution_[iColumn]=lower_[iColumn];
02386 setStatus(iColumn,atLowerBound);
02387 } else if (fabs(solution_[iColumn]-upper_[iColumn])
02388 <=primalTolerance_) {
02389 solution_[iColumn]=upper_[iColumn];
02390 setStatus(iColumn,atUpperBound);
02391 } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
02392 setStatus(iColumn,isFree);
02393 break;
02394 } else {
02395 break;
02396 }
02397 }
02398 }
02399 }
02400 firstFree_ = iColumn;
02401 if (firstFree_==numberRows_+numberColumns_)
02402 firstFree_=-1;
02403 return returnValue;
02404 } else {
02405 return -1;
02406 }
02407 }
02408 void
02409 ClpSimplexPrimal::clearAll()
02410 {
02411 int number=rowArray_[1]->getNumElements();
02412 int * which=rowArray_[1]->getIndices();
02413
02414 int iIndex;
02415 for (iIndex=0;iIndex<number;iIndex++) {
02416
02417 int iRow = which[iIndex];
02418 clearActive(iRow);
02419 }
02420 rowArray_[1]->clear();
02421 }
02422