Changeset 1825
- Timestamp:
- Nov 20, 2011 11:02:57 AM (9 years ago)
- Location:
- trunk/Clp/src
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Clp/src/ClpEventHandler.cpp
r1665 r1825 63 63 return 0; // say normal exit 64 64 } 65 /* This can do whatever it likes. Return code -1 means no action. 66 This passes in something 67 */ 68 int 69 ClpEventHandler::eventWithInfo(Event whichEvent, void * info) 70 { 71 return -1; 72 } 65 73 /* set model. */ 66 74 void -
trunk/Clp/src/ClpEventHandler.hpp
r1790 r1825 65 65 endOfCreateRim, 66 66 slightlyInfeasible, 67 modifyMatrixInMiniPresolve, 68 moreMiniPresolve, 69 modifyMatrixInMiniPostsolve, 67 70 noTheta // At end (because no pivot) 68 71 }; … … 78 81 */ 79 82 virtual int event(Event whichEvent); 83 /** This can do whatever it likes. Return code -1 means no action. 84 This passes in something 85 */ 86 virtual int eventWithInfo(Event whichEvent, void * info) ; 80 87 //@} 81 88 -
trunk/Clp/src/ClpModel.cpp
r1792 r1825 1583 1583 setRowScale(NULL); 1584 1584 setColumnScale(NULL); 1585 } 1586 // Deletes rows AND columns (does not reallocate) 1587 void 1588 ClpModel::deleteRowsAndColumns(int numberRows, const int * whichRows, 1589 int numberColumns, const int * whichColumns) 1590 { 1591 if (!numberColumns) { 1592 deleteRows(numberRows,whichRows); 1593 } else if (!numberRows) { 1594 deleteColumns(numberColumns,whichColumns); 1595 } else { 1596 whatsChanged_ &= ~511; // all changed 1597 bool doStatus = status_!=NULL; 1598 int numberTotal=numberRows_+numberColumns_; 1599 int * backRows = new int [numberTotal]; 1600 int * backColumns = backRows+numberRows_; 1601 memset(backRows,0,numberTotal*sizeof(int)); 1602 int newNumberColumns=0; 1603 for (int i=0;i<numberColumns;i++) { 1604 int iColumn=whichColumns[i]; 1605 if (iColumn>=0&&iColumn<numberColumns_) 1606 backColumns[iColumn]=-1; 1607 } 1608 assert (objective_->type()==1); 1609 double * obj = objective(); 1610 for (int i=0;i<numberColumns_;i++) { 1611 if (!backColumns[i]) { 1612 columnActivity_[newNumberColumns] = columnActivity_[i]; 1613 reducedCost_[newNumberColumns] = reducedCost_[i]; 1614 obj[newNumberColumns] = obj[i]; 1615 columnLower_[newNumberColumns] = columnLower_[i]; 1616 columnUpper_[newNumberColumns] = columnUpper_[i]; 1617 if (doStatus) 1618 status_[newNumberColumns] = status_[i]; 1619 backColumns[i]=newNumberColumns++; 1620 } 1621 } 1622 integerType_ = deleteChar(integerType_, numberColumns_, 1623 numberColumns, whichColumns, newNumberColumns, true); 1624 #ifndef CLP_NO_STD 1625 // Now works if which out of order 1626 if (lengthNames_) { 1627 for (int i=0;i<numberColumns_;i++) { 1628 int iColumn=backColumns[i]; 1629 if (iColumn) 1630 columnNames_[iColumn] = columnNames_[i]; 1631 } 1632 columnNames_.erase(columnNames_.begin() + newNumberColumns, columnNames_.end()); 1633 } 1634 #endif 1635 int newNumberRows=0; 1636 assert (!rowObjective_); 1637 unsigned char * status2 = status_ + numberColumns_; 1638 for (int i=0;i<numberRows;i++) { 1639 int iRow=whichRows[i]; 1640 if (iRow>=0&&iRow<numberRows_) 1641 backRows[iRow]=-1; 1642 } 1643 for (int i=0;i<numberRows_;i++) { 1644 if (!backRows[i]) { 1645 rowActivity_[newNumberRows] = rowActivity_[i]; 1646 dual_[newNumberRows] = dual_[i]; 1647 rowLower_[newNumberRows] = rowLower_[i]; 1648 rowUpper_[newNumberRows] = rowUpper_[i]; 1649 if (doStatus) 1650 status2[newNumberRows] = status2[i]; 1651 backRows[i]=newNumberRows++; 1652 } 1653 } 1654 #ifndef CLP_NO_STD 1655 // Now works if which out of order 1656 if (lengthNames_) { 1657 for (int i=0;i<numberRows_;i++) { 1658 int iRow=backRows[i]; 1659 if (iRow) 1660 rowNames_[iRow] = rowNames_[i]; 1661 } 1662 rowNames_.erase(rowNames_.begin() + newNumberRows, rowNames_.end()); 1663 } 1664 #endif 1665 // possible matrix is not full 1666 ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *>(matrix_); 1667 CoinPackedMatrix * matrix = clpMatrix ? clpMatrix->matrix() : NULL; 1668 if (matrix_->getNumCols() < numberColumns_) { 1669 assert (matrix); 1670 CoinBigIndex nel=matrix->getNumElements(); 1671 int n=matrix->getNumCols(); 1672 matrix->reserve(numberColumns_,nel); 1673 CoinBigIndex * columnStart = matrix->getMutableVectorStarts(); 1674 int * columnLength = matrix->getMutableVectorLengths(); 1675 for (int i=n;i<numberColumns_;i++) { 1676 columnStart[i]=nel; 1677 columnLength[i]=0; 1678 } 1679 } 1680 if (matrix) { 1681 matrix->setExtraMajor(0.1); 1682 //CoinPackedMatrix temp(*matrix); 1683 matrix->setExtraGap(0.0); 1684 matrix->setExtraMajor(0.0); 1685 int * row = matrix->getMutableIndices(); 1686 CoinBigIndex * columnStart = matrix->getMutableVectorStarts(); 1687 int * columnLength = matrix->getMutableVectorLengths(); 1688 double * element = matrix->getMutableElements(); 1689 newNumberColumns=0; 1690 CoinBigIndex n=0; 1691 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 1692 if (backColumns[iColumn]>=0) { 1693 CoinBigIndex start = columnStart[iColumn]; 1694 int nSave=n; 1695 columnStart[newNumberColumns]=n; 1696 for (CoinBigIndex j=start;j<start+columnLength[iColumn];j++) { 1697 int iRow=row[j]; 1698 iRow = backRows[iRow]; 1699 if (iRow>=0) { 1700 row[n]=iRow; 1701 element[n++]=element[j]; 1702 } 1703 } 1704 columnLength[newNumberColumns++]=n-nSave; 1705 } 1706 } 1707 columnStart[newNumberColumns]=n; 1708 matrix->setNumElements(n); 1709 matrix->setMajorDim(newNumberColumns); 1710 matrix->setMinorDim(newNumberRows); 1711 clpMatrix->setNumberActiveColumns(newNumberColumns); 1712 //temp.deleteRows(numberRows, whichRows); 1713 //temp.deleteCols(numberColumns, whichColumns); 1714 //assert(matrix->isEquivalent2(temp)); 1715 //*matrix=temp; 1716 } else { 1717 matrix_->deleteRows(numberRows, whichRows); 1718 matrix_->deleteCols(numberColumns, whichColumns); 1719 } 1720 numberColumns_ = newNumberColumns; 1721 numberRows_ = newNumberRows; 1722 delete [] backRows; 1723 // set state back to unknown 1724 problemStatus_ = -1; 1725 secondaryStatus_ = 0; 1726 delete [] ray_; 1727 ray_ = NULL; 1728 if (savedRowScale_ != rowScale_) { 1729 delete [] rowScale_; 1730 delete [] columnScale_; 1731 } 1732 rowScale_ = NULL; 1733 columnScale_ = NULL; 1734 delete scaledMatrix_; 1735 scaledMatrix_ = NULL; 1736 delete rowCopy_; 1737 rowCopy_ = NULL; 1738 } 1585 1739 } 1586 1740 // Add one row -
trunk/Clp/src/ClpModel.hpp
r1780 r1825 184 184 /// Deletes columns 185 185 void deleteColumns(int number, const int * which); 186 /// Deletes rows AND columns (keeps old sizes) 187 void deleteRowsAndColumns(int numberRows, const int * whichRows, 188 int numberColumns, const int * whichColumns); 186 189 /// Add one column 187 190 void addColumn(int numberInColumn, -
trunk/Clp/src/ClpObjective.hpp
r1665 r1825 96 96 //@{ 97 97 /// Returns type (above 63 is extra information) 98 inline int type() {98 inline int type() const { 99 99 return type_; 100 } 101 /// Sets type (above 63 is extra information) 102 inline void setType(int value) { 103 type_ = value; 100 104 } 101 105 /// Whether activated -
trunk/Clp/src/ClpPackedMatrix.hpp
r1722 r1825 334 334 flags_ = (matrix_->hasGaps()) ? (flags_ | 2) : (flags_ & (~2)); 335 335 } 336 /// number of active columns (normally same as number of columns) 337 inline int numberActiveColumns() const 338 { return numberActiveColumns_;} 339 /// Set number of active columns (normally same as number of columns) 340 inline void setNumberActiveColumns(int value) 341 { numberActiveColumns_ = value;} 336 342 //@} 337 343 -
trunk/Clp/src/ClpSimplex.hpp
r1780 r1825 373 373 but only if > threshold */ 374 374 void removeSuperBasicSlacks(int threshold=0); 375 /** Mini presolve (faster) 376 Char arrays must be numberRows and numberColumns long 377 on entry second part must be filled in as follows - 378 0 - possible 379 >0 - take out and do something (depending on value - TBD) 380 -1 row/column can't vanish but can have entries removed/changed 381 -2 don't touch at all 382 on exit <=0 ones will be in presolved problem 383 struct will be created and will be long enough 384 (information on length etc in first entry) 385 user must delete struct 386 */ 387 ClpSimplex * miniPresolve(char * rowType, char * columnType,void ** info); 388 /// After mini presolve 389 void miniPostsolve(const ClpSimplex * presolvedModel,void * info); 375 390 /** Write the basis in MPS format to the specified file. 376 391 If writeValues true writes values of structurals -
trunk/Clp/src/ClpSimplexOther.cpp
r1817 r1825 28 28 #include <string> 29 29 #include <stdio.h> 30 #include <iostream> 30 #include <iostream> 31 #ifdef HAS_CILK 32 #include <cilk/cilk.h> 33 #else 34 #define cilk_for for 35 #define cilk_spawn 36 #define cilk_sync 37 #endif 38 #ifdef INT_IS_8 39 #define COIN_ANY_BITS_PER_INT 64 40 #define COIN_ANY_SHIFT_PER_INT 6 41 #define COIN_ANY_MASK_PER_INT 0x3f 42 #else 43 #define COIN_ANY_BITS_PER_INT 32 44 #define COIN_ANY_SHIFT_PER_INT 5 45 #define COIN_ANY_MASK_PER_INT 0x1f 46 #endif 31 47 /* Dual ranging. 32 48 This computes increase/decrease in cost for each given variable and corresponding … … 95 111 int iRow = backPivot[iSequence]; 96 112 assert (iRow >= 0); 113 #ifndef COIN_FAC_NEW 97 114 double plusOne = 1.0; 98 115 rowArray_[0]->createPacked(1, &iRow, &plusOne); 116 #else 117 rowArray_[0]->createOneUnpackedElement( iRow, 1.0); 118 #endif 99 119 factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]); 100 120 // put row of tableau in rowArray[0] and columnArray[0] 101 121 matrix_->transposeTimes(this, -1.0, 102 122 rowArray_[0], columnArray_[1], columnArray_[0]); 123 #if COIN_FAC_NEW 124 assert (!rowArray_[0]->packedMode()); 125 #endif 103 126 double alphaIncrease; 104 127 double alphaDecrease; … … 115 138 } else { 116 139 int number = rowArray_[0]->getNumElements(); 140 #if COIN_FAC_NEW 141 const int * index = rowArray_[0]->getIndices(); 142 #endif 117 143 double scale2 = 0.0; 118 144 int j; 119 145 for (j = 0; j < number; j++) { 146 #ifndef COIN_FAC_NEW 120 147 scale2 += arrayX[j] * arrayX[j]; 148 #else 149 int iRow=index[j]; 150 scale2 += arrayX[iRow] * arrayX[iRow]; 151 #endif 121 152 } 122 153 scale2 = 1.0 / sqrt(scale2); … … 271 302 int iSequence = which[i]; 272 303 int iSequence2 = iSequence + addSequence; 304 #ifndef COIN_FAC_NEW 273 305 double alpha = work[i]; 306 #else 307 double alpha = !addSequence ? work[i] : work[iSequence]; 308 #endif 274 309 if (fabs(alpha) < acceptablePivot) 275 310 continue; … … 381 416 // Non trivial 382 417 // Other bound is ignored 418 #ifndef COIN_FAC_NEW 383 419 unpackPacked(rowArray_[1], iSequence); 420 #else 421 unpack(rowArray_[1], iSequence); 422 #endif 384 423 factorization_->updateColumn(rowArray_[2], rowArray_[1]); 385 424 // Get extra rows … … 451 490 { 452 491 // Other bound is ignored 492 #ifndef COIN_FAC_NEW 453 493 unpackPacked(rowArray_[1], iSequence); 494 #else 495 unpack(rowArray_[1], iSequence); 496 #endif 454 497 factorization_->updateColumn(rowArray_[2], rowArray_[1]); 455 498 // Get extra rows … … 467 510 468 511 int iRow = which[iIndex]; 512 #ifndef COIN_FAC_NEW 469 513 double alpha = work[iIndex] * way; 514 #else 515 double alpha = work[iRow] * way; 516 #endif 470 517 int iPivot = pivotVariable_[iRow]; 471 518 if (iPivot == whichOther) { … … 543 590 544 591 int iRow = which[iIndex]; 592 #ifndef COIN_FAC_NEW 545 593 double alpha = work[iIndex] * way; 594 #else 595 double alpha = work[iRow] * way; 596 #endif 546 597 int iPivot = pivotVariable_[iRow]; 547 598 double oldValue = solution_[iPivot]; … … 2829 2880 // To mark as odd 2830 2881 char * markDone = reinterpret_cast<char *>(lowerList+numberTotal); 2831 memset(markDone,0,numberTotal);2882 //memset(markDone,0,numberTotal); 2832 2883 int * backwardBasic = lowerList+2*numberTotal; 2833 2884 parametricsData paramData; … … 2968 3019 #endif 2969 3020 if (!returnCode) { 3021 assert (objective_->type()==1); 3022 objective_->setType(2); // in case matrix empty 2970 3023 returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0); 3024 objective_->setType(1); 2971 3025 if (!returnCode) { 2972 3026 double saveDualBound=dualBound_; … … 3050 3104 dualRowPivot_->setModel(this); 3051 3105 #endif 3052 for (int i=0;i<numberRows_+numberColumns_;i++)3053 3106 //for (int i=0;i<numberRows_+numberColumns_;i++) 3107 //setFakeBound(i, noFake); 3054 3108 // Now do parametrics 3055 3109 handler_->message(CLP_PARAMETRICS_STATS, messages_) … … 3085 3139 dualBound_ = saveDualBound; 3086 3140 //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data); 3141 } else { 3142 // check if empty 3143 //if (!numberRows_||!matrix_->getNumElements()) { 3144 // success 3145 #ifdef CLP_USER_DRIVEN 3146 //theta_ = endingTheta; 3147 //eventHandler_->event(ClpEventHandler::theta); 3148 #endif 3149 //} 3087 3150 } 3088 3151 if (problemStatus_==2) { … … 3238 3301 double * saveDuals = NULL; 3239 3302 problemStatus_=-1; 3303 ClpObjective * saveObjective = objective_; 3240 3304 reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data); 3305 if (saveObjective!=objective_) { 3306 delete objective_; 3307 objective_=saveObjective; 3308 } 3241 3309 int pass=100; 3242 3310 double moved=0.0; 3243 3311 while(sumPrimalInfeasibilities_) { 3312 //printf("INFEAS pass %d %d %g\n",100-pass,numberPrimalInfeasibilities_, 3313 // sumPrimalInfeasibilities_); 3244 3314 pass--; 3245 3315 if (!pass) … … 3528 3598 for (int i=0;i<n;i++) { 3529 3599 int iSequence = lowerList[i]; 3530 lower_[iSequence] += change * lowerChange[iSequence]; 3600 double newValue = lower_[iSequence] + change * lowerChange[iSequence]; 3601 lower_[iSequence] = newValue; 3531 3602 if(getStatus(iSequence)==atLowerBound) { 3532 solution_[iSequence] = lower_[iSequence];3603 solution_[iSequence] = newValue; 3533 3604 } 3605 #if 0 3606 // may have to adjust other bound 3607 double otherValue = upper_[iSequence]; 3608 if (otherValue-newValue<dualBound_) { 3609 //originalBound(iSequence,useTheta,lowerChange,upperChange); 3610 //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence); 3611 //ClpTraceDebug (fabs(lower_[iSequence]-newValue)<1.0e-5); 3612 } 3613 #endif 3534 3614 } 3535 3615 n=upperList[-1]; 3536 3616 for (int i=0;i<n;i++) { 3537 3617 int iSequence = upperList[i]; 3538 upper_[iSequence] += change * upperChange[iSequence]; 3618 double newValue = upper_[iSequence] + change * upperChange[iSequence]; 3619 upper_[iSequence] = newValue; 3539 3620 if(getStatus(iSequence)==atUpperBound|| 3540 3621 getStatus(iSequence)==isFixed) { 3541 solution_[iSequence] = upper_[iSequence]; 3622 solution_[iSequence] = newValue; 3623 } 3624 // may have to adjust other bound 3625 double otherValue = lower_[iSequence]; 3626 if (newValue-otherValue<dualBound_) { 3627 //originalBound(iSequence,useTheta,lowerChange,upperChange); 3628 //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence); 3629 //ClpTraceDebug (fabs(upper_[iSequence]-newValue)<1.0e-5); 3542 3630 } 3543 3631 } … … 3593 3681 // update the incoming column 3594 3682 double btranAlpha = -alpha_ * directionOut_; // for check 3683 #ifndef COIN_FAC_NEW 3595 3684 unpackPacked(rowArray_[1]); 3685 #else 3686 unpack(rowArray_[1]); 3687 #endif 3596 3688 // and update dual weights (can do in parallel - with extra array) 3597 3689 rowArray_[2]->clear(); … … 3773 3865 int * which = rowArray_[1]->getIndices(); 3774 3866 assert (!rowArray_[4]->packedMode()); 3867 #ifndef COIN_FAC_NEW 3775 3868 assert (rowArray_[1]->packedMode()); 3869 #else 3870 assert (!rowArray_[1]->packedMode()); 3871 #endif 3776 3872 double pivotValue = rowArray_[4]->denseVector()[pivotRow_]; 3777 3873 double multiplier = -pivotValue/alpha_; … … 3779 3875 for (int i = 0; i < number; i++) { 3780 3876 int iRow = which[i]; 3877 #ifndef COIN_FAC_NEW 3781 3878 rowArray_[4]->quickAddNonZero(iRow,multiplier*work[i]); 3879 #else 3880 rowArray_[4]->quickAddNonZero(iRow,multiplier*work[iRow]); 3881 #endif 3782 3882 } 3783 3883 } … … 3855 3955 solution_[sequenceOut_] = valueOut_; 3856 3956 int whatNext = housekeeping(objectiveChange); 3957 reinterpret_cast<ClpSimplexDual *>(this)->originalBound(sequenceIn_); 3857 3958 assert (backwardBasic[sequenceOut_]==pivotRow_); 3858 3959 backwardBasic[sequenceOut_]=-1; … … 4147 4248 // create as packed 4148 4249 double direction = directionOut_; 4250 #ifndef COIN_FAC_NEW 4149 4251 rowArray_[0]->createPacked(1, &pivotRow_, &direction); 4252 #else 4253 rowArray_[0]->createOneUnpackedElement(pivotRow_, direction); 4254 #endif 4150 4255 factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]); 4151 4256 // put row of tableau in rowArray[0] and columnArray[0] … … 4165 4270 // Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none) 4166 4271 int 4167 ClpSimplexOther::nextTheta(int type, double maxTheta, parametricsData & paramData,4272 ClpSimplexOther::nextTheta(int /*type*/, double maxTheta, parametricsData & paramData, 4168 4273 const double * /*changeObjective*/) 4169 4274 { … … 4173 4278 const int * upperList = paramData.upperList; 4174 4279 int iSequence; 4175 theta_ = maxTheta;4176 4280 bool toLower = false; 4177 assert (type==1);4281 //assert (type==1); 4178 4282 // may need to decide based on model? 4179 4283 bool needFullUpdate = rowArray_[4]->getNumElements()==0; … … 4375 4479 } 4376 4480 } 4377 pivotRow_ = -1;4378 4481 const int * index = rowArray_[4]->getIndices(); 4379 4482 int number = rowArray_[4]->getNumElements(); 4380 char * markDone = paramData.markDone; 4483 int * markDone = reinterpret_cast<int *>(paramData.markDone); 4484 int nToZero=(numberRows_+numberColumns_+COIN_ANY_BITS_PER_INT-1)>>COIN_ANY_SHIFT_PER_INT; 4485 memset(markDone,0,nToZero*sizeof(int)); 4381 4486 const int * backwardBasic = paramData.backwardBasic; 4382 4487 // first ones with alpha 4383 for (int i=0;i<number;i++) { 4488 double theta1=maxTheta; 4489 double theta2=maxTheta; 4490 //bool toLower2=true; 4491 int pivotRow1=-1; 4492 int pivotRow2=-1; 4493 cilk_for (int i=0;i<number;i++) { 4384 4494 int iPivot=index[i]; 4385 4495 iSequence = pivotVariable_[iPivot]; 4386 assert(!markDone[iSequence]); 4387 markDone[iSequence]=1; 4496 //assert(!markDone[iSequence]); 4497 int word = iSequence >> COIN_ANY_SHIFT_PER_INT; 4498 int bit = iSequence & COIN_ANY_MASK_PER_INT; 4499 markDone[word] |= ( 1 << bit ); 4500 //markDone[iSequence]=1; 4388 4501 // solution value will be sol - theta*alpha 4389 4502 // bounds will be bounds + change *theta … … 4392 4505 double thetaCoefficientLower = lowerChange[iSequence] + alpha; 4393 4506 double thetaCoefficientUpper = upperChange[iSequence] + alpha; 4507 //#define NO_TEST 4508 #ifndef NO_TEST 4394 4509 if (thetaCoefficientLower > 1.0e-8) { 4510 #endif 4395 4511 double currentLower = lower_[iSequence]; 4396 4512 ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_); 4397 4513 double gap=currentSolution-currentLower; 4398 if (thetaCoefficientLower*theta_>gap) { 4399 theta_ = gap/thetaCoefficientLower; 4400 toLower=true; 4401 pivotRow_=iPivot; 4402 } 4514 if (thetaCoefficientLower*theta1>gap) { 4515 theta1 = gap/thetaCoefficientLower; 4516 //toLower=true; 4517 pivotRow1=iPivot; 4518 } 4519 #ifndef NO_TEST 4403 4520 } 4404 4521 if (thetaCoefficientUpper < -1.0e-8) { 4522 #endif 4405 4523 double currentUpper = upper_[iSequence]; 4406 4524 ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_); 4407 double gap=currentSolution-currentUpper; //negative 4408 if (thetaCoefficientUpper*theta_<gap) { 4409 theta_ = gap/thetaCoefficientUpper; 4410 toLower=false; 4411 pivotRow_=iPivot; 4412 } 4413 } 4525 double gap2=currentSolution-currentUpper; //negative 4526 if (thetaCoefficientUpper*theta2<gap2) { 4527 theta2 = gap2/thetaCoefficientUpper; 4528 //toLower=false; 4529 pivotRow2=iPivot; 4530 } 4531 #ifndef NO_TEST 4532 } 4533 #endif 4414 4534 } 4415 4535 // now others … … 4417 4537 for (int i=0;i<nLook;i++) { 4418 4538 int iSequence = lowerList[i]; 4419 if (getColumnStatus(iSequence)==basic&&!markDone[iSequence]) { 4539 int word = iSequence >> COIN_ANY_SHIFT_PER_INT; 4540 int bit = iSequence & COIN_ANY_MASK_PER_INT; 4541 if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) { 4420 4542 double currentSolution = solution_[iSequence]; 4421 4543 double currentLower = lower_[iSequence]; … … 4424 4546 if (thetaCoefficient > 0.0) { 4425 4547 double gap=currentSolution-currentLower; 4426 if (thetaCoefficient*theta _>gap) {4427 theta _= gap/thetaCoefficient;4428 toLower=true;4429 pivotRow _= backwardBasic[iSequence];4548 if (thetaCoefficient*theta1>gap) { 4549 theta1 = gap/thetaCoefficient; 4550 //toLower=true; 4551 pivotRow1 = backwardBasic[iSequence]; 4430 4552 } 4431 4553 } … … 4435 4557 for (int i=0;i<nLook;i++) { 4436 4558 int iSequence = upperList[i]; 4437 if (getColumnStatus(iSequence)==basic&&!markDone[iSequence]) { 4559 int word = iSequence >> COIN_ANY_SHIFT_PER_INT; 4560 int bit = iSequence & COIN_ANY_MASK_PER_INT; 4561 if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) { 4438 4562 double currentSolution = solution_[iSequence]; 4439 4563 double currentUpper = upper_[iSequence]; … … 4442 4566 if (thetaCoefficient < 0) { 4443 4567 double gap=currentSolution-currentUpper; //negative 4444 if (thetaCoefficient*theta_<gap) { 4445 theta_ = gap/thetaCoefficient; 4446 toLower=false; 4447 pivotRow_ = backwardBasic[iSequence]; 4448 } 4449 } 4450 } 4568 if (thetaCoefficient*theta2<gap) { 4569 theta2 = gap/thetaCoefficient; 4570 //toLower=false; 4571 pivotRow2 = backwardBasic[iSequence]; 4572 } 4573 } 4574 } 4575 } 4576 if (theta2<theta1) { 4577 theta_=theta2; 4578 toLower=false; 4579 pivotRow_=pivotRow2; 4580 } else { 4581 theta_=theta1; 4582 toLower=true; 4583 pivotRow_=pivotRow1; 4451 4584 } 4452 4585 theta_ = CoinMax(theta_,0.0); … … 4456 4589 int iPivot = index[iRow]; 4457 4590 iSequence = pivotVariable_[iPivot]; 4458 markDone[iSequence]=0;4591 //markDone[iSequence]=0; 4459 4592 // solution value will be sol - theta*alpha 4460 4593 double alpha = array[iPivot]; 4461 4594 solution_[iSequence] -= theta_ * alpha; 4462 4595 } 4596 #if 0 4463 4597 } else { 4464 4598 for (int iRow = 0; iRow < number; iRow++) { … … 4467 4601 markDone[iSequence]=0; 4468 4602 } 4603 #endif 4469 4604 } 4470 4605 #if 0 … … 6234 6369 // unpack column 6235 6370 assert(!vec[0]->getNumElements()); 6371 #ifndef COIN_FAC_NEW 6236 6372 unpackPacked(vec[0]); 6373 #else 6374 unpack(vec[0]); 6375 #endif 6237 6376 // update 6238 6377 assert(!vec[1]->getNumElements()); 6239 6378 factorization_->updateColumnFT(vec[1], vec[0]); 6379 const double * array = vec[0]->denseVector(); 6380 #ifndef COIN_FAC_NEW 6240 6381 // Find alpha 6241 6382 const int * indices = vec[0]->getIndices(); 6242 const double * array = vec[0]->denseVector();6243 6383 int n=vec[0]->getNumElements(); 6244 6384 alpha_=0.0; … … 6249 6389 } 6250 6390 } 6391 #else 6392 alpha_ = array[pivotRow_]; 6393 #endif 6251 6394 int updateStatus=2; 6252 6395 if (fabs(alpha_)>1.0e-7) … … 6406 6549 if (!(state&1)) { 6407 6550 // update the incoming column 6551 #ifndef COIN_FAC_NEW 6408 6552 unpackPacked(rowArray_[1]); 6553 #else 6554 unpack(rowArray_[1]); 6555 #endif 6409 6556 factorization_->updateColumnFT(rowArray_[2], rowArray_[1]); 6410 6557 } … … 6427 6574 assert (!rowArray_[0]->getNumElements()); 6428 6575 #endif 6576 #ifndef COIN_FAC_NEW 6429 6577 rowArray_[0]->createPacked(1, &pivotRow_, &direction); 6578 #else 6579 rowArray_[0]->createOneUnpackedElement(pivotRow_, direction); 6580 #endif 6430 6581 factorization_->updateColumnTranspose(rowArray_[2], rowArray_[0]); 6431 6582 rowArray_[3]->clear(); … … 6459 6610 for (int i = 0; i < number; i++) { 6460 6611 int iRow = which[i]; 6612 #ifndef COIN_FAC_NEW 6461 6613 double alpha = work[i]; 6614 #else 6615 double alpha = work[iRow]; 6616 #endif 6462 6617 int iPivot = pivotVariable_[iRow]; 6463 6618 dualIn_ -= alpha * cost_[iPivot]; … … 6482 6637 number = rowArray_[0]->getNumElements(); 6483 6638 element = rowArray_[0]->denseVector(); 6639 #ifndef COIN_FAC_NEW 6484 6640 assert (rowArray_[0]->packedMode()); 6485 6641 for (i = 0; i < number; i++) { … … 6489 6645 element[i] = 0.0; 6490 6646 } 6647 #else 6648 assert (!rowArray_[0]->packedMode()); 6649 for (i = 0; i < number; i++) { 6650 int iSequence = index[i]; 6651 dj_[iSequence+numberColumns_] += multiplier*element[iSequence]; 6652 dual_[iSequence] = dj_[iSequence+numberColumns_]; 6653 element[iSequence] = 0.0; 6654 } 6655 #endif 6491 6656 rowArray_[0]->setNumElements(0); 6492 6657 double oldCost = cost_[sequenceOut_]; … … 6565 6730 double btranAlpha = -alpha_ * directionOut_; // for check 6566 6731 rowArray_[1]->clear(); 6732 #ifndef COIN_FAC_NEW 6567 6733 unpackPacked(rowArray_[1]); 6734 #else 6735 unpack(rowArray_[1]); 6736 #endif 6568 6737 // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]); 6569 6738 // and update dual weights (can do in parallel - with extra array) … … 7202 7371 delete [] whichRows; 7203 7372 } 7373 /* 7374 1 (and 4) redundant (and 8 is user) 7375 2 sub 7376 11 movable column 7377 13 empty (or initially fixed) column 7378 14 doubleton 7379 */ 7380 typedef struct { 7381 double oldRowLower; 7382 double oldRowUpper; 7383 int row; 7384 int lengthRow; 7385 } clpPresolveInfo1_4_8; 7386 // can be used instead of 1_4_8 7387 typedef struct { 7388 double oldRowLower; 7389 double oldRowUpper; 7390 int row; 7391 int lengthRow; 7392 double * rowLowerX; 7393 double * rowUpperX; 7394 double * tempElement; 7395 int * tempIndex; 7396 int otherRow; 7397 } clpPresolveInfo8; 7398 typedef struct { 7399 double oldRowLower; 7400 double oldRowUpper; 7401 double oldColumnLower; 7402 double oldColumnUpper; 7403 double coefficient; 7404 // 2 is upper 7405 double oldRowLower2; 7406 double oldRowUpper2; 7407 double coefficient2; 7408 int row; 7409 int row2; 7410 int column; 7411 } clpPresolveInfo2; 7412 typedef struct { 7413 double oldColumnLower; 7414 double oldColumnUpper; 7415 double fixedTo; 7416 int column; 7417 int lengthColumn; 7418 } clpPresolveInfo11; 7419 typedef struct { 7420 double oldColumnLower; 7421 double oldColumnUpper; 7422 int column; 7423 } clpPresolveInfo13; 7424 typedef struct { 7425 double oldColumnLower; 7426 double oldColumnUpper; 7427 double oldColumnLower2; 7428 double oldColumnUpper2; 7429 double oldObjective2; 7430 double value1; 7431 double rhs; 7432 int type; 7433 int row; 7434 int column; 7435 int column2; 7436 int lengthColumn2; 7437 } clpPresolveInfo14; 7438 typedef struct { 7439 int infoOffset; 7440 int type; 7441 } clpPresolveInfo; 7442 typedef struct { 7443 int numberEntries; 7444 int maximumEntries; 7445 int numberInitial; 7446 clpPresolveInfo * start; 7447 } listInfo; 7448 typedef struct { 7449 char * putStuff; 7450 char * startStuff; 7451 CoinBigIndex maxStuff; 7452 } saveInfo; 7453 typedef struct { 7454 double * elements; 7455 int * indices; 7456 char * startStuff; 7457 } restoreInfo; 7458 // struct must match in handler 7459 typedef struct { 7460 ClpSimplex * model; 7461 CoinPackedMatrix * rowCopy; 7462 char * rowType; 7463 char * columnType; 7464 saveInfo * stuff; 7465 clpPresolveInfo * info; 7466 int * nActions; 7467 } clpPresolveMore; 7468 void ClpCopyToMiniSave(saveInfo & where, const char * info, unsigned int sizeInfo,int numberElements, 7469 const int * indices, const double * elements) 7470 { 7471 char * put = where.putStuff; 7472 int n = numberElements*(sizeof(int)+sizeof(double))+sizeInfo; 7473 if (n+(put-where.startStuff)>where.maxStuff) { 7474 where.maxStuff += CoinMax(where.maxStuff/2 + 10000, 2*n); 7475 char * temp = new char[where.maxStuff]; 7476 int k = put-where.startStuff; 7477 memcpy(temp,where.startStuff,k); 7478 delete [] where.startStuff; 7479 where.startStuff=temp; 7480 put = temp+k; 7481 } 7482 memcpy(put,info,sizeInfo); 7483 put += sizeInfo; 7484 memcpy(put,indices,numberElements*sizeof(int)); 7485 put += numberElements*sizeof(int); 7486 memcpy(put,elements,numberElements*sizeof(double)); 7487 put += numberElements*sizeof(double); 7488 where.putStuff=put; 7489 } 7490 static void copyFromSave(restoreInfo & where, clpPresolveInfo & info, void * thisInfoX) 7491 { 7492 char * get = where.startStuff+info.infoOffset; 7493 int type = info.type; 7494 int n=0; 7495 switch(type) { 7496 case 1: 7497 case 4: 7498 // redundant 7499 { 7500 clpPresolveInfo1_4_8 thisInfo; 7501 memcpy(&thisInfo,get,sizeof(clpPresolveInfo1_4_8)); 7502 memcpy(thisInfoX,get,sizeof(clpPresolveInfo1_4_8)); 7503 get += sizeof(clpPresolveInfo1_4_8); 7504 n = thisInfo.lengthRow; 7505 } 7506 break; 7507 case 8: 7508 case 9: 7509 // redundant 7510 { 7511 clpPresolveInfo8 thisInfo; 7512 memcpy(&thisInfo,get,sizeof(clpPresolveInfo8)); 7513 memcpy(thisInfoX,get,sizeof(clpPresolveInfo8)); 7514 get += sizeof(clpPresolveInfo8); 7515 n = thisInfo.lengthRow; 7516 } 7517 break; 7518 case 2: 7519 // sub 7520 { 7521 clpPresolveInfo2 thisInfo; 7522 memcpy(&thisInfo,get,sizeof(clpPresolveInfo2)); 7523 memcpy(thisInfoX,get,sizeof(clpPresolveInfo2)); 7524 get += sizeof(clpPresolveInfo2); 7525 } 7526 break; 7527 case 11: 7528 // movable column 7529 { 7530 clpPresolveInfo11 thisInfo; 7531 memcpy(&thisInfo,get,sizeof(clpPresolveInfo11)); 7532 memcpy(thisInfoX,get,sizeof(clpPresolveInfo11)); 7533 get += sizeof(clpPresolveInfo11); 7534 n = thisInfo.lengthColumn; 7535 } 7536 break; 7537 case 13: 7538 // empty (or initially fixed) column 7539 { 7540 clpPresolveInfo13 thisInfo; 7541 memcpy(&thisInfo,get,sizeof(clpPresolveInfo13)); 7542 memcpy(thisInfoX,get,sizeof(clpPresolveInfo13)); 7543 get += sizeof(clpPresolveInfo13); 7544 } 7545 break; 7546 case 14: 7547 // doubleton 7548 { 7549 clpPresolveInfo14 thisInfo; 7550 memcpy(&thisInfo,get,sizeof(clpPresolveInfo14)); 7551 memcpy(thisInfoX,get,sizeof(clpPresolveInfo14)); 7552 get += sizeof(clpPresolveInfo14); 7553 n = thisInfo.lengthColumn2; 7554 } 7555 break; 7556 } 7557 if (n) { 7558 memcpy(where.indices,get,n*sizeof(int)); 7559 get += n*sizeof(int); 7560 memcpy(where.elements,get,n*sizeof(double)); 7561 } 7562 } 7563 #define DEBUG_SOME 0 7564 // need more space 7565 static 7566 void moveAround(int numberColumns,CoinBigIndex numberElementsOriginal, 7567 int iColumn,int lengthNeeded, 7568 int * forward,int * backward, 7569 CoinBigIndex * columnStart,int * columnLength, 7570 int * row,double * element) 7571 { 7572 // we only get here if can't fit so if iColumn is last one need shuffle 7573 int last=backward[numberColumns]; 7574 bool needCompaction=false; 7575 CoinBigIndex lastElement=columnStart[numberColumns]; 7576 //assert(lastElement==2*(numberElementsOriginal+numberColumns)); 7577 // save length 7578 int length=columnLength[iColumn]; 7579 if (iColumn!=last) { 7580 CoinBigIndex put=columnStart[last]+columnLength[last]+3; 7581 if (put+lengthNeeded<=lastElement) { 7582 // copy 7583 CoinBigIndex start = columnStart[iColumn]; 7584 columnStart[iColumn]=put; 7585 memcpy(element+put,element+start,length*sizeof(double)); 7586 memcpy(row+put,row+start,length*sizeof(int)); 7587 // forward backward 7588 int iLast=backward[iColumn]; 7589 int iNext=forward[iColumn]; 7590 forward[iLast]=iNext; 7591 backward[iNext]=iLast; 7592 forward[last]=iColumn; 7593 backward[iColumn]=last; 7594 forward[iColumn]=numberColumns; 7595 backward[numberColumns]=iColumn; 7596 } else { 7597 needCompaction=true; 7598 } 7599 } else { 7600 needCompaction=true; 7601 } 7602 if (needCompaction) { 7603 printf("compacting\n"); 7604 // size is lastElement+numberElementsOriginal 7605 #ifndef NDEBUG 7606 CoinBigIndex total=lengthNeeded-columnLength[iColumn]; 7607 for (int i=0;i<numberColumns;i++) 7608 total += columnLength[i]; 7609 assert (total<=numberElementsOriginal+lengthNeeded); 7610 #endif 7611 CoinBigIndex put=lastElement; 7612 for (int i=0;i<numberColumns;i++) { 7613 CoinBigIndex start = columnStart[i]; 7614 columnStart[i]=put; 7615 int n=columnLength[i]; 7616 memcpy(element+put,element+start,n*sizeof(double)); 7617 memcpy(row+put,row+start,n*sizeof(int)); 7618 put += n; 7619 } 7620 // replace length (may mean copying uninitialized) 7621 columnLength[iColumn]=lengthNeeded; 7622 int spare = (2*lastElement-put-(lengthNeeded-length)-numberElementsOriginal)/numberColumns; 7623 assert (spare>=0); 7624 // copy back 7625 put=0; 7626 for (int i=0;i<numberColumns;i++) { 7627 CoinBigIndex start = columnStart[i]; 7628 columnStart[i]=put; 7629 int n=columnLength[i]; 7630 memcpy(element+put,element+start,n*sizeof(double)); 7631 memcpy(row+put,row+start,n*sizeof(int)); 7632 put += n+spare; 7633 } 7634 assert (put<=lastElement); 7635 columnLength[iColumn]=length; 7636 // redo forward,backward 7637 for (int i=-1;i<numberColumns;i++) 7638 forward[i]=i+1; 7639 forward[numberColumns]=-1; 7640 for (int i=0;i<=numberColumns;i++) 7641 backward[i]=i-1; 7642 backward[-1]=-1; 7643 } 7644 //abort(); 7645 #if DEBUG_SOME > 0 7646 printf("moved column %d\n",iColumn); 7647 #endif 7648 } 7649 #if DEBUG_SOME > 0 7650 #ifndef NDEBUG 7651 static void checkBasis(ClpSimplex * model,char * rowType, char * columnType) 7652 { 7653 int numberRows=model->numberRows(); 7654 int nRowBasic=0; 7655 int nRows=0; 7656 for (int i=0;i<numberRows;i++) { 7657 if (rowType[i]<=0||rowType[i]==55) { 7658 nRows++; 7659 if(model->getRowStatus(i)==ClpSimplex::basic) 7660 nRowBasic++; 7661 } 7662 } 7663 int numberColumns=model->numberColumns(); 7664 int nColumnBasic=0; 7665 for (int i=0;i<numberColumns;i++) { 7666 if ((columnType[i]<11||columnType[i]==55)&&model->getColumnStatus(i)==ClpSimplex::basic) 7667 nColumnBasic++; 7668 } 7669 ClpTraceDebug (nRowBasic+nColumnBasic==nRows); 7670 } 7671 #endif 7672 #endif 7673 #if DEBUG_SOME > 0 7674 static int xxxxxx=2999999; 7675 #endif 7676 /* Mini presolve (faster) 7677 Char arrays must be numberRows and numberColumns long 7678 on entry second part must be filled in as follows - 7679 0 - possible 7680 >0 - take out and do something (depending on value - TBD) 7681 1 - redundant row 7682 2 - sub 7683 11 - column can be moved to bound 7684 4 - row redundant (duplicate) 7685 13 - empty (or initially fixed) column 7686 14 - == row (also column deleted by row) 7687 3 - column altered by a 14 7688 5 - temporary marker for truly redundant sub row 7689 8 - other 7690 -1 row/column can't vanish but can have entries removed/changed 7691 -2 don't touch at all 7692 on exit <=0 ones will be in presolved problem 7693 struct will be created and will be long enough 7694 (information on length etc in first entry) 7695 user must delete struct 7696 */ 7697 ClpSimplex * 7698 ClpSimplex::miniPresolve(char * rowType, char * columnType,void ** infoOut) 7699 { 7700 // Big enough structure 7701 int numberTotal=numberRows_+numberColumns_; 7702 CoinBigIndex lastElement = matrix_->getNumElements(); 7703 int maxInfoStuff = 5*lastElement*sizeof(double)+numberTotal*sizeof(clpPresolveInfo2); 7704 clpPresolveInfo * infoA = new clpPresolveInfo[numberTotal]; 7705 char * startStuff = new char [maxInfoStuff]; 7706 memset(infoA,'B',numberTotal*sizeof(clpPresolveInfo)); 7707 memset(startStuff,'B',maxInfoStuff); 7708 int nActions=0; 7709 int * whichRows = new int [2*numberRows_+numberColumns_]; 7710 int * whichColumns = whichRows + numberRows_; 7711 int * whichRows2 = whichColumns + numberColumns_; 7712 double * array = new double [numberRows_]; 7713 memset(array,0,numberRows_*sizeof(double)); 7714 // New model (put in modification to increase size of matrix) and pack 7715 bool needExtension=numberColumns_>matrix_->getNumCols(); 7716 if (needExtension) { 7717 matrix()->reserve(numberColumns_,lastElement,true); 7718 CoinBigIndex * columnStart = matrix()->getMutableVectorStarts(); 7719 for (int i=numberColumns_;i>=0;i--) { 7720 if (columnStart[i]==0) 7721 columnStart[i]=lastElement; 7722 else 7723 break; 7724 } 7725 assert (lastElement==columnStart[numberColumns_]); 7726 } 7727 #define TWOFER 7728 #ifdef TWOFER 7729 ClpSimplex * newModel = NULL; 7730 CoinBigIndex lastPossible=3*lastElement; 7731 CoinBigIndex lastGood=2*lastElement; 7732 clpPresolveMore moreInfo; 7733 moreInfo.model=NULL; 7734 moreInfo.rowType=rowType; 7735 moreInfo.columnType=columnType; 7736 int addColumns = eventHandler_->eventWithInfo(ClpEventHandler::modifyMatrixInMiniPresolve,&moreInfo); 7737 if (moreInfo.model) { 7738 newModel = moreInfo.model; 7739 } else { 7740 newModel = new ClpSimplex(*this); 7741 newModel->matrix()->reserve(numberColumns_+addColumns,lastPossible,true); 7742 } 7743 #else 7744 ClpSimplex * newModel = new ClpSimplex(*this); 7745 //newModel->matrix()->reserve(numberColumns_,lastElement,true); 7746 #endif 7747 double * rowLower = newModel->rowLower(); 7748 double * rowUpper = newModel->rowUpper(); 7749 //double * rowActivity = newModel->primalRowSolution(); 7750 //unsigned char * rowStatus = newModel->statusArray()+numberColumns_; 7751 double * columnLower = newModel->columnLower(); 7752 double * columnUpper = newModel->columnUpper(); 7753 //double * columnActivity = newModel->primalColumnSolution(); 7754 //unsigned char * columnStatus = newModel->statusArray(); 7755 // Take out marked stuff 7756 saveInfo stuff; 7757 stuff.putStuff=startStuff; 7758 stuff.startStuff=startStuff; 7759 stuff.maxStuff=maxInfoStuff; 7760 CoinPackedMatrix * matrix = newModel->matrix(); 7761 int * row = matrix->getMutableIndices(); 7762 CoinBigIndex * columnStart = matrix->getMutableVectorStarts(); 7763 int * columnLength = matrix->getMutableVectorLengths(); 7764 double * element = matrix->getMutableElements(); 7765 // get row copy 7766 CoinPackedMatrix rowCopy = *matrix; 7767 rowCopy.reverseOrdering(); 7768 int * column = rowCopy.getMutableIndices(); 7769 CoinBigIndex * rowStart = rowCopy.getMutableVectorStarts(); 7770 double * elementByRow = rowCopy.getMutableElements(); 7771 int * rowLength = rowCopy.getMutableVectorLengths(); 7772 for (int iRow=0;iRow<numberRows_;iRow++) { 7773 if (rowType[iRow]>0) { 7774 clpPresolveInfo1_4_8 thisInfo; 7775 thisInfo.row=iRow; 7776 thisInfo.oldRowLower=(rowLower_[iRow]>-1.0e30) ? rowLower_[iRow]-rowLower[iRow] : rowLower[iRow]; 7777 thisInfo.oldRowUpper=(rowUpper_[iRow]<1.0e30) ? rowUpper_[iRow]-rowUpper[iRow] : rowUpper[iRow]; 7778 int n=rowLength[iRow]; 7779 CoinBigIndex start=rowStart[iRow]; 7780 thisInfo.lengthRow=n; 7781 //thisInfo.column=-1; 7782 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 7783 infoA[nActions].type=4; //rowType[iRow]; 7784 nActions++; 7785 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo1_4_8), 7786 n,column+start,elementByRow+start); 7787 } 7788 } 7789 CoinBigIndex put=0; 7790 bool anyDeleted=false; 7791 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 7792 CoinBigIndex start=columnStart[iColumn]; 7793 int length = columnLength[iColumn]; 7794 if (columnType[iColumn]>0||(columnType[iColumn]==0&& 7795 (!length||columnLower_[iColumn]==columnUpper_[iColumn]))) { 7796 clpPresolveInfo13 thisInfo; 7797 //thisInfo.row=-1; 7798 thisInfo.oldColumnLower=columnLower[iColumn]; 7799 thisInfo.oldColumnUpper=columnUpper[iColumn]; 7800 thisInfo.column=iColumn; 7801 CoinBigIndex start=columnStart[iColumn]; 7802 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 7803 infoA[nActions].type=(columnType[iColumn]>0) ? columnType[iColumn] : 13; 7804 nActions++; 7805 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo13), 7806 0,NULL,NULL); 7807 columnType[iColumn]=13; 7808 if (length) { 7809 double solValue=columnLower[iColumn]; 7810 if (solValue) { 7811 for (int j=start;j<start+length;j++) { 7812 int iRow=row[j]; 7813 double value = element[j]*solValue; 7814 double lower = rowLower[iRow]; 7815 if (lower>-1.0e20) 7816 rowLower[iRow]=lower-value; 7817 double upper = rowUpper[iRow]; 7818 if (upper<1.0e20) 7819 rowUpper[iRow]=upper-value; 7820 } 7821 } 7822 anyDeleted=true; 7823 length=0; 7824 } 7825 } 7826 columnStart[iColumn]=put; 7827 for (CoinBigIndex i=start;i<start+length;i++) { 7828 int iRow=row[i]; 7829 if (rowType[iRow]<=0) { 7830 row[put]=iRow; 7831 element[put++]=element[i]; 7832 } 7833 } 7834 columnLength[iColumn]=put-columnStart[iColumn]; 7835 } 7836 int numberInitial=nActions; 7837 columnStart[numberColumns_]=put; 7838 matrix->setNumElements(put); 7839 // get row copy if changed 7840 if (anyDeleted) { 7841 rowCopy = *matrix; 7842 rowCopy.reverseOrdering(); 7843 column = rowCopy.getMutableIndices(); 7844 rowStart = rowCopy.getMutableVectorStarts(); 7845 elementByRow = rowCopy.getMutableElements(); 7846 rowLength = rowCopy.getMutableVectorLengths(); 7847 } 7848 double * objective = newModel->objective(); 7849 double offset = objectiveOffset(); 7850 int numberRowsLook=0; 7851 #ifdef TWOFER 7852 bool orderedMatrix=true; 7853 #endif 7854 int nChanged = 1; 7855 bool feasible=true; 7856 for (int iRow=0;iRow<numberRows_;iRow++) { 7857 #if DEBUG_SOME>0 7858 #ifndef NDEBUG 7859 checkBasis(newModel, rowType, columnType); 7860 #endif 7861 #endif 7862 int nPoss=2; 7863 #if DEBUG_SOME > 0 7864 xxxxxx--; 7865 if (xxxxxx<=0) 7866 nPoss=1; 7867 if (xxxxxx<-1000) 7868 nPoss=-1; 7869 if (xxxxxx==1) { 7870 printf("bad\n"); 7871 } 7872 #endif 7873 if (rowLength[iRow]<=nPoss&&!rowType[iRow]) { 7874 if (rowLength[iRow]<=1) { 7875 if (rowLength[iRow]==1) { 7876 whichRows[numberRowsLook++]=iRow; 7877 } else { 7878 #if DEBUG_SOME > 0 7879 printf("Dropping null row %d (status %d) - nActions %d\n", 7880 iRow,getRowStatus(iRow),nActions); 7881 #endif 7882 if (rowLower[iRow]>0.0||rowUpper[iRow]<0.0) { 7883 feasible=false; 7884 nChanged=-1; 7885 numberRowsLook=0; 7886 break; 7887 } 7888 rowType[iRow]=1; 7889 clpPresolveInfo1_4_8 thisInfo; 7890 thisInfo.oldRowLower=(rowLower_[iRow]>-1.0e30) ? rowLower_[iRow]-rowLower[iRow] : rowLower[iRow]; 7891 thisInfo.oldRowUpper=(rowUpper_[iRow]<1.0e30) ? rowUpper_[iRow]-rowUpper[iRow] : rowUpper[iRow]; 7892 thisInfo.row=iRow; 7893 int n=rowLength[iRow]; 7894 CoinBigIndex start=rowStart[iRow]; 7895 thisInfo.lengthRow=n; 7896 //thisInfo.column=-1; 7897 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 7898 infoA[nActions].type=1; 7899 nActions++; 7900 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo1_4_8), 7901 n,column+start,elementByRow+start); 7902 } 7903 #ifdef TWOFER 7904 } else if (rowLower[iRow]==rowUpper[iRow]) { 7905 whichRows[numberRowsLook++]=iRow; 7906 CoinBigIndex start = rowStart[iRow]; 7907 int iColumn1 = column[start]; 7908 double value1 = elementByRow[start]; 7909 int iColumn2 = column[start+1]; 7910 double value2 = elementByRow[start+1]; 7911 bool swap=false; 7912 double ratio = fabs(value1/value2); 7913 if (ratio<0.001||ratio>1000.0) { 7914 if (fabs(value1)<fabs(value2)) { 7915 swap=true; 7916 } 7917 } else if (columnLength[iColumn1]<columnLength[iColumn2]) { 7918 swap=true; 7919 } 7920 if(swap) { 7921 iColumn1 = iColumn2; 7922 value1 = value2; 7923 iColumn2 = column[start]; 7924 value2 = elementByRow[start]; 7925 } 7926 // column bounds 7927 double dropLower = columnLower[iColumn2]; 7928 double dropUpper = columnUpper[iColumn2]; 7929 double newLower; 7930 double newUpper; 7931 double rhs = rowLower[iRow]/value1; 7932 double multiplier = value2/value1; 7933 if (multiplier>0.0) { 7934 newLower = (dropUpper<1.0e30) ? rhs - multiplier*dropUpper : -COIN_DBL_MAX; 7935 newUpper = (dropLower>-1.0e30) ? rhs - multiplier*dropLower : COIN_DBL_MAX; 7936 } else { 7937 newUpper = (dropUpper<1.0e30) ? rhs - multiplier*dropUpper : COIN_DBL_MAX; 7938 newLower = (dropLower>-1.0e30) ? rhs - multiplier*dropLower : -COIN_DBL_MAX; 7939 } 7940 //columnType[iColumn1]=3; 7941 columnType[iColumn2]=14; 7942 rowType[iRow]=14; 7943 rhs = rowLower[iRow]/value2; 7944 multiplier = value1/value2; 7945 clpPresolveInfo14 thisInfo; 7946 thisInfo.oldColumnLower=columnLower[iColumn1]; 7947 thisInfo.oldColumnUpper=columnUpper[iColumn1]; 7948 thisInfo.oldColumnLower2=columnLower[iColumn2]; 7949 thisInfo.oldColumnUpper2=columnUpper[iColumn2]; 7950 thisInfo.oldObjective2=objective[iColumn2]; 7951 thisInfo.value1=value1; 7952 thisInfo.rhs=rowLower[iRow]; 7953 thisInfo.row=iRow; 7954 thisInfo.column=iColumn1; 7955 thisInfo.column2=iColumn2; 7956 int nel=columnLength[iColumn2]; 7957 CoinBigIndex startCol=columnStart[iColumn2]; 7958 thisInfo.lengthColumn2=nel; 7959 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 7960 infoA[nActions].type=14; 7961 nActions++; 7962 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo14), 7963 nel,row+startCol,element+startCol); 7964 newLower = CoinMax(newLower,columnLower[iColumn1]); 7965 newUpper = CoinMin(newUpper,columnUpper[iColumn1]); 7966 if (newLower>newUpper+primalTolerance_) { 7967 feasible=false; 7968 nChanged=-1; 7969 numberRowsLook=0; 7970 break; 7971 } 7972 columnLower[iColumn1]=newLower; 7973 columnUpper[iColumn1]=newUpper; 7974 #if DEBUG_SOME > 0 7975 printf("Dropping doubleton row %d (status %d) keeping column %d (status %d) dropping %d (status %d) (mult,rhs %g %g) - nActions %d\n", 7976 iRow,getRowStatus(iRow),iColumn1,getColumnStatus(iColumn1),iColumn2,getColumnStatus(iColumn2),multiplier,rhs,nActions); 7977 #endif 7978 objective[iColumn1] -= objective[iColumn2]*multiplier; 7979 offset -= rowLower[iRow]*(objective[iColumn2]*multiplier); 7980 bool needDrop=false; 7981 if (newModel->getRowStatus(iRow)!=basic) { 7982 if (newModel->getColumnStatus(iColumn2)!=basic) { 7983 // On way back may as well have column basic 7984 newModel->setColumnStatus(iColumn2,basic); 7985 // need to drop basic 7986 if (newModel->getColumnStatus(iColumn1)==basic) { 7987 //setColumnStatus(iColumn1,superBasic); 7988 newModel->setColumnStatus(iColumn1,superBasic); 7989 } else { 7990 // not much we can do 7991 #if DEBUG_SOME > 0 7992 printf("dropping but no basic a\n"); 7993 #endif 7994 } 7995 } else { 7996 // looks good 7997 } 7998 } else { 7999 if (newModel->getColumnStatus(iColumn2)!=basic) { 8000 // looks good 8001 } else { 8002 // need to keep a basic 8003 if (newModel->getColumnStatus(iColumn1)!=basic) { 8004 //setColumnStatus(iColumn2,superBasic); 8005 //setColumnStatus(iColumn1,basic); 8006 newModel->setColumnStatus(iColumn1,basic); 8007 } else { 8008 // not much we can do 8009 #if DEBUG_SOME > 0 8010 printf("dropping but all basic a\n"); 8011 #endif 8012 needDrop=true; 8013 //setColumnStatus(iColumn2,superBasic); 8014 } 8015 } 8016 } 8017 int n=0; 8018 start = columnStart[iColumn1]; 8019 for (int i=start;i<start+columnLength[iColumn1]; 8020 i++) { 8021 int jRow=row[i]; 8022 if (jRow!=iRow) { 8023 array[jRow]=element[i]; 8024 whichRows2[n++]=jRow; 8025 } 8026 } 8027 rowLength[iRow]=0; 8028 start = columnStart[iColumn2]; 8029 for (int i=start;i<start+columnLength[iColumn2]; 8030 i++) { 8031 int jRow=row[i]; 8032 if (jRow!=iRow) { 8033 double value = array[jRow]; 8034 double valueNew = value -multiplier*element[i]; 8035 double rhsMod = rhs*element[i]; 8036 if (rowLower[jRow]>-1.0e30) 8037 rowLower[jRow] -= rhsMod; 8038 if (rowUpper[jRow]<1.0e30) 8039 rowUpper[jRow] -= rhsMod; 8040 if (!value) { 8041 array[jRow]=valueNew; 8042 whichRows2[n++]=jRow; 8043 } else { 8044 if (!valueNew) 8045 valueNew=1.0e-100; 8046 array[jRow]=valueNew; 8047 } 8048 } 8049 } 8050 columnLength[iColumn2]=0; 8051 start = columnStart[iColumn1]; 8052 if (n>columnLength[iColumn1]) { 8053 orderedMatrix=false; 8054 if (lastElement+n>lastGood) { 8055 // pack down 8056 CoinBigIndex put=lastElement; 8057 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8058 CoinBigIndex start = columnStart[iColumn]; 8059 columnStart[iColumn]=put-lastElement; 8060 int length=columnLength[iColumn]; 8061 for (CoinBigIndex j=start;j<start+length;j++) { 8062 row[put]=row[j]; 8063 element[put++]=element[j]; 8064 } 8065 } 8066 int numberElements = put-lastElement; 8067 columnStart[numberColumns_]=numberElements; 8068 memcpy(row,row+lastElement,numberElements*sizeof(int)); 8069 memcpy(element,element+lastElement,numberElements*sizeof(double)); 8070 lastElement=numberElements; 8071 } 8072 start=lastElement; 8073 columnStart[iColumn1]=start; 8074 } 8075 CoinBigIndex put = start; 8076 for (int i=0;i<n;i++) { 8077 int jRow = whichRows2[i]; 8078 double value = array[jRow]; 8079 //#define FUNNY_CHECK 8080 #ifdef FUNNY_CHECK 8081 if (rowType[jRow]==-2) 8082 printf("iColumn1 %d iColumn2 %d row %d value %g\n", 8083 iColumn1,iColumn2,jRow,value); 8084 #endif 8085 array[jRow]=0.0; 8086 if (fabs(value)<1.0e-13) 8087 value=0.0; 8088 if (value) { 8089 row[put] = jRow; 8090 element[put++]=value; 8091 if(needDrop&&newModel->getRowStatus(jRow)!=basic) { 8092 newModel->setRowStatus(jRow,basic); 8093 needDrop=false; 8094 } 8095 } 8096 // take out of row copy 8097 int startR = rowStart[jRow]; 8098 int putR=startR; 8099 for (int i=startR;i<startR+rowLength[jRow];i++) { 8100 int iColumn = column[i]; 8101 if (iColumn!=iColumn1&&iColumn!=iColumn2) { 8102 column[putR]=iColumn; 8103 elementByRow[putR++]=elementByRow[i]; 8104 } else if (value) { 8105 column[putR]=iColumn1; 8106 elementByRow[putR++]=value; 8107 value=0.0; 8108 } 8109 } 8110 int rowLength2=putR-startR; 8111 if (rowLength2<=1&&rowLength[jRow]>1&&jRow<iRow) { 8112 // may be interesting 8113 whichRows[numberRowsLook++]=jRow; 8114 } 8115 rowLength[jRow]=rowLength2; 8116 } 8117 columnLength[iColumn1]=put-start; 8118 lastElement=CoinMax(lastElement,put); 8119 #endif 8120 } 8121 } 8122 } 8123 #if DEBUG_SOME>0 8124 if (xxxxxx<-1000) 8125 nChanged=-1; 8126 #endif 8127 while (nChanged>0) { 8128 nChanged=0; 8129 int numberColumnsLook=0; 8130 for (int i=0;i<numberRowsLook;i++) { 8131 #if DEBUG_SOME>0 8132 #ifndef NDEBUG 8133 checkBasis(newModel, rowType, columnType); 8134 #endif 8135 #endif 8136 int iRow = whichRows[i]; 8137 if (rowLength[iRow]==1) { 8138 //rowType[iRow]=55; 8139 int iColumn = column[rowStart[iRow]]; 8140 if (/*columnType[iColumn]==14||*/columnType[iColumn]<-1) 8141 continue; 8142 if (!columnType[iColumn]) { 8143 columnType[iColumn]=55; 8144 whichColumns[numberColumnsLook++]=iColumn; 8145 nChanged++; 8146 } 8147 #if 0 8148 } else if (rowLength[iRow]==2) { 8149 if (rowLower[iRow]==rowUpper[iRow]) { 8150 CoinBigIndex start = rowStart[iRow]; 8151 int iColumn1 = column[start]; 8152 int iColumn2 = column[start+1]; 8153 if (!columnType[iColumn1]&&!columnType[iColumn2]) { 8154 if (fabs(elementByRow[start])<fabs(elementByRow[start+1])) { 8155 iColumn1 = iColumn2; 8156 iColumn2 = column[start]; 8157 } 8158 // iColumn2 will be deleted 8159 columnType[iColumn1]=56; 8160 columnType[iColumn2]=56; 8161 } 8162 } else { 8163 printf("why here in miniPresolve row %d\n",iRow); 8164 } 8165 #endif 8166 } 8167 } 8168 // mark one row as ub and take out row in column copy 8169 numberRowsLook=0; 8170 for (int iLook=0;iLook<numberColumnsLook;iLook++) { 8171 nChanged++; 8172 int iColumn = whichColumns[iLook]; 8173 if (columnType[iColumn]!=55&&columnType[iColumn]>10) 8174 continue; 8175 if (columnType[iColumn]==55) 8176 columnType[iColumn]=0; 8177 CoinBigIndex start=columnStart[iColumn]; 8178 int jRowLower=-1; 8179 double newLower=columnLower[iColumn]; 8180 int jRowUpper=-1; 8181 double newUpper=columnUpper[iColumn]; 8182 double coefficientLower=0.0; 8183 double coefficientUpper=0.0; 8184 for (CoinBigIndex i=start;i<start+columnLength[iColumn]; 8185 i++) { 8186 int iRow=row[i]; 8187 if (rowLength[iRow]==1&&rowType[iRow]==0) { 8188 rowType[iRow]=1; 8189 assert(columnType[iColumn]>=0); 8190 // adjust bounds 8191 double value = elementByRow[rowStart[iRow]]; 8192 double lower = newLower; 8193 double upper = newUpper; 8194 if (value>0.0) { 8195 if (rowUpper[iRow]<1.0e30) 8196 upper = rowUpper[iRow]/value; 8197 if (rowLower[iRow]>-1.0e30) 8198 lower = rowLower[iRow]/value; 8199 } else { 8200 if (rowUpper[iRow]<1.0e30) 8201 lower = rowUpper[iRow]/value; 8202 if (rowLower[iRow]>-1.0e30) 8203 upper = rowLower[iRow]/value; 8204 } 8205 if (lower>newLower+primalTolerance_) { 8206 if (lower>newUpper+primalTolerance_) { 8207 feasible=false; 8208 nChanged=-1; 8209 numberColumnsLook=0; 8210 break; 8211 } else if (lower>newUpper-primalTolerance_) { 8212 newLower=newUpper; 8213 } else { 8214 newLower = CoinMax(lower,newLower); 8215 } 8216 jRowLower=iRow; 8217 coefficientLower=value; 8218 } 8219 if (upper<newUpper-primalTolerance_) { 8220 if (upper<newLower-primalTolerance_) { 8221 feasible=false; 8222 nChanged=-1; 8223 numberColumnsLook=0; 8224 break; 8225 } else if (upper<newLower+primalTolerance_) { 8226 newUpper = newLower; 8227 } else { 8228 newUpper = CoinMin(upper,newUpper); 8229 } 8230 jRowUpper=iRow; 8231 coefficientUpper=value; 8232 } 8233 } 8234 } 8235 int put=start; 8236 int iFlag=0; 8237 int nonFree=0; 8238 int numberNonBasicSlacksOut=0; 8239 if (jRowLower>=0||jRowUpper>=0) { 8240 clpPresolveInfo2 thisInfo; 8241 if (jRowLower>=0) { 8242 thisInfo.oldRowLower=(rowLower_[jRowLower]>-1.0e30) ? rowLower_[jRowLower]-rowLower[jRowLower] : rowLower[jRowLower]; 8243 thisInfo.oldRowUpper=(rowUpper_[jRowLower]<1.0e30) ? rowUpper_[jRowLower]-rowUpper[jRowLower] : rowUpper[jRowLower]; 8244 thisInfo.row=jRowLower; 8245 thisInfo.coefficient=coefficientLower; 8246 } else { 8247 thisInfo.row=-1; 8248 #ifndef NDEBUG 8249 thisInfo.oldRowLower=COIN_DBL_MAX; 8250 thisInfo.oldRowUpper=-COIN_DBL_MAX; 8251 thisInfo.coefficient=0.0; 8252 #endif 8253 } 8254 if (jRowUpper>=0&&jRowLower!=jRowUpper) { 8255 thisInfo.oldRowLower2=(rowLower_[jRowUpper]>-1.0e30) ? rowLower_[jRowUpper]-rowLower[jRowUpper] : rowLower[jRowUpper]; 8256 thisInfo.oldRowUpper2=(rowUpper_[jRowUpper]<1.0e30) ? rowUpper_[jRowUpper]-rowUpper[jRowUpper] : rowUpper[jRowUpper]; 8257 thisInfo.row2=jRowUpper; 8258 thisInfo.coefficient2=coefficientUpper; 8259 } else { 8260 thisInfo.row2=-1; 8261 #ifndef NDEBUG 8262 thisInfo.oldRowLower2=COIN_DBL_MAX; 8263 thisInfo.oldRowUpper2=-COIN_DBL_MAX; 8264 thisInfo.coefficient2=0.0; 8265 #endif 8266 } 8267 thisInfo.oldColumnLower=columnLower[iColumn]; 8268 thisInfo.oldColumnUpper=columnUpper[iColumn]; 8269 columnLower[iColumn]=newLower; 8270 columnUpper[iColumn]=newUpper; 8271 thisInfo.column=iColumn; 8272 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 8273 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo2), 8274 0,NULL,NULL); 8275 infoA[nActions].type=2; 8276 nActions++; 8277 } 8278 for (CoinBigIndex i=start;i<start+columnLength[iColumn]; 8279 i++) { 8280 int iRow=row[i]; 8281 if (rowLength[iRow]==1&&rowType[iRow]>=0&&rowType[iRow]!=5) { 8282 #if DEBUG_SOME > 0 8283 printf("Dropping singleton row %d (status %d) because of column %d (status %d) - jRow lower/upper %d/%d - nActions %d\n", 8284 iRow,getRowStatus(iRow),iColumn, 8285 getColumnStatus(iColumn),jRowLower,jRowUpper,nActions); 8286 #endif 8287 if (newModel->getRowStatus(iRow)!=basic) 8288 numberNonBasicSlacksOut++; 8289 rowType[iRow]=1; 8290 if (iRow!=jRowLower&&iRow!=jRowUpper) { 8291 // mark as redundant 8292 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 8293 clpPresolveInfo1_4_8 thisInfo; 8294 thisInfo.oldRowLower=(rowLower_[iRow]>-1.0e30) ? rowLower_[iRow]-rowLower[iRow] : rowLower[iRow]; 8295 thisInfo.oldRowUpper=(rowUpper_[iRow]<1.0e30) ? rowUpper_[iRow]-rowUpper[iRow] : rowUpper[iRow]; 8296 thisInfo.row=iRow; 8297 int n=rowLength[iRow]; 8298 CoinBigIndex start=rowStart[iRow]; 8299 thisInfo.lengthRow=n; 8300 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo1_4_8), 8301 n,column+start,elementByRow+start); 8302 infoA[nActions].type=1; 8303 nActions++; 8304 } 8305 rowLength[iRow]=0; 8306 } else if (rowType[iRow]<=0) { 8307 row[put]=iRow; 8308 double value = element[i]; 8309 element[put++]=value; 8310 if (rowType[iRow]>=0&&iFlag<3) { 8311 assert(rowType[iRow]==0); 8312 double lower = rowLower[iRow]; 8313 double upper = rowUpper[iRow]; 8314 if (-1.0e20 < lower && upper < 1.0e20) { 8315 // bounded - we lose 8316 iFlag=-1; 8317 //break; 8318 } else if (-1.0e20 < lower || upper < 1.0e20) { 8319 nonFree++; 8320 } 8321 // see what this particular row says 8322 // jFlag == 2 ==> up is towards feasibility 8323 int jFlag = (value > 0.0 8324 ? (upper > 1.0e20 ? 2 : 1) 8325 : (lower < -1.0e20 ? 2 : 1)); 8326 8327 if (iFlag) { 8328 // check that it agrees with iFlag. 8329 if (iFlag!=jFlag) { 8330 iFlag=-1; 8331 } 8332 } else { 8333 // first row -- initialize iFlag 8334 iFlag=jFlag; 8335 } 8336 } else if (rowType[iRow]<0) { 8337 iFlag=-1; // be safe 8338 } 8339 } 8340 } 8341 // Do we need to switch status of iColumn? 8342 if (numberNonBasicSlacksOut>0) { 8343 // make iColumn non basic if possible 8344 if (newModel->getColumnStatus(iColumn)==basic) { 8345 newModel->setColumnStatus(iColumn,superBasic); 8346 } 8347 } 8348 double cost = objective[iColumn]*optimizationDirection_; 8349 int length = put-columnStart[iColumn]; 8350 if (!length) { 8351 if (!cost) { 8352 // put to closest to zero 8353 if (fabs(columnLower[iColumn])<fabs(columnUpper[iColumn])) 8354 iFlag=1; 8355 else 8356 iFlag=2; 8357 } else if (cost>0.0) { 8358 iFlag=1; 8359 } else { 8360 iFlag=2; 8361 } 8362 } else { 8363 if (cost>0.0&&iFlag==2) 8364 iFlag=-1; 8365 else if (cost<0.0&&iFlag==1) 8366 iFlag=-1; 8367 } 8368 columnLength[iColumn]=length; 8369 if (iFlag>0&&nonFree) { 8370 double newValue; 8371 if (iFlag==2) { 8372 // fix to upper 8373 newValue =CoinMin(columnUpper[iColumn],1.0e20); 8374 } else { 8375 // fix to lower 8376 newValue =CoinMax(columnLower[iColumn],-1.0e20); 8377 } 8378 columnActivity_[iColumn]=newValue; 8379 #if DEBUG_SOME > 0 8380 if (newModel->getColumnStatus(iColumn)== 8381 basic) { 8382 // ? move basic back onto sub if can? 8383 iFlag += 2; 8384 } 8385 printf("Dropping movable column %d - iFlag %d - jRow lower/upper %d/%d - nActions %d\n", 8386 iColumn,iFlag,jRowLower,jRowUpper,nActions); 8387 #endif 8388 columnType[iColumn]=11; 8389 if (newModel->getColumnStatus(iColumn)== 8390 basic) { 8391 // need to put status somewhere else 8392 int shortestNumber=numberColumns_; 8393 int shortest=-1; 8394 for (int j=start;j<start+length;j++) { 8395 int iRow=row[j]; 8396 if (rowLength[iRow]<shortestNumber&& 8397 newModel->getRowStatus(iRow)!= 8398 basic) { 8399 shortest=iRow; 8400 shortestNumber = rowLength[iRow]; 8401 } 8402 } 8403 if (shortest>=0) { 8404 // make basic 8405 newModel->setRowStatus(shortest,basic); 8406 newModel->setColumnStatus(iColumn,superBasic); 8407 } else { 8408 // put on a column 8409 shortestNumber=numberColumns_; 8410 shortest=-1; 8411 for (int j=start;j<start+length;j++) { 8412 int iRow=row[j]; 8413 if (rowLength[iRow]<shortestNumber) { 8414 int start = rowStart[iRow]; 8415 for (int i=start;i<start+rowLength[iRow];i++) { 8416 int jColumn = column[i]; 8417 if (iColumn!=jColumn&& 8418 newModel->getColumnStatus(jColumn)!= 8419 basic) { 8420 shortest=jColumn; 8421 shortestNumber = rowLength[iRow]; 8422 } 8423 } 8424 } 8425 } 8426 if (shortest>=0) { 8427 // make basic 8428 newModel->setColumnStatus(shortest,basic); 8429 } else { 8430 #if DEBUG_SOME > 0 8431 printf("what now - dropping - basic\n"); 8432 #endif 8433 } 8434 } 8435 } 8436 clpPresolveInfo11 thisInfo; 8437 thisInfo.oldColumnLower=columnLower[iColumn]; 8438 thisInfo.oldColumnUpper=columnUpper[iColumn]; 8439 thisInfo.fixedTo=newValue; 8440 columnLower[iColumn]=newValue; 8441 columnUpper[iColumn]=newValue; 8442 thisInfo.column=iColumn; 8443 int n=columnLength[iColumn]; 8444 CoinBigIndex start=columnStart[iColumn]; 8445 thisInfo.lengthColumn=n; 8446 infoA[nActions].infoOffset=stuff.putStuff-startStuff; 8447 infoA[nActions].type=11; 8448 nActions++; 8449 ClpCopyToMiniSave(stuff,reinterpret_cast<char *>(&thisInfo),sizeof(clpPresolveInfo11), 8450 n,row+start,element+start); 8451 // adjust rhs and take out of rows 8452 columnLength[iColumn]=0; 8453 nChanged++; 8454 for (int j=start;j<start+length;j++) { 8455 int iRow=row[j]; 8456 double value = element[j]*newValue; 8457 double lower = rowLower[iRow]; 8458 if (lower>-1.0e20) 8459 rowLower[iRow]=lower-value; 8460 double upper = rowUpper[iRow]; 8461 if (upper<1.0e20) 8462 rowUpper[iRow]=upper-value; 8463 // take out of row copy (and put on list) 8464 assert (rowType[iRow]<=0&&rowType[iRow]>-2); 8465 //if (!rowType[iRow]) { 8466 //rowType[iRow]=-1; 8467 whichRows[numberRowsLook++]=iRow; 8468 //} 8469 int start = rowStart[iRow]; 8470 int put=start; 8471 for (int i=start;i<start+rowLength[iRow];i++) { 8472 int jColumn = column[i]; 8473 if (iColumn!=jColumn) { 8474 column[put]=jColumn; 8475 elementByRow[put++]=elementByRow[i]; 8476 } 8477 } 8478 rowLength[iRow]=put-start; 8479 } 8480 } 8481 } 8482 } 8483 if (nActions&&feasible) { 8484 clpPresolveMore moreInfo; 8485 moreInfo.model=newModel; 8486 moreInfo.rowCopy=&rowCopy; 8487 moreInfo.rowType=rowType; 8488 moreInfo.columnType=columnType; 8489 moreInfo.stuff=&stuff; 8490 moreInfo.info=infoA; 8491 moreInfo.nActions=&nActions; 8492 eventHandler_->eventWithInfo(ClpEventHandler::moreMiniPresolve,&moreInfo); 8493 newModel->setObjectiveOffset(offset); 8494 int nChar2 = nActions*sizeof(clpPresolveInfo)+(stuff.putStuff-startStuff); 8495 clpPresolveInfo * infoData = reinterpret_cast<clpPresolveInfo *>(new char[nChar2]); 8496 memcpy(infoData,infoA,nActions*sizeof(clpPresolveInfo)); 8497 char * info2 = reinterpret_cast<char *>(infoData+nActions); 8498 memcpy(info2,startStuff,stuff.putStuff-startStuff); 8499 listInfo * infoNew = new listInfo; 8500 infoNew->numberEntries=nActions; 8501 infoNew->maximumEntries=nActions; 8502 infoNew->start=infoData; 8503 infoNew->numberInitial=numberInitial; 8504 *infoOut=infoNew; 8505 int nRows=0; 8506 for (int iRow=0;iRow<numberRows_;iRow++) { 8507 if (rowType[iRow]>0) 8508 whichRows[nRows++]=iRow; 8509 } 8510 int nColumns=0; 8511 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8512 if (columnType[iColumn]>10) 8513 whichColumns[nColumns++]=iColumn; 8514 } 8515 #ifdef FUNNY_CHECK 8516 { 8517 CoinPackedMatrix rowCopy2 = *this->matrix(); 8518 rowCopy2.reverseOrdering(); 8519 int * column2 = rowCopy2.getMutableIndices(); 8520 CoinBigIndex * rowStart2 = rowCopy2.getMutableVectorStarts(); 8521 double * elementByRow2 = rowCopy2.getMutableElements(); 8522 int * rowLength2 = rowCopy2.getMutableVectorLengths(); 8523 printf("Odd rows\n"); 8524 for (int iRow=0;iRow<numberRows_;iRow++) { 8525 if (rowType[iRow]==-2) { 8526 CoinBigIndex start = rowStart[iRow]; 8527 int length=rowLength[iRow]; 8528 printf("Odd row %d -> ",iRow); 8529 for (CoinBigIndex j=start;j<start+length;j++) { 8530 int iColumn=column[j]; 8531 printf("(%d %g) ",iColumn,elementByRow[j]); 8532 } 8533 printf("Original "); 8534 start = rowStart2[iRow]; 8535 length=rowLength2[iRow]; 8536 for (CoinBigIndex j=start;j<start+length;j++) { 8537 int iColumn=column2[j]; 8538 printf("(%d %g) ",iColumn,elementByRow2[j]); 8539 } 8540 printf(">= %g/%g\n",rowLower[iRow],rowLower_[iRow]); 8541 } 8542 } 8543 printf("now columns\n"); 8544 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8545 CoinBigIndex start = columnStart[iColumn]; 8546 int length=columnLength[iColumn]; 8547 for (CoinBigIndex j=start;j<start+length;j++) { 8548 int iRow=row[j]; 8549 if (rowType[iRow]==-2) 8550 printf("Column %d row %d value %g\n", 8551 iColumn,iRow,element[j]); 8552 } 8553 } 8554 } 8555 #endif 8556 #ifdef TWOFER 8557 if (!orderedMatrix) { 8558 //lastElement; 8559 // move up 8560 CoinBigIndex put=lastElement; 8561 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8562 CoinBigIndex start = columnStart[iColumn]; 8563 columnStart[iColumn]=put-lastElement; 8564 int length=columnLength[iColumn]; 8565 for (CoinBigIndex j=start;j<start+length;j++) { 8566 row[put]=row[j]; 8567 element[put++]=element[j]; 8568 } 8569 } 8570 int numberElements = put-lastElement; 8571 columnStart[numberColumns_]=numberElements; 8572 memcpy(row,row+lastElement,numberElements*sizeof(int)); 8573 memcpy(element,element+lastElement,numberElements*sizeof(double)); 8574 } 8575 #endif 8576 // could just shrink - would be faster 8577 #if DEBUG_SOME > 0 8578 printf("%d Row types and lookup\n",nRows); 8579 int nBNew=0; 8580 int iNew=0; 8581 for (int iRow=0;iRow<numberRows_;iRow++) { 8582 int type=rowType[iRow]; 8583 char xNew='O'; 8584 if (type<=0) { 8585 xNew='N'; 8586 if (newModel->getRowStatus(iRow)==basic) { 8587 nBNew++; 8588 xNew='B'; 8589 } 8590 printf("%d -> %d type %d - new status %c\n",iRow,iNew,rowType[iRow],xNew); 8591 iNew++; 8592 } 8593 } 8594 printf("%d Column types and lookup\n",nColumns); 8595 iNew=0; 8596 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8597 int type=columnType[iColumn]; 8598 char xNew='O'; 8599 if (type<=10) { 8600 xNew='N'; 8601 if (newModel->getColumnStatus(iColumn)==basic) { 8602 nBNew++; 8603 xNew='B'; 8604 } 8605 printf("%d -> %d type %d - new status %c\n",iColumn,iNew,columnType[iColumn],xNew); 8606 iNew++; 8607 } 8608 } 8609 printf("Deleting %d rows (now %d) and %d columns (%d basic)\n",nRows,numberRows_-nRows, 8610 nColumns,nBNew); 8611 #else 8612 #if DEBUG_SOME >0 8613 printf("Deleting %d rows (now %d) and %d columns\n",nRows,numberRows_-nRows, 8614 nColumns); 8615 #endif 8616 #endif 8617 #if 0 8618 newModel->deleteRows(nRows,whichRows); 8619 newModel->deleteColumns(nColumns,whichColumns); 8620 #else 8621 newModel->deleteRowsAndColumns(nRows,whichRows,nColumns,whichColumns); 8622 #endif 8623 } else { 8624 delete newModel; 8625 newModel=NULL; 8626 *infoOut=NULL; 8627 } 8628 delete [] whichRows; 8629 delete [] infoA; 8630 delete [] startStuff; 8631 delete [] array; 8632 return newModel; 8633 } 8634 // After mini presolve 8635 void 8636 ClpSimplex::miniPostsolve(const ClpSimplex * presolvedModel, void * infoIn) 8637 { 8638 int numberTotal=numberRows_+numberColumns_; 8639 listInfo * infoX = reinterpret_cast<listInfo *>(infoIn); 8640 int nActions=infoX->numberEntries; 8641 #if DEBUG_SOME > 0 8642 #ifndef NDEBUG 8643 int numberInitial=infoX->numberInitial; 8644 #endif 8645 #endif 8646 clpPresolveInfo * infoA = infoX->start; 8647 char * startStuff = reinterpret_cast<char *>(infoA+nActions); 8648 // move status and solution across 8649 int numberColumns2=presolvedModel->numberColumns(); 8650 const double * solution2 = presolvedModel->primalColumnSolution(); 8651 unsigned char * rowStatus2 = presolvedModel->statusArray()+numberColumns2; 8652 unsigned char * columnStatus2 = presolvedModel->statusArray(); 8653 unsigned char * rowStatus = status_+numberColumns_; 8654 unsigned char * columnStatus = status_; 8655 char * rowType = new char [numberTotal]; 8656 memset(rowType,0,numberTotal); 8657 char * columnType = rowType+numberRows_; 8658 double * rowLowerX = new double [3*numberRows_+3*numberColumns_+CoinMax(numberRows_,numberColumns_)]; 8659 double * rowUpperX = rowLowerX+numberRows_; 8660 double * columnLowerX = rowUpperX+numberRows_; 8661 double * columnUpperX = columnLowerX+numberColumns_; 8662 double * objectiveX = columnUpperX+numberColumns_; 8663 double * tempElement = objectiveX+numberColumns_; 8664 double * array = tempElement+CoinMax(numberRows_,numberColumns_); 8665 memset(array,0,numberRows_*sizeof(double)); 8666 int * tempIndex = new int [CoinMax(numberRows_,numberColumns_)+4+2*numberColumns_+numberRows_]; 8667 int * forward = tempIndex+CoinMax(numberRows_,numberColumns_)+1; 8668 int * backward = forward + numberColumns_+2; 8669 int * whichRows2 = backward + numberColumns_+1; 8670 for (int i=-1;i<numberColumns_;i++) 8671 forward[i]=i+1; 8672 forward[numberColumns_]=-1; 8673 for (int i=0;i<=numberColumns_;i++) 8674 backward[i]=i-1; 8675 backward[-1]=-1; 8676 restoreInfo restore; 8677 restore.elements=tempElement; 8678 restore.indices=tempIndex; 8679 restore.startStuff=startStuff; 8680 #ifndef NDEBUG 8681 for (int i=0;i<numberRows_;i++) { 8682 rowLowerX[i]=COIN_DBL_MAX; 8683 rowUpperX[i]=-COIN_DBL_MAX; 8684 } 8685 for (int i=0;i<numberColumns_;i++) { 8686 columnLowerX[i]=COIN_DBL_MAX; 8687 columnUpperX[i]=-COIN_DBL_MAX; 8688 } 8689 #endif 8690 memset(rowActivity_,0,numberRows_*sizeof(double)); 8691 memset(dual_,0,numberRows_*sizeof(double)); 8692 // Get presolved contribution 8693 const int * row = matrix()->getIndices(); 8694 const CoinBigIndex * columnStart = matrix()->getVectorStarts(); 8695 const int * columnLength = matrix()->getVectorLengths(); 8696 const double * element = matrix()->getElements(); 8697 // forwards so dropped column at end 8698 for (int i=0;i<nActions;i++) { 8699 int type=infoA[i].type; 8700 char * getStuff = startStuff+infoA[i].infoOffset; 8701 int iRow=-1; 8702 int iColumn=-1; 8703 switch (type) { 8704 case 1: 8705 case 4: 8706 case 8: 8707 case 9: 8708 // redundant 8709 { 8710 clpPresolveInfo1_4_8 thisInfo; 8711 memcpy(&thisInfo,getStuff,sizeof(clpPresolveInfo1_4_8)); 8712 iRow = thisInfo.row; 8713 rowType[iRow]=static_cast<char>(type); 8714 } 8715 break; 8716 case 2: 8717 // sub 8718 { 8719 clpPresolveInfo2 thisInfo; 8720 memcpy(&thisInfo,getStuff,sizeof(clpPresolveInfo2)); 8721 iRow = thisInfo.row; 8722 if (iRow>=0) 8723 rowType[iRow]=2; 8724 iRow = thisInfo.row2; 8725 if (iRow>=0) 8726 rowType[iRow]=2; 8727 iColumn = thisInfo.column; 8728 columnType[iColumn]=2; 8729 } 8730 break; 8731 case 11: 8732 // movable column 8733 { 8734 clpPresolveInfo11 thisInfo; 8735 memcpy(&thisInfo,getStuff,sizeof(clpPresolveInfo11)); 8736 iColumn = thisInfo.column; 8737 columnType[iColumn]=11; 8738 } 8739 break; 8740 case 13: 8741 // empty column 8742 { 8743 clpPresolveInfo13 thisInfo; 8744 memcpy(&thisInfo,getStuff,sizeof(clpPresolveInfo13)); 8745 iColumn = thisInfo.column; 8746 columnType[iColumn]=13; 8747 } 8748 break; 8749 case 14: 8750 // doubleton 8751 { 8752 clpPresolveInfo14 thisInfo; 8753 memcpy(&thisInfo,getStuff,sizeof(clpPresolveInfo14)); 8754 iRow = thisInfo.row; 8755 iColumn = thisInfo.column; 8756 int iColumn2 = thisInfo.column2; 8757 columnType[iColumn]=3; 8758 columnType[iColumn2]=14; 8759 rowType[iRow]=3; 8760 } 8761 break; 8762 } 8763 #if DEBUG_SOME > 0 8764 printf("Action %d type %d row %d column %d\n", 8765 i,type,iRow,iColumn); 8766 #endif 8767 } 8768 int iGet=0; 8769 const double * rowLowerY = presolvedModel->rowLower(); 8770 const double * rowUpperY = presolvedModel->rowUpper(); 8771 for (int iRow=0;iRow<numberRows_;iRow++) { 8772 if (!rowType[iRow]) { 8773 rowStatus[iRow]=rowStatus2[iGet]; 8774 rowActivity_[iRow]=presolvedModel->rowActivity_[iGet]; 8775 dual_[iRow]=presolvedModel->dual_[iGet]; 8776 rowLowerX[iRow]=rowLowerY[iGet]; 8777 rowUpperX[iRow]=rowUpperY[iGet]; 8778 tempIndex[iGet]=iRow; 8779 iGet++; 8780 } else { 8781 setRowStatus(iRow,basic); 8782 } 8783 } 8784 assert (iGet==presolvedModel->numberRows()); 8785 CoinPackedMatrix matrixX; 8786 int numberElementsOriginal=matrix_->getNumElements(); 8787 const int * rowY = presolvedModel->matrix()->getIndices(); 8788 const CoinBigIndex * columnStartY = presolvedModel->matrix()->getVectorStarts(); 8789 const int * columnLengthY = presolvedModel->matrix()->getVectorLengths(); 8790 const double * elementY = presolvedModel->matrix()->getElements(); 8791 iGet=0; 8792 CoinBigIndex put=0; 8793 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8794 if (columnType[iColumn]<11) { 8795 put += CoinMax(columnLength[iColumn],columnLengthY[iGet]); 8796 iGet++; 8797 } else { 8798 put += columnLength[iColumn]; 8799 } 8800 } 8801 int lastElement=2*(put+numberColumns_)+1000; 8802 int spare = (lastElement-put)/numberColumns_; 8803 //printf("spare %d\n",spare); 8804 matrixX.reserve(numberColumns_,lastElement+2*numberElementsOriginal); 8805 int * rowX = matrixX.getMutableIndices(); 8806 CoinBigIndex * columnStartX = matrixX.getMutableVectorStarts(); 8807 int * columnLengthX = matrixX.getMutableVectorLengths(); 8808 double * elementX = matrixX.getMutableElements(); 8809 const double * columnLowerY = presolvedModel->columnLower(); 8810 const double * columnUpperY = presolvedModel->columnUpper(); 8811 iGet=0; 8812 put=0; 8813 memcpy(objectiveX,this->objective(),numberColumns_*sizeof(double)); 8814 const double * objectiveY=presolvedModel->objective(); 8815 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8816 columnStartX[iColumn]=put; 8817 if (columnType[iColumn]<11) { 8818 CoinBigIndex next=put+columnLengthY[iGet]; 8819 if (spare>=0) 8820 next += CoinMax(columnLength[iColumn]-columnLengthY[iGet],0)+spare; 8821 columnStatus[iColumn]=columnStatus2[iGet]; 8822 columnActivity_[iColumn]=solution2[iGet]; 8823 columnLowerX[iColumn]=columnLowerY[iGet]; 8824 columnUpperX[iColumn]=columnUpperY[iGet]; 8825 columnLengthX[iColumn]=columnLengthY[iGet]; 8826 objectiveX[iColumn]=objectiveY[iGet]; 8827 for (CoinBigIndex j=columnStartY[iGet];j<columnStartY[iGet]+columnLengthY[iGet];j++) { 8828 rowX[put]=tempIndex[rowY[j]]; 8829 elementX[put++]=elementY[j]; 8830 } 8831 columnActivity_[iColumn]=presolvedModel->columnActivity_[iGet]; 8832 iGet++; 8833 columnType[iColumn]=0; 8834 put = next; 8835 } else { 8836 put += CoinMax(columnLength[iColumn]+spare,0); 8837 columnActivity_[iColumn]=0.0; 8838 columnType[iColumn]=1; 8839 columnLengthX[iColumn]=0; 8840 setColumnStatus(iColumn,superBasic); 8841 } 8842 } 8843 assert (put<=lastElement); 8844 columnStartX[numberColumns_]=lastElement+numberElementsOriginal; 8845 assert (put<=lastElement); 8846 assert (iGet==numberColumns2); 8847 matrixX.times(columnActivity_,rowActivity_); 8848 if (optimizationDirection_<0) { 8849 for (int i=0;i<numberColumns_;i++) 8850 objectiveX[i]=-objectiveX[i]; 8851 } 8852 #if 0 8853 // get row copy 8854 CoinPackedMatrix rowCopy = *matrix(); 8855 rowCopy.reverseOrdering(); 8856 const int * column = rowCopy.getIndices(); 8857 const CoinBigIndex * rowStart = rowCopy.getVectorStarts(); 8858 const int * rowLength = rowCopy.getVectorLengths(); 8859 const double * elementByRow = rowCopy.getElements(); 8860 #endif 8861 #if 1 //DEBUG_SOME > 0 8862 double * tempRhs=new double[numberRows_]; 8863 #endif 8864 #ifndef NDEBUG 8865 bool checkMatrixAtEnd=true; 8866 #endif 8867 for (int i=nActions-1;i>=0;i--) { 8868 #if DEBUG_SOME > 0 8869 #if 1 //ndef NDEBUG 8870 memset(tempRhs,0,numberRows_*sizeof(double)); 8871 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8872 //if (columnActivity_[iColumn]<columnLower_[iColumn]-1.0e-5|| 8873 // columnActivity_[iColumn]>columnUpper_[iColumn]+1.0e-5) 8874 //printf ("Bad column %d %g <= %g <= %g\n",iColumn, 8875 // columnLower_[iColumn],columnActivity_[iColumn],columnUpper_[iColumn]); 8876 double value=columnActivity_[iColumn]; 8877 for (CoinBigIndex j=columnStartX[iColumn]; 8878 j<columnStartX[iColumn]+columnLengthX[iColumn];j++) { 8879 int jRow=rowX[j]; 8880 tempRhs[jRow] += value*elementX[j]; 8881 } 8882 } 8883 int nBadRow=0; 8884 for (int iRow=0;iRow<numberRows_;iRow++) { 8885 if (rowLowerX[iRow]==COIN_DBL_MAX) 8886 continue; 8887 if (fabs(tempRhs[iRow]-rowActivity_[iRow])>1.0e-4) 8888 printf("Row %d temprhs %g rowActivity_ %g\n",iRow,tempRhs[iRow], 8889 rowActivity_[iRow]); 8890 if (tempRhs[iRow]<rowLowerX[iRow]-1.0e-4|| 8891 tempRhs[iRow]>rowUpperX[iRow]+1.0e-4) { 8892 printf("Row %d %g %g %g\n",iRow,rowLowerX[iRow], 8893 tempRhs[iRow],rowUpperX[iRow]); 8894 nBadRow++; 8895 } 8896 } 8897 assert (!nBadRow); 8898 #endif 8899 #ifndef NDEBUG 8900 if (i>=numberInitial&&true) { 8901 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 8902 if (columnLengthX[iColumn]) { 8903 double djValue = objectiveX[iColumn]; 8904 for (CoinBigIndex j=columnStartX[iColumn]; 8905 j<columnStartX[iColumn]+columnLengthX[iColumn];j++) { 8906 int jRow=rowX[j]; 8907 djValue -= dual_[jRow]*elementX[j]; 8908 } 8909 char inOut=columnType[iColumn] ? 'O' : 'I'; 8910 if (djValue>10.0*dualTolerance_&& 8911 columnActivity_[iColumn]>columnLowerX[iColumn]+primalTolerance_) 8912 printf("bad dj on column %d lb dj %g in/out %c\n",iColumn,djValue,inOut); 8913 else if (djValue<-10.0*dualTolerance_&& 8914 columnActivity_[iColumn]<columnUpperX[iColumn]-primalTolerance_) 8915 printf("bad dj on column %d ub dj %g in/out %c\n",iColumn,djValue,inOut); 8916 } 8917 } 8918 } 8919 #endif 8920 #endif 8921 int type=infoA[i].type; 8922 int iRow=-1; 8923 int iColumn=-1; 8924 #if DEBUG_SOME>0 8925 printf("Action %d type %d\n",i,type); 8926 #ifndef NDEBUG 8927 checkBasis(this, rowType, columnType); 8928 #endif 8929 #endif 8930 switch (type) { 8931 case 1: 8932 case 4: 8933 case 8: 8934 case 9: 8935 // redundant 8936 { 8937 #ifndef NDEBUG 8938 if (type==4) 8939 checkMatrixAtEnd=false; 8940 #endif 8941 clpPresolveInfo8 thisInfo; 8942 copyFromSave(restore,infoA[i],&thisInfo); 8943 iRow = thisInfo.row; 8944 // Insert row and modify RHS 8945 #ifdef CLP_USER_DRIVEN 8946 if (type>=8) { 8947 thisInfo.tempElement=tempElement; 8948 thisInfo.tempIndex=tempIndex; 8949 thisInfo.rowLowerX=rowLowerX; 8950 thisInfo.rowUpperX=rowUpperX; 8951 eventHandler_->eventWithInfo(ClpEventHandler::modifyMatrixInMiniPostsolve,&thisInfo); 8952 } else { 8953 #endif 8954 if (rowLower_[iRow]>-1.0e30) 8955 rowLowerX[iRow]=rowLower_[iRow]-thisInfo.oldRowLower; 8956 else 8957 rowLowerX[iRow]=thisInfo.oldRowLower; 8958 if (rowUpper_[iRow]<1.0e30) 8959 rowUpperX[iRow]=rowUpper_[iRow]-thisInfo.oldRowUpper; 8960 else 8961 rowUpperX[iRow]=thisInfo.oldRowUpper; 8962 #ifdef CLP_USER_DRIVEN 8963 } 8964 #endif 8965 int n=thisInfo.lengthRow; 8966 double rhs=0.0; 8967 for (int i=0;i<n;i++) { 8968 int iColumn = tempIndex[i]; 8969 double value = tempElement[i]; 8970 rhs += value*columnActivity_[iColumn]; 8971 int length=columnLengthX[iColumn]; 8972 CoinBigIndex start=columnStartX[iColumn]; 8973 int nextColumn=forward[iColumn]; 8974 CoinBigIndex startNext=columnStartX[nextColumn]; 8975 if (start+length==startNext) { 8976 // need more 8977 moveAround(numberColumns_,numberElementsOriginal, 8978 iColumn,length+1, 8979 forward,backward,columnStartX,columnLengthX, 8980 rowX,elementX); 8981 start=columnStartX[iColumn]; 8982 } 8983 columnLengthX[iColumn]++; 8984 rowX[start+length]=iRow; 8985 elementX[start+length]=value; 8986 } 8987 rowActivity_[iRow]=rhs; 8988 rowType[iRow]=0; 8989 assert (getRowStatus(iRow)==basic); 8990 } 8991 break; 8992 case 2: 8993 // sub 8994 { 8995 clpPresolveInfo2 thisInfo; 8996 copyFromSave(restore,infoA[i],&thisInfo); 8997 iColumn = thisInfo.column; 8998 columnType[iColumn]=0; 8999 int jRowLower=thisInfo.row; 9000 int jRowUpper=thisInfo.row2; 9001 int length=columnLengthX[iColumn]; 9002 CoinBigIndex start=columnStartX[iColumn]; 9003 int nextColumn=forward[iColumn]; 9004 CoinBigIndex startNext=columnStartX[nextColumn]; 9005 if (start+length==startNext) { 9006 // need more 9007 moveAround(numberColumns_,numberElementsOriginal, 9008 iColumn,length+(jRowLower<0||jRowUpper<0) ? 1 : 2, 9009 forward,backward,columnStartX,columnLengthX, 9010 rowX,elementX); 9011 start=columnStartX[iColumn]; 9012 } 9013 // modify 9014 double valueLower=thisInfo.coefficient; 9015 double valueUpper=thisInfo.coefficient2; 9016 if (jRowLower>=0) { 9017 rowType[jRowLower]=0; 9018 if (rowLower_[jRowLower]>-1.0e30) 9019 rowLowerX[jRowLower]=rowLower_[jRowLower]-thisInfo.oldRowLower; 9020 else 9021 rowLowerX[jRowLower]=thisInfo.oldRowLower; 9022 if (rowUpper_[jRowLower]<1.0e30) 9023 rowUpperX[jRowLower]=rowUpper_[jRowLower]-thisInfo.oldRowUpper; 9024 else 9025 rowUpperX[jRowLower]=thisInfo.oldRowUpper; 9026 columnLengthX[iColumn]++; 9027 rowX[start+length]=jRowLower; 9028 elementX[start+length]=valueLower; 9029 rowActivity_[jRowLower]=valueLower*columnActivity_[iColumn]; 9030 // start off with slack basic 9031 setRowStatus(jRowLower,basic); 9032 length++; 9033 } 9034 if (jRowUpper>=0) { 9035 rowType[jRowUpper]=0; 9036 if (rowLower_[jRowUpper]>-1.0e30) 9037 rowLowerX[jRowUpper]=rowLower_[jRowUpper]-thisInfo.oldRowLower2; 9038 else 9039 rowLowerX[jRowUpper]=thisInfo.oldRowLower2; 9040 if (rowUpper_[jRowUpper]<1.0e30) 9041 rowUpperX[jRowUpper]=rowUpper_[jRowUpper]-thisInfo.oldRowUpper2; 9042 else 9043 rowUpperX[jRowUpper]=thisInfo.oldRowUpper2; 9044 columnLengthX[iColumn]++; 9045 rowX[start+length]=jRowUpper; 9046 elementX[start+length]=valueUpper; 9047 rowActivity_[jRowUpper]=valueUpper*columnActivity_[iColumn]; 9048 // start off with slack basic 9049 setRowStatus(jRowUpper,basic); 9050 length++; 9051 } 9052 double djValue = objectiveX[iColumn]; 9053 for (CoinBigIndex j=columnStartX[iColumn]; 9054 j<columnStartX[iColumn]+columnLengthX[iColumn];j++) { 9055 int jRow=rowX[j]; 9056 djValue -= dual_[jRow]*elementX[j]; 9057 } 9058 if (thisInfo.oldColumnLower==thisInfo.oldColumnUpper) 9059 djValue=0.0; 9060 if (getColumnStatus(iColumn)!=basic) { 9061 setColumnStatus(iColumn,superBasic); 9062 double solValue = columnActivity_[iColumn]; 9063 bool hasToBeBasic=false; 9064 if (getColumnStatus(iColumn)!=isFree) { 9065 if (solValue>thisInfo.oldColumnLower+primalTolerance_&& 9066 solValue<thisInfo.oldColumnUpper-primalTolerance_) 9067 hasToBeBasic=true; 9068 } else { 9069 hasToBeBasic=fabs(djValue)>dualTolerance_; 9070 } 9071 if (solValue<=thisInfo.oldColumnLower+primalTolerance_&&djValue<-dualTolerance_) { 9072 #if DEBUG_SOME > 1 9073 printf("odd column %d at lb of %g has dj of %g\n", 9074 iColumn,thisInfo.oldColumnLower,djValue); 9075 #endif 9076 hasToBeBasic=true; 9077 } else if (solValue>=thisInfo.oldColumnUpper-primalTolerance_&&djValue>dualTolerance_) { 9078 #if DEBUG_SOME > 1 9079 printf("odd column %d at ub of %g has dj of %g\n", 9080 iColumn,thisInfo.oldColumnUpper,djValue); 9081 #endif 9082 hasToBeBasic=true; 9083 } 9084 if (hasToBeBasic) { 9085 if (jRowLower>=0&&jRowUpper>=0) { 9086 // choose 9087 #if DEBUG_SOME > 1 9088 printf("lower row %d %g <= %g <= %g\n", 9089 jRowLower,rowLowerX[jRowLower],rowActivity_[jRowLower], 9090 rowUpperX[jRowLower]); 9091 #endif 9092 double awayLower = CoinMin(rowActivity_[jRowLower]-rowLowerX[jRowLower], 9093 rowUpperX[jRowLower]-rowActivity_[jRowLower]); 9094 #if DEBUG_SOME > 1 9095 printf("upper row %d %g <= %g <= %g\n", 9096 jRowUpper,rowLowerX[jRowUpper],rowActivity_[jRowUpper], 9097 rowUpperX[jRowUpper]); 9098 #endif 9099 double awayUpper = CoinMin(rowActivity_[jRowUpper]-rowLowerX[jRowUpper], 9100 rowUpperX[jRowUpper]-rowActivity_[jRowUpper]); 9101 if (awayLower>awayUpper) 9102 jRowLower=-1; 9103 } 9104 if (jRowLower<0) { 9105 setRowStatus(jRowUpper,superBasic); 9106 setColumnStatus(iColumn,basic); 9107 dual_[jRowUpper]=djValue/valueUpper; 9108 } else { 9109 setRowStatus(jRowLower,superBasic); 9110 setColumnStatus(iColumn,basic); 9111 dual_[jRowLower]=djValue/valueLower; 9112 } 9113 } 9114 } 9115 columnLowerX[iColumn]=thisInfo.oldColumnLower; 9116 columnUpperX[iColumn]=thisInfo.oldColumnUpper; 9117 } 9118 break; 9119 case 11: 9120 // movable column 9121 { 9122 clpPresolveInfo11 thisInfo; 9123 copyFromSave(restore,infoA[i],&thisInfo); 9124 iColumn = thisInfo.column; 9125 columnType[iColumn]=0; 9126 assert (!columnLengthX[iColumn]); 9127 int newLength=thisInfo.lengthColumn; 9128 CoinBigIndex start=columnStartX[iColumn]; 9129 int nextColumn=forward[iColumn]; 9130 CoinBigIndex startNext=columnStartX[nextColumn]; 9131 if (start+newLength>startNext) { 9132 // need more 9133 moveAround(numberColumns_,numberElementsOriginal, 9134 iColumn,newLength, 9135 forward,backward,columnStartX,columnLengthX, 9136 rowX,elementX); 9137 start=columnStartX[iColumn]; 9138 } 9139 columnLengthX[iColumn]=newLength; 9140 memcpy(rowX+start,tempIndex,newLength*sizeof(int)); 9141 memcpy(elementX+start,tempElement,newLength*sizeof(double)); 9142 double solValue=thisInfo.fixedTo; 9143 columnActivity_[iColumn]=solValue; 9144 if (solValue) { 9145 for (int j=columnStartX[iColumn]; 9146 j<columnStartX[iColumn]+columnLengthX[iColumn];j++) { 9147 int jRow=rowX[j]; 9148 double value=elementX[j]*solValue; 9149 rowActivity_[jRow] += value; 9150 if (rowLowerX[jRow]>-1.0e30) 9151 rowLowerX[jRow] += value; 9152 if (rowUpperX[jRow]<1.0e30) 9153 rowUpperX[jRow] += value; 9154 } 9155 } 9156 double djValue = objectiveX[iColumn]; 9157 for (CoinBigIndex j=columnStartX[iColumn]; 9158 j<columnStartX[iColumn]+columnLengthX[iColumn];j++) { 9159 int jRow=rowX[j]; 9160 djValue -= dual_[jRow]*elementX[j]; 9161 } 9162 if (thisInfo.oldColumnLower==thisInfo.oldColumnUpper) 9163 djValue=0.0; 9164 if (getColumnStatus(iColumn)!=basic) { 9165 setColumnStatus(iColumn,superBasic); 9166 bool hasToBeBasic=false; 9167 if (getColumnStatus(iColumn)!=isFree) { 9168 if (solValue>thisInfo.oldColumnLower+primalTolerance_&& 9169 solValue<thisInfo.oldColumnUpper-primalTolerance_) 9170 hasToBeBasic=true; 9171 } else { 9172 hasToBeBasic=fabs(djValue)>dualTolerance_; 9173 } 9174 if (solValue<=thisInfo.oldColumnLower+primalTolerance_&&djValue<-dualTolerance_) { 9175 #if DEBUG_SOME > 1 9176 printf("odd column %d at lb of %g has dj of %g\n", 9177 iColumn,thisInfo.oldColumnLower,djValue); 9178 #endif 9179 hasToBeBasic=true; 9180 } else if (solValue>=thisInfo.oldColumnUpper-primalTolerance_&&djValue>dualTolerance_) { 9181 #if DEBUG_SOME > 1 9182 printf("odd column %d at ub of %g has dj of %g\n", 9183 iColumn,thisInfo.oldColumnUpper,djValue); 9184 #endif 9185 hasToBeBasic=true; 9186 } 9187 if (hasToBeBasic) { 9188 abort(); 9189 //setRowStatus(iRow,superBasic); 9190 setColumnStatus(iColumn,basic); 9191 //dual_[iRow]=djValue/value; 9192 } 9193 } 9194 columnLowerX[iColumn]=thisInfo.oldColumnLower; 9195 columnUpperX[iColumn]=thisInfo.oldColumnUpper; 9196 } 9197 break; 9198 case 13: 9199 // empty (or initially fixed) column 9200 { 9201 clpPresolveInfo13 thisInfo; 9202 copyFromSave(restore,infoA[i],&thisInfo); 9203 iColumn = thisInfo.column; 9204 columnType[iColumn]=0; 9205 columnLowerX[iColumn]=thisInfo.oldColumnLower; 9206 columnUpperX[iColumn]=thisInfo.oldColumnUpper; 9207 assert (!columnLengthX[iColumn]); 9208 int newLength=columnLength[iColumn]; 9209 CoinBigIndex startOriginal = columnStart[iColumn]; 9210 double solValue=columnLower_[iColumn]; 9211 columnActivity_[iColumn]=solValue; 9212 if (newLength&&solValue) { 9213 for (int j=startOriginal;j<startOriginal+newLength;j++) { 9214 int iRow=row[j]; 9215 double value = element[j]*solValue; 9216 rowActivity_[iRow] += value; 9217 double lower = rowLowerX[iRow]; 9218 if (lower>-1.0e20) 9219 rowLowerX[iRow]=lower+value; 9220 double upper = rowUpperX[iRow]; 9221 if (upper<1.0e20) 9222 rowUpperX[iRow]=upper+value; 9223 } 9224 } 9225 #ifndef NDEBUG 9226 // copy from original 9227 CoinBigIndex start=columnStartX[iColumn]; 9228 int nextColumn=forward[iColumn]; 9229 CoinBigIndex startNext=columnStartX[nextColumn]; 9230 if (start+newLength>startNext) { 9231 // need more 9232 moveAround(numberColumns_,numberElementsOriginal, 9233 iColumn,newLength, 9234 forward,backward,columnStartX,columnLengthX, 9235 rowX,elementX); 9236 start=columnStartX[iColumn]; 9237 } 9238 columnLengthX[iColumn]=newLength; 9239 memcpy(rowX+start,row+startOriginal,newLength*sizeof(int)); 9240 memcpy(elementX+start,element+startOriginal,newLength*sizeof(double)); 9241 #endif 9242 } 9243 break; 9244 case 14: 9245 // doubleton 9246 { 9247 clpPresolveInfo14 thisInfo; 9248 copyFromSave(restore,infoA[i],&thisInfo); 9249 iRow = thisInfo.row; 9250 rowType[iRow]=0; 9251 int iColumn1 = thisInfo.column; 9252 int iColumn2 = thisInfo.column2; 9253 columnType[iColumn2]=0; 9254 double value1=thisInfo.value1; 9255 int length=columnLengthX[iColumn1]; 9256 CoinBigIndex start=columnStartX[iColumn1]; 9257 for (int i=0;i<length;i++) { 9258 int jRow=rowX[i+start]; 9259 array[jRow]=elementX[i+start]; 9260 whichRows2[i]=jRow; 9261 } 9262 int n=length; 9263 length=columnLengthX[iColumn2]; 9264 assert (!length); 9265 int newLength=thisInfo.lengthColumn2; 9266 start=columnStartX[iColumn2]; 9267 int nextColumn=forward[iColumn2]; 9268 CoinBigIndex startNext=columnStartX[nextColumn]; 9269 if (start+newLength>startNext) { 9270 // need more 9271 moveAround(numberColumns_,numberElementsOriginal, 9272 iColumn2,newLength, 9273 forward,backward,columnStartX,columnLengthX, 9274 rowX,elementX); 9275 start=columnStartX[iColumn2]; 9276 } 9277 columnLengthX[iColumn2]=newLength; 9278 memcpy(rowX+start,tempIndex,newLength*sizeof(int)); 9279 memcpy(elementX+start,tempElement,newLength*sizeof(double)); 9280 double value2=0.0; 9281 for (int i=0;i<newLength;i++) { 9282 int jRow=tempIndex[i]; 9283 double value=tempElement[i]; 9284 if (jRow==iRow) { 9285 value2=value; 9286 break; 9287 } 9288 } 9289 assert (value2); 9290 double rhs=thisInfo.rhs; 9291 double multiplier1 = rhs/value2; 9292 double multiplier = value1/value2; 9293 for (int i=0;i<newLength;i++) { 9294 int jRow=tempIndex[i]; 9295 double value = array[jRow]; 9296 double valueNew = value + multiplier*tempElement[i]; 9297 double rhsMod = multiplier1*tempElement[i]; 9298 if (rowLowerX[jRow]>-1.0e30) 9299 rowLowerX[jRow] += rhsMod; 9300 if (rowUpperX[jRow]<1.0e30) 9301 rowUpperX[jRow] += rhsMod; 9302 rowActivity_[jRow] += rhsMod; 9303 if (!value) { 9304 array[jRow]=valueNew; 9305 whichRows2[n++]=jRow; 9306 } else { 9307 if (!valueNew) 9308 valueNew=1.0e-100; 9309 array[jRow]=valueNew; 9310 } 9311 } 9312 length=0; 9313 for (int i=0;i<n;i++) { 9314 int jRow = whichRows2[i]; 9315 double value = array[jRow]; 9316 array[jRow]=0.0; 9317 if (fabs(value)<1.0e-13) 9318 value=0.0; 9319 if (value) { 9320 tempIndex[length] = jRow; 9321 tempElement[length++]=value; 9322 } 9323 } 9324 start=columnStartX[iColumn1]; 9325 nextColumn=forward[iColumn1]; 9326 startNext=columnStartX[nextColumn]; 9327 if (start+length>startNext) { 9328 // need more 9329 moveAround(numberColumns_,numberElementsOriginal, 9330 iColumn1,length, 9331 forward,backward,columnStartX,columnLengthX, 9332 rowX,elementX); 9333 start=columnStartX[iColumn1]; 9334 } 9335 columnLengthX[iColumn1]=length; 9336 memcpy(rowX+start,tempIndex,length*sizeof(int)); 9337 memcpy(elementX+start,tempElement,length*sizeof(double)); 9338 objectiveX[iColumn2]=thisInfo.oldObjective2; 9339 objectiveX[iColumn1] += objectiveX[iColumn2]*multiplier; 9340 //offset += rhs*(objectiveX[iColumn2]*multiplier); 9341 double value = multiplier1-columnActivity_[iColumn1]*multiplier; 9342 #if DEBUG_SOME > 0 9343 printf("type14 column2 %d rhs %g value1 %g - value %g\n", 9344 iColumn2,thisInfo.rhs,thisInfo.value1,value); 9345 #endif 9346 columnActivity_[iColumn2]=value; 9347 double djValue1 = objectiveX[iColumn1]; 9348 for (CoinBigIndex j=columnStartX[iColumn1]; 9349 j<columnStartX[iColumn1]+columnLengthX[iColumn1];j++) { 9350 int jRow=rowX[j]; 9351 djValue1 -= dual_[jRow]*elementX[j]; 9352 } 9353 bool fixed = (thisInfo.oldColumnLower==thisInfo.oldColumnUpper); 9354 columnLowerX[iColumn1]=thisInfo.oldColumnLower; 9355 columnUpperX[iColumn1]=thisInfo.oldColumnUpper; 9356 double djValue2 = objectiveX[iColumn2]; 9357 for (CoinBigIndex j=columnStartX[iColumn2]; 9358 j<columnStartX[iColumn2]+columnLengthX[iColumn2];j++) { 9359 int jRow=rowX[j]; 9360 djValue2 -= dual_[jRow]*elementX[j]; 9361 } 9362 bool fixed2 = (thisInfo.oldColumnLower2==thisInfo.oldColumnUpper2); 9363 columnLowerX[iColumn2]=thisInfo.oldColumnLower2; 9364 columnUpperX[iColumn2]=thisInfo.oldColumnUpper2; 9365 if (getColumnStatus(iColumn1)!=basic) { 9366 setColumnStatus(iColumn1,superBasic); 9367 double solValue = columnActivity_[iColumn1]; 9368 double lowerValue = columnLowerX[iColumn1]; 9369 double upperValue = columnUpperX[iColumn1]; 9370 setColumnStatus(iColumn2,superBasic); 9371 double solValue2 = columnActivity_[iColumn2]; 9372 double lowerValue2 = columnLowerX[iColumn2]; 9373 double upperValue2 = columnUpperX[iColumn2]; 9374 int choice=-1; 9375 double best=COIN_DBL_MAX; 9376 for (int iTry=0;iTry<3;iTry++) { 9377 double pi=0.0; 9378 switch (iTry) { 9379 // leave 9380 case 0: 9381 break; 9382 // iColumn1 basic 9383 case 1: 9384 pi=djValue1/value1; 9385 break; 9386 // iColumn2 basic 9387 case 2: 9388 pi=djValue2/value2; 9389 break; 9390 } 9391 // djs 9392 double dj1 = djValue1-value1*pi;; 9393 double bad1=0.0; 9394 if (!fixed) { 9395 if (dj1<-dualTolerance_&&solValue<=lowerValue+primalTolerance_) 9396 bad1=-dj1; 9397 else if (dj1>dualTolerance_&&solValue>=upperValue-primalTolerance_) 9398 bad1=dj1; 9399 } 9400 double dj2 = djValue2-value2*pi; 9401 double bad2=0.0; 9402 if (!fixed2) { 9403 if (dj2<-dualTolerance_&&solValue2<=lowerValue2+primalTolerance_) 9404 bad2=-dj2; 9405 else if (dj2>dualTolerance_&&solValue2>=upperValue2-primalTolerance_) 9406 bad2=dj2; 9407 } 9408 if (CoinMax(bad1,bad2)<best) { 9409 best=CoinMax(bad1,bad2); 9410 choice=iTry; 9411 } 9412 } 9413 // but values override 9414 if (solValue>lowerValue+primalTolerance_&& 9415 solValue<upperValue-primalTolerance_) { 9416 if (getColumnStatus(iColumn1)!=isFree|| 9417 fabs(djValue1)>dualTolerance_) 9418 choice=1; 9419 } else if (solValue2>lowerValue2+primalTolerance_&& 9420 solValue2<upperValue2-primalTolerance_) { 9421 if (getColumnStatus(iColumn2)!=isFree|| 9422 fabs(djValue2)>dualTolerance_) 9423 choice=2; 9424 } 9425 if (choice==1) { 9426 // iColumn1 in basis 9427 setRowStatus(iRow,superBasic); 9428 setColumnStatus(iColumn1,basic); 9429 dual_[iRow]=djValue1/value1; 9430 } else if (choice==2) { 9431 // iColumn2 in basis 9432 setRowStatus(iRow,superBasic); 9433 setColumnStatus(iColumn2,basic); 9434 dual_[iRow]=djValue2/value2; 9435 } 9436 } else { 9437 if (fabs(djValue1)>10.0*dualTolerance_) { 9438 // iColumn2 in basis 9439 setRowStatus(iRow,superBasic); 9440 setColumnStatus(iColumn2,basic); 9441 dual_[iRow]=djValue2/value2; 9442 } else { 9443 // do we need iColumn2 in basis 9444 double solValue2 = columnActivity_[iColumn2]; 9445 double lowerValue2 = columnLowerX[iColumn2]; 9446 double upperValue2 = columnUpperX[iColumn2]; 9447 if (solValue2>lowerValue2+primalTolerance_&& 9448 solValue2<upperValue2-primalTolerance_&& 9449 getColumnStatus(iColumn2)!=isFree) { 9450 // iColumn2 in basis 9451 setRowStatus(iRow,superBasic); 9452 setColumnStatus(iColumn2,basic); 9453 dual_[iRow]=djValue2/value2; 9454 } 9455 } 9456 } 9457 rowLowerX[iRow]=rhs; 9458 rowUpperX[iRow]=rhs; 9459 //abort(); 9460 } 9461 break; 9462 } 9463 } 9464 #ifndef NDEBUG 9465 #ifndef CLP_USER_DRIVEN 9466 // otherwise user may be fudging rhs 9467 for (int iRow=0;iRow<numberRows_;iRow++) { 9468 #ifndef CLP_USER_DRIVEN 9469 assert (fabs(rowLower_[iRow]-rowLowerX[iRow])<1.0e-4); 9470 assert (fabs(rowUpper_[iRow]-rowUpperX[iRow])<1.0e-4); 9471 #else 9472 if (fabs(rowLower_[iRow]-rowLowerX[iRow])>1.0e-4|| 9473 fabs(rowUpper_[iRow]-rowUpperX[iRow])>1.0e-4) 9474 printf("USER row %d lower,x %g %g upper,x %g %g\n",iRow,rowLower_[iRow],rowLowerX[iRow], 9475 rowUpper_[iRow],rowUpperX[iRow]); 9476 9477 #endif 9478 } 9479 double * objective = this->objective(); 9480 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 9481 assert (columnLower_[iColumn]==columnLowerX[iColumn]); 9482 assert (columnUpper_[iColumn]==columnUpperX[iColumn]); 9483 assert (fabs(objective[iColumn]-objectiveX[iColumn]*optimizationDirection_)<1.0e-3); 9484 } 9485 #endif 9486 if (checkMatrixAtEnd) { 9487 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 9488 assert (columnLength[iColumn]==columnLengthX[iColumn]); 9489 int length=columnLengthX[iColumn]; 9490 CoinBigIndex startX=columnStartX[iColumn]; 9491 CoinSort_2(rowX+startX,rowX+startX+length,elementX+startX); 9492 CoinBigIndex start=columnStart[iColumn]; 9493 memcpy(tempIndex,row+start,length*sizeof(int)); 9494 memcpy(tempElement,element+start,length*sizeof(double)); 9495 CoinSort_2(tempIndex,tempIndex+length,tempElement); 9496 for (int i=0;i<length;i++) { 9497 assert (rowX[i+startX]==tempIndex[i]); 9498 assert (fabs(elementX[i+startX]-tempElement[i])<1.0e-5); 9499 } 9500 } 9501 } 9502 #endif 9503 delete [] rowType; 9504 delete [] tempIndex; 9505 // use temp matrix and do this every action 9506 #if DEBUG_SOME>0 9507 matrix()->times(columnActivity_,rowActivity_); 9508 for (int i=0;i<numberColumns_;i++) 9509 assert (columnActivity_[i]>=columnLower_[i]-1.0e-5&& 9510 columnActivity_[i]<=columnUpper_[i]+1.0e-5); 9511 int nBad=0; 9512 for (int i=0;i<numberRows_;i++) { 9513 if (rowActivity_[i]<rowLower_[i]-1.0e-5|| 9514 rowActivity_[i]>rowUpper_[i]+1.0e-5) { 9515 printf("Row %d %g %g %g\n",i,rowLower_[i], 9516 rowActivity_[i],rowUpper_[i]); 9517 nBad++; 9518 } 9519 } 9520 ClpTraceDebug (!nBad); 9521 #endif 9522 #if 1 //DEBUG_SOME > 0 9523 delete [] tempRhs; 9524 #endif 9525 delete infoX; 9526 delete [] infoA; 9527 delete [] rowLowerX; 9528 } -
trunk/Clp/src/ClpSimplexOther.hpp
r1785 r1825 101 101 const double * changeLowerBound, const double * changeUpperBound, 102 102 const double * changeLowerRhs, const double * changeUpperRhs); 103 int parametricsObj(double startingTheta, double & endingTheta, 104 const double * changeObjective); 103 105 /// Finds best possible pivot 104 106 double bestPivot(bool justColumns=false); … … 129 131 int parametricsLoop(parametricsData & paramData, 130 132 ClpDataSave & data,bool canSkipFactorization=false); 133 int parametricsObjLoop(parametricsData & paramData, 134 ClpDataSave & data,bool canSkipFactorization=false); 131 135 /** Refactorizes if necessary 132 136 Checks if finished. Updates status. … … 137 141 */ 138 142 void statusOfProblemInParametrics(int type, ClpDataSave & saveData); 143 void statusOfProblemInParametricsObj(int type, ClpDataSave & saveData); 139 144 /** This has the flow between re-factorizations 140 145 … … 155 160 int nextTheta(int type, double maxTheta, parametricsData & paramData, 156 161 const double * changeObjective); 162 int whileIteratingObj(parametricsData & paramData); 163 int nextThetaObj(double maxTheta, parametricsData & paramData); 157 164 /// Restores bound to original bound 158 165 void originalBound(int iSequence, double theta, const double * changeLower,
Note: See TracChangeset
for help on using the changeset viewer.