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/ClpGubMatrix.cpp

    r2271 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>
     
    2827// Default Constructor
    2928//-------------------------------------------------------------------
    30 ClpGubMatrix::ClpGubMatrix ()
    31      : ClpPackedMatrix(),
    32        sumDualInfeasibilities_(0.0),
    33        sumPrimalInfeasibilities_(0.0),
    34        sumOfRelaxedDualInfeasibilities_(0.0),
    35        sumOfRelaxedPrimalInfeasibilities_(0.0),
    36        infeasibilityWeight_(0.0),
    37        start_(NULL),
    38        end_(NULL),
    39        lower_(NULL),
    40        upper_(NULL),
    41        status_(NULL),
    42        saveStatus_(NULL),
    43        savedKeyVariable_(NULL),
    44        backward_(NULL),
    45        backToPivotRow_(NULL),
    46        changeCost_(NULL),
    47        keyVariable_(NULL),
    48        next_(NULL),
    49        toIndex_(NULL),
    50        fromIndex_(NULL),
    51        model_(NULL),
    52        numberDualInfeasibilities_(0),
    53        numberPrimalInfeasibilities_(0),
    54        noCheck_(-1),
    55        numberSets_(0),
    56        saveNumber_(0),
    57        possiblePivotKey_(0),
    58        gubSlackIn_(-1),
    59        firstGub_(0),
    60        lastGub_(0),
    61       gubType_(0)
     29ClpGubMatrix::ClpGubMatrix()
     30  : ClpPackedMatrix()
     31  , sumDualInfeasibilities_(0.0)
     32  , sumPrimalInfeasibilities_(0.0)
     33  , sumOfRelaxedDualInfeasibilities_(0.0)
     34  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     35  , infeasibilityWeight_(0.0)
     36  , start_(NULL)
     37  , end_(NULL)
     38  , lower_(NULL)
     39  , upper_(NULL)
     40  , status_(NULL)
     41  , saveStatus_(NULL)
     42  , savedKeyVariable_(NULL)
     43  , backward_(NULL)
     44  , backToPivotRow_(NULL)
     45  , changeCost_(NULL)
     46  , keyVariable_(NULL)
     47  , next_(NULL)
     48  , toIndex_(NULL)
     49  , fromIndex_(NULL)
     50  , model_(NULL)
     51  , numberDualInfeasibilities_(0)
     52  , numberPrimalInfeasibilities_(0)
     53  , noCheck_(-1)
     54  , numberSets_(0)
     55  , saveNumber_(0)
     56  , possiblePivotKey_(0)
     57  , gubSlackIn_(-1)
     58  , firstGub_(0)
     59  , lastGub_(0)
     60  , gubType_(0)
    6261{
    63      setType(16);
     62  setType(16);
    6463}
    6564
     
    6766// Copy constructor
    6867//-------------------------------------------------------------------
    69 ClpGubMatrix::ClpGubMatrix (const ClpGubMatrix & rhs)
    70      : ClpPackedMatrix(rhs)
     68ClpGubMatrix::ClpGubMatrix(const ClpGubMatrix &rhs)
     69  : ClpPackedMatrix(rhs)
    7170{
    72      numberSets_ = rhs.numberSets_;
    73      saveNumber_ = rhs.saveNumber_;
    74      possiblePivotKey_ = rhs.possiblePivotKey_;
    75      gubSlackIn_ = rhs.gubSlackIn_;
    76      start_ = ClpCopyOfArray(rhs.start_, numberSets_);
    77      end_ = ClpCopyOfArray(rhs.end_, numberSets_);
    78      lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
    79      upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
    80      status_ = ClpCopyOfArray(rhs.status_, numberSets_);
    81      saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
    82      savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
    83      int numberColumns = getNumCols();
    84      backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
    85      backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
    86      changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
    87      fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
    88      keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
    89      // find longest set
    90      int * longest = new int[numberSets_];
    91      CoinZeroN(longest, numberSets_);
    92      int j;
    93      for (j = 0; j < numberColumns; j++) {
    94           int iSet = backward_[j];
    95           if (iSet >= 0)
    96                longest[iSet]++;
    97      }
    98      int length = 0;
    99      for (j = 0; j < numberSets_; j++)
    100           length = CoinMax(length, longest[j]);
    101      next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
    102      toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
    103      sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
    104      sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    105      sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    106      sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    107      infeasibilityWeight_ = rhs.infeasibilityWeight_;
    108      numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    109      numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    110      noCheck_ = rhs.noCheck_;
    111      firstGub_ = rhs.firstGub_;
    112      lastGub_ = rhs.lastGub_;
    113      gubType_ = rhs.gubType_;
    114      model_ = rhs.model_;
     71  numberSets_ = rhs.numberSets_;
     72  saveNumber_ = rhs.saveNumber_;
     73  possiblePivotKey_ = rhs.possiblePivotKey_;
     74  gubSlackIn_ = rhs.gubSlackIn_;
     75  start_ = ClpCopyOfArray(rhs.start_, numberSets_);
     76  end_ = ClpCopyOfArray(rhs.end_, numberSets_);
     77  lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
     78  upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
     79  status_ = ClpCopyOfArray(rhs.status_, numberSets_);
     80  saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
     81  savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
     82  int numberColumns = getNumCols();
     83  backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
     84  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
     85  changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
     86  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
     87  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
     88  // find longest set
     89  int *longest = new int[numberSets_];
     90  CoinZeroN(longest, numberSets_);
     91  int j;
     92  for (j = 0; j < numberColumns; j++) {
     93    int iSet = backward_[j];
     94    if (iSet >= 0)
     95      longest[iSet]++;
     96  }
     97  int length = 0;
     98  for (j = 0; j < numberSets_; j++)
     99    length = CoinMax(length, longest[j]);
     100  next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
     101  toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
     102  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
     103  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
     104  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
     105  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     106  infeasibilityWeight_ = rhs.infeasibilityWeight_;
     107  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
     108  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
     109  noCheck_ = rhs.noCheck_;
     110  firstGub_ = rhs.firstGub_;
     111  lastGub_ = rhs.lastGub_;
     112  gubType_ = rhs.gubType_;
     113  model_ = rhs.model_;
    115114}
    116115
     
    118117// assign matrix (for space reasons)
    119118//-------------------------------------------------------------------
    120 ClpGubMatrix::ClpGubMatrix (CoinPackedMatrix * rhs)
    121      : ClpPackedMatrix(rhs),
    122        sumDualInfeasibilities_(0.0),
    123        sumPrimalInfeasibilities_(0.0),
    124        sumOfRelaxedDualInfeasibilities_(0.0),
    125        sumOfRelaxedPrimalInfeasibilities_(0.0),
    126        infeasibilityWeight_(0.0),
    127        start_(NULL),
    128        end_(NULL),
    129        lower_(NULL),
    130        upper_(NULL),
    131        status_(NULL),
    132        saveStatus_(NULL),
    133        savedKeyVariable_(NULL),
    134        backward_(NULL),
    135        backToPivotRow_(NULL),
    136        changeCost_(NULL),
    137        keyVariable_(NULL),
    138        next_(NULL),
    139        toIndex_(NULL),
    140        fromIndex_(NULL),
    141        model_(NULL),
    142        numberDualInfeasibilities_(0),
    143        numberPrimalInfeasibilities_(0),
    144        noCheck_(-1),
    145        numberSets_(0),
    146        saveNumber_(0),
    147        possiblePivotKey_(0),
    148        gubSlackIn_(-1),
    149        firstGub_(0),
    150        lastGub_(0),
    151       gubType_(0)
     119ClpGubMatrix::ClpGubMatrix(CoinPackedMatrix *rhs)
     120  : ClpPackedMatrix(rhs)
     121  , sumDualInfeasibilities_(0.0)
     122  , sumPrimalInfeasibilities_(0.0)
     123  , sumOfRelaxedDualInfeasibilities_(0.0)
     124  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     125  , infeasibilityWeight_(0.0)
     126  , start_(NULL)
     127  , end_(NULL)
     128  , lower_(NULL)
     129  , upper_(NULL)
     130  , status_(NULL)
     131  , saveStatus_(NULL)
     132  , savedKeyVariable_(NULL)
     133  , backward_(NULL)
     134  , backToPivotRow_(NULL)
     135  , changeCost_(NULL)
     136  , keyVariable_(NULL)
     137  , next_(NULL)
     138  , toIndex_(NULL)
     139  , fromIndex_(NULL)
     140  , model_(NULL)
     141  , numberDualInfeasibilities_(0)
     142  , numberPrimalInfeasibilities_(0)
     143  , noCheck_(-1)
     144  , numberSets_(0)
     145  , saveNumber_(0)
     146  , possiblePivotKey_(0)
     147  , gubSlackIn_(-1)
     148  , firstGub_(0)
     149  , lastGub_(0)
     150  , gubType_(0)
    152151{
    153      setType(16);
     152  setType(16);
    154153}
    155154
    156155/* This takes over ownership (for space reasons) and is the
    157156   real constructor*/
    158 ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix * matrix, int numberSets,
    159                            const int * start, const int * end,
    160                            const double * lower, const double * upper,
    161                            const unsigned char * status)
    162      : ClpPackedMatrix(matrix->matrix()),
    163        sumDualInfeasibilities_(0.0),
    164        sumPrimalInfeasibilities_(0.0),
    165        sumOfRelaxedDualInfeasibilities_(0.0),
    166        sumOfRelaxedPrimalInfeasibilities_(0.0),
    167        numberDualInfeasibilities_(0),
    168        numberPrimalInfeasibilities_(0),
    169        saveNumber_(0),
    170        possiblePivotKey_(0),
    171       gubSlackIn_(-1)
     157ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix *matrix, int numberSets,
     158  const int *start, const int *end,
     159  const double *lower, const double *upper,
     160  const unsigned char *status)
     161  : ClpPackedMatrix(matrix->matrix())
     162  , sumDualInfeasibilities_(0.0)
     163  , sumPrimalInfeasibilities_(0.0)
     164  , sumOfRelaxedDualInfeasibilities_(0.0)
     165  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     166  , numberDualInfeasibilities_(0)
     167  , numberPrimalInfeasibilities_(0)
     168  , saveNumber_(0)
     169  , possiblePivotKey_(0)
     170  , gubSlackIn_(-1)
    172171{
    173      model_ = NULL;
    174      numberSets_ = numberSets;
    175      start_ = ClpCopyOfArray(start, numberSets_);
    176      end_ = ClpCopyOfArray(end, numberSets_);
    177      lower_ = ClpCopyOfArray(lower, numberSets_);
    178      upper_ = ClpCopyOfArray(upper, numberSets_);
    179      // Check valid and ordered
    180      int last = -1;
    181      int numberColumns = matrix_->getNumCols();
    182      int numberRows = matrix_->getNumRows();
    183      backward_ = new int[numberColumns];
    184      backToPivotRow_ = new int[numberColumns];
    185      changeCost_ = new double [numberRows+numberSets_];
    186      keyVariable_ = new int[numberSets_];
    187      // signal to need new ordering
    188      next_ = NULL;
    189      for (int iColumn = 0; iColumn < numberColumns; iColumn++)
    190           backward_[iColumn] = -1;
     172  model_ = NULL;
     173  numberSets_ = numberSets;
     174  start_ = ClpCopyOfArray(start, numberSets_);
     175  end_ = ClpCopyOfArray(end, numberSets_);
     176  lower_ = ClpCopyOfArray(lower, numberSets_);
     177  upper_ = ClpCopyOfArray(upper, numberSets_);
     178  // Check valid and ordered
     179  int last = -1;
     180  int numberColumns = matrix_->getNumCols();
     181  int numberRows = matrix_->getNumRows();
     182  backward_ = new int[numberColumns];
     183  backToPivotRow_ = new int[numberColumns];
     184  changeCost_ = new double[numberRows + numberSets_];
     185  keyVariable_ = new int[numberSets_];
     186  // signal to need new ordering
     187  next_ = NULL;
     188  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
     189    backward_[iColumn] = -1;
    191190
    192      int iSet;
    193      for (iSet = 0; iSet < numberSets_; iSet++) {
    194           // set key variable as slack
    195           keyVariable_[iSet] = iSet + numberColumns;
    196           if (start_[iSet] < 0 || start_[iSet] >= numberColumns)
    197                throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
    198           if (end_[iSet] < 0 || end_[iSet] > numberColumns)
    199                throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
    200           if (end_[iSet] <= start_[iSet])
    201                throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix");
    202           if (start_[iSet] < last)
    203                throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix");
    204           last = end_[iSet];
    205           int j;
    206           for (j = start_[iSet]; j < end_[iSet]; j++)
    207                backward_[j] = iSet;
    208      }
    209      // Find type of gub
    210      firstGub_ = numberColumns + 1;
    211      lastGub_ = -1;
    212      int i;
    213      for (i = 0; i < numberColumns; i++) {
    214           if (backward_[i] >= 0) {
    215                firstGub_ = CoinMin(firstGub_, i);
    216                lastGub_ = CoinMax(lastGub_, i);
    217           }
    218      }
    219      gubType_ = 0;
    220      // adjust lastGub_
    221      if (lastGub_ > 0)
    222           lastGub_++;
    223      for (i = firstGub_; i < lastGub_; i++) {
    224           if (backward_[i] < 0) {
    225                gubType_ = 1;
    226                printf("interior non gub %d\n", i);
    227                break;
    228           }
    229      }
    230      if (status) {
    231           status_ = ClpCopyOfArray(status, numberSets_);
    232      } else {
    233           status_ = new unsigned char [numberSets_];
    234           memset(status_, 0, numberSets_);
    235           int i;
    236           for (i = 0; i < numberSets_; i++) {
    237                // make slack key
    238                setStatus(i, ClpSimplex::basic);
    239           }
    240      }
    241      saveStatus_ = new unsigned char [numberSets_];
    242      memset(saveStatus_, 0, numberSets_);
    243      savedKeyVariable_ = new int [numberSets_];
    244      memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
    245      noCheck_ = -1;
    246      infeasibilityWeight_ = 0.0;
     191  int iSet;
     192  for (iSet = 0; iSet < numberSets_; iSet++) {
     193    // set key variable as slack
     194    keyVariable_[iSet] = iSet + numberColumns;
     195    if (start_[iSet] < 0 || start_[iSet] >= numberColumns)
     196      throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
     197    if (end_[iSet] < 0 || end_[iSet] > numberColumns)
     198      throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
     199    if (end_[iSet] <= start_[iSet])
     200      throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix");
     201    if (start_[iSet] < last)
     202      throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix");
     203    last = end_[iSet];
     204    int j;
     205    for (j = start_[iSet]; j < end_[iSet]; j++)
     206      backward_[j] = iSet;
     207  }
     208  // Find type of gub
     209  firstGub_ = numberColumns + 1;
     210  lastGub_ = -1;
     211  int i;
     212  for (i = 0; i < numberColumns; i++) {
     213    if (backward_[i] >= 0) {
     214      firstGub_ = CoinMin(firstGub_, i);
     215      lastGub_ = CoinMax(lastGub_, i);
     216    }
     217  }
     218  gubType_ = 0;
     219  // adjust lastGub_
     220  if (lastGub_ > 0)
     221    lastGub_++;
     222  for (i = firstGub_; i < lastGub_; i++) {
     223    if (backward_[i] < 0) {
     224      gubType_ = 1;
     225      printf("interior non gub %d\n", i);
     226      break;
     227    }
     228  }
     229  if (status) {
     230    status_ = ClpCopyOfArray(status, numberSets_);
     231  } else {
     232    status_ = new unsigned char[numberSets_];
     233    memset(status_, 0, numberSets_);
     234    int i;
     235    for (i = 0; i < numberSets_; i++) {
     236      // make slack key
     237      setStatus(i, ClpSimplex::basic);
     238    }
     239  }
     240  saveStatus_ = new unsigned char[numberSets_];
     241  memset(saveStatus_, 0, numberSets_);
     242  savedKeyVariable_ = new int[numberSets_];
     243  memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
     244  noCheck_ = -1;
     245  infeasibilityWeight_ = 0.0;
    247246}
    248247
    249 ClpGubMatrix::ClpGubMatrix (const CoinPackedMatrix & rhs)
    250      : ClpPackedMatrix(rhs),
    251        sumDualInfeasibilities_(0.0),
    252        sumPrimalInfeasibilities_(0.0),
    253        sumOfRelaxedDualInfeasibilities_(0.0),
    254        sumOfRelaxedPrimalInfeasibilities_(0.0),
    255        infeasibilityWeight_(0.0),
    256        start_(NULL),
    257        end_(NULL),
    258        lower_(NULL),
    259        upper_(NULL),
    260        status_(NULL),
    261        saveStatus_(NULL),
    262        savedKeyVariable_(NULL),
    263        backward_(NULL),
    264        backToPivotRow_(NULL),
    265        changeCost_(NULL),
    266        keyVariable_(NULL),
    267        next_(NULL),
    268        toIndex_(NULL),
    269        fromIndex_(NULL),
    270        model_(NULL),
    271        numberDualInfeasibilities_(0),
    272        numberPrimalInfeasibilities_(0),
    273        noCheck_(-1),
    274        numberSets_(0),
    275        saveNumber_(0),
    276        possiblePivotKey_(0),
    277        gubSlackIn_(-1),
    278        firstGub_(0),
    279        lastGub_(0),
    280       gubType_(0)
     248ClpGubMatrix::ClpGubMatrix(const CoinPackedMatrix &rhs)
     249  : ClpPackedMatrix(rhs)
     250  , sumDualInfeasibilities_(0.0)
     251  , sumPrimalInfeasibilities_(0.0)
     252  , sumOfRelaxedDualInfeasibilities_(0.0)
     253  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     254  , infeasibilityWeight_(0.0)
     255  , start_(NULL)
     256  , end_(NULL)
     257  , lower_(NULL)
     258  , upper_(NULL)
     259  , status_(NULL)
     260  , saveStatus_(NULL)
     261  , savedKeyVariable_(NULL)
     262  , backward_(NULL)
     263  , backToPivotRow_(NULL)
     264  , changeCost_(NULL)
     265  , keyVariable_(NULL)
     266  , next_(NULL)
     267  , toIndex_(NULL)
     268  , fromIndex_(NULL)
     269  , model_(NULL)
     270  , numberDualInfeasibilities_(0)
     271  , numberPrimalInfeasibilities_(0)
     272  , noCheck_(-1)
     273  , numberSets_(0)
     274  , saveNumber_(0)
     275  , possiblePivotKey_(0)
     276  , gubSlackIn_(-1)
     277  , firstGub_(0)
     278  , lastGub_(0)
     279  , gubType_(0)
    281280{
    282      setType(16);
    283 
     281  setType(16);
    284282}
    285283
     
    287285// Destructor
    288286//-------------------------------------------------------------------
    289 ClpGubMatrix::~ClpGubMatrix ()
     287ClpGubMatrix::~ClpGubMatrix()
    290288{
    291      delete [] start_;
    292      delete [] end_;
    293      delete [] lower_;
    294      delete [] upper_;
    295      delete [] status_;
    296      delete [] saveStatus_;
    297      delete [] savedKeyVariable_;
    298      delete [] backward_;
    299      delete [] backToPivotRow_;
    300      delete [] changeCost_;
    301      delete [] keyVariable_;
    302      delete [] next_;
    303      delete [] toIndex_;
    304      delete [] fromIndex_;
     289  delete[] start_;
     290  delete[] end_;
     291  delete[] lower_;
     292  delete[] upper_;
     293  delete[] status_;
     294  delete[] saveStatus_;
     295  delete[] savedKeyVariable_;
     296  delete[] backward_;
     297  delete[] backToPivotRow_;
     298  delete[] changeCost_;
     299  delete[] keyVariable_;
     300  delete[] next_;
     301  delete[] toIndex_;
     302  delete[] fromIndex_;
    305303}
    306304
     
    309307//-------------------------------------------------------------------
    310308ClpGubMatrix &
    311 ClpGubMatrix::operator=(const ClpGubMatrix& rhs)
     309ClpGubMatrix::operator=(const ClpGubMatrix &rhs)
    312310{
    313      if (this != &rhs) {
    314           ClpPackedMatrix::operator=(rhs);
    315           delete [] start_;
    316           delete [] end_;
    317           delete [] lower_;
    318           delete [] upper_;
    319           delete [] status_;
    320           delete [] saveStatus_;
    321           delete [] savedKeyVariable_;
    322           delete [] backward_;
    323           delete [] backToPivotRow_;
    324           delete [] changeCost_;
    325           delete [] keyVariable_;
    326           delete [] next_;
    327           delete [] toIndex_;
    328           delete [] fromIndex_;
    329           numberSets_ = rhs.numberSets_;
    330           saveNumber_ = rhs.saveNumber_;
    331           possiblePivotKey_ = rhs.possiblePivotKey_;
    332           gubSlackIn_ = rhs.gubSlackIn_;
    333           start_ = ClpCopyOfArray(rhs.start_, numberSets_);
    334           end_ = ClpCopyOfArray(rhs.end_, numberSets_);
    335           lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
    336           upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
    337           status_ = ClpCopyOfArray(rhs.status_, numberSets_);
    338           saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
    339           savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
    340           int numberColumns = getNumCols();
    341           backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
    342           backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
    343           changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
    344           fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
    345           keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
    346           // find longest set
    347           int * longest = new int[numberSets_];
    348           CoinZeroN(longest, numberSets_);
    349           int j;
    350           for (j = 0; j < numberColumns; j++) {
    351                int iSet = backward_[j];
    352                if (iSet >= 0)
    353                     longest[iSet]++;
    354           }
    355           int length = 0;
    356           for (j = 0; j < numberSets_; j++)
    357                length = CoinMax(length, longest[j]);
    358           next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
    359           toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
    360           sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
    361           sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    362           sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    363           sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    364           infeasibilityWeight_ = rhs.infeasibilityWeight_;
    365           numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    366           numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    367           noCheck_ = rhs.noCheck_;
    368           firstGub_ = rhs.firstGub_;
    369           lastGub_ = rhs.lastGub_;
    370           gubType_ = rhs.gubType_;
    371           model_ = rhs.model_;
    372      }
    373      return *this;
     311  if (this != &rhs) {
     312    ClpPackedMatrix::operator=(rhs);
     313    delete[] start_;
     314    delete[] end_;
     315    delete[] lower_;
     316    delete[] upper_;
     317    delete[] status_;
     318    delete[] saveStatus_;
     319    delete[] savedKeyVariable_;
     320    delete[] backward_;
     321    delete[] backToPivotRow_;
     322    delete[] changeCost_;
     323    delete[] keyVariable_;
     324    delete[] next_;
     325    delete[] toIndex_;
     326    delete[] fromIndex_;
     327    numberSets_ = rhs.numberSets_;
     328    saveNumber_ = rhs.saveNumber_;
     329    possiblePivotKey_ = rhs.possiblePivotKey_;
     330    gubSlackIn_ = rhs.gubSlackIn_;
     331    start_ = ClpCopyOfArray(rhs.start_, numberSets_);
     332    end_ = ClpCopyOfArray(rhs.end_, numberSets_);
     333    lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
     334    upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
     335    status_ = ClpCopyOfArray(rhs.status_, numberSets_);
     336    saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
     337    savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
     338    int numberColumns = getNumCols();
     339    backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
     340    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
     341    changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
     342    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
     343    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
     344    // find longest set
     345    int *longest = new int[numberSets_];
     346    CoinZeroN(longest, numberSets_);
     347    int j;
     348    for (j = 0; j < numberColumns; j++) {
     349      int iSet = backward_[j];
     350      if (iSet >= 0)
     351        longest[iSet]++;
     352    }
     353    int length = 0;
     354    for (j = 0; j < numberSets_; j++)
     355      length = CoinMax(length, longest[j]);
     356    next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
     357    toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
     358    sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
     359    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
     360    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
     361    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     362    infeasibilityWeight_ = rhs.infeasibilityWeight_;
     363    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
     364    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
     365    noCheck_ = rhs.noCheck_;
     366    firstGub_ = rhs.firstGub_;
     367    lastGub_ = rhs.lastGub_;
     368    gubType_ = rhs.gubType_;
     369    model_ = rhs.model_;
     370  }
     371  return *this;
    374372}
    375373//-------------------------------------------------------------------
    376374// Clone
    377375//-------------------------------------------------------------------
    378 ClpMatrixBase * ClpGubMatrix::clone() const
     376ClpMatrixBase *ClpGubMatrix::clone() const
    379377{
    380      return new ClpGubMatrix(*this);
     378  return new ClpGubMatrix(*this);
    381379}
    382380/* Subset clone (without gaps).  Duplicates are allowed
    383381   and order is as given */
    384382ClpMatrixBase *
    385 ClpGubMatrix::subsetClone (int numberRows, const int * whichRows,
    386                            int numberColumns,
    387                            const int * whichColumns) const
     383ClpGubMatrix::subsetClone(int numberRows, const int *whichRows,
     384  int numberColumns,
     385  const int *whichColumns) const
    388386{
    389      return new ClpGubMatrix(*this, numberRows, whichRows,
    390                              numberColumns, whichColumns);
     387  return new ClpGubMatrix(*this, numberRows, whichRows,
     388    numberColumns, whichColumns);
    391389}
    392390/* Returns a new matrix in reverse order without gaps
     
    395393ClpGubMatrix::reverseOrderedCopy() const
    396394{
    397      return NULL;
     395  return NULL;
    398396}
    399 int
    400 ClpGubMatrix::hiddenRows() const
     397int ClpGubMatrix::hiddenRows() const
    401398{
    402      return numberSets_;
     399  return numberSets_;
    403400}
    404401/* Subset constructor (without gaps).  Duplicates are allowed
    405402   and order is as given */
    406 ClpGubMatrix::ClpGubMatrix (
    407      const ClpGubMatrix & rhs,
    408      int numberRows, const int * whichRows,
    409      int numberColumns, const int * whichColumns)
    410      : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
     403ClpGubMatrix::ClpGubMatrix(
     404  const ClpGubMatrix &rhs,
     405  int numberRows, const int *whichRows,
     406  int numberColumns, const int *whichColumns)
     407  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
    411408{
    412      // Assuming no gub rows deleted
    413      // We also assume all sets in same order
    414      // Get array with backward pointers
    415      int numberColumnsOld = rhs.matrix_->getNumCols();
    416      int * array = new int [ numberColumnsOld];
    417      int i;
    418      for (i = 0; i < numberColumnsOld; i++)
    419           array[i] = -1;
    420      for (int iSet = 0; iSet < numberSets_; iSet++) {
    421           for (int j = start_[iSet]; j < end_[iSet]; j++)
    422                array[j] = iSet;
    423      }
    424      numberSets_ = -1;
    425      int lastSet = -1;
    426      bool inSet = false;
    427      for (i = 0; i < numberColumns; i++) {
    428           int iColumn = whichColumns[i];
    429           int iSet = array[iColumn];
    430           if (iSet < 0) {
    431                inSet = false;
    432           } else {
    433                if (!inSet) {
    434                     // start of new set but check okay
    435                     if (iSet <= lastSet)
    436                          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
    437                     lastSet = iSet;
    438                     numberSets_++;
    439                     start_[numberSets_] = i;
    440                     end_[numberSets_] = i + 1;
    441                     lower_[numberSets_] = lower_[iSet];
    442                     upper_[numberSets_] = upper_[iSet];
    443                     inSet = true;
    444                } else {
    445                     if (iSet < lastSet) {
    446                          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
    447                     } else if (iSet == lastSet) {
    448                          end_[numberSets_] = i + 1;
    449                     } else {
    450                          // new set
    451                          lastSet = iSet;
    452                          numberSets_++;
    453                          start_[numberSets_] = i;
    454                          end_[numberSets_] = i + 1;
    455                          lower_[numberSets_] = lower_[iSet];
    456                          upper_[numberSets_] = upper_[iSet];
    457                     }
    458                }
    459           }
    460      }
    461      delete [] array;
    462      numberSets_++; // adjust
    463      // Find type of gub
    464      firstGub_ = numberColumns + 1;
    465      lastGub_ = -1;
    466      for (i = 0; i < numberColumns; i++) {
    467           if (backward_[i] >= 0) {
    468                firstGub_ = CoinMin(firstGub_, i);
    469                lastGub_ = CoinMax(lastGub_, i);
    470           }
    471      }
    472      if (lastGub_ > 0)
    473           lastGub_++;
    474      gubType_ = 0;
    475      for (i = firstGub_; i < lastGub_; i++) {
    476           if (backward_[i] < 0) {
    477                gubType_ = 1;
    478                break;
    479           }
    480      }
     409  // Assuming no gub rows deleted
     410  // We also assume all sets in same order
     411  // Get array with backward pointers
     412  int numberColumnsOld = rhs.matrix_->getNumCols();
     413  int *array = new int[numberColumnsOld];
     414  int i;
     415  for (i = 0; i < numberColumnsOld; i++)
     416    array[i] = -1;
     417  for (int iSet = 0; iSet < numberSets_; iSet++) {
     418    for (int j = start_[iSet]; j < end_[iSet]; j++)
     419      array[j] = iSet;
     420  }
     421  numberSets_ = -1;
     422  int lastSet = -1;
     423  bool inSet = false;
     424  for (i = 0; i < numberColumns; i++) {
     425    int iColumn = whichColumns[i];
     426    int iSet = array[iColumn];
     427    if (iSet < 0) {
     428      inSet = false;
     429    } else {
     430      if (!inSet) {
     431        // start of new set but check okay
     432        if (iSet <= lastSet)
     433          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
     434        lastSet = iSet;
     435        numberSets_++;
     436        start_[numberSets_] = i;
     437        end_[numberSets_] = i + 1;
     438        lower_[numberSets_] = lower_[iSet];
     439        upper_[numberSets_] = upper_[iSet];
     440        inSet = true;
     441      } else {
     442        if (iSet < lastSet) {
     443          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
     444        } else if (iSet == lastSet) {
     445          end_[numberSets_] = i + 1;
     446        } else {
     447          // new set
     448          lastSet = iSet;
     449          numberSets_++;
     450          start_[numberSets_] = i;
     451          end_[numberSets_] = i + 1;
     452          lower_[numberSets_] = lower_[iSet];
     453          upper_[numberSets_] = upper_[iSet];
     454        }
     455      }
     456    }
     457  }
     458  delete[] array;
     459  numberSets_++; // adjust
     460  // Find type of gub
     461  firstGub_ = numberColumns + 1;
     462  lastGub_ = -1;
     463  for (i = 0; i < numberColumns; i++) {
     464    if (backward_[i] >= 0) {
     465      firstGub_ = CoinMin(firstGub_, i);
     466      lastGub_ = CoinMax(lastGub_, i);
     467    }
     468  }
     469  if (lastGub_ > 0)
     470    lastGub_++;
     471  gubType_ = 0;
     472  for (i = firstGub_; i < lastGub_; i++) {
     473    if (backward_[i] < 0) {
     474      gubType_ = 1;
     475      break;
     476    }
     477  }
    481478
    482      // Make sure key is feasible if only key in set
     479  // Make sure key is feasible if only key in set
    483480}
    484 ClpGubMatrix::ClpGubMatrix (
    485      const CoinPackedMatrix & rhs,
    486      int numberRows, const int * whichRows,
    487      int numberColumns, const int * whichColumns)
    488      : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns),
    489        sumDualInfeasibilities_(0.0),
    490        sumPrimalInfeasibilities_(0.0),
    491        sumOfRelaxedDualInfeasibilities_(0.0),
    492        sumOfRelaxedPrimalInfeasibilities_(0.0),
    493        start_(NULL),
    494        end_(NULL),
    495        lower_(NULL),
    496        upper_(NULL),
    497        backward_(NULL),
    498        backToPivotRow_(NULL),
    499        changeCost_(NULL),
    500        keyVariable_(NULL),
    501        next_(NULL),
    502        toIndex_(NULL),
    503        fromIndex_(NULL),
    504        numberDualInfeasibilities_(0),
    505        numberPrimalInfeasibilities_(0),
    506        numberSets_(0),
    507        saveNumber_(0),
    508        possiblePivotKey_(0),
    509        gubSlackIn_(-1),
    510        firstGub_(0),
    511        lastGub_(0),
    512       gubType_(0)
     481ClpGubMatrix::ClpGubMatrix(
     482  const CoinPackedMatrix &rhs,
     483  int numberRows, const int *whichRows,
     484  int numberColumns, const int *whichColumns)
     485  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
     486  , sumDualInfeasibilities_(0.0)
     487  , sumPrimalInfeasibilities_(0.0)
     488  , sumOfRelaxedDualInfeasibilities_(0.0)
     489  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     490  , start_(NULL)
     491  , end_(NULL)
     492  , lower_(NULL)
     493  , upper_(NULL)
     494  , backward_(NULL)
     495  , backToPivotRow_(NULL)
     496  , changeCost_(NULL)
     497  , keyVariable_(NULL)
     498  , next_(NULL)
     499  , toIndex_(NULL)
     500  , fromIndex_(NULL)
     501  , numberDualInfeasibilities_(0)
     502  , numberPrimalInfeasibilities_(0)
     503  , numberSets_(0)
     504  , saveNumber_(0)
     505  , possiblePivotKey_(0)
     506  , gubSlackIn_(-1)
     507  , firstGub_(0)
     508  , lastGub_(0)
     509  , gubType_(0)
    513510{
    514      setType(16);
     511  setType(16);
    515512}
    516513/* Return <code>x * A + y</code> in <code>z</code>.
    517514        Squashes small elements and knows about ClpSimplex */
    518 void
    519 ClpGubMatrix::transposeTimes(const ClpSimplex * model, double scalar,
    520                              const CoinIndexedVector * rowArray,
    521                              CoinIndexedVector * y,
    522                              CoinIndexedVector * columnArray) const
     515void ClpGubMatrix::transposeTimes(const ClpSimplex *model, double scalar,
     516  const CoinIndexedVector *rowArray,
     517  CoinIndexedVector *y,
     518  CoinIndexedVector *columnArray) const
    523519{
    524      columnArray->clear();
    525      double * pi = rowArray->denseVector();
    526      int numberNonZero = 0;
    527      int * index = columnArray->getIndices();
    528      double * array = columnArray->denseVector();
    529      int numberInRowArray = rowArray->getNumElements();
    530      // maybe I need one in OsiSimplex
    531      double zeroTolerance = model->zeroTolerance();
    532      int numberRows = model->numberRows();
    533      ClpPackedMatrix* rowCopy =
    534           dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
    535      bool packed = rowArray->packedMode();
    536      double factor = 0.3;
    537      // We may not want to do by row if there may be cache problems
    538      int numberColumns = model->numberColumns();
    539      // It would be nice to find L2 cache size - for moment 512K
    540      // Be slightly optimistic
    541      if (numberColumns * sizeof(double) > 1000000) {
    542           if (numberRows * 10 < numberColumns)
    543                factor = 0.1;
    544           else if (numberRows * 4 < numberColumns)
    545                factor = 0.15;
    546           else if (numberRows * 2 < numberColumns)
    547                factor = 0.2;
    548           //if (model->numberIterations()%50==0)
    549           //printf("%d nonzero\n",numberInRowArray);
    550      }
    551      // reduce for gub
    552      factor *= 0.5;
    553      assert (!y->getNumElements());
    554      if (numberInRowArray > factor * numberRows || !rowCopy) {
    555           // do by column
    556           int iColumn;
    557           // get matrix data pointers
    558           const int * row = matrix_->getIndices();
    559           const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    560           const int * columnLength = matrix_->getVectorLengths();
    561           const double * elementByColumn = matrix_->getElements();
    562           const double * rowScale = model->rowScale();
    563           int numberColumns = model->numberColumns();
    564           int iSet = -1;
    565           double djMod = 0.0;
    566           if (packed) {
    567                // need to expand pi into y
    568                assert(y->capacity() >= numberRows);
    569                double * piOld = pi;
    570                pi = y->denseVector();
    571                const int * whichRow = rowArray->getIndices();
    572                int i;
    573                if (!rowScale) {
    574                     // modify pi so can collapse to one loop
    575                     for (i = 0; i < numberInRowArray; i++) {
    576                          int iRow = whichRow[i];
    577                          pi[iRow] = scalar * piOld[i];
    578                     }
    579                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    580                          if (backward_[iColumn] != iSet) {
    581                               // get pi on gub row
    582                               iSet = backward_[iColumn];
    583                               if (iSet >= 0) {
    584                                    int iBasic = keyVariable_[iSet];
    585                                    if (iBasic < numberColumns) {
    586                                         // get dj without
    587                                         assert (model->getStatus(iBasic) == ClpSimplex::basic);
    588                                         djMod = 0.0;
    589                                         for (CoinBigIndex j = columnStart[iBasic];
    590                                                   j < columnStart[iBasic] + columnLength[iBasic]; j++) {
    591                                              int jRow = row[j];
    592                                              djMod -= pi[jRow] * elementByColumn[j];
    593                                         }
    594                                    } else {
    595                                         djMod = 0.0;
    596                                    }
    597                               } else {
    598                                    djMod = 0.0;
    599                               }
    600                          }
    601                          double value = -djMod;
    602                          CoinBigIndex j;
    603                          for (j = columnStart[iColumn];
    604                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    605                               int iRow = row[j];
    606                               value += pi[iRow] * elementByColumn[j];
    607                          }
    608                          if (fabs(value) > zeroTolerance) {
    609                               array[numberNonZero] = value;
    610                               index[numberNonZero++] = iColumn;
    611                          }
    612                     }
    613                } else {
    614                     // scaled
    615                     // modify pi so can collapse to one loop
    616                     for (i = 0; i < numberInRowArray; i++) {
    617                          int iRow = whichRow[i];
    618                          pi[iRow] = scalar * piOld[i] * rowScale[iRow];
    619                     }
    620                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    621                          if (backward_[iColumn] != iSet) {
    622                               // get pi on gub row
    623                               iSet = backward_[iColumn];
    624                               if (iSet >= 0) {
    625                                    int iBasic = keyVariable_[iSet];
    626                                    if (iBasic < numberColumns) {
    627                                         // get dj without
    628                                         assert (model->getStatus(iBasic) == ClpSimplex::basic);
    629                                         djMod = 0.0;
    630                                         // scaled
    631                                         for (CoinBigIndex j = columnStart[iBasic];
    632                                                   j < columnStart[iBasic] + columnLength[iBasic]; j++) {
    633                                              int jRow = row[j];
    634                                              djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow];
    635                                         }
    636                                    } else {
    637                                         djMod = 0.0;
    638                                    }
    639                               } else {
    640                                    djMod = 0.0;
    641                               }
    642                          }
    643                          double value = -djMod;
    644                          CoinBigIndex j;
    645                          const double * columnScale = model->columnScale();
    646                          for (j = columnStart[iColumn];
    647                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    648                               int iRow = row[j];
    649                               value += pi[iRow] * elementByColumn[j];
    650                          }
    651                          value *= columnScale[iColumn];
    652                          if (fabs(value) > zeroTolerance) {
    653                               array[numberNonZero] = value;
    654                               index[numberNonZero++] = iColumn;
    655                          }
    656                     }
    657                }
    658                // zero out
    659                for (i = 0; i < numberInRowArray; i++) {
    660                     int iRow = whichRow[i];
    661                     pi[iRow] = 0.0;
    662                }
    663           } else {
    664                // code later
    665                assert (packed);
    666                if (!rowScale) {
    667                     if (scalar == -1.0) {
    668                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    669                               double value = 0.0;
    670                               CoinBigIndex j;
    671                               for (j = columnStart[iColumn];
    672                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    673                                    int iRow = row[j];
    674                                    value += pi[iRow] * elementByColumn[j];
    675                               }
    676                               if (fabs(value) > zeroTolerance) {
    677                                    index[numberNonZero++] = iColumn;
    678                                    array[iColumn] = -value;
    679                               }
    680                          }
    681                     } else if (scalar == 1.0) {
    682                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    683                               double value = 0.0;
    684                               CoinBigIndex j;
    685                               for (j = columnStart[iColumn];
    686                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    687                                    int iRow = row[j];
    688                                    value += pi[iRow] * elementByColumn[j];
    689                               }
    690                               if (fabs(value) > zeroTolerance) {
    691                                    index[numberNonZero++] = iColumn;
    692                                    array[iColumn] = value;
    693                               }
    694                          }
    695                     } else {
    696                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    697                               double value = 0.0;
    698                               CoinBigIndex j;
    699                               for (j = columnStart[iColumn];
    700                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    701                                    int iRow = row[j];
    702                                    value += pi[iRow] * elementByColumn[j];
    703                               }
    704                               value *= scalar;
    705                               if (fabs(value) > zeroTolerance) {
    706                                    index[numberNonZero++] = iColumn;
    707                                    array[iColumn] = value;
    708                               }
    709                          }
    710                     }
    711                } else {
    712                     // scaled
    713                     if (scalar == -1.0) {
    714                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    715                               double value = 0.0;
    716                               CoinBigIndex j;
    717                               const double * columnScale = model->columnScale();
    718                               for (j = columnStart[iColumn];
    719                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    720                                    int iRow = row[j];
    721                                    value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    722                               }
    723                               value *= columnScale[iColumn];
    724                               if (fabs(value) > zeroTolerance) {
    725                                    index[numberNonZero++] = iColumn;
    726                                    array[iColumn] = -value;
    727                               }
    728                          }
    729                     } else if (scalar == 1.0) {
    730                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    731                               double value = 0.0;
    732                               CoinBigIndex j;
    733                               const double * columnScale = model->columnScale();
    734                               for (j = columnStart[iColumn];
    735                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    736                                    int iRow = row[j];
    737                                    value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    738                               }
    739                               value *= columnScale[iColumn];
    740                               if (fabs(value) > zeroTolerance) {
    741                                    index[numberNonZero++] = iColumn;
    742                                    array[iColumn] = value;
    743                               }
    744                          }
    745                     } else {
    746                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    747                               double value = 0.0;
    748                               CoinBigIndex j;
    749                               const double * columnScale = model->columnScale();
    750                               for (j = columnStart[iColumn];
    751                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    752                                    int iRow = row[j];
    753                                    value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    754                               }
    755                               value *= scalar * columnScale[iColumn];
    756                               if (fabs(value) > zeroTolerance) {
    757                                    index[numberNonZero++] = iColumn;
    758                                    array[iColumn] = value;
    759                               }
    760                          }
    761                     }
    762                }
    763           }
    764           columnArray->setNumElements(numberNonZero);
    765           y->setNumElements(0);
    766      } else {
    767           // do by row
    768           transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    769      }
    770      if (packed)
    771           columnArray->setPackedMode(true);
    772      if (0) {
    773           columnArray->checkClean();
    774           int numberNonZero = columnArray->getNumElements();;
    775           int * index = columnArray->getIndices();
    776           double * array = columnArray->denseVector();
    777           int i;
    778           for (i = 0; i < numberNonZero; i++) {
    779                int j = index[i];
    780                double value;
    781                if (packed)
    782                     value = array[i];
    783                else
    784                     value = array[j];
    785                printf("Ti %d %d %g\n", i, j, value);
    786           }
    787      }
     520  columnArray->clear();
     521  double *pi = rowArray->denseVector();
     522  int numberNonZero = 0;
     523  int *index = columnArray->getIndices();
     524  double *array = columnArray->denseVector();
     525  int numberInRowArray = rowArray->getNumElements();
     526  // maybe I need one in OsiSimplex
     527  double zeroTolerance = model->zeroTolerance();
     528  int numberRows = model->numberRows();
     529  ClpPackedMatrix *rowCopy = dynamic_cast< ClpPackedMatrix * >(model->rowCopy());
     530  bool packed = rowArray->packedMode();
     531  double factor = 0.3;
     532  // We may not want to do by row if there may be cache problems
     533  int numberColumns = model->numberColumns();
     534  // It would be nice to find L2 cache size - for moment 512K
     535  // Be slightly optimistic
     536  if (numberColumns * sizeof(double) > 1000000) {
     537    if (numberRows * 10 < numberColumns)
     538      factor = 0.1;
     539    else if (numberRows * 4 < numberColumns)
     540      factor = 0.15;
     541    else if (numberRows * 2 < numberColumns)
     542      factor = 0.2;
     543    //if (model->numberIterations()%50==0)
     544    //printf("%d nonzero\n",numberInRowArray);
     545  }
     546  // reduce for gub
     547  factor *= 0.5;
     548  assert(!y->getNumElements());
     549  if (numberInRowArray > factor * numberRows || !rowCopy) {
     550    // do by column
     551    int iColumn;
     552    // get matrix data pointers
     553    const int *row = matrix_->getIndices();
     554    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     555    const int *columnLength = matrix_->getVectorLengths();
     556    const double *elementByColumn = matrix_->getElements();
     557    const double *rowScale = model->rowScale();
     558    int numberColumns = model->numberColumns();
     559    int iSet = -1;
     560    double djMod = 0.0;
     561    if (packed) {
     562      // need to expand pi into y
     563      assert(y->capacity() >= numberRows);
     564      double *piOld = pi;
     565      pi = y->denseVector();
     566      const int *whichRow = rowArray->getIndices();
     567      int i;
     568      if (!rowScale) {
     569        // modify pi so can collapse to one loop
     570        for (i = 0; i < numberInRowArray; i++) {
     571          int iRow = whichRow[i];
     572          pi[iRow] = scalar * piOld[i];
     573        }
     574        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     575          if (backward_[iColumn] != iSet) {
     576            // get pi on gub row
     577            iSet = backward_[iColumn];
     578            if (iSet >= 0) {
     579              int iBasic = keyVariable_[iSet];
     580              if (iBasic < numberColumns) {
     581                // get dj without
     582                assert(model->getStatus(iBasic) == ClpSimplex::basic);
     583                djMod = 0.0;
     584                for (CoinBigIndex j = columnStart[iBasic];
     585                     j < columnStart[iBasic] + columnLength[iBasic]; j++) {
     586                  int jRow = row[j];
     587                  djMod -= pi[jRow] * elementByColumn[j];
     588                }
     589              } else {
     590                djMod = 0.0;
     591              }
     592            } else {
     593              djMod = 0.0;
     594            }
     595          }
     596          double value = -djMod;
     597          CoinBigIndex j;
     598          for (j = columnStart[iColumn];
     599               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     600            int iRow = row[j];
     601            value += pi[iRow] * elementByColumn[j];
     602          }
     603          if (fabs(value) > zeroTolerance) {
     604            array[numberNonZero] = value;
     605            index[numberNonZero++] = iColumn;
     606          }
     607        }
     608      } else {
     609        // scaled
     610        // modify pi so can collapse to one loop
     611        for (i = 0; i < numberInRowArray; i++) {
     612          int iRow = whichRow[i];
     613          pi[iRow] = scalar * piOld[i] * rowScale[iRow];
     614        }
     615        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     616          if (backward_[iColumn] != iSet) {
     617            // get pi on gub row
     618            iSet = backward_[iColumn];
     619            if (iSet >= 0) {
     620              int iBasic = keyVariable_[iSet];
     621              if (iBasic < numberColumns) {
     622                // get dj without
     623                assert(model->getStatus(iBasic) == ClpSimplex::basic);
     624                djMod = 0.0;
     625                // scaled
     626                for (CoinBigIndex j = columnStart[iBasic];
     627                     j < columnStart[iBasic] + columnLength[iBasic]; j++) {
     628                  int jRow = row[j];
     629                  djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow];
     630                }
     631              } else {
     632                djMod = 0.0;
     633              }
     634            } else {
     635              djMod = 0.0;
     636            }
     637          }
     638          double value = -djMod;
     639          CoinBigIndex j;
     640          const double *columnScale = model->columnScale();
     641          for (j = columnStart[iColumn];
     642               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     643            int iRow = row[j];
     644            value += pi[iRow] * elementByColumn[j];
     645          }
     646          value *= columnScale[iColumn];
     647          if (fabs(value) > zeroTolerance) {
     648            array[numberNonZero] = value;
     649            index[numberNonZero++] = iColumn;
     650          }
     651        }
     652      }
     653      // zero out
     654      for (i = 0; i < numberInRowArray; i++) {
     655        int iRow = whichRow[i];
     656        pi[iRow] = 0.0;
     657      }
     658    } else {
     659      // code later
     660      assert(packed);
     661      if (!rowScale) {
     662        if (scalar == -1.0) {
     663          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     664            double value = 0.0;
     665            CoinBigIndex j;
     666            for (j = columnStart[iColumn];
     667                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     668              int iRow = row[j];
     669              value += pi[iRow] * elementByColumn[j];
     670            }
     671            if (fabs(value) > zeroTolerance) {
     672              index[numberNonZero++] = iColumn;
     673              array[iColumn] = -value;
     674            }
     675          }
     676        } else if (scalar == 1.0) {
     677          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     678            double value = 0.0;
     679            CoinBigIndex j;
     680            for (j = columnStart[iColumn];
     681                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     682              int iRow = row[j];
     683              value += pi[iRow] * elementByColumn[j];
     684            }
     685            if (fabs(value) > zeroTolerance) {
     686              index[numberNonZero++] = iColumn;
     687              array[iColumn] = value;
     688            }
     689          }
     690        } else {
     691          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     692            double value = 0.0;
     693            CoinBigIndex j;
     694            for (j = columnStart[iColumn];
     695                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     696              int iRow = row[j];
     697              value += pi[iRow] * elementByColumn[j];
     698            }
     699            value *= scalar;
     700            if (fabs(value) > zeroTolerance) {
     701              index[numberNonZero++] = iColumn;
     702              array[iColumn] = value;
     703            }
     704          }
     705        }
     706      } else {
     707        // scaled
     708        if (scalar == -1.0) {
     709          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     710            double value = 0.0;
     711            CoinBigIndex j;
     712            const double *columnScale = model->columnScale();
     713            for (j = columnStart[iColumn];
     714                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     715              int iRow = row[j];
     716              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     717            }
     718            value *= columnScale[iColumn];
     719            if (fabs(value) > zeroTolerance) {
     720              index[numberNonZero++] = iColumn;
     721              array[iColumn] = -value;
     722            }
     723          }
     724        } else if (scalar == 1.0) {
     725          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     726            double value = 0.0;
     727            CoinBigIndex j;
     728            const double *columnScale = model->columnScale();
     729            for (j = columnStart[iColumn];
     730                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     731              int iRow = row[j];
     732              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     733            }
     734            value *= columnScale[iColumn];
     735            if (fabs(value) > zeroTolerance) {
     736              index[numberNonZero++] = iColumn;
     737              array[iColumn] = value;
     738            }
     739          }
     740        } else {
     741          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     742            double value = 0.0;
     743            CoinBigIndex j;
     744            const double *columnScale = model->columnScale();
     745            for (j = columnStart[iColumn];
     746                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     747              int iRow = row[j];
     748              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     749            }
     750            value *= scalar * columnScale[iColumn];
     751            if (fabs(value) > zeroTolerance) {
     752              index[numberNonZero++] = iColumn;
     753              array[iColumn] = value;
     754            }
     755          }
     756        }
     757      }
     758    }
     759    columnArray->setNumElements(numberNonZero);
     760    y->setNumElements(0);
     761  } else {
     762    // do by row
     763    transposeTimesByRow(model, scalar, rowArray, y, columnArray);
     764  }
     765  if (packed)
     766    columnArray->setPackedMode(true);
     767  if (0) {
     768    columnArray->checkClean();
     769    int numberNonZero = columnArray->getNumElements();
     770    ;
     771    int *index = columnArray->getIndices();
     772    double *array = columnArray->denseVector();
     773    int i;
     774    for (i = 0; i < numberNonZero; i++) {
     775      int j = index[i];
     776      double value;
     777      if (packed)
     778        value = array[i];
     779      else
     780        value = array[j];
     781      printf("Ti %d %d %g\n", i, j, value);
     782    }
     783  }
    788784}
    789785/* Return <code>x * A + y</code> in <code>z</code>.
    790786        Squashes small elements and knows about ClpSimplex */
    791 void
    792 ClpGubMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
    793                                   const CoinIndexedVector * rowArray,
    794                                   CoinIndexedVector * y,
    795                                   CoinIndexedVector * columnArray) const
     787void ClpGubMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar,
     788  const CoinIndexedVector *rowArray,
     789  CoinIndexedVector *y,
     790  CoinIndexedVector *columnArray) const
    796791{
    797      // Do packed part
    798      ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    799      if (numberSets_) {
    800           /* what we need to do is do by row as normal but get list of sets touched
     792  // Do packed part
     793  ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray);
     794  if (numberSets_) {
     795    /* what we need to do is do by row as normal but get list of sets touched
    801796             and then update those ones */
    802           abort();
    803      }
     797    abort();
     798  }
    804799}
    805800/* Return <code>x *A in <code>z</code> but
    806801   just for indices in y. */
    807 void
    808 ClpGubMatrix::subsetTransposeTimes(const ClpSimplex * model,
    809                                    const CoinIndexedVector * rowArray,
    810                                    const CoinIndexedVector * y,
    811                                    CoinIndexedVector * columnArray) const
     802void ClpGubMatrix::subsetTransposeTimes(const ClpSimplex *model,
     803  const CoinIndexedVector *rowArray,
     804  const CoinIndexedVector *y,
     805  CoinIndexedVector *columnArray) const
    812806{
    813      columnArray->clear();
    814      double * pi = rowArray->denseVector();
    815      double * array = columnArray->denseVector();
    816      int jColumn;
    817      // get matrix data pointers
    818      const int * row = matrix_->getIndices();
    819      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    820      const int * columnLength = matrix_->getVectorLengths();
    821      const double * elementByColumn = matrix_->getElements();
    822      const double * rowScale = model->rowScale();
    823      int numberToDo = y->getNumElements();
    824      const int * which = y->getIndices();
    825      assert (!rowArray->packedMode());
    826      columnArray->setPacked();
    827      int numberTouched = 0;
    828      if (!rowScale) {
    829           for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    830                int iColumn = which[jColumn];
    831                double value = 0.0;
    832                CoinBigIndex j;
    833                for (j = columnStart[iColumn];
    834                          j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    835                     int iRow = row[j];
    836                     value += pi[iRow] * elementByColumn[j];
    837                }
    838                array[jColumn] = value;
    839                if (value) {
    840                     int iSet = backward_[iColumn];
    841                     if (iSet >= 0) {
    842                          int iBasic = keyVariable_[iSet];
    843                          if (iBasic == iColumn) {
    844                               toIndex_[iSet] = jColumn;
    845                               fromIndex_[numberTouched++] = iSet;
    846                          }
    847                     }
    848                }
    849           }
    850      } else {
    851           // scaled
    852           for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    853                int iColumn = which[jColumn];
    854                double value = 0.0;
    855                CoinBigIndex j;
    856                const double * columnScale = model->columnScale();
    857                for (j = columnStart[iColumn];
    858                          j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    859                     int iRow = row[j];
    860                     value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    861                }
    862                value *= columnScale[iColumn];
    863                array[jColumn] = value;
    864                if (value) {
    865                     int iSet = backward_[iColumn];
    866                     if (iSet >= 0) {
    867                          int iBasic = keyVariable_[iSet];
    868                          if (iBasic == iColumn) {
    869                               toIndex_[iSet] = jColumn;
    870                               fromIndex_[numberTouched++] = iSet;
    871                          }
    872                     }
    873                }
    874           }
    875      }
    876      // adjust djs
    877      for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    878           int iColumn = which[jColumn];
    879           int iSet = backward_[iColumn];
    880           if (iSet >= 0) {
    881                int kColumn = toIndex_[iSet];
    882                if (kColumn >= 0)
    883                     array[jColumn] -= array[kColumn];
    884           }
    885      }
    886      // and clear basic
    887      for (int j = 0; j < numberTouched; j++) {
    888           int iSet = fromIndex_[j];
    889           int kColumn = toIndex_[iSet];
    890           toIndex_[iSet] = -1;
    891           array[kColumn] = 0.0;
    892      }
     807  columnArray->clear();
     808  double *pi = rowArray->denseVector();
     809  double *array = columnArray->denseVector();
     810  int jColumn;
     811  // get matrix data pointers
     812  const int *row = matrix_->getIndices();
     813  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     814  const int *columnLength = matrix_->getVectorLengths();
     815  const double *elementByColumn = matrix_->getElements();
     816  const double *rowScale = model->rowScale();
     817  int numberToDo = y->getNumElements();
     818  const int *which = y->getIndices();
     819  assert(!rowArray->packedMode());
     820  columnArray->setPacked();
     821  int numberTouched = 0;
     822  if (!rowScale) {
     823    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     824      int iColumn = which[jColumn];
     825      double value = 0.0;
     826      CoinBigIndex j;
     827      for (j = columnStart[iColumn];
     828           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     829        int iRow = row[j];
     830        value += pi[iRow] * elementByColumn[j];
     831      }
     832      array[jColumn] = value;
     833      if (value) {
     834        int iSet = backward_[iColumn];
     835        if (iSet >= 0) {
     836          int iBasic = keyVariable_[iSet];
     837          if (iBasic == iColumn) {
     838            toIndex_[iSet] = jColumn;
     839            fromIndex_[numberTouched++] = iSet;
     840          }
     841        }
     842      }
     843    }
     844  } else {
     845    // scaled
     846    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     847      int iColumn = which[jColumn];
     848      double value = 0.0;
     849      CoinBigIndex j;
     850      const double *columnScale = model->columnScale();
     851      for (j = columnStart[iColumn];
     852           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     853        int iRow = row[j];
     854        value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     855      }
     856      value *= columnScale[iColumn];
     857      array[jColumn] = value;
     858      if (value) {
     859        int iSet = backward_[iColumn];
     860        if (iSet >= 0) {
     861          int iBasic = keyVariable_[iSet];
     862          if (iBasic == iColumn) {
     863            toIndex_[iSet] = jColumn;
     864            fromIndex_[numberTouched++] = iSet;
     865          }
     866        }
     867      }
     868    }
     869  }
     870  // adjust djs
     871  for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     872    int iColumn = which[jColumn];
     873    int iSet = backward_[iColumn];
     874    if (iSet >= 0) {
     875      int kColumn = toIndex_[iSet];
     876      if (kColumn >= 0)
     877        array[jColumn] -= array[kColumn];
     878    }
     879  }
     880  // and clear basic
     881  for (int j = 0; j < numberTouched; j++) {
     882    int iSet = fromIndex_[j];
     883    int kColumn = toIndex_[iSet];
     884    toIndex_[iSet] = -1;
     885    array[kColumn] = 0.0;
     886  }
    893887}
    894888/// returns number of elements in column part of basis,
    895 int
    896 ClpGubMatrix::countBasis(const int * whichColumn,
    897                          int & numberColumnBasic)
     889int ClpGubMatrix::countBasis(const int *whichColumn,
     890  int &numberColumnBasic)
    898891{
    899      int i;
    900      int numberColumns = getNumCols();
    901      const int * columnLength = matrix_->getVectorLengths();
    902      int numberRows = getNumRows();
    903      int numberBasic = 0;
    904      CoinBigIndex numberElements = 0;
    905      int lastSet = -1;
    906      int key = -1;
    907      int keyLength = -1;
    908      double * work = new double[numberRows];
    909      CoinZeroN(work, numberRows);
    910      char * mark = new char[numberRows];
    911      CoinZeroN(mark, numberRows);
    912      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    913      const int * row = matrix_->getIndices();
    914      const double * elementByColumn = matrix_->getElements();
    915      //ClpGubDynamicMatrix* gubx =
    916      //dynamic_cast< ClpGubDynamicMatrix*>(this);
    917      //int * id = gubx->id();
    918      // just count
    919      for (i = 0; i < numberColumnBasic; i++) {
    920           int iColumn = whichColumn[i];
    921           int iSet = backward_[iColumn];
    922           int length = columnLength[iColumn];
    923           if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
    924                numberElements += length;
    925                numberBasic++;
    926                //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length);
     892  int i;
     893  int numberColumns = getNumCols();
     894  const int *columnLength = matrix_->getVectorLengths();
     895  int numberRows = getNumRows();
     896  int numberBasic = 0;
     897  CoinBigIndex numberElements = 0;
     898  int lastSet = -1;
     899  int key = -1;
     900  int keyLength = -1;
     901  double *work = new double[numberRows];
     902  CoinZeroN(work, numberRows);
     903  char *mark = new char[numberRows];
     904  CoinZeroN(mark, numberRows);
     905  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     906  const int *row = matrix_->getIndices();
     907  const double *elementByColumn = matrix_->getElements();
     908  //ClpGubDynamicMatrix* gubx =
     909  //dynamic_cast< ClpGubDynamicMatrix*>(this);
     910  //int * id = gubx->id();
     911  // just count
     912  for (i = 0; i < numberColumnBasic; i++) {
     913    int iColumn = whichColumn[i];
     914    int iSet = backward_[iColumn];
     915    int length = columnLength[iColumn];
     916    if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
     917      numberElements += length;
     918      numberBasic++;
     919      //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length);
     920    } else {
     921      // in gub set
     922      if (iColumn != keyVariable_[iSet]) {
     923        numberBasic++;
     924        CoinBigIndex j;
     925        // not key
     926        if (lastSet < iSet) {
     927          // erase work
     928          if (key >= 0) {
     929            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
     930              work[row[j]] = 0.0;
     931          }
     932          key = keyVariable_[iSet];
     933          lastSet = iSet;
     934          keyLength = columnLength[key];
     935          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
     936            work[row[j]] = elementByColumn[j];
     937        }
     938        int extra = keyLength;
     939        for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
     940          int iRow = row[j];
     941          double keyValue = work[iRow];
     942          double value = elementByColumn[j];
     943          if (!keyValue) {
     944            if (fabs(value) > 1.0e-20)
     945              extra++;
    927946          } else {
    928                // in gub set
    929                if (iColumn != keyVariable_[iSet]) {
    930                     numberBasic++;
    931                     CoinBigIndex j;
    932                     // not key
    933                     if (lastSet < iSet) {
    934                          // erase work
    935                          if (key >= 0) {
    936                               for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
    937                                    work[row[j]] = 0.0;
    938                          }
    939                          key = keyVariable_[iSet];
    940                          lastSet = iSet;
    941                          keyLength = columnLength[key];
    942                          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
    943                               work[row[j]] = elementByColumn[j];
    944                     }
    945                     int extra = keyLength;
    946                     for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
    947                          int iRow = row[j];
    948                          double keyValue = work[iRow];
    949                          double value = elementByColumn[j];
    950                          if (!keyValue) {
    951                               if (fabs(value) > 1.0e-20)
    952                                    extra++;
    953                          } else {
    954                               value -= keyValue;
    955                               if (fabs(value) <= 1.0e-20)
    956                                    extra--;
    957                          }
    958                     }
    959                     numberElements += extra;
    960                     //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra);
    961                }
    962           }
    963      }
    964      delete [] work;
    965      delete [] mark;
    966      // update number of column basic
    967      numberColumnBasic = numberBasic;
     947            value -= keyValue;
     948            if (fabs(value) <= 1.0e-20)
     949              extra--;
     950          }
     951        }
     952        numberElements += extra;
     953        //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra);
     954      }
     955    }
     956  }
     957  delete[] work;
     958  delete[] mark;
     959  // update number of column basic
     960  numberColumnBasic = numberBasic;
    968961#if COIN_BIG_INDEX
    969      if (numberElements>COIN_INT_MAX) {
    970        printf("Factorization too large\n");
    971        abort();
    972      }
     962  if (numberElements > COIN_INT_MAX) {
     963    printf("Factorization too large\n");
     964    abort();
     965  }
    973966#endif
    974      return static_cast<int>(numberElements);
     967  return static_cast< int >(numberElements);
    975968}
    976 void
    977 ClpGubMatrix::fillBasis(ClpSimplex * model,
    978                         const int * whichColumn,
    979                         int & numberColumnBasic,
    980                         int * indexRowU, int * start,
    981                         int * rowCount, int * columnCount,
    982                         CoinFactorizationDouble * elementU)
     969void ClpGubMatrix::fillBasis(ClpSimplex *model,
     970  const int *whichColumn,
     971  int &numberColumnBasic,
     972  int *indexRowU, int *start,
     973  int *rowCount, int *columnCount,
     974  CoinFactorizationDouble *elementU)
    983975{
    984      int i;
    985      int numberColumns = getNumCols();
    986      const int * columnLength = matrix_->getVectorLengths();
    987      int numberRows = getNumRows();
    988      assert (next_ || !elementU) ;
    989      CoinBigIndex numberElements = start[0];
    990      int lastSet = -1;
    991      int key = -1;
    992      int keyLength = -1;
    993      double * work = new double[numberRows];
    994      CoinZeroN(work, numberRows);
    995      char * mark = new char[numberRows];
    996      CoinZeroN(mark, numberRows);
    997      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    998      const int * row = matrix_->getIndices();
    999      const double * elementByColumn = matrix_->getElements();
    1000      const double * rowScale = model->rowScale();
    1001      int numberBasic = 0;
    1002      if (0) {
    1003           printf("%d basiccolumns\n", numberColumnBasic);
    1004           int i;
    1005           for (i = 0; i < numberSets_; i++) {
    1006                int k = keyVariable_[i];
    1007                if (k < numberColumns) {
    1008                     printf("key %d on set %d, %d elements\n", k, i, columnStart[k+1] - columnStart[k]);
    1009                     for (int j = columnStart[k]; j < columnStart[k+1]; j++)
    1010                          printf("row %d el %g\n", row[j], elementByColumn[j]);
    1011                } else {
    1012                     printf("slack key on set %d\n", i);
    1013                }
    1014           }
    1015      }
    1016      // fill
    1017      if (!rowScale) {
    1018           // no scaling
    1019           for (i = 0; i < numberColumnBasic; i++) {
    1020                int iColumn = whichColumn[i];
    1021                int iSet = backward_[iColumn];
    1022                int length = columnLength[iColumn];
    1023                if (0) {
    1024                     int k = iColumn;
    1025                     printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k+1] - columnStart[k]);
    1026                     for (int j = columnStart[k]; j < columnStart[k+1]; j++)
    1027                          printf("row %d el %g\n", row[j], elementByColumn[j]);
    1028                }
    1029                CoinBigIndex j;
    1030                if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
    1031                     for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1032                          double value = elementByColumn[j];
    1033                          if (fabs(value) > 1.0e-20) {
    1034                               int iRow = row[j];
    1035                               indexRowU[numberElements] = iRow;
    1036                               rowCount[iRow]++;
    1037                               elementU[numberElements++] = value;
    1038                          }
    1039                     }
    1040                     // end of column
    1041                     columnCount[numberBasic] = numberElements - start[numberBasic];
    1042                     numberBasic++;
    1043                     start[numberBasic] = numberElements;
    1044                } else {
    1045                     // in gub set
    1046                     if (iColumn != keyVariable_[iSet]) {
    1047                          // not key
    1048                          if (lastSet != iSet) {
    1049                               // erase work
    1050                               if (key >= 0) {
    1051                                    for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
    1052                                         int iRow = row[j];
    1053                                         work[iRow] = 0.0;
    1054                                         mark[iRow] = 0;
    1055                                    }
    1056                               }
    1057                               key = keyVariable_[iSet];
    1058                               lastSet = iSet;
    1059                               keyLength = columnLength[key];
    1060                               for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
    1061                                    int iRow = row[j];
    1062                                    work[iRow] = elementByColumn[j];
    1063                                    mark[iRow] = 1;
    1064                               }
    1065                          }
    1066                          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
    1067                               int iRow = row[j];
    1068                               double value = elementByColumn[j];
    1069                               if (mark[iRow]) {
    1070                                    mark[iRow] = 0;
    1071                                    double keyValue = work[iRow];
    1072                                    value -= keyValue;
    1073                               }
    1074                               if (fabs(value) > 1.0e-20) {
    1075                                    indexRowU[numberElements] = iRow;
    1076                                    rowCount[iRow]++;
    1077                                    elementU[numberElements++] = value;
    1078                               }
    1079                          }
    1080                          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
    1081                               int iRow = row[j];
    1082                               if (mark[iRow]) {
    1083                                    double value = -work[iRow];
    1084                                    if (fabs(value) > 1.0e-20) {
    1085                                         indexRowU[numberElements] = iRow;
    1086                                         rowCount[iRow]++;
    1087                                         elementU[numberElements++] = value;
    1088                                    }
    1089                               } else {
    1090                                    // just put back mark
    1091                                    mark[iRow] = 1;
    1092                               }
    1093                          }
    1094                          // end of column
    1095                          columnCount[numberBasic] = numberElements - start[numberBasic];
    1096                          numberBasic++;
    1097                          start[numberBasic] = numberElements;
    1098                     }
    1099                }
    1100           }
    1101      } else {
    1102           // scaling
    1103           const double * columnScale = model->columnScale();
    1104           for (i = 0; i < numberColumnBasic; i++) {
    1105                int iColumn = whichColumn[i];
    1106                int iSet = backward_[iColumn];
    1107                int length = columnLength[iColumn];
    1108                CoinBigIndex j;
    1109                if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
    1110                     double scale = columnScale[iColumn];
    1111                     for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1112                          int iRow = row[j];
    1113                          double value = elementByColumn[j] * scale * rowScale[iRow];
    1114                          if (fabs(value) > 1.0e-20) {
    1115                               indexRowU[numberElements] = iRow;
    1116                               rowCount[iRow]++;
    1117                               elementU[numberElements++] = value;
    1118                          }
    1119                     }
    1120                     // end of column
    1121                     columnCount[numberBasic] = numberElements - start[numberBasic];
    1122                     numberBasic++;
    1123                     start[numberBasic] = numberElements;
    1124                } else {
    1125                     // in gub set
    1126                     if (iColumn != keyVariable_[iSet]) {
    1127                          double scale = columnScale[iColumn];
    1128                          // not key
    1129                          if (lastSet < iSet) {
    1130                               // erase work
    1131                               if (key >= 0) {
    1132                                    for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
    1133                                         int iRow = row[j];
    1134                                         work[iRow] = 0.0;
    1135                                         mark[iRow] = 0;
    1136                                    }
    1137                               }
    1138                               key = keyVariable_[iSet];
    1139                               lastSet = iSet;
    1140                               keyLength = columnLength[key];
    1141                               double scale = columnScale[key];
    1142                               for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
    1143                                    int iRow = row[j];
    1144                                    work[iRow] = elementByColumn[j] * scale * rowScale[iRow];
    1145                                    mark[iRow] = 1;
    1146                               }
    1147                          }
    1148                          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
    1149                               int iRow = row[j];
    1150                               double value = elementByColumn[j] * scale * rowScale[iRow];
    1151                               if (mark[iRow]) {
    1152                                    mark[iRow] = 0;
    1153                                    double keyValue = work[iRow];
    1154                                    value -= keyValue;
    1155                               }
    1156                               if (fabs(value) > 1.0e-20) {
    1157                                    indexRowU[numberElements] = iRow;
    1158                                    rowCount[iRow]++;
    1159                                    elementU[numberElements++] = value;
    1160                               }
    1161                          }
    1162                          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
    1163                               int iRow = row[j];
    1164                               if (mark[iRow]) {
    1165                                    double value = -work[iRow];
    1166                                    if (fabs(value) > 1.0e-20) {
    1167                                         indexRowU[numberElements] = iRow;
    1168                                         rowCount[iRow]++;
    1169                                         elementU[numberElements++] = value;
    1170                                    }
    1171                               } else {
    1172                                    // just put back mark
    1173                                    mark[iRow] = 1;
    1174                               }
    1175                          }
    1176                          // end of column
    1177                          columnCount[numberBasic] = numberElements - start[numberBasic];
    1178                          numberBasic++;
    1179                          start[numberBasic] = numberElements;
    1180                     }
    1181                }
    1182           }
    1183      }
    1184      delete [] work;
    1185      delete [] mark;
    1186      // update number of column basic
    1187      numberColumnBasic = numberBasic;
     976  int i;
     977  int numberColumns = getNumCols();
     978  const int *columnLength = matrix_->getVectorLengths();
     979  int numberRows = getNumRows();
     980  assert(next_ || !elementU);
     981  CoinBigIndex numberElements = start[0];
     982  int lastSet = -1;
     983  int key = -1;
     984  int keyLength = -1;
     985  double *work = new double[numberRows];
     986  CoinZeroN(work, numberRows);
     987  char *mark = new char[numberRows];
     988  CoinZeroN(mark, numberRows);
     989  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     990  const int *row = matrix_->getIndices();
     991  const double *elementByColumn = matrix_->getElements();
     992  const double *rowScale = model->rowScale();
     993  int numberBasic = 0;
     994  if (0) {
     995    printf("%d basiccolumns\n", numberColumnBasic);
     996    int i;
     997    for (i = 0; i < numberSets_; i++) {
     998      int k = keyVariable_[i];
     999      if (k < numberColumns) {
     1000        printf("key %d on set %d, %d elements\n", k, i, columnStart[k + 1] - columnStart[k]);
     1001        for (int j = columnStart[k]; j < columnStart[k + 1]; j++)
     1002          printf("row %d el %g\n", row[j], elementByColumn[j]);
     1003      } else {
     1004        printf("slack key on set %d\n", i);
     1005      }
     1006    }
     1007  }
     1008  // fill
     1009  if (!rowScale) {
     1010    // no scaling
     1011    for (i = 0; i < numberColumnBasic; i++) {
     1012      int iColumn = whichColumn[i];
     1013      int iSet = backward_[iColumn];
     1014      int length = columnLength[iColumn];
     1015      if (0) {
     1016        int k = iColumn;
     1017        printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k + 1] - columnStart[k]);
     1018        for (int j = columnStart[k]; j < columnStart[k + 1]; j++)
     1019          printf("row %d el %g\n", row[j], elementByColumn[j]);
     1020      }
     1021      CoinBigIndex j;
     1022      if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
     1023        for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1024          double value = elementByColumn[j];
     1025          if (fabs(value) > 1.0e-20) {
     1026            int iRow = row[j];
     1027            indexRowU[numberElements] = iRow;
     1028            rowCount[iRow]++;
     1029            elementU[numberElements++] = value;
     1030          }
     1031        }
     1032        // end of column
     1033        columnCount[numberBasic] = numberElements - start[numberBasic];
     1034        numberBasic++;
     1035        start[numberBasic] = numberElements;
     1036      } else {
     1037        // in gub set
     1038        if (iColumn != keyVariable_[iSet]) {
     1039          // not key
     1040          if (lastSet != iSet) {
     1041            // erase work
     1042            if (key >= 0) {
     1043              for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
     1044                int iRow = row[j];
     1045                work[iRow] = 0.0;
     1046                mark[iRow] = 0;
     1047              }
     1048            }
     1049            key = keyVariable_[iSet];
     1050            lastSet = iSet;
     1051            keyLength = columnLength[key];
     1052            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
     1053              int iRow = row[j];
     1054              work[iRow] = elementByColumn[j];
     1055              mark[iRow] = 1;
     1056            }
     1057          }
     1058          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
     1059            int iRow = row[j];
     1060            double value = elementByColumn[j];
     1061            if (mark[iRow]) {
     1062              mark[iRow] = 0;
     1063              double keyValue = work[iRow];
     1064              value -= keyValue;
     1065            }
     1066            if (fabs(value) > 1.0e-20) {
     1067              indexRowU[numberElements] = iRow;
     1068              rowCount[iRow]++;
     1069              elementU[numberElements++] = value;
     1070            }
     1071          }
     1072          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
     1073            int iRow = row[j];
     1074            if (mark[iRow]) {
     1075              double value = -work[iRow];
     1076              if (fabs(value) > 1.0e-20) {
     1077                indexRowU[numberElements] = iRow;
     1078                rowCount[iRow]++;
     1079                elementU[numberElements++] = value;
     1080              }
     1081            } else {
     1082              // just put back mark
     1083              mark[iRow] = 1;
     1084            }
     1085          }
     1086          // end of column
     1087          columnCount[numberBasic] = numberElements - start[numberBasic];
     1088          numberBasic++;
     1089          start[numberBasic] = numberElements;
     1090        }
     1091      }
     1092    }
     1093  } else {
     1094    // scaling
     1095    const double *columnScale = model->columnScale();
     1096    for (i = 0; i < numberColumnBasic; i++) {
     1097      int iColumn = whichColumn[i];
     1098      int iSet = backward_[iColumn];
     1099      int length = columnLength[iColumn];
     1100      CoinBigIndex j;
     1101      if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
     1102        double scale = columnScale[iColumn];
     1103        for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1104          int iRow = row[j];
     1105          double value = elementByColumn[j] * scale * rowScale[iRow];
     1106          if (fabs(value) > 1.0e-20) {
     1107            indexRowU[numberElements] = iRow;
     1108            rowCount[iRow]++;
     1109            elementU[numberElements++] = value;
     1110          }
     1111        }
     1112        // end of column
     1113        columnCount[numberBasic] = numberElements - start[numberBasic];
     1114        numberBasic++;
     1115        start[numberBasic] = numberElements;
     1116      } else {
     1117        // in gub set
     1118        if (iColumn != keyVariable_[iSet]) {
     1119          double scale = columnScale[iColumn];
     1120          // not key
     1121          if (lastSet < iSet) {
     1122            // erase work
     1123            if (key >= 0) {
     1124              for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
     1125                int iRow = row[j];
     1126                work[iRow] = 0.0;
     1127                mark[iRow] = 0;
     1128              }
     1129            }
     1130            key = keyVariable_[iSet];
     1131            lastSet = iSet;
     1132            keyLength = columnLength[key];
     1133            double scale = columnScale[key];
     1134            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
     1135              int iRow = row[j];
     1136              work[iRow] = elementByColumn[j] * scale * rowScale[iRow];
     1137              mark[iRow] = 1;
     1138            }
     1139          }
     1140          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
     1141            int iRow = row[j];
     1142            double value = elementByColumn[j] * scale * rowScale[iRow];
     1143            if (mark[iRow]) {
     1144              mark[iRow] = 0;
     1145              double keyValue = work[iRow];
     1146              value -= keyValue;
     1147            }
     1148            if (fabs(value) > 1.0e-20) {
     1149              indexRowU[numberElements] = iRow;
     1150              rowCount[iRow]++;
     1151              elementU[numberElements++] = value;
     1152            }
     1153          }
     1154          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
     1155            int iRow = row[j];
     1156            if (mark[iRow]) {
     1157              double value = -work[iRow];
     1158              if (fabs(value) > 1.0e-20) {
     1159                indexRowU[numberElements] = iRow;
     1160                rowCount[iRow]++;
     1161                elementU[numberElements++] = value;
     1162              }
     1163            } else {
     1164              // just put back mark
     1165              mark[iRow] = 1;
     1166            }
     1167          }
     1168          // end of column
     1169          columnCount[numberBasic] = numberElements - start[numberBasic];
     1170          numberBasic++;
     1171          start[numberBasic] = numberElements;
     1172        }
     1173      }
     1174    }
     1175  }
     1176  delete[] work;
     1177  delete[] mark;
     1178  // update number of column basic
     1179  numberColumnBasic = numberBasic;
    11881180}
    11891181/* Unpacks a column into an CoinIndexedvector
    11901182 */
    1191 void
    1192 ClpGubMatrix::unpack(const ClpSimplex * model, CoinIndexedVector * rowArray,
    1193                      int iColumn) const
     1183void ClpGubMatrix::unpack(const ClpSimplex *model, CoinIndexedVector *rowArray,
     1184  int iColumn) const
    11941185{
    1195      assert (iColumn < model->numberColumns());
    1196      // Do packed part
    1197      ClpPackedMatrix::unpack(model, rowArray, iColumn);
    1198      int iSet = backward_[iColumn];
    1199      if (iSet >= 0) {
    1200           int iBasic = keyVariable_[iSet];
    1201           if (iBasic < model->numberColumns()) {
    1202                add(model, rowArray, iBasic, -1.0);
    1203           }
    1204      }
     1186  assert(iColumn < model->numberColumns());
     1187  // Do packed part
     1188  ClpPackedMatrix::unpack(model, rowArray, iColumn);
     1189  int iSet = backward_[iColumn];
     1190  if (iSet >= 0) {
     1191    int iBasic = keyVariable_[iSet];
     1192    if (iBasic < model->numberColumns()) {
     1193      add(model, rowArray, iBasic, -1.0);
     1194    }
     1195  }
    12051196}
    12061197/* Unpacks a column into a CoinIndexedVector
     
    12081199Note that model is NOT const.  Bounds and objective could
    12091200be modified if doing column generation (just for this variable) */
    1210 void
    1211 ClpGubMatrix::unpackPacked(ClpSimplex * model,
    1212                            CoinIndexedVector * rowArray,
    1213                            int iColumn) const
     1201void ClpGubMatrix::unpackPacked(ClpSimplex *model,
     1202  CoinIndexedVector *rowArray,
     1203  int iColumn) const
    12141204{
    1215      int numberColumns = model->numberColumns();
    1216      if (iColumn < numberColumns) {
    1217           // Do packed part
    1218           ClpPackedMatrix::unpackPacked(model, rowArray, iColumn);
    1219           int iSet = backward_[iColumn];
    1220           if (iSet >= 0) {
    1221                // columns are in order
    1222                int iBasic = keyVariable_[iSet];
    1223                if (iBasic < numberColumns) {
    1224                     int number = rowArray->getNumElements();
    1225                     const double * rowScale = model->rowScale();
    1226                     const int * row = matrix_->getIndices();
    1227                     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1228                     const int * columnLength = matrix_->getVectorLengths();
    1229                     const double * elementByColumn = matrix_->getElements();
    1230                     double * array = rowArray->denseVector();
    1231                     int * index = rowArray->getIndices();
    1232                     CoinBigIndex i;
    1233                     int numberOld = number;
    1234                     int lastIndex = 0;
    1235                     int next = index[lastIndex];
    1236                     if (!rowScale) {
    1237                          for (i = columnStart[iBasic];
    1238                                    i < columnStart[iBasic] + columnLength[iBasic]; i++) {
    1239                               int iRow = row[i];
    1240                               while (iRow > next) {
    1241                                    lastIndex++;
    1242                                    if (lastIndex == numberOld)
    1243                                         next = matrix_->getNumRows();
    1244                                    else
    1245                                         next = index[lastIndex];
    1246                               }
    1247                               if (iRow < next) {
    1248                                    array[number] = -elementByColumn[i];
    1249                                    index[number++] = iRow;
    1250                               } else {
    1251                                    assert (iRow == next);
    1252                                    array[lastIndex] -= elementByColumn[i];
    1253                                    if (!array[lastIndex])
    1254                                         array[lastIndex] = 1.0e-100;
    1255                               }
    1256                          }
    1257                     } else {
    1258                          // apply scaling
    1259                          double scale = model->columnScale()[iBasic];
    1260                          for (i = columnStart[iBasic];
    1261                                    i < columnStart[iBasic] + columnLength[iBasic]; i++) {
    1262                               int iRow = row[i];
    1263                               while (iRow > next) {
    1264                                    lastIndex++;
    1265                                    if (lastIndex == numberOld)
    1266                                         next = matrix_->getNumRows();
    1267                                    else
    1268                                         next = index[lastIndex];
    1269                               }
    1270                               if (iRow < next) {
    1271                                    array[number] = -elementByColumn[i] * scale * rowScale[iRow];
    1272                                    index[number++] = iRow;
    1273                               } else {
    1274                                    assert (iRow == next);
    1275                                    array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow];
    1276                                    if (!array[lastIndex])
    1277                                         array[lastIndex] = 1.0e-100;
    1278                               }
    1279                          }
    1280                     }
    1281                     rowArray->setNumElements(number);
    1282                }
    1283           }
    1284      } else {
    1285           // key slack entering
    1286           int iBasic = keyVariable_[gubSlackIn_];
    1287           assert (iBasic < numberColumns);
    1288           int number = 0;
    1289           const double * rowScale = model->rowScale();
    1290           const int * row = matrix_->getIndices();
    1291           const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1292           const int * columnLength = matrix_->getVectorLengths();
    1293           const double * elementByColumn = matrix_->getElements();
    1294           double * array = rowArray->denseVector();
    1295           int * index = rowArray->getIndices();
    1296           CoinBigIndex i;
    1297           if (!rowScale) {
    1298                for (i = columnStart[iBasic];
    1299                          i < columnStart[iBasic] + columnLength[iBasic]; i++) {
    1300                     int iRow = row[i];
    1301                     array[number] = elementByColumn[i];
    1302                     index[number++] = iRow;
    1303                }
    1304           } else {
    1305                // apply scaling
    1306                double scale = model->columnScale()[iBasic];
    1307                for (i = columnStart[iBasic];
    1308                          i < columnStart[iBasic] + columnLength[iBasic]; i++) {
    1309                     int iRow = row[i];
    1310                     array[number] = elementByColumn[i] * scale * rowScale[iRow];
    1311                     index[number++] = iRow;
    1312                }
    1313           }
    1314           rowArray->setNumElements(number);
    1315           rowArray->setPacked();
    1316      }
     1205  int numberColumns = model->numberColumns();
     1206  if (iColumn < numberColumns) {
     1207    // Do packed part
     1208    ClpPackedMatrix::unpackPacked(model, rowArray, iColumn);
     1209    int iSet = backward_[iColumn];
     1210    if (iSet >= 0) {
     1211      // columns are in order
     1212      int iBasic = keyVariable_[iSet];
     1213      if (iBasic < numberColumns) {
     1214        int number = rowArray->getNumElements();
     1215        const double *rowScale = model->rowScale();
     1216        const int *row = matrix_->getIndices();
     1217        const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     1218        const int *columnLength = matrix_->getVectorLengths();
     1219        const double *elementByColumn = matrix_->getElements();
     1220        double *array = rowArray->denseVector();
     1221        int *index = rowArray->getIndices();
     1222        CoinBigIndex i;
     1223        int numberOld = number;
     1224        int lastIndex = 0;
     1225        int next = index[lastIndex];
     1226        if (!rowScale) {
     1227          for (i = columnStart[iBasic];
     1228               i < columnStart[iBasic] + columnLength[iBasic]; i++) {
     1229            int iRow = row[i];
     1230            while (iRow > next) {
     1231              lastIndex++;
     1232              if (lastIndex == numberOld)
     1233                next = matrix_->getNumRows();
     1234              else
     1235                next = index[lastIndex];
     1236            }
     1237            if (iRow < next) {
     1238              array[number] = -elementByColumn[i];
     1239              index[number++] = iRow;
     1240            } else {
     1241              assert(iRow == next);
     1242              array[lastIndex] -= elementByColumn[i];
     1243              if (!array[lastIndex])
     1244                array[lastIndex] = 1.0e-100;
     1245            }
     1246          }
     1247        } else {
     1248          // apply scaling
     1249          double scale = model->columnScale()[iBasic];
     1250          for (i = columnStart[iBasic];
     1251               i < columnStart[iBasic] + columnLength[iBasic]; i++) {
     1252            int iRow = row[i];
     1253            while (iRow > next) {
     1254              lastIndex++;
     1255              if (lastIndex == numberOld)
     1256                next = matrix_->getNumRows();
     1257              else
     1258                next = index[lastIndex];
     1259            }
     1260            if (iRow < next) {
     1261              array[number] = -elementByColumn[i] * scale * rowScale[iRow];
     1262              index[number++] = iRow;
     1263            } else {
     1264              assert(iRow == next);
     1265              array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow];
     1266              if (!array[lastIndex])
     1267                array[lastIndex] = 1.0e-100;
     1268            }
     1269          }
     1270        }
     1271        rowArray->setNumElements(number);
     1272      }
     1273    }
     1274  } else {
     1275    // key slack entering
     1276    int iBasic = keyVariable_[gubSlackIn_];
     1277    assert(iBasic < numberColumns);
     1278    int number = 0;
     1279    const double *rowScale = model->rowScale();
     1280    const int *row = matrix_->getIndices();
     1281    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     1282    const int *columnLength = matrix_->getVectorLengths();
     1283    const double *elementByColumn = matrix_->getElements();
     1284    double *array = rowArray->denseVector();
     1285    int *index = rowArray->getIndices();
     1286    CoinBigIndex i;
     1287    if (!rowScale) {
     1288      for (i = columnStart[iBasic];
     1289           i < columnStart[iBasic] + columnLength[iBasic]; i++) {
     1290        int iRow = row[i];
     1291        array[number] = elementByColumn[i];
     1292        index[number++] = iRow;
     1293      }
     1294    } else {
     1295      // apply scaling
     1296      double scale = model->columnScale()[iBasic];
     1297      for (i = columnStart[iBasic];
     1298           i < columnStart[iBasic] + columnLength[iBasic]; i++) {
     1299        int iRow = row[i];
     1300        array[number] = elementByColumn[i] * scale * rowScale[iRow];
     1301        index[number++] = iRow;
     1302      }
     1303    }
     1304    rowArray->setNumElements(number);
     1305    rowArray->setPacked();
     1306  }
    13171307}
    13181308/* Adds multiple of a column into an CoinIndexedvector
    13191309      You can use quickAdd to add to vector */
    1320 void
    1321 ClpGubMatrix::add(const ClpSimplex * model, CoinIndexedVector * rowArray,
    1322                   int iColumn, double multiplier) const
     1310void ClpGubMatrix::add(const ClpSimplex *model, CoinIndexedVector *rowArray,
     1311  int iColumn, double multiplier) const
    13231312{
    1324      assert (iColumn < model->numberColumns());
    1325      // Do packed part
    1326      ClpPackedMatrix::add(model, rowArray, iColumn, multiplier);
    1327      int iSet = backward_[iColumn];
    1328      if (iSet >= 0 && iColumn != keyVariable_[iSet]) {
    1329           ClpPackedMatrix::add(model, rowArray, keyVariable_[iSet], -multiplier);
    1330      }
     1313  assert(iColumn < model->numberColumns());
     1314  // Do packed part
     1315  ClpPackedMatrix::add(model, rowArray, iColumn, multiplier);
     1316  int iSet = backward_[iColumn];
     1317  if (iSet >= 0 && iColumn != keyVariable_[iSet]) {
     1318    ClpPackedMatrix::add(model, rowArray, keyVariable_[iSet], -multiplier);
     1319  }
    13311320}
    13321321/* Adds multiple of a column into an array */
    1333 void
    1334 ClpGubMatrix::add(const ClpSimplex * model, double * array,
    1335                   int iColumn, double multiplier) const
     1322void ClpGubMatrix::add(const ClpSimplex *model, double *array,
     1323  int iColumn, double multiplier) const
    13361324{
    1337      assert (iColumn < model->numberColumns());
    1338      // Do packed part
    1339      ClpPackedMatrix::add(model, array, iColumn, multiplier);
    1340      if (iColumn < model->numberColumns()) {
    1341           int iSet = backward_[iColumn];
    1342           if (iSet >= 0 && iColumn != keyVariable_[iSet] && keyVariable_[iSet] < model->numberColumns()) {
    1343                ClpPackedMatrix::add(model, array, keyVariable_[iSet], -multiplier);
    1344           }
    1345      }
     1325  assert(iColumn < model->numberColumns());
     1326  // Do packed part
     1327  ClpPackedMatrix::add(model, array, iColumn, multiplier);
     1328  if (iColumn < model->numberColumns()) {
     1329    int iSet = backward_[iColumn];
     1330    if (iSet >= 0 && iColumn != keyVariable_[iSet] && keyVariable_[iSet] < model->numberColumns()) {
     1331      ClpPackedMatrix::add(model, array, keyVariable_[iSet], -multiplier);
     1332    }
     1333  }
    13461334}
    13471335// Partial pricing
    1348 void
    1349 ClpGubMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    1350                              int & bestSequence, int & numberWanted)
     1336void ClpGubMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
     1337  int &bestSequence, int &numberWanted)
    13511338{
    1352      numberWanted = currentWanted_;
    1353      if (numberSets_) {
    1354           // Do packed part before gub
    1355           int numberColumns = matrix_->getNumCols();
    1356           double ratio = static_cast<double> (firstGub_) /
    1357                          static_cast<double> (numberColumns);
    1358           ClpPackedMatrix::partialPricing(model, startFraction * ratio,
    1359                                           endFraction * ratio, bestSequence, numberWanted);
    1360           if (numberWanted || minimumGoodReducedCosts_ < -1) {
    1361                // do gub
    1362                const double * element = matrix_->getElements();
    1363                const int * row = matrix_->getIndices();
    1364                const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    1365                const int * length = matrix_->getVectorLengths();
    1366                const double * rowScale = model->rowScale();
    1367                const double * columnScale = model->columnScale();
    1368                int iSequence;
    1369                CoinBigIndex j;
    1370                double tolerance = model->currentDualTolerance();
    1371                double * reducedCost = model->djRegion();
    1372                const double * duals = model->dualRowSolution();
    1373                const double * cost = model->costRegion();
    1374                double bestDj;
    1375                int numberColumns = model->numberColumns();
    1376                int numberRows = model->numberRows();
    1377                if (bestSequence >= 0)
    1378                     bestDj = fabs(this->reducedCost(model, bestSequence));
    1379                else
    1380                     bestDj = tolerance;
    1381                int sequenceOut = model->sequenceOut();
    1382                int saveSequence = bestSequence;
    1383                int startG = firstGub_ + static_cast<int> (startFraction * (lastGub_ - firstGub_));
    1384                int endG = firstGub_ + static_cast<int> (endFraction * (lastGub_ - firstGub_));
    1385                endG = CoinMin(lastGub_, endG + 1);
    1386                // If nothing found yet can go all the way to end
    1387                int endAll = endG;
    1388                if (bestSequence < 0 && !startG)
    1389                     endAll = lastGub_;
    1390                int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
    1391                int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_;
    1392                int nSets = 0;
    1393                int iSet = -1;
    1394                double djMod = 0.0;
    1395                double infeasibilityCost = model->infeasibilityCost();
    1396                if (rowScale) {
    1397                     double bestDjMod = 0.0;
    1398                     // scaled
    1399                     for (iSequence = startG; iSequence < endAll; iSequence++) {
    1400                          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
    1401                               // give up
    1402                               numberWanted = 0;
    1403                               break;
    1404                          } else if (iSequence == endG && bestSequence >= 0) {
    1405                               break;
    1406                          }
    1407                          if (backward_[iSequence] != iSet) {
    1408                               // get pi on gub row
    1409                               iSet = backward_[iSequence];
    1410                               if (iSet >= 0) {
    1411                                    nSets++;
    1412                                    int iBasic = keyVariable_[iSet];
    1413                                    if (iBasic >= numberColumns) {
    1414                                         djMod = - weight(iSet) * infeasibilityCost;
    1415                                    } else {
    1416                                         // get dj without
    1417                                         assert (model->getStatus(iBasic) == ClpSimplex::basic);
    1418                                         djMod = 0.0;
    1419                                         // scaled
    1420                                         for (j = startColumn[iBasic];
    1421                                                   j < startColumn[iBasic] + length[iBasic]; j++) {
    1422                                              int jRow = row[j];
    1423                                              djMod -= duals[jRow] * element[j] * rowScale[jRow];
    1424                                         }
    1425                                         // allow for scaling
    1426                                         djMod +=  cost[iBasic] / columnScale[iBasic];
    1427                                         // See if gub slack possible - dj is djMod
    1428                                         if (getStatus(iSet) == ClpSimplex::atLowerBound) {
    1429                                              double value = -djMod;
    1430                                              if (value > tolerance) {
    1431                                                   numberWanted--;
    1432                                                   if (value > bestDj) {
    1433                                                        // check flagged variable and correct dj
    1434                                                        if (!flagged(iSet)) {
    1435                                                             bestDj = value;
    1436                                                             bestSequence = numberRows + numberColumns + iSet;
    1437                                                             bestDjMod = djMod;
    1438                                                        } else {
    1439                                                             // just to make sure we don't exit before got something
    1440                                                             numberWanted++;
    1441                                                             abort();
    1442                                                        }
    1443                                                   }
    1444                                              }
    1445                                         } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
    1446                                              double value = djMod;
    1447                                              if (value > tolerance) {
    1448                                                   numberWanted--;
    1449                                                   if (value > bestDj) {
    1450                                                        // check flagged variable and correct dj
    1451                                                        if (!flagged(iSet)) {
    1452                                                             bestDj = value;
    1453                                                             bestSequence = numberRows + numberColumns + iSet;
    1454                                                             bestDjMod = djMod;
    1455                                                        } else {
    1456                                                             // just to make sure we don't exit before got something
    1457                                                             numberWanted++;
    1458                                                             abort();
    1459                                                        }
    1460                                                   }
    1461                                              }
    1462                                         }
    1463                                    }
    1464                               } else {
    1465                                    // not in set
    1466                                    djMod = 0.0;
    1467                               }
    1468                          }
    1469                          if (iSequence != sequenceOut) {
    1470                               double value;
    1471                               ClpSimplex::Status status = model->getStatus(iSequence);
     1339  numberWanted = currentWanted_;
     1340  if (numberSets_) {
     1341    // Do packed part before gub
     1342    int numberColumns = matrix_->getNumCols();
     1343    double ratio = static_cast< double >(firstGub_) / static_cast< double >(numberColumns);
     1344    ClpPackedMatrix::partialPricing(model, startFraction * ratio,
     1345      endFraction * ratio, bestSequence, numberWanted);
     1346    if (numberWanted || minimumGoodReducedCosts_ < -1) {
     1347      // do gub
     1348      const double *element = matrix_->getElements();
     1349      const int *row = matrix_->getIndices();
     1350      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
     1351      const int *length = matrix_->getVectorLengths();
     1352      const double *rowScale = model->rowScale();
     1353      const double *columnScale = model->columnScale();
     1354      int iSequence;
     1355      CoinBigIndex j;
     1356      double tolerance = model->currentDualTolerance();
     1357      double *reducedCost = model->djRegion();
     1358      const double *duals = model->dualRowSolution();
     1359      const double *cost = model->costRegion();
     1360      double bestDj;
     1361      int numberColumns = model->numberColumns();
     1362      int numberRows = model->numberRows();
     1363      if (bestSequence >= 0)
     1364        bestDj = fabs(this->reducedCost(model, bestSequence));
     1365      else
     1366        bestDj = tolerance;
     1367      int sequenceOut = model->sequenceOut();
     1368      int saveSequence = bestSequence;
     1369      int startG = firstGub_ + static_cast< int >(startFraction * (lastGub_ - firstGub_));
     1370      int endG = firstGub_ + static_cast< int >(endFraction * (lastGub_ - firstGub_));
     1371      endG = CoinMin(lastGub_, endG + 1);
     1372      // If nothing found yet can go all the way to end
     1373      int endAll = endG;
     1374      if (bestSequence < 0 && !startG)
     1375        endAll = lastGub_;
     1376      int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
     1377      int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_;
     1378      int nSets = 0;
     1379      int iSet = -1;
     1380      double djMod = 0.0;
     1381      double infeasibilityCost = model->infeasibilityCost();
     1382      if (rowScale) {
     1383        double bestDjMod = 0.0;
     1384        // scaled
     1385        for (iSequence = startG; iSequence < endAll; iSequence++) {
     1386          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
     1387            // give up
     1388            numberWanted = 0;
     1389            break;
     1390          } else if (iSequence == endG && bestSequence >= 0) {
     1391            break;
     1392          }
     1393          if (backward_[iSequence] != iSet) {
     1394            // get pi on gub row
     1395            iSet = backward_[iSequence];
     1396            if (iSet >= 0) {
     1397              nSets++;
     1398              int iBasic = keyVariable_[iSet];
     1399              if (iBasic >= numberColumns) {
     1400                djMod = -weight(iSet) * infeasibilityCost;
     1401              } else {
     1402                // get dj without
     1403                assert(model->getStatus(iBasic) == ClpSimplex::basic);
     1404                djMod = 0.0;
     1405                // scaled
     1406                for (j = startColumn[iBasic];
     1407                     j < startColumn[iBasic] + length[iBasic]; j++) {
     1408                  int jRow = row[j];
     1409                  djMod -= duals[jRow] * element[j] * rowScale[jRow];
     1410                }
     1411                // allow for scaling
     1412                djMod += cost[iBasic] / columnScale[iBasic];
     1413                // See if gub slack possible - dj is djMod
     1414                if (getStatus(iSet) == ClpSimplex::atLowerBound) {
     1415                  double value = -djMod;
     1416                  if (value > tolerance) {
     1417                    numberWanted--;
     1418                    if (value > bestDj) {
     1419                      // check flagged variable and correct dj
     1420                      if (!flagged(iSet)) {
     1421                        bestDj = value;
     1422                        bestSequence = numberRows + numberColumns + iSet;
     1423                        bestDjMod = djMod;
     1424                      } else {
     1425                        // just to make sure we don't exit before got something
     1426                        numberWanted++;
     1427                        abort();
     1428                      }
     1429                    }
     1430                  }
     1431                } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
     1432                  double value = djMod;
     1433                  if (value > tolerance) {
     1434                    numberWanted--;
     1435                    if (value > bestDj) {
     1436                      // check flagged variable and correct dj
     1437                      if (!flagged(iSet)) {
     1438                        bestDj = value;
     1439                        bestSequence = numberRows + numberColumns + iSet;
     1440                        bestDjMod = djMod;
     1441                      } else {
     1442                        // just to make sure we don't exit before got something
     1443                        numberWanted++;
     1444                        abort();
     1445                      }
     1446                    }
     1447                  }
     1448                }
     1449              }
     1450            } else {
     1451              // not in set
     1452              djMod = 0.0;
     1453            }
     1454          }
     1455          if (iSequence != sequenceOut) {
     1456            double value;
     1457            ClpSimplex::Status status = model->getStatus(iSequence);
    14721458
    1473                               switch(status) {
     1459            switch (status) {
    14741460
    1475                               case ClpSimplex::basic:
    1476                               case ClpSimplex::isFixed:
    1477                                    break;
    1478                               case ClpSimplex::isFree:
    1479                               case ClpSimplex::superBasic:
    1480                                    value = -djMod;
    1481                                    // scaled
    1482                                    for (j = startColumn[iSequence];
    1483                                              j < startColumn[iSequence] + length[iSequence]; j++) {
    1484                                         int jRow = row[j];
    1485                                         value -= duals[jRow] * element[j] * rowScale[jRow];
    1486                                    }
    1487                                    value = fabs(cost[iSequence] + value * columnScale[iSequence]);
    1488                                    if (value > FREE_ACCEPT * tolerance) {
    1489                                         numberWanted--;
    1490                                         // we are going to bias towards free (but only if reasonable)
    1491                                         value *= FREE_BIAS;
    1492                                         if (value > bestDj) {
    1493                                              // check flagged variable and correct dj
    1494                                              if (!model->flagged(iSequence)) {
    1495                                                   bestDj = value;
    1496                                                   bestSequence = iSequence;
    1497                                                   bestDjMod = djMod;
    1498                                              } else {
    1499                                                   // just to make sure we don't exit before got something
    1500                                                   numberWanted++;
    1501                                              }
    1502                                         }
    1503                                    }
    1504                                    break;
    1505                               case ClpSimplex::atUpperBound:
    1506                                    value = -djMod;
    1507                                    // scaled
    1508                                    for (j = startColumn[iSequence];
    1509                                              j < startColumn[iSequence] + length[iSequence]; j++) {
    1510                                         int jRow = row[j];
    1511                                         value -= duals[jRow] * element[j] * rowScale[jRow];
    1512                                    }
    1513                                    value = cost[iSequence] + value * columnScale[iSequence];
    1514                                    if (value > tolerance) {
    1515                                         numberWanted--;
    1516                                         if (value > bestDj) {
    1517                                              // check flagged variable and correct dj
    1518                                              if (!model->flagged(iSequence)) {
    1519                                                   bestDj = value;
    1520                                                   bestSequence = iSequence;
    1521                                                   bestDjMod = djMod;
    1522                                              } else {
    1523                                                   // just to make sure we don't exit before got something
    1524                                                   numberWanted++;
    1525                                              }
    1526                                         }
    1527                                    }
    1528                                    break;
    1529                               case ClpSimplex::atLowerBound:
    1530                                    value = -djMod;
    1531                                    // scaled
    1532                                    for (j = startColumn[iSequence];
    1533                                              j < startColumn[iSequence] + length[iSequence]; j++) {
    1534                                         int jRow = row[j];
    1535                                         value -= duals[jRow] * element[j] * rowScale[jRow];
    1536                                    }
    1537                                    value = -(cost[iSequence] + value * columnScale[iSequence]);
    1538                                    if (value > tolerance) {
    1539                                         numberWanted--;
    1540                                         if (value > bestDj) {
    1541                                              // check flagged variable and correct dj
    1542                                              if (!model->flagged(iSequence)) {
    1543                                                   bestDj = value;
    1544                                                   bestSequence = iSequence;
    1545                                                   bestDjMod = djMod;
    1546                                              } else {
    1547                                                   // just to make sure we don't exit before got something
    1548                                                   numberWanted++;
    1549                                              }
    1550                                         }
    1551                                    }
    1552                                    break;
    1553                               }
    1554                          }
    1555                          if (!numberWanted)
    1556                               break;
     1461            case ClpSimplex::basic:
     1462            case ClpSimplex::isFixed:
     1463              break;
     1464            case ClpSimplex::isFree:
     1465            case ClpSimplex::superBasic:
     1466              value = -djMod;
     1467              // scaled
     1468              for (j = startColumn[iSequence];
     1469                   j < startColumn[iSequence] + length[iSequence]; j++) {
     1470                int jRow = row[j];
     1471                value -= duals[jRow] * element[j] * rowScale[jRow];
     1472              }
     1473              value = fabs(cost[iSequence] + value * columnScale[iSequence]);
     1474              if (value > FREE_ACCEPT * tolerance) {
     1475                numberWanted--;
     1476                // we are going to bias towards free (but only if reasonable)
     1477                value *= FREE_BIAS;
     1478                if (value > bestDj) {
     1479                  // check flagged variable and correct dj
     1480                  if (!model->flagged(iSequence)) {
     1481                    bestDj = value;
     1482                    bestSequence = iSequence;
     1483                    bestDjMod = djMod;
     1484                  } else {
     1485                    // just to make sure we don't exit before got something
     1486                    numberWanted++;
     1487                  }
     1488                }
     1489              }
     1490              break;
     1491            case ClpSimplex::atUpperBound:
     1492              value = -djMod;
     1493              // scaled
     1494              for (j = startColumn[iSequence];
     1495                   j < startColumn[iSequence] + length[iSequence]; j++) {
     1496                int jRow = row[j];
     1497                value -= duals[jRow] * element[j] * rowScale[jRow];
     1498              }
     1499              value = cost[iSequence] + value * columnScale[iSequence];
     1500              if (value > tolerance) {
     1501                numberWanted--;
     1502                if (value > bestDj) {
     1503                  // check flagged variable and correct dj
     1504                  if (!model->flagged(iSequence)) {
     1505                    bestDj = value;
     1506                    bestSequence = iSequence;
     1507                    bestDjMod = djMod;
     1508                  } else {
     1509                    // just to make sure we don't exit before got something
     1510                    numberWanted++;
     1511                  }
     1512                }
     1513              }
     1514              break;
     1515            case ClpSimplex::atLowerBound:
     1516              value = -djMod;
     1517              // scaled
     1518              for (j = startColumn[iSequence];
     1519                   j < startColumn[iSequence] + length[iSequence]; j++) {
     1520                int jRow = row[j];
     1521                value -= duals[jRow] * element[j] * rowScale[jRow];
     1522              }
     1523              value = -(cost[iSequence] + value * columnScale[iSequence]);
     1524              if (value > tolerance) {
     1525                numberWanted--;
     1526                if (value > bestDj) {
     1527                  // check flagged variable and correct dj
     1528                  if (!model->flagged(iSequence)) {
     1529                    bestDj = value;
     1530                    bestSequence = iSequence;
     1531                    bestDjMod = djMod;
     1532                  } else {
     1533                    // just to make sure we don't exit before got something
     1534                    numberWanted++;
     1535                  }
     1536                }
     1537              }
     1538              break;
     1539            }
     1540          }
     1541          if (!numberWanted)
     1542            break;
     1543        }
     1544        if (bestSequence != saveSequence) {
     1545          if (bestSequence < numberRows + numberColumns) {
     1546            // recompute dj
     1547            double value = bestDjMod;
     1548            // scaled
     1549            for (j = startColumn[bestSequence];
     1550                 j < startColumn[bestSequence] + length[bestSequence]; j++) {
     1551              int jRow = row[j];
     1552              value -= duals[jRow] * element[j] * rowScale[jRow];
     1553            }
     1554            reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
     1555            gubSlackIn_ = -1;
     1556          } else {
     1557            // slack - make last column
     1558            gubSlackIn_ = bestSequence - numberRows - numberColumns;
     1559            bestSequence = numberColumns + 2 * numberRows;
     1560            reducedCost[bestSequence] = bestDjMod;
     1561            model->setStatus(bestSequence, getStatus(gubSlackIn_));
     1562            if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
     1563              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
     1564            else
     1565              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
     1566            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
     1567            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
     1568            model->costRegion()[bestSequence] = 0.0;
     1569          }
     1570          savedBestSequence_ = bestSequence;
     1571          savedBestDj_ = reducedCost[savedBestSequence_];
     1572        }
     1573      } else {
     1574        double bestDjMod = 0.0;
     1575        //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
     1576        //     startG,endG,numberWanted);
     1577        for (iSequence = startG; iSequence < endG; iSequence++) {
     1578          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
     1579            // give up
     1580            numberWanted = 0;
     1581            break;
     1582          } else if (iSequence == endG && bestSequence >= 0) {
     1583            break;
     1584          }
     1585          if (backward_[iSequence] != iSet) {
     1586            // get pi on gub row
     1587            iSet = backward_[iSequence];
     1588            if (iSet >= 0) {
     1589              nSets++;
     1590              int iBasic = keyVariable_[iSet];
     1591              if (iBasic >= numberColumns) {
     1592                djMod = -weight(iSet) * infeasibilityCost;
     1593              } else {
     1594                // get dj without
     1595                assert(model->getStatus(iBasic) == ClpSimplex::basic);
     1596                djMod = 0.0;
     1597
     1598                for (j = startColumn[iBasic];
     1599                     j < startColumn[iBasic] + length[iBasic]; j++) {
     1600                  int jRow = row[j];
     1601                  djMod -= duals[jRow] * element[j];
     1602                }
     1603                djMod += cost[iBasic];
     1604                // See if gub slack possible - dj is djMod
     1605                if (getStatus(iSet) == ClpSimplex::atLowerBound) {
     1606                  double value = -djMod;
     1607                  if (value > tolerance) {
     1608                    numberWanted--;
     1609                    if (value > bestDj) {
     1610                      // check flagged variable and correct dj
     1611                      if (!flagged(iSet)) {
     1612                        bestDj = value;
     1613                        bestSequence = numberRows + numberColumns + iSet;
     1614                        bestDjMod = djMod;
     1615                      } else {
     1616                        // just to make sure we don't exit before got something
     1617                        numberWanted++;
     1618                        abort();
     1619                      }
    15571620                    }
    1558                     if (bestSequence != saveSequence) {
    1559                          if (bestSequence < numberRows + numberColumns) {
    1560                               // recompute dj
    1561                               double value = bestDjMod;
    1562                               // scaled
    1563                               for (j = startColumn[bestSequence];
    1564                                         j < startColumn[bestSequence] + length[bestSequence]; j++) {
    1565                                    int jRow = row[j];
    1566                                    value -= duals[jRow] * element[j] * rowScale[jRow];
    1567                               }
    1568                               reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
    1569                               gubSlackIn_ = -1;
    1570                          } else {
    1571                               // slack - make last column
    1572                               gubSlackIn_ = bestSequence - numberRows - numberColumns;
    1573                               bestSequence = numberColumns + 2 * numberRows;
    1574                               reducedCost[bestSequence] = bestDjMod;
    1575                               model->setStatus(bestSequence, getStatus(gubSlackIn_));
    1576                               if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
    1577                                    model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
    1578                               else
    1579                                    model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
    1580                               model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
    1581                               model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
    1582                               model->costRegion()[bestSequence] = 0.0;
    1583                          }
    1584                          savedBestSequence_ = bestSequence;
    1585                          savedBestDj_ = reducedCost[savedBestSequence_];
     1621                  }
     1622                } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
     1623                  double value = djMod;
     1624                  if (value > tolerance) {
     1625                    numberWanted--;
     1626                    if (value > bestDj) {
     1627                      // check flagged variable and correct dj
     1628                      if (!flagged(iSet)) {
     1629                        bestDj = value;
     1630                        bestSequence = numberRows + numberColumns + iSet;
     1631                        bestDjMod = djMod;
     1632                      } else {
     1633                        // just to make sure we don't exit before got something
     1634                        numberWanted++;
     1635                        abort();
     1636                      }
    15861637                    }
    1587                } else {
    1588                     double bestDjMod = 0.0;
    1589                     //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
    1590                     //     startG,endG,numberWanted);
    1591                     for (iSequence = startG; iSequence < endG; iSequence++) {
    1592                          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
    1593                               // give up
    1594                               numberWanted = 0;
    1595                               break;
    1596                          } else if (iSequence == endG && bestSequence >= 0) {
    1597                               break;
    1598                          }
    1599                          if (backward_[iSequence] != iSet) {
    1600                               // get pi on gub row
    1601                               iSet = backward_[iSequence];
    1602                               if (iSet >= 0) {
    1603                                    nSets++;
    1604                                    int iBasic = keyVariable_[iSet];
    1605                                    if (iBasic >= numberColumns) {
    1606                                         djMod = - weight(iSet) * infeasibilityCost;
    1607                                    } else {
    1608                                         // get dj without
    1609                                         assert (model->getStatus(iBasic) == ClpSimplex::basic);
    1610                                         djMod = 0.0;
     1638                  }
     1639                }
     1640              }
     1641            } else {
     1642              // not in set
     1643              djMod = 0.0;
     1644            }
     1645          }
     1646          if (iSequence != sequenceOut) {
     1647            double value;
     1648            ClpSimplex::Status status = model->getStatus(iSequence);
    16111649
    1612                                         for (j = startColumn[iBasic];
    1613                                                   j < startColumn[iBasic] + length[iBasic]; j++) {
    1614                                              int jRow = row[j];
    1615                                              djMod -= duals[jRow] * element[j];
    1616                                         }
    1617                                         djMod += cost[iBasic];
    1618                                         // See if gub slack possible - dj is djMod
    1619                                         if (getStatus(iSet) == ClpSimplex::atLowerBound) {
    1620                                              double value = -djMod;
    1621                                              if (value > tolerance) {
    1622                                                   numberWanted--;
    1623                                                   if (value > bestDj) {
    1624                                                        // check flagged variable and correct dj
    1625                                                        if (!flagged(iSet)) {
    1626                                                             bestDj = value;
    1627                                                             bestSequence = numberRows + numberColumns + iSet;
    1628                                                             bestDjMod = djMod;
    1629                                                        } else {
    1630                                                             // just to make sure we don't exit before got something
    1631                                                             numberWanted++;
    1632                                                             abort();
    1633                                                        }
    1634                                                   }
    1635                                              }
    1636                                         } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
    1637                                              double value = djMod;
    1638                                              if (value > tolerance) {
    1639                                                   numberWanted--;
    1640                                                   if (value > bestDj) {
    1641                                                        // check flagged variable and correct dj
    1642                                                        if (!flagged(iSet)) {
    1643                                                             bestDj = value;
    1644                                                             bestSequence = numberRows + numberColumns + iSet;
    1645                                                             bestDjMod = djMod;
    1646                                                        } else {
    1647                                                             // just to make sure we don't exit before got something
    1648                                                             numberWanted++;
    1649                                                             abort();
    1650                                                        }
    1651                                                   }
    1652                                              }
    1653                                         }
    1654                                    }
    1655                               } else {
    1656                                    // not in set
    1657                                    djMod = 0.0;
    1658                               }
    1659                          }
    1660                          if (iSequence != sequenceOut) {
    1661                               double value;
    1662                               ClpSimplex::Status status = model->getStatus(iSequence);
     1650            switch (status) {
    16631651
    1664                               switch(status) {
    1665 
    1666                               case ClpSimplex::basic:
    1667                               case ClpSimplex::isFixed:
    1668                                    break;
    1669                               case ClpSimplex::isFree:
    1670                               case ClpSimplex::superBasic:
    1671                                    value = cost[iSequence] - djMod;
    1672                                    for (j = startColumn[iSequence];
    1673                                              j < startColumn[iSequence] + length[iSequence]; j++) {
    1674                                         int jRow = row[j];
    1675                                         value -= duals[jRow] * element[j];
    1676                                    }
    1677                                    value = fabs(value);
    1678                                    if (value > FREE_ACCEPT * tolerance) {
    1679                                         numberWanted--;
    1680                                         // we are going to bias towards free (but only if reasonable)
    1681                                         value *= FREE_BIAS;
    1682                                         if (value > bestDj) {
    1683                                              // check flagged variable and correct dj
    1684                                              if (!model->flagged(iSequence)) {
    1685                                                   bestDj = value;
    1686                                                   bestSequence = iSequence;
    1687                                                   bestDjMod = djMod;
    1688                                              } else {
    1689                                                   // just to make sure we don't exit before got something
    1690                                                   numberWanted++;
    1691                                              }
    1692                                         }
    1693                                    }
    1694                                    break;
    1695                               case ClpSimplex::atUpperBound:
    1696                                    value = cost[iSequence] - djMod;
    1697                                    for (j = startColumn[iSequence];
    1698                                              j < startColumn[iSequence] + length[iSequence]; j++) {
    1699                                         int jRow = row[j];
    1700                                         value -= duals[jRow] * element[j];
    1701                                    }
    1702                                    if (value > tolerance) {
    1703                                         numberWanted--;
    1704                                         if (value > bestDj) {
    1705                                              // check flagged variable and correct dj
    1706                                              if (!model->flagged(iSequence)) {
    1707                                                   bestDj = value;
    1708                                                   bestSequence = iSequence;
    1709                                                   bestDjMod = djMod;
    1710                                              } else {
    1711                                                   // just to make sure we don't exit before got something
    1712                                                   numberWanted++;
    1713                                              }
    1714                                         }
    1715                                    }
    1716                                    break;
    1717                               case ClpSimplex::atLowerBound:
    1718                                    value = cost[iSequence] - djMod;
    1719                                    for (j = startColumn[iSequence];
    1720                                              j < startColumn[iSequence] + length[iSequence]; j++) {
    1721                                         int jRow = row[j];
    1722                                         value -= duals[jRow] * element[j];
    1723                                    }
    1724                                    value = -value;
    1725                                    if (value > tolerance) {
    1726                                         numberWanted--;
    1727                                         if (value > bestDj) {
    1728                                              // check flagged variable and correct dj
    1729                                              if (!model->flagged(iSequence)) {
    1730                                                   bestDj = value;
    1731                                                   bestSequence = iSequence;
    1732                                                   bestDjMod = djMod;
    1733                                              } else {
    1734                                                   // just to make sure we don't exit before got something
    1735                                                   numberWanted++;
    1736                                              }
    1737                                         }
    1738                                    }
    1739                                    break;
    1740                               }
    1741                          }
    1742                          if (!numberWanted)
    1743                               break;
    1744                     }
    1745                     if (bestSequence != saveSequence) {
    1746                          if (bestSequence < numberRows + numberColumns) {
    1747                               // recompute dj
    1748                               double value = cost[bestSequence] - bestDjMod;
    1749                               for (j = startColumn[bestSequence];
    1750                                         j < startColumn[bestSequence] + length[bestSequence]; j++) {
    1751                                    int jRow = row[j];
    1752                                    value -= duals[jRow] * element[j];
    1753                               }
    1754                               //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
    1755                               reducedCost[bestSequence] = value;
    1756                               gubSlackIn_ = -1;
    1757                          } else {
    1758                               // slack - make last column
    1759                               gubSlackIn_ = bestSequence - numberRows - numberColumns;
    1760                               bestSequence = numberColumns + 2 * numberRows;
    1761                               reducedCost[bestSequence] = bestDjMod;
    1762                               //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
    1763                               model->setStatus(bestSequence, getStatus(gubSlackIn_));
    1764                               if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
    1765                                    model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
    1766                               else
    1767                                    model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
    1768                               model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
    1769                               model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
    1770                               model->costRegion()[bestSequence] = 0.0;
    1771                          }
    1772                     }
    1773                }
    1774                // See if may be finished
    1775                if (startG == firstGub_ && bestSequence < 0)
    1776                     infeasibilityWeight_ = model_->infeasibilityCost();
    1777                else if (bestSequence >= 0)
    1778                     infeasibilityWeight_ = -1.0;
    1779           }
    1780           if (numberWanted) {
    1781                // Do packed part after gub
    1782                double offset = static_cast<double> (lastGub_) /
    1783                                static_cast<double> (numberColumns);
    1784                double ratio = static_cast<double> (numberColumns) /
    1785                               static_cast<double> (numberColumns) - offset;
    1786                double start2 = offset + ratio * startFraction;
    1787                double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6);
    1788                ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted);
    1789           }
    1790      } else {
    1791           // no gub
    1792           ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
    1793      }
    1794      if (bestSequence >= 0)
    1795           infeasibilityWeight_ = -1.0; // not optimal
    1796      currentWanted_ = numberWanted;
     1652            case ClpSimplex::basic:
     1653            case ClpSimplex::isFixed:
     1654              break;
     1655            case ClpSimplex::isFree:
     1656            case ClpSimplex::superBasic:
     1657              value = cost[iSequence] - djMod;
     1658              for (j = startColumn[iSequence];
     1659                   j < startColumn[iSequence] + length[iSequence]; j++) {
     1660                int jRow = row[j];
     1661                value -= duals[jRow] * element[j];
     1662              }
     1663              value = fabs(value);
     1664              if (value > FREE_ACCEPT * tolerance) {
     1665                numberWanted--;
     1666                // we are going to bias towards free (but only if reasonable)
     1667                value *= FREE_BIAS;
     1668                if (value > bestDj) {
     1669                  // check flagged variable and correct dj
     1670                  if (!model->flagged(iSequence)) {
     1671                    bestDj = value;
     1672                    bestSequence = iSequence;
     1673                    bestDjMod = djMod;
     1674                  } else {
     1675                    // just to make sure we don't exit before got something
     1676                    numberWanted++;
     1677                  }
     1678                }
     1679              }
     1680              break;
     1681            case ClpSimplex::atUpperBound:
     1682              value = cost[iSequence] - djMod;
     1683              for (j = startColumn[iSequence];
     1684                   j < startColumn[iSequence] + length[iSequence]; j++) {
     1685                int jRow = row[j];
     1686                value -= duals[jRow] * element[j];
     1687              }
     1688              if (value > tolerance) {
     1689                numberWanted--;
     1690                if (value > bestDj) {
     1691                  // check flagged variable and correct dj
     1692                  if (!model->flagged(iSequence)) {
     1693                    bestDj = value;
     1694                    bestSequence = iSequence;
     1695                    bestDjMod = djMod;
     1696                  } else {
     1697                    // just to make sure we don't exit before got something
     1698                    numberWanted++;
     1699                  }
     1700                }
     1701              }
     1702              break;
     1703            case ClpSimplex::atLowerBound:
     1704              value = cost[iSequence] - djMod;
     1705              for (j = startColumn[iSequence];
     1706                   j < startColumn[iSequence] + length[iSequence]; j++) {
     1707                int jRow = row[j];
     1708                value -= duals[jRow] * element[j];
     1709              }
     1710              value = -value;
     1711              if (value > tolerance) {
     1712                numberWanted--;
     1713                if (value > bestDj) {
     1714                  // check flagged variable and correct dj
     1715                  if (!model->flagged(iSequence)) {
     1716                    bestDj = value;
     1717                    bestSequence = iSequence;
     1718                    bestDjMod = djMod;
     1719                  } else {
     1720                    // just to make sure we don't exit before got something
     1721                    numberWanted++;
     1722                  }
     1723                }
     1724              }
     1725              break;
     1726            }
     1727          }
     1728          if (!numberWanted)
     1729            break;
     1730        }
     1731        if (bestSequence != saveSequence) {
     1732          if (bestSequence < numberRows + numberColumns) {
     1733            // recompute dj
     1734            double value = cost[bestSequence] - bestDjMod;
     1735            for (j = startColumn[bestSequence];
     1736                 j < startColumn[bestSequence] + length[bestSequence]; j++) {
     1737              int jRow = row[j];
     1738              value -= duals[jRow] * element[j];
     1739            }
     1740            //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
     1741            reducedCost[bestSequence] = value;
     1742            gubSlackIn_ = -1;
     1743          } else {
     1744            // slack - make last column
     1745            gubSlackIn_ = bestSequence - numberRows - numberColumns;
     1746            bestSequence = numberColumns + 2 * numberRows;
     1747            reducedCost[bestSequence] = bestDjMod;
     1748            //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
     1749            model->setStatus(bestSequence, getStatus(gubSlackIn_));
     1750            if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
     1751              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
     1752            else
     1753              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
     1754            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
     1755            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
     1756            model->costRegion()[bestSequence] = 0.0;
     1757          }
     1758        }
     1759      }
     1760      // See if may be finished
     1761      if (startG == firstGub_ && bestSequence < 0)
     1762        infeasibilityWeight_ = model_->infeasibilityCost();
     1763      else if (bestSequence >= 0)
     1764        infeasibilityWeight_ = -1.0;
     1765    }
     1766    if (numberWanted) {
     1767      // Do packed part after gub
     1768      double offset = static_cast< double >(lastGub_) / static_cast< double >(numberColumns);
     1769      double ratio = static_cast< double >(numberColumns) / static_cast< double >(numberColumns) - offset;
     1770      double start2 = offset + ratio * startFraction;
     1771      double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6);
     1772      ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted);
     1773    }
     1774  } else {
     1775    // no gub
     1776    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
     1777  }
     1778  if (bestSequence >= 0)
     1779    infeasibilityWeight_ = -1.0; // not optimal
     1780  currentWanted_ = numberWanted;
    17971781}
    17981782/* expands an updated column to allow for extra rows which the main
    17991783   solver does not know about and returns number added.
    18001784*/
    1801 int
    1802 ClpGubMatrix::extendUpdated(ClpSimplex * model, CoinIndexedVector * update, int mode)
     1785int ClpGubMatrix::extendUpdated(ClpSimplex *model, CoinIndexedVector *update, int mode)
    18031786{
    1804      // I think we only need to bother about sets with two in basis or incoming set
    1805      int number = update->getNumElements();
    1806      double * array = update->denseVector();
    1807      int * index = update->getIndices();
    1808      int i;
    1809      assert (!number || update->packedMode());
    1810      int * pivotVariable = model->pivotVariable();
    1811      int numberRows = model->numberRows();
    1812      int numberColumns = model->numberColumns();
    1813      int numberTotal = numberRows + numberColumns;
    1814      int sequenceIn = model->sequenceIn();
    1815      int returnCode = 0;
    1816      int iSetIn;
    1817      if (sequenceIn < numberColumns) {
    1818           iSetIn = backward_[sequenceIn];
    1819           gubSlackIn_ = -1; // in case set
    1820      } else if (sequenceIn < numberRows + numberColumns) {
    1821           iSetIn = -1;
    1822           gubSlackIn_ = -1; // in case set
    1823      } else {
    1824           iSetIn = gubSlackIn_;
    1825      }
    1826      double * lower = model->lowerRegion();
    1827      double * upper = model->upperRegion();
    1828      double * cost = model->costRegion();
    1829      double * solution = model->solutionRegion();
    1830      int number2 = number;
    1831      if (!mode) {
    1832           double primalTolerance = model->primalTolerance();
    1833           double infeasibilityCost = model->infeasibilityCost();
    1834           // extend
    1835           saveNumber_ = number;
    1836           for (i = 0; i < number; i++) {
    1837                int iRow = index[i];
    1838                int iPivot = pivotVariable[iRow];
    1839                if (iPivot < numberColumns) {
    1840                     int iSet = backward_[iPivot];
    1841                     if (iSet >= 0) {
    1842                          // two (or more) in set
    1843                          int iIndex = toIndex_[iSet];
    1844                          double otherValue = array[i];
    1845                          double value;
    1846                          if (iIndex < 0) {
    1847                               toIndex_[iSet] = number2;
    1848                               int iNew = number2 - number;
    1849                               fromIndex_[number2-number] = iSet;
    1850                               iIndex = number2;
    1851                               index[number2] = numberRows + iNew;
    1852                               // do key stuff
    1853                               int iKey = keyVariable_[iSet];
    1854                               if (iKey < numberColumns) {
    1855                                    // Save current cost of key
    1856                                    changeCost_[number2-number] = cost[iKey];
    1857                                    if (iSet != iSetIn)
    1858                                         value = 0.0;
    1859                                    else if (iSetIn != gubSlackIn_)
    1860                                         value = 1.0;
    1861                                    else
    1862                                         value = -1.0;
    1863                                    pivotVariable[numberRows+iNew] = iKey;
    1864                                    // Do I need to recompute?
    1865                                    double sol;
    1866                                    assert (getStatus(iSet) != ClpSimplex::basic);
    1867                                    if (getStatus(iSet) == ClpSimplex::atLowerBound)
    1868                                         sol = lower_[iSet];
    1869                                    else
    1870                                         sol = upper_[iSet];
    1871                                    if ((gubType_ & 8) != 0) {
    1872                                         int iColumn = next_[iKey];
    1873                                         // sum all non-key variables
    1874                                         while(iColumn >= 0) {
    1875                                              sol -= solution[iColumn];
    1876                                              iColumn = next_[iColumn];
    1877                                         }
    1878                                    } else {
    1879                                         int stop = -(iKey + 1);
    1880                                         int iColumn = next_[iKey];
    1881                                         // sum all non-key variables
    1882                                         while(iColumn != stop) {
    1883                                              if (iColumn < 0)
    1884                                                   iColumn = -iColumn - 1;
    1885                                              sol -= solution[iColumn];
    1886                                              iColumn = next_[iColumn];
    1887                                         }
    1888                                    }
    1889                                    solution[iKey] = sol;
    1890                                    if (model->algorithm() > 0)
    1891                                         model->nonLinearCost()->setOne(iKey, sol);
    1892                                    //assert (fabs(sol-solution[iKey])<1.0e-3);
    1893                               } else {
    1894                                    // gub slack is basic
    1895                                    // Save current cost of key
    1896                                    changeCost_[number2-number] = -weight(iSet) * infeasibilityCost;
    1897                                    otherValue = - otherValue; //allow for - sign on slack
    1898                                    if (iSet != iSetIn)
    1899                                         value = 0.0;
    1900                                    else
    1901                                         value = -1.0;
    1902                                    pivotVariable[numberRows+iNew] = iNew + numberTotal;
    1903                                    model->djRegion()[iNew+numberTotal] = 0.0;
    1904                                    double sol = 0.0;
    1905                                    if ((gubType_ & 8) != 0) {
    1906                                         int iColumn = next_[iKey];
    1907                                         // sum all non-key variables
    1908                                         while(iColumn >= 0) {
    1909                                              sol += solution[iColumn];
    1910                                              iColumn = next_[iColumn];
    1911                                         }
    1912                                    } else {
    1913                                         int stop = -(iKey + 1);
    1914                                         int iColumn = next_[iKey];
    1915                                         // sum all non-key variables
    1916                                         while(iColumn != stop) {
    1917                                              if (iColumn < 0)
    1918                                                   iColumn = -iColumn - 1;
    1919                                              sol += solution[iColumn];
    1920                                              iColumn = next_[iColumn];
    1921                                         }
    1922                                    }
    1923                                    solution[iNew+numberTotal] = sol;
    1924                                    // and do cost in nonLinearCost
    1925                                    if (model->algorithm() > 0)
    1926                                         model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]);
    1927                                    if (sol > upper_[iSet] + primalTolerance) {
    1928                                         setAbove(iSet);
    1929                                         lower[iNew+numberTotal] = upper_[iSet];
    1930                                         upper[iNew+numberTotal] = COIN_DBL_MAX;
    1931                                    } else if (sol < lower_[iSet] - primalTolerance) {
    1932                                         setBelow(iSet);
    1933                                         lower[iNew+numberTotal] = -COIN_DBL_MAX;
    1934                                         upper[iNew+numberTotal] = lower_[iSet];
    1935                                    } else {
    1936                                         setFeasible(iSet);
    1937                                         lower[iNew+numberTotal] = lower_[iSet];
    1938                                         upper[iNew+numberTotal] = upper_[iSet];
    1939                                    }
    1940                                    cost[iNew+numberTotal] = weight(iSet) * infeasibilityCost;
    1941                               }
    1942                               number2++;
    1943                          } else {
    1944                               value = array[iIndex];
    1945                               int iKey = keyVariable_[iSet];
    1946                               if (iKey >= numberColumns)
    1947                                    otherValue = - otherValue; //allow for - sign on slack
    1948                          }
    1949                          value -= otherValue;
    1950                          array[iIndex] = value;
    1951                     }
    1952                }
    1953           }
    1954           if (iSetIn >= 0 && toIndex_[iSetIn] < 0) {
    1955                // Do incoming
    1956                update->setPacked(); // just in case no elements
    1957                toIndex_[iSetIn] = number2;
    1958                int iNew = number2 - number;
    1959                fromIndex_[number2-number] = iSetIn;
    1960                // Save current cost of key
    1961                double currentCost;
    1962                int key = keyVariable_[iSetIn];
    1963                if (key < numberColumns)
    1964                     currentCost = cost[key];
    1965                else
    1966                     currentCost = -weight(iSetIn) * infeasibilityCost;
    1967                changeCost_[number2-number] = currentCost;
    1968                index[number2] = numberRows + iNew;
    1969                // do key stuff
    1970                int iKey = keyVariable_[iSetIn];
    1971                if (iKey < numberColumns) {
    1972                     if (gubSlackIn_ < 0)
    1973                          array[number2] = 1.0;
    1974                     else
    1975                          array[number2] = -1.0;
    1976                     pivotVariable[numberRows+iNew] = iKey;
    1977                     // Do I need to recompute?
    1978                     double sol;
    1979                     assert (getStatus(iSetIn) != ClpSimplex::basic);
    1980                     if (getStatus(iSetIn) == ClpSimplex::atLowerBound)
    1981                          sol = lower_[iSetIn];
    1982                     else
    1983                          sol = upper_[iSetIn];
    1984                     if ((gubType_ & 8) != 0) {
    1985                          int iColumn = next_[iKey];
    1986                          // sum all non-key variables
    1987                          while(iColumn >= 0) {
    1988                               sol -= solution[iColumn];
    1989                               iColumn = next_[iColumn];
    1990                          }
    1991                     } else {
    1992                          // bounds exist - sum over all except key
    1993                          int stop = -(iKey + 1);
    1994                          int iColumn = next_[iKey];
    1995                          // sum all non-key variables
    1996                          while(iColumn != stop) {
    1997                               if (iColumn < 0)
    1998                                    iColumn = -iColumn - 1;
    1999                               sol -= solution[iColumn];
    2000                               iColumn = next_[iColumn];
    2001                          }
    2002                     }
    2003                     solution[iKey] = sol;
    2004                     if (model->algorithm() > 0)
    2005                          model->nonLinearCost()->setOne(iKey, sol);
    2006                     //assert (fabs(sol-solution[iKey])<1.0e-3);
    2007                } else {
    2008                     // gub slack is basic
    2009                     array[number2] = -1.0;
    2010                     pivotVariable[numberRows+iNew] = iNew + numberTotal;
    2011                     model->djRegion()[iNew+numberTotal] = 0.0;
    2012                     double sol = 0.0;
    2013                     if ((gubType_ & 8) != 0) {
    2014                          int iColumn = next_[iKey];
    2015                          // sum all non-key variables
    2016                          while(iColumn >= 0) {
    2017                               sol += solution[iColumn];
    2018                               iColumn = next_[iColumn];
    2019                          }
    2020                     } else {
    2021                          // bounds exist - sum over all except key
    2022                          int stop = -(iKey + 1);
    2023                          int iColumn = next_[iKey];
    2024                          // sum all non-key variables
    2025                          while(iColumn != stop) {
    2026                               if (iColumn < 0)
    2027                                    iColumn = -iColumn - 1;
    2028                               sol += solution[iColumn];
    2029                               iColumn = next_[iColumn];
    2030                          }
    2031                     }
    2032                     solution[iNew+numberTotal] = sol;
    2033                     // and do cost in nonLinearCost
    2034                     if (model->algorithm() > 0)
    2035                          model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]);
    2036                     if (sol > upper_[iSetIn] + primalTolerance) {
    2037                          setAbove(iSetIn);
    2038                          lower[iNew+numberTotal] = upper_[iSetIn];
    2039                          upper[iNew+numberTotal] = COIN_DBL_MAX;
    2040                     } else if (sol < lower_[iSetIn] - primalTolerance) {
    2041                          setBelow(iSetIn);
    2042                          lower[iNew+numberTotal] = -COIN_DBL_MAX;
    2043                          upper[iNew+numberTotal] = lower_[iSetIn];
    2044                     } else {
    2045                          setFeasible(iSetIn);
    2046                          lower[iNew+numberTotal] = lower_[iSetIn];
    2047                          upper[iNew+numberTotal] = upper_[iSetIn];
    2048                     }
    2049                     cost[iNew+numberTotal] = weight(iSetIn) * infeasibilityCost;
    2050                }
    2051                number2++;
    2052           }
    2053           // mark end
    2054           fromIndex_[number2-number] = -1;
    2055           returnCode = number2 - number;
    2056           // make sure lower_ upper_ adjusted
    2057           synchronize(model, 9);
    2058      } else {
    2059           // take off?
    2060           if (number > saveNumber_) {
    2061                // clear
    2062                double theta = model->theta();
    2063                double * solution = model->solutionRegion();
    2064                for (i = saveNumber_; i < number; i++) {
    2065                     int iRow = index[i];
    2066                     int iColumn = pivotVariable[iRow];
     1787  // I think we only need to bother about sets with two in basis or incoming set
     1788  int number = update->getNumElements();
     1789  double *array = update->denseVector();
     1790  int *index = update->getIndices();
     1791  int i;
     1792  assert(!number || update->packedMode());
     1793  int *pivotVariable = model->pivotVariable();
     1794  int numberRows = model->numberRows();
     1795  int numberColumns = model->numberColumns();
     1796  int numberTotal = numberRows + numberColumns;
     1797  int sequenceIn = model->sequenceIn();
     1798  int returnCode = 0;
     1799  int iSetIn;
     1800  if (sequenceIn < numberColumns) {
     1801    iSetIn = backward_[sequenceIn];
     1802    gubSlackIn_ = -1; // in case set
     1803  } else if (sequenceIn < numberRows + numberColumns) {
     1804    iSetIn = -1;
     1805    gubSlackIn_ = -1; // in case set
     1806  } else {
     1807    iSetIn = gubSlackIn_;
     1808  }
     1809  double *lower = model->lowerRegion();
     1810  double *upper = model->upperRegion();
     1811  double *cost = model->costRegion();
     1812  double *solution = model->solutionRegion();
     1813  int number2 = number;
     1814  if (!mode) {
     1815    double primalTolerance = model->primalTolerance();
     1816    double infeasibilityCost = model->infeasibilityCost();
     1817    // extend
     1818    saveNumber_ = number;
     1819    for (i = 0; i < number; i++) {
     1820      int iRow = index[i];
     1821      int iPivot = pivotVariable[iRow];
     1822      if (iPivot < numberColumns) {
     1823        int iSet = backward_[iPivot];
     1824        if (iSet >= 0) {
     1825          // two (or more) in set
     1826          int iIndex = toIndex_[iSet];
     1827          double otherValue = array[i];
     1828          double value;
     1829          if (iIndex < 0) {
     1830            toIndex_[iSet] = number2;
     1831            int iNew = number2 - number;
     1832            fromIndex_[number2 - number] = iSet;
     1833            iIndex = number2;
     1834            index[number2] = numberRows + iNew;
     1835            // do key stuff
     1836            int iKey = keyVariable_[iSet];
     1837            if (iKey < numberColumns) {
     1838              // Save current cost of key
     1839              changeCost_[number2 - number] = cost[iKey];
     1840              if (iSet != iSetIn)
     1841                value = 0.0;
     1842              else if (iSetIn != gubSlackIn_)
     1843                value = 1.0;
     1844              else
     1845                value = -1.0;
     1846              pivotVariable[numberRows + iNew] = iKey;
     1847              // Do I need to recompute?
     1848              double sol;
     1849              assert(getStatus(iSet) != ClpSimplex::basic);
     1850              if (getStatus(iSet) == ClpSimplex::atLowerBound)
     1851                sol = lower_[iSet];
     1852              else
     1853                sol = upper_[iSet];
     1854              if ((gubType_ & 8) != 0) {
     1855                int iColumn = next_[iKey];
     1856                // sum all non-key variables
     1857                while (iColumn >= 0) {
     1858                  sol -= solution[iColumn];
     1859                  iColumn = next_[iColumn];
     1860                }
     1861              } else {
     1862                int stop = -(iKey + 1);
     1863                int iColumn = next_[iKey];
     1864                // sum all non-key variables
     1865                while (iColumn != stop) {
     1866                  if (iColumn < 0)
     1867                    iColumn = -iColumn - 1;
     1868                  sol -= solution[iColumn];
     1869                  iColumn = next_[iColumn];
     1870                }
     1871              }
     1872              solution[iKey] = sol;
     1873              if (model->algorithm() > 0)
     1874                model->nonLinearCost()->setOne(iKey, sol);
     1875              //assert (fabs(sol-solution[iKey])<1.0e-3);
     1876            } else {
     1877              // gub slack is basic
     1878              // Save current cost of key
     1879              changeCost_[number2 - number] = -weight(iSet) * infeasibilityCost;
     1880              otherValue = -otherValue; //allow for - sign on slack
     1881              if (iSet != iSetIn)
     1882                value = 0.0;
     1883              else
     1884                value = -1.0;
     1885              pivotVariable[numberRows + iNew] = iNew + numberTotal;
     1886              model->djRegion()[iNew + numberTotal] = 0.0;
     1887              double sol = 0.0;
     1888              if ((gubType_ & 8) != 0) {
     1889                int iColumn = next_[iKey];
     1890                // sum all non-key variables
     1891                while (iColumn >= 0) {
     1892                  sol += solution[iColumn];
     1893                  iColumn = next_[iColumn];
     1894                }
     1895              } else {
     1896                int stop = -(iKey + 1);
     1897                int iColumn = next_[iKey];
     1898                // sum all non-key variables
     1899                while (iColumn != stop) {
     1900                  if (iColumn < 0)
     1901                    iColumn = -iColumn - 1;
     1902                  sol += solution[iColumn];
     1903                  iColumn = next_[iColumn];
     1904                }
     1905              }
     1906              solution[iNew + numberTotal] = sol;
     1907              // and do cost in nonLinearCost
     1908              if (model->algorithm() > 0)
     1909                model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]);
     1910              if (sol > upper_[iSet] + primalTolerance) {
     1911                setAbove(iSet);
     1912                lower[iNew + numberTotal] = upper_[iSet];
     1913                upper[iNew + numberTotal] = COIN_DBL_MAX;
     1914              } else if (sol < lower_[iSet] - primalTolerance) {
     1915                setBelow(iSet);
     1916                lower[iNew + numberTotal] = -COIN_DBL_MAX;
     1917                upper[iNew + numberTotal] = lower_[iSet];
     1918              } else {
     1919                setFeasible(iSet);
     1920                lower[iNew + numberTotal] = lower_[iSet];
     1921                upper[iNew + numberTotal] = upper_[iSet];
     1922              }
     1923              cost[iNew + numberTotal] = weight(iSet) * infeasibilityCost;
     1924            }
     1925            number2++;
     1926          } else {
     1927            value = array[iIndex];
     1928            int iKey = keyVariable_[iSet];
     1929            if (iKey >= numberColumns)
     1930              otherValue = -otherValue; //allow for - sign on slack
     1931          }
     1932          value -= otherValue;
     1933          array[iIndex] = value;
     1934        }
     1935      }
     1936    }
     1937    if (iSetIn >= 0 && toIndex_[iSetIn] < 0) {
     1938      // Do incoming
     1939      update->setPacked(); // just in case no elements
     1940      toIndex_[iSetIn] = number2;
     1941      int iNew = number2 - number;
     1942      fromIndex_[number2 - number] = iSetIn;
     1943      // Save current cost of key
     1944      double currentCost;
     1945      int key = keyVariable_[iSetIn];
     1946      if (key < numberColumns)
     1947        currentCost = cost[key];
     1948      else
     1949        currentCost = -weight(iSetIn) * infeasibilityCost;
     1950      changeCost_[number2 - number] = currentCost;
     1951      index[number2] = numberRows + iNew;
     1952      // do key stuff
     1953      int iKey = keyVariable_[iSetIn];
     1954      if (iKey < numberColumns) {
     1955        if (gubSlackIn_ < 0)
     1956          array[number2] = 1.0;
     1957        else
     1958          array[number2] = -1.0;
     1959        pivotVariable[numberRows + iNew] = iKey;
     1960        // Do I need to recompute?
     1961        double sol;
     1962        assert(getStatus(iSetIn) != ClpSimplex::basic);
     1963        if (getStatus(iSetIn) == ClpSimplex::atLowerBound)
     1964          sol = lower_[iSetIn];
     1965        else
     1966          sol = upper_[iSetIn];
     1967        if ((gubType_ & 8) != 0) {
     1968          int iColumn = next_[iKey];
     1969          // sum all non-key variables
     1970          while (iColumn >= 0) {
     1971            sol -= solution[iColumn];
     1972            iColumn = next_[iColumn];
     1973          }
     1974        } else {
     1975          // bounds exist - sum over all except key
     1976          int stop = -(iKey + 1);
     1977          int iColumn = next_[iKey];
     1978          // sum all non-key variables
     1979          while (iColumn != stop) {
     1980            if (iColumn < 0)
     1981              iColumn = -iColumn - 1;
     1982            sol -= solution[iColumn];
     1983            iColumn = next_[iColumn];
     1984          }
     1985        }
     1986        solution[iKey] = sol;
     1987        if (model->algorithm() > 0)
     1988          model->nonLinearCost()->setOne(iKey, sol);
     1989        //assert (fabs(sol-solution[iKey])<1.0e-3);
     1990      } else {
     1991        // gub slack is basic
     1992        array[number2] = -1.0;
     1993        pivotVariable[numberRows + iNew] = iNew + numberTotal;
     1994        model->djRegion()[iNew + numberTotal] = 0.0;
     1995        double sol = 0.0;
     1996        if ((gubType_ & 8) != 0) {
     1997          int iColumn = next_[iKey];
     1998          // sum all non-key variables
     1999          while (iColumn >= 0) {
     2000            sol += solution[iColumn];
     2001            iColumn = next_[iColumn];
     2002          }
     2003        } else {
     2004          // bounds exist - sum over all except key
     2005          int stop = -(iKey + 1);
     2006          int iColumn = next_[iKey];
     2007          // sum all non-key variables
     2008          while (iColumn != stop) {
     2009            if (iColumn < 0)
     2010              iColumn = -iColumn - 1;
     2011            sol += solution[iColumn];
     2012            iColumn = next_[iColumn];
     2013          }
     2014        }
     2015        solution[iNew + numberTotal] = sol;
     2016        // and do cost in nonLinearCost
     2017        if (model->algorithm() > 0)
     2018          model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]);
     2019        if (sol > upper_[iSetIn] + primalTolerance) {
     2020          setAbove(iSetIn);
     2021          lower[iNew + numberTotal] = upper_[iSetIn];
     2022          upper[iNew + numberTotal] = COIN_DBL_MAX;
     2023        } else if (sol < lower_[iSetIn] - primalTolerance) {
     2024          setBelow(iSetIn);
     2025          lower[iNew + numberTotal] = -COIN_DBL_MAX;
     2026          upper[iNew + numberTotal] = lower_[iSetIn];
     2027        } else {
     2028          setFeasible(iSetIn);
     2029          lower[iNew + numberTotal] = lower_[iSetIn];
     2030          upper[iNew + numberTotal] = upper_[iSetIn];
     2031        }
     2032        cost[iNew + numberTotal] = weight(iSetIn) * infeasibilityCost;
     2033      }
     2034      number2++;
     2035    }
     2036    // mark end
     2037    fromIndex_[number2 - number] = -1;
     2038    returnCode = number2 - number;
     2039    // make sure lower_ upper_ adjusted
     2040    synchronize(model, 9);
     2041  } else {
     2042    // take off?
     2043    if (number > saveNumber_) {
     2044      // clear
     2045      double theta = model->theta();
     2046      double *solution = model->solutionRegion();
     2047      for (i = saveNumber_; i < number; i++) {
     2048        int iRow = index[i];
     2049        int iColumn = pivotVariable[iRow];
    20672050#ifdef CLP_DEBUG_PRINT
    2068                     printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
    2069                            iColumn, fromIndex_[i-saveNumber_], lower[iColumn], upper[iColumn], array[i],
    2070                            solution[iColumn], solution[iColumn] - model->theta()*array[i], model->theta());
     2051        printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
     2052          iColumn, fromIndex_[i - saveNumber_], lower[iColumn], upper[iColumn], array[i],
     2053          solution[iColumn], solution[iColumn] - model->theta() * array[i], model->theta());
    20712054#endif
    2072                     double value = array[i];
    2073                     array[i] = 0.0;
    2074                     int iSet = fromIndex_[i-saveNumber_];
    2075                     toIndex_[iSet] = -1;
    2076                     if (iSet == iSetIn && iColumn < numberColumns) {
    2077                          // update as may need value
    2078                          solution[iColumn] -= theta * value;
    2079                     }
    2080                }
    2081           }
     2055        double value = array[i];
     2056        array[i] = 0.0;
     2057        int iSet = fromIndex_[i - saveNumber_];
     2058        toIndex_[iSet] = -1;
     2059        if (iSet == iSetIn && iColumn < numberColumns) {
     2060          // update as may need value
     2061          solution[iColumn] -= theta * value;
     2062        }
     2063      }
     2064    }
    20822065#ifdef CLP_DEBUG
    2083           for (i = 0; i < numberSets_; i++)
    2084                assert(toIndex_[i] == -1);
     2066    for (i = 0; i < numberSets_; i++)
     2067      assert(toIndex_[i] == -1);
    20852068#endif
    2086           number2 = saveNumber_;
    2087      }
    2088      update->setNumElements(number2);
    2089      return returnCode;
     2069    number2 = saveNumber_;
     2070  }
     2071  update->setNumElements(number2);
     2072  return returnCode;
    20902073}
    20912074/*
     
    20942077     Remember to update here when settled down
    20952078*/
    2096 void
    2097 ClpGubMatrix::primalExpanded(ClpSimplex * model, int mode)
     2079void ClpGubMatrix::primalExpanded(ClpSimplex *model, int mode)
    20982080{
    2099      int numberColumns = model->numberColumns();
    2100      switch (mode) {
    2101           // If key variable then slot in gub rhs so will get correct contribution
    2102      case 0: {
    2103           int i;
    2104           double * solution = model->solutionRegion();
    2105           ClpSimplex::Status iStatus;
    2106           for (i = 0; i < numberSets_; i++) {
    2107                int iColumn = keyVariable_[i];
    2108                if (iColumn < numberColumns) {
    2109                     // key is structural - where is slack
    2110                     iStatus = getStatus(i);
    2111                     assert (iStatus != ClpSimplex::basic);
    2112                     if (iStatus == ClpSimplex::atLowerBound)
    2113                          solution[iColumn] = lower_[i];
    2114                     else
    2115                          solution[iColumn] = upper_[i];
    2116                }
    2117           }
    2118      }
    2119      break;
    2120      // Compute values of key variables
    2121      case 1: {
    2122           int i;
    2123           double * solution = model->solutionRegion();
    2124           //const int * columnLength = matrix_->getVectorLengths();
    2125           //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    2126           //const int * row = matrix_->getIndices();
    2127           //const double * elementByColumn = matrix_->getElements();
    2128           //int * pivotVariable = model->pivotVariable();
    2129           sumPrimalInfeasibilities_ = 0.0;
    2130           numberPrimalInfeasibilities_ = 0;
    2131           double primalTolerance = model->primalTolerance();
    2132           double relaxedTolerance = primalTolerance;
    2133           // we can't really trust infeasibilities if there is primal error
    2134           double error = CoinMin(1.0e-2, model->largestPrimalError());
    2135           // allow tolerance at least slightly bigger than standard
    2136           relaxedTolerance = relaxedTolerance +  error;
    2137           // but we will be using difference
    2138           relaxedTolerance -= primalTolerance;
    2139           sumOfRelaxedPrimalInfeasibilities_ = 0.0;
    2140           for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds)
    2141                int kColumn = keyVariable_[i];
    2142                double value = 0.0;
    2143                if ((gubType_ & 8) != 0) {
    2144                     int iColumn = next_[kColumn];
    2145                     // sum all non-key variables
    2146                     while(iColumn >= 0) {
    2147                          value += solution[iColumn];
    2148                          iColumn = next_[iColumn];
    2149                     }
    2150                } else {
    2151                     // bounds exist - sum over all except key
    2152                     int stop = -(kColumn + 1);
    2153                     int iColumn = next_[kColumn];
    2154                     // sum all non-key variables
    2155                     while(iColumn != stop) {
    2156                          if (iColumn < 0)
    2157                               iColumn = -iColumn - 1;
    2158                          value += solution[iColumn];
    2159                          iColumn = next_[iColumn];
    2160                     }
    2161                }
    2162                if (kColumn < numberColumns) {
    2163                     // make sure key is basic - so will be skipped in values pass
    2164                     model->setStatus(kColumn, ClpSimplex::basic);
    2165                     // feasibility will be done later
    2166                     assert (getStatus(i) != ClpSimplex::basic);
    2167                     if (getStatus(i) == ClpSimplex::atUpperBound)
    2168                          solution[kColumn] = upper_[i] - value;
    2169                     else
    2170                          solution[kColumn] = lower_[i] - value;
    2171                     //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
    2172                } else {
    2173                     // slack is key
    2174                     assert (getStatus(i) == ClpSimplex::basic);
    2175                     double infeasibility = 0.0;
    2176                     if (value > upper_[i] + primalTolerance) {
    2177                          infeasibility = value - upper_[i] - primalTolerance;
    2178                          setAbove(i);
    2179                     } else if (value < lower_[i] - primalTolerance) {
    2180                          infeasibility = lower_[i] - value - primalTolerance ;
    2181                          setBelow(i);
    2182                     } else {
    2183                          setFeasible(i);
    2184                     }
    2185                     //printf("Value of key slack for set %d is %g\n",i,value);
    2186                     if (infeasibility > 0.0) {
    2187                          sumPrimalInfeasibilities_ += infeasibility;
    2188                          if (infeasibility > relaxedTolerance)
    2189                               sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
    2190                          numberPrimalInfeasibilities_ ++;
    2191                     }
    2192                }
    2193           }
    2194      }
    2195      break;
    2196      // Report on infeasibilities of key variables
    2197      case 2: {
    2198           model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() +
    2199                                              sumPrimalInfeasibilities_);
    2200           model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() +
    2201                                                 numberPrimalInfeasibilities_);
    2202           model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() +
    2203                     sumOfRelaxedPrimalInfeasibilities_);
    2204      }
    2205      break;
    2206      }
     2081  int numberColumns = model->numberColumns();
     2082  switch (mode) {
     2083    // If key variable then slot in gub rhs so will get correct contribution
     2084  case 0: {
     2085    int i;
     2086    double *solution = model->solutionRegion();
     2087    ClpSimplex::Status iStatus;
     2088    for (i = 0; i < numberSets_; i++) {
     2089      int iColumn = keyVariable_[i];
     2090      if (iColumn < numberColumns) {
     2091        // key is structural - where is slack
     2092        iStatus = getStatus(i);
     2093        assert(iStatus != ClpSimplex::basic);
     2094        if (iStatus == ClpSimplex::atLowerBound)
     2095          solution[iColumn] = lower_[i];
     2096        else
     2097          solution[iColumn] = upper_[i];
     2098      }
     2099    }
     2100  } break;
     2101  // Compute values of key variables
     2102  case 1: {
     2103    int i;
     2104    double *solution = model->solutionRegion();
     2105    //const int * columnLength = matrix_->getVectorLengths();
     2106    //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     2107    //const int * row = matrix_->getIndices();
     2108    //const double * elementByColumn = matrix_->getElements();
     2109    //int * pivotVariable = model->pivotVariable();
     2110    sumPrimalInfeasibilities_ = 0.0;
     2111    numberPrimalInfeasibilities_ = 0;
     2112    double primalTolerance = model->primalTolerance();
     2113    double relaxedTolerance = primalTolerance;
     2114    // we can't really trust infeasibilities if there is primal error
     2115    double error = CoinMin(1.0e-2, model->largestPrimalError());
     2116    // allow tolerance at least slightly bigger than standard
     2117    relaxedTolerance = relaxedTolerance + error;
     2118    // but we will be using difference
     2119    relaxedTolerance -= primalTolerance;
     2120    sumOfRelaxedPrimalInfeasibilities_ = 0.0;
     2121    for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds)
     2122      int kColumn = keyVariable_[i];
     2123      double value = 0.0;
     2124      if ((gubType_ & 8) != 0) {
     2125        int iColumn = next_[kColumn];
     2126        // sum all non-key variables
     2127        while (iColumn >= 0) {
     2128          value += solution[iColumn];
     2129          iColumn = next_[iColumn];
     2130        }
     2131      } else {
     2132        // bounds exist - sum over all except key
     2133        int stop = -(kColumn + 1);
     2134        int iColumn = next_[kColumn];
     2135        // sum all non-key variables
     2136        while (iColumn != stop) {
     2137          if (iColumn < 0)
     2138            iColumn = -iColumn - 1;
     2139          value += solution[iColumn];
     2140          iColumn = next_[iColumn];
     2141        }
     2142      }
     2143      if (kColumn < numberColumns) {
     2144        // make sure key is basic - so will be skipped in values pass
     2145        model->setStatus(kColumn, ClpSimplex::basic);
     2146        // feasibility will be done later
     2147        assert(getStatus(i) != ClpSimplex::basic);
     2148        if (getStatus(i) == ClpSimplex::atUpperBound)
     2149          solution[kColumn] = upper_[i] - value;
     2150        else
     2151          solution[kColumn] = lower_[i] - value;
     2152        //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
     2153      } else {
     2154        // slack is key
     2155        assert(getStatus(i) == ClpSimplex::basic);
     2156        double infeasibility = 0.0;
     2157        if (value > upper_[i] + primalTolerance) {
     2158          infeasibility = value - upper_[i] - primalTolerance;
     2159          setAbove(i);
     2160        } else if (value < lower_[i] - primalTolerance) {
     2161          infeasibility = lower_[i] - value - primalTolerance;
     2162          setBelow(i);
     2163        } else {
     2164          setFeasible(i);
     2165        }
     2166        //printf("Value of key slack for set %d is %g\n",i,value);
     2167        if (infeasibility > 0.0) {
     2168          sumPrimalInfeasibilities_ += infeasibility;
     2169          if (infeasibility > relaxedTolerance)
     2170            sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
     2171          numberPrimalInfeasibilities_++;
     2172        }
     2173      }
     2174    }
     2175  } break;
     2176  // Report on infeasibilities of key variables
     2177  case 2: {
     2178    model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() + sumPrimalInfeasibilities_);
     2179    model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() + numberPrimalInfeasibilities_);
     2180    model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() + sumOfRelaxedPrimalInfeasibilities_);
     2181  } break;
     2182  }
    22072183}
    22082184/*
     
    22112187     Remember to update here when settled down
    22122188*/
    2213 void
    2214 ClpGubMatrix::dualExpanded(ClpSimplex * model,
    2215                            CoinIndexedVector * array,
    2216                            double * /*other*/, int mode)
     2189void ClpGubMatrix::dualExpanded(ClpSimplex *model,
     2190  CoinIndexedVector *array,
     2191  double * /*other*/, int mode)
    22172192{
    2218      switch (mode) {
    2219           // modify costs before transposeUpdate
    2220      case 0: {
    2221           int i;
    2222           double * cost = model->costRegion();
    2223           // not dual values yet
    2224           //assert (!other);
    2225           //double * work = array->denseVector();
    2226           double infeasibilityCost = model->infeasibilityCost();
    2227           int * pivotVariable = model->pivotVariable();
    2228           int numberRows = model->numberRows();
    2229           int numberColumns = model->numberColumns();
    2230           for (i = 0; i < numberRows; i++) {
    2231                int iPivot = pivotVariable[i];
    2232                if (iPivot < numberColumns) {
    2233                     int iSet = backward_[iPivot];
    2234                     if (iSet >= 0) {
    2235                          int kColumn = keyVariable_[iSet];
    2236                          double costValue;
    2237                          if (kColumn < numberColumns) {
    2238                               // structural has cost
    2239                               costValue = cost[kColumn];
    2240                          } else {
    2241                               // slack is key
    2242                               assert (getStatus(iSet) == ClpSimplex::basic);
    2243                               // negative as -1.0 for slack
    2244                               costValue = -weight(iSet) * infeasibilityCost;
    2245                          }
    2246                          array->add(i, -costValue); // was work[i]-costValue
    2247                     }
    2248                }
    2249           }
    2250      }
    2251      break;
    2252      // create duals for key variables (without check on dual infeasible)
    2253      case 1: {
    2254           // If key slack then dual 0.0 (if feasible)
    2255           // dj for key is zero so that defines dual on set
    2256           int i;
    2257           double * dj = model->djRegion();
    2258           int numberColumns = model->numberColumns();
    2259           double infeasibilityCost = model->infeasibilityCost();
    2260           for (i = 0; i < numberSets_; i++) {
    2261                int kColumn = keyVariable_[i];
    2262                if (kColumn < numberColumns) {
    2263                     // dj without set
    2264                     double value = dj[kColumn];
    2265                     // Now subtract out from all
    2266                     dj[kColumn] = 0.0;
    2267                     int iColumn = next_[kColumn];
    2268                     // modify all non-key variables
    2269                     while(iColumn >= 0) {
    2270                          dj[iColumn] -= value;
    2271                          iColumn = next_[iColumn];
    2272                     }
    2273                } else {
    2274                     // slack key - may not be feasible
    2275                     assert (getStatus(i) == ClpSimplex::basic);
    2276                     // negative as -1.0 for slack
    2277                     double value = -weight(i) * infeasibilityCost;
    2278                     if (value) {
    2279                          int iColumn = next_[kColumn];
    2280                          // modify all non-key variables basic
    2281                          while(iColumn >= 0) {
    2282                               dj[iColumn] -= value;
    2283                               iColumn = next_[iColumn];
    2284                          }
    2285                     }
    2286                }
    2287           }
    2288      }
    2289      break;
    2290      // as 1 but check slacks and compute djs
    2291      case 2: {
    2292           // If key slack then dual 0.0
    2293           // If not then slack could be dual infeasible
    2294           // dj for key is zero so that defines dual on set
    2295           int i;
    2296           // make sure fromIndex will not confuse pricing
    2297           fromIndex_[0] = -1;
     2193  switch (mode) {
     2194    // modify costs before transposeUpdate
     2195  case 0: {
     2196    int i;
     2197    double *cost = model->costRegion();
     2198    // not dual values yet
     2199    //assert (!other);
     2200    //double * work = array->denseVector();
     2201    double infeasibilityCost = model->infeasibilityCost();
     2202    int *pivotVariable = model->pivotVariable();
     2203    int numberRows = model->numberRows();
     2204    int numberColumns = model->numberColumns();
     2205    for (i = 0; i < numberRows; i++) {
     2206      int iPivot = pivotVariable[i];
     2207      if (iPivot < numberColumns) {
     2208        int iSet = backward_[iPivot];
     2209        if (iSet >= 0) {
     2210          int kColumn = keyVariable_[iSet];
     2211          double costValue;
     2212          if (kColumn < numberColumns) {
     2213            // structural has cost
     2214            costValue = cost[kColumn];
     2215          } else {
     2216            // slack is key
     2217            assert(getStatus(iSet) == ClpSimplex::basic);
     2218            // negative as -1.0 for slack
     2219            costValue = -weight(iSet) * infeasibilityCost;
     2220          }
     2221          array->add(i, -costValue); // was work[i]-costValue
     2222        }
     2223      }
     2224    }
     2225  } break;
     2226  // create duals for key variables (without check on dual infeasible)
     2227  case 1: {
     2228    // If key slack then dual 0.0 (if feasible)
     2229    // dj for key is zero so that defines dual on set
     2230    int i;
     2231    double *dj = model->djRegion();
     2232    int numberColumns = model->numberColumns();
     2233    double infeasibilityCost = model->infeasibilityCost();
     2234    for (i = 0; i < numberSets_; i++) {
     2235      int kColumn = keyVariable_[i];
     2236      if (kColumn < numberColumns) {
     2237        // dj without set
     2238        double value = dj[kColumn];
     2239        // Now subtract out from all
     2240        dj[kColumn] = 0.0;
     2241        int iColumn = next_[kColumn];
     2242        // modify all non-key variables
     2243        while (iColumn >= 0) {
     2244          dj[iColumn] -= value;
     2245          iColumn = next_[iColumn];
     2246        }
     2247      } else {
     2248        // slack key - may not be feasible
     2249        assert(getStatus(i) == ClpSimplex::basic);
     2250        // negative as -1.0 for slack
     2251        double value = -weight(i) * infeasibilityCost;
     2252        if (value) {
     2253          int iColumn = next_[kColumn];
     2254          // modify all non-key variables basic
     2255          while (iColumn >= 0) {
     2256            dj[iColumn] -= value;
     2257            iColumn = next_[iColumn];
     2258          }
     2259        }
     2260      }
     2261    }
     2262  } break;
     2263  // as 1 but check slacks and compute djs
     2264  case 2: {
     2265    // If key slack then dual 0.0
     2266    // If not then slack could be dual infeasible
     2267    // dj for key is zero so that defines dual on set
     2268    int i;
     2269    // make sure fromIndex will not confuse pricing
     2270    fromIndex_[0] = -1;
     2271    possiblePivotKey_ = -1;
     2272    // Create array
     2273    int numberColumns = model->numberColumns();
     2274    int *pivotVariable = model->pivotVariable();
     2275    int numberRows = model->numberRows();
     2276    for (i = 0; i < numberRows; i++) {
     2277      int iPivot = pivotVariable[i];
     2278      if (iPivot < numberColumns)
     2279        backToPivotRow_[iPivot] = i;
     2280    }
     2281    if (noCheck_ >= 0) {
     2282      if (infeasibilityWeight_ != model->infeasibilityCost()) {
     2283        // don't bother checking
     2284        sumDualInfeasibilities_ = 100.0;
     2285        numberDualInfeasibilities_ = 1;
     2286        sumOfRelaxedDualInfeasibilities_ = 100.0;
     2287        return;
     2288      }
     2289    }
     2290    double *dj = model->djRegion();
     2291    double *dual = model->dualRowSolution();
     2292    double *cost = model->costRegion();
     2293    ClpSimplex::Status iStatus;
     2294    const int *columnLength = matrix_->getVectorLengths();
     2295    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     2296    const int *row = matrix_->getIndices();
     2297    const double *elementByColumn = matrix_->getElements();
     2298    double infeasibilityCost = model->infeasibilityCost();
     2299    sumDualInfeasibilities_ = 0.0;
     2300    numberDualInfeasibilities_ = 0;
     2301    double dualTolerance = model->dualTolerance();
     2302    double relaxedTolerance = dualTolerance;
     2303    // we can't really trust infeasibilities if there is dual error
     2304    double error = CoinMin(1.0e-2, model->largestDualError());
     2305    // allow tolerance at least slightly bigger than standard
     2306    relaxedTolerance = relaxedTolerance + error;
     2307    // but we will be using difference
     2308    relaxedTolerance -= dualTolerance;
     2309    sumOfRelaxedDualInfeasibilities_ = 0.0;
     2310    for (i = 0; i < numberSets_; i++) {
     2311      int kColumn = keyVariable_[i];
     2312      if (kColumn < numberColumns) {
     2313        // dj without set
     2314        double value = cost[kColumn];
     2315        for (CoinBigIndex j = columnStart[kColumn];
     2316             j < columnStart[kColumn] + columnLength[kColumn]; j++) {
     2317          int iRow = row[j];
     2318          value -= dual[iRow] * elementByColumn[j];
     2319        }
     2320        // Now subtract out from all
     2321        dj[kColumn] -= value;
     2322        int stop = -(kColumn + 1);
     2323        kColumn = next_[kColumn];
     2324        while (kColumn != stop) {
     2325          if (kColumn < 0)
     2326            kColumn = -kColumn - 1;
     2327          double djValue = dj[kColumn] - value;
     2328          dj[kColumn] = djValue;
     2329          ;
     2330          double infeasibility = 0.0;
     2331          iStatus = model->getStatus(kColumn);
     2332          if (iStatus == ClpSimplex::atLowerBound) {
     2333            if (djValue < -dualTolerance)
     2334              infeasibility = -djValue - dualTolerance;
     2335          } else if (iStatus == ClpSimplex::atUpperBound) {
     2336            // at upper bound
     2337            if (djValue > dualTolerance)
     2338              infeasibility = djValue - dualTolerance;
     2339          }
     2340          if (infeasibility > 0.0) {
     2341            sumDualInfeasibilities_ += infeasibility;
     2342            if (infeasibility > relaxedTolerance)
     2343              sumOfRelaxedDualInfeasibilities_ += infeasibility;
     2344            numberDualInfeasibilities_++;
     2345          }
     2346          kColumn = next_[kColumn];
     2347        }
     2348        // check slack
     2349        iStatus = getStatus(i);
     2350        assert(iStatus != ClpSimplex::basic);
     2351        double infeasibility = 0.0;
     2352        // dj of slack is -(-1.0)value
     2353        if (iStatus == ClpSimplex::atLowerBound) {
     2354          if (value < -dualTolerance)
     2355            infeasibility = -value - dualTolerance;
     2356        } else if (iStatus == ClpSimplex::atUpperBound) {
     2357          // at upper bound
     2358          if (value > dualTolerance)
     2359            infeasibility = value - dualTolerance;
     2360        }
     2361        if (infeasibility > 0.0) {
     2362          sumDualInfeasibilities_ += infeasibility;
     2363          if (infeasibility > relaxedTolerance)
     2364            sumOfRelaxedDualInfeasibilities_ += infeasibility;
     2365          numberDualInfeasibilities_++;
     2366        }
     2367      } else {
     2368        // slack key - may not be feasible
     2369        assert(getStatus(i) == ClpSimplex::basic);
     2370        // negative as -1.0 for slack
     2371        double value = -weight(i) * infeasibilityCost;
     2372        if (value) {
     2373          // Now subtract out from all
     2374          int kColumn = i + numberColumns;
     2375          int stop = -(kColumn + 1);
     2376          kColumn = next_[kColumn];
     2377          while (kColumn != stop) {
     2378            if (kColumn < 0)
     2379              kColumn = -kColumn - 1;
     2380            double djValue = dj[kColumn] - value;
     2381            dj[kColumn] = djValue;
     2382            ;
     2383            double infeasibility = 0.0;
     2384            iStatus = model->getStatus(kColumn);
     2385            if (iStatus == ClpSimplex::atLowerBound) {
     2386              if (djValue < -dualTolerance)
     2387                infeasibility = -djValue - dualTolerance;
     2388            } else if (iStatus == ClpSimplex::atUpperBound) {
     2389              // at upper bound
     2390              if (djValue > dualTolerance)
     2391                infeasibility = djValue - dualTolerance;
     2392            }
     2393            if (infeasibility > 0.0) {
     2394              sumDualInfeasibilities_ += infeasibility;
     2395              if (infeasibility > relaxedTolerance)
     2396                sumOfRelaxedDualInfeasibilities_ += infeasibility;
     2397              numberDualInfeasibilities_++;
     2398            }
     2399            kColumn = next_[kColumn];
     2400          }
     2401        }
     2402      }
     2403    }
     2404    // and get statistics for column generation
     2405    synchronize(model, 4);
     2406    infeasibilityWeight_ = -1.0;
     2407  } break;
     2408  // Report on infeasibilities of key variables
     2409  case 3: {
     2410    model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + sumDualInfeasibilities_);
     2411    model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + numberDualInfeasibilities_);
     2412    model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_);
     2413  } break;
     2414  // modify costs before transposeUpdate for partial pricing
     2415  case 4: {
     2416    // First compute new costs etc for interesting gubs
     2417    int iLook = 0;
     2418    int iSet = fromIndex_[0];
     2419    double primalTolerance = model->primalTolerance();
     2420    const double *cost = model->costRegion();
     2421    double *solution = model->solutionRegion();
     2422    double infeasibilityCost = model->infeasibilityCost();
     2423    int numberColumns = model->numberColumns();
     2424    int numberChanged = 0;
     2425    int *pivotVariable = model->pivotVariable();
     2426    while (iSet >= 0) {
     2427      int key = keyVariable_[iSet];
     2428      double value = 0.0;
     2429      // sum over all except key
     2430      if ((gubType_ & 8) != 0) {
     2431        int iColumn = next_[key];
     2432        // sum all non-key variables
     2433        while (iColumn >= 0) {
     2434          value += solution[iColumn];
     2435          iColumn = next_[iColumn];
     2436        }
     2437      } else {
     2438        // bounds exist - sum over all except key
     2439        int stop = -(key + 1);
     2440        int iColumn = next_[key];
     2441        // sum all non-key variables
     2442        while (iColumn != stop) {
     2443          if (iColumn < 0)
     2444            iColumn = -iColumn - 1;
     2445          value += solution[iColumn];
     2446          iColumn = next_[iColumn];
     2447        }
     2448      }
     2449      double costChange;
     2450      double oldCost = changeCost_[iLook];
     2451      if (key < numberColumns) {
     2452        assert(getStatus(iSet) != ClpSimplex::basic);
     2453        double sol;
     2454        if (getStatus(iSet) == ClpSimplex::atUpperBound)
     2455          sol = upper_[iSet] - value;
     2456        else
     2457          sol = lower_[iSet] - value;
     2458        solution[key] = sol;
     2459        // fix up cost
     2460        model->nonLinearCost()->setOne(key, sol);
     2461#ifdef CLP_DEBUG_PRINT
     2462        printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol,
     2463          cost[key], oldCost);
     2464#endif
     2465        costChange = cost[key] - oldCost;
     2466      } else {
     2467        // slack is key
     2468        if (value > upper_[iSet] + primalTolerance) {
     2469          setAbove(iSet);
     2470        } else if (value < lower_[iSet] - primalTolerance) {
     2471          setBelow(iSet);
     2472        } else {
     2473          setFeasible(iSet);
     2474        }
     2475        // negative as -1.0 for slack
     2476        costChange = -weight(iSet) * infeasibilityCost - oldCost;
     2477#ifdef CLP_DEBUG_PRINT
     2478        printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value,
     2479          weight(iSet) * infeasibilityCost, oldCost);
     2480#endif
     2481      }
     2482      if (costChange) {
     2483        fromIndex_[numberChanged] = iSet;
     2484        toIndex_[iSet] = numberChanged;
     2485        changeCost_[numberChanged++] = costChange;
     2486      }
     2487      iSet = fromIndex_[++iLook];
     2488    }
     2489    if (numberChanged || possiblePivotKey_ >= 0) {
     2490      // first do those in list already
     2491      int number = array->getNumElements();
     2492      array->setPacked();
     2493      int i;
     2494      double *work = array->denseVector();
     2495      int *which = array->getIndices();
     2496      for (i = 0; i < number; i++) {
     2497        int iRow = which[i];
     2498        int iPivot = pivotVariable[iRow];
     2499        if (iPivot < numberColumns) {
     2500          int iSet = backward_[iPivot];
     2501          if (iSet >= 0 && toIndex_[iSet] >= 0) {
     2502            double newValue = work[i] + changeCost_[toIndex_[iSet]];
     2503            if (!newValue)
     2504              newValue = 1.0e-100;
     2505            work[i] = newValue;
     2506            // mark as done
     2507            backward_[iPivot] = -1;
     2508          }
     2509        }
     2510        if (possiblePivotKey_ == iRow) {
     2511          double newValue = work[i] - model->dualIn();
     2512          if (!newValue)
     2513            newValue = 1.0e-100;
     2514          work[i] = newValue;
    22982515          possiblePivotKey_ = -1;
    2299           // Create array
    2300           int numberColumns = model->numberColumns();
    2301           int * pivotVariable = model->pivotVariable();
    2302           int numberRows = model->numberRows();
    2303           for (i = 0; i < numberRows; i++) {
    2304                int iPivot = pivotVariable[i];
    2305                if (iPivot < numberColumns)
    2306                     backToPivotRow_[iPivot] = i;
    2307           }
    2308           if (noCheck_ >= 0) {
    2309                if (infeasibilityWeight_ != model->infeasibilityCost()) {
    2310                     // don't bother checking
    2311                     sumDualInfeasibilities_ = 100.0;
    2312                     numberDualInfeasibilities_ = 1;
    2313                     sumOfRelaxedDualInfeasibilities_ = 100.0;
    2314                     return;
    2315                }
    2316           }
    2317           double * dj = model->djRegion();
    2318           double * dual = model->dualRowSolution();
    2319           double * cost = model->costRegion();
    2320           ClpSimplex::Status iStatus;
    2321           const int * columnLength = matrix_->getVectorLengths();
    2322           const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    2323           const int * row = matrix_->getIndices();
    2324           const double * elementByColumn = matrix_->getElements();
    2325           double infeasibilityCost = model->infeasibilityCost();
    2326           sumDualInfeasibilities_ = 0.0;
    2327           numberDualInfeasibilities_ = 0;
    2328           double dualTolerance = model->dualTolerance();
    2329           double relaxedTolerance = dualTolerance;
    2330           // we can't really trust infeasibilities if there is dual error
    2331           double error = CoinMin(1.0e-2, model->largestDualError());
    2332           // allow tolerance at least slightly bigger than standard
    2333           relaxedTolerance = relaxedTolerance +  error;
    2334           // but we will be using difference
    2335           relaxedTolerance -= dualTolerance;
    2336           sumOfRelaxedDualInfeasibilities_ = 0.0;
    2337           for (i = 0; i < numberSets_; i++) {
    2338                int kColumn = keyVariable_[i];
    2339                if (kColumn < numberColumns) {
    2340                     // dj without set
    2341                     double value = cost[kColumn];
    2342                     for (CoinBigIndex j = columnStart[kColumn];
    2343                               j < columnStart[kColumn] + columnLength[kColumn]; j++) {
    2344                          int iRow = row[j];
    2345                          value -= dual[iRow] * elementByColumn[j];
    2346                     }
    2347                     // Now subtract out from all
    2348                     dj[kColumn] -= value;
    2349                     int stop = -(kColumn + 1);
    2350                     kColumn = next_[kColumn];
    2351                     while (kColumn != stop) {
    2352                          if (kColumn < 0)
    2353                               kColumn = -kColumn - 1;
    2354                          double djValue = dj[kColumn] - value;
    2355                          dj[kColumn] = djValue;;
    2356                          double infeasibility = 0.0;
    2357                          iStatus = model->getStatus(kColumn);
    2358                          if (iStatus == ClpSimplex::atLowerBound) {
    2359                               if (djValue < -dualTolerance)
    2360                                    infeasibility = -djValue - dualTolerance;
    2361                          } else if (iStatus == ClpSimplex::atUpperBound) {
    2362                               // at upper bound
    2363                               if (djValue > dualTolerance)
    2364                                    infeasibility = djValue - dualTolerance;
    2365                          }
    2366                          if (infeasibility > 0.0) {
    2367                               sumDualInfeasibilities_ += infeasibility;
    2368                               if (infeasibility > relaxedTolerance)
    2369                                    sumOfRelaxedDualInfeasibilities_ += infeasibility;
    2370                               numberDualInfeasibilities_ ++;
    2371                          }
    2372                          kColumn = next_[kColumn];
    2373                     }
    2374                     // check slack
    2375                     iStatus = getStatus(i);
    2376                     assert (iStatus != ClpSimplex::basic);
    2377                     double infeasibility = 0.0;
    2378                     // dj of slack is -(-1.0)value
    2379                     if (iStatus == ClpSimplex::atLowerBound) {
    2380                          if (value < -dualTolerance)
    2381                               infeasibility = -value - dualTolerance;
    2382                     } else if (iStatus == ClpSimplex::atUpperBound) {
    2383                          // at upper bound
    2384                          if (value > dualTolerance)
    2385                               infeasibility = value - dualTolerance;
    2386                     }
    2387                     if (infeasibility > 0.0) {
    2388                          sumDualInfeasibilities_ += infeasibility;
    2389                          if (infeasibility > relaxedTolerance)
    2390                               sumOfRelaxedDualInfeasibilities_ += infeasibility;
    2391                          numberDualInfeasibilities_ ++;
    2392                     }
    2393                } else {
    2394                     // slack key - may not be feasible
    2395                     assert (getStatus(i) == ClpSimplex::basic);
    2396                     // negative as -1.0 for slack
    2397                     double value = -weight(i) * infeasibilityCost;
    2398                     if (value) {
    2399                          // Now subtract out from all
    2400                          int kColumn = i + numberColumns;
    2401                          int stop = -(kColumn + 1);
    2402                          kColumn = next_[kColumn];
    2403                          while (kColumn != stop) {
    2404                               if (kColumn < 0)
    2405                                    kColumn = -kColumn - 1;
    2406                               double djValue = dj[kColumn] - value;
    2407                               dj[kColumn] = djValue;;
    2408                               double infeasibility = 0.0;
    2409                               iStatus = model->getStatus(kColumn);
    2410                               if (iStatus == ClpSimplex::atLowerBound) {
    2411                                    if (djValue < -dualTolerance)
    2412                                         infeasibility = -djValue - dualTolerance;
    2413                               } else if (iStatus == ClpSimplex::atUpperBound) {
    2414                                    // at upper bound
    2415                                    if (djValue > dualTolerance)
    2416                                         infeasibility = djValue - dualTolerance;
    2417                               }
    2418                               if (infeasibility > 0.0) {
    2419                                    sumDualInfeasibilities_ += infeasibility;
    2420                                    if (infeasibility > relaxedTolerance)
    2421                                         sumOfRelaxedDualInfeasibilities_ += infeasibility;
    2422                                    numberDualInfeasibilities_ ++;
    2423                               }
    2424                               kColumn = next_[kColumn];
    2425                          }
    2426                     }
    2427                }
    2428           }
    2429           // and get statistics for column generation
    2430           synchronize(model, 4);
    2431           infeasibilityWeight_ = -1.0;
    2432      }
    2433      break;
    2434      // Report on infeasibilities of key variables
    2435      case 3: {
    2436           model->setSumDualInfeasibilities(model->sumDualInfeasibilities() +
    2437                                            sumDualInfeasibilities_);
    2438           model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() +
    2439                                               numberDualInfeasibilities_);
    2440           model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() +
    2441                     sumOfRelaxedDualInfeasibilities_);
    2442      }
    2443      break;
    2444      // modify costs before transposeUpdate for partial pricing
    2445      case 4: {
    2446           // First compute new costs etc for interesting gubs
    2447           int iLook = 0;
    2448           int iSet = fromIndex_[0];
    2449           double primalTolerance = model->primalTolerance();
    2450           const double * cost = model->costRegion();
    2451           double * solution = model->solutionRegion();
    2452           double infeasibilityCost = model->infeasibilityCost();
    2453           int numberColumns = model->numberColumns();
    2454           int numberChanged = 0;
    2455           int * pivotVariable = model->pivotVariable();
    2456           while (iSet >= 0) {
    2457                int key = keyVariable_[iSet];
    2458                double value = 0.0;
    2459                // sum over all except key
    2460                if ((gubType_ & 8) != 0) {
    2461                     int iColumn = next_[key];
    2462                     // sum all non-key variables
    2463                     while(iColumn >= 0) {
    2464                          value += solution[iColumn];
    2465                          iColumn = next_[iColumn];
    2466                     }
    2467                } else {
    2468                     // bounds exist - sum over all except key
    2469                     int stop = -(key + 1);
    2470                     int iColumn = next_[key];
    2471                     // sum all non-key variables
    2472                     while(iColumn != stop) {
    2473                          if (iColumn < 0)
    2474                               iColumn = -iColumn - 1;
    2475                          value += solution[iColumn];
    2476                          iColumn = next_[iColumn];
    2477                     }
    2478                }
    2479                double costChange;
    2480                double oldCost = changeCost_[iLook];
    2481                if (key < numberColumns) {
    2482                     assert (getStatus(iSet) != ClpSimplex::basic);
    2483                     double sol;
    2484                     if (getStatus(iSet) == ClpSimplex::atUpperBound)
    2485                          sol = upper_[iSet] - value;
    2486                     else
    2487                          sol = lower_[iSet] - value;
    2488                     solution[key] = sol;
    2489                     // fix up cost
    2490                     model->nonLinearCost()->setOne(key, sol);
    2491 #ifdef CLP_DEBUG_PRINT
    2492                     printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol,
    2493                            cost[key], oldCost);
    2494 #endif
    2495                     costChange = cost[key] - oldCost;
    2496                } else {
    2497                     // slack is key
    2498                     if (value > upper_[iSet] + primalTolerance) {
    2499                          setAbove(iSet);
    2500                     } else if (value < lower_[iSet] - primalTolerance) {
    2501                          setBelow(iSet);
    2502                     } else {
    2503                          setFeasible(iSet);
    2504                     }
    2505                     // negative as -1.0 for slack
    2506                     costChange = -weight(iSet) * infeasibilityCost - oldCost;
    2507 #ifdef CLP_DEBUG_PRINT
    2508                     printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value,
    2509                            weight(iSet)*infeasibilityCost, oldCost);
    2510 #endif
    2511                }
    2512                if (costChange) {
    2513                     fromIndex_[numberChanged] = iSet;
    2514                     toIndex_[iSet] = numberChanged;
    2515                     changeCost_[numberChanged++] = costChange;
    2516                }
    2517                iSet = fromIndex_[++iLook];
    2518           }
    2519           if (numberChanged || possiblePivotKey_ >= 0) {
    2520                // first do those in list already
    2521                int number = array->getNumElements();
    2522                array->setPacked();
    2523                int i;
    2524                double * work = array->denseVector();
    2525                int * which = array->getIndices();
    2526                for (i = 0; i < number; i++) {
    2527                     int iRow = which[i];
    2528                     int iPivot = pivotVariable[iRow];
    2529                     if (iPivot < numberColumns) {
    2530                          int iSet = backward_[iPivot];
    2531                          if (iSet >= 0 && toIndex_[iSet] >= 0) {
    2532                               double newValue = work[i] + changeCost_[toIndex_[iSet]];
    2533                               if (!newValue)
    2534                                    newValue = 1.0e-100;
    2535                               work[i] = newValue;
    2536                               // mark as done
    2537                               backward_[iPivot] = -1;
    2538                          }
    2539                     }
    2540                     if (possiblePivotKey_ == iRow) {
    2541                          double newValue = work[i] - model->dualIn();
    2542                          if (!newValue)
    2543                               newValue = 1.0e-100;
    2544                          work[i] = newValue;
    2545                          possiblePivotKey_ = -1;
    2546                     }
    2547                }
    2548                // now do rest and clean up
    2549                for (i = 0; i < numberChanged; i++) {
    2550                     int iSet = fromIndex_[i];
    2551                     int key = keyVariable_[iSet];
    2552                     int iColumn = next_[key];
    2553                     double change = changeCost_[i];
    2554                     while (iColumn >= 0) {
    2555                          if (backward_[iColumn] >= 0) {
    2556                               int iRow = backToPivotRow_[iColumn];
    2557                               assert (iRow >= 0);
    2558                               work[number] = change;
    2559                               if (possiblePivotKey_ == iRow) {
    2560                                    double newValue = work[number] - model->dualIn();
    2561                                    if (!newValue)
    2562                                         newValue = 1.0e-100;
    2563                                    work[number] = newValue;
    2564                                    possiblePivotKey_ = -1;
    2565                               }
    2566                               which[number++] = iRow;
    2567                          } else {
    2568                               // reset
    2569                               backward_[iColumn] = iSet;
    2570                          }
    2571                          iColumn = next_[iColumn];
    2572                     }
    2573                     toIndex_[iSet] = -1;
    2574                }
    2575                if (possiblePivotKey_ >= 0) {
    2576                     work[number] = -model->dualIn();
    2577                     which[number++] = possiblePivotKey_;
    2578                     possiblePivotKey_ = -1;
    2579                }
    2580                fromIndex_[0] = -1;
    2581                array->setNumElements(number);
    2582           }
    2583      }
    2584      break;
    2585      }
     2516        }
     2517      }
     2518      // now do rest and clean up
     2519      for (i = 0; i < numberChanged; i++) {
     2520        int iSet = fromIndex_[i];
     2521        int key = keyVariable_[iSet];
     2522        int iColumn = next_[key];
     2523        double change = changeCost_[i];
     2524        while (iColumn >= 0) {
     2525          if (backward_[iColumn] >= 0) {
     2526            int iRow = backToPivotRow_[iColumn];
     2527            assert(iRow >= 0);
     2528            work[number] = change;
     2529            if (possiblePivotKey_ == iRow) {
     2530              double newValue = work[number] - model->dualIn();
     2531              if (!newValue)
     2532                newValue = 1.0e-100;
     2533              work[number] = newValue;
     2534              possiblePivotKey_ = -1;
     2535            }
     2536            which[number++] = iRow;
     2537          } else {
     2538            // reset
     2539            backward_[iColumn] = iSet;
     2540          }
     2541          iColumn = next_[iColumn];
     2542        }
     2543        toIndex_[iSet] = -1;
     2544      }
     2545      if (possiblePivotKey_ >= 0) {
     2546        work[number] = -model->dualIn();
     2547        which[number++] = possiblePivotKey_;
     2548        possiblePivotKey_ = -1;
     2549      }
     2550      fromIndex_[0] = -1;
     2551      array->setNumElements(number);
     2552    }
     2553  } break;
     2554  }
    25862555}
    25872556// This is local to Gub to allow synchronization when status is good
    2588 int
    2589 ClpGubMatrix::synchronize(ClpSimplex *, int)
     2557int ClpGubMatrix::synchronize(ClpSimplex *, int)
    25902558{
    2591      return 0;
     2559  return 0;
    25922560}
    25932561/*
     
    25962564     Remember to update here when settled down
    25972565*/
    2598 int
    2599 ClpGubMatrix::generalExpanded(ClpSimplex * model, int mode, int &number)
     2566int ClpGubMatrix::generalExpanded(ClpSimplex *model, int mode, int &number)
    26002567{
    2601      int returnCode = 0;
    2602      int numberColumns = model->numberColumns();
    2603      switch (mode) {
    2604           // Fill in pivotVariable but not for key variables
    2605      case 0: {
    2606           if (!next_ ) {
    2607                // do ordering
    2608                assert (!rhsOffset_);
    2609                // create and do gub crash
    2610                useEffectiveRhs(model, false);
    2611           }
    2612           int i;
    2613           int numberBasic = number;
    2614           // Use different array so can build from true pivotVariable_
    2615           //int * pivotVariable = model->pivotVariable();
    2616           int * pivotVariable = model->rowArray(0)->getIndices();
    2617           for (i = 0; i < numberColumns; i++) {
    2618                if (model->getColumnStatus(i) == ClpSimplex::basic) {
    2619                     int iSet = backward_[i];
    2620                     if (iSet < 0 || i != keyVariable_[iSet])
    2621                          pivotVariable[numberBasic++] = i;
    2622                }
    2623           }
    2624           number = numberBasic;
    2625           if (model->numberIterations())
    2626                assert (number == model->numberRows());
    2627      }
    2628      break;
    2629      // Make all key variables basic
    2630      case 1: {
    2631           int i;
    2632           for (i = 0; i < numberSets_; i++) {
    2633                int iColumn = keyVariable_[i];
    2634                if (iColumn < numberColumns)
    2635                     model->setColumnStatus(iColumn, ClpSimplex::basic);
    2636           }
    2637      }
    2638      break;
    2639      // Do initial extra rows + maximum basic
    2640      case 2: {
    2641           returnCode = getNumRows() + 1;
    2642           number = model->numberRows() + numberSets_;
    2643      }
    2644      break;
    2645      // Before normal replaceColumn
    2646      case 3: {
    2647           int sequenceIn = model->sequenceIn();
    2648           int sequenceOut = model->sequenceOut();
    2649           int numberColumns = model->numberColumns();
    2650           int numberRows = model->numberRows();
    2651           int pivotRow = model->pivotRow();
    2652           if (gubSlackIn_ >= 0)
    2653                assert (sequenceIn > numberRows + numberColumns);
    2654           if (sequenceIn == sequenceOut)
    2655                return -1;
    2656           int iSetIn = -1;
    2657           int iSetOut = -1;
    2658           if (sequenceOut < numberColumns) {
    2659                iSetOut = backward_[sequenceOut];
    2660           } else if (sequenceOut >= numberRows + numberColumns) {
    2661                assert (pivotRow >= numberRows);
    2662                int iExtra = pivotRow - numberRows;
    2663                assert (iExtra >= 0);
    2664                if (iSetOut < 0)
    2665                     iSetOut = fromIndex_[iExtra];
    2666                else
    2667                     assert(iSetOut == fromIndex_[iExtra]);
    2668           }
    2669           if (sequenceIn < numberColumns) {
    2670                iSetIn = backward_[sequenceIn];
    2671           } else if (gubSlackIn_ >= 0) {
    2672                iSetIn = gubSlackIn_;
    2673           }
    2674           possiblePivotKey_ = -1;
    2675           number = 0; // say do ordinary
    2676           int * pivotVariable = model->pivotVariable();
    2677           if (pivotRow >= numberRows) {
    2678                int iExtra = pivotRow - numberRows;
    2679                //const int * length = matrix_->getVectorLengths();
     2568  int returnCode = 0;
     2569  int numberColumns = model->numberColumns();
     2570  switch (mode) {
     2571    // Fill in pivotVariable but not for key variables
     2572  case 0: {
     2573    if (!next_) {
     2574      // do ordering
     2575      assert(!rhsOffset_);
     2576      // create and do gub crash
     2577      useEffectiveRhs(model, false);
     2578    }
     2579    int i;
     2580    int numberBasic = number;
     2581    // Use different array so can build from true pivotVariable_
     2582    //int * pivotVariable = model->pivotVariable();
     2583    int *pivotVariable = model->rowArray(0)->getIndices();
     2584    for (i = 0; i < numberColumns; i++) {
     2585      if (model->getColumnStatus(i) == ClpSimplex::basic) {
     2586        int iSet = backward_[i];
     2587        if (iSet < 0 || i != keyVariable_[iSet])
     2588          pivotVariable[numberBasic++] = i;
     2589      }
     2590    }
     2591    number = numberBasic;
     2592    if (model->numberIterations())
     2593      assert(number == model->numberRows());
     2594  } break;
     2595  // Make all key variables basic
     2596  case 1: {
     2597    int i;
     2598    for (i = 0; i < numberSets_; i++) {
     2599      int iColumn = keyVariable_[i];
     2600      if (iColumn < numberColumns)
     2601        model->setColumnStatus(iColumn, ClpSimplex::basic);
     2602    }
     2603  } break;
     2604  // Do initial extra rows + maximum basic
     2605  case 2: {
     2606    returnCode = getNumRows() + 1;
     2607    number = model->numberRows() + numberSets_;
     2608  } break;
     2609  // Before normal replaceColumn
     2610  case 3: {
     2611    int sequenceIn = model->sequenceIn();
     2612    int sequenceOut = model->sequenceOut();
     2613    int numberColumns = model->numberColumns();
     2614    int numberRows = model->numberRows();
     2615    int pivotRow = model->pivotRow();
     2616    if (gubSlackIn_ >= 0)
     2617      assert(sequenceIn > numberRows + numberColumns);
     2618    if (sequenceIn == sequenceOut)
     2619      return -1;
     2620    int iSetIn = -1;
     2621    int iSetOut = -1;
     2622    if (sequenceOut < numberColumns) {
     2623      iSetOut = backward_[sequenceOut];
     2624    } else if (sequenceOut >= numberRows + numberColumns) {
     2625      assert(pivotRow >= numberRows);
     2626      int iExtra = pivotRow - numberRows;
     2627      assert(iExtra >= 0);
     2628      if (iSetOut < 0)
     2629        iSetOut = fromIndex_[iExtra];
     2630      else
     2631        assert(iSetOut == fromIndex_[iExtra]);
     2632    }
     2633    if (sequenceIn < numberColumns) {
     2634      iSetIn = backward_[sequenceIn];
     2635    } else if (gubSlackIn_ >= 0) {
     2636      iSetIn = gubSlackIn_;
     2637    }
     2638    possiblePivotKey_ = -1;
     2639    number = 0; // say do ordinary
     2640    int *pivotVariable = model->pivotVariable();
     2641    if (pivotRow >= numberRows) {
     2642      int iExtra = pivotRow - numberRows;
     2643      //const int * length = matrix_->getVectorLengths();
    26802644
    2681                assert (sequenceOut >= numberRows + numberColumns ||
    2682                        sequenceOut == keyVariable_[iSetOut]);
    2683                int incomingColumn = sequenceIn; // to be used in updates
    2684                if (iSetIn != iSetOut) {
    2685                     // We need to find a possible pivot for incoming
    2686                     // look through rowArray_[1]
    2687                     int n = model->rowArray(1)->getNumElements();
    2688                     int * which = model->rowArray(1)->getIndices();
    2689                     double * array = model->rowArray(1)->denseVector();
    2690                     double bestAlpha = 1.0e-5;
    2691                     //int shortest=numberRows+1;
    2692                     for (int i = 0; i < n; i++) {
    2693                          int iRow = which[i];
    2694                          int iPivot = pivotVariable[iRow];
    2695                          if (iPivot < numberColumns && backward_[iPivot] == iSetOut) {
    2696                               if (fabs(array[i]) > fabs(bestAlpha)) {
    2697                                    bestAlpha = array[i];
    2698                                    possiblePivotKey_ = iRow;
    2699                               }
    2700                          }
    2701                     }
    2702                     assert (possiblePivotKey_ >= 0); // could set returnCode=4
    2703                     number = 1;
    2704                     if (sequenceIn >= numberRows + numberColumns) {
    2705                          number = 3;
    2706                          // need swap as gub slack in and must become key
    2707                          // is this best way
    2708                          int key = keyVariable_[iSetIn];
    2709                          assert (key < numberColumns);
    2710                          // check other basic
    2711                          int iColumn = next_[key];
    2712                          // set new key to be used by unpack
    2713                          keyVariable_[iSetIn] = iSetIn + numberColumns;
    2714                          // change cost in changeCost
    2715                          {
    2716                               int iLook = 0;
    2717                               int iSet = fromIndex_[0];
    2718                               while (iSet >= 0) {
    2719                                    if (iSet == iSetIn) {
    2720                                         changeCost_[iLook] = 0.0;
    2721                                         break;
    2722                                    }
    2723                                    iSet = fromIndex_[++iLook];
    2724                               }
    2725                          }
    2726                          while (iColumn >= 0) {
    2727                               if (iColumn != sequenceOut) {
    2728                                    // need partial ftran and skip accuracy check in replaceColumn
     2645      assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
     2646      int incomingColumn = sequenceIn; // to be used in updates
     2647      if (iSetIn != iSetOut) {
     2648        // We need to find a possible pivot for incoming
     2649        // look through rowArray_[1]
     2650        int n = model->rowArray(1)->getNumElements();
     2651        int *which = model->rowArray(1)->getIndices();
     2652        double *array = model->rowArray(1)->denseVector();
     2653        double bestAlpha = 1.0e-5;
     2654        //int shortest=numberRows+1;
     2655        for (int i = 0; i < n; i++) {
     2656          int iRow = which[i];
     2657          int iPivot = pivotVariable[iRow];
     2658          if (iPivot < numberColumns && backward_[iPivot] == iSetOut) {
     2659            if (fabs(array[i]) > fabs(bestAlpha)) {
     2660              bestAlpha = array[i];
     2661              possiblePivotKey_ = iRow;
     2662            }
     2663          }
     2664        }
     2665        assert(possiblePivotKey_ >= 0); // could set returnCode=4
     2666        number = 1;
     2667        if (sequenceIn >= numberRows + numberColumns) {
     2668          number = 3;
     2669          // need swap as gub slack in and must become key
     2670          // is this best way
     2671          int key = keyVariable_[iSetIn];
     2672          assert(key < numberColumns);
     2673          // check other basic
     2674          int iColumn = next_[key];
     2675          // set new key to be used by unpack
     2676          keyVariable_[iSetIn] = iSetIn + numberColumns;
     2677          // change cost in changeCost
     2678          {
     2679            int iLook = 0;
     2680            int iSet = fromIndex_[0];
     2681            while (iSet >= 0) {
     2682        &nbs