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 //-------------------------------------------------------------------
CoinStructuredModel()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 */
CoinStructuredModel(const char * fileName,int decomposeType,int maxBlocks)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 //-------------------------------------------------------------------
CoinStructuredModel(const CoinStructuredModel & rhs)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 //-------------------------------------------------------------------
~CoinStructuredModel()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 *
clone() const113 CoinStructuredModel::clone() const
114 {
115 return new CoinStructuredModel(*this);
116 }
117
118 //----------------------------------------------------------------
119 // Assignment operator
120 //-------------------------------------------------------------------
121 CoinStructuredModel &
operator =(const CoinStructuredModel & rhs)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 }
sameValues(const double * a,const double * b,int n)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 }
sameValues(const int * a,const int * b,int n)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 }
sameValues(const CoinModel * a,const CoinModel * b,bool doRows)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
addRowBlock(int numberRows,const std::string & name)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
rowBlock(const std::string & name) const235 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
addColumnBlock(int numberColumns,const std::string & name)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
columnBlock(const std::string & name) const264 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
numberElements() const277 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 *
coinBlock(int i) const287 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 *
coinModelBlock(CoinModelBlockInfo & info)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
setCoinModel(CoinModel * block,int iBlock)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
refresh(int iBlock)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
fillInfo(CoinModelBlockInfo & info,const CoinModel * block)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
fillInfo(CoinModelBlockInfo & info,const CoinStructuredModel * block)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
addBlock(const std::string & rowBlock,const std::string & columnBlock,const CoinBaseModel & block)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
addBlock(const std::string & rowBlock,const std::string & columnBlock,const CoinPackedMatrix & matrix,const double * rowLower,const double * rowUpper,const double * columnLower,const double * columnUpper,const double * objective)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
addBlock(const std::string & rowBlock,const std::string & columnBlock,CoinBaseModel * block)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
addBlock(const CoinBaseModel & block)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
decompose(const CoinPackedMatrix & matrix,const double * rowLower,const double * rowUpper,const double * columnLower,const double * columnUpper,const double * objective,int type,int maxBlocks,double objectiveOffset)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
decompose(const CoinModel & coinModel,int type,int maxBlocks)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 *
block(int row,int column) const1808 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 *
coinBlock(int row,int column) const1824 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
block(int row,int column,const double * & rowLower,const double * & rowUpper,const double * & columnLower,const double * & columnUpper,const double * & objective) const1842 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
blockIndex(int row,int column) const1878 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