Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpGubDynamicMatrix.cpp

    r1732 r2385  
    33// Corporation and others.  All Rights Reserved.
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    5 
    65
    76#include <cstdio>
     
    2726// Default Constructor
    2827//-------------------------------------------------------------------
    29 ClpGubDynamicMatrix::ClpGubDynamicMatrix ()
    30      : ClpGubMatrix(),
    31        objectiveOffset_(0.0),
    32        startColumn_(NULL),
    33        row_(NULL),
    34        element_(NULL),
    35        cost_(NULL),
    36        fullStart_(NULL),
    37        id_(NULL),
    38        dynamicStatus_(NULL),
    39        lowerColumn_(NULL),
    40        upperColumn_(NULL),
    41        lowerSet_(NULL),
    42        upperSet_(NULL),
    43        numberGubColumns_(0),
    44        firstAvailable_(0),
    45        savedFirstAvailable_(0),
    46        firstDynamic_(0),
    47        lastDynamic_(0),
    48       numberElements_(0)
     28ClpGubDynamicMatrix::ClpGubDynamicMatrix()
     29  : ClpGubMatrix()
     30  , objectiveOffset_(0.0)
     31  , startColumn_(NULL)
     32  , row_(NULL)
     33  , element_(NULL)
     34  , cost_(NULL)
     35  , fullStart_(NULL)
     36  , id_(NULL)
     37  , dynamicStatus_(NULL)
     38  , lowerColumn_(NULL)
     39  , upperColumn_(NULL)
     40  , lowerSet_(NULL)
     41  , upperSet_(NULL)
     42  , numberGubColumns_(0)
     43  , firstAvailable_(0)
     44  , savedFirstAvailable_(0)
     45  , firstDynamic_(0)
     46  , lastDynamic_(0)
     47  , numberElements_(0)
    4948{
    50      setType(13);
     49  setType(13);
    5150}
    5251
     
    5453// Copy constructor
    5554//-------------------------------------------------------------------
    56 ClpGubDynamicMatrix::ClpGubDynamicMatrix (const ClpGubDynamicMatrix & rhs)
    57      : ClpGubMatrix(rhs)
     55ClpGubDynamicMatrix::ClpGubDynamicMatrix(const ClpGubDynamicMatrix &rhs)
     56  : ClpGubMatrix(rhs)
    5857{
    59      objectiveOffset_ = rhs.objectiveOffset_;
    60      numberGubColumns_ = rhs.numberGubColumns_;
    61      firstAvailable_ = rhs.firstAvailable_;
    62      savedFirstAvailable_ = rhs.savedFirstAvailable_;
    63      firstDynamic_ = rhs.firstDynamic_;
    64      lastDynamic_ = rhs.lastDynamic_;
    65      numberElements_ = rhs.numberElements_;
    66      startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
    67      CoinBigIndex numberElements = startColumn_[numberGubColumns_];
    68      row_ = ClpCopyOfArray(rhs.row_, numberElements);;
    69      element_ = ClpCopyOfArray(rhs.element_, numberElements);;
    70      cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
    71      fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
    72      id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
    73      lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
    74      upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
    75      dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
    76      lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
    77      upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
     58  objectiveOffset_ = rhs.objectiveOffset_;
     59  numberGubColumns_ = rhs.numberGubColumns_;
     60  firstAvailable_ = rhs.firstAvailable_;
     61  savedFirstAvailable_ = rhs.savedFirstAvailable_;
     62  firstDynamic_ = rhs.firstDynamic_;
     63  lastDynamic_ = rhs.lastDynamic_;
     64  numberElements_ = rhs.numberElements_;
     65  startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
     66  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
     67  row_ = ClpCopyOfArray(rhs.row_, numberElements);
     68  ;
     69  element_ = ClpCopyOfArray(rhs.element_, numberElements);
     70  ;
     71  cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
     72  fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
     73  id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
     74  lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
     75  upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
     76  dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
     77  lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
     78  upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
    7879}
    7980
    8081/* This is the real constructor*/
    81 ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex * model, int numberSets,
    82           int numberGubColumns, const int * starts,
    83           const double * lower, const double * upper,
    84           const CoinBigIndex * startColumn, const int * row,
    85           const double * element, const double * cost,
    86           const double * lowerColumn, const double * upperColumn,
    87           const unsigned char * status)
    88      : ClpGubMatrix()
     82ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex *model, int numberSets,
     83  int numberGubColumns, const int *starts,
     84  const double *lower, const double *upper,
     85  const CoinBigIndex *startColumn, const int *row,
     86  const double *element, const double *cost,
     87  const double *lowerColumn, const double *upperColumn,
     88  const unsigned char *status)
     89  : ClpGubMatrix()
    8990{
    90      objectiveOffset_ = model->objectiveOffset();
    91      model_ = model;
    92      numberSets_ = numberSets;
    93      numberGubColumns_ = numberGubColumns;
    94      fullStart_ = ClpCopyOfArray(starts, numberSets_ + 1);
    95      lower_ = ClpCopyOfArray(lower, numberSets_);
    96      upper_ = ClpCopyOfArray(upper, numberSets_);
    97      int numberColumns = model->numberColumns();
    98      int numberRows = model->numberRows();
    99      // Number of columns needed
    100      int numberGubInSmall = numberSets_ + numberRows + 2 * model->factorizationFrequency() + 2;
    101      // for small problems this could be too big
    102      //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
    103      int numberNeeded = numberGubInSmall + numberColumns;
    104      firstAvailable_ = numberColumns;
    105      savedFirstAvailable_ = numberColumns;
    106      firstDynamic_ = numberColumns;
    107      lastDynamic_ = numberNeeded;
    108      startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
    109      CoinBigIndex numberElements = startColumn_[numberGubColumns_];
    110      row_ = ClpCopyOfArray(row, numberElements);
    111      element_ = new double[numberElements];
    112      CoinBigIndex i;
    113      for (i = 0; i < numberElements; i++)
    114           element_[i] = element[i];
    115      cost_ = new double[numberGubColumns_];
    116      for (i = 0; i < numberGubColumns_; i++) {
    117           cost_[i] = cost[i];
    118           // need sorted
    119           CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]);
    120      }
    121      if (lowerColumn) {
    122           lowerColumn_ = new double[numberGubColumns_];
    123           for (i = 0; i < numberGubColumns_; i++)
    124                lowerColumn_[i] = lowerColumn[i];
    125      } else {
    126           lowerColumn_ = NULL;
    127      }
    128      if (upperColumn) {
    129           upperColumn_ = new double[numberGubColumns_];
    130           for (i = 0; i < numberGubColumns_; i++)
    131                upperColumn_[i] = upperColumn[i];
    132      } else {
    133           upperColumn_ = NULL;
    134      }
    135      if (upperColumn || lowerColumn) {
    136           lowerSet_ = new double[numberSets_];
    137           for (i = 0; i < numberSets_; i++) {
    138                if (lower[i] > -1.0e20)
    139                     lowerSet_[i] = lower[i];
    140                else
    141                     lowerSet_[i] = -1.0e30;
    142           }
    143           upperSet_ = new double[numberSets_];
    144           for (i = 0; i < numberSets_; i++) {
    145                if (upper[i] < 1.0e20)
    146                     upperSet_[i] = upper[i];
    147                else
    148                     upperSet_[i] = 1.0e30;
    149           }
    150      } else {
    151           lowerSet_ = NULL;
    152           upperSet_ = NULL;
    153      }
    154      start_ = NULL;
    155      end_ = NULL;
    156      dynamicStatus_ = NULL;
    157      id_ = new int[numberGubInSmall];
    158      for (i = 0; i < numberGubInSmall; i++)
    159           id_[i] = -1;
    160      ClpPackedMatrix* originalMatrixA =
    161           dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
    162      assert (originalMatrixA);
    163      CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
    164      originalMatrixA->setMatrixNull(); // so can be deleted safely
    165      // guess how much space needed
    166      double guess = originalMatrix->getNumElements() + 10;
    167      guess /= static_cast<double> (numberColumns);
    168      guess *= 2 * numberGubColumns_;
    169      numberElements_ = static_cast<int> (CoinMin(guess, 10000000.0));
    170      numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
    171      matrix_ = originalMatrix;
    172      flags_ &= ~1;
    173      // resize model (matrix stays same)
    174      model->resize(numberRows, numberNeeded);
    175      if (upperColumn_) {
    176           // set all upper bounds so we have enough space
    177           double * columnUpper = model->columnUpper();
    178           for(i = firstDynamic_; i < lastDynamic_; i++)
    179                columnUpper[i] = 1.0e10;
    180      }
    181      // resize matrix
    182      // extra 1 is so can keep number of elements handy
    183      originalMatrix->reserve(numberNeeded, numberElements_, true);
    184      originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
    185      originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
    186      // redo number of columns
    187      numberColumns = matrix_->getNumCols();
    188      backward_ = new int[numberNeeded];
    189      backToPivotRow_ = new int[numberNeeded];
    190      // We know a bit better
    191      delete [] changeCost_;
    192      changeCost_ = new double [numberRows+numberSets_];
    193      keyVariable_ = new int[numberSets_];
    194      // signal to need new ordering
    195      next_ = NULL;
    196      for (int iColumn = 0; iColumn < numberNeeded; iColumn++)
    197           backward_[iColumn] = -1;
     91  objectiveOffset_ = model->objectiveOffset();
     92  model_ = model;
     93  numberSets_ = numberSets;
     94  numberGubColumns_ = numberGubColumns;
     95  fullStart_ = ClpCopyOfArray(starts, numberSets_ + 1);
     96  lower_ = ClpCopyOfArray(lower, numberSets_);
     97  upper_ = ClpCopyOfArray(upper, numberSets_);
     98  int numberColumns = model->numberColumns();
     99  int numberRows = model->numberRows();
     100  // Number of columns needed
     101  int numberGubInSmall = numberSets_ + numberRows + 2 * model->factorizationFrequency() + 2;
     102  // for small problems this could be too big
     103  //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
     104  int numberNeeded = numberGubInSmall + numberColumns;
     105  firstAvailable_ = numberColumns;
     106  savedFirstAvailable_ = numberColumns;
     107  firstDynamic_ = numberColumns;
     108  lastDynamic_ = numberNeeded;
     109  startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
     110  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
     111  row_ = ClpCopyOfArray(row, numberElements);
     112  element_ = new double[numberElements];
     113  CoinBigIndex i;
     114  for (i = 0; i < numberElements; i++)
     115    element_[i] = element[i];
     116  cost_ = new double[numberGubColumns_];
     117  for (i = 0; i < numberGubColumns_; i++) {
     118    cost_[i] = cost[i];
     119    // need sorted
     120    CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i + 1], element_ + startColumn_[i]);
     121  }
     122  if (lowerColumn) {
     123    lowerColumn_ = new double[numberGubColumns_];
     124    for (i = 0; i < numberGubColumns_; i++)
     125      lowerColumn_[i] = lowerColumn[i];
     126  } else {
     127    lowerColumn_ = NULL;
     128  }
     129  if (upperColumn) {
     130    upperColumn_ = new double[numberGubColumns_];
     131    for (i = 0; i < numberGubColumns_; i++)
     132      upperColumn_[i] = upperColumn[i];
     133  } else {
     134    upperColumn_ = NULL;
     135  }
     136  if (upperColumn || lowerColumn) {
     137    lowerSet_ = new double[numberSets_];
     138    for (i = 0; i < numberSets_; i++) {
     139      if (lower[i] > -1.0e20)
     140        lowerSet_[i] = lower[i];
     141      else
     142        lowerSet_[i] = -1.0e30;
     143    }
     144    upperSet_ = new double[numberSets_];
     145    for (i = 0; i < numberSets_; i++) {
     146      if (upper[i] < 1.0e20)
     147        upperSet_[i] = upper[i];
     148      else
     149        upperSet_[i] = 1.0e30;
     150    }
     151  } else {
     152    lowerSet_ = NULL;
     153    upperSet_ = NULL;
     154  }
     155  start_ = NULL;
     156  end_ = NULL;
     157  dynamicStatus_ = NULL;
     158  id_ = new int[numberGubInSmall];
     159  for (i = 0; i < numberGubInSmall; i++)
     160    id_[i] = -1;
     161  ClpPackedMatrix *originalMatrixA = dynamic_cast< ClpPackedMatrix * >(model->clpMatrix());
     162  assert(originalMatrixA);
     163  CoinPackedMatrix *originalMatrix = originalMatrixA->getPackedMatrix();
     164  originalMatrixA->setMatrixNull(); // so can be deleted safely
     165  // guess how much space needed
     166  double guess = originalMatrix->getNumElements() + 10;
     167  guess /= static_cast< double >(numberColumns);
     168  guess *= 2 * numberGubColumns_;
     169  numberElements_ = static_cast< int >(CoinMin(guess, 10000000.0));
     170  numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
     171  matrix_ = originalMatrix;
     172  flags_ &= ~1;
     173  // resize model (matrix stays same)
     174  model->resize(numberRows, numberNeeded);
     175  if (upperColumn_) {
     176    // set all upper bounds so we have enough space
     177    double *columnUpper = model->columnUpper();
     178    for (i = firstDynamic_; i < lastDynamic_; i++)
     179      columnUpper[i] = 1.0e10;
     180  }
     181  // resize matrix
     182  // extra 1 is so can keep number of elements handy
     183  originalMatrix->reserve(numberNeeded, numberElements_, true);
     184  originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
     185  originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
     186  // redo number of columns
     187  numberColumns = matrix_->getNumCols();
     188  backward_ = new int[numberNeeded];
     189  backToPivotRow_ = new int[numberNeeded];
     190  // We know a bit better
     191  delete[] changeCost_;
     192  changeCost_ = new double[numberRows + numberSets_];
     193  keyVariable_ = new int[numberSets_];
     194  // signal to need new ordering
     195  next_ = NULL;
     196  for (int iColumn = 0; iColumn < numberNeeded; iColumn++)
     197    backward_[iColumn] = -1;
    198198
    199      firstGub_ = firstDynamic_;
    200      lastGub_ = lastDynamic_;
    201      if (!lowerColumn_ && !upperColumn_)
    202           gubType_ = 8;
    203      if (status) {
    204           status_ = ClpCopyOfArray(status, numberSets_);
    205      } else {
    206           status_ = new unsigned char [numberSets_];
    207           memset(status_, 0, numberSets_);
    208           int i;
    209           for (i = 0; i < numberSets_; i++) {
    210                // make slack key
    211                setStatus(i, ClpSimplex::basic);
    212           }
    213      }
    214      saveStatus_ = new unsigned char [numberSets_];
    215      memset(saveStatus_, 0, numberSets_);
    216      savedKeyVariable_ = new int [numberSets_];
    217      memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
     199  firstGub_ = firstDynamic_;
     200  lastGub_ = lastDynamic_;
     201  if (!lowerColumn_ && !upperColumn_)
     202    gubType_ = 8;
     203  if (status) {
     204    status_ = ClpCopyOfArray(status, numberSets_);
     205  } else {
     206    status_ = new unsigned char[numberSets_];
     207    memset(status_, 0, numberSets_);
     208    int i;
     209    for (i = 0; i < numberSets_; i++) {
     210      // make slack key
     211      setStatus(i, ClpSimplex::basic);
     212    }
     213  }
     214  saveStatus_ = new unsigned char[numberSets_];
     215  memset(saveStatus_, 0, numberSets_);
     216  savedKeyVariable_ = new int[numberSets_];
     217  memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
    218218}
    219219
     
    221221// Destructor
    222222//-------------------------------------------------------------------
    223 ClpGubDynamicMatrix::~ClpGubDynamicMatrix ()
     223ClpGubDynamicMatrix::~ClpGubDynamicMatrix()
    224224{
    225      delete [] startColumn_;
    226      delete [] row_;
    227      delete [] element_;
    228      delete [] cost_;
    229      delete [] fullStart_;
    230      delete [] id_;
    231      delete [] dynamicStatus_;
    232      delete [] lowerColumn_;
    233      delete [] upperColumn_;
    234      delete [] lowerSet_;
    235      delete [] upperSet_;
     225  delete[] startColumn_;
     226  delete[] row_;
     227  delete[] element_;
     228  delete[] cost_;
     229  delete[] fullStart_;
     230  delete[] id_;
     231  delete[] dynamicStatus_;
     232  delete[] lowerColumn_;
     233  delete[] upperColumn_;
     234  delete[] lowerSet_;
     235  delete[] upperSet_;
    236236}
    237237
     
    240240//-------------------------------------------------------------------
    241241ClpGubDynamicMatrix &
    242 ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix& rhs)
     242ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix &rhs)
    243243{
    244      if (this != &rhs) {
    245           ClpGubMatrix::operator=(rhs);
    246           delete [] startColumn_;
    247           delete [] row_;
    248           delete [] element_;
    249           delete [] cost_;
    250           delete [] fullStart_;
    251           delete [] id_;
    252           delete [] dynamicStatus_;
    253           delete [] lowerColumn_;
    254           delete [] upperColumn_;
    255           delete [] lowerSet_;
    256           delete [] upperSet_;
    257           objectiveOffset_ = rhs.objectiveOffset_;
    258           numberGubColumns_ = rhs.numberGubColumns_;
    259           firstAvailable_ = rhs.firstAvailable_;
    260           savedFirstAvailable_ = rhs.savedFirstAvailable_;
    261           firstDynamic_ = rhs.firstDynamic_;
    262           lastDynamic_ = rhs.lastDynamic_;
    263           numberElements_ = rhs.numberElements_;
    264           startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
    265           int numberElements = startColumn_[numberGubColumns_];
    266           row_ = ClpCopyOfArray(rhs.row_, numberElements);;
    267           element_ = ClpCopyOfArray(rhs.element_, numberElements);;
    268           cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
    269           fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
    270           id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
    271           lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
    272           upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
    273           dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
    274           lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
    275           upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
    276      }
    277      return *this;
     244  if (this != &rhs) {
     245    ClpGubMatrix::operator=(rhs);
     246    delete[] startColumn_;
     247    delete[] row_;
     248    delete[] element_;
     249    delete[] cost_;
     250    delete[] fullStart_;
     251    delete[] id_;
     252    delete[] dynamicStatus_;
     253    delete[] lowerColumn_;
     254    delete[] upperColumn_;
     255    delete[] lowerSet_;
     256    delete[] upperSet_;
     257    objectiveOffset_ = rhs.objectiveOffset_;
     258    numberGubColumns_ = rhs.numberGubColumns_;
     259    firstAvailable_ = rhs.firstAvailable_;
     260    savedFirstAvailable_ = rhs.savedFirstAvailable_;
     261    firstDynamic_ = rhs.firstDynamic_;
     262    lastDynamic_ = rhs.lastDynamic_;
     263    numberElements_ = rhs.numberElements_;
     264    startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
     265    int numberElements = startColumn_[numberGubColumns_];
     266    row_ = ClpCopyOfArray(rhs.row_, numberElements);
     267    ;
     268    element_ = ClpCopyOfArray(rhs.element_, numberElements);
     269    ;
     270    cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
     271    fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
     272    id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
     273    lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
     274    upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
     275    dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
     276    lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
     277    upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
     278  }
     279  return *this;
    278280}
    279281//-------------------------------------------------------------------
    280282// Clone
    281283//-------------------------------------------------------------------
    282 ClpMatrixBase * ClpGubDynamicMatrix::clone() const
     284ClpMatrixBase *ClpGubDynamicMatrix::clone() const
    283285{
    284      return new ClpGubDynamicMatrix(*this);
     286  return new ClpGubDynamicMatrix(*this);
    285287}
    286288// Partial pricing
    287 void
    288 ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    289                                     int & bestSequence, int & numberWanted)
     289void ClpGubDynamicMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
     290  int &bestSequence, int &numberWanted)
    290291{
    291      assert(!model->rowScale());
    292      numberWanted = currentWanted_;
    293      if (!numberSets_) {
    294           // no gub
    295           ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
    296           return;
    297      } else {
    298           // and do some proportion of full set
    299           int startG2 = static_cast<int> (startFraction * numberSets_);
    300           int endG2 = static_cast<int> (endFraction * numberSets_ + 0.1);
    301           endG2 = CoinMin(endG2, numberSets_);
    302           //printf("gub price - set start %d end %d\n",
    303           //   startG2,endG2);
    304           double tolerance = model->currentDualTolerance();
    305           double * reducedCost = model->djRegion();
    306           const double * duals = model->dualRowSolution();
    307           double * cost = model->costRegion();
    308           double bestDj;
    309           int numberRows = model->numberRows();
    310           int numberColumns = lastDynamic_;
    311           // If nothing found yet can go all the way to end
    312           int endAll = endG2;
    313           if (bestSequence < 0 && !startG2)
    314                endAll = numberSets_;
    315           if (bestSequence >= 0)
    316                bestDj = fabs(reducedCost[bestSequence]);
    317           else
    318                bestDj = tolerance;
    319           int saveSequence = bestSequence;
    320           double djMod = 0.0;
    321           double infeasibilityCost = model->infeasibilityCost();
    322           double bestDjMod = 0.0;
    323           //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
    324           //     startG2,endG2,numberWanted);
    325           int bestType = -1;
    326           int bestSet = -1;
    327           const double * element = matrix_->getElements();
    328           const int * row = matrix_->getIndices();
    329           const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    330           int * length = matrix_->getMutableVectorLengths();
     292  assert(!model->rowScale());
     293  numberWanted = currentWanted_;
     294  if (!numberSets_) {
     295    // no gub
     296    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
     297    return;
     298  } else {
     299    // and do some proportion of full set
     300    int startG2 = static_cast< int >(startFraction * numberSets_);
     301    int endG2 = static_cast< int >(endFraction * numberSets_ + 0.1);
     302    endG2 = CoinMin(endG2, numberSets_);
     303    //printf("gub price - set start %d end %d\n",
     304    //   startG2,endG2);
     305    double tolerance = model->currentDualTolerance();
     306    double *reducedCost = model->djRegion();
     307    const double *duals = model->dualRowSolution();
     308    double *cost = model->costRegion();
     309    double bestDj;
     310    int numberRows = model->numberRows();
     311    int numberColumns = lastDynamic_;
     312    // If nothing found yet can go all the way to end
     313    int endAll = endG2;
     314    if (bestSequence < 0 && !startG2)
     315      endAll = numberSets_;
     316    if (bestSequence >= 0)
     317      bestDj = fabs(reducedCost[bestSequence]);
     318    else
     319      bestDj = tolerance;
     320    int saveSequence = bestSequence;
     321    double djMod = 0.0;
     322    double infeasibilityCost = model->infeasibilityCost();
     323    double bestDjMod = 0.0;
     324    //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
     325    //     startG2,endG2,numberWanted);
     326    int bestType = -1;
     327    int bestSet = -1;
     328    const double *element = matrix_->getElements();
     329    const int *row = matrix_->getIndices();
     330    const CoinBigIndex *startColumn = matrix_->getVectorStarts();
     331    int *length = matrix_->getMutableVectorLengths();
    331332#if 0
    332333          // make sure first available is clean (in case last iteration rejected)
     
    341342#endif
    342343#ifdef CLP_DEBUG
    343           {
    344                for (int i = firstDynamic_; i < firstAvailable_; i++) {
    345                     assert (getDynamicStatus(id_[i-firstDynamic_]) == inSmall);
    346                }
    347           }
     344    {
     345      for (int i = firstDynamic_; i < firstAvailable_; i++) {
     346        assert(getDynamicStatus(id_[i - firstDynamic_]) == inSmall);
     347      }
     348    }
    348349#endif
    349           int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
    350           int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
    351           for (int iSet = startG2; iSet < endAll; iSet++) {
    352                if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
    353                     // give up
    354                     numberWanted = 0;
    355                     break;
    356                } else if (iSet == endG2 && bestSequence >= 0) {
    357                     break;
    358                }
    359                CoinBigIndex j;
    360                int iBasic = keyVariable_[iSet];
    361                if (iBasic >= numberColumns) {
    362                     djMod = - weight(iSet) * infeasibilityCost;
    363                } else {
    364                     // get dj without
    365                     assert (model->getStatus(iBasic) == ClpSimplex::basic);
    366                     djMod = 0.0;
     350    int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
     351    int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
     352    for (int iSet = startG2; iSet < endAll; iSet++) {
     353      if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
     354        // give up
     355        numberWanted = 0;
     356        break;
     357      } else if (iSet == endG2 && bestSequence >= 0) {
     358        break;
     359      }
     360      CoinBigIndex j;
     361      int iBasic = keyVariable_[iSet];
     362      if (iBasic >= numberColumns) {
     363        djMod = -weight(iSet) * infeasibilityCost;
     364      } else {
     365        // get dj without
     366        assert(model->getStatus(iBasic) == ClpSimplex::basic);
     367        djMod = 0.0;
    367368
    368                     for (j = startColumn[iBasic];
    369                               j < startColumn[iBasic] + length[iBasic]; j++) {
    370                          int jRow = row[j];
    371                          djMod -= duals[jRow] * element[j];
    372                     }
    373                     djMod += cost[iBasic];
    374                     // See if gub slack possible - dj is djMod
    375                     if (getStatus(iSet) == ClpSimplex::atLowerBound) {
    376                          double value = -djMod;
    377                          if (value > tolerance) {
    378                               numberWanted--;
    379                               if (value > bestDj) {
    380                                    // check flagged variable and correct dj
    381                                    if (!flagged(iSet)) {
    382                                         bestDj = value;
    383                                         bestSequence = numberRows + numberColumns + iSet;
    384                                         bestDjMod = djMod;
    385                                         bestType = 0;
    386                                         bestSet = iSet;
    387                                    } else {
    388                                         // just to make sure we don't exit before got something
    389                                         numberWanted++;
    390                                         abort();
    391                                    }
    392                               }
    393                          }
    394                     } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
    395                          double value = djMod;
    396                          if (value > tolerance) {
    397                               numberWanted--;
    398                               if (value > bestDj) {
    399                                    // check flagged variable and correct dj
    400                                    if (!flagged(iSet)) {
    401                                         bestDj = value;
    402                                         bestSequence = numberRows + numberColumns + iSet;
    403                                         bestDjMod = djMod;
    404                                         bestType = 0;
    405                                         bestSet = iSet;
    406                                    } else {
    407                                         // just to make sure we don't exit before got something
    408                                         numberWanted++;
    409                                         abort();
    410                                    }
    411                               }
    412                          }
    413                     }
    414                }
    415                for (int iSequence = fullStart_[iSet]; iSequence < fullStart_[iSet+1]; iSequence++) {
    416                     DynamicStatus status = getDynamicStatus(iSequence);
    417                     if (status != inSmall) {
    418                          double value = cost_[iSequence] - djMod;
    419                          for (j = startColumn_[iSequence];
    420                                    j < startColumn_[iSequence+1]; j++) {
    421                               int jRow = row_[j];
    422                               value -= duals[jRow] * element_[j];
    423                          }
    424                          // change sign if at lower bound
    425                          if (status == atLowerBound)
    426                               value = -value;
    427                          if (value > tolerance) {
    428                               numberWanted--;
    429                               if (value > bestDj) {
    430                                    // check flagged variable and correct dj
    431                                    if (!flagged(iSequence)) {
    432                                         bestDj = value;
    433                                         bestSequence = iSequence;
    434                                         bestDjMod = djMod;
    435                                         bestType = 1;
    436                                         bestSet = iSet;
    437                                    } else {
    438                                         // just to make sure we don't exit before got something
    439                                         numberWanted++;
    440                                    }
    441                               }
    442                          }
    443                     }
    444                }
    445                if (numberWanted <= 0) {
    446                     numberWanted = 0;
    447                     break;
    448                }
    449           }
    450           // Do packed part before gub and small gub - but lightly
    451           int saveMinNeg = minimumGoodReducedCosts_;
    452           int saveSequence2 = bestSequence;
    453           if (bestSequence >= 0)
    454                minimumGoodReducedCosts_ = -2;
    455           int saveLast = lastGub_;
    456           lastGub_ = firstAvailable_;
    457           currentWanted_ = numberWanted;
    458           ClpGubMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
    459           minimumGoodReducedCosts_ = saveMinNeg;
    460           lastGub_ = saveLast;
    461           if (bestSequence != saveSequence2) {
    462                bestType = -1; // in normal or small gub part
    463                saveSequence = bestSequence;
    464           }
    465           if (bestSequence != saveSequence || bestType >= 0) {
    466                double * lowerColumn = model->lowerRegion();
    467                double * upperColumn = model->upperRegion();
    468                double * solution = model->solutionRegion();
    469                if (bestType > 0) {
    470                     // recompute dj and create
    471                     double value = cost_[bestSequence] - bestDjMod;
    472                     for (CoinBigIndex jBigIndex = startColumn_[bestSequence];
    473                               jBigIndex < startColumn_[bestSequence+1]; jBigIndex++) {
    474                          int jRow = row_[jBigIndex];
    475                          value -= duals[jRow] * element_[jBigIndex];
    476                     }
    477                     double * element = matrix_->getMutableElements();
    478                     int * row = matrix_->getMutableIndices();
    479                     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
    480                     int * length = matrix_->getMutableVectorLengths();
    481                     CoinBigIndex numberElements = startColumn[firstAvailable_];
    482                     int numberThis = startColumn_[bestSequence+1] - startColumn_[bestSequence];
    483                     if (numberElements + numberThis > numberElements_) {
    484                          // need to redo
    485                          numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
    486                          matrix_->reserve(numberColumns, numberElements_);
    487                          element = matrix_->getMutableElements();
    488                          row = matrix_->getMutableIndices();
    489                          // these probably okay but be safe
    490                          startColumn = matrix_->getMutableVectorStarts();
    491                          length = matrix_->getMutableVectorLengths();
    492                     }
    493                     // already set startColumn[firstAvailable_]=numberElements;
    494                     length[firstAvailable_] = numberThis;
    495                     model->costRegion()[firstAvailable_] = cost_[bestSequence];
    496                     CoinBigIndex base = startColumn_[bestSequence];
    497                     for (int j = 0; j < numberThis; j++) {
    498                          row[numberElements] = row_[base+j];
    499                          element[numberElements++] = element_[base+j];
    500                     }
    501                     id_[firstAvailable_-firstDynamic_] = bestSequence;
    502                     //printf("best %d\n",bestSequence);
    503                     backward_[firstAvailable_] = bestSet;
    504                     model->solutionRegion()[firstAvailable_] = 0.0;
    505                     if (!lowerColumn_ && !upperColumn_) {
    506                          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
    507                          lowerColumn[firstAvailable_] = 0.0;
    508                          upperColumn[firstAvailable_] = COIN_DBL_MAX;
    509                     } else {
    510                          DynamicStatus status = getDynamicStatus(bestSequence);
    511                          if (lowerColumn_)
    512                               lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
    513                          else
    514                               lowerColumn[firstAvailable_] = 0.0;
    515                          if (upperColumn_)
    516                               upperColumn[firstAvailable_] = upperColumn_[bestSequence];
    517                          else
    518                               upperColumn[firstAvailable_] = COIN_DBL_MAX;
    519                          if (status == atLowerBound) {
    520                               solution[firstAvailable_] = lowerColumn[firstAvailable_];
    521                               model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
    522                          } else {
    523                               solution[firstAvailable_] = upperColumn[firstAvailable_];
    524                               model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
    525                          }
    526                     }
    527                     model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
    528                                                    lowerColumn[firstAvailable_],
    529                                                    upperColumn[firstAvailable_], cost_[bestSequence]);
    530                     bestSequence = firstAvailable_;
    531                     // firstAvailable_ only updated if good pivot (in updatePivot)
    532                     startColumn[firstAvailable_+1] = numberElements;
    533                     //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
    534                     reducedCost[bestSequence] = value;
    535                     gubSlackIn_ = -1;
    536                } else {
    537                     // slack - make last column
    538                     gubSlackIn_ = bestSequence - numberRows - numberColumns;
    539                     bestSequence = numberColumns + 2 * numberRows;
    540                     reducedCost[bestSequence] = bestDjMod;
    541                     //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
    542                     model->setStatus(bestSequence, getStatus(gubSlackIn_));
    543                     if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
    544                          solution[bestSequence] = upper_[gubSlackIn_];
    545                     else
    546                          solution[bestSequence] = lower_[gubSlackIn_];
    547                     lowerColumn[bestSequence] = lower_[gubSlackIn_];
    548                     upperColumn[bestSequence] = upper_[gubSlackIn_];
    549                     model->costRegion()[bestSequence] = 0.0;
    550                     model->nonLinearCost()->setOne(bestSequence, solution[bestSequence], lowerColumn[bestSequence],
    551                                                    upperColumn[bestSequence], 0.0);
    552                }
    553                savedBestSequence_ = bestSequence;
    554                savedBestDj_ = reducedCost[savedBestSequence_];
    555           }
    556           // See if may be finished
    557           if (!startG2 && bestSequence < 0)
    558                infeasibilityWeight_ = model_->infeasibilityCost();
    559           else if (bestSequence >= 0)
    560                infeasibilityWeight_ = -1.0;
    561      }
    562      currentWanted_ = numberWanted;
     369        for (j = startColumn[iBasic];
     370             j < startColumn[iBasic] + length[iBasic]; j++) {
     371          int jRow = row[j];
     372          djMod -= duals[jRow] * element[j];
     373        }
     374        djMod += cost[iBasic];
     375        // See if gub slack possible - dj is djMod
     376        if (getStatus(iSet) == ClpSimplex::atLowerBound) {
     377          double value = -djMod;
     378          if (value > tolerance) {
     379            numberWanted--;
     380            if (value > bestDj) {
     381              // check flagged variable and correct dj
     382              if (!flagged(iSet)) {
     383                bestDj = value;
     384                bestSequence = numberRows + numberColumns + iSet;
     385                bestDjMod = djMod;
     386                bestType = 0;
     387                bestSet = iSet;
     388              } else {
     389                // just to make sure we don't exit before got something
     390                numberWanted++;
     391                abort();
     392              }
     393            }
     394          }
     395        } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
     396          double value = djMod;
     397          if (value > tolerance) {
     398            numberWanted--;
     399            if (value > bestDj) {
     400              // check flagged variable and correct dj
     401              if (!flagged(iSet)) {
     402                bestDj = value;
     403                bestSequence = numberRows + numberColumns + iSet;
     404                bestDjMod = djMod;
     405                bestType = 0;
     406                bestSet = iSet;
     407              } else {
     408                // just to make sure we don't exit before got something
     409                numberWanted++;
     410                abort();
     411              }
     412            }
     413          }
     414        }
     415      }
     416      for (int iSequence = fullStart_[iSet]; iSequence < fullStart_[iSet + 1]; iSequence++) {
     417        DynamicStatus status = getDynamicStatus(iSequence);
     418        if (status != inSmall) {
     419          double value = cost_[iSequence] - djMod;
     420          for (j = startColumn_[iSequence];
     421               j < startColumn_[iSequence + 1]; j++) {
     422            int jRow = row_[j];
     423            value -= duals[jRow] * element_[j];
     424          }
     425          // change sign if at lower bound
     426          if (status == atLowerBound)
     427            value = -value;
     428          if (value > tolerance) {
     429            numberWanted--;
     430            if (value > bestDj) {
     431              // check flagged variable and correct dj
     432              if (!flagged(iSequence)) {
     433                bestDj = value;
     434                bestSequence = iSequence;
     435                bestDjMod = djMod;
     436                bestType = 1;
     437                bestSet = iSet;
     438              } else {
     439                // just to make sure we don't exit before got something
     440                numberWanted++;
     441              }
     442            }
     443          }
     444        }
     445      }
     446      if (numberWanted <= 0) {
     447        numberWanted = 0;
     448        break;
     449      }
     450    }
     451    // Do packed part before gub and small gub - but lightly
     452    int saveMinNeg = minimumGoodReducedCosts_;
     453    int saveSequence2 = bestSequence;
     454    if (bestSequence >= 0)
     455      minimumGoodReducedCosts_ = -2;
     456    int saveLast = lastGub_;
     457    lastGub_ = firstAvailable_;
     458    currentWanted_ = numberWanted;
     459    ClpGubMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
     460    minimumGoodReducedCosts_ = saveMinNeg;
     461    lastGub_ = saveLast;
     462    if (bestSequence != saveSequence2) {
     463      bestType = -1; // in normal or small gub part
     464      saveSequence = bestSequence;
     465    }
     466    if (bestSequence != saveSequence || bestType >= 0) {
     467      double *lowerColumn = model->lowerRegion();
     468      double *upperColumn = model->upperRegion();
     469      double *solution = model->solutionRegion();
     470      if (bestType > 0) {
     471        // recompute dj and create
     472        double value = cost_[bestSequence] - bestDjMod;
     473        for (CoinBigIndex jBigIndex = startColumn_[bestSequence];
     474             jBigIndex < startColumn_[bestSequence + 1]; jBigIndex++) {
     475          int jRow = row_[jBigIndex];
     476          value -= duals[jRow] * element_[jBigIndex];
     477        }
     478        double *element = matrix_->getMutableElements();
     479        int *row = matrix_->getMutableIndices();
     480        CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
     481        int *length = matrix_->getMutableVectorLengths();
     482        CoinBigIndex numberElements = startColumn[firstAvailable_];
     483        int numberThis = startColumn_[bestSequence + 1] - startColumn_[bestSequence];
     484        if (numberElements + numberThis > numberElements_) {
     485          // need to redo
     486          numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
     487          matrix_->reserve(numberColumns, numberElements_);
     488          element = matrix_->getMutableElements();
     489          row = matrix_->getMutableIndices();
     490          // these probably okay but be safe
     491          startColumn = matrix_->getMutableVectorStarts();
     492          length = matrix_->getMutableVectorLengths();
     493        }
     494        // already set startColumn[firstAvailable_]=numberElements;
     495        length[firstAvailable_] = numberThis;
     496        model->costRegion()[firstAvailable_] = cost_[bestSequence];
     497        CoinBigIndex base = startColumn_[bestSequence];
     498        for (int j = 0; j < numberThis; j++) {
     499          row[numberElements] = row_[base + j];
     500          element[numberElements++] = element_[base + j];
     501        }
     502        id_[firstAvailable_ - firstDynamic_] = bestSequence;
     503        //printf("best %d\n",bestSequence);
     504        backward_[firstAvailable_] = bestSet;
     505        model->solutionRegion()[firstAvailable_] = 0.0;
     506        if (!lowerColumn_ && !upperColumn_) {
     507          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
     508          lowerColumn[firstAvailable_] = 0.0;
     509          upperColumn[firstAvailable_] = COIN_DBL_MAX;
     510        } else {
     511          DynamicStatus status = getDynamicStatus(bestSequence);
     512          if (lowerColumn_)
     513            lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
     514          else
     515            lowerColumn[firstAvailable_] = 0.0;
     516          if (upperColumn_)
     517            upperColumn[firstAvailable_] = upperColumn_[bestSequence];
     518          else
     519            upperColumn[firstAvailable_] = COIN_DBL_MAX;
     520          if (status == atLowerBound) {
     521            solution[firstAvailable_] = lowerColumn[firstAvailable_];
     522            model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
     523          } else {
     524            solution[firstAvailable_] = upperColumn[firstAvailable_];
     525            model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
     526          }
     527        }
     528        model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
     529          lowerColumn[firstAvailable_],
     530          upperColumn[firstAvailable_], cost_[bestSequence]);
     531        bestSequence = firstAvailable_;
     532        // firstAvailable_ only updated if good pivot (in updatePivot)
     533        startColumn[firstAvailable_ + 1] = numberElements;
     534        //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
     535        reducedCost[bestSequence] = value;
     536        gubSlackIn_ = -1;
     537      } else {
     538        // slack - make last column
     539        gubSlackIn_ = bestSequence - numberRows - numberColumns;
     540        bestSequence = numberColumns + 2 * numberRows;
     541        reducedCost[bestSequence] = bestDjMod;
     542        //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
     543        model->setStatus(bestSequence, getStatus(gubSlackIn_));
     544        if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
     545          solution[bestSequence] = upper_[gubSlackIn_];
     546        else
     547          solution[bestSequence] = lower_[gubSlackIn_];
     548        lowerColumn[bestSequence] = lower_[gubSlackIn_];
     549        upperColumn[bestSequence] = upper_[gubSlackIn_];
     550        model->costRegion()[bestSequence] = 0.0;
     551        model->nonLinearCost()->setOne(bestSequence, solution[bestSequence], lowerColumn[bestSequence],
     552          upperColumn[bestSequence], 0.0);
     553      }
     554      savedBestSequence_ = bestSequence;
     555      savedBestDj_ = reducedCost[savedBestSequence_];
     556    }
     557    // See if may be finished
     558    if (!startG2 && bestSequence < 0)
     559      infeasibilityWeight_ = model_->infeasibilityCost();
     560    else if (bestSequence >= 0)
     561      infeasibilityWeight_ = -1.0;
     562  }
     563  currentWanted_ = numberWanted;
    563564}
    564565// This is local to Gub to allow synchronization when status is good
    565 int
    566 ClpGubDynamicMatrix::synchronize(ClpSimplex * model, int mode)
     566int ClpGubDynamicMatrix::synchronize(ClpSimplex *model, int mode)
    567567{
    568      int returnNumber = 0;
    569      switch (mode) {
    570      case 0: {
     568  int returnNumber = 0;
     569  switch (mode) {
     570  case 0: {
    571571#ifdef CLP_DEBUG
    572           {
    573                for (int i = 0; i < numberSets_; i++)
    574                     assert(toIndex_[i] == -1);
    575           }
     572    {
     573      for (int i = 0; i < numberSets_; i++)
     574        assert(toIndex_[i] == -1);
     575    }
    576576#endif
    577           // lookup array
    578           int * lookup = new int[lastDynamic_];
    579           int iColumn;
    580           int numberColumns = model->numberColumns();
    581           double * element = matrix_->getMutableElements();
    582           int * row = matrix_->getMutableIndices();
    583           CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
    584           int * length = matrix_->getMutableVectorLengths();
    585           double * cost = model->costRegion();
    586           double * lowerColumn = model->lowerRegion();
    587           double * upperColumn = model->upperRegion();
    588           int * pivotVariable = model->pivotVariable();
    589           CoinBigIndex numberElements = startColumn[firstDynamic_];
    590           // first just do lookup and basic stuff
    591           int currentNumber = firstAvailable_;
    592           firstAvailable_ = firstDynamic_;
    593           int numberToDo = 0;
    594           double objectiveChange = 0.0;
    595           double * solution = model->solutionRegion();
    596           for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
    597                int iSet = backward_[iColumn];
    598                if (toIndex_[iSet] < 0) {
    599                     toIndex_[iSet] = 0;
    600                     fromIndex_[numberToDo++] = iSet;
    601                }
    602                if (model->getStatus(iColumn) == ClpSimplex::basic || iColumn == keyVariable_[iSet]) {
    603                     lookup[iColumn] = firstAvailable_;
    604                     if (iColumn != keyVariable_[iSet]) {
    605                          int iPivot = backToPivotRow_[iColumn];
    606                          backToPivotRow_[firstAvailable_] = iPivot;
    607                          pivotVariable[iPivot] = firstAvailable_;
    608                     }
    609                     firstAvailable_++;
    610                } else {
    611                     int jColumn = id_[iColumn-firstDynamic_];
    612                     setDynamicStatus(jColumn, atLowerBound);
    613                     if (lowerColumn_ || upperColumn_) {
    614                          if (model->getStatus(iColumn) == ClpSimplex::atUpperBound)
    615                               setDynamicStatus(jColumn, atUpperBound);
    616                          // treat solution as if exactly at a bound
    617                          double value = solution[iColumn];
    618                          if (fabs(value - lowerColumn[iColumn]) < fabs(value - upperColumn[iColumn]))
    619                               value = lowerColumn[iColumn];
    620                          else
    621                               value = upperColumn[iColumn];
    622                          objectiveChange += cost[iColumn] * value;
    623                          // redo lower and upper on sets
    624                          double shift = value;
    625                          if (lowerSet_[iSet] > -1.0e20)
    626                               lower_[iSet] = lowerSet_[iSet] - shift;
    627                          if (upperSet_[iSet] < 1.0e20)
    628                               upper_[iSet] = upperSet_[iSet] - shift;
    629                     }
    630                     lookup[iColumn] = -1;
    631                }
    632           }
    633           model->setObjectiveOffset(model->objectiveOffset() + objectiveChange);
    634           firstAvailable_ = firstDynamic_;
    635           for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
    636                if (lookup[iColumn] >= 0) {
    637                     // move
    638                     int jColumn = id_[iColumn-firstDynamic_];
    639                     id_[firstAvailable_-firstDynamic_] = jColumn;
    640                     int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn];
    641                     length[firstAvailable_] = numberThis;
    642                     cost[firstAvailable_] = cost[iColumn];
    643                     lowerColumn[firstAvailable_] = lowerColumn[iColumn];
    644                     upperColumn[firstAvailable_] = upperColumn[iColumn];
    645                     double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
    646                     double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
    647                     if (originalUpper > 1.0e30)
    648                          originalUpper = COIN_DBL_MAX;
    649                     model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn],
    650                                                    originalLower, originalUpper,
    651                                                    cost_[jColumn]);
    652                     CoinBigIndex base = startColumn_[jColumn];
    653                     for (int j = 0; j < numberThis; j++) {
    654                          row[numberElements] = row_[base+j];
    655                          element[numberElements++] = element_[base+j];
    656                     }
    657                     model->setStatus(firstAvailable_, model->getStatus(iColumn));
    658                     backward_[firstAvailable_] = backward_[iColumn];
    659                     solution[firstAvailable_] = solution[iColumn];
    660                     firstAvailable_++;
    661                     startColumn[firstAvailable_] = numberElements;
    662                }
    663           }
    664           // clean up next_
    665           int * temp = new int [firstAvailable_];
    666           for (int jSet = 0; jSet < numberToDo; jSet++) {
    667                int iSet = fromIndex_[jSet];
    668                toIndex_[iSet] = -1;
    669                int last = keyVariable_[iSet];
    670                int j = next_[last];
    671                bool setTemp = true;
    672                if (last < lastDynamic_) {
    673                     last = lookup[last];
    674                     assert (last >= 0);
    675                     keyVariable_[iSet] = last;
    676                } else if (j >= 0) {
    677                     int newJ = lookup[j];
    678                     assert (newJ >= 0);
    679                     j = next_[j];
    680                     next_[last] = newJ;
    681                     last = newJ;
    682                } else {
    683                     next_[last] = -(iSet + numberColumns + 1);
    684                     setTemp = false;
    685                }
    686                while (j >= 0) {
    687                     int newJ = lookup[j];
    688                     assert (newJ >= 0);
    689                     temp[last] = newJ;
    690                     last = newJ;
    691                     j = next_[j];
    692                }
    693                if (setTemp)
    694                     temp[last] = -(keyVariable_[iSet] + 1);
    695                if (lowerSet_) {
    696                     // we only need to get lower_ and upper_ correct
    697                     double shift = 0.0;
    698                     for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
    699                          if (getDynamicStatus(j) == atUpperBound)
    700                               shift += upperColumn_[j];
    701                          else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
    702                               shift += lowerColumn_[j];
    703                     if (lowerSet_[iSet] > -1.0e20)
    704                          lower_[iSet] = lowerSet_[iSet] - shift;
    705                     if (upperSet_[iSet] < 1.0e20)
    706                          upper_[iSet] = upperSet_[iSet] - shift;
    707                }
    708           }
    709           // move to next_
    710           CoinMemcpyN(temp + firstDynamic_, (firstAvailable_ - firstDynamic_), next_ + firstDynamic_);
    711           // if odd iterations may be one out so adjust currentNumber
    712           currentNumber = CoinMin(currentNumber + 1, lastDynamic_);
    713           // zero solution
    714           CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
    715           // zero cost
    716           CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
    717           // zero lengths
    718           CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
    719           for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
    720                model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
    721                model->setStatus(iColumn, ClpSimplex::atLowerBound);
    722                backward_[iColumn] = -1;
    723           }
    724           delete [] lookup;
    725           delete [] temp;
    726           // make sure fromIndex clean
    727           fromIndex_[0] = -1;
    728           //#define CLP_DEBUG
     577    // lookup array
     578    int *lookup = new int[lastDynamic_];
     579    int iColumn;
     580    int numberColumns = model->numberColumns();
     581    double *element = matrix_->getMutableElements();
     582    int *row = matrix_->getMutableIndices();
     583    CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
     584    int *length = matrix_->getMutableVectorLengths();
     585    double *cost = model->costRegion();
     586    double *lowerColumn = model->lowerRegion();
     587    double *upperColumn = model->upperRegion();
     588    int *pivotVariable = model->pivotVariable();
     589    CoinBigIndex numberElements = startColumn[firstDynamic_];
     590    // first just do lookup and basic stuff
     591    int currentNumber = firstAvailable_;
     592    firstAvailable_ = firstDynamic_;
     593    int numberToDo = 0;
     594    double objectiveChange = 0.0;
     595    double *solution = model->solutionRegion();
     596    for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
     597      int iSet = backward_[iColumn];
     598      if (toIndex_[iSet] < 0) {
     599        toIndex_[iSet] = 0;
     600        fromIndex_[numberToDo++] = iSet;
     601      }
     602      if (model->getStatus(iColumn) == ClpSimplex::basic || iColumn == keyVariable_[iSet]) {
     603        lookup[iColumn] = firstAvailable_;
     604        if (iColumn != keyVariable_[iSet]) {
     605          int iPivot = backToPivotRow_[iColumn];
     606          backToPivotRow_[firstAvailable_] = iPivot;
     607          pivotVariable[iPivot] = firstAvailable_;
     608        }
     609        firstAvailable_++;
     610      } else {
     611        int jColumn = id_[iColumn - firstDynamic_];
     612        setDynamicStatus(jColumn, atLowerBound);
     613        if (lowerColumn_ || upperColumn_) {
     614          if (model->getStatus(iColumn) == ClpSimplex::atUpperBound)
     615            setDynamicStatus(jColumn, atUpperBound);
     616          // treat solution as if exactly at a bound
     617          double value = solution[iColumn];
     618          if (fabs(value - lowerColumn[iColumn]) < fabs(value - upperColumn[iColumn]))
     619            value = lowerColumn[iColumn];
     620          else
     621            value = upperColumn[iColumn];
     622          objectiveChange += cost[iColumn] * value;
     623          // redo lower and upper on sets
     624          double shift = value;
     625          if (lowerSet_[iSet] > -1.0e20)
     626            lower_[iSet] = lowerSet_[iSet] - shift;
     627          if (upperSet_[iSet] < 1.0e20)
     628            upper_[iSet] = upperSet_[iSet] - shift;
     629        }
     630        lookup[iColumn] = -1;
     631      }
     632    }
     633    model->setObjectiveOffset(model->objectiveOffset() + objectiveChange);
     634    firstAvailable_ = firstDynamic_;
     635    for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
     636      if (lookup[iColumn] >= 0) {
     637        // move
     638        int jColumn = id_[iColumn - firstDynamic_];
     639        id_[firstAvailable_ - firstDynamic_] = jColumn;
     640        int numberThis = startColumn_[jColumn + 1] - startColumn_[jColumn];
     641        length[firstAvailable_] = numberThis;
     642        cost[firstAvailable_] = cost[iColumn];
     643        lowerColumn[firstAvailable_] = lowerColumn[iColumn];
     644        upperColumn[firstAvailable_] = upperColumn[iColumn];
     645        double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
     646        double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
     647        if (originalUpper > 1.0e30)
     648          originalUpper = COIN_DBL_MAX;
     649        model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn],
     650          originalLower, originalUpper,
     651          cost_[jColumn]);
     652        CoinBigIndex base = startColumn_[jColumn];
     653        for (int j = 0; j < numberThis; j++) {
     654          row[numberElements] = row_[base + j];
     655          element[numberElements++] = element_[base + j];
     656        }
     657        model->setStatus(firstAvailable_, model->getStatus(iColumn));
     658        backward_[firstAvailable_] = backward_[iColumn];
     659        solution[firstAvailable_] = solution[iColumn];
     660        firstAvailable_++;
     661        startColumn[firstAvailable_] = numberElements;
     662      }
     663    }
     664    // clean up next_
     665    int *temp = new int[firstAvailable_];
     666    for (int jSet = 0; jSet < numberToDo; jSet++) {
     667      int iSet = fromIndex_[jSet];
     668      toIndex_[iSet] = -1;
     669      int last = keyVariable_[iSet];
     670      int j = next_[last];
     671      bool setTemp = true;
     672      if (last < lastDynamic_) {
     673        last = lookup[last];
     674        assert(last >= 0);
     675        keyVariable_[iSet] = last;
     676      } else if (j >= 0) {
     677        int newJ = lookup[j];
     678        assert(newJ >= 0);
     679        j = next_[j];
     680        next_[last] = newJ;
     681        last = newJ;
     682      } else {
     683        next_[last] = -(iSet + numberColumns + 1);
     684        setTemp = false;
     685      }
     686      while (j >= 0) {
     687        int newJ = lookup[j];
     688        assert(newJ >= 0);
     689        temp[last] = newJ;
     690        last = newJ;
     691        j = next_[j];
     692      }
     693      if (setTemp)
     694        temp[last] = -(keyVariable_[iSet] + 1);
     695      if (lowerSet_) {
     696        // we only need to get lower_ and upper_ correct
     697        double shift = 0.0;
     698        for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
     699          if (getDynamicStatus(j) == atUpperBound)
     700            shift += upperColumn_[j];
     701          else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
     702            shift += lowerColumn_[j];
     703        if (lowerSet_[iSet] > -1.0e20)
     704          lower_[iSet] = lowerSet_[iSet] - shift;
     705        if (upperSet_[iSet] < 1.0e20)
     706          upper_[iSet] = upperSet_[iSet] - shift;
     707      }
     708    }
     709    // move to next_
     710    CoinMemcpyN(temp + firstDynamic_, (firstAvailable_ - firstDynamic_), next_ + firstDynamic_);
     711    // if odd iterations may be one out so adjust currentNumber
     712    currentNumber = CoinMin(currentNumber + 1, lastDynamic_);
     713    // zero solution
     714    CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
     715    // zero cost
     716    CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
     717    // zero lengths
     718    CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
     719    for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
     720      model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
     721      model->setStatus(iColumn, ClpSimplex::atLowerBound);
     722      backward_[iColumn] = -1;
     723    }
     724    delete[] lookup;
     725    delete[] temp;
     726    // make sure fromIndex clean
     727    fromIndex_[0] = -1;
     728    //#define CLP_DEBUG
    729729#ifdef CLP_DEBUG
    730           // debug
    731           {
    732                int i;
    733                int numberRows = model->numberRows();
    734                char * xxxx = new char[numberColumns];
    735                memset(xxxx, 0, numberColumns);
    736                for (i = 0; i < numberRows; i++) {
    737                     int iPivot = pivotVariable[i];
    738                     assert (model->getStatus(iPivot) == ClpSimplex::basic);
    739                     if (iPivot < numberColumns && backward_[iPivot] >= 0)
    740                          xxxx[iPivot] = 1;
    741                }
    742                for (i = 0; i < numberSets_; i++) {
    743                     int key = keyVariable_[i];
    744                     int iColumn = next_[key];
    745                     int k = 0;
    746                     while(iColumn >= 0) {
    747                          k++;
    748                          assert (k < 100);
    749                          assert (backward_[iColumn] == i);
    750                          iColumn = next_[iColumn];
    751                     }
    752                     int stop = -(key + 1);
    753                     while (iColumn != stop) {
    754                          assert (iColumn < 0);
    755                          iColumn = -iColumn - 1;
    756                          k++;
    757                          assert (k < 100);
    758                          assert (backward_[iColumn] == i);
    759                          iColumn = next_[iColumn];
    760                     }
    761                     iColumn = next_[key];
    762                     while (iColumn >= 0) {
    763                          assert (xxxx[iColumn]);
    764                          xxxx[iColumn] = 0;
    765                          iColumn = next_[iColumn];
    766                     }
    767                }
    768                for (i = 0; i < numberColumns; i++) {
    769                     if (i < numberColumns && backward_[i] >= 0) {
    770                          assert (!xxxx[i] || i == keyVariable_[backward_[i]]);
    771                     }
    772                }
    773                delete [] xxxx;
    774           }
    775           {
    776                for (int i = 0; i < numberSets_; i++)
    777                     assert(toIndex_[i] == -1);
    778           }
     730    // debug
     731    {
     732      int i;
     733      int numberRows = model->numberRows();
     734      char *xxxx = new char[numberColumns];
     735      memset(xxxx, 0, numberColumns);
     736      for (i = 0; i < numberRows; i++) {
     737        int iPivot = pivotVariable[i];
     738        assert(model->getStatus(iPivot) == ClpSimplex::basic);
     739        if (iPivot < numberColumns && backward_[iPivot] >= 0)
     740          xxxx[iPivot] = 1;
     741      }
     742      for (i = 0; i < numberSets_; i++) {
     743        int key = keyVariable_[i];
     744        int iColumn = next_[key];
     745        int k = 0;
     746        while (iColumn >= 0) {
     747          k++;
     748          assert(k < 100);
     749          assert(backward_[iColumn] == i);
     750          iColumn = next_[iColumn];
     751        }
     752        int stop = -(key + 1);
     753        while (iColumn != stop) {
     754          assert(iColumn < 0);
     755          iColumn = -iColumn - 1;
     756          k++;
     757          assert(k < 100);
     758          assert(backward_[iColumn] == i);
     759          iColumn = next_[iColumn];
     760        }
     761        iColumn = next_[key];
     762        while (iColumn >= 0) {
     763          assert(xxxx[iColumn]);
     764          xxxx[iColumn] = 0;
     765          iColumn = next_[iColumn];
     766        }
     767      }
     768      for (i = 0; i < numberColumns; i++) {
     769        if (i < numberColumns && backward_[i] >= 0) {
     770          assert(!xxxx[i] || i == keyVariable_[backward_[i]]);
     771        }
     772      }
     773      delete[] xxxx;
     774    }
     775    {
     776      for (int i = 0; i < numberSets_; i++)
     777        assert(toIndex_[i] == -1);
     778    }
    779779#endif
    780           savedFirstAvailable_ = firstAvailable_;
    781      }
    782      break;
    783      // flag a variable
    784      case 1: {
    785           // id will be sitting at firstAvailable
    786           int sequence = id_[firstAvailable_-firstDynamic_];
    787           assert (!flagged(sequence));
    788           setFlagged(sequence);
    789           model->clearFlagged(firstAvailable_);
    790      }
    791      break;
    792      // unflag all variables
    793      case 2: {
    794           for (int i = 0; i < numberGubColumns_; i++) {
    795                if (flagged(i)) {
    796                     unsetFlagged(i);
    797                     returnNumber++;
    798                }
    799           }
    800      }
    801      break;
    802      //  just reset costs and bounds (primal)
    803      case 3: {
    804           double * cost = model->costRegion();
    805           double * solution = model->solutionRegion();
    806           double * lowerColumn = model->columnLower();
    807           double * upperColumn = model->columnUpper();
    808           for (int i = firstDynamic_; i < firstAvailable_; i++) {
    809                int jColumn = id_[i-firstDynamic_];
    810                cost[i] = cost_[jColumn];
    811                if (!lowerColumn_ && !upperColumn_) {
    812                     lowerColumn[i] = 0.0;
    813                     upperColumn[i] = COIN_DBL_MAX;
    814                }  else {
    815                     if (lowerColumn_)
    816                          lowerColumn[i] = lowerColumn_[jColumn];
    817                     else
    818                          lowerColumn[i] = 0.0;
    819                     if (upperColumn_)
    820                          upperColumn[i] = upperColumn_[jColumn];
    821                     else
    822                          upperColumn[i] = COIN_DBL_MAX;
    823                }
    824                if (model->nonLinearCost())
    825                     model->nonLinearCost()->setOne(i, solution[i],
    826                                                    lowerColumn[i],
    827                                                    upperColumn[i], cost_[jColumn]);
    828           }
    829           if (!model->numberIterations() && rhsOffset_) {
    830                lastRefresh_ = - refreshFrequency_; // force refresh
    831           }
    832      }
    833      break;
    834      // and get statistics for column generation
    835      case 4: {
    836           // In theory we should subtract out ones we have done but ....
    837           // If key slack then dual 0.0
    838           // If not then slack could be dual infeasible
    839           // dj for key is zero so that defines dual on set
    840           int i;
    841           int numberColumns = model->numberColumns();
    842           double * dual = model->dualRowSolution();
    843           double infeasibilityCost = model->infeasibilityCost();
    844           double dualTolerance = model->dualTolerance();
    845           double relaxedTolerance = dualTolerance;
    846           // we can't really trust infeasibilities if there is dual error
    847           double error = CoinMin(1.0e-2, model->largestDualError());
    848           // allow tolerance at least slightly bigger than standard
    849           relaxedTolerance = relaxedTolerance +  error;
    850           // but we will be using difference
    851           relaxedTolerance -= dualTolerance;
    852           double objectiveOffset = 0.0;
    853           for (i = 0; i < numberSets_; i++) {
    854                int kColumn = keyVariable_[i];
    855                double value = 0.0;
    856                if (kColumn < numberColumns) {
    857                     kColumn = id_[kColumn-firstDynamic_];
    858                     // dj without set
    859                     value = cost_[kColumn];
    860                     for (CoinBigIndex j = startColumn_[kColumn];
    861                               j < startColumn_[kColumn+1]; j++) {
    862                          int iRow = row_[j];
    863                          value -= dual[iRow] * element_[j];
    864                     }
    865                     double infeasibility = 0.0;
    866                     if (getStatus(i) == ClpSimplex::atLowerBound) {
    867                          if (-value > dualTolerance)
    868                               infeasibility = -value - dualTolerance;
    869                     } else if (getStatus(i) == ClpSimplex::atUpperBound) {
    870                          if (value > dualTolerance)
    871                               infeasibility = value - dualTolerance;
    872                     }
    873                     if (infeasibility > 0.0) {
    874                          sumDualInfeasibilities_ += infeasibility;
    875                          if (infeasibility > relaxedTolerance)
    876                               sumOfRelaxedDualInfeasibilities_ += infeasibility;
    877                          numberDualInfeasibilities_ ++;
    878                     }
    879                } else {
    880                     // slack key - may not be feasible
    881                     assert (getStatus(i) == ClpSimplex::basic);
    882                     // negative as -1.0 for slack
    883                     value = -weight(i) * infeasibilityCost;
    884                }
    885                // Now subtract out from all
    886                for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i+1]; k++) {
    887                     if (getDynamicStatus(k) != inSmall) {
    888                          double djValue = cost_[k] - value;
    889                          for (CoinBigIndex j = startColumn_[k];
    890                                    j < startColumn_[k+1]; j++) {
    891                               int iRow = row_[j];
    892                               djValue -= dual[iRow] * element_[j];
    893                          }
    894                          double infeasibility = 0.0;
    895                          double shift = 0.0;
    896                          if (getDynamicStatus(k) == atLowerBound) {
    897                               if (lowerColumn_)
    898                                    shift = lowerColumn_[k];
    899                               if (djValue < -dualTolerance)
    900                                    infeasibility = -djValue - dualTolerance;
    901                          } else {
    902                               // at upper bound
    903                               shift = upperColumn_[k];
    904                               if (djValue > dualTolerance)
    905                                    infeasibility = djValue - dualTolerance;
    906                          }
    907                          objectiveOffset += shift * cost_[k];
    908                          if (infeasibility > 0.0) {
    909                               sumDualInfeasibilities_ += infeasibility;
    910                               if (infeasibility > relaxedTolerance)
    911                                    sumOfRelaxedDualInfeasibilities_ += infeasibility;
    912                               numberDualInfeasibilities_ ++;
    913                          }
    914                     }
    915                }
    916           }
    917           model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
    918      }
    919      break;
    920      // see if time to re-factorize
    921      case 5: {
    922           if (firstAvailable_ > numberSets_ + model->numberRows() + model->factorizationFrequency())
    923                returnNumber = 4;
    924      }
    925      break;
    926      // return 1 if there may be changing bounds on variable (column generation)
    927      case 6: {
    928           returnNumber = (lowerColumn_ != NULL || upperColumn_ != NULL) ? 1 : 0;
     780    savedFirstAvailable_ = firstAvailable_;
     781  } break;
     782  // flag a variable
     783  case 1: {
     784    // id will be sitting at firstAvailable
     785    int sequence = id_[firstAvailable_ - firstDynamic_];
     786    assert(!flagged(sequence));
     787    setFlagged(sequence);
     788    model->clearFlagged(firstAvailable_);
     789  } break;
     790  // unflag all variables
     791  case 2: {
     792    for (int i = 0; i < numberGubColumns_; i++) {
     793      if (flagged(i)) {
     794        unsetFlagged(i);
     795        returnNumber++;
     796      }
     797    }
     798  } break;
     799  //  just reset costs and bounds (primal)
     800  case 3: {
     801    double *cost = model->costRegion();
     802    double *solution = model->solutionRegion();
     803    double *lowerColumn = model->columnLower();
     804    double *upperColumn = model->columnUpper();
     805    for (int i = firstDynamic_; i < firstAvailable_; i++) {
     806      int jColumn = id_[i - firstDynamic_];
     807      cost[i] = cost_[jColumn];
     808      if (!lowerColumn_ && !upperColumn_) {
     809        lowerColumn[i] = 0.0;
     810        upperColumn[i] = COIN_DBL_MAX;
     811      } else {
     812        if (lowerColumn_)
     813          lowerColumn[i] = lowerColumn_[jColumn];
     814        else
     815          lowerColumn[i] = 0.0;
     816        if (upperColumn_)
     817          upperColumn[i] = upperColumn_[jColumn];
     818        else
     819          upperColumn[i] = COIN_DBL_MAX;
     820      }
     821      if (model->nonLinearCost())
     822        model->nonLinearCost()->setOne(i, solution[i],
     823          lowerColumn[i],
     824          upperColumn[i], cost_[jColumn]);
     825    }
     826    if (!model->numberIterations() && rhsOffset_) {
     827      lastRefresh_ = -refreshFrequency_; // force refresh
     828    }
     829  } break;
     830  // and get statistics for column generation
     831  case 4: {
     832    // In theory we should subtract out ones we have done but ....
     833    // If key slack then dual 0.0
     834    // If not then slack could be dual infeasible
     835    // dj for key is zero so that defines dual on set
     836    int i;
     837    int numberColumns = model->numberColumns();
     838    double *dual = model->dualRowSolution();
     839    double infeasibilityCost = model->infeasibilityCost();
     840    double dualTolerance = model->dualTolerance();
     841    double relaxedTolerance = dualTolerance;
     842    // we can't really trust infeasibilities if there is dual error
     843    double error = CoinMin(1.0e-2, model->largestDualError());
     844    // allow tolerance at least slightly bigger than standard
     845    relaxedTolerance = relaxedTolerance + error;
     846    // but we will be using difference
     847    relaxedTolerance -= dualTolerance;
     848    double objectiveOffset = 0.0;
     849    for (i = 0; i < numberSets_; i++) {
     850      int kColumn = keyVariable_[i];
     851      double value = 0.0;
     852      if (kColumn < numberColumns) {
     853        kColumn = id_[kColumn - firstDynamic_];
     854        // dj without set
     855        value = cost_[kColumn];
     856        for (CoinBigIndex j = startColumn_[kColumn];
     857             j < startColumn_[kColumn + 1]; j++) {
     858          int iRow = row_[j];
     859          value -= dual[iRow] * element_[j];
     860        }
     861        double infeasibility = 0.0;
     862        if (getStatus(i) == ClpSimplex::atLowerBound) {
     863          if (-value > dualTolerance)
     864            infeasibility = -value - dualTolerance;
     865        } else if (getStatus(i) == ClpSimplex::atUpperBound) {
     866          if (value > dualTolerance)
     867            infeasibility = value - dualTolerance;
     868        }
     869        if (infeasibility > 0.0) {
     870          sumDualInfeasibilities_ += infeasibility;
     871          if (infeasibility > relaxedTolerance)
     872            sumOfRelaxedDualInfeasibilities_ += infeasibility;
     873          numberDualInfeasibilities_++;
     874        }
     875      } else {
     876        // slack key - may not be feasible
     877        assert(getStatus(i) == ClpSimplex::basic);
     878        // negative as -1.0 for slack
     879        value = -weight(i) * infeasibilityCost;
     880      }
     881      // Now subtract out from all
     882      for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i + 1]; k++) {
     883        if (getDynamicStatus(k) != inSmall) {
     884          double djValue = cost_[k] - value;
     885          for (CoinBigIndex j = startColumn_[k];
     886               j < startColumn_[k + 1]; j++) {
     887            int iRow = row_[j];
     888            djValue -= dual[iRow] * element_[j];
     889          }
     890          double infeasibility = 0.0;
     891          double shift = 0.0;
     892          if (getDynamicStatus(k) == atLowerBound) {
     893            if (lowerColumn_)
     894              shift = lowerColumn_[k];
     895            if (djValue < -dualTolerance)
     896              infeasibility = -djValue - dualTolerance;
     897          } else {
     898            // at upper bound
     899            shift = upperColumn_[k];
     900            if (djValue > dualTolerance)
     901              infeasibility = djValue - dualTolerance;
     902          }
     903          objectiveOffset += shift * cost_[k];
     904          if (infeasibility > 0.0) {
     905            sumDualInfeasibilities_ += infeasibility;
     906            if (infeasibility > relaxedTolerance)
     907              sumOfRelaxedDualInfeasibilities_ += infeasibility;
     908            numberDualInfeasibilities_++;
     909          }
     910        }
     911      }
     912    }
     913    model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
     914  } break;
     915  // see if time to re-factorize
     916  case 5: {
     917    if (firstAvailable_ > numberSets_ + model->numberRows() + model->factorizationFrequency())
     918      returnNumber = 4;
     919  } break;
     920  // return 1 if there may be changing bounds on variable (column generation)
     921  case 6: {
     922    returnNumber = (lowerColumn_ != NULL || upperColumn_ != NULL) ? 1 : 0;
    929923#if 0
    930924          if (!returnNumber) {
     
    938932          }
    939933#endif
    940      }
    941      break;
    942      // restore firstAvailable_
    943      case 7: {
    944           int iColumn;
    945           int * length = matrix_->getMutableVectorLengths();
    946           double * cost = model->costRegion();
    947           double * solution = model->solutionRegion();
    948           int currentNumber = firstAvailable_;
    949           firstAvailable_ = savedFirstAvailable_;
    950           // zero solution
    951           CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
    952           // zero cost
    953           CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
    954           // zero lengths
    955           CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
    956           for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
    957                model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
    958                model->setStatus(iColumn, ClpSimplex::atLowerBound);
    959                backward_[iColumn] = -1;
    960           }
    961      }
    962      break;
    963      // make sure set is clean
    964      case 8: {
    965           int sequenceIn = model->sequenceIn();
    966           if (sequenceIn < model->numberColumns()) {
    967                int iSet = backward_[sequenceIn];
    968                if (iSet >= 0 && lowerSet_) {
    969                     // we only need to get lower_ and upper_ correct
    970                     double shift = 0.0;
    971                     for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
    972                          if (getDynamicStatus(j) == atUpperBound)
    973                               shift += upperColumn_[j];
    974                          else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
    975                               shift += lowerColumn_[j];
    976                     if (lowerSet_[iSet] > -1.0e20)
    977                          lower_[iSet] = lowerSet_[iSet] - shift;
    978                     if (upperSet_[iSet] < 1.0e20)
    979                          upper_[iSet] = upperSet_[iSet] - shift;
    980                }
    981                if (sequenceIn == firstAvailable_) {
    982                     // not really in small problem
    983                     int iBig = id_[sequenceIn-firstDynamic_];
    984                     if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound)
    985                          setDynamicStatus(iBig, atLowerBound);
    986                     else
    987                          setDynamicStatus(iBig, atUpperBound);
    988                }
    989           }
    990      }
    991      break;
    992      // adjust lower,upper
    993      case 9: {
    994           int sequenceIn = model->sequenceIn();
    995           if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_ && lowerSet_) {
    996                int iSet = backward_[sequenceIn];
    997                assert (iSet >= 0);
    998                int inBig = id_[sequenceIn-firstDynamic_];
    999                const double * solution = model->solutionRegion();
    1000                setDynamicStatus(inBig, inSmall);
    1001                if (lowerSet_[iSet] > -1.0e20)
    1002                     lower_[iSet] += solution[sequenceIn];
    1003                if (upperSet_[iSet] < 1.0e20)
    1004                     upper_[iSet] += solution[sequenceIn];
    1005                model->setObjectiveOffset(model->objectiveOffset() -
    1006                                          solution[sequenceIn]*cost_[inBig]);
    1007           }
    1008      }
    1009      }
    1010      return returnNumber;
     934  } break;
     935  // restore firstAvailable_
     936  case 7: {
     937    int iColumn;
     938    int *length = matrix_->getMutableVectorLengths();
     939    double *cost = model->costRegion();
     940    double *solution = model->solutionRegion();
     941    int currentNumber = firstAvailable_;
     942    firstAvailable_ = savedFirstAvailable_;
     943    // zero solution
     944    CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
     945    // zero cost
     946    CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
     947    // zero lengths
     948    CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
     949    for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
     950      model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
     951      model->setStatus(iColumn, ClpSimplex::atLowerBound);
     952      backward_[iColumn] = -1;
     953    }
     954  } break;
     955  // make sure set is clean
     956  case 8: {
     957    int sequenceIn = model->sequenceIn();
     958    if (sequenceIn < model->numberColumns()) {
     959      int iSet = backward_[sequenceIn];
     960      if (iSet >= 0 && lowerSet_) {
     961        // we only need to get lower_ and upper_ correct
     962        double shift = 0.0;
     963        for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
     964          if (getDynamicStatus(j) == atUpperBound)
     965            shift += upperColumn_[j];
     966          else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
     967            shift += lowerColumn_[j];
     968        if (lowerSet_[iSet] > -1.0e20)
     969          lower_[iSet] = lowerSet_[iSet] - shift;
     970        if (upperSet_[iSet] < 1.0e20)
     971          upper_[iSet] = upperSet_[iSet] - shift;
     972      }
     973      if (sequenceIn == firstAvailable_) {
     974        // not really in small problem
     975        int iBig = id_[sequenceIn - firstDynamic_];
     976        if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound)
     977          setDynamicStatus(iBig, atLowerBound);
     978        else
     979          setDynamicStatus(iBig, atUpperBound);
     980      }
     981    }
     982  } break;
     983  // adjust lower,upper
     984  case 9: {
     985    int sequenceIn = model->sequenceIn();
     986    if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_ && lowerSet_) {
     987      int iSet = backward_[sequenceIn];
     988      assert(iSet >= 0);
     989      int inBig = id_[sequenceIn - firstDynamic_];
     990      const double *solution = model->solutionRegion();
     991      setDynamicStatus(inBig, inSmall);
     992      if (lowerSet_[iSet] > -1.0e20)
     993        lower_[iSet] += solution[sequenceIn];
     994      if (upperSet_[iSet] < 1.0e20)
     995        upper_[iSet] += solution[sequenceIn];
     996      model->setObjectiveOffset(model->objectiveOffset() - solution[sequenceIn] * cost_[inBig]);
     997    }
     998  }
     999  }
     1000  return returnNumber;
    10111001}
    10121002// Add a new variable to a set
    1013 void
    1014 ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
     1003void ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
    10151004{
    1016      int last = keyVariable_[iSet];
    1017      int j = next_[last];
    1018      while (j >= 0) {
    1019           last = j;
    1020           j = next_[j];
    1021      }
    1022      next_[last] = -(sequence + 1);
    1023      next_[sequence] = j;
     1005  int last = keyVariable_[iSet];
     1006  int j = next_[last];
     1007  while (j >= 0) {
     1008    last = j;
     1009    j = next_[j];
     1010  }
     1011  next_[last] = -(sequence + 1);
     1012  next_[sequence] = j;
    10241013}
    10251014// Sets up an effective RHS and does gub crash if needed
    1026 void
    1027 ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
     1015void ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex *model, bool cheapest)
    10281016{
    1029      // Do basis - cheapest or slack if feasible (unless cheapest set)
    1030      int longestSet = 0;
    1031      int iSet;
    1032      for (iSet = 0; iSet < numberSets_; iSet++)
    1033           longestSet = CoinMax(longestSet, fullStart_[iSet+1] - fullStart_[iSet]);
     1017  // Do basis - cheapest or slack if feasible (unless cheapest set)
     1018  int longestSet = 0;
     1019  int iSet;
     1020  for (iSet = 0; iSet < numberSets_; iSet++)
     1021    longestSet = CoinMax(longestSet, fullStart_[iSet + 1] - fullStart_[iSet]);
    10341022
    1035      double * upper = new double[longestSet+1];
    1036      double * cost = new double[longestSet+1];
    1037      double * lower = new double[longestSet+1];
    1038      double * solution = new double[longestSet+1];
    1039      assert (!next_);
    1040      delete [] next_;
    1041      int numberColumns = model->numberColumns();
    1042      next_ = new int[numberColumns+numberSets_+CoinMax(2*longestSet, lastDynamic_-firstDynamic_)];
    1043      char * mark = new char[numberColumns];
    1044      memset(mark, 0, numberColumns);
    1045      for (int iColumn = 0; iColumn < numberColumns; iColumn++)
    1046           next_[iColumn] = COIN_INT_MAX;
    1047      int i;
    1048      int * keys = new int[numberSets_];
    1049      int * back = new int[numberGubColumns_];
    1050      CoinFillN(back, numberGubColumns_, -1);
    1051      for (i = 0; i < numberSets_; i++)
    1052           keys[i] = COIN_INT_MAX;
    1053      delete [] dynamicStatus_;
    1054      dynamicStatus_ = new unsigned char [numberGubColumns_];
    1055      memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
    1056      for (i = 0; i < numberGubColumns_; i++)
    1057           setDynamicStatus(i, atLowerBound);
    1058      // set up chains
    1059      for (i = firstDynamic_; i < lastDynamic_; i++) {
    1060           if (id_[i-firstDynamic_] >= 0) {
    1061                if (model->getStatus(i) == ClpSimplex::basic)
    1062                     mark[i] = 1;
    1063                int iSet = backward_[i];
    1064                assert (iSet >= 0);
    1065                int iNext = keys[iSet];
    1066                next_[i] = iNext;
    1067                keys[iSet] = i;
    1068                back[id_[i-firstDynamic_]] = i;
     1023  double *upper = new double[longestSet + 1];
     1024  double *cost = new double[longestSet + 1];
     1025  double *lower = new double[longestSet + 1];
     1026  double *solution = new double[longestSet + 1];
     1027  assert(!next_);
     1028  delete[] next_;
     1029  int numberColumns = model->numberColumns();
     1030  next_ = new int[numberColumns + numberSets_ + CoinMax(2 * longestSet, lastDynamic_ - firstDynamic_)];
     1031  char *mark = new char[numberColumns];
     1032  memset(mark, 0, numberColumns);
     1033  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
     1034    next_[iColumn] = COIN_INT_MAX;
     1035  int i;
     1036  int *keys = new int[numberSets_];
     1037  int *back = new int[numberGubColumns_];
     1038  CoinFillN(back, numberGubColumns_, -1);
     1039  for (i = 0; i < numberSets_; i++)
     1040    keys[i] = COIN_INT_MAX;
     1041  delete[] dynamicStatus_;
     1042  dynamicStatus_ = new unsigned char[numberGubColumns_];
     1043  memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
     1044  for (i = 0; i < numberGubColumns_; i++)
     1045    setDynamicStatus(i, atLowerBound);
     1046  // set up chains
     1047  for (i = firstDynamic_; i < lastDynamic_; i++) {
     1048    if (id_[i - firstDynamic_] >= 0) {
     1049      if (model->getStatus(i) == ClpSimplex::basic)
     1050        mark[i] = 1;
     1051      int iSet = backward_[i];
     1052      assert(iSet >= 0);
     1053      int iNext = keys[iSet];
     1054      next_[i] = iNext;
     1055      keys[iSet] = i;
     1056      back[id_[i - firstDynamic_]] = i;
     1057    } else {
     1058      model->setStatus(i, ClpSimplex::atLowerBound);
     1059      backward_[i] = -1;
     1060    }
     1061  }
     1062  double *columnSolution = model->solutionRegion();
     1063  int numberRows = getNumRows();
     1064  toIndex_ = new int[numberSets_];
     1065  for (iSet = 0; iSet < numberSets_; iSet++)
     1066    toIndex_[iSet] = -1;
     1067  fromIndex_ = new int[numberRows + numberSets_];
     1068  double tolerance = model->primalTolerance();
     1069  double *element = matrix_->getMutableElements();
     1070  int *row = matrix_->getMutableIndices();
     1071  CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
     1072  int *length = matrix_->getMutableVectorLengths();
     1073  double objectiveOffset = 0.0;
     1074  for (iSet = 0; iSet < numberSets_; iSet++) {
     1075    int j;
     1076    int numberBasic = 0;
     1077    int iBasic = -1;
     1078    int iStart = fullStart_[iSet];
     1079    int iEnd = fullStart_[iSet + 1];
     1080    // find one with smallest length
     1081    int smallest = numberRows + 1;
     1082    double value = 0.0;
     1083    j = keys[iSet];
     1084    while (j != COIN_INT_MAX) {
     1085      if (model->getStatus(j) == ClpSimplex::basic) {
     1086        if (length[j] < smallest) {
     1087          smallest = length[j];
     1088          iBasic = j;
     1089        }
     1090        numberBasic++;
     1091      }
     1092      value += columnSolution[j];
     1093      j = next_[j];
     1094    }
     1095    bool done = false;
     1096    if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
     1097      if (getStatus(iSet) == ClpSimplex::basic)
     1098        iBasic = iSet + numberColumns; // slack key - use
     1099      done = true;
     1100    } else if (numberBasic == 1) {
     1101      // see if can be key
     1102      double thisSolution = columnSolution[iBasic];
     1103      if (thisSolution < 0.0) {
     1104        value -= thisSolution;
     1105        thisSolution = 0.0;
     1106        columnSolution[iBasic] = thisSolution;
     1107      }
     1108      // try setting slack to a bound
     1109      assert(upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
     1110      double cost1 = COIN_DBL_MAX;
     1111      int whichBound = -1;
     1112      if (upper_[iSet] < 1.0e20) {
     1113        // try slack at ub
     1114        double newBasic = thisSolution + upper_[iSet] - value;
     1115        if (newBasic >= -tolerance) {
     1116          // can go
     1117          whichBound = 1;
     1118          cost1 = newBasic * cost_[iBasic];
     1119          // But if exact then may be good solution
     1120          if (fabs(upper_[iSet] - value) < tolerance)
     1121            cost1 = -COIN_DBL_MAX;
     1122        }
     1123      }
     1124      if (lower_[iSet] > -1.0e20) {
     1125        // try slack at lb
     1126        double newBasic = thisSolution + lower_[iSet] - value;
     1127        if (newBasic >= -tolerance) {
     1128          // can go but is it cheaper
     1129          double cost2 = newBasic * cost_[iBasic];
     1130          // But if exact then may be good solution
     1131          if (fabs(lower_[iSet] - value) < tolerance)
     1132            cost2 = -COIN_DBL_MAX;
     1133          if (cost2 < cost1)
     1134            whichBound = 0;
     1135        }
     1136      }
     1137      if (whichBound != -1) {
     1138        // key
     1139        done = true;
     1140        if (whichBound) {
     1141          // slack to upper
     1142          columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
     1143          setStatus(iSet, ClpSimplex::atUpperBound);
     1144        } else {
     1145          // slack to lower
     1146          columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
     1147          setStatus(iSet, ClpSimplex::atLowerBound);
     1148        }
     1149      }
     1150    }
     1151    if (!done) {
     1152      if (!cheapest) {
     1153        // see if slack can be key
     1154        if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
     1155          done = true;
     1156          setStatus(iSet, ClpSimplex::basic);
     1157          iBasic = iSet + numberColumns;
     1158        }
     1159      }
     1160      if (!done) {
     1161        // set non basic if there was one
     1162        if (iBasic >= 0)
     1163          model->setStatus(iBasic, ClpSimplex::atLowerBound);
     1164        // find cheapest
     1165        int numberInSet = iEnd - iStart;
     1166        if (!lowerColumn_) {
     1167          CoinZeroN(lower, numberInSet);
     1168        } else {
     1169          for (int j = 0; j < numberInSet; j++)
     1170            lower[j] = lowerColumn_[j + iStart];
     1171        }
     1172        if (!upperColumn_) {
     1173          CoinFillN(upper, numberInSet, COIN_DBL_MAX);
     1174        } else {
     1175          for (int j = 0; j < numberInSet; j++)
     1176            upper[j] = upperColumn_[j + iStart];
     1177        }
     1178        CoinFillN(solution, numberInSet, 0.0);
     1179        // and slack
     1180        iBasic = numberInSet;
     1181        solution[iBasic] = -value;
     1182        lower[iBasic] = -upper_[iSet];
     1183        upper[iBasic] = -lower_[iSet];
     1184        int kphase;
     1185        if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
     1186          // feasible
     1187          kphase = 1;
     1188          cost[iBasic] = 0.0;
     1189          for (int j = 0; j < numberInSet; j++)
     1190            cost[j] = cost_[j + iStart];
     1191        } else {
     1192          // infeasible
     1193          kphase = 0;
     1194          // remember bounds are flipped so opposite to natural
     1195          if (value < lower_[iSet] - tolerance)
     1196            cost[iBasic] = 1.0;
     1197          else
     1198            cost[iBasic] = -1.0;
     1199          CoinZeroN(cost, numberInSet);
     1200        }
     1201        double dualTolerance = model->dualTolerance();
     1202        for (int iphase = kphase; iphase < 2; iphase++) {
     1203          if (iphase) {
     1204            cost[numberInSet] = 0.0;
     1205            for (int j = 0; j < numberInSet; j++)
     1206              cost[j] = cost_[j + iStart];
     1207          }
     1208          // now do one row lp
     1209          bool improve = true;
     1210          while (improve) {
     1211            improve = false;
     1212            double dual = cost[iBasic];
     1213            int chosen = -1;
     1214            double best = dualTolerance;
     1215            int way = 0;
     1216            for (int i = 0; i <= numberInSet; i++) {
     1217              double dj = cost[i] - dual;
     1218              double improvement = 0.0;
     1219              if (iphase || i < numberInSet)
     1220                assert(solution[i] >= lower[i] && solution[i] <= upper[i]);
     1221              if (dj > dualTolerance)
     1222                improvement = dj * (solution[i] - lower[i]);
     1223              else if (dj < -dualTolerance)
     1224                improvement = dj * (solution[i] - upper[i]);
     1225              if (improvement > best) {
     1226                best = improvement;
     1227                chosen = i;
     1228                if (dj < 0.0) {
     1229                  way = 1;
     1230                } else {
     1231                  way = -1;
     1232                }
     1233              }
     1234            }
     1235            if (chosen >= 0) {
     1236              improve = true;
     1237              // now see how far
     1238              if (way > 0) {
     1239                // incoming increasing so basic decreasing
     1240                // if phase 0 then go to nearest bound
     1241                double distance = upper[chosen] - solution[chosen];
     1242                double basicDistance;
     1243                if (!iphase) {
     1244                  assert(iBasic == numberInSet);
     1245                  assert(solution[iBasic] > upper[iBasic]);
     1246                  basicDistance = solution[iBasic] - upper[iBasic];
     1247                } else {
     1248                  basicDistance = solution[iBasic] - lower[iBasic];
     1249                }
     1250                // need extra coding for unbounded
     1251                assert(CoinMin(distance, basicDistance) < 1.0e20);
     1252                if (distance > basicDistance) {
     1253                  // incoming becomes basic
     1254                  solution[chosen] += basicDistance;
     1255                  if (!iphase)
     1256                    solution[iBasic] = upper[iBasic];
     1257                  else
     1258                    solution[iBasic] = lower[iBasic];
     1259                  iBasic = chosen;
     1260                } else {
     1261                  // flip
     1262                  solution[chosen] = upper[chosen];
     1263                  solution[iBasic] -= distance;
     1264                }
     1265              } else {
     1266                // incoming decreasing so basic increasing
     1267                // if phase 0 then go to nearest bound
     1268                double distance = solution[chosen] - lower[chosen];
     1269                double basicDistance;
     1270                if (!iphase) {
     1271                  assert(iBasic == numberInSet);
     1272                  assert(solution[iBasic] < lower[iBasic]);
     1273                  basicDistance = lower[iBasic] - solution[iBasic];
     1274                } else {
     1275                  basicDistance = upper[iBasic] - solution[iBasic];
     1276                }
     1277                // need extra coding for unbounded - for now just exit
     1278                if (CoinMin(distance, basicDistance) > 1.0e20) {
     1279                  printf("unbounded on set %d\n", iSet);
     1280                  iphase = 1;
     1281                  iBasic = numberInSet;
     1282                  break;
     1283                }
     1284                if (distance > basicDistance) {
     1285                  // incoming becomes basic
     1286                  solution[chosen] -= basicDistance;
     1287                  if (!iphase)
     1288                    solution[iBasic] = lower[iBasic];
     1289                  else
     1290                    solution[iBasic] = upper[iBasic];
     1291                  iBasic = chosen;
     1292                } else {
     1293                  // flip
     1294                  solution[chosen] = lower[chosen];
     1295                  solution[iBasic] += distance;
     1296                }
     1297              }
     1298              if (!iphase) {
     1299                if (iBasic < numberInSet)
     1300                  break; // feasible
     1301                else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic])
     1302                  break; // feasible (on flip)
     1303              }
     1304            }
     1305          }
     1306        }
     1307        // do solution i.e. bounds
     1308        if (lowerColumn_ || upperColumn_) {
     1309          for (int j = 0; j < numberInSet; j++) {
     1310            if (j != iBasic) {
     1311              objectiveOffset += solution[j] * cost[j];
     1312              if (lowerColumn_ && upperColumn_) {
     1313                if (fabs(solution[j] - lowerColumn_[j + iStart]) > fabs(solution[j] - upperColumn_[j + iStart]))
     1314                  setDynamicStatus(j + iStart, atUpperBound);
     1315              } else if (upperColumn_ && solution[j] > 0.0) {
     1316                setDynamicStatus(j + iStart, atUpperBound);
     1317              } else {
     1318                setDynamicStatus(j + iStart, atLowerBound);
     1319              }
     1320            }
     1321          }
     1322        }
     1323        // convert iBasic back and do bounds
     1324        if (iBasic == numberInSet) {
     1325          // slack basic
     1326          setStatus(iSet, ClpSimplex::basic);
     1327          iBasic = iSet + numberColumns;
     1328        } else {
     1329          iBasic += fullStart_[iSet];
     1330          if (back[iBasic] >= 0) {
     1331            // exists
     1332            iBasic = back[iBasic];
    10691333          } else {
    1070                model->setStatus(i, ClpSimplex::atLowerBound);
    1071                backward_[i] = -1;
    1072           }
    1073      }
    1074      double * columnSolution = model->solutionRegion();
    1075      int numberRows = getNumRows();
    1076      toIndex_ = new int[numberSets_];
    1077      for (iSet = 0; iSet < numberSets_; iSet++)
    1078           toIndex_[iSet] = -1;
    1079      fromIndex_ = new int [numberRows+numberSets_];
    1080      double tolerance = model->primalTolerance();
    1081      double * element =  matrix_->getMutableElements();
    1082      int * row = matrix_->getMutableIndices();
    1083      CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
    1084      int * length = matrix_->getMutableVectorLengths();
    1085      double objectiveOffset = 0.0;
    1086      for (iSet = 0; iSet < numberSets_; iSet++) {
    1087           int j;
    1088           int numberBasic = 0;
    1089           int iBasic = -1;
    1090           int iStart = fullStart_[iSet];
    1091           int iEnd = fullStart_[iSet+1];
    1092           // find one with smallest length
    1093           int smallest = numberRows + 1;
    1094           double value = 0.0;
    1095           j = keys[iSet];
    1096           while (j != COIN_INT_MAX) {
    1097                if (model->getStatus(j) == ClpSimplex::basic) {
    1098                     if (length[j] < smallest) {
    1099                          smallest = length[j];
    1100                          iBasic = j;
    1101                     }
    1102                     numberBasic++;
    1103                }
    1104                value += columnSolution[j];
    1105                j = next_[j];
    1106           }
    1107           bool done = false;
    1108           if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
    1109                if (getStatus(iSet) == ClpSimplex::basic)
    1110                     iBasic = iSet + numberColumns; // slack key - use
    1111                done = true;
    1112           } else if (numberBasic == 1) {
    1113                // see if can be key
    1114                double thisSolution = columnSolution[iBasic];
    1115                if (thisSolution < 0.0) {
    1116                     value -= thisSolution;
    1117                     thisSolution = 0.0;
    1118                     columnSolution[iBasic] = thisSolution;
    1119                }
    1120                // try setting slack to a bound
    1121                assert (upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
    1122                double cost1 = COIN_DBL_MAX;
    1123                int whichBound = -1;
    1124                if (upper_[iSet] < 1.0e20) {
    1125                     // try slack at ub
    1126                     double newBasic = thisSolution + upper_[iSet] - value;
    1127                     if (newBasic >= -tolerance) {
    1128                          // can go
    1129                          whichBound = 1;
    1130                          cost1 = newBasic * cost_[iBasic];
    1131                          // But if exact then may be good solution
    1132                          if (fabs(upper_[iSet] - value) < tolerance)
    1133                               cost1 = -COIN_DBL_MAX;
    1134                     }
    1135                }
    1136                if (lower_[iSet] > -1.0e20) {
    1137                     // try slack at lb
    1138                     double newBasic = thisSolution + lower_[iSet] - value;
    1139                     if (newBasic >= -tolerance) {
    1140                          // can go but is it cheaper
    1141                          double cost2 = newBasic * cost_[iBasic];
    1142                          // But if exact then may be good solution
    1143                          if (fabs(lower_[iSet] - value) < tolerance)
    1144                               cost2 = -COIN_DBL_MAX;
    1145                          if (cost2 < cost1)
    1146                               whichBound = 0;
    1147                     }
    1148                }
    1149                if (whichBound != -1) {
    1150                     // key
    1151                     done = true;
    1152                     if (whichBound) {
    1153                          // slack to upper
    1154                          columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
    1155                          setStatus(iSet, ClpSimplex::atUpperBound);
    1156                     } else {
    1157                          // slack to lower
    1158                          columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
    1159                          setStatus(iSet, ClpSimplex::atLowerBound);
    1160                     }
    1161                }
    1162           }
    1163           if (!done) {
    1164                if (!cheapest) {
    1165                     // see if slack can be key
    1166                     if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
    1167                          done = true;
    1168                          setStatus(iSet, ClpSimplex::basic);
    1169                          iBasic = iSet + numberColumns;
    1170                     }
    1171                }
    1172                if (!done) {
    1173                     // set non basic if there was one
    1174                     if (iBasic >= 0)
    1175                          model->setStatus(iBasic, ClpSimplex::atLowerBound);
    1176                     // find cheapest
    1177                     int numberInSet = iEnd - iStart;
    1178                     if (!lowerColumn_) {
    1179                          CoinZeroN(lower, numberInSet);
    1180                     } else {
    1181                          for (int j = 0; j < numberInSet; j++)
    1182                               lower[j] = lowerColumn_[j+iStart];
    1183                     }
    1184                     if (!upperColumn_) {
    1185                          CoinFillN(upper, numberInSet, COIN_DBL_MAX);
    1186                     } else {
    1187                          for (int j = 0; j < numberInSet; j++)
    1188                               upper[j] = upperColumn_[j+iStart];
    1189                     }
    1190                     CoinFillN(solution, numberInSet, 0.0);
    1191                     // and slack
    1192                     iBasic = numberInSet;
    1193                     solution[iBasic] = -value;
    1194                     lower[iBasic] = -upper_[iSet];
    1195                     upper[iBasic] = -lower_[iSet];
    1196                     int kphase;
    1197                     if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
    1198                          // feasible
    1199                          kphase = 1;
    1200                          cost[iBasic] = 0.0;
    1201                          for (int j = 0; j < numberInSet; j++)
    1202                               cost[j] = cost_[j+iStart];
    1203                     } else {
    1204                          // infeasible
    1205                          kphase = 0;
    1206                          // remember bounds are flipped so opposite to natural
    1207                          if (value < lower_[iSet] - tolerance)
    1208                               cost[iBasic] = 1.0;
    1209                          else
    1210                               cost[iBasic] = -1.0;
    1211                          CoinZeroN(cost, numberInSet);
    1212                     }
    1213                     double dualTolerance = model->dualTolerance();
    1214                     for (int iphase = kphase; iphase < 2; iphase++) {
    1215                          if (iphase) {
    1216                               cost[numberInSet] = 0.0;
    1217                               for (int j = 0; j < numberInSet; j++)
    1218                                    cost[j] = cost_[j+iStart];
    1219                          }
    1220                          // now do one row lp
    1221                          bool improve = true;
    1222                          while (improve) {
    1223                               improve = false;
    1224                               double dual = cost[iBasic];
    1225                               int chosen = -1;
    1226                               double best = dualTolerance;
    1227                               int way = 0;
    1228                               for (int i = 0; i <= numberInSet; i++) {
    1229                                    double dj = cost[i] - dual;
    1230                                    double improvement = 0.0;
    1231                                    if (iphase || i < numberInSet)
    1232                                         assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
    1233                                    if (dj > dualTolerance)
    1234                                         improvement = dj * (solution[i] - lower[i]);
    1235                                    else if (dj < -dualTolerance)
    1236                                         improvement = dj * (solution[i] - upper[i]);
    1237                                    if (improvement > best) {
    1238                                         best = improvement;
    1239                                         chosen = i;
    1240                                         if (dj < 0.0) {
    1241                                              way = 1;
    1242                                         } else {
    1243                                              way = -1;
    1244                                         }
    1245                                    }
    1246                               }
    1247                               if (chosen >= 0) {
    1248                                    improve = true;
    1249                                    // now see how far
    1250                                    if (way > 0) {
    1251                                         // incoming increasing so basic decreasing
    1252                                         // if phase 0 then go to nearest bound
    1253                                         double distance = upper[chosen] - solution[chosen];
    1254                                         double basicDistance;
    1255                                         if (!iphase) {
    1256                                              assert (iBasic == numberInSet);
    1257                                              assert (solution[iBasic] > upper[iBasic]);
    1258                                              basicDistance = solution[iBasic] - upper[iBasic];
    1259                                         } else {
    1260                                              basicDistance = solution[iBasic] - lower[iBasic];
    1261                                         }
    1262                                         // need extra coding for unbounded
    1263                                         assert (CoinMin(distance, basicDistance) < 1.0e20);
    1264                                         if (distance > basicDistance) {
    1265                                              // incoming becomes basic
    1266                                              solution[chosen] += basicDistance;
    1267                                              if (!iphase)
    1268                                                   solution[iBasic] = upper[iBasic];
    1269                                              else
    1270                                                   solution[iBasic] = lower[iBasic];
    1271                                              iBasic = chosen;
    1272                                         } else {
    1273                                              // flip
    1274                                              solution[chosen] = upper[chosen];
    1275                                              solution[iBasic] -= distance;
    1276                                         }
    1277                                    } else {
    1278                                         // incoming decreasing so basic increasing
    1279                                         // if phase 0 then go to nearest bound
    1280                                         double distance = solution[chosen] - lower[chosen];
    1281                                         double basicDistance;
    1282                                         if (!iphase) {
    1283                                              assert (iBasic == numberInSet);
    1284                                              assert (solution[iBasic] < lower[iBasic]);
    1285                                              basicDistance = lower[iBasic] - solution[iBasic];
    1286                                         } else {
    1287                                              basicDistance = upper[iBasic] - solution[iBasic];
    1288                                         }
    1289                                         // need extra coding for unbounded - for now just exit
    1290                                         if (CoinMin(distance, basicDistance) > 1.0e20) {
    1291                                              printf("unbounded on set %d\n", iSet);
    1292                                              iphase = 1;
    1293                                              iBasic = numberInSet;
    1294                                              break;
    1295                                         }
    1296                                         if (distance > basicDistance) {
    1297                                              // incoming becomes basic
    1298                                              solution[chosen] -= basicDistance;
    1299                                              if (!iphase)
    1300                                                   solution[iBasic] = lower[iBasic];
    1301                                              else
    1302                                                   solution[iBasic] = upper[iBasic];
    1303                                              iBasic = chosen;
    1304                                         } else {
    1305                                              // flip
    1306                                              solution[chosen] = lower[chosen];
    1307                                              solution[iBasic] += distance;
    1308                                         }
    1309                                    }
    1310                                    if (!iphase) {
    1311                                         if(iBasic < numberInSet)
    1312                                              break; // feasible
    1313                                         else if (solution[iBasic] >= lower[iBasic] &&
    1314                                                   solution[iBasic] <= upper[iBasic])
    1315                                              break; // feasible (on flip)
    1316                                    }
    1317                               }
    1318                          }
    1319                     }
    1320                     // do solution i.e. bounds
    1321                     if (lowerColumn_ || upperColumn_) {
    1322                          for (int j = 0; j < numberInSet; j++) {
    1323                               if (j != iBasic) {
    1324                                    objectiveOffset += solution[j] * cost[j];
    1325                                    if (lowerColumn_ && upperColumn_) {
    1326                                         if (fabs(solution[j] - lowerColumn_[j+iStart]) >
    1327                                                   fabs(solution[j] - upperColumn_[j+iStart]))
    1328                                              setDynamicStatus(j + iStart, atUpperBound);
    1329                                    } else if (upperColumn_ && solution[j] > 0.0) {
    1330                                         setDynamicStatus(j + iStart, atUpperBound);
    1331                                    } else {
    1332                                         setDynamicStatus(j + iStart, atLowerBound);
    1333                                    }
    1334                               }
    1335                          }
    1336                     }
    1337                     // convert iBasic back and do bounds
    1338                     if (iBasic == numberInSet) {
    1339                          // slack basic
    1340                          setStatus(iSet, ClpSimplex::basic);
    1341                          iBasic = iSet + numberColumns;
    1342                     } else {
    1343                          iBasic += fullStart_[iSet];
    1344                          if (back[iBasic] >= 0) {
    1345                               // exists
    1346                               iBasic = back[iBasic];
    1347                          } else {
    1348                               // create
    1349                               CoinBigIndex numberElements = startColumn[firstAvailable_];
    1350                               int numberThis = startColumn_[iBasic+1] - startColumn_[iBasic];
    1351                               if (numberElements + numberThis > numberElements_) {
    1352                                    // need to redo
    1353                                    numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
    1354                                    matrix_->reserve(numberColumns, numberElements_);
    1355                                    element =  matrix_->getMutableElements();
    1356                                    row = matrix_->getMutableIndices();
    1357                                    // these probably okay but be safe
    1358                                    startColumn = matrix_->getMutableVectorStarts();
    1359                                    length = matrix_->getMutableVectorLengths();
    1360                               }
    1361                               length[firstAvailable_] = numberThis;
    1362                               model->costRegion()[firstAvailable_] = cost_[iBasic];
    1363                               if (lowerColumn_)
    1364                                    model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
    1365                               else
    1366                                    model->lowerRegion()[firstAvailable_] = 0.0;
    1367                               if (upperColumn_)
    1368                                    model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
    1369                               else
    1370                                    model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
    1371                               columnSolution[firstAvailable_] = solution[iBasic-fullStart_[iSet]];
    1372                               CoinBigIndex base = startColumn_[iBasic];
    1373                               for (int j = 0; j < numberThis; j++) {
    1374                                    row[numberElements] = row_[base+j];
    1375                                    element[numberElements++] = element_[base+j];
    1376                               }
    1377                               // already set startColumn[firstAvailable_]=numberElements;
    1378                               id_[firstAvailable_-firstDynamic_] = iBasic;
    1379                               setDynamicStatus(iBasic, inSmall);
    1380                               backward_[firstAvailable_] = iSet;
    1381                               iBasic = firstAvailable_;
    1382                               firstAvailable_++;
    1383                               startColumn[firstAvailable_] = numberElements;
    1384                          }
    1385                          model->setStatus(iBasic, ClpSimplex::basic);
    1386                          // remember bounds flipped
    1387                          if (upper[numberInSet] == lower[numberInSet])
    1388                               setStatus(iSet, ClpSimplex::isFixed);
    1389                          else if (solution[numberInSet] == upper[numberInSet])
    1390                               setStatus(iSet, ClpSimplex::atLowerBound);
    1391                          else if (solution[numberInSet] == lower[numberInSet])
    1392                               setStatus(iSet, ClpSimplex::atUpperBound);
    1393                          else
    1394                               abort();
    1395                     }
    1396                     for (j = iStart; j < iEnd; j++) {
    1397                          int iBack = back[j];
    1398                          if (iBack >= 0) {
    1399                               if (model->getStatus(iBack) != ClpSimplex::basic) {
    1400                                    int inSet = j - iStart;
    1401                                    columnSolution[iBack] = solution[inSet];
    1402                                    if (upper[inSet] == lower[inSet])
    1403                                         model->setStatus(iBack, ClpSimplex::isFixed);
    1404                                    else if (solution[inSet] == upper[inSet])
    1405                                         model->setStatus(iBack, ClpSimplex::atUpperBound);
    1406                                    else if (solution[inSet] == lower[inSet])
    1407                                         model->setStatus(iBack, ClpSimplex::atLowerBound);
    1408                               }
    1409                          }
    1410                     }
    1411                }
    1412           }
    1413           keyVariable_[iSet] = iBasic;
    1414      }
    1415      model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
    1416      delete [] lower;
    1417      delete [] solution;
    1418      delete [] upper;
    1419      delete [] cost;
    1420      // make sure matrix is in good shape
    1421      matrix_->orderMatrix();
    1422      // create effective rhs
    1423      delete [] rhsOffset_;
    1424      rhsOffset_ = new double[numberRows];
    1425      // and redo chains
    1426      memset(mark, 0, numberColumns);
    1427      for (int iColumnX = 0; iColumnX < firstAvailable_; iColumnX++)
    1428           next_[iColumnX] = COIN_INT_MAX;
    1429      for (i = 0; i < numberSets_; i++) {
    1430           keys[i] = COIN_INT_MAX;
    1431           int iKey = keyVariable_[i];
    1432           if (iKey < numberColumns)
    1433                model->setStatus(iKey, ClpSimplex::basic);
    1434      }
    1435      // set up chains
    1436      for (i = 0; i < firstAvailable_; i++) {
    1437           if (model->getStatus(i) == ClpSimplex::basic)
    1438                mark[i] = 1;
    1439           int iSet = backward_[i];
    1440           if (iSet >= 0) {
    1441                int iNext = keys[iSet];
    1442                next_[i] = iNext;
    1443                keys[iSet] = i;
    1444           }
    1445      }
    1446      for (i = 0; i < numberSets_; i++) {
    1447           if (keys[i] != COIN_INT_MAX) {
    1448                // something in set
    1449                int j;
    1450                if (getStatus(i) != ClpSimplex::basic) {
    1451                     // make sure fixed if it is
    1452                     if (upper_[i] == lower_[i])
    1453                          setStatus(i, ClpSimplex::isFixed);
    1454                     // slack not key - choose one with smallest length
    1455                     int smallest = numberRows + 1;
    1456                     int key = -1;
    1457                     j = keys[i];
    1458                     while (1) {
    1459                          if (mark[j] && length[j] < smallest) {
    1460                               key = j;
    1461                               smallest = length[j];
    1462                          }
    1463                          if (next_[j] != COIN_INT_MAX) {
    1464                               j = next_[j];
    1465                          } else {
    1466                               // correct end
    1467                               next_[j] = -(keys[i] + 1);
    1468                               break;
    1469                          }
    1470                     }
    1471                     if (key >= 0) {
    1472                          keyVariable_[i] = key;
    1473                     } else {
    1474                          // nothing basic - make slack key
    1475                          //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
    1476                          // fudge to avoid const problem
    1477                          status_[i] = 1;
    1478                     }
    1479                } else {
    1480                     // slack key
    1481                     keyVariable_[i] = numberColumns + i;
    1482                     int j;
    1483                     double sol = 0.0;
    1484                     j = keys[i];
    1485                     while (1) {
    1486                          sol += columnSolution[j];
    1487                          if (next_[j] != COIN_INT_MAX) {
    1488                               j = next_[j];
    1489                          } else {
    1490                               // correct end
    1491                               next_[j] = -(keys[i] + 1);
    1492                               break;
    1493                          }
    1494                     }
    1495                     if (sol > upper_[i] + tolerance) {
    1496                          setAbove(i);
    1497                     } else if (sol < lower_[i] - tolerance) {
    1498                          setBelow(i);
    1499                     } else {
    1500                          setFeasible(i);
    1501                     }
    1502                }
    1503                // Create next_
    1504                int key = keyVariable_[i];
    1505                redoSet(model, key, keys[i], i);
     1334            // create
     1335            CoinBigIndex numberElements = startColumn[firstAvailable_];
     1336            int numberThis = startColumn_[iBasic + 1] - startColumn_[iBasic];
     1337            if (numberElements + numberThis > numberElements_) {
     1338              // need to redo
     1339              numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
     1340              matrix_->reserve(numberColumns, numberElements_);
     1341              element = matrix_->getMutableElements();
     1342              row = matrix_->getMutableIndices();
     1343              // these probably okay but be safe
     1344              startColumn = matrix_->getMutableVectorStarts();
     1345              length = matrix_->getMutableVectorLengths();
     1346            }
     1347            length[firstAvailable_] = numberThis;
     1348            model->costRegion()[firstAvailable_] = cost_[iBasic];
     1349            if (lowerColumn_)
     1350              model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
     1351            else
     1352              model->lowerRegion()[firstAvailable_] = 0.0;
     1353            if (upperColumn_)
     1354              model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
     1355            else
     1356              model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
     1357            columnSolution[firstAvailable_] = solution[iBasic - fullStart_[iSet]];
     1358            CoinBigIndex base = startColumn_[iBasic];
     1359            for (int j = 0; j < numberThis; j++) {
     1360              row[numberElements] = row_[base + j];
     1361              element[numberElements++] = element_[base + j];
     1362            }
     1363            // already set startColumn[firstAvailable_]=numberElements;
     1364            id_[firstAvailable_ - firstDynamic_] = iBasic;
     1365            setDynamicStatus(iBasic, inSmall);
     1366            backward_[firstAvailable_] = iSet;
     1367            iBasic = firstAvailable_;
     1368            firstAvailable_++;
     1369            startColumn[firstAvailable_] = numberElements;
     1370          }
     1371          model->setStatus(iBasic, ClpSimplex::basic);
     1372          // remember bounds flipped
     1373          if (upper[numberInSet] == lower[numberInSet])
     1374            setStatus(iSet, ClpSimplex::isFixed);
     1375          else if (solution[numberInSet] == upper[numberInSet])
     1376            setStatus(iSet, ClpSimplex::atLowerBound);
     1377          else if (solution[numberInSet] == lower[numberInSet])
     1378            setStatus(iSet, ClpSimplex::atUpperBound);
     1379          else
     1380            abort();
     1381        }
     1382        for (j = iStart; j < iEnd; j++) {
     1383          int iBack = back[j];
     1384          if (iBack >= 0) {
     1385            if (model->getStatus(iBack) != ClpSimplex::basic) {
     1386              int inSet = j - iStart;
     1387              columnSolution[iBack] = solution[inSet];
     1388              if (upper[inSet] == lower[inSet])
     1389                model->setStatus(iBack, ClpSimplex::isFixed);
     1390              else if (solution[inSet] == upper[inSet])
     1391                model->setStatus(iBack, ClpSimplex::atUpperBound);
     1392              else if (solution[inSet] == lower[inSet])
     1393                model->setStatus(iBack, ClpSimplex::atLowerBound);
     1394            }
     1395          }
     1396        }
     1397      }
     1398    }
     1399    keyVariable_[iSet] = iBasic;
     1400  }
     1401  model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
     1402  delete[] lower;
     1403  delete[] solution;
     1404  delete[] upper;
     1405  delete[] cost;
     1406  // make sure matrix is in good shape
     1407  matrix_->orderMatrix();
     1408  // create effective rhs
     1409  delete[] rhsOffset_;
     1410  rhsOffset_ = new double[numberRows];
     1411  // and redo chains
     1412  memset(mark, 0, numberColumns);
     1413  for (int iColumnX = 0; iColumnX < firstAvailable_; iColumnX++)
     1414    next_[iColumnX] = COIN_INT_MAX;
     1415  for (i = 0; i < numberSets_; i++) {
     1416    keys[i] = COIN_INT_MAX;
     1417    int iKey = keyVariable_[i];
     1418    if (iKey < numberColumns)
     1419      model->setStatus(iKey, ClpSimplex::basic);
     1420  }
     1421  // set up chains
     1422  for (i = 0; i < firstAvailable_; i++) {
     1423    if (model->getStatus(i) == ClpSimplex::basic)
     1424      mark[i] = 1;
     1425    int iSet = backward_[i];
     1426    if (iSet >= 0) {
     1427      int iNext = keys[iSet];
     1428      next_[i] = iNext;
     1429      keys[iSet] = i;
     1430    }
     1431  }
     1432  for (i = 0; i < numberSets_; i++) {
     1433    if (keys[i] != COIN_INT_MAX) {
     1434      // something in set
     1435      int j;
     1436      if (getStatus(i) != ClpSimplex::basic) {
     1437        // make sure fixed if it is
     1438        if (upper_[i] == lower_[i])
     1439          setStatus(i, ClpSimplex::isFixed);
     1440        // slack not key - choose one with smallest length
     1441        int smallest = numberRows + 1;
     1442        int key = -1;
     1443        j = keys[i];
     1444        while (1) {
     1445          if (mark[j] && length[j] < smallest) {
     1446            key = j;
     1447            smallest = length[j];
     1448          }
     1449          if (next_[j] != COIN_INT_MAX) {
     1450            j = next_[j];
    15061451          } else {
    1507                // nothing in set!
    1508                next_[i+numberColumns] = -(i + numberColumns + 1);
    1509                keyVariable_[i] = numberColumns + i;
    1510                double sol = 0.0;
    1511                if (sol > upper_[i] + tolerance) {
    1512                     setAbove(i);
    1513                } else if (sol < lower_[i] - tolerance) {
    1514                     setBelow(i);
    1515                } else {
    1516                     setFeasible(i);
    1517                }
    1518           }
    1519      }
    1520      delete [] keys;
    1521      delete [] mark;
    1522      delete [] back;
    1523      rhsOffset(model, true);
     1452            // correct end
     1453            next_[j] = -(keys[i] + 1);
     1454            break;
     1455          }
     1456        }
     1457        if (key >= 0) {
     1458          keyVariable_[i] = key;
     1459        } else {
     1460          // nothing basic - make slack key
     1461          //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
     1462          // fudge to avoid const problem
     1463          status_[i] = 1;
     1464        }
     1465      } else {
     1466        // slack key
     1467        keyVariable_[i] = numberColumns + i;
     1468        int j;
     1469        double sol = 0.0;
     1470        j = keys[i];
     1471        while (1) {
     1472          sol += columnSolution[j];
     1473          if (next_[j] != COIN_INT_MAX) {
     1474            j = next_[j];
     1475          } else {
     1476            // correct end
     1477            next_[j] = -(keys[i] + 1);
     1478            break;
     1479          }
     1480        }
     1481        if (sol > upper_[i] + tolerance) {
     1482          setAbove(i);
     1483        } else if (sol < lower_[i] - tolerance) {
     1484          setBelow(i);
     1485        } else {
     1486          setFeasible(i);
     1487        }
     1488      }
     1489      // Create next_
     1490      int key = keyVariable_[i];
     1491      redoSet(model, key, keys[i], i);
     1492    } else {
     1493      // nothing in set!
     1494      next_[i + numberColumns] = -(i + numberColumns + 1);
     1495      keyVariable_[i] = numberColumns + i;
     1496      double sol = 0.0;
     1497      if (sol > upper_[i] + tolerance) {
     1498        setAbove(i);
     1499      } else if (sol < lower_[i] - tolerance) {
     1500        setBelow(i);
     1501      } else {
     1502        setFeasible(i);
     1503      }
     1504    }
     1505  }
     1506  delete[] keys;
     1507  delete[] mark;
     1508  delete[] back;
     1509  rhsOffset(model, true);
    15241510}
    15251511/* Returns effective RHS if it is being used.  This is used for long problems
     
    15271513   expensive.  This may re-compute */
    15281514double *
    1529 ClpGubDynamicMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh,
    1530                                bool
     1515ClpGubDynamicMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh,
     1516  bool
    15311517#ifdef CLP_DEBUG
    1532                                check
     1518    check
    15331519#endif
    1534                               )
     1520)
    15351521{
    1536      //forceRefresh=true;
    1537      //check=false;
     1522  //forceRefresh=true;
     1523  //check=false;
    15381524#ifdef CLP_DEBUG
    1539      double * saveE = NULL;
    1540      if (rhsOffset_ && check) {
    1541           int numberRows = model->numberRows();
    1542           saveE = new double[numberRows];
    1543      }
     1525  double *saveE = NULL;
     1526  if (rhsOffset_ && check) {
     1527    int numberRows = model->numberRows();
     1528    saveE = new double[numberRows];
     1529  }
    15441530#endif
    1545      if (rhsOffset_) {
     1531  if (rhsOffset_) {
    15461532#ifdef CLP_DEBUG
    1547           if (check) {
    1548                // no need - but check anyway
    1549                int numberRows = model->numberRows();
    1550                double * rhs = new double[numberRows];
    1551                int numberColumns = model->numberColumns();
    1552                int iRow;
    1553                CoinZeroN(rhs, numberRows);
    1554                // do ones at bounds before gub
    1555                const double * smallSolution = model->solutionRegion();
    1556                const double * element = matrix_->getElements();
    1557                const int * row = matrix_->getIndices();
    1558                const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    1559                const int * length = matrix_->getVectorLengths();
    1560                int iColumn;
    1561                for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
    1562                     if (model->getStatus(iColumn) != ClpSimplex::basic) {
    1563                          double value = smallSolution[iColumn];
    1564                          for (CoinBigIndex j = startColumn[iColumn];
    1565                                    j < startColumn[iColumn] + length[iColumn]; j++) {
    1566                               int jRow = row[j];
    1567                               rhs[jRow] -= value * element[j];
    1568                          }
    1569                     }
    1570                }
    1571                if (lowerColumn_ || upperColumn_) {
    1572                     double * solution = new double [numberGubColumns_];
    1573                     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
    1574                          double value = 0.0;
    1575                          if(getDynamicStatus(iColumn) == atUpperBound)
    1576                               value = upperColumn_[iColumn];
    1577                          else if (lowerColumn_)
    1578                               value = lowerColumn_[iColumn];
    1579                          solution[iColumn] = value;
    1580                     }
    1581                     // ones at bounds in small and gub
    1582                     for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
    1583                          int jFull = id_[iColumn-firstDynamic_];
    1584                          solution[jFull] = smallSolution[iColumn];
    1585                     }
    1586                     // zero all basic in small model
    1587                     int * pivotVariable = model->pivotVariable();
    1588                     for (iRow = 0; iRow < numberRows; iRow++) {
    1589                          int iColumn = pivotVariable[iRow];
    1590                          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
    1591                               int iSequence = id_[iColumn-firstDynamic_];
    1592                               solution[iSequence] = 0.0;
    1593                          }
    1594                     }
    1595                     // and now compute value to use for key
    1596                     ClpSimplex::Status iStatus;
    1597                     for (int iSet = 0; iSet < numberSets_; iSet++) {
    1598                          iColumn = keyVariable_[iSet];
    1599                          if (iColumn < numberColumns) {
    1600                               int iSequence = id_[iColumn-firstDynamic_];
    1601                               solution[iSequence] = 0.0;
    1602                               double b = 0.0;
    1603                               // key is structural - where is slack
    1604                               iStatus = getStatus(iSet);
    1605                               assert (iStatus != ClpSimplex::basic);
    1606                               if (iStatus == ClpSimplex::atLowerBound)
    1607                                    b = lowerSet_[iSet];
    1608                               else
    1609                                    b = upperSet_[iSet];
    1610                               // subtract out others at bounds
    1611                               for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
    1612                                    b -= solution[j];
    1613                               solution[iSequence] = b;
    1614                          }
    1615                     }
    1616                     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
    1617                          double value = solution[iColumn];
    1618                          if (value) {
    1619                               for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
    1620                                    int iRow = row_[j];
    1621                                    rhs[iRow] -= element_[j] * value;
    1622                               }
    1623                          }
    1624                     }
    1625                     // now do lower and upper bounds on sets
    1626                     for (int iSet = 0; iSet < numberSets_; iSet++) {
    1627                          iColumn = keyVariable_[iSet];
    1628                          double shift = 0.0;
    1629                          for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++) {
    1630                               if (getDynamicStatus(j) != inSmall && j != iColumn) {
    1631                                    if (getDynamicStatus(j) == atLowerBound) {
    1632                                         if (lowerColumn_)
    1633                                              shift += lowerColumn_[j];
    1634                                    } else {
    1635                                         shift += upperColumn_[j];
    1636                                    }
    1637                               }
    1638                          }
    1639                          if (lowerSet_[iSet] > -1.0e20)
    1640                               assert(fabs(lower_[iSet] - (lowerSet_[iSet] - shift)) < 1.0e-3);
    1641                          if (upperSet_[iSet] < 1.0e20)
    1642                               assert(fabs(upper_[iSet] - ( upperSet_[iSet] - shift)) < 1.0e-3);
    1643                     }
    1644                     delete [] solution;
    1645                } else {
    1646                     // no bounds
    1647                     ClpSimplex::Status iStatus;
    1648                     for (int iSet = 0; iSet < numberSets_; iSet++) {
    1649                          int iColumn = keyVariable_[iSet];
    1650                          if (iColumn < numberColumns) {
    1651                               int iSequence = id_[iColumn-firstDynamic_];
    1652                               double b = 0.0;
    1653                               // key is structural - where is slack
    1654                               iStatus = getStatus(iSet);
    1655                               assert (iStatus != ClpSimplex::basic);
    1656                               if (iStatus == ClpSimplex::atLowerBound)
    1657                                    b = lower_[iSet];
    1658                               else
    1659                                    b = upper_[iSet];
    1660                               if (b) {
    1661                                    for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) {
    1662                                         int iRow = row_[j];
    1663                                         rhs[iRow] -= element_[j] * b;
    1664                                    }
    1665                               }
    1666                          }
    1667                     }
    1668                }
    1669                for (iRow = 0; iRow < numberRows; iRow++) {
    1670                     if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
    1671                          printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
    1672                }
    1673                CoinMemcpyN(rhs, numberRows, saveE);
    1674                delete [] rhs;
    1675           }
     1533    if (check) {
     1534      // no need - but check anyway
     1535      int numberRows = model->numberRows();
     1536      double *rhs = new double[numberRows];
     1537      int numberColumns = model->numberColumns();
     1538      int iRow;
     1539      CoinZeroN(rhs, numberRows);
     1540      // do ones at bounds before gub
     1541      const double *smallSolution = model->solutionRegion();
     1542      const double *element = matrix_->getElements();
     1543      const int *row = matrix_->getIndices();
     1544      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
     1545      const int *length = matrix_->getVectorLengths();
     1546      int iColumn;
     1547      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
     1548        if (model->getStatus(iColumn) != ClpSimplex::basic) {
     1549          double value = smallSolution[iColumn];
     1550          for (CoinBigIndex j = startColumn[iColumn];
     1551               j < startColumn[iColumn] + length[iColumn]; j++) {
     1552            int jRow = row[j];
     1553            rhs[jRow] -= value * element[j];
     1554          }
     1555        }
     1556      }
     1557      if (lowerColumn_ || upperColumn_) {
     1558        double *solution = new double[numberGubColumns_];
     1559        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
     1560          double value = 0.0;
     1561          if (getDynamicStatus(iColumn) == atUpperBound)
     1562            value = upperColumn_[iColumn];
     1563          else if (lowerColumn_)
     1564            value = lowerColumn_[iColumn];
     1565          solution[iColumn] = value;
     1566        }
     1567        // ones at bounds in small and gub
     1568        for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
     1569          int jFull = id_[iColumn - firstDynamic_];
     1570          solution[jFull] = smallSolution[iColumn];
     1571        }
     1572        // zero all basic in small model
     1573        int *pivotVariable = model->pivotVariable();
     1574        for (iRow = 0; iRow < numberRows; iRow++) {
     1575          int iColumn = pivotVariable[iRow];
     1576          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
     1577            int iSequence = id_[iColumn - firstDynamic_];
     1578            solution[iSequence] = 0.0;
     1579          }
     1580        }
     1581        // and now compute value to use for key
     1582        ClpSimplex::Status iStatus;
     1583        for (int iSet = 0; iSet < numberSets_; iSet++) {
     1584          iColumn = keyVariable_[iSet];
     1585          if (iColumn < numberColumns) {
     1586            int iSequence = id_[iColumn - firstDynamic_];
     1587            solution[iSequence] = 0.0;
     1588            double b = 0.0;
     1589            // key is structural - where is slack
     1590            iStatus = getStatus(iSet);
     1591            assert(iStatus != ClpSimplex::basic);
     1592            if (iStatus == ClpSimplex::atLowerBound)
     1593              b = lowerSet_[iSet];
     1594            else
     1595              b = upperSet_[iSet];
     1596            // subtract out others at bounds
     1597            for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
     1598              b -= solution[j];
     1599            solution[iSequence] = b;
     1600          }
     1601        }
     1602        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
     1603          double value = solution[iColumn];
     1604          if (value) {
     1605            for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
     1606              int iRow = row_[j];
     1607              rhs[iRow] -= element_[j] * value;
     1608            }
     1609          }
     1610        }
     1611        // now do lower and upper bounds on sets
     1612        for (int iSet = 0; iSet < numberSets_; iSet++) {
     1613          iColumn = keyVariable_[iSet];
     1614          double shift = 0.0;
     1615          for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++) {
     1616            if (getDynamicStatus(j) != inSmall && j != iColumn) {
     1617              if (getDynamicStatus(j) == atLowerBound) {
     1618                if (lowerColumn_)
     1619                  shift += lowerColumn_[j];
     1620              } else {
     1621                shift += upperColumn_[j];
     1622              }
     1623            }
     1624          }
     1625          if (lowerSet_[iSet] > -1.0e20)
     1626            assert(fabs(lower_[iSet] - (lowerSet_[iSet] - shift)) < 1.0e-3);
     1627          if (upperSet_[iSet] < 1.0e20)
     1628            assert(fabs(upper_[iSet] - (upperSet_[iSet] - shift)) < 1.0e-3);
     1629        }
     1630        delete[] solution;
     1631      } else {
     1632        // no bounds
     1633        ClpSimplex::Status iStatus;
     1634        for (int iSet = 0; iSet < numberSets_; iSet++) {
     1635          int iColumn = keyVariable_[iSet];
     1636          if (iColumn < numberColumns) {
     1637            int iSequence = id_[iColumn - firstDynamic_];
     1638            double b = 0.0;
     1639            // key is structural - where is slack
     1640            iStatus = getStatus(iSet);
     1641            assert(iStatus != ClpSimplex::basic);
     1642            if (iStatus == ClpSimplex::atLowerBound)
     1643              b = lower_[iSet];
     1644            else
     1645              b = upper_[iSet];
     1646            if (b) {
     1647              for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence + 1]; j++) {
     1648                int iRow = row_[j];
     1649                rhs[iRow] -= element_[j] * b;
     1650              }
     1651            }
     1652          }
     1653        }
     1654      }
     1655      for (iRow = 0; iRow < numberRows; iRow++) {
     1656        if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
     1657          printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
     1658      }
     1659      CoinMemcpyN(rhs, numberRows, saveE);
     1660      delete[] rhs;
     1661    }
    16761662#endif
    1677           if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
    1678                                lastRefresh_ + refreshFrequency_)) {
    1679                int numberRows = model->numberRows();
    1680                int numberColumns = model->numberColumns();
    1681                int iRow;
    1682                CoinZeroN(rhsOffset_, numberRows);
    1683                // do ones at bounds before gub
    1684                const double * smallSolution = model->solutionRegion();
    1685                const double * element = matrix_->getElements();
    1686                const int * row = matrix_->getIndices();
    1687                const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    1688                const int * length = matrix_->getVectorLengths();
    1689                int iColumn;
    1690                for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
    1691                     if (model->getStatus(iColumn) != ClpSimplex::basic) {
    1692                          double value = smallSolution[iColumn];
    1693                          for (CoinBigIndex j = startColumn[iColumn];
    1694                                    j < startColumn[iColumn] + length[iColumn]; j++) {
    1695                               int jRow = row[j];
    1696                               rhsOffset_[jRow] -= value * element[j];
    1697                          }
    1698                     }
    1699                }
    1700                if (lowerColumn_ || upperColumn_) {
    1701                     double * solution = new double [numberGubColumns_];
    1702                     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
    1703                          double value = 0.0;
    1704                          if(getDynamicStatus(iColumn) == atUpperBound)
    1705                               value = upperColumn_[iColumn];
    1706                          else if (lowerColumn_)
    1707                               value = lowerColumn_[iColumn];
    1708                          solution[iColumn] = value;
    1709                     }
    1710                     // ones in gub and in small problem
    1711                     for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
    1712                          int jFull = id_[iColumn-firstDynamic_];
    1713                          solution[jFull] = smallSolution[iColumn];
    1714                     }
    1715                     // zero all basic in small model
    1716                     int * pivotVariable = model->pivotVariable();
    1717                     for (iRow = 0; iRow < numberRows; iRow++) {
    1718                          int iColumn = pivotVariable[iRow];
    1719                          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
    1720                               int iSequence = id_[iColumn-firstDynamic_];
    1721                               solution[iSequence] = 0.0;
    1722                          }
    1723                     }
    1724                     // and now compute value to use for key
    1725                     ClpSimplex::Status iStatus;
    1726                     int iSet;
    1727                     for ( iSet = 0; iSet < numberSets_; iSet++) {
    1728                          iColumn = keyVariable_[iSet];
    1729                          if (iColumn < numberColumns) {
    1730                               int iSequence = id_[iColumn-firstDynamic_];
    1731                               solution[iSequence] = 0.0;
    1732                               double b = 0.0;
    1733                               // key is structural - where is slack
    1734                               iStatus = getStatus(iSet);
    1735                               assert (iStatus != ClpSimplex::basic);
    1736                               if (iStatus == ClpSimplex::atLowerBound)
    1737                                    b = lowerSet_[iSet];
    1738                               else
    1739                                    b = upperSet_[iSet];
    1740                               // subtract out others at bounds
    1741                               for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
    1742                                    b -= solution[j];
    1743                               solution[iSequence] = b;
    1744                          }
    1745                     }
    1746                     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
    1747                          double value = solution[iColumn];
    1748                          if (value) {
    1749                               for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
    1750                                    int iRow = row_[j];
    1751                                    rhsOffset_[iRow] -= element_[j] * value;
    1752                               }
    1753                          }
    1754                     }
    1755                     // now do lower and upper bounds on sets
    1756                     // and offset
    1757                     double objectiveOffset = 0.0;
    1758                     for ( iSet = 0; iSet < numberSets_; iSet++) {
    1759                          iColumn = keyVariable_[iSet];
    1760                          double shift = 0.0;
    1761                          for (CoinBigIndex j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++) {
    1762                               if (getDynamicStatus(j) != inSmall) {
    1763                                    double value = 0.0;
    1764                                    if (getDynamicStatus(j) == atLowerBound) {
    1765                                         if (lowerColumn_)
    1766                                              value = lowerColumn_[j];
    1767                                    } else {
    1768                                         value = upperColumn_[j];
    1769                                    }
    1770                                    if (j != iColumn)
    1771                                         shift += value;
    1772                                    objectiveOffset += value * cost_[j];
    1773                               }
    1774                          }
    1775                          if (lowerSet_[iSet] > -1.0e20)
    1776                               lower_[iSet] = lowerSet_[iSet] - shift;
    1777                          if (upperSet_[iSet] < 1.0e20)
    1778                               upper_[iSet] = upperSet_[iSet] - shift;
    1779                     }
    1780                     delete [] solution;
    1781                     model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
    1782                } else {
    1783                     // no bounds
    1784                     ClpSimplex::Status iStatus;
    1785                     for (int iSet = 0; iSet < numberSets_; iSet++) {
    1786                          int iColumn = keyVariable_[iSet];
    1787                          if (iColumn < numberColumns) {
    1788                               int iSequence = id_[iColumn-firstDynamic_];
    1789                               double b = 0.0;
    1790                               // key is structural - where is slack
    1791                               iStatus = getStatus(iSet);
    1792                               assert (iStatus != ClpSimplex::basic);
    1793                               if (iStatus == ClpSimplex::atLowerBound)
    1794                                    b = lower_[iSet];
    1795                               else
    1796                                    b = upper_[iSet];
    1797                               if (b) {
    1798                                    for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) {
    1799                                         int iRow = row_[j];
    1800                                         rhsOffset_[iRow] -= element_[j] * b;
    1801                                    }
    1802                               }
    1803                          }
    1804                     }
    1805                }
     1663    if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) {
     1664      int numberRows = model->numberRows();
     1665      int numberColumns = model->numberColumns();
     1666      int iRow;
     1667      CoinZeroN(rhsOffset_, numberRows);
     1668      // do ones at bounds before gub
     1669      const double *smallSolution = model->solutionRegion();
     1670      const double *element = matrix_->getElements();
     1671      const int *row = matrix_->getIndices();
     1672      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
     1673      const int *length = matrix_->getVectorLengths();
     1674      int iColumn;
     1675      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
     1676        if (model->getStatus(iColumn) != ClpSimplex::basic) {
     1677          double value = smallSolution[iColumn];
     1678          for (CoinBigIndex j = startColumn[iColumn];
     1679               j < startColumn[iColumn] + length[iColumn]; j++) {
     1680            int jRow = row[j];
     1681            rhsOffset_[jRow] -= value * element[j];
     1682          }
     1683        }
     1684      }
     1685      if (lowerColumn_ || upperColumn_) {
     1686        double *solution = new double[numberGubColumns_];
     1687        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
     1688          double value = 0.0;
     1689          if (getDynamicStatus(iColumn) == atUpperBound)
     1690            value = upperColumn_[iColumn];
     1691          else if (lowerColumn_)
     1692            value = lowerColumn_[iColumn];
     1693          solution[iColumn] = value;
     1694        }
     1695        // ones in gub and in small problem
     1696        for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
     1697          int jFull = id_[iColumn - firstDynamic_];
     1698          solution[jFull] = smallSolution[iColumn];
     1699        }
     1700        // zero all basic in small model
     1701        int *pivotVariable = model->pivotVariable();
     1702        for (iRow = 0; iRow < numberRows; iRow++) {
     1703          int iColumn = pivotVariable[iRow];
     1704          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
     1705            int iSequence = id_[iColumn - firstDynamic_];
     1706            solution[iSequence] = 0.0;
     1707          }
     1708        }
     1709        // and now compute value to use for key
     1710        ClpSimplex::Status iStatus;
     1711        int iSet;
     1712        for (iSet = 0; iSet < numberSets_; iSet++) {
     1713          iColumn = keyVariable_[iSet];
     1714          if (iColumn < numberColumns) {
     1715            int iSequence = id_[iColumn - firstDynamic_];
     1716            solution[iSequence] = 0.0;
     1717            double b = 0.0;
     1718            // key is structural - where is slack
     1719            iStatus = getStatus(iSet);
     1720            assert(iStatus != ClpSimplex::basic);
     1721            if (iStatus == ClpSimplex::atLowerBound)
     1722              b = lowerSet_[iSet];
     1723            else
     1724              b = upperSet_[iSet];
     1725            // subtract out others at bounds
     1726            for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
     1727              b -= solution[j];
     1728            solution[iSequence] = b;
     1729          }
     1730        }
     1731        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
     1732          double value = solution[iColumn];
     1733          if (value) {
     1734            for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
     1735              int iRow = row_[j];
     1736              rhsOffset_[iRow] -= element_[j] * value;
     1737            }
     1738          }
     1739        }
     1740        // now do lower and upper bounds on sets
     1741        // and offset
     1742        double objectiveOffset = 0.0;
     1743        for (iSet = 0; iSet < numberSets_; iSet++) {
     1744          iColumn = keyVariable_[iSet];
     1745          double shift = 0.0;
     1746          for (CoinBigIndex j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++) {
     1747            if (getDynamicStatus(j) != inSmall) {
     1748              double value = 0.0;
     1749              if (getDynamicStatus(j) == atLowerBound) {
     1750                if (lowerColumn_)
     1751                  value = lowerColumn_[j];
     1752              } else {
     1753                value = upperColumn_[j];
     1754              }
     1755              if (j != iColumn)
     1756                shift += value;
     1757              objectiveOffset += value * cost_[j];
     1758            }
     1759          }
     1760          if (lowerSet_[iSet] > -1.0e20)
     1761            lower_[iSet] = lowerSet_[iSet] - shift;
     1762          if (upperSet_[iSet] < 1.0e20)
     1763            upper_[iSet] = upperSet_[iSet] - shift;
     1764        }
     1765        delete[] solution;
     1766        model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
     1767      } else {
     1768        // no bounds
     1769        ClpSimplex::Status iStatus;
     1770        for (int iSet = 0; iSet < numberSets_; iSet++) {
     1771          int iColumn = keyVariable_[iSet];
     1772          if (iColumn < numberColumns) {
     1773            int iSequence = id_[iColumn - firstDynamic_];
     1774            double b = 0.0;
     1775            // key is structural - where is slack
     1776            iStatus = getStatus(iSet);
     1777            assert(iStatus != ClpSimplex::basic);
     1778            if (iStatus == ClpSimplex::atLowerBound)
     1779              b = lower_[iSet];
     1780            else
     1781              b = upper_[iSet];
     1782            if (b) {
     1783              for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence + 1]; j++) {
     1784                int iRow = row_[j];
     1785                rhsOffset_[iRow] -= element_[j] * b;
     1786              }
     1787            }
     1788          }
     1789        }
     1790      }
    18061791#ifdef CLP_DEBUG
    1807                if (saveE) {
    1808                     for (iRow = 0; iRow < numberRows; iRow++) {
    1809                          if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
    1810                               printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
    1811                     }
    1812                     delete [] saveE;
    1813                }
     1792      if (saveE) {
     1793        for (iRow = 0; iRow < numberRows; iRow++) {
     1794          if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
     1795            printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
     1796        }
     1797        delete[] saveE;
     1798      }
    18141799#endif
    1815                lastRefresh_ = model->numberIterations();
    1816           }
    1817      }
    1818      return rhsOffset_;
     1800      lastRefresh_ = model->numberIterations();
     1801    }
     1802  }
     1803  return rhsOffset_;
    18191804}
    18201805/*
    18211806  update information for a pivot (and effective rhs)
    18221807*/
    1823 int
    1824 ClpGubDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue)
     1808int ClpGubDynamicMatrix::updatePivot(ClpSimplex *model, double oldInValue, double oldOutValue)
    18251809{
    18261810
    1827      // now update working model
    1828      int sequenceIn = model->sequenceIn();
    1829      int sequenceOut = model->sequenceOut();
    1830      bool doPrinting = (model->messageHandler()->logLevel() == 63);
    1831      bool print = false;
    1832      int iSet;
    1833      int trueIn = -1;
    1834      int trueOut = -1;
    1835      int numberRows = model->numberRows();
    1836      int numberColumns = model->numberColumns();
    1837      if (sequenceIn == firstAvailable_) {
    1838           if (doPrinting)
    1839                printf("New variable ");
    1840           if (sequenceIn != sequenceOut) {
    1841                insertNonBasic(firstAvailable_, backward_[firstAvailable_]);
    1842                setDynamicStatus(id_[sequenceIn-firstDynamic_], inSmall);
    1843                firstAvailable_++;
    1844           } else {
    1845                int bigSequence = id_[sequenceIn-firstDynamic_];
    1846                if (model->getStatus(sequenceIn) == ClpSimplex::atUpperBound)
    1847                     setDynamicStatus(bigSequence, atUpperBound);
    1848                else
    1849                     setDynamicStatus(bigSequence, atLowerBound);
    1850           }
    1851           synchronize(model, 8);
    1852      }
    1853      if (sequenceIn < lastDynamic_) {
    1854           iSet = backward_[sequenceIn];
    1855           if (iSet >= 0) {
    1856                int bigSequence = id_[sequenceIn-firstDynamic_];
    1857                trueIn = bigSequence + numberRows + numberColumns + numberSets_;
    1858                if (doPrinting)
    1859                     printf(" incoming set %d big seq %d", iSet, bigSequence);
    1860                print = true;
    1861           }
    1862      } else if (sequenceIn >= numberRows + numberColumns) {
    1863           trueIn = numberRows + numberColumns + gubSlackIn_;
    1864      }
    1865      if (sequenceOut < lastDynamic_) {
    1866           iSet = backward_[sequenceOut];
    1867           if (iSet >= 0) {
    1868                int bigSequence = id_[sequenceOut-firstDynamic_];
    1869                trueOut = bigSequence + firstDynamic_;
    1870                if (getDynamicStatus(bigSequence) != inSmall) {
    1871                     if (model->getStatus(sequenceOut) == ClpSimplex::atUpperBound)
    1872                          setDynamicStatus(bigSequence, atUpperBound);
    1873                     else
    1874                          setDynamicStatus(bigSequence, atLowerBound);
    1875                }
    1876                if (doPrinting)
    1877                     printf(" ,outgoing set %d big seq %d,", iSet, bigSequence);
    1878                print = true;
    1879                model->setSequenceIn(sequenceOut);
    1880                synchronize(model, 8);
    1881                model->setSequenceIn(sequenceIn);
    1882           }
    1883      }
    1884      if (print && doPrinting)
    1885           printf("\n");
    1886      ClpGubMatrix::updatePivot(model, oldInValue, oldOutValue);
    1887      // Redo true in and out
    1888      if (trueIn >= 0)
    1889           trueSequenceIn_ = trueIn;
    1890      if (trueOut >= 0)
    1891           trueSequenceOut_ = trueOut;
    1892      if (doPrinting && 0) {
    1893           for (int i = 0; i < numberSets_; i++) {
    1894                printf("set %d key %d lower %g upper %g\n", i, keyVariable_[i], lower_[i], upper_[i]);
    1895                for (int j = fullStart_[i]; j < fullStart_[i+1]; j++)
    1896                     if (getDynamicStatus(j) == atUpperBound) {
    1897                          bool print = true;
    1898                          for (int k = firstDynamic_; k < firstAvailable_; k++) {
    1899                               if (id_[k-firstDynamic_] == j)
    1900                                    print = false;
    1901                               if (id_[k-firstDynamic_] == j)
    1902                                    assert(getDynamicStatus(j) == inSmall);
    1903                          }
    1904                          if (print)
    1905                               printf("variable %d at ub\n", j);
    1906                     }
    1907           }
    1908      }
     1811  // now update working model
     1812  int sequenceIn = model->sequenceIn();
     1813  int sequenceOut = model->sequenceOut();
     1814  bool doPrinting = (model->messageHandler()->logLevel() == 63);
     1815  bool print = false;
     1816  int iSet;
     1817  int trueIn = -1;
     1818  int trueOut = -1;
     1819  int numberRows = model->numberRows();
     1820  int numberColumns = model->numberColumns();
     1821  if (sequenceIn == firstAvailable_) {
     1822    if (doPrinting)
     1823      printf("New variable ");
     1824    if (sequenceIn != sequenceOut) {
     1825      insertNonBasic(firstAvailable_, backward_[firstAvailable_]);
     1826      setDynamicStatus(id_[sequenceIn - firstDynamic_], inSmall);
     1827      firstAvailable_++;
     1828    } else {
     1829      int bigSequence = id_[sequenceIn - firstDynamic_];
     1830      if (model->getStatus(sequenceIn) == ClpSimplex::atUpperBound)
     1831        setDynamicStatus(bigSequence, atUpperBound);
     1832      else
     1833        setDynamicStatus(bigSequence, atLowerBound);
     1834    }
     1835    synchronize(model, 8);
     1836  }
     1837  if (sequenceIn < lastDynamic_) {
     1838    iSet = backward_[sequenceIn];
     1839    if (iSet >= 0) {
     1840      int bigSequence = id_[sequenceIn - firstDynamic_];
     1841      trueIn = bigSequence + numberRows + numberColumns + numberSets_;
     1842      if (doPrinting)
     1843        printf(" incoming set %d big seq %d", iSet, bigSequence);
     1844      print = true;
     1845    }
     1846  } else if (sequenceIn >= numberRows + numberColumns) {
     1847    trueIn = numberRows + numberColumns + gubSlackIn_;
     1848  }
     1849  if (sequenceOut < lastDynamic_) {
     1850    iSet = backward_[sequenceOut];
     1851    if (iSet >= 0) {
     1852      int bigSequence = id_[sequenceOut - firstDynamic_];
     1853      trueOut = bigSequence + firstDynamic_;
     1854      if (getDynamicStatus(bigSequence) != inSmall) {
     1855        if (model->getStatus(sequenceOut) == ClpSimplex::atUpperBound)
     1856          setDynamicStatus(bigSequence, atUpperBound);
     1857        else
     1858          setDynamicStatus(bigSequence, atLowerBound);
     1859      }
     1860      if (doPrinting)
     1861        printf(" ,outgoing set %d big seq %d,", iSet, bigSequence);
     1862      print = true;
     1863      model->setSequenceIn(sequenceOut);
     1864      synchronize(model, 8);
     1865      model->setSequenceIn(sequenceIn);
     1866    }
     1867  }
     1868  if (print && doPrinting)
     1869    printf("\n");
     1870  ClpGubMatrix::updatePivot(model, oldInValue, oldOutValue);
     1871  // Redo true in and out
     1872  if (trueIn >= 0)
     1873    trueSequenceIn_ = trueIn;
     1874  if (trueOut >= 0)
     1875    trueSequenceOut_ = trueOut;
     1876  if (doPrinting && 0) {
     1877    for (int i = 0; i < numberSets_; i++) {
     1878      printf("set %d key %d lower %g upper %g\n", i, keyVariable_[i], lower_[i], upper_[i]);
     1879      for (int j = fullStart_[i]; j < fullStart_[i + 1]; j++)
     1880        if (getDynamicStatus(j) == atUpperBound) {
     1881          bool print = true;
     1882          for (int k = firstDynamic_; k < firstAvailable_; k++) {
     1883            if (id_[k - firstDynamic_] == j)
     1884              print = false;
     1885            if (id_[k - firstDynamic_] == j)
     1886              assert(getDynamicStatus(j) == inSmall);
     1887          }
     1888          if (print)
     1889            printf("variable %d at ub\n", j);
     1890        }
     1891    }
     1892  }
    19091893#ifdef CLP_DEBUG
    1910      char * inSmall = new char [numberGubColumns_];
    1911      memset(inSmall, 0, numberGubColumns_);
    1912      for (int i = 0; i < numberGubColumns_; i++)
    1913           if (getDynamicStatus(i) == ClpGubDynamicMatrix::inSmall)
    1914                inSmall[i] = 1;
    1915      for (int i = firstDynamic_; i < firstAvailable_; i++) {
    1916           int k = id_[i-firstDynamic_];
    1917           inSmall[k] = 0;
    1918      }
    1919      for (int i = 0; i < numberGubColumns_; i++)
    1920           assert (!inSmall[i]);
    1921      delete [] inSmall;
     1894  char *inSmall = new char[numberGubColumns_];
     1895  memset(inSmall, 0, numberGubColumns_);
     1896  for (int i = 0; i < numberGubColumns_; i++)
     1897    if (getDynamicStatus(i) == ClpGubDynamicMatrix::inSmall)
     1898      inSmall[i] = 1;
     1899  for (int i = firstDynamic_; i < firstAvailable_; i++) {
     1900    int k = id_[i - firstDynamic_];
     1901    inSmall[k] = 0;
     1902  }
     1903  for (int i = 0; i < numberGubColumns_; i++)
     1904    assert(!inSmall[i]);
     1905  delete[] inSmall;
    19221906#endif
    1923      return 0;
     1907  return 0;
    19241908}
    1925 void
    1926 ClpGubDynamicMatrix::times(double scalar,
    1927                            const double * x, double * y) const
     1909void ClpGubDynamicMatrix::times(double scalar,
     1910  const double *x, double *y) const
    19281911{
    1929      if (model_->specialOptions() != 16) {
    1930           ClpPackedMatrix::times(scalar, x, y);
    1931      } else {
    1932           int iRow;
    1933           int numberColumns = model_->numberColumns();
    1934           int numberRows = model_->numberRows();
    1935           const double * element = matrix_->getElements();
    1936           const int * row = matrix_->getIndices();
    1937           const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    1938           const int * length = matrix_->getVectorLengths();
    1939           int * pivotVariable = model_->pivotVariable();
    1940           int numberToDo = 0;
    1941           for (iRow = 0; iRow < numberRows; iRow++) {
    1942                y[iRow] -= scalar * rhsOffset_[iRow];
    1943                int iColumn = pivotVariable[iRow];
    1944                if (iColumn < numberColumns) {
    1945                     int iSet = backward_[iColumn];
    1946                     if (iSet >= 0 && toIndex_[iSet] < 0) {
    1947                          toIndex_[iSet] = 0;
    1948                          fromIndex_[numberToDo++] = iSet;
    1949                     }
    1950                     CoinBigIndex j;
    1951                     double value = scalar * x[iColumn];
    1952                     if (value) {
    1953                          for (j = startColumn[iColumn];
    1954                                    j < startColumn[iColumn] + length[iColumn]; j++) {
    1955                               int jRow = row[j];
    1956                               y[jRow] += value * element[j];
    1957                          }
    1958                     }
    1959                }
    1960           }
    1961           // and gubs which are interacting
    1962           for (int jSet = 0; jSet < numberToDo; jSet++) {
    1963                int iSet = fromIndex_[jSet];
    1964                toIndex_[iSet] = -1;
    1965                int iKey = keyVariable_[iSet];
    1966                if (iKey < numberColumns) {
    1967                     double valueKey;
    1968                     if (getStatus(iSet) == ClpSimplex::atLowerBound)
    1969                          valueKey = lower_[iSet];
    1970                     else
    1971                          valueKey = upper_[iSet];
    1972                     double value = scalar * (x[iKey] - valueKey);
    1973                     if (value) {
    1974                          for (CoinBigIndex j = startColumn[iKey];
    1975                                    j < startColumn[iKey] + length[iKey]; j++) {
    1976                               int jRow = row[j];
    1977                               y[jRow] += value * element[j];
    1978                          }
    1979                     }
    1980                }
    1981           }
    1982      }
     1912  if (model_->specialOptions() != 16) {
     1913    ClpPackedMatrix::times(scalar, x, y);
     1914  } else {
     1915    int iRow;
     1916    int numberColumns = model_->numberColumns();
     1917    int numberRows = model_->numberRows();
     1918    const double *element = matrix_->getElements();
     1919    const int *row = matrix_->getIndices();
     1920    const CoinBigIndex *startColumn = matrix_->getVectorStarts();
     1921    const int *length = matrix_->getVectorLengths();
     1922    int *pivotVariable = model_->pivotVariable();
     1923    int numberToDo = 0;
     1924    for (iRow = 0; iRow < numberRows; iRow++) {
     1925      y[iRow] -= scalar * rhsOffset_[iRow];
     1926      int iColumn = pivotVariable[iRow];
     1927      if (iColumn < numberColumns) {
     1928        int iSet = backward_[iColumn];
     1929        if (iSet >= 0 && toIndex_[iSet] < 0) {
     1930          toIndex_[iSet] = 0;
     1931          fromIndex_[numberToDo++] = iSet;
     1932        }
     1933        CoinBigIndex j;
     1934        double value = scalar * x[iColumn];
     1935        if (value) {
     1936          for (j = startColumn[iColumn];
     1937               j < startColumn[iColumn] + length[iColumn]; j++) {
     1938            int jRow = row[j];
     1939            y[jRow] += value * element[j];
     1940          }
     1941        }
     1942      }
     1943    }
     1944    // and gubs which are interacting
     1945    for (int jSet = 0; jSet < numberToDo; jSet++) {
     1946      int iSet = fromIndex_[jSet];
     1947      toIndex_[iSet] = -1;
     1948      int iKey = keyVariable_[iSet];
     1949      if (iKey < numberColumns) {
     1950        double valueKey;
     1951        if (getStatus(iSet) == ClpSimplex::atLowerBound)
     1952          valueKey = lower_[iSet];
     1953        else
     1954          valueKey = upper_[iSet];
     1955        double value = scalar * (x[iKey] - valueKey);
     1956        if (value) {
     1957          for (CoinBigIndex j = startColumn[iKey];
     1958               j < startColumn[iKey] + length[iKey]; j++) {
     1959            int jRow = row[j];
     1960            y[jRow] += value * element[j];
     1961          }
     1962        }
     1963      }
     1964    }
     1965  }
    19831966}
    19841967/* Just for debug - may be extended to other matrix types later.
    19851968   Returns number and sum of primal infeasibilities.
    19861969*/
    1987 int
    1988 ClpGubDynamicMatrix::checkFeasible(ClpSimplex * /*model*/, double & sum) const
     1970int ClpGubDynamicMatrix::checkFeasible(ClpSimplex * /*model*/, double &sum) const
    19891971{
    1990      int numberRows = model_->numberRows();
    1991      double * rhs = new double[numberRows];
    1992      int numberColumns = model_->numberColumns();
    1993      int iRow;
    1994      CoinZeroN(rhs, numberRows);
    1995      // do ones at bounds before gub
    1996      const double * smallSolution = model_->solutionRegion();
    1997      const double * element = matrix_->getElements();
    1998      const int * row = matrix_->getIndices();
    1999      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    2000      const int * length = matrix_->getVectorLengths();
    2001      int iColumn;
    2002      int numberInfeasible = 0;
    2003      const double * rowLower = model_->rowLower();
    2004      const double * rowUpper = model_->rowUpper();
    2005      sum = 0.0;
    2006      for (iRow = 0; iRow < numberRows; iRow++) {
    2007           double value = smallSolution[numberColumns+iRow];
    2008           if (value < rowLower[iRow] - 1.0e-5 ||
    2009                     value > rowUpper[iRow] + 1.0e-5) {
    2010                //printf("row %d %g %g %g\n",
    2011                //     iRow,rowLower[iRow],value,rowUpper[iRow]);
    2012                numberInfeasible++;
    2013                sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
    2014           }
    2015           rhs[iRow] = value;
    2016      }
    2017      const double * columnLower = model_->columnLower();
    2018      const double * columnUpper = model_->columnUpper();
    2019      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
    2020           double value = smallSolution[iColumn];
    2021           if (value < columnLower[iColumn] - 1.0e-5 ||
    2022                     value > columnUpper[iColumn] + 1.0e-5) {
    2023                //printf("column %d %g %g %g\n",
    2024                //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
    2025                numberInfeasible++;
    2026                sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
    2027           }
    2028           for (CoinBigIndex j = startColumn[iColumn];
    2029                     j < startColumn[iColumn] + length[iColumn]; j++) {
    2030                int jRow = row[j];
    2031                rhs[jRow] -= value * element[j];
    2032           }
    2033      }
    2034      double * solution = new double [numberGubColumns_];
    2035      for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
    2036           double value = 0.0;
    2037           if(getDynamicStatus(iColumn) == atUpperBound)
    2038                value = upperColumn_[iColumn];
    2039           else if (lowerColumn_)
    2040                value = lowerColumn_[iColumn];
    2041           solution[iColumn] = value;
    2042      }
    2043      // ones in small and gub
    2044      for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
    2045           int jFull = id_[iColumn-firstDynamic_];
    2046           solution[jFull] = smallSolution[iColumn];
    2047      }
    2048      // fill in all basic in small model
    2049      int * pivotVariable = model_->pivotVariable();
    2050      for (iRow = 0; iRow < numberRows; iRow++) {
    2051           int iColumn = pivotVariable[iRow];
    2052           if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
    2053                int iSequence = id_[iColumn-firstDynamic_];
    2054                solution[iSequence] = smallSolution[iColumn];
    2055           }
    2056      }
    2057      // and now compute value to use for key
    2058      ClpSimplex::Status iStatus;
    2059      for (int iSet = 0; iSet < numberSets_; iSet++) {
    2060           iColumn = keyVariable_[iSet];
    2061           if (iColumn < numberColumns) {
    2062                int iSequence = id_[iColumn-firstDynamic_];
    2063                solution[iSequence] = 0.0;
    2064                double b = 0.0;
    2065                // key is structural - where is slack
    2066                iStatus = getStatus(iSet);
    2067                assert (iStatus != ClpSimplex::basic);
    2068                if (iStatus == ClpSimplex::atLowerBound)
    2069                     b = lower_[iSet];
    2070                else
    2071                     b = upper_[iSet];
    2072                // subtract out others at bounds
    2073                for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
    2074                     b -= solution[j];
    2075                solution[iSequence] = b;
    2076           }
    2077      }
    2078      for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
    2079           double value = solution[iColumn];
    2080           if ((lowerColumn_ && value < lowerColumn_[iColumn] - 1.0e-5) ||
    2081                     (!lowerColumn_ && value < -1.0e-5) ||
    2082                     (upperColumn_ && value > upperColumn_[iColumn] + 1.0e-5)) {
    2083                //printf("column %d %g %g %g\n",
    2084                //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
    2085                numberInfeasible++;
    2086           }
    2087           if (value) {
    2088                for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
    2089                     int iRow = row_[j];
    2090                     rhs[iRow] -= element_[j] * value;
    2091                }
    2092           }
    2093      }
    2094      for (iRow = 0; iRow < numberRows; iRow++) {
    2095           if (fabs(rhs[iRow]) > 1.0e-5)
    2096                printf("rhs mismatch %d %g\n", iRow, rhs[iRow]);
    2097      }
    2098      delete [] solution;
    2099      delete [] rhs;
    2100      return numberInfeasible;
     1972  int numberRows = model_->numberRows();
     1973  double *rhs = new double[numberRows];
     1974  int numberColumns = model_->numberColumns();
     1975  int iRow;
     1976  CoinZeroN(rhs, numberRows);
     1977  // do ones at bounds before gub
     1978  const double *smallSolution = model_->solutionRegion();
     1979  const double *element = matrix_->getElements();
     1980  const int *row = matrix_->getIndices();
     1981  const CoinBigIndex *startColumn = matrix_->getVectorStarts();
     1982  const int *length = matrix_->getVectorLengths();
     1983  int iColumn;
     1984  int numberInfeasible = 0;
     1985  const double *rowLower = model_->rowLower();
     1986  const double *rowUpper = model_->rowUpper();
     1987  sum = 0.0;
     1988  for (iRow = 0; iRow < numberRows; iRow++) {
     1989    double value = smallSolution[numberColumns + iRow];
     1990    if (value < rowLower[iRow] - 1.0e-5 || value > rowUpper[iRow] + 1.0e-5) {
     1991      //printf("row %d %g %g %g\n",
     1992      //     iRow,rowLower[iRow],value,rowUpper[iRow]);
     1993      numberInfeasible++;
     1994      sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
     1995    }
     1996    rhs[iRow] = value;
     1997  }
     1998  const double *columnLower = model_->columnLower();
     1999  const double *columnUpper = model_->columnUpper();
     2000  for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
     2001    double value = smallSolution[iColumn];
     2002    if (value < columnLower[iColumn] - 1.0e-5 || value > columnUpper[iColumn] + 1.0e-5) {
     2003      //printf("column %d %g %g %g\n",
     2004      //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
     2005      numberInfeasible++;
     2006      sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
     2007    }
     2008    for (CoinBigIndex j = startColumn[iColumn];
     2009         j < startColumn[iColumn] + length[iColumn]; j++) {
     2010      int jRow = row[j];
     2011      rhs[jRow] -= value * element[j];
     2012    }
     2013  }
     2014  double *solution = new double[numberGubColumns_];
     2015  for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
     2016    double value = 0.0;
     2017    if (getDynamicStatus(iColumn) == atUpperBound)
     2018      value = upperColumn_[iColumn];
     2019    else if (lowerColumn_)
     2020      value = lowerColumn_[iColumn];
     2021    solution[iColumn] = value;
     2022  }
     2023  // ones in small and gub
     2024  for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
     2025    int jFull = id_[iColumn - firstDynamic_];
     2026    solution[jFull] = smallSolution[iColumn];
     2027  }
     2028  // fill in all basic in small model
     2029  int *pivotVariable = model_->pivotVariable();
     2030  for (iRow = 0; iRow < numberRows; iRow++) {
     2031    int iColumn = pivotVariable[iRow];
     2032    if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
     2033      int iSequence = id_[iColumn - firstDynamic_];
     2034      solution[iSequence] = smallSolution[iColumn];
     2035    }
     2036  }
     2037  // and now compute value to use for key
     2038  ClpSimplex::Status iStatus;
     2039  for (int iSet = 0; iSet < numberSets_; iSet++) {
     2040    iColumn = keyVariable_[iSet];
     2041    if (iColumn < numberColumns) {
     2042      int iSequence = id_[iColumn - firstDynamic_];
     2043      solution[iSequence] = 0.0;
     2044      double b = 0.0;
     2045      // key is structural - where is slack
     2046      iStatus = getStatus(iSet);
     2047      assert(iStatus != ClpSimplex::basic);
     2048      if (iStatus == ClpSimplex::atLowerBound)
     2049        b = lower_[iSet];
     2050      else
     2051        b = upper_[iSet];
     2052      // subtract out others at bounds
     2053      for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
     2054        b -= solution[j];
     2055      solution[iSequence] = b;
     2056    }
     2057  }
     2058  for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
     2059    double value = solution[iColumn];
     2060    if ((lowerColumn_ && value < lowerColumn_[iColumn] - 1.0e-5) || (!lowerColumn_ && value < -1.0e-5) || (upperColumn_ && value > upperColumn_[iColumn] + 1.0e-5)) {
     2061      //printf("column %d %g %g %g\n",
     2062      //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
     2063      numberInfeasible++;
     2064    }
     2065    if (value) {
     2066      for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
     2067        int iRow = row_[j];
     2068        rhs[iRow] -= element_[j] * value;
     2069      }
     2070    }
     2071  }
     2072  for (iRow = 0; iRow < numberRows; iRow++) {
     2073    if (fabs(rhs[iRow]) > 1.0e-5)
     2074      printf("rhs mismatch %d %g\n", iRow, rhs[iRow]);
     2075  }
     2076  delete[] solution;
     2077  delete[] rhs;
     2078  return numberInfeasible;
    21012079}
    21022080// Cleans data after setWarmStart
    2103 void
    2104 ClpGubDynamicMatrix::cleanData(ClpSimplex * model)
     2081void ClpGubDynamicMatrix::cleanData(ClpSimplex *model)
    21052082{
    2106      // and redo chains
    2107      int numberColumns = model->numberColumns();
    2108      int iColumn;
    2109      // do backward
    2110      int * mark = new int [numberGubColumns_];
    2111      for (iColumn = 0; iColumn < numberGubColumns_; iColumn++)
    2112           mark[iColumn] = -1;
    2113      int i;
    2114      for (i = 0; i < firstDynamic_; i++) {
    2115           assert (backward_[i] == -1);
    2116           next_[i] = -1;
    2117      }
    2118      for (i = firstDynamic_; i < firstAvailable_; i++) {
    2119           iColumn = id_[i-firstDynamic_];
    2120           mark[iColumn] = i;
    2121      }
    2122      for (i = 0; i < numberSets_; i++) {
    2123           int iKey = keyVariable_[i];
    2124           int lastNext = -1;
    2125           int firstNext = -1;
    2126           for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i+1]; k++) {
    2127                iColumn = mark[k];
    2128                if (iColumn >= 0) {
    2129                     if (iColumn != iKey) {
    2130                          if (lastNext >= 0)
    2131                               next_[lastNext] = iColumn;
    2132                          else
    2133                               firstNext = iColumn;
    2134                          lastNext = iColumn;
    2135                     }
    2136                     backward_[iColumn] = i;
    2137                }
    2138           }
    2139           setFeasible(i);
    2140           if (firstNext >= 0) {
    2141                // others
    2142                next_[iKey] = firstNext;
    2143                next_[lastNext] = -(iKey + 1);
    2144           } else if (iKey < numberColumns) {
    2145                next_[iKey] = -(iKey + 1);
    2146           }
    2147      }
    2148      delete [] mark;
    2149      // fill matrix
    2150      double * element = matrix_->getMutableElements();
    2151      int * row = matrix_->getMutableIndices();
    2152      CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
    2153      int * length = matrix_->getMutableVectorLengths();
    2154      CoinBigIndex numberElements = startColumn[firstDynamic_];
    2155      for (i = firstDynamic_; i < firstAvailable_; i++) {
    2156           int iColumn = id_[i-firstDynamic_];
    2157           int numberThis = startColumn_[iColumn+1] - startColumn_[iColumn];
    2158           length[i] = numberThis;
    2159           for (CoinBigIndex jBigIndex = startColumn_[iColumn];
    2160                     jBigIndex < startColumn_[iColumn+1]; jBigIndex++) {
    2161                row[numberElements] = row_[jBigIndex];
    2162                element[numberElements++] = element_[jBigIndex];
    2163           }
    2164           startColumn[i+1] = numberElements;
    2165      }
     2083  // and redo chains
     2084  int numberColumns = model->numberColumns();
     2085  int iColumn;
     2086  // do backward
     2087  int *mark = new int[numberGubColumns_];
     2088  for (iColumn = 0; iColumn < numberGubColumns_; iColumn++)
     2089    mark[iColumn] = -1;
     2090  int i;
     2091  for (i = 0; i < firstDynamic_; i++) {
     2092    assert(backward_[i] == -1);
     2093    next_[i] = -1;
     2094  }
     2095  for (i = firstDynamic_; i < firstAvailable_; i++) {
     2096    iColumn = id_[i - firstDynamic_];
     2097    mark[iColumn] = i;
     2098  }
     2099  for (i = 0; i < numberSets_; i++) {
     2100    int iKey = keyVariable_[i];
     2101    int lastNext = -1;
     2102    int firstNext = -1;
     2103    for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i + 1]; k++) {
     2104      iColumn = mark[k];
     2105      if (iColumn >= 0) {
     2106        if (iColumn != iKey) {
     2107          if (lastNext >= 0)
     2108            next_[lastNext] = iColumn;
     2109          else
     2110            firstNext = iColumn;
     2111          lastNext = iColumn;
     2112        }
     2113        backward_[iColumn] = i;
     2114      }
     2115    }
     2116    setFeasible(i);
     2117    if (firstNext >= 0) {
     2118      // others
     2119      next_[iKey] = firstNext;
     2120      next_[lastNext] = -(iKey + 1);
     2121    } else if (iKey < numberColumns) {
     2122      next_[iKey] = -(iKey + 1);
     2123    }
     2124  }
     2125  delete[] mark;
     2126  // fill matrix
     2127  double *element = matrix_->getMutableElements();
     2128  int *row = matrix_->getMutableIndices();
     2129  CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
     2130  int *length = matrix_->getMutableVectorLengths();
     2131  CoinBigIndex numberElements = startColumn[firstDynamic_];
     2132  for (i = firstDynamic_; i < firstAvailable_; i++) {
     2133    int iColumn = id_[i - firstDynamic_];
     2134    int numberThis = startColumn_[iColumn + 1] - startColumn_[iColumn];
     2135    length[i] = numberThis;
     2136    for (CoinBigIndex jBigIndex = startColumn_[iColumn];
     2137         jBigIndex < startColumn_[iColumn + 1]; jBigIndex++) {
     2138      row[numberElements] = row_[jBigIndex];
     2139      element[numberElements++] = element_[jBigIndex];
     2140    }
     2141    startColumn[i + 1] = numberElements;
     2142  }
    21662143}
     2144
     2145/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     2146*/
Note: See TracChangeset for help on using the changeset viewer.