- Timestamp:
- Jan 20, 2011 12:00:28 PM (9 years ago)
- Location:
- trunk/Clp/src
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Clp/src/ClpDynamicExampleMatrix.cpp
r1665 r1676 655 655 } 656 656 } 657 ClpDynamicMatrix::createVariable(model, bestSequence );657 ClpDynamicMatrix::createVariable(model, bestSequence/*, bestSequence2*/); 658 658 // clear for next iteration 659 659 savedBestSequence_ = -1; -
trunk/Clp/src/ClpDynamicMatrix.cpp
r1665 r1676 94 94 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_); 95 95 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_); 96 status_ = ClpCopyOfArray(rhs.status_, numberSets_);96 status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_); 97 97 model_ = rhs.model_; 98 98 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; … … 110 110 maximumGubColumns_ = rhs.maximumGubColumns_; 111 111 maximumElements_ = rhs.maximumElements_; 112 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_ );112 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1); 113 113 next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_); 114 114 startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1); … … 119 119 columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_); 120 120 columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_); 121 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, maximumGubColumns_);121 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_); 122 122 } 123 123 … … 143 143 else 144 144 maximumElements_ = 0; 145 startSet_ = new int [numberSets_ ];145 startSet_ = new int [numberSets_+1]; 146 146 next_ = new int [maximumGubColumns_]; 147 147 // fill in startSet and next … … 157 157 } 158 158 } 159 startSet_[numberSets_] = starts[numberSets_]; 159 160 int numberColumns = model->numberColumns(); 160 161 int numberRows = model->numberRows(); … … 165 166 int frequency = model->factorizationFrequency(); 166 167 int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4; 168 // But we may have two per row + one for incoming (make it two) 169 numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall); 167 170 // for small problems this could be too big 168 171 //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_); … … 256 259 keyVariable_ = new int[numberSets_]; 257 260 if (status) { 258 status_ = ClpCopyOfArray(status, numberSets_);261 status_ = ClpCopyOfArray(status, 2*numberSets_); 259 262 assert (dynamicStatus); 260 dynamicStatus_ = ClpCopyOfArray(dynamicStatus, numberGubColumns_);263 dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_); 261 264 } else { 262 status_ = new unsigned char [ numberSets_];265 status_ = new unsigned char [2*numberSets_]; 263 266 memset(status_, 0, numberSets_); 264 267 int i; … … 267 270 setStatus(i, ClpSimplex::basic); 268 271 } 269 dynamicStatus_ = new unsigned char [ numberGubColumns_];272 dynamicStatus_ = new unsigned char [2*numberGubColumns_]; 270 273 memset(dynamicStatus_, 0, numberGubColumns_); // for clarity 271 274 for (i = 0; i < numberGubColumns_; i++) … … 353 356 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_); 354 357 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_); 355 status_ = ClpCopyOfArray(rhs.status_, numberSets_);358 status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_); 356 359 model_ = rhs.model_; 357 360 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; … … 369 372 maximumGubColumns_ = rhs.maximumGubColumns_; 370 373 maximumElements_ = rhs.maximumElements_; 371 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_ );374 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1); 372 375 next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_); 373 376 startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1); … … 378 381 columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_); 379 382 columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_); 380 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, maximumGubColumns_);383 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_); 381 384 } 382 385 return *this; … … 533 536 // check flagged variable and correct dj 534 537 if (!flagged(iSequence)) { 538 if (false/*status == atLowerBound 539 &&keyValue(iSet)<1.0e-7*/) { 540 // can't come in because 541 // of ones at ub 542 numberWanted++; 543 } else { 544 535 545 bestDj = value; 536 546 bestSequence = structuralOffset + iSequence; 537 547 bestDjMod = djMod; 538 548 bestSet = iSet; 549 } 539 550 } else { 540 551 // just to make sure we don't exit before got something … … 576 587 ) 577 588 { 578 // forceRefresh=true;589 // forceRefresh=true;printf("take out forceRefresh\n"); 579 590 if (!model_->numberIterations()) 580 591 forceRefresh = true; … … 735 746 } else if (getDynamicStatus(j) == atUpperBound) { 736 747 value = columnUpper_[j]; 748 assert (value<1.0e30); 737 749 } else if (getDynamicStatus(j) == soloKey) { 738 750 value = keyValue(iSet); 739 751 } 740 752 objectiveOffset += value * cost_[j]; 741 753 } … … 753 765 for (iSet = 0; iSet < numberSets_; iSet++) { 754 766 int kRow = toIndex_[iSet]; 755 if (kRow >= 0) 756 kRow += numberStaticRows_; 757 int j = startSet_[iSet]; 758 while (j >= 0) { 759 double value = solution[j]; 760 if (value) { 761 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) { 762 int iRow = row_[k]; 763 rhsOffset_[iRow] -= element_[k] * value; 764 } 765 if (kRow >= 0) 766 rhsOffset_[kRow] -= value; 767 } 768 j = next_[j]; //onto next in set 769 } 767 if (kRow >= 0) 768 kRow += numberStaticRows_; 769 int j = startSet_[iSet]; 770 while (j >= 0) { 771 //? use DynamicStatus status = getDynamicStatus(j); 772 double value = solution[j]; 773 if (value) { 774 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) { 775 int iRow = row_[k]; 776 rhsOffset_[iRow] -= element_[k] * value; 777 } 778 if (kRow >= 0) 779 rhsOffset_[kRow] -= value; 780 } 781 j = next_[j]; //onto next in set 782 } 770 783 } 771 784 delete [] solution; … … 1073 1086 // save status 1074 1087 case 5: { 1088 memcpy(status_+numberSets_,status_,numberSets_); 1089 memcpy(dynamicStatus_+maximumGubColumns_, 1090 dynamicStatus_,maximumGubColumns_); 1075 1091 } 1076 1092 break; 1077 1093 // restore status 1078 1094 case 6: { 1079 // maybe mismatch - redo all in small model 1080 //for (int i=firstDynamic_;i<firstAvailable_;i++) { 1081 //int jColumn = id_[i-firstDynamic_]; 1082 //setDynamicStatus(jColumn,inSmall); 1083 //} 1095 memcpy(status_,status_+numberSets_,numberSets_); 1096 memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_, 1097 maximumGubColumns_); 1098 initialProblem(); 1084 1099 } 1085 1100 break; … … 1172 1187 // first flag 1173 1188 if (number >= firstDynamic_ && number < lastDynamic_) { 1174 assert (number == model->sequenceIn()); 1175 // id will be sitting at firstAvailable 1176 int sequence = id_[firstAvailable_-firstDynamic_]; 1177 setFlagged(sequence); 1178 model->clearFlagged(firstAvailable_); 1179 } 1189 int sequence = id_[number-firstDynamic_]; 1190 setFlagged(sequence); 1191 //model->clearFlagged(number); 1192 } else if (number>=model_->numberColumns()+numberStaticRows_) { 1193 // slack 1194 int iSet = fromIndex_[number-model_->numberColumns()- 1195 numberStaticRows_]; 1196 setFlaggedSlack(iSet); 1197 //model->clearFlagged(number); 1198 } 1180 1199 } 1181 1200 case 11: { 1182 int sequenceIn = model->sequenceIn();1183 if ( sequenceIn >= firstDynamic_ && sequenceIn< lastDynamic_) {1184 1201 //int sequenceIn = model->sequenceIn(); 1202 if (number >= firstDynamic_ && number < lastDynamic_) { 1203 //assert (number == model->sequenceIn()); 1185 1204 // take out variable (but leave key) 1186 1205 double * cost = model->costRegion(); … … 1190 1209 int * length = matrix_->getMutableVectorLengths(); 1191 1210 #ifndef NDEBUG 1192 if (length[ sequenceIn]) {1211 if (length[number]) { 1193 1212 int * row = matrix_->getMutableIndices(); 1194 1213 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); 1195 int iRow = row[startColumn[ sequenceIn] + length[sequenceIn] - 1];1214 int iRow = row[startColumn[number] + length[number] - 1]; 1196 1215 int which = iRow - numberStaticRows_; 1197 1216 assert (which >= 0); … … 1210 1229 1211 1230 // not really in small problem 1212 int iBig = id_[ sequenceIn-firstDynamic_];1213 if (model->getStatus( sequenceIn) == ClpSimplex::atLowerBound) {1231 int iBig = id_[number-firstDynamic_]; 1232 if (model->getStatus(number) == ClpSimplex::atLowerBound) { 1214 1233 setDynamicStatus(iBig, atLowerBound); 1215 1234 if (columnLower_) 1216 modifyOffset( sequenceIn, columnLower_[iBig]);1235 modifyOffset(number, columnLower_[iBig]); 1217 1236 } else { 1218 1237 setDynamicStatus(iBig, atUpperBound); 1219 modifyOffset(sequenceIn, columnUpper_[iBig]); 1220 } 1238 modifyOffset(number, columnUpper_[iBig]); 1239 } 1240 } else if (number>=model_->numberColumns()+numberStaticRows_) { 1241 // slack 1242 int iSet = fromIndex_[number-model_->numberColumns()- 1243 numberStaticRows_]; 1244 printf("what now - set %d\n",iSet); 1221 1245 } 1222 1246 } … … 1302 1326 double value = solution[iColumn]; 1303 1327 bool toLowerBound = true; 1328 assert (jColumn>=0);assert (iColumn>=0); 1304 1329 if (columnUpper_) { 1305 1330 if (!columnLower_) { 1306 1331 if (fabs(value - columnUpper_[jColumn]) < fabs(value)) 1307 1332 toLowerBound = false; 1308 } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[ iColumn])) {1333 } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) { 1309 1334 toLowerBound = false; 1310 1335 } … … 1509 1534 n++; 1510 1535 } 1536 int k=0; 1537 for (int j=startSet_[i];j<startSet_[i+1];j++) { 1538 if (getDynamicStatus(j)==soloKey) 1539 k++; 1540 } 1541 assert (k<2); 1511 1542 } 1512 1543 assert (n == numberSets_); … … 1597 1628 j = next_[j]; //onto next in set 1598 1629 } 1599 } 1630 #if 0 1631 // slack must be feasible 1632 double oldValue=value; 1633 value = CoinMax(value,lowerSet_[iSet]); 1634 value = CoinMin(value,upperSet_[iSet]); 1635 if (value!=oldValue) 1636 printf("using %g (not %g) for slack on set %d (%g,%g)\n", 1637 value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]); 1638 #endif 1639 } 1600 1640 } 1601 1641 return value; … … 1676 1716 model->setStatus(iSequence, getStatus(savedBestSet_)); 1677 1717 model->djRegion()[iSequence] = savedBestGubDual_; 1678 if (getStatus(savedBestSet_) == ClpSimplex::atUpperBound) 1679 solution[iSequence] = columnUpper[iSequence]; 1680 else 1681 solution[iSequence] = columnLower[iSequence]; 1682 // create variable and pivot in 1718 solution[iSequence] = valueOfKey; 1719 // create variable and pivot in 1683 1720 int key = keyVariable_[savedBestSet_]; 1684 1721 setDynamicStatus(key, inSmall); … … 1770 1807 model->setStatus(iSequence, ClpSimplex::basic); 1771 1808 model->djRegion()[iSequence] = 0.0; 1772 solution[iSequence] = shift;1809 solution[iSequence] = valueOfKey+shift; 1773 1810 rhsOffset_[newRow] = -shift; // sign? 1774 1811 } … … 1816 1853 //printf("best %d\n",bestSequence2); 1817 1854 model->solutionRegion()[firstAvailable_] = 0.0; 1855 model->clearFlagged(firstAvailable_); 1818 1856 if (!columnLower_ && !columnUpper_) { 1819 1857 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound); … … 2109 2147 int numberActive = 0; 2110 2148 int whichKey = -1; 2111 if (getStatus(iSet) == ClpSimplex::basic) 2149 if (getStatus(iSet) == ClpSimplex::basic) { 2112 2150 whichKey = maximumGubColumns_ + iSet; 2113 else 2151 numberActive = 1; 2152 } else { 2114 2153 whichKey = -1; 2154 } 2115 2155 int j = startSet_[iSet]; 2116 2156 while (j >= 0) { 2117 2157 assert (getDynamicStatus(j) != soloKey || whichKey == -1); 2118 if (getDynamicStatus(j) == inSmall) 2158 if (getDynamicStatus(j) == inSmall) { 2119 2159 numberActive++; 2120 else if (getDynamicStatus(j) == soloKey)2160 } else if (getDynamicStatus(j) == soloKey) { 2121 2161 whichKey = j; 2162 numberActive++; 2163 } 2122 2164 j = next_[j]; //onto next in set 2123 2165 } 2124 if (getStatus(iSet) == ClpSimplex::basic && numberActive) 2125 numberActive++; 2126 if (numberActive) { 2127 assert (numberActive > 1); 2166 if (numberActive > 1) { 2128 2167 int iRow = numberActiveSets_ + numberStaticRows_; 2129 2168 rowSolution[iRow] = 0.0; … … 2189 2228 else 2190 2229 columnUpper[firstAvailable_] = COIN_DBL_MAX; 2191 if (status == atLowerBound) {2230 if (status != atUpperBound) { 2192 2231 solution[firstAvailable_] = columnLower[firstAvailable_]; 2193 2232 } else { … … 2203 2242 toIndex_[iSet] = numberActiveSets_; 2204 2243 fromIndex_[numberActiveSets_++] = iSet; 2244 } else { 2245 // solo key 2246 bool needKey; 2247 if (numberActive) { 2248 if (whichKey<maximumGubColumns_) { 2249 // structural - assume ok 2250 needKey = false; 2251 } else { 2252 // slack 2253 keyVariable_[iSet] = maximumGubColumns_ + iSet; 2254 double value = keyValue(iSet); 2255 if (value<lowerSet_[iSet]-1.0e-8|| 2256 value>upperSet_[iSet]+1.0e-8) 2257 needKey=true; 2258 } 2259 } else { 2260 needKey = true; 2261 } 2262 if (needKey) { 2263 // all to lb then up some (slack/null if possible) 2264 int length=99999999; 2265 int which=-1; 2266 double sum=0.0; 2267 for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) { 2268 setDynamicStatus(iColumn,atLowerBound); 2269 sum += columnLower_[iColumn]; 2270 if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) { 2271 which=iColumn; 2272 length=startColumn_[iColumn+1]-startColumn_[iColumn]; 2273 } 2274 } 2275 if (sum>lowerSet_[iSet]-1.0e-8) { 2276 // slack can be basic 2277 setStatus(iSet,ClpSimplex::basic); 2278 keyVariable_[iSet] = maximumGubColumns_ + iSet; 2279 } else { 2280 // use shortest 2281 setDynamicStatus(which,soloKey); 2282 keyVariable_[iSet] = which; 2283 setStatus(iSet,ClpSimplex::atLowerBound); 2284 } 2285 } 2205 2286 } 2206 2287 assert (toIndex_[iSet] >= 0 || whichKey >= 0); 2207 2288 keyVariable_[iSet] = whichKey; 2208 2289 } 2290 numberActiveColumns_ = firstAvailable_; 2209 2291 return; 2292 } 2293 // Writes out model (without names) 2294 void 2295 ClpDynamicMatrix::writeMps(const char * name) 2296 { 2297 int numberTotalRows = numberStaticRows_+numberSets_; 2298 int numberTotalColumns = firstDynamic_+numberGubColumns_; 2299 // over estimate 2300 int numberElements = getNumElements()+startColumn_[numberGubColumns_] 2301 + numberGubColumns_; 2302 double * columnLower = new double [numberTotalColumns]; 2303 double * columnUpper = new double [numberTotalColumns]; 2304 double * cost = new double [numberTotalColumns]; 2305 double * rowLower = new double [numberTotalRows]; 2306 double * rowUpper = new double [numberTotalRows]; 2307 CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1]; 2308 int * row = new int [numberElements]; 2309 double * element = new double [numberElements]; 2310 // Fill in 2311 const CoinBigIndex * startA = getVectorStarts(); 2312 const int * lengthA = getVectorLengths(); 2313 const int * rowA = getIndices(); 2314 const double * elementA = getElements(); 2315 const double * columnLowerA = model_->columnLower(); 2316 const double * columnUpperA = model_->columnUpper(); 2317 const double * costA = model_->objective(); 2318 const double * rowLowerA = model_->rowLower(); 2319 const double * rowUpperA = model_->rowUpper(); 2320 start[0]=0; 2321 numberElements=0; 2322 for (int i=0;i<firstDynamic_;i++) { 2323 columnLower[i] = columnLowerA[i]; 2324 columnUpper[i] = columnUpperA[i]; 2325 cost[i] = costA[i]; 2326 for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) { 2327 row[numberElements] = rowA[j]; 2328 element[numberElements++]=elementA[j]; 2329 } 2330 start[i+1]=numberElements; 2331 } 2332 for (int i=0;i<numberStaticRows_;i++) { 2333 rowLower[i] = rowLowerA[i]; 2334 rowUpper[i] = rowUpperA[i]; 2335 } 2336 int putC=firstDynamic_; 2337 int putR=numberStaticRows_; 2338 for (int i=0;i<numberSets_;i++) { 2339 rowLower[putR]=lowerSet_[i]; 2340 rowUpper[putR]=upperSet_[i]; 2341 for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) { 2342 columnLower[putC]=columnLower_[k]; 2343 columnUpper[putC]=columnUpper_[k]; 2344 cost[putC]=cost_[k]; 2345 putC++; 2346 for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) { 2347 row[numberElements] = row_[j]; 2348 element[numberElements++]=element_[j]; 2349 } 2350 row[numberElements] = putR; 2351 element[numberElements++]=1.0; 2352 start[putC]=numberElements; 2353 } 2354 putR++; 2355 } 2356 2357 assert (putR==numberTotalRows); 2358 assert (putC==numberTotalColumns); 2359 ClpSimplex modelOut; 2360 modelOut.loadProblem(numberTotalColumns,numberTotalRows, 2361 start,row,element, 2362 columnLower,columnUpper,cost, 2363 rowLower,rowUpper); 2364 modelOut.writeMps(name); 2365 delete [] columnLower; 2366 delete [] columnUpper; 2367 delete [] cost; 2368 delete [] rowLower; 2369 delete [] rowUpper; 2370 delete [] start; 2371 delete [] row; 2372 delete [] element; 2210 2373 } 2211 2374 // Adds in a column to gub structure (called from descendant) … … 2236 2399 if (columnUpper_ && upper != columnUpper_[j]) 2237 2400 odd = true; 2238 if (odd) 2401 if (odd) { 2239 2402 printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n", 2240 2403 cost, lower, upper, cost_[j], 2241 2404 columnLower_ ? columnLower_[j] : 0.0, 2242 2405 columnUpper_ ? columnUpper_[j] : 1.0e100); 2243 else 2406 } else { 2407 setDynamicStatus(j, status); 2244 2408 return j; 2409 } 2245 2410 } 2246 2411 } … … 2259 2424 for (i = 0; i < numberGubColumns_; i++) { 2260 2425 CoinBigIndex end = startColumn_[i+1]; 2426 // what about ubs if column generation? 2261 2427 if (getDynamicStatus(i) != atLowerBound) { 2262 2428 // keep in -
trunk/Clp/src/ClpDynamicMatrix.hpp
r1665 r1676 95 95 /// Does gub crash 96 96 void gubCrash(); 97 /// Writes out model (without names) 98 void writeMps(const char * name); 97 99 /// Populates initial matrix from dynamic status 98 100 void initialProblem(); … … 166 168 st_byte = static_cast<unsigned char>(st_byte | status); 167 169 } 170 /// Whether flagged slack 171 inline bool flaggedSlack(int i) const { 172 return (status_[i] & 8) != 0; 173 } 174 inline void setFlaggedSlack(int i) { 175 status_[i] = static_cast<unsigned char>(status_[i] | 8); 176 } 177 inline void unsetFlaggedSlack(int i) { 178 status_[i] = static_cast<unsigned char>(status_[i] & ~8); 179 } 168 180 /// Number of sets (dynamic rows) 169 181 inline int numberSets() const { 170 182 return numberSets_; 171 183 } 184 /// Number of possible gub variables 185 inline int numberGubEntries() const 186 { return startSet_[numberSets_];} 187 /// Sets 188 inline int * startSets() const 189 { return startSet_;} 172 190 /// Whether flagged 173 191 inline bool flagged(int i) const { -
trunk/Clp/src/ClpNonLinearCost.cpp
r1665 r1676 608 608 } 609 609 } 610 //#define PRINT_DETAIL7 2 611 #if PRINT_DETAIL7>1 612 printf("NL %d sol %g bounds %g %g\n", 613 iSequence,solution[iSequence], 614 lowerValue,upperValue); 615 #endif 610 616 switch(status) { 611 617 … … 625 631 printf("nonlincostb %d %g %g %g\n", 626 632 iSequence, lowerValue, solution[iSequence], lower_[iRange+2]); 633 #endif 634 #if PRINT_DETAIL7 635 printf("**NL %d sol %g below %g\n", 636 iSequence,solution[iSequence], 637 lowerValue); 627 638 #endif 628 639 sumInfeasibilities_ += value; … … 642 653 printf("nonlincostu %d %g %g %g\n", 643 654 iSequence, lower_[iRange-1], solution[iSequence], upperValue); 655 #endif 656 #if PRINT_DETAIL7 657 printf("**NL %d sol %g above %g\n", 658 iSequence,solution[iSequence], 659 upperValue); 644 660 #endif 645 661 sumInfeasibilities_ += value; -
trunk/Clp/src/ClpSimplex.cpp
r1665 r1676 1933 1933 int cycle = progress_.cycle(in, out, 1934 1934 directionIn_, directionOut_); 1935 if (cycle > 0 && objective_->type() < 2 ) {1935 if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) { 1936 1936 //if (cycle>0) { 1937 1937 if (handler_->logLevel() >= 63) … … 1945 1945 forceFactorization_ = CoinMax(1, cycle - off[extra]); 1946 1946 } else { 1947 // need to reject something 1947 /* need to reject something 1948 should be better if don't reject incoming 1949 as it is in basis */ 1948 1950 int iSequence; 1949 if (algorithm_ > 0)1950 1951 else1951 //if (algorithm_ > 0) 1952 // iSequence = sequenceIn_; 1953 //else 1952 1954 iSequence = sequenceOut_; 1953 1955 char x = isColumn(iSequence) ? 'C' : 'R'; … … 1993 1995 forceFactorization_ = -1; //off 1994 1996 return 1; 1995 } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) ) {1997 } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) { 1996 1998 double random = randomNumberGenerator_.randomDouble(); 1997 1999 int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots); … … 8250 8252 whatsChanged_ &= ~0xffff; 8251 8253 } 8254 double saveObjValue = objectiveValue_; 8252 8255 deleteRim(getRidOfData); 8256 if (matrix_->type()>=15) 8257 objectiveValue_ = saveObjValue; 8253 8258 // Skip message if changing algorithms 8254 8259 if (problemStatus_ != 10) { -
trunk/Clp/src/ClpSimplexOther.cpp
r1665 r1676 16 16 #include "ClpFactorization.hpp" 17 17 #include "ClpDualRowDantzig.hpp" 18 #include "ClpDynamicMatrix.hpp" 18 19 #include "CoinPackedMatrix.hpp" 19 20 #include "CoinIndexedVector.hpp" … … 4068 4069 checkSolutionInternal(); 4069 4070 } 4071 // Returns gub version of model or NULL 4072 ClpSimplex * 4073 ClpSimplexOther::gubVersion(int * whichRows, int * whichColumns, 4074 int neededGub, 4075 int factorizationFrequency) 4076 { 4077 // find gub 4078 int numberRows = this->numberRows(); 4079 int numberColumns = this->numberColumns(); 4080 int iRow, iColumn; 4081 int * columnIsGub = new int [numberColumns]; 4082 const double * columnLower = this->columnLower(); 4083 const double * columnUpper = this->columnUpper(); 4084 int numberFixed=0; 4085 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 4086 if (columnUpper[iColumn] == columnLower[iColumn]) { 4087 columnIsGub[iColumn]=-2; 4088 numberFixed++; 4089 } else if (columnLower[iColumn]>=0) { 4090 columnIsGub[iColumn]=-1; 4091 } else { 4092 columnIsGub[iColumn]=-3; 4093 } 4094 } 4095 CoinPackedMatrix * matrix = this->matrix(); 4096 // get row copy 4097 CoinPackedMatrix rowCopy = *matrix; 4098 rowCopy.reverseOrdering(); 4099 const int * column = rowCopy.getIndices(); 4100 const int * rowLength = rowCopy.getVectorLengths(); 4101 const CoinBigIndex * rowStart = rowCopy.getVectorStarts(); 4102 const double * element = rowCopy.getElements(); 4103 int numberNonGub = 0; 4104 int numberEmpty = numberRows; 4105 int * rowIsGub = new int [numberRows]; 4106 int smallestGubRow=-1; 4107 int count=numberColumns+1; 4108 double * rowLower = this->rowLower(); 4109 double * rowUpper = this->rowUpper(); 4110 // make sure we can get rid of upper bounds 4111 double * fixedRow = new double [numberRows]; 4112 for (iRow = 0 ; iRow < numberRows ; iRow++) { 4113 double sumFixed=0.0; 4114 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) { 4115 int iColumn = column[j]; 4116 double value = columnLower[iColumn]; 4117 if (value) 4118 sumFixed += element[j] * value; 4119 } 4120 fixedRow[iRow]=rowUpper[iRow]-sumFixed; 4121 } 4122 for (iRow = numberRows - 1; iRow >= 0; iRow--) { 4123 bool gubRow = true; 4124 int numberInRow=0; 4125 double sumFixed=0.0; 4126 double gap = fixedRow[iRow]-1.0e-12; 4127 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) { 4128 int iColumn = column[j]; 4129 if (columnIsGub[iColumn]!=-2) { 4130 if (element[j] != 1.0||columnIsGub[iColumn]==-3|| 4131 columnUpper[iColumn]-columnLower[iColumn]<gap) { 4132 gubRow = false; 4133 break; 4134 } else { 4135 numberInRow++; 4136 if (columnIsGub[iColumn] >= 0) { 4137 gubRow = false; 4138 break; 4139 } 4140 } 4141 } else { 4142 sumFixed += columnLower[iColumn]*element[j]; 4143 } 4144 } 4145 if (!gubRow) { 4146 whichRows[numberNonGub++] = iRow; 4147 rowIsGub[iRow] = -1; 4148 } else if (numberInRow) { 4149 if (numberInRow<count) { 4150 count = numberInRow; 4151 smallestGubRow=iRow; 4152 } 4153 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) { 4154 int iColumn = column[j]; 4155 if (columnIsGub[iColumn]!=-2) 4156 columnIsGub[iColumn] = iRow; 4157 } 4158 rowIsGub[iRow] = 0; 4159 } else { 4160 // empty row! 4161 whichRows[--numberEmpty] = iRow; 4162 rowIsGub[iRow] = -2; 4163 if (sumFixed>rowUpper[iRow]+1.0e-4|| 4164 sumFixed<rowLower[iRow]-1.0e-4) { 4165 fprintf(stderr,"******** No infeasible empty rows - please!\n"); 4166 abort(); 4167 } 4168 } 4169 } 4170 delete [] fixedRow; 4171 char message[100]; 4172 int numberGub = numberEmpty - numberNonGub; 4173 if (numberGub >= neededGub) { 4174 sprintf(message,"%d gub rows", numberGub); 4175 handler_->message(CLP_GENERAL, messages_) 4176 << message << CoinMessageEol; 4177 int numberNormal = 0; 4178 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 4179 if (columnIsGub[iColumn] < 0 && columnIsGub[iColumn] !=-2) { 4180 whichColumns[numberNormal++] = iColumn; 4181 } 4182 } 4183 if (!numberNormal) { 4184 sprintf(message,"Putting back one gub row to make non-empty"); 4185 handler_->message(CLP_GENERAL, messages_) 4186 << message << CoinMessageEol; 4187 rowIsGub[smallestGubRow]=-1; 4188 whichRows[numberNonGub++] = smallestGubRow; 4189 for (int j = rowStart[smallestGubRow]; 4190 j < rowStart[smallestGubRow] + rowLength[smallestGubRow]; j++) { 4191 int iColumn = column[j]; 4192 if (columnIsGub[iColumn]>=0) { 4193 columnIsGub[iColumn]=-4; 4194 whichColumns[numberNormal++] = iColumn; 4195 } 4196 } 4197 } 4198 std::sort(whichRows,whichRows+numberNonGub); 4199 std::sort(whichColumns,whichColumns+numberNormal); 4200 double * lower = CoinCopyOfArray(this->rowLower(),numberRows); 4201 double * upper = CoinCopyOfArray(this->rowUpper(),numberRows); 4202 // leave empty rows at end 4203 numberEmpty = numberRows-numberEmpty; 4204 const int * row = matrix->getIndices(); 4205 const int * columnLength = matrix->getVectorLengths(); 4206 const CoinBigIndex * columnStart = matrix->getVectorStarts(); 4207 const double * elementByColumn = matrix->getElements(); 4208 // Fixed at end 4209 int put2 = numberColumns-numberFixed; 4210 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 4211 if (columnIsGub[iColumn] ==-2) { 4212 whichColumns[put2++] = iColumn; 4213 double value = columnLower[iColumn]; 4214 for (int j = columnStart[iColumn]; 4215 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 4216 int iRow = row[j]; 4217 if (lower[iRow]>-1.0e20) 4218 lower[iRow] -= value*element[j]; 4219 if (upper[iRow]<1.0e20) 4220 upper[iRow] -= value*element[j]; 4221 } 4222 } 4223 } 4224 int put = numberNormal; 4225 ClpSimplex * model2 = 4226 new ClpSimplex(this, numberNonGub, whichRows , numberNormal, whichColumns); 4227 // scale 4228 double * scaleArray = new double [numberRows]; 4229 for (int i=0;i<numberRows;i++) { 4230 scaleArray[i]=1.0; 4231 if (rowIsGub[i]==-1) { 4232 double largest = 1.0e-30; 4233 double smallest = 1.0e30; 4234 for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) { 4235 int iColumn = column[j]; 4236 if (columnIsGub[iColumn]!=-2) { 4237 double value =fabs(element[j]); 4238 largest = CoinMax(value,largest); 4239 smallest = CoinMin(value,smallest); 4240 } 4241 } 4242 double scale = CoinMax(0.001,1.0/sqrt(largest*smallest)); 4243 scaleArray[i]=scale; 4244 if (lower[i]>-1.0e30) 4245 lower[i] *= scale; 4246 if (upper[i]<1.0e30) 4247 upper[i] *= scale; 4248 } 4249 } 4250 // scale partial matrix 4251 { 4252 CoinPackedMatrix * matrix = model2->matrix(); 4253 const int * row = matrix->getIndices(); 4254 const int * columnLength = matrix->getVectorLengths(); 4255 const CoinBigIndex * columnStart = matrix->getVectorStarts(); 4256 double * element = matrix->getMutableElements(); 4257 for (int i=0;i<numberNormal;i++) { 4258 for (int j = columnStart[i]; 4259 j < columnStart[i] + columnLength[i]; j++) { 4260 int iRow = row[j]; 4261 iRow = whichRows[iRow]; 4262 double scaleBy = scaleArray[iRow]; 4263 element[j] *= scaleBy; 4264 } 4265 } 4266 } 4267 // adjust rhs 4268 double * rowLower = model2->rowLower(); 4269 double * rowUpper = model2->rowUpper(); 4270 for (int i=0;i<numberNonGub;i++) { 4271 int iRow = whichRows[i]; 4272 rowLower[i] = lower[iRow]; 4273 rowUpper[i] = upper[iRow]; 4274 } 4275 int numberGubColumns = numberColumns - put - numberFixed; 4276 CoinBigIndex numberElements=0; 4277 int * temp1 = new int [numberRows+1]; 4278 // get counts 4279 memset(temp1,0,numberRows*sizeof(int)); 4280 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 4281 int iGub = columnIsGub[iColumn]; 4282 if (iGub>=0) { 4283 numberElements += columnLength[iColumn]-1; 4284 temp1[iGub]++; 4285 } 4286 } 4287 /* Optional but means can eventually simplify coding 4288 could even add in fixed slacks to deal with 4289 singularities - but should not be necessary */ 4290 int numberSlacks=0; 4291 for (int i = 0; i < numberRows; i++) { 4292 if (rowIsGub[i]>=0) { 4293 if (lower[i]<upper[i]) { 4294 numberSlacks++; 4295 temp1[i]++; 4296 } 4297 } 4298 } 4299 int * gubStart = new int [numberGub+1]; 4300 numberGub=0; 4301 gubStart[0]=0; 4302 for (int i = 0; i < numberRows; i++) { 4303 if (rowIsGub[i]>=0) { 4304 rowIsGub[i]=numberGub; 4305 gubStart[numberGub+1]=gubStart[numberGub]+temp1[i]; 4306 temp1[numberGub]=0; 4307 lower[numberGub]=lower[i]; 4308 upper[numberGub]=upper[i]; 4309 whichRows[numberNonGub+numberGub]=i; 4310 numberGub++; 4311 } 4312 } 4313 int numberGubColumnsPlus = numberGubColumns + numberSlacks; 4314 double * lowerColumn2 = new double [numberGubColumnsPlus]; 4315 CoinFillN(lowerColumn2, numberGubColumnsPlus, 0.0); 4316 double * upperColumn2 = new double [numberGubColumnsPlus]; 4317 CoinFillN(upperColumn2, numberGubColumnsPlus, COIN_DBL_MAX); 4318 int * start2 = new int[numberGubColumnsPlus+1]; 4319 int * row2 = new int[numberElements]; 4320 double * element2 = new double[numberElements]; 4321 double * cost2 = new double [numberGubColumnsPlus]; 4322 CoinFillN(cost2, numberGubColumnsPlus, 0.0); 4323 const double * cost = this->objective(); 4324 put = numberNormal; 4325 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 4326 int iGub = columnIsGub[iColumn]; 4327 if (iGub>=0) { 4328 // TEMP 4329 //this->setColUpper(iColumn,COIN_DBL_MAX); 4330 iGub = rowIsGub[iGub]; 4331 assert (iGub>=0); 4332 int kPut = put+gubStart[iGub]+temp1[iGub]; 4333 temp1[iGub]++; 4334 whichColumns[kPut]=iColumn; 4335 } 4336 } 4337 for (int i = 0; i < numberRows; i++) { 4338 if (rowIsGub[i]>=0) { 4339 int iGub = rowIsGub[i]; 4340 if (lower[iGub]<upper[iGub]) { 4341 int kPut = put+gubStart[iGub]+temp1[iGub]; 4342 temp1[iGub]++; 4343 whichColumns[kPut]=iGub+numberColumns; 4344 } 4345 } 4346 } 4347 //this->primal(1); // TEMP 4348 // redo rowIsGub to give lookup 4349 for (int i=0;i<numberRows;i++) 4350 rowIsGub[i]=-1; 4351 for (int i=0;i<numberNonGub;i++) 4352 rowIsGub[whichRows[i]]=i; 4353 start2[0]=0; 4354 numberElements = 0; 4355 for (int i=0;i<numberGubColumnsPlus;i++) { 4356 int iColumn = whichColumns[put++]; 4357 if (iColumn<numberColumns) { 4358 cost2[i] = cost[iColumn]; 4359 lowerColumn2[i] = columnLower[iColumn]; 4360 upperColumn2[i] = columnUpper[iColumn]; 4361 upperColumn2[i] = COIN_DBL_MAX; 4362 for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { 4363 int iRow = row[j]; 4364 double scaleBy = scaleArray[iRow]; 4365 iRow = rowIsGub[iRow]; 4366 if (iRow >= 0) { 4367 row2[numberElements] = iRow; 4368 element2[numberElements++] = elementByColumn[j]*scaleBy; 4369 } 4370 } 4371 } else { 4372 // slack 4373 int iGub = iColumn-numberColumns; 4374 double slack = upper[iGub]-lower[iGub]; 4375 assert (upper[iGub]<1.0e20); 4376 lower[iGub]=upper[iGub]; 4377 cost2[i] = 0; 4378 lowerColumn2[i] = 0; 4379 upperColumn2[i] = slack; 4380 upperColumn2[i] = COIN_DBL_MAX; 4381 } 4382 start2[i+1] = numberElements; 4383 } 4384 // clean up bounds on variables 4385 for (int iSet=0;iSet<numberGub;iSet++) { 4386 double lowerValue=0.0; 4387 for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) { 4388 lowerValue += lowerColumn2[i]; 4389 } 4390 assert (lowerValue<upper[iSet]+1.0e-6); 4391 double gap = CoinMax(0.0,upper[iSet]-lowerValue); 4392 for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) { 4393 if (upperColumn2[i]<1.0e30) { 4394 upperColumn2[i] = CoinMin(upperColumn2[i], 4395 lowerColumn2[i]+gap); 4396 } 4397 } 4398 } 4399 sprintf(message,"** Before adding matrix there are %d rows and %d columns", 4400 model2->numberRows(), model2->numberColumns()); 4401 handler_->message(CLP_GENERAL, messages_) 4402 << message << CoinMessageEol; 4403 delete [] scaleArray; 4404 delete [] temp1; 4405 model2->setFactorizationFrequency(factorizationFrequency); 4406 ClpDynamicMatrix * newMatrix = 4407 new ClpDynamicMatrix(model2, numberGub, 4408 numberGubColumnsPlus, gubStart, 4409 lower, upper, 4410 start2, row2, element2, cost2, 4411 lowerColumn2, upperColumn2); 4412 delete [] gubStart; 4413 delete [] lowerColumn2; 4414 delete [] upperColumn2; 4415 delete [] start2; 4416 delete [] row2; 4417 delete [] element2; 4418 delete [] cost2; 4419 delete [] lower; 4420 delete [] upper; 4421 model2->replaceMatrix(newMatrix,true); 4422 #ifdef EVERY_ITERATION 4423 { 4424 ClpDynamicMatrix * gubMatrix = 4425 dynamic_cast< ClpDynamicMatrix*>(model2->clpMatrix()); 4426 assert(gubMatrix); 4427 gubMatrix->writeMps("gub.mps"); 4428 } 4429 #endif 4430 delete [] columnIsGub; 4431 delete [] rowIsGub; 4432 newMatrix->switchOffCheck(); 4433 #ifdef EVERY_ITERATION 4434 newMatrix->setRefreshFrequency(1/*000*/); 4435 #else 4436 newMatrix->setRefreshFrequency(1000); 4437 #endif 4438 sprintf(message, 4439 "** While after adding matrix there are %d rows and %d columns", 4440 model2->numberRows(), model2->numberColumns()); 4441 handler_->message(CLP_GENERAL, messages_) 4442 << message << CoinMessageEol; 4443 model2->setSpecialOptions(4); // exactly to bound 4444 // Scaling off (done by hand) 4445 model2->scaling(0); 4446 return model2; 4447 } else { 4448 delete [] columnIsGub; 4449 delete [] rowIsGub; 4450 return NULL; 4451 } 4452 } 4453 // Sets basis from original 4454 void 4455 ClpSimplexOther::setGubBasis(const ClpSimplex &original,const int * whichRows, 4456 const int * whichColumns) 4457 { 4458 ClpDynamicMatrix * gubMatrix = 4459 dynamic_cast< ClpDynamicMatrix*>(clpMatrix()); 4460 assert(gubMatrix); 4461 int numberGubColumns = gubMatrix->numberGubColumns(); 4462 int numberNormal = gubMatrix->firstDynamic(); 4463 //int lastOdd = gubMatrix->firstAvailable(); 4464 //int numberTotalColumns = numberNormal + numberGubColumns; 4465 //assert (numberTotalColumns==numberColumns+numberSlacks); 4466 int numberRows = original.numberRows(); 4467 int numberColumns = original.numberColumns(); 4468 int * columnIsGub = new int [numberColumns]; 4469 int numberNonGub = gubMatrix->numberStaticRows(); 4470 //assert (firstOdd==numberNormal); 4471 double * solution = primalColumnSolution(); 4472 const double * originalSolution = original.primalColumnSolution(); 4473 int numberSets = gubMatrix->numberSets(); 4474 const int * startSet = gubMatrix->startSets(); 4475 const CoinBigIndex * startColumn = gubMatrix->startColumn(); 4476 const double * columnLower = gubMatrix->columnLower(); 4477 for (int i=0;i<numberSets;i++) { 4478 for (int j=startSet[i];j<startSet[i+1];j++) { 4479 int iColumn = whichColumns[j+numberNormal]; 4480 if (iColumn<numberColumns) 4481 columnIsGub[iColumn] = whichRows[numberNonGub+i]; 4482 } 4483 } 4484 int * numberKey = new int [numberRows]; 4485 memset(numberKey,0,numberRows*sizeof(int)); 4486 for (int i=0;i<numberGubColumns;i++) { 4487 int iOrig = whichColumns[i+numberNormal]; 4488 if (iOrig<numberColumns) { 4489 if (original.getColumnStatus(iOrig)==ClpSimplex::basic) { 4490 int iRow = columnIsGub[iOrig]; 4491 assert (iRow>=0); 4492 numberKey[iRow]++; 4493 } 4494 } else { 4495 // Set slack 4496 int iSet = iOrig - numberColumns; 4497 int iRow = whichRows[iSet+numberNonGub]; 4498 if (original.getRowStatus(iRow)==ClpSimplex::basic) 4499 numberKey[iRow]++; 4500 } 4501 } 4502 /* Before going into cleanMatrix we need 4503 gub status set (inSmall just means basic and active) 4504 row status set 4505 */ 4506 for (int i = 0; i < numberSets; i++) { 4507 gubMatrix->setStatus(i,ClpSimplex::isFixed); 4508 } 4509 for (int i = 0; i < numberGubColumns; i++) { 4510 int iOrig = whichColumns[i+numberNormal]; 4511 if (iOrig<numberColumns) { 4512 ClpSimplex::Status status = original.getColumnStatus(iOrig); 4513 if (status==ClpSimplex::atUpperBound) { 4514 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atUpperBound); 4515 } else if (status==ClpSimplex::atLowerBound) { 4516 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound); 4517 } else if (status==ClpSimplex::basic) { 4518 int iRow = columnIsGub[iOrig]; 4519 assert (iRow>=0); 4520 assert(numberKey[iRow]); 4521 if (numberKey[iRow]==1) 4522 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey); 4523 else 4524 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall); 4525 } 4526 } else { 4527 // slack 4528 int iSet = iOrig - numberColumns; 4529 int iRow = whichRows[iSet+numberNonGub]; 4530 if (original.getRowStatus(iRow)==ClpSimplex::basic) { 4531 assert(numberKey[iRow]); 4532 if (numberKey[iRow]==1) 4533 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey); 4534 else 4535 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall); 4536 } else { 4537 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound); 4538 } 4539 } 4540 } 4541 // deal with sets without key 4542 for (int i = 0; i < numberSets; i++) { 4543 int iRow = whichRows[numberNonGub+i]; 4544 if (!numberKey[iRow]) { 4545 if (original.getRowStatus(iRow)==ClpSimplex::basic) { 4546 gubMatrix->setStatus(i,ClpSimplex::basic); 4547 } else { 4548 // If not at lb make key otherwise one with smallest number els 4549 double largest=0.0; 4550 int fewest=numberRows+1; 4551 int chosen=-1; 4552 for (int j=startSet[i];j<startSet[i+1];j++) { 4553 int length=startColumn[j+1]-startColumn[j]; 4554 int iOrig = whichColumns[j+numberNormal]; 4555 double value; 4556 if (iOrig<numberColumns) { 4557 value=originalSolution[iOrig]-columnLower[j]; 4558 } else { 4559 // slack - take value as 0.0 as will win on length 4560 value=0.0; 4561 } 4562 if (value>largest+1.0e-8) { 4563 largest=value; 4564 fewest=length; 4565 chosen=j; 4566 } else if (fabs(value-largest)<=1.0e-8&&length<fewest) { 4567 largest=value; 4568 fewest=length; 4569 chosen=j; 4570 } 4571 } 4572 assert(chosen>=0); 4573 // set as key 4574 for (int j=startSet[i];j<startSet[i+1];j++) { 4575 if (j!=chosen) 4576 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound); 4577 else 4578 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::soloKey); 4579 } 4580 } 4581 } 4582 } 4583 for (int i = 0; i < numberNormal; i++) { 4584 int iOrig = whichColumns[i]; 4585 setColumnStatus(i,original.getColumnStatus(iOrig)); 4586 solution[i]=originalSolution[iOrig]; 4587 } 4588 for (int i = 0; i < numberNonGub; i++) { 4589 int iOrig = whichRows[i]; 4590 setRowStatus(i,original.getRowStatus(iOrig)); 4591 } 4592 // Fill in current matrix 4593 gubMatrix->initialProblem(); 4594 delete [] numberKey; 4595 delete [] columnIsGub; 4596 } 4597 // Restores basis to original 4598 void 4599 ClpSimplexOther::getGubBasis(ClpSimplex &original,const int * whichRows, 4600 const int * whichColumns) const 4601 { 4602 ClpDynamicMatrix * gubMatrix = 4603 dynamic_cast< ClpDynamicMatrix*>(clpMatrix()); 4604 assert(gubMatrix); 4605 int numberGubColumns = gubMatrix->numberGubColumns(); 4606 int numberNormal = gubMatrix->firstDynamic(); 4607 //int lastOdd = gubMatrix->firstAvailable(); 4608 //int numberRows = original.numberRows(); 4609 int numberColumns = original.numberColumns(); 4610 int numberNonGub = gubMatrix->numberStaticRows(); 4611 //assert (firstOdd==numberNormal); 4612 double * solution = primalColumnSolution(); 4613 double * originalSolution = original.primalColumnSolution(); 4614 int numberSets = gubMatrix->numberSets(); 4615 const double * cost = original.objective(); 4616 int lastOdd = gubMatrix->firstAvailable(); 4617 //assert (numberTotalColumns==numberColumns+numberSlacks); 4618 int numberRows = original.numberRows(); 4619 //int numberStaticRows = gubMatrix->numberStaticRows(); 4620 const int * startSet = gubMatrix->startSets(); 4621 unsigned char * status = original.statusArray(); 4622 unsigned char * rowStatus = status+numberColumns; 4623 //assert (firstOdd==numberNormal); 4624 for (int i=0;i<numberSets;i++) { 4625 int iRow = whichRows[i+numberNonGub]; 4626 original.setRowStatus(iRow,ClpSimplex::atLowerBound); 4627 } 4628 const int * id = gubMatrix->id(); 4629 const double * columnLower = gubMatrix->columnLower(); 4630 const double * columnUpper = gubMatrix->columnUpper(); 4631 for (int i = 0; i < numberGubColumns; i++) { 4632 int iOrig = whichColumns[i+numberNormal]; 4633 if (iOrig<numberColumns) { 4634 if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) { 4635 originalSolution[iOrig] = columnUpper[i]; 4636 status[iOrig] = 2; 4637 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound && columnLower) { 4638 originalSolution[iOrig] = columnLower[i]; 4639 status[iOrig] = 3; 4640 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) { 4641 int iSet = gubMatrix->whichSet(i); 4642 originalSolution[iOrig] = gubMatrix->keyValue(iSet); 4643 status[iOrig] = 1; 4644 } else { 4645 originalSolution[iOrig] = 0.0; 4646 status[iOrig] = 4; 4647 } 4648 } else { 4649 // slack 4650 int iSet = iOrig - numberColumns; 4651 int iRow = whichRows[iSet+numberNonGub]; 4652 if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) { 4653 original.setRowStatus(iRow,ClpSimplex::atLowerBound); 4654 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound) { 4655 original.setRowStatus(iRow,ClpSimplex::atUpperBound); 4656 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) { 4657 original.setRowStatus(iRow,ClpSimplex::basic); 4658 } 4659 } 4660 } 4661 for (int i = 0; i < numberNormal; i++) { 4662 int iOrig = whichColumns[i]; 4663 ClpSimplex::Status thisStatus = getStatus(i); 4664 if (thisStatus == ClpSimplex::basic) 4665 status[iOrig] = 1; 4666 else if (thisStatus == ClpSimplex::atLowerBound) 4667 status[iOrig] = 3; 4668 else if (thisStatus == ClpSimplex::atUpperBound) 4669 status[iOrig] = 2; 4670 else if (thisStatus == ClpSimplex::isFixed) 4671 status[iOrig] = 5; 4672 else 4673 abort(); 4674 originalSolution[iOrig] = solution[i]; 4675 } 4676 for (int i = numberNormal; i < lastOdd; i++) { 4677 int iOrig = whichColumns[id[i-numberNormal] + numberNormal]; 4678 if (iOrig<numberColumns) { 4679 ClpSimplex::Status thisStatus = getStatus(i); 4680 if (thisStatus == ClpSimplex::basic) 4681 status[iOrig] = 1; 4682 else if (thisStatus == ClpSimplex::atLowerBound) 4683 status[iOrig] = 3; 4684 else if (thisStatus == ClpSimplex::atUpperBound) 4685 status[iOrig] = 2; 4686 else if (thisStatus == ClpSimplex::isFixed) 4687 status[iOrig] = 5; 4688 else 4689 abort(); 4690 originalSolution[iOrig] = solution[i]; 4691 } else { 4692 // slack (basic probably) 4693 int iSet = iOrig - numberColumns; 4694 int iRow = whichRows[iSet+numberNonGub]; 4695 ClpSimplex::Status thisStatus = getStatus(i); 4696 if (thisStatus == ClpSimplex::atLowerBound) 4697 thisStatus = ClpSimplex::atUpperBound; 4698 else if (thisStatus == ClpSimplex::atUpperBound) 4699 thisStatus = ClpSimplex::atLowerBound; 4700 original.setRowStatus(iRow,thisStatus); 4701 } 4702 } 4703 for (int i = 0; i < numberNonGub; i++) { 4704 int iOrig = whichRows[i]; 4705 ClpSimplex::Status thisStatus = getRowStatus(i); 4706 if (thisStatus == ClpSimplex::basic) 4707 rowStatus[iOrig] = 1; 4708 else if (thisStatus == ClpSimplex::atLowerBound) 4709 rowStatus[iOrig] = 3; 4710 else if (thisStatus == ClpSimplex::atUpperBound) 4711 rowStatus[iOrig] = 2; 4712 else if (thisStatus == ClpSimplex::isFixed) 4713 rowStatus[iOrig] = 5; 4714 else 4715 abort(); 4716 } 4717 int * numberKey = new int [numberRows]; 4718 memset(numberKey,0,numberRows*sizeof(int)); 4719 for (int i=0;i<numberSets;i++) { 4720 int iRow = whichRows[i+numberNonGub]; 4721 for (int j=startSet[i];j<startSet[i+1];j++) { 4722 int iOrig = whichColumns[j+numberNormal]; 4723 if (iOrig<numberColumns) { 4724 if (original.getColumnStatus(iOrig)==ClpSimplex::basic) { 4725 numberKey[iRow]++; 4726 } 4727 } else { 4728 // slack 4729 if (original.getRowStatus(iRow)==ClpSimplex::basic) 4730 numberKey[iRow]++; 4731 } 4732 } 4733 } 4734 for (int i=0;i<numberSets;i++) { 4735 int iRow = whichRows[i+numberNonGub]; 4736 if (!numberKey[iRow]) { 4737 original.setRowStatus(iRow,ClpSimplex::basic); 4738 } 4739 } 4740 delete [] numberKey; 4741 double objValue = 0.0; 4742 for (int i = 0; i < numberColumns; i++) 4743 objValue += cost[i] * originalSolution[i]; 4744 //printf("objective value is %g\n", objValue); 4745 } -
trunk/Clp/src/ClpSimplexOther.hpp
r1665 r1676 190 190 const int * whichRows, const int * whichColumns, 191 191 int nBound); 192 /** Returns gub version of model or NULL 193 whichRows has to be numberRows 194 whichColumns has to be numberRows+numberColumns */ 195 ClpSimplex * gubVersion(int * whichRows, int * whichColumns, 196 int neededGub, 197 int factorizationFrequency=50); 198 /// Sets basis from original 199 void setGubBasis(const ClpSimplex &original,const int * whichRows, 200 const int * whichColumns); 201 /// Restores basis to original 202 void getGubBasis(ClpSimplex &original,const int * whichRows, 203 const int * whichColumns) const; 192 204 /// Quick try at cleaning up duals if postsolve gets wrong 193 205 void cleanupAfterPostsolve(); -
trunk/Clp/src/ClpSimplexPrimal.cpp
r1671 r1676 1092 1092 objectiveWork_ = cost_; 1093 1093 rowObjectiveWork_ = cost_ + numberColumns_; 1094 if (n )1094 if (n||matrix_->type()>=15) 1095 1095 gutsOfSolution(NULL, NULL, ifValuesPass != 0); 1096 1096 } … … 1118 1118 && firstFree_ < 0) { 1119 1119 if (handler_->logLevel() == 63) 1120 printf("lastobj %g this %g force %d 1120 printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_); 1121 1121 int maxFactor = factorization_->maximumPivots(); 1122 1122 if (maxFactor > 10) { … … 1130 1130 && firstFree_ < 0) { 1131 1131 if (handler_->logLevel() == 63) 1132 printf("lastobj3 %g this3 %g `force %d", lastObj3, thisObj, forceFactorization_);1132 printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_); 1133 1133 int maxFactor = factorization_->maximumPivots(); 1134 1134 if (maxFactor > 10) { -
trunk/Clp/src/ClpSolve.cpp
r1665 r1676 3025 3025 if (matched == (1 << (CLP_PROGRESS - 1))) 3026 3026 numberMatched = 0; 3027 if (numberMatched ) {3027 if (numberMatched && model_->clpMatrix()->type() < 15) { 3028 3028 model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages()) 3029 3029 << numberMatched
Note: See TracChangeset
for help on using the changeset viewer.