| 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 | |