1 | /* $Id: CoinStructuredModel.cpp 1373 2011-01-03 23:57:44Z lou $ */ |
2 | // Copyright (C) 2008, International Business Machines |
3 | // Corporation and others. All Rights Reserved. |
4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
5 | |
6 | |
7 | #include "CoinUtilsConfig.h" |
8 | #include "CoinHelperFunctions.hpp" |
9 | #include "CoinStructuredModel.hpp" |
10 | #include "CoinSort.hpp" |
11 | #include "CoinMpsIO.hpp" |
12 | #include "CoinFloatEqual.hpp" |
13 | |
14 | //############################################################################# |
15 | // Constructors / Destructor / Assignment |
16 | //############################################################################# |
17 | |
18 | //------------------------------------------------------------------- |
19 | // Default Constructor |
20 | //------------------------------------------------------------------- |
21 | CoinStructuredModel::CoinStructuredModel () |
22 | : CoinBaseModel(), |
23 | numberRowBlocks_(0), |
24 | numberColumnBlocks_(0), |
25 | numberElementBlocks_(0), |
26 | maximumElementBlocks_(0), |
27 | blocks_(NULL), |
28 | coinModelBlocks_(NULL), |
29 | blockType_(NULL) |
30 | { |
31 | } |
32 | /* Read a problem in MPS or GAMS format from the given filename. |
33 | */ |
34 | CoinStructuredModel::CoinStructuredModel(const char *fileName, |
35 | int decomposeType, |
36 | int maxBlocks) |
37 | : CoinBaseModel(), |
38 | numberRowBlocks_(0), |
39 | numberColumnBlocks_(0), |
40 | numberElementBlocks_(0), |
41 | maximumElementBlocks_(0), |
42 | blocks_(NULL), |
43 | coinModelBlocks_(NULL), |
44 | blockType_(NULL) |
45 | { |
46 | CoinModel coinModel(fileName,false); |
47 | if (coinModel.numberRows()) { |
48 | problemName_ = coinModel.getProblemName(); |
49 | optimizationDirection_ = coinModel.optimizationDirection(); |
50 | objectiveOffset_ = coinModel.objectiveOffset(); |
51 | if (!decomposeType) { |
52 | addBlock("row_master" ,"column_master" ,coinModel); |
53 | } else { |
54 | const CoinPackedMatrix * matrix = coinModel.packedMatrix(); |
55 | if (!matrix) |
56 | coinModel.convertMatrix(); |
57 | decompose(coinModel,decomposeType,maxBlocks); |
58 | //addBlock("row_master","column_master",coinModel); |
59 | } |
60 | } |
61 | } |
62 | |
63 | //------------------------------------------------------------------- |
64 | // Copy constructor |
65 | //------------------------------------------------------------------- |
66 | CoinStructuredModel::CoinStructuredModel (const CoinStructuredModel & rhs) |
67 | : CoinBaseModel(rhs), |
68 | numberRowBlocks_(rhs.numberRowBlocks_), |
69 | numberColumnBlocks_(rhs.numberColumnBlocks_), |
70 | numberElementBlocks_(rhs.numberElementBlocks_), |
71 | maximumElementBlocks_(rhs.maximumElementBlocks_) |
72 | { |
73 | if (maximumElementBlocks_) { |
74 | blocks_ = CoinCopyOfArray(rhs.blocks_,maximumElementBlocks_); |
75 | for (int i=0;i<numberElementBlocks_;i++) |
76 | blocks_[i]= rhs.blocks_[i]->clone(); |
77 | blockType_ = CoinCopyOfArray(rhs.blockType_,maximumElementBlocks_); |
78 | if (rhs.coinModelBlocks_) { |
79 | coinModelBlocks_ = CoinCopyOfArray(rhs.coinModelBlocks_, |
80 | maximumElementBlocks_); |
81 | for (int i=0;i<numberElementBlocks_;i++) |
82 | coinModelBlocks_[i]= new CoinModel(*rhs.coinModelBlocks_[i]); |
83 | } else { |
84 | coinModelBlocks_ = NULL; |
85 | } |
86 | } else { |
87 | blocks_ = NULL; |
88 | blockType_ = NULL; |
89 | coinModelBlocks_ = NULL; |
90 | } |
91 | rowBlockNames_ = rhs.rowBlockNames_; |
92 | columnBlockNames_ = rhs.columnBlockNames_; |
93 | } |
94 | |
95 | //------------------------------------------------------------------- |
96 | // Destructor |
97 | //------------------------------------------------------------------- |
98 | CoinStructuredModel::~CoinStructuredModel () |
99 | { |
100 | for (int i=0;i<numberElementBlocks_;i++) |
101 | delete blocks_[i]; |
102 | delete [] blocks_; |
103 | delete [] blockType_; |
104 | if (coinModelBlocks_) { |
105 | for (int i=0;i<numberElementBlocks_;i++) |
106 | delete coinModelBlocks_[i]; |
107 | delete [] coinModelBlocks_; |
108 | } |
109 | } |
110 | |
111 | // Clone |
112 | CoinBaseModel * |
113 | CoinStructuredModel::clone() const |
114 | { |
115 | return new CoinStructuredModel(*this); |
116 | } |
117 | |
118 | //---------------------------------------------------------------- |
119 | // Assignment operator |
120 | //------------------------------------------------------------------- |
121 | CoinStructuredModel & |
122 | CoinStructuredModel::operator=(const CoinStructuredModel& rhs) |
123 | { |
124 | if (this != &rhs) { |
125 | CoinBaseModel::operator=(rhs); |
126 | for (int i=0;i<numberElementBlocks_;i++) |
127 | delete blocks_[i]; |
128 | delete [] blocks_; |
129 | delete [] blockType_; |
130 | if (coinModelBlocks_) { |
131 | for (int i=0;i<numberElementBlocks_;i++) |
132 | delete coinModelBlocks_[i]; |
133 | delete [] coinModelBlocks_; |
134 | } |
135 | numberRowBlocks_ = rhs.numberRowBlocks_; |
136 | numberColumnBlocks_ = rhs.numberColumnBlocks_; |
137 | numberElementBlocks_ = rhs.numberElementBlocks_; |
138 | maximumElementBlocks_ = rhs.maximumElementBlocks_; |
139 | if (maximumElementBlocks_) { |
140 | blocks_ = CoinCopyOfArray(rhs.blocks_,maximumElementBlocks_); |
141 | for (int i=0;i<numberElementBlocks_;i++) |
142 | blocks_[i]= rhs.blocks_[i]->clone(); |
143 | blockType_ = CoinCopyOfArray(rhs.blockType_,maximumElementBlocks_); |
144 | if (rhs.coinModelBlocks_) { |
145 | coinModelBlocks_ = CoinCopyOfArray(rhs.coinModelBlocks_, |
146 | maximumElementBlocks_); |
147 | for (int i=0;i<numberElementBlocks_;i++) |
148 | coinModelBlocks_[i]= new CoinModel(*rhs.coinModelBlocks_[i]); |
149 | } else { |
150 | coinModelBlocks_ = NULL; |
151 | } |
152 | } else { |
153 | blocks_ = NULL; |
154 | blockType_ = NULL; |
155 | coinModelBlocks_ = NULL; |
156 | } |
157 | rowBlockNames_ = rhs.rowBlockNames_; |
158 | columnBlockNames_ = rhs.columnBlockNames_; |
159 | } |
160 | return *this; |
161 | } |
162 | static bool sameValues(const double * a, const double *b, int n) |
163 | { |
164 | int i; |
165 | for ( i=0;i<n;i++) { |
166 | if (a[i]!=b[i]) |
167 | break; |
168 | } |
169 | return (i==n); |
170 | } |
171 | static bool sameValues(const int * a, const int *b, int n) |
172 | { |
173 | int i; |
174 | for ( i=0;i<n;i++) { |
175 | if (a[i]!=b[i]) |
176 | break; |
177 | } |
178 | return (i==n); |
179 | } |
180 | static bool sameValues(const CoinModel * a, const CoinModel *b, bool doRows) |
181 | { |
182 | int i=0; |
183 | int n=0; |
184 | if (doRows) { |
185 | n = a->numberRows(); |
186 | for ( i=0;i<n;i++) { |
187 | const char * aName = a->getRowName(i); |
188 | const char * bName = b->getRowName(i); |
189 | bool good=true; |
190 | if (aName) { |
191 | if (!bName||strcmp(aName,bName)) |
192 | good=false; |
193 | } else if (bName) { |
194 | good=false; |
195 | } |
196 | if (!good) |
197 | break; |
198 | } |
199 | } else { |
200 | n = a->numberColumns(); |
201 | for ( i=0;i<n;i++) { |
202 | const char * aName = a->getColumnName(i); |
203 | const char * bName = b->getColumnName(i); |
204 | bool good=true; |
205 | if (aName) { |
206 | if (!bName||strcmp(aName,bName)) |
207 | good=false; |
208 | } else if (bName) { |
209 | good=false; |
210 | } |
211 | if (!good) |
212 | break; |
213 | } |
214 | } |
215 | return (i==n); |
216 | } |
217 | // Add a row block name and number of rows |
218 | int |
219 | CoinStructuredModel::addRowBlock(int numberRows,const std::string &name) |
220 | { |
221 | int iRowBlock; |
222 | for (iRowBlock=0;iRowBlock<numberRowBlocks_;iRowBlock++) { |
223 | if (name==rowBlockNames_[iRowBlock]) |
224 | break; |
225 | } |
226 | if (iRowBlock==numberRowBlocks_) { |
227 | rowBlockNames_.push_back(name); |
228 | numberRowBlocks_++; |
229 | numberRows_ += numberRows; |
230 | } |
231 | return iRowBlock; |
232 | } |
233 | // Return a row block index given a row block name |
234 | int |
235 | CoinStructuredModel::rowBlock(const std::string &name) const |
236 | { |
237 | int iRowBlock; |
238 | for (iRowBlock=0;iRowBlock<numberRowBlocks_;iRowBlock++) { |
239 | if (name==rowBlockNames_[iRowBlock]) |
240 | break; |
241 | } |
242 | if (iRowBlock==numberRowBlocks_) |
243 | iRowBlock=-1; |
244 | return iRowBlock; |
245 | } |
246 | // Add a column block name and number of columns |
247 | int |
248 | CoinStructuredModel::addColumnBlock(int numberColumns,const std::string &name) |
249 | { |
250 | int iColumnBlock; |
251 | for (iColumnBlock=0;iColumnBlock<numberColumnBlocks_;iColumnBlock++) { |
252 | if (name==columnBlockNames_[iColumnBlock]) |
253 | break; |
254 | } |
255 | if (iColumnBlock==numberColumnBlocks_) { |
256 | columnBlockNames_.push_back(name); |
257 | numberColumnBlocks_++; |
258 | numberColumns_ += numberColumns; |
259 | } |
260 | return iColumnBlock; |
261 | } |
262 | // Return a column block index given a column block name |
263 | int |
264 | CoinStructuredModel::columnBlock(const std::string &name) const |
265 | { |
266 | int iColumnBlock; |
267 | for (iColumnBlock=0;iColumnBlock<numberColumnBlocks_;iColumnBlock++) { |
268 | if (name==columnBlockNames_[iColumnBlock]) |
269 | break; |
270 | } |
271 | if (iColumnBlock==numberColumnBlocks_) |
272 | iColumnBlock=-1; |
273 | return iColumnBlock; |
274 | } |
275 | // Return number of elements |
276 | CoinBigIndex |
277 | CoinStructuredModel::numberElements() const |
278 | { |
279 | CoinBigIndex numberElements=0; |
280 | for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) { |
281 | numberElements += blocks_[iBlock]->numberElements(); |
282 | } |
283 | return numberElements; |
284 | } |
285 | // Return i'th block as CoinModel (or NULL) |
286 | CoinModel * |
287 | CoinStructuredModel::coinBlock(int i) const |
288 | { |
289 | CoinModel * block = dynamic_cast<CoinModel *>(blocks_[i]); |
290 | if (block) |
291 | return block; |
292 | else if (coinModelBlocks_) |
293 | return coinModelBlocks_[i]; |
294 | else |
295 | return NULL; |
296 | } |
297 | |
298 | /* Return block as a CoinModel block |
299 | and fill in info structure and update counts |
300 | */ |
301 | CoinModel * |
302 | CoinStructuredModel::coinModelBlock(CoinModelBlockInfo & info) |
303 | { |
304 | // CoinStructuredModel |
305 | int numberRows=this->numberRows(); |
306 | int numberColumns =this-> numberColumns(); |
307 | int numberRowBlocks=this->numberRowBlocks(); |
308 | int numberColumnBlocks =this-> numberColumnBlocks(); |
309 | int numberElementBlocks =this-> numberElementBlocks(); |
310 | CoinBigIndex numberElements=this->numberElements(); |
311 | // See what is needed |
312 | double * rowLower = NULL; |
313 | double * rowUpper = NULL; |
314 | double * columnLower = NULL; |
315 | double * columnUpper = NULL; |
316 | double * objective = NULL; |
317 | int * integerType = NULL; |
318 | info = CoinModelBlockInfo(); |
319 | CoinModel ** blocks = new CoinModel * [numberElementBlocks]; |
320 | for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) { |
321 | CoinModelBlockInfo thisInfo=blockType_[iBlock]; |
322 | CoinStructuredModel * subModel = dynamic_cast<CoinStructuredModel *>(blocks_[iBlock]); |
323 | CoinModel * thisBlock; |
324 | if (subModel) { |
325 | thisBlock = subModel->coinModelBlock(thisInfo); |
326 | fillInfo(thisInfo,subModel); |
327 | setCoinModel(thisBlock,iBlock); |
328 | } else { |
329 | thisBlock = dynamic_cast<CoinModel *>(blocks_[iBlock]); |
330 | assert (thisBlock); |
331 | fillInfo(thisInfo,thisBlock); |
332 | } |
333 | blocks[iBlock]=thisBlock; |
334 | if (thisInfo.rhs&&!info.rhs) { |
335 | info.rhs=1; |
336 | rowLower = new double [numberRows]; |
337 | rowUpper = new double [numberRows]; |
338 | CoinFillN(rowLower,numberRows,-COIN_DBL_MAX); |
339 | CoinFillN(rowUpper,numberRows,COIN_DBL_MAX); |
340 | } |
341 | if (thisInfo.bounds&&!info.bounds) { |
342 | info.bounds=1; |
343 | columnLower = new double [numberColumns]; |
344 | columnUpper = new double [numberColumns]; |
345 | objective = new double [numberColumns]; |
346 | CoinFillN(columnLower,numberColumns,0.0); |
347 | CoinFillN(columnUpper,numberColumns,COIN_DBL_MAX); |
348 | CoinFillN(objective,numberColumns,0.0); |
349 | } |
350 | if (thisInfo.integer&&!info.integer) { |
351 | info.integer=1; |
352 | integerType = new int [numberColumns]; |
353 | CoinFillN(integerType,numberColumns,0); |
354 | } |
355 | if (thisInfo.rowName&&!info.rowName) { |
356 | info.rowName=1; |
357 | } |
358 | if (thisInfo.columnName&&!info.columnName) { |
359 | info.columnName=1; |
360 | } |
361 | } |
362 | // Space for elements |
363 | int * row = new int[numberElements]; |
364 | int * column = new int[numberElements]; |
365 | double * element = new double[numberElements]; |
366 | numberElements=0; |
367 | // Bases for blocks |
368 | int * rowBase = new int[numberRowBlocks]; |
369 | CoinFillN(rowBase,numberRowBlocks,-1); |
370 | CoinModelBlockInfo * rowBlockInfo = |
371 | new CoinModelBlockInfo [numberRowBlocks]; |
372 | int * columnBase = new int[numberColumnBlocks]; |
373 | CoinFillN(columnBase,numberColumnBlocks,-1); |
374 | CoinModelBlockInfo * columnBlockInfo = |
375 | new CoinModelBlockInfo [numberColumnBlocks]; |
376 | for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) { |
377 | int iRowBlock = rowBlock(blocks[iBlock]->getRowBlock()); |
378 | assert (iRowBlock>=0&&iRowBlock<numberRowBlocks); |
379 | if (rowBase[iRowBlock]==-1) |
380 | rowBase[iRowBlock]=blocks[iBlock]->numberRows(); |
381 | else |
382 | assert (rowBase[iRowBlock]==blocks[iBlock]->numberRows()); |
383 | int iColumnBlock = columnBlock(blocks[iBlock]->getColumnBlock()); |
384 | assert (iColumnBlock>=0&&iColumnBlock<numberColumnBlocks); |
385 | if (columnBase[iColumnBlock]==-1) |
386 | columnBase[iColumnBlock]=blocks[iBlock]->numberColumns(); |
387 | else |
388 | assert (columnBase[iColumnBlock]==blocks[iBlock]->numberColumns()); |
389 | } |
390 | int n=0; |
391 | for (int iBlock=0;iBlock<numberRowBlocks;iBlock++) { |
392 | int k = rowBase[iBlock]; |
393 | rowBase[iBlock]=n; |
394 | assert (k>=0); |
395 | n+=k; |
396 | } |
397 | assert (n==numberRows); |
398 | n=0; |
399 | for (int iBlock=0;iBlock<numberColumnBlocks;iBlock++) { |
400 | int k = columnBase[iBlock]; |
401 | columnBase[iBlock]=n; |
402 | assert (k>=0); |
403 | n+=k; |
404 | } |
405 | assert (n==numberColumns); |
406 | for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) { |
407 | CoinModelBlockInfo thisInfo=blockType_[iBlock]; |
408 | CoinModel * thisBlock = blocks[iBlock]; |
409 | int iRowBlock = rowBlock(blocks[iBlock]->getRowBlock()); |
410 | int iRowBase = rowBase[iRowBlock]; |
411 | int nRows = thisBlock->numberRows(); |
412 | // could check if identical before error |
413 | if (thisInfo.rhs) { |
414 | assert (!rowBlockInfo[iRowBlock].rhs); |
415 | rowBlockInfo[iRowBlock].rhs=1; |
416 | memcpy(rowLower+iRowBase,thisBlock->rowLowerArray(), |
417 | nRows*sizeof(double)); |
418 | memcpy(rowUpper+iRowBase,thisBlock->rowUpperArray(), |
419 | nRows*sizeof(double)); |
420 | } |
421 | int iColumnBlock = columnBlock(blocks[iBlock]->getColumnBlock()); |
422 | int iColumnBase = columnBase[iColumnBlock]; |
423 | int nColumns = thisBlock->numberColumns(); |
424 | if (thisInfo.bounds) { |
425 | assert (!columnBlockInfo[iColumnBlock].bounds); |
426 | columnBlockInfo[iColumnBlock].bounds=1; |
427 | memcpy(columnLower+iColumnBase,thisBlock->columnLowerArray(), |
428 | nColumns*sizeof(double)); |
429 | memcpy(columnUpper+iColumnBase,thisBlock->columnUpperArray(), |
430 | nColumns*sizeof(double)); |
431 | memcpy(objective+iColumnBase,thisBlock->objectiveArray(), |
432 | nColumns*sizeof(double)); |
433 | } |
434 | if (thisInfo.integer) { |
435 | assert (!columnBlockInfo[iColumnBlock].integer); |
436 | columnBlockInfo[iColumnBlock].integer=1; |
437 | memcpy(integerType+iColumnBase,thisBlock->integerTypeArray(), |
438 | nColumns*sizeof(int)); |
439 | } |
440 | const CoinPackedMatrix * elementBlock = thisBlock->packedMatrix(); |
441 | // get matrix data pointers |
442 | const int * row2 = elementBlock->getIndices(); |
443 | const CoinBigIndex * columnStart = elementBlock->getVectorStarts(); |
444 | const double * elementByColumn = elementBlock->getElements(); |
445 | const int * columnLength = elementBlock->getVectorLengths(); |
446 | int n = elementBlock->getNumCols(); |
447 | assert (elementBlock->isColOrdered()); |
448 | for (int iColumn=0;iColumn<n;iColumn++) { |
449 | CoinBigIndex j; |
450 | int jColumn = iColumn+iColumnBase; |
451 | for (j=columnStart[iColumn]; |
452 | j<columnStart[iColumn]+columnLength[iColumn];j++) { |
453 | row[numberElements]=row2[j]+iRowBase; |
454 | column[numberElements]=jColumn; |
455 | element[numberElements++]=elementByColumn[j]; |
456 | } |
457 | } |
458 | } |
459 | delete [] rowBlockInfo; |
460 | delete [] columnBlockInfo; |
461 | CoinPackedMatrix matrix(true,row,column,element,numberElements); |
462 | if (numberElements) |
463 | info.matrix=1; |
464 | delete [] row; |
465 | delete [] column; |
466 | delete [] element; |
467 | CoinModel * block = new CoinModel(numberRows, numberColumns, |
468 | &matrix, |
469 | rowLower, rowUpper, |
470 | columnLower, columnUpper, objective); |
471 | delete [] rowLower; |
472 | delete [] rowUpper; |
473 | delete [] columnLower; |
474 | delete [] columnUpper; |
475 | delete [] objective; |
476 | // Do integers if wanted |
477 | if (integerType) { |
478 | for (int iColumn=0;iColumn<numberColumns;iColumn++) { |
479 | block->setColumnIsInteger(iColumn,integerType[iColumn]!=0); |
480 | } |
481 | } |
482 | delete [] integerType; |
483 | block->setObjectiveOffset(objectiveOffset()); |
484 | if (info.rowName||info.columnName) { |
485 | for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) { |
486 | CoinModelBlockInfo thisInfo; |
487 | CoinModel * thisBlock = blocks[iBlock]; |
488 | int iRowBlock = rowBlock(thisBlock->getRowBlock()); |
489 | int iRowBase = rowBase[iRowBlock]; |
490 | if (thisInfo.rowName) { |
491 | int numberItems = thisBlock->rowNames()->numberItems(); |
492 | assert( thisBlock->numberRows()>=numberItems); |
493 | if (numberItems) { |
494 | const char *const * rowNames=thisBlock->rowNames()->names(); |
495 | for (int i=0;i<numberItems;i++) { |
496 | std::string name = rowNames[i]; |
497 | block->setRowName(i+iRowBase,name.c_str()); |
498 | } |
499 | } |
500 | } |
501 | int iColumnBlock = columnBlock(thisBlock->getColumnBlock()); |
502 | int iColumnBase = columnBase[iColumnBlock]; |
503 | if (thisInfo.columnName) { |
504 | int numberItems = thisBlock->columnNames()->numberItems(); |
505 | assert( thisBlock->numberColumns()>=numberItems); |
506 | if (numberItems) { |
507 | const char *const * columnNames=thisBlock->columnNames()->names(); |
508 | for (int i=0;i<numberItems;i++) { |
509 | std::string name = columnNames[i]; |
510 | block->setColumnName(i+iColumnBase,name.c_str()); |
511 | } |
512 | } |
513 | } |
514 | } |
515 | } |
516 | delete [] rowBase; |
517 | delete [] columnBase; |
518 | for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) { |
519 | if (static_cast<CoinBaseModel *> (blocks[iBlock])!= |
520 | static_cast<CoinBaseModel *> (blocks_[iBlock])) |
521 | delete blocks[iBlock]; |
522 | } |
523 | delete [] blocks; |
524 | return block; |
525 | } |
526 | // Sets given block into coinModelBlocks_ |
527 | void |
528 | CoinStructuredModel::setCoinModel(CoinModel * block, int iBlock) |
529 | { |
530 | if (!coinModelBlocks_) { |
531 | coinModelBlocks_ = new CoinModel * [maximumElementBlocks_]; |
532 | CoinZeroN(coinModelBlocks_,maximumElementBlocks_); |
533 | } |
534 | delete coinModelBlocks_[iBlock]; |
535 | coinModelBlocks_[iBlock]=block; |
536 | } |
537 | // Refresh info in blockType_ |
538 | void |
539 | CoinStructuredModel::refresh(int iBlock) |
540 | { |
541 | fillInfo(blockType_[iBlock],coinBlock(iBlock)); |
542 | } |
543 | /* Fill in info structure and update counts |
544 | Returns number of inconsistencies on border |
545 | */ |
546 | int |
547 | CoinStructuredModel::fillInfo(CoinModelBlockInfo & info, |
548 | const CoinModel * block) |
549 | { |
550 | int whatsSet = block->whatIsSet(); |
551 | info.matrix = static_cast<char>(((whatsSet&1)!=0) ? 1 : 0); |
552 | info.rhs = static_cast<char>(((whatsSet&2)!=0) ? 1 : 0); |
553 | info.rowName = static_cast<char>(((whatsSet&4)!=0) ? 1 : 0); |
554 | info.integer = static_cast<char>(((whatsSet&32)!=0) ? 1 : 0); |
555 | info.bounds = static_cast<char>(((whatsSet&8)!=0) ? 1 : 0); |
556 | info.columnName = static_cast<char>(((whatsSet&16)!=0) ? 1 : 0); |
557 | int numberRows = block->numberRows(); |
558 | int numberColumns = block->numberColumns(); |
559 | // Which block |
560 | int iRowBlock=addRowBlock(numberRows,block->getRowBlock()); |
561 | info.rowBlock=iRowBlock; |
562 | int iColumnBlock=addColumnBlock(numberColumns,block->getColumnBlock()); |
563 | info.columnBlock=iColumnBlock; |
564 | int numberErrors=0; |
565 | CoinModelBlockInfo sumInfo=blockType_[numberElementBlocks_-1]; |
566 | int iRhs=(sumInfo.rhs) ? numberElementBlocks_-1 : -1; |
567 | int iRowName=(sumInfo.rowName) ? numberElementBlocks_-1 : -1; |
568 | int iBounds=(sumInfo.bounds) ? numberElementBlocks_-1 : -1; |
569 | int iColumnName=(sumInfo.columnName) ? numberElementBlocks_-1 : -1; |
570 | int iInteger=(sumInfo.integer) ? numberElementBlocks_-1 : -1; |
571 | for (int i=0;i<numberElementBlocks_-1;i++) { |
572 | if (iRowBlock==blockType_[i].rowBlock) { |
573 | if (numberRows!=blocks_[i]->numberRows()) |
574 | numberErrors+=1000; |
575 | if (blockType_[i].rhs) { |
576 | if (iRhs<0) { |
577 | iRhs=i; |
578 | } else { |
579 | // check |
580 | const double * a = static_cast<CoinModel *>(blocks_[iRhs])->rowLowerArray(); |
581 | const double * b = static_cast<CoinModel *>(blocks_[i])->rowLowerArray(); |
582 | if (!sameValues(a,b,numberRows)) |
583 | numberErrors++; |
584 | a = static_cast<CoinModel *>(blocks_[iRhs])->rowUpperArray(); |
585 | b = static_cast<CoinModel *>(blocks_[i])->rowUpperArray(); |
586 | if (!sameValues(a,b,numberRows)) |
587 | numberErrors++; |
588 | } |
589 | } |
590 | if (blockType_[i].rowName) { |
591 | if (iRowName<0) { |
592 | iRowName=i; |
593 | } else { |
594 | // check |
595 | if (!sameValues(static_cast<CoinModel *>(blocks_[iRowName]), |
596 | static_cast<CoinModel *>(blocks_[i]),true)) |
597 | numberErrors++; |
598 | } |
599 | } |
600 | } |
601 | if (iColumnBlock==blockType_[i].columnBlock) { |
602 | if (numberColumns!=blocks_[i]->numberColumns()) |
603 | numberErrors+=1000; |
604 | if (blockType_[i].bounds) { |
605 | if (iBounds<0) { |
606 | iBounds=i; |
607 | } else { |
608 | // check |
609 | const double * a = static_cast<CoinModel *>(blocks_[iBounds])->columnLowerArray(); |
610 | const double * b = static_cast<CoinModel *>(blocks_[i])->columnLowerArray(); |
611 | if (!sameValues(a,b,numberColumns)) |
612 | numberErrors++; |
613 | a = static_cast<CoinModel *>(blocks_[iBounds])->columnUpperArray(); |
614 | b = static_cast<CoinModel *>(blocks_[i])->columnUpperArray(); |
615 | if (!sameValues(a,b,numberColumns)) |
616 | numberErrors++; |
617 | a = static_cast<CoinModel *>(blocks_[iBounds])->objectiveArray(); |
618 | b = static_cast<CoinModel *>(blocks_[i])->objectiveArray(); |
619 | if (!sameValues(a,b,numberColumns)) |
620 | numberErrors++; |
621 | } |
622 | } |
623 | if (blockType_[i].columnName) { |
624 | if (iColumnName<0) { |
625 | iColumnName=i; |
626 | } else { |
627 | // check |
628 | if (!sameValues(static_cast<CoinModel *>(blocks_[iColumnName]), |
629 | static_cast<CoinModel *>(blocks_[i]),false)) |
630 | numberErrors++; |
631 | } |
632 | } |
633 | if (blockType_[i].integer) { |
634 | if (iInteger<0) { |
635 | iInteger=i; |
636 | } else { |
637 | // check |
638 | const int * a = static_cast<CoinModel *>(blocks_[iInteger])->integerTypeArray(); |
639 | const int * b = static_cast<CoinModel *>(blocks_[i])->integerTypeArray(); |
640 | if (!sameValues(a,b,numberColumns)) |
641 | numberErrors++; |
642 | } |
643 | } |
644 | } |
645 | } |
646 | return numberErrors; |
647 | } |
648 | /* Fill in info structure and update counts |
649 | */ |
650 | void |
651 | CoinStructuredModel::fillInfo(CoinModelBlockInfo & info, |
652 | const CoinStructuredModel * block) |
653 | { |
654 | int numberRows = block->numberRows(); |
655 | int numberColumns = block->numberColumns(); |
656 | // Which block |
657 | int iRowBlock=addRowBlock(numberRows,block->rowBlockName_); |
658 | info.rowBlock=iRowBlock; |
659 | int iColumnBlock=addColumnBlock(numberColumns,block->columnBlockName_); |
660 | info.columnBlock=iColumnBlock; |
661 | } |
662 | /* add a block from a CoinModel without names*/ |
663 | int |
664 | CoinStructuredModel::addBlock(const std::string & rowBlock, |
665 | const std::string & columnBlock, |
666 | const CoinBaseModel & block) |
667 | { |
668 | CoinBaseModel * block2 = block.clone(); |
669 | return addBlock(rowBlock,columnBlock,block2); |
670 | } |
671 | /* add a block using names |
672 | */ |
673 | int |
674 | CoinStructuredModel::addBlock(const std::string & rowBlock, |
675 | const std::string & columnBlock, |
676 | const CoinPackedMatrix & matrix, |
677 | const double * rowLower, const double * rowUpper, |
678 | const double * columnLower, const double * columnUpper, |
679 | const double * objective) |
680 | { |
681 | CoinModel * block = new CoinModel(); |
682 | block->loadBlock(matrix,columnLower,columnUpper,objective, |
683 | rowLower,rowUpper); |
684 | return addBlock(rowBlock,columnBlock,block); |
685 | } |
686 | /* add a block from a CoinModel without names*/ |
687 | int |
688 | CoinStructuredModel::addBlock(const std::string & rowBlock, |
689 | const std::string & columnBlock, |
690 | CoinBaseModel * block) |
691 | { |
692 | if (numberElementBlocks_==maximumElementBlocks_) { |
693 | maximumElementBlocks_ = 3*(maximumElementBlocks_+10)/2; |
694 | CoinBaseModel ** temp = new CoinBaseModel * [maximumElementBlocks_]; |
695 | memcpy(temp,blocks_,numberElementBlocks_*sizeof(CoinBaseModel *)); |
696 | delete [] blocks_; |
697 | blocks_ = temp; |
698 | CoinModelBlockInfo * temp2 = new CoinModelBlockInfo [maximumElementBlocks_]; |
699 | memcpy(temp2,blockType_,numberElementBlocks_*sizeof(CoinModelBlockInfo)); |
700 | delete [] blockType_; |
701 | blockType_ = temp2; |
702 | if (coinModelBlocks_) { |
703 | CoinModel ** temp = new CoinModel * [maximumElementBlocks_]; |
704 | CoinZeroN(temp,maximumElementBlocks_); |
705 | memcpy(temp,coinModelBlocks_,numberElementBlocks_*sizeof(CoinModel *)); |
706 | delete [] coinModelBlocks_; |
707 | coinModelBlocks_ = temp; |
708 | } |
709 | } |
710 | blocks_[numberElementBlocks_++]=block; |
711 | block->setRowBlock(rowBlock); |
712 | block->setColumnBlock(columnBlock); |
713 | int numberErrors=0; |
714 | CoinModel * coinBlock = dynamic_cast<CoinModel *>(block); |
715 | if (coinBlock) { |
716 | // Convert matrix |
717 | if (coinBlock->type()!=3) |
718 | coinBlock->convertMatrix(); |
719 | numberErrors=fillInfo(blockType_[numberElementBlocks_-1],coinBlock); |
720 | } else { |
721 | CoinStructuredModel * subModel = dynamic_cast<CoinStructuredModel *>(block); |
722 | assert (subModel); |
723 | CoinModel * blockX = |
724 | subModel->coinModelBlock(blockType_[numberElementBlocks_-1]); |
725 | fillInfo(blockType_[numberElementBlocks_-1],subModel); |
726 | setCoinModel(blockX,numberElementBlocks_-1); |
727 | } |
728 | return numberErrors; |
729 | } |
730 | |
731 | /* add a block from a CoinModel with names*/ |
732 | int |
733 | CoinStructuredModel::addBlock(const CoinBaseModel & block) |
734 | { |
735 | |
736 | //inline const std::string & getRowBlock() const |
737 | //abort(); |
738 | return addBlock(block.getRowBlock(),block.getColumnBlock(), |
739 | block); |
740 | } |
741 | /* Decompose a model specified as arrays + CoinPackedMatrix |
742 | 1 - try D-W |
743 | 2 - try Benders |
744 | 3 - try Staircase |
745 | Returns number of blocks or zero if no structure |
746 | */ |
747 | int |
748 | CoinStructuredModel::decompose(const CoinPackedMatrix & matrix, |
749 | const double * rowLower, const double * rowUpper, |
750 | const double * columnLower, const double * columnUpper, |
751 | const double * objective, int type,int maxBlocks, |
752 | double objectiveOffset) |
753 | { |
754 | setObjectiveOffset(objectiveOffset); |
755 | int numberBlocks=0; |
756 | if (type==1) { |
757 | // Try master at top and bottom |
758 | bool goodDW=true; |
759 | // get row copy |
760 | CoinPackedMatrix rowCopy = matrix; |
761 | rowCopy.reverseOrdering(); |
762 | const int * row = matrix.getIndices(); |
763 | const int * columnLength = matrix.getVectorLengths(); |
764 | const CoinBigIndex * columnStart = matrix.getVectorStarts(); |
765 | //const double * elementByColumn = matrix.getElements(); |
766 | const int * column = rowCopy.getIndices(); |
767 | const int * rowLength = rowCopy.getVectorLengths(); |
768 | const CoinBigIndex * rowStart = rowCopy.getVectorStarts(); |
769 | //const double * elementByRow = rowCopy.getElements(); |
770 | int numberRows = matrix.getNumRows(); |
771 | int * rowBlock = new int[numberRows+1]; |
772 | int iRow; |
773 | // Row counts (maybe look at long rows for master) |
774 | CoinZeroN(rowBlock,numberRows+1); |
775 | for (iRow=0;iRow<numberRows;iRow++) { |
776 | int length = rowLength[iRow]; |
777 | rowBlock[length]++; |
778 | } |
779 | for (iRow=0;iRow<numberRows+1;iRow++) { |
780 | if (rowBlock[iRow]) |
781 | printf("%d rows have %d elements\n" ,rowBlock[iRow],iRow); |
782 | } |
783 | bool newWay=true; |
784 | // to say if column looked at |
785 | int numberColumns = matrix.getNumCols(); |
786 | int * columnBlock = new int[numberColumns]; |
787 | int iColumn; |
788 | int * whichRow = new int [numberRows]; |
789 | int * whichColumn = new int [numberColumns]; |
790 | int * stack = new int [numberRows]; |
791 | if (newWay) { |
792 | //double best2[3]={COIN_DBL_MAX,COIN_DBL_MAX,COIN_DBL_MAX}; |
793 | double best2[3]={0.0,0.0,0.0}; |
794 | int row2[3]={-1,-1,-1}; |
795 | // try forward and backward and sorted |
796 | for (int iWay=0;iWay<3;iWay++) { |
797 | if (iWay==0) { |
798 | // forwards |
799 | for (int i=0;i<numberRows;i++) |
800 | stack[i]=i; |
801 | } else if (iWay==1) { |
802 | // backwards |
803 | for (int i=0;i<numberRows;i++) |
804 | stack[i]=numberRows-1-i; |
805 | } else { |
806 | // sparsest first |
807 | for (int i=0;i<numberRows;i++) { |
808 | rowBlock[i]=rowLength[i]; |
809 | stack[i]=i; |
810 | } |
811 | CoinSort_2(rowBlock,rowBlock+numberRows,stack); |
812 | } |
813 | CoinFillN(rowBlock,numberRows,-1); |
814 | rowBlock[numberRows]=0; |
815 | CoinFillN(whichRow,numberRows,-1); |
816 | CoinFillN(columnBlock,numberColumns,-1); |
817 | CoinFillN(whichColumn,numberColumns,-1); |
818 | int numberMarkedColumns = 0; |
819 | numberBlocks=0; |
820 | int bestRow = -1; |
821 | int maximumInBlock = 0; |
822 | int rowsDone=0; |
823 | int checkAfterRows = (5*numberRows)/10+1; |
824 | #define OSL_WAY |
825 | #ifdef OSL_WAY |
826 | double best = COIN_DBL_MAX; |
827 | int bestRowsDone=-1; |
828 | #else |
829 | double best = 0.0; //COIN_DBL_MAX; |
830 | #endif |
831 | int numberGoodBlocks=0; |
832 | |
833 | for (int kRow=0;kRow<numberRows;kRow++) { |
834 | iRow = stack[kRow]; |
835 | CoinBigIndex start = rowStart[iRow]; |
836 | CoinBigIndex end = start+rowLength[iRow]; |
837 | int iBlock=-1; |
838 | for (CoinBigIndex j=start;j<end;j++) { |
839 | int iColumn = column[j]; |
840 | if (columnBlock[iColumn]>=0) { |
841 | // already marked |
842 | if (iBlock<0) { |
843 | iBlock = columnBlock[iColumn]; |
844 | } else if (iBlock != columnBlock[iColumn]) { |
845 | // join two blocks |
846 | int jBlock = columnBlock[iColumn]; |
847 | numberGoodBlocks--; |
848 | // Increase count of iBlock |
849 | rowBlock[iBlock] += rowBlock[jBlock]; |
850 | rowBlock[jBlock]=0; |
851 | // First column of block jBlock |
852 | int jColumn = whichRow[jBlock]; |
853 | while (jColumn>=0) { |
854 | columnBlock[jColumn]=iBlock; |
855 | iColumn = jColumn; |
856 | jColumn = whichColumn[jColumn]; |
857 | } |
858 | whichColumn[iColumn] = whichRow[iBlock]; |
859 | whichRow[iBlock] = whichRow[jBlock]; |
860 | whichRow[jBlock]=-1; |
861 | } |
862 | } |
863 | } |
864 | int n=end-start; |
865 | // If not in block - then start one |
866 | if (iBlock<0) { |
867 | // unless null row |
868 | if (n) { |
869 | iBlock = numberBlocks; |
870 | numberBlocks++; |
871 | numberGoodBlocks++; |
872 | int jColumn = column[start]; |
873 | columnBlock[jColumn]=iBlock; |
874 | whichRow[iBlock]=jColumn; |
875 | numberMarkedColumns += n; |
876 | rowBlock[iBlock] = n; |
877 | for (CoinBigIndex j=start+1;j<end;j++) { |
878 | int iColumn = column[j]; |
879 | columnBlock[iColumn]=iBlock; |
880 | whichColumn[jColumn]=iColumn; |
881 | jColumn=iColumn; |
882 | } |
883 | } else { |
884 | rowBlock[numberRows]++; |
885 | } |
886 | } else { |
887 | // add all to this block if not already in |
888 | int jColumn = whichRow[iBlock]; |
889 | for (CoinBigIndex j=start;j<end;j++) { |
890 | int iColumn = column[j]; |
891 | if (columnBlock[iColumn]<0) { |
892 | numberMarkedColumns++; |
893 | rowBlock[iBlock]++; |
894 | columnBlock[iColumn]=iBlock; |
895 | whichColumn[iColumn]=jColumn; |
896 | jColumn=iColumn; |
897 | } |
898 | } |
899 | whichRow[iBlock]=jColumn; |
900 | } |
901 | #if 0 |
902 | { |
903 | int nn=0; |
904 | int * temp = new int [numberRows]; |
905 | CoinZeroN(temp,numberRows); |
906 | for (int i=0;i<numberColumns;i++) { |
907 | int iBlock = columnBlock[i]; |
908 | if (iBlock>=0) { |
909 | nn++; |
910 | assert (iBlock<numberBlocks); |
911 | temp[iBlock]++; |
912 | } |
913 | } |
914 | for (int i=0;i<numberBlocks;i++) |
915 | assert (temp[i]==rowBlock[i]); |
916 | assert (nn==numberMarkedColumns); |
917 | for (int i=0;i<numberBlocks;i++) { |
918 | // First column of block i |
919 | int jColumn = whichRow[i]; |
920 | int n=0; |
921 | while (jColumn>=0) { |
922 | n++; |
923 | jColumn = whichColumn[jColumn]; |
924 | } |
925 | assert (n==temp[i]); |
926 | } |
927 | delete [] temp; |
928 | } |
929 | #endif |
930 | rowsDone++; |
931 | if (iBlock>=0) |
932 | maximumInBlock = CoinMax(maximumInBlock,rowBlock[iBlock]); |
933 | if (rowsDone>=checkAfterRows) { |
934 | assert (numberGoodBlocks>0); |
935 | double averageSize = static_cast<double>(numberMarkedColumns)/ |
936 | static_cast<double>(numberGoodBlocks); |
937 | #ifndef OSL_WAY |
938 | double notionalBlocks = static_cast<double>(numberMarkedColumns)/ |
939 | averageSize; |
940 | if (maximumInBlock<3*averageSize&&numberGoodBlocks>2) { |
941 | if(best*(numberRows-rowsDone) < notionalBlocks) { |
942 | best = notionalBlocks/ |
943 | static_cast<double> (numberRows-rowsDone); |
944 | bestRow = kRow; |
945 | } |
946 | } |
947 | #else |
948 | if (maximumInBlock*10<numberColumns*11&&numberGoodBlocks>1) { |
949 | double test = maximumInBlock + 0.0*averageSize; |
950 | if(best*static_cast<double>(rowsDone) > test) { |
951 | best = test/static_cast<double> (rowsDone); |
952 | bestRow = kRow; |
953 | bestRowsDone=rowsDone; |
954 | } |
955 | } |
956 | #endif |
957 | } |
958 | } |
959 | #ifndef OSL_WAY |
960 | best2[iWay]=best; |
961 | #else |
962 | if (bestRowsDone<numberRows) |
963 | best2[iWay]=-(numberRows-bestRowsDone); |
964 | else |
965 | best2[iWay]=-numberRows; |
966 | #endif |
967 | row2[iWay]=bestRow; |
968 | } |
969 | // mark rows |
970 | int nMaster; |
971 | CoinFillN(rowBlock,numberRows,-2); |
972 | if (best2[2]<best2[0]||best2[2]<best2[1]) { |
973 | int iRow1; |
974 | int iRow2; |
975 | if (best2[0]>best2[1]) { |
976 | // Bottom rows in master |
977 | iRow1=row2[0]+1; |
978 | iRow2=numberRows; |
979 | } else { |
980 | // Top rows in master |
981 | iRow1=0; |
982 | iRow2=numberRows-row2[1]; |
983 | } |
984 | nMaster = iRow2-iRow1; |
985 | CoinFillN(rowBlock+iRow1,nMaster,-1); |
986 | } else { |
987 | // sorted |
988 | // Bottom rows in master (in order) |
989 | int iRow1=row2[2]+1; |
990 | nMaster = numberRows-iRow1; |
991 | for (int i=iRow1;i<numberRows;i++) |
992 | rowBlock[stack[i]]=-1; |
993 | } |
994 | if (nMaster*2>numberRows) { |
995 | goodDW=false; |
996 | printf("%d rows out of %d would be in master - no good\n" , |
997 | nMaster,numberRows); |
998 | delete [] rowBlock; |
999 | delete [] columnBlock; |
1000 | delete [] whichRow; |
1001 | delete [] whichColumn; |
1002 | delete [] stack; |
1003 | CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper, |
1004 | columnLower,columnUpper,objective); |
1005 | model.setObjectiveOffset(objectiveOffset); |
1006 | addBlock("row_master" ,"column_master" ,model); |
1007 | return 0; |
1008 | } |
1009 | } else { |
1010 | for (iRow=0;iRow<numberRows;iRow++) |
1011 | rowBlock[iRow]=-2; |
1012 | // these are master rows |
1013 | if (numberRows==105127) { |
1014 | // ken-18 |
1015 | for (iRow=104976;iRow<numberRows;iRow++) |
1016 | rowBlock[iRow]=-1; |
1017 | } else if (numberRows==2426) { |
1018 | // ken-7 |
1019 | for (iRow=2401;iRow<numberRows;iRow++) |
1020 | rowBlock[iRow]=-1; |
1021 | } else if (numberRows==810) { |
1022 | for (iRow=81;iRow<84;iRow++) |
1023 | rowBlock[iRow]=-1; |
1024 | } else if (numberRows==5418) { |
1025 | for (iRow=564;iRow<603;iRow++) |
1026 | rowBlock[iRow]=-1; |
1027 | } else if (numberRows==10280) { |
1028 | // osa-60 |
1029 | for (iRow=10198;iRow<10280;iRow++) |
1030 | rowBlock[iRow]=-1; |
1031 | } else if (numberRows==1503) { |
1032 | // degen3 |
1033 | for (iRow=0;iRow<561;iRow++) |
1034 | rowBlock[iRow]=-1; |
1035 | } else if (numberRows==929) { |
1036 | // czprob |
1037 | for (iRow=0;iRow<39;iRow++) |
1038 | rowBlock[iRow]=-1; |
1039 | } else if (numberRows==33874) { |
1040 | // pds-20 |
1041 | for (iRow=31427;iRow<33874;iRow++) |
1042 | rowBlock[iRow]=-1; |
1043 | } else if (numberRows==24902) { |
1044 | // allgrade |
1045 | int kRow=818; |
1046 | for (iRow=0;iRow<kRow;iRow++) |
1047 | rowBlock[iRow]=-1; |
1048 | } |
1049 | } |
1050 | numberBlocks=0; |
1051 | CoinFillN(columnBlock,numberColumns,-2); |
1052 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1053 | int kstart = columnStart[iColumn]; |
1054 | int kend = columnStart[iColumn]+columnLength[iColumn]; |
1055 | if (columnBlock[iColumn]==-2) { |
1056 | // column not allocated |
1057 | int j; |
1058 | int nstack=0; |
1059 | for (j=kstart;j<kend;j++) { |
1060 | int iRow= row[j]; |
1061 | if (rowBlock[iRow]!=-1) { |
1062 | assert(rowBlock[iRow]==-2); |
1063 | rowBlock[iRow]=numberBlocks; // mark |
1064 | stack[nstack++] = iRow; |
1065 | } |
1066 | } |
1067 | if (nstack) { |
1068 | // new block - put all connected in |
1069 | numberBlocks++; |
1070 | columnBlock[iColumn]=numberBlocks-1; |
1071 | while (nstack) { |
1072 | int iRow = stack[--nstack]; |
1073 | int k; |
1074 | for (k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) { |
1075 | int iColumn = column[k]; |
1076 | int kkstart = columnStart[iColumn]; |
1077 | int kkend = kkstart + columnLength[iColumn]; |
1078 | if (columnBlock[iColumn]==-2) { |
1079 | columnBlock[iColumn]=numberBlocks-1; // mark |
1080 | // column not allocated |
1081 | int jj; |
1082 | for (jj=kkstart;jj<kkend;jj++) { |
1083 | int jRow= row[jj]; |
1084 | if (rowBlock[jRow]==-2) { |
1085 | rowBlock[jRow]=numberBlocks-1; |
1086 | stack[nstack++]=jRow; |
1087 | } |
1088 | } |
1089 | } else { |
1090 | assert (columnBlock[iColumn]==numberBlocks-1); |
1091 | } |
1092 | } |
1093 | } |
1094 | } else { |
1095 | // Only in master |
1096 | columnBlock[iColumn]=-1; |
1097 | } |
1098 | } |
1099 | } |
1100 | delete [] stack; |
1101 | int numberMasterRows=0; |
1102 | for (iRow=0;iRow<numberRows;iRow++) { |
1103 | int iBlock = rowBlock[iRow]; |
1104 | if (iBlock==-1) |
1105 | numberMasterRows++; |
1106 | } |
1107 | int numberMasterColumns=0; |
1108 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1109 | int iBlock = columnBlock[iColumn]; |
1110 | if (iBlock==-1) |
1111 | numberMasterColumns++; |
1112 | } |
1113 | if (numberBlocks<=maxBlocks) |
1114 | printf("%d blocks found - %d rows, %d columns in master\n" , |
1115 | numberBlocks,numberMasterRows,numberMasterColumns); |
1116 | else |
1117 | printf("%d blocks found (reduced to %d) - %d rows, %d columns in master\n" , |
1118 | numberBlocks,maxBlocks,numberMasterRows,numberMasterColumns); |
1119 | if (numberBlocks) { |
1120 | if (numberBlocks>maxBlocks) { |
1121 | int iBlock; |
1122 | for (iRow=0;iRow<numberRows;iRow++) { |
1123 | iBlock = rowBlock[iRow]; |
1124 | if (iBlock>=0) |
1125 | rowBlock[iRow] = iBlock%maxBlocks; |
1126 | } |
1127 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1128 | iBlock = columnBlock[iColumn]; |
1129 | if (iBlock>=0) |
1130 | columnBlock[iColumn] = iBlock%maxBlocks; |
1131 | } |
1132 | numberBlocks=maxBlocks; |
1133 | } |
1134 | } |
1135 | // make up problems |
1136 | // Create all sub problems |
1137 | // Space for creating |
1138 | double * obj = new double [numberColumns]; |
1139 | double * columnLo = new double [numberColumns]; |
1140 | double * columnUp = new double [numberColumns]; |
1141 | double * rowLo = new double [numberRows]; |
1142 | double * rowUp = new double [numberRows]; |
1143 | // Counts |
1144 | int * rowCount = reinterpret_cast<int *>(rowLo); |
1145 | CoinZeroN(rowCount,numberBlocks); |
1146 | for (int i=0;i<numberRows;i++) { |
1147 | int iBlock=rowBlock[i]; |
1148 | if (iBlock>=0) |
1149 | rowCount[iBlock]++; |
1150 | } |
1151 | // allocate empty rows |
1152 | for (int i=0;i<numberRows;i++) { |
1153 | int iBlock=rowBlock[i]; |
1154 | if (iBlock==-2) { |
1155 | // find block with smallest count |
1156 | int iSmall=-1; |
1157 | int smallest = numberRows; |
1158 | for (int j=0;j<numberBlocks;j++) { |
1159 | if (rowCount[j]<smallest) { |
1160 | iSmall=j; |
1161 | smallest=rowCount[j]; |
1162 | } |
1163 | } |
1164 | rowBlock[i]=iSmall; |
1165 | rowCount[iSmall]++; |
1166 | } |
1167 | } |
1168 | int * columnCount = reinterpret_cast<int *>(rowUp); |
1169 | CoinZeroN(columnCount,numberBlocks); |
1170 | for (int i=0;i<numberColumns;i++) { |
1171 | int iBlock=columnBlock[i]; |
1172 | if (iBlock>=0) |
1173 | columnCount[iBlock]++; |
1174 | } |
1175 | int maximumSize=0; |
1176 | for (int i=0;i<numberBlocks;i++) { |
1177 | printf("Block %d has %d rows and %d columns\n" ,i, |
1178 | rowCount[i],columnCount[i]); |
1179 | int k=2*rowCount[i]+columnCount[i]; |
1180 | maximumSize = CoinMax(maximumSize,k); |
1181 | } |
1182 | if (maximumSize*10>4*(2*numberRows+numberColumns)) { |
1183 | // No good |
1184 | printf("Doesn't look good\n" ); |
1185 | delete [] rowBlock; |
1186 | delete [] columnBlock; |
1187 | delete [] whichRow; |
1188 | delete [] whichColumn; |
1189 | delete [] obj ; |
1190 | delete [] columnLo ; |
1191 | delete [] columnUp ; |
1192 | delete [] rowLo ; |
1193 | delete [] rowUp ; |
1194 | CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper, |
1195 | columnLower,columnUpper,objective); |
1196 | model.setObjectiveOffset(objectiveOffset); |
1197 | addBlock("row_master" ,"column_master" ,model); |
1198 | return 0; |
1199 | } |
1200 | // Name for master so at top |
1201 | addRowBlock(numberMasterRows,"row_master" ); |
1202 | // Arrays |
1203 | // get full matrix |
1204 | CoinPackedMatrix fullMatrix = matrix; |
1205 | int numberRow2,numberColumn2; |
1206 | int iBlock; |
1207 | for (iBlock=0;iBlock<numberBlocks;iBlock++) { |
1208 | char rowName[20]; |
1209 | sprintf(rowName,"row_%d" ,iBlock); |
1210 | char columnName[20]; |
1211 | sprintf(columnName,"column_%d" ,iBlock); |
1212 | numberRow2=0; |
1213 | numberColumn2=0; |
1214 | for (iRow=0;iRow<numberRows;iRow++) { |
1215 | if (iBlock==rowBlock[iRow]) { |
1216 | rowLo[numberRow2]=rowLower[iRow]; |
1217 | rowUp[numberRow2]=rowUpper[iRow]; |
1218 | whichRow[numberRow2++]=iRow; |
1219 | } |
1220 | } |
1221 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1222 | if (iBlock==columnBlock[iColumn]) { |
1223 | obj[numberColumn2]=objective[iColumn]; |
1224 | columnLo[numberColumn2]=columnLower[iColumn]; |
1225 | columnUp[numberColumn2]=columnUpper[iColumn]; |
1226 | whichColumn[numberColumn2++]=iColumn; |
1227 | } |
1228 | } |
1229 | // Diagonal block |
1230 | CoinPackedMatrix mat(fullMatrix, |
1231 | numberRow2,whichRow, |
1232 | numberColumn2,whichColumn); |
1233 | // make sure correct dimensions |
1234 | mat.setDimensions(numberRow2,numberColumn2); |
1235 | CoinModel * block = new CoinModel(numberRow2,numberColumn2,&mat, |
1236 | rowLo,rowUp,NULL,NULL,NULL); |
1237 | block->setOriginalIndices(whichRow,whichColumn); |
1238 | addBlock(rowName,columnName,block); // takes ownership |
1239 | // and top block |
1240 | numberRow2=0; |
1241 | // get top matrix |
1242 | for (iRow=0;iRow<numberRows;iRow++) { |
1243 | int iBlock = rowBlock[iRow]; |
1244 | if (iBlock==-1) { |
1245 | whichRow[numberRow2++]=iRow; |
1246 | } |
1247 | } |
1248 | CoinPackedMatrix top(fullMatrix, |
1249 | numberRow2,whichRow, |
1250 | numberColumn2,whichColumn); |
1251 | // make sure correct dimensions |
1252 | top.setDimensions(numberRow2,numberColumn2); |
1253 | block = new CoinModel(numberMasterRows,numberColumn2,&top, |
1254 | NULL,NULL,columnLo,columnUp,obj); |
1255 | block->setOriginalIndices(whichRow,whichColumn); |
1256 | addBlock("row_master" ,columnName,block); // takes ownership |
1257 | } |
1258 | // and master |
1259 | numberRow2=0; |
1260 | numberColumn2=0; |
1261 | for (iRow=0;iRow<numberRows;iRow++) { |
1262 | int iBlock = rowBlock[iRow]; |
1263 | if (iBlock==-1) { |
1264 | rowLo[numberRow2]=rowLower[iRow]; |
1265 | rowUp[numberRow2]=rowUpper[iRow]; |
1266 | whichRow[numberRow2++]=iRow; |
1267 | } |
1268 | } |
1269 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1270 | int iBlock = columnBlock[iColumn]; |
1271 | if (iBlock<0) { |
1272 | obj[numberColumn2]=objective[iColumn]; |
1273 | columnLo[numberColumn2]=columnLower[iColumn]; |
1274 | columnUp[numberColumn2]=columnUpper[iColumn]; |
1275 | whichColumn[numberColumn2++]=iColumn; |
1276 | } |
1277 | } |
1278 | delete [] rowBlock; |
1279 | delete [] columnBlock; |
1280 | CoinPackedMatrix top(fullMatrix, |
1281 | numberRow2,whichRow, |
1282 | numberColumn2,whichColumn); |
1283 | // make sure correct dimensions |
1284 | top.setDimensions(numberRow2,numberColumn2); |
1285 | CoinModel * block = new CoinModel(numberRow2,numberColumn2,&top, |
1286 | rowLo,rowUp, |
1287 | columnLo,columnUp,obj); |
1288 | block->setOriginalIndices(whichRow,whichColumn); |
1289 | addBlock("row_master" ,"column_master" ,block); // takes ownership |
1290 | delete [] whichRow; |
1291 | delete [] whichColumn; |
1292 | delete [] obj ; |
1293 | delete [] columnLo ; |
1294 | delete [] columnUp ; |
1295 | delete [] rowLo ; |
1296 | delete [] rowUp ; |
1297 | } else if (type==2) { |
1298 | // Try master at beginning and end |
1299 | bool goodBenders=true; |
1300 | // get row copy |
1301 | CoinPackedMatrix rowCopy = matrix; |
1302 | rowCopy.reverseOrdering(); |
1303 | const int * row = matrix.getIndices(); |
1304 | const int * columnLength = matrix.getVectorLengths(); |
1305 | const CoinBigIndex * columnStart = matrix.getVectorStarts(); |
1306 | //const double * elementByColumn = matrix.getElements(); |
1307 | const int * column = rowCopy.getIndices(); |
1308 | const int * rowLength = rowCopy.getVectorLengths(); |
1309 | const CoinBigIndex * rowStart = rowCopy.getVectorStarts(); |
1310 | //const double * elementByRow = rowCopy.getElements(); |
1311 | int numberColumns = matrix.getNumCols(); |
1312 | int * columnBlock = new int[numberColumns+1]; |
1313 | int iColumn; |
1314 | // Column counts (maybe look at long columns for master) |
1315 | CoinZeroN(columnBlock,numberColumns+1); |
1316 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1317 | int length = columnLength[iColumn]; |
1318 | columnBlock[length]++; |
1319 | } |
1320 | for (iColumn=0;iColumn<numberColumns+1;iColumn++) { |
1321 | if (columnBlock[iColumn]) |
1322 | printf("%d columns have %d elements\n" ,columnBlock[iColumn],iColumn); |
1323 | } |
1324 | bool newWay=false; |
1325 | // to say if row looked at |
1326 | int numberRows = matrix.getNumRows(); |
1327 | int * rowBlock = new int[numberRows]; |
1328 | int iRow; |
1329 | int * whichRow = new int [numberRows]; |
1330 | int * whichColumn = new int [numberColumns]; |
1331 | int * stack = new int [numberColumns]; |
1332 | if (newWay) { |
1333 | //double best2[3]={COIN_DBL_MAX,COIN_DBL_MAX,COIN_DBL_MAX}; |
1334 | double best2[3]={0.0,0.0,0.0}; |
1335 | int column2[3]={-1,-1,-1}; |
1336 | // try forward and backward and sorted |
1337 | for (int iWay=0;iWay<3;iWay++) { |
1338 | if (iWay==0) { |
1339 | // forwards |
1340 | for (int i=0;i<numberColumns;i++) |
1341 | stack[i]=i; |
1342 | } else if (iWay==1) { |
1343 | // backwards |
1344 | for (int i=0;i<numberColumns;i++) |
1345 | stack[i]=numberColumns-1-i; |
1346 | } else { |
1347 | // sparsest first |
1348 | for (int i=0;i<numberColumns;i++) { |
1349 | columnBlock[i]=columnLength[i]; |
1350 | stack[i]=i; |
1351 | } |
1352 | CoinSort_2(columnBlock,columnBlock+numberColumns,stack); |
1353 | } |
1354 | CoinFillN(columnBlock,numberColumns,-1); |
1355 | columnBlock[numberColumns]=0; |
1356 | CoinFillN(whichColumn,numberColumns,-1); |
1357 | CoinFillN(rowBlock,numberRows,-1); |
1358 | CoinFillN(whichRow,numberRows,-1); |
1359 | int numberMarkedRows = 0; |
1360 | numberBlocks=0; |
1361 | int bestColumn = -1; |
1362 | int maximumInBlock = 0; |
1363 | int columnsDone=0; |
1364 | int checkAfterColumns = (5*numberColumns)/10+1; |
1365 | #ifdef OSL_WAY |
1366 | double best = COIN_DBL_MAX; |
1367 | int bestColumnsDone=-1; |
1368 | #else |
1369 | double best = 0.0; //COIN_DBL_MAX; |
1370 | #endif |
1371 | int numberGoodBlocks=0; |
1372 | |
1373 | for (int kColumn=0;kColumn<numberColumns;kColumn++) { |
1374 | iColumn = stack[kColumn]; |
1375 | CoinBigIndex start = columnStart[iColumn]; |
1376 | CoinBigIndex end = start+columnLength[iColumn]; |
1377 | int iBlock=-1; |
1378 | for (CoinBigIndex j=start;j<end;j++) { |
1379 | int iRow = row[j]; |
1380 | if (rowBlock[iRow]>=0) { |
1381 | // already marked |
1382 | if (iBlock<0) { |
1383 | iBlock = rowBlock[iRow]; |
1384 | } else if (iBlock != rowBlock[iRow]) { |
1385 | // join two blocks |
1386 | int jBlock = rowBlock[iRow]; |
1387 | numberGoodBlocks--; |
1388 | // Increase count of iBlock |
1389 | columnBlock[iBlock] += columnBlock[jBlock]; |
1390 | columnBlock[jBlock]=0; |
1391 | // First row of block jBlock |
1392 | int jRow = whichColumn[jBlock]; |
1393 | while (jRow>=0) { |
1394 | rowBlock[jRow]=iBlock; |
1395 | iRow = jRow; |
1396 | jRow = whichRow[jRow]; |
1397 | } |
1398 | whichRow[iRow] = whichColumn[iBlock]; |
1399 | whichColumn[iBlock] = whichColumn[jBlock]; |
1400 | whichColumn[jBlock]=-1; |
1401 | } |
1402 | } |
1403 | } |
1404 | int n=end-start; |
1405 | // If not in block - then start one |
1406 | if (iBlock<0) { |
1407 | // unless null column |
1408 | if (n) { |
1409 | iBlock = numberBlocks; |
1410 | numberBlocks++; |
1411 | numberGoodBlocks++; |
1412 | int jRow = row[start]; |
1413 | rowBlock[jRow]=iBlock; |
1414 | whichColumn[iBlock]=jRow; |
1415 | numberMarkedRows += n; |
1416 | columnBlock[iBlock] = n; |
1417 | for (CoinBigIndex j=start+1;j<end;j++) { |
1418 | int iRow = row[j]; |
1419 | rowBlock[iRow]=iBlock; |
1420 | whichRow[jRow]=iRow; |
1421 | jRow=iRow; |
1422 | } |
1423 | } else { |
1424 | columnBlock[numberColumns]++; |
1425 | } |
1426 | } else { |
1427 | // add all to this block if not already in |
1428 | int jRow = whichColumn[iBlock]; |
1429 | for (CoinBigIndex j=start;j<end;j++) { |
1430 | int iRow = row[j]; |
1431 | if (rowBlock[iRow]<0) { |
1432 | numberMarkedRows++; |
1433 | columnBlock[iBlock]++; |
1434 | rowBlock[iRow]=iBlock; |
1435 | whichRow[iRow]=jRow; |
1436 | jRow=iRow; |
1437 | } |
1438 | } |
1439 | whichColumn[iBlock]=jRow; |
1440 | } |
1441 | columnsDone++; |
1442 | if (iBlock>=0) |
1443 | maximumInBlock = CoinMax(maximumInBlock,columnBlock[iBlock]); |
1444 | if (columnsDone>=checkAfterColumns) { |
1445 | assert (numberGoodBlocks>0); |
1446 | double averageSize = static_cast<double>(numberMarkedRows)/ |
1447 | static_cast<double>(numberGoodBlocks); |
1448 | #ifndef OSL_WAY |
1449 | double notionalBlocks = static_cast<double>(numberMarkedRows)/ |
1450 | averageSize; |
1451 | if (maximumInBlock<3*averageSize&&numberGoodBlocks>2) { |
1452 | if(best*(numberColumns-columnsDone) < notionalBlocks) { |
1453 | best = notionalBlocks/ |
1454 | static_cast<double> (numberColumns-columnsDone); |
1455 | bestColumn = kColumn; |
1456 | } |
1457 | } |
1458 | #else |
1459 | if (maximumInBlock*10<numberRows*11&&numberGoodBlocks>1) { |
1460 | double test = maximumInBlock + 0.0*averageSize; |
1461 | if(best*static_cast<double>(columnsDone) > test) { |
1462 | best = test/static_cast<double> (columnsDone); |
1463 | bestColumn = kColumn; |
1464 | bestColumnsDone=columnsDone; |
1465 | } |
1466 | } |
1467 | #endif |
1468 | } |
1469 | } |
1470 | #ifndef OSL_WAY |
1471 | best2[iWay]=best; |
1472 | #else |
1473 | if (bestColumnsDone<numberColumns) |
1474 | best2[iWay]=-(numberColumns-bestColumnsDone); |
1475 | else |
1476 | best2[iWay]=-numberColumns; |
1477 | #endif |
1478 | column2[iWay]=bestColumn; |
1479 | } |
1480 | // mark columns |
1481 | int nMaster; |
1482 | CoinFillN(columnBlock,numberColumns,-2); |
1483 | if (best2[2]<best2[0]||best2[2]<best2[1]) { |
1484 | int iColumn1; |
1485 | int iColumn2; |
1486 | if (best2[0]>best2[1]) { |
1487 | // End columns in master |
1488 | iColumn1=column2[0]+1; |
1489 | iColumn2=numberColumns; |
1490 | } else { |
1491 | // Beginning columns in master |
1492 | iColumn1=0; |
1493 | iColumn2=numberColumns-column2[1]; |
1494 | } |
1495 | nMaster = iColumn2-iColumn1; |
1496 | CoinFillN(columnBlock+iColumn1,nMaster,-1); |
1497 | } else { |
1498 | // sorted |
1499 | // End columns in master (in order) |
1500 | int iColumn1=column2[2]+1; |
1501 | nMaster = numberColumns-iColumn1; |
1502 | for (int i=iColumn1;i<numberColumns;i++) |
1503 | columnBlock[stack[i]]=-1; |
1504 | } |
1505 | if (nMaster*2>numberColumns) { |
1506 | goodBenders=false; |
1507 | printf("%d columns out of %d would be in master - no good\n" , |
1508 | nMaster,numberColumns); |
1509 | delete [] rowBlock; |
1510 | delete [] columnBlock; |
1511 | delete [] whichRow; |
1512 | delete [] whichColumn; |
1513 | delete [] stack; |
1514 | CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper, |
1515 | columnLower,columnUpper,objective); |
1516 | model.setObjectiveOffset(objectiveOffset); |
1517 | addBlock("row_master" ,"column_master" ,model); |
1518 | return 0; |
1519 | } |
1520 | } else { |
1521 | for (iColumn=0;iColumn<numberColumns;iColumn++) |
1522 | columnBlock[iColumn]=-2; |
1523 | // these are master columns |
1524 | if (numberColumns==2426) { |
1525 | // ken-7 dual |
1526 | for (iColumn=2401;iColumn<numberColumns;iColumn++) |
1527 | columnBlock[iColumn]=-1; |
1528 | } |
1529 | } |
1530 | numberBlocks=0; |
1531 | CoinFillN(rowBlock,numberRows,-2); |
1532 | for (iRow=0;iRow<numberRows;iRow++) { |
1533 | int kstart = rowStart[iRow]; |
1534 | int kend = rowStart[iRow]+rowLength[iRow]; |
1535 | if (rowBlock[iRow]==-2) { |
1536 | // row not allocated |
1537 | int j; |
1538 | int nstack=0; |
1539 | for (j=kstart;j<kend;j++) { |
1540 | int iColumn= column[j]; |
1541 | if (columnBlock[iColumn]!=-1) { |
1542 | assert(columnBlock[iColumn]==-2); |
1543 | columnBlock[iColumn]=numberBlocks; // mark |
1544 | stack[nstack++] = iColumn; |
1545 | } |
1546 | } |
1547 | if (nstack) { |
1548 | // new block - put all connected in |
1549 | numberBlocks++; |
1550 | rowBlock[iRow]=numberBlocks-1; |
1551 | while (nstack) { |
1552 | int iColumn = stack[--nstack]; |
1553 | int k; |
1554 | for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];k++) { |
1555 | int iRow = row[k]; |
1556 | int kkstart = rowStart[iRow]; |
1557 | int kkend = kkstart + rowLength[iRow]; |
1558 | if (rowBlock[iRow]==-2) { |
1559 | rowBlock[iRow]=numberBlocks-1; // mark |
1560 | // row not allocated |
1561 | int jj; |
1562 | for (jj=kkstart;jj<kkend;jj++) { |
1563 | int jColumn= column[jj]; |
1564 | if (columnBlock[jColumn]==-2) { |
1565 | columnBlock[jColumn]=numberBlocks-1; |
1566 | stack[nstack++]=jColumn; |
1567 | } |
1568 | } |
1569 | } else { |
1570 | assert (rowBlock[iRow]==numberBlocks-1); |
1571 | } |
1572 | } |
1573 | } |
1574 | } else { |
1575 | // Only in master |
1576 | rowBlock[iRow]=-1; |
1577 | } |
1578 | } |
1579 | } |
1580 | delete [] stack; |
1581 | int numberMasterColumns=0; |
1582 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1583 | int iBlock = columnBlock[iColumn]; |
1584 | if (iBlock==-1) |
1585 | numberMasterColumns++; |
1586 | } |
1587 | int numberMasterRows=0; |
1588 | for (iRow=0;iRow<numberRows;iRow++) { |
1589 | int iBlock = rowBlock[iRow]; |
1590 | if (iBlock==-1) |
1591 | numberMasterRows++; |
1592 | } |
1593 | if (numberBlocks<=maxBlocks) |
1594 | printf("%d blocks found - %d columns, %d rows in master\n" , |
1595 | numberBlocks,numberMasterColumns,numberMasterRows); |
1596 | else |
1597 | printf("%d blocks found (reduced to %d) - %d columns, %d rows in master\n" , |
1598 | numberBlocks,maxBlocks,numberMasterColumns,numberMasterRows); |
1599 | if (numberBlocks) { |
1600 | if (numberBlocks>maxBlocks) { |
1601 | int iBlock; |
1602 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1603 | iBlock = columnBlock[iColumn]; |
1604 | if (iBlock>=0) |
1605 | columnBlock[iColumn] = iBlock%maxBlocks; |
1606 | } |
1607 | for (iRow=0;iRow<numberRows;iRow++) { |
1608 | iBlock = rowBlock[iRow]; |
1609 | if (iBlock>=0) |
1610 | rowBlock[iRow] = iBlock%maxBlocks; |
1611 | } |
1612 | numberBlocks=maxBlocks; |
1613 | } |
1614 | } |
1615 | // make up problems |
1616 | // Create all sub problems |
1617 | // Space for creating |
1618 | double * obj = new double [numberColumns]; |
1619 | double * rowLo = new double [numberRows]; |
1620 | double * rowUp = new double [numberRows]; |
1621 | double * columnLo = new double [numberColumns]; |
1622 | double * columnUp = new double [numberColumns]; |
1623 | // Counts |
1624 | int * columnCount = reinterpret_cast<int *>(columnLo); |
1625 | CoinZeroN(columnCount,numberBlocks); |
1626 | for (int i=0;i<numberColumns;i++) { |
1627 | int iBlock=columnBlock[i]; |
1628 | if (iBlock>=0) |
1629 | columnCount[iBlock]++; |
1630 | } |
1631 | // allocate empty columns |
1632 | for (int i=0;i<numberColumns;i++) { |
1633 | int iBlock=columnBlock[i]; |
1634 | if (iBlock==-2) { |
1635 | // find block with smallest count |
1636 | int iSmall=-1; |
1637 | int smallest = numberColumns; |
1638 | for (int j=0;j<numberBlocks;j++) { |
1639 | if (columnCount[j]<smallest) { |
1640 | iSmall=j; |
1641 | smallest=columnCount[j]; |
1642 | } |
1643 | } |
1644 | columnBlock[i]=iSmall; |
1645 | columnCount[iSmall]++; |
1646 | } |
1647 | } |
1648 | int * rowCount = reinterpret_cast<int *>(columnUp); |
1649 | CoinZeroN(rowCount,numberBlocks); |
1650 | for (int i=0;i<numberRows;i++) { |
1651 | int iBlock=rowBlock[i]; |
1652 | if (iBlock>=0) |
1653 | rowCount[iBlock]++; |
1654 | } |
1655 | int maximumSize=0; |
1656 | for (int i=0;i<numberBlocks;i++) { |
1657 | printf("Block %d has %d columns and %d rows\n" ,i, |
1658 | columnCount[i],rowCount[i]); |
1659 | int k=2*columnCount[i]+rowCount[i]; |
1660 | maximumSize = CoinMax(maximumSize,k); |
1661 | } |
1662 | if (maximumSize*10>4*(2*numberColumns+numberRows)) { |
1663 | // No good |
1664 | printf("Doesn't look good\n" ); |
1665 | delete [] rowBlock; |
1666 | delete [] columnBlock; |
1667 | delete [] whichRow; |
1668 | delete [] whichColumn; |
1669 | delete [] obj ; |
1670 | delete [] columnLo ; |
1671 | delete [] columnUp ; |
1672 | delete [] rowLo ; |
1673 | delete [] rowUp ; |
1674 | CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper, |
1675 | columnLower,columnUpper,objective); |
1676 | model.setObjectiveOffset(objectiveOffset); |
1677 | addBlock("row_master" ,"column_master" ,model); |
1678 | return 0; |
1679 | } |
1680 | // Name for master so at beginning |
1681 | addColumnBlock(numberMasterColumns,"column_master" ); |
1682 | // Arrays |
1683 | // get full matrix |
1684 | CoinPackedMatrix fullMatrix = matrix; |
1685 | int numberRow2,numberColumn2; |
1686 | int iBlock; |
1687 | for (iBlock=0;iBlock<numberBlocks;iBlock++) { |
1688 | char rowName[20]; |
1689 | sprintf(rowName,"row_%d" ,iBlock); |
1690 | char columnName[20]; |
1691 | sprintf(columnName,"column_%d" ,iBlock); |
1692 | numberRow2=0; |
1693 | numberColumn2=0; |
1694 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1695 | if (iBlock==columnBlock[iColumn]) { |
1696 | obj[numberColumn2]=objective[iColumn]; |
1697 | columnLo[numberColumn2]=columnLower[iColumn]; |
1698 | columnUp[numberColumn2]=columnUpper[iColumn]; |
1699 | whichColumn[numberColumn2++]=iColumn; |
1700 | } |
1701 | } |
1702 | for (iRow=0;iRow<numberRows;iRow++) { |
1703 | if (iBlock==rowBlock[iRow]) { |
1704 | rowLo[numberRow2]=rowLower[iRow]; |
1705 | rowUp[numberRow2]=rowUpper[iRow]; |
1706 | whichRow[numberRow2++]=iRow; |
1707 | } |
1708 | } |
1709 | // Diagonal block |
1710 | CoinPackedMatrix mat(fullMatrix, |
1711 | numberRow2,whichRow, |
1712 | numberColumn2,whichColumn); |
1713 | // make sure correct dimensions |
1714 | mat.setDimensions(numberRow2,numberColumn2); |
1715 | CoinModel * block = new CoinModel(numberRow2,numberColumn2,&mat, |
1716 | rowLo,rowUp,columnLo,columnUp,obj); |
1717 | block->setOriginalIndices(whichRow,whichColumn); |
1718 | addBlock(rowName,columnName,block); // takes ownership |
1719 | // and beginning block |
1720 | numberColumn2=0; |
1721 | // get beginning matrix |
1722 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1723 | int iBlock = columnBlock[iColumn]; |
1724 | if (iBlock==-1) { |
1725 | whichColumn[numberColumn2++]=iColumn; |
1726 | } |
1727 | } |
1728 | CoinPackedMatrix beginning(fullMatrix, |
1729 | numberRow2,whichRow, |
1730 | numberColumn2,whichColumn); |
1731 | // make sure correct dimensions ********* |
1732 | beginning.setDimensions(numberRow2,numberColumn2); |
1733 | block = new CoinModel(numberRow2,numberMasterColumns,&beginning, |
1734 | NULL,NULL,NULL,NULL,NULL); |
1735 | block->setOriginalIndices(whichRow,whichColumn); |
1736 | addBlock(rowName,"column_master" ,block); // takes ownership |
1737 | } |
1738 | // and master |
1739 | numberRow2=0; |
1740 | numberColumn2=0; |
1741 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
1742 | int iBlock = columnBlock[iColumn]; |
1743 | if (iBlock==-1) { |
1744 | obj[numberColumn2]=objective[iColumn]; |
1745 | columnLo[numberColumn2]=columnLower[iColumn]; |
1746 | columnUp[numberColumn2]=columnUpper[iColumn]; |
1747 | whichColumn[numberColumn2++]=iColumn; |
1748 | } |
1749 | } |
1750 | for (iRow=0;iRow<numberRows;iRow++) { |
1751 | int iBlock = rowBlock[iRow]; |
1752 | if (iBlock<0) { |
1753 | rowLo[numberRow2]=rowLower[iRow]; |
1754 | rowUp[numberRow2]=rowUpper[iRow]; |
1755 | whichRow[numberRow2++]=iRow; |
1756 | } |
1757 | } |
1758 | delete [] rowBlock; |
1759 | delete [] columnBlock; |
1760 | CoinPackedMatrix beginning(fullMatrix, |
1761 | numberRow2,whichRow, |
1762 | numberColumn2,whichColumn); |
1763 | // make sure correct dimensions |
1764 | beginning.setDimensions(numberRow2,numberColumn2); |
1765 | CoinModel * block = new CoinModel(numberRow2,numberColumn2,&beginning, |
1766 | rowLo,rowUp, |
1767 | columnLo,columnUp,obj); |
1768 | block->setOriginalIndices(whichRow,whichColumn); |
1769 | addBlock("row_master" ,"column_master" ,block); // takes ownership |
1770 | delete [] whichRow; |
1771 | delete [] whichColumn; |
1772 | delete [] obj ; |
1773 | delete [] columnLo ; |
1774 | delete [] columnUp ; |
1775 | delete [] rowLo ; |
1776 | delete [] rowUp ; |
1777 | } else { |
1778 | abort(); |
1779 | } |
1780 | return numberBlocks; |
1781 | } |
1782 | /* Decompose a CoinModel |
1783 | 1 - try D-W |
1784 | 2 - try Benders |
1785 | 3 - try Staircase |
1786 | Returns number of blocks or zero if no structure |
1787 | */ |
1788 | int |
1789 | CoinStructuredModel::decompose(const CoinModel & coinModel, int type, |
1790 | int maxBlocks) |
1791 | { |
1792 | const CoinPackedMatrix * matrix = coinModel.packedMatrix(); |
1793 | assert (matrix!=NULL); |
1794 | // Arrays |
1795 | const double * objective = coinModel.objectiveArray(); |
1796 | const double * columnLower = coinModel.columnLowerArray(); |
1797 | const double * columnUpper = coinModel.columnUpperArray(); |
1798 | const double * rowLower = coinModel.rowLowerArray(); |
1799 | const double * rowUpper = coinModel.rowUpperArray(); |
1800 | return decompose(*matrix, |
1801 | rowLower, rowUpper, |
1802 | columnLower, columnUpper, |
1803 | objective, type,maxBlocks, |
1804 | coinModel.objectiveOffset()); |
1805 | } |
1806 | // Return block corresponding to row and column |
1807 | const CoinBaseModel * |
1808 | CoinStructuredModel::block(int row,int column) const |
1809 | { |
1810 | const CoinBaseModel * block = NULL; |
1811 | if (blockType_) { |
1812 | for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) { |
1813 | if (blockType_[iBlock].rowBlock==row&& |
1814 | blockType_[iBlock].columnBlock==column) { |
1815 | block = blocks_[iBlock]; |
1816 | break; |
1817 | } |
1818 | } |
1819 | } |
1820 | return block; |
1821 | } |
1822 | // Return block corresponding to row and column as CoinModel |
1823 | const CoinBaseModel * |
1824 | CoinStructuredModel::coinBlock(int row,int column) const |
1825 | { |
1826 | const CoinModel * block = NULL; |
1827 | if (blockType_) { |
1828 | for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) { |
1829 | if (blockType_[iBlock].rowBlock==row&& |
1830 | blockType_[iBlock].columnBlock==column) { |
1831 | block = dynamic_cast<CoinModel *>(blocks_[iBlock]); |
1832 | assert (block); |
1833 | break; |
1834 | } |
1835 | } |
1836 | } |
1837 | return block; |
1838 | } |
1839 | /* Fill pointers corresponding to row and column. |
1840 | False if any missing */ |
1841 | CoinModelBlockInfo |
1842 | CoinStructuredModel::block(int row,int column, |
1843 | const double * & rowLower, const double * & rowUpper, |
1844 | const double * & columnLower, const double * & columnUpper, |
1845 | const double * & objective) const |
1846 | { |
1847 | CoinModelBlockInfo info; |
1848 | //memset(&info,0,sizeof(info)); |
1849 | rowLower=NULL; |
1850 | rowUpper=NULL; |
1851 | columnLower=NULL; |
1852 | columnUpper=NULL; |
1853 | objective=NULL; |
1854 | if (blockType_) { |
1855 | for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) { |
1856 | CoinModel * thisBlock = coinBlock(iBlock); |
1857 | if (blockType_[iBlock].rowBlock==row) { |
1858 | if (blockType_[iBlock].rhs) { |
1859 | info.rhs=1; |
1860 | rowLower = thisBlock->rowLowerArray(); |
1861 | rowUpper = thisBlock->rowUpperArray(); |
1862 | } |
1863 | } |
1864 | if (blockType_[iBlock].columnBlock==column) { |
1865 | if (blockType_[iBlock].bounds) { |
1866 | info.bounds=1; |
1867 | columnLower = thisBlock->columnLowerArray(); |
1868 | columnUpper = thisBlock->columnUpperArray(); |
1869 | objective = thisBlock->objectiveArray(); |
1870 | } |
1871 | } |
1872 | } |
1873 | } |
1874 | return info; |
1875 | } |
1876 | // Return block number corresponding to row and column |
1877 | int |
1878 | CoinStructuredModel::blockIndex(int row,int column) const |
1879 | { |
1880 | int block=-1; |
1881 | if (blockType_) { |
1882 | for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) { |
1883 | if (blockType_[iBlock].rowBlock==row&& |
1884 | blockType_[iBlock].columnBlock==column) { |
1885 | block = iBlock; |
1886 | break; |
1887 | } |
1888 | } |
1889 | } |
1890 | return block; |
1891 | } |
1892 | |