Changeset 1552 for trunk/Clp/examples/decompose.cpp
 Timestamp:
 May 24, 2010 9:03:59 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/examples/decompose.cpp
r1370 r1552 9 9 int main (int argc, const char *argv[]) 10 10 { 11 ClpSimplex model;12 int status;13 // Keep names14 if (argc<2) {15 //status=model.readMps("/home/forrest/data/ken18.mps.gz",true);16 status=model.readMps("../../Data/Netlib/czprob.mps",true);17 } else {18 status=model.readMps(argv[1],true);19 }20 if (status)21 exit(10);22 /*23 This driver does a simple Dantzig Wolfe decomposition24 */25 // Get master rows in some mysterious way26 int numberRows = model.numberRows();27 int * rowBlock = new int[numberRows];28 int iRow;29 for (iRow=0;iRow<numberRows;iRow++)30 rowBlock[iRow]=2;31 // these are master rows32 if (numberRows==105127) {33 // ken1834 for (iRow=104976;iRow<numberRows;iRow++)35 rowBlock[iRow]=1;36 } else if (numberRows==2426) {37 // ken738 for (iRow=2401;iRow<numberRows;iRow++)39 rowBlock[iRow]=1;40 } else if (numberRows==810) {41 for (iRow=81;iRow<84;iRow++)42 rowBlock[iRow]=1;43 } else if (numberRows==5418) {44 for (iRow=564;iRow<603;iRow++)45 rowBlock[iRow]=1;46 } else if (numberRows==10280) {47 // osa6048 for (iRow=10198;iRow<10280;iRow++)49 rowBlock[iRow]=1;50 } else if (numberRows==1503) {51 // degen352 for (iRow=0;iRow<561;iRow++)53 rowBlock[iRow]=1;54 } else if (numberRows==929) {55 // czprob56 for (iRow=0;iRow<39;iRow++)57 rowBlock[iRow]=1;58 }59 CoinPackedMatrix * matrix = model.matrix();60 // get row copy61 CoinPackedMatrix rowCopy = *matrix;62 rowCopy.reverseOrdering();63 const int * row = matrix>getIndices();64 const int * columnLength = matrix>getVectorLengths();65 const CoinBigIndex * columnStart = matrix>getVectorStarts();66 //const double * elementByColumn = matrix>getElements();67 const int * column = rowCopy.getIndices();68 const int * rowLength = rowCopy.getVectorLengths();69 const CoinBigIndex * rowStart = rowCopy.getVectorStarts();70 //const double * elementByRow = rowCopy.getElements();71 int numberBlocks = 0;72 int * stack = new int [numberRows];73 // to say if column looked at74 int numberColumns = model.numberColumns();75 int * columnBlock = new int[numberColumns];76 int iColumn;77 for (iColumn=0;iColumn<numberColumns;iColumn++)78 columnBlock[iColumn]=2;79 for (iColumn=0;iColumn<numberColumns;iColumn++) {80 int kstart = columnStart[iColumn];81 int kend = columnStart[iColumn]+columnLength[iColumn];82 if (columnBlock[iColumn]==2) {83 // column not allocated84 int j;85 int nstack=0;86 for (j=kstart;j<kend;j++) {87 int iRow= row[j];88 if (rowBlock[iRow]!=1) {89 assert(rowBlock[iRow]==2);90 rowBlock[iRow]=numberBlocks; // mark91 92 93 }94 if (nstack) {95 96 97 columnBlock[iColumn]=numberBlocks1;98 99 100 101 for (k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {102 103 104 105 if (columnBlock[iColumn]==2) {106 columnBlock[iColumn]=numberBlocks1; // mark107 108 109 for (jj=kkstart;jj<kkend;jj++) {110 int jRow= row[jj];111 if (rowBlock[jRow]==2) {112 rowBlock[jRow]=numberBlocks1;113 stack[nstack++]=jRow;114 115 116 117 assert (columnBlock[iColumn]==numberBlocks1);118 119 120 121 } else {122 123 columnBlock[iColumn]=1;124 }125 }126 }127 printf("%d blocks found\n",numberBlocks);128 if (numberBlocks>50) {129 int iBlock;130 for (iRow=0;iRow<numberRows;iRow++) {131 iBlock = rowBlock[iRow];132 if (iBlock>=0)133 rowBlock[iRow] = iBlock%50;134 }135 for (iColumn=0;iColumn<numberColumns;iColumn++) {136 iBlock = columnBlock[iColumn];137 if (iBlock>=0)138 columnBlock[iColumn] = iBlock%50;139 }140 numberBlocks=50;141 }142 delete [] stack;143 // make up problems144 CoinPackedMatrix * top = new CoinPackedMatrix [numberBlocks];145 ClpSimplex * sub = new ClpSimplex [numberBlocks];146 ClpSimplex master;147 // Create all sub problems148 // Could do much faster  but do that later149 int * whichRow = new int [numberRows];150 int * whichColumn = new int [numberColumns];151 // get top matrix152 CoinPackedMatrix topMatrix = *model.matrix();153 int numberRow2,numberColumn2;154 numberRow2=0;155 for (iRow=0;iRow<numberRows;iRow++)156 if (rowBlock[iRow]>=0)157 whichRow[numberRow2++]=iRow;158 topMatrix.deleteRows(numberRow2,whichRow);159 int iBlock;160 for (iBlock=0;iBlock<numberBlocks;iBlock++) {161 numberRow2=0;162 numberColumn2=0;163 for (iRow=0;iRow<numberRows;iRow++)164 if (iBlock==rowBlock[iRow])165 whichRow[numberRow2++]=iRow;166 for (iColumn=0;iColumn<numberColumns;iColumn++)167 if (iBlock==columnBlock[iColumn])168 whichColumn[numberColumn2++]=iColumn;169 sub[iBlock]=ClpSimplex(&model,numberRow2,whichRow,170 numberColumn2,whichColumn);11 ClpSimplex model; 12 int status; 13 // Keep names 14 if (argc < 2) { 15 //status=model.readMps("/home/forrest/data/ken18.mps.gz",true); 16 status = model.readMps("../../Data/Netlib/czprob.mps", true); 17 } else { 18 status = model.readMps(argv[1], true); 19 } 20 if (status) 21 exit(10); 22 /* 23 This driver does a simple Dantzig Wolfe decomposition 24 */ 25 // Get master rows in some mysterious way 26 int numberRows = model.numberRows(); 27 int * rowBlock = new int[numberRows]; 28 int iRow; 29 for (iRow = 0; iRow < numberRows; iRow++) 30 rowBlock[iRow] = 2; 31 // these are master rows 32 if (numberRows == 105127) { 33 // ken18 34 for (iRow = 104976; iRow < numberRows; iRow++) 35 rowBlock[iRow] = 1; 36 } else if (numberRows == 2426) { 37 // ken7 38 for (iRow = 2401; iRow < numberRows; iRow++) 39 rowBlock[iRow] = 1; 40 } else if (numberRows == 810) { 41 for (iRow = 81; iRow < 84; iRow++) 42 rowBlock[iRow] = 1; 43 } else if (numberRows == 5418) { 44 for (iRow = 564; iRow < 603; iRow++) 45 rowBlock[iRow] = 1; 46 } else if (numberRows == 10280) { 47 // osa60 48 for (iRow = 10198; iRow < 10280; iRow++) 49 rowBlock[iRow] = 1; 50 } else if (numberRows == 1503) { 51 // degen3 52 for (iRow = 0; iRow < 561; iRow++) 53 rowBlock[iRow] = 1; 54 } else if (numberRows == 929) { 55 // czprob 56 for (iRow = 0; iRow < 39; iRow++) 57 rowBlock[iRow] = 1; 58 } 59 CoinPackedMatrix * matrix = model.matrix(); 60 // get row copy 61 CoinPackedMatrix rowCopy = *matrix; 62 rowCopy.reverseOrdering(); 63 const int * row = matrix>getIndices(); 64 const int * columnLength = matrix>getVectorLengths(); 65 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 66 //const double * elementByColumn = matrix>getElements(); 67 const int * column = rowCopy.getIndices(); 68 const int * rowLength = rowCopy.getVectorLengths(); 69 const CoinBigIndex * rowStart = rowCopy.getVectorStarts(); 70 //const double * elementByRow = rowCopy.getElements(); 71 int numberBlocks = 0; 72 int * stack = new int [numberRows]; 73 // to say if column looked at 74 int numberColumns = model.numberColumns(); 75 int * columnBlock = new int[numberColumns]; 76 int iColumn; 77 for (iColumn = 0; iColumn < numberColumns; iColumn++) 78 columnBlock[iColumn] = 2; 79 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 80 int kstart = columnStart[iColumn]; 81 int kend = columnStart[iColumn] + columnLength[iColumn]; 82 if (columnBlock[iColumn] == 2) { 83 // column not allocated 84 int j; 85 int nstack = 0; 86 for (j = kstart; j < kend; j++) { 87 int iRow = row[j]; 88 if (rowBlock[iRow] != 1) { 89 assert(rowBlock[iRow] == 2); 90 rowBlock[iRow] = numberBlocks; // mark 91 stack[nstack++] = iRow; 92 } 93 } 94 if (nstack) { 95 // new block  put all connected in 96 numberBlocks++; 97 columnBlock[iColumn] = numberBlocks  1; 98 while (nstack) { 99 int iRow = stack[nstack]; 100 int k; 101 for (k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) { 102 int iColumn = column[k]; 103 int kkstart = columnStart[iColumn]; 104 int kkend = kkstart + columnLength[iColumn]; 105 if (columnBlock[iColumn] == 2) { 106 columnBlock[iColumn] = numberBlocks  1; // mark 107 // column not allocated 108 int jj; 109 for (jj = kkstart; jj < kkend; jj++) { 110 int jRow = row[jj]; 111 if (rowBlock[jRow] == 2) { 112 rowBlock[jRow] = numberBlocks  1; 113 stack[nstack++] = jRow; 114 } 115 } 116 } else { 117 assert (columnBlock[iColumn] == numberBlocks  1); 118 } 119 } 120 } 121 } else { 122 // Only in master 123 columnBlock[iColumn] = 1; 124 } 125 } 126 } 127 printf("%d blocks found\n", numberBlocks); 128 if (numberBlocks > 50) { 129 int iBlock; 130 for (iRow = 0; iRow < numberRows; iRow++) { 131 iBlock = rowBlock[iRow]; 132 if (iBlock >= 0) 133 rowBlock[iRow] = iBlock % 50; 134 } 135 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 136 iBlock = columnBlock[iColumn]; 137 if (iBlock >= 0) 138 columnBlock[iColumn] = iBlock % 50; 139 } 140 numberBlocks = 50; 141 } 142 delete [] stack; 143 // make up problems 144 CoinPackedMatrix * top = new CoinPackedMatrix [numberBlocks]; 145 ClpSimplex * sub = new ClpSimplex [numberBlocks]; 146 ClpSimplex master; 147 // Create all sub problems 148 // Could do much faster  but do that later 149 int * whichRow = new int [numberRows]; 150 int * whichColumn = new int [numberColumns]; 151 // get top matrix 152 CoinPackedMatrix topMatrix = *model.matrix(); 153 int numberRow2, numberColumn2; 154 numberRow2 = 0; 155 for (iRow = 0; iRow < numberRows; iRow++) 156 if (rowBlock[iRow] >= 0) 157 whichRow[numberRow2++] = iRow; 158 topMatrix.deleteRows(numberRow2, whichRow); 159 int iBlock; 160 for (iBlock = 0; iBlock < numberBlocks; iBlock++) { 161 numberRow2 = 0; 162 numberColumn2 = 0; 163 for (iRow = 0; iRow < numberRows; iRow++) 164 if (iBlock == rowBlock[iRow]) 165 whichRow[numberRow2++] = iRow; 166 for (iColumn = 0; iColumn < numberColumns; iColumn++) 167 if (iBlock == columnBlock[iColumn]) 168 whichColumn[numberColumn2++] = iColumn; 169 sub[iBlock] = ClpSimplex(&model, numberRow2, whichRow, 170 numberColumn2, whichColumn); 171 171 #if 0 172 // temp173 double * upper = sub[iBlock].columnUpper();174 for (iColumn=0;iColumn<numberColumn2;iColumn++)175 upper[iColumn]=100.0;172 // temp 173 double * upper = sub[iBlock].columnUpper(); 174 for (iColumn = 0; iColumn < numberColumn2; iColumn++) 175 upper[iColumn] = 100.0; 176 176 #endif 177 // and top matrix178 CoinPackedMatrix matrix = topMatrix;179 // and delete bits180 numberColumn2=0;181 for (iColumn=0;iColumn<numberColumns;iColumn++)182 if (iBlock!=columnBlock[iColumn])183 whichColumn[numberColumn2++]=iColumn;184 matrix.deleteCols(numberColumn2,whichColumn);185 top[iBlock]=matrix;186 }187 // and master188 numberRow2=0;189 numberColumn2=0;190 for (iRow=0;iRow<numberRows;iRow++)191 if (rowBlock[iRow]<0)192 whichRow[numberRow2++]=iRow;193 for (iColumn=0;iColumn<numberColumns;iColumn++)194 if (columnBlock[iColumn]==1)195 whichColumn[numberColumn2++]=iColumn;196 ClpModel masterModel(&model,numberRow2,whichRow,197 numberColumn2,whichColumn);198 master=ClpSimplex(masterModel);199 delete [] whichRow;200 delete [] whichColumn;177 // and top matrix 178 CoinPackedMatrix matrix = topMatrix; 179 // and delete bits 180 numberColumn2 = 0; 181 for (iColumn = 0; iColumn < numberColumns; iColumn++) 182 if (iBlock != columnBlock[iColumn]) 183 whichColumn[numberColumn2++] = iColumn; 184 matrix.deleteCols(numberColumn2, whichColumn); 185 top[iBlock] = matrix; 186 } 187 // and master 188 numberRow2 = 0; 189 numberColumn2 = 0; 190 for (iRow = 0; iRow < numberRows; iRow++) 191 if (rowBlock[iRow] < 0) 192 whichRow[numberRow2++] = iRow; 193 for (iColumn = 0; iColumn < numberColumns; iColumn++) 194 if (columnBlock[iColumn] == 1) 195 whichColumn[numberColumn2++] = iColumn; 196 ClpModel masterModel(&model, numberRow2, whichRow, 197 numberColumn2, whichColumn); 198 master = ClpSimplex(masterModel); 199 delete [] whichRow; 200 delete [] whichColumn; 201 201 202 // Overkill in terms of space203 int numberMasterRows = master.numberRows();204 int * columnAdd = new int[numberBlocks+1];205 int * rowAdd = new int[numberBlocks*(numberMasterRows+1)];206 double * elementAdd = new double[numberBlocks*(numberMasterRows+1)];207 double * objective = new double[numberBlocks];208 int maxPass=500;209 int iPass;210 double lastObjective=1.0e31;211 // Create convexity rows for proposals212 int numberMasterColumns = master.numberColumns();213 master.resize(numberMasterRows+numberBlocks,numberMasterColumns);214 // Arrays to say which block and when created215 int maximumColumns = 2*numberMasterRows+10*numberBlocks;216 int * whichBlock = new int[maximumColumns];217 int * when = new int[maximumColumns];218 int numberColumnsGenerated=numberBlocks;219 // fill in rhs and add in artificials220 {221 double * rowLower = master.rowLower();222 double * rowUpper = master.rowUpper();223 int iBlock;224 columnAdd[0] = 0;225 for (iBlock=0;iBlock<numberBlocks;iBlock++) {226 int iRow = iBlock + numberMasterRows;;227 rowLower[iRow]=1.0;228 rowUpper[iRow]=1.0;229 rowAdd[iBlock] = iRow;230 elementAdd[iBlock] = 1.0;231 objective[iBlock] = 1.0e9;232 columnAdd[iBlock+1] = iBlock+1;233 when[iBlock]=1;234 whichBlock[iBlock] = iBlock;235 }236 master.addColumns(numberBlocks,NULL,NULL,objective,237 238 }239 // and resize matrix to double check clp will be happy240 //master.matrix()>setDimensions(numberMasterRows+numberBlocks,241 // numberMasterColumns+numberBlocks);202 // Overkill in terms of space 203 int numberMasterRows = master.numberRows(); 204 int * columnAdd = new int[numberBlocks+1]; 205 int * rowAdd = new int[numberBlocks*(numberMasterRows+1)]; 206 double * elementAdd = new double[numberBlocks*(numberMasterRows+1)]; 207 double * objective = new double[numberBlocks]; 208 int maxPass = 500; 209 int iPass; 210 double lastObjective = 1.0e31; 211 // Create convexity rows for proposals 212 int numberMasterColumns = master.numberColumns(); 213 master.resize(numberMasterRows + numberBlocks, numberMasterColumns); 214 // Arrays to say which block and when created 215 int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks; 216 int * whichBlock = new int[maximumColumns]; 217 int * when = new int[maximumColumns]; 218 int numberColumnsGenerated = numberBlocks; 219 // fill in rhs and add in artificials 220 { 221 double * rowLower = master.rowLower(); 222 double * rowUpper = master.rowUpper(); 223 int iBlock; 224 columnAdd[0] = 0; 225 for (iBlock = 0; iBlock < numberBlocks; iBlock++) { 226 int iRow = iBlock + numberMasterRows;; 227 rowLower[iRow] = 1.0; 228 rowUpper[iRow] = 1.0; 229 rowAdd[iBlock] = iRow; 230 elementAdd[iBlock] = 1.0; 231 objective[iBlock] = 1.0e9; 232 columnAdd[iBlock+1] = iBlock + 1; 233 when[iBlock] = 1; 234 whichBlock[iBlock] = iBlock; 235 } 236 master.addColumns(numberBlocks, NULL, NULL, objective, 237 columnAdd, rowAdd, elementAdd); 238 } 239 // and resize matrix to double check clp will be happy 240 //master.matrix()>setDimensions(numberMasterRows+numberBlocks, 241 // numberMasterColumns+numberBlocks); 242 242 243 for (iPass=0;iPass<maxPass;iPass++) {244 printf("Start of pass %d\n",iPass);245 // Solve master  may be infeasible246 master.scaling(false);247 if (0) {248 master.writeMps("yy.mps");249 }243 for (iPass = 0; iPass < maxPass; iPass++) { 244 printf("Start of pass %d\n", iPass); 245 // Solve master  may be infeasible 246 master.scaling(false); 247 if (0) { 248 master.writeMps("yy.mps"); 249 } 250 250 251 master.primal();252 int problemStatus = master.status(); // do here as can change (delcols)253 if (master.numberIterations()==0&&iPass)254 break; // finished255 if (master.objectiveValue()>lastObjective1.0e7&&iPass>555)256 break; // finished257 lastObjective = master.objectiveValue();258 // mark basic ones and delete if necessary259 int iColumn;260 numberColumnsGenerated=master.numberColumns()numberMasterColumns;261 for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {262 if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)263 when[iColumn]=iPass;264 }265 if (numberColumnsGenerated+numberBlocks>maximumColumns) {266 // delete267 int numberKeep=0;268 int numberDelete=0;269 int * whichDelete = new int[numberColumns];270 for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {271 if (when[iColumn]>iPass7) {272 273 when[numberKeep]=when[iColumn];274 whichBlock[numberKeep++]=whichBlock[iColumn];275 276 277 whichDelete[numberDelete++]=iColumn+numberMasterColumns;278 279 }280 numberColumnsGenerated = numberDelete;281 master.deleteColumns(numberDelete,whichDelete);282 delete [] whichDelete;283 }284 const double * dual=NULL;285 bool deleteDual=false;286 if (problemStatus==0) {287 dual = master.dualRowSolution();288 } else if (problemStatus==1) {289 // could do composite objective290 dual = master.infeasibilityRay();291 deleteDual = true;292 printf("The sum of infeasibilities is %g\n",293 294 } else if (!master.numberColumns()) {295 assert(!iPass);296 dual = master.dualRowSolution();297 memset(master.dualRowSolution(),298 0,(numberMasterRows+numberBlocks)*sizeof(double));299 } else {300 abort();301 }302 // Create objective for sub problems and solve303 columnAdd[0]=0;304 int numberProposals=0;305 for (iBlock=0;iBlock<numberBlocks;iBlock++) {306 int numberColumns2 = sub[iBlock].numberColumns();307 double * saveObj = new double [numberColumns2];308 double * objective2 = sub[iBlock].objective();309 memcpy(saveObj,objective2,numberColumns2*sizeof(double));310 // new objective311 top[iBlock].transposeTimes(dual,objective2);312 int i;313 if (problemStatus==0) {314 for (i=0;i<numberColumns2;i++)315 objective2[i] = saveObj[i]objective2[i];316 } else {317 for (i=0;i<numberColumns2;i++)318 319 }320 sub[iBlock].primal();321 memcpy(objective2,saveObj,numberColumns2*sizeof(double));322 // get proposal323 if (sub[iBlock].numberIterations()!iPass) {324 double objValue =0.0;325 326 327 328 329 top[iBlock].times(solution,elementAdd+start);330 for (i=0;i<numberColumns2;i++)331 objValue += solution[i]*saveObj[i];332 333 int number=start;334 335 if (problemStatus) 336 dj=0.0;337 double smallest=1.0e100;338 double largest=0.0;339 for (i=0;i<numberMasterRows;i++) {340 341 if (fabs(value)>1.0e15) {342 dj = dual[i]*value;343 smallest = CoinMin(smallest,fabs(value));344 largest = CoinMax(largest,fabs(value));345 rowAdd[number]=i;346 elementAdd[number++]=value;347 348 349 350 351 rowAdd[number]=numberMasterRows+iBlock;352 elementAdd[number++]=1.0;353 354 355 356 iBlock,smallest,largest,dj);357 if (dj<1.0e6!iPass) {358 359 objective[numberProposals]=objValue;360 columnAdd[++numberProposals]=number;361 when[numberColumnsGenerated]=iPass;362 whichBlock[numberColumnsGenerated++]=iBlock;363 364 365 366 367 top[iBlock].times(solution,elementAdd+start);368 for (i=0;i<numberColumns2;i++)369 objValue += solution[i]*saveObj[i];370 371 int number=start;372 373 double smallest=1.0e100;374 double largest=0.0;375 for (i=0;i<numberMasterRows;i++) {376 377 if (fabs(value)>1.0e15) {378 dj = dual[i]*value;379 smallest = CoinMin(smallest,fabs(value));380 largest = CoinMax(largest,fabs(value));381 rowAdd[number]=i;382 elementAdd[number++]=value;383 384 385 386 387 388 iBlock,smallest,largest,dj);389 if (dj<1.0e6) {390 391 objective[numberProposals]=objValue;392 columnAdd[++numberProposals]=number;393 when[numberColumnsGenerated]=iPass;394 whichBlock[numberColumnsGenerated++]=iBlock;395 396 397 398 399 }400 delete [] saveObj;401 }402 if (deleteDual)403 delete [] dual;404 if (numberProposals)405 master.addColumns(numberProposals,NULL,NULL,objective,406 columnAdd,rowAdd,elementAdd);407 }408 // now put back a good solution409 double * lower = new double[numberMasterRows+numberBlocks];410 double * upper = new double[numberMasterRows+numberBlocks];411 numberColumnsGenerated += numberMasterColumns;412 double * sol = new double[numberColumnsGenerated];413 const double * solution = master.primalColumnSolution();414 const double * masterLower = master.rowLower();415 const double * masterUpper = master.rowUpper();416 double * fullSolution = model.primalColumnSolution();417 const double * fullLower = model.columnLower();418 const double * fullUpper = model.columnUpper();419 const double * rowSolution = master.primalRowSolution();420 double * fullRowSolution = model.primalRowSolution();421 int kRow=0;422 for (iRow=0;iRow<numberRows;iRow++) {423 if (rowBlock[iRow]==1) {424 model.setRowStatus(iRow,master.getRowStatus(kRow));425 fullRowSolution[iRow]=rowSolution[kRow++];426 }427 }428 int kColumn=0;429 for (iColumn=0;iColumn<numberColumns;iColumn++) {430 if (columnBlock[iColumn]==1) {431 model.setStatus(iColumn,master.getStatus(kColumn));432 fullSolution[iColumn]=solution[kColumn++];433 }434 }435 for (iBlock=0;iBlock<numberBlocks;iBlock++) {436 // convert top bit to by rows437 top[iBlock].reverseOrdering();438 // zero solution439 memset(sol,0,numberColumnsGenerated*sizeof(double));440 int i;441 for (i=numberMasterColumns;i<numberColumnsGenerated;i++)442 if (whichBlock[inumberMasterColumns]==iBlock)443 444 memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));445 master.times(1.0,sol,lower);446 for (iRow=0;iRow<numberMasterRows;iRow++) {447 double value = lower[iRow];448 if (masterUpper[iRow]<1.0e20)449 450 else451 upper[iRow]=COIN_DBL_MAX;452 if (masterLower[iRow]>1.0e20)453 454 else455 lower[iRow]=COIN_DBL_MAX;456 }457 sub[iBlock].addRows(numberMasterRows,lower,upper,458 459 460 461 462 sub[iBlock].primal();463 const double * subSolution = sub[iBlock].primalColumnSolution();464 const double * subRowSolution = sub[iBlock].primalRowSolution();465 // move solution466 kColumn=0;467 for (iColumn=0;iColumn<numberColumns;iColumn++) {468 if (columnBlock[iColumn]==iBlock) {469 model.setStatus(iColumn,sub[iBlock].getStatus(kColumn));470 fullSolution[iColumn]=subSolution[kColumn++];471 }472 }473 assert(kColumn==sub[iBlock].numberColumns());474 kRow=0;475 for (iRow=0;iRow<numberRows;iRow++) {476 if (rowBlock[iRow]==iBlock) {477 model.setRowStatus(iRow,sub[iBlock].getRowStatus(kRow));478 fullRowSolution[iRow]=subRowSolution[kRow++];479 }480 }481 assert(kRow == sub[iBlock].numberRows()numberMasterRows);482 }483 for (iColumn=0;iColumn<numberColumns;iColumn++) {484 if (fullSolution[iColumn]<fullUpper[iColumn]1.0e8&&485 fullSolution[iColumn]>fullLower[iColumn]+1.0e8) {486 assert(model.getStatus(iColumn)==ClpSimplex::basic);487 } else if (fullSolution[iColumn]>=fullUpper[iColumn]1.0e8) {488 // may help to make rest non basic489 model.setStatus(iColumn,ClpSimplex::atUpperBound);490 } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e8) {491 // may help to make rest non basic492 model.setStatus(iColumn,ClpSimplex::atLowerBound);493 }494 }495 for (iRow=0;iRow<numberRows;iRow++)496 model.setRowStatus(iRow,ClpSimplex::superBasic);497 model.primal(1);498 delete [] sol;499 delete [] lower;500 delete [] upper;501 delete [] whichBlock;502 delete [] when;503 delete [] columnAdd;504 delete [] rowAdd;505 delete [] elementAdd;506 delete [] objective;507 delete [] top;508 delete [] sub;509 delete [] rowBlock;510 delete [] columnBlock;511 return 0;512 } 251 master.primal(); 252 int problemStatus = master.status(); // do here as can change (delcols) 253 if (master.numberIterations() == 0 && iPass) 254 break; // finished 255 if (master.objectiveValue() > lastObjective  1.0e7 && iPass > 555) 256 break; // finished 257 lastObjective = master.objectiveValue(); 258 // mark basic ones and delete if necessary 259 int iColumn; 260 numberColumnsGenerated = master.numberColumns()  numberMasterColumns; 261 for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) { 262 if (master.getStatus(iColumn + numberMasterColumns) == ClpSimplex::basic) 263 when[iColumn] = iPass; 264 } 265 if (numberColumnsGenerated + numberBlocks > maximumColumns) { 266 // delete 267 int numberKeep = 0; 268 int numberDelete = 0; 269 int * whichDelete = new int[numberColumns]; 270 for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) { 271 if (when[iColumn] > iPass  7) { 272 // keep 273 when[numberKeep] = when[iColumn]; 274 whichBlock[numberKeep++] = whichBlock[iColumn]; 275 } else { 276 // delete 277 whichDelete[numberDelete++] = iColumn + numberMasterColumns; 278 } 279 } 280 numberColumnsGenerated = numberDelete; 281 master.deleteColumns(numberDelete, whichDelete); 282 delete [] whichDelete; 283 } 284 const double * dual = NULL; 285 bool deleteDual = false; 286 if (problemStatus == 0) { 287 dual = master.dualRowSolution(); 288 } else if (problemStatus == 1) { 289 // could do composite objective 290 dual = master.infeasibilityRay(); 291 deleteDual = true; 292 printf("The sum of infeasibilities is %g\n", 293 master.sumPrimalInfeasibilities()); 294 } else if (!master.numberColumns()) { 295 assert(!iPass); 296 dual = master.dualRowSolution(); 297 memset(master.dualRowSolution(), 298 0, (numberMasterRows + numberBlocks)*sizeof(double)); 299 } else { 300 abort(); 301 } 302 // Create objective for sub problems and solve 303 columnAdd[0] = 0; 304 int numberProposals = 0; 305 for (iBlock = 0; iBlock < numberBlocks; iBlock++) { 306 int numberColumns2 = sub[iBlock].numberColumns(); 307 double * saveObj = new double [numberColumns2]; 308 double * objective2 = sub[iBlock].objective(); 309 memcpy(saveObj, objective2, numberColumns2 * sizeof(double)); 310 // new objective 311 top[iBlock].transposeTimes(dual, objective2); 312 int i; 313 if (problemStatus == 0) { 314 for (i = 0; i < numberColumns2; i++) 315 objective2[i] = saveObj[i]  objective2[i]; 316 } else { 317 for (i = 0; i < numberColumns2; i++) 318 objective2[i] = objective2[i]; 319 } 320 sub[iBlock].primal(); 321 memcpy(objective2, saveObj, numberColumns2 * sizeof(double)); 322 // get proposal 323 if (sub[iBlock].numberIterations()  !iPass) { 324 double objValue = 0.0; 325 int start = columnAdd[numberProposals]; 326 // proposal 327 if (sub[iBlock].isProvenOptimal()) { 328 const double * solution = sub[iBlock].primalColumnSolution(); 329 top[iBlock].times(solution, elementAdd + start); 330 for (i = 0; i < numberColumns2; i++) 331 objValue += solution[i] * saveObj[i]; 332 // See if good dj and pack down 333 int number = start; 334 double dj = objValue; 335 if (problemStatus) 336 dj = 0.0; 337 double smallest = 1.0e100; 338 double largest = 0.0; 339 for (i = 0; i < numberMasterRows; i++) { 340 double value = elementAdd[start+i]; 341 if (fabs(value) > 1.0e15) { 342 dj = dual[i] * value; 343 smallest = CoinMin(smallest, fabs(value)); 344 largest = CoinMax(largest, fabs(value)); 345 rowAdd[number] = i; 346 elementAdd[number++] = value; 347 } 348 } 349 // and convexity 350 dj = dual[numberMasterRows+iBlock]; 351 rowAdd[number] = numberMasterRows + iBlock; 352 elementAdd[number++] = 1.0; 353 // if elements large then scale? 354 //if (largest>1.0e8smallest<1.0e8) 355 printf("For subproblem %d smallest  %g, largest %g  dj %g\n", 356 iBlock, smallest, largest, dj); 357 if (dj < 1.0e6  !iPass) { 358 // take 359 objective[numberProposals] = objValue; 360 columnAdd[++numberProposals] = number; 361 when[numberColumnsGenerated] = iPass; 362 whichBlock[numberColumnsGenerated++] = iBlock; 363 } 364 } else if (sub[iBlock].isProvenDualInfeasible()) { 365 // use ray 366 const double * solution = sub[iBlock].unboundedRay(); 367 top[iBlock].times(solution, elementAdd + start); 368 for (i = 0; i < numberColumns2; i++) 369 objValue += solution[i] * saveObj[i]; 370 // See if good dj and pack down 371 int number = start; 372 double dj = objValue; 373 double smallest = 1.0e100; 374 double largest = 0.0; 375 for (i = 0; i < numberMasterRows; i++) { 376 double value = elementAdd[start+i]; 377 if (fabs(value) > 1.0e15) { 378 dj = dual[i] * value; 379 smallest = CoinMin(smallest, fabs(value)); 380 largest = CoinMax(largest, fabs(value)); 381 rowAdd[number] = i; 382 elementAdd[number++] = value; 383 } 384 } 385 // if elements large or small then scale? 386 //if (largest>1.0e8smallest<1.0e8) 387 printf("For subproblem ray %d smallest  %g, largest %g  dj %g\n", 388 iBlock, smallest, largest, dj); 389 if (dj < 1.0e6) { 390 // take 391 objective[numberProposals] = objValue; 392 columnAdd[++numberProposals] = number; 393 when[numberColumnsGenerated] = iPass; 394 whichBlock[numberColumnsGenerated++] = iBlock; 395 } 396 } else { 397 abort(); 398 } 399 } 400 delete [] saveObj; 401 } 402 if (deleteDual) 403 delete [] dual; 404 if (numberProposals) 405 master.addColumns(numberProposals, NULL, NULL, objective, 406 columnAdd, rowAdd, elementAdd); 407 } 408 // now put back a good solution 409 double * lower = new double[numberMasterRows+numberBlocks]; 410 double * upper = new double[numberMasterRows+numberBlocks]; 411 numberColumnsGenerated += numberMasterColumns; 412 double * sol = new double[numberColumnsGenerated]; 413 const double * solution = master.primalColumnSolution(); 414 const double * masterLower = master.rowLower(); 415 const double * masterUpper = master.rowUpper(); 416 double * fullSolution = model.primalColumnSolution(); 417 const double * fullLower = model.columnLower(); 418 const double * fullUpper = model.columnUpper(); 419 const double * rowSolution = master.primalRowSolution(); 420 double * fullRowSolution = model.primalRowSolution(); 421 int kRow = 0; 422 for (iRow = 0; iRow < numberRows; iRow++) { 423 if (rowBlock[iRow] == 1) { 424 model.setRowStatus(iRow, master.getRowStatus(kRow)); 425 fullRowSolution[iRow] = rowSolution[kRow++]; 426 } 427 } 428 int kColumn = 0; 429 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 430 if (columnBlock[iColumn] == 1) { 431 model.setStatus(iColumn, master.getStatus(kColumn)); 432 fullSolution[iColumn] = solution[kColumn++]; 433 } 434 } 435 for (iBlock = 0; iBlock < numberBlocks; iBlock++) { 436 // convert top bit to by rows 437 top[iBlock].reverseOrdering(); 438 // zero solution 439 memset(sol, 0, numberColumnsGenerated * sizeof(double)); 440 int i; 441 for (i = numberMasterColumns; i < numberColumnsGenerated; i++) 442 if (whichBlock[inumberMasterColumns] == iBlock) 443 sol[i] = solution[i]; 444 memset(lower, 0, (numberMasterRows + numberBlocks)*sizeof(double)); 445 master.times(1.0, sol, lower); 446 for (iRow = 0; iRow < numberMasterRows; iRow++) { 447 double value = lower[iRow]; 448 if (masterUpper[iRow] < 1.0e20) 449 upper[iRow] = value; 450 else 451 upper[iRow] = COIN_DBL_MAX; 452 if (masterLower[iRow] > 1.0e20) 453 lower[iRow] = value; 454 else 455 lower[iRow] = COIN_DBL_MAX; 456 } 457 sub[iBlock].addRows(numberMasterRows, lower, upper, 458 top[iBlock].getVectorStarts(), 459 top[iBlock].getVectorLengths(), 460 top[iBlock].getIndices(), 461 top[iBlock].getElements()); 462 sub[iBlock].primal(); 463 const double * subSolution = sub[iBlock].primalColumnSolution(); 464 const double * subRowSolution = sub[iBlock].primalRowSolution(); 465 // move solution 466 kColumn = 0; 467 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 468 if (columnBlock[iColumn] == iBlock) { 469 model.setStatus(iColumn, sub[iBlock].getStatus(kColumn)); 470 fullSolution[iColumn] = subSolution[kColumn++]; 471 } 472 } 473 assert(kColumn == sub[iBlock].numberColumns()); 474 kRow = 0; 475 for (iRow = 0; iRow < numberRows; iRow++) { 476 if (rowBlock[iRow] == iBlock) { 477 model.setRowStatus(iRow, sub[iBlock].getRowStatus(kRow)); 478 fullRowSolution[iRow] = subRowSolution[kRow++]; 479 } 480 } 481 assert(kRow == sub[iBlock].numberRows()  numberMasterRows); 482 } 483 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 484 if (fullSolution[iColumn] < fullUpper[iColumn]  1.0e8 && 485 fullSolution[iColumn] > fullLower[iColumn] + 1.0e8) { 486 assert(model.getStatus(iColumn) == ClpSimplex::basic); 487 } else if (fullSolution[iColumn] >= fullUpper[iColumn]  1.0e8) { 488 // may help to make rest non basic 489 model.setStatus(iColumn, ClpSimplex::atUpperBound); 490 } else if (fullSolution[iColumn] <= fullLower[iColumn] + 1.0e8) { 491 // may help to make rest non basic 492 model.setStatus(iColumn, ClpSimplex::atLowerBound); 493 } 494 } 495 for (iRow = 0; iRow < numberRows; iRow++) 496 model.setRowStatus(iRow, ClpSimplex::superBasic); 497 model.primal(1); 498 delete [] sol; 499 delete [] lower; 500 delete [] upper; 501 delete [] whichBlock; 502 delete [] when; 503 delete [] columnAdd; 504 delete [] rowAdd; 505 delete [] elementAdd; 506 delete [] objective; 507 delete [] top; 508 delete [] sub; 509 delete [] rowBlock; 510 delete [] columnBlock; 511 return 0; 512 }
Note: See TracChangeset
for help on using the changeset viewer.