/* $Id: CoinStructuredModel.cpp 1373 2011-01-03 23:57:44Z lou $ */
// Copyright (C) 2008, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).


#include "CoinUtilsConfig.h"
#include "CoinHelperFunctions.hpp"
#include "CoinStructuredModel.hpp"
#include "CoinSort.hpp"
#include "CoinMpsIO.hpp"
#include "CoinFloatEqual.hpp"

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CoinStructuredModel::CoinStructuredModel ()
  :  CoinBaseModel(),
     numberRowBlocks_(0),
     numberColumnBlocks_(0),
     numberElementBlocks_(0),
     maximumElementBlocks_(0),
     blocks_(NULL),
     coinModelBlocks_(NULL),
     blockType_(NULL)
{
}
/* Read a problem in MPS or GAMS format from the given filename.
 */
CoinStructuredModel::CoinStructuredModel(const char *fileName,
					 int decomposeType,
					 int maxBlocks)
  :  CoinBaseModel(),
     numberRowBlocks_(0),
     numberColumnBlocks_(0),
     numberElementBlocks_(0),
     maximumElementBlocks_(0),
     blocks_(NULL),
     coinModelBlocks_(NULL),
     blockType_(NULL)
{
  CoinModel coinModel(fileName,false);
  if (coinModel.numberRows()) {
    problemName_ = coinModel.getProblemName();
    optimizationDirection_ = coinModel.optimizationDirection();
    objectiveOffset_ = coinModel.objectiveOffset();
    if (!decomposeType) {
      addBlock("row_master","column_master",coinModel);
    } else {
      const CoinPackedMatrix * matrix = coinModel.packedMatrix();
      if (!matrix)
	coinModel.convertMatrix();
      decompose(coinModel,decomposeType,maxBlocks);
      //addBlock("row_master","column_master",coinModel);
    }
  }
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CoinStructuredModel::CoinStructuredModel (const CoinStructuredModel & rhs)
  : CoinBaseModel(rhs),
    numberRowBlocks_(rhs.numberRowBlocks_),
    numberColumnBlocks_(rhs.numberColumnBlocks_),
    numberElementBlocks_(rhs.numberElementBlocks_),
    maximumElementBlocks_(rhs.maximumElementBlocks_)
{
  if (maximumElementBlocks_) {
    blocks_ = CoinCopyOfArray(rhs.blocks_,maximumElementBlocks_);
    for (int i=0;i<numberElementBlocks_;i++)
      blocks_[i]= rhs.blocks_[i]->clone();
    blockType_ = CoinCopyOfArray(rhs.blockType_,maximumElementBlocks_);
    if (rhs.coinModelBlocks_) {
      coinModelBlocks_ = CoinCopyOfArray(rhs.coinModelBlocks_,
					 maximumElementBlocks_);
      for (int i=0;i<numberElementBlocks_;i++)
	coinModelBlocks_[i]= new CoinModel(*rhs.coinModelBlocks_[i]);
    } else {
      coinModelBlocks_ = NULL;
    }
  } else {
    blocks_ = NULL;
    blockType_ = NULL;
    coinModelBlocks_ = NULL;
  }
  rowBlockNames_ = rhs.rowBlockNames_;
  columnBlockNames_ = rhs.columnBlockNames_;
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
CoinStructuredModel::~CoinStructuredModel ()
{
  for (int i=0;i<numberElementBlocks_;i++)
    delete blocks_[i];
  delete [] blocks_;
  delete [] blockType_;
  if (coinModelBlocks_) {
    for (int i=0;i<numberElementBlocks_;i++)
      delete coinModelBlocks_[i];
    delete [] coinModelBlocks_;
  }
}

// Clone
CoinBaseModel *
CoinStructuredModel::clone() const
{
  return new CoinStructuredModel(*this);
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CoinStructuredModel &
CoinStructuredModel::operator=(const CoinStructuredModel& rhs)
{
  if (this != &rhs) {
    CoinBaseModel::operator=(rhs);
    for (int i=0;i<numberElementBlocks_;i++)
      delete blocks_[i];
    delete [] blocks_;
    delete [] blockType_;
    if (coinModelBlocks_) {
      for (int i=0;i<numberElementBlocks_;i++)
	delete coinModelBlocks_[i];
      delete [] coinModelBlocks_;
    }
    numberRowBlocks_ = rhs.numberRowBlocks_;
    numberColumnBlocks_ = rhs.numberColumnBlocks_;
    numberElementBlocks_ = rhs.numberElementBlocks_;
    maximumElementBlocks_ = rhs.maximumElementBlocks_;
    if (maximumElementBlocks_) {
      blocks_ = CoinCopyOfArray(rhs.blocks_,maximumElementBlocks_);
      for (int i=0;i<numberElementBlocks_;i++)
	blocks_[i]= rhs.blocks_[i]->clone();
      blockType_ = CoinCopyOfArray(rhs.blockType_,maximumElementBlocks_);
      if (rhs.coinModelBlocks_) {
	coinModelBlocks_ = CoinCopyOfArray(rhs.coinModelBlocks_,
					   maximumElementBlocks_);
	for (int i=0;i<numberElementBlocks_;i++)
	  coinModelBlocks_[i]= new CoinModel(*rhs.coinModelBlocks_[i]);
      } else {
	coinModelBlocks_ = NULL;
      }
    } else {
      blocks_ = NULL;
      blockType_ = NULL;
      coinModelBlocks_ = NULL;
    }
    rowBlockNames_ = rhs.rowBlockNames_;
    columnBlockNames_ = rhs.columnBlockNames_;
  }
  return *this;
}
static bool sameValues(const double * a, const double *b, int n)
{
  int i;
  for ( i=0;i<n;i++) {
    if (a[i]!=b[i])
      break;
  }
  return (i==n);
}
static bool sameValues(const int * a, const int *b, int n)
{
  int i;
  for ( i=0;i<n;i++) {
    if (a[i]!=b[i])
      break;
  }
  return (i==n);
}
static bool sameValues(const CoinModel * a, const CoinModel *b, bool doRows)
{
  int i=0;
  int n=0;
  if (doRows) {
    n = a->numberRows();
    for ( i=0;i<n;i++) {
      const char * aName = a->getRowName(i);
      const char * bName = b->getRowName(i);
      bool good=true;
      if (aName) {
	if (!bName||strcmp(aName,bName))
	  good=false;
      } else if (bName) {
	good=false;
      }
      if (!good)
	break;
    }
  } else {
    n = a->numberColumns();
    for ( i=0;i<n;i++) {
      const char * aName = a->getColumnName(i);
      const char * bName = b->getColumnName(i);
      bool good=true;
      if (aName) {
	if (!bName||strcmp(aName,bName))
	  good=false;
      } else if (bName) {
	good=false;
      }
      if (!good)
	break;
    }
  }
  return (i==n);
}
// Add a row block name and number of rows
int
CoinStructuredModel::addRowBlock(int numberRows,const std::string &name)
{
  int iRowBlock;
  for (iRowBlock=0;iRowBlock<numberRowBlocks_;iRowBlock++) {
    if (name==rowBlockNames_[iRowBlock])
      break;
  }
  if (iRowBlock==numberRowBlocks_) {
    rowBlockNames_.push_back(name);
    numberRowBlocks_++;
    numberRows_ += numberRows;
  }
  return iRowBlock;
}
// Return a row block index given a row block name
int
CoinStructuredModel::rowBlock(const std::string &name) const
{
  int iRowBlock;
  for (iRowBlock=0;iRowBlock<numberRowBlocks_;iRowBlock++) {
    if (name==rowBlockNames_[iRowBlock])
      break;
  }
  if (iRowBlock==numberRowBlocks_)
    iRowBlock=-1;
  return iRowBlock;
}
// Add a column block name and number of columns
int
CoinStructuredModel::addColumnBlock(int numberColumns,const std::string &name)
{
  int iColumnBlock;
  for (iColumnBlock=0;iColumnBlock<numberColumnBlocks_;iColumnBlock++) {
    if (name==columnBlockNames_[iColumnBlock])
      break;
  }
  if (iColumnBlock==numberColumnBlocks_) {
    columnBlockNames_.push_back(name);
    numberColumnBlocks_++;
    numberColumns_ += numberColumns;
  }
  return iColumnBlock;
}
// Return a column block index given a column block name
int
CoinStructuredModel::columnBlock(const std::string &name) const
{
  int iColumnBlock;
  for (iColumnBlock=0;iColumnBlock<numberColumnBlocks_;iColumnBlock++) {
    if (name==columnBlockNames_[iColumnBlock])
      break;
  }
  if (iColumnBlock==numberColumnBlocks_)
    iColumnBlock=-1;
  return iColumnBlock;
}
// Return number of elements
CoinBigIndex
CoinStructuredModel::numberElements() const
{
  CoinBigIndex numberElements=0;
  for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) {
    numberElements += blocks_[iBlock]->numberElements();
  }
  return numberElements;
}
// Return i'th block as CoinModel (or NULL)
CoinModel *
CoinStructuredModel::coinBlock(int i) const
{
  CoinModel * block = dynamic_cast<CoinModel *>(blocks_[i]);
  if (block)
    return block;
  else if (coinModelBlocks_)
    return coinModelBlocks_[i];
  else
    return NULL;
}

/* Return block as a CoinModel block
   and fill in info structure and update counts
*/
CoinModel *
CoinStructuredModel::coinModelBlock(CoinModelBlockInfo & info)
{
  // CoinStructuredModel
  int numberRows=this->numberRows();
  int numberColumns =this-> numberColumns();
  int numberRowBlocks=this->numberRowBlocks();
  int numberColumnBlocks =this-> numberColumnBlocks();
  int numberElementBlocks =this-> numberElementBlocks();
  CoinBigIndex numberElements=this->numberElements();
  // See what is needed
  double * rowLower = NULL;
  double * rowUpper = NULL;
  double * columnLower = NULL;
  double * columnUpper = NULL;
  double * objective = NULL;
  int * integerType = NULL;
  info = CoinModelBlockInfo();
  CoinModel ** blocks = new CoinModel * [numberElementBlocks];
  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
    CoinModelBlockInfo thisInfo=blockType_[iBlock];
    CoinStructuredModel * subModel = dynamic_cast<CoinStructuredModel *>(blocks_[iBlock]);
    CoinModel * thisBlock;
    if (subModel) {
      thisBlock = subModel->coinModelBlock(thisInfo);
      fillInfo(thisInfo,subModel);
      setCoinModel(thisBlock,iBlock);
    } else {
      thisBlock = dynamic_cast<CoinModel *>(blocks_[iBlock]);
      assert (thisBlock);
      fillInfo(thisInfo,thisBlock);
    }
    blocks[iBlock]=thisBlock;
    if (thisInfo.rhs&&!info.rhs) {
      info.rhs=1;
      rowLower = new double [numberRows];
      rowUpper = new double [numberRows];
      CoinFillN(rowLower,numberRows,-COIN_DBL_MAX);
      CoinFillN(rowUpper,numberRows,COIN_DBL_MAX);
    }
    if (thisInfo.bounds&&!info.bounds) {
      info.bounds=1;
      columnLower = new double [numberColumns];
      columnUpper = new double [numberColumns];
      objective = new double [numberColumns];
      CoinFillN(columnLower,numberColumns,0.0);
      CoinFillN(columnUpper,numberColumns,COIN_DBL_MAX);
      CoinFillN(objective,numberColumns,0.0);
    }
    if (thisInfo.integer&&!info.integer) {
      info.integer=1;
      integerType = new int [numberColumns];
      CoinFillN(integerType,numberColumns,0);
    }
    if (thisInfo.rowName&&!info.rowName) {
      info.rowName=1;
    }
    if (thisInfo.columnName&&!info.columnName) {
      info.columnName=1;
    }
  }
  // Space for elements
  int * row = new int[numberElements];
  int * column = new int[numberElements];
  double * element = new double[numberElements];
  numberElements=0;
  // Bases for blocks
  int * rowBase = new int[numberRowBlocks];
  CoinFillN(rowBase,numberRowBlocks,-1);
  CoinModelBlockInfo * rowBlockInfo =
    new CoinModelBlockInfo [numberRowBlocks];
  int * columnBase = new int[numberColumnBlocks];
  CoinFillN(columnBase,numberColumnBlocks,-1);
  CoinModelBlockInfo * columnBlockInfo =
    new CoinModelBlockInfo [numberColumnBlocks];
  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
    int iRowBlock = rowBlock(blocks[iBlock]->getRowBlock());
    assert (iRowBlock>=0&&iRowBlock<numberRowBlocks);
    if (rowBase[iRowBlock]==-1)
      rowBase[iRowBlock]=blocks[iBlock]->numberRows();
    else
      assert (rowBase[iRowBlock]==blocks[iBlock]->numberRows());
    int iColumnBlock = columnBlock(blocks[iBlock]->getColumnBlock());
    assert (iColumnBlock>=0&&iColumnBlock<numberColumnBlocks);
    if (columnBase[iColumnBlock]==-1)
      columnBase[iColumnBlock]=blocks[iBlock]->numberColumns();
    else
      assert (columnBase[iColumnBlock]==blocks[iBlock]->numberColumns());
  }
  int n=0;
  for (int iBlock=0;iBlock<numberRowBlocks;iBlock++) {
    int k = rowBase[iBlock];
    rowBase[iBlock]=n;
    assert (k>=0);
    n+=k;
  }
  assert (n==numberRows);
  n=0;
  for (int iBlock=0;iBlock<numberColumnBlocks;iBlock++) {
    int k = columnBase[iBlock];
    columnBase[iBlock]=n;
    assert (k>=0);
    n+=k;
  }
  assert (n==numberColumns);
  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
    CoinModelBlockInfo thisInfo=blockType_[iBlock];
    CoinModel * thisBlock = blocks[iBlock];
    int iRowBlock = rowBlock(blocks[iBlock]->getRowBlock());
    int iRowBase = rowBase[iRowBlock];
    int nRows = thisBlock->numberRows();
    // could check if identical before error
    if (thisInfo.rhs) {
      assert (!rowBlockInfo[iRowBlock].rhs);
      rowBlockInfo[iRowBlock].rhs=1;
      memcpy(rowLower+iRowBase,thisBlock->rowLowerArray(),
	     nRows*sizeof(double));
      memcpy(rowUpper+iRowBase,thisBlock->rowUpperArray(),
	     nRows*sizeof(double));
    }
    int iColumnBlock = columnBlock(blocks[iBlock]->getColumnBlock());
    int iColumnBase = columnBase[iColumnBlock];
    int nColumns = thisBlock->numberColumns();
    if (thisInfo.bounds) {
      assert (!columnBlockInfo[iColumnBlock].bounds);
      columnBlockInfo[iColumnBlock].bounds=1;
      memcpy(columnLower+iColumnBase,thisBlock->columnLowerArray(),
	     nColumns*sizeof(double));
      memcpy(columnUpper+iColumnBase,thisBlock->columnUpperArray(),
	     nColumns*sizeof(double));
      memcpy(objective+iColumnBase,thisBlock->objectiveArray(),
	     nColumns*sizeof(double));
    }
    if (thisInfo.integer) {
      assert (!columnBlockInfo[iColumnBlock].integer);
      columnBlockInfo[iColumnBlock].integer=1;
      memcpy(integerType+iColumnBase,thisBlock->integerTypeArray(),
	     nColumns*sizeof(int));
    }
    const CoinPackedMatrix * elementBlock = thisBlock->packedMatrix();
    // get matrix data pointers
    const int * row2 = elementBlock->getIndices();
    const CoinBigIndex * columnStart = elementBlock->getVectorStarts();
    const double * elementByColumn = elementBlock->getElements();
    const int * columnLength = elementBlock->getVectorLengths();
    int n = elementBlock->getNumCols();
    assert (elementBlock->isColOrdered());
    for (int iColumn=0;iColumn<n;iColumn++) {
      CoinBigIndex j;
      int jColumn = iColumn+iColumnBase;
      for (j=columnStart[iColumn];
	   j<columnStart[iColumn]+columnLength[iColumn];j++) {
	row[numberElements]=row2[j]+iRowBase;
	column[numberElements]=jColumn;
	element[numberElements++]=elementByColumn[j];
	}
    }
  }
  delete [] rowBlockInfo;
  delete [] columnBlockInfo;
  CoinPackedMatrix matrix(true,row,column,element,numberElements);
  if (numberElements)
    info.matrix=1;
  delete [] row;
  delete [] column;
  delete [] element;
  CoinModel * block = new CoinModel(numberRows, numberColumns,
			&matrix,
			rowLower, rowUpper,
			columnLower, columnUpper, objective);
  delete [] rowLower;
  delete [] rowUpper;
  delete [] columnLower;
  delete [] columnUpper;
  delete [] objective;
  // Do integers if wanted
  if (integerType) {
    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      block->setColumnIsInteger(iColumn,integerType[iColumn]!=0);
    }
  }
  delete [] integerType;
  block->setObjectiveOffset(objectiveOffset());
  if (info.rowName||info.columnName) {
    for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
      CoinModelBlockInfo thisInfo;
      CoinModel * thisBlock = blocks[iBlock];
      int iRowBlock = rowBlock(thisBlock->getRowBlock());
      int iRowBase = rowBase[iRowBlock];
      if (thisInfo.rowName) {
	int numberItems = thisBlock->rowNames()->numberItems();
	assert( thisBlock->numberRows()>=numberItems);
	if (numberItems) {
	  const char *const * rowNames=thisBlock->rowNames()->names();
	  for (int i=0;i<numberItems;i++) {
	    std::string name = rowNames[i];
	    block->setRowName(i+iRowBase,name.c_str());
	  }
	}
      }
      int iColumnBlock = columnBlock(thisBlock->getColumnBlock());
      int iColumnBase = columnBase[iColumnBlock];
      if (thisInfo.columnName) {
	int numberItems = thisBlock->columnNames()->numberItems();
	assert( thisBlock->numberColumns()>=numberItems);
	if (numberItems) {
	  const char *const * columnNames=thisBlock->columnNames()->names();
	  for (int i=0;i<numberItems;i++) {
	    std::string name = columnNames[i];
	    block->setColumnName(i+iColumnBase,name.c_str());
	  }
	}
      }
    }
  }
  delete [] rowBase;
  delete [] columnBase;
  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
    if (static_cast<CoinBaseModel *> (blocks[iBlock])!=
	static_cast<CoinBaseModel *> (blocks_[iBlock]))
      delete blocks[iBlock];
  }
  delete [] blocks;
  return block;
}
// Sets given block into coinModelBlocks_
void
CoinStructuredModel::setCoinModel(CoinModel * block, int iBlock)
{
 if (!coinModelBlocks_) {
    coinModelBlocks_ = new CoinModel * [maximumElementBlocks_];
    CoinZeroN(coinModelBlocks_,maximumElementBlocks_);
  }
  delete coinModelBlocks_[iBlock];
  coinModelBlocks_[iBlock]=block;
 }
// Refresh info in blockType_
void
CoinStructuredModel::refresh(int iBlock)
{
  fillInfo(blockType_[iBlock],coinBlock(iBlock));
}
/* Fill in info structure and update counts
   Returns number of inconsistencies on border
*/
int
CoinStructuredModel::fillInfo(CoinModelBlockInfo & info,
			      const CoinModel * block)
{
  int whatsSet = block->whatIsSet();
  info.matrix = static_cast<char>(((whatsSet&1)!=0) ? 1 : 0);
  info.rhs = static_cast<char>(((whatsSet&2)!=0) ? 1 : 0);
  info.rowName = static_cast<char>(((whatsSet&4)!=0) ? 1 : 0);
  info.integer = static_cast<char>(((whatsSet&32)!=0) ? 1 : 0);
  info.bounds = static_cast<char>(((whatsSet&8)!=0) ? 1 : 0);
  info.columnName = static_cast<char>(((whatsSet&16)!=0) ? 1 : 0);
  int numberRows = block->numberRows();
  int numberColumns = block->numberColumns();
  // Which block
  int iRowBlock=addRowBlock(numberRows,block->getRowBlock());
  info.rowBlock=iRowBlock;
  int iColumnBlock=addColumnBlock(numberColumns,block->getColumnBlock());
  info.columnBlock=iColumnBlock;
  int numberErrors=0;
  CoinModelBlockInfo sumInfo=blockType_[numberElementBlocks_-1];
  int iRhs=(sumInfo.rhs) ? numberElementBlocks_-1 : -1;
  int iRowName=(sumInfo.rowName) ? numberElementBlocks_-1 : -1;
  int iBounds=(sumInfo.bounds) ? numberElementBlocks_-1 : -1;
  int iColumnName=(sumInfo.columnName) ? numberElementBlocks_-1 : -1;
  int iInteger=(sumInfo.integer) ? numberElementBlocks_-1 : -1;
  for (int i=0;i<numberElementBlocks_-1;i++) {
    if (iRowBlock==blockType_[i].rowBlock) {
      if (numberRows!=blocks_[i]->numberRows())
	numberErrors+=1000;
      if (blockType_[i].rhs) {
	if (iRhs<0) {
	  iRhs=i;
	} else {
	  // check
	  const double * a = static_cast<CoinModel *>(blocks_[iRhs])->rowLowerArray();
	  const double * b = static_cast<CoinModel *>(blocks_[i])->rowLowerArray();
	  if (!sameValues(a,b,numberRows))
	    numberErrors++;
	  a = static_cast<CoinModel *>(blocks_[iRhs])->rowUpperArray();
	  b = static_cast<CoinModel *>(blocks_[i])->rowUpperArray();
	  if (!sameValues(a,b,numberRows))
	    numberErrors++;
	}
      }
      if (blockType_[i].rowName) {
	if (iRowName<0) {
	  iRowName=i;
	} else {
	  // check
	  if (!sameValues(static_cast<CoinModel *>(blocks_[iRowName]),
			  static_cast<CoinModel *>(blocks_[i]),true))
	    numberErrors++;
	}
      }
    }
    if (iColumnBlock==blockType_[i].columnBlock) {
      if (numberColumns!=blocks_[i]->numberColumns())
	numberErrors+=1000;
      if (blockType_[i].bounds) {
	if (iBounds<0) {
	  iBounds=i;
	} else {
	  // check
	  const double * a = static_cast<CoinModel *>(blocks_[iBounds])->columnLowerArray();
	  const double * b = static_cast<CoinModel *>(blocks_[i])->columnLowerArray();
	  if (!sameValues(a,b,numberColumns))
	    numberErrors++;
	  a = static_cast<CoinModel *>(blocks_[iBounds])->columnUpperArray();
	  b = static_cast<CoinModel *>(blocks_[i])->columnUpperArray();
	  if (!sameValues(a,b,numberColumns))
	    numberErrors++;
	  a = static_cast<CoinModel *>(blocks_[iBounds])->objectiveArray();
	  b = static_cast<CoinModel *>(blocks_[i])->objectiveArray();
	  if (!sameValues(a,b,numberColumns))
	    numberErrors++;
	}
      }
      if (blockType_[i].columnName) {
	if (iColumnName<0) {
	  iColumnName=i;
	} else {
	  // check
	  if (!sameValues(static_cast<CoinModel *>(blocks_[iColumnName]),
			  static_cast<CoinModel *>(blocks_[i]),false))
	    numberErrors++;
	}
      }
      if (blockType_[i].integer) {
	if (iInteger<0) {
	  iInteger=i;
	} else {
	  // check
	  const int * a = static_cast<CoinModel *>(blocks_[iInteger])->integerTypeArray();
	  const int * b = static_cast<CoinModel *>(blocks_[i])->integerTypeArray();
	  if (!sameValues(a,b,numberColumns))
	    numberErrors++;
	}
      }
    }
  }
  return numberErrors;
}
/* Fill in info structure and update counts
*/
void
CoinStructuredModel::fillInfo(CoinModelBlockInfo & info,
			      const CoinStructuredModel * block)
{
  int numberRows = block->numberRows();
  int numberColumns = block->numberColumns();
  // Which block
  int iRowBlock=addRowBlock(numberRows,block->rowBlockName_);
  info.rowBlock=iRowBlock;
  int iColumnBlock=addColumnBlock(numberColumns,block->columnBlockName_);
  info.columnBlock=iColumnBlock;
}
/* add a block from a CoinModel without names*/
int
CoinStructuredModel::addBlock(const std::string & rowBlock,
			      const std::string & columnBlock,
			      const CoinBaseModel & block)
{
  CoinBaseModel * block2 = block.clone();
  return addBlock(rowBlock,columnBlock,block2);
}
/* add a block using names
 */
int
CoinStructuredModel::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)
{
  CoinModel * block = new CoinModel();
  block->loadBlock(matrix,columnLower,columnUpper,objective,
		   rowLower,rowUpper);
  return addBlock(rowBlock,columnBlock,block);
}
/* add a block from a CoinModel without names*/
int
CoinStructuredModel::addBlock(const std::string & rowBlock,
			      const std::string & columnBlock,
			      CoinBaseModel * block)
{
  if (numberElementBlocks_==maximumElementBlocks_) {
    maximumElementBlocks_ = 3*(maximumElementBlocks_+10)/2;
    CoinBaseModel ** temp = new CoinBaseModel * [maximumElementBlocks_];
    memcpy(temp,blocks_,numberElementBlocks_*sizeof(CoinBaseModel *));
    delete [] blocks_;
    blocks_ = temp;
    CoinModelBlockInfo * temp2 = new CoinModelBlockInfo [maximumElementBlocks_];
    memcpy(temp2,blockType_,numberElementBlocks_*sizeof(CoinModelBlockInfo));
    delete [] blockType_;
    blockType_ = temp2;
    if (coinModelBlocks_) {
      CoinModel ** temp = new CoinModel * [maximumElementBlocks_];
      CoinZeroN(temp,maximumElementBlocks_);
      memcpy(temp,coinModelBlocks_,numberElementBlocks_*sizeof(CoinModel *));
      delete [] coinModelBlocks_;
      coinModelBlocks_ = temp;
    }
  }
  blocks_[numberElementBlocks_++]=block;
  block->setRowBlock(rowBlock);
  block->setColumnBlock(columnBlock);
  int numberErrors=0;
  CoinModel * coinBlock = dynamic_cast<CoinModel *>(block);
  if (coinBlock) {
    // Convert matrix
    if (coinBlock->type()!=3)
      coinBlock->convertMatrix();
    numberErrors=fillInfo(blockType_[numberElementBlocks_-1],coinBlock);
  } else {
    CoinStructuredModel * subModel = dynamic_cast<CoinStructuredModel *>(block);
    assert (subModel);
    CoinModel * blockX =
      subModel->coinModelBlock(blockType_[numberElementBlocks_-1]);
    fillInfo(blockType_[numberElementBlocks_-1],subModel);
    setCoinModel(blockX,numberElementBlocks_-1);
  }
  return numberErrors;
}

/* add a block from a CoinModel with names*/
int
CoinStructuredModel::addBlock(const CoinBaseModel & block)
{

  //inline const std::string & getRowBlock() const
  //abort();
  return addBlock(block.getRowBlock(),block.getColumnBlock(),
		  block);
}
/* Decompose a model specified as arrays + CoinPackedMatrix
   1 - try D-W
   2 - try Benders
   3 - try Staircase
   Returns number of blocks or zero if no structure
*/
int
CoinStructuredModel::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)
{
  setObjectiveOffset(objectiveOffset);
  int numberBlocks=0;
  if (type==1) {
    // Try master at top and bottom
    bool goodDW=true;
    // get row copy
    CoinPackedMatrix rowCopy = matrix;
    rowCopy.reverseOrdering();
    const int * row = matrix.getIndices();
    const int * columnLength = matrix.getVectorLengths();
    const CoinBigIndex * columnStart = matrix.getVectorStarts();
    //const double * elementByColumn = matrix.getElements();
    const int * column = rowCopy.getIndices();
    const int * rowLength = rowCopy.getVectorLengths();
    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
    //const double * elementByRow = rowCopy.getElements();
    int numberRows = matrix.getNumRows();
    int * rowBlock = new int[numberRows+1];
    int iRow;
    // Row counts (maybe look at long rows for master)
    CoinZeroN(rowBlock,numberRows+1);
    for (iRow=0;iRow<numberRows;iRow++) {
      int length = rowLength[iRow];
      rowBlock[length]++;
    }
    for (iRow=0;iRow<numberRows+1;iRow++) {
      if (rowBlock[iRow])
	printf("%d rows have %d elements\n",rowBlock[iRow],iRow);
    }
    bool newWay=true;
    // to say if column looked at
    int numberColumns = matrix.getNumCols();
    int * columnBlock = new int[numberColumns];
    int iColumn;
    int * whichRow = new int [numberRows];
    int * whichColumn = new int [numberColumns];
    int * stack = new int [numberRows];
    if (newWay) {
      //double best2[3]={COIN_DBL_MAX,COIN_DBL_MAX,COIN_DBL_MAX};
      double best2[3]={0.0,0.0,0.0};
      int row2[3]={-1,-1,-1};
      // try forward and backward and sorted
      for (int iWay=0;iWay<3;iWay++) {
	if (iWay==0) {
	  // forwards
	  for (int i=0;i<numberRows;i++)
	    stack[i]=i;
	} else if (iWay==1) {
	  // backwards
	  for (int i=0;i<numberRows;i++)
	    stack[i]=numberRows-1-i;
	} else {
	  // sparsest first
	  for (int i=0;i<numberRows;i++) {
	    rowBlock[i]=rowLength[i];
	    stack[i]=i;
	  }
	  CoinSort_2(rowBlock,rowBlock+numberRows,stack);
	}
	CoinFillN(rowBlock,numberRows,-1);
	rowBlock[numberRows]=0;
	CoinFillN(whichRow,numberRows,-1);
	CoinFillN(columnBlock,numberColumns,-1);
	CoinFillN(whichColumn,numberColumns,-1);
	int numberMarkedColumns = 0;
	numberBlocks=0;
	int bestRow = -1;
	int maximumInBlock = 0;
	int rowsDone=0;
	int checkAfterRows = (5*numberRows)/10+1;
#define OSL_WAY
#ifdef OSL_WAY
	double best = COIN_DBL_MAX;
	int bestRowsDone=-1;
#else
	double best = 0.0; //COIN_DBL_MAX;
#endif
	int numberGoodBlocks=0;

	for (int kRow=0;kRow<numberRows;kRow++) {
	  iRow = stack[kRow];
	  CoinBigIndex start = rowStart[iRow];
	  CoinBigIndex end = start+rowLength[iRow];
	  int iBlock=-1;
	  for (CoinBigIndex j=start;j<end;j++) {
	    int iColumn = column[j];
	    if (columnBlock[iColumn]>=0) {
	      // already marked
	      if (iBlock<0) {
		iBlock = columnBlock[iColumn];
	      } else if (iBlock != columnBlock[iColumn]) {
		// join two blocks
		int jBlock = columnBlock[iColumn];
		numberGoodBlocks--;
		// Increase count of iBlock
		rowBlock[iBlock] += rowBlock[jBlock];
		rowBlock[jBlock]=0;
		// First column of block jBlock
		int jColumn = whichRow[jBlock];
		while (jColumn>=0) {
		  columnBlock[jColumn]=iBlock;
		  iColumn = jColumn;
		  jColumn = whichColumn[jColumn];
		}
		whichColumn[iColumn] = whichRow[iBlock];
		whichRow[iBlock] = whichRow[jBlock];
		whichRow[jBlock]=-1;
	      }
	    }
	  }
	  int n=end-start;
	  // If not in block - then start one
	  if (iBlock<0) {
	    // unless null row
	    if (n) {
	      iBlock = numberBlocks;
	      numberBlocks++;
	      numberGoodBlocks++;
	      int jColumn = column[start];
	      columnBlock[jColumn]=iBlock;
	      whichRow[iBlock]=jColumn;
	      numberMarkedColumns += n;
	      rowBlock[iBlock] = n;
	      for (CoinBigIndex j=start+1;j<end;j++) {
		int iColumn = column[j];
		columnBlock[iColumn]=iBlock;
		whichColumn[jColumn]=iColumn;
		jColumn=iColumn;
	      }
	    } else {
	      rowBlock[numberRows]++;
	    }
	  } else {
	    // add all to this block if not already in
	    int jColumn = whichRow[iBlock];
	    for (CoinBigIndex j=start;j<end;j++) {
	      int iColumn = column[j];
	      if (columnBlock[iColumn]<0) {
		numberMarkedColumns++;
		rowBlock[iBlock]++;
		columnBlock[iColumn]=iBlock;
		whichColumn[iColumn]=jColumn;
		jColumn=iColumn;
	      }
	    }
	    whichRow[iBlock]=jColumn;
	  }
#if 0
	  {
	    int nn=0;
	    int * temp = new int [numberRows];
	    CoinZeroN(temp,numberRows);
	    for (int i=0;i<numberColumns;i++) {
	      int iBlock = columnBlock[i];
	      if (iBlock>=0) {
		nn++;
		assert (iBlock<numberBlocks);
		temp[iBlock]++;
	      }
	    }
	    for (int i=0;i<numberBlocks;i++)
	      assert (temp[i]==rowBlock[i]);
	    assert (nn==numberMarkedColumns);
	    for (int i=0;i<numberBlocks;i++) {
	      // First column of block i
	      int jColumn = whichRow[i];
	      int n=0;
	      while (jColumn>=0) {
		n++;
		jColumn = whichColumn[jColumn];
	      }
	      assert (n==temp[i]);
	    }
	    delete [] temp;
	  }
#endif
	  rowsDone++;
	  if (iBlock>=0)
	    maximumInBlock = CoinMax(maximumInBlock,rowBlock[iBlock]);
	  if (rowsDone>=checkAfterRows) {
	    assert (numberGoodBlocks>0);
	    double averageSize = static_cast<double>(numberMarkedColumns)/
	      static_cast<double>(numberGoodBlocks);
#ifndef OSL_WAY
	    double notionalBlocks = static_cast<double>(numberMarkedColumns)/
	      averageSize;
	    if (maximumInBlock<3*averageSize&&numberGoodBlocks>2) {
	      if(best*(numberRows-rowsDone) < notionalBlocks) {
		best = notionalBlocks/
		  static_cast<double> (numberRows-rowsDone);
		bestRow = kRow;
	      }
	    }
#else
	    if (maximumInBlock*10<numberColumns*11&&numberGoodBlocks>1) {
	      double test = maximumInBlock + 0.0*averageSize;
	      if(best*static_cast<double>(rowsDone) > test) {
		best = test/static_cast<double> (rowsDone);
		bestRow = kRow;
		bestRowsDone=rowsDone;
	      }
	    }
#endif
	  }
	}
#ifndef OSL_WAY
	best2[iWay]=best;
#else
	if (bestRowsDone<numberRows)
	  best2[iWay]=-(numberRows-bestRowsDone);
	else
	  best2[iWay]=-numberRows;
#endif
	row2[iWay]=bestRow;
      }
      // mark rows
      int nMaster;
      CoinFillN(rowBlock,numberRows,-2);
      if (best2[2]<best2[0]||best2[2]<best2[1]) {
	int iRow1;
	int iRow2;
	if (best2[0]>best2[1]) {
	  // Bottom rows in master
	  iRow1=row2[0]+1;
	  iRow2=numberRows;
	} else {
	  // Top rows in master
	  iRow1=0;
	  iRow2=numberRows-row2[1];
	}
	nMaster = iRow2-iRow1;
	CoinFillN(rowBlock+iRow1,nMaster,-1);
      } else {
	// sorted
	// Bottom rows in master (in order)
	int iRow1=row2[2]+1;
	nMaster = numberRows-iRow1;
	for (int i=iRow1;i<numberRows;i++)
	  rowBlock[stack[i]]=-1;
      }
      if (nMaster*2>numberRows) {
	goodDW=false;
	printf("%d rows out of %d would be in master - no good\n",
	       nMaster,numberRows);
	delete [] rowBlock;
	delete [] columnBlock;
	delete [] whichRow;
	delete [] whichColumn;
	delete [] stack;
	CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper,
			columnLower,columnUpper,objective);
	model.setObjectiveOffset(objectiveOffset);
	addBlock("row_master","column_master",model);
	return 0;
      }
    } else {
      for (iRow=0;iRow<numberRows;iRow++)
	rowBlock[iRow]=-2;
      // these are master rows
      if (numberRows==105127) {
	// ken-18
	for (iRow=104976;iRow<numberRows;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==2426) {
	// ken-7
	for (iRow=2401;iRow<numberRows;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==810) {
	for (iRow=81;iRow<84;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==5418) {
	for (iRow=564;iRow<603;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==10280) {
	// osa-60
	for (iRow=10198;iRow<10280;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==1503) {
	// degen3
	for (iRow=0;iRow<561;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==929) {
	// czprob
	for (iRow=0;iRow<39;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==33874) {
	// pds-20
	for (iRow=31427;iRow<33874;iRow++)
	  rowBlock[iRow]=-1;
      } else if (numberRows==24902) {
	// allgrade
	int kRow=818;
	for (iRow=0;iRow<kRow;iRow++)
	  rowBlock[iRow]=-1;
      }
    }
    numberBlocks=0;
    CoinFillN(columnBlock,numberColumns,-2);
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int kstart = columnStart[iColumn];
      int kend = columnStart[iColumn]+columnLength[iColumn];
      if (columnBlock[iColumn]==-2) {
	// column not allocated
	int j;
	int nstack=0;
	for (j=kstart;j<kend;j++) {
	  int iRow= row[j];
	  if (rowBlock[iRow]!=-1) {
	    assert(rowBlock[iRow]==-2);
	    rowBlock[iRow]=numberBlocks; // mark
	    stack[nstack++] = iRow;
	  }
	}
	if (nstack) {
	  // new block - put all connected in
	  numberBlocks++;
	  columnBlock[iColumn]=numberBlocks-1;
	  while (nstack) {
	    int iRow = stack[--nstack];
	    int k;
	    for (k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {
	      int iColumn = column[k];
	      int kkstart = columnStart[iColumn];
	      int kkend = kkstart + columnLength[iColumn];
	      if (columnBlock[iColumn]==-2) {
		columnBlock[iColumn]=numberBlocks-1; // mark
		// column not allocated
		int jj;
		for (jj=kkstart;jj<kkend;jj++) {
		  int jRow= row[jj];
		  if (rowBlock[jRow]==-2) {
		    rowBlock[jRow]=numberBlocks-1;
		    stack[nstack++]=jRow;
		  }
		}
	      } else {
		assert (columnBlock[iColumn]==numberBlocks-1);
	      }
	    }
	  }
	} else {
	  // Only in master
	  columnBlock[iColumn]=-1;
	}
      }
    }
    delete [] stack;
    int numberMasterRows=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock==-1)
	numberMasterRows++;
    }
    int numberMasterColumns=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock==-1)
	numberMasterColumns++;
    }
    if (numberBlocks<=maxBlocks)
      printf("%d blocks found - %d rows, %d columns in master\n",
	     numberBlocks,numberMasterRows,numberMasterColumns);
    else
      printf("%d blocks found (reduced to %d) - %d rows, %d columns in master\n",
	     numberBlocks,maxBlocks,numberMasterRows,numberMasterColumns);
    if (numberBlocks) {
      if (numberBlocks>maxBlocks) {
	int iBlock;
	for (iRow=0;iRow<numberRows;iRow++) {
	  iBlock = rowBlock[iRow];
	  if (iBlock>=0)
	    rowBlock[iRow] = iBlock%maxBlocks;
	}
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  iBlock = columnBlock[iColumn];
	  if (iBlock>=0)
	    columnBlock[iColumn] = iBlock%maxBlocks;
	}
	numberBlocks=maxBlocks;
      }
    }
    // make up problems
    // Create all sub problems
    // Space for creating
    double * obj = new double [numberColumns];
    double * columnLo = new double [numberColumns];
    double * columnUp = new double [numberColumns];
    double * rowLo = new double [numberRows];
    double * rowUp = new double [numberRows];
    // Counts
    int * rowCount = reinterpret_cast<int *>(rowLo);
    CoinZeroN(rowCount,numberBlocks);
    for (int i=0;i<numberRows;i++) {
      int iBlock=rowBlock[i];
      if (iBlock>=0)
	rowCount[iBlock]++;
    }
    // allocate empty rows
    for (int i=0;i<numberRows;i++) {
      int iBlock=rowBlock[i];
      if (iBlock==-2) {
	// find block with smallest count
	int iSmall=-1;
	int smallest = numberRows;
	for (int j=0;j<numberBlocks;j++) {
	  if (rowCount[j]<smallest) {
	    iSmall=j;
	    smallest=rowCount[j];
	  }
	}
	rowBlock[i]=iSmall;
	rowCount[iSmall]++;
      }
    }
    int * columnCount = reinterpret_cast<int *>(rowUp);
    CoinZeroN(columnCount,numberBlocks);
    for (int i=0;i<numberColumns;i++) {
      int iBlock=columnBlock[i];
      if (iBlock>=0)
	columnCount[iBlock]++;
    }
    int maximumSize=0;
    for (int i=0;i<numberBlocks;i++) {
      printf("Block %d has %d rows and %d columns\n",i,
	     rowCount[i],columnCount[i]);
      int k=2*rowCount[i]+columnCount[i];
      maximumSize = CoinMax(maximumSize,k);
    }
    if (maximumSize*10>4*(2*numberRows+numberColumns)) {
      // No good
      printf("Doesn't look good\n");
      delete [] rowBlock;
      delete [] columnBlock;
      delete [] whichRow;
      delete [] whichColumn;
      delete [] obj ;
      delete [] columnLo ;
      delete [] columnUp ;
      delete [] rowLo ;
      delete [] rowUp ;
      CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper,
		      columnLower,columnUpper,objective);
      model.setObjectiveOffset(objectiveOffset);
      addBlock("row_master","column_master",model);
      return 0;
    }
    // Name for master so at top
    addRowBlock(numberMasterRows,"row_master");
    // Arrays
    // get full matrix
    CoinPackedMatrix fullMatrix = matrix;
    int numberRow2,numberColumn2;
    int iBlock;
    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
      char rowName[20];
      sprintf(rowName,"row_%d",iBlock);
      char columnName[20];
      sprintf(columnName,"column_%d",iBlock);
      numberRow2=0;
      numberColumn2=0;
      for (iRow=0;iRow<numberRows;iRow++) {
	if (iBlock==rowBlock[iRow]) {
	  rowLo[numberRow2]=rowLower[iRow];
	  rowUp[numberRow2]=rowUpper[iRow];
	  whichRow[numberRow2++]=iRow;
	}
      }
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	if (iBlock==columnBlock[iColumn]) {
	  obj[numberColumn2]=objective[iColumn];
	  columnLo[numberColumn2]=columnLower[iColumn];
	  columnUp[numberColumn2]=columnUpper[iColumn];
	  whichColumn[numberColumn2++]=iColumn;
	}
      }
      // Diagonal block
      CoinPackedMatrix mat(fullMatrix,
			   numberRow2,whichRow,
			   numberColumn2,whichColumn);
      // make sure correct dimensions
      mat.setDimensions(numberRow2,numberColumn2);
      CoinModel * block = new CoinModel(numberRow2,numberColumn2,&mat,
					rowLo,rowUp,NULL,NULL,NULL);
      block->setOriginalIndices(whichRow,whichColumn);
      addBlock(rowName,columnName,block); // takes ownership
      // and top block
      numberRow2=0;
      // get top matrix
      for (iRow=0;iRow<numberRows;iRow++) {
	int iBlock = rowBlock[iRow];
	if (iBlock==-1) {
	  whichRow[numberRow2++]=iRow;
	}
      }
      CoinPackedMatrix top(fullMatrix,
			   numberRow2,whichRow,
			   numberColumn2,whichColumn);
      // make sure correct dimensions
      top.setDimensions(numberRow2,numberColumn2);
      block = new CoinModel(numberMasterRows,numberColumn2,&top,
			    NULL,NULL,columnLo,columnUp,obj);
      block->setOriginalIndices(whichRow,whichColumn);
      addBlock("row_master",columnName,block); // takes ownership
    }
    // and master
    numberRow2=0;
    numberColumn2=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock==-1) {
	rowLo[numberRow2]=rowLower[iRow];
	rowUp[numberRow2]=rowUpper[iRow];
	whichRow[numberRow2++]=iRow;
      }
    }
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock<0) {
	obj[numberColumn2]=objective[iColumn];
	columnLo[numberColumn2]=columnLower[iColumn];
	columnUp[numberColumn2]=columnUpper[iColumn];
	whichColumn[numberColumn2++]=iColumn;
      }
    }
    delete [] rowBlock;
    delete [] columnBlock;
    CoinPackedMatrix top(fullMatrix,
			 numberRow2,whichRow,
			 numberColumn2,whichColumn);
    // make sure correct dimensions
    top.setDimensions(numberRow2,numberColumn2);
    CoinModel * block = new CoinModel(numberRow2,numberColumn2,&top,
				      rowLo,rowUp,
				      columnLo,columnUp,obj);
    block->setOriginalIndices(whichRow,whichColumn);
    addBlock("row_master","column_master",block); // takes ownership
    delete [] whichRow;
    delete [] whichColumn;
    delete [] obj ;
    delete [] columnLo ;
    delete [] columnUp ;
    delete [] rowLo ;
    delete [] rowUp ;
  } else if (type==2) {
    // Try master at beginning and end
    bool goodBenders=true;
    // get row copy
    CoinPackedMatrix rowCopy = matrix;
    rowCopy.reverseOrdering();
    const int * row = matrix.getIndices();
    const int * columnLength = matrix.getVectorLengths();
    const CoinBigIndex * columnStart = matrix.getVectorStarts();
    //const double * elementByColumn = matrix.getElements();
    const int * column = rowCopy.getIndices();
    const int * rowLength = rowCopy.getVectorLengths();
    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
    //const double * elementByRow = rowCopy.getElements();
    int numberColumns = matrix.getNumCols();
    int * columnBlock = new int[numberColumns+1];
    int iColumn;
    // Column counts (maybe look at long columns for master)
    CoinZeroN(columnBlock,numberColumns+1);
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int length = columnLength[iColumn];
      columnBlock[length]++;
    }
    for (iColumn=0;iColumn<numberColumns+1;iColumn++) {
      if (columnBlock[iColumn])
	printf("%d columns have %d elements\n",columnBlock[iColumn],iColumn);
    }
    bool newWay=false;
    // to say if row looked at
    int numberRows = matrix.getNumRows();
    int * rowBlock = new int[numberRows];
    int iRow;
    int * whichRow = new int [numberRows];
    int * whichColumn = new int [numberColumns];
    int * stack = new int [numberColumns];
    if (newWay) {
      //double best2[3]={COIN_DBL_MAX,COIN_DBL_MAX,COIN_DBL_MAX};
      double best2[3]={0.0,0.0,0.0};
      int column2[3]={-1,-1,-1};
      // try forward and backward and sorted
      for (int iWay=0;iWay<3;iWay++) {
	if (iWay==0) {
	  // forwards
	  for (int i=0;i<numberColumns;i++)
	    stack[i]=i;
	} else if (iWay==1) {
	  // backwards
	  for (int i=0;i<numberColumns;i++)
	    stack[i]=numberColumns-1-i;
	} else {
	  // sparsest first
	  for (int i=0;i<numberColumns;i++) {
	    columnBlock[i]=columnLength[i];
	    stack[i]=i;
	  }
	  CoinSort_2(columnBlock,columnBlock+numberColumns,stack);
	}
	CoinFillN(columnBlock,numberColumns,-1);
	columnBlock[numberColumns]=0;
	CoinFillN(whichColumn,numberColumns,-1);
	CoinFillN(rowBlock,numberRows,-1);
	CoinFillN(whichRow,numberRows,-1);
	int numberMarkedRows = 0;
	numberBlocks=0;
	int bestColumn = -1;
	int maximumInBlock = 0;
	int columnsDone=0;
	int checkAfterColumns = (5*numberColumns)/10+1;
#ifdef OSL_WAY
	double best = COIN_DBL_MAX;
	int bestColumnsDone=-1;
#else
	double best = 0.0; //COIN_DBL_MAX;
#endif
	int numberGoodBlocks=0;

	for (int kColumn=0;kColumn<numberColumns;kColumn++) {
	  iColumn = stack[kColumn];
	  CoinBigIndex start = columnStart[iColumn];
	  CoinBigIndex end = start+columnLength[iColumn];
	  int iBlock=-1;
	  for (CoinBigIndex j=start;j<end;j++) {
	    int iRow = row[j];
	    if (rowBlock[iRow]>=0) {
	      // already marked
	      if (iBlock<0) {
		iBlock = rowBlock[iRow];
	      } else if (iBlock != rowBlock[iRow]) {
		// join two blocks
		int jBlock = rowBlock[iRow];
		numberGoodBlocks--;
		// Increase count of iBlock
		columnBlock[iBlock] += columnBlock[jBlock];
		columnBlock[jBlock]=0;
		// First row of block jBlock
		int jRow = whichColumn[jBlock];
		while (jRow>=0) {
		  rowBlock[jRow]=iBlock;
		  iRow = jRow;
		  jRow = whichRow[jRow];
		}
		whichRow[iRow] = whichColumn[iBlock];
		whichColumn[iBlock] = whichColumn[jBlock];
		whichColumn[jBlock]=-1;
	      }
	    }
	  }
	  int n=end-start;
	  // If not in block - then start one
	  if (iBlock<0) {
	    // unless null column
	    if (n) {
	      iBlock = numberBlocks;
	      numberBlocks++;
	      numberGoodBlocks++;
	      int jRow = row[start];
	      rowBlock[jRow]=iBlock;
	      whichColumn[iBlock]=jRow;
	      numberMarkedRows += n;
	      columnBlock[iBlock] = n;
	      for (CoinBigIndex j=start+1;j<end;j++) {
		int iRow = row[j];
		rowBlock[iRow]=iBlock;
		whichRow[jRow]=iRow;
		jRow=iRow;
	      }
	    } else {
	      columnBlock[numberColumns]++;
	    }
	  } else {
	    // add all to this block if not already in
	    int jRow = whichColumn[iBlock];
	    for (CoinBigIndex j=start;j<end;j++) {
	      int iRow = row[j];
	      if (rowBlock[iRow]<0) {
		numberMarkedRows++;
		columnBlock[iBlock]++;
		rowBlock[iRow]=iBlock;
		whichRow[iRow]=jRow;
		jRow=iRow;
	      }
	    }
	    whichColumn[iBlock]=jRow;
	  }
	  columnsDone++;
	  if (iBlock>=0)
	    maximumInBlock = CoinMax(maximumInBlock,columnBlock[iBlock]);
	  if (columnsDone>=checkAfterColumns) {
	    assert (numberGoodBlocks>0);
	    double averageSize = static_cast<double>(numberMarkedRows)/
	      static_cast<double>(numberGoodBlocks);
#ifndef OSL_WAY
	    double notionalBlocks = static_cast<double>(numberMarkedRows)/
	      averageSize;
	    if (maximumInBlock<3*averageSize&&numberGoodBlocks>2) {
	      if(best*(numberColumns-columnsDone) < notionalBlocks) {
		best = notionalBlocks/
		  static_cast<double> (numberColumns-columnsDone);
		bestColumn = kColumn;
	      }
	    }
#else
	    if (maximumInBlock*10<numberRows*11&&numberGoodBlocks>1) {
	      double test = maximumInBlock + 0.0*averageSize;
	      if(best*static_cast<double>(columnsDone) > test) {
		best = test/static_cast<double> (columnsDone);
		bestColumn = kColumn;
		bestColumnsDone=columnsDone;
	      }
	    }
#endif
	  }
	}
#ifndef OSL_WAY
	best2[iWay]=best;
#else
	if (bestColumnsDone<numberColumns)
	  best2[iWay]=-(numberColumns-bestColumnsDone);
	else
	  best2[iWay]=-numberColumns;
#endif
	column2[iWay]=bestColumn;
      }
      // mark columns
      int nMaster;
      CoinFillN(columnBlock,numberColumns,-2);
      if (best2[2]<best2[0]||best2[2]<best2[1]) {
	int iColumn1;
	int iColumn2;
	if (best2[0]>best2[1]) {
	  // End columns in master
	  iColumn1=column2[0]+1;
	  iColumn2=numberColumns;
	} else {
	  // Beginning columns in master
	  iColumn1=0;
	  iColumn2=numberColumns-column2[1];
	}
	nMaster = iColumn2-iColumn1;
	CoinFillN(columnBlock+iColumn1,nMaster,-1);
      } else {
	// sorted
	// End columns in master (in order)
	int iColumn1=column2[2]+1;
	nMaster = numberColumns-iColumn1;
	for (int i=iColumn1;i<numberColumns;i++)
	  columnBlock[stack[i]]=-1;
      }
      if (nMaster*2>numberColumns) {
	goodBenders=false;
	printf("%d columns out of %d would be in master - no good\n",
	       nMaster,numberColumns);
	delete [] rowBlock;
	delete [] columnBlock;
	delete [] whichRow;
	delete [] whichColumn;
	delete [] stack;
	CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper,
			columnLower,columnUpper,objective);
	model.setObjectiveOffset(objectiveOffset);
	addBlock("row_master","column_master",model);
	return 0;
      }
    } else {
      for (iColumn=0;iColumn<numberColumns;iColumn++)
	columnBlock[iColumn]=-2;
      // these are master columns
      if (numberColumns==2426) {
	// ken-7 dual
	for (iColumn=2401;iColumn<numberColumns;iColumn++)
	  columnBlock[iColumn]=-1;
      }
    }
    numberBlocks=0;
    CoinFillN(rowBlock,numberRows,-2);
    for (iRow=0;iRow<numberRows;iRow++) {
      int kstart = rowStart[iRow];
      int kend = rowStart[iRow]+rowLength[iRow];
      if (rowBlock[iRow]==-2) {
	// row not allocated
	int j;
	int nstack=0;
	for (j=kstart;j<kend;j++) {
	  int iColumn= column[j];
	  if (columnBlock[iColumn]!=-1) {
	    assert(columnBlock[iColumn]==-2);
	    columnBlock[iColumn]=numberBlocks; // mark
	    stack[nstack++] = iColumn;
	  }
	}
	if (nstack) {
	  // new block - put all connected in
	  numberBlocks++;
	  rowBlock[iRow]=numberBlocks-1;
	  while (nstack) {
	    int iColumn = stack[--nstack];
	    int k;
	    for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];k++) {
	      int iRow = row[k];
	      int kkstart = rowStart[iRow];
	      int kkend = kkstart + rowLength[iRow];
	      if (rowBlock[iRow]==-2) {
		rowBlock[iRow]=numberBlocks-1; // mark
		// row not allocated
		int jj;
		for (jj=kkstart;jj<kkend;jj++) {
		  int jColumn= column[jj];
		  if (columnBlock[jColumn]==-2) {
		    columnBlock[jColumn]=numberBlocks-1;
		    stack[nstack++]=jColumn;
		  }
		}
	      } else {
		assert (rowBlock[iRow]==numberBlocks-1);
	      }
	    }
	  }
	} else {
	  // Only in master
	  rowBlock[iRow]=-1;
	}
      }
    }
    delete [] stack;
    int numberMasterColumns=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock==-1)
	numberMasterColumns++;
    }
    int numberMasterRows=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock==-1)
	numberMasterRows++;
    }
    if (numberBlocks<=maxBlocks)
      printf("%d blocks found - %d columns, %d rows in master\n",
	     numberBlocks,numberMasterColumns,numberMasterRows);
    else
      printf("%d blocks found (reduced to %d) - %d columns, %d rows in master\n",
	     numberBlocks,maxBlocks,numberMasterColumns,numberMasterRows);
    if (numberBlocks) {
      if (numberBlocks>maxBlocks) {
	int iBlock;
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  iBlock = columnBlock[iColumn];
	  if (iBlock>=0)
	    columnBlock[iColumn] = iBlock%maxBlocks;
	}
	for (iRow=0;iRow<numberRows;iRow++) {
	  iBlock = rowBlock[iRow];
	  if (iBlock>=0)
	    rowBlock[iRow] = iBlock%maxBlocks;
	}
	numberBlocks=maxBlocks;
      }
    }
    // make up problems
    // Create all sub problems
    // Space for creating
    double * obj = new double [numberColumns];
    double * rowLo = new double [numberRows];
    double * rowUp = new double [numberRows];
    double * columnLo = new double [numberColumns];
    double * columnUp = new double [numberColumns];
    // Counts
    int * columnCount = reinterpret_cast<int *>(columnLo);
    CoinZeroN(columnCount,numberBlocks);
    for (int i=0;i<numberColumns;i++) {
      int iBlock=columnBlock[i];
      if (iBlock>=0)
	columnCount[iBlock]++;
    }
    // allocate empty columns
    for (int i=0;i<numberColumns;i++) {
      int iBlock=columnBlock[i];
      if (iBlock==-2) {
	// find block with smallest count
	int iSmall=-1;
	int smallest = numberColumns;
	for (int j=0;j<numberBlocks;j++) {
	  if (columnCount[j]<smallest) {
	    iSmall=j;
	    smallest=columnCount[j];
	  }
	}
	columnBlock[i]=iSmall;
	columnCount[iSmall]++;
      }
    }
    int * rowCount = reinterpret_cast<int *>(columnUp);
    CoinZeroN(rowCount,numberBlocks);
    for (int i=0;i<numberRows;i++) {
      int iBlock=rowBlock[i];
      if (iBlock>=0)
	rowCount[iBlock]++;
    }
    int maximumSize=0;
    for (int i=0;i<numberBlocks;i++) {
      printf("Block %d has %d columns and %d rows\n",i,
	     columnCount[i],rowCount[i]);
      int k=2*columnCount[i]+rowCount[i];
      maximumSize = CoinMax(maximumSize,k);
    }
    if (maximumSize*10>4*(2*numberColumns+numberRows)) {
      // No good
      printf("Doesn't look good\n");
      delete [] rowBlock;
      delete [] columnBlock;
      delete [] whichRow;
      delete [] whichColumn;
      delete [] obj ;
      delete [] columnLo ;
      delete [] columnUp ;
      delete [] rowLo ;
      delete [] rowUp ;
      CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper,
		      columnLower,columnUpper,objective);
      model.setObjectiveOffset(objectiveOffset);
      addBlock("row_master","column_master",model);
      return 0;
    }
    // Name for master so at beginning
    addColumnBlock(numberMasterColumns,"column_master");
    // Arrays
    // get full matrix
    CoinPackedMatrix fullMatrix = matrix;
    int numberRow2,numberColumn2;
    int iBlock;
    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
      char rowName[20];
      sprintf(rowName,"row_%d",iBlock);
      char columnName[20];
      sprintf(columnName,"column_%d",iBlock);
      numberRow2=0;
      numberColumn2=0;
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	if (iBlock==columnBlock[iColumn]) {
	  obj[numberColumn2]=objective[iColumn];
	  columnLo[numberColumn2]=columnLower[iColumn];
	  columnUp[numberColumn2]=columnUpper[iColumn];
	  whichColumn[numberColumn2++]=iColumn;
	}
      }
      for (iRow=0;iRow<numberRows;iRow++) {
	if (iBlock==rowBlock[iRow]) {
	  rowLo[numberRow2]=rowLower[iRow];
	  rowUp[numberRow2]=rowUpper[iRow];
	  whichRow[numberRow2++]=iRow;
	}
      }
      // Diagonal block
      CoinPackedMatrix mat(fullMatrix,
			   numberRow2,whichRow,
			   numberColumn2,whichColumn);
      // make sure correct dimensions
      mat.setDimensions(numberRow2,numberColumn2);
      CoinModel * block = new CoinModel(numberRow2,numberColumn2,&mat,
					rowLo,rowUp,columnLo,columnUp,obj);
      block->setOriginalIndices(whichRow,whichColumn);
      addBlock(rowName,columnName,block); // takes ownership
      // and beginning block
      numberColumn2=0;
      // get beginning matrix
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	int iBlock = columnBlock[iColumn];
	if (iBlock==-1) {
	  whichColumn[numberColumn2++]=iColumn;
	}
      }
      CoinPackedMatrix beginning(fullMatrix,
			   numberRow2,whichRow,
			   numberColumn2,whichColumn);
      // make sure correct dimensions *********
      beginning.setDimensions(numberRow2,numberColumn2);
      block = new CoinModel(numberRow2,numberMasterColumns,&beginning,
			    NULL,NULL,NULL,NULL,NULL);
      block->setOriginalIndices(whichRow,whichColumn);
      addBlock(rowName,"column_master",block); // takes ownership
    }
    // and master
    numberRow2=0;
    numberColumn2=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock==-1) {
	obj[numberColumn2]=objective[iColumn];
	columnLo[numberColumn2]=columnLower[iColumn];
	columnUp[numberColumn2]=columnUpper[iColumn];
	whichColumn[numberColumn2++]=iColumn;
      }
    }
    for (iRow=0;iRow<numberRows;iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock<0) {
	rowLo[numberRow2]=rowLower[iRow];
	rowUp[numberRow2]=rowUpper[iRow];
	whichRow[numberRow2++]=iRow;
      }
    }
    delete [] rowBlock;
    delete [] columnBlock;
    CoinPackedMatrix beginning(fullMatrix,
			 numberRow2,whichRow,
			 numberColumn2,whichColumn);
    // make sure correct dimensions
    beginning.setDimensions(numberRow2,numberColumn2);
    CoinModel * block = new CoinModel(numberRow2,numberColumn2,&beginning,
				      rowLo,rowUp,
				      columnLo,columnUp,obj);
    block->setOriginalIndices(whichRow,whichColumn);
    addBlock("row_master","column_master",block); // takes ownership
    delete [] whichRow;
    delete [] whichColumn;
    delete [] obj ;
    delete [] columnLo ;
    delete [] columnUp ;
    delete [] rowLo ;
    delete [] rowUp ;
  } else {
    abort();
  }
  return numberBlocks;
}
/* Decompose a CoinModel
   1 - try D-W
   2 - try Benders
   3 - try Staircase
   Returns number of blocks or zero if no structure
*/
int
CoinStructuredModel::decompose(const CoinModel & coinModel, int type,
			       int maxBlocks)
{
  const CoinPackedMatrix * matrix = coinModel.packedMatrix();
  assert (matrix!=NULL);
  // Arrays
  const double * objective = coinModel.objectiveArray();
  const double * columnLower = coinModel.columnLowerArray();
  const double * columnUpper = coinModel.columnUpperArray();
  const double * rowLower = coinModel.rowLowerArray();
  const double * rowUpper = coinModel.rowUpperArray();
  return decompose(*matrix,
		   rowLower,  rowUpper,
		   columnLower,  columnUpper,
		   objective, type,maxBlocks,
		   coinModel.objectiveOffset());
}
// Return block corresponding to row and column
const CoinBaseModel *
CoinStructuredModel::block(int row,int column) const
{
  const CoinBaseModel * block = NULL;
  if (blockType_) {
    for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) {
      if (blockType_[iBlock].rowBlock==row&&
	  blockType_[iBlock].columnBlock==column) {
	block = blocks_[iBlock];
	break;
      }
    }
  }
  return block;
}
// Return block corresponding to row and column as CoinModel
const CoinBaseModel *
CoinStructuredModel::coinBlock(int row,int column) const
{
  const CoinModel * block = NULL;
  if (blockType_) {
    for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) {
      if (blockType_[iBlock].rowBlock==row&&
	  blockType_[iBlock].columnBlock==column) {
	block = dynamic_cast<CoinModel *>(blocks_[iBlock]);
	assert (block);
	break;
      }
    }
  }
  return block;
}
/* Fill pointers corresponding to row and column.
   False if any missing */
CoinModelBlockInfo
CoinStructuredModel::block(int row,int column,
			   const double * & rowLower, const double * & rowUpper,
			   const double * & columnLower, const double * & columnUpper,
			   const double * & objective) const
{
  CoinModelBlockInfo info;
  //memset(&info,0,sizeof(info));
  rowLower=NULL;
  rowUpper=NULL;
  columnLower=NULL;
  columnUpper=NULL;
  objective=NULL;
  if (blockType_) {
    for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) {
      CoinModel * thisBlock = coinBlock(iBlock);
      if (blockType_[iBlock].rowBlock==row) {
	if (blockType_[iBlock].rhs) {
	  info.rhs=1;
	  rowLower = thisBlock->rowLowerArray();
	  rowUpper = thisBlock->rowUpperArray();
	}
      }
      if (blockType_[iBlock].columnBlock==column) {
	if (blockType_[iBlock].bounds) {
	  info.bounds=1;
	  columnLower = thisBlock->columnLowerArray();
	  columnUpper = thisBlock->columnUpperArray();
	  objective = thisBlock->objectiveArray();
	}
      }
    }
  }
  return info;
}
// Return block number corresponding to row and column
int
CoinStructuredModel::blockIndex(int row,int column) const
{
  int block=-1;
  if (blockType_) {
    for (int iBlock=0;iBlock<numberElementBlocks_;iBlock++) {
      if (blockType_[iBlock].rowBlock==row&&
	  blockType_[iBlock].columnBlock==column) {
	block = iBlock;
	break;
      }
    }
  }
  return block;
}