Changeset 213 for branches/pre/ClpInterior.cpp
 Timestamp:
 Oct 2, 2003 2:06:08 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/pre/ClpInterior.cpp
r212 r213 11 11 #include "CoinHelperFunctions.hpp" 12 12 #include "ClpInterior.hpp" 13 #include "ClpFactorization.hpp"14 13 #include "ClpPackedMatrix.hpp" 15 14 #include "CoinIndexedVector.hpp" 16 #include "ClpDualRowDantzig.hpp"17 #include "ClpDualRowSteepest.hpp"18 #include "ClpPrimalColumnDantzig.hpp"19 #include "ClpPrimalColumnSteepest.hpp"20 #include "ClpNonLinearCost.hpp"21 15 #include "ClpMessage.hpp" 22 16 #include "ClpLinearObjective.hpp" … … 31 25 32 26 ClpModel(), 33 columnPrimalInfeasibility_(0.0),34 rowPrimalInfeasibility_(0.0),35 columnPrimalSequence_(2),36 rowPrimalSequence_(2),37 columnDualInfeasibility_(0.0),38 rowDualInfeasibility_(0.0),39 columnDualSequence_(2),40 rowDualSequence_(2),41 primalToleranceToGetOptimal_(1.0),42 remainingDualInfeasibility_(0.0),43 largeValue_(1.0e15),44 27 largestPrimalError_(0.0), 45 28 largestDualError_(0.0), 46 largestSolutionError_(0.0),47 dualBound_(1.0e10),48 alpha_(0.0),49 theta_(0.0),50 lowerIn_(0.0),51 valueIn_(0.0),52 upperIn_(0.0),53 dualIn_(0.0),54 lowerOut_(1),55 valueOut_(1),56 upperOut_(1),57 dualOut_(1),58 dualTolerance_(0.0),59 primalTolerance_(0.0),60 29 sumDualInfeasibilities_(0.0), 61 30 sumPrimalInfeasibilities_(0.0), 62 infeasibilityCost_(1.0e10),63 sumOfRelaxedDualInfeasibilities_(0.0),64 sumOfRelaxedPrimalInfeasibilities_(0.0),65 31 lower_(NULL), 66 32 rowLowerWork_(NULL), … … 69 35 rowUpperWork_(NULL), 70 36 columnUpperWork_(NULL), 71 cost_(NULL), 72 rowObjectiveWork_(NULL), 73 objectiveWork_(NULL), 74 sequenceIn_(1), 75 directionIn_(1), 76 sequenceOut_(1), 77 directionOut_(1), 78 pivotRow_(1), 79 lastGoodIteration_(100), 80 dj_(NULL), 81 rowReducedCost_(NULL), 82 reducedCostWork_(NULL), 83 solution_(NULL), 84 rowActivityWork_(NULL), 85 columnActivityWork_(NULL), 86 numberDualInfeasibilities_(0), 87 numberDualInfeasibilitiesWithoutFree_(0), 88 numberPrimalInfeasibilities_(100), 89 numberRefinements_(0), 90 pivotVariable_(NULL), 91 factorization_(NULL), 92 rowScale_(NULL), 93 savedSolution_(NULL), 94 columnScale_(NULL), 95 scalingFlag_(3), 96 numberTimesOptimal_(0), 97 changeMade_(1), 98 algorithm_(0), 99 forceFactorization_(1), 100 perturbation_(100), 101 nonLinearCost_(NULL), 102 specialOptions_(0), 103 lastBadIteration_(999999), 104 numberFake_(0), 105 progressFlag_(0), 106 firstFree_(1), 107 incomingInfeasibility_(1.0), 108 allowedInfeasibility_(10.0), 109 progress_(NULL) 110 { 111 int i; 112 for (i=0;i<6;i++) { 113 rowArray_[i]=NULL; 114 columnArray_[i]=NULL; 115 } 116 saveStatus_=NULL; 117 // get an empty factorization so we can set tolerances etc 118 factorization_ = new ClpFactorization(); 119 // Say sparse 120 factorization_>sparseThreshold(1); 121 // say Steepest pricing 122 dualRowPivot_ = new ClpDualRowSteepest(); 123 // say Steepest pricing 124 primalColumnPivot_ = new ClpPrimalColumnSteepest(); 125 solveType_=1; // say simplex based life form 126 37 cost_(NULL) 38 { 39 solveType_=2; // say interior based life form 127 40 } 128 41 … … 134 47 : ClpModel(rhs, numberRows, whichRow, 135 48 numberColumns,whichColumn,dropNames,dropIntegers), 136 columnPrimalInfeasibility_(0.0),137 rowPrimalInfeasibility_(0.0),138 columnPrimalSequence_(2),139 rowPrimalSequence_(2),140 columnDualInfeasibility_(0.0),141 rowDualInfeasibility_(0.0),142 columnDualSequence_(2),143 rowDualSequence_(2),144 primalToleranceToGetOptimal_(1.0),145 remainingDualInfeasibility_(0.0),146 largeValue_(1.0e15),147 49 largestPrimalError_(0.0), 148 50 largestDualError_(0.0), 149 largestSolutionError_(0.0),150 dualBound_(1.0e10),151 alpha_(0.0),152 theta_(0.0),153 lowerIn_(0.0),154 valueIn_(0.0),155 upperIn_(0.0),156 dualIn_(0.0),157 lowerOut_(1),158 valueOut_(1),159 upperOut_(1),160 dualOut_(1),161 dualTolerance_(0.0),162 primalTolerance_(0.0),163 51 sumDualInfeasibilities_(0.0), 164 52 sumPrimalInfeasibilities_(0.0), 165 infeasibilityCost_(1.0e10),166 sumOfRelaxedDualInfeasibilities_(0.0),167 sumOfRelaxedPrimalInfeasibilities_(0.0),168 53 lower_(NULL), 169 54 rowLowerWork_(NULL), … … 172 57 rowUpperWork_(NULL), 173 58 columnUpperWork_(NULL), 174 cost_(NULL), 175 rowObjectiveWork_(NULL), 176 objectiveWork_(NULL), 177 sequenceIn_(1), 178 directionIn_(1), 179 sequenceOut_(1), 180 directionOut_(1), 181 pivotRow_(1), 182 lastGoodIteration_(100), 183 dj_(NULL), 184 rowReducedCost_(NULL), 185 reducedCostWork_(NULL), 186 solution_(NULL), 187 rowActivityWork_(NULL), 188 columnActivityWork_(NULL), 189 numberDualInfeasibilities_(0), 190 numberDualInfeasibilitiesWithoutFree_(0), 191 numberPrimalInfeasibilities_(100), 192 numberRefinements_(0), 193 pivotVariable_(NULL), 194 factorization_(NULL), 195 rowScale_(NULL), 196 savedSolution_(NULL), 197 columnScale_(NULL), 198 scalingFlag_(3), 199 numberTimesOptimal_(0), 200 changeMade_(1), 201 algorithm_(0), 202 forceFactorization_(1), 203 perturbation_(100), 204 nonLinearCost_(NULL), 205 specialOptions_(0), 206 lastBadIteration_(999999), 207 numberFake_(0), 208 progressFlag_(0), 209 firstFree_(1), 210 incomingInfeasibility_(1.0), 211 allowedInfeasibility_(10.0), 212 progress_(NULL) 213 { 214 int i; 215 for (i=0;i<6;i++) { 216 rowArray_[i]=NULL; 217 columnArray_[i]=NULL; 218 } 219 saveStatus_=NULL; 220 // get an empty factorization so we can set tolerances etc 221 factorization_ = new ClpFactorization(); 222 // say Steepest pricing 223 dualRowPivot_ = new ClpDualRowSteepest(); 224 // say Steepest pricing 225 primalColumnPivot_ = new ClpPrimalColumnSteepest(); 226 solveType_=1; // say simplex based life form 59 cost_(NULL) 60 { 61 solveType_=2; // say interior based life form 227 62 228 63 } … … 232 67 ClpInterior::~ClpInterior () 233 68 { 234 gutsOfDelete(0); 235 delete nonLinearCost_; 69 gutsOfDelete(); 236 70 } 237 71 //############################################################################# 238 void ClpInterior::setLargeValue( double value)239 {240 if (value>0.0&&value<COIN_DBL_MAX)241 largeValue_=value;242 }243 int244 ClpInterior::gutsOfSolution ( double * givenDuals,245 const double * givenPrimals,246 bool valuesPass)247 {248 249 250 // if values pass, save values of basic variables251 double * save = NULL;252 double oldValue=0.0;253 if (valuesPass) {254 assert(algorithm_>0); // only primal at present255 assert(nonLinearCost_);256 int iRow;257 checkPrimalSolution( rowActivityWork_, columnActivityWork_);258 // get correct bounds on all variables259 nonLinearCost_>checkInfeasibilities(primalTolerance_);260 oldValue = nonLinearCost_>largestInfeasibility();261 save = new double[numberRows_];262 for (iRow=0;iRow<numberRows_;iRow++) {263 int iPivot=pivotVariable_[iRow];264 save[iRow] = solution_[iPivot];265 }266 }267 // do work268 computePrimals(rowActivityWork_, columnActivityWork_);269 // If necessary  override results270 if (givenPrimals) {271 memcpy(columnActivityWork_,givenPrimals,numberColumns_*sizeof(double));272 memset(rowActivityWork_,0,numberRows_*sizeof(double));273 times(1.0,columnActivityWork_,rowActivityWork_);274 }275 double objectiveModification = 0.0;276 if (algorithm_>0&&nonLinearCost_!=NULL) {277 // primal algorithm278 // get correct bounds on all variables279 // If 4 bit set  Force outgoing variables to exact bound (primal)280 if ((specialOptions_&4)==0)281 nonLinearCost_>checkInfeasibilities(primalTolerance_);282 else283 nonLinearCost_>checkInfeasibilities(0.0);284 objectiveModification += nonLinearCost_>changeInCost();285 if (nonLinearCost_>numberInfeasibilities())286 handler_>message(CLP_SIMPLEX_NONLINEAR,messages_)287 <<nonLinearCost_>changeInCost()288 <<nonLinearCost_>numberInfeasibilities()289 <<CoinMessageEol;290 }291 if (valuesPass) {292 #ifdef CLP_DEBUG293 std::cout<<"Largest given infeasibility "<<oldValue294 <<" now "<<nonLinearCost_>largestInfeasibility()<<std::endl;295 #endif296 int numberOut=0;297 if (oldValue<incomingInfeasibility_298 &&nonLinearCost_>largestInfeasibility()>299 max(incomingInfeasibility_,allowedInfeasibility_)300 largestPrimalError_>1.0e3) {301 printf("Original largest infeas %g, now %g, primalError %g\n",302 oldValue,nonLinearCost_>largestInfeasibility(),303 largestPrimalError_);304 // throw out up to 1000 structurals305 int iRow;306 int * sort = new int[numberRows_];307 // first put back solution and store difference308 for (iRow=0;iRow<numberRows_;iRow++) {309 int iPivot=pivotVariable_[iRow];310 double difference = fabs(solution_[iPivot]save[iRow]);311 solution_[iPivot]=save[iRow];312 save[iRow]=difference;313 }314 for (iRow=0;iRow<numberRows_;iRow++) {315 int iPivot=pivotVariable_[iRow];316 317 if (iPivot<numberColumns_) {318 // column319 double difference= save[iRow];320 if (difference>1.0e4) {321 sort[numberOut]=iPivot;322 save[numberOut++]=difference;323 }324 }325 }326 CoinSort_2(save, save + numberOut, sort,327 CoinFirstGreater_2<double, int>());328 numberOut = min(1000,numberOut);329 for (iRow=0;iRow<numberOut;iRow++) {330 int iColumn=sort[iRow];331 setColumnStatus(iColumn,superBasic);332 333 }334 delete [] sort;335 }336 delete [] save;337 if (numberOut)338 return numberOut;339 }340 341 computeDuals(givenDuals);342 343 // now check solutions344 checkPrimalSolution( rowActivityWork_, columnActivityWork_);345 objectiveValue_ += objectiveModification;346 checkDualSolution();347 if (handler_>logLevel()>3(largestPrimalError_>1.0e2348 largestDualError_>1.0e2))349 handler_>message(CLP_SIMPLEX_ACCURACY,messages_)350 <<largestPrimalError_351 <<largestDualError_352 <<CoinMessageEol;353 // Switch off false values pass indicator354 if (!valuesPass&&algorithm_>0)355 firstFree_ = 1;356 return 0;357 }358 void359 ClpInterior::computePrimals ( const double * rowActivities,360 const double * columnActivities)361 {362 363 //work space364 CoinIndexedVector * workSpace = rowArray_[0];365 366 CoinIndexedVector arrayVector;367 arrayVector.reserve(numberRows_+1);368 CoinIndexedVector previousVector;369 previousVector.reserve(numberRows_+1);370 371 // accumulate non basic stuff372 373 int iRow;374 // order is this way for scaling375 // Use whole matrix every time to make it easier for ClpMatrixBase376 // So zero out basic377 if (columnActivities!=columnActivityWork_)378 ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);379 if (rowActivities!=rowActivityWork_)380 ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);381 for (iRow=0;iRow<numberRows_;iRow++) {382 int iPivot=pivotVariable_[iRow];383 solution_[iPivot] = 0.0;384 }385 double * array = arrayVector.denseVector();386 times(1.0,columnActivityWork_,array);387 388 int * index = arrayVector.getIndices();389 int number=0;390 for (iRow=0;iRow<numberRows_;iRow++) {391 double value = array[iRow] + rowActivityWork_[iRow];392 if (value) {393 array[iRow]=value;394 index[number++]=iRow;395 } else {396 array[iRow]=0.0;397 }398 }399 arrayVector.setNumElements(number);400 401 // Ftran adjusted RHS and iterate to improve accuracy402 double lastError=COIN_DBL_MAX;403 int iRefine;404 double * work = workSpace>denseVector();405 CoinIndexedVector * thisVector = &arrayVector;406 CoinIndexedVector * lastVector = &previousVector;407 factorization_>updateColumn(workSpace,thisVector);408 bool goodSolution=true;409 for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {410 411 int numberIn = thisVector>getNumElements();412 int * indexIn = thisVector>getIndices();413 double * arrayIn = thisVector>denseVector();414 // put solution in correct place415 int j;416 for (j=0;j<numberIn;j++) {417 iRow = indexIn[j];418 int iPivot=pivotVariable_[iRow];419 solution_[iPivot] = arrayIn[iRow];420 }421 // check Ax == b (for all)422 times(1.0,columnActivityWork_,work);423 largestPrimalError_=0.0;424 double multiplier = 131072.0;425 for (iRow=0;iRow<numberRows_;iRow++) {426 double value = work[iRow] + rowActivityWork_[iRow];427 work[iRow] = value*multiplier;428 if (fabs(value)>largestPrimalError_) {429 largestPrimalError_=fabs(value);430 }431 }432 if (largestPrimalError_>=lastError) {433 // restore434 CoinIndexedVector * temp = thisVector;435 thisVector = lastVector;436 lastVector=temp;437 goodSolution=false;438 break;439 }440 if (iRefine<numberRefinements_&&largestPrimalError_>1.0e10) {441 // try and make better442 // save this443 CoinIndexedVector * temp = thisVector;444 thisVector = lastVector;445 lastVector=temp;446 int * indexOut = thisVector>getIndices();447 int number=0;448 array = thisVector>denseVector();449 thisVector>clear();450 for (iRow=0;iRow<numberRows_;iRow++) {451 double value = work[iRow];452 if (value) {453 array[iRow]=value;454 indexOut[number++]=iRow;455 work[iRow]=0.0;456 }457 }458 thisVector>setNumElements(number);459 lastError=largestPrimalError_;460 factorization_>updateColumn(workSpace,thisVector);461 multiplier = 1.0/multiplier;462 double * previous = lastVector>denseVector();463 number=0;464 for (iRow=0;iRow<numberRows_;iRow++) {465 double value = previous[iRow] + multiplier*array[iRow];466 if (value) {467 array[iRow]=value;468 indexOut[number++]=iRow;469 } else {470 array[iRow]=0.0;471 }472 }473 thisVector>setNumElements(number);474 } else {475 break;476 }477 }478 479 // solution as accurate as we are going to get480 ClpFillN(work,numberRows_,0.0);481 if (!goodSolution) {482 array = thisVector>denseVector();483 // put solution in correct place484 for (iRow=0;iRow<numberRows_;iRow++) {485 int iPivot=pivotVariable_[iRow];486 solution_[iPivot] = array[iRow];487 }488 }489 }490 // now dual side491 void492 ClpInterior::computeDuals(double * givenDjs)493 {494 //work space495 CoinIndexedVector * workSpace = rowArray_[0];496 497 CoinIndexedVector arrayVector;498 arrayVector.reserve(numberRows_+1);499 CoinIndexedVector previousVector;500 previousVector.reserve(numberRows_+1);501 502 503 int iRow;504 #ifdef CLP_DEBUG505 workSpace>checkClear();506 #endif507 double * array = arrayVector.denseVector();508 int * index = arrayVector.getIndices();509 int number=0;510 if (!givenDjs) {511 for (iRow=0;iRow<numberRows_;iRow++) {512 int iPivot=pivotVariable_[iRow];513 double value = cost_[iPivot];514 if (value) {515 array[iRow]=value;516 index[number++]=iRow;517 }518 }519 } else {520 // dual values pass  djs may not be zero521 for (iRow=0;iRow<numberRows_;iRow++) {522 int iPivot=pivotVariable_[iRow];523 // make sure zero if done524 if (!pivoted(iPivot))525 givenDjs[iPivot]=0.0;526 double value =cost_[iPivot]givenDjs[iPivot];527 if (value) {528 array[iRow]=value;529 index[number++]=iRow;530 }531 }532 }533 arrayVector.setNumElements(number);534 535 // Btran basic costs and get as accurate as possible536 double lastError=COIN_DBL_MAX;537 int iRefine;538 double * work = workSpace>denseVector();539 CoinIndexedVector * thisVector = &arrayVector;540 CoinIndexedVector * lastVector = &previousVector;541 factorization_>updateColumnTranspose(workSpace,thisVector);542 543 for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {544 // check basic reduced costs zero545 largestDualError_=0.0;546 // would be faster to do just for basic but this reduces code547 ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);548 transposeTimes(1.0,array,reducedCostWork_);549 if (!givenDjs) {550 for (iRow=0;iRow<numberRows_;iRow++) {551 int iPivot=pivotVariable_[iRow];552 double value;553 if (iPivot>=numberColumns_) {554 // slack555 value = rowObjectiveWork_[iPivotnumberColumns_]556 + array[iPivotnumberColumns_];557 } else {558 // column559 value = reducedCostWork_[iPivot];560 }561 work[iRow]=value;562 if (fabs(value)>largestDualError_) {563 largestDualError_=fabs(value);564 }565 }566 } else {567 for (iRow=0;iRow<numberRows_;iRow++) {568 int iPivot=pivotVariable_[iRow];569 if (iPivot>=numberColumns_) {570 // slack571 work[iRow] = rowObjectiveWork_[iPivotnumberColumns_]572 + array[iPivotnumberColumns_]givenDjs[iPivot];573 } else {574 // column575 work[iRow] = reducedCostWork_[iPivot] givenDjs[iPivot];576 }577 if (fabs(work[iRow])>largestDualError_) {578 largestDualError_=fabs(work[iRow]);579 //assert (largestDualError_<1.0e7);580 //if (largestDualError_>1.0e7)581 //printf("large dual error %g\n",largestDualError_);582 }583 }584 }585 if (largestDualError_>=lastError) {586 // restore587 CoinIndexedVector * temp = thisVector;588 thisVector = lastVector;589 lastVector=temp;590 break;591 }592 if (iRefine<numberRefinements_&&largestDualError_>1.0e10593 &&!givenDjs) {594 // try and make better595 // save this596 CoinIndexedVector * temp = thisVector;597 thisVector = lastVector;598 lastVector=temp;599 int * indexOut = thisVector>getIndices();600 int number=0;601 array = thisVector>denseVector();602 thisVector>clear();603 double multiplier = 131072.0;604 for (iRow=0;iRow<numberRows_;iRow++) {605 double value = multiplier*work[iRow];606 if (value) {607 array[iRow]=value;608 indexOut[number++]=iRow;609 work[iRow]=0.0;610 }611 work[iRow]=0.0;612 }613 thisVector>setNumElements(number);614 lastError=largestDualError_;615 factorization_>updateColumnTranspose(workSpace,thisVector);616 multiplier = 1.0/multiplier;617 double * previous = lastVector>denseVector();618 number=0;619 for (iRow=0;iRow<numberRows_;iRow++) {620 double value = previous[iRow] + multiplier*array[iRow];621 if (value) {622 array[iRow]=value;623 indexOut[number++]=iRow;624 } else {625 array[iRow]=0.0;626 }627 }628 thisVector>setNumElements(number);629 } else {630 break;631 }632 }633 ClpFillN(work,numberRows_,0.0);634 // now look at dual solution635 array = thisVector>denseVector();636 for (iRow=0;iRow<numberRows_;iRow++) {637 // slack638 double value = array[iRow];639 dual_[iRow]=value;640 value += rowObjectiveWork_[iRow];641 rowReducedCost_[iRow]=value;642 }643 ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);644 transposeTimes(1.0,dual_,reducedCostWork_);645 // If necessary  override results646 if (givenDjs) {647 // restore accurate duals648 memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));649 }650 651 }652 /* Given an existing factorization computes and checks653 primal and dual solutions. Uses input arrays for variables at654 bounds. Returns feasibility states */655 int ClpInterior::getSolution ( const double * rowActivities,656 const double * columnActivities)657 {658 if (!factorization_>status()) {659 // put in standard form660 createRim(7+8+16+32);661 // do work662 gutsOfSolution ( NULL,NULL);663 // release extra memory664 deleteRim(0);665 }666 return factorization_>status();667 }668 /* Given an existing factorization computes and checks669 primal and dual solutions. Uses current problem arrays for670 bounds. Returns feasibility states */671 int ClpInterior::getSolution ( )672 {673 double * rowActivities = new double[numberRows_];674 double * columnActivities = new double[numberColumns_];675 ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);676 ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);677 int status = getSolution( rowActivities, columnActivities);678 delete [] rowActivities;679 delete [] columnActivities;680 return status;681 }682 // Factorizes using current basis. This is for external use683 // Return codes are as from ClpFactorization684 int ClpInterior::factorize ()685 {686 // put in standard form687 createRim(7+8+16+32,false);688 // do work689 int status = internalFactorize(1);690 // release extra memory691 deleteRim(0);692 693 return status;694 }695 // Clean up status696 void697 ClpInterior::cleanStatus()698 {699 int iRow,iColumn;700 int numberBasic=0;701 // make row activities correct702 memset(rowActivityWork_,0,numberRows_*sizeof(double));703 times(1.0,columnActivityWork_,rowActivityWork_);704 if (!status_)705 createStatus();706 for (iRow=0;iRow<numberRows_;iRow++) {707 if (getRowStatus(iRow)==basic)708 numberBasic++;709 else {710 setRowStatus(iRow,superBasic);711 // but put to bound if close712 if (fabs(rowActivityWork_[iRow]rowLowerWork_[iRow])713 <=primalTolerance_) {714 rowActivityWork_[iRow]=rowLowerWork_[iRow];715 setRowStatus(iRow,atLowerBound);716 } else if (fabs(rowActivityWork_[iRow]rowUpperWork_[iRow])717 <=primalTolerance_) {718 rowActivityWork_[iRow]=rowUpperWork_[iRow];719 setRowStatus(iRow,atUpperBound);720 }721 }722 }723 for (iColumn=0;iColumn<numberColumns_;iColumn++) {724 if (getColumnStatus(iColumn)==basic) {725 if (numberBasic==numberRows_) {726 // take out of basis727 setColumnStatus(iColumn,superBasic);728 // but put to bound if close729 if (fabs(columnActivityWork_[iColumn]columnLowerWork_[iColumn])730 <=primalTolerance_) {731 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];732 setColumnStatus(iColumn,atLowerBound);733 } else if (fabs(columnActivityWork_[iColumn]734 columnUpperWork_[iColumn])735 <=primalTolerance_) {736 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];737 setColumnStatus(iColumn,atUpperBound);738 }739 } else740 numberBasic++;741 } else {742 setColumnStatus(iColumn,superBasic);743 // but put to bound if close744 if (fabs(columnActivityWork_[iColumn]columnLowerWork_[iColumn])745 <=primalTolerance_) {746 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];747 setColumnStatus(iColumn,atLowerBound);748 } else if (fabs(columnActivityWork_[iColumn]749 columnUpperWork_[iColumn])750 <=primalTolerance_) {751 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];752 setColumnStatus(iColumn,atUpperBound);753 }754 }755 }756 }757 758 /* Factorizes using current basis.759 solveType  1 iterating, 0 initial, 1 external760  2 then iterating but can throw out of basis761 If 10 added then in primal values pass762 */763 /* Return codes are as from ClpFactorization unless initial factorization764 when total number of singularities is returned */765 int ClpInterior::internalFactorize ( int solveType)766 {767 768 int iRow,iColumn;769 int totalSlacks=numberRows_;770 if (!status_)771 createStatus();772 773 bool valuesPass=false;774 if (solveType>=10) {775 valuesPass=true;776 solveType = 10;777 }778 #ifdef CLP_DEBUG779 if (solveType>0) {780 int numberFreeIn=0,numberFreeOut=0;781 double biggestDj=0.0;782 for (iColumn=0;iColumn<numberColumns_;iColumn++) {783 switch(getColumnStatus(iColumn)) {784 785 case basic:786 if (columnLower_[iColumn]<largeValue_787 &&columnUpper_[iColumn]>largeValue_)788 numberFreeIn++;789 break;790 default:791 if (columnLower_[iColumn]<largeValue_792 &&columnUpper_[iColumn]>largeValue_) {793 numberFreeOut++;794 biggestDj = max(fabs(dj_[iColumn]),biggestDj);795 }796 break;797 }798 }799 if (numberFreeIn+numberFreeOut)800 printf("%d in basis, %d out  largest dj %g\n",801 numberFreeIn,numberFreeOut,biggestDj);802 }803 #endif804 if (solveType<=0) {805 // Make sure everything is clean806 for (iRow=0;iRow<numberRows_;iRow++) {807 if(getRowStatus(iRow)==isFixed) {808 // double check fixed809 if (rowUpperWork_[iRow]>rowLowerWork_[iRow])810 setRowStatus(iRow,atLowerBound);811 } else if (getRowStatus(iRow)==isFree) {812 // may not be free after all813 if (rowLowerWork_[iRow]>largeValue_rowUpperWork_[iRow]<largeValue_)814 setRowStatus(iRow,superBasic);815 }816 }817 for (iColumn=0;iColumn<numberColumns_;iColumn++) {818 if(getColumnStatus(iColumn)==isFixed) {819 // double check fixed820 if (columnUpperWork_[iColumn]>columnLowerWork_[iColumn])821 setColumnStatus(iColumn,atLowerBound);822 } else if (getColumnStatus(iColumn)==isFree) {823 // may not be free after all824 if (columnLowerWork_[iColumn]>largeValue_columnUpperWork_[iColumn]<largeValue_)825 setColumnStatus(iColumn,superBasic);826 }827 }828 if (!valuesPass) {829 // not values pass so set to bounds830 bool allSlack=true;831 if (status_) {832 for (iRow=0;iRow<numberRows_;iRow++) {833 if (getRowStatus(iRow)!=basic) {834 allSlack=false;835 break;836 }837 }838 }839 if (!allSlack) {840 // set values from warm start (if sensible)841 int numberBasic=0;842 for (iRow=0;iRow<numberRows_;iRow++) {843 switch(getRowStatus(iRow)) {844 845 case basic:846 numberBasic++;847 break;848 case atUpperBound:849 rowActivityWork_[iRow]=rowUpperWork_[iRow];850 if (rowActivityWork_[iRow]>largeValue_) {851 if (rowLowerWork_[iRow]>largeValue_) {852 rowActivityWork_[iRow]=rowLowerWork_[iRow];853 setRowStatus(iRow,atLowerBound);854 } else {855 // say free856 setRowStatus(iRow,isFree);857 rowActivityWork_[iRow]=0.0;858 }859 }860 break;861 case ClpInterior::isFixed:862 case atLowerBound:863 rowActivityWork_[iRow]=rowLowerWork_[iRow];864 if (rowActivityWork_[iRow]<largeValue_) {865 if (rowUpperWork_[iRow]<largeValue_) {866 rowActivityWork_[iRow]=rowUpperWork_[iRow];867 setRowStatus(iRow,atUpperBound);868 } else {869 // say free870 setRowStatus(iRow,isFree);871 rowActivityWork_[iRow]=0.0;872 }873 }874 break;875 case isFree:876 break;877 // not really free  fall through to superbasic878 case superBasic:879 if (rowUpperWork_[iRow]>largeValue_) {880 if (rowLowerWork_[iRow]>largeValue_) {881 rowActivityWork_[iRow]=rowLowerWork_[iRow];882 setRowStatus(iRow,atLowerBound);883 } else {884 // say free885 setRowStatus(iRow,isFree);886 rowActivityWork_[iRow]=0.0;887 }888 } else {889 if (rowLowerWork_[iRow]>largeValue_) {890 // set to nearest891 if (fabs(rowActivityWork_[iRow]rowLowerWork_[iRow])892 <fabs(rowActivityWork_[iRow]rowLowerWork_[iRow])) {893 rowActivityWork_[iRow]=rowLowerWork_[iRow];894 setRowStatus(iRow,atLowerBound);895 } else {896 rowActivityWork_[iRow]=rowUpperWork_[iRow];897 setRowStatus(iRow,atUpperBound);898 }899 } else {900 rowActivityWork_[iRow]=rowUpperWork_[iRow];901 setRowStatus(iRow,atUpperBound);902 }903 }904 break;905 }906 }907 totalSlacks=numberBasic;908 909 for (iColumn=0;iColumn<numberColumns_;iColumn++) {910 switch(getColumnStatus(iColumn)) {911 912 case basic:913 if (numberBasic==numberRows_) {914 // take out of basis915 if (columnLowerWork_[iColumn]>largeValue_) {916 if (columnActivityWork_[iColumn]columnLowerWork_[iColumn]<917 columnUpperWork_[iColumn]columnActivityWork_[iColumn]) {918 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];919 setColumnStatus(iColumn,atLowerBound);920 } else {921 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];922 setColumnStatus(iColumn,atUpperBound);923 }924 } else if (columnUpperWork_[iColumn]<largeValue_) {925 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];926 setColumnStatus(iColumn,atUpperBound);927 } else {928 columnActivityWork_[iColumn]=0.0;929 setColumnStatus(iColumn,isFree);930 }931 } else {932 numberBasic++;933 }934 break;935 case atUpperBound:936 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];937 if (columnActivityWork_[iColumn]>largeValue_) {938 if (columnLowerWork_[iColumn]<largeValue_) {939 columnActivityWork_[iColumn]=0.0;940 setColumnStatus(iColumn,isFree);941 } else {942 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];943 setColumnStatus(iColumn,atLowerBound);944 }945 }946 break;947 case isFixed:948 case atLowerBound:949 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];950 if (columnActivityWork_[iColumn]<largeValue_) {951 if (columnUpperWork_[iColumn]>largeValue_) {952 columnActivityWork_[iColumn]=0.0;953 setColumnStatus(iColumn,isFree);954 } else {955 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];956 setColumnStatus(iColumn,atUpperBound);957 }958 }959 break;960 case isFree:961 break;962 // not really free  fall through to superbasic963 case superBasic:964 if (columnUpperWork_[iColumn]>largeValue_) {965 if (columnLowerWork_[iColumn]>largeValue_) {966 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];967 setColumnStatus(iColumn,atLowerBound);968 } else {969 // say free970 setColumnStatus(iColumn,isFree);971 columnActivityWork_[iColumn]=0.0;972 }973 } else {974 if (columnLowerWork_[iColumn]>largeValue_) {975 // set to nearest976 if (fabs(columnActivityWork_[iColumn]columnLowerWork_[iColumn])977 <fabs(columnActivityWork_[iColumn]columnLowerWork_[iColumn])) {978 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];979 setColumnStatus(iColumn,atLowerBound);980 } else {981 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];982 setColumnStatus(iColumn,atUpperBound);983 }984 } else {985 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];986 setColumnStatus(iColumn,atUpperBound);987 }988 }989 break;990 }991 }992 } else {993 // all slack basis994 int numberBasic=0;995 if (!status_) {996 createStatus();997 }998 for (iRow=0;iRow<numberRows_;iRow++) {999 double lower=rowLowerWork_[iRow];1000 double upper=rowUpperWork_[iRow];1001 if (lower>largeValue_upper<largeValue_) {1002 if (fabs(lower)<=fabs(upper)) {1003 rowActivityWork_[iRow]=lower;1004 } else {1005 rowActivityWork_[iRow]=upper;1006 }1007 } else {1008 rowActivityWork_[iRow]=0.0;1009 }1010 setRowStatus(iRow,basic);1011 numberBasic++;1012 }1013 for (iColumn=0;iColumn<numberColumns_;iColumn++) {1014 double lower=columnLowerWork_[iColumn];1015 double upper=columnUpperWork_[iColumn];1016 double big_bound = largeValue_;1017 if (lower>big_boundupper<big_bound) {1018 if ((getColumnStatus(iColumn)==atLowerBound&&1019 columnActivityWork_[iColumn]==lower)1020 (getColumnStatus(iColumn)==atUpperBound&&1021 columnActivityWork_[iColumn]==upper)) {1022 // status looks plausible1023 } else {1024 // set to sensible1025 if (fabs(lower)<=fabs(upper)) {1026 setColumnStatus(iColumn,atLowerBound);1027 columnActivityWork_[iColumn]=lower;1028 } else {1029 setColumnStatus(iColumn,atUpperBound);1030 columnActivityWork_[iColumn]=upper;1031 }1032 }1033 } else {1034 setColumnStatus(iColumn,isFree);1035 columnActivityWork_[iColumn]=0.0;1036 }1037 }1038 }1039 } else {1040 // values pass has less coding1041 // make row activities correct and clean basis a bit1042 cleanStatus();1043 if (status_) {1044 int numberBasic=0;1045 for (iRow=0;iRow<numberRows_;iRow++) {1046 if (getRowStatus(iRow)==basic)1047 numberBasic++;1048 }1049 totalSlacks=numberBasic;1050 for (iColumn=0;iColumn<numberColumns_;iColumn++) {1051 if (getColumnStatus(iColumn)==basic)1052 numberBasic++;1053 }1054 } else {1055 // all slack basis1056 int numberBasic=0;1057 if (!status_) {1058 createStatus();1059 }1060 for (iRow=0;iRow<numberRows_;iRow++) {1061 setRowStatus(iRow,basic);1062 numberBasic++;1063 }1064 for (iColumn=0;iColumn<numberColumns_;iColumn++) {1065 setColumnStatus(iColumn,superBasic);1066 // but put to bound if close1067 if (fabs(columnActivityWork_[iColumn]columnLowerWork_[iColumn])1068 <=primalTolerance_) {1069 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];1070 setColumnStatus(iColumn,atLowerBound);1071 } else if (fabs(columnActivityWork_[iColumn]1072 columnUpperWork_[iColumn])1073 <=primalTolerance_) {1074 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];1075 setColumnStatus(iColumn,atUpperBound);1076 }1077 }1078 }1079 }1080 numberRefinements_=1;1081 // set fixed if they are1082 for (iRow=0;iRow<numberRows_;iRow++) {1083 if (getRowStatus(iRow)!=basic ) {1084 if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {1085 rowActivityWork_[iRow]=rowLowerWork_[iRow];1086 setRowStatus(iRow,isFixed);1087 }1088 }1089 }1090 for (iColumn=0;iColumn<numberColumns_;iColumn++) {1091 if (getColumnStatus(iColumn)!=basic ) {1092 if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {1093 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];1094 setColumnStatus(iColumn,isFixed);1095 }1096 }1097 }1098 }1099 if (0) {1100 static int k=0;1101 printf("start basis\n");1102 int i;1103 for (i=0;i<numberRows_;i++)1104 printf ("xx %d %d\n",i,pivotVariable_[i]);1105 for (i=0;i<numberRows_+numberColumns_;i++)1106 if (getColumnStatus(i)==basic)1107 printf ("yy %d basic\n",i);1108 if (k>20)1109 exit(0);1110 k++;1111 }1112 int status = factorization_>factorize(this, solveType,valuesPass);1113 if (status) {1114 handler_>message(CLP_SIMPLEX_BADFACTOR,messages_)1115 <<status1116 <<CoinMessageEol;1117 return 1;1118 } else if (!solveType) {1119 // Initial basis  return number of singularities1120 int numberSlacks=0;1121 for (iRow=0;iRow<numberRows_;iRow++) {1122 if (getRowStatus(iRow) == basic)1123 numberSlacks++;1124 }1125 status= max(numberSlackstotalSlacks,0);1126 }1127 1128 // sparse methods1129 //if (factorization_>sparseThreshold()) {1130 // get default value1131 factorization_>sparseThreshold(0);1132 factorization_>goSparse();1133 //}1134 1135 return status;1136 }1137 72 /* 1138 This does basis housekeeping and does values for in/out variables. 1139 Can also decide to refactorize 73 This does housekeeping 1140 74 */ 1141 75 int 1142 ClpInterior::housekeeping( double objectiveChange)76 ClpInterior::housekeeping() 1143 77 { 1144 78 numberIterations_++; 1145 changeMade_++; // something has happened 1146 // incoming variable 1147 handler_>message(CLP_SIMPLEX_HOUSE1,messages_) 1148 <<directionOut_ 1149 <<directionIn_<<theta_ 1150 <<dualOut_<<dualIn_<<alpha_ 1151 <<CoinMessageEol; 1152 if (getStatus(sequenceIn_)==isFree) { 1153 handler_>message(CLP_SIMPLEX_FREEIN,messages_) 1154 <<sequenceIn_ 1155 <<CoinMessageEol; 1156 } 1157 // change of incoming 1158 char rowcol[]={'R','C'}; 1159 if (pivotRow_>=0) 1160 pivotVariable_[pivotRow_]=sequenceIn(); 1161 if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<1.0e20) 1162 progressFlag_ = 2; // making real progress 1163 solution_[sequenceIn_]=valueIn_; 1164 if (upper_[sequenceOut_]lower_[sequenceOut_]<1.0e12) 1165 progressFlag_ = 1; // making real progress 1166 if (sequenceIn_!=sequenceOut_) { 1167 //assert( getStatus(sequenceOut_)== basic); 1168 setStatus(sequenceIn_,basic); 1169 if (upper_[sequenceOut_]lower_[sequenceOut_]>0) { 1170 // As Nonlinear costs may have moved bounds (to more feasible) 1171 // Redo using value 1172 if (fabs(valueOut_lower_[sequenceOut_])<fabs(valueOut_upper_[sequenceOut_])) { 1173 // going to lower 1174 setStatus(sequenceOut_,atLowerBound); 1175 } else { 1176 // going to upper 1177 setStatus(sequenceOut_,atUpperBound); 1178 } 1179 } else { 1180 // fixed 1181 setStatus(sequenceOut_,isFixed); 1182 } 1183 solution_[sequenceOut_]=valueOut_; 1184 } else { 1185 // flip from bound to bound 1186 // As Nonlinear costs may have moved bounds (to more feasible) 1187 // Redo using value 1188 if (fabs(valueIn_lower_[sequenceIn_])<fabs(valueIn_upper_[sequenceIn_])) { 1189 // as if from upper bound 1190 setStatus(sequenceIn_, atLowerBound); 1191 } else { 1192 // as if from lower bound 1193 setStatus(sequenceIn_, atUpperBound); 1194 } 1195 } 1196 objectiveValue_ += objectiveChange; 1197 handler_>message(CLP_SIMPLEX_HOUSE2,messages_) 1198 <<numberIterations_<<objectiveValue() 1199 <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_) 1200 <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_); 1201 handler_>printing(algorithm_<0)<<theta_<<dualOut_; 1202 handler_>printing(algorithm_>0)<<dualIn_<<theta_; 1203 handler_>message()<<CoinMessageEol; 1204 if (hitMaximumIterations()) 1205 return 2; 1206 #if 0 1207 // check for small cycles 1208 int cycle=progress_>cycle(sequenceIn_,sequenceOut_, 1209 directionIn_,directionOut_); 1210 if (cycle>0) { 1211 printf("Cycle of %d\n",cycle); 1212 if (factorization_>pivots()>cycle) 1213 forceFactorization_=cycle1; 1214 else 1215 forceFactorization_=1; 1216 return 1; 1217 } 1218 #endif 1219 // only time to refactorize if one before real time 1220 // this is so user won't be surprised that maximumPivots has exact meaning 1221 if (factorization_>pivots()==factorization_>maximumPivots()) { 1222 return 1; 1223 } else { 1224 if (forceFactorization_>0&& 1225 factorization_>pivots()==forceFactorization_) { 1226 // relax 1227 forceFactorization_ = (3+5*forceFactorization_)/4; 1228 if (forceFactorization_>factorization_>maximumPivots()) 1229 forceFactorization_ = 1; //off 1230 return 1; 1231 } else { 1232 // carry on iterating 1233 return 0; 1234 } 1235 } 79 return 0; 1236 80 } 1237 81 // Copy constructor. 1238 82 ClpInterior::ClpInterior(const ClpInterior &rhs) : 1239 83 ClpModel(rhs), 1240 columnPrimalInfeasibility_(0.0),1241 rowPrimalInfeasibility_(0.0),1242 columnPrimalSequence_(2),1243 rowPrimalSequence_(2),1244 columnDualInfeasibility_(0.0),1245 rowDualInfeasibility_(0.0),1246 columnDualSequence_(2),1247 rowDualSequence_(2),1248 primalToleranceToGetOptimal_(1.0),1249 remainingDualInfeasibility_(0.0),1250 largeValue_(1.0e15),1251 84 largestPrimalError_(0.0), 1252 85 largestDualError_(0.0), 1253 largestSolutionError_(0.0),1254 dualBound_(1.0e10),1255 alpha_(0.0),1256 theta_(0.0),1257 lowerIn_(0.0),1258 valueIn_(0.0),1259 upperIn_(0.0),1260 dualIn_(0.0),1261 lowerOut_(1),1262 valueOut_(1),1263 upperOut_(1),1264 dualOut_(1),1265 dualTolerance_(0.0),1266 primalTolerance_(0.0),1267 86 sumDualInfeasibilities_(0.0), 1268 87 sumPrimalInfeasibilities_(0.0), 1269 infeasibilityCost_(1.0e10),1270 sumOfRelaxedDualInfeasibilities_(0.0),1271 sumOfRelaxedPrimalInfeasibilities_(0.0),1272 88 lower_(NULL), 1273 89 rowLowerWork_(NULL), … … 1276 92 rowUpperWork_(NULL), 1277 93 columnUpperWork_(NULL), 1278 cost_(NULL), 1279 rowObjectiveWork_(NULL), 1280 objectiveWork_(NULL), 1281 sequenceIn_(1), 1282 directionIn_(1), 1283 sequenceOut_(1), 1284 directionOut_(1), 1285 pivotRow_(1), 1286 lastGoodIteration_(100), 1287 dj_(NULL), 1288 rowReducedCost_(NULL), 1289 reducedCostWork_(NULL), 1290 solution_(NULL), 1291 rowActivityWork_(NULL), 1292 columnActivityWork_(NULL), 1293 numberDualInfeasibilities_(0), 1294 numberDualInfeasibilitiesWithoutFree_(0), 1295 numberPrimalInfeasibilities_(100), 1296 numberRefinements_(0), 1297 pivotVariable_(NULL), 1298 factorization_(NULL), 1299 rowScale_(NULL), 1300 savedSolution_(NULL), 1301 columnScale_(NULL), 1302 scalingFlag_(3), 1303 numberTimesOptimal_(0), 1304 changeMade_(1), 1305 algorithm_(0), 1306 forceFactorization_(1), 1307 perturbation_(100), 1308 nonLinearCost_(NULL), 1309 specialOptions_(0), 1310 lastBadIteration_(999999), 1311 numberFake_(0), 1312 progressFlag_(0), 1313 firstFree_(1), 1314 incomingInfeasibility_(1.0), 1315 allowedInfeasibility_(10.0), 1316 progress_(NULL) 1317 { 1318 int i; 1319 for (i=0;i<6;i++) { 1320 rowArray_[i]=NULL; 1321 columnArray_[i]=NULL; 1322 } 1323 saveStatus_=NULL; 1324 factorization_ = NULL; 1325 dualRowPivot_ = NULL; 1326 primalColumnPivot_ = NULL; 1327 gutsOfDelete(0); 1328 specialOptions_ =0; 1329 delete nonLinearCost_; 1330 nonLinearCost_ = NULL; 94 cost_(NULL) 95 { 96 gutsOfDelete(); 1331 97 gutsOfCopy(rhs); 1332 solveType_= 1; // say simplexbased life form98 solveType_=2; // say interior based life form 1333 99 } 1334 100 // Copy constructor from model 1335 101 ClpInterior::ClpInterior(const ClpModel &rhs) : 1336 102 ClpModel(rhs), 1337 columnPrimalInfeasibility_(0.0),1338 rowPrimalInfeasibility_(0.0),1339 columnPrimalSequence_(2),1340 rowPrimalSequence_(2),1341 columnDualInfeasibility_(0.0),1342 rowDualInfeasibility_(0.0),1343 columnDualSequence_(2),1344 rowDualSequence_(2),1345 primalToleranceToGetOptimal_(1.0),1346 remainingDualInfeasibility_(0.0),1347 largeValue_(1.0e15),1348 103 largestPrimalError_(0.0), 1349 104 largestDualError_(0.0), 1350 largestSolutionError_(0.0),1351 dualBound_(1.0e10),1352 alpha_(0.0),1353 theta_(0.0),1354 lowerIn_(0.0),1355 valueIn_(0.0),1356 upperIn_(0.0),1357 dualIn_(0.0),1358 lowerOut_(1),1359 valueOut_(1),1360 upperOut_(1),1361 dualOut_(1),1362 dualTolerance_(0.0),1363 primalTolerance_(0.0),1364 105 sumDualInfeasibilities_(0.0), 1365 106 sumPrimalInfeasibilities_(0.0), 1366 infeasibilityCost_(1.0e10),1367 sumOfRelaxedDualInfeasibilities_(0.0),1368 sumOfRelaxedPrimalInfeasibilities_(0.0),1369 107 lower_(NULL), 1370 108 rowLowerWork_(NULL), … … 1373 111 rowUpperWork_(NULL), 1374 112 columnUpperWork_(NULL), 1375 cost_(NULL), 1376 rowObjectiveWork_(NULL), 1377 objectiveWork_(NULL), 1378 sequenceIn_(1), 1379 directionIn_(1), 1380 sequenceOut_(1), 1381 directionOut_(1), 1382 pivotRow_(1), 1383 lastGoodIteration_(100), 1384 dj_(NULL), 1385 rowReducedCost_(NULL), 1386 reducedCostWork_(NULL), 1387 solution_(NULL), 1388 rowActivityWork_(NULL), 1389 columnActivityWork_(NULL), 1390 numberDualInfeasibilities_(0), 1391 numberDualInfeasibilitiesWithoutFree_(0), 1392 numberPrimalInfeasibilities_(100), 1393 numberRefinements_(0), 1394 pivotVariable_(NULL), 1395 factorization_(NULL), 1396 rowScale_(NULL), 1397 savedSolution_(NULL), 1398 columnScale_(NULL), 1399 scalingFlag_(3), 1400 numberTimesOptimal_(0), 1401 changeMade_(1), 1402 algorithm_(0), 1403 forceFactorization_(1), 1404 perturbation_(100), 1405 nonLinearCost_(NULL), 1406 specialOptions_(0), 1407 lastBadIteration_(999999), 1408 numberFake_(0), 1409 progressFlag_(0), 1410 firstFree_(1), 1411 incomingInfeasibility_(1.0), 1412 allowedInfeasibility_(10.0), 1413 progress_(NULL) 1414 { 1415 int i; 1416 for (i=0;i<6;i++) { 1417 rowArray_[i]=NULL; 1418 columnArray_[i]=NULL; 1419 } 1420 saveStatus_=NULL; 1421 // get an empty factorization so we can set tolerances etc 1422 factorization_ = new ClpFactorization(); 1423 // say Steepest pricing 1424 dualRowPivot_ = new ClpDualRowSteepest(); 1425 // say Steepest pricing 1426 primalColumnPivot_ = new ClpPrimalColumnSteepest(); 1427 solveType_=1; // say simplex based life form 1428 113 cost_(NULL) 114 { 115 solveType_=2; // say interior based life form 1429 116 } 1430 117 // Assignment operator. This copies the data … … 1433 120 { 1434 121 if (this != &rhs) { 1435 gutsOfDelete(0); 1436 specialOptions_=0; 1437 delete nonLinearCost_; 1438 nonLinearCost_ = NULL; 122 gutsOfDelete(); 1439 123 ClpModel::operator=(rhs); 1440 124 gutsOfCopy(rhs); … … 1452 136 columnUpperWork_ = upper_; 1453 137 //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_)); 1454 cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_)); 1455 objectiveWork_ = cost_; 1456 rowObjectiveWork_ = cost_+numberColumns_; 1457 dj_ = ClpCopyOfArray(rhs.dj_,numberRows_+numberColumns_); 1458 if (dj_) { 1459 reducedCostWork_ = dj_; 1460 rowReducedCost_ = dj_+numberColumns_; 1461 } 1462 solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_); 1463 if (solution_) { 1464 columnActivityWork_ = solution_; 1465 rowActivityWork_ = solution_+numberColumns_; 1466 } 1467 if (rhs.pivotVariable_) { 1468 pivotVariable_ = new int[numberRows_]; 1469 ClpDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_); 1470 } else { 1471 pivotVariable_=NULL; 1472 } 1473 if (rhs.factorization_) { 1474 factorization_ = new ClpFactorization(*rhs.factorization_); 1475 } else { 1476 factorization_=NULL; 1477 } 1478 rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_); 1479 savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_); 1480 columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_); 1481 int i; 1482 for (i=0;i<6;i++) { 1483 rowArray_[i]=NULL; 1484 if (rhs.rowArray_[i]) 1485 rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]); 1486 columnArray_[i]=NULL; 1487 if (rhs.columnArray_[i]) 1488 columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]); 1489 } 1490 if (rhs.saveStatus_) { 1491 saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_); 1492 } 1493 columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_; 1494 columnPrimalSequence_ = rhs.columnPrimalSequence_; 1495 rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_; 1496 rowPrimalSequence_ = rhs.rowPrimalSequence_; 1497 columnDualInfeasibility_ = rhs.columnDualInfeasibility_; 1498 columnDualSequence_ = rhs.columnDualSequence_; 1499 rowDualInfeasibility_ = rhs.rowDualInfeasibility_; 1500 rowDualSequence_ = rhs.rowDualSequence_; 1501 primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_; 1502 remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_; 1503 largeValue_ = rhs.largeValue_; 138 cost_ = ClpCopyOfArray(rhs.cost_,numberColumns_); 1504 139 largestPrimalError_ = rhs.largestPrimalError_; 1505 140 largestDualError_ = rhs.largestDualError_; 1506 largestSolutionError_ = rhs.largestSolutionError_;1507 dualBound_ = rhs.dualBound_;1508 alpha_ = rhs.alpha_;1509 theta_ = rhs.theta_;1510 lowerIn_ = rhs.lowerIn_;1511 valueIn_ = rhs.valueIn_;1512 upperIn_ = rhs.upperIn_;1513 dualIn_ = rhs.dualIn_;1514 sequenceIn_ = rhs.sequenceIn_;1515 directionIn_ = rhs.directionIn_;1516 lowerOut_ = rhs.lowerOut_;1517 valueOut_ = rhs.valueOut_;1518 upperOut_ = rhs.upperOut_;1519 dualOut_ = rhs.dualOut_;1520 sequenceOut_ = rhs.sequenceOut_;1521 directionOut_ = rhs.directionOut_;1522 pivotRow_ = rhs.pivotRow_;1523 lastGoodIteration_ = rhs.lastGoodIteration_;1524 numberRefinements_ = rhs.numberRefinements_;1525 dualTolerance_ = rhs.dualTolerance_;1526 primalTolerance_ = rhs.primalTolerance_;1527 141 sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; 1528 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;1529 numberDualInfeasibilitiesWithoutFree_ =1530 rhs.numberDualInfeasibilitiesWithoutFree_;1531 142 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; 1532 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;1533 dualRowPivot_ = rhs.dualRowPivot_>clone(true);1534 primalColumnPivot_ = rhs.primalColumnPivot_>clone(true);1535 scalingFlag_ = rhs.scalingFlag_;1536 numberTimesOptimal_ = rhs.numberTimesOptimal_;1537 changeMade_ = rhs.changeMade_;1538 algorithm_ = rhs.algorithm_;1539 forceFactorization_ = rhs.forceFactorization_;1540 perturbation_ = rhs.perturbation_;1541 infeasibilityCost_ = rhs.infeasibilityCost_;1542 specialOptions_ = rhs.specialOptions_;1543 lastBadIteration_ = rhs.lastBadIteration_;1544 numberFake_ = rhs.numberFake_;1545 progressFlag_ = rhs.progressFlag_;1546 firstFree_ = rhs.firstFree_;1547 incomingInfeasibility_ = rhs.incomingInfeasibility_;1548 allowedInfeasibility_ = rhs.allowedInfeasibility_;1549 if (rhs.progress_)1550 progress_ = new ClpInteriorProgress(*rhs.progress_);1551 else1552 progress_=NULL;1553 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;1554 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;1555 if (rhs.nonLinearCost_!=NULL)1556 nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);1557 else1558 nonLinearCost_=NULL;1559 143 solveType_=rhs.solveType_; 1560 144 } 1561 145 // type == 0 do everything, most + pivot data, 2 factorization data as well 1562 146 void 1563 ClpInterior::gutsOfDelete( int type)147 ClpInterior::gutsOfDelete() 1564 148 { 1565 149 delete [] lower_; … … 1573 157 delete [] cost_; 1574 158 cost_=NULL; 1575 objectiveWork_=NULL;1576 rowObjectiveWork_=NULL;1577 delete [] dj_;1578 dj_=NULL;1579 reducedCostWork_=NULL;1580 rowReducedCost_=NULL;1581 delete [] solution_;1582 solution_=NULL;1583 rowActivityWork_=NULL;1584 columnActivityWork_=NULL;1585 delete [] rowScale_;1586 rowScale_ = NULL;1587 delete [] savedSolution_;1588 savedSolution_ = NULL;1589 delete [] columnScale_;1590 columnScale_ = NULL;1591 if ((specialOptions_&2)==0) {1592 delete nonLinearCost_;1593 nonLinearCost_ = NULL;1594 }1595 int i;1596 for (i=0;i<6;i++) {1597 delete rowArray_[i];1598 rowArray_[i]=NULL;1599 delete columnArray_[i];1600 columnArray_[i]=NULL;1601 }1602 delete rowCopy_;1603 rowCopy_=NULL;1604 delete [] saveStatus_;1605 saveStatus_=NULL;1606 if (!type) {1607 // delete everything1608 delete factorization_;1609 factorization_ = NULL;1610 delete [] pivotVariable_;1611 pivotVariable_=NULL;1612 delete dualRowPivot_;1613 dualRowPivot_ = NULL;1614 delete primalColumnPivot_;1615 primalColumnPivot_ = NULL;1616 delete progress_;1617 progress_=NULL;1618 } else {1619 // delete any size information in methods1620 if (type>1) {1621 factorization_>clearArrays();1622 delete [] pivotVariable_;1623 pivotVariable_=NULL;1624 }1625 dualRowPivot_>clearArrays();1626 primalColumnPivot_>clearArrays();1627 }1628 }1629 // This sets largest infeasibility and most infeasible1630 void1631 ClpInterior::checkPrimalSolution(const double * rowActivities,1632 const double * columnActivities)1633 {1634 double * solution;1635 int iRow,iColumn;1636 1637 objectiveValue_ = 0.0;1638 // now look at primal solution1639 columnPrimalInfeasibility_=0.0;1640 columnPrimalSequence_=1;1641 rowPrimalInfeasibility_=0.0;1642 rowPrimalSequence_=1;1643 largestSolutionError_=0.0;1644 solution = rowActivityWork_;1645 sumPrimalInfeasibilities_=0.0;1646 numberPrimalInfeasibilities_=0;1647 double primalTolerance = primalTolerance_;1648 double relaxedTolerance=primalTolerance_;1649 // we can't really trust infeasibilities if there is primal error1650 double error = min(1.0e3,largestPrimalError_);1651 // allow tolerance at least slightly bigger than standard1652 relaxedTolerance = relaxedTolerance + error;1653 sumOfRelaxedPrimalInfeasibilities_ = 0.0;1654 1655 for (iRow=0;iRow<numberRows_;iRow++) {1656 //assert (fabs(solution[iRow])<1.0e15getRowStatus(iRow) == basic);1657 double infeasibility=0.0;1658 objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];1659 if (solution[iRow]>rowUpperWork_[iRow]) {1660 infeasibility=solution[iRow]rowUpperWork_[iRow];1661 } else if (solution[iRow]<rowLowerWork_[iRow]) {1662 infeasibility=rowLowerWork_[iRow]solution[iRow];1663 }1664 if (infeasibility>primalTolerance) {1665 sumPrimalInfeasibilities_ += infeasibilityprimalTolerance_;1666 if (infeasibility>relaxedTolerance)1667 sumOfRelaxedPrimalInfeasibilities_ += infeasibilityrelaxedTolerance;1668 numberPrimalInfeasibilities_ ++;1669 }1670 if (infeasibility>rowPrimalInfeasibility_) {1671 rowPrimalInfeasibility_=infeasibility;1672 rowPrimalSequence_=iRow;1673 }1674 infeasibility = fabs(rowActivities[iRow]solution[iRow]);1675 if (infeasibility>largestSolutionError_)1676 largestSolutionError_=infeasibility;1677 }1678 solution = columnActivityWork_;1679 for (iColumn=0;iColumn<numberColumns_;iColumn++) {1680 //assert (fabs(solution[iColumn])<1.0e15getColumnStatus(iColumn) == basic);1681 double infeasibility=0.0;1682 objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];1683 if (solution[iColumn]>columnUpperWork_[iColumn]) {1684 infeasibility=solution[iColumn]columnUpperWork_[iColumn];1685 } else if (solution[iColumn]<columnLowerWork_[iColumn]) {1686 infeasibility=columnLowerWork_[iColumn]solution[iColumn];1687 }1688 if (infeasibility>columnPrimalInfeasibility_) {1689 columnPrimalInfeasibility_=infeasibility;1690 columnPrimalSequence_=iColumn;1691 }1692 if (infeasibility>primalTolerance) {1693 sumPrimalInfeasibilities_ += infeasibilityprimalTolerance_;1694 if (infeasibility>relaxedTolerance)1695 sumOfRelaxedPrimalInfeasibilities_ += infeasibilityrelaxedTolerance;1696 numberPrimalInfeasibilities_ ++;1697 }1698 infeasibility = fabs(columnActivities[iColumn]solution[iColumn]);1699 if (infeasibility>largestSolutionError_)1700 largestSolutionError_=infeasibility;1701 }1702 }1703 void1704 ClpInterior::checkDualSolution()1705 {1706 1707 int iRow,iColumn;1708 sumDualInfeasibilities_=0.0;1709 numberDualInfeasibilities_=0;1710 numberDualInfeasibilitiesWithoutFree_=0;1711 columnDualInfeasibility_=0.0;1712 columnDualSequence_=1;1713 rowDualInfeasibility_=0.0;1714 rowDualSequence_=1;1715 int firstFreePrimal = 1;1716 int firstFreeDual = 1;1717 int numberSuperBasicWithDj=0;1718 primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,1719 columnPrimalInfeasibility_);1720 remainingDualInfeasibility_=0.0;1721 double relaxedTolerance=dualTolerance_;1722 // we can't really trust infeasibilities if there is dual error1723 double error = min(1.0e3,largestDualError_);1724 // allow tolerance at least slightly bigger than standard1725 relaxedTolerance = relaxedTolerance + error;1726 sumOfRelaxedDualInfeasibilities_ = 0.0;1727 1728 for (iColumn=0;iColumn<numberColumns_;iColumn++) {1729 if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {1730 // not basic1731 double distanceUp = columnUpperWork_[iColumn]1732 columnActivityWork_[iColumn];1733 double distanceDown = columnActivityWork_[iColumn] 1734 columnLowerWork_[iColumn];1735 if (distanceUp>primalTolerance_) {1736 double value = reducedCostWork_[iColumn];1737 // Check if "free"1738 if (distanceDown>primalTolerance_) {1739 if (fabs(value)>1.0e2*relaxedTolerance) {1740 numberSuperBasicWithDj++;1741 if (firstFreeDual<0)1742 firstFreeDual = iColumn;1743 }1744 if (firstFreePrimal<0)1745 firstFreePrimal = iColumn;1746 }1747 // should not be negative1748 if (value<0.0) {1749 value =  value;1750 if (value>columnDualInfeasibility_) {1751 columnDualInfeasibility_=value;1752 columnDualSequence_=iColumn;1753 }1754 if (value>dualTolerance_) {1755 if (getColumnStatus(iColumn) != isFree) {1756 numberDualInfeasibilitiesWithoutFree_ ++;1757 sumDualInfeasibilities_ += valuedualTolerance_;1758 if (value>relaxedTolerance)1759 sumOfRelaxedDualInfeasibilities_ += valuerelaxedTolerance;1760 numberDualInfeasibilities_ ++;1761 } else {1762 // free so relax a lot1763 value *= 0.01;1764 if (value>dualTolerance_) {1765 sumDualInfeasibilities_ += valuedualTolerance_;1766 if (value>relaxedTolerance)1767 sumOfRelaxedDualInfeasibilities_ += valuerelaxedTolerance;1768 numberDualInfeasibilities_ ++;1769 }1770 }1771 // maybe we can make feasible by increasing tolerance1772 if (distanceUp<largeValue_) {1773 if (distanceUp>primalToleranceToGetOptimal_)1774 primalToleranceToGetOptimal_=distanceUp;1775 } else {1776 //gap too big for any tolerance1777 remainingDualInfeasibility_=1778 max(remainingDualInfeasibility_,value);1779 }1780 }1781 }1782 }1783 if (distanceDown>primalTolerance_) {1784 double value = reducedCostWork_[iColumn];1785 // should not be positive1786 if (value>0.0) {1787 if (value>columnDualInfeasibility_) {1788 columnDualInfeasibility_=value;1789 columnDualSequence_=iColumn;1790 }1791 if (value>dualTolerance_) {1792 sumDualInfeasibilities_ += valuedualTolerance_;1793 if (value>relaxedTolerance)1794 sumOfRelaxedDualInfeasibilities_ += valuerelaxedTolerance;1795 numberDualInfeasibilities_ ++;1796 if (getColumnStatus(iColumn) != isFree)1797 numberDualInfeasibilitiesWithoutFree_ ++;1798 // maybe we can make feasible by increasing tolerance1799 if (distanceDown<largeValue_&&1800 distanceDown>primalToleranceToGetOptimal_)1801 primalToleranceToGetOptimal_=distanceDown;1802 }1803 }1804 }1805 }1806 }1807 for (iRow=0;iRow<numberRows_;iRow++) {1808 if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {1809 // not basic1810 double distanceUp = rowUpperWork_[iRow]rowActivityWork_[iRow];1811 double distanceDown = rowActivityWork_[iRow] rowLowerWork_[iRow];1812 if (distanceUp>primalTolerance_) {1813 double value = rowReducedCost_[iRow];1814 // Check if "free"1815 if (distanceDown>primalTolerance_) {1816 if (fabs(value)>1.0e2*relaxedTolerance) {1817 numberSuperBasicWithDj++;1818 if (firstFreeDual<0)1819 firstFreeDual = iRow+numberColumns_;1820 }1821 if (firstFreePrimal<0)1822 firstFreePrimal = iRow+numberColumns_;1823 }1824 // should not be negative1825 if (value<0.0) {1826 value =  value;1827 if (value>rowDualInfeasibility_) {1828 rowDualInfeasibility_=value;1829 rowDualSequence_=iRow;1830 }1831 if (value>dualTolerance_) {1832 sumDualInfeasibilities_ += valuedualTolerance_;1833 if (value>relaxedTolerance)1834 sumOfRelaxedDualInfeasibilities_ += valuerelaxedTolerance;1835 numberDualInfeasibilities_ ++;1836 if (getRowStatus(iRow) != isFree)1837 numberDualInfeasibilitiesWithoutFree_ ++;1838 // maybe we can make feasible by increasing tolerance1839 if (distanceUp<largeValue_) {1840 if (distanceUp>primalToleranceToGetOptimal_)1841 primalToleranceToGetOptimal_=distanceUp;1842 } else {1843 //gap too big for any tolerance1844 remainingDualInfeasibility_=1845 max(remainingDualInfeasibility_,value);1846 }1847 }1848 }1849 }1850 if (distanceDown>primalTolerance_) {1851 double value = rowReducedCost_[iRow];1852 // should not be positive1853 if (value>0.0) {1854 if (value>rowDualInfeasibility_) {1855 rowDualInfeasibility_=value;1856 rowDualSequence_=iRow;1857 }1858 if (value>dualTolerance_) {1859 sumDualInfeasibilities_ += valuedualTolerance_;1860 if (value>relaxedTolerance)1861 sumOfRelaxedDualInfeasibilities_ += valuerelaxedTolerance;1862 numberDualInfeasibilities_ ++;1863 if (getRowStatus(iRow) != isFree)1864 numberDualInfeasibilitiesWithoutFree_ ++;1865 // maybe we can make feasible by increasing tolerance1866 if (distanceDown<largeValue_&&1867 distanceDown>primalToleranceToGetOptimal_)1868 primalToleranceToGetOptimal_=distanceDown;1869 }1870 }1871 }1872 }1873 }1874 if (algorithm_<0&&firstFreeDual>=0) {1875 // dual1876 firstFree_ = firstFreeDual;1877 } else if (numberSuperBasicWithDj1878 (progress_&&progress_>lastIterationNumber(0)<=0)) {1879 firstFree_=firstFreePrimal;1880 }1881 }1882 /*1883 Unpacks one column of the matrix into indexed array1884 */1885 void1886 ClpInterior::unpack(CoinIndexedVector * rowArray) const1887 {1888 rowArray>clear();1889 if (sequenceIn_>=numberColumns_) {1890 //slack1891 rowArray>insert(sequenceIn_numberColumns_,1.0);1892 } else {1893 // column1894 matrix_>unpack(this,rowArray,sequenceIn_);1895 }1896 }1897 void1898 ClpInterior::unpack(CoinIndexedVector * rowArray,int sequence) const1899 {1900 rowArray>clear();1901 if (sequence>=numberColumns_) {1902 //slack1903 rowArray>insert(sequencenumberColumns_,1.0);1904 } else {1905 // column1906 matrix_>unpack(this,rowArray,sequence);1907 }1908 }1909 /*1910 Unpacks one column of the matrix into indexed array1911 */1912 void1913 ClpInterior::unpackPacked(CoinIndexedVector * rowArray)1914 {1915 rowArray>clear();1916 if (sequenceIn_>=numberColumns_) {1917 //slack1918 int * index = rowArray>getIndices();1919 double * array = rowArray>denseVector();1920 array[0]=1.0;1921 index[0]=sequenceIn_numberColumns_;1922 rowArray>setNumElements(1);1923 rowArray>setPackedMode(true);1924 } else {1925 // column1926 matrix_>unpackPacked(this,rowArray,sequenceIn_);1927 }1928 }1929 void1930 ClpInterior::unpackPacked(CoinIndexedVector * rowArray,int sequence)1931 {1932 rowArray>clear();1933 if (sequence>=numberColumns_) {1934 //slack1935 int * index = rowArray>getIndices();1936 double * array = rowArray>denseVector();1937 array[0]=1.0;1938 index[0]=sequencenumberColumns_;1939 rowArray>setNumElements(1);1940 rowArray>setPackedMode(true);1941 } else {1942 // column1943 matrix_>unpackPacked(this,rowArray,sequence);1944 }1945 159 } 1946 160 bool 1947 ClpInterior::create Rim(int what,bool makeRowCopy)161 ClpInterior::createWorkingData() 1948 162 { 1949 163 bool goodMatrix=true; 1950 int saveLevel=handler_>logLevel(); 1951 if (problemStatus_==10) { 1952 handler_>setLogLevel(0); // switch off messages 1953 } else if (factorization_) { 1954 // match up factorization messages 1955 if (handler_>logLevel()<3) 1956 factorization_>messageLevel(0); 1957 else 1958 factorization_>messageLevel(3); 1959 } 1960 int i; 1961 if ((what&(16+32))!=0) { 1962 // move information to work arrays 1963 double direction = optimizationDirection_; 1964 // direction is actually scale out not scale in 1965 if (direction) 1966 direction = 1.0/direction; 1967 if (direction!=1.0) { 1968 // reverse all dual signs 1969 for (i=0;i<numberColumns_;i++) 1970 reducedCost_[i] *= direction; 1971 for (i=0;i<numberRows_;i++) 1972 dual_[i] *= direction; 1973 } 1974 // row reduced costs 1975 if (!dj_) { 1976 dj_ = new double[numberRows_+numberColumns_]; 1977 reducedCostWork_ = dj_; 1978 rowReducedCost_ = dj_+numberColumns_; 1979 memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double)); 1980 memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double)); 1981 } 1982 if (!solution_(what&32)!=0) { 1983 if (!solution_) 1984 solution_ = new double[numberRows_+numberColumns_]; 1985 columnActivityWork_ = solution_; 1986 rowActivityWork_ = solution_+numberColumns_; 1987 memcpy(columnActivityWork_,columnActivity_, 1988 numberColumns_*sizeof(double)); 1989 memcpy(rowActivityWork_,rowActivity_, 1990 numberRows_*sizeof(double)); 1991 } 1992 } 1993 if ((what&16)!=0) { 1994 //check matrix 1995 if (!matrix_>allElementsInRange(this,smallElement_,1.0e20)) { 1996 problemStatus_=4; 1997 goodMatrix= false; 1998 } 1999 if (makeRowCopy) { 2000 delete rowCopy_; 2001 // may return NULL if can't give row copy 2002 rowCopy_ = matrix_>reverseOrderedCopy(); 2003 #ifdef TAKEOUT 2004 { 2005 2006 ClpPackedMatrix* rowCopy = 2007 dynamic_cast< ClpPackedMatrix*>(rowCopy_); 2008 const int * column = rowCopy>getIndices(); 2009 const CoinBigIndex * rowStart = rowCopy>getVectorStarts(); 2010 const double * element = rowCopy>getElements(); 2011 int i; 2012 for (i=133;i<numberRows_;i++) { 2013 if (rowStart[i+1]rowStart[i]==10rowStart[i+1]rowStart[i]==15) 2014 printf("Row %d has %d elements\n",i,rowStart[i+1]rowStart[i]); 2015 } 2016 } 2017 #endif 2018 } 2019 } 2020 if ((what&4)!=0) { 2021 delete [] cost_; 2022 // extra copy with original costs 2023 int nTotal = numberRows_+numberColumns_; 2024 //cost_ = new double[2*nTotal]; 2025 cost_ = new double[nTotal]; 2026 objectiveWork_ = cost_; 2027 rowObjectiveWork_ = cost_+numberColumns_; 2028 memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double)); 2029 if (rowObjective_) 2030 memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double)); 2031 else 2032 memset(rowObjectiveWork_,0,numberRows_*sizeof(double)); 2033 // and initialize changes to zero 2034 //memset(cost_+nTotal,0,nTotal*sizeof(double)); 2035 } 2036 if ((what&1)!=0) { 2037 delete [] lower_; 2038 delete [] upper_; 2039 lower_ = new double[numberColumns_+numberRows_]; 2040 upper_ = new double[numberColumns_+numberRows_]; 2041 rowLowerWork_ = lower_+numberColumns_; 2042 columnLowerWork_ = lower_; 2043 rowUpperWork_ = upper_+numberColumns_; 2044 columnUpperWork_ = upper_; 2045 memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double)); 2046 memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double)); 2047 memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double)); 2048 memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double)); 2049 // clean up any mismatches on infinity 2050 for (i=0;i<numberColumns_;i++) { 2051 if (columnLowerWork_[i]<1.0e30) 2052 columnLowerWork_[i] = COIN_DBL_MAX; 2053 if (columnUpperWork_[i]>1.0e30) 2054 columnUpperWork_[i] = COIN_DBL_MAX; 2055 } 2056 // clean up any mismatches on infinity 2057 for (i=0;i<numberRows_;i++) { 2058 if (rowLowerWork_[i]<1.0e30) 2059 rowLowerWork_[i] = COIN_DBL_MAX; 2060 if (rowUpperWork_[i]>1.0e30) 2061 rowUpperWork_[i] = COIN_DBL_MAX; 2062 } 2063 } 2064 // do scaling if needed 2065 if (scalingFlag_>0&&!rowScale_) { 2066 if (matrix_>scale(this)) 2067 scalingFlag_=scalingFlag_; // not scaled after all 2068 } 2069 if ((what&4)!=0) { 2070 double direction = optimizationDirection_; 2071 // direction is actually scale out not scale in 2072 if (direction) 2073 direction = 1.0/direction; 2074 // but also scale by scale factors 2075 // not really sure about row scaling 2076 if (!rowScale_) { 2077 if (direction!=1.0) { 2078 for (i=0;i<numberRows_;i++) 2079 rowObjectiveWork_[i] *= direction; 2080 for (i=0;i<numberColumns_;i++) 2081 objectiveWork_[i] *= direction; 2082 } 2083 } else { 2084 for (i=0;i<numberRows_;i++) 2085 rowObjectiveWork_[i] *= direction/rowScale_[i]; 2086 for (i=0;i<numberColumns_;i++) 2087 objectiveWork_[i] *= direction*columnScale_[i]; 2088 } 2089 } 2090 if ((what&(1+32))!=0&&rowScale_) { 2091 for (i=0;i<numberColumns_;i++) { 2092 double multiplier = 1.0/columnScale_[i]; 2093 if (columnLowerWork_[i]>1.0e50) 2094 columnLowerWork_[i] *= multiplier; 2095 if (columnUpperWork_[i]<1.0e50) 2096 columnUpperWork_[i] *= multiplier; 2097 2098 } 2099 for (i=0;i<numberRows_;i++) { 2100 if (rowLowerWork_[i]>1.0e50) 2101 rowLowerWork_[i] *= rowScale_[i]; 2102 if (rowUpperWork_[i]<1.0e50) 2103 rowUpperWork_[i] *= rowScale_[i]; 2104 } 2105 } 2106 if ((what&(8+32))!=0&&rowScale_) { 2107 // on entry 2108 for (i=0;i<numberColumns_;i++) { 2109 columnActivityWork_[i] /= columnScale_[i]; 2110 reducedCostWork_[i] *= columnScale_[i]; 2111 } 2112 for (i=0;i<numberRows_;i++) { 2113 rowActivityWork_[i] *= rowScale_[i]; 2114 dual_[i] /= rowScale_[i]; 2115 rowReducedCost_[i] = dual_[i]; 2116 } 2117 } 2118 2119 if ((what&16)!=0) { 2120 // check rim of problem okay 2121 if (!sanityCheck()) 2122 goodMatrix=false; 2123 } 2124 // we need to treat matrix as if each element by rowScaleIn and columnScaleout?? 2125 // maybe we need to move scales to SimplexModel for factorization? 2126 if ((what&8)!=0&&!pivotVariable_) { 2127 pivotVariable_=new int[numberRows_]; 2128 } 2129 2130 if ((what&16)!=0&&!rowArray_[2]) { 2131 // get some arrays 2132 int iRow,iColumn; 2133 // these are "indexed" arrays so we always know where nonzeros are 2134 /********************************************************** 2135 rowArray_[3] is long enough for rows+columns 2136 *********************************************************/ 2137 for (iRow=0;iRow<4;iRow++) { 2138 if (!rowArray_[iRow]) { 2139 rowArray_[iRow]=new CoinIndexedVector(); 2140 int length =numberRows_+factorization_>maximumPivots(); 2141 if (iRow==3) 2142 length += numberColumns_; 2143 rowArray_[iRow]>reserve(length); 2144 } 2145 } 2146 2147 for (iColumn=0;iColumn<2;iColumn++) { 2148 if (!columnArray_[iColumn]) { 2149 columnArray_[iColumn]=new CoinIndexedVector(); 2150 columnArray_[iColumn]>reserve(numberColumns_); 2151 } 2152 } 2153 } 2154 double primalTolerance=dblParam_[ClpPrimalTolerance]; 2155 if ((what&1)!=0) { 2156 // fix any variables with tiny gaps 2157 for (i=0;i<numberColumns_;i++) { 2158 if (columnUpperWork_[i]columnLowerWork_[i]<=primalTolerance) { 2159 if (columnLowerWork_[i]>=0.0) { 2160 columnUpperWork_[i] = columnLowerWork_[i]; 2161 } else if (columnUpperWork_[i]<=0.0) { 2162 columnLowerWork_[i] = columnUpperWork_[i]; 2163 } else { 2164 columnUpperWork_[i] = 0.0; 2165 columnLowerWork_[i] = 0.0; 2166 } 2167 } 2168 } 2169 for (i=0;i<numberRows_;i++) { 2170 if (rowUpperWork_[i]rowLowerWork_[i]<=primalTolerance) { 2171 if (rowLowerWork_[i]>=0.0) { 2172 rowUpperWork_[i] = rowLowerWork_[i]; 2173 } else if (rowUpperWork_[i]<=0.0) { 2174 rowLowerWork_[i] = rowUpperWork_[i]; 2175 } else { 2176 rowUpperWork_[i] = 0.0; 2177 rowLowerWork_[i] = 0.0; 2178 } 2179 } 2180 } 2181 } 2182 if (problemStatus_==10) { 2183 problemStatus_=1; 2184 handler_>setLogLevel(saveLevel); // switch back messages 2185 } 164 //check matrix 165 if (!matrix_>allElementsInRange(this,1.0e12,1.0e20)) { 166 problemStatus_=4; 167 goodMatrix= false; 168 } 169 // check rim of problem okay 170 if (!sanityCheck()) 171 goodMatrix=false; 2186 172 return goodMatrix; 2187 173 } 2188 174 void 2189 ClpInterior::deleteRim(int getRidOfFactorizationData) 2190 { 2191 int i; 2192 if (problemStatus_!=1&&problemStatus_!=2) { 2193 delete [] ray_; 2194 ray_=NULL; 2195 } 2196 // ray may be null if in branch and bound 2197 if (rowScale_) { 2198 for (i=0;i<numberColumns_;i++) { 2199 columnActivity_[i] = columnActivityWork_[i]*columnScale_[i]; 2200 reducedCost_[i] = reducedCostWork_[i]/columnScale_[i]; 2201 } 2202 for (i=0;i<numberRows_;i++) { 2203 rowActivity_[i] = rowActivityWork_[i]/rowScale_[i]; 2204 dual_[i] *= rowScale_[i]; 2205 } 2206 if (problemStatus_==2) { 2207 for (i=0;i<numberColumns_;i++) { 2208 ray_[i] *= columnScale_[i]; 2209 } 2210 } else if (problemStatus_==1&&ray_) { 2211 for (i=0;i<numberRows_;i++) { 2212 ray_[i] *= rowScale_[i]; 2213 } 2214 } 2215 } else { 2216 for (i=0;i<numberColumns_;i++) { 2217 columnActivity_[i] = columnActivityWork_[i]; 2218 reducedCost_[i] = reducedCostWork_[i]; 2219 } 2220 for (i=0;i<numberRows_;i++) { 2221 rowActivity_[i] = rowActivityWork_[i]; 2222 } 2223 } 2224 if (optimizationDirection_!=1.0) { 2225 // and modify all dual signs 2226 for (i=0;i<numberColumns_;i++) 2227 reducedCost_[i] *= optimizationDirection_; 2228 for (i=0;i<numberRows_;i++) 2229 dual_[i] *= optimizationDirection_; 2230 } 2231 // scaling may have been turned off 2232 scalingFlag_ = abs(scalingFlag_); 2233 if(getRidOfFactorizationData>=0) 2234 gutsOfDelete(getRidOfFactorizationData+1); 2235 } 2236 void 2237 ClpInterior::setDualBound(double value) 2238 { 2239 if (value>0.0) 2240 dualBound_=value; 2241 } 2242 void 2243 ClpInterior::setInfeasibilityCost(double value) 2244 { 2245 if (value>0.0) 2246 infeasibilityCost_=value; 2247 } 2248 void ClpInterior::setNumberRefinements( int value) 2249 { 2250 if (value>=0&&value<10) 2251 numberRefinements_=value; 2252 } 2253 // Sets row pivot choice algorithm in dual 2254 void 2255 ClpInterior::setDualRowPivotAlgorithm(ClpDualRowPivot & choice) 2256 { 2257 delete dualRowPivot_; 2258 dualRowPivot_ = choice.clone(true); 2259 } 2260 // Sets row pivot choice algorithm in dual 2261 void 2262 ClpInterior::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice) 2263 { 2264 delete primalColumnPivot_; 2265 primalColumnPivot_ = choice.clone(true); 2266 } 2267 // Sets or unsets scaling, 0 off, 1 on, 2 dynamic(later) 2268 void 2269 ClpInterior::scaling(int mode) 2270 { 2271 if (mode>0&&mode<4) { 2272 scalingFlag_=mode; 2273 } else if (!mode) { 2274 scalingFlag_=0; 2275 delete [] rowScale_; 2276 rowScale_ = NULL; 2277 delete [] columnScale_; 2278 columnScale_ = NULL; 2279 } 2280 } 2281 // Passes in factorization 2282 void 2283 ClpInterior::setFactorization( ClpFactorization & factorization) 2284 { 2285 delete factorization_; 2286 factorization_= new ClpFactorization(factorization); 2287 } 2288 void 2289 ClpInterior::times(double scalar, 2290 const double * x, double * y) const 2291 { 2292 if (rowScale_) 2293 matrix_>times(scalar,x,y,rowScale_,columnScale_); 2294 else 2295 matrix_>times(scalar,x,y); 2296 } 2297 void 2298 ClpInterior::transposeTimes(double scalar, 2299 const double * x, double * y) const 2300 { 2301 if (rowScale_) 2302 matrix_>transposeTimes(scalar,x,y,rowScale_,columnScale_); 2303 else 2304 matrix_>transposeTimes(scalar,x,y); 2305 } 2306 /* Perturbation: 2307 50 to +50  perturb by this power of ten (6 sounds good) 2308 100  auto perturb if takes too long (1.0e6 largest nonzero) 2309 101  we are perturbed 2310 102  don't try perturbing again 2311 default is 100 2312 */ 2313 void 2314 ClpInterior::setPerturbation(int value) 2315 { 2316 if(value<=100&&value >=1000) { 2317 perturbation_=value; 2318 } 2319 } 2320 // Sparsity on or off 2321 bool 2322 ClpInterior::sparseFactorization() const 2323 { 2324 return factorization_>sparseThreshold()!=0; 2325 } 2326 void 2327 ClpInterior::setSparseFactorization(bool value) 2328 { 2329 if (value) { 2330 if (!factorization_>sparseThreshold()) 2331 factorization_>goSparse(); 2332 } else { 2333 factorization_>sparseThreshold(0); 2334 } 2335 } 2336 /* Tightens primal bounds to make dual faster. Unless 2337 fixed, bounds are slightly looser than they could be. 2338 This is to make dual go faster and is probably not needed 2339 with a presolve. Returns nonzero if problem infeasible 2340 2341 Fudge for branch and bound  put bounds on columns of factor * 2342 largest value (at continuous)  should improve stability 2343 in branch and bound on infeasible branches (0.0 is off) 2344 */ 2345 int 2346 ClpInterior::tightenPrimalBounds(double factor) 2347 { 2348 2349 // Get a row copy in standard format 2350 CoinPackedMatrix copy; 2351 copy.reverseOrderedCopyOf(*matrix()); 2352 // get matrix data pointers 2353 const int * column = copy.getIndices(); 2354 const CoinBigIndex * rowStart = copy.getVectorStarts(); 2355 const int * rowLength = copy.getVectorLengths(); 2356 const double * element = copy.getElements(); 2357 int numberChanged=1,iPass=0; 2358 double large = largeValue(); // treat bounds > this as infinite 2359 int numberInfeasible=0; 2360 int totalTightened = 0; 2361 2362 double tolerance = primalTolerance(); 2363 2364 2365 // Save column bounds 2366 double * saveLower = new double [numberColumns_]; 2367 memcpy(saveLower,columnLower_,numberColumns_*sizeof(double)); 2368 double * saveUpper = new double [numberColumns_]; 2369 memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double)); 2370 2371 int iRow, iColumn; 2372 2373 // If wanted  tighten column bounds using solution 2374 if (factor) { 2375 assert (factor>1.0); 2376 double largest=0.0; 2377 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 2378 if (columnUpper_[iColumn]columnLower_[iColumn]>tolerance) { 2379 largest = max(largest,fabs(columnActivity_[iColumn])); 2380 } 2381 } 2382 largest *= factor; 2383 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 2384 if (columnUpper_[iColumn]columnLower_[iColumn]>tolerance) { 2385 columnUpper_[iColumn] = min(columnUpper_[iColumn],largest); 2386 columnLower_[iColumn] = max(columnLower_[iColumn],largest); 2387 } 2388 } 2389 } 2390 #define MAXPASS 10 2391 2392 // Loop round seeing if we can tighten bounds 2393 // Would be faster to have a stack of possible rows 2394 // and we put altered rows back on stack 2395 int numberCheck=1; 2396 while(numberChanged>numberCheck) { 2397 2398 numberChanged = 0; // Bounds tightened this pass 2399 2400 if (iPass==MAXPASS) break; 2401 iPass++; 2402 2403 for (iRow = 0; iRow < numberRows_; iRow++) { 2404 2405 if (rowLower_[iRow]>largerowUpper_[iRow]<large) { 2406 2407 // possible row 2408 int infiniteUpper = 0; 2409 int infiniteLower = 0; 2410 double maximumUp = 0.0; 2411 double maximumDown = 0.0; 2412 double newBound; 2413 CoinBigIndex rStart = rowStart[iRow]; 2414 CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow]; 2415 CoinBigIndex j; 2416 // Compute possible lower and upper ranges 2417 2418 for (j = rStart; j < rEnd; ++j) { 2419 double value=element[j]; 2420 iColumn = column[j]; 2421 if (value > 0.0) { 2422 if (columnUpper_[iColumn] >= large) { 2423 ++infiniteUpper; 2424 } else { 2425 maximumUp += columnUpper_[iColumn] * value; 2426 } 2427 if (columnLower_[iColumn] <= large) { 2428 ++infiniteLower; 2429 } else { 2430 maximumDown += columnLower_[iColumn] * value; 2431 } 2432 } else if (value<0.0) { 2433 if (columnUpper_[iColumn] >= large) { 2434 ++infiniteLower; 2435 } else { 2436 maximumDown += columnUpper_[iColumn] * value; 2437 } 2438 if (columnLower_[iColumn] <= large) { 2439 ++infiniteUpper; 2440 } else { 2441 maximumUp += columnLower_[iColumn] * value; 2442 } 2443 } 2444 } 2445 double maxUp = maximumUp+infiniteUpper*1.0e31; 2446 double maxDown = maximumDowninfiniteLower*1.0e31; 2447 if (maxUp <= rowUpper_[iRow] + tolerance && 2448 maxDown >= rowLower_[iRow]  tolerance) { 2449 2450 // Row is redundant  make totally free 2451 // NO  can't do this for postsolve 2452 // rowLower_[iRow]=COIN_DBL_MAX; 2453 // rowUpper_[iRow]=COIN_DBL_MAX; 2454 //printf("Redundant row in presolveX %d\n",iRow); 2455 2456 } else { 2457 if (maxUp < rowLower_[iRow] tolerance  2458 maxDown > rowUpper_[iRow]+tolerance) { 2459 // problem is infeasible  exit at once 2460 numberInfeasible++; 2461 break; 2462 } 2463 double lower = rowLower_[iRow]; 2464 double upper = rowUpper_[iRow]; 2465 for (j = rStart; j < rEnd; ++j) { 2466 double value=element[j]; 2467 iColumn = column[j]; 2468 double nowLower = columnLower_[iColumn]; 2469 double nowUpper = columnUpper_[iColumn]; 2470 if (value > 0.0) { 2471 // positive value 2472 if (lower>large) { 2473 if (!infiniteUpper) { 2474 assert(nowUpper < large); 2475 newBound = nowUpper + 2476 (lower  maximumUp) / value; 2477 } else if (infiniteUpper==1&&nowUpper>large) { 2478 newBound = (lower maximumUp) / value; 2479 } else { 2480 newBound = COIN_DBL_MAX; 2481 } 2482 if (newBound > nowLower + 1.0e12) { 2483 // Tighten the lower bound 2484 columnLower_[iColumn] = newBound; 2485 numberChanged++; 2486 // check infeasible (relaxed) 2487 if (nowUpper  newBound < 2488 100.0*tolerance) { 2489 numberInfeasible++; 2490 } 2491 // adjust 2492 double now; 2493 if (nowLower<large) { 2494 now=0.0; 2495 infiniteLower; 2496 } else { 2497 now = nowLower; 2498 } 2499 maximumDown += (newBoundnow) * value; 2500 nowLower = newBound; 2501 } 2502 } 2503 if (upper <large) { 2504 if (!infiniteLower) { 2505 assert(nowLower > large); 2506 newBound = nowLower + 2507 (upper  maximumDown) / value; 2508 } else if (infiniteLower==1&&nowLower<large) { 2509 newBound = (upper  maximumDown) / value; 2510 } else { 2511 newBound = COIN_DBL_MAX; 2512 } 2513 if (newBound < nowUpper  1.0e12) { 2514 // Tighten the upper bound 2515 columnUpper_[iColumn] = newBound; 2516 numberChanged++; 2517 // check infeasible (relaxed) 2518 if (newBound  nowLower < 2519 100.0*tolerance) { 2520 numberInfeasible++; 2521 } 2522 // adjust 2523 double now; 2524 if (nowUpper>large) { 2525 now=0.0; 2526 infiniteUpper; 2527 } else { 2528 now = nowUpper; 2529 } 2530 maximumUp += (newBoundnow) * value; 2531 nowUpper = newBound; 2532 } 2533 } 2534 } else { 2535 // negative value 2536 if (lower>large) { 2537 if (!infiniteUpper) { 2538 assert(nowLower > large); 2539 newBound = nowLower + 2540 (lower  maximumUp) / value; 2541 } else if (infiniteUpper==1&&nowLower<large) { 2542 newBound = (lower maximumUp) / value; 2543 } else { 2544 newBound = COIN_DBL_MAX; 2545 } 2546 if (newBound < nowUpper  1.0e12) { 2547 // Tighten the upper bound 2548 columnUpper_[iColumn] = newBound; 2549 numberChanged++; 2550 // check infeasible (relaxed) 2551 if (newBound  nowLower < 2552 100.0*tolerance) { 2553 numberInfeasible++; 2554 } 2555 // adjust 2556 double now; 2557 if (nowUpper>large) { 2558 now=0.0; 2559 infiniteLower; 2560 } else { 2561 now = nowUpper; 2562 } 2563 maximumDown += (newBoundnow) * value; 2564 nowUpper = newBound; 2565 } 2566 } 2567 if (upper <large) { 2568 if (!infiniteLower) { 2569 assert(nowUpper < large); 2570 newBound = nowUpper + 2571 (upper  maximumDown) / value; 2572 } else if (infiniteLower==1&&nowUpper>large) { 2573 newBound = (upper  maximumDown) / value; 2574 } else { 2575 newBound = COIN_DBL_MAX; 2576 } 2577 if (newBound > nowLower + 1.0e12) { 2578 // Tighten the lower bound 2579 columnLower_[iColumn] = newBound; 2580 numberChanged++; 2581 // check infeasible (relaxed) 2582 if (nowUpper  newBound < 2583 100.0*tolerance) { 2584 numberInfeasible++; 2585 } 2586 // adjust 2587 double now; 2588 if (nowLower<large) { 2589 now=0.0; 2590 infiniteUpper; 2591 } else { 2592 now = nowLower; 2593 } 2594 maximumUp += (newBoundnow) * value; 2595 nowLower = newBound; 2596 } 2597 } 2598 } 2599 } 2600 } 2601 } 2602 } 2603 totalTightened += numberChanged; 2604 if (iPass==1) 2605 numberCheck=numberChanged>>4; 2606 if (numberInfeasible) break; 2607 } 2608 if (!numberInfeasible) { 2609 handler_>message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_) 2610 <<totalTightened 2611 <<CoinMessageEol; 2612 // Set bounds slightly loose 2613 double useTolerance = 1.0e3; 2614 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 2615 if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) { 2616 if (columnUpper_[iColumn]columnLower_[iColumn]<useTolerance+1.0e8) { 2617 // relax enough so will have correct dj 2618 #if 1 2619 columnLower_[iColumn]=max(saveLower[iColumn], 2620 columnLower_[iColumn]100.0*useTolerance); 2621 columnUpper_[iColumn]=min(saveUpper[iColumn], 2622 columnUpper_[iColumn]+100.0*useTolerance); 2623 #else 2624 if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) { 2625 if (columnUpper_[iColumn] 100.0*useTolerance>saveLower[iColumn]) { 2626 columnLower_[iColumn]=columnUpper_[iColumn]100.0*useTolerance; 2627 } else { 2628 columnLower_[iColumn]=saveLower[iColumn]; 2629 columnUpper_[iColumn]=min(saveUpper[iColumn], 2630 saveLower[iColumn]+100.0*useTolerance); 2631 } 2632 } else { 2633 if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) { 2634 columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance; 2635 } else { 2636 columnUpper_[iColumn]=saveUpper[iColumn]; 2637 columnLower_[iColumn]=max(saveLower[iColumn], 2638 saveUpper[iColumn]100.0*useTolerance); 2639 } 2640 } 2641 #endif 2642 } else { 2643 if (columnUpper_[iColumn]<saveUpper[iColumn]) { 2644 // relax a bit 2645 columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*useTolerance, 2646 saveUpper[iColumn]); 2647 } 2648 if (columnLower_[iColumn]>saveLower[iColumn]) { 2649 // relax a bit 2650 columnLower_[iColumn] = max(columnLower_[iColumn]100.0*useTolerance, 2651 saveLower[iColumn]); 2652 } 2653 } 2654 } 2655 } 2656 } else { 2657 handler_>message(CLP_SIMPLEX_INFEASIBILITIES,messages_) 2658 <<numberInfeasible 2659 <<CoinMessageEol; 2660 // restore column bounds 2661 memcpy(columnLower_,saveLower,numberColumns_*sizeof(double)); 2662 memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double)); 2663 } 2664 delete [] saveLower; 2665 delete [] saveUpper; 2666 return (numberInfeasible); 2667 } 2668 // dual 2669 #include "ClpInteriorDual.hpp" 2670 #include "ClpInteriorPrimal.hpp" 2671 int ClpInterior::dual (int ifValuesPass ) 2672 { 2673 int returnCode = ((ClpInteriorDual *) this)>dual(ifValuesPass); 2674 if (problemStatus_==10) { 2675 //printf("Cleaning up with primal\n"); 2676 int savePerturbation = perturbation_; 2677 perturbation_=100; 2678 bool denseFactorization = initialDenseFactorization(); 2679 // It will be safe to allow dense 2680 setInitialDenseFactorization(true); 2681 returnCode = ((ClpInteriorPrimal *) this)>primal(1); 2682 setInitialDenseFactorization(denseFactorization); 2683 perturbation_=savePerturbation; 2684 } 2685 return returnCode; 2686 } 2687 // primal 2688 int ClpInterior::primal (int ifValuesPass ) 2689 { 2690 int returnCode = ((ClpInteriorPrimal *) this)>primal(ifValuesPass); 2691 if (problemStatus_==10) { 2692 //printf("Cleaning up with dual\n"); 2693 int savePerturbation = perturbation_; 2694 perturbation_=100; 2695 bool denseFactorization = initialDenseFactorization(); 2696 // It will be safe to allow dense 2697 setInitialDenseFactorization(true); 2698 returnCode = ((ClpInteriorDual *) this)>dual(0); 2699 setInitialDenseFactorization(denseFactorization); 2700 perturbation_=savePerturbation; 2701 } 2702 return returnCode; 2703 } 2704 #include "ClpInteriorPrimalQuadratic.hpp" 2705 /* Solves quadratic problem using SLP  may be used as crash 2706 for other algorithms when number of iterations small 2707 */ 2708 int 2709 ClpInterior::quadraticSLP(int numberPasses, double deltaTolerance) 2710 { 2711 return ((ClpInteriorPrimalQuadratic *) this)>primalSLP(numberPasses,deltaTolerance); 2712 } 2713 // Solves quadratic using Dantzig's algorithm  primal 2714 int 2715 ClpInterior::quadraticPrimal(int phase) 2716 { 2717 return ((ClpInteriorPrimalQuadratic *) this)>primalQuadratic(phase); 2718 } 2719 /* For strong branching. On input lower and upper are new bounds 2720 while on output they are objective function values (>1.0e50 infeasible). 2721 Return code is 0 if nothing interesting, 1 if infeasible both 2722 ways and +1 if infeasible one way (check values to see which one(s)) 2723 */ 2724 int ClpInterior::strongBranching(int numberVariables,const int * variables, 2725 double * newLower, double * newUpper, 2726 double ** outputSolution, 2727 int * outputStatus, int * outputIterations, 2728 bool stopOnFirstInfeasible, 2729 bool alwaysFinish) 2730 { 2731 return ((ClpInteriorDual *) this)>strongBranching(numberVariables,variables, 2732 newLower, newUpper,outputSolution, 2733 outputStatus, outputIterations, 2734 stopOnFirstInfeasible, 2735 alwaysFinish); 2736 } 2737 /* Borrow model. This is so we dont have to copy large amounts 2738 of data around. It assumes a derived class wants to overwrite 2739 an empty model with a real one  while it does an algorithm. 2740 This is same as ClpModel one, but sets scaling on etc. */ 2741 void 2742 ClpInterior::borrowModel(ClpModel & otherModel) 2743 { 2744 ClpModel::borrowModel(otherModel); 2745 createStatus(); 2746 ClpDualRowSteepest steep1; 2747 setDualRowPivotAlgorithm(steep1); 2748 ClpPrimalColumnSteepest steep2; 2749 setPrimalColumnPivotAlgorithm(steep2); 2750 } 2751 typedef struct { 2752 double optimizationDirection; 2753 double dblParam[ClpLastDblParam]; 2754 double objectiveValue; 2755 double dualBound; 2756 double dualTolerance; 2757 double primalTolerance; 2758 double sumDualInfeasibilities; 2759 double sumPrimalInfeasibilities; 2760 double infeasibilityCost; 2761 int numberRows; 2762 int numberColumns; 2763 int intParam[ClpLastIntParam]; 2764 int numberIterations; 2765 int problemStatus; 2766 int maximumIterations; 2767 int lengthNames; 2768 int numberDualInfeasibilities; 2769 int numberDualInfeasibilitiesWithoutFree; 2770 int numberPrimalInfeasibilities; 2771 int numberRefinements; 2772 int scalingFlag; 2773 int algorithm; 2774 unsigned int specialOptions; 2775 int dualPivotChoice; 2776 int primalPivotChoice; 2777 int matrixStorageChoice; 2778 } Clp_scalars; 2779 2780 int outDoubleArray(double * array, int length, FILE * fp) 2781 { 2782 int numberWritten; 2783 if (array&&length) { 2784 numberWritten = fwrite(&length,sizeof(int),1,fp); 2785 if (numberWritten!=1) 2786 return 1; 2787 numberWritten = fwrite(array,sizeof(double),length,fp); 2788 if (numberWritten!=length) 2789 return 1; 2790 } else { 2791 length = 0; 2792 numberWritten = fwrite(&length,sizeof(int),1,fp); 2793 if (numberWritten!=1) 2794 return 1; 2795 } 2796 return 0; 2797 } 2798 // Save model to file, returns 0 if success 2799 int 2800 ClpInterior::saveModel(const char * fileName) 2801 { 2802 FILE * fp = fopen(fileName,"wb"); 2803 if (fp) { 2804 Clp_scalars scalars; 2805 int i; 2806 CoinBigIndex numberWritten; 2807 // Fill in scalars 2808 scalars.optimizationDirection = optimizationDirection_; 2809 memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double)); 2810 scalars.objectiveValue = objectiveValue_; 2811 scalars.dualBound = dualBound_; 2812 scalars.dualTolerance = dualTolerance_; 2813 scalars.primalTolerance = primalTolerance_; 2814 scalars.sumDualInfeasibilities = sumDualInfeasibilities_; 2815 scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_; 2816 scalars.infeasibilityCost = infeasibilityCost_; 2817 scalars.numberRows = numberRows_; 2818 scalars.numberColumns = numberColumns_; 2819 memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double)); 2820 scalars.numberIterations = numberIterations_; 2821 scalars.problemStatus = problemStatus_; 2822 scalars.maximumIterations = maximumIterations(); 2823 scalars.lengthNames = lengthNames_; 2824 scalars.numberDualInfeasibilities = numberDualInfeasibilities_; 2825 scalars.numberDualInfeasibilitiesWithoutFree 2826 = numberDualInfeasibilitiesWithoutFree_; 2827 scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_; 2828 scalars.numberRefinements = numberRefinements_; 2829 scalars.scalingFlag = scalingFlag_; 2830 scalars.algorithm = algorithm_; 2831 scalars.specialOptions = specialOptions_; 2832 scalars.dualPivotChoice = dualRowPivot_>type(); 2833 scalars.primalPivotChoice = primalColumnPivot_>type(); 2834 scalars.matrixStorageChoice = matrix_>type(); 2835 2836 // put out scalars 2837 numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp); 2838 if (numberWritten!=1) 2839 return 1; 2840 // strings 2841 CoinBigIndex length; 2842 for (i=0;i<ClpLastStrParam;i++) { 2843 length = strParam_[i].size(); 2844 numberWritten = fwrite(&length,sizeof(int),1,fp); 2845 if (numberWritten!=1) 2846 return 1; 2847 if (length) { 2848 numberWritten = fwrite(strParam_[i].c_str(),length,1,fp); 2849 if (numberWritten!=1) 2850 return 1; 2851 } 2852 } 2853 // arrays  in no particular order 2854 if (outDoubleArray(rowActivity_,numberRows_,fp)) 2855 return 1; 2856 if (outDoubleArray(columnActivity_,numberColumns_,fp)) 2857 return 1; 2858 if (outDoubleArray(dual_,numberRows_,fp)) 2859 return 1; 2860 if (outDoubleArray(reducedCost_,numberColumns_,fp)) 2861 return 1; 2862 if (outDoubleArray(rowLower_,numberRows_,fp)) 2863 return 1; 2864 if (outDoubleArray(rowUpper_,numberRows_,fp)) 2865 return 1; 2866 if (outDoubleArray(objective(),numberColumns_,fp)) 2867 return 1; 2868 if (outDoubleArray(rowObjective_,numberRows_,fp)) 2869 return 1; 2870 if (outDoubleArray(columnLower_,numberColumns_,fp)) 2871 return 1; 2872 if (outDoubleArray(columnUpper_,numberColumns_,fp)) 2873 return 1; 2874 if (ray_) { 2875 if (problemStatus_==1) 2876 if (outDoubleArray(ray_,numberRows_,fp)) 2877 return 1; 2878 else if (problemStatus_==2) 2879 if (outDoubleArray(ray_,numberColumns_,fp)) 2880 return 1; 2881 else 2882 if (outDoubleArray(NULL,0,fp)) 2883 return 1; 2884 } else { 2885 if (outDoubleArray(NULL,0,fp)) 2886 return 1; 2887 } 2888 if (status_&&(numberRows_+numberColumns_)>0) { 2889 length = numberRows_+numberColumns_; 2890 numberWritten = fwrite(&length,sizeof(int),1,fp); 2891 if (numberWritten!=1) 2892 return 1; 2893 numberWritten = fwrite(status_,sizeof(char),length, fp); 2894 if (numberWritten!=length) 2895 return 1; 2896 } else { 2897 length = 0; 2898 numberWritten = fwrite(&length,sizeof(int),1,fp); 2899 if (numberWritten!=1) 2900 return 1; 2901 } 2902 if (lengthNames_) { 2903 char * array = 2904 new char[max(numberRows_,numberColumns_)*(lengthNames_+1)]; 2905 char * put = array; 2906 assert (numberRows_ == (int) rowNames_.size()); 2907 for (i=0;i<numberRows_;i++) { 2908 assert((int)rowNames_[i].size()<=lengthNames_); 2909 strcpy(put,rowNames_[i].c_str()); 2910 put += lengthNames_+1; 2911 } 2912 numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp); 2913 if (numberWritten!=numberRows_) 2914 return 1; 2915 put=array; 2916 assert (numberColumns_ == (int) columnNames_.size()); 2917 for (i=0;i<numberColumns_;i++) { 2918 assert((int) columnNames_[i].size()<=lengthNames_); 2919 strcpy(put,columnNames_[i].c_str()); 2920 put += lengthNames_+1; 2921 } 2922 numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp); 2923 if (numberWritten!=numberColumns_) 2924 return 1; 2925 delete [] array; 2926 } 2927 // just standard type at present 2928 assert (matrix_>type()==1); 2929 assert (matrix_>getNumCols() == numberColumns_); 2930 assert (matrix_>getNumRows() == numberRows_); 2931 // we are going to save with gaps 2932 length = matrix_>getVectorStarts()[numberColumns_1] 2933 + matrix_>getVectorLengths()[numberColumns_1]; 2934 numberWritten = fwrite(&length,sizeof(int),1,fp); 2935 if (numberWritten!=1) 2936 return 1; 2937 numberWritten = fwrite(matrix_>getElements(), 2938 sizeof(double),length,fp); 2939 if (numberWritten!=length) 2940 return 1; 2941 numberWritten = fwrite(matrix_>getIndices(), 2942 sizeof(int),length,fp); 2943 if (numberWritten!=length) 2944 return 1; 2945 numberWritten = fwrite(matrix_>getVectorStarts(), 2946 sizeof(int),numberColumns_+1,fp); 2947 if (numberWritten!=numberColumns_+1) 2948 return 1; 2949 numberWritten = fwrite(matrix_>getVectorLengths(), 2950 sizeof(int),numberColumns_,fp); 2951 if (numberWritten!=numberColumns_) 2952 return 1; 2953 // finished 2954 fclose(fp); 2955 return 0; 2956 } else { 2957 return 1; 2958 } 2959 } 2960 2961 int inDoubleArray(double * &array, int length, FILE * fp) 2962 { 2963 int numberRead; 2964 int length2; 2965 numberRead = fread(&length2,sizeof(int),1,fp); 2966 if (numberRead!=1) 2967 return 1; 2968 if (length2) { 2969 // lengths must match 2970 if (length!=length2) 2971 return 2; 2972 array = new double[length]; 2973 numberRead = fread(array,sizeof(double),length,fp); 2974 if (numberRead!=length) 2975 return 1; 2976 } 2977 return 0; 2978 } 2979 /* Restore model from file, returns 0 if success, 2980 deletes current model */ 2981 int 2982 ClpInterior::restoreModel(const char * fileName) 2983 { 2984 FILE * fp = fopen(fileName,"rb"); 2985 if (fp) { 2986 // Get rid of current model 2987 ClpModel::gutsOfDelete(); 2988 gutsOfDelete(0); 2989 int i; 2990 for (i=0;i<6;i++) { 2991 rowArray_[i]=NULL; 2992 columnArray_[i]=NULL; 2993 } 2994 // get an empty factorization so we can set tolerances etc 2995 factorization_ = new ClpFactorization(); 2996 // Say sparse 2997 factorization_>sparseThreshold(1); 2998 Clp_scalars scalars; 2999 CoinBigIndex numberRead; 3000 3001 // get scalars 3002 numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp); 3003 if (numberRead!=1) 3004 return 1; 3005 // Fill in scalars 3006 optimizationDirection_ = scalars.optimizationDirection; 3007 memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double)); 3008 objectiveValue_ = scalars.objectiveValue; 3009 dualBound_ = scalars.dualBound; 3010 dualTolerance_ = scalars.dualTolerance; 3011 primalTolerance_ = scalars.primalTolerance; 3012 sumDualInfeasibilities_ = scalars.sumDualInfeasibilities; 3013 sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities; 3014 infeasibilityCost_ = scalars.infeasibilityCost; 3015 numberRows_ = scalars.numberRows; 3016 numberColumns_ = scalars.numberColumns; 3017 memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double)); 3018 numberIterations_ = scalars.numberIterations; 3019 problemStatus_ = scalars.problemStatus; 3020 setMaximumIterations(scalars.maximumIterations); 3021 lengthNames_ = scalars.lengthNames; 3022 numberDualInfeasibilities_ = scalars.numberDualInfeasibilities; 3023 numberDualInfeasibilitiesWithoutFree_ 3024 = scalars.numberDualInfeasibilitiesWithoutFree; 3025 numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities; 3026 numberRefinements_ = scalars.numberRefinements; 3027 scalingFlag_ = scalars.scalingFlag; 3028 algorithm_ = scalars.algorithm; 3029 specialOptions_ = scalars.specialOptions; 3030 // strings 3031 CoinBigIndex length; 3032 for (i=0;i<ClpLastStrParam;i++) { 3033 numberRead = fread(&length,sizeof(int),1,fp); 3034 if (numberRead!=1) 3035 return 1; 3036 if (length) { 3037 char * array = new char[length+1]; 3038 numberRead = fread(array,length,1,fp); 3039 if (numberRead!=1) 3040 return 1; 3041 array[length]='\0'; 3042 strParam_[i]=array; 3043 delete [] array; 3044 } 3045 } 3046 // arrays  in no particular order 3047 if (inDoubleArray(rowActivity_,numberRows_,fp)) 3048 return 1; 3049 if (inDoubleArray(columnActivity_,numberColumns_,fp)) 3050 return 1; 3051 if (inDoubleArray(dual_,numberRows_,fp)) 3052 return 1; 3053 if (inDoubleArray(reducedCost_,numberColumns_,fp)) 3054 return 1; 3055 if (inDoubleArray(rowLower_,numberRows_,fp)) 3056 return 1; 3057 if (inDoubleArray(rowUpper_,numberRows_,fp)) 3058 return 1; 3059 double * objective; 3060 if (inDoubleArray(objective,numberColumns_,fp)) 3061 return 1; 3062 delete objective_; 3063 objective_ = new ClpLinearObjective(objective,numberColumns_); 3064 delete [] objective; 3065 if (inDoubleArray(rowObjective_,numberRows_,fp)) 3066 return 1; 3067 if (inDoubleArray(columnLower_,numberColumns_,fp)) 3068 return 1; 3069 if (inDoubleArray(columnUpper_,numberColumns_,fp)) 3070 return 1; 3071 if (problemStatus_==1) { 3072 if (inDoubleArray(ray_,numberRows_,fp)) 3073 return 1; 3074 } else if (problemStatus_==2) { 3075 if (inDoubleArray(ray_,numberColumns_,fp)) 3076 return 1; 3077 } else { 3078 // ray should be null 3079 numberRead = fread(&length,sizeof(int),1,fp); 3080 if (numberRead!=1) 3081 return 1; 3082 if (length) 3083 return 2; 3084 } 3085 delete [] status_; 3086 status_=NULL; 3087 // status region 3088 numberRead = fread(&length,sizeof(int),1,fp); 3089 if (numberRead!=1) 3090 return 1; 3091 if (length) { 3092 if (length!=numberRows_+numberColumns_) 3093 return 1; 3094 status_ = new char unsigned[length]; 3095 numberRead = fread(status_,sizeof(char),length, fp); 3096 if (numberRead!=length) 3097 return 1; 3098 } 3099 if (lengthNames_) { 3100 char * array = 3101 new char[max(numberRows_,numberColumns_)*(lengthNames_+1)]; 3102 char * get = array; 3103 numberRead = fread(array,lengthNames_+1,numberRows_,fp); 3104 if (numberRead!=numberRows_) 3105 return 1; 3106 rowNames_ = std::vector<std::string> (); 3107 rowNames_.resize(numberRows_); 3108 for (i=0;i<numberRows_;i++) { 3109 rowNames_[i]=get; 3110 get += lengthNames_+1; 3111 } 3112 get = array; 3113 numberRead = fread(array,lengthNames_+1,numberColumns_,fp); 3114 if (numberRead!=numberColumns_) 3115 return 1; 3116 columnNames_ = std::vector<std::string> (); 3117 columnNames_.resize(numberColumns_); 3118 for (i=0;i<numberColumns_;i++) { 3119 columnNames_[i]=get; 3120 get += lengthNames_+1; 3121 } 3122 delete [] array; 3123 } 3124 // Pivot choices 3125 assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3); 3126 delete dualRowPivot_; 3127 switch ((scalars.dualPivotChoice&63)) { 3128 default: 3129 printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63); 3130 case 1: 3131 // Dantzig 3132 dualRowPivot_ = new ClpDualRowDantzig(); 3133 break; 3134 case 2: 3135 // Steepest  use mode 3136 dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6); 3137 break; 3138 } 3139 assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3); 3140 delete primalColumnPivot_; 3141 switch ((scalars.primalPivotChoice&63)) { 3142 default: 3143 printf("Need another primalPivot case %d\n", 3144 scalars.primalPivotChoice&63); 3145 case 1: 3146 // Dantzig 3147 primalColumnPivot_ = new ClpPrimalColumnDantzig(); 3148 break; 3149 case 2: 3150 // Steepest  use mode 3151 primalColumnPivot_ 3152 = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6); 3153 break; 3154 } 3155 assert(scalars.matrixStorageChoice==1); 3156 delete matrix_; 3157 // get arrays 3158 numberRead = fread(&length,sizeof(int),1,fp); 3159 if (numberRead!=1) 3160 return 1; 3161 double * elements = new double[length]; 3162 int * indices = new int[length]; 3163 CoinBigIndex * starts = new CoinBigIndex[numberColumns_]; 3164 int * lengths = new int[numberColumns_]; 3165 numberRead = fread(elements, sizeof(double),length,fp); 3166 if (numberRead!=length) 3167 return 1; 3168 numberRead = fread(indices, sizeof(int),length,fp); 3169 if (numberRead!=length) 3170 return 1; 3171 numberRead = fread(starts, sizeof(int),numberColumns_+1,fp); 3172 if (numberRead!=numberColumns_+1) 3173 return 1; 3174 numberRead = fread(lengths, sizeof(int),numberColumns_,fp); 3175 if (numberRead!=numberColumns_) 3176 return 1; 3177 // assign matrix 3178 CoinPackedMatrix * matrix = new CoinPackedMatrix(); 3179 matrix>assignMatrix(true, numberRows_, numberColumns_, 3180 length, elements, indices, starts, lengths); 3181 // and transfer to Clp 3182 matrix_ = new ClpPackedMatrix(matrix); 3183 // finished 3184 fclose(fp); 3185 return 0; 3186 } else { 3187 return 1; 3188 } 3189 return 0; 3190 } 3191 // value of incoming variable (in Dual) 3192 double 3193 ClpInterior::valueIncomingDual() const 3194 { 3195 // Need value of incoming for list of infeasibilities as may be infeasible 3196 double valueIncoming = (dualOut_/alpha_)*directionOut_; 3197 if (directionIn_==1) 3198 valueIncoming = upperIn_valueIncoming; 3199 else 3200 valueIncoming = lowerIn_valueIncoming; 3201 return valueIncoming; 175 ClpInterior::deleteWorkingData() 176 { 3202 177 } 3203 178 // Sanity check on input data  returns true if okay … … 3229 204 largestObj=0.0; 3230 205 // If bounds are too close  fix 3231 double fixTolerance = 1.1*primalTolerance _;206 double fixTolerance = 1.1*primalTolerance(); 3232 207 for (i=numberColumns_;i<numberColumns_+numberRows_;i++) { 3233 208 double value; … … 3244 219 } 3245 220 value=upper_[i]lower_[i]; 3246 if (value<primalTolerance _) {221 if (value<primalTolerance()) { 3247 222 numberBad++; 3248 223 if (firstBad<0) … … 3296 271 } 3297 272 value=upper_[i]lower_[i]; 3298 if (value<primalTolerance _) {273 if (value<primalTolerance()) { 3299 274 numberBad++; 3300 275 if (firstBad<0) … … 3349 324 return true; 3350 325 } 3351 // Set up status array (for OsiClp)3352 void3353 ClpInterior::createStatus()3354 {3355 if(!status_)3356 status_ = new unsigned char [numberColumns_+numberRows_];3357 memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));3358 int i;3359 // set column status to one nearest zero3360 for (i=0;i<numberColumns_;i++) {3361 #if 03362 if (columnLower_[i]>=0.0) {3363 setColumnStatus(i,atLowerBound);3364 } else if (columnUpper_[i]<=0.0) {3365 setColumnStatus(i,atUpperBound);3366 } else if (columnLower_[i]<1.0e20&&columnUpper_[i]>1.0e20) {3367 // free3368 setColumnStatus(i,isFree);3369 } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {3370 setColumnStatus(i,atLowerBound);3371 } else {3372 setColumnStatus(i,atUpperBound);3373 }3374 #else3375 setColumnStatus(i,atLowerBound);3376 #endif3377 }3378 for (i=0;i<numberRows_;i++) {3379 setRowStatus(i,basic);3380 }3381 }3382 /* Loads a problem (the constraints on the3383 rows are given by lower and upper bounds). If a pointer is 0 then the3384 following values are the default:3385 <ul>3386 <li> <code>colub</code>: all columns have upper bound infinity3387 <li> <code>collb</code>: all columns have lower bound 03388 <li> <code>rowub</code>: all rows have upper bound infinity3389 <li> <code>rowlb</code>: all rows have lower bound infinity3390 <li> <code>obj</code>: all variables have 0 objective coefficient3391 </ul>3392 */3393 void3394 ClpInterior::loadProblem ( const ClpMatrixBase& matrix,3395 const double* collb, const double* colub,3396 const double* obj,3397 const double* rowlb, const double* rowub,3398 const double * rowObjective)3399 {3400 ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,3401 rowObjective);3402 createStatus();3403 }3404 void3405 ClpInterior::loadProblem ( const CoinPackedMatrix& matrix,3406 const double* collb, const double* colub,3407 const double* obj,3408 const double* rowlb, const double* rowub,3409 const double * rowObjective)3410 {3411 ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,3412 rowObjective);3413 createStatus();3414 }3415 3416 /* Just like the other loadProblem() method except that the matrix is3417 given in a standard column major ordered format (without gaps). */3418 void3419 ClpInterior::loadProblem ( const int numcols, const int numrows,3420 const CoinBigIndex* start, const int* index,3421 const double* value,3422 const double* collb, const double* colub,3423 const double* obj,3424 const double* rowlb, const double* rowub,3425 const double * rowObjective)3426 {3427 ClpModel::loadProblem(numcols, numrows, start, index, value,3428 collb, colub, obj, rowlb, rowub,3429 rowObjective);3430 createStatus();3431 }3432 void3433 ClpInterior::loadProblem ( const int numcols, const int numrows,3434 const CoinBigIndex* start, const int* index,3435 const double* value,const int * length,3436 const double* collb, const double* colub,3437 const double* obj,3438 const double* rowlb, const double* rowub,3439 const double * rowObjective)3440 {3441 ClpModel::loadProblem(numcols, numrows, start, index, value, length,3442 collb, colub, obj, rowlb, rowub,3443 rowObjective);3444 createStatus();3445 }3446 // Read an mps file from the given filename3447 int3448 ClpInterior::readMps(const char *filename,3449 bool keepNames,3450 bool ignoreErrors)3451 {3452 int status = ClpModel::readMps(filename,keepNames,ignoreErrors);3453 createStatus();3454 return status;3455 }3456 // Just check solution (for external use)3457 void3458 ClpInterior::checkSolution()3459 {3460 // put in standard form3461 createRim(7+8+16);3462 dualTolerance_=dblParam_[ClpDualTolerance];3463 primalTolerance_=dblParam_[ClpPrimalTolerance];3464 checkPrimalSolution( rowActivityWork_, columnActivityWork_);3465 checkDualSolution();3466 if (!numberDualInfeasibilities_&&3467 !numberPrimalInfeasibilities_)3468 problemStatus_=0;3469 else3470 problemStatus_=1;3471 #ifdef CLP_DEBUG3472 int i;3473 double value=0.0;3474 for (i=0;i<numberRows_+numberColumns_;i++)3475 value += dj_[i]*solution_[i];3476 printf("dual value %g, primal %g\n",value,objectiveValue());3477 #endif3478 // release extra memory3479 deleteRim(0);3480 }3481 /* Crash  at present just aimed at dual, returns3482 2 if dual preferred and crash basis created3483 1 if dual preferred and all slack basis preferred3484 0 if basis going in was not all slack3485 1 if primal preferred and all slack basis preferred3486 2 if primal preferred and crash basis created.3487 3488 if gap between bounds <="gap" variables can be flipped3489 3490 If "pivot" is3491 0 No pivoting (so will just be choice of algorithm)3492 1 Simple pivoting e.g. gub3493 2 Mini iterations3494 */3495 int3496 ClpInterior::crash(double gap,int pivot)3497 {3498 assert(!rowObjective_); // not coded3499 int iColumn;3500 int numberBad=0;3501 int numberBasic=0;3502 double dualTolerance=dblParam_[ClpDualTolerance];3503 //double primalTolerance=dblParam_[ClpPrimalTolerance];3504 3505 // If no basis then make all slack one3506 if (!status_)3507 createStatus();3508 3509 for (iColumn=0;iColumn<numberColumns_;iColumn++) {3510 if (getColumnStatus(iColumn)==basic)3511 numberBasic++;3512 }3513 if (numberBasic) {3514 // not all slack3515 return 0;3516 } else {3517 double * dj = new double [numberColumns_];3518 double * solution = columnActivity_;3519 const double * linearObjective = objective();3520 //double objectiveValue=0.0;3521 int iColumn;3522 double direction = optimizationDirection_;3523 // direction is actually scale out not scale in3524 if (direction)3525 direction = 1.0/direction;3526 for (iColumn=0;iColumn<numberColumns_;iColumn++)3527 dj[iColumn] = direction*linearObjective[iColumn];3528 for (iColumn=0;iColumn<numberColumns_;iColumn++) {3529 // assume natural place is closest to zero3530 double lowerBound = columnLower_[iColumn];3531 double upperBound = columnUpper_[iColumn];3532 if (lowerBound>1.0e20upperBound<1.0e20) {3533 bool atLower;3534 if (fabs(upperBound)<fabs(lowerBound)) {3535 atLower=false;3536 setColumnStatus(iColumn,atUpperBound);3537 solution[iColumn]=upperBound;3538 } else {3539 atLower=true;3540 setColumnStatus(iColumn,atLowerBound);3541 solution[iColumn]=lowerBound;3542 }3543 if (dj[iColumn]<0.0) {3544 // should be at upper bound3545 if (atLower) {3546 // can we flip3547 if (upperBoundlowerBound<=gap) {3548 columnActivity_[iColumn]=upperBound;3549 setColumnStatus(iColumn,atUpperBound);3550 } else if (dj[iColumn]<dualTolerance) {3551 numberBad++;3552 }3553 }3554 } else if (dj[iColumn]>0.0) {3555 // should be at lower bound3556 if (!atLower) {3557 // can we flip3558 if (upperBoundlowerBound<=gap) {3559 columnActivity_[iColumn]=lowerBound;3560 setColumnStatus(iColumn,atLowerBound);3561 } else if (dj[iColumn]>dualTolerance) {3562 numberBad++;3563 }3564 }3565 }3566 } else {3567 // free3568 setColumnStatus(iColumn,isFree);3569 if (fabs(dj[iColumn])>dualTolerance)3570 numberBad++;3571 }3572 }3573 if (numberBadpivot) {3574 if (!pivot) {3575 delete [] dj;3576 return 1;3577 } else {3578 // see if can be made dual feasible with gubs etc3579 double * pi = new double[numberRows_];3580 memset (pi,0,numberRows_*sizeof(double));3581 int * way = new int[numberColumns_];3582 int numberIn = 0;3583 3584 // Get column copy3585 CoinPackedMatrix * columnCopy = matrix();3586 // Get a row copy in standard format3587 CoinPackedMatrix copy;3588 copy.reverseOrderedCopyOf(*columnCopy);3589 // get matrix data pointers3590 const int * column = copy.getIndices();3591 const CoinBigIndex * rowStart = copy.getVectorStarts();3592 const int * rowLength = copy.getVectorLengths();3593 const double * elementByRow = copy.getElements();3594 //const int * row = columnCopy>getIndices();3595 //const CoinBigIndex * columnStart = columnCopy>getVectorStarts();3596 //const int * columnLength = columnCopy>getVectorLengths();3597 //const double * element = columnCopy>getElements();3598 3599 3600 // if equality row and bounds mean artificial in basis bad3601 // then do anyway3602 3603 for (iColumn=0;iColumn<numberColumns_;iColumn++) {3604 //  if we want to reduce dj, + if we want to increase3605 int thisWay = 100;3606 double lowerBound = columnLower_[iColumn];3607 double upperBound = columnUpper_[iColumn];3608 if (upperBound>lowerBound) {3609 switch(getColumnStatus(iColumn)) {3610 3611 case basic:3612 thisWay=0;3613 case ClpInterior::isFixed:3614 break;3615 case isFree:3616 case superBasic:3617 if (dj[iColumn]<dualTolerance)3618 thisWay = 1;3619 else if (dj[iColumn]>dualTolerance)3620 thisWay = 1;3621 else3622 thisWay =0;3623 break;3624 case atUpperBound:3625 if (dj[iColumn]>dualTolerance)3626 thisWay = 1;3627 else if (dj[iColumn]<dualTolerance)3628 thisWay = 3;3629 else3630 thisWay = 2;3631 break;3632 case atLowerBound:3633 if (dj[iColumn]<dualTolerance)3634 thisWay = 1;3635 else if (dj[iColumn]>dualTolerance)3636 thisWay = 3;3637 else3638 thisWay = 2;3639 break;3640 }3641 }3642 way[iColumn] = thisWay;3643 }3644 /*if (!numberBad)3645 printf("Was dual feasible before passes  rows %d\n",3646 numberRows_);*/3647 int lastNumberIn = 100000;3648 int numberPasses=5;3649 while (numberIn>lastNumberIn+numberRows_/100) {3650 lastNumberIn = numberIn;3651 // we need to maximize chance of doing good3652 int iRow;3653 for (iRow=0;iRow<numberRows_;iRow++) {3654 double lowerBound = rowLower_[iRow];3655 double upperBound = rowUpper_[iRow];3656 if (getRowStatus(iRow)==basic) {3657 // see if we can find a column to pivot on3658 int j;3659 // down is amount pi can go down3660 double maximumDown = COIN_DBL_MAX;3661 double maximumUp = COIN_DBL_MAX;3662 double minimumDown =0.0;3663 double minimumUp =0.0;3664 int iUp=1;3665 int iDown=1;3666 int iUpB=1;3667 int iDownB=1;3668 if (lowerBound<1.0e20)3669 maximumUp = 1.0;3670 if (upperBound>1.0e20)3671 maximumDown = 1.0;3672 for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {3673 int iColumn = column[j];3674 double value = elementByRow[j];3675 double djValue = dj[iColumn];3676 /* way 3677 3  okay at upper bound with negative dj3678 2  marginal at upper bound with zero dj  can only decrease3679 1  bad at upper bound3680 0  we can never pivot on this row3681 1  bad at lower bound3682 2  marginal at lower bound with zero dj  can only increase3683 3  okay at lower bound with positive dj3684 100  fine we can just ignore3685 */3686 if (way[iColumn]!=100) {3687 switch(way[iColumn]) {3688 3689 case 3:3690 if (value>0.0) {3691 if (maximumDown*value>djValue) {3692 maximumDown = djValue/value;3693 iDown = iColumn;3694 }3695 } else {3696 if (maximumUp*value>djValue) {3697 maximumUp = djValue/value;3698 iUp = iColumn;3699 }3700 }3701 break;3702 case 2:3703 if (value>0.0) {3704 maximumDown = 0.0;3705 } else {3706 maximumUp = 0.0;3707 }3708 break;3709 case 1:3710 // see if could be satisfied3711 // dj value > 03712 if (value>0.0) {3713 maximumDown=0.0;3714 if (maximumUp*value<djValuedualTolerance) {3715 maximumUp = 0.0; // would improve but not enough3716 } else {3717 if (minimumUp*value<djValue) {3718 minimumUp = djValue/value;3719 iUpB = iColumn;3720 }3721 }3722 } else {3723 maximumUp=0.0;3724 if (maximumDown*value<djValuedualTolerance) {3725 maximumDown = 0.0; // would improve but not enough3726 } else {3727 if (minimumDown*value<djValue) {3728 minimumDown = djValue/value;3729 iDownB = iColumn;3730 }3731 }3732 }3733 3734 break;3735 case 0:3736 maximumDown = 1.0;3737 maximumUp=1.0;3738 break;3739 case 1:3740 // see if could be satisfied3741 // dj value < 03742 if (value>0.0) {3743 maximumUp=0.0;3744 if (maximumDown*value<djValuedualTolerance) {3745 maximumDown = 0.0; // would improve but not enough3746 } else {3747 if (minimumDown*value<djValue) {3748 minimumDown = djValue/value;3749 iDownB = iColumn;3750 }3751 }3752 } else {3753 maximumDown=0.0;3754 if (maximumUp*value<djValuedualTolerance) {3755 maximumUp = 0.0; // would improve but not enough3756 } else {3757 if (minimumUp*value<djValue) {3758 minimumUp = djValue/value;3759 iUpB = iColumn;3760 }3761 }3762 }3763 3764 break;3765 case 2:3766 if (value>0.0) {3767 maximumUp = 0.0;3768 } else {3769 maximumDown = 0.0;3770 }3771 3772 break;3773 case 3:3774 if (value>0.0) {3775 if (maximumUp*value>djValue) {3776 maximumUp = djValue/value;3777 iUp = iColumn;3778 }3779 } else {3780 if (maximumDown*value>djValue) {3781 maximumDown = djValue/value;3782 iDown = iColumn;3783 }3784 }3785 3786 break;3787 default:3788 break;3789 }3790 }3791 }3792 if (iUpB>=0)3793 iUp=iUpB;3794 if (maximumUp<=dualTolerancemaximumUp<minimumUp)3795 iUp=1;3796 if (iDownB>=0)3797 iDown=iDownB;3798 if (maximumDown<=dualTolerancemaximumDown<minimumDown)3799 iDown=1;3800 if (iUp>=0iDown>=0) {3801 // do something3802 if (iUp>=0&&iDown>=0) {3803 if (maximumDown>maximumUp)3804 iUp=1;3805 }3806 double change;3807 int kColumn;3808 if (iUp>=0) {3809 kColumn=iUp;3810 change=maximumUp;3811 // just do minimum if was dual infeasible3812 // ? only if maximum large?3813 if (minimumUp>0.0)3814 change=minimumUp;3815 setRowStatus(iRow,atUpperBound);3816 } else {3817 kColumn=iDown;3818 change=maximumDown;3819 // just do minimum if was dual infeasible3820 // ? only if maximum large?3821 if (minimumDown>0.0)3822 change=minimumDown;3823 setRowStatus(iRow,atLowerBound);3824 }3825 assert (fabs(change)<1.0e20);3826 setColumnStatus(kColumn,basic);3827 numberIn++;3828 pi[iRow]=change;3829 for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {3830 int iColumn = column[j];3831 double value = elementByRow[j];3832 double djValue = dj[iColumn]change*value;3833 dj[iColumn]=djValue;3834 if (abs(way[iColumn])==1) {3835 numberBad;3836 /*if (!numberBad)3837 printf("Became dual feasible at row %d out of %d\n",3838 iRow, numberRows_);*/3839 lastNumberIn=1000000;3840 }3841 int thisWay = 100;3842 double lowerBound = columnLower_[iColumn];3843 double upperBound = columnUpper_[iColumn];3844 if (upperBound>lowerBound) {3845 switch(getColumnStatus(iColumn)) {3846 3847 case basic:3848 thisWay=0;3849 case isFixed:3850 break;3851 case isFree:3852 case superBasic:3853 if (djValue<dualTolerance)3854 thisWay = 1;3855 else if (djValue>dualTolerance)3856 thisWay = 1;3857 else3858 { thisWay =0; abort();}3859 break;3860 case atUpperBound:3861 if (djValue>dualTolerance)3862 { thisWay =1; abort();}3863 else if (djValue<dualTolerance)3864 thisWay = 3;3865 else3866 thisWay = 2;3867 break;3868 case atLowerBound:3869 if (djValue<dualTolerance)3870 { thisWay =1; abort();}3871 else if (djValue>dualTolerance)3872 thisWay = 3;3873 else3874 thisWay = 2;3875 break;3876 }3877 }3878 way[iColumn] = thisWay;3879 }3880 }3881 }3882 }3883 if (numberIn==lastNumberInnumberBadpivot<2)3884 break;3885 if (!(numberPasses))3886 break;3887 //printf("%d put in so far\n",numberIn);3888 }3889 // last attempt to flip3890 for (iColumn=0;iColumn<numberColumns_;iColumn++) {3891 double lowerBound = columnLower_[iColumn];3892 double upperBound = columnUpper_[iColumn];3893 if (upperBoundlowerBound<=gap&&upperBound>lowerBound) {3894 double djValue=dj[iColumn];3895 switch(getColumnStatus(iColumn)) {3896 3897 case basic:3898 case ClpInterior::isFixed:3899 break;3900 case isFree:3901 case superBasic:3902 break;3903 case atUpperBound:3904 if (djValue>dualTolerance) {3905 setColumnStatus(iColumn,atUpperBound);3906 solution[iColumn]=upperBound;3907 }3908 break;3909 case atLowerBound:3910 if (djValue<dualTolerance) {3911 setColumnStatus(iColumn,atUpperBound);3912 solution[iColumn]=upperBound;3913 }3914 break;3915 }3916 }3917 }3918 delete [] pi;3919 delete [] dj;3920 delete [] way;3921 handler_>message(CLP_CRASH,messages_)3922 <<numberIn3923 <<numberBad3924 <<CoinMessageEol;3925 return 1;3926 }3927 } else {3928 delete [] dj;3929 return 1;3930 }3931 }3932 }3933 /* Pivot in a variable and out a variable. Returns 0 if okay,3934 1 if inaccuracy forced refactorization, 1 if would be singular.3935 Also updates primal/dual infeasibilities.3936 Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.3937 */3938 int ClpInterior::pivot()3939 {3940 // scaling not allowed3941 assert (!scalingFlag_);3942 // assume In_ and Out_ are correct and directionOut_ set3943 // (or In_ if flip3944 lowerIn_ = lower_[sequenceIn_];3945 valueIn_ = solution_[sequenceIn_];3946 upperIn_ = upper_[sequenceIn_];3947 dualIn_ = dj_[sequenceIn_];3948 if (sequenceOut_>=0&&sequenceIn_!=sequenceIn_) {3949 assert (pivotRow_>=0&&pivotRow_<numberRows_);3950 assert (pivotVariable_[pivotRow_]==sequenceOut_);3951 lowerOut_ = lower_[sequenceOut_];3952 valueOut_ = solution_[sequenceOut_];3953 upperOut_ = upper_[sequenceOut_];3954 // for now assume primal is feasible (or in dual)3955 dualOut_ = dj_[sequenceOut_];3956 assert(fabs(dualOut_)<1.0e6);3957 } else {3958 assert (pivotRow_<0);3959 }3960 bool roundAgain = true;3961 int returnCode=0;3962 while (roundAgain) {3963 roundAgain=false;3964 unpack(rowArray_[1]);3965 factorization_>updateColumnFT(rowArray_[2],rowArray_[1]);3966 // we are going to subtract movement from current basic3967 double movement;3968 // see where incoming will go to3969 if (sequenceOut_<0sequenceIn_==sequenceOut_) {3970 // flip so go to bound3971 movement = ((directionIn_>0) ? upperIn_ : lowerIn_)  valueIn_;3972 } else {3973 // get where outgoing needs to get to3974 double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;3975 // solutionOut_  movement*alpha_ == outValue3976 movement = (outValuevalueOut_)/alpha_;3977 // set directionIn_ correctly3978 directionIn_ = (movement>0) ? 1 :1;3979 }3980 // update primal solution3981 {3982 int i;3983 int * index = rowArray_[1]>getIndices();3984 int number = rowArray_[1]>getNumElements();3985 double * element = rowArray_[1]>denseVector();3986 for (i=0;i<number;i++) {3987 int ii = index[i];3988 // get column3989 ii = pivotVariable_[ii];3990 solution_[ii] = movement*element[i];3991 element[i]=0.0;3992 }3993 // see where something went to3994 if (sequenceOut_<0) {3995 if (directionIn_<0) {3996 assert (fabs(solution_[sequenceIn_]upperIn_)<1.0e7);3997 solution_[sequenceIn_]=upperIn_;3998 } else {3999 assert (fabs(solution_[sequenceIn_]lowerIn_)<1.0e7);4000 solution_[sequenceIn_]=lowerIn_;4001 }4002 } else {4003 if (directionOut_<0) {4004 assert (fabs(solution_[sequenceOut_]upperOut_)<1.0e7);4005 solution_[sequenceOut_]=upperOut_;4006 } else {4007 assert (fabs(solution_[sequenceOut_]lowerOut_)<1.0e7);4008 solution_[sequenceOut_]=lowerOut_;4009 }4010 solution_[sequenceIn_]=valueIn_+movement;4011 }4012 }4013 double objectiveChange = dualIn_*movement;4014 // update duals4015 if (pivotRow_>=0) {4016 alpha_ = rowArray_[1]>denseVector()[pivotRow_];4017 assert (fabs(alpha_)>1.0e8);4018 double multiplier = dualIn_/alpha_;4019 rowArray_[0]>insert(pivotRow_,multiplier);4020 factorization_>updateColumnTranspose(rowArray_[2],rowArray_[0]);4021 // put row of tableau in rowArray[0] and columnArray[0]4022 matrix_>transposeTimes(this,1.0,4023 rowArray_[0],columnArray_[1],columnArray_[0]);4024 // update column djs4025 int i;4026 int * index = columnArray_[0]>getIndices();4027 int number = columnArray_[0]>getNumElements();4028 double * element = columnArray_[0]>denseVector();4029 for (i=0;i<number;i++) {4030 int ii = index[i];4031 dj_[ii] += element[ii];4032 element[ii]=0.0;4033 }4034 columnArray_[0]>setNumElements(0);4035 // and row djs4036 index = rowArray_[0]>getIndices();4037 number = rowArray_[0]>getNumElements();4038 element = rowArray_[0]>denseVector();4039 for (i=0;i<number;i++) {4040 int ii = index[i];4041 dj_[ii+numberColumns_] += element[ii];4042 dual_[ii] = dj_[ii+numberColumns_];4043 element[ii]=0.0;4044 }4045 rowArray_[0]>setNumElements(0);4046 // check incoming4047 assert (fabs(dj_[sequenceIn_])<1.0e6);4048 }4049 4050 // if stable replace in basis4051 int updateStatus = factorization_>replaceColumn(rowArray_[2],4052 pivotRow_,4053 alpha_);4054 bool takePivot=true;4055 // if no pivots, bad update but reasonable alpha  take and invert4056 if (updateStatus==2&&4057 lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e5)4058 updateStatus=4;4059 if (updateStatus==1updateStatus==4) {4060 // slight error4061 if (factorization_>pivots()>5updateStatus==4) {4062 returnCode=1;4063 }4064 } else if (updateStatus==2) {4065 // major error4066 rowArray_[1]>clear();4067 takePivot=false;4068 if (factorization_>pivots()) {4069 // refactorize here4070 statusOfProblem();4071 roundAgain=true;4072 } else {4073 returnCode=1;4074 }4075 } else if (updateStatus==3) {4076 // out of memory4077 // increase space if not many iterations4078 if (factorization_>pivots()<4079 0.5*factorization_>maximumPivots()&&4080 factorization_>pivots()<200)4081 factorization_>areaFactor(4082 factorization_>areaFactor() * 1.1);4083 returnCode =1; // factorize now4084 }4085 if (takePivot) {4086 int save = algorithm_;4087 // make simple so always primal4088 algorithm_=1;4089 housekeeping(objectiveChange);4090 algorithm_=save;4091 }4092 }4093 if (returnCode == 1) {4094 // refactorize here4095 statusOfProblem();4096 } else {4097 // just for now  recompute anyway4098 gutsOfSolution(NULL,NULL);4099 }4100 return returnCode;4101 }4102 4103 /* Pivot in a variable and choose an outgoing one. Assumes primal4104 feasible  will not go through a bound. Returns step length in theta4105 Returns ray in ray_ (or NULL if no pivot)4106 Return codes as before but 1 means no acceptable pivot4107 */4108 int ClpInterior::primalPivotResult()4109 {4110 assert (sequenceIn_>=0);4111 valueIn_=solution_[sequenceIn_];4112 lowerIn_=lower_[sequenceIn_];4113 upperIn_=upper_[sequenceIn_];4114 dualIn_=dj_[sequenceIn_];4115 4116 int returnCode = ((ClpInteriorPrimal *) this)>pivotResult();4117 if (returnCode<0&&returnCode>4) {4118 return 0;4119 } else {4120 printf("Return code of %d from ClpInteriorPrimal::pivotResult\n",4121 returnCode);4122 return 1;4123 }4124 }4125 4126 /* Pivot out a variable and choose an incoing one. Assumes dual4127 feasible  will not go through a reduced cost.4128 Returns step length in theta4129 Returns ray in ray_ (or NULL if no pivot)4130 Return codes as before but 1 means no acceptable pivot4131 */4132 int4133 ClpInterior::dualPivotResult()4134 {4135 return ((ClpInteriorDual *) this)>pivotResult();4136 }4137 // Factorization frequency4138 int4139 ClpInterior::factorizationFrequency() const4140 {4141 if (factorization_)4142 return factorization_>maximumPivots();4143 else4144 return 1;4145 }4146 void4147 ClpInterior::setFactorizationFrequency(int value)4148 {4149 if (factorization_)4150 factorization_>maximumPivots(value);4151 }4152 // Common bits of coding for dual and primal4153 int4154 ClpInterior::startup(int ifValuesPass)4155 {4156 // sanity check4157 assert (numberRows_==matrix_>getNumRows());4158 assert (numberColumns_==matrix_>getNumCols());4159 // for moment all arrays must exist4160 assert(columnLower_);4161 assert(columnUpper_);4162 assert(rowLower_);4163 assert(rowUpper_);4164 pivotRow_=1;4165 sequenceIn_=1;4166 sequenceOut_=1;4167 secondaryStatus_=0;4168 4169 primalTolerance_=dblParam_[ClpPrimalTolerance];4170 dualTolerance_=dblParam_[ClpDualTolerance];4171 if (problemStatus_!=10)4172 numberIterations_=0;4173 4174 // put in standard form (and make row copy)4175 // create modifiable copies of model rim and do optional scaling4176 bool goodMatrix=createRim(7+8+16,true);4177 4178 if (goodMatrix) {4179 // Model looks okay4180 // Do initial factorization4181 // and set certain stuff4182 // We can either set increasing rows so ...IsBasic gives pivot row4183 // or we can just increment iBasic one by one4184 // for now let ...iBasic give pivot row4185 factorization_>increasingRows(2);4186 // row activities have negative sign4187 factorization_>slackValue(1.0);4188 factorization_>zeroTolerance(1.0e13);4189 // Switch off dense (unless special option set)4190 int saveThreshold = factorization_>denseThreshold();4191 factorization_>setDenseThreshold(0);4192 // If values pass then perturb (otherwise may be optimal so leave a bit)4193 if (ifValuesPass) {4194 // do perturbation if asked for4195 4196 if (perturbation_<100) {4197 if (algorithm_>0) {4198 ((ClpInteriorPrimal *) this)>perturb(0);4199 } else if (algorithm_<0) {4200 ((ClpInteriorDual *) this)>perturb();4201 }4202 }4203 }4204 // for primal we will change bounds using infeasibilityCost_4205 if (nonLinearCost_==NULL&&algorithm_>0) {4206 // get a valid nonlinear cost function4207 delete nonLinearCost_;4208 nonLinearCost_= new ClpNonLinearCost(this);4209 }4210 4211 // loop round to clean up solution if values pass4212 int numberThrownOut = 1;4213 int totalNumberThrownOut=0;4214 while(numberThrownOut) {4215 int status = internalFactorize(0+10*ifValuesPass);4216 if (status<0)4217 return 1; // some error4218 else4219 numberThrownOut = status;4220 4221 // for this we need clean basis so it is after factorize4222 if (!numberThrownOut)4223 numberThrownOut = gutsOfSolution( NULL,NULL,4224 ifValuesPass!=0);4225 totalNumberThrownOut+= numberThrownOut;4226 4227 }4228 4229 if (totalNumberThrownOut)4230 handler_>message(CLP_SINGULARITIES,messages_)4231 <<totalNumberThrownOut4232 <<CoinMessageEol;4233 // Switch back dense4234 factorization_>setDenseThreshold(saveThreshold);4235 4236 problemStatus_ = 1;4237 4238 // number of times we have declared optimality4239 numberTimesOptimal_=0;4240 4241 return 0;4242 } else {4243 // bad matrix4244 return 2;4245 }4246 4247 }4248 4249 4250 void4251 ClpInterior::finish()4252 {4253 // Get rid of some arrays and empty factorization4254 deleteRim();4255 // Skip message if changing algorithms4256 if (problemStatus_!=10) {4257 assert(problemStatus_>=0&&problemStatus_<5);4258 handler_>message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)4259 <<objectiveValue()4260 <<CoinMessageEol;4261 }4262 factorization_>relaxAccuracyCheck(1.0);4263 // get rid of any network stuff  could do more4264 factorization_>cleanUp();4265 }4266 // Save data4267 ClpDataSave4268 ClpInterior::saveData()4269 {4270 ClpDataSave saved;4271 saved.dualBound_ = dualBound_;4272 saved.infeasibilityCost_ = infeasibilityCost_;4273 saved.sparseThreshold_ = factorization_>sparseThreshold();4274 saved.perturbation_ = perturbation_;4275 // Progress indicator4276 delete progress_;4277 progress_ = new ClpInteriorProgress (this);4278 return saved;4279 }4280 // Restore data4281 void4282 ClpInterior::restoreData(ClpDataSave saved)4283 {4284 factorization_>sparseThreshold(saved.sparseThreshold_);4285 perturbation_ = saved.perturbation_;4286 infeasibilityCost_ = saved.infeasibilityCost_;4287 dualBound_ = saved.dualBound_;4288 delete progress_;4289 progress_=NULL;4290 }4291 /* Factorizes and returns true if optimal. Used by user */4292 bool4293 ClpInterior::statusOfProblem()4294 {4295 // is factorization okay?4296 assert (internalFactorize(1)==0);4297 // put back original costs and then check4298 // also move to work arrays4299 createRim(4+32);4300 //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));4301 //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));4302 gutsOfSolution(NULL,NULL);4303 //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));4304 //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));4305 //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));4306 deleteRim(1);4307 return (primalFeasible()&&dualFeasible());4308 }4309 /* Return model  updates any scalars */4310 void4311 ClpInterior::returnModel(ClpInterior & otherModel)4312 {4313 ClpModel::returnModel(otherModel);4314 otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;4315 otherModel.columnPrimalSequence_ = columnPrimalSequence_;4316 otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;4317 otherModel.rowPrimalSequence_ = rowPrimalSequence_;4318 otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;4319 otherModel.columnDualSequence_ = columnDualSequence_;4320 otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;4321 otherModel.rowDualSequence_ = rowDualSequence_;4322 otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;4323 otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;4324 otherModel.largestPrimalError_ = largestPrimalError_;4325 otherModel.largestDualError_ = largestDualError_;4326 otherModel.largestSolutionError_ = largestSolutionError_;4327 otherModel.alpha_ = alpha_;4328 otherModel.theta_ = theta_;4329 otherModel.lowerIn_ = lowerIn_;4330 otherModel.valueIn_ = valueIn_;4331 otherModel.upperIn_ = upperIn_;4332 otherModel.dualIn_ = dualIn_;4333 otherModel.sequenceIn_ = sequenceIn_;4334 otherModel.directionIn_ = directionIn_;4335 otherModel.lowerOut_ = lowerOut_;4336 otherModel.valueOut_ = valueOut_;4337 otherModel.upperOut_ = upperOut_;4338 otherModel.dualOut_ = dualOut_;4339 otherModel.sequenceOut_ = sequenceOut_;4340 otherModel.directionOut_ = directionOut_;4341 otherModel.pivotRow_ = pivotRow_;4342 otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;4343 otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;4344 otherModel.numberDualInfeasibilitiesWithoutFree_ =4345 numberDualInfeasibilitiesWithoutFree_;4346 otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;4347 otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;4348 otherModel.numberTimesOptimal_ = numberTimesOptimal_;4349 otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;4350 otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;4351 }4352 /* Constructs a non linear cost from list of nonlinearities (columns only)4353 First lower of each column is taken as real lower4354 Last lower is taken as real upper and cost ignored4355 4356 Returns nonzero if bad data e.g. lowers not monotonic4357 */4358 int4359 ClpInterior::createPiecewiseLinearCosts(const int * starts,4360 const double * lower, const double * gradient)4361 {4362 delete nonLinearCost_;4363 // Set up feasible bounds and check monotonicity4364 int iColumn;4365 int returnCode=0;4366 4367 for (iColumn=0;iColumn<numberColumns_;iColumn++) {4368 int iIndex = starts[iColumn];4369 int end = starts[iColumn+1]1;4370 columnLower_[iColumn] = lower[iIndex];4371 columnUpper_[iColumn] = lower[end];4372 double value = columnLower_[iColumn];4373 iIndex++;4374 for (;iIndex<end;iIndex++) {4375 if (lower[iIndex]<value)4376 returnCode++; // not monotonic4377 value=lower[iIndex];4378 }4379 }4380 nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);4381 specialOptions_ = 2; // say keep4382 return returnCode;4383 }4384 /* For advanced use. When doing iterative solves things can get4385 nasty so on values pass if incoming solution has largest4386 infeasibility < incomingInfeasibility throw out variables4387 from basis until largest infeasibility < allowedInfeasibility4388 or incoming largest infeasibility.4389 If allowedInfeasibility>= incomingInfeasibility this is4390 always possible altough you may end up with an all slack basis.4391 4392 Defaults are 1.0,10.04393 */4394 void4395 ClpInterior::setValuesPassAction(float incomingInfeasibility,4396 float allowedInfeasibility)4397 {4398 incomingInfeasibility_=incomingInfeasibility;4399 allowedInfeasibility_=allowedInfeasibility;4400 assert(incomingInfeasibility_>=0.0);4401 assert(allowedInfeasibility_>=incomingInfeasibility_);4402 }4403 //#############################################################################4404 4405 ClpInteriorProgress::ClpInteriorProgress ()4406 {4407 int i;4408 for (i=0;i<CLP_PROGRESS;i++) {4409 objective_[i] = COIN_DBL_MAX;4410 infeasibility_[i] = 1.0; // set to an impossible value4411 numberInfeasibilities_[i]=1;4412 iterationNumber_[i]=1;4413 }4414 for (i=0;i<CLP_CYCLE;i++) {4415 in_[i]=1;4416 out_[i]=1;4417 way_[i]=0;4418 }4419 numberTimes_ = 0;4420 numberBadTimes_ = 0;4421 model_ = NULL;4422 }4423 4424 4425 //4426 4427 ClpInteriorProgress::~ClpInteriorProgress ()4428 {4429 }4430 // Copy constructor.4431 ClpInteriorProgress::ClpInteriorProgress(const ClpInteriorProgress &rhs)4432 {4433 int i;4434 for (i=0;i<CLP_PROGRESS;i++) {4435 objective_[i] = rhs.objective_[i];4436 infeasibility_[i] = rhs.infeasibility_[i];4437 numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];4438 iterationNumber_[i]=rhs.iterationNumber_[i];4439 }4440 for (i=0;i<CLP_CYCLE;i++) {4441 in_[i]=rhs.in_[i];4442 out_[i]=rhs.out_[i];4443 way_[i]=rhs.way_[i];4444 }4445 numberTimes_ = rhs.numberTimes_;4446 numberBadTimes_ = rhs.numberBadTimes_;4447 model_ = rhs.model_;4448 }4449 // Copy constructor.from model4450 ClpInteriorProgress::ClpInteriorProgress(ClpInterior * model)4451 {4452 model_ = model;4453 int i;4454 for (i=0;i<CLP_PROGRESS;i++) {4455 if (model_>algorithm()>=0)4456 objective_[i] = COIN_DBL_MAX;4457 else4458 objective_[i] = COIN_DBL_MAX;4459 infeasibility_[i] = 1.0; // set to an impossible value4460 numberInfeasibilities_[i]=1;4461 iterationNumber_[i]=1;4462 }4463 for (i=0;i<CLP_CYCLE;i++) {4464 in_[i]=1;4465 out_[i]=1;4466 way_[i]=0;4467 }4468 numberTimes_ = 0;4469 numberBadTimes_ = 0;4470 }4471 // Assignment operator. This copies the data4472 ClpInteriorProgress &4473 ClpInteriorProgress::operator=(const ClpInteriorProgress & rhs)4474 {4475 if (this != &rhs) {4476 int i;4477 for (i=0;i<CLP_PROGRESS;i++) {4478 objective_[i] = rhs.objective_[i];4479 infeasibility_[i] = rhs.infeasibility_[i];4480 numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];4481 iterationNumber_[i]=rhs.iterationNumber_[i];4482 }4483 for (i=0;i<CLP_CYCLE;i++) {4484 in_[i]=rhs.in_[i];4485 out_[i]=rhs.out_[i];4486 way_[i]=rhs.way_[i];4487 }4488 numberTimes_ = rhs.numberTimes_;4489 numberBadTimes_ = rhs.numberBadTimes_;4490 model_ = rhs.model_;4491 }4492 return *this;4493 }4494 // Seems to be something odd about exact comparison of doubles on linux4495 static bool equalDouble(double value1, double value2)4496 {4497 unsigned int *i1 = (unsigned int *) &value1;4498 unsigned int *i2 = (unsigned int *) &value2;4499 if (sizeof(unsigned int)*2==sizeof(double))4500 return (i1[0]==i2[0]&&i1[1]==i2[1]);4501 else4502 return (i1[0]==i2[0]);4503 }4504 int4505 ClpInteriorProgress::looping()4506 {4507 if (!model_)4508 return 1;4509 double objective = model_>rawObjectiveValue();4510 double infeasibility;4511 int numberInfeasibilities;4512 int iterationNumber = model_>numberIterations();4513 if (model_>algorithm()<0) {4514 // dual4515 infeasibility = model_>sumPrimalInfeasibilities();4516 numberInfeasibilities = model_>numberPrimalInfeasibilities();4517 } else {4518 //primal4519 infeasibility = model_>sumDualInfeasibilities();4520 numberInfeasibilities = model_>numberDualInfeasibilities();4521 }4522 int i;4523 int numberMatched=0;4524 int matched=0;4525 int nsame=0;4526 for (i=0;i<CLP_PROGRESS;i++) {4527 bool matchedOnObjective = equalDouble(objective,objective_[i]);4528 bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);4529 bool matchedOnInfeasibilities =4530 (numberInfeasibilities==numberInfeasibilities_[i]);4531 4532 if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {4533 matched = (1<<i);4534 // Check not same iteration4535 if (iterationNumber!=iterationNumber_[i]) {4536 numberMatched++;4537 // here mainly to get over compiler bug?4538 if (model_>messageHandler()>logLevel()>10)4539 printf("%d %d %d %d %d loop check\n",i,numberMatched,4540 matchedOnObjective, matchedOnInfeasibility,4541 matchedOnInfeasibilities);4542 } else {4543 // stuck but code should notice4544 nsame++;4545 }4546 }4547 if (i) {4548 objective_[i1] = objective_[i];4549 infeasibility_[i1] = infeasibility_[i];4550 numberInfeasibilities_[i1]=numberInfeasibilities_[i];4551 iterationNumber_[i1]=iterationNumber_[i];4552 }4553 }4554 objective_[CLP_PROGRESS1] = objective;4555 infeasibility_[CLP_PROGRESS1] = infeasibility;4556 numberInfeasibilities_[CLP_PROGRESS1]=numberInfeasibilities;4557 iterationNumber_[CLP_PROGRESS1]=iterationNumber;4558 if (nsame==CLP_PROGRESS)4559 numberMatched=CLP_PROGRESS; // really stuck4560 if (model_>progressFlag())4561 numberMatched=0;4562 numberTimes_++;4563 if (numberTimes_<10)4564 numberMatched=0;4565 // skip if just last time as may be checking something4566 if (matched==(1<<(CLP_PROGRESS1)))4567 numberMatched=0;4568 if (numberMatched) {4569 model_>messageHandler()>message(CLP_POSSIBLELOOP,model_>messages())4570 <<numberMatched4571 <<matched4572 <<numberTimes_4573 <<CoinMessageEol;4574 numberBadTimes_++;4575 if (numberBadTimes_<10) {4576 // make factorize every iteration4577 model_>forceFactorization(1);4578 if (model_>algorithm()<0) {4579 // dual  change tolerance4580 model_>setCurrentDualTolerance(model_>currentDualTolerance()*1.05);4581 // if infeasible increase dual bound4582 if (model_>dualBound()<1.0e19) {4583 model_>setDualBound(model_>dualBound()*1.1);4584 }4585 } else {4586 // primal  change tolerance model_>setCurrentPrimalTolerance(model_>currentPrimalTolerance()*1.05);4587 // if infeasible increase infeasibility cost4588 if (model_>nonLinearCost()>numberInfeasibilities()&&4589 model_>infeasibilityCost()<1.0e19) {4590 model_>setInfeasibilityCost(model_>infeasibilityCost()*1.1);4591 }4592 }4593 return 2;4594 } else {4595 // look at solution and maybe declare victory4596 if (infeasibility<1.0e4) {4597 return 0;4598 } else {4599 model_>messageHandler()>message(CLP_LOOP,model_>messages())4600 <<CoinMessageEol;4601 #ifndef NDEBUG4602 abort();4603 #endif4604 return 3;4605 }4606 }4607 }4608 return 1;4609 }4610 // Returns previous objective (if 1)  current if (0)4611 double4612 ClpInteriorProgress::lastObjective(int back) const4613 {4614 return objective_[CLP_PROGRESS1back];4615 }4616 // Modify objective e.g. if dual infeasible in dual4617 void4618 ClpInteriorProgress::modifyObjective(double value)4619 {4620 objective_[CLP_PROGRESS1]=value;4621 }4622 // Returns previous iteration number (if 1)  current if (0)4623 int4624 ClpInteriorProgress::lastIterationNumber(int back) const4625 {4626 return iterationNumber_[CLP_PROGRESS1back];4627 }4628 // Start check at beginning of whileIterating4629 void4630 ClpInteriorProgress::startCheck()4631 {4632 int i;4633 for (i=0;i<CLP_CYCLE;i++) {4634 in_[i]=1;4635 out_[i]=1;4636 way_[i]=0;4637 }4638 }4639 // Returns cycle length in whileIterating4640 int4641 ClpInteriorProgress::cycle(int in, int out,int wayIn,int wayOut)4642 {4643 int i;4644 int matched=0;4645 return 0;4646 // first see if in matches any out4647 for (i=1;i<CLP_CYCLE;i++) {4648 if (in==out_[i]) {4649 // even if flip then suspicious4650 matched=1;4651 break;4652 }4653 }4654 if (!matched) {4655 // can't be cycle4656 for (i=0;i<CLP_CYCLE1;i++) {4657 in_[i]=in_[i+1];4658 out_[i]=out_[i+1];4659 way_[i]=way_[i+1];4660 }4661 } else {4662 // possible cycle4663 matched=0;4664 for (i=0;i<CLP_CYCLE1;i++) {4665 int k;4666 char wayThis = way_[i];4667 int inThis = in_[i];4668 int outThis = out_[i];4669 for(k=i+1;k<CLP_CYCLE;k++) {4670 if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {4671 matched=ki;4672 }4673 }4674 in_[i]=in_[i+1];4675 out_[i]=out_[i+1];4676 way_[i]=way_[i+1];4677 }4678 }4679 char way = 1wayIn+4*(1wayOut);4680 in_[CLP_CYCLE1]=in;4681 out_[CLP_CYCLE1]=out;4682 way_[CLP_CYCLE1]=way;4683 return matched;4684 }4685 // Allow initial dense factorization4686 void4687 ClpInterior::setInitialDenseFactorization(bool onOff)4688 {4689 if (onOff)4690 specialOptions_ = 8;4691 else4692 specialOptions_ &= ~8;4693 }4694 bool4695 ClpInterior::initialDenseFactorization() const4696 {4697 return (specialOptions_&8)!=0;4698 }
Note: See TracChangeset
for help on using the changeset viewer.